MaterialEffects/src/EnergyLossUpdator.cc

0001 #include "TrackingTools/MaterialEffects/interface/EnergyLossUpdator.h"
0002 #include "DataFormats/GeometrySurface/interface/MediumProperties.h"
0003
0004 #include "DataFormats/Math/interface/approx_exp.h"
0005 #include "DataFormats/Math/interface/approx_log.h"
0006
0007 void oldComputeBetheBloch(const LocalVector& localP, const MediumProperties& materialConstants, double mass);
0008 void oldComputeElectrons(const LocalVector& localP,
0009                          const MediumProperties& materialConstants,
0010                          const PropagationDirection propDir);
0011
0012 //
0013 // Computation of contribution of energy loss to momentum and covariance
0014 // matrix of local parameters based on Bethe-Bloch. For electrons
0015 // contribution of radiation acc. to Bethe & Heitler.
0016 //
0017 void EnergyLossUpdator::compute(const TrajectoryStateOnSurface& TSoS,
0018                                 const PropagationDirection propDir,
0019                                 Effect& effect) const {
0020   //
0021   // Get surface
0022   //
0023   const Surface& surface = TSoS.surface();
0024   //
0025   //
0026   // Now get information on medium
0027   //
0028   if (surface.mediumProperties().isValid()) {
0029     //
0030     // Bethe-Bloch
0031     //
0032     if (mass() > 0.001)
0033       computeBetheBloch(TSoS.localMomentum(), surface.mediumProperties(), effect);
0034     //
0035     // Special treatment for electrons (currently rather crude
0036     // distinction using mass)
0037     //
0038     else
0039       computeElectrons(TSoS.localMomentum(), surface.mediumProperties(), propDir, effect);
0040     if (propDir != alongMomentum)
0041       effect.deltaP *= -1.;
0042   }
0043 }
0044 //
0045 // Computation of energy loss according to Bethe-Bloch
0046 //
0047 void EnergyLossUpdator::computeBetheBloch(const LocalVector& localP,
0048                                           const MediumProperties& materialConstants,
0049                                           Effect& effect) const {
0050   //
0051   // calculate absolute momentum and correction to path length from angle
0052   // of incidence
0053   //
0054
0055   typedef float Float;
0056
0057   Float p2 = localP.mag2();
0058   Float xf = std::abs(std::sqrt(p2) / localP.z());
0059
0060   // constants
0061   const Float m2 = mass() * mass();  // use mass hypothesis from constructor
0062
0063   constexpr Float emass = 0.511e-3;
0064   constexpr Float poti = 16.e-9 * 10.75;                 // = 16 eV * Z**0.9, for Si Z=14
0065   const Float eplasma = 28.816e-9 * sqrt(2.33 * 0.498);  // 28.816 eV * sqrt(rho*(Z/A)) for Si
0066   const Float delta0 = 2 * log(eplasma / poti) - 1.;
0067
0068   // calculate general physics things
0069   Float im2 = Float(1.) / m2;
0070   Float e2 = p2 + m2;
0071   Float e = std::sqrt(e2);
0072   Float beta2 = p2 / e2;
0073   Float eta2 = p2 * im2;
0074   Float ratio2 = (emass * emass) * im2;
0075   Float emax = Float(2.) * emass * eta2 / (Float(1.) + Float(2.) * emass * e * im2 + ratio2);
0076
0077   Float xi = materialConstants.xi() * xf;
0078   xi /= beta2;
0079
0080   Float dEdx = xi * (unsafe_logf<2>(Float(2.) * emass * emax / (poti * poti)) - Float(2.) * (beta2)-delta0);
0081
0082   Float dEdx2 = xi * emax * (Float(1.) - Float(0.5) * beta2);
0083   Float dP = dEdx / std::sqrt(beta2);
0084   Float sigp2 = dEdx2 / (beta2 * p2 * p2);
0085   effect.deltaP += -dP;
0086   using namespace materialEffect;
0087   effect.deltaCov[elos] += sigp2;
0088
0089   // std::cout << "pion new " <<  theDeltaP << " " << theDeltaCov(0,0) << std::endl;
0090   // oldComputeBetheBloch (localP, materialConstants, mass());
0091 }
0092 //
0093 // Computation of energy loss for electrons
0094 //
0095 void EnergyLossUpdator::computeElectrons(const LocalVector& localP,
0096                                          const MediumProperties& materialConstants,
0097                                          const PropagationDirection propDir,
0098                                          Effect& effect) const {
0099   //
0100   // calculate absolute momentum and correction to path length from angle
0101   // of incidence
0102   //
0103   float p2 = localP.mag2();
0104   float p = std::sqrt(p2);
0105   float normalisedPath = std::abs(p / localP.z()) * materialConstants.radLen();
0106   //
0107   // Energy loss and variance according to Bethe and Heitler, see also
0108   // Comp. Phys. Comm. 79 (1994) 157.
0109   //
0110   const float l3ol2 = std::log(3.) / std::log(2.);
0111   float z = unsafe_expf<3>(-normalisedPath);
0112   float varz = unsafe_expf<3>(-normalisedPath * l3ol2) - z * z;
0113   // exp(-2*normalisedPath);
0114
0115   if (propDir == oppositeToMomentum) {
0116     //
0117     // for backward propagation: delta(1/p) is linear in z=p_outside/p_inside
0118     // convert to obtain equivalent delta(p). Sign of deltaP is corrected
0119     // in method compute -> deltaP<0 at this place!!!
0120     //
0121     effect.deltaP += -p * (1.f / z - 1.f);
0122     using namespace materialEffect;
0123     effect.deltaCov[elos] += varz / p2;
0124   } else {
0125     //
0126     // for forward propagation: calculate in p (linear in 1/z=p_inside/p_outside),
0127     // then convert sig(p) to sig(1/p).
0128     //
0129     effect.deltaP += p * (z - 1.f);
0130     //    float f = 1/p/z/z;
0131     // patch to ensure consistency between for- and backward propagation
0132     float f2 = 1.f / (p2 * z * z);
0133     using namespace materialEffect;
0134     effect.deltaCov[elos] += f2 * varz;
0135   }
0136
0137   // std::cout << "electron new " <<  theDeltaP << " " << theDeltaCov(0,0) << std::endl;
0138   // oldComputeElectrons (localP, materialConstants, propDir);
0139 }
0140
0141 ///
0142
0143 void oldComputeBetheBloch(const LocalVector& localP, const MediumProperties& materialConstants, double mass) {
0144   double theDeltaP = 0., theDeltaCov = 0;
0145
0146   //
0147   // calculate absolute momentum and correction to path length from angle
0148   // of incidence
0149   //
0150   double p = localP.mag();
0151   double xf = fabs(p / localP.z());
0152   // constants
0153   const double m = mass;  // use mass hypothesis from constructor
0154
0155   const double emass = 0.511e-3;
0156   const double poti = 16.e-9 * 10.75;                     // = 16 eV * Z**0.9, for Si Z=14
0157   const double eplasma = 28.816e-9 * sqrt(2.33 * 0.498);  // 28.816 eV * sqrt(rho*(Z/A)) for Si
0158   const double delta0 = 2 * log(eplasma / poti) - 1.;
0159
0160   // calculate general physics things
0161   double e = sqrt(p * p + m * m);
0162   double beta = p / e;
0163   double gamma = e / m;
0164   double eta2 = beta * gamma;
0165   eta2 *= eta2;
0166   //  double lnEta2 = log(eta2);
0167   double ratio = emass / m;
0168   double emax = 2. * emass * eta2 / (1. + 2. * ratio * gamma + ratio * ratio);
0169   // double delta = delta0 + lnEta2;
0170
0171   // calculate the mean and sigma of energy loss
0172   // xi = d[g/cm2] * 0.307075MeV/(g/cm2) * Z/A * 1/2
0173   double xi = materialConstants.xi() * xf;
0174   xi /= (beta * beta);
0175
0176   //   double dEdx = xi*(log(2.*emass*eta2*emax/(poti*poti)) - 2.*(beta*beta));
0177   //double dEdx = xi*(log(2.*emass*emax/(poti*poti))+lnEta2 - 2.*(beta*beta) - delta);
0178
0179   double dEdx = xi * (log(2. * emass * emax / (poti * poti)) - 2. * (beta * beta) - delta0);
0180   double dEdx2 = xi * emax * (1. - 0.5 * (beta * beta));
0181   double dP = dEdx / beta;
0182   double sigp2 = dEdx2 * e * e / (p * p * p * p * p * p);
0183   theDeltaP += -dP;
0184   theDeltaCov += sigp2;
0185
0186   std::cout << "pion old " << theDeltaP << " " << theDeltaCov << std::endl;
0187 }
0188
0189 void oldComputeElectrons(const LocalVector& localP,
0190                          const MediumProperties& materialConstants,
0191                          const PropagationDirection propDir) {
0192   double theDeltaP = 0., theDeltaCov = 0;
0193
0194   //
0195   // calculate absolute momentum and correction to path length from angle
0196   // of incidence
0197   //
0198   double p = localP.mag();
0199   double normalisedPath = fabs(p / localP.z()) * materialConstants.radLen();
0200   //
0201   // Energy loss and variance according to Bethe and Heitler, see also
0202   // Comp. Phys. Comm. 79 (1994) 157.
0203   //
0204   double z = exp(-normalisedPath);
0205   double varz = exp(-normalisedPath * log(3.) / log(2.)) - z * z;
0206   // exp(-2*normalisedPath);
0207
0208   if (propDir == oppositeToMomentum) {
0209     //
0210     // for backward propagation: delta(1/p) is linear in z=p_outside/p_inside
0211     // convert to obtain equivalent delta(p). Sign of deltaP is corrected
0212     // in method compute -> deltaP<0 at this place!!!
0213     //
0214     theDeltaP += -p * (1 / z - 1);
0215     theDeltaCov += varz / (p * p);
0216   } else {
0217     //
0218     // for forward propagation: calculate in p (linear in 1/z=p_inside/p_outside),
0219     // then convert sig(p) to sig(1/p).
0220     //
0221     theDeltaP += p * (z - 1);
0222     //    double f = 1/p/z/z;
0223     // patch to ensure consistency between for- and backward propagation
0224     double f = 1. / (p * z);
0225     theDeltaCov += f * f * varz;
0226   }
0227
0228   std::cout << "electron old " << theDeltaP << " " << theDeltaCov << std::endl;
0229 }