MultipleScatteringUpdator.cc File Reference

#include "TrackingTools/MaterialEffects/interface/MultipleScatteringUpdator.h"
#include "DataFormats/GeometrySurface/interface/MediumProperties.h"
#include "DataFormats/Math/interface/approx_log.h"

Go to the source code of this file.

Functions
void	oldMUcompute (const TrajectoryStateOnSurface &TSoS, const PropagationDirection propDir, double mass, double ptmin)
void	oldMUcompute (const TrajectoryStateOnSurface &TSoS, const PropagationDirection propDir, float mass, float thePtMin)

Function Documentation

void oldMUcompute	(	const TrajectoryStateOnSurface &	TSoS,
		const PropagationDirection	propDir,
		double	mass,
		double	ptmin
	)

void oldMUcompute	(	const TrajectoryStateOnSurface &	TSoS,
		const PropagationDirection	propDir,
		float	mass,
		float	thePtMin
	)

Definition at line 101 of file MultipleScatteringUpdator.cc.

References a, zMuMuMuonUserData::beta, GetRecoTauVFromDQM_MC_cff::cl, gather_cfg::cout, ztail::d, MillePedeFileConverter_cfg::e, TrajectoryStateOnSurface::globalMomentum(), MediumProperties::isValid(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localMomentum(), dqm-mbProfile::log, visualization-live-secondInstance_cfg::m, PV3DBase< T, PVType, FrameType >::mag(), MaterialEffectsUpdator::mass(), LocalTrajectoryError::matrix(), Surface::mediumProperties(), AlCaHLTBitMon_ParallelJobs::p, PV3DBase< T, PVType, FrameType >::perp(), ALCARECOTkAlMinBias_cff::pMin, MediumProperties::radLen(), fastSimProducer_cff::radLen, mathSSE::sqrt(), TrajectoryStateOnSurface::surface(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                                                                        {
  //
  // Get surface
  //
  const Surface& surface = TSoS.surface();
  //
  //
  // Now get information on medium
  //
  if (surface.mediumProperties().isValid()) {
    // Momentum vector
    LocalVector d = TSoS.localMomentum();
    float p = d.mag();
    d *= 1. / p;
    // MediumProperties mp(0.02, .5e-4);
    const MediumProperties& mp = surface.mediumProperties();
    float xf = 1. / fabs(d.z());  // increase of path due to angle of incidence
    // calculate general physics things
    const float amscon = 1.8496e-4;  // (13.6MeV)**2
    const float m = mass;            // use mass hypothesis from constructor
    float e = sqrt(p * p + m * m);
    float beta = p / e;
    // calculate the multiple scattering angle
    float radLen = mp.radLen() * xf;  // effective rad. length
    float sigt2 = 0.;                 // sigma(alpha)**2
    if (radLen > 0) {
      // Calculated rms scattering angle squared.
      float a = (1. + 0.038 * log(radLen)) / (beta * p);
      a *= a;
      sigt2 = amscon * radLen * a;
      if (thePtMin > 0) {
#ifdef DBG_MSU
        std::cout << "Original rms scattering = " << sqrt(sigt2);
#endif
        // Inflate estimated rms scattering angle, to take into account
        // that 1/p is not known precisely.
        AlgebraicSymMatrix55 covMatrix = TSoS.localError().matrix();
        float error2_QoverP = covMatrix(0, 0);
        // Formula valid for ultra-relativistic particles.
        //      sigt2 *= (1. + (p*p) * error2_QoverP);
        // Exact formula
        sigt2 *= (1. + (p * p) * error2_QoverP * (1. + 5 * m * m / (e * e) + 3 * m * m * beta * beta * error2_QoverP));
#ifdef DBG_MSU
        std::cout << " new = " << sqrt(sigt2);
#endif
        // Convert Pt constraint to P constraint, neglecting uncertainty in
        // track angle.
        float pMin = thePtMin * (TSoS.globalMomentum().mag() / TSoS.globalMomentum().perp());
        // Use P constraint to calculate rms maximum scattering angle.
        float betaMin = pMin / sqrt(pMin * pMin + m * m);
        float a_max = (1. + 0.038 * log(radLen)) / (betaMin * pMin);
        a_max *= a_max;
        float sigt2_max = amscon * radLen * a_max;
        if (sigt2 > sigt2_max)
          sigt2 = sigt2_max;
#ifdef DBG_MSU
        std::cout << " after P constraint (" << pMin << ") = " << sqrt(sigt2);
        std::cout << " for track with 1/p=" << 1 / p << "+-" << sqrt(error2_QoverP) << std::endl;
#endif
      }
    }
    float sl = d.perp();
    float cl = d.z();
    float cf = d.x() / sl;
    float sf = d.y() / sl;
    // Create update (transformation of independant variations
    //   on angle in orthogonal planes to local parameters.
    std::cout << "old " << sigt2 * (sf * sf * cl * cl + cf * cf) / (cl * cl * cl * cl) << " "
              << sigt2 * (cf * sf * sl * sl) / (cl * cl * cl * cl) << " "
              << sigt2 * (cf * cf * cl * cl + sf * sf) / (cl * cl * cl * cl) << std::endl;
  }
}