Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TpcEndCapSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TpcEndCapSteppingAction.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 //____________________________________________________________________________..
44  : PHG4SteppingAction(detector->GetName())
45  , m_Detector(detector)
46  , m_Params(parameters)
47  , m_ActiveFlag(m_Params->get_int_param("active"))
48  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
49 {
50 }
51 
52 //____________________________________________________________________________..
54 {
55  // if the last hit was a zero energie deposit hit, it is just reset
56  // and the memory is still allocated, so we need to delete it here
57  // if the last hit was saved, hit is a nullptr pointer which are
58  // legal to delete (it results in a no operation)
59  delete m_Hit;
60 }
61 
62 //____________________________________________________________________________..
63 // This is the implementation of the G4 UserSteppingAction
64 bool PHG4TpcEndCapSteppingAction::UserSteppingAction(const G4Step *aStep, bool /*was_used*/)
65 {
66  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
67  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
68  // get volume of the current step
69  G4VPhysicalVolume *volume = touch->GetVolume();
70  // IsInDetector(volume) returns
71  // == 0 outside of detector
72  // > 0 for hits in active volume
73  // < 0 for hits in passive material
74  int whichactive = m_Detector->IsInDetector(volume);
75  if (!whichactive)
76  {
77  return false;
78  }
79  // collect energy and track length step by step
80  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
81  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
82  const G4Track *aTrack = aStep->GetTrack();
83  // if this detector stops everything, just put all kinetic energy into edep
84  if (m_BlackHoleFlag)
85  {
86  edep = aTrack->GetKineticEnergy() / GeV;
87  G4Track *killtrack = const_cast<G4Track *>(aTrack);
88  killtrack->SetTrackStatus(fStopAndKill);
89  }
90  // we use here only one detector in this simple example
91  // if you deal with multiple detectors in this stepping action
92  // the detector id can be used to distinguish between them
93  // hits can easily be analyzed later according to their detector id
94  int detector_id = 0; // we use here only one detector in this simple example
95  bool geantino = false;
96  // the check for the pdg code speeds things up, I do not want to make
97  // an expensive string compare for every track when we know
98  // geantino or chargedgeantino has pid=0
99  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
100  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") !=
101  std::string::npos) // this also accounts for "chargedgeantino"
102  {
103  geantino = true;
104  }
105  G4StepPoint *prePoint = aStep->GetPreStepPoint();
106  G4StepPoint *postPoint = aStep->GetPostStepPoint();
107 
108  // Here we have to decide if we need to create a new hit. Normally this should
109  // only be neccessary if a G4 Track enters a new volume or is freshly created
110  // For this we look at the step status of the prePoint (beginning of the G4 Step).
111  // This should be either fGeomBoundary (G4 Track crosses into volume) or
112  // fUndefined (G4 Track newly created)
113  // Sadly over the years with different G4 versions we have observed cases where
114  // G4 produces "impossible hits" which we try to catch here
115  // These errors were always rare and it is not clear if they still exist but we
116  // still check for them for safety. We can reproduce G4 runs identically (if given
117  // the sequence of random number seeds you find in the log), the printouts help
118  // us giving the G4 support information about those failures
119  //
120  switch (prePoint->GetStepStatus())
121  {
122  case fPostStepDoItProc:
123  if (m_SavePostStepStatus != fGeomBoundary)
124  {
125  // this is the okay case, fPostStepDoItProc called in a volume, not first thing inside
126  // a new volume, just proceed here
127  break;
128  }
129  else
130  {
131  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
132  // this is still with us
133  std::cout << GetName() << ": New Hit for " << std::endl;
134  std::cout << "prestep status: "
135  << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
136  << ", poststep status: "
137  << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
138  << ", last pre step status: "
140  << ", last post step status: "
142  std::cout << "last track: " << m_SaveTrackId
143  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
144  std::cout << "phys pre vol: " << volume->GetName()
145  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
146  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
147  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
148  }
149  [[fallthrough]];
150  // These are the normal cases
151  case fGeomBoundary:
152  case fUndefined:
153  if (!m_Hit)
154  {
155  m_Hit = new PHG4Hitv1();
156  }
157  m_Hit->set_layer(detector_id);
158  // here we set the entrance values in cm
159  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
160  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
161  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
162  // time in ns
163  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
164  // set the track ID
165  m_Hit->set_trkid(aTrack->GetTrackID());
166  m_SaveTrackId = aTrack->GetTrackID();
167  // set the initial energy deposit
168  m_EdepSum = 0;
169  // implement your own here://
170  // add the properties you are interested in via set_XXX methods
171  // you can find existing set methods in $OFFLINE_MAIN/include/g4main/PHG4Hit.h
172  // this is initialization of your value. This is not needed you can just set the final
173  // value at the last step in this volume later one
174  if (whichactive > 0)
175  {
176  m_EionSum = 0; // assuming the ionization energy is only needed for active
177  // volumes (scintillators)
178  m_Hit->set_eion(0);
180  }
181  else
182  {
183  std::cout << "implement stuff for whichactive < 0 (inactive volumes)" << std::endl;
184  gSystem->Exit(1);
185  }
186  // this is for the tracking of the truth info
187  if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
188  {
189  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
190  {
191  m_Hit->set_trkid(pp->GetUserTrackId());
192  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
193  }
194  }
195  break;
196  default:
197  break;
198  }
199 
200  // This section is called for every step
201  // some sanity checks for inconsistencies (aka bugs) we have seen over the years
202  // check if this hit was created, if not print out last post step status
203  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
204  {
205  std::cout << GetName() << ": hit was not created" << std::endl;
206  std::cout << "prestep status: "
207  << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
208  << ", poststep status: "
209  << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
210  << ", last pre step status: "
212  << ", last post step status: "
214  std::cout << "last track: " << m_SaveTrackId
215  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
216  std::cout << "phys pre vol: " << volume->GetName()
217  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
218  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
219  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
220  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
221  gSystem->Exit(1);
222  }
223  // check if track id matches the initial one when the hit was created
224  if (aTrack->GetTrackID() != m_SaveTrackId)
225  {
226  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
227  std::cout << "saved track: " << m_SaveTrackId
228  << ", current trackid: " << aTrack->GetTrackID()
229  << ", prestep status: " << prePoint->GetStepStatus()
230  << ", previous post step status: " << m_SavePostStepStatus << std::endl;
231  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
232  gSystem->Exit(1);
233  }
234 
235  // We need to cache a few things from one step to the next
236  // to identify impossible hits and subsequent debugging printout
237  m_SavePreStepStatus = prePoint->GetStepStatus();
238  m_SavePostStepStatus = postPoint->GetStepStatus();
239  m_SaveVolPre = volume;
240  m_SaveVolPost = touchpost->GetVolume();
241  // here we just update the exit values, it will be overwritten
242  // for every step until we leave the volume or the particle
243  // ceases to exist
244  // sum up the energy to get total deposited
245  m_EdepSum += edep;
246  if (whichactive > 0)
247  {
248  m_EionSum += eion;
249  }
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);
272  if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
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  if (whichactive > 0)
286  {
287  m_Hit->set_eion(-1);
288  }
289  }
290  else
291  {
293  }
294  //implement your own here://
295  // what you set here will be saved in the output
296  if (whichactive > 0)
297  {
299  }
300  m_SaveHitContainer->AddHit(detector_id, m_Hit);
301  // ownership has been transferred to container, set to null
302  // so we will create a new hit for the next track
303  m_Hit = nullptr;
304  }
305  else
306  {
307  // if this hit has no energy deposit, just reset it for reuse
308  // this means we have to delete it in the dtor. If this was
309  // the last hit we processed the memory is still allocated
310  m_Hit->Reset();
311  }
312  }
313  // return true to indicate the hit was used
314  return true;
315 }
316 
317 //____________________________________________________________________________..
319 {
320  // now look for the map and grab a pointer to it.
321  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
322  // if we do not find the node we need to make it.
323  if (!m_HitContainer)
324  {
325  std::cout << "PHG4TpcEndCapSteppingAction::SetTopNode - unable to find "
326  << m_HitNodeName << std::endl;
327  gSystem->Exit(1);
328  }
329 }
330 
332 {
333  if (type == "G4HIT")
334  {
336  return;
337  }
338  std::cout << "Invalid output hit node type " << type << std::endl;
339  gSystem->Exit(1);
340  return;
341 }