ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4ZDCSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4ZDCSteppingAction.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 // -1/June/2021 First ZDC design (Crystal + FoCal style) Shima Shimizu
20 // Attach Layer and channel ID to G4Hits for TTree/Ntuple production
21 // -14/Dec/2021 Second ZDC design. Modifications for digitization.
22 
23 #include "EICG4ZDCSteppingAction.h"
24 
25 #include "EICG4ZDCDetector.h"
26 #include "EICG4ZDCdetid.h"
27 
28 #include <phparameter/PHParameters.h>
29 
31 
32 #include <g4main/PHG4Hit.h>
34 #include <g4main/PHG4Hitv1.h>
35 #include <g4main/PHG4Shower.h>
38 
39 #include <phool/getClass.h>
40 
41 #include <TSystem.h>
42 
43 #include <Geant4/G4ParticleDefinition.hh>
44 #include <Geant4/G4ReferenceCountedHandle.hh>
45 #include <Geant4/G4Step.hh>
46 #include <Geant4/G4StepPoint.hh>
47 #include <Geant4/G4StepStatus.hh>
48 #include <Geant4/G4String.hh>
49 #include <Geant4/G4SystemOfUnits.hh>
50 #include <Geant4/G4ThreeVector.hh>
51 #include <Geant4/G4TouchableHandle.hh>
52 #include <Geant4/G4Track.hh>
53 #include <Geant4/G4TrackStatus.hh>
54 #include <Geant4/G4Types.hh>
55 #include <Geant4/G4VPhysicalVolume.hh>
56 #include <Geant4/G4VTouchable.hh>
57 #include <Geant4/G4VUserTrackInformation.hh>
58 
59 #include <cmath>
60 #include <iostream>
61 #include <string>
62 
63 class PHCompositeNode;
64 
65 //____________________________________________________________________________..
67  : PHG4SteppingAction(detector->GetName())
68  , m_Detector(detector)
69  , m_Params(parameters)
70  , m_HitContainer(nullptr)
71  , m_Hit(nullptr)
72  , m_SaveHitContainer(nullptr)
73  , m_SaveVolPre(nullptr)
74  , m_SaveVolPost(nullptr)
75  , m_SaveShower(nullptr)
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_EdepSum(0)
82  , m_EionSum(0)
83  , m_LightYield(0)
84 {
85 
86 }
87 
88 //____________________________________________________________________________..
90 {
91  // if the last hit was a zero energie deposit hit, it is just reset
92  // and the memory is still allocated, so we need to delete it here
93  // if the last hit was saved, hit is a nullptr pointer which are
94  // legal to delete (it results in a no operation)
95  delete m_Hit;
96 }
97 
98 //____________________________________________________________________________..
99 // This is the implementation of the G4 UserSteppingAction
100 bool EICG4ZDCSteppingAction::UserSteppingAction(const G4Step *aStep,bool was_used)
101 {
103  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
104  // get volume of the current step
105  G4VPhysicalVolume *volume = touch->GetVolume();
106 
107  // IsInDetector(volume) returns
108  // == 0 outside of detector
109  // > 0 for hits in active volume
110  // < 0 for hits in passive material
111  int whichactive = m_Detector->IsInDetector(volume);
112  if (!whichactive)
113  {
114  return false;
115  }
116 
117 
118  // collect energy and track length step by step
119  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
120  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
121  G4double light_yield = 0;
122  if(whichactive>0) light_yield = GetVisibleEnergyDeposition(aStep);
123  const G4Track *aTrack = aStep->GetTrack();
124  // if this detector stops everything, just put all kinetic energy into edep
125  if (m_BlackHoleFlag)
126  {
127  edep = aTrack->GetKineticEnergy() / GeV;
128  G4Track *killtrack = const_cast<G4Track *>(aTrack);
129  killtrack->SetTrackStatus(fStopAndKill);
130  }
131  // we use here only one detector in this simple example
132  // if you deal with multiple detectors in this stepping action
133  // the detector id can be used to distinguish between them
134  // hits can easily be analyzed later according to their detector id
135 
136  /*------------------------------------------------*/
137  /*--Here check the ZDC detector ID and layers-----*/
138  /*------------------------------------------------*/
139 
140  int detflag = -1;
141  if(whichactive>0) detflag = m_Detector->GetActiveVolumeInfo(volume);
142  else if (whichactive<0) detflag = m_Detector->GetAbsorberVolumeInfo(volume);
143  int detector_id = detflag%100;
144  int detector_layer = (detflag%10000)/100;
145  int detector_nlyrbox =(detflag%1000000)/10000;
146  int detector_system= (detflag/1000000)*1000000;
147 
148  int xid = touch->GetCopyNumber(1);
149  int yid = touch->GetCopyNumber();
150 
151  int layer_id = -1;
152 
153  if(whichactive>0){
154  if(detector_system == ZDCID::CrystalTower){
155  int zid = touch->GetCopyNumber(2);
156  if(detector_id == ZDCID::SI_PIXEL) layer_id = zid * 2;
157  else if(detector_id==ZDCID::Crystal) layer_id = zid * 2 +1;
158 
159  }else if(detector_system == ZDCID::EMLayer){
160  if(detector_id == ZDCID::SI_PIXEL) layer_id = detector_layer + touch->GetCopyNumber(3);
161  if(detector_id == ZDCID::SI_PAD) {
162  int boxid = touch->GetCopyNumber(4);
163  int zid =touch->GetCopyNumber(3);
164  int nlyr = detector_nlyrbox +1;
165  layer_id = detector_layer + zid + boxid * nlyr;
166  }
167  }else if (detector_system == ZDCID::HCPadLayer){
168  if(detector_id == ZDCID::SI_PAD) layer_id = detector_layer + touch->GetCopyNumber(3);
169  }else if (detector_system == ZDCID::HCSciLayer){
170  if(detector_id == ZDCID::Scintillator){
171  int boxid = touch->GetCopyNumber(4);
172  int zid = touch->GetCopyNumber(3);
173  int nlyr = detector_nlyrbox;
174  layer_id = detector_layer + zid + boxid *nlyr;
175  }
176  }
177  }
178 
179  if(whichactive<0){
180  layer_id = detector_layer;
181  if(detector_system == ZDCID::HCSciLayer){
182  int boxid = touch->GetCopyNumber(2);
183  layer_id = detector_layer + boxid * detector_nlyrbox;
184  }
185  }
186 
187  bool geantino = false;
188  // the check for the pdg code speeds things up, I do not want to make
189  // an expensive string compare for every track when we know
190  // geantino or chargedgeantino has pid=0
191  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
192  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") !=
193  std::string::npos) // this also accounts for "chargedgeantino"
194  {
195  geantino = true;
196  }
197  G4StepPoint *prePoint = aStep->GetPreStepPoint();
198  G4StepPoint *postPoint = aStep->GetPostStepPoint();
199 
200 // Here we have to decide if we need to create a new hit. Normally this should
201 // only be neccessary if a G4 Track enters a new volume or is freshly created
202 // For this we look at the step status of the prePoint (beginning of the G4 Step).
203 // This should be either fGeomBoundary (G4 Track crosses into volume) or
204 // fUndefined (G4 Track newly created)
205 // Sadly over the years with different G4 versions we have observed cases where
206 // G4 produces "impossible hits" which we try to catch here
207 // These errors were always rare and it is not clear if they still exist but we
208 // still check for them for safety. We can reproduce G4 runs identically (if given
209 // the sequence of random number seeds you find in the log), the printouts help
210 // us giving the G4 support information about those failures
211 //
212  switch (prePoint->GetStepStatus()){
213  case fPostStepDoItProc:
215  // this is the okay case, fPostStepDoItProc called in a volume, not first thing inside
216  // a new volume, just proceed here
217  break;
218  } else {
219  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
220  // this is still with us
221  std::cout << GetName() << ": New Hit for " << std::endl;
222  std::cout << "prestep status: "
224  << ", poststep status: "
226  << ", last pre step status: "
228  << ", last post step status: "
230  std::cout << "last track: " << m_SaveTrackId
231  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
232  std::cout << "phys pre vol: " << volume->GetName()
233  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
234  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
235  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
236  }
237 // These are the normal cases
238  case fGeomBoundary:
239  case fUndefined:
240  if (!m_Hit) {
241  m_Hit = new PHG4Hitv1();
242  }
243  // here we set the entrance values in cm
244  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
245  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
246  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
247  // time in ns
248  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
249  // set the track ID
250  m_Hit->set_trkid(aTrack->GetTrackID());
251  m_SaveTrackId = aTrack->GetTrackID();
252  // set the initial energy deposit
253  m_EdepSum = 0;
254  // implement your own here://
255  // add the properties you are interested in via set_XXX methods
256  // you can find existing set methods in $OFFLINE_MAIN/include/g4main/PHG4Hit.h
257  // this is initialization of your value. This is not needed you can just set the final
258  // value at the last step in this volume later one
259 
260  if (whichactive > 0)
261  {
262  m_EionSum = 0; // assuming the ionization energy is only needed for active
263  // volumes (scintillators)
264  m_LightYield = 0;
265  m_Hit->set_eion(0);
266  m_Hit->set_light_yield(0);;
268  }
269  else
270  {
272  // std::cout << "implement stuff for whichactive < 0 (inactive volumes)" << std::endl;
273  // gSystem->Exit(1);
274  }
275  // this is for the tracking of the truth info
277  {
278  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
279  {
280  m_Hit->set_trkid(pp->GetUserTrackId());
281  m_Hit->set_shower_id(pp->GetShower()->get_id());
282  m_SaveShower=pp->GetShower();
283  }
284  }
285  break;
286  default:
287  break;
288  }
289 
290  // This section is called for every step
291  // some sanity checks for inconsistencies (aka bugs) we have seen over the years
292  // check if this hit was created, if not print out last post step status
293  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
294  {
295  std::cout << GetName() << ": hit was not created" << std::endl;
296  std::cout << "prestep status: "
298  << ", poststep status: "
300  << ", last pre step status: "
302  << ", last post step status: "
304  std::cout << "last track: " << m_SaveTrackId
305  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
306  std::cout << "phys pre vol: " << volume->GetName()
307  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
308  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
309  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
310  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
311  gSystem->Exit(1);
312  }
313  // check if track id matches the initial one when the hit was created
314  if (aTrack->GetTrackID() != m_SaveTrackId){
315  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
316  std::cout << "saved track: " << m_SaveTrackId
317  << ", current trackid: " << aTrack->GetTrackID()
318  << ", prestep status: " << prePoint->GetStepStatus()
319  << ", previous post step status: " << m_SavePostStepStatus << std::endl;
320  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
321  gSystem->Exit(1);
322  }
323 
324 // We need to cache a few things from one step to the next
325 // to identify impossible hits and subsequent debugging printout
326  m_SavePreStepStatus = prePoint->GetStepStatus();
327  m_SavePostStepStatus = postPoint->GetStepStatus();
329  m_SaveVolPost = touchpost->GetVolume();
330  // here we just update the exit values, it will be overwritten
331  // for every step until we leave the volume or the particle
332  // ceases to exist
333  // sum up the energy to get total deposited
334  m_EdepSum += edep;
335  if (whichactive > 0) {
336  m_EionSum += eion;
337  m_LightYield += light_yield;
338  }
339 
340  // if any of these conditions is true this is the last step in
341  // this volume and we need to save the hit
342  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
343  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
344  // (happens when your detector goes outside world volume)
345  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
346  // aTrack->GetTrackStatus() == fStopAndKill is also set)
347  // aTrack->GetTrackStatus() == fStopAndKill: track ends
348 
349  if (postPoint->GetStepStatus() == fGeomBoundary ||
350  postPoint->GetStepStatus() == fWorldBoundary ||
351  postPoint->GetStepStatus() == fAtRestDoItProc ||
352  aTrack->GetTrackStatus() == fStopAndKill){
353  // save only hits with energy deposit (or geantino)
354  if (m_EdepSum > 0 || geantino){
355  // update values at exit coordinates and set keep flag
356  // of track to keep
357  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
358  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
359  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
360  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
361  if (G4VUserTrackInformation *p = aTrack->GetUserInformation()){
362  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p)){
363  pp->SetKeep(1); // we want to keep the track
364  }
365  }
366 
367  if (geantino){
368  //implement your own here://
369  // if you want to do something special for geantinos (normally you do not)
370  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way
371  // geantinos survive the g4hit compression
372  if (whichactive > 0){
373  m_Hit->set_eion(-1);
374  m_Hit->set_light_yield(-1);
375  }
376  } else {
378  }
379  //implement your own here://
380  // what you set here will be saved in the output
381  m_Hit->set_layer(layer_id);
382  m_Hit->set_index_i(xid);
383  m_Hit->set_index_j(yid);
384  m_Hit->set_hit_type(detector_id);
387 
388  m_SaveHitContainer->AddHit(detector_id, m_Hit);
389  if (m_SaveShower)
390  {
392  }
393 
394  // ownership has been transferred to container, set to null
395  // so we will create a new hit for the next track
396  m_Hit = nullptr;
397 
398  }else{
399  // if this hit has no energy deposit, just reset it for reuse
400  // this means we have to delete it in the dtor. If this was
401  // the last hit we processed the memory is still allocated
402  m_Hit->Reset();
403  }
404  }
405  // return true to indicate the hit was used
406  return true;
407 }
408 
409 //____________________________________________________________________________..
411 {
412  std::string hitnodename = "G4HIT_" + m_Detector->GetName();
413  // now look for the map and grab a pointer to it.
414  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
415  // if we do not find the node we need to make it.
416  if (!m_HitContainer)
417  {
418  std::cout << "EICG4ZDCSteppingAction::SetTopNode - unable to find "
419  << hitnodename << std::endl;
420  }
421 }