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