ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4InnerHcalSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4InnerHcalSteppingAction.cc
2 
4 #include "PHG4StepStatusDecode.h"
5 
6 #include <phparameter/PHParameters.h>
7 
8 #include <g4main/PHG4Hit.h>
10 #include <g4main/PHG4Hitv1.h>
11 #include <g4main/PHG4Shower.h>
12 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
14 
15 #include <phool/getClass.h>
16 
17 #include <TSystem.h>
18 
19 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
20 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
21 #include <Geant4/G4Step.hh>
22 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
23 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
24 #include <Geant4/G4String.hh> // for G4String
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 <cmath> // for isfinite
36 #include <iostream>
37 #include <string> // for operator<<, operator+
38 #include <utility> // for pair
39 
40 class PHCompositeNode;
41 
42 using namespace std;
43 //____________________________________________________________________________..
45  : PHG4SteppingAction(detector->GetName())
46  , m_Detector(detector)
47  , m_Hits(nullptr)
48  , m_AbsorberHits(nullptr)
49  , m_Hit(nullptr)
50  , m_Params(parameters)
51  , m_SaveHitContainer(nullptr)
52  , m_SaveShower(nullptr)
53  , m_SaveVolPre(nullptr)
54  , m_SaveVolPost(nullptr)
55  , m_SaveTrackId(-1)
56  , m_SavePreStepStatus(-1)
57  , m_SavePostStepStatus(-1)
58  , m_IsActive(m_Params->get_int_param("active"))
59  , m_IsBlackHole(m_Params->get_int_param("blackhole"))
60  , m_LightScintModel(m_Params->get_int_param("light_scint_model"))
61 {
62  SetLightCorrection(m_Params->get_double_param("light_balance_inner_radius") * cm,
63  m_Params->get_double_param("light_balance_inner_corr"),
64  m_Params->get_double_param("light_balance_outer_radius") * cm,
65  m_Params->get_double_param("light_balance_outer_corr"));
66 }
67 
69 {
70  // if the last hit was a zero energie deposit hit, it is just reset
71  // and the memory is still allocated, so we need to delete it here
72  // if the last hit was saved, hit is a nullptr pointer which are
73  // legal to delete (it results in a no operation)
74  delete m_Hit;
75 }
76 //____________________________________________________________________________..
78 {
81  // get volume of the current step
83 
84  // m_Detector->IsInInnerHcal(volume)
85  // returns
86  // 0 is outside of InnerHcal
87  // 1 is inside scintillator
88  // -1 is steel absorber
89 
90  int whichactive = m_Detector->IsInInnerHcal(volume);
91 
92  if (!whichactive)
93  {
94  return false;
95  }
96  int layer_id = -1;
97  int tower_id = -1;
98  if (whichactive > 0) // scintillator
99  {
100  pair<int, int> layer_tower = m_Detector->GetLayerTowerId(volume);
101  layer_id = layer_tower.first;
102  tower_id = layer_tower.second;
103  // cout << "name " << volume->GetName() << ", mid: " << layer_id
104  // << ", twr: " << tower_id << endl;
105  }
106  else
107  {
108  layer_id = touch->GetCopyNumber(); // steel plate id
109  }
110  // collect energy and track length step by step
111  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
112  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
113  G4double light_yield = 0;
114  const G4Track* aTrack = aStep->GetTrack();
115 
116  // if this block stops everything, just put all kinetic energy into edep
117  if (m_IsBlackHole)
118  {
119  edep = aTrack->GetKineticEnergy() / GeV;
120  G4Track* killtrack = const_cast<G4Track*>(aTrack);
121  killtrack->SetTrackStatus(fStopAndKill);
122  }
123 
124  // make sure we are in a volume
125  if (m_IsActive)
126  {
127  bool geantino = false;
128 
129  // the check for the pdg code speeds things up, I do not want to make
130  // an expensive string compare for every track when we know
131  // geantino or chargedgeantino has pid=0
132  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
133  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
134  {
135  geantino = true;
136  }
137  G4StepPoint* prePoint = aStep->GetPreStepPoint();
138  G4StepPoint* postPoint = aStep->GetPostStepPoint();
139  // cout << "track id " << aTrack->GetTrackID() << endl;
140  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
141  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
142  switch (prePoint->GetStepStatus())
143  {
144  case fPostStepDoItProc:
146  {
147  break;
148  }
149  else
150  {
151  cout << GetName() << ": New Hit for " << endl;
152  cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
153  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
154  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
155  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
156  cout << "last track: " << m_SaveTrackId
157  << ", current trackid: " << aTrack->GetTrackID() << endl;
158  cout << "phys pre vol: " << volume->GetName()
159  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
160  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
161  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
162  }
163  [[fallthrough]];
164  case fGeomBoundary:
165  case fUndefined:
166  // if previous hit was saved, hit pointer was set to nullptr
167  // and we have to make a new one
168  if (!m_Hit)
169  {
170  m_Hit = new PHG4Hitv1();
171  }
172  //here we set the entrance values in cm
173  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
174  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
175  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
176  // time in ns
177  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
178  //set and save the track ID
179  m_Hit->set_trkid(aTrack->GetTrackID());
180  m_SaveTrackId = aTrack->GetTrackID();
181  //set the initial energy deposit
182  m_Hit->set_edep(0);
183  if (whichactive > 0) // return of IsInInnerHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
184  {
185  m_Hit->set_scint_id(tower_id); // the slat id
186  m_Hit->set_eion(0); // only implemented for v5 otherwise empty
187  m_Hit->set_light_yield(0); // for scintillator only, initialize light yields
188  // Now save the container we want to add this hit to
190  }
191  else
192  {
194  }
196  {
197  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
198  {
199  m_Hit->set_trkid(pp->GetUserTrackId());
200  m_Hit->set_shower_id(pp->GetShower()->get_id());
201  m_SaveShower = pp->GetShower();
202  }
203  }
204  break;
205  default:
206  break;
207  }
208  // some sanity checks for inconsistencies
209  // check if this hit was created, if not print out last post step status
210  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
211  {
212  cout << GetName() << ": hit was not created" << endl;
213  cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
214  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
215  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
216  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
217  cout << "last track: " << m_SaveTrackId
218  << ", current trackid: " << aTrack->GetTrackID() << endl;
219  cout << "phys pre vol: " << volume->GetName()
220  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
221  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
222  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
223  gSystem->Exit(1);
224  }
225  // check if track id matches the initial one when the hit was created
226  if (aTrack->GetTrackID() != m_SaveTrackId)
227  {
228  cout << GetName() << ": hits do not belong to the same track" << endl;
229  cout << "saved track: " << m_SaveTrackId
230  << ", current trackid: " << aTrack->GetTrackID()
231  << ", prestep status: " << prePoint->GetStepStatus()
232  << ", previous post step status: " << m_SavePostStepStatus
233  << endl;
234 
235  gSystem->Exit(1);
236  }
237  m_SavePreStepStatus = prePoint->GetStepStatus();
238  m_SavePostStepStatus = postPoint->GetStepStatus();
240  m_SaveVolPost = touchpost->GetVolume();
241  // here we just update the exit values, it will be overwritten
242  // for every step until we leave the volume or the particle
243  // ceases to exist
244  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
245  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
246  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
247 
248  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
249 
250  //sum up the energy to get total deposited
252  if (whichactive > 0) // return of IsInInnerHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
253  {
254  m_Hit->set_eion(m_Hit->get_eion() + eion);
255  light_yield = eion;
256  if (m_LightScintModel)
257  {
258  light_yield = GetVisibleEnergyDeposition(aStep); // for scintillator only, calculate light yields
259  }
260  light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
261  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
262  }
263  if (geantino)
264  {
265  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
266  m_Hit->set_eion(-1);
267  }
268  if (edep > 0)
269  {
271  {
272  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
273  {
274  pp->SetKeep(1); // we want to keep the track
275  }
276  }
277  }
278 
279  // if any of these conditions is true this is the last step in
280  // this volume and we need to save the hit
281  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
282  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
283  // (happens when your detector goes outside world volume)
284  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
285  // aTrack->GetTrackStatus() == fStopAndKill is also set)
286  // aTrack->GetTrackStatus() == fStopAndKill: track ends
287  if (postPoint->GetStepStatus() == fGeomBoundary ||
288  postPoint->GetStepStatus() == fWorldBoundary ||
289  postPoint->GetStepStatus() == fAtRestDoItProc ||
290  aTrack->GetTrackStatus() == fStopAndKill)
291  {
292  // save only hits with energy deposit (or -1 for geantino)
293  if (m_Hit->get_edep())
294  {
295  m_SaveHitContainer->AddHit(layer_id, m_Hit);
296  if (m_SaveShower)
297  {
299  }
300  // ownership has been transferred to container, set to null
301  // so we will create a new hit for the next track
302  m_Hit = nullptr;
303  }
304  else
305  {
306  // if this hit has no energy deposit, just reset it for reuse
307  // this means we have to delete it in the dtor. If this was
308  // the last hit we processed the memory is still allocated
309  m_Hit->Reset();
310  }
311  }
312  // return true to indicate the hit was used
313  return true;
314  }
315  else
316  {
317  return false;
318  }
319 }
320 
321 //____________________________________________________________________________..
323 {
324  string hitnodename;
325  string absorbernodename;
326  if (m_Detector->SuperDetector() != "NONE")
327  {
328  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
329  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
330  }
331  else
332  {
333  hitnodename = "G4HIT_" + m_Detector->GetName();
334  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
335  }
336 
337  //now look for the map and grab a pointer to it.
338  m_Hits = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
339  m_AbsorberHits = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
340 
341  // if we do not find the node it's messed up.
342  if (!m_Hits)
343  {
344  std::cout << "PHG4InnerHcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
345  }
346  if (!m_AbsorberHits)
347  {
348  if (Verbosity() > 1)
349  {
350  cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
351  }
352  }
353 }