Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4Prototype2InnerHcalSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4Prototype2InnerHcalSteppingAction.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
13 
14 #include <phool/getClass.h>
15 
16 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
17 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
18 #include <Geant4/G4Step.hh>
19 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
20 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
21 #include <Geant4/G4String.hh> // for G4String
22 #include <Geant4/G4SystemOfUnits.hh>
23 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
24 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
25 #include <Geant4/G4Track.hh> // for G4Track
26 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
27 #include <Geant4/G4Types.hh> // for G4double
28 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
29 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
30 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
31 
32 #include <cmath> // for isfinite, sqrt
33 #include <iostream>
34 #include <string> // for string, operator+, ope...
35 
36 class PHCompositeNode;
37 
38 //#define TESTSINGLESLAT
39 #ifdef TESTSINGLESLAT
40 static const int nrow = 4;
41 static const int nslat = 9;
42 #endif
43 
44 using namespace std;
45 //____________________________________________________________________________..
47  : PHG4SteppingAction(detector->GetName())
48  , m_Detector(detector)
49  , m_HitContainer(nullptr)
50  , m_AbsorberHitContainer(nullptr)
51  , m_Hit(nullptr)
52  , m_Params(parameters)
53  , m_SaveHitContainer(nullptr)
54  , m_SaveShower(nullptr)
55  , m_AbsorberTruthFlag(m_Params->get_int_param("absorbertruth"))
56  , m_IsActiveFlag(m_Params->get_int_param("active"))
57  , m_IsBlackHoleFlag(m_Params->get_int_param("blackhole"))
58  , m_LightScintModelFlag(m_Params->get_int_param("light_scint_model"))
59  , m_LightBalanceInnerCorr(m_Params->get_double_param("light_balance_inner_corr"))
60  , m_LightBalanceInnerRadius(m_Params->get_double_param("light_balance_inner_radius") * cm)
61  , m_LightBalanceOuterCorr(m_Params->get_double_param("light_balance_outer_corr"))
62  , m_LightBalanceOuterRadius(m_Params->get_double_param("light_balance_outer_radius") * cm)
63 {
64 }
65 
67 {
68  // if the last hit was a zero energie deposit hit, it is just reset
69  // and the memory is still allocated, so we need to delete it here
70  // if the last hit was saved, hit is a nullptr pointer which are
71  // legal to delete (it results in a no operation)
72  delete m_Hit;
73 }
74 
75 //____________________________________________________________________________..
77 {
78  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
79  // get volume of the current step
80  G4VPhysicalVolume* volume = touch->GetVolume();
81 
82  // m_Detector->IsInPrototype2InnerHcal(volume)
83  // returns
84  // 0 is outside of Prototype2InnerHcal
85  // 1 is inside scintillator
86  // -1 is steel absorber
87 
88  int whichactive = m_Detector->IsInPrototype2InnerHcal(volume);
89 
90  if (!whichactive)
91  {
92  return false;
93  }
94  int row_id = -1;
95  int slat_id = -1;
96  if (whichactive > 0) // scintillator
97  {
98  slat_id = volume->GetCopyNo();
99  // the row id comes from saved info in the detector construction
100  G4VPhysicalVolume* mothervolume = touch->GetVolume(1);
101  row_id = m_Detector->get_scinti_row_id(mothervolume->GetName());
102  // cout << "mother volume: " << mothervolume->GetName()
103  // << ", volume name " << volume->GetName() << ", row: " << row_id
104  // << ", column: " << slat_id << endl;
105 #ifdef TESTSINGLESLAT
106  if (row_id != nrow)
107  {
108  return false;
109  }
110  if (slat_id != nslat)
111  {
112  return false;
113  }
114 #endif
115  }
116  else
117  {
118  // the row id comes from saved info in the detector construction
119  row_id = m_Detector->get_steel_plate_id(volume->GetName());
120 #ifdef TESTSINGLESLAT
121  return false;
122 #endif
123  }
124  // collect energy and track length step by step
125  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
126  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
127  G4double light_yield = 0;
128  const G4Track* aTrack = aStep->GetTrack();
129 
130  // if this block stops everything, just put all kinetic energy into edep
131  if (m_IsBlackHoleFlag)
132  {
133  edep = aTrack->GetKineticEnergy() / GeV;
134  G4Track* killtrack = const_cast<G4Track*>(aTrack);
135  killtrack->SetTrackStatus(fStopAndKill);
136  }
137  int layer_id = m_Detector->get_Layer();
138  // make sure we are in a volume
139  if (m_IsActiveFlag)
140  {
141  bool geantino = false;
142 
143  // the check for the pdg code speeds things up, I do not want to make
144  // an expensive string compare for every track when we know
145  // geantino or chargedgeantino has pid=0
146  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
147  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
148  {
149  geantino = true;
150  }
151  G4StepPoint* prePoint = aStep->GetPreStepPoint();
152  G4StepPoint* postPoint = aStep->GetPostStepPoint();
153  // cout << "track id " << aTrack->GetTrackID() << endl;
154  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
155  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
156  switch (prePoint->GetStepStatus())
157  {
158  case fGeomBoundary:
159  case fUndefined:
160  if (!m_Hit)
161  {
162  m_Hit = new PHG4Hitv1();
163  }
164  m_Hit->set_row(row_id);
165  if (whichactive > 0) // only for scintillators
166  {
167  m_Hit->set_scint_id(slat_id); // the slat id in the mother volume (or steel plate id), the column
168  }
169  //here we set the entrance values in cm
170  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
171  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
172  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
173  // time in ns
174  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
175  //set the track ID
176  m_Hit->set_trkid(aTrack->GetTrackID());
177  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
178  {
179  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
180  {
181  m_Hit->set_trkid(pp->GetUserTrackId());
182  m_Hit->set_shower_id(pp->GetShower()->get_id());
183  }
184  }
185 
186  //set the initial energy deposit
187  m_Hit->set_edep(0);
188 
189  m_Hit->set_hit_type(0);
190  if ((aTrack->GetParticleDefinition()->GetParticleName().find("e+") != string::npos) ||
191  (aTrack->GetParticleDefinition()->GetParticleName().find("e-") != string::npos))
192  m_Hit->set_hit_type(1);
193 
194  if (whichactive > 0) // return of IsInPrototype2InnerHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
195  {
197  m_Hit->set_light_yield(0); // for scintillator only, initialize light yields
198  m_Hit->set_eion(0);
199  }
200  else
201  {
203  }
204 
205  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
206  {
207  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
208  {
209  m_Hit->set_trkid(pp->GetUserTrackId());
210  m_Hit->set_shower_id(pp->GetShower()->get_id());
211  m_SaveShower = pp->GetShower();
212  }
213  }
214  break;
215  default:
216  break;
217  }
218  // here we just update the exit values, it will be overwritten
219  // for every step until we leave the volume or the particle
220  // ceases to exist
221  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
222  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
223  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
224 
225  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
226 
227  if (whichactive > 0) // return of IsInPrototype2InnerHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
228  {
229  m_Hit->set_eion(m_Hit->get_eion() + eion);
230  light_yield = eion;
232  {
233  light_yield = GetVisibleEnergyDeposition(aStep); // for scintillator only, calculate light yields
234  }
235  if (isfinite(m_LightBalanceOuterRadius) &&
236  isfinite(m_LightBalanceInnerRadius) &&
237  isfinite(m_LightBalanceOuterCorr) &&
238  isfinite(m_LightBalanceInnerCorr))
239  {
240  double r = sqrt(postPoint->GetPosition().x() * postPoint->GetPosition().x() + postPoint->GetPosition().y() * postPoint->GetPosition().y());
241  double cor = GetLightCorrection(r);
242  light_yield = light_yield * cor;
243  }
244  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
245  }
246 
247  //sum up the energy to get total deposited
248  m_Hit->set_edep(m_Hit->get_edep() + edep);
249  if (geantino)
250  {
251  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
252  if (whichactive > 0) // add light yield for scintillators
253  {
254  m_Hit->set_light_yield(-1);
255  m_Hit->set_eion(-1);
256  }
257  }
258  if (edep > 0 && (whichactive > 0 || m_AbsorberTruthFlag > 0))
259  {
260  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
261  {
262  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
263  {
264  pp->SetKeep(1); // we want to keep the track
265  }
266  }
267  }
268  // if any of these conditions is true this is the last step in
269  // this volume and we need to save the hit
270  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
271  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
272  // (happens when your detector goes outside world volume)
273  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
274  // aTrack->GetTrackStatus() == fStopAndKill is also set)
275  // aTrack->GetTrackStatus() == fStopAndKill: track ends
276  if (postPoint->GetStepStatus() == fGeomBoundary ||
277  postPoint->GetStepStatus() == fWorldBoundary ||
278  postPoint->GetStepStatus() == fAtRestDoItProc ||
279  aTrack->GetTrackStatus() == fStopAndKill)
280  {
281  // save only hits with energy deposit (or -1 for geantino)
282  if (m_Hit->get_edep())
283  {
284  m_SaveHitContainer->AddHit(layer_id, m_Hit);
285  if (m_SaveShower)
286  {
288  }
289  // ownership has been transferred to container, set to null
290  // so we will create a new hit for the next track
291  m_Hit = nullptr;
292  }
293  else
294  {
295  // if this hit has no energy deposit, just reset it for reuse
296  // this means we have to delete it in the dtor. If this was
297  // the last hit we processed the memory is still allocated
298  m_Hit->Reset();
299  }
300  }
301 
302  // m_Hit->identify();
303  // return true to indicate the hit was used
304  return true;
305  }
306  else
307  {
308  return false;
309  }
310 }
311 
312 //____________________________________________________________________________..
314 {
315  string hitnodename;
316  string absorbernodename;
317  if (m_Detector->SuperDetector() != "NONE")
318  {
319  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
320  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
321  }
322  else
323  {
324  hitnodename = "G4HIT_" + m_Detector->GetName();
325  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
326  }
327 
328  //now look for the map and grab a pointer to it.
329  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
330  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
331 
332  // if we do not find the node it's messed up.
333  if (!m_HitContainer)
334  {
335  std::cout << "PHG4Prototype2InnerHcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
336  }
337  if (!m_AbsorberHitContainer)
338  {
339  if (Verbosity() > 1)
340  {
341  cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
342  }
343  }
344 }
345 
346 double
348 {
351  double value = m * r + b;
352  if (value > 1.0) return 1.0;
353  if (value < 0.0) return 0.0;
354 
355  return value;
356 }