ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4PSTOFSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4PSTOFSteppingAction.cc
2 #include "PHG4PSTOFDetector.h"
3 #include "PHG4StepStatusDecode.h"
4 
5 #include <g4main/PHG4Hit.h>
7 #include <g4main/PHG4Hitv1.h>
8 #include <g4main/PHG4Shower.h>
9 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
11 
12 #include <phool/getClass.h>
13 #include <phool/phool.h> // for PHWHERE
14 
15 #include <TSystem.h>
16 
17 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
18 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
19 #include <Geant4/G4Step.hh>
20 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
21 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
22 #include <Geant4/G4String.hh> // for G4String
23 #include <Geant4/G4SystemOfUnits.hh>
24 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
25 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
26 #include <Geant4/G4Track.hh> // for G4Track
27 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
28 #include <Geant4/G4Types.hh> // for G4double
29 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
30 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
31 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
32 
33 #include <cmath> // for isfinite
34 #include <iostream>
35 #include <string> // for operator<<, string
36 
37 class PHCompositeNode;
38 
39 using namespace std;
40 //____________________________________________________________________________..
42  : PHG4SteppingAction(detector->GetName())
43  , detector_(detector)
44  , hits_(nullptr)
45  , hit(nullptr)
46  , savehitcontainer(nullptr)
47  , savevolpre(nullptr)
48  , savevolpost(nullptr)
49  , savetrackid(-1)
50  , saveprestepstatus(-1)
51  , savepoststepstatus(-1)
52  , edepsum(0)
53  , eionsum(0)
54 {}
55 
57 {
58  // if the last hit was a zero energie deposit hit, it is just reset
59  // and the memory is still allocated, so we need to delete it here
60  // if the last hit was saved, hit is a nullptr pointer which are
61  // legal to delete (it results in a no operation)
62  delete hit;
63 }
64 
65 //____________________________________________________________________________..
66 bool PHG4PSTOFSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
67 {
70  // get volume of the current step
72 // IsInPSTOF(volume) returns
73 // == 0 outside of pstof
74 // > 0 for hits in active volume
75 // < 0 for hits in passive material
76  int whichactive = detector_->IsInPSTOF(volume);
77  if (!whichactive)
78  {
79  return false;
80  }
81 
82  // collect energy and track length step by step
84  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
85  const G4Track* aTrack = aStep->GetTrack();
86 
87  int layer_id = 0; // what the heck is this?
88  bool geantino = false;
89  // the check for the pdg code speeds things up, I do not want to make
90  // an expensive string compare for every track when we know
91  // geantino or chargedgeantino has pid=0
92  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
93  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
94  {
95  geantino = true;
96  }
97  G4StepPoint* prePoint = aStep->GetPreStepPoint();
98  G4StepPoint* postPoint = aStep->GetPostStepPoint();
99  // cout << "track id " << aTrack->GetTrackID() << endl;
100  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
101  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
102 
103  layer_id = touch->GetCopyNumber();
104  if (layer_id != whichactive)
105  {
106  cout << PHWHERE << " inconsistency between G4 copy number: "
107  << layer_id << " and module id from detector: "
108  << whichactive << endl;
109  gSystem->Exit(1);
110  }
111 
112  switch (prePoint->GetStepStatus())
113  {
114  case fPostStepDoItProc:
116  {
117  break;
118  }
119  else
120  {
121  cout << GetName() << ": New Hit for " << endl;
122 cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
123 << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
124  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(saveprestepstatus)
125  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(savepoststepstatus) << endl;
126  cout << "last track: " << savetrackid
127  << ", current trackid: " << aTrack->GetTrackID() << endl;
128  cout << "phys pre vol: " << volume->GetName()
129  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
130  cout << " previous phys pre vol: " << savevolpre->GetName()
131  << " previous phys post vol: " << savevolpost->GetName() << endl;
132  }
133  [[fallthrough]];
134  case fGeomBoundary:
135  case fUndefined:
136  if (!hit)
137  {
138  hit = new PHG4Hitv1();
139  }
140  hit->set_layer(layer_id);
141  //here we set the entrance values in cm
142  hit->set_x(0, prePoint->GetPosition().x() / cm);
143  hit->set_y(0, prePoint->GetPosition().y() / cm);
144  hit->set_z(0, prePoint->GetPosition().z() / cm);
145  // time in ns
146  hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
147  //set the track ID
148  hit->set_trkid(aTrack->GetTrackID());
149  savetrackid = aTrack->GetTrackID();
151  {
152  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
153  {
154  hit->set_trkid(pp->GetUserTrackId());
155  }
156  }
157  //set the initial energy deposit
158  edepsum = 0;
159  if (whichactive > 0)
160  {
161  eionsum = 0;
162  hit->set_eion(0);
164  }
165  else
166  {
167  cout << "implement stuff for whichactive < 0" << endl;
168  gSystem->Exit(1);
169  }
171  {
172  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
173  {
174  hit->set_trkid(pp->GetUserTrackId());
175  pp->GetShower()->add_g4hit_id(savehitcontainer->GetID(), hit->get_hit_id());
176  }
177  }
178 
179  break;
180  default:
181  break;
182  }
183 
184  // some sanity checks for inconsistencies
185  // check if this hit was created, if not print out last post step status
186  if (!hit || !isfinite(hit->get_x(0)))
187  {
188  cout << GetName() << ": hit was not created" << endl;
189  cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
190  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
191  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(saveprestepstatus)
192  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(savepoststepstatus) << endl;
193  cout << "last track: " << savetrackid
194  << ", current trackid: " << aTrack->GetTrackID() << endl;
195  cout << "phys pre vol: " << volume->GetName()
196  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
197  cout << " previous phys pre vol: " << savevolpre->GetName()
198  << " previous phys post vol: " << savevolpost->GetName() << endl;
199  gSystem->Exit(1);
200  }
201  // check if track id matches the initial one when the hit was created
202  if (aTrack->GetTrackID() != savetrackid)
203  {
204  cout << GetName() << ": hits do not belong to the same track" << endl;
205  cout << "saved track: " << savetrackid
206  << ", current trackid: " << aTrack->GetTrackID()
207  << ", prestep status: " << prePoint->GetStepStatus()
208  << ", previous post step status: " << savepoststepstatus
209  << endl;
210 
211  gSystem->Exit(1);
212  }
213  saveprestepstatus = prePoint->GetStepStatus();
214  savepoststepstatus = postPoint->GetStepStatus();
215  savevolpre = volume;
216  savevolpost = touchpost->GetVolume();
217 
218  // here we just update the exit values, it will be overwritten
219  // for every step until we leave the volume or the particle
220  // ceases to exist
221  //sum up the energy to get total deposited
222  edepsum += edep;
223  if (whichactive > 0)
224  {
225  eionsum += eion;
226  }
227  // if any of these conditions is true this is the last step in
228  // this volume and we need to save the hit
229  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
230  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
231  // (happens when your detector goes outside world volume)
232  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
233  // aTrack->GetTrackStatus() == fStopAndKill is also set)
234  // aTrack->GetTrackStatus() == fStopAndKill: track ends
235  if (postPoint->GetStepStatus() == fGeomBoundary ||
236  postPoint->GetStepStatus() == fWorldBoundary ||
237  postPoint->GetStepStatus() == fAtRestDoItProc ||
238  aTrack->GetTrackStatus() == fStopAndKill)
239  {
240  // save only hits with energy deposit (or geantino)
241  if (edepsum > 0 || geantino)
242  {
243 // update values at exit coordinates and set keep flag
244 // of track to keep
245  hit->set_x(1, postPoint->GetPosition().x() / cm);
246  hit->set_y(1, postPoint->GetPosition().y() / cm);
247  hit->set_z(1, postPoint->GetPosition().z() / cm);
248 
249  hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
251  {
252  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
253  {
254  pp->SetKeep(1); // we want to keep the track
255  }
256  }
257  if (geantino)
258  {
259  hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
260  if (whichactive > 0)
261  {
262  hit->set_eion(-1);
263  }
264  }
265  else
266  {
267  hit->set_edep(edepsum);
268  }
269  if (whichactive > 0)
270  {
271  hit->set_eion(eionsum);
272  }
273  savehitcontainer->AddHit(layer_id, hit);
274  // ownership has been transferred to container, set to null
275  // so we will create a new hit for the next track
276  hit = nullptr;
277  }
278  else
279  {
280  // if this hit has no energy deposit, just reset it for reuse
281  // this means we have to delete it in the dtor. If this was
282  // the last hit we processed the memory is still allocated
283  hit->Reset();
284  }
285  }
286  // return true to indicate the hit was used
287  return true;
288 }
289 
290 //____________________________________________________________________________..
292 {
293  string hitnodename;
294  if (detector_->SuperDetector() != "NONE")
295  {
296  hitnodename = "G4HIT_" + detector_->SuperDetector();
297  }
298  else
299  {
300  hitnodename = "G4HIT_" + detector_->GetName();
301  }
302 
303  //now look for the map and grab a pointer to it.
304  hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
305 
306  // if we do not find the node we need to make it.
307  if (!hits_)
308  {
309  std::cout << "PHG4PSTOFSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
310  }
311 }