Belle II Software  release-05-01-25
EnergyLossEstimator.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Oliver Frost *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #include <tracking/trackFindingCDC/eventdata/utils/EnergyLossEstimator.h>
11 
12 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCBFieldUtil.h>
13 
14 #include <tracking/trackFindingCDC/numerics/ESign.h>
15 
16 #include <geometry/GeometryManager.h>
17 
18 #include <framework/gearbox/Unit.h>
19 #include <framework/logging/Logger.h>
20 
21 #include "G4ThreeVector.hh"
22 #include "G4Navigator.hh"
23 #include "G4VPhysicalVolume.hh"
24 #include "G4Material.hh"
25 
26 using namespace Belle2;
27 using namespace TrackFindingCDC;
28 
30 {
31 
32  // Look up material properties from the Belle2 geometry
33  G4Navigator g4Nav;
34  G4VPhysicalVolume* g4World = geometry::GeometryManager::getInstance().getTopVolume();
35  g4Nav.SetWorldVolume(g4World);
36 
37  // Choose position well inside the CDC
38  double posX = 0;
39  double posY = 50;
40  double posZ = 0;
41 
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();
45 
46  double A = 0;
47  double Z = 0;
48 
49  if (g4Mat->GetNumberOfElements() == 1) {
50  A = g4Mat->GetA() * CLHEP::mole / CLHEP::g;
51  Z = g4Mat->GetZ();
52  } else {
53  // Calculate weight-averaged A, Z
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 << ")");
60  Z += elementZ * frac;
61  A += elementA * frac;
62  }
63  }
64 
65  // Make sure to translate to Belle units from the Geant4 / CLHEP units
66  // Z has correct units (is unitless)
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;
70 
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);
76 
77  const double eDensity = Z * density / A;
78  B2RESULT("Received electron density " << eDensity);
79 
80  const double bZ = CDCBFieldUtil::getBFieldZ();
81 
82  return EnergyLossEstimator(eDensity, mEE, bZ);
83 }
84 
85 double EnergyLossEstimator::getMass(int pdgCode)
86 {
87  if (std::abs(pdgCode) == 11) {
88  const double electronMass = 0.511 * Unit::MeV;
89  return electronMass;
90  } else if (std::abs(pdgCode) == 13) {
91  const double muonMass = 195.658 * Unit::MeV;
92  return muonMass;
93  } else if (std::abs(pdgCode) == 321) {
94  const double kaonMass = 0.493677;
95  return kaonMass;
96  } else if (std::abs(pdgCode) == 211) {
97  const double pionMass = 0.13957;
98  return pionMass;
99  } else if (std::abs(pdgCode) == 2212) {
100  const double protonMass = 0.938272;
101  return protonMass;
102  }
103  return NAN;
104 }
105 
107 {
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);
118  }
119  return 0;
120 }
121 
122 EnergyLossEstimator::EnergyLossEstimator(double eDensity, double I, double bZ)
123  : m_eDensity(eDensity)
124  , m_I(I)
125  , m_bZ(bZ)
126 {
127 }
128 
129 double EnergyLossEstimator::getBetheStoppingPower(double p, int pdgCode) const
130 {
131  static const double eMass = getMass(11);
132 
133  const double M = getMass(pdgCode);
134 
135  const double E = std::sqrt(p * p + M * M);
136  const double gamma = E / M;
137  const double beta = p / E;
138 
139  const double beta2 = beta * beta;
140  const double gamma2 = gamma * gamma;
141 
142  const double Wmax = 2 * eMass * beta2 * gamma2 / (1 + 2 * gamma * eMass / M);
143  const double I2 = m_I * m_I;
144 
145  static const double K = 0.307075 * Unit::MeV * Unit::cm2;
146  const double dEdx = K * m_eDensity / beta2 *
147  (1.0 / 2.0 * std::log(2 * eMass * beta2 * gamma2 * Wmax / I2) - beta2);
148  return dEdx;
149 }
150 
151 double EnergyLossEstimator::getEnergyLoss(double p, int pdgCode, double arcLength) const
152 {
153  const double dEdx = getBetheStoppingPower(p, pdgCode);
154  const double eLoss = arcLength * dEdx;
155  return eLoss;
156 }
157 
158 double EnergyLossEstimator::getMomentumLossFactor(double p, int pdgCode, double arcLength) const
159 {
160  const double eLoss = getEnergyLoss(p, pdgCode, arcLength);
161  return (p - eLoss) / p;
162 
163  const double mass = getMass(pdgCode);
164 
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);
168  return pAfter / p;
169 }
170 
171 double EnergyLossEstimator::getLossDist2D(double pt, int pdgCode, double arcLength2D) const
172 {
173  const double eLoss = getEnergyLoss(pt, pdgCode, arcLength2D);
174 
175  const int q = getCharge(pdgCode);
176  const double radius = q * CDCBFieldUtil::absMom2DToBendRadius(pt, m_bZ);;
177  return radius * eLoss / (pt - eLoss);
178 
179  // const double pFactor = getMomentumLossFactor(p, pdgCode, arcLength2D);
180  // return radius * (1 / pFactor - 1);
181 }
Belle2::geometry::GeometryManager::getTopVolume
G4VPhysicalVolume * getTopVolume()
Return a pointer to the top volume.
Definition: GeometryManager.h:59
Belle2::TrackFindingCDC::EnergyLossEstimator
Helper struct to provide consistent energy loss estimation throughout the CDC track finding.
Definition: EnergyLossEstimator.h:32
Belle2::geometry::GeometryManager::getInstance
static GeometryManager & getInstance()
Return a reference to the instance.
Definition: GeometryManager.cc:98
Belle2::TrackFindingCDC::EnergyLossEstimator::m_bZ
double m_bZ
B field to be used for the distance translation.
Definition: EnergyLossEstimator.h:91
Belle2::TrackFindingCDC::EnergyLossEstimator::getLossDist2D
double getLossDist2D(double pt, int pdgCode, double arcLength2D) const
Calculates a correction term for the two dimensional distance undoing the energy loss after the given...
Definition: EnergyLossEstimator.cc:171
Belle2::Unit::MeV
static const double MeV
[megaelectronvolt]
Definition: Unit.h:124
Belle2::TrackFindingCDC::EnergyLossEstimator::getBetheStoppingPower
double getBetheStoppingPower(double p, int pdgCode) const
Stopping power aka energy loss / arc length.
Definition: EnergyLossEstimator.cc:129
Belle2::TrackFindingCDC::EnergyLossEstimator::forCDC
static EnergyLossEstimator forCDC()
Create an energy loss estimator with the material properties of the CDC.
Definition: EnergyLossEstimator.cc:29
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TrackFindingCDC::EnergyLossEstimator::getMomentumLossFactor
double getMomentumLossFactor(double p, int pdgCode, double arcLength) const
Calculates a factor applicable scaling the current momentum to the momentum after traveling given arc...
Definition: EnergyLossEstimator.cc:158
Belle2::TrackFindingCDC::EnergyLossEstimator::getCharge
static int getCharge(int pdgCode)
Lookup the charge for the given pdg code.
Definition: EnergyLossEstimator.cc:106
Belle2::TrackFindingCDC::CDCBFieldUtil::getBFieldZ
static double getBFieldZ()
Getter for the signed magnetic field stength in z direction at the origin ( in Tesla )
Definition: CDCBFieldUtil.cc:40
Belle2::TrackFindingCDC::EnergyLossEstimator::EnergyLossEstimator
EnergyLossEstimator(double eDensity, double I, double bZ=NAN)
Constructor from the material properties.
Definition: EnergyLossEstimator.cc:122
Belle2::TrackFindingCDC::EnergyLossEstimator::getMass
static double getMass(int pdgCode)
Lookup the mass for the given pdg code.
Definition: EnergyLossEstimator.cc:85
Belle2::Unit::cm2
static const double cm2
[square centimeters]
Definition: Unit.h:88
Belle2::TrackFindingCDC::EnergyLossEstimator::getEnergyLoss
double getEnergyLoss(double p, int pdgCode, double arcLength) const
Calculates the total energy loss after travelling the given distance.
Definition: EnergyLossEstimator.cc:151
Belle2::TrackFindingCDC::EnergyLossEstimator::m_eDensity
double m_eDensity
Electron density in mol / cm^3.
Definition: EnergyLossEstimator.h:85
Belle2::TrackFindingCDC::CDCBFieldUtil::absMom2DToBendRadius
static double absMom2DToBendRadius(double absMom2D, double bZ)
Conversion helper for momenta to two dimensional (absolute) bend radius.
Definition: CDCBFieldUtil.cc:92
Belle2::TrackFindingCDC::EnergyLossEstimator::m_I
double m_I
Mean excitation energy in GeV.
Definition: EnergyLossEstimator.h:88