ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4BlockSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4BlockSteppingAction.cc
2 
3 #include "PHG4BlockDetector.h"
4 #include "PHG4StepStatusDecode.h"
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
14 
15 #include <phool/getClass.h>
16 
17 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
18 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
19 #include <Geant4/G4Step.hh>
20 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
21 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fPostSt...
22 #include <Geant4/G4String.hh> // for G4String
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 <cmath> // for isfinite
34 #include <cstdlib> // for exit
35 #include <iostream>
36 #include <string> // for operator<<, string
37 
38 class PHCompositeNode;
39 
40 using namespace std;
41 //____________________________________________________________________________..
43  : PHG4SteppingAction(detector->GetName())
44  , m_Detector(detector)
45  , m_Params(parameters)
46  , m_HitContainer(nullptr)
47  , m_Hit(nullptr)
48  , m_SaveShower(nullptr)
49  , m_SaveVolPre(nullptr)
50  , m_SaveVolPost(nullptr)
51  , m_SaveTrackId(-1)
52  , m_SavePreStepStatus(-1)
53  , m_SavePostStepStatus(-1)
54  , m_ActiveFlag(m_Params->get_int_param("active"))
55  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
56  , m_UseG4StepsFlag(m_Params->get_int_param("use_g4steps"))
57 {
58 }
59 
61 {
62  // if the last hit was a zero energie deposit hit, it is just reset
63  // and the memory is still allocated, so we need to delete it here
64  // if the last hit was saved, hit is a nullptr pointer which are
65  // legal to delete (it results in a no operation)
66  delete m_Hit;
67 }
68 
69 //____________________________________________________________________________..
71 {
74  // get volume of the current step
76 
77  if (!m_Detector->IsInBlock(volume))
78  {
79  return false;
80  }
81 
82  // collect energy and track length step by step
84  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
85  const G4Track* aTrack = aStep->GetTrack();
86 
87  // if this block stops everything, just put all kinetic energy into edep
88  if (m_BlackHoleFlag)
89  {
90  edep = aTrack->GetKineticEnergy() / GeV;
91  G4Track* killtrack = const_cast<G4Track*>(aTrack);
92  killtrack->SetTrackStatus(fStopAndKill);
93  }
94 
95  // make sure we are in a volume
96  if (m_ActiveFlag)
97  {
98  int layer_id = m_Detector->get_Layer();
99  bool geantino = false;
100  // the check for the pdg code speeds things up, I do not want to make
101  // an expensive string compare for every track when we know
102  // geantino or chargedgeantino has pid=0
103  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
104  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
105  {
106  geantino = true;
107  }
108  G4StepPoint* prePoint = aStep->GetPreStepPoint();
109  G4StepPoint* postPoint = aStep->GetPostStepPoint();
110  // cout << "track id " << aTrack->GetTrackID() << endl;
111  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
112  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
113  int prepointstatus = prePoint->GetStepStatus();
114  if (prepointstatus == fGeomBoundary ||
115  prepointstatus == fUndefined ||
116  (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary) ||
117  m_UseG4StepsFlag > 0)
118  {
119  if (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary)
120  {
121  cout << GetName() << ": New Hit for " << endl;
122  cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
123  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
124  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
125  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
126  cout << "last track: " << m_SaveTrackId
127  << ", current trackid: " << aTrack->GetTrackID() << endl;
128  cout << "phys pre vol: " << volume->GetName()
129  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
130  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
131  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
132  }
133  if (!m_Hit)
134  {
135  m_Hit = new PHG4Hitv1();
136  }
137  m_Hit->set_layer(layer_id);
138  //here we set the entrance values in cm
139  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
140  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
141  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
142  // time in ns
143  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
144  //set the track ID
145  m_Hit->set_trkid(aTrack->GetTrackID());
146  m_SaveTrackId = aTrack->GetTrackID();
147 
148  //set the initial energy deposit
149  m_Hit->set_edep(0);
150  if (!geantino && !m_BlackHoleFlag && m_ActiveFlag)
151  {
152  m_Hit->set_eion(0);
153  }
155  {
156  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
157  {
158  m_Hit->set_trkid(pp->GetUserTrackId());
159  m_Hit->set_shower_id(pp->GetShower()->get_id());
160  m_SaveShower = pp->GetShower();
161  }
162  }
163  }
164 
165  // some sanity checks for inconsistencies
166  // check if this hit was created, if not print out last post step status
167  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
168  {
169  cout << GetName() << ": hit was not created" << endl;
170  cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
171  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
172  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
173  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
174  cout << "last track: " << m_SaveTrackId
175  << ", current trackid: " << aTrack->GetTrackID() << endl;
176  cout << "phys pre vol: " << volume->GetName()
177  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
178  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
179  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
180  exit(1);
181  }
182  m_SavePostStepStatus = postPoint->GetStepStatus();
183  // check if track id matches the initial one when the hit was created
184  if (aTrack->GetTrackID() != m_SaveTrackId)
185  {
186  cout << "hits do not belong to the same track" << endl;
187  cout << "saved track: " << m_SaveTrackId
188  << ", current trackid: " << aTrack->GetTrackID()
189  << endl;
190  exit(1);
191  }
192  m_SavePreStepStatus = prePoint->GetStepStatus();
193  m_SavePostStepStatus = postPoint->GetStepStatus();
195  m_SaveVolPost = touchpost->GetVolume();
196 
197  // here we just update the exit values, it will be overwritten
198  // for every step until we leave the volume or the particle
199  // ceases to exist
200  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
201  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
202  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
203 
204  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
205  //sum up the energy to get total deposited
207  // update ionization energy only for active volumes, not for black holes or geantinos
208  // if the hit is created without eion, get_eion() returns NAN
209  // if you see NANs check the creation of the hit
210  if (m_ActiveFlag && !m_BlackHoleFlag && !geantino)
211  {
212  m_Hit->set_eion(m_Hit->get_eion() + eion);
213  }
214  if (geantino)
215  {
216  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
217  }
218  if (edep > 0)
219  {
221  {
222  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
223  {
224  pp->SetKeep(1); // we want to keep the track
225  }
226  }
227  }
228  // if any of these conditions is true this is the last step in
229  // this volume and we need to save the hit
230  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
231  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
232  // (happens when your detector goes outside world volume)
233  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
234  // aTrack->GetTrackStatus() == fStopAndKill is also set)
235  // aTrack->GetTrackStatus() == fStopAndKill: track ends
236  if (postPoint->GetStepStatus() == fGeomBoundary ||
237  postPoint->GetStepStatus() == fWorldBoundary ||
238  postPoint->GetStepStatus() == fAtRestDoItProc ||
239  aTrack->GetTrackStatus() == fStopAndKill ||
240  m_UseG4StepsFlag > 0)
241  {
242  // save only hits with energy deposit (or -1 for geantino)
243  if (m_Hit->get_edep())
244  {
245  m_HitContainer->AddHit(layer_id, m_Hit);
246  if (m_SaveShower)
247  {
249  }
250  // ownership has been transferred to container, set to null
251  // so we will create a new hit for the next track
252  m_Hit = nullptr;
253  }
254  else
255  {
256  // if this hit has no energy deposit, just reset it for reuse
257  // this means we have to delete it in the dtor. If this was
258  // the last hit we processed the memory is still allocated
259  m_Hit->Reset();
260  }
261  }
262  // return true to indicate the hit was used
263  return true;
264  }
265  else
266  {
267  return false;
268  }
269 }
270 
271 //____________________________________________________________________________..
273 {
274  string hitnodename;
275  if (m_Detector->SuperDetector() != "NONE")
276  {
277  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
278  }
279  else
280  {
281  hitnodename = "G4HIT_" + m_Detector->GetName();
282  }
283 
284  //now look for the map and grab a pointer to it.
285  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
286 
287  // if we do not find the node we need to make it.
288  if (!m_HitContainer)
289  {
290  std::cout << "PHG4BlockSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
291  }
292 }