ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4CylinderStripSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4CylinderStripSteppingAction.cc
4 
5 #include <phparameter/PHParameters.h>
6 
7 #include <g4main/PHG4Hit.h>
9 #include <g4main/PHG4Hitv1.h>
10 #include <g4main/PHG4Shower.h>
11 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
12 
14 
15 #include <phool/getClass.h>
16 
17 #include <TSystem.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 using namespace std;
46 //____________________________________________________________________________..
48  : PHG4SteppingAction(detector->GetName())
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->IsInTileC(volume)) && (!m_Detector->IsInTileZ(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 ((!isfinite(m_Tmin) && !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") != string::npos)
125  {
126  geantino = true;
127  }
128  G4StepPoint* prePoint = aStep->GetPreStepPoint();
129  G4StepPoint* postPoint = aStep->GetPostStepPoint();
130  // cout << "time prepoint: " << prePoint->GetGlobalTime()/ns << endl;
131  // cout << "time postpoint: " << postPoint->GetGlobalTime()/ns << endl;
132  // cout << "kinetic energy: " << aTrack->GetKineticEnergy()/GeV << endl;
133  // G4ParticleDefinition* def = aTrack->GetDefinition();
134  // cout << "Particle: " << def->GetParticleName() << 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  cout << GetName() << ": New Hit for " << endl;
144  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) << endl;
148  cout << "last track: " << m_SaveTrackId
149  << ", current trackid: " << aTrack->GetTrackID() << endl;
150  cout << "phys pre vol: " << volume->GetName()
151  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
152  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
153  << " previous phys post vol: " << m_SaveVolPost->GetName() << 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  if (m_Detector->IsInTileC(volume)) m_Hit->set_hit_type(1);
177  if (m_Detector->IsInTileZ(volume)) m_Hit->set_hit_type(2);
178  m_SaveTrackId = aTrack->GetTrackID();
179  //set the initial energy deposit
180  m_Hit->set_edep(0);
181  if (!geantino && !m_BlackHoleFlag)
182  {
183  m_Hit->set_eion(0);
184  }
186  {
188  }
190  {
191  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
192  {
193  m_Hit->set_trkid(pp->GetUserTrackId());
194  m_Hit->set_shower_id(pp->GetShower()->get_id());
195  m_SaveShower = pp->GetShower();
196  }
197  }
198 
199  if (m_Hit->get_z(0) * cm > m_Zmax || m_Hit->get_z(0) * cm < m_Zmin)
200  {
201  boost::io::ios_precision_saver ips(cout);
202  cout << m_Detector->SuperDetector() << std::setprecision(9)
203  << "PHG4CylinderStripSteppingAction: Entry hit z " << m_Hit->get_z(0) * cm
204  << " outside acceptance, zmin " << m_Zmin
205  << ", zmax " << m_Zmax << ", layer: " << layer_id << endl;
206  }
207  }
208  // here we just update the exit values, it will be overwritten
209  // for every step until we leave the volume or the particle
210  // ceases to exist
211  // some sanity checks for inconsistencies
212  // check if this hit was created, if not print out last post step status
213  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
214  {
215  cout << GetName() << ": hit was not created" << endl;
216  cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
217  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
218  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
219  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
220  cout << "last track: " << m_SaveTrackId
221  << ", current trackid: " << aTrack->GetTrackID() << endl;
222  cout << "phys pre vol: " << volume->GetName()
223  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
224  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
225  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
226  gSystem->Exit(1);
227  }
228  m_SavePostStepStatus = postPoint->GetStepStatus();
229  // check if track id matches the initial one when the hit was created
230  if (aTrack->GetTrackID() != m_SaveTrackId)
231  {
232  cout << "hits do not belong to the same track" << endl;
233  cout << "saved track: " << m_SaveTrackId
234  << ", current trackid: " << aTrack->GetTrackID()
235  << endl;
236  gSystem->Exit(1);
237  }
238  m_SavePreStepStatus = prePoint->GetStepStatus();
239  m_SavePostStepStatus = postPoint->GetStepStatus();
241  m_SaveVolPost = touchpost->GetVolume();
242 
243  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
244  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
245  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
246 
247  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
248  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
249  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
250 
251  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
252  //sum up the energy to get total deposited
254  if (m_Hit->get_z(1) * cm > m_Zmax || m_Hit->get_z(1) * cm < m_Zmin)
255  {
256  cout << m_Detector->SuperDetector() << std::setprecision(9)
257  << " PHG4CylinderStripSteppingAction: Exit hit z " << m_Hit->get_z(1) * cm
258  << " outside acceptance zmin " << m_Zmin
259  << ", zmax " << m_Zmax << ", layer: " << layer_id << endl;
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  }
274  {
275  double light_yield = GetVisibleEnergyDeposition(aStep);
276  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
277  }
278  if (edep > 0)
279  {
281  {
282  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
283  {
284  pp->SetKeep(1); // we want to keep the track
285  }
286  }
287  }
288  // if any of these conditions is true this is the last step in
289  // this volume and we need to save the hit
290  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
291  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
292  // (happens when your detector goes outside world volume)
293  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
294  // aTrack->GetTrackStatus() == fStopAndKill is also set)
295  // aTrack->GetTrackStatus() == fStopAndKill: track ends
296  if (postPoint->GetStepStatus() == fGeomBoundary ||
297  postPoint->GetStepStatus() == fWorldBoundary ||
298  postPoint->GetStepStatus() == fAtRestDoItProc ||
299  aTrack->GetTrackStatus() == fStopAndKill ||
300  m_UseG4StepsFlag > 0)
301  {
302  // save only hits with energy deposit (or -1 for geantino)
303  if (m_Hit->get_edep())
304  {
305  m_HitContainer->AddHit(layer_id, m_Hit);
306  if (m_SaveShower)
307  {
309  }
310  // ownership has been transferred to container, set to null
311  // so we will create a new hit for the next track
312  m_Hit = nullptr;
313  }
314  else
315  {
316  // if this hit has no energy deposit, just reset it for reuse
317  // this means we have to delete it in the dtor. If this was
318  // the last hit we processed the memory is still allocated
319  m_Hit->Reset();
320  }
321  }
322  // return true to indicate the hit was used
323  return true;
324  }
325  else
326  {
327  return false;
328  }
329 }
330 
331 //____________________________________________________________________________..
333 {
334  string hitnodename;
335  if (m_Detector->SuperDetector() != "NONE")
336  {
337  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
338  }
339  else
340  {
341  hitnodename = "G4HIT_" + m_Detector->GetName();
342  }
343 
344  //now look for the map and grab a pointer to it.
345  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
346 
347  // if we do not find the node we need to make it.
349  {
350  std::cout << "PHG4CylinderStripSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
351  }
352 }