/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/TrackingTools/MaterialEffects/src/MultipleScatteringUpdator.cc File Reference

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

Go to the source code of this file.

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

---

Function Documentation

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

Definition at line 108 of file MultipleScatteringUpdator.cc.

References a, beta, gather_cfg::cout, ExpressReco_HICollisions_FallBack::e, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localMomentum(), funct::log(), m, PV3DBase< T, PVType, FrameType >::mag(), LocalTrajectoryError::matrix(), Surface::mediumProperties(), L1TEmulatorMonitor_cff::p, PV3DBase< T, PVType, FrameType >::perp(), ExpressReco_HICollisions_FallBack::pMin, MediumProperties::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()) {
    // Momentum vector
    LocalVector d = TSoS.localMomentum();
    double p = d.mag();
    d *= 1./p;
    // MediumProperties mp(0.02, .5e-4);
    const MediumProperties& mp = *surface.mediumProperties();
    double xf = 1./fabs(d.z());         // increase of path due to angle of incidence
    // calculate general physics things
    const double amscon = 1.8496e-4;    // (13.6MeV)**2
    const double m = mass;            // use mass hypothesis from constructor
    double e     = sqrt(p*p + m*m);
    double beta  = p/e;
    // calculate the multiple scattering angle
    double radLen = mp.radLen()*xf;     // effective rad. length
    double sigt2 = 0.;                  // sigma(alpha)**2
    if (radLen > 0) {
      // Calculated rms scattering angle squared.
      double 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();
        double 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.
        double pMin = thePtMin*(TSoS.globalMomentum().mag()/TSoS.globalMomentum().perp());       
        // Use P constraint to calculate rms maximum scattering angle.
        double betaMin = pMin/sqrt(pMin * pMin + m*m);
        double a_max = (1. + 0.038*log(radLen))/(betaMin * pMin); a_max *= a_max;
        double 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
      }
    }
    double sl = d.perp();
    double cl = d.z();
    double cf = d.x()/sl;
    double 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;
  }
}