ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4dRICHSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4dRICHSteppingAction.cc
2 #include "EICG4dRICHDetector.h"
3 #include "EICG4dRICHHit.h"
4 
5 #include <phparameter/PHParameters.h>
6 
7 #include <fun4all/Fun4AllBase.h>
8 
10 
11 #include <g4main/PHG4Hit.h>
13 #include <g4main/PHG4Hitv1.h>
14 #include <g4main/PHG4Shower.h>
17 
18 #include <phool/getClass.h>
19 
20 #include <TSystem.h>
21 
22 #include <G4ParticleDefinition.hh>
24 #include <G4Step.hh>
25 #include <G4StepStatus.hh>
26 #include <G4SystemOfUnits.hh>
27 #include <G4ThreeVector.hh>
28 #include <G4TouchableHandle.hh>
29 #include <G4TrackStatus.hh>
30 #include <G4Types.hh>
31 #include <G4VPhysicalVolume.hh>
32 #include <G4VProcess.hh>
33 #include <G4VTouchable.hh>
35 
36 #include <cmath>
37 #include <iostream>
38 #include <string>
39 
40 class PHCompositeNode;
41 
42 //____________________________________________________________________________..
44  const PHParameters *parameters)
45  : PHG4SteppingAction(detector->GetName())
46  , m_Detector(detector)
47  , m_Params(parameters)
48  , m_HitContainer(nullptr)
49  , m_Hit(nullptr)
50  , m_SaveHitContainer(nullptr)
51  , m_SaveVolPre(nullptr)
52  , m_SaveVolPost(nullptr)
53  , m_SaveTrackId(-1)
54  , m_SavePreStepStatus(-1)
55  , m_SavePostStepStatus(-1)
56  , m_ActiveFlag(m_Params->get_int_param("active"))
57  , m_EdepSum(0)
58  , m_EionSum(0)
59  , hitType(-1)
60  , hitSubtype(-1)
61 {
62  // hit type strings
63  hitTypeStr[hEntrance] = "entrance";
64  hitTypeStr[hExit] = "exit";
65  hitTypeStr[hPSST] = "psst";
66  hitTypeStr[hIgnore] = "ignore";
67  // hit subtype strings
68  // - entrances
69  hitSubtypeStr[entPrimary] = "primary";
70  hitSubtypeStr[entSecondary] = "secondary";
71  hitSubtypeStr[entPostStep] = "postStep";
72  // - exits
73  hitSubtypeStr[exPrimary] = "primary";
74  hitSubtypeStr[exSecondary] = "secondary";
75  // - photosensor hits
76  hitSubtypeStr[psOptical] = "optical";
77  hitSubtypeStr[psGamma] = "gamma";
78  hitSubtypeStr[psOther] = "other";
79  // - unknown
80  hitSubtypeStr[subtypeUnknown] = "unknown";
81 }
82 
83 //____________________________________________________________________________..
85 {
86  // if the last hit was a zero energy deposit hit, it is
87  // just reset and the memory is still allocated, so we
88  // need to delete it here if the last hit was saved, hit
89  // is a nullptr pointer which are legal to delete (it
90  // results in a no operation)
91  delete m_Hit;
92 }
93 
94 //____________________________________________________________________________..
95 // This is the implementation of the G4 UserSteppingAction
97  bool was_used)
98 {
100  {
101  std::cout << "[>>>>>] call EICG4dRICHSteppingAction::UserSteppingAction" << std::endl;
102  }
103 
104  // get touchables, points, volumes
106  G4TouchableHandle postTouch = aStep->GetPostStepPoint()->GetTouchableHandle();
107  G4VPhysicalVolume *preVol = preTouch->GetVolume();
108  G4VPhysicalVolume *postVol = postTouch->GetVolume();
109  G4StepPoint *prePoint = aStep->GetPreStepPoint();
110  G4StepPoint *postPoint = aStep->GetPostStepPoint();
111 
112  // skip this step, if leaving the world (postVol will be nullptr)
113  if (postPoint->GetStepStatus() == fWorldBoundary)
114  {
115  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "... skip this step (leaving world)" << std::endl;
116  if (m_Hit) m_Hit->Reset();
117  return false;
118  }
119 
120  // get volume names
121  G4String prePointVolName = prePoint->GetPhysicalVolume()->GetName();
122  G4String postPointVolName = postPoint->GetPhysicalVolume()->GetName();
123  G4String preTouchVolName = preVol->GetName();
124  G4String postTouchVolName = postVol->GetName();
125 
126  // get track
127  const G4Track *aTrack = aStep->GetTrack();
128  G4String particleName = aTrack->GetParticleDefinition()->GetParticleName();
130  {
131  std::cout << "[-] track ID=" << aTrack->GetTrackID()
132  << ", particle=" << particleName << std::endl;
133  }
134 
135  // IsInDetector(preVol) returns
136  // == 0 outside of detector
137  // > 0 for hits in active volume
138  // < 0 for hits in passive material
139  int whichactive = m_Detector->IsInDetector(preVol);
141  {
142  std::cout << "[_] step preVol=" << preTouchVolName
143  << ", postVol=" << postTouchVolName << ", whichactive=" << whichactive
144  << std::endl;
145  }
146 
147  // reset hit classifiers
148  //hitType = -1;
149  //hitSubtype = -1;
150 
151  // classify hit type
152  if (prePointVolName.contains("dRICHpetal") && postPointVolName.contains("dRICHpsst"))
153  {
154  hitType = hPSST;
155  }
156  else if (prePointVolName.contains("World") && postPointVolName.contains("dRICHvessel"))
157  {
158  hitType = hEntrance;
159  }
160  else if (prePointVolName.contains("dRICHvessel") && postPointVolName.contains("World"))
161  {
162  hitType = hExit;
163  }
164  else
165  {
166  hitType = hIgnore;
167  }
168 
170  {
171  std::cout << "[__] step is ENTERING vessel" << std::endl;
172  }
173 
174  // skip this step, if it's outside the detector, and not an entrance
175  // or exit of the vessel
176  if (!whichactive && hitType != hEntrance && hitType != hExit)
177  {
178  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "... skip this step" << std::endl;
179  return false;
180  }
181 
182  // get step energy // TODO: do we need `eion`?
183  G4double edep = 0;
184  G4double eion = 0;
185  if (hitType != hEntrance)
186  {
187  edep = aStep->GetTotalEnergyDeposit() / GeV;
188  eion = (aStep->GetTotalEnergyDeposit() -
189  aStep->GetNonIonizingEnergyDeposit()) /
190  GeV;
191  }
192  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "[_] step edep=" << edep << ", eion=" << eion << std::endl;
193 
194  // Here we have to decide if we need to create a new hit. Normally this
195  // should only be neccessary if a G4 Track enters a new volume or is freshly
196  // created For this we look at the step status of the prePoint (beginning of
197  // the G4 Step). This should be either fGeomBoundary (G4 Track crosses into
198  // volume) or fUndefined (G4 Track newly created) Sadly over the years with
199  // different G4 versions we have observed cases where G4 produces "impossible
200  // hits" which we try to catch here These errors were always rare and it is
201  // not clear if they still exist but we still check for them for safety. We
202  // can reproduce G4 runs identically (if given the sequence of random number
203  // seeds you find in the log), the printouts help us giving the G4 support
204  // information about those failures
205  switch (prePoint->GetStepStatus())
206  {
207  // --- abnormal cases
208  case fPostStepDoItProc: // step defined by PostStepDoItVector
210  {
211  // this is the okay case, fPostStepDoItProc called in a volume, not first
212  // thing inside a new volume, just proceed here
214  std::cout << "[__] first step in a new volume" << std::endl;
215  }
216  else
217  {
218  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "[ + ] step was defined by PostStepDoItVector" << std::endl;
219  if (hitType != hEntrance)
220  {
221  // this is an impossible G4 Step print out diagnostic to help debug, not
222  // sure if this is still with us
223  std::cerr << "ERROR: impossible G4 Step" << std::endl;
224  std::cout << GetName() << ": New Hit for " << std::endl;
225  std::cout << "prestep status: "
227  << ", poststep status: "
229  << ", last pre step status: "
231  << ", last post step status: "
233  << std::endl;
234  std::cout << "last track: " << m_SaveTrackId
235  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
236  std::cout << "phys pre vol: " << preTouchVolName
237  << " post vol : " << postTouch->GetVolume()->GetName() << std::endl;
238  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
239  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
240  }
241  else
242  {
244  }
245  }
246 
247  // if this step is incident on the vessel, and we have not yet created a
248  // hit, create one
249  if (hitType == hEntrance)
250  {
251  m_Hit = nullptr; // kill any leftover hit
252  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "[++++] NEW hit (entrance)" << std::endl;
253  m_Hit = new EICG4dRICHHit();
254  this->InitHit(prePoint, aTrack, true);
255  }
256  break;
257 
258  // --- normal cases
259  case fGeomBoundary:
260  case fUndefined:
261  default:
262 
263  // do nothing if not geometry boundary, not undefined, and not entrance
264  if (prePoint->GetStepStatus() != fGeomBoundary &&
265  prePoint->GetStepStatus() != fUndefined && hitType != hEntrance)
266  {
267  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "[+] prepoint status ignored" << std::endl;
268  break;
269  }
270 
271  // create new hit
272  if (!m_Hit)
273  {
274  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "[++++] NEW hit" << std::endl;
275  m_Hit = new EICG4dRICHHit();
276  this->InitHit(prePoint, aTrack, true);
277  }
278  else
279  {
280  // if hit already exists, then Reset() has likely just been
281  // called; in this case, initialize, but don't reset accumulators
282  if (m_Hit->get_trkid() != aTrack->GetTrackID())
283  this->InitHit(prePoint, aTrack, true);
284  else
285  this->InitHit(prePoint, aTrack, false);
286  }
287 
288  // print info
290  {
291  std::cout << "[+] prepoint status=";
292  if (prePoint->GetStepStatus() == fGeomBoundary)
293  std::cout << "fGeomBoundary";
294  else if (prePoint->GetStepStatus() == fUndefined)
295  std::cout << "fUndefined";
296  else
297  std::cout << "UNKNOWN!";
298  std::cout << std::endl;
299  if (aTrack->GetTrackID() > 1 && aTrack->GetCreatorProcess())
300  {
301  std::cout << "[-] secondary track, creator process="
302  << aTrack->GetCreatorProcess()->GetProcessName();
303  }
304  else
305  {
306  std::cout << "[-] primary track, particle=" << particleName;
307  }
308  std::cout << std::endl;
309  }
310 
311  // tracking of the truth info // TODO: not used yet?
313  {
314  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
315  {
316  m_Hit->set_trkid(pp->GetUserTrackId());
317  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
318  }
319  }
320  break;
321  }
322 
323  // This section is called for every step -------------------------------------
324  // - some sanity checks for inconsistencies (aka bugs) we have seen over the years
325  // - check if this hit was created, if not print out last post step status
326  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
327  {
328  std::cout << "m_Hit @" << static_cast<void *>(m_Hit) << " isfinite = " << std::isfinite(m_Hit->get_x(0)) << std::endl;
329  std::cout << GetName() << ": hit was not created" << std::endl;
330  std::cout << "prestep status: "
332  << ", poststep status: "
334  << ", last pre step status: "
336  << ", last post step status: "
338  std::cout << "last track: " << m_SaveTrackId
339  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
340  std::cout << "phys pre vol: " << preTouchVolName
341  << " post vol : " << postTouch->GetVolume()->GetName() << std::endl;
342  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
343  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
344  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
345  gSystem->Exit(1);
346  }
347  // check if track id matches the initial one when the hit was created
348  if (aTrack->GetTrackID() != m_SaveTrackId)
349  {
350  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
351  std::cout << "saved track: " << m_SaveTrackId
352  << ", current trackid: " << aTrack->GetTrackID()
353  << ", prestep status: " << prePoint->GetStepStatus()
354  << ", previous post step status: " << m_SavePostStepStatus << std::endl;
355  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
356  gSystem->Exit(1);
357  }
358  // We need to cache a few things from one step to the next
359  // to identify impossible hits and subsequent debugging printout
360  m_SavePreStepStatus = prePoint->GetStepStatus();
361  m_SavePostStepStatus = postPoint->GetStepStatus();
362  m_SaveVolPre = preVol;
363  m_SaveVolPost = postVol;
364 
365  // update accumulators
366  m_EdepSum += edep;
367  if (whichactive > 0)
368  {
369  m_EionSum += eion;
370  }
372  {
373  std::cout << "[_] m_EdepSum=" << m_EdepSum << ", m_EionSum=" << m_EionSum << std::endl;
374  }
375 
376  // get petal number // TODO: does not work for entrance/exit hits
377  // since they are identified by world<->vessel
378  // crossings; the vessel has no petal number;
379  // note that there is currently a gap between
380  // the vessel and the petal volumes, which may
381  // be unrealistic
382  int petal = hitType == hEntrance ? m_Detector->GetPetal(postVol)
383  : m_Detector->GetPetal(preVol);
384 
385  // save the hit ---------------------------------------------
386  // if any of these conditions is true, this is the last step in
387  // this volume and we consider saving the hit
388  if (postPoint->GetStepStatus() == fGeomBoundary /*left volume*/
389  || postPoint->GetStepStatus() == fWorldBoundary /*left world*/
390  || postPoint->GetStepStatus() == fAtRestDoItProc /*track stops*/
391  || aTrack->GetTrackStatus() == fStopAndKill) /*track ends*/
392  {
394  {
395  std::cout << "[---+] last step in the volume (pre=" << prePointVolName << ", post=" << postPointVolName << ")" << std::endl;
396  }
397 
398  // hits to keep +++++++++++++++++++++++
399  if (hitType != hIgnore)
400  {
401  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "[-+] " << hitTypeStr[hitType] << " hit, KEEP!" << std::endl;
402 
403  // classify hit subtype
404  switch (hitType)
405  {
406  case hEntrance:
407  if (hitSubtype != entPostStep)
408  {
409  if (aTrack->GetTrackID() == 1)
411  else
413  }
414  break;
415  case hExit:
416  if (aTrack->GetTrackID() == 1)
418  else
420  break;
421  case hPSST:
422  if (particleName == "opticalphoton")
424  else if (particleName == "gamma")
426  else
428  break;
429  default:
431  }
432  if (hitSubtype == -1) hitSubtype = subtypeUnknown;
433 
434  // set hit vars
437  m_Hit->set_petal(petal);
438  m_Hit->set_psst(m_Detector->GetPSST(postVol));
440  m_Hit->set_particle_name(particleName);
441 
442  switch (hitType)
443  {
444  case hEntrance:
445  if (hitSubtype == entPostStep)
446  m_Hit->set_process("postStep");
447  else if (hitSubtype == entPrimary)
448  m_Hit->set_process("primary");
449  else if (aTrack->GetCreatorProcess())
451  else m_Hit->set_process("primary");
452  break;
453  case hExit:
454  m_Hit->set_process("exitProcess");
455  break;
456  default:
457  if (aTrack->GetCreatorProcess())
459  else
460  m_Hit->set_process("unknown");
461  }
462 
463  m_Hit->set_parent_id(aTrack->GetParentID());
464  m_Hit->set_position(1, postPoint->GetPosition() / cm);
465  m_Hit->set_momentum(aTrack->GetMomentum() / GeV);
466  m_Hit->set_momentum_dir(aTrack->GetMomentumDirection() /*unitless*/);
469  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
470 
471  // tracking of the truth info // TODO: not used yet?
473  {
474  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
475  {
476  pp->SetKeep(1); // we want to keep the track
477  }
478  }
479 
480  // total accumulators
484  {
485  std::cout << "[-+] m_EdepSum=" << m_EdepSum << ", m_EionSum=" << m_EionSum << std::endl;
486  }
487 
488  // transfer ownership to container
490  m_Hit = nullptr; // so that next track will create new hit
491  }
492  else
493  {
494  // do not save this hit ++++++++++++++++++++++++++++
495  // - reset hit object for reuse
496  // - if last hit overall, memory still allocated
497  // - does not reset local vars, such as m_EdepSum
498  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "[-+] not keeping this hit" << std::endl;
499  m_Hit->Reset();
500  }
501  }
502 
503  // return true to indicate the hit was used
504  return true;
505 }
506 
507 //____________________________________________________________________________..
509  const G4Track *aTrack_,
510  bool resetAccumulators)
511 {
512  // set some entrance attributes, and track ID
513  m_Hit->set_position(0, prePoint_->GetPosition() / cm);
514  m_Hit->set_t(0, prePoint_->GetGlobalTime() / nanosecond);
515  m_Hit->set_trkid(aTrack_->GetTrackID());
516  m_SaveTrackId = aTrack_->GetTrackID();
518  // initializate accumulators (e.g., for total energy deposition)
519  if (resetAccumulators)
520  {
521  m_EdepSum = 0;
522  m_EionSum = 0;
523  }
524 }
525 
526 //____________________________________________________________________________..
528 {
529  std::string hitnodename = "G4HIT_" + m_Detector->GetName();
530  // now look for the map and grab a pointer to it.
531  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
532  // if we do not find the node we need to make it.
533  if (!m_HitContainer)
534  {
535  std::cout << "EICG4dRICHSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
536  }
537 }