ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ForwardEcalSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ForwardEcalSteppingAction.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->IsInForwardEcal(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->IsInForwardEcal(volume);
74 
75  if (!whichactive)
76  {
77  return false;
78  }
79 
80  int layer_id = m_Detector->get_Layer();
81  int towertype = m_Detector->get_TowerType();
82  unsigned int icopy = touch->GetVolume(3)->GetCopyNo();
83  if (towertype != 2)
84  icopy = touch->GetVolume(1)->GetCopyNo();
85  int idx_j = icopy >> 16;
86  int idx_k = icopy & 0xFFFF;
87 
88  if (Verbosity() > 2)
89  std::cout << "\t" << icopy << "\t idx_j =" << idx_j << ", idx_k =" << idx_k << "\t id:" << icopy << "\t type: " << towertype << std::endl;
90 
91  /* Get energy deposited by this step */
92  double edep = aStep->GetTotalEnergyDeposit() / GeV;
93  double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
94  double light_yield = 0;
95 
96  /* Get pointer to associated Geant4 track */
97  const G4Track* aTrack = aStep->GetTrack();
98 
99  // if this block stops everything, just put all kinetic energy into edep
100  if (m_BlackHoleFlag)
101  {
102  edep = aTrack->GetKineticEnergy() / GeV;
103  G4Track* killtrack = const_cast<G4Track*>(aTrack);
104  killtrack->SetTrackStatus(fStopAndKill);
105  }
106 
107  /* Make sure we are in a volume */
108  if (m_ActiveFlag)
109  {
110  /* Check if particle is 'geantino' */
111  bool geantino = false;
112  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
113  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
114  {
115  geantino = true;
116  }
117 
118  /* Get Geant4 pre- and post-step points */
119  G4StepPoint* prePoint = aStep->GetPreStepPoint();
120  G4StepPoint* postPoint = aStep->GetPostStepPoint();
121 
122  switch (prePoint->GetStepStatus())
123  {
124  case fGeomBoundary:
125  case fUndefined:
126  if (!m_Hit)
127  {
128  m_Hit = new PHG4Hitv1();
129  }
130 
131  /* Set hit location (space point) */
132  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
133  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
134  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
135 
136  /* Set hit time */
137  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
138 
139  /* Set hit location (tower index) */
140  m_Hit->set_index_j(idx_j);
141  m_Hit->set_index_k(idx_k);
142 
143  //set the track ID
144  m_Hit->set_trkid(aTrack->GetTrackID());
145  /* set intial energy deposit */
146  m_Hit->set_edep(0);
147 
148  /* Now add the hit to the hit collection */
149  // here we do things which are different between scintillator and absorber hits
150  if (whichactive > 0)
151  {
153  m_Hit->set_eion(0);
154  m_Hit->set_light_yield(0); // for scintillator only, initialize light yields
155  }
156  else
157  {
159  }
160  // here we set what is common for scintillator and absorber hits
162  {
163  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
164  {
165  m_Hit->set_trkid(pp->GetUserTrackId());
166  m_Hit->set_shower_id(pp->GetShower()->get_id());
167  m_CurrentShower = pp->GetShower();
168  }
169  }
170  break;
171  default:
172  break;
173  }
174 
175  if (whichactive > 0)
176  {
177  light_yield = GetVisibleEnergyDeposition(aStep); // for scintillator only, calculate light yields
178  static bool once = true;
179  if (once && edep > 0)
180  {
181  once = false;
182 
183  if (Verbosity() > 0)
184  {
185  std::cout << "PHG4ForwardEcalSteppingAction::UserSteppingAction::"
186  //
187  << m_Detector->GetName() << " - "
188  << " use scintillating light model at each Geant4 steps. "
189  << "First step: "
190  << "Material = "
191  << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
192  << ", "
193  << "Birk Constant = "
195  << ","
196  << "edep = " << edep << ", "
197  << "eion = " << eion
198  << ", "
199  << "light_yield = " << light_yield << std::endl;
200  }
201  }
202  }
203 
204  /* Update exit values- will be overwritten with every step until
205  * we leave the volume or the particle ceases to exist */
206  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
207  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
208  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
209 
210  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
211 
212  /* sum up the energy to get total deposited */
214  if (whichactive > 0)
215  {
216  m_Hit->set_eion(m_Hit->get_eion() + eion);
217  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
218  }
219 
220  if (geantino)
221  {
222  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
223  if (whichactive > 0)
224  {
225  m_Hit->set_eion(-1);
226  m_Hit->set_light_yield(-1);
227  }
228  }
229  if (edep > 0 && (whichactive > 0 || m_AbsorberTruthFlag > 0))
230  {
232  {
233  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
234  {
235  pp->SetKeep(1); // we want to keep the track
236  }
237  }
238  }
239  // if any of these conditions is true this is the last step in
240  // this volume and we need to save the hit
241  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
242  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
243  // (not sure if this will ever be the case)
244  // aTrack->GetTrackStatus() == fStopAndKill: track ends
245  if (postPoint->GetStepStatus() == fGeomBoundary ||
246  postPoint->GetStepStatus() == fWorldBoundary ||
247  postPoint->GetStepStatus() == fAtRestDoItProc ||
248  aTrack->GetTrackStatus() == fStopAndKill)
249  {
250  // save only hits with energy deposit (or -1 for geantino)
251  if (m_Hit->get_edep())
252  {
253  m_CurrentHitContainer->AddHit(layer_id, m_Hit);
254  if (m_CurrentShower)
255  {
257  }
258  // ownership has been transferred to container, set to null
259  // so we will create a new hit for the next track
260  m_Hit = nullptr;
261  }
262  else
263  {
264  // if this hit has no energy deposit, just reset it for reuse
265  // this means we have to delete it in the dtor. If this was
266  // the last hit we processed the memory is still allocated
267  m_Hit->Reset();
268  }
269  }
270  return true;
271  }
272  else
273  {
274  return false;
275  }
276 }
277 
278 //____________________________________________________________________________..
280 {
281  std::string hitnodename;
282  std::string absorbernodename;
283 
284  if (m_Detector->SuperDetector() != "NONE")
285  {
286  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
287  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
288  }
289  else
290  {
291  hitnodename = "G4HIT_" + m_Detector->GetName();
292  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
293  }
294 
295  //now look for the map and grab a pointer to it.
296  m_SignalHitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
297  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
298 
299  // if we do not find the node it's messed up.
301  {
302  std::cout << "PHG4ForwardEcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
303  gSystem->Exit(1);
304  }
305  // this is perfectly fine if absorber hits are disabled
306  if (!m_AbsorberHitContainer)
307  {
308  if (Verbosity() > 0)
309  {
310  std::cout << "PHG4ForwardEcalSteppingAction::SetTopNode - unable to find " << absorbernodename << std::endl;
311  }
312  }
313 }