ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TRDSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TRDSteppingAction.cc
2 #include "PHG4TRDDetector.h"
3 #include "PHG4TRDSubsystem.h"
4 
5 //#include "PHG4StepStatusDecode.h"
6 
7 #include <phparameter/PHParameters.h>
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 <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
19 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
20 #include <Geant4/G4Step.hh>
21 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
22 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fPostSt...
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 <boost/io/ios_state.hpp>
34 
35 #include <cmath> // for isfinite, copysign
36 #include <cstdlib> // for exit
37 #include <iomanip>
38 #include <iostream>
39 #include <string> // for operator<<, char_traits
40 
41 class PHCompositeNode;
42 
43 using namespace std;
44 
45 //____________________________________________________________________________..
46 //PHG4TRDSteppingAction::PHG4TRDSteppingAction(PHG4TRDSubsystem *subsys, PHG4TRDDetector* detector, const PHParameters* parameters)
48  : PHG4SteppingAction(detector->GetName())
49  , m_Detector(detector)
50  , m_Params(parameters)
51  , m_HitContainer(nullptr)
52  , m_ActiveGasHits(nullptr)
53  , m_Hit(nullptr)
54  , m_SaveHitContainer(nullptr)
55  , m_SaveShower(nullptr)
56  , m_SaveVolPre(nullptr)
57  , m_SaveVolPost(nullptr)
58  , m_SaveTrackId(-1)
59  , m_SavePreStepStatus(-1)
60  , m_SavePostStepStatus(-1)
61  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
62  , m_ActiveFlag(m_Params->get_int_param("active"))
63  , m_UseG4StepsFlag(m_Params->get_int_param("use_g4steps"))
64  , m_Zmin(m_Params->get_double_param("PosZ") * cm - m_Params->get_double_param("ThicknessZ") * cm / 2.)
65  , m_Zmax(m_Params->get_double_param("PosZ") * cm + m_Params->get_double_param("ThicknessZ") * cm / 2.)
66  //, m_EdepSum(0)
67 {
68  m_Zmin -= copysign(m_Zmin, 1. / 1e6 * cm);
69  m_Zmax += copysign(m_Zmax, 1. / 1e6 * cm);
70 }
71 
73 {
74  // if the last hit was a zero energie deposit hit, it is just reset
75  // and the memory is still allocated, so we need to delete it here
76  // if the last hit was saved, hit is a nullptr pointer which are
77  // legal to delete (it results in a no operation)
78  delete m_Hit;
79 }
80 
81 //____________________________________________________________________________..
83 {
84  // get volume of the current step
87 
89  // G4 just calls UserSteppingAction for every step (and we loop then over all our
90  // steppingactions. First we have to check if we are actually in our volume
91  if (!m_Detector->IsInTRD(volume))
92  {
93  return false;
94  }
95 
96  int whichactive = m_Detector->IsInTRD(volume);
97 
98  // collect energy and track length step by step
100 
101  const G4Track* aTrack = aStep->GetTrack();
102  // if this volume stops everything, just put all kinetic energy into edep
103  if (m_BlackHoleFlag)
104  {
105  edep = aTrack->GetKineticEnergy() / GeV;
106  G4Track* killtrack = const_cast<G4Track*>(aTrack);
107  killtrack->SetTrackStatus(fStopAndKill);
108  }
109  int layer_id = m_Detector->get_Layer();
110 
111  if (m_ActiveFlag)
112  {
113  bool geantino = false;
114  // the check for the pdg code speeds things up, I do not want to make
115  // an expensive string compare for every track when we know
116  // geantino or chargedgeantino has pid=0
117  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
118  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
119  {
120  geantino = true;
121  }
122  G4StepPoint* prePoint = aStep->GetPreStepPoint();
123  G4StepPoint* postPoint = aStep->GetPostStepPoint();
124  // std::cout << "time prepoint: " << prePoint->GetGlobalTime()/ns << std::endl;
125  // std::cout << "time postpoint: " << postPoint->GetGlobalTime()/ns << std::endl;
126  // std::cout << "kinetic energy: " << aTrack->GetKineticEnergy()/GeV << std::endl;
127  // G4ParticleDefinition* def = aTrack->GetDefinition();
128  // std::cout << "Particle: " << def->GetParticleName() << std::endl;
129  int prepointstatus = prePoint->GetStepStatus();
130  if (prepointstatus == fGeomBoundary ||
131  prepointstatus == fUndefined ||
132  (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary) ||
133  m_UseG4StepsFlag > 0)
134  {
135  if (!m_Hit)
136  {
137  m_Hit = new PHG4Hitv1();
138  }
139 
140  m_Hit->set_layer((unsigned int) layer_id);
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 
151  // time in ns
152  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
153  //set and save the track ID
154  m_Hit->set_trkid(aTrack->GetTrackID());
155  m_SaveTrackId = aTrack->GetTrackID();
156  //set the initial energy deposit
157  m_Hit->set_edep(0);
158 
159  //cout << " which active : " << whichactive << " layer ID :" << layer_id << " hit id : " << aTrack->GetTrackID() <<" vol name : " << volume->GetName() <<" z : " << prePoint->GetPosition().z()/cm << " end z : " << postPoint->GetPosition().z()/cm << endl;
160  /*
161  if(whichactive > 0 )
162  {
163  m_SaveHitContainer = m_AbsorberHits;
164  }
165  else
166  {
167  m_SaveHitContainer = m_HitContainer;
168  }
169  */
170  if (whichactive < 0)
171  {
173  }
174  else
175  {
176  //m_SaveHitContainer = m_AbsorberHits;
178  }
179 
180  if (!geantino && !m_BlackHoleFlag)
181  {
182  m_Hit->set_eion(0);
183  }
185  {
186  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
187  {
188  m_Hit->set_trkid(pp->GetUserTrackId());
189  m_Hit->set_shower_id(pp->GetShower()->get_id());
190  m_SaveShower = pp->GetShower();
191  }
192  }
193  /*
194  if (!hasMotherSubsystem() && (m_Hit->get_z(0) * cm > m_Zmax || m_Hit->get_z(0) * cm < m_Zmin))
195  {
196  boost::io::ios_precision_saver ips(std::cout);
197  std::cout << m_Detector->SuperDetector() << std::setprecision(9)
198  << " PHG4TRDSteppingAction: Entry hit z " << m_Hit->get_z(0) * cm
199  << " outside acceptance, zmin " << m_Zmin
200  << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
201  }
202  */
203  } // END ..... || m_UseG4StepsFlag > 0)
204  /*
205  // here we just update the exit values, it will be overwritten
206  // for every step until we leave the volume or the particle
207  // ceases to exist
208  // some sanity checks for inconsistencies
209  // check if this hit was created, if not print out last post step status
210  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
211  {
212  std::cout << GetName() << ": hit was not created" << std::endl;
213  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
214  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
215  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
216  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
217  std::cout << "last track: " << m_SaveTrackId
218  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
219  std::cout << "phys pre vol: " << volume->GetName()
220  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
221  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
222  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
223  exit(1);
224  }
225  */
226  m_SavePostStepStatus = postPoint->GetStepStatus();
227  // check if track id matches the initial one when the hit was created
228  if (aTrack->GetTrackID() != m_SaveTrackId)
229  {
230  std::cout << "hits do not belong to the same track" << std::endl;
231  std::cout << "saved track: " << m_SaveTrackId
232  << ", current trackid: " << aTrack->GetTrackID()
233  << std::endl;
234  exit(1);
235  }
236  m_SavePreStepStatus = prePoint->GetStepStatus();
237  m_SavePostStepStatus = postPoint->GetStepStatus();
239  m_SaveVolPost = touchpost->GetVolume();
240 
241  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
242  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
243  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
244 
245  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
246  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
247  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
248 
249  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
250  //sum up the energy to get total deposited
252  /*
253  if (!hasMotherSubsystem() && (m_Hit->get_z(1) * cm > m_Zmax || m_Hit->get_z(1) * cm < m_Zmin))
254  {
255  std::cout << m_Detector->SuperDetector() << std::setprecision(9)
256  << " PHG4TRDSteppingAction: Exit hit z " << m_Hit->get_z(1) * cm
257  << " outside acceptance zmin " << m_Zmin
258  << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
259  }
260  */
261  if (geantino)
262  {
263  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
264  }
265  else
266  {
267  if (!m_BlackHoleFlag)
268  {
269  double eion = edep - aStep->GetNonIonizingEnergyDeposit() / GeV;
270  m_Hit->set_eion(m_Hit->get_eion() + eion);
271  }
272  }
273 
274  // if any of these conditions is true this is the last step in
275  // this volume and we need to save the hit
276  if (postPoint->GetStepStatus() == fGeomBoundary ||
277  postPoint->GetStepStatus() == fWorldBoundary ||
278  postPoint->GetStepStatus() == fAtRestDoItProc ||
279  aTrack->GetTrackStatus() == fStopAndKill ||
280  m_UseG4StepsFlag > 0)
281  {
282  // save only hits with energy deposit (or -1 for geantino) or if save all hits flag is set
283  if (m_Hit->get_edep())
284  {
285  m_SaveHitContainer->AddHit(layer_id, m_Hit);
286 
287  if (m_SaveShower)
288  {
290  }
291  // ownership has been transferred to container, set to null
292  // so we will create a new hit for the next track
293  m_Hit = nullptr;
294  }
295  else
296  {
297  // if this hit has no energy deposit, just reset it for reuse
298  // this means we have to delete it in the dtor. If this was
299  // the last hit we processed the memory is still allocated
300  m_Hit->Reset();
301  }
302  }
303  // return true to indicate the hit was used
304  return true;
305  } // END Acitve flag condition
306  else
307  {
308  return false;
309  }
310 }
311 
312 //____________________________________________________________________________..
314 {
315  /*
316  // Node Name is passed down from PHG4TRDSubsystem
317  //now look for the map and grab a pointer to it.
318  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
319 
320  // if we do not find the node we need to scream.
321  if (!m_HitContainer && !m_BlackHoleFlag)
322  {
323  std::cout << "PHG4TRDSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
324  }
325  */
326 
327  string hitnodename;
328  string gasnodename;
329  if (m_Detector->SuperDetector() != "NONE")
330  {
331  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
332  gasnodename = "G4HIT_ACTIVEGAS_" + m_Detector->SuperDetector();
333  }
334  else
335  {
336  hitnodename = "G4HIT_" + m_Detector->GetName();
337  gasnodename = "G4HIT_ACTIVEGAS_" + m_Detector->SuperDetector();
338  }
339 
340  //now look for the map and grab a pointer to it.
341  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
342  m_ActiveGasHits = findNode::getClass<PHG4HitContainer>(topNode, gasnodename.c_str());
343 
344  // if we do not find the node it's messed up.
345  if (!m_HitContainer)
346  {
347  std::cout << "PHG4TRDSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
348  }
349  // if (!m_AbsorberHits)
350  if (!m_ActiveGasHits)
351  {
352  if (Verbosity() > 1)
353  {
354  cout << "PHG4TRDSteppingAction::SetTopNode - unable to find " << gasnodename << endl;
355  }
356  }
357 }
358 
359 /*
360 bool PHG4TRDSteppingAction::hasMotherSubsystem() const
361 {
362  if (m_Subsystem->GetMotherSubsystem())
363  {
364  return true;
365  }
366  return false;
367 }
368 */