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