Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TpcSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TpcSteppingAction.cc
2 #include "PHG4TpcDetector.h"
3 
5 
6 #include <phparameter/PHParameters.h>
7 
8 #include <g4main/PHG4Hit.h>
10 #include <g4main/PHG4Hitv1.h>
11 #include <g4main/PHG4Shower.h>
12 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
13 
15 
16 #include <phool/getClass.h>
17 
18 #include <TSystem.h>
19 
20 #include <Geant4/G4ParticleDefinition.hh>
21 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
22 #include <Geant4/G4Step.hh>
23 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
24 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fPostSt...
25 #include <Geant4/G4String.hh> // for G4String
26 #include <Geant4/G4SystemOfUnits.hh>
27 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
28 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
29 #include <Geant4/G4Track.hh> // for G4Track
30 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
31 #include <Geant4/G4Types.hh> // for G4double
32 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
33 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
34 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
35 
36 #include <cmath> // for isfinite
37 #include <cstdlib> // for exit
38 #include <iostream>
39 #include <string> // for operator<<, operator+
40 
41 class PHCompositeNode;
42 
43 //____________________________________________________________________________..
45  : PHG4SteppingAction(detector->GetName())
46  , m_Detector(detector)
47  , m_Params(parameters)
48  , m_IsBlackHoleFlag(m_Params->get_int_param("blackhole"))
49 {
50  if (std::isfinite(m_Params->get_double_param("steplimits")))
51  {
52  m_UseG4StepsFlag = 1;
53  }
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 bool PHG4TpcSteppingAction::UserSteppingAction(const G4Step* aStep, bool)
67 {
68  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
69  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
70  // get volume of the current step
71  G4VPhysicalVolume* volume = touch->GetVolume();
72 
73  // m_Detector->IsInTpc(volume)
74  // returns
75  // 0 is outside of Tpc
76  // 1 is inside tpc gas
77  // <0 is in tpc support structure (cage, endcaps,...)
78 
79  int whichactive = m_Detector->IsInTpc(volume);
80  if (!whichactive)
81  {
82  return false;
83  }
84  // collect energy and track length step by step
85  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
86  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
87  const G4Track* aTrack = aStep->GetTrack();
88 
89  // if this block stops everything, just put all kinetic energy into edep
91  {
92  edep = aTrack->GetKineticEnergy() / GeV;
93  G4Track* killtrack = const_cast<G4Track*>(aTrack);
94  killtrack->SetTrackStatus(fStopAndKill);
95  }
96 
97  bool geantino = false;
98 
99  // the check for the pdg code speeds things up, I do not want to make
100  // an expensive string compare for every track when we know
101  // geantino or chargedgeantino has pid=0
102  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
103  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
104  {
105  geantino = true;
106  }
107  G4StepPoint* prePoint = aStep->GetPreStepPoint();
108  G4StepPoint* postPoint = aStep->GetPostStepPoint();
109  int prepointstatus = prePoint->GetStepStatus();
110 
111  // std::cout << "track id " << aTrack->GetTrackID() << std::endl;
112  // std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
113  // std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
114 
115  if ((m_UseG4StepsFlag > 0 && whichactive > 0) ||
116  prepointstatus == fGeomBoundary ||
117  prepointstatus == fUndefined ||
118  (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary))
119  {
120  unsigned int layer_id = 99; // no layer number for the hit, use a non-existent one for now, replace it later
121  // this is for debugging weird occurances we have occasionally
122  if (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary)
123  {
124  std::cout << GetName() << ": New Hit for " << std::endl;
125  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
126  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
127  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
128  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
129  std::cout << "last track: " << m_SaveTrackId
130  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
131  std::cout << "phys pre vol: " << volume->GetName()
132  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
133  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
134  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
135  }
136  // if previous hit was saved, hit pointer was set to nullptr
137  // and we have to make a new one
138  if (!m_Hit)
139  {
140  m_Hit = new PHG4Hitv1();
141  }
142  m_Hit->set_layer(layer_id);
143  //here we set the entrance values in cm
144  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
145  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
146  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
147 
148  // momentum
149  m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
150  m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
151  m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
152 
153  // time in ns
154  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
155  //set and save the track ID
156  m_Hit->set_trkid(aTrack->GetTrackID());
157  m_SaveTrackId = aTrack->GetTrackID();
158  //set the initial energy deposit
159  m_Hit->set_edep(0);
160  if (whichactive > 0) // return of IsInTpcDetector, > 0 hit in tpc gas volume, < 0 hit in support structures
161  {
162  m_Hit->set_eion(0);
163  // Now save the container we want to add this hit to
165  }
166  else
167  {
169  }
170  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
171  {
172  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
173  {
174  m_Hit->set_trkid(pp->GetUserTrackId());
175  m_Hit->set_shower_id(pp->GetShower()->get_id());
176  m_Shower = pp->GetShower();
177  }
178  }
179  // some sanity checks for inconsistencies
180  // check if this hit was created, if not print out last post step status
181  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
182  {
183  std::cout << GetName() << ": hit was not created" << std::endl;
184  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
185  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
186  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
187  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
188  std::cout << "last track: " << m_SaveTrackId
189  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
190  std::cout << "phys pre vol: " << volume->GetName()
191  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
192  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
193  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
194  exit(1);
195  }
196  // check if track id matches the initial one when the hit was created
197  if (aTrack->GetTrackID() != m_SaveTrackId)
198  {
199  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
200  std::cout << "saved track: " << m_SaveTrackId
201  << ", current trackid: " << aTrack->GetTrackID()
202  << ", prestep status: " << prePoint->GetStepStatus()
203  << ", previous post step status: " << m_SavePostStepStatus
204  << std::endl;
205 
206  exit(1);
207  }
208  m_SavePreStepStatus = prePoint->GetStepStatus();
209  m_SavePostStepStatus = postPoint->GetStepStatus();
210  m_SaveVolPre = volume;
211  m_SaveVolPost = touchpost->GetVolume();
212  // here we just update the exit values, it will be overwritten
213  // for every step until we leave the volume or the particle
214  // ceases to exist
215  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
216  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
217  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
218 
219  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
220  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
221  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
222 
223  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
224 
225  //sum up the energy to get total deposited
226  m_Hit->set_edep(m_Hit->get_edep() + edep);
227  if (whichactive > 0)
228  {
229  m_Hit->set_eion(m_Hit->get_eion() + eion);
230  }
231  if (geantino)
232  {
233  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
234  if (whichactive > 0)
235  {
236  m_Hit->set_eion(-1);
237  }
238  }
239  if (edep > 0)
240  {
241  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
242  {
243  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
244  {
245  pp->SetKeep(1); // we want to keep the track
246  }
247  }
248  }
249 
250  // if any of these conditions is true this is the last step in
251  // this volume and we need to save the hit
252  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
253  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
254  // (happens when your detector goes outside world volume)
255  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
256  // aTrack->GetTrackStatus() == fStopAndKill is also set)
257  // aTrack->GetTrackStatus() == fStopAndKill: track ends
258  if ((m_UseG4StepsFlag > 0 && whichactive > 0) ||
259  postPoint->GetStepStatus() == fGeomBoundary ||
260  postPoint->GetStepStatus() == fWorldBoundary ||
261  postPoint->GetStepStatus() == fAtRestDoItProc ||
262  aTrack->GetTrackStatus() == fStopAndKill)
263  {
264  // save only hits with energy deposit (or -1 for geantino)
265  if (m_Hit->get_edep())
266  {
267  m_CurrentHitContainer->AddHit(layer_id, m_Hit);
268  if (m_Shower)
269  {
271  }
272 
273  double rin = sqrt(m_Hit->get_x(0) * m_Hit->get_x(0) + m_Hit->get_y(0) * m_Hit->get_y(0));
274  double rout = sqrt(m_Hit->get_x(1) * m_Hit->get_x(1) + m_Hit->get_y(1) * m_Hit->get_y(1));
275  if (Verbosity() > 10)
276  if ((rin > 69.0 && rin < 70.125) || (rout > 69.0 && rout < 70.125))
277  {
278  std::cout << "Added Tpc g4hit with rin, rout = " << rin << " " << rout
279  << " g4hitid " << m_Hit->get_hit_id() << std::endl;
280  std::cout << " xin " << m_Hit->get_x(0)
281  << " yin " << m_Hit->get_y(0)
282  << " zin " << m_Hit->get_z(0)
283  << " rin " << rin
284  << std::endl;
285  std::cout << " xout " << m_Hit->get_x(1)
286  << " yout " << m_Hit->get_y(1)
287  << " zout " << m_Hit->get_z(1)
288  << " rout " << rout
289  << std::endl;
290  std::cout << " xav " << (m_Hit->get_x(1) + m_Hit->get_x(0)) / 2.0
291  << " yav " << (m_Hit->get_y(1) + m_Hit->get_y(0)) / 2.0
292  << " zav " << (m_Hit->get_z(1) + m_Hit->get_z(0)) / 2.0
293  << " rav " << (rout + rin) / 2.0
294  << std::endl;
295  }
296 
297  // ownership has been transferred to container, set to null
298  // so we will create a new hit for the next track
299  m_Hit = nullptr;
300  }
301  else
302  {
303  // if this hit has no energy deposit, just reset it for reuse
304  // this means we have to delete it in the dtor. If this was
305  // the last hit we processed the memory is still allocated
306  m_Hit->Reset();
307  }
308  }
309  // return true to indicate the hit was used
310  return true;
311  }
312  else
313  {
314  return false;
315  }
316 }
317 
318 //____________________________________________________________________________..
320 {
321  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
322  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
323 
324  // if we do not find the node it's messed up.
325  if (!m_HitContainer)
326  {
327  std::cout << "PHG4TpcSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
328  }
329  if (!m_AbsorberHitContainer)
330  {
331  if (Verbosity() > 1)
332  {
333  std::cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
334  }
335  }
336 }
337 
339 {
340  if (type == "G4HIT")
341  {
343  return;
344  }
345  else if (type == "G4HIT_ABSORBER")
346  {
348  return;
349  }
350  std::cout << "Invalid output hit node type " << type << std::endl;
351  gSystem->Exit(1);
352  return;
353 }