ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4JLeicBeamLineMagnetSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4JLeicBeamLineMagnetSteppingAction.cc
3 
5 
6 #include <g4main/PHG4Hit.h>
8 #include <g4main/PHG4Hitv1.h>
9 #include <g4main/PHG4Shower.h>
10 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
12 
13 #include <phool/getClass.h>
14 
15 #include <TSystem.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, fAtRest...
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 <iostream>
35 #include <string> // for operator<<, string
36 
37 class PHCompositeNode;
38 
39 using namespace std;
40 //____________________________________________________________________________..
42  : PHG4SteppingAction(detector->GetName())
43  , m_Detector(detector)
44  , m_HitContainer(nullptr)
45  , m_AbsorberHitContainer(nullptr)
46  , m_Hit(nullptr)
47  , m_SaveHitContainer(nullptr)
48  , m_SaveVolPre(nullptr)
49  , m_SaveVolPost(nullptr)
50  , m_SaveTrackId(-1)
51  , m_SavePreStepStatus(-1)
52  , m_SavePostStepStatus(-1)
53  , m_EdepSum(0)
54  , m_EionSum(0)
55 {
56 }
57 
59 {
60  // if the last hit was a zero energie deposit hit, it is just reset
61  // and the memory is still allocated, so we need to delete it here
62  // if the last hit was saved, hit is a nullptr pointer which are
63  // legal to delete (it results in a no operation)
64  delete m_Hit;
65 }
66 
67 //____________________________________________________________________________..
69 {
72  // get volume of the current step
74  // IsInVTX(volume) returns
75  // == 0 outside of pstof
76  // > 0 for hits in active volume
77  // < 0 for hits in passive material
78  int whichactive = m_Detector->IsInBeamLineMagnet(volume);
79  if (!whichactive)
80  {
81  return false;
82  }
83 
84  // collect energy and track length step by step
86  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
87  const G4Track* aTrack = aStep->GetTrack();
88 
89  int layer_id = whichactive - 1; // layer id is +1 in detector description
90  bool geantino = false;
91  // the check for the pdg code speeds things up, I do not want to make
92  // an expensive string compare for every track when we know
93  // geantino or chargedgeantino has pid=0
94  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
95  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
96  {
97  geantino = true;
98  }
99  G4StepPoint* prePoint = aStep->GetPreStepPoint();
100  G4StepPoint* postPoint = aStep->GetPostStepPoint();
101  // cout << "track id " << aTrack->GetTrackID() << endl;
102  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
103  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
104 
105  switch (prePoint->GetStepStatus())
106  {
107  case fPostStepDoItProc:
109  {
110  break;
111  }
112  else
113  {
114  cout << GetName() << ": New Hit for " << endl;
115  cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
116  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
117  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
118  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
119  cout << "last track: " << m_SaveTrackId
120  << ", current trackid: " << aTrack->GetTrackID() << endl;
121  cout << "phys pre vol: " << volume->GetName()
122  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
123  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
124  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
125  }
126  case fGeomBoundary:
127  case fUndefined:
128  if (!m_Hit)
129  {
130  m_Hit = new PHG4Hitv1();
131  }
132  m_Hit->set_layer(layer_id);
133  //here we set the entrance values in cm
134  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
135  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
136  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
137  // time in ns
138  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
139  //set the track ID
140  m_Hit->set_trkid(aTrack->GetTrackID());
141  m_SaveTrackId = aTrack->GetTrackID();
143  {
144  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
145  {
146  m_Hit->set_trkid(pp->GetUserTrackId());
147  }
148  }
149  //set the initial energy deposit
150  m_EdepSum = 0;
151  if (whichactive > 0)
152  {
153  m_EionSum = 0;
154  m_Hit->set_eion(0);
156  }
157  else
158  {
160  }
162  {
163  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
164  {
165  m_Hit->set_trkid(pp->GetUserTrackId());
166  pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
167  }
168  }
169 
170  break;
171  default:
172  break;
173  }
174 
175  // some sanity checks for inconsistencies
176  // check if this hit was created, if not print out last post step status
177  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
178  {
179  cout << GetName() << ": hit was not created" << endl;
180  cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
181  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
182  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
183  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
184  cout << "last track: " << m_SaveTrackId
185  << ", current trackid: " << aTrack->GetTrackID() << endl;
186  cout << "phys pre vol: " << volume->GetName()
187  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
188  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
189  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
190  gSystem->Exit(1);
191  }
192  // check if track id matches the initial one when the hit was created
193  if (aTrack->GetTrackID() != m_SaveTrackId)
194  {
195  cout << GetName() << ": hits do not belong to the same track" << endl;
196  cout << "saved track: " << m_SaveTrackId
197  << ", current trackid: " << aTrack->GetTrackID()
198  << ", prestep status: " << prePoint->GetStepStatus()
199  << ", previous post step status: " << m_SavePostStepStatus
200  << endl;
201 
202  gSystem->Exit(1);
203  }
204  m_SavePreStepStatus = prePoint->GetStepStatus();
205  m_SavePostStepStatus = postPoint->GetStepStatus();
207  m_SaveVolPost = touchpost->GetVolume();
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  //sum up the energy to get total deposited
213  m_EdepSum += edep;
214  if (whichactive > 0)
215  {
216  m_EionSum += eion;
217  }
218  // if any of these conditions is true this is the last step in
219  // this volume and we need to save the hit
220  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
221  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
222  // (happens when your detector goes outside world volume)
223  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
224  // aTrack->GetTrackStatus() == fStopAndKill is also set)
225  // aTrack->GetTrackStatus() == fStopAndKill: track ends
226  if (postPoint->GetStepStatus() == fGeomBoundary ||
227  postPoint->GetStepStatus() == fWorldBoundary ||
228  postPoint->GetStepStatus() == fAtRestDoItProc ||
229  aTrack->GetTrackStatus() == fStopAndKill)
230  {
231  // save only hits with energy deposit (or geantino)
232  if (m_EdepSum > 0 || geantino)
233  {
234  // update values at exit coordinates and set keep flag
235  // of track to keep
236  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
237  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
238  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
239 
240  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
242  {
243  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
244  {
245  pp->SetKeep(1); // we want to keep the track
246  }
247  }
248  if (geantino)
249  {
250  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
251  if (whichactive > 0)
252  {
253  m_Hit->set_eion(-1);
254  }
255  }
256  else
257  {
259  }
260  if (whichactive > 0)
261  {
263  }
264  m_SaveHitContainer->AddHit(layer_id, m_Hit);
265  // ownership has been transferred to container, set to null
266  // so we will create a new hit for the next track
267  m_Hit = nullptr;
268  }
269  else
270  {
271  // if this hit has no energy deposit, just reset it for reuse
272  // this means we have to delete it in the dtor. If this was
273  // the last hit we processed the memory is still allocated
274  m_Hit->Reset();
275  }
276  }
277  // return true to indicate the hit was used
278  return true;
279 }
280 
281 //____________________________________________________________________________..
283 {
284  string hitnodename;
285  string absorbernodename;
286  if (m_Detector->SuperDetector() != "NONE")
287  {
288  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
289  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
290  }
291  else
292  {
293  hitnodename = "G4HIT_" + m_Detector->GetName();
294  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
295  }
296 
297  //now look for the map and grab a pointer to it.
298  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
299  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
300 
301  // if we do not find the node we need to make it.
302  if (!m_HitContainer)
303  {
304  std::cout << "G4JLeicBeamLineMagnetSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
305  }
306  if (!m_AbsorberHitContainer)
307  {
308  if (Verbosity() > 1)
309  {
310  cout << "G4JLeicBeamLineMagnetSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
311  }
312  }
313 }