GsfTracking/src/GsfMultipleScatteringUpdator.cc

0001 #include "TrackingTools/GsfTracking/interface/GsfMultipleScatteringUpdator.h"
0002
0003 #include "DataFormats/GeometrySurface/interface/MediumProperties.h"
0004
0005 // #include "CommonDet/DetUtilities/interface/DetExceptions.h"
0006
0007 // #include "Utilities/Notification/interface/Verbose.h"
0008
0009 void GsfMultipleScatteringUpdator::compute(const TrajectoryStateOnSurface& TSoS,
0010                                            const PropagationDirection propDir,
0011                                            Effect effects[]) const {
0012   //
0013   // Get surface and check presence of medium properties
0014   //
0015   const Surface& surface = TSoS.surface();
0016   //
0017   // calculate components
0018   //
0019   if (surface.mediumProperties().isValid()) {
0020     LocalVector pvec = TSoS.localMomentum();
0021     float p = TSoS.localMomentum().mag();
0022     pvec *= 1. / p;
0023     // thickness in radiation lengths
0024     float rl = surface.mediumProperties().radLen() / fabs(pvec.z());
0025     // auxiliary variables for modified X0
0026     constexpr float z = 14;  // atomic number of silicon
0027     const float logz = log(z);
0028     const float h = (z + 1) / z * log(287 * sqrt(z)) / log(159 * pow(z, -1. / 3.));
0029     float beta2 = 1. / (1. + mass() * mass() / p / p);
0030     // reduced thickness
0031     float dp1 = rl / beta2 / h;
0032     float logdp1 = log(dp1);
0033     float logdp2 = 2. / 3. * logz + logdp1;
0034     // weights
0035     float w2;
0036     if (logdp2 < log(0.5))
0037       w2 = 0.05283 + 0.0077 * logdp2 + 0.00069 * logdp2 * logdp2;
0038     else
0039       w2 = -0.01517 + 0.1151 * logdp2 - 0.00653 * logdp2 * logdp2;
0040     float w1 = 1. - w2;
0041     effects[0].weight *= w1;
0042     effects[1].weight *= w2;
0043     // reduced variances
0044     float var1 = 0.8510 + 0.03314 * logdp1 - 0.001825 * logdp1 * logdp1;
0045     float var2 = (1. - w1 * var1) / w2;
0046     for (int ic = 0; ic < 2; ic++) {
0047       // choose component and multiply with total variance
0048       float var = ic == 0 ? var1 : var2;
0049       var *= 225.e-6 * dp1 / p / p;
0050       AlgebraicSymMatrix55 cov;
0051       // transform from orthogonal planes containing the
0052       // momentum vector to local parameters
0053       float sl = pvec.perp();
0054       float cl = pvec.z();
0055       float cf = pvec.x() / sl;
0056       float sf = pvec.y() / sl;
0057       using namespace materialEffect;
0058       effects[ic].deltaCov[msxx] += var * (sf * sf * cl * cl + cf * cf) / (cl * cl * cl * cl);
0059       effects[ic].deltaCov[msxy] += var * (cf * sf * sl * sl) / (cl * cl * cl * cl);
0060       effects[ic].deltaCov[msyy] += var * (cf * cf * cl * cl + sf * sf) / (cl * cl * cl * cl);
0061     }
0062   }
0063 }