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