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