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