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