ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4RPSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4RPSteppingAction.cc
1 //____________________________________________________________________________..
2 //
3 //
4 //____________________________________________________________________________..
5 
7 
8 #include "EICG4RPDetector.h"
9 #include "EICG4RPSubsystem.h"
10 
11 #include <phparameter/PHParameters.h>
12 
14 
15 #include <g4main/PHG4Hit.h>
17 #include <g4main/PHG4Hitv1.h>
18 #include <g4main/PHG4Shower.h>
21 
22 #include <phool/getClass.h>
23 
24 #include <TSystem.h>
25 
26 #include <Geant4/G4NavigationHistory.hh>
27 #include <Geant4/G4ParticleDefinition.hh>
28 #include <Geant4/G4ReferenceCountedHandle.hh>
29 #include <Geant4/G4Step.hh>
30 #include <Geant4/G4StepPoint.hh>
31 #include <Geant4/G4StepStatus.hh>
32 #include <Geant4/G4String.hh>
33 #include <Geant4/G4SystemOfUnits.hh>
34 #include <Geant4/G4ThreeVector.hh>
35 #include <Geant4/G4TouchableHandle.hh>
36 #include <Geant4/G4Track.hh>
37 #include <Geant4/G4TrackStatus.hh>
38 #include <Geant4/G4Types.hh>
39 #include <Geant4/G4VPhysicalVolume.hh>
40 #include <Geant4/G4VTouchable.hh>
41 #include <Geant4/G4VUserTrackInformation.hh>
42 
43 #include <cmath>
44 #include <cstdlib>
45 #include <iomanip>
46 #include <iostream>
47 #include <string>
48 
49 class PHCompositeNode;
50 
51 //____________________________________________________________________________..
53  : PHG4SteppingAction(detector->GetName())
54  , m_Subsystem(subsys)
55  , m_Detector(detector)
56  , m_Params(parameters)
57  , m_HitContainer(nullptr)
58  , m_Hit(nullptr)
59  , m_SaveShower(nullptr)
60  , m_SaveVolPre(nullptr)
61  , m_SaveVolPost(nullptr)
62  , m_SaveLightYieldFlag(m_Params->get_int_param("lightyield"))
63  , m_SaveTrackId(-1)
64  , m_SavePreStepStatus(-1)
65  , m_SavePostStepStatus(-1)
66  , m_ActiveFlag(m_Params->get_int_param("active"))
67  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
68  , m_UseG4StepsFlag(m_Params->get_int_param("use_g4steps"))
69  , m_Tmin(m_Params->get_double_param("tmin") * ns)
70  , m_Tmax(m_Params->get_double_param("tmax") * ns)
71  , m_EdepSum(0)
72  , m_EabsSum(0)
73  , m_EionSum(0)
74 {
75 }
76 
77 //____________________________________________________________________________..
79 {
80  // if the last hit was a zero energie deposit hit, it is just reset
81  // and the memory is still allocated, so we need to delete it here
82  // if the last hit was saved, hit is a nullptr pointer which are
83  // legal to delete (it results in a no operation)
84  delete m_Hit;
85 }
86 
87 //____________________________________________________________________________..
88 // This is the implementation of the G4 UserSteppingAction
89 bool EICG4RPSteppingAction::UserSteppingAction(const G4Step *aStep, bool was_used)
90 {
93  // get volume of the current step
95  // IsInDetector(volume) returns
96  // == 0 outside of detector
97  // > 0 for hits in active volume
98  // < 0 for hits in passive material
99  int activeMaterial = m_Detector->IsInDetector(volume);
100  int virtualMaterial = m_Detector->IsInVirtualDetector(volume);
101 
102  if ( ! activeMaterial && ! virtualMaterial )
103  {
104  return false;
105  }
106 
107  // collect energy and track length step by step
108  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
109  // G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
110  const G4Track *aTrack = aStep->GetTrack();
111  // if this detector stops everything, just put all kinetic energy into edep
112  if (m_BlackHoleFlag)
113  {
114  if ((!std::isfinite(m_Tmin) && !std::isfinite(m_Tmax)) ||
115  aTrack->GetGlobalTime() < m_Tmin ||
116  aTrack->GetGlobalTime() > m_Tmax)
117  {
118  edep = aTrack->GetKineticEnergy() / GeV;
119  G4Track *killtrack = const_cast<G4Track *>(aTrack);
120  killtrack->SetTrackStatus(fStopAndKill);
121  }
122  }
123  // we use here only one detector in this simple example
124  // if you deal with multiple detectors in this stepping action
125  // the detector id can be used to distinguish between them
126  // hits can easily be analyzed later according to their detector id
127  int layer_id = m_Detector->get_Layer();
128 
129  if (!m_ActiveFlag)
130  {
131  return false;
132  }
133  bool geantino = false;
134  // the check for the pdg code speeds things up, I do not want to make
135  // an expensive string compare for every track when we know
136  // geantino or chargedgeantino has pid=0
137  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
138  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") !=
139  std::string::npos) // this also accounts for "chargedgeantino"
140  {
141  geantino = true;
142  }
143  G4StepPoint *prePoint = aStep->GetPreStepPoint();
144  G4StepPoint *postPoint = aStep->GetPostStepPoint();
145 
146  // Here we have to decide if we need to create a new hit. Normally this should
147  // only be neccessary if a G4 Track enters a new volume or is freshly created
148  // For this we look at the step status of the prePoint (beginning of the G4 Step).
149  // This should be either fGeomBoundary (G4 Track crosses into volume) or
150  // fUndefined (G4 Track newly created)
151  // Sadly over the years with different G4 versions we have observed cases where
152  // G4 produces "impossible hits" which we try to catch here
153  // These errors were always rare and it is not clear if they still exist but we
154  // still check for them for safety. We can reproduce G4 runs identically (if given
155  // the sequence of random number seeds you find in the log), the printouts help
156  // us giving the G4 support information about those failures
157  //
158  switch (prePoint->GetStepStatus())
159  {
160  case fPostStepDoItProc:
162  {
163  // this is the okay case, fPostStepDoItProc called in a volume, not first thing inside
164  // a new volume, just proceed here
165  break;
166  }
167  else
168  {
169  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
170  // this is still with us
171  std::cout << GetName() << ": New Hit for " << std::endl;
172  std::cout << "prestep status: "
174  << ", poststep status: "
176  << ", last pre step status: "
178  << ", last post step status: "
180  std::cout << "last track: " << m_SaveTrackId
181  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
182  std::cout << "phys pre vol: " << volume->GetName()
183  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
184  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
185  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
186  }
187  break;
188  // These are the normal cases
189  case fGeomBoundary:
190  case fUndefined:
191  if (!m_Hit)
192  {
193  m_Hit = new PHG4Hitv1();
194  }
195  m_Hit->set_layer((unsigned int) layer_id);
196  // here we set the entrance values in cm
197  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
198  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
199  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
200 
201  m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
202  m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
203  m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
204 
205  // time in ns
206  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
207  // set the track ID
208  m_Hit->set_trkid(aTrack->GetTrackID());
209  m_SaveTrackId = aTrack->GetTrackID();
210  // set the initial energy deposit
211  m_EdepSum = 0;
212  m_EabsSum = 0;
213 
214  m_Hit->set_edep(0);
215  if (!geantino && !m_BlackHoleFlag)
216  {
217  m_Hit->set_eion(0);
218  }
220  {
222  }
223  // implement your own here://
224  // add the properties you are interested in via set_XXX methods
225  // you can find existing set methods in $OFFLINE_MAIN/include/g4main/PHG4Hit.h
226  // this is initialization of your value. This is not needed you can just set the final
227  // value at the last step in this volume later one
228  /* if (whichactive > 0)
229  {
230  m_EionSum = 0; // assuming the ionization energy is only needed for active
231  // volumes (scintillators)
232  m_Hit->set_eion(0);
233  m_SaveHitContainer = m_HitContainer;
234  }
235  else
236  {
237  std::cout << "implement stuff for whichactive < 0 (inactive volumes)" << std::endl;
238  gSystem->Exit(1);
239  }
240  */
241  // this is for the tracking of the truth info
243  {
244  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
245  {
246  m_Hit->set_trkid(pp->GetUserTrackId());
247  m_Hit->set_shower_id(pp->GetShower()->get_id());
248  m_SaveShower = pp->GetShower();
249  }
250  }
251 
252  break;
253  default:
254  break;
255  }
256 
257  // This section is called for every step
258  // some sanity checks for inconsistencies (aka bugs) we have seen over the years
259  // check if this hit was created, if not print out last post step status
260  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
261  {
262  std::cout << GetName() << ": hit was not created" << std::endl;
263  std::cout << "prestep status: "
265  << ", poststep status: "
267  << ", last pre step status: "
269  << ", last post step status: "
271  std::cout << "last track: " << m_SaveTrackId
272  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
273  std::cout << "phys pre vol: " << volume->GetName()
274  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
275  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
276  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
277  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
278  gSystem->Exit(1);
279  }
280  m_SavePostStepStatus = postPoint->GetStepStatus();
281  // check if track id matches the initial one when the hit was created
282  if (aTrack->GetTrackID() != m_SaveTrackId)
283  {
284  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
285  std::cout << "saved track: " << m_SaveTrackId
286  << ", current trackid: " << aTrack->GetTrackID()
287  << ", prestep status: " << prePoint->GetStepStatus()
288  << ", previous post step status: " << m_SavePostStepStatus << std::endl;
289  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
290  gSystem->Exit(1);
291  }
292 
293  // We need to cache a few things from one step to the next
294  // to identify impossible hits and subsequent debugging printout
295  m_SavePreStepStatus = prePoint->GetStepStatus();
296  m_SavePostStepStatus = postPoint->GetStepStatus();
298  m_SaveVolPost = touchpost->GetVolume();
299 
300  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
301  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
302  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
303 
304  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
305  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
306  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
307 
308  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
309 
310  //sum up the energy to get total deposited
312 
313  if (geantino)
314  {
315  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
316  }
317  else
318  {
319  if (!m_BlackHoleFlag)
320  {
321  double eion = edep - aStep->GetNonIonizingEnergyDeposit() / GeV;
322  m_Hit->set_eion(m_Hit->get_eion() + eion);
323  }
324  }
326  {
327  double light_yield = GetVisibleEnergyDeposition(aStep);
328  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
329  }
330  if (edep > 0 || m_SaveAllHitsFlag)
331  {
333  {
334  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
335  {
336  pp->SetKeep(1); // we want to keep the track
337  }
338  }
339  }
340  // here we just update the exit values, it will be overwritten
341  // for every step until we leave the volume or the particle
342  // ceases to exist
343  // sum up the energy to get total deposited
349  // if any of these conditions is true this is the last step in
350  // this volume and we need to save the hit
351  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
352  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
353  // (happens when your detector goes outside world volume)
354  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
355  // aTrack->GetTrackStatus() == fStopAndKill is also set)
356  // aTrack->GetTrackStatus() == fStopAndKill: track ends
357  if (postPoint->GetStepStatus() == fGeomBoundary ||
358  postPoint->GetStepStatus() == fWorldBoundary ||
359  postPoint->GetStepStatus() == fAtRestDoItProc ||
360  aTrack->GetTrackStatus() == fStopAndKill ||
361  m_UseG4StepsFlag > 0)
362  {
363  // save only hits with energy deposit (or geantino)
364  if (m_Hit->get_edep() || m_SaveAllHitsFlag)
365  {
366  // update values at exit coordinates and set keep flag
367  // of track to keep
368  m_Hit->set_layer(layer_id);
369  m_Hit->set_hit_type( 0 );
371 
372  if( activeMaterial )
373  {
374  m_HitContainer->AddHit(layer_id, m_Hit);
375  }
376  if( virtualMaterial )
377  {
378  m_HitContainerVirt->AddHit(layer_id, m_Hit);
379  }
380 
381  if (m_SaveShower)
382  {
384  }
385  m_Hit = nullptr;
386  }
387  else
388  {
389  // if this hit has no energy deposit, just reset it for reuse
390  // this means we have to delete it in the dtor. If this was
391  // the last hit we processed the memory is still allocated
392  m_Hit->Reset();
393  }
394  }
395  // return true to indicate the hit was used
396  return true;
397 }
398 
399 //____________________________________________________________________________..
401 {
402  //std::string hitnodename = "G4HIT_" + m_Detector->GetName();
403  // std::cout << " ---> !!! hitnodename: " << hitnodename << std::endl;
404  // now look for the map and grab a pointer to it.
405  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
406  m_HitContainerVirt = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeNameVirt);
407  // if we do not find the node we need to make it.
408  if (!m_HitContainer)
409  {
410  std::cout << "EICG4RPSteppingAction::SetTopNode - unable to find "
411  << m_HitNodeName << std::endl;
412  }
413 }
415 {
417  {
418  return true;
419  }
420  return false;
421 }