ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AllSiliconTrackerSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AllSiliconTrackerSteppingAction.cc
2 
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_AbsorberHitContainer(nullptr)
51  , m_Hit(nullptr)
52  , m_SaveHitContainer(nullptr)
53  , m_SaveVolPre(nullptr)
54  , m_SaveVolPost(nullptr)
55  , m_SaveTrackId(-1)
56  , m_SavePreStepStatus(-1)
57  , m_SavePostStepStatus(-1)
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
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  int detector_id = m_Detector->get_detid(volume, whichactive); // the detector id is coded into the IsInDetector return
101  // cout << "Name: " << volume->GetName() << endl;
102  // cout << "det id: " << whichactive << endl;
103  bool geantino = false;
104  // the check for the pdg code speeds things up, I do not want to make
105  // an expensive string compare for every track when we know
106  // geantino or chargedgeantino has pid=0
107  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
108  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") !=
109  string::npos) // this also accounts for "chargedgeantino"
110  {
111  geantino = true;
112  }
113  G4StepPoint *prePoint = aStep->GetPreStepPoint();
114  G4StepPoint *postPoint = aStep->GetPostStepPoint();
115 
116  // Here we have to decide if we need to create a new hit. Normally this should
117  // only be neccessary if a G4 Track enters a new volume or is freshly created
118  // For this we look at the step status of the prePoint (beginning of the G4 Step).
119  // This should be either fGeomBoundary (G4 Track crosses into volume) or
120  // fUndefined (G4 Track newly created)
121  // Sadly over the years with different G4 versions we have observed cases where
122  // G4 produces "impossible hits" which we try to catch here
123  // These errors were always rare and it is not clear if they still exist but we
124  // still check for them for safety. We can reproduce G4 runs identically (if given
125  // the sequence of random number seeds you find in the log), the printouts help
126  // us giving the G4 support information about those failures
127  //
128  switch (prePoint->GetStepStatus())
129  {
130  case fPostStepDoItProc:
132  {
133  // this is the okay case, fPostStepDoItProc called in a volume, not first thing inside
134  // a new volume, just proceed here
135  break;
136  }
137  else
138  {
139  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
140  // this is still with us
141  cout << GetName() << ": New Hit for " << endl;
142  cout << "prestep status: "
144  << ", poststep status: "
146  << ", last pre step status: "
148  << ", last post step status: "
150  cout << "last track: " << m_SaveTrackId
151  << ", current trackid: " << aTrack->GetTrackID() << endl;
152  cout << "phys pre vol: " << volume->GetName()
153  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
154  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
155  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
156  }
157  // These are the normal cases
158  case fGeomBoundary:
159  case fUndefined:
160  if (!m_Hit)
161  {
162  m_Hit = new PHG4Hitv1();
163  }
164  m_Hit->set_layer(detector_id);
165  // here we set the entrance values in cm
166  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
167  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
168  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
169  // time in ns
170  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
171  // set the track ID
172  m_Hit->set_trkid(aTrack->GetTrackID());
173  m_SaveTrackId = aTrack->GetTrackID();
174  // set the initial energy deposit
175  m_EdepSum = 0;
176  // implement your own here://
177  // add the properties you are interested in via set_XXX methods
178  // you can find existing set methods in $OFFLINE_MAIN/include/g4main/PHG4Hit.h
179  // this is initialization of your value. This is not needed you can just set the final
180  // value at the last step in this volume later one
181  if (whichactive > 0)
182  {
184  }
185  else
186  {
187  // all absorber hits go into one node
189  }
190  // this is for the tracking of the truth info
192  {
193  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
194  {
195  m_Hit->set_trkid(pp->GetUserTrackId());
196  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
197  }
198  }
199  break;
200  default:
201  break;
202  }
203 
204  // This section is called for every step
205  // some sanity checks for inconsistencies (aka bugs) we have seen over the years
206  // check if this hit was created, if not print out last post step status
207  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
208  {
209  cout << GetName() << ": hit was not created" << endl;
210  cout << "prestep status: "
212  << ", poststep status: "
214  << ", last pre step status: "
216  << ", last post step status: "
218  cout << "last track: " << m_SaveTrackId
219  << ", current trackid: " << aTrack->GetTrackID() << endl;
220  cout << "phys pre vol: " << volume->GetName()
221  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
222  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
223  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
224  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
225  gSystem->Exit(1);
226  }
227  // check if track id matches the initial one when the hit was created
228  if (aTrack->GetTrackID() != m_SaveTrackId)
229  {
230  cout << GetName() << ": hits do not belong to the same track" << endl;
231  cout << "saved track: " << m_SaveTrackId
232  << ", current trackid: " << aTrack->GetTrackID()
233  << ", prestep status: " << prePoint->GetStepStatus()
234  << ", previous post step status: " << m_SavePostStepStatus << endl;
235  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
236  gSystem->Exit(1);
237  }
238 
239  // We need to cache a few things from one step to the next
240  // to identify impossible hits and subsequent debugging printout
241  m_SavePreStepStatus = prePoint->GetStepStatus();
242  m_SavePostStepStatus = postPoint->GetStepStatus();
244  m_SaveVolPost = touchpost->GetVolume();
245  // here we just update the exit values, it will be overwritten
246  // for every step until we leave the volume or the particle
247  // ceases to exist
248  // sum up the energy to get total deposited
249  m_EdepSum += edep;
250  // if any of these conditions is true this is the last step in
251  // this volume and we need to save the hit
252  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
253  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
254  // (happens when your detector goes outside world volume)
255  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
256  // aTrack->GetTrackStatus() == fStopAndKill is also set)
257  // aTrack->GetTrackStatus() == fStopAndKill: track ends
258  if (postPoint->GetStepStatus() == fGeomBoundary ||
259  postPoint->GetStepStatus() == fWorldBoundary ||
260  postPoint->GetStepStatus() == fAtRestDoItProc ||
261  aTrack->GetTrackStatus() == fStopAndKill)
262  {
263  // save only hits with energy deposit (or geantino)
264  if (m_EdepSum > 0 || geantino)
265  {
266  // update values at exit coordinates and set keep flag
267  // of track to keep
268  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
269  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
270  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
271  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
273  {
274  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
275  {
276  pp->SetKeep(1); // we want to keep the track
277  }
278  }
279  if (geantino)
280  {
281  //implement your own here://
282  // if you want to do something special for geantinos (normally you do not)
283  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way
284  // geantinos survive the g4hit compression
285  }
286  else
287  {
289  }
290  m_SaveHitContainer->AddHit(detector_id, m_Hit);
291  // ownership has been transferred to container, set to null
292  // so we will create a new hit for the next track
293  m_Hit = nullptr;
294  }
295  else
296  {
297  // if this hit has no energy deposit, just reset it for reuse
298  // this means we have to delete it in the dtor. If this was
299  // the last hit we processed the memory is still allocated
300  m_Hit->Reset();
301  }
302  }
303  // return true to indicate the hit was used
304  return true;
305 }