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