Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4SteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4SteppingAction.cc
1 
9 #include "PHG4SteppingAction.h"
10 
11 #include "PHG4Hit.h"
12 
13 #include <Geant4/G4AffineTransform.hh> // for G4AffineTransform
14 #include <Geant4/G4EmSaturation.hh>
15 #include <Geant4/G4LossTableManager.hh>
16 #include <Geant4/G4Material.hh>
17 #include <Geant4/G4MaterialPropertiesTable.hh> // for G4MaterialProperties...
18 #include <Geant4/G4NavigationHistory.hh>
19 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHa...
20 #include <Geant4/G4Step.hh>
21 #include <Geant4/G4StepPoint.hh>
22 #include <Geant4/G4String.hh> // for G4String
23 #include <Geant4/G4SystemOfUnits.hh>
24 #include <Geant4/G4ThreeVector.hh>
25 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
26 #include <Geant4/G4Track.hh>
27 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
28 
29 #include <algorithm>
30 #include <cassert>
31 #include <cmath> // for isfinite, NAN, sqrt
32 #include <iostream>
33 
35  : m_Verbosity(i)
36  , m_LightBalanceInnerRadius(NAN)
37  , m_LightBalanceInnerCorr(NAN)
38  , m_LightBalanceOuterRadius(NAN)
39  , m_LightBalanceOuterCorr(NAN)
40  , m_Name(name)
41 {
42 }
43 
44 double
46 {
47  // pirated from G4Scintillation::PostStepDoIt()
48 
49  double light_yield = 0;
50 
51  const G4Track* aTrack = step->GetTrack();
52  const G4Material* aMaterial = aTrack->GetMaterial();
53  G4MaterialPropertiesTable* aMaterialPropertiesTable =
54  aMaterial->GetMaterialPropertiesTable();
55  if (!aMaterialPropertiesTable)
56  {
57  const std::string &mname(aMaterial->GetName());
58 
59  std::set<std::string>::const_iterator it =
61 
63  {
65 
66  std::cout << "PHG4SteppingAction::GetScintLightYield - WARNING - "
67  << "can not find Material Properties Table for material " << mname
68  << ", will assume it do NOT scintillate. "
69  << "Please ignore this warning if you do not expect scintillation light from "
70  << mname << std::endl;
71  }
72 
73  return 0.;
74  }
75 
76  if (aMaterialPropertiesTable->ConstPropertyExists("SCINTILLATIONYIELD"))
77  {
78  light_yield = aMaterialPropertiesTable->GetConstProperty("SCINTILLATIONYIELD") * GetVisibleEnergyDeposition(step) * GeV;
79 
80  return light_yield;
81  } // if (aMaterialPropertiesTable->ConstPropertyExists("SCINTILLATIONYIELD"))
82  else
83  {
84  const std::string &mname(aMaterial->GetName());
85 
86  std::set<std::string>::const_iterator it =
88 
90  {
92 
93  std::cout << "PHG4SteppingAction::GetScintLightYield - WARNING - "
94  << "can not find scintillation light yield for material " << mname
95  << ", will assume it do NOT scintillate. "
96  << "Please ignore this warning if you do not expect scintillation light from "
97  << mname << std::endl;
98  }
99 
100  return 0.;
101  }
102 
103  return light_yield;
104 }
105 
106 double
108 {
109  G4EmSaturation* emSaturation = G4LossTableManager::Instance()->EmSaturation();
110  if (emSaturation)
111  {
112  if (m_Verbosity)
113  {
114  emSaturation->SetVerbose(m_Verbosity);
115  }
116  double visen = emSaturation->VisibleEnergyDepositionAtAStep(step) / GeV;
117  return visen;
118  }
119  else
120  {
121  std::cout
122  << "PHG4SteppingAction::GetScintLightYield - ERROR - can NOT initialize G4EmSaturation!"
123  << std::endl;
124 
125  return 0.;
126  }
127 }
128 
129 void PHG4SteppingAction::StoreLocalCoordinate(PHG4Hit* hit, const G4Step* aStep,
130  const bool do_prepoint, const bool do_postpoint)
131 {
132  assert(hit);
133  assert(aStep);
134 
135  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
136  const G4TouchableHandle& theTouchable = preStepPoint->GetTouchableHandle();
137 
138  if (do_prepoint)
139  {
140  const G4ThreeVector& worldPosition = preStepPoint->GetPosition();
141  G4ThreeVector localPosition =
142  theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
143 
144  hit->set_local_x(0, localPosition.x() / cm);
145  hit->set_local_y(0, localPosition.y() / cm);
146  hit->set_local_z(0, localPosition.z() / cm);
147  }
148  if (do_postpoint)
149  {
150  G4StepPoint* postPoint = aStep->GetPostStepPoint();
151 
152  const G4ThreeVector& worldPosition = postPoint->GetPosition();
153  G4ThreeVector localPosition =
154  theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
155 
156  hit->set_local_x(1, localPosition.x() / cm);
157  hit->set_local_y(1, localPosition.y() / cm);
158  hit->set_local_z(1, localPosition.z() / cm);
159  }
160 }
161 
162 void PHG4SteppingAction::SetLightCorrection(const double inner_radius, const double inner_corr, const double outer_radius, const double outer_corr)
163 {
165  m_LightBalanceInnerCorr = inner_corr;
167  m_LightBalanceOuterCorr = outer_corr;
168  return;
169 }
170 
171 double PHG4SteppingAction::GetLightCorrection(const double xpos, const double ypos) const
172 {
173  double r = sqrt(xpos * xpos + ypos * ypos);
174  double correction = GetLightCorrection(r);
175  return correction;
176 }
177 
178 double PHG4SteppingAction::GetLightCorrection(const double r) const
179 {
180  double correction = 1.;
181  if (ValidCorrection())
182  {
185  correction = slope * r + b;
186  correction = std::min(correction, 1.);
187  correction = std::max(correction, 0.);
188  }
189  return correction;
190 }
191 
193 {
194  if (std::isfinite(m_LightBalanceOuterRadius) &&
195  std::isfinite(m_LightBalanceInnerRadius) &&
196  std::isfinite(m_LightBalanceOuterCorr) &&
197  std::isfinite(m_LightBalanceInnerCorr))
198  {
199  return true;
200  }
201  return false;
202 }