Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Example02SteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Example02SteppingAction.cc
2 
3 #include "G4Example02Detector.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") != string::npos) // this also accounts for "chargedgeantino"
113  {
114  geantino = true;
115  }
116  G4StepPoint *prePoint = aStep->GetPreStepPoint();
117  G4StepPoint *postPoint = aStep->GetPostStepPoint();
118  // cout << "track id " << aTrack->GetTrackID() << endl;
119  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
120  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
121 
122  switch (prePoint->GetStepStatus())
123  {
124  case fPostStepDoItProc:
125  if (m_SavePostStepStatus != fGeomBoundary)
126  {
127  break;
128  }
129  else
130  {
131  // this is bad from G4 print out diagnostic to help debug, not sure if
132  // this is still with us
133  cout << GetName() << ": New Hit for " << endl;
134  cout << "prestep status: "
135  << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
136  << ", poststep status: "
137  << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
138  << ", last pre step status: "
140  << ", last post step status: "
142  cout << "last track: " << m_SaveTrackId
143  << ", current trackid: " << aTrack->GetTrackID() << endl;
144  cout << "phys pre vol: " << volume->GetName()
145  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
146  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
147  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
148  }
149  case fGeomBoundary:
150  case fUndefined:
151  if (!m_Hit)
152  {
153  m_Hit = new PHG4Hitv1();
154  }
155  m_Hit->set_layer(detector_id);
156  // here we set the entrance values in cm
157  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
158  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
159  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
160  // time in ns
161  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
162  // set the track ID
163  m_Hit->set_trkid(aTrack->GetTrackID());
164  m_SaveTrackId = aTrack->GetTrackID();
165  // set the initial energy deposit
166  m_EdepSum = 0;
167  if (whichactive > 0)
168  {
169  m_EionSum = 0; // assuming the ionization energy is only needed for active
170  // volumes (scintillators)
171  m_Hit->set_eion(0);
173  }
174  else
175  {
176  cout << "implement stuff for whichactive < 0 (inactive volumes)" << endl;
177  gSystem->Exit(1);
178  }
179  // this is for the tracking of the truth info
180  if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
181  {
182  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
183  {
184  m_Hit->set_trkid(pp->GetUserTrackId());
185  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(),
186  m_Hit->get_hit_id());
187  }
188  }
189 
190  break;
191  default:
192  break;
193  }
194 
195  // some sanity checks for inconsistencies (aka bugs)
196  // check if this hit was created, if not print out last post step status
197  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
198  {
199  cout << GetName() << ": hit was not created" << endl;
200  cout << "prestep status: "
201  << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
202  << ", poststep status: "
203  << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
204  << ", last pre step status: "
206  << ", last post step status: "
208  cout << "last track: " << m_SaveTrackId
209  << ", current trackid: " << aTrack->GetTrackID() << endl;
210  cout << "phys pre vol: " << volume->GetName()
211  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
212  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
213  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
214  gSystem->Exit(1);
215  }
216  // check if track id matches the initial one when the hit was created
217  if (aTrack->GetTrackID() != m_SaveTrackId)
218  {
219  cout << GetName() << ": hits do not belong to the same track" << endl;
220  cout << "saved track: " << m_SaveTrackId
221  << ", current trackid: " << aTrack->GetTrackID()
222  << ", prestep status: " << prePoint->GetStepStatus()
223  << ", previous post step status: " << m_SavePostStepStatus << endl;
224 
225  gSystem->Exit(1);
226  }
227  m_SavePreStepStatus = prePoint->GetStepStatus();
228  m_SavePostStepStatus = postPoint->GetStepStatus();
229  m_SaveVolPre = volume;
230  m_SaveVolPost = touchpost->GetVolume();
231 
232  // here we just update the exit values, it will be overwritten
233  // for every step until we leave the volume or the particle
234  // ceases to exist
235  // sum up the energy to get total deposited
236  m_EdepSum += edep;
237  if (whichactive > 0)
238  {
239  m_EionSum += eion;
240  }
241  // if any of these conditions is true this is the last step in
242  // this volume and we need to save the hit
243  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
244  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
245  // (happens when your detector goes outside world volume)
246  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
247  // aTrack->GetTrackStatus() == fStopAndKill is also set)
248  // aTrack->GetTrackStatus() == fStopAndKill: track ends
249  if (postPoint->GetStepStatus() == fGeomBoundary ||
250  postPoint->GetStepStatus() == fWorldBoundary ||
251  postPoint->GetStepStatus() == fAtRestDoItProc ||
252  aTrack->GetTrackStatus() == fStopAndKill)
253  {
254  // save only hits with energy deposit (or geantino)
255  if (m_EdepSum > 0 || geantino)
256  {
257  // update values at exit coordinates and set keep flag
258  // of track to keep
259  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
260  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
261  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
262 
263  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
264  if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
265  {
266  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
267  {
268  pp->SetKeep(1); // we want to keep the track
269  }
270  }
271  if (geantino)
272  {
273  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way
274  // geantinos survive the g4hit compression
275  if (whichactive > 0)
276  {
277  m_Hit->set_eion(-1);
278  }
279  }
280  else
281  {
283  }
284  if (whichactive > 0)
285  {
287  }
288  m_SaveHitContainer->AddHit(detector_id, m_Hit);
289  // ownership has been transferred to container, set to null
290  // so we will create a new hit for the next track
291  m_Hit = nullptr;
292  }
293  else
294  {
295  // if this hit has no energy deposit, just reset it for reuse
296  // this means we have to delete it in the dtor. If this was
297  // the last hit we processed the memory is still allocated
298  m_Hit->Reset();
299  }
300  }
301  // return true to indicate the hit was used
302  return true;
303 }
304 
305 //____________________________________________________________________________..
307 {
308  string hitnodename = "G4HIT_" + m_Detector->GetName();
309 
310  // now look for the map and grab a pointer to it.
311  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
312 
313  // if we do not find the node we need to make it.
314  if (!m_HitContainer)
315  {
316  std::cout << "G4Example02SteppingAction::SetTopNode - unable to find "
317  << hitnodename << std::endl;
318  }
319 }