ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4CrystalCalorimeterSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4CrystalCalorimeterSteppingAction.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);
205  }
206 
207  if (geantino)
208  {
209  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
210  if (whichactive > 0)
211  {
212  m_Hit->set_eion(-1);
213  m_Hit->set_light_yield(-1);
214  }
215  }
216  if (edep > 0)
217  {
219  {
220  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
221  {
222  pp->SetKeep(1); // we want to keep the track
223  }
224  }
225  }
226  // if any of these conditions is true this is the last step in
227  // this volume and we need to save the hit
228  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
229  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
230  // (not sure if this will ever be the case)
231  // aTrack->GetTrackStatus() == fStopAndKill: track ends
232  if (postPoint->GetStepStatus() == fGeomBoundary ||
233  postPoint->GetStepStatus() == fWorldBoundary ||
234  postPoint->GetStepStatus() == fAtRestDoItProc ||
235  aTrack->GetTrackStatus() == fStopAndKill)
236  {
237  // save only hits with energy deposit (or -1 for geantino)
238  if (m_Hit->get_edep())
239  {
240  m_SaveHitContainer->AddHit(layer_id, m_Hit);
241  if (m_SaveShower)
242  {
244  }
245  // ownership has been transferred to container, set to null
246  // so we will create a new hit for the next track
247  m_Hit = nullptr;
248  }
249  else
250  {
251  // if this hit has no energy deposit, just reset it for reuse
252  // this means we have to delete it in the dtor. If this was
253  // the last hit we processed the memory is still allocated
254  m_Hit->Reset();
255  }
256  }
257  return true;
258  }
259  else
260  {
261  return false;
262  }
263 }
264 
265 //____________________________________________________________________________..
267 {
268  std::string hitnodename;
269  std::string absorbernodename;
270 
271  if (m_Detector->SuperDetector() != "NONE")
272  {
273  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
274  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
275  }
276  else
277  {
278  hitnodename = "G4HIT_" + m_Detector->GetName();
279  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
280  }
281 
282  //now look for the map and grab a pointer to it.
283  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
284  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
285 
286  // if we do not find the node it's messed up.
287  if (!m_HitContainer)
288  {
289  std::cout << "PHG4CrystalCalorimeterSteppingAction::SetInterfacePointers - unable to find " << hitnodename << std::endl;
290  gSystem->Exit(1);
291  }
292  if (!m_AbsorberHitContainer)
293  {
294  if (Verbosity() > 0)
295  {
296  std::cout << "PHG4CrystalCalorimeterSteppingAction::SetInterfacePointers - unable to find " << absorbernodename << std::endl;
297  }
298  }
299 }