ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ECAPToFSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ECAPToFSteppingAction.cc
2 #include "PHG4ECAPToFDetector.h"
3 #include "PHG4ECAPToFSubsystem.h"
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
13 
14 #include <phool/getClass.h>
15 
16 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
17 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
18 #include <Geant4/G4Step.hh>
19 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
20 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fPostSt...
21 #include <Geant4/G4SystemOfUnits.hh>
22 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
23 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
24 #include <Geant4/G4Track.hh> // for G4Track
25 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
26 #include <Geant4/G4Types.hh> // for G4double
27 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
28 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
29 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
30 
31 #include <boost/io/ios_state.hpp>
32 
33 #include <cmath> // for isfinite, copysign
34 #include <cstdlib> // for exit
35 #include <iomanip>
36 #include <iostream>
37 #include <string> // for operator<<, char_traits
38 
39 class PHCompositeNode;
40 
42  : PHG4SteppingAction(detector->GetName())
43  , m_Detector(detector)
44  , m_Params(parameters)
45  , m_HitContainer(nullptr)
46  //, m_AbsorberHits(nullptr)
47  , m_ActiveGasHits(nullptr)
48  , m_Hit(nullptr)
49  , m_SaveHitContainer(nullptr)
50  , m_SaveShower(nullptr)
51  , m_SaveVolPre(nullptr)
52  , m_SaveVolPost(nullptr)
53  , m_SaveTrackId(-1)
54  , m_SavePreStepStatus(-1)
55  , m_SavePostStepStatus(-1)
56  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
57  , m_ActiveFlag(m_Params->get_int_param("active"))
58  , m_UseG4StepsFlag(m_Params->get_int_param("use_g4steps"))
59  , m_Zmin(m_Params->get_double_param("z_begin") * cm)
60  , m_Zmax(m_Params->get_double_param("z_begin") * cm + 3.0)
61 {
62  m_Zmin -= copysign(m_Zmin, 1. / 1e6 * cm);
63  m_Zmax += copysign(m_Zmax, 1. / 1e6 * cm);
64 }
65 
67 {
68  delete m_Hit;
69 }
70 
72 {
73  // get volume of the current step
76 
78  // G4 just calls UserSteppingAction for every step (and we loop then over all our
79  // steppingactions. First we have to check if we are actually in our volume
80  if (!m_Detector->IsInToF(volume))
81  {
82  return false;
83  }
84 
85  int whichactive = m_Detector->IsInToF(volume);
86 
87  // collect energy and track length step by step
89 
90  const G4Track* aTrack = aStep->GetTrack();
91  // if this volume stops everything, just put all kinetic energy into edep
92  if (m_BlackHoleFlag)
93  {
94  edep = aTrack->GetKineticEnergy() / GeV;
95  G4Track* killtrack = const_cast<G4Track*>(aTrack);
96  killtrack->SetTrackStatus(fStopAndKill);
97  }
98  int layer_id = m_Detector->get_Layer();
99 
100  if (m_ActiveFlag)
101  {
102  bool geantino = false;
103  // the check for the pdg code speeds things up, I do not want to make
104  // an expensive string compare for every track when we know
105  // geantino or chargedgeantino has pid=0
106  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
107  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
108  {
109  geantino = true;
110  }
111  G4StepPoint* prePoint = aStep->GetPreStepPoint();
112  G4StepPoint* postPoint = aStep->GetPostStepPoint();
113  // std::cout << "time prepoint: " << prePoint->GetGlobalTime()/ns << std::endl;
114  // std::cout << "time postpoint: " << postPoint->GetGlobalTime()/ns << std::endl;
115  // std::cout << "kinetic energy: " << aTrack->GetKineticEnergy()/GeV << std::endl;
116  // G4ParticleDefinition* def = aTrack->GetDefinition();
117  // std::cout << "Particle: " << def->GetParticleName() << std::endl;
118  int prepointstatus = prePoint->GetStepStatus();
119  if (prepointstatus == fGeomBoundary ||
120  prepointstatus == fUndefined ||
121  (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary) ||
122  m_UseG4StepsFlag > 0)
123  {
124  if (!m_Hit)
125  {
126  m_Hit = new PHG4Hitv1();
127  }
128 
129  m_Hit->set_layer((unsigned int) layer_id);
130 
131  //here we set the entrance values in cm
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  m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
137  m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
138  m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
139 
140  // time in ns
141  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
142  //set and save the track ID
143  m_Hit->set_trkid(aTrack->GetTrackID());
144  m_SaveTrackId = aTrack->GetTrackID();
145  //set the initial energy deposit
146  m_Hit->set_edep(0);
147 
148  if (whichactive < 0)
149  {
151  }
152  else
153  {
154  //m_SaveHitContainer = m_AbsorberHits;
156  }
157 
158  if (!geantino && !m_BlackHoleFlag)
159  {
160  m_Hit->set_eion(0);
161  }
163  {
164  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
165  {
166  m_Hit->set_trkid(pp->GetUserTrackId());
167  m_Hit->set_shower_id(pp->GetShower()->get_id());
168  m_SaveShower = pp->GetShower();
169  }
170  }
171 
172  } // END ..... || m_UseG4StepsFlag > 0)
173 
174  m_SavePostStepStatus = postPoint->GetStepStatus();
175  // check if track id matches the initial one when the hit was created
176  if (aTrack->GetTrackID() != m_SaveTrackId)
177  {
178  std::cout << "hits do not belong to the same track" << std::endl;
179  std::cout << "saved track: " << m_SaveTrackId
180  << ", current trackid: " << aTrack->GetTrackID()
181  << std::endl;
182  exit(1);
183  }
184  m_SavePreStepStatus = prePoint->GetStepStatus();
185  m_SavePostStepStatus = postPoint->GetStepStatus();
187  m_SaveVolPost = touchpost->GetVolume();
188 
189  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
190  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
191  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
192 
193  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
194  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
195  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
196 
197  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
198  //sum up the energy to get total deposited
200 
201  if (geantino)
202  {
203  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
204  }
205  else
206  {
207  if (!m_BlackHoleFlag)
208  {
209  double eion = edep - aStep->GetNonIonizingEnergyDeposit() / GeV;
210  m_Hit->set_eion(m_Hit->get_eion() + eion);
211  }
212  }
213 
214  if (edep > 0 || m_SaveAllHitsFlag)
215  {
217  {
218  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
219  {
220  pp->SetKeep(1); // we want to keep the track
221  }
222  }
223  }
224 
225  // if any of these conditions is true this is the last step in
226  // this volume and we need to save the hit
227  if (postPoint->GetStepStatus() == fGeomBoundary ||
228  postPoint->GetStepStatus() == fWorldBoundary ||
229  postPoint->GetStepStatus() == fAtRestDoItProc ||
230  aTrack->GetTrackStatus() == fStopAndKill ||
231  m_UseG4StepsFlag > 0)
232  {
233  // save only hits with energy deposit (or -1 for geantino) or if save all hits flag is set
234  if (m_Hit->get_edep())
235  {
236  m_SaveHitContainer->AddHit(layer_id, m_Hit);
237 
238  if (m_SaveShower)
239  {
241  }
242  // ownership has been transferred to container, set to null
243  // so we will create a new hit for the next track
244  m_Hit = nullptr;
245  }
246  else
247  {
248  // if this hit has no energy deposit, just reset it for reuse
249  // this means we have to delete it in the dtor. If this was
250  // the last hit we processed the memory is still allocated
251  m_Hit->Reset();
252  }
253  }
254  // return true to indicate the hit was used
255  return true;
256  } // END Acitve flag condition
257 
258  else
259  {
260  return false;
261  }
262 }
263 
265 {
266  std::string hitnodename;
267  //string absorbernodename;
268  std::string gasnodename;
269 
270  if (m_Detector->SuperDetector() != "NONE")
271  {
272  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
273  gasnodename = "G4HIT_ACTIVEGAS_" + m_Detector->SuperDetector();
274  }
275  else
276  {
277  hitnodename = "G4HIT_" + m_Detector->GetName();
278  gasnodename = "G4HIT_ACTIVEGAS_" + m_Detector->SuperDetector();
279  }
280 
281  //now look for the map and grab a pointer to it.
282  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
283  m_ActiveGasHits = findNode::getClass<PHG4HitContainer>(topNode, gasnodename.c_str());
284 
285  // if we do not find the node it's messed up.
286  if (!m_HitContainer)
287  {
288  std::cout << "PHG4ECAPToFSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
289  }
290  // if (!m_AbsorberHits)
291  if (!m_ActiveGasHits)
292  {
293  if (Verbosity() > 1)
294  {
295  std::cout << "PHG4ECAPToFSteppingAction::SetTopNode - unable to find " << gasnodename << std::endl;
296  }
297  }
298 }