ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4CylinderSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4CylinderSteppingAction.cc
2 #include "PHG4CylinderDetector.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
14 
16 
17 #include <phool/getClass.h>
18 
19 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
20 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
21 #include <Geant4/G4Step.hh>
22 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
23 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fPostSt...
24 #include <Geant4/G4String.hh> // for G4String
25 #include <Geant4/G4SystemOfUnits.hh>
26 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
27 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
28 #include <Geant4/G4Track.hh> // for G4Track
29 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
30 #include <Geant4/G4Types.hh> // for G4double
31 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
32 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
33 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
34 
35 #include <boost/io/ios_state.hpp>
36 
37 #include <cmath> // for isfinite, copysign
38 #include <cstdlib> // for exit
39 #include <iomanip>
40 #include <iostream>
41 #include <string> // for operator<<, char_traits
42 
43 class PHCompositeNode;
44 
45 //____________________________________________________________________________..
47  : PHG4SteppingAction(detector->GetName())
48  , m_Subsystem(subsys)
49  , m_Detector(detector)
50  , m_Params(parameters)
51  , m_HitContainer(nullptr)
52  , m_Hit(nullptr)
53  , m_SaveShower(nullptr)
54  , m_SaveVolPre(nullptr)
55  , m_SaveVolPost(nullptr)
56  , m_SaveLightYieldFlag(m_Params->get_int_param("lightyield"))
57  , m_SaveTrackId(-1)
58  , m_SavePreStepStatus(-1)
59  , m_SavePostStepStatus(-1)
60  , m_ActiveFlag(m_Params->get_int_param("active"))
61  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
62  , m_UseG4StepsFlag(m_Params->get_int_param("use_g4steps"))
63  , m_Zmin(m_Params->get_double_param("place_z") * cm - m_Params->get_double_param("length") * cm / 2.)
64  , m_Zmax(m_Params->get_double_param("place_z") * cm + m_Params->get_double_param("length") * cm / 2.)
65  , m_Tmin(m_Params->get_double_param("tmin") * ns)
66  , m_Tmax(m_Params->get_double_param("tmax") * ns)
67 {
68  // G4 seems to have issues in the um range
69  m_Zmin -= copysign(m_Zmin, 1. / 1e6 * cm);
70  m_Zmax += copysign(m_Zmax, 1. / 1e6 * cm);
71 }
72 
74 {
75  // if the last hit was a zero energie deposit hit, it is just reset
76  // and the memory is still allocated, so we need to delete it here
77  // if the last hit was saved, hit is a nullptr pointer which are
78  // legal to delete (it results in a no operation)
79  delete m_Hit;
80 }
81 
82 //____________________________________________________________________________..
84 {
85  // get volume of the current step
88 
90  // G4 just calls UserSteppingAction for every step (and we loop then over all our
91  // steppingactions. First we have to check if we are actually in our volume
92  if (!m_Detector->IsInCylinder(volume))
93  {
94  return false;
95  }
96 
97  // collect energy and track length step by step
99 
100  const G4Track* aTrack = aStep->GetTrack();
101 
102  // if this cylinder stops everything, just put all kinetic energy into edep
103  if (m_BlackHoleFlag)
104  {
105  if ((!std::isfinite(m_Tmin) && !std::isfinite(m_Tmax)) ||
106  aTrack->GetGlobalTime() < m_Tmin ||
107  aTrack->GetGlobalTime() > m_Tmax)
108  {
109  edep = aTrack->GetKineticEnergy() / GeV;
110  G4Track* killtrack = const_cast<G4Track*>(aTrack);
111  killtrack->SetTrackStatus(fStopAndKill);
112  }
113  }
114 
115  int layer_id = m_Detector->get_Layer();
116  // test if we are active
117  if (m_ActiveFlag)
118  {
119  bool geantino = false;
120  // the check for the pdg code speeds things up, I do not want to make
121  // an expensive string compare for every track when we know
122  // geantino or chargedgeantino has pid=0
123  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
124  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
125  {
126  geantino = true;
127  }
128  G4StepPoint* prePoint = aStep->GetPreStepPoint();
129  G4StepPoint* postPoint = aStep->GetPostStepPoint();
130  // std::cout << "time prepoint: " << prePoint->GetGlobalTime()/ns << std::endl;
131  // std::cout << "time postpoint: " << postPoint->GetGlobalTime()/ns << std::endl;
132  // std::cout << "kinetic energy: " << aTrack->GetKineticEnergy()/GeV << std::endl;
133  // G4ParticleDefinition* def = aTrack->GetDefinition();
134  // std::cout << "Particle: " << def->GetParticleName() << std::endl;
135  int prepointstatus = prePoint->GetStepStatus();
136  if (prepointstatus == fGeomBoundary ||
137  prepointstatus == fUndefined ||
138  (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary) ||
139  m_UseG4StepsFlag > 0)
140  {
141  if (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary)
142  {
143  std::cout << GetName() << ": New Hit for " << std::endl;
144  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
145  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
146  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
147  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
148  std::cout << "last track: " << m_SaveTrackId
149  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
150  std::cout << "phys pre vol: " << volume->GetName()
151  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
152  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
153  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
154  }
155 
156  if (!m_Hit)
157  {
158  m_Hit = new PHG4Hitv1();
159  }
160 
161  m_Hit->set_layer((unsigned int) layer_id);
162 
163  //here we set the entrance values in cm
164  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
165  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
166  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
167 
168  m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
169  m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
170  m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
171 
172  // time in ns
173  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
174  //set and save the track ID
175  m_Hit->set_trkid(aTrack->GetTrackID());
176  m_SaveTrackId = aTrack->GetTrackID();
177  //set the initial energy deposit
178  m_Hit->set_edep(0);
179  if (!geantino && !m_BlackHoleFlag)
180  {
181  m_Hit->set_eion(0);
182  }
184  {
186  }
188  {
189  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
190  {
191  m_Hit->set_trkid(pp->GetUserTrackId());
192  m_Hit->set_shower_id(pp->GetShower()->get_id());
193  m_SaveShower = pp->GetShower();
194  }
195  }
196 
197  if (!hasMotherSubsystem() && (m_Hit->get_z(0) * cm > m_Zmax || m_Hit->get_z(0) * cm < m_Zmin))
198  {
199  boost::io::ios_precision_saver ips(std::cout);
200  std::cout << m_Detector->SuperDetector() << std::setprecision(9)
201  << " PHG4CylinderSteppingAction: Entry hit z " << m_Hit->get_z(0) * cm
202  << " outside acceptance, zmin " << m_Zmin
203  << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
204  }
205  }
206  // here we just update the exit values, it will be overwritten
207  // for every step until we leave the volume or the particle
208  // ceases to exist
209  // some sanity checks for inconsistencies
210  // check if this hit was created, if not print out last post step status
211  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
212  {
213  std::cout << GetName() << ": hit was not created" << std::endl;
214  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
215  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
216  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
217  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
218  std::cout << "last track: " << m_SaveTrackId
219  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
220  std::cout << "phys pre vol: " << volume->GetName()
221  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
222  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
223  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
224  exit(1);
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  if (!hasMotherSubsystem() && (m_Hit->get_z(1) * cm > m_Zmax || m_Hit->get_z(1) * cm < m_Zmin))
253  {
254  std::cout << m_Detector->SuperDetector() << std::setprecision(9)
255  << " PHG4CylinderSteppingAction: Exit hit z " << m_Hit->get_z(1) * cm
256  << " outside acceptance zmin " << m_Zmin
257  << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
258  }
259  if (geantino)
260  {
261  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
262  }
263  else
264  {
265  if (!m_BlackHoleFlag)
266  {
267  double eion = edep - aStep->GetNonIonizingEnergyDeposit() / GeV;
268  m_Hit->set_eion(m_Hit->get_eion() + eion);
269  }
270  }
272  {
273  double light_yield = GetVisibleEnergyDeposition(aStep);
274  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
275  }
276  if (edep > 0 || m_SaveAllHitsFlag)
277  {
279  {
280  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
281  {
282  pp->SetKeep(1); // we want to keep the track
283  }
284  }
285  }
286  // if any of these conditions is true this is the last step in
287  // this volume and we need to save the hit
288  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
289  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
290  // (happens when your detector goes outside world volume)
291  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
292  // aTrack->GetTrackStatus() == fStopAndKill is also set)
293  // aTrack->GetTrackStatus() == fStopAndKill: track ends
294  if (postPoint->GetStepStatus() == fGeomBoundary ||
295  postPoint->GetStepStatus() == fWorldBoundary ||
296  postPoint->GetStepStatus() == fAtRestDoItProc ||
297  aTrack->GetTrackStatus() == fStopAndKill ||
298  m_UseG4StepsFlag > 0)
299  {
300  // save only hits with energy deposit (or -1 for geantino) or if save all hits flag is set
301  if (m_Hit->get_edep() || m_SaveAllHitsFlag)
302  {
303  m_HitContainer->AddHit(layer_id, m_Hit);
304  if (m_SaveShower)
305  {
307  }
308  // ownership has been transferred to container, set to null
309  // so we will create a new hit for the next track
310  m_Hit = nullptr;
311  }
312  else
313  {
314  // if this hit has no energy deposit, just reset it for reuse
315  // this means we have to delete it in the dtor. If this was
316  // the last hit we processed the memory is still allocated
317  m_Hit->Reset();
318  }
319  }
320  // return true to indicate the hit was used
321  return true;
322  }
323  else
324  {
325  return false;
326  }
327 }
328 
329 //____________________________________________________________________________..
331 {
332  // Node Name is passed down from PHG4CylinderSubsystem
333  //now look for the map and grab a pointer to it.
334  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
335 
336  // if we do not find the node we need to scream.
338  {
339  std::cout << "PHG4CylinderSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
340  }
341 }
342 
344 {
346  {
347  return true;
348  }
349  return false;
350 }