Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BeamLineMagnetSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BeamLineMagnetSteppingAction.cc
2 
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 <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, fAtRest...
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 <cmath> // for isfinite
36 #include <iostream>
37 #include <string> // for operator<<, string
38 
39 class PHCompositeNode;
40 
41 //____________________________________________________________________________..
43  : PHG4SteppingAction(detector->GetName())
44  , m_Detector(detector)
45  , m_Params(parameters)
46  , m_ActiveFlag(m_Params->get_int_param("active"))
47  , m_AbsorberActiveFlag(m_Params->get_int_param("absorberactive"))
48  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
49 {
50 }
51 
53 {
54  // if the last hit was a zero energie deposit hit, it is just reset
55  // and the memory is still allocated, so we need to delete it here
56  // if the last hit was saved, hit is a nullptr pointer which are
57  // legal to delete (it results in a no operation)
58  delete m_Hit;
59 }
60 
61 //____________________________________________________________________________..
62 bool BeamLineMagnetSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
63 {
64  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
65  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
66  // get volume of the current step
67  G4VPhysicalVolume* volume = touch->GetVolume();
68  // IsInVTX(volume) returns
69  // == 0 outside of pstof
70  // > 0 for hits in active volume
71  // < 0 for hits in passive material
72  int whichactive = m_Detector->IsInBeamLineMagnet(volume);
73  if (!whichactive)
74  {
75  return false;
76  }
77 
78  // collect energy and track length step by step
79  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
80  const G4Track* aTrack = aStep->GetTrack();
81 
82  // if this block stops everything ONLY if in the magnet iron, just put all kinetic energy into edep
83  if (whichactive == -1)
84  {
85  if (m_BlackHoleFlag)
86  {
87  edep = aTrack->GetKineticEnergy() / GeV;
88  G4Track* killtrack = const_cast<G4Track*>(aTrack);
89  killtrack->SetTrackStatus(fStopAndKill);
90  }
92  {
93  return false;
94  }
95  }
96  if (whichactive == -2)
97  {
98  edep = aTrack->GetKineticEnergy() / GeV;
99  G4Track* killtrack = const_cast<G4Track*>(aTrack);
100  killtrack->SetTrackStatus(fStopAndKill);
101  }
102  if (!m_ActiveFlag)
103  {
104  return false;
105  }
106 
107  int magnet_id = volume->GetCopyNo(); // magnet id is stored in copy number
108  bool geantino = false;
109  // the check for the pdg code speeds things up, I do not want to make
110  // an expensive string compare for every track when we know
111  // geantino or chargedgeantino has pid=0
112  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
113  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
114  {
115  geantino = true;
116  }
117  G4StepPoint* prePoint = aStep->GetPreStepPoint();
118  G4StepPoint* postPoint = aStep->GetPostStepPoint();
119  // std::cout << "track id " << aTrack->GetTrackID() << std::endl;
120  // std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
121  // std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
122 
123  switch (prePoint->GetStepStatus())
124  {
125  case fPostStepDoItProc:
126  if (m_SavePostStepStatus != fGeomBoundary)
127  {
128  break;
129  }
130  else
131  {
132  std::cout << GetName() << ": New Hit for " << std::endl;
133  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
134  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
135  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
136  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
137  std::cout << "last track: " << m_SaveTrackId
138  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
139  std::cout << "phys pre vol: " << volume->GetName()
140  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
141  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
142  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
143  }
144  [[fallthrough]];
145  case fGeomBoundary:
146  case fUndefined:
147  if (!m_Hit)
148  {
149  m_Hit = new PHG4Hitv1();
150  }
151  m_Hit->set_layer(magnet_id);
152  //here we set the entrance values in cm
153  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
154  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
155  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
156  // time in ns
157  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
158  //set the track ID
159  m_Hit->set_trkid(aTrack->GetTrackID());
160  m_SaveTrackId = aTrack->GetTrackID();
161  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
162  {
163  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
164  {
165  m_Hit->set_trkid(pp->GetUserTrackId());
166  }
167  }
168  //set the initial energy deposit
169  m_EdepSum = 0;
170  if (whichactive > 0)
171  {
173  }
174  else
175  {
177  }
178  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
179  {
180  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
181  {
182  m_Hit->set_trkid(pp->GetUserTrackId());
183  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
184  }
185  }
186 
187  break;
188  default:
189  break;
190  }
191 
192  // some sanity checks for inconsistencies
193  // check if this hit was created, if not print out last post step status
194  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
195  {
196  std::cout << GetName() << ": hit was not created" << std::endl;
197  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
198  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
199  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
200  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
201  std::cout << "last track: " << m_SaveTrackId
202  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
203  std::cout << "phys pre vol: " << volume->GetName()
204  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
205  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
206  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
207  gSystem->Exit(1);
208  }
209  // check if track id matches the initial one when the hit was created
210  if (aTrack->GetTrackID() != m_SaveTrackId)
211  {
212  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
213  std::cout << "saved track: " << m_SaveTrackId
214  << ", current trackid: " << aTrack->GetTrackID()
215  << ", prestep status: " << prePoint->GetStepStatus()
216  << ", previous post step status: " << m_SavePostStepStatus
217  << std::endl;
218 
219  gSystem->Exit(1);
220  }
221  m_SavePreStepStatus = prePoint->GetStepStatus();
222  m_SavePostStepStatus = postPoint->GetStepStatus();
223  m_SaveVolPre = volume;
224  m_SaveVolPost = touchpost->GetVolume();
225 
226  // here we just update the exit values, it will be overwritten
227  // for every step until we leave the volume or the particle
228  // ceases to exist
229  //sum up the energy to get total deposited
230  m_EdepSum += edep;
231  // if any of these conditions is true this is the last step in
232  // this volume and we need to save the hit
233  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
234  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
235  // (happens when your detector goes outside world volume)
236  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
237  // aTrack->GetTrackStatus() == fStopAndKill is also set)
238  // aTrack->GetTrackStatus() == fStopAndKill: track ends
239  if (postPoint->GetStepStatus() == fGeomBoundary ||
240  postPoint->GetStepStatus() == fWorldBoundary ||
241  postPoint->GetStepStatus() == fAtRestDoItProc ||
242  aTrack->GetTrackStatus() == fStopAndKill)
243  {
244  // save only hits with energy deposit (or geantino)
245  if (m_EdepSum > 0 || geantino)
246  {
247  // update values at exit coordinates and set keep flag
248  // of track to keep
249  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
250  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
251  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
252 
253  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
254  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
255  {
256  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
257  {
258  pp->SetKeep(1); // we want to keep the track
259  }
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  {
268  }
269  m_SaveHitContainer->AddHit(magnet_id, m_Hit);
270  // ownership has been transferred to container, set to null
271  // so we will create a new hit for the next track
272  m_Hit = nullptr;
273  }
274  else
275  {
276  // if this hit has no energy deposit, just reset it for reuse
277  // this means we have to delete it in the dtor. If this was
278  // the last hit we processed the memory is still allocated
279  m_Hit->Reset();
280  }
281  }
282  // return true to indicate the hit was used
283  return true;
284 }
285 
286 //____________________________________________________________________________..
288 {
289  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
290  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
291 
292  // If this is a black hole but otherwise not an active detector there is
293  // no hitnode, otherwise this will crash in the user stepping action
294  if (!m_HitContainer)
295  {
296  if (!m_BlackHoleFlag)
297  {
298  std::cout << "BeamLineMagnetSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
299  }
300  }
301  if (!m_AbsorberHitContainer)
302  {
303  if (Verbosity() > 1)
304  {
305  std::cout << "BeamLineMagnetSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
306  }
307  }
308 }
309 
311 {
312  if (type == "G4HIT")
313  {
315  return;
316  }
317  else if (type == "G4HIT_ABSORBER")
318  {
320  return;
321  }
322  std::cout << "Invalid output hit node type " << type << std::endl;
323  gSystem->Exit(1);
324  return;
325 }