ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EicFRichSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EicFRichSteppingAction.cc
1 //____________________________________________________________________________..
2 //
3 // This is a working template for the Stepping Action which needs to be implemented
4 // for active detectors. Most of the code is error handling and access to the G4 objects
5 // and our data structures. It does not need any adjustment. The only thing you need to
6 // do is to add the properties of the G4Hits you want to save for later analysis
7 // This needs to be done in 2 places, G4Hits are generated when a G4 track enters a new
8 // volume (or is created). Here you give it an initial value. When the G4 track leaves
9 // the volume the final value needs to be set.
10 // The places to do this is marked by //implement your own here//
11 //
12 // As guidance you can look at the total (integrated over all steps in a volume) energy
13 // deposit which should always be saved.
14 // Additionally the total ionization energy is saved - this can be removed if you are not
15 // interested in this. Naturally you may want remove these comments in your version
16 //
17 //____________________________________________________________________________..
18 
19 #include "EicFRichSteppingAction.h"
20 
21 #include "EicFRichDetector.h"
22 
23 #include <phparameter/PHParameters.h>
24 
26 
27 #include <g4main/PHG4Hit.h>
29 #include <g4main/PHG4Hitv1.h>
30 #include <g4main/PHG4Shower.h>
33 
34 #include <phool/getClass.h>
35 
36 #include <TSystem.h>
37 
38 #include <Geant4/G4ParticleDefinition.hh>
39 #include <Geant4/G4ReferenceCountedHandle.hh>
40 #include <Geant4/G4Step.hh>
41 #include <Geant4/G4StepPoint.hh>
42 #include <Geant4/G4StepStatus.hh>
43 #include <Geant4/G4String.hh>
44 #include <Geant4/G4SystemOfUnits.hh>
45 #include <Geant4/G4ThreeVector.hh>
46 #include <Geant4/G4TouchableHandle.hh>
47 #include <Geant4/G4Track.hh>
48 #include <Geant4/G4TrackStatus.hh>
49 #include <Geant4/G4Types.hh>
50 #include <Geant4/G4VPhysicalVolume.hh>
51 #include <Geant4/G4VTouchable.hh>
52 #include <Geant4/G4VUserTrackInformation.hh>
53 
54 #include <cmath>
55 #include <iostream>
56 #include <string>
57 
58 class PHCompositeNode;
59 
60 using namespace std;
61 
62 //____________________________________________________________________________..
64  : PHG4SteppingAction(detector->GetName())
65  , m_Detector(detector)
66  , m_Params(parameters)
67  , m_HitContainer(nullptr)
68  , m_Hit(nullptr)
69  , m_SaveHitContainer(nullptr)
70  , m_SaveVolPre(nullptr)
71  , m_SaveVolPost(nullptr)
72  , m_SaveTrackId(-1)
73  , m_SavePreStepStatus(-1)
74  , m_SavePostStepStatus(-1)
75  , m_ActiveFlag(m_Params->get_int_param("active"))
76  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
77  , m_EdepSum(0)
78  , m_EionSum(0)
79 {
80 }
81 
82 //____________________________________________________________________________..
84 {
85  // if the last hit was a zero energie deposit hit, it is just reset
86  // and the memory is still allocated, so we need to delete it here
87  // if the last hit was saved, hit is a nullptr pointer which are
88  // legal to delete (it results in a no operation)
89  delete m_Hit;
90 }
91 
92 //____________________________________________________________________________..
93 // This is the implementation of the G4 UserSteppingAction
94 bool EicFRichSteppingAction::UserSteppingAction(const G4Step *aStep,bool was_used)
95 {
98  // get volume of the current step
100  // IsInDetector(volume) returns
101  // == 0 outside of detector
102  // > 0 for hits in active volume
103  // < 0 for hits in passive material
104  int whichactive = m_Detector->IsInDetector(volume);
105  if (!whichactive)
106  {
107  return false;
108  }
109  // collect energy and track length step by step
110  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
111  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
112  const G4Track *aTrack = aStep->GetTrack();
113  // if this detector stops everything, just put all kinetic energy into edep
114  if (m_BlackHoleFlag)
115  {
116  edep = aTrack->GetKineticEnergy() / GeV;
117  G4Track *killtrack = const_cast<G4Track *>(aTrack);
118  killtrack->SetTrackStatus(fStopAndKill);
119  }
120  // we use here only one detector in this simple example
121  // if you deal with multiple detectors in this stepping action
122  // the detector id can be used to distinguish between them
123  // hits can easily be analyzed later according to their detector id
124  int detector_id = 0; // we use here only one detector in this simple example
125  bool geantino = false;
126  // the check for the pdg code speeds things up, I do not want to make
127  // an expensive string compare for every track when we know
128  // geantino or chargedgeantino has pid=0
129  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
130  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") !=
131  string::npos) // this also accounts for "chargedgeantino"
132  {
133  geantino = true;
134  }
135  G4StepPoint *prePoint = aStep->GetPreStepPoint();
136  G4StepPoint *postPoint = aStep->GetPostStepPoint();
137 
138 // Here we have to decide if we need to create a new hit. Normally this should
139 // only be neccessary if a G4 Track enters a new volume or is freshly created
140 // For this we look at the step status of the prePoint (beginning of the G4 Step).
141 // This should be either fGeomBoundary (G4 Track crosses into volume) or
142 // fUndefined (G4 Track newly created)
143 // Sadly over the years with different G4 versions we have observed cases where
144 // G4 produces "impossible hits" which we try to catch here
145 // These errors were always rare and it is not clear if they still exist but we
146 // still check for them for safety. We can reproduce G4 runs identically (if given
147 // the sequence of random number seeds you find in the log), the printouts help
148 // us giving the G4 support information about those failures
149 //
150  switch (prePoint->GetStepStatus())
151  {
152  case fPostStepDoItProc:
154  {
155  // this is the okay case, fPostStepDoItProc called in a volume, not first thing inside
156  // a new volume, just proceed here
157  break;
158  }
159  else
160  {
161  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
162  // this is still with us
163  cout << GetName() << ": New Hit for " << endl;
164  cout << "prestep status: "
166  << ", poststep status: "
168  << ", last pre step status: "
170  << ", last post step status: "
172  cout << "last track: " << m_SaveTrackId
173  << ", current trackid: " << aTrack->GetTrackID() << endl;
174  cout << "phys pre vol: " << volume->GetName()
175  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
176  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
177  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
178  }
179 // These are the normal cases
180  case fGeomBoundary:
181  case fUndefined:
182  if (!m_Hit)
183  {
184  m_Hit = new PHG4Hitv1();
185  }
186  m_Hit->set_layer(detector_id);
187  // here we set the entrance values in cm
188  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
189  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
190  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
191  // time in ns
192  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
193  // set the track ID
194  m_Hit->set_trkid(aTrack->GetTrackID());
195  m_SaveTrackId = aTrack->GetTrackID();
196  // set the initial energy deposit
197  m_EdepSum = 0;
198  // implement your own here://
199  // add the properties you are interested in via set_XXX methods
200  // you can find existing set methods in $OFFLINE_MAIN/include/g4main/PHG4Hit.h
201  // this is initialization of your value. This is not needed you can just set the final
202  // value at the last step in this volume later one
203  if (whichactive > 0)
204  {
205  m_EionSum = 0; // assuming the ionization energy is only needed for active
206  // volumes (scintillators)
207  m_Hit->set_eion(0);
209  }
210  else
211  {
212  cout << "implement stuff for whichactive < 0 (inactive volumes)" << endl;
213  gSystem->Exit(1);
214  }
215  // this is for the tracking of the truth info
217  {
218  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
219  {
220  m_Hit->set_trkid(pp->GetUserTrackId());
221  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
222  }
223  }
224  break;
225  default:
226  break;
227  }
228 
229  // This section is called for every step
230  // some sanity checks for inconsistencies (aka bugs) we have seen over the years
231  // check if this hit was created, if not print out last post step status
232  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
233  {
234  cout << GetName() << ": hit was not created" << endl;
235  cout << "prestep status: "
237  << ", poststep status: "
239  << ", last pre step status: "
241  << ", last post step status: "
243  cout << "last track: " << m_SaveTrackId
244  << ", current trackid: " << aTrack->GetTrackID() << endl;
245  cout << "phys pre vol: " << volume->GetName()
246  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
247  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
248  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
249  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
250  gSystem->Exit(1);
251  }
252  // check if track id matches the initial one when the hit was created
253  if (aTrack->GetTrackID() != m_SaveTrackId)
254  {
255  cout << GetName() << ": hits do not belong to the same track" << endl;
256  cout << "saved track: " << m_SaveTrackId
257  << ", current trackid: " << aTrack->GetTrackID()
258  << ", prestep status: " << prePoint->GetStepStatus()
259  << ", previous post step status: " << m_SavePostStepStatus << endl;
260  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
261  gSystem->Exit(1);
262  }
263 
264 // We need to cache a few things from one step to the next
265 // to identify impossible hits and subsequent debugging printout
266  m_SavePreStepStatus = prePoint->GetStepStatus();
267  m_SavePostStepStatus = postPoint->GetStepStatus();
269  m_SaveVolPost = touchpost->GetVolume();
270  // here we just update the exit values, it will be overwritten
271  // for every step until we leave the volume or the particle
272  // ceases to exist
273  // sum up the energy to get total deposited
274  m_EdepSum += edep;
275  if (whichactive > 0)
276  {
277  m_EionSum += eion;
278  }
279  // if any of these conditions is true this is the last step in
280  // this volume and we need to save the hit
281  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
282  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
283  // (happens when your detector goes outside world volume)
284  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
285  // aTrack->GetTrackStatus() == fStopAndKill is also set)
286  // aTrack->GetTrackStatus() == fStopAndKill: track ends
287  if (postPoint->GetStepStatus() == fGeomBoundary ||
288  postPoint->GetStepStatus() == fWorldBoundary ||
289  postPoint->GetStepStatus() == fAtRestDoItProc ||
290  aTrack->GetTrackStatus() == fStopAndKill)
291  {
292  // save only hits with energy deposit (or geantino)
293  if (m_EdepSum > 0 || geantino)
294  {
295  // update values at exit coordinates and set keep flag
296  // of track to keep
297  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
298  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
299  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
300  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
302  {
303  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
304  {
305  pp->SetKeep(1); // we want to keep the track
306  }
307  }
308  if (geantino)
309  {
310  //implement your own here://
311  // if you want to do something special for geantinos (normally you do not)
312  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way
313  // geantinos survive the g4hit compression
314  if (whichactive > 0)
315  {
316  m_Hit->set_eion(-1);
317  }
318  }
319  else
320  {
322  }
323  //implement your own here://
324  // what you set here will be saved in the output
325  if (whichactive > 0)
326  {
328  }
329  m_SaveHitContainer->AddHit(detector_id, m_Hit);
330  // ownership has been transferred to container, set to null
331  // so we will create a new hit for the next track
332  m_Hit = nullptr;
333  }
334  else
335  {
336  // if this hit has no energy deposit, just reset it for reuse
337  // this means we have to delete it in the dtor. If this was
338  // the last hit we processed the memory is still allocated
339  m_Hit->Reset();
340  }
341  }
342  // return true to indicate the hit was used
343  return true;
344 }
345 
346 //____________________________________________________________________________..
348 {
349  string hitnodename = "G4HIT_" + m_Detector->GetName();
350  // now look for the map and grab a pointer to it.
351  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
352  // if we do not find the node we need to make it.
353  if (!m_HitContainer)
354  {
355  std::cout << "EicFRichSteppingAction::SetTopNode - unable to find "
356  << hitnodename << std::endl;
357  }
358 }