Belle II Software  release-05-02-19
TrackingAction.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Guofu Cao, Andreas Moll, Marko Staric *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <simulation/kernel/TrackingAction.h>
12 #include <simulation/kernel/UserInfo.h>
13 #include <simulation/kernel/SensitiveDetectorBase.h>
14 #include <framework/logging/Logger.h>
15 #include <framework/gearbox/Unit.h>
16 
17 #include <G4TrackingManager.hh>
18 #include <G4Track.hh>
19 #include <G4ParticleDefinition.hh>
20 #include <G4ParticleTypes.hh>
21 #include <G4EmProcessSubType.hh>
22 #include <G4DecayProcessType.hh>
23 #include <G4Event.hh>
24 #include <G4ParticleTable.hh>
25 
26 using namespace Belle2;
27 using namespace Belle2::Simulation;
28 
30  G4UserTrackingAction(), m_mcParticleGraph(mcParticleGraph),
31  m_ignoreOpticalPhotons(false),
32  m_ignoreSecondaries(false), m_secondariesEnergyCut(0.0),
33  m_ignoreBremsstrahlungPhotons(false), m_bremsstrahlungPhotonsEnergyCut(0.0),
34  m_ignorePairConversions(false), m_pairConversionsEnergyCut(0.0),
35  m_storeTrajectories(false), m_distanceTolerance(0),
36  m_storeMCTrajectories(), m_relMCTrajectories(StoreArray<MCParticle>(), m_storeMCTrajectories)
37 {
38 }
39 
41 {
42 }
43 
44 void TrackingAction::setStoreTrajectories(int store, double distanceTolerance)
45 {
46  m_storeTrajectories = store;
47  m_distanceTolerance = distanceTolerance;
48  if (store) {
49  // registration of store arrays and relations
50  m_storeMCTrajectories.registerInDataStore();
51  StoreArray<MCParticle> mcParticles;
53 
54  // additional registration of MCParticle relation
55  // (note: m_relMCTrajectories is already defined by TrackingAction::TrackingAction)
57  }
58 }
59 
60 void TrackingAction::PreUserTrackingAction(const G4Track* track)
61 {
62  //We only want to do the following for new tracks, not for suspended and reactivated ones"
63  if (track->GetCurrentStepNumber() > 0) return;
64 
65  const G4DynamicParticle* dynamicParticle = track->GetDynamicParticle();
66 
67  try {
68  //Check if the dynamic particle has a primary particle attached.
69  //If the answer is yes, the UserInfo of the primary particle is recycled as the UserInfo of the track.
70  const bool isPrimary = dynamicParticle->GetPrimaryParticle() != nullptr;
71  bool neutral_llp = false;
72  if (isPrimary) {
73  const G4PrimaryParticle* primaryParticle = dynamicParticle->GetPrimaryParticle();
74  if (primaryParticle->GetUserInformation() != nullptr) {
75  const_cast<G4Track*>(track)->SetUserInformation(new TrackInfo(ParticleInfo::getInfo(*primaryParticle)));
76  //check for neutral long-lived primary particle
77  if (primaryParticle->GetParticleDefinition() == G4ParticleTable::GetParticleTable()->FindParticle("LongLivedNeutralParticle")) {
78  neutral_llp = true;
79  }
80  } else {
81  B2WARNING(track->GetDefinition()->GetPDGEncoding() << " has no MCParticle user information !");
82  }
83  }
84 
85  //Get particle of current track
86  MCParticleGraph::GraphParticle& currParticle = TrackInfo::getInfo(*track);
87 
88  G4ThreeVector dpMom = dynamicParticle->GetMomentum() / CLHEP::MeV * Unit::MeV;
89  currParticle.setTrackID(track->GetTrackID());
90  // If the particle is a primary particle with simulated parents update the
91  // production vertex and time to the real values from Geant4. As convention
92  // we set NaN as production time in the MCParticleGenerator if the particle
93  // wasn't directly placed with a primary vertex in Geant4.
94  if (!isPrimary or std::isnan(currParticle.getProductionTime())) {
95  G4ThreeVector trVtxPos = track->GetVertexPosition() / CLHEP::mm * Unit::mm;
96  currParticle.setProductionTime(track->GetGlobalTime()); // Time does not need a conversion factor
97  currParticle.setProductionVertex(trVtxPos.x(), trVtxPos.y(), trVtxPos.z());
98  }
99 
100  //Set the Values of the particle which are already known
101  if (!neutral_llp) {
102  currParticle.setPDG(dynamicParticle->GetPDGcode());
103  currParticle.setMass(dynamicParticle->GetMass() / CLHEP::MeV * Unit::MeV);
104  currParticle.setEnergy(dynamicParticle->GetTotalEnergy() / CLHEP::MeV * Unit::MeV);
105  currParticle.setMomentum(dpMom.x(), dpMom.y(), dpMom.z());
106  }
107 
108  //Primary or secondary particle?
109  if (isPrimary) {
110  //Primary particle
111  currParticle.setSecondaryPhysicsProcess(0);
112  currParticle.setIgnore(false); //Store the generator info in the MCParticles block.
113 
114  } else if (track->GetCreatorProcess() != nullptr) {
115  //Secondary particle
116  const int& processSubType = track->GetCreatorProcess()->GetProcessSubType();
117  currParticle.setSecondaryPhysicsProcess(processSubType); //Store the physics process type.
118 
119  //Decay-in-flight
120  if (processSubType >= static_cast<int>(DECAY) && processSubType <= static_cast<int>(DECAY_External))
121  currParticle.setIgnore(false); //Store the generator info in the MCParticles block.
122 
123  } else {
124  //Unknown origin. This could be a bug originated from Geant4.
125  currParticle.setSecondaryPhysicsProcess(-1);
126  }
127 
128  //Either we store all the trajectories
129  if (m_storeTrajectories > 2 ||
130  //Or only the primary ones
131  (m_storeTrajectories == 1 && isPrimary) ||
132  //Or all except optical photons
133  (m_storeTrajectories == 2 && track->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition())) {
134  TrackInfo& info = dynamic_cast<TrackInfo&>(*track->GetUserInformation());
135  m_relMCTrajectories.add(track->GetTrackID(), m_storeMCTrajectories.getEntries());
136  info.setTrajectory(m_storeMCTrajectories.appendNew());
137  }
138 
139  } catch (CouldNotFindUserInfo& exc) {
140  B2FATAL(exc.what());
141  }
142 }
143 
144 
145 void TrackingAction::PostUserTrackingAction(const G4Track* track)
146 {
147  G4StepPoint* postStep = track->GetStep()->GetPostStepPoint();
148 
149  // Get particle of the current track
150  try {
151  MCParticleGraph::GraphParticle& currParticle = TrackInfo::getInfo(*track);
152 
153  // Add particle and decay Information to all the secondaries
154  for (G4Track* daughterTrack : *fpTrackingManager->GimmeSecondaries()) {
155 
156  // Add the particle to the particle graph and as UserInfo to the track
157  // if it is a secondary particle created by Geant4.
158  if (daughterTrack->GetDynamicParticle()->GetPrimaryParticle() == NULL &&
159  daughterTrack->GetUserInformation() == NULL) {
161  const_cast<G4Track*>(daughterTrack)->SetUserInformation(new TrackInfo(daughterParticle));
162 
163  currParticle.decaysInto(daughterParticle); //Add the decay history
164 
165  // Optical photons and secondaries: steering of output to MCParticles
166  if (daughterTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) {
167 
168  // Optical photons
169  if (m_ignoreOpticalPhotons) daughterParticle.setIgnore();
170  // to apply quantum efficiency only once, if optical photon is a daugher of optical photon
171  if (track->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) {
172  TrackInfo* currInfo = dynamic_cast<TrackInfo*>(track->GetUserInformation());
173  TrackInfo* daughterInfo = dynamic_cast<TrackInfo*>(daughterTrack->GetUserInformation());
174  daughterInfo->setStatus(currInfo->getStatus());
175  daughterInfo->setFraction(currInfo->getFraction());
176  }
177 
178  } else if (daughterTrack->GetCreatorProcess()->GetProcessSubType() == fBremsstrahlung) {
179 
180  // Bremsstrahlung photons
181  // Do not store the generator info in the final MCParticles block
182  // if the ignore flag is set or its energy is too low in [MeV].
183  if (m_ignoreBremsstrahlungPhotons || daughterTrack->GetKineticEnergy() < m_bremsstrahlungPhotonsEnergyCut)
184  daughterParticle.setIgnore();
185 
186  } else if (daughterTrack->GetCreatorProcess()->GetProcessSubType() == fGammaConversion) {
187 
188  // e+ or e- created by gamma conversion to pairs
189  // Do not store the generator info in the final MCParticles block
190  // if the ignore flag is set or kinetic energy is too low in [MeV].
191  if (m_ignorePairConversions || daughterTrack->GetKineticEnergy() < m_pairConversionsEnergyCut)
192  daughterParticle.setIgnore();
193 
194  } else {
195 
196  // Do not store the generator info in the final MCParticles block
197  // if the ignore flag is set or its kinetic energy is too low in [MeV].
198  if (m_ignoreSecondaries || daughterTrack->GetKineticEnergy() < m_secondariesEnergyCut)
199  daughterParticle.setIgnore();
200 
201  //B2INFO("Secondary Physics Process: " << daughterTrack->GetCreatorProcess()->GetProcessSubType());
202  }
203 
204  }
205  }
206  //If the track is just suspended we can return here: the rest should be filled once the track is done
207  if (track->GetTrackStatus() == fSuspend) return;
208 
209  //Check if particle left the detector/simulation volume.
210  if (postStep->GetStepStatus() == fWorldBoundary) {
212  }
213 
214  //Check if particle was stopped in the detector/simulation volume
215  if (track->GetKineticEnergy() <= 0.0) {
217  }
218 
219  //Set the values for the particle
220  G4ThreeVector decVtx = postStep->GetPosition() / CLHEP::mm * Unit::mm;
221  currParticle.setDecayVertex(decVtx.x(), decVtx.y(), decVtx.z());
222  currParticle.setDecayTime(postStep->GetGlobalTime()); // Time does not need a conversion factor
223  currParticle.setValidVertex(true);
224 
225  //Check if we can remove some points from the trajectory
226  if (m_storeTrajectories) {
227  MCParticleTrajectory* tr = dynamic_cast<TrackInfo*>(track->GetUserInformation())->getTrajectory();
228  if (tr) {
230  }
231  }
232  } catch (CouldNotFindUserInfo& exc) {
233  B2FATAL(exc.what());
234  }
235 }
Belle2::MCParticleGraph::GraphParticle::setTrackID
void setTrackID(int trackID)
Set the track ID for the particle.
Definition: MCParticleGraph.h:158
Belle2::Simulation::TrackingAction::m_storeMCTrajectories
StoreArray< MCParticleTrajectory > m_storeMCTrajectories
Store array for the Trajectories.
Definition: TrackingAction.h:140
Belle2::MCParticle::c_StoppedInDetector
@ c_StoppedInDetector
bit 3: Particle was stopped in the detector (the simulation volume).
Definition: MCParticle.h:64
Belle2::Simulation::TrackingAction::m_ignoreSecondaries
bool m_ignoreSecondaries
do not store secondaries in MCParticles
Definition: TrackingAction.h:128
Belle2::StoreArray::registerRelationTo
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:150
Belle2::RelationArray::c_negativeWeight
@ c_negativeWeight
Flip the sign of the weight to become negative if the original element got re-attributed.
Definition: RelationArray.h:89
Belle2::Simulation::TrackingAction::setStoreTrajectories
void setStoreTrajectories(int store, double distanceTolerance)
Sets the trajectory option to enable storing of the simulated particle trajectories.
Definition: TrackingAction.cc:44
Belle2::Simulation::UserInfo::getInfo
static Payload getInfo(Carrier &obj)
Static function to just return UserInformation attached to the obj of type Carrier.
Definition: UserInfo.h:110
Belle2::MCParticleGraph::GraphParticle::setIgnore
void setIgnore(bool ignore=true)
Set or remove the ignore flag.
Definition: MCParticleGraph.h:143
Belle2::Simulation::TrackingAction::m_ignoreBremsstrahlungPhotons
bool m_ignoreBremsstrahlungPhotons
do not store bremsstrahlung photons in MCParticles
Definition: TrackingAction.h:130
Belle2::Simulation::TrackingAction::TrackingAction
TrackingAction(MCParticleGraph &mcParticleGraph)
Constructor.
Definition: TrackingAction.cc:29
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2::Simulation::UserInfo
UserInfo class which is used to attach additional information to Geant4 particles and tracks.
Definition: UserInfo.h:46
Belle2::MCParticle::setDecayVertex
void setDecayVertex(const TVector3 &vertex)
Set decay vertex.
Definition: MCParticle.h:452
Belle2::Simulation::UserInfo::setStatus
void setStatus(int status)
Set status of optical photon (used for performance speed-ups)
Definition: UserInfo.h:87
Belle2::Simulation::SensitiveDetectorBase::registerMCParticleRelation
static void registerMCParticleRelation(const std::string &name, RelationArray::EConsolidationAction ignoreAction=RelationArray::c_negativeWeight)
Register an relation involving MCParticles.
Definition: SensitiveDetectorBase.cc:24
Belle2::Simulation::TrackingAction::m_bremsstrahlungPhotonsEnergyCut
double m_bremsstrahlungPhotonsEnergyCut
kinetic energy cut for stored bremsstrahlung photons [MeV]
Definition: TrackingAction.h:131
Belle2::Simulation::TrackingAction::PreUserTrackingAction
void PreUserTrackingAction(const G4Track *track)
Checks if the particle associated to the track is already in the MCParticle list.
Definition: TrackingAction.cc:60
Belle2::MCParticle::c_LeftDetector
@ c_LeftDetector
bit 2: Particle left the detector (the simulation volume).
Definition: MCParticle.h:62
Belle2::Simulation::UserInfo::getStatus
int getStatus()
Get status of optical photon (used for performance speed-ups) 0 initial 1 prescaled in StackingAction...
Definition: UserInfo.h:69
Belle2::Simulation::TrackingAction::~TrackingAction
virtual ~TrackingAction()
Destructor.
Definition: TrackingAction.cc:40
Belle2::Unit::MeV
static const double MeV
[megaelectronvolt]
Definition: Unit.h:124
Belle2::MCParticleTrajectory
Class to save the full simulated trajectory of a particle.
Definition: MCParticleTrajectory.h:32
Belle2::MCParticleTrajectory::simplify
void simplify(float distanceTolerance)
Simplify the trajectory using the Ramer-Douglas-Peuker algorithm.
Definition: MCParticleTrajectory.cc:35
Belle2::MCParticle::setProductionTime
void setProductionTime(float time)
Set production time.
Definition: MCParticle.h:392
Belle2::Simulation::TrackingAction::m_secondariesEnergyCut
double m_secondariesEnergyCut
kinetic energy cut for stored secondaries [MeV]
Definition: TrackingAction.h:129
Belle2::MCParticle::setSecondaryPhysicsProcess
void setSecondaryPhysicsProcess(int physicsProcess)
Sets the physics process type of a secondary particle.
Definition: MCParticle.h:473
Belle2::Simulation::UserInfo::setFraction
void setFraction(double fraction)
Store optical photon propagation fraction (used for performance speed-ups)
Definition: UserInfo.h:93
Belle2::Simulation::TrackingAction::m_relMCTrajectories
RelationArray m_relMCTrajectories
RelationArry for the relation between MCParticles and Trajectories.
Definition: TrackingAction.h:142
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::MCParticle::setPDG
void setPDG(int pdg)
Set PDG code of the particle.
Definition: MCParticle.h:343
Belle2::Simulation::UserInfo::getFraction
double getFraction()
Get optical photon propagation fraction (used for performance speed-ups) status=0: fraction=1 status=...
Definition: UserInfo.h:78
Belle2::Simulation::TrackingAction::m_pairConversionsEnergyCut
double m_pairConversionsEnergyCut
kinetic energy cut for stored e+ or e- from pair conversions [MeV]
Definition: TrackingAction.h:133
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::RelationArray::add
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
Definition: RelationArray.h:291
Belle2::MCParticle::setProductionVertex
void setProductionVertex(const TVector3 &vertex)
Set production vertex position.
Definition: MCParticle.h:404
Belle2::MCParticle::addStatus
void addStatus(unsigned short int bitmask)
Add bitmask to current status.
Definition: MCParticle.h:361
Belle2::Simulation::TrackingAction::m_distanceTolerance
double m_distanceTolerance
distance tolerance to merge trajectory points
Definition: TrackingAction.h:137
Belle2::Simulation::TrackingAction::m_mcParticleGraph
MCParticleGraph & m_mcParticleGraph
< Reference to the MCParticle graph which is updated by the tracking action.
Definition: TrackingAction.h:125
Belle2::Simulation::TrackingAction::PostUserTrackingAction
void PostUserTrackingAction(const G4Track *track)
Updates the data of the MCParticle associated with the Geant4 track.
Definition: TrackingAction.cc:145
Belle2::MCParticle::setEnergy
void setEnergy(float energy)
Set energy.
Definition: MCParticle.h:380
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
Belle2::Simulation::TrackingAction::m_ignorePairConversions
bool m_ignorePairConversions
do not store e+ or e- from pair conversions in MCparticles
Definition: TrackingAction.h:132
Belle2::MCParticle::setMass
void setMass(float mass)
Set particle mass.
Definition: MCParticle.h:374
Belle2::Simulation::TrackingAction::m_ignoreOpticalPhotons
bool m_ignoreOpticalPhotons
do not store optical photons in MCParticles
Definition: TrackingAction.h:127
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::MCParticleGraph::GraphParticle::decaysInto
void decaysInto(GraphParticle &daughter)
Tells the graph that this particle decays into daughter.
Definition: MCParticleGraph.h:110
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::MCParticle::setDecayTime
void setDecayTime(float time)
Set decay time.
Definition: MCParticle.h:398
Belle2::Simulation::TrackingAction::m_storeTrajectories
int m_storeTrajectories
Store trajectories for 0=none, 1=primary or 2=all particles.
Definition: TrackingAction.h:136
Belle2::MCParticle::setMomentum
void setMomentum(const TVector3 &momentum)
Set particle momentum.
Definition: MCParticle.h:425
Belle2::MCParticle::setValidVertex
void setValidVertex(bool valid)
Set indication wether vertex and time information is valid or just default.
Definition: MCParticle.h:386
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86
Belle2::MCParticle::getProductionTime
float getProductionTime() const
Return production time in ns.
Definition: MCParticle.h:170