Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4MicromegasSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4MicromegasSteppingAction.cc
1 
7 
9 
10 #include <phparameter/PHParameters.h>
11 
13 
14 #include <g4main/PHG4Hit.h>
16 #include <g4main/PHG4Hitv1.h>
17 #include <g4main/PHG4Shower.h>
20 
21 #include <phool/getClass.h>
22 
23 #include <TSystem.h>
24 
25 #include <Geant4/G4ParticleDefinition.hh>
26 #include <Geant4/G4ReferenceCountedHandle.hh>
27 #include <Geant4/G4Step.hh>
28 #include <Geant4/G4StepPoint.hh>
29 #include <Geant4/G4StepStatus.hh>
30 #include <Geant4/G4String.hh>
31 #include <Geant4/G4SystemOfUnits.hh>
32 #include <Geant4/G4ThreeVector.hh>
33 #include <Geant4/G4TouchableHandle.hh>
34 #include <Geant4/G4Track.hh>
35 #include <Geant4/G4TrackStatus.hh>
36 #include <Geant4/G4Types.hh>
37 #include <Geant4/G4VPhysicalVolume.hh>
38 #include <Geant4/G4VTouchable.hh>
39 #include <Geant4/G4VUserTrackInformation.hh>
40 
41 #include <cmath>
42 #include <iostream>
43 #include <string>
44 
45 class PHCompositeNode;
46 
47 //____________________________________________________________________________..
49  : PHG4SteppingAction(detector->GetName())
50  , m_Detector(detector)
51  , m_Params(parameters)
52  , m_ActiveFlag(m_Params->get_int_param("active"))
53  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
54 {
55 }
56 
57 //____________________________________________________________________________..
58 // This is the implementation of the G4 UserSteppingAction
59 bool PHG4MicromegasSteppingAction::UserSteppingAction(const G4Step *aStep, bool /*was_used*/)
60 {
61  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
62  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
63 
64  // get volume of the current step
65  G4VPhysicalVolume *volume = touch->GetVolume();
66  int whichactive = m_Detector->IsInDetector(volume);
67  if (whichactive == 0)
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  // if this detector stops everything, just put all kinetic energy into edep
78  if (m_BlackHoleFlag)
79  {
80  edep = aTrack->GetKineticEnergy() / GeV;
81  auto killtrack = const_cast<G4Track *>(aTrack);
82  killtrack->SetTrackStatus(fStopAndKill);
83  }
84 
85  bool geantino =
86  aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
87  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos;
88 
89  G4StepPoint *prePoint = aStep->GetPreStepPoint();
90  G4StepPoint *postPoint = aStep->GetPostStepPoint();
91 
92  // Here we have to decide if we need to create a new hit. Normally this should
93  // only be neccessary if a G4 Track enters a new volume or is freshly created
94  // For this we look at the step status of the prePoint (beginning of the G4 Step).
95  // This should be either fGeomBoundary (G4 Track crosses into volume) or
96  // fUndefined (G4 Track newly created)
97  // Sadly over the years with different G4 versions we have observed cases where
98  // G4 produces "impossible hits" which we try to catch here
99  // These errors were always rare and it is not clear if they still exist but we
100  // still check for them for safety. We can reproduce G4 runs identically (if given
101  // the sequence of random number seeds you find in the log), the printouts help
102  // us giving the G4 support information about those failures
103  switch (prePoint->GetStepStatus())
104  {
105  case fPostStepDoItProc:
106  if (m_SavePostStepStatus == fGeomBoundary)
107  {
108  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
109  // this is still with us
110  std::cout << "PHG4MicromegasSteppingAction::UserSteppingAction - " << GetName() << ": New Hit for " << std::endl;
111  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
112  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
113  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
114  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus)
115  << std::endl;
116 
117  std::cout << "last track: " << m_SaveTrackId << ", current trackid: " << aTrack->GetTrackID() << std::endl;
118  std::cout << "phys pre vol: " << volume->GetName() << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
119  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName() << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
120  }
121  break;
122 
123  // These are the normal cases
124  case fGeomBoundary:
125  case fUndefined:
126  {
127  if (!m_hit)
128  {
129  m_hit.reset(new PHG4Hitv1());
130  }
131 
132  if (whichactive > 0)
133  {
134  // assign layer
135  m_hit->set_layer(m_Detector->get_layer(volume));
136 
137  const auto tileid = m_Detector->get_tileid(volume);
138  m_hit->set_property(PHG4Hit::prop_index_i, tileid);
139  // momentum
140  m_hit->set_px(0, prePoint->GetMomentum().x() / GeV);
141  m_hit->set_py(0, prePoint->GetMomentum().y() / GeV);
142  m_hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
143 
144  m_hit->set_eion(0);
146  }
147  else
148  {
150  }
151  // here we set the entrance values in cm
152  m_hit->set_x(0, prePoint->GetPosition().x() / cm);
153  m_hit->set_y(0, prePoint->GetPosition().y() / cm);
154  m_hit->set_z(0, prePoint->GetPosition().z() / cm);
155 
156  // time in ns
157  m_hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
158 
159  // set the track ID
160  m_hit->set_trkid(aTrack->GetTrackID());
161  m_SaveTrackId = aTrack->GetTrackID();
162 
163  // reset the initial energy deposit
164  m_EdepSum = 0;
165  m_EionSum = 0;
166  m_hit->set_edep(0);
167 
168  // this is for the tracking of the truth info
169  if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
170  {
171  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
172  {
173  m_hit->set_trkid(pp->GetUserTrackId());
174  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_hit->get_hit_id());
175  }
176  }
177  break;
178  }
179 
180  default:
181  break;
182  }
183 
184  if (!m_hit || !std::isfinite(m_hit->get_x(0)))
185  {
186  std::cout << GetName() << ": hit was not created" << std::endl;
187  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
188  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
189  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
190  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus)
191  << std::endl;
192  std::cout << "last track: " << m_SaveTrackId << ", current trackid: " << aTrack->GetTrackID() << std::endl;
193  std::cout << "phys pre vol: " << volume->GetName() << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
194  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName() << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
195 
196  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
197  gSystem->Exit(1);
198  }
199 
200  // check if track id matches the initial one when the hit was created
201  if (aTrack->GetTrackID() != m_SaveTrackId)
202  {
203  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
204  std::cout << "saved track: " << m_SaveTrackId
205  << ", current trackid: " << aTrack->GetTrackID()
206  << ", prestep status: " << prePoint->GetStepStatus()
207  << ", previous post step status: " << m_SavePostStepStatus << std::endl;
208  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
209  gSystem->Exit(1);
210  }
211 
212  // We need to cache a few things from one step to the next
213  // to identify impossible hits and subsequent debugging printout
214  m_SavePreStepStatus = prePoint->GetStepStatus();
215  m_SavePostStepStatus = postPoint->GetStepStatus();
216  m_SaveVolPre = volume;
217  m_SaveVolPost = touchpost->GetVolume();
218 
219  // here we just update the exit values, it will be overwritten
220  // for every step until we leave the volume or the particle
221  // ceases to exist
222  // sum up the energy to get total deposited
223  m_EdepSum += edep;
224  m_EionSum += eion;
225 
226  // if any of these conditions is true this is the last step in
227  // this volume and we need to save the hit
228  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
229  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
230  // (happens when your detector goes outside world volume)
231  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
232  // aTrack->GetTrackStatus() == fStopAndKill is also set)
233  // aTrack->GetTrackStatus() == fStopAndKill: track ends
234  if (postPoint->GetStepStatus() == fGeomBoundary ||
235  postPoint->GetStepStatus() == fWorldBoundary ||
236  postPoint->GetStepStatus() == fAtRestDoItProc ||
237  aTrack->GetTrackStatus() == fStopAndKill)
238  {
239  // save only hits with energy deposit (or geantino)
240  if (m_EdepSum > 0 || geantino)
241  {
242  // update values at exit coordinates and set keep flag
243  // of track to keep
244  m_hit->set_x(1, postPoint->GetPosition().x() / cm);
245  m_hit->set_y(1, postPoint->GetPosition().y() / cm);
246  m_hit->set_z(1, postPoint->GetPosition().z() / cm);
247 
248  m_hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
249  if (whichactive > 0)
250  {
251  m_hit->set_px(1, postPoint->GetMomentum().x() / GeV);
252  m_hit->set_py(1, postPoint->GetMomentum().y() / GeV);
253  m_hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
254  }
255 
256  if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
257  {
258  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
259  {
260  pp->SetKeep(1); // we want to keep the track
261  }
262  }
263 
264  if (geantino)
265  {
266  m_hit->set_edep(-1);
267  }
268  else
269  {
270  m_hit->set_edep(m_EdepSum);
271  }
272  if (whichactive > 0)
273  {
274  if (geantino)
275  {
276  m_hit->set_eion(-1);
277  }
278  else
279  {
280  m_hit->set_eion(m_EionSum);
281  }
282  }
283 
284  // add in container
285  m_SaveHitContainer->AddHit(m_hit->get_layer(), m_hit.get());
286 
287  // ownership is transferred to container
288  // so we release the hit
289  m_hit.release();
290  }
291  else
292  {
293  m_hit->Reset();
294  }
295  }
296  // return true to indicate the hit was used
297  return true;
298 }
299 
300 //____________________________________________________________________________..
302 {
303  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
304  m_SupportHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_SupportNodeName);
305 
306  if (!m_HitContainer)
307  {
308  if (!m_Params->get_int_param("blackhole")) // not messed up if we have a black hole
309  {
310  std::cout << "PHG4MicromegasSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
311  gSystem->Exit(1);
312  }
313  }
314  // this is perfectly fine if support hits are disabled
315  if (!m_SupportHitContainer)
316  {
317  if (Verbosity() > 0)
318  {
319  std::cout << "PHG4MicromegasSteppingAction::SetTopNode - unable to find " << m_SupportNodeName << std::endl;
320  }
321  }
322 }
323 
325 {
326  if (type == "G4HIT")
327  {
329  return;
330  }
331  else if (type == "G4HIT_SUPPORT")
332  {
334  return;
335  }
336  std::cout << "Invalid output hit node type " << type << std::endl;
337  gSystem->Exit(1);
338  return;
339 }