StateManager


   // StateManager.h

#
ifndef __ScoreKeeper_StateManager_h__ #define __ScoreKeeper_StateManager_h__ #include "Track.h" #include "StackManager.h" #include "Callbacks.h" #ifdef __ENABLE_STARMC__   #include "TGiant3.h"   #include "St_geant_Maker/St_geant_Maker.h"   #include "StarMCCallbacks.h" #endif //#define __ENABLE_G4__ #ifdef __ENABLE_G4__   #include "G4Event.hh"   #include "G4Track.hh"   #include "G4EventManager.hh"   #include "G4Callbacks.h" #endif #include <iostream> // Use compensated sum for accumulated quantities. The second element of accumulated // arrays will be the compensation term // https://en.wikipedia.org/wiki/Kahan_summation_algorithm namespace ScoreKeeper { struct State_t { // event int idevent; // track identity int pdgid; // track enumeration, to be defined by stacking mechanism int idtrack; int idparent; // state flags int flags; // current state of track on geometry int nstep; double step; double astep[2]; // accumulated step size double x, y, z, px, py, pz; double tof; // time of flight (from prod vtx to here) double etotal; // total energy double atof[2]; // time of flight since collion (accumlated over parent tracks) // physical response of detector to passage of particle double destep; double adestep[2]; // accumulated destep }; template<typename MC> class StateManager {   public: typedef MC MonteCarlo_t; // Update the current tracking state and return a reference // to the state data State_t& Update() { // Previous state _prev = _state; // Obtain state flags _state.flags = getStateFlags ( _mc ); _state.idevent = getIdEvent ( _mc ); // _state.idtrack = getIdTrack ( _mc ); _state.pdgid = getTrackPdgId ( _mc ); // _state.idparent = getIdParent ( _mc ); _state.nstep = getStepNumber ( _mc ); _state.step = getStepSize ( _mc ); _state.x = getx( _mc ); _state.y = gety( _mc ); _state.z = getz( _mc ); _state.tof = gett( _mc ); _state.px = getpx( _mc ); _state.py = getpy( _mc ); _state.pz = getpz( _mc ); _state.etotal = getTotalEnergy( _mc ); _state.destep = getEnergyDeposit( _mc ); kahanSum( _state.astep, _state.step ); kahanSum( _state.atof, _state.tof ); kahanSum( _state.adestep, _state.destep ); return _state; } // Signal a new event, clear the state and increment event id. // (Event ID may be overwritten by underlying MC) State_t& NewEvent() { State_t old = _state; State_t s = {}; _state = s; _state.idevent = old.idevent + 1; return _state; } // Signal a new track, clear all and increment track id // (Track ID may be overwritten by underlying MC) State_t& NewTrack() { State_t old = _state; State_t s = {}; _state = s; _state.idtrack = old.idtrack + 1; return _state; }; State_t& Current(){ return _state; } State_t& Previous(){ return _prev; }   private:   protected: State_t _state; // current state State_t _prev; // previous state MonteCarlo_t _mc; // MC type // Performs kahan compensated summation void kahanSum( double* x, const double in ) { double& sum = x[0]; double& c = x[1]; double y = in - c; // On first call, c = 0 double t = sum + y; // Alas, sum is big, y small, so low-order digits of y are lost. c = (t - sum) - y; // (t - sum) cancels the high-order part of y; subtracting y recovers negative (low part of y) sum = t; }; }; std::ostream& operator << (std::ostream &out, const State_t& s ) { out << "-- mc state --------------------------------------------------------" << std::endl << "event = " << s.idevent << std::endl << "track = " << s.idtrack << std::endl << "pdgid = " << s.pdgid << std::endl << "flags = " << s.flags << std::endl << "astep = " << s.astep[0] << std::endl << "atof = " << s.atof[0] << std::endl << "adest = " << s.adestep[0] << std::endl; return out; } }; #endif // G4Callbacks.h #ifndef __ScoreKeeper_Geant4Callbacks_h__ #define __ScoreKeeper_Geant4Callbacks_h__ #include "Callbacks.h" #include <iostream> namespace ScoreKeeper { #ifdef __ENABLEG4__ template<> int getStateFlags<Geant4>( Geant4 ){ return 0; } template<> int getIdEvent <Geant4>( Geant4 ){ int id = 0; auto* mgr = G4EventManager::GetEventManager(); auto* evt = mgr->GetConstCurrentEvent(); id = evt->GetEventID(); return id; } #endif }; #endif