GsfMultipleScatteringUpdator.cc

Go to the documentation of this file.

#include "TrackingTools/GsfTracking/interface/GsfMultipleScatteringUpdator.h"

#include "DataFormats/GeometrySurface/interface/MediumProperties.h"

// #include "CommonDet/DetUtilities/interface/DetExceptions.h"

// #include "Utilities/Notification/interface/Verbose.h"

void GsfMultipleScatteringUpdator::compute(const TrajectoryStateOnSurface& TSoS,
                                           const PropagationDirection propDir,
                                           Effect effects[]) const {
  //
  // Get surface and check presence of medium properties
  //
  const Surface& surface = TSoS.surface();
  //
  // calculate components
  //
  if (surface.mediumProperties().isValid()) {
    LocalVector pvec = TSoS.localMomentum();
    float p = TSoS.localMomentum().mag();
    pvec *= 1. / p;
    // thickness in radiation lengths
    float rl = surface.mediumProperties().radLen() / fabs(pvec.z());
    // auxiliary variables for modified X0
    constexpr float z = 14;  // atomic number of silicon
    const float logz = log(z);
    const float h = (z + 1) / z * log(287 * sqrt(z)) / log(159 * pow(z, -1. / 3.));
    float beta2 = 1. / (1. + mass() * mass() / p / p);
    // reduced thickness
    float dp1 = rl / beta2 / h;
    float logdp1 = log(dp1);
    float logdp2 = 2. / 3. * logz + logdp1;
    // weights
    float w2;
    if (logdp2 < log(0.5))
      w2 = 0.05283 + 0.0077 * logdp2 + 0.00069 * logdp2 * logdp2;
    else
      w2 = -0.01517 + 0.1151 * logdp2 - 0.00653 * logdp2 * logdp2;
    float w1 = 1. - w2;
    effects[0].weight *= w1;
    effects[1].weight *= w2;
    // reduced variances
    float var1 = 0.8510 + 0.03314 * logdp1 - 0.001825 * logdp1 * logdp1;
    float var2 = (1. - w1 * var1) / w2;
    for (int ic = 0; ic < 2; ic++) {
      // choose component and multiply with total variance
      float var = ic == 0 ? var1 : var2;
      var *= 225.e-6 * dp1 / p / p;
      AlgebraicSymMatrix55 cov;
      // transform from orthogonal planes containing the
      // momentum vector to local parameters
      float sl = pvec.perp();
      float cl = pvec.z();
      float cf = pvec.x() / sl;
      float sf = pvec.y() / sl;
      using namespace materialEffect;
      effects[ic].deltaCov[msxx] += var * (sf * sf * cl * cl + cf * cf) / (cl * cl * cl * cl);
      effects[ic].deltaCov[msxy] += var * (cf * sf * sl * sl) / (cl * cl * cl * cl);
      effects[ic].deltaCov[msyy] += var * (cf * cf * cl * cl + sf * sf) / (cl * cl * cl * cl);
    }
  }
}