ECCE @ EIC 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 // This is the implementation of the G4 UserSteppingAction
58 bool PHG4MicromegasSteppingAction::UserSteppingAction(const G4Step *aStep,bool /*was_used*/)
59 {
62 
63  // get volume of the current step
65  if (!m_Detector->IsInDetector(volume)) return false;
66 
67  // collect energy and track length step by step
69  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
70  const G4Track *aTrack = aStep->GetTrack();
71 
72  // if this detector stops everything, just put all kinetic energy into edep
73  if (m_BlackHoleFlag)
74  {
75  edep = aTrack->GetKineticEnergy() / GeV;
76  auto killtrack = const_cast<G4Track *>(aTrack);
77  killtrack->SetTrackStatus(fStopAndKill);
78  }
79 
80  bool geantino =
81  aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
82  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos;
83 
84  G4StepPoint *prePoint = aStep->GetPreStepPoint();
85  G4StepPoint *postPoint = aStep->GetPostStepPoint();
86 
87  // Here we have to decide if we need to create a new hit. Normally this should
88  // only be neccessary if a G4 Track enters a new volume or is freshly created
89  // For this we look at the step status of the prePoint (beginning of the G4 Step).
90  // This should be either fGeomBoundary (G4 Track crosses into volume) or
91  // fUndefined (G4 Track newly created)
92  // Sadly over the years with different G4 versions we have observed cases where
93  // G4 produces "impossible hits" which we try to catch here
94  // These errors were always rare and it is not clear if they still exist but we
95  // still check for them for safety. We can reproduce G4 runs identically (if given
96  // the sequence of random number seeds you find in the log), the printouts help
97  // us giving the G4 support information about those failures
98  switch (prePoint->GetStepStatus())
99  {
100 
101  case fPostStepDoItProc:
103  {
104 
105  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
106  // this is still with us
107  std::cout << "PHG4MicromegasSteppingAction::UserSteppingAction - " << GetName() << ": New Hit for " << std::endl;
108  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
109  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
110  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
111  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus)
112  << std::endl;
113 
114  std::cout << "last track: " << m_SaveTrackId << ", current trackid: " << aTrack->GetTrackID() << std::endl;
115  std::cout << "phys pre vol: " << volume->GetName() << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
116  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName() << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
117  }
118  break;
119 
120  // These are the normal cases
121  case fGeomBoundary:
122  case fUndefined:
123  {
124  if (!m_hit) m_hit.reset( new PHG4Hitv1() );
125 
126  // assign layer
127  m_hit->set_layer(m_Detector->get_layer(volume));
128 
129  // here we set the entrance values in cm
130  m_hit->set_x(0, prePoint->GetPosition().x() / cm);
131  m_hit->set_y(0, prePoint->GetPosition().y() / cm);
132  m_hit->set_z(0, prePoint->GetPosition().z() / cm);
133 
134  // momentum
135  m_hit->set_px(0, prePoint->GetMomentum().x() / GeV);
136  m_hit->set_py(0, prePoint->GetMomentum().y() / GeV);
137  m_hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
138 
139  // time in ns
140  m_hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
141 
142  // set the track ID
143  m_hit->set_trkid(aTrack->GetTrackID());
144  m_SaveTrackId = aTrack->GetTrackID();
145 
146  // reset the initial energy deposit
147  m_EdepSum = 0;
148  m_EionSum = 0;
149  m_hit->set_edep(0);
150  m_hit->set_eion(0);
152 
153  // this is for the tracking of the truth info
155  {
156  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
157  {
158  m_hit->set_trkid(pp->GetUserTrackId());
159  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_hit->get_hit_id());
160  }
161  }
162  break;
163  }
164 
165  default: break;
166 
167  }
168 
169  if (!m_hit || !std::isfinite(m_hit->get_x(0)))
170  {
171  std::cout << GetName() << ": hit was not created" << std::endl;
172  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
173  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
174  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
175  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus)
176  << std::endl;
177  std::cout << "last track: " << m_SaveTrackId << ", current trackid: " << aTrack->GetTrackID() << std::endl;
178  std::cout << "phys pre vol: " << volume->GetName() << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
179  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName() << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
180 
181  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
182  gSystem->Exit(1);
183  }
184 
185  // check if track id matches the initial one when the hit was created
186  if (aTrack->GetTrackID() != m_SaveTrackId)
187  {
188  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
189  std::cout << "saved track: " << m_SaveTrackId
190  << ", current trackid: " << aTrack->GetTrackID()
191  << ", prestep status: " << prePoint->GetStepStatus()
192  << ", previous post step status: " << m_SavePostStepStatus << std::endl;
193  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
194  gSystem->Exit(1);
195  }
196 
197  // We need to cache a few things from one step to the next
198  // to identify impossible hits and subsequent debugging printout
199  m_SavePreStepStatus = prePoint->GetStepStatus();
200  m_SavePostStepStatus = postPoint->GetStepStatus();
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  m_EionSum += eion;
210 
211  // if any of these conditions is true this is the last step in
212  // this volume and we need to save the hit
213  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
214  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
215  // (happens when your detector goes outside world volume)
216  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
217  // aTrack->GetTrackStatus() == fStopAndKill is also set)
218  // aTrack->GetTrackStatus() == fStopAndKill: track ends
219  if (postPoint->GetStepStatus() == fGeomBoundary ||
220  postPoint->GetStepStatus() == fWorldBoundary ||
221  postPoint->GetStepStatus() == fAtRestDoItProc ||
222  aTrack->GetTrackStatus() == fStopAndKill)
223  {
224  // save only hits with energy deposit (or geantino)
225  if (m_EdepSum > 0 || geantino)
226  {
227  // update values at exit coordinates and set keep flag
228  // of track to keep
229  m_hit->set_x(1, postPoint->GetPosition().x() / cm);
230  m_hit->set_y(1, postPoint->GetPosition().y() / cm);
231  m_hit->set_z(1, postPoint->GetPosition().z() / cm);
232 
233  m_hit->set_px(1, postPoint->GetMomentum().x() / GeV);
234  m_hit->set_py(1, postPoint->GetMomentum().y() / GeV);
235  m_hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
236 
237  m_hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
239  {
240  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
241  {
242  pp->SetKeep(1); // we want to keep the track
243  }
244  }
245  if (geantino)
246  {
247  m_hit->set_edep(-1);
248  m_hit->set_eion(-1);
249  }
250  else
251  {
252  m_hit->set_edep(m_EdepSum);
253  m_hit->set_eion(m_EionSum);
254  }
255 
256 
257  // add in container
258  m_SaveHitContainer->AddHit(m_hit->get_layer(), m_hit.get());
259 
260  // ownership is transferred to container
261  // so we release the hit
262  m_hit.release();
263 
264  } else {
265 
266  m_hit->Reset();
267 
268  }
269  }
270  // return true to indicate the hit was used
271  return true;
272 }
273 
274 //____________________________________________________________________________..
276 {
277  std::string hitnodename = "G4HIT_" + m_Detector->SuperDetector();
278  m_hitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
279  if (!m_hitContainer) std::cout << "PHG4MicromegasSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
280 }