ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4HcalSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4HcalSteppingAction.cc
2 #include "PHG4HcalDetector.h"
3 
4 #include <g4main/PHG4Hit.h>
6 #include <g4main/PHG4Hitv1.h>
7 #include <g4main/PHG4Shower.h>
8 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
9 
11 
12 #include <phool/getClass.h>
13 
14 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
15 #include <Geant4/G4Step.hh>
16 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
17 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRestD...
18 #include <Geant4/G4String.hh> // for G4String
19 #include <Geant4/G4SystemOfUnits.hh> // for cm, GeV, nanosecond
20 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
21 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
22 #include <Geant4/G4Track.hh> // for G4Track
23 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
24 #include <Geant4/G4Types.hh> // for G4double
25 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
26 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
27 
28 #include <iostream>
29 #include <string> // for char_traits, operator<<
30 
31 class G4VPhysicalVolume;
32 class PHCompositeNode;
33 
34 using namespace std;
35 //____________________________________________________________________________..
37  : PHG4SteppingAction(detector->GetName())
38  , detector_(detector)
39 {
40 }
41 
42 //____________________________________________________________________________..
44 {
45  // if the last hit was a zero energie deposit hit, it is just reset
46  // and the memory is still allocated, so we need to delete it here
47  // if the last hit was saved, hit is a nullptr pointer which are
48  // legal to delete (it results in a no operation)
49  delete m_Hit;
50 }
51 
52 //____________________________________________________________________________..
54 {
55  // get volume of the current step
57 
58  // collect energy and track length step by step
60  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
61  G4double light_yield = 0;
62 
63  const G4Track* aTrack = aStep->GetTrack();
64 
65  int layer_id = detector_->get_Layer();
66  // make sure we are in a volume
67  // IsInCylinderActive returns the number of the scintillator
68  // slat which has fired
69  int isactive = detector_->IsInCylinderActive(volume);
70  if (isactive > PHG4HcalDetector::INACTIVE)
71  {
72  bool geantino = false;
73  // the check for the pdg code speeds things up, I do not want to make
74  // an expensive string compare for every track when we know
75  // geantino or chargedgeantino has pid=0
76  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
77  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
78  {
79  geantino = true;
80  }
81  G4StepPoint* prePoint = aStep->GetPreStepPoint();
82  G4StepPoint* postPoint = aStep->GetPostStepPoint();
83  // cout << "track id " << aTrack->GetTrackID() << endl;
84  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
85  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
86  switch (prePoint->GetStepStatus())
87  {
88  case fGeomBoundary:
89  case fUndefined:
90 
91  if (!m_Hit)
92  {
93  m_Hit = new PHG4Hitv1();
94  }
95  m_Hit->set_layer((unsigned int) layer_id);
96  m_Hit->set_scint_id(isactive); // isactive contains the scintillator slat id
97  //here we set the entrance values in cm
98  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
99  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
100  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
101 
102  // time in ns
103  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
104  //set the track ID
105  m_Hit->set_trkid(aTrack->GetTrackID());
107  {
108  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
109  {
110  m_Hit->set_trkid(pp->GetUserTrackId());
111  m_Hit->set_shower_id(pp->GetShower()->get_id());
112  m_SaveShower = pp->GetShower();
113  }
114  }
115 
116  //set the initial energy deposit
117  m_Hit->set_edep(0);
118  m_Hit->set_eion(0); // only implemented for v5 otherwise empty
120  // m_Hit->print();
121  // Now add the hit
122  if (isactive >= 0) // the slat ids start with zero
123  {
125  // unsigned int shift_layer_id = layer_id << (phg4hitdefs::keybits - 3);
126  }
127  else
128  {
130  }
131  // m_SaveHitContainer->AddHit(layer_id, m_Hit);
132 
134  {
135  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
136  {
137  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
138  }
139  }
140  if (m_Hit->get_z(0) > zmax || m_Hit->get_z(0) < zmin)
141  {
142  cout << "PHG4HcalSteppingAction: hit outside acceptance, layer: " << layer_id << endl;
143  m_Hit->identify();
144  }
145  break;
146  default:
147  break;
148  }
149  // here we just update the exit values, it will be overwritten
150  // for every step until we leave the volume or the particle
151  // ceases to exist
152  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
153  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
154  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
155 
156  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
157 
158  if (isactive >= 0) // the slat ids start with zero
159  {
160  if (light_scint_model_)
161  {
162  light_yield = GetVisibleEnergyDeposition(aStep);
163  }
164  else
165  {
166  light_yield = eion;
167  }
168 
169  if (ValidCorrection())
170  {
171  light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x() / cm, (postPoint->GetPosition().y() / cm));
172  }
173  }
174 
175  //sum up the energy to get total deposited
177  m_Hit->set_eion(m_Hit->get_eion() + eion);
178  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
179  m_Hit->set_path_length(aTrack->GetTrackLength() / cm);
180 
181  if (m_Hit->get_z(1) > zmax || m_Hit->get_z(1) < zmin)
182  {
183  cout << "PHG4HcalSteppingAction: hit outside acceptance zmin " << zmin << ", zmax " << zmax << " at exit" << endl;
184  m_Hit->identify();
185  }
186  if (geantino)
187  {
188  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
189  m_Hit->set_eion(-1);
190  }
191  if (edep > 0)
192  {
194  {
195  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
196  {
197  pp->SetKeep(1); // we want to keep the track
198  }
199  }
200  }
201  // if any of these conditions is true this is the last step in
202  // this volume and we need to save the hit
203  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
204  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
205  // (happens when your detector goes outside world volume)
206  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
207  // aTrack->GetTrackStatus() == fStopAndKill is also set)
208  // aTrack->GetTrackStatus() == fStopAndKill: track ends
209  if (postPoint->GetStepStatus() == fGeomBoundary ||
210  postPoint->GetStepStatus() == fWorldBoundary ||
211  postPoint->GetStepStatus() == fAtRestDoItProc ||
212  aTrack->GetTrackStatus() == fStopAndKill)
213  {
214  // save only hits with energy deposit (or -1 for geantino)
215  if (m_Hit->get_edep())
216  {
217  m_SaveHitContainer->AddHit(layer_id, m_Hit);
218  if (m_SaveShower)
219  {
221  }
222  // ownership has been transferred to container, set to null
223  // so we will create a new hit for the next track
224  m_Hit = nullptr;
225  }
226  else
227  {
228  // if this hit has no energy deposit, just reset it for reuse
229  // this means we have to delete it in the dtor. If this was
230  // the last hit we processed the memory is still allocated
231  m_Hit->Reset();
232  }
233  }
234  // return true to indicate the hit was used
235  // hit->identify();
236  // return true to indicate the hit was used
237  return true;
238  }
239  else
240  {
241  return false;
242  }
243 }
244 
245 //____________________________________________________________________________..
247 {
248  string hitnodename;
249  string absorbernodename;
250  if (detector_->SuperDetector() != "NONE")
251  {
252  hitnodename = "G4HIT_" + detector_->SuperDetector();
253  absorbernodename = "G4HIT_ABSORBER_" + detector_->SuperDetector();
254  }
255  else
256  {
257  hitnodename = "G4HIT_" + detector_->GetName();
258  absorbernodename = "G4HIT_ABSORBER_" + detector_->GetName();
259  }
260 
261  //now look for the map and grab a pointer to it.
262  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
263  m_AbsorberHits = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
264  // if we do not find the node we need to make it.
265  if (!m_HitContainer)
266  {
267  std::cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
268  }
269  if (!m_AbsorberHits)
270  {
271  if (Verbosity() > 0)
272  {
273  std::cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << std::endl;
274  }
275  }
276 }