MaterialEffects/src/MultipleScatteringUpdator.cc

0001 #include "TrackingTools/MaterialEffects/interface/MultipleScatteringUpdator.h"
0002 #include "DataFormats/GeometrySurface/interface/MediumProperties.h"
0003
0004 #include "DataFormats/Math/interface/approx_log.h"
0005
0006 void oldMUcompute(const TrajectoryStateOnSurface& TSoS, const PropagationDirection propDir, double mass, double ptmin);
0007
0008 //#define DBG_MSU
0009
0010 //
0011 // Computation of contribution of multiple scatterning to covariance matrix
0012 //   of local parameters based on Highland formula for sigma(alpha) in plane.
0013 //
0014 void MultipleScatteringUpdator::compute(const TrajectoryStateOnSurface& TSoS,
0015                                         const PropagationDirection propDir,
0016                                         Effect& effect) const {
0017   //
0018   // Get surface
0019   //
0020   const Surface& surface = TSoS.surface();
0021   //
0022   //
0023   // Now get information on medium
0024   //
0025   const MediumProperties& mp = surface.mediumProperties();
0026   if UNLIKELY (mp.radLen() == 0)
0027     return;
0028
0029   // Momentum vector
0030   LocalVector d = TSoS.localMomentum();
0031   float p2 = d.mag2();
0032   d *= 1.f / sqrt(p2);
0033   float xf = 1.f / std::abs(d.z());  // increase of path due to angle of incidence
0034   // calculate general physics things
0035   constexpr float amscon = 1.8496e-4;  // (13.6MeV)**2
0036   const float m2 = mass() * mass();    // use mass hypothesis from constructor
0037   float e2 = p2 + m2;
0038   float beta2 = p2 / e2;
0039   // calculate the multiple scattering angle
0040   float radLen = mp.radLen() * xf;  // effective rad. length
0041   float sigt2 = 0.;                 // sigma(alpha)**2
0042
0043   // Calculated rms scattering angle squared.
0044   float fact = 1.f + 0.038f * unsafe_logf<2>(radLen);
0045   fact *= fact;
0046   float a = fact / (beta2 * p2);
0047   sigt2 = amscon * radLen * a;
0048
0049   if (thePtMin > 0) {
0050 #ifdef DBG_MSU
0051     std::cout << "Original rms scattering = " << sqrt(sigt2);
0052 #endif
0053     // Inflate estimated rms scattering angle, to take into account
0054     // that 1/p is not known precisely.
0055     AlgebraicSymMatrix55 const& covMatrix = TSoS.localError().matrix();
0056     float error2_QoverP = covMatrix(0, 0);
0057     // Formula valid for ultra-relativistic particles.
0058     //      sigt2 *= (1. + p2 * error2_QoverP);
0059     // Exact formula
0060     sigt2 *= 1.f + error2_QoverP * (p2 + m2 * beta2 * (5.f + 3.f * p2 * error2_QoverP));
0061 #ifdef DBG_MSU
0062     std::cout << " new = " << sqrt(sigt2);
0063 #endif
0064     // Convert Pt constraint to P constraint, neglecting uncertainty in
0065     // track angle.
0066     float pMin2 = thePtMin * thePtMin * (p2 / TSoS.globalMomentum().perp2());
0067     // Use P constraint to calculate rms maximum scattering angle.
0068     float betaMin2 = pMin2 / (pMin2 + m2);
0069     float a_max = fact / (betaMin2 * pMin2);
0070     float sigt2_max = amscon * radLen * a_max;
0071     if (sigt2 > sigt2_max)
0072       sigt2 = sigt2_max;
0073 #ifdef DBG_MSU
0074     std::cout << " after P constraint (" << pMin << ") = " << sqrt(sigt2);
0075     std::cout << " for track with 1/p=" << 1 / p << "+-" << sqrt(error2_QoverP) << std::endl;
0076 #endif
0077   }
0078
0079   float isl2 = 1.f / d.perp2();
0080   float cl2 = (d.z() * d.z());
0081   float cf2 = (d.x() * d.x()) * isl2;
0082   float sf2 = (d.y() * d.y()) * isl2;
0083   // Create update (transformation of independant variations
0084   //   on angle in orthogonal planes to local parameters.
0085   float den = 1.f / (cl2 * cl2);
0086   using namespace materialEffect;
0087   effect.deltaCov[msxx] += (den * sigt2) * (sf2 * cl2 + cf2);
0088   effect.deltaCov[msxy] += (den * sigt2) * (d.x() * d.y());
0089   effect.deltaCov[msyy] += (den * sigt2) * (cf2 * cl2 + sf2);
0090
0091   /*
0092     std::cout << "new " <<  theDeltaCov(1,1) << " " <<  theDeltaCov(1,2)  << " " <<  theDeltaCov(2,2) << std::endl;
0093     oldMUcompute(TSoS,propDir, mass(), thePtMin);
0094   */
0095 }
0096
0097 //
0098 // Computation of contribution of multiple scatterning to covariance matrix
0099 //   of local parameters based on Highland formula for sigma(alpha) in plane.
0100 //
0101 void oldMUcompute(const TrajectoryStateOnSurface& TSoS, const PropagationDirection propDir, float mass, float thePtMin) {
0102   //
0103   // Get surface
0104   //
0105   const Surface& surface = TSoS.surface();
0106   //
0107   //
0108   // Now get information on medium
0109   //
0110   if (surface.mediumProperties().isValid()) {
0111     // Momentum vector
0112     LocalVector d = TSoS.localMomentum();
0113     float p = d.mag();
0114     d *= 1. / p;
0115     // MediumProperties mp(0.02, .5e-4);
0116     const MediumProperties& mp = surface.mediumProperties();
0117     float xf = 1. / fabs(d.z());  // increase of path due to angle of incidence
0118     // calculate general physics things
0119     const float amscon = 1.8496e-4;  // (13.6MeV)**2
0120     const float m = mass;            // use mass hypothesis from constructor
0121     float e = sqrt(p * p + m * m);
0122     float beta = p / e;
0123     // calculate the multiple scattering angle
0124     float radLen = mp.radLen() * xf;  // effective rad. length
0125     float sigt2 = 0.;                 // sigma(alpha)**2
0126     if (radLen > 0) {
0127       // Calculated rms scattering angle squared.
0128       float a = (1. + 0.038 * log(radLen)) / (beta * p);
0129       a *= a;
0130       sigt2 = amscon * radLen * a;
0131       if (thePtMin > 0) {
0132 #ifdef DBG_MSU
0133         std::cout << "Original rms scattering = " << sqrt(sigt2);
0134 #endif
0135         // Inflate estimated rms scattering angle, to take into account
0136         // that 1/p is not known precisely.
0137         AlgebraicSymMatrix55 covMatrix = TSoS.localError().matrix();
0138         float error2_QoverP = covMatrix(0, 0);
0139         // Formula valid for ultra-relativistic particles.
0140         //      sigt2 *= (1. + (p*p) * error2_QoverP);
0141         // Exact formula
0142         sigt2 *= (1. + (p * p) * error2_QoverP * (1. + 5 * m * m / (e * e) + 3 * m * m * beta * beta * error2_QoverP));
0143 #ifdef DBG_MSU
0144         std::cout << " new = " << sqrt(sigt2);
0145 #endif
0146         // Convert Pt constraint to P constraint, neglecting uncertainty in
0147         // track angle.
0148         float pMin = thePtMin * (TSoS.globalMomentum().mag() / TSoS.globalMomentum().perp());
0149         // Use P constraint to calculate rms maximum scattering angle.
0150         float betaMin = pMin / sqrt(pMin * pMin + m * m);
0151         float a_max = (1. + 0.038 * log(radLen)) / (betaMin * pMin);
0152         a_max *= a_max;
0153         float sigt2_max = amscon * radLen * a_max;
0154         if (sigt2 > sigt2_max)
0155           sigt2 = sigt2_max;
0156 #ifdef DBG_MSU
0157         std::cout << " after P constraint (" << pMin << ") = " << sqrt(sigt2);
0158         std::cout << " for track with 1/p=" << 1 / p << "+-" << sqrt(error2_QoverP) << std::endl;
0159 #endif
0160       }
0161     }
0162     float sl = d.perp();
0163     float cl = d.z();
0164     float cf = d.x() / sl;
0165     float sf = d.y() / sl;
0166     // Create update (transformation of independant variations
0167     //   on angle in orthogonal planes to local parameters.
0168     std::cout << "old " << sigt2 * (sf * sf * cl * cl + cf * cf) / (cl * cl * cl * cl) << " "
0169               << sigt2 * (cf * sf * sl * sl) / (cl * cl * cl * cl) << " "
0170               << sigt2 * (cf * cf * cl * cl + sf * sf) / (cl * cl * cl * cl) << std::endl;
0171   }
0172 }