ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4EPDSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4EPDSteppingAction.cc
1 /* vim: set sw=2 ft=cpp: */
2 
4 
5 #include "PHG4EPDDetector.h"
6 
8 
9 #include <phool/getClass.h>
10 
11 #include <g4main/PHG4Hit.h>
13 #include <g4main/PHG4Hitv1.h>
14 #include <g4main/PHG4Shower.h>
17 
18 #include <TSystem.h>
19 
20 #include <Geant4/G4ParticleDefinition.hh>
21 #include <Geant4/G4ReferenceCountedHandle.hh>
22 #include <Geant4/G4Step.hh>
23 #include <Geant4/G4StepPoint.hh>
24 #include <Geant4/G4StepStatus.hh>
25 #include <Geant4/G4SystemOfUnits.hh>
26 #include <Geant4/G4ThreeVector.hh>
27 #include <Geant4/G4TouchableHandle.hh>
28 #include <Geant4/G4Track.hh>
29 #include <Geant4/G4TrackStatus.hh>
30 #include <Geant4/G4Types.hh>
31 #include <Geant4/G4VTouchable.hh>
32 #include <Geant4/G4VUserTrackInformation.hh>
33 
34 #include <cstdint> // for int32_t
35 #include <iostream>
36 #include <string>
37 
38 class G4VPhysicalVolume;
39 
41  const PHParameters*)
42  : PHG4SteppingAction(detector->GetName())
43  , m_Detector(detector)
44 {
45 }
46 
48 {
49  delete m_Hit;
50 }
51 
53 {
54  G4StepPoint* prestep = aStep->GetPreStepPoint();
55  G4TouchableHandle prehandle = prestep->GetTouchableHandle();
56 
57  G4VPhysicalVolume* volume = prehandle->GetVolume();
58 
59  // m_Detector->IsInDetector(volume)
60  // returns
61  // 0 is outside of EPD
62  // 1 is inside scintillator
63  // -1 is inside support structure (dead material)
64  int whichactive = m_Detector->IsInDetector(volume);
65  if (!whichactive)
66  {
67  return false;
68  }
69 
70  G4double deposit = aStep->GetTotalEnergyDeposit() / GeV;
71  G4double ionising = deposit - aStep->GetNonIonizingEnergyDeposit() / GeV;
72  G4double light_yield = GetVisibleEnergyDeposition(aStep) / GeV;
73 
74  // G4StepStatus prestatus = prestep->GetStepStatus();
75 
76  int32_t tile_id = m_Detector->module_id_for(volume);
77 
78  G4Track const* track = aStep->GetTrack();
79 
81 
82  bool geantino = (particle->GetPDGEncoding() == 0 && particle->GetParticleName().find("geantino") != std::string::npos);
83 
84  G4StepPoint* prePoint = aStep->GetPreStepPoint();
85  G4StepPoint* postPoint = aStep->GetPostStepPoint();
86 
87  switch (prePoint->GetStepStatus())
88  {
89  case fPostStepDoItProc:
91  {
92  break;
93  }
94  else
95  {
96  std::cout << GetName() << ": New Hit for " << std::endl;
97  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
98  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
99  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
100  }
101  [[fallthrough]];
102  case fGeomBoundary:
103  case fUndefined:
104  if (m_Hit == nullptr)
105  {
106  m_Hit = new PHG4Hitv1();
107  }
108 
109  // only for active columes (scintillators)
110  if (whichactive > 0)
111  {
112  m_Hit->set_scint_id(tile_id);
113  m_Hit->set_eion(0);
115  }
116  m_Hit->set_x(0, prestep->GetPosition().x() / cm);
117  m_Hit->set_y(0, prestep->GetPosition().y() / cm);
118  m_Hit->set_z(0, prestep->GetPosition().z() / cm);
119  m_Hit->set_t(0, prestep->GetGlobalTime() / nanosecond);
120 
121  m_Hit->set_trkid(track->GetTrackID());
122 
123  if (PHG4TrackUserInfoV1* userinfo = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation()))
124  {
125  m_Hit->set_trkid(userinfo->GetUserTrackId());
126 
127  userinfo->GetShower()->add_g4hit_id(m_HitContainer->GetID(), m_Hit->get_hit_id());
128  }
129 
130  m_Hit->set_edep(0);
131  break;
132  default:
133  break;
134  }
135 
136  G4StepPoint* poststep = aStep->GetPostStepPoint();
137  const G4ThreeVector postpos = poststep->GetPosition();
138  m_SavePostStepStatus = postPoint->GetStepStatus();
139  m_Hit->set_edep(m_Hit->get_edep() + deposit);
140  if (whichactive > 0)
141  {
142  m_Hit->set_eion(m_Hit->get_eion() + ionising);
143  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield * GetLightCorrection(postpos.x(), postpos.y()));
144  }
145 
146  if (postPoint->GetStepStatus() != fGeomBoundary &&
147  postPoint->GetStepStatus() != fWorldBoundary &&
148  postPoint->GetStepStatus() != fAtRestDoItProc &&
149  track->GetTrackStatus() != fStopAndKill)
150  {
151  return true;
152  }
153  if (m_Hit->get_edep() <= 0 && !geantino)
154  {
155  m_Hit->Reset();
156 
157  return true;
158  }
159 
160  m_Hit->set_x(1, poststep->GetPosition().x() / cm);
161  m_Hit->set_y(1, poststep->GetPosition().y() / cm);
162  m_Hit->set_z(1, poststep->GetPosition().z() / cm);
163  m_Hit->set_t(1, poststep->GetGlobalTime() / nanosecond);
164 
165  PHG4TrackUserInfoV1* userinfo = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation());
166 
167  if (userinfo != nullptr)
168  {
169  userinfo->SetKeep(1);
170  }
171 
172  if (geantino)
173  {
174  m_Hit->set_edep(-1.);
175  m_Hit->set_eion(-1.);
176  m_Hit->set_light_yield(-1.);
177  }
178 
179  m_HitContainer->AddHit(tile_id, m_Hit);
180 
181  m_Hit = nullptr;
182 
183  return true;
184 }
185 
187 {
188  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
189 
190  m_SupportHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_SupportNodeName);
191  // if we do not find the node it's messed up.
192  if (!m_HitContainer)
193  {
194  std::cout << "PHG4EPDSteppingAction::SetTopNode - unable to find hit node " << m_HitNodeName << std::endl;
195  gSystem->Exit(1);
196  }
197  // this is perfectly fine if support hits are disabled
198  if (!m_SupportHitContainer)
199  {
200  if (Verbosity() > 0)
201  {
202  std::cout << "PHG4EPDSteppingAction::SetTopNode - unable to find support hit node " << m_SupportNodeName << std::endl;
203  }
204  }
205 }
206 
207 void PHG4EPDSteppingAction::SetHitNodeName(const std::string& type, const std::string& name)
208 {
209  if (type == "G4HIT")
210  {
212  return;
213  }
214  else if (type == "G4HIT_SUPPORT")
215  {
217  return;
218  }
219  std::cout << "Invalid output hit node type " << type << std::endl;
220  gSystem->Exit(1);
221  return;
222 }