ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4B0SteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4B0SteppingAction.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 "EICG4B0SteppingAction.h"
20 
21 #include "EICG4B0Detector.h"
22 #include "EICG4B0Subsystem.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 
37 #include <TSystem.h>
38 
39 #include <Geant4/G4NavigationHistory.hh>
40 #include <Geant4/G4ParticleDefinition.hh>
41 #include <Geant4/G4ReferenceCountedHandle.hh>
42 #include <Geant4/G4Step.hh>
43 #include <Geant4/G4StepPoint.hh>
44 #include <Geant4/G4StepStatus.hh>
45 #include <Geant4/G4String.hh>
46 #include <Geant4/G4SystemOfUnits.hh>
47 #include <Geant4/G4ThreeVector.hh>
48 #include <Geant4/G4TouchableHandle.hh>
49 #include <Geant4/G4Track.hh>
50 #include <Geant4/G4TrackStatus.hh>
51 #include <Geant4/G4Types.hh>
52 #include <Geant4/G4VPhysicalVolume.hh>
53 #include <Geant4/G4VTouchable.hh>
54 #include <Geant4/G4VUserTrackInformation.hh>
55 
56 #include <cmath>
57 #include <cstdlib>
58 #include <iomanip>
59 #include <iostream>
60 #include <string>
61 
62 class PHCompositeNode;
63 
64 //____________________________________________________________________________..
66  : PHG4SteppingAction(detector->GetName())
67  , m_Subsystem(subsys)
68  , m_Detector(detector)
69  , m_Params(parameters)
70  , m_HitContainer(nullptr)
71  , m_Hit(nullptr)
72  , m_SaveShower(nullptr)
73  , m_SaveVolPre(nullptr)
74  , m_SaveVolPost(nullptr)
75  // , m_SaveLightYieldFlag(m_Params->get_int_param("lightyield"))
76  , m_SaveTrackId(-1)
77  , m_SavePreStepStatus(-1)
78  , m_SavePostStepStatus(-1)
79  , m_ActiveFlag(m_Params->get_int_param("active"))
80  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
81  , m_UseG4StepsFlag(m_Params->get_int_param("use_g4steps"))
82  , m_Zmin(m_Params->get_double_param("place_z") * cm - m_Params->get_double_param("length") * cm / 2.)
83  , m_Zmax(m_Params->get_double_param("place_z") * cm + m_Params->get_double_param("length") * cm / 2.)
84  , m_Tmin(m_Params->get_double_param("tmin") * ns)
85  , m_Tmax(m_Params->get_double_param("tmax") * ns)
86  , m_EdepSum(0)
87  , m_EabsSum(0)
88  , _towerdivision(0.0)
89  , _tower_size(2.0)
90  , _readout_size(2.0)
91  , _detector_size(20)
92 {
93  // G4 seems to have issues in the um range
94  m_Zmin -= copysign(m_Zmin, 1. / 1e6 * cm);
95  m_Zmax += copysign(m_Zmax, 1. / 1e6 * cm);
96 }
97 
98 //____________________________________________________________________________..
100 {
101  // if the last hit was a zero energie deposit hit, it is just reset
102  // and the memory is still allocated, so we need to delete it here
103  // if the last hit was saved, hit is a nullptr pointer which are
104  // legal to delete (it results in a no operation)
105  delete m_Hit;
106 }
107 
108 //____________________________________________________________________________..
109 // This is the implementation of the G4 UserSteppingAction
110 bool EICG4B0SteppingAction::UserSteppingAction(const G4Step *aStep, bool was_used)
111 {
113  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
114  // get volume of the current step
115  G4VPhysicalVolume *volume = touch->GetVolume();
116  // IsInDetector(volume) returns
117  // == 0 outside of detector
118  // > 0 for hits in active volume
119  // < 0 for hits in passive material
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  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
127  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
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 = 0;
150  if (m_Params->get_string_param("material") == "G4_Cu")
151  layer_type = 0;
152  else
153  layer_type = 1;
154  //if (layer_id % 2 == 1) //Proper implementation for several layer configurations
155  /*
156  if (m_Params->get_string_param("material") == "G4_Cu")
157  {
158  return false;
159  }
160 */
161  if (!m_ActiveFlag)
162  {
163  return false;
164  }
165  bool geantino = false;
166  // the check for the pdg code speeds things up, I do not want to make
167  // an expensive string compare for every track when we know
168  // geantino or chargedgeantino has pid=0
169  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
170  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") !=
171  std::string::npos) // this also accounts for "chargedgeantino"
172  {
173  geantino = true;
174  }
175  G4StepPoint *prePoint = aStep->GetPreStepPoint();
176  G4StepPoint *postPoint = aStep->GetPostStepPoint();
177 
178  // Here we have to decide if we need to create a new hit. Normally this should
179  // only be neccessary if a G4 Track enters a new volume or is freshly created
180  // For this we look at the step status of the prePoint (beginning of the G4 Step).
181  // This should be either fGeomBoundary (G4 Track crosses into volume) or
182  // fUndefined (G4 Track newly created)
183  // Sadly over the years with different G4 versions we have observed cases where
184  // G4 produces "impossible hits" which we try to catch here
185  // These errors were always rare and it is not clear if they still exist but we
186  // still check for them for safety. We can reproduce G4 runs identically (if given
187  // the sequence of random number seeds you find in the log), the printouts help
188  // us giving the G4 support information about those failures
189  //
190  switch (prePoint->GetStepStatus())
191  {
192  case fPostStepDoItProc:
194  {
195  // this is the okay case, fPostStepDoItProc called in a volume, not first thing inside
196  // a new volume, just proceed here
197  break;
198  }
199  else
200  {
201  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
202  // this is still with us
203  std::cout << GetName() << ": New Hit for " << std::endl;
204  std::cout << "prestep status: "
206  << ", poststep status: "
208  << ", last pre step status: "
210  << ", last post step status: "
212  std::cout << "last track: " << m_SaveTrackId
213  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
214  std::cout << "phys pre vol: " << volume->GetName()
215  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
216  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
217  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
218  }
219  break;
220  // These are the normal cases
221  case fGeomBoundary:
222  case fUndefined:
223  if (!m_Hit)
224  {
225  m_Hit = new PHG4Hitv1();
226  }
227  m_Hit->set_layer((unsigned int) layer_id);
228  // here we set the entrance values in cm
229  FindTowerIndexFromPosition(prePoint, idx_j, idx_k);
230  //std::cout << "B0 Hits: "<<postPoint->GetPosition().x() / cm<<" "<<postPoint->GetPosition().y() / cm<<" "<<edep<<std::endl;
231 
232  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
233  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
234  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
235 
236  m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
237  m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
238  m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
239 
240  // time in ns
241  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
242  // set the track ID
243  m_Hit->set_trkid(aTrack->GetTrackID());
244  m_SaveTrackId = aTrack->GetTrackID();
245  // set the initial energy deposit
246  m_EdepSum = 0;
247  m_EabsSum = 0;
248 
249  m_Hit->set_edep(0);
250  if (!geantino && !m_BlackHoleFlag)
251  {
252  m_Hit->set_eion(0);
253  }
254  if (layer_type)
255  {
256  // m_SaveHitContainer=m_HitContainer;
257  m_Hit->set_eion(0);
258  // if (m_SaveLightYieldFlag)
259  // {
260  // m_Hit->set_light_yield(0);
261  m_Hit->set_index_j(idx_j);
262  m_Hit->set_index_k(idx_k);
263  // }
264  }
265  /* else{
266  m_SaveHitContainer=m_AbsorberHitContainer;
267  }*/
268  // implement your own here://
269  // add the properties you are interested in via set_XXX methods
270  // you can find existing set methods in $OFFLINE_MAIN/include/g4main/PHG4Hit.h
271  // this is initialization of your value. This is not needed you can just set the final
272  // value at the last step in this volume later one
273  /* if (whichactive > 0)
274  {
275  m_EionSum = 0; // assuming the ionization energy is only needed for active
276  // volumes (scintillators)
277  m_Hit->set_eion(0);
278  m_SaveHitContainer = m_HitContainer;
279  }
280  else
281  {
282  std::cout << "implement stuff for whichactive < 0 (inactive volumes)" << std::endl;
283  gSystem->Exit(1);
284  }
285 */
286  // this is for the tracking of the truth info
288  {
289  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
290  {
291  m_Hit->set_trkid(pp->GetUserTrackId());
292  m_Hit->set_shower_id(pp->GetShower()->get_id());
293  m_SaveShower = pp->GetShower();
294  // pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
295  }
296  }
297  if (!hasMotherSubsystem() && (m_Hit->get_z(0) * cm > m_Zmax || m_Hit->get_z(0) * cm < m_Zmin))
298  {
299  std::cout << m_Detector->SuperDetector() << std::setprecision(9)
300  << " PHG4CylinderSteppingAction: Entry hit z " << m_Hit->get_z(0) * cm
301  << " outside acceptance, zmin " << m_Zmin
302  << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
303  }
304  break;
305  default:
306  break;
307  }
308 
309  // This section is called for every step
310  // some sanity checks for inconsistencies (aka bugs) we have seen over the years
311  // check if this hit was created, if not print out last post step status
312  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
313  {
314  std::cout << GetName() << ": hit was not created" << std::endl;
315  std::cout << "prestep status: "
317  << ", poststep status: "
319  << ", last pre step status: "
321  << ", last post step status: "
323  std::cout << "last track: " << m_SaveTrackId
324  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
325  std::cout << "phys pre vol: " << volume->GetName()
326  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
327  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
328  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
329  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
330  gSystem->Exit(1);
331  }
332  m_SavePostStepStatus = postPoint->GetStepStatus();
333  // check if track id matches the initial one when the hit was created
334  if (aTrack->GetTrackID() != m_SaveTrackId)
335  {
336  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
337  std::cout << "saved track: " << m_SaveTrackId
338  << ", current trackid: " << aTrack->GetTrackID()
339  << ", prestep status: " << prePoint->GetStepStatus()
340  << ", previous post step status: " << m_SavePostStepStatus << std::endl;
341  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
342  gSystem->Exit(1);
343  }
344 
345  // We need to cache a few things from one step to the next
346  // to identify impossible hits and subsequent debugging printout
347  m_SavePreStepStatus = prePoint->GetStepStatus();
348  m_SavePostStepStatus = postPoint->GetStepStatus();
350  m_SaveVolPost = touchpost->GetVolume();
351 
352  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
353  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
354  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
355 
356  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
357  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
358  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
359 
360  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
361  //sum up the energy to get total deposited
363  if (layer_type)
364  {
365  m_Hit->set_eion(m_Hit->get_eion() + eion);
367  }
368  if (!layer_type) m_EabsSum += edep;
369  m_EdepSum += edep;
370  if (!hasMotherSubsystem() && (m_Hit->get_z(1) * cm > m_Zmax || m_Hit->get_z(1) * cm < m_Zmin))
371  {
372  std::cout << m_Detector->SuperDetector() << std::setprecision(9)
373  << " PHG4CylinderSteppingAction: Exit hit z " << m_Hit->get_z(1) * cm
374  << " outside acceptance zmin " << m_Zmin
375  << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
376  }
377  if (geantino)
378  {
379  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
380  if (layer_type)
381  {
382  m_Hit->set_eion(-1);
383  // m_Hit->set_light_yield(-1);
384  }
385  }
386  else
387  {
388  if (!m_BlackHoleFlag)
389  {
390  double eion_2 = edep - aStep->GetNonIonizingEnergyDeposit() / GeV;
391  m_Hit->set_eion(m_Hit->get_eion() + eion_2);
392  }
393  }
394 
395  /* if (m_SaveLightYieldFlag)
396  {
397  double light_yield = GetVisibleEnergyDeposition(aStep) / GeV;
398  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
399  }
400 */
401  if (edep > 0 || m_SaveAllHitsFlag)
402  {
404  {
405  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
406  {
407  pp->SetKeep(1); // we want to keep the track
408  }
409  }
410  }
411  // here we just update the exit values, it will be overwritten
412  // for every step until we leave the volume or the particle
413  // ceases to exist
414  // sum up the energy to get total deposited
420  // if any of these conditions is true this is the last step in
421  // this volume and we need to save the hit
422  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
423  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
424  // (happens when your detector goes outside world volume)
425  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
426  // aTrack->GetTrackStatus() == fStopAndKill is also set)
427  // aTrack->GetTrackStatus() == fStopAndKill: track ends
428  if (postPoint->GetStepStatus() == fGeomBoundary ||
429  postPoint->GetStepStatus() == fWorldBoundary ||
430  postPoint->GetStepStatus() == fAtRestDoItProc ||
431  aTrack->GetTrackStatus() == fStopAndKill ||
432  m_UseG4StepsFlag > 0)
433  {
434  // save only hits with energy deposit (or geantino)
435  if (m_Hit->get_edep() || m_SaveAllHitsFlag)
436  {
437  m_Hit->set_layer(layer_id);
438  m_Hit->set_hit_type(layer_type);
439  // m_Hit->set_eion(m_EionSum);
440  // m_Hit->set_edep(m_EdepSum);
441  m_HitContainer->AddHit(layer_id, m_Hit);
442  if (m_SaveShower)
443  {
445  }
446  m_Hit = nullptr;
447  }
448  else
449  {
450  // if this hit has no energy deposit, just reset it for reuse
451  // this means we have to delete it in the dtor. If this was
452  // the last hit we processed the memory is still allocated
453  m_Hit->Reset();
454  }
455  }
456  // return true to indicate the hit was used
457  return true;
458 }
459 
460 //____________________________________________________________________________..
462 {
463  int j_0 = 0; //The j and k indices for the scintillator / tower
464  int k_0 = 0; //The j and k indices for the scintillator / tower
465  int towersize = 2;
466  if (_towerdivision == 0.)
467  {
468  int maxsubtow = (int) ((_tower_size) / (_readout_size));
469  _towerdivision = (_tower_size - (maxsubtow * _readout_size)) / maxsubtow;
471  }
472  // j_0 = (int) ( ( _detector_size + ( prePoint->GetPosition().x() ) ) / _towerdivision ); //TODO DRCALO TOWER SIZE
473  // k_0 = (int) ( ( _detector_size + ( prePoint->GetPosition().y() ) ) / _towerdivision ); //TODO DRCALO TOWER SIZE
474 
475  j_0 = (int) ((36.6 + (prePoint->GetPosition().x() / cm)) / towersize); //TODO DRCALO TOWER SIZE
476  k_0 = (int) ((22 + (prePoint->GetPosition().y() / cm)) / towersize); //TODO DRCALO TOWER SIZE
477  j = (j_0 * 1);
478  k = (k_0 * 1);
479 
480  return 0;
481 }
482 //____________________________________________________________________________..
484 {
485  //std::string hitnodename = "G4HIT_" + m_Detector->GetName();
486  // std::cout << " ---> !!! hitnodename: " << hitnodename << std::endl;
487  // now look for the map and grab a pointer to it.
488  // std::string hitnodename;
489  // std::string absorbernodename;
490  // if (m_Detector->SuperDetector()!= "NONE")
491  // {
492  // hitnodename = "G4HIT_" + m_Detector->SuperDetector();
493  // absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
494  // }
495  // else
496  // {
497  // hitnodename = "G4HIT_" + m_Detector->GetName();
498  // absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
499  // }
500  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
501  // m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
502  // if we do not find the node we need to make it.
503  if (!m_HitContainer)
504  {
505  std::cout << "EICG4B0SteppingAction::SetTopNode - unable to find "
506  << m_HitNodeName << std::endl;
507  gSystem->Exit(1);
508  }
509  /* if (!m_AbsorberHitContainer)
510  {
511  if (Verbosity() > 0)
512  {
513  std::cout << "EICG4B0SteppingAction::SetInterfacePointers - unable to find " << absorbernodename << std::endl;
514  std::cout <<"Are you running Realistic B0 ???"<<std::endl;
515  }
516  }*/
517 }
519 {
521  {
522  return true;
523  }
524  return false;
525 }