ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ForwardHcalSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ForwardHcalSteppingAction.cc
2 
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_SupportTruthFlag(parameters->get_int_param("supportactive"))
47  , m_BlackHoleFlag(parameters->get_int_param("blackhole"))
48 {
49 }
50 
52 {
53  // if the last hit was a zero energie deposit hit, it is just reset
54  // and the memory is still allocated, so we need to delete it here
55  // if the last hit was saved, hit is a nullptr pointer which are
56  // legal to delete (it results in a no operation)
57  delete m_Hit;
58 }
59 
60 //____________________________________________________________________________..
62 {
65 
66  // m_Detector->IsInForwardHcal(volume)
67  // returns
68  // 0 is outside of Forward HCAL
69  // 1 is inside scintillator
70  // -1 is inside absorber (dead material)
71  // -2 is inside the support (or other volume)
72 
73  int whichactive = m_Detector->IsInForwardHcal(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(1)->GetCopyNo();
82  int idx_j = icopy >> 16;
83  int idx_k = icopy & 0xFFFF;
84 
85  /* Get energy deposited by this step */
86  double edep = aStep->GetTotalEnergyDeposit() / GeV;
87  double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
88 
89  /* Get pointer to associated Geant4 track */
90  const G4Track* aTrack = aStep->GetTrack();
91 
92  // if this block stops everything, just put all kinetic energy into edep
93  if (m_BlackHoleFlag)
94  {
95  edep = aTrack->GetKineticEnergy() / GeV;
96  G4Track* killtrack = const_cast<G4Track*>(aTrack);
97  killtrack->SetTrackStatus(fStopAndKill);
98  }
99 
100  /* Make sure we are in a volume */
101  if (m_ActiveFlag)
102  {
103  /* Check if particle is 'geantino' */
104  bool geantino = false;
105  double light_yield = 0;
106  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
107  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
108  {
109  geantino = true;
110  }
111 
112  /* Get Geant4 pre- and post-step points */
113  G4StepPoint* prePoint = aStep->GetPreStepPoint();
114  G4StepPoint* postPoint = aStep->GetPostStepPoint();
115 
116  switch (prePoint->GetStepStatus())
117  {
118  case fGeomBoundary:
119  case fUndefined:
120  if (!m_Hit)
121  {
122  m_Hit = new PHG4Hitv1();
123  }
124 
125  /* Set hit location (space point) */
126  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
127  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
128  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
129 
130  /* Set hit time */
131  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
132 
133  // set the tower index
134  m_Hit->set_index_j(idx_j);
135  m_Hit->set_index_k(idx_k);
136 
137  //set the track ID
138  m_Hit->set_trkid(aTrack->GetTrackID());
139  /* set intial energy deposit */
140  m_Hit->set_edep(0);
141  /* Set hit location (tower index) */
142 
143  /* Now add the hit to the hit collection */
144  if (whichactive > 0)
145  {
146  m_Hit->set_eion(0);
149  }
150  else
151  {
152  if (whichactive == -1)
153  {
155  }
156  else
157  {
159  }
160  }
162  {
163  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
164  {
165  m_Hit->set_trkid(pp->GetUserTrackId());
166  m_Hit->set_shower_id(pp->GetShower()->get_id());
167  m_SaveShower = pp->GetShower();
168  }
169  }
170  break;
171  default:
172  break;
173  }
174 
175  if (whichactive > 0)
176  {
177  light_yield = GetVisibleEnergyDeposition(aStep); // for scintillator only, calculate light yields
178  static bool once = true;
179  if (once && edep > 0)
180  {
181  once = false;
182 
183  if (Verbosity() > 0)
184  {
185  std::cout << "PHG4ForwardHcalSteppingAction::UserSteppingAction::"
186  //
187  << m_Detector->GetName() << " - "
188  << " use scintillating light model at each Geant4 steps. "
189  << "First step: "
190  << "Material = "
191  << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
192  << ", "
193  << "Birk Constant = "
195  << ","
196  << "edep = " << edep << ", "
197  << "eion = " << eion
198  << ", "
199  << "light_yield = " << light_yield << std::endl;
200  }
201  }
202  }
203 
204  /* Update exit values- will be overwritten with every step until
205  * we leave the volume or the particle ceases to exist */
206  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
207  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
208  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
209 
210  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
211 
212  /* sum up the energy to get total deposited */
214  if (whichactive > 0)
215  {
216  m_Hit->set_eion(m_Hit->get_eion() + eion);
217  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
218  }
219 
220  if (geantino)
221  {
222  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
223  if (whichactive > 0)
224  {
225  m_Hit->set_eion(-1);
226  m_Hit->set_light_yield(-1);
227  }
228  }
229  if (edep > 0 && (whichactive > 0 ||
230  (whichactive == -1 && m_AbsorberTruthFlag > 0) ||
231  (whichactive < -1 && m_SupportTruthFlag > 0)))
232  {
234  {
235  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
236  {
237  pp->SetKeep(1); // we want to keep the track
238  }
239  }
240  }
241  // if any of these conditions is true this is the last step in
242  // this volume and we need to save the hit
243  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
244  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
245  // (not sure if this will ever be the case)
246  // aTrack->GetTrackStatus() == fStopAndKill: track ends
247  if (postPoint->GetStepStatus() == fGeomBoundary ||
248  postPoint->GetStepStatus() == fWorldBoundary ||
249  postPoint->GetStepStatus() == fAtRestDoItProc ||
250  aTrack->GetTrackStatus() == fStopAndKill)
251  {
252  // save only hits with energy deposit (or -1 for geantino)
253  if (m_Hit->get_edep())
254  {
255  m_SaveHitContainer->AddHit(layer_id, m_Hit);
256  if (m_SaveShower)
257  {
259  }
260  // ownership has been transferred to container, set to null
261  // so we will create a new hit for the next track
262  m_Hit = nullptr;
263  }
264  else
265  {
266  // if this hit has no energy deposit, just reset it for reuse
267  // this means we have to delete it in the dtor. If this was
268  // the last hit we processed the memory is still allocated
269  m_Hit->Reset();
270  }
271  }
272 
273  return true;
274  }
275  else
276  {
277  return false;
278  }
279 }
280 
281 //____________________________________________________________________________..
283 {
284  //now look for the map and grab a pointer to it.
285  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
286  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
287  m_SupportHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_SupportNodeName);
288  // if we do not find the node it's messed up.
289  if (!m_HitContainer)
290  {
291  std::cout << "PHG4ForwardHcalSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
292  }
293  if (!m_AbsorberHitContainer)
294  {
295  if (Verbosity() > 0)
296  {
297  std::cout << "PHG4ForwardHcalSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
298  }
299  }
300  if (!m_SupportHitContainer)
301  {
302  if (Verbosity() > 0)
303  {
304  std::cout << "PHG4ForwardHcalSteppingAction::SetTopNode - unable to find " << m_SupportNodeName << std::endl;
305  }
306  }
307 }