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