ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LBLVtxSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LBLVtxSteppingAction.cc
2 
3 #include "G4LBLVtxDetector.h"
4 
5 #include <phparameter/PHParameters.h>
6 
8 
9 #include <g4main/PHG4Hit.h>
11 #include <g4main/PHG4Hitv1.h>
12 #include <g4main/PHG4Shower.h>
13 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
15 
16 #include <phool/getClass.h>
17 
18 #include <TSystem.h>
19 
20 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
21 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
22 #include <Geant4/G4Step.hh>
23 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
24 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
25 #include <Geant4/G4String.hh> // for G4String
26 #include <Geant4/G4SystemOfUnits.hh>
27 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
28 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
29 #include <Geant4/G4Track.hh> // for G4Track
30 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
31 #include <Geant4/G4Types.hh> // for G4double
32 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
33 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
34 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
35 
36 #include <cmath> // for isfinite
37 #include <iostream>
38 #include <string> // for operator<<, string
39 
40 class PHCompositeNode;
41 
42 using namespace std;
43 //____________________________________________________________________________..
45  G4LBLVtxDetector *detector, const PHParameters *parameters)
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_ActiveFlag(m_Params->get_int_param("active"))
59  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
60  , m_EdepSum(0)
61  , m_EionSum(0)
62 {
63 }
64 
66 {
67  // if the last hit was a zero energie deposit hit, it is just reset
68  // and the memory is still allocated, so we need to delete it here
69  // if the last hit was saved, hit is a nullptr pointer which are
70  // legal to delete (it results in a no operation)
71  delete m_Hit;
72 }
73 
74 //____________________________________________________________________________..
76  bool was_used)
77 {
80  // get volume of the current step
82  // IsInDetector(volume) returns
83  // == 0 outside of detector
84  // > 0 for hits in active volume
85  // < 0 for hits in passive material
86  int whichactive = m_Detector->IsInDetector(volume);
87  if (!whichactive)
88  {
89  return false;
90  }
91 
92  // collect energy and track length step by step
94  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
95  const G4Track *aTrack = aStep->GetTrack();
96 
97  // if this block stops everything, just put all kinetic energy into edep
98  if (m_BlackHoleFlag)
99  {
100  edep = aTrack->GetKineticEnergy() / GeV;
101  G4Track *killtrack = const_cast<G4Track *>(aTrack);
102  killtrack->SetTrackStatus(fStopAndKill);
103  }
104 
105  int detector_id = 0; // we use here only one detector in this simple example
106  bool geantino = false;
107  // the check for the pdg code speeds things up, I do not want to make
108  // an expensive string compare for every track when we know
109  // geantino or chargedgeantino has pid=0
110  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
111  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") !=
112  string::npos) // this also accounts for "chargedgeantino"
113  {
114  geantino = true;
115  }
116  G4StepPoint *prePoint = aStep->GetPreStepPoint();
117  G4StepPoint *postPoint = aStep->GetPostStepPoint();
118  // cout << "track id " << aTrack->GetTrackID() << endl;
119  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
120  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
121 
122  switch (prePoint->GetStepStatus())
123  {
124  case fPostStepDoItProc:
126  {
127  break;
128  }
129  else
130  {
131  // this is bad from G4 print out diagnostic to help debug, not sure if
132  // this is still with us
133  cout << GetName() << ": New Hit for " << endl;
134  cout << "prestep status: "
136  << ", poststep status: "
138  << ", last pre step status: "
140  << ", last post step status: "
142  cout << "last track: " << m_SaveTrackId
143  << ", current trackid: " << aTrack->GetTrackID() << endl;
144  cout << "phys pre vol: " << volume->GetName()
145  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
146  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
147  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
148  }
149  case fGeomBoundary:
150  case fUndefined:
151  if (!m_Hit)
152  {
153  m_Hit = new PHG4Hitv1();
154  }
155  m_Hit->set_layer(detector_id);
156  // here we set the entrance values in cm
157  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
158  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
159  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
160  // time in ns
161  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
162  // set the track ID
163  m_Hit->set_trkid(aTrack->GetTrackID());
164  m_SaveTrackId = aTrack->GetTrackID();
165  // set the initial energy deposit
166  m_EdepSum = 0;
167  if (whichactive > 0)
168  {
169  m_EionSum = 0; // assuming the ionization energy is only needed for active
170  // volumes (silicon)
171  m_Hit->set_eion(0);
173  }
174  else
175  {
176  // for the absorber we want only the energy deposition
178  }
179  // this is for the tracking of the truth info
181  {
182  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
183  {
184  m_Hit->set_trkid(pp->GetUserTrackId());
185  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(),
186  m_Hit->get_hit_id());
187  }
188  }
189 
190  break;
191  default:
192  break;
193  }
194 
195  // some sanity checks for inconsistencies (aka bugs)
196  // check if this hit was created, if not print out last post step status
197  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
198  {
199  cout << GetName() << ": hit was not created" << endl;
200  cout << "prestep status: "
202  << ", poststep status: "
204  << ", last pre step status: "
206  << ", last post step status: "
208  cout << "last track: " << m_SaveTrackId
209  << ", current trackid: " << aTrack->GetTrackID() << endl;
210  cout << "phys pre vol: " << volume->GetName()
211  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
212  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
213  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
214  gSystem->Exit(1);
215  }
216  // check if track id matches the initial one when the hit was created
217  if (aTrack->GetTrackID() != m_SaveTrackId)
218  {
219  cout << GetName() << ": hits do not belong to the same track" << endl;
220  cout << "saved track: " << m_SaveTrackId
221  << ", current trackid: " << aTrack->GetTrackID()
222  << ", prestep status: " << prePoint->GetStepStatus()
223  << ", previous post step status: " << m_SavePostStepStatus << endl;
224 
225  gSystem->Exit(1);
226  }
227  m_SavePreStepStatus = prePoint->GetStepStatus();
228  m_SavePostStepStatus = postPoint->GetStepStatus();
230  m_SaveVolPost = touchpost->GetVolume();
231 
232  // here we just update the exit values, it will be overwritten
233  // for every step until we leave the volume or the particle
234  // ceases to exist
235  // sum up the energy to get total deposited
236  m_EdepSum += edep;
237  if (whichactive > 0)
238  {
239  m_EionSum += eion;
240  }
241  // if any of these conditions is true this is the last step in
242  // this volume and we need to save the hit
243  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
244  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
245  // (happens when your detector goes outside world volume)
246  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
247  // aTrack->GetTrackStatus() == fStopAndKill is also set)
248  // aTrack->GetTrackStatus() == fStopAndKill: track ends
249  if (postPoint->GetStepStatus() == fGeomBoundary ||
250  postPoint->GetStepStatus() == fWorldBoundary ||
251  postPoint->GetStepStatus() == fAtRestDoItProc ||
252  aTrack->GetTrackStatus() == fStopAndKill)
253  {
254  // save only hits with energy deposit (or geantino)
255  if (m_EdepSum > 0 || geantino)
256  {
257  // update values at exit coordinates and set keep flag
258  // of track to keep
259  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
260  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
261  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
262 
263  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
265  {
266  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
267  {
268  pp->SetKeep(1); // we want to keep the track
269  }
270  }
271  if (geantino)
272  {
273  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way
274  // geantinos survive the g4hit compression
275  if (whichactive > 0)
276  {
277  m_Hit->set_eion(-1);
278  }
279  }
280  else
281  {
283  }
284  if (whichactive > 0)
285  {
287  }
288  m_SaveHitContainer->AddHit(detector_id, m_Hit);
289  // ownership has been transferred to container, set to null
290  // so we will create a new hit for the next track
291  m_Hit = nullptr;
292  }
293  else
294  {
295  // if this hit has no energy deposit, just reset it for reuse
296  // this means we have to delete it in the dtor. If this was
297  // the last hit we processed the memory is still allocated
298  m_Hit->Reset();
299  }
300  }
301  // return true to indicate the hit was used
302  return true;
303 }
304 
305 //____________________________________________________________________________..
307 {
308  string hitnodename;
309  string absorbernodename;
310  if (m_Detector->SuperDetector() != "NONE")
311  {
312  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
313  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
314  }
315  else
316  {
317  hitnodename = "G4HIT_" + m_Detector->GetName();
318  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
319  }
320 
321  //now look for the map and grab a pointer to it.
322  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
323  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
324 
325  // if we do not find the node it's messed up.
326  if (!m_HitContainer)
327  {
328  cout << "G4LBLVtxSteppingAction::SetTopNode - unable to find " << hitnodename << endl;
329  }
330  // if absorber hits are turned off we will not find this node
331  // so only print out if verbosity is enabled.
332  if (!m_AbsorberHitContainer)
333  {
334  if (Verbosity() > 1)
335  {
336  cout << "G4LBLVtxSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
337  }
338  }
339 }