Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4BbcSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4BbcSteppingAction.cc
2 
3 #include "PHG4BbcDetector.h"
4 #include "PHG4StepStatusDecode.h"
5 
6 #include <phparameter/PHParameters.h>
7 
8 #include <g4main/PHG4Hit.h>
10 #include <g4main/PHG4Hitv1.h>
11 #include <g4main/PHG4Shower.h>
12 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
14 
15 #include <phool/getClass.h>
16 
17 #include <TSystem.h>
18 
19 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
20 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
21 #include <Geant4/G4Step.hh>
22 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
23 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
24 #include <Geant4/G4String.hh> // for G4String
25 #include <Geant4/G4SystemOfUnits.hh>
26 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
27 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
28 #include <Geant4/G4Track.hh> // for G4Track
29 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
30 #include <Geant4/G4Types.hh> // for G4double
31 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
32 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
33 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
34 
35 #include <cmath> // for isfinite
36 #include <iostream>
37 #include <string> // for operator<<, string
38 
39 class PHCompositeNode;
40 
41 //____________________________________________________________________________..
43  : PHG4SteppingAction(detector->GetName())
44  , m_Detector(detector)
45  , m_Params(parameters)
46  , m_ActiveFlag(m_Params->get_int_param("active"))
47  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
48  , m_SupportFlag(m_Params->get_int_param("supportactive"))
49 {
50 }
51 
53 {
54  // if the last hit was a zero energie deposit hit, it is just reset
55  // and the memory is still allocated, so we need to delete it here
56  // if the last hit was saved, hit is a nullptr pointer which are
57  // legal to delete (it results in a no operation)
58  delete m_Hit;
59 }
60 
61 //____________________________________________________________________________..
62 bool PHG4BbcSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
63 {
64  //std::cout << PHWHERE << " In PHG4BbcSteppingAction::UserSteppingAction()" << std::endl;
65 
66  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
67  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
68  // get volume of the current step
69  G4VPhysicalVolume* volume = touch->GetVolume();
70  // IsInBbc(volume) returns
71  // == 0 outside of bbc
72  // > 0 for hits in active volume
73  // < 0 for hits in passive material
74  int whichactive = m_Detector->IsInBbc(volume);
75  if (!whichactive)
76  {
77  return false;
78  }
79 
80  //std::cout << PHWHERE << " Found Hit" << std::endl;
81 
82  // collect energy and track length step by step
83  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
84  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
85  G4double steplen = aStep->GetStepLength() / cm;
86  const G4Track* aTrack = aStep->GetTrack();
87 
88  // check steplength vs tracklength
89  /*
90  G4double tracklen = aTrack->GetTrackLength() / cm;
91  if ( tracklen != steplen )
92  {
93  std::cout << "YYY " << tracklen << "\t" << steplen << std::endl;
94  }
95  */
96 
97  // if this detector stops everything, just put all kinetic energy into edep
98  if (m_BlackHoleFlag)
99  {
100  edep = aTrack->GetKineticEnergy() / GeV;
101  G4Track* killtrack = const_cast<G4Track*>(aTrack);
102  killtrack->SetTrackStatus(fStopAndKill);
103  if (!m_ActiveFlag)
104  {
105  return false;
106  }
107  }
108  if (whichactive < 0 && !m_SupportFlag)
109  {
110  return false;
111  }
112  bool geantino = false;
113  // the check for the pdg code speeds things up, I do not want to make
114  // an expensive string compare for every track when we know
115  // geantino or chargedgeantino has pid=0
116  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
117  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
118  {
119  geantino = true;
120  }
121 
122  G4StepPoint* prePoint = aStep->GetPreStepPoint();
123  G4StepPoint* postPoint = aStep->GetPostStepPoint();
124  // std::cout << "track id " << aTrack->GetTrackID() << std::endl;
125  // std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
126  // std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
127 
128  //int detector_id = touch->GetCopyNumber(); // not used
129  int tube_id = touch->GetCopyNumber(1);
130 
131  // Create a new hit if a G4 Track enters a new volume or is freshly created
132  // For this we look at the step status of the prePoint (beginning of the G4 Step).
133  // This should be either fGeomBoundary (G4 Track crosses into volume) or
134  // fUndefined (G4 Track newly created)
135 
136  switch (prePoint->GetStepStatus())
137  {
138  case fPostStepDoItProc:
139  if (m_SavePostStepStatus != fGeomBoundary)
140  {
141  // this is the okay case, fPostStepDoItProc called in a volume, not first thing inside
142  // a new volume, just proceed here
143  break;
144  }
145  else
146  {
147  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
148  // this is still with us
149  std::cout << GetName() << ": New Hit for " << std::endl;
150  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
151  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
152  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
153  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
154  std::cout << "last track: " << m_SaveTrackId
155  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
156  std::cout << "phys pre vol: " << volume->GetName()
157  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
158  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
159  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
160  }
161  [[fallthrough]];
162  // These are the normal cases
163  case fGeomBoundary:
164  case fUndefined:
165  if (!m_Hit)
166  {
167  m_Hit = new PHG4Hitv1();
168  }
169  m_Hit->set_layer(tube_id);
170  m_Hit->set_scint_id(tube_id);
171 
172  //here we set the entrance values in cm
173  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
174  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
175  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
176  // time in ns
177  m_Hit->set_t(0, prePoint->GetGlobalTime() / ns);
178  //set the track ID
179  m_Hit->set_trkid(aTrack->GetTrackID());
180  m_SaveTrackId = aTrack->GetTrackID();
181 
182  //set the initial energy deposit
183  m_EdepSum = 0;
184  if (whichactive > 0)
185  {
186  m_EionSum = 0;
187  m_Hit->set_eion(0);
188 
189  m_PathLen = 0.;
191 
193  }
194  else
195  {
197  }
198 
199  // this is for the tracking of the truth info
200  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
201  {
202  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
203  {
204  m_Hit->set_trkid(pp->GetUserTrackId());
205  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
206  }
207  }
208 
209  break;
210  default:
211  break;
212  }
213 
214  // some sanity checks for inconsistencies
215  // check if this hit was created, if not print out last post step status
216  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
217  {
218  std::cout << GetName() << ": hit was not created" << std::endl;
219  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
220  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
221  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
222  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
223  std::cout << "last track: " << m_SaveTrackId
224  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
225  std::cout << "phys pre vol: " << volume->GetName()
226  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
227  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
228  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
229  gSystem->Exit(1);
230  }
231 
232  // check if track id matches the initial one when the hit was created
233  if (aTrack->GetTrackID() != m_SaveTrackId)
234  {
235  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
236  std::cout << "saved track: " << m_SaveTrackId
237  << ", current trackid: " << aTrack->GetTrackID()
238  << ", prestep status: " << prePoint->GetStepStatus()
239  << ", previous post step status: " << m_SavePostStepStatus
240  << std::endl;
241 
242  gSystem->Exit(1);
243  }
244 
245  // We need to cache a few things from one step to the next
246  // to identify impossible hits and subsequent debugging printout
247  m_SavePreStepStatus = prePoint->GetStepStatus();
248  m_SavePostStepStatus = postPoint->GetStepStatus();
249  m_SaveVolPre = volume;
250  m_SaveVolPost = touchpost->GetVolume();
251 
252  // here we just update the exit values, it will be overwritten
253  // for every step until we leave the volume or the particle
254  // ceases to exist
255 
256  // Sum up the energies and lengths to get totals
257  m_EdepSum += edep;
258  if (whichactive > 0)
259  {
260  m_EionSum += eion;
261  m_PathLen += steplen;
262  }
263 
264  // if any of these conditions is true this is the last step in
265  // this volume and we need to save the hit
266  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
267  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
268  // (happens when your detector goes outside world volume)
269  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
270  // aTrack->GetTrackStatus() == fStopAndKill is also set)
271  // aTrack->GetTrackStatus() == fStopAndKill: track ends
272  if (postPoint->GetStepStatus() == fGeomBoundary ||
273  postPoint->GetStepStatus() == fWorldBoundary ||
274  postPoint->GetStepStatus() == fAtRestDoItProc ||
275  aTrack->GetTrackStatus() == fStopAndKill)
276  {
277  // save only hits with energy deposit (or geantino)
278  if (m_EdepSum > 0 || geantino)
279  {
280  // update values at exit coordinates and set keep flag
281  // of track to keep
282  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
283  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
284  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
285 
286  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
287  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
288  {
289  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
290  {
291  pp->SetKeep(1); // we want to keep the track
292  }
293  }
294  if (geantino)
295  {
296  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
297  if (whichactive > 0)
298  {
299  m_Hit->set_eion(-1);
300  m_Hit->set_path_length(-1);
301  }
302  }
303  else
304  {
306  }
307  if (whichactive > 0)
308  {
311  }
312  m_SaveHitContainer->AddHit(tube_id, m_Hit);
313  // ownership has been transferred to container, set to null
314  // so we will create a new hit for the next track
315  m_Hit = nullptr;
316  }
317  else
318  {
319  // if this hit has no energy deposit, just reset it for reuse
320  // this means we have to delete it in the dtor. If this was
321  // the last hit we processed the memory is still allocated
322  m_Hit->Reset();
323  }
324  }
325  // return true to indicate the hit was used
326  return true;
327 }
328 
329 //____________________________________________________________________________..
331 {
332  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
333  m_SupportHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_SupportNodeName);
334  // if we do not find the node it's messed up.
335  if (!m_HitContainer)
336  {
337  if (!m_Params->get_int_param("blackhole")) // not messed up if we have a black hole
338  {
339  std::cout << "PHG4BbcSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
340  gSystem->Exit(1);
341  }
342  }
343  // this is perfectly fine if support hits are disabled
344  if (!m_SupportHitContainer)
345  {
346  if (Verbosity() > 0)
347  {
348  std::cout << "PHG4BbcSteppingAction::SetTopNode - unable to find " << m_SupportNodeName << std::endl;
349  }
350  }
351 }
352 
354 {
355  if (type == "G4HIT")
356  {
358  return;
359  }
360  else if (type == "G4HIT_SUPPORT")
361  {
363  return;
364  }
365  std::cout << "Invalid output hit node type " << type << std::endl;
366  gSystem->Exit(1);
367  return;
368 }