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