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