ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4HybridHomogeneousCalorimeterSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4HybridHomogeneousCalorimeterSteppingAction.cc
2 
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
13 
15 
16 #include <phool/getClass.h>
17 
18 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
19 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
20 #include <Geant4/G4Step.hh>
21 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
22 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
23 #include <Geant4/G4SystemOfUnits.hh>
24 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
25 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
26 #include <Geant4/G4Track.hh> // for G4Track
27 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
28 #include <Geant4/G4Types.hh> // for G4double
29 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
30 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
31 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
32 
33 #include <TSystem.h>
34 
35 #include <iostream>
36 #include <string> // for basic_string, string
37 
38 class PHCompositeNode;
39 
40 //____________________________________________________________________________..
42  : PHG4SteppingAction(detector->GetName())
43  , m_Detector(detector)
44  , m_ActiveFlag(parameters->get_int_param("active"))
45  , m_BlackHoleFlag(parameters->get_int_param("blackhole"))
46 {
47 }
48 
50 {
51  // if the last hit was a zero energie deposit hit, it is just reset
52  // and the memory is still allocated, so we need to delete it here
53  // if the last hit was saved, hit is a nullptr pointer which are
54  // legal to delete (it results in a no operation)
55  delete m_Hit;
56 }
57 
58 //____________________________________________________________________________..
60 {
63 
64  // m_Detector->IsInCrystalCalorimeter(volume)
65  // returns
66  // 0 is outside of Crystal Calorimeter
67  // 1 is inside scintillator scrystal
68  // -1 is absorber (dead material)
69 
70  int whichactive = m_Detector->IsInCrystalCalorimeter(volume);
71 
72  if (!whichactive)
73  {
74  return false;
75  }
76 
77  int layer_id = m_Detector->get_DetectorId();
78  int idx_j = -1;
79  int idx_k = -1;
80 
81  if (whichactive > 0) // in crystal
82  {
83  // Find indices of crystal containing this step. The indices are
84  // coded into the copy number of the physical volume and we extract them from there
86  {
87  int j[3];
88  int k[3];
89  for (int i = 0; i < 3; i++)
90  {
91  unsigned int icopy = touch->GetVolume(i)->GetCopyNo();
92  j[i] = icopy >> 16;
93  k[i] = icopy & 0xFFFF;
94  }
95  idx_j = j[0] + j[1] * 2 + j[2] * 4;
96  idx_k = k[0] + k[1] * 2 + k[2] * 4;
97  }
99  {
100  unsigned int icopy = touch->GetVolume(1)->GetCopyNo();
101  idx_j = icopy >> 16;
102  idx_k = icopy & 0xFFFF;
103  }
104  }
105 
106  /* Get energy deposited by this step */
107  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
108  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
109  // G4double light_yield = 0;
110 
111  /* Get pointer to associated Geant4 track */
112  const G4Track* aTrack = aStep->GetTrack();
113 
114  // if this block stops everything, just put all kinetic energy into edep
115  if (m_BlackHoleFlag)
116  {
117  edep = aTrack->GetKineticEnergy() / GeV;
118  G4Track* killtrack = const_cast<G4Track*>(aTrack);
119  killtrack->SetTrackStatus(fStopAndKill);
120  }
121 
122  /* Make sure we are in a volume */
123  if (m_ActiveFlag)
124  {
125  /* Check if particle is 'geantino' */
126  bool geantino = false;
127  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
128  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
129  {
130  geantino = true;
131  }
132 
133  /* Get Geant4 pre- and post-step points */
134  G4StepPoint* prePoint = aStep->GetPreStepPoint();
135  G4StepPoint* postPoint = aStep->GetPostStepPoint();
136 
137  switch (prePoint->GetStepStatus())
138  {
139  case fGeomBoundary:
140  case fUndefined:
141  if (!m_Hit)
142  {
143  m_Hit = new PHG4Hitv1();
144  }
145 
146  /* Set hit location (space point)*/
147  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
148  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
149  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
150 
151  /* Set hit time */
152  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
153 
154  /* Set the track ID */
155  m_Hit->set_trkid(aTrack->GetTrackID());
156 
157  /* Set intial energy deposit */
158  m_Hit->set_edep(0);
159 
160  /* Now add the hit to the hit collection */
161  // here we do things which are different between scintillator and absorber hits
162  if (whichactive > 0)
163  {
165  m_Hit->set_light_yield(0); // for scintillator only, initialize light yields
166  m_Hit->set_eion(0);
167  /* Set hit location (tower index) */
168  m_Hit->set_index_j(idx_j);
169  m_Hit->set_index_k(idx_k);
170  }
171  else
172  {
174  }
175 
176  // here we set what is common for scintillator and absorber hits
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  m_SaveShower = pp->GetShower();
184  }
185  }
186  break;
187  default:
188  break;
189  }
190 
191  /* Update exit values- will be overwritten with every step until
192  * we leave the volume or the particle ceases to exist */
193  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
194  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
195  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
196 
197  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
198 
199  /* sum up the energy to get total deposited */
201  if (whichactive > 0)
202  {
203  m_Hit->set_eion(m_Hit->get_eion() + eion);
204 
205  // using the G4 scintillation light production model
206  double light_yield = GetVisibleEnergyDeposition(aStep);
207  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
208  }
209 
210  if (geantino)
211  {
212  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
213  if (whichactive > 0)
214  {
215  m_Hit->set_eion(-1);
216  m_Hit->set_light_yield(-1);
217  }
218  }
219  if (edep > 0)
220  {
222  {
223  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
224  {
225  pp->SetKeep(1); // we want to keep the track
226  }
227  }
228  }
229  // if any of these conditions is true this is the last step in
230  // this volume and we need to save the hit
231  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
232  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
233  // (not sure if this will ever be the case)
234  // aTrack->GetTrackStatus() == fStopAndKill: track ends
235  if (postPoint->GetStepStatus() == fGeomBoundary ||
236  postPoint->GetStepStatus() == fWorldBoundary ||
237  postPoint->GetStepStatus() == fAtRestDoItProc ||
238  aTrack->GetTrackStatus() == fStopAndKill)
239  {
240  // save only hits with energy deposit (or -1 for geantino)
241  if (m_Hit->get_edep())
242  {
243  m_SaveHitContainer->AddHit(layer_id, m_Hit);
244  if (m_SaveShower)
245  {
247  }
248  // ownership has been transferred to container, set to null
249  // so we will create a new hit for the next track
250  m_Hit = nullptr;
251  }
252  else
253  {
254  // if this hit has no energy deposit, just reset it for reuse
255  // this means we have to delete it in the dtor. If this was
256  // the last hit we processed the memory is still allocated
257  m_Hit->Reset();
258  }
259  }
260  return true;
261  }
262  else
263  {
264  return false;
265  }
266 }
267 
268 //____________________________________________________________________________..
270 {
271  std::string hitnodename;
272  std::string absorbernodename;
273 
274  if (m_Detector->SuperDetector() != "NONE")
275  {
276  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
277  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
278  }
279  else
280  {
281  hitnodename = "G4HIT_" + m_Detector->GetName();
282  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
283  }
284 
285  //now look for the map and grab a pointer to it.
286  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
287  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
288 
289  // if we do not find the node it's messed up.
290  if (!m_HitContainer)
291  {
292  std::cout << "PHG4HybridHomogeneousCalorimeterSteppingAction::SetInterfacePointers - unable to find " << hitnodename << std::endl;
293  gSystem->Exit(1);
294  }
295  if (!m_AbsorberHitContainer)
296  {
297  if (Verbosity() > 0)
298  {
299  std::cout << "PHG4HybridHomogeneousCalorimeterSteppingAction::SetInterfacePointers - unable to find " << absorbernodename << std::endl;
300  }
301  }
302 }