ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4BwdSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4BwdSteppingAction.cc
1 //____________________________________________________________________________..
2 //
3 // This is a working template for the Stepping Action which needs to be implemented
4 // for active detectors. Most of the code is error handling and access to the G4 objects
5 // and our data structures. It does not need any adjustment. The only thing you need to
6 // do is to add the properties of the G4Hits you want to save for later analysis
7 // This needs to be done in 2 places, G4Hits are generated when a G4 track enters a new
8 // volume (or is created). Here you give it an initial value. When the G4 track leaves
9 // the volume the final value needs to be set.
10 // The places to do this is marked by //implement your own here//
11 //
12 // As guidance you can look at the total (integrated over all steps in a volume) energy
13 // deposit which should always be saved.
14 // Additionally the total ionization energy is saved - this can be removed if you are not
15 // interested in this. Naturally you may want remove these comments in your version
16 //
17 //____________________________________________________________________________..
18 
19 #include "EICG4BwdSteppingAction.h"
20 
21 #include "EICG4BwdDetector.h"
22 #include "EICG4BwdSubsystem.h"
23 
24 #include <phparameter/PHParameters.h>
25 
27 
28 #include <g4main/PHG4Hit.h>
30 #include <g4main/PHG4Hitv1.h>
31 #include <g4main/PHG4Shower.h>
34 
35 #include <phool/getClass.h>
36 #include <phool/PHCompositeNode.h>
37 
38 #include <TSystem.h>
39 
40 #include <Geant4/G4NavigationHistory.hh>
41 #include <Geant4/G4ParticleDefinition.hh>
42 #include <Geant4/G4ReferenceCountedHandle.hh>
43 #include <Geant4/G4Step.hh>
44 #include <Geant4/G4StepPoint.hh>
45 #include <Geant4/G4StepStatus.hh>
46 #include <Geant4/G4String.hh>
47 #include <Geant4/G4SystemOfUnits.hh>
48 #include <Geant4/G4ThreeVector.hh>
49 #include <Geant4/G4TouchableHandle.hh>
50 #include <Geant4/G4Track.hh>
51 #include <Geant4/G4TrackStatus.hh>
52 #include <Geant4/G4Types.hh>
53 #include <Geant4/G4VPhysicalVolume.hh>
54 #include <Geant4/G4VTouchable.hh>
55 #include <Geant4/G4VUserTrackInformation.hh>
56 
57 #include <cmath>
58 #include <iostream>
59 #include <string>
60 #include <cstdlib>
61 #include <iomanip>
62 
63 class PHCompositeNode;
64 
65 //____________________________________________________________________________..
67  : PHG4SteppingAction(detector->GetName())
68  , m_Subsystem(subsys)
69  , m_Detector(detector)
70  , m_Params(parameters)
71  , m_HitContainer(nullptr)
72  , m_Hit(nullptr)
73  , m_SaveShower(nullptr)
74  , m_SaveVolPre(nullptr)
75  , m_SaveVolPost(nullptr)
76  , m_SaveLightYieldFlag(m_Params->get_int_param("lightyield"))
77  , m_SaveTrackId(-1)
78  , m_SavePreStepStatus(-1)
79  , m_SavePostStepStatus(-1)
80  , m_ActiveFlag(m_Params->get_int_param("active"))
81  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
82  , m_UseG4StepsFlag(m_Params->get_int_param("use_g4steps"))
83  , m_Zmin(m_Params->get_double_param("place_z") * cm - m_Params->get_double_param("length") * cm / 2.)
84  , m_Zmax(m_Params->get_double_param("place_z") * cm + m_Params->get_double_param("length") * cm / 2.)
85  , m_Tmin(m_Params->get_double_param("tmin") * ns)
86  , m_Tmax(m_Params->get_double_param("tmax") * ns)
87  , m_EdepSum(0)
88  , m_EabsSum(0)
89  , m_EionSum(0)
90 {
91 // G4 seems to have issues in the um range
92  m_Zmin -= copysign(m_Zmin, 1. / 1e6 * cm);
93  m_Zmax += copysign(m_Zmax, 1. / 1e6 * cm);
94 }
95 
96 //____________________________________________________________________________..
98 {
99  // if the last hit was a zero energie deposit hit, it is just reset
100  // and the memory is still allocated, so we need to delete it here
101  // if the last hit was saved, hit is a nullptr pointer which are
102  // legal to delete (it results in a no operation)
103  delete m_Hit;
104 }
105 
106 //____________________________________________________________________________..
107 // This is the implementation of the G4 UserSteppingAction
108 bool EICG4BwdSteppingAction::UserSteppingAction(const G4Step *aStep, bool was_used)
109 {
111  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
112  // get volume of the current step
113  G4VPhysicalVolume *volume = touch->GetVolume();
114  // IsInDetector(volume) returns
115  // == 0 outside of detector
116  // > 0 for hits in active volume
117  // < 0 for hits in passive material
118 
119 
120  int whichactive = m_Detector->IsInDetector(volume);
121  if (!whichactive)
122  {
123  return false;
124  }
125  // collect energy and track length step by step
126 
127  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
128  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
129  G4double light_yield = GetVisibleEnergyDeposition(aStep);
130  const G4Track *aTrack = aStep->GetTrack();
131  // if this detector stops everything, just put all kinetic energy into edep
132  if (m_BlackHoleFlag)
133  {
134  if ((!std::isfinite(m_Tmin) && !std::isfinite(m_Tmax)) ||
135  aTrack->GetGlobalTime() < m_Tmin ||
136  aTrack->GetGlobalTime() > m_Tmax)
137  {
138  edep = aTrack->GetKineticEnergy() / GeV;
139  G4Track *killtrack = const_cast<G4Track *>(aTrack);
140  killtrack->SetTrackStatus(fStopAndKill);
141  }
142  }
143  // we use here only one detector in this simple example
144  // if you deal with multiple detectors in this stepping action
145  // the detector id can be used to distinguish between them
146  // hits can easily be analyzed later according to their detector id
147  int layer_id = m_Detector->get_Layer();
148  int idx_j = -1;
149  int idx_k = -1;
150 // deadhits for dead material
151  int layer_type = 1;
152  if (!m_ActiveFlag)
153  {
154  return false;
155  }
156  bool geantino = false;
157  // the check for the pdg code speeds things up, I do not want to make
158  // an expensive string compare for every track when we know
159  // geantino or chargedgeantino has pid=0
160  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
161  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") !=
162  std::string::npos) // this also accounts for "chargedgeantino"
163  {
164  geantino = true;
165  }
166  G4StepPoint *prePoint = aStep->GetPreStepPoint();
167  G4StepPoint *postPoint = aStep->GetPostStepPoint();
168 
169  // Here we have to decide if we need to create a new hit. Normally this should
170  // only be neccessary if a G4 Track enters a new volume or is freshly created
171  // For this we look at the step status of the prePoint (beginning of the G4 Step).
172  // This should be either fGeomBoundary (G4 Track crosses into volume) or
173  // fUndefined (G4 Track newly created)
174  // Sadly over the years with different G4 versions we have observed cases where
175  // G4 produces "impossible hits" which we try to catch here
176  // These errors were always rare and it is not clear if they still exist but we
177  // still check for them for safety. We can reproduce G4 runs identically (if given
178  // the sequence of random number seeds you find in the log), the printouts help
179  // us giving the G4 support information about those failures
180  //
181  switch (prePoint->GetStepStatus())
182  {
183  case fPostStepDoItProc:
185  {
186  // this is the okay case, fPostStepDoItProc called in a volume, not first thing inside
187  // a new volume, just proceed here
188  break;
189  }
190  else
191  {
192  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
193  // this is still with us
194  std::cout << GetName() << ": New Hit for " << std::endl;
195  std::cout << "prestep status: "
197  << ", poststep status: "
199  << ", last pre step status: "
201  << ", last post step status: "
203  std::cout << "last track: " << m_SaveTrackId
204  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
205  std::cout << "phys pre vol: " << volume->GetName()
206  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
207  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
208  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
209  }
210  break;
211  // These are the normal cases
212  case fGeomBoundary:
213  case fUndefined:
214  if (!m_Hit)
215  {
216  m_Hit = new PHG4Hitv1();
217  }
218  m_Hit->set_layer((unsigned int) layer_id);
219  // here we set the entrance values in cm
220  FindTowerIndexFromPosition(prePoint, idx_j, idx_k);
221  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
222  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
223  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
224 
225  m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
226  m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
227  m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
228 
229  // time in ns
230  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
231  // set the track ID
232  m_Hit->set_trkid(aTrack->GetTrackID());
233  m_SaveTrackId = aTrack->GetTrackID();
234  // set the initial energy deposit
235  m_EdepSum = 0;
236  m_EabsSum = 0;
237 
238  m_Hit->set_edep(0);
239  if (!geantino && !m_BlackHoleFlag)
240  {
241  m_Hit->set_eion(0);
242  }
244  {
246  m_Hit->set_index_j(idx_j);
247  m_Hit->set_index_k(idx_k);
248  }
249  // implement your own here://
250  // add the properties you are interested in via set_XXX methods
251  // you can find existing set methods in $OFFLINE_MAIN/include/g4main/PHG4Hit.h
252  // this is initialization of your value. This is not needed you can just set the final
253  // value at the last step in this volume later one
254  // this is for the tracking of the truth info
256  {
257  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
258  {
259  m_Hit->set_trkid(pp->GetUserTrackId());
260  m_Hit->set_shower_id(pp->GetShower()->get_id());
261  m_SaveShower = pp->GetShower();
262 // pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
263  }
264  }
265  if (!hasMotherSubsystem() && (m_Hit->get_z(0) * cm > m_Zmax || m_Hit->get_z(0) * cm < m_Zmin))
266  {
267  std::cout << m_Detector->SuperDetector() << std::setprecision(9)
268  << " PHG4CylinderSteppingAction: Entry hit z " << m_Hit->get_z(0) * cm
269  << " outside acceptance, zmin " << m_Zmin
270  << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
271  }
272  break;
273  default:
274  break;
275  }
276 
277  // This section is called for every step
278  // some sanity checks for inconsistencies (aka bugs) we have seen over the years
279  // check if this hit was created, if not print out last post step status
280  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
281  {
282  std::cout << GetName() << ": hit was not created" << std::endl;
283  std::cout << "prestep status: "
285  << ", poststep status: "
287  << ", last pre step status: "
289  << ", last post step status: "
291  std::cout << "last track: " << m_SaveTrackId
292  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
293  std::cout << "phys pre vol: " << volume->GetName()
294  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
295  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
296  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
297  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
298  gSystem->Exit(1);
299  }
300  m_SavePostStepStatus = postPoint->GetStepStatus();
301  // check if track id matches the initial one when the hit was created
302  if (aTrack->GetTrackID() != m_SaveTrackId)
303  {
304  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
305  std::cout << "saved track: " << m_SaveTrackId
306  << ", current trackid: " << aTrack->GetTrackID()
307  << ", prestep status: " << prePoint->GetStepStatus()
308  << ", previous post step status: " << m_SavePostStepStatus << std::endl;
309  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
310  gSystem->Exit(1);
311  }
312 
313  // We need to cache a few things from one step to the next
314  // to identify impossible hits and subsequent debugging printout
315  m_SavePreStepStatus = prePoint->GetStepStatus();
316  m_SavePostStepStatus = postPoint->GetStepStatus();
318  m_SaveVolPost = touchpost->GetVolume();
319 
320  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
321  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
322  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
323  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
324  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
325  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
326 
327  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
328  //sum up the energy to get total deposited
330  m_EdepSum += edep;
331  if (!hasMotherSubsystem() && (m_Hit->get_z(1) * cm > m_Zmax || m_Hit->get_z(1) * cm < m_Zmin))
332  {
333  std::cout << m_Detector->SuperDetector() << std::setprecision(9)
334  << " PHG4CylinderSteppingAction: Exit hit z " << m_Hit->get_z(1) * cm
335  << " outside acceptance zmin " << m_Zmin
336  << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
337  }
338  if (geantino)
339  {
340  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
341  if(layer_type){
342  m_Hit->set_eion(-1);
343  m_Hit->set_light_yield(-1);
344  }
345  }
346  else
347  {
348  if (!m_BlackHoleFlag)
349  {
350  eion = edep - aStep->GetNonIonizingEnergyDeposit() / GeV;
351  m_Hit->set_eion(m_Hit->get_eion() + eion);
353  m_Hit->set_light_yield(m_Hit->get_light_yield()+light_yield);
354  }
355  }
356 
357  if (edep > 0 || m_SaveAllHitsFlag)
358  {
360  {
361  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
362  {
363  pp->SetKeep(1); // we want to keep the track
364  }
365  }
366  }
367  // here we just update the exit values, it will be overwritten
368  // for every step until we leave the volume or the particle
369  // ceases to exist
370  // sum up the energy to get total deposited
376  // if any of these conditions is true this is the last step in
377  // this volume and we need to save the hit
378  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
379  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
380  // (happens when your detector goes outside world volume)
381  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
382  // aTrack->GetTrackStatus() == fStopAndKill is also set)
383  // aTrack->GetTrackStatus() == fStopAndKill: track ends
384  if (postPoint->GetStepStatus() == fGeomBoundary ||
385  postPoint->GetStepStatus() == fWorldBoundary ||
386  postPoint->GetStepStatus() == fAtRestDoItProc ||
387  aTrack->GetTrackStatus() == fStopAndKill ||
388  m_UseG4StepsFlag > 0)
389  {
390  // save only hits with energy deposit (or geantino)
391  if (m_Hit->get_edep() || m_SaveAllHitsFlag)
392  {
393  m_Hit->set_layer(layer_id);
394  m_Hit->set_hit_type(layer_type);
395  m_HitContainer->AddHit(layer_id, m_Hit);
396  if (m_SaveShower)
397  {
399  }
400  m_Hit = nullptr;
401  }
402  else
403  {
404  // if this hit has no energy deposit, just reset it for reuse
405  // this means we have to delete it in the dtor. If this was
406  // the last hit we processed the memory is still allocated
407  m_Hit->Reset();
408  }
409  }
410  // return true to indicate the hit was used
411  return true;
412 }
413 
414 //____________________________________________________________________________..
416 {
417  int j_0 = 0; //The j and k indices for the scintillator / tower
418  int k_0 = 0; //The j and k indices for the scintillator / tower
419  double towersize = m_Params->get_double_param("tower_size");
420  double width = m_Params->get_double_param("width");
421  double height = m_Params->get_double_param("height");
422  double det_x_pos = m_Params->get_double_param("global_x");
423  double det_y_pos = m_Params->get_double_param("global_y");
424  j_0 = (int) ( ( width - det_x_pos + towersize + ( prePoint->GetPosition().x() / cm ) ) / towersize ); //TODO DRCALO TOWER SIZE
425  k_0 = (int) ( ( height - det_y_pos + towersize + ( prePoint->GetPosition().y() / cm ) ) / towersize ); //TODO DRCALO TOWER SIZE
426  j = (j_0 * 1);
427  k = (k_0 * 1);
428 
429  return 0;
430 }
431 //____________________________________________________________________________..
433 {
434  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
435  if (!m_HitContainer)
436  {
437  std::cout << "EICG4BwdSteppingAction::SetTopNode - unable to find "
438  << m_HitNodeName << std::endl;
439  gSystem->Exit(1);
440  }
441 }
443 {
445  {
446  return true;
447  }
448  return false;
449 }