ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4LFHcalSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4LFHcalSteppingAction.cc
2 
3 #include "PHG4LFHcalDetector.h"
4 
5 #include <phparameter/PHParameters.h>
6 
7 #include <g4main/PHG4Hit.h>
9 #include <g4main/PHG4Hitv1.h>
10 #include <g4main/PHG4Shower.h>
11 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
12 
14 
15 #include <phool/getClass.h>
16 
17 #include <Geant4/G4IonisParamMat.hh> // for G4IonisParamMat
18 #include <Geant4/G4Material.hh> // for G4Material
19 #include <Geant4/G4MaterialCutsCouple.hh>
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/G4SystemOfUnits.hh>
26 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
27 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
28 #include <Geant4/G4Track.hh> // for G4Track
29 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
30 #include <Geant4/G4Types.hh> // for G4double
31 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
32 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
33 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
34 
35 #include <iostream>
36 #include <string> // for basic_string, operator+
37 
38 class PHCompositeNode;
39 
40 //____________________________________________________________________________..
42  : PHG4SteppingAction(detector->GetName())
43  , m_Detector(detector)
44  , m_ActiveFlag(parameters->get_int_param("active"))
45  , m_AbsorberTruthFlag(parameters->get_int_param("absorberactive"))
46  , m_BlackHoleFlag(parameters->get_int_param("blackhole"))
47  , m_NlayersPerTowerSeg(parameters->get_int_param("nlayerspertowerseg"))
48  , m_zdepthtailcatcher(parameters->get_double_param("zdepthcatcheroffset"))
49  , m_NLayerTailCatcherOffset(parameters->get_int_param("nLayerOffsetTailcatcher"))
50 {
51 }
52 
54 {
55  // if the last hit was a zero energie deposit hit, it is just reset
56  // and the memory is still allocated, so we need to delete it here
57  // if the last hit was saved, hit is a nullptr pointer which are
58  // legal to delete (it results in a no operation)
59  delete m_Hit;
60 }
61 
62 //____________________________________________________________________________..
64 {
67  // m_Detector->IsInLFHcal(volume)
68  // returns
69  // 0 is outside of Forward HCAL
70  // 1 is inside scintillator
71  // -1 is inside absorber (dead material)
72 
73  int whichactive = m_Detector->IsInLFHcal(volume);
74 
75  if (!whichactive)
76  {
77  return false;
78  }
79 
80  int layer_id = m_Detector->get_Layer();
81  unsigned int icopy = touch->GetVolume(3)->GetCopyNo();
82  unsigned int layer = touch->GetVolume(1)->GetCopyNo();
83  int idx_j = icopy >> 16;
84  int idx_k = icopy & 0xFFFF;
85  int idx_l = (int) (layer / m_NlayersPerTowerSeg);
86 
87  if (Verbosity() > 2)
88  std::cout << "\t" << icopy << "\t idx_j =" << idx_j << ", idx_k =" << idx_k << ", idx_l =" << idx_l << "\t id:" << icopy << "\t layer:" << layer << std::endl;
89 
90  /* Get energy deposited by this step */
91  double edep = aStep->GetTotalEnergyDeposit() / GeV;
92  double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
93 
94  /* Get pointer to associated Geant4 track */
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  /* Make sure we are in a volume */
106  if (m_ActiveFlag)
107  {
108  /* Check if particle is 'geantino' */
109  bool geantino = false;
110  double light_yield = 0;
111  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
112  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
113  {
114  geantino = true;
115  }
116 
117  /* Get Geant4 pre- and post-step points */
118  G4StepPoint* prePoint = aStep->GetPreStepPoint();
119  G4StepPoint* postPoint = aStep->GetPostStepPoint();
120  // std::cout << "\t" << icopy << "\t idx_j =" << idx_j << ", idx_k =" << idx_k << ", idx_l =" << idx_l <<"\t id:" << icopy << "\t layer_id:" << layer_id << "\t layer:" << layer << "\tz: " << printf("%3.0f",prePoint->GetPosition().z() / cm) << "\ttco: " << m_zdepthtailcatcher << "\tnlayoffs: " << m_NLayerTailCatcherOffset;
121  // if(prePoint->GetPosition().z() / cm > m_zdepthtailcatcher){
122  // idx_l = (int)((layer + m_NLayerTailCatcherOffset)/ m_NlayersPerTowerSeg);
123  // std::cout << "\tidx_l fixed: " << idx_l << std::endl;
124  // }
125  // else { std::cout <<std::endl;}
126  // std::cout << "x: " << prePoint->GetPosition().x() << "\ty: " << prePoint->GetPosition().y() << "\tz: " << prePoint->GetPosition().z() << std::endl;
127 
128  switch (prePoint->GetStepStatus())
129  {
130  case fGeomBoundary:
131  case fUndefined:
132  if (!m_Hit)
133  {
134  m_Hit = new PHG4Hitv1();
135  }
136 
137  /* Set hit location (space point) */
138  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
139  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
140  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
141 
142  /* Set hit time */
143  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
144 
145  // set the tower index
146  m_Hit->set_index_j(idx_j);
147  m_Hit->set_index_k(idx_k);
148  m_Hit->set_index_l(idx_l);
149 
150  //set the track ID
151  m_Hit->set_trkid(aTrack->GetTrackID());
152  /* set intial energy deposit */
153  m_Hit->set_edep(0);
154  /* Set hit location (tower index) */
155 
156  /* Now add the hit to the hit collection */
157  if (whichactive > 0)
158  {
159  m_Hit->set_eion(0);
162  }
163  else
164  {
166  }
168  {
169  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
170  {
171  m_Hit->set_trkid(pp->GetUserTrackId());
172  m_Hit->set_shower_id(pp->GetShower()->get_id());
173  m_SaveShower = pp->GetShower();
174  }
175  }
176  break;
177  default:
178  break;
179  }
180 
181  if (whichactive > 0)
182  {
183  light_yield = GetVisibleEnergyDeposition(aStep); // for scintillator only, calculate light yields
184  static bool once = true;
185  if (once && edep > 0)
186  {
187  once = false;
188 
189  if (Verbosity() > 0)
190  {
191  std::cout << "PHG4LFHcalSteppingAction::UserSteppingAction::"
192  //
193  << m_Detector->GetName() << " - "
194  << " use scintillating light model at each Geant4 steps. "
195  << "First step: "
196  << "Material = "
197  << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
198  << ", "
199  << "Birk Constant = "
201  << ","
202  << "edep = " << edep << ", "
203  << "eion = " << eion
204  << ", "
205  << "light_yield = " << light_yield << std::endl;
206  }
207  }
208  }
209 
210  /* Update exit values- will be overwritten with every step until
211  * we leave the volume or the particle ceases to exist */
212  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
213  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
214  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
215 
216  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
217 
218  /* sum up the energy to get total deposited */
220  if (whichactive > 0)
221  {
222  m_Hit->set_eion(m_Hit->get_eion() + eion);
223  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
224  }
225 
226  if (geantino)
227  {
228  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
229  if (whichactive > 0)
230  {
231  m_Hit->set_eion(-1);
232  m_Hit->set_light_yield(-1);
233  }
234  }
235  if (edep > 0 && (whichactive > 0 || m_AbsorberTruthFlag > 0))
236  {
238  {
239  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
240  {
241  pp->SetKeep(1); // we want to keep the track
242  }
243  }
244  }
245  // if any of these conditions is true this is the last step in
246  // this volume and we need to save the hit
247  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
248  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
249  // (not sure if this will ever be the case)
250  // aTrack->GetTrackStatus() == fStopAndKill: track ends
251  if (postPoint->GetStepStatus() == fGeomBoundary ||
252  postPoint->GetStepStatus() == fWorldBoundary ||
253  postPoint->GetStepStatus() == fAtRestDoItProc ||
254  aTrack->GetTrackStatus() == fStopAndKill)
255  {
256  // save only hits with energy deposit (or -1 for geantino)
257  if (m_Hit->get_edep())
258  {
259  m_SaveHitContainer->AddHit(layer_id, m_Hit);
260  if (m_SaveShower)
261  {
263  }
264  // ownership has been transferred to container, set to null
265  // so we will create a new hit for the next track
266  m_Hit = nullptr;
267  }
268  else
269  {
270  // if this hit has no energy deposit, just reset it for reuse
271  // this means we have to delete it in the dtor. If this was
272  // the last hit we processed the memory is still allocated
273  m_Hit->Reset();
274  }
275  }
276 
277  return true;
278  }
279  else
280  {
281  return false;
282  }
283 }
284 
285 //____________________________________________________________________________..
287 {
288  std::string hitnodename;
289  std::string absorbernodename;
290 
291  if (m_Detector->SuperDetector() != "NONE")
292  {
293  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
294  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
295  }
296  else
297  {
298  hitnodename = "G4HIT_" + m_Detector->GetName();
299  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
300  }
301 
302  //now look for the map and grab a pointer to it.
303  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
304  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
305 
306  // if we do not find the node it's messed up.
307  if (!m_HitContainer)
308  {
309  std::cout << "PHG4LFHcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
310  }
311  if (!m_AbsorberHitContainer)
312  {
313  if (Verbosity() > 0)
314  {
315  std::cout << "PHG4LFHcalSteppingAction::SetTopNode - unable to find " << absorbernodename << std::endl;
316  }
317  }
318 }