ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EicDircSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EicDircSteppingAction.cc
2 #include "G4EicDircDetector.h"
3 #include "PrtHit.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 #include <TVector3.h>
20 
21 #include <Geant4/G4NavigationHistory.hh>
22 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
23 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
24 #include <Geant4/G4Step.hh>
25 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
26 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
27 #include <Geant4/G4String.hh> // for G4String
28 #include <Geant4/G4SystemOfUnits.hh>
29 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
30 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
31 #include <Geant4/G4Track.hh> // for G4Track
32 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
33 #include <Geant4/G4TransportationManager.hh>
34 #include <Geant4/G4Types.hh> // for G4double
35 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
36 #include <Geant4/G4VProcess.hh>
37 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
38 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
39 #include <Geant4/Randomize.hh>
40 
41 #include <cmath> // for isfinite
42 #include <iostream>
43 #include <string> // for operator<<, string
44 
45 class PHCompositeNode;
46 //____________________________________________________________________________..
48  G4EicDircDetector* detector, const PHParameters* parameters)
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 
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 //____________________________________________________________________________..
68  bool was_used)
69 {
72  // get volume of the current step
74  //G4VPhysicalVolume* volume_post = touchpost->GetVolume();
75  // IsInDetector(volume) returns
76  // == 0 outside of detector
77  // > 0 for hits in active volume
78  // < 0 for hits in passive material
79  int whichactive_int = m_Detector->IsInDetector(volume);
80  //int whichactive_int_post = m_Detector->IsInDetector(volume_post);
81  bool whichactive = (whichactive_int > 0 && whichactive_int < 12);
82 
83  //int whichactive = m_Detector->IsInDetector(volume);
84  if (!whichactive)
85  {
86  return false;
87  }
88 
89  // collect energy and track length step by step
90  const G4Track* aTrack = aStep->GetTrack();
92  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
93 
94  // if this block stops everything, just put all kinetic energy into edep
95  if (m_BlackHoleFlag)
96  {
97  edep = aTrack->GetKineticEnergy() / GeV;
98  G4Track* killtrack = const_cast<G4Track*>(aTrack);
99  killtrack->SetTrackStatus(fStopAndKill);
100  }
101 
102  // make sure we are in a volume
103  if (m_ActiveFlag)
104  {
105  bool geantino = false;
106  // the check for the pdg code speeds things up, I do not want to make
107  // an expensive string compare for every track when we know
108  // geantino or chargedgeantino has pid=0
109  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
110  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
111  {
112  geantino = true;
113  }
114  G4StepPoint* prePoint = aStep->GetPreStepPoint();
115  G4StepPoint* postPoint = aStep->GetPostStepPoint();
116  // std::cout << "track id " << aTrack->GetTrackID() << std::endl;
117  // std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
118  // std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
119  int prepointstatus = prePoint->GetStepStatus();
120  if (prepointstatus == fGeomBoundary ||
121  prepointstatus == fUndefined ||
122  (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary))
123  {
124  if (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary)
125  {
126  std::cout << GetName() << ": New Hit for " << std::endl;
127  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
128  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
129  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
130  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
131  std::cout << "last track: " << m_SaveTrackId
132  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
133  std::cout << "phys pre vol: " << volume->GetName()
134  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
135  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
136  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
137  }
138  if (!m_Hit)
139  {
140  m_Hit = new PHG4Hitv1();
141  }
142  //here we set the entrance values in cm
143  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
144  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
145  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
146 
147  m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
148  m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
149  m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
150  // time in ns
151  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
152  //set the track ID
153  m_Hit->set_trkid(aTrack->GetTrackID());
154  m_SaveTrackId = aTrack->GetTrackID();
155 
156  //set the initial energy deposit
157  m_Hit->set_edep(0);
158  if (!geantino && !m_BlackHoleFlag && m_ActiveFlag)
159  {
160  m_Hit->set_eion(0);
161  }
163  {
164  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
165  {
166  m_Hit->set_trkid(pp->GetUserTrackId());
167  m_Hit->set_shower_id(pp->GetShower()->get_id());
168  }
169  }
170  }
171 
172  // some sanity checks for inconsistencies
173  // check if this hit was created, if not print out last post step status
174  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
175  {
176  std::cout << GetName() << ": hit was not created" << std::endl;
177  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
178  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
179  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
180  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
181  std::cout << "last track: " << m_SaveTrackId
182  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
183  std::cout << "phys pre vol: " << volume->GetName()
184  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
185  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
186  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
187  exit(1);
188  }
189  m_SavePostStepStatus = postPoint->GetStepStatus();
190  // check if track id matches the initial one when the hit was created
191  if (aTrack->GetTrackID() != m_SaveTrackId)
192  {
193  std::cout << "hits do not belong to the same track" << std::endl;
194  std::cout << "saved track: " << m_SaveTrackId
195  << ", current trackid: " << aTrack->GetTrackID()
196  << std::endl;
197  exit(1);
198  }
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  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
208  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
209  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
210 
211  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
212  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
213  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
214 
215  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
216  //sum up the energy to get total deposited
218  // update ionization energy only for active volumes, not for black holes or geantinos
219  // if the hit is created without eion, get_eion() returns NAN
220  // if you see NANs check the creation of the hit
221  if (m_ActiveFlag && !m_BlackHoleFlag && !geantino)
222  {
223  m_Hit->set_eion(m_Hit->get_eion() + eion);
224  }
225  if (geantino)
226  {
227  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
228  }
229  if (edep > 0)
230  {
232  {
233  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
234  {
235  pp->SetKeep(1); // we want to keep the track
236  }
237  }
238  }
239  // if any of these conditions is true this is the last step in
240  // this volume and we need to save the hit
241  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
242  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
243  // (happens when your detector goes outside world volume)
244  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
245  // aTrack->GetTrackStatus() == fStopAndKill is also set)
246  // aTrack->GetTrackStatus() == fStopAndKill: track ends
247  if (postPoint->GetStepStatus() == fGeomBoundary ||
248  postPoint->GetStepStatus() == fWorldBoundary ||
249  postPoint->GetStepStatus() == fAtRestDoItProc ||
250  aTrack->GetTrackStatus() == fStopAndKill)
251  {
252  // save only hits with energy deposit (or -1 for geantino)
253  if (m_Hit->get_edep())
254  {
256  // ownership has been transferred to container, set to null
257  // so we will create a new hit for the next track
258  m_Hit = nullptr;
259  }
260  else
261  {
262  // if this hit has no energy deposit, just reset it for reuse
263  // this means we have to delete it in the dtor. If this was
264  // the last hit we processed the memory is still allocated
265  m_Hit->Reset();
266  }
267  }
268  // return true to indicate the hit was used
269  return true;
270  }
271  else
272  {
273  return false;
274  }
275 }
276 
277 //____________________________________________________________________________..
279 {
280  if (!m_HitNodeName.empty())
281  {
282  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
283  }
284  if (!m_AbsorberNodeName.empty())
285  {
286  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
288  {
289  if (Verbosity() > 0)
290  {
291  std::cout << "G4EicDircSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
292  }
293  }
294  }
295  if (!m_SupportNodeName.empty())
296  {
297  m_SupportHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_SupportNodeName);
299  {
300  if (Verbosity() > 0)
301  {
302  std::cout << "G4EicDircSteppingAction::SetTopNode - unable to find " << m_SupportNodeName << std::endl;
303  }
304  }
305  }
306  if (!m_HitContainer)
307  {
308  std::cout << "G4EicDircSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
309  }
310 }