10 #include <tracking/trackFindingCDC/eventdata/utils/EnergyLossEstimator.h>
12 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCBFieldUtil.h>
14 #include <tracking/trackFindingCDC/numerics/ESign.h>
16 #include <geometry/GeometryManager.h>
18 #include <framework/gearbox/Unit.h>
19 #include <framework/logging/Logger.h>
21 #include "G4ThreeVector.hh"
22 #include "G4Navigator.hh"
23 #include "G4VPhysicalVolume.hh"
24 #include "G4Material.hh"
27 using namespace TrackFindingCDC;
35 g4Nav.SetWorldVolume(g4World);
42 G4ThreeVector g4Pos(posX * CLHEP::cm, posY * CLHEP::cm, posZ * CLHEP::cm);
43 const G4VPhysicalVolume* g4Volume = g4Nav.LocateGlobalPointAndSetup(g4Pos);
44 const G4Material* g4Mat = g4Volume->GetLogicalVolume()->GetMaterial();
49 if (g4Mat->GetNumberOfElements() == 1) {
50 A = g4Mat->GetA() * CLHEP::mole / CLHEP::g;
54 for (
unsigned i = 0; i < g4Mat->GetNumberOfElements(); ++i) {
55 const G4Element* element = (*g4Mat->GetElementVector())[i];
56 const double elementA = element->GetA() * CLHEP::mole / CLHEP::g;
57 const double elementZ = element->GetZ();
58 const double frac = g4Mat->GetFractionVector()[i];
59 B2RESULT(
"Part " << i <<
" Z=" << elementZ <<
" A=" << elementA <<
" (" << frac <<
")");
67 const double density = g4Mat->GetDensity() / CLHEP::g * CLHEP::cm3;
68 const double radiationLength = g4Mat->GetRadlen() / CLHEP::cm;
69 const double mEE = g4Mat->GetIonisation()->GetMeanExcitationEnergy() / CLHEP::GeV;
71 B2RESULT(
"Received Z " << Z);
72 B2RESULT(
"Received A " << A);
73 B2RESULT(
"Received density " << density);
74 B2RESULT(
"Received radiation length " << radiationLength);
75 B2RESULT(
"Received mean excitation energy " << mEE);
77 const double eDensity = Z * density / A;
78 B2RESULT(
"Received electron density " << eDensity);
87 if (std::abs(pdgCode) == 11) {
88 const double electronMass = 0.511 *
Unit::MeV;
90 }
else if (std::abs(pdgCode) == 13) {
91 const double muonMass = 195.658 *
Unit::MeV;
93 }
else if (std::abs(pdgCode) == 321) {
94 const double kaonMass = 0.493677;
96 }
else if (std::abs(pdgCode) == 211) {
97 const double pionMass = 0.13957;
99 }
else if (std::abs(pdgCode) == 2212) {
100 const double protonMass = 0.938272;
108 if (std::abs(pdgCode) == 11) {
109 return -sign(pdgCode);
110 }
else if (std::abs(pdgCode) == 13) {
111 return -sign(pdgCode);
112 }
else if (std::abs(pdgCode) == 321) {
113 return sign(pdgCode);
114 }
else if (std::abs(pdgCode) == 211) {
115 return sign(pdgCode);
116 }
else if (std::abs(pdgCode) == 2212) {
117 return sign(pdgCode);
123 : m_eDensity(eDensity)
131 static const double eMass =
getMass(11);
133 const double M =
getMass(pdgCode);
135 const double E = std::sqrt(p * p + M * M);
136 const double gamma = E / M;
137 const double beta = p / E;
139 const double beta2 = beta * beta;
140 const double gamma2 = gamma * gamma;
142 const double Wmax = 2 * eMass * beta2 * gamma2 / (1 + 2 * gamma * eMass / M);
143 const double I2 =
m_I *
m_I;
147 (1.0 / 2.0 * std::log(2 * eMass * beta2 * gamma2 * Wmax / I2) - beta2);
154 const double eLoss = arcLength * dEdx;
161 return (p - eLoss) / p;
163 const double mass =
getMass(pdgCode);
165 const double eBefore = std::sqrt(p * p + mass * mass);
166 const double eAfter = eBefore - eLoss;
167 const double pAfter = std::sqrt(eAfter * eAfter - mass * mass);
173 const double eLoss =
getEnergyLoss(pt, pdgCode, arcLength2D);
177 return radius * eLoss / (pt - eLoss);