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