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