EnergyLossUpdator.cc File Reference

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

Go to the source code of this file.

Functions
void	oldComputeBetheBloch (const LocalVector &localP, const MediumProperties &materialConstants, double mass)
void	oldComputeElectrons (const LocalVector &localP, const MediumProperties &materialConstants, const PropagationDirection propDir)

Function Documentation

void oldComputeBetheBloch	(	const LocalVector &	localP,
		const MediumProperties &	materialConstants,
		double	mass
	)

Definition at line 161 of file EnergyLossUpdator.cc.

References pfBoostedDoubleSVAK8TagInfos_cfi::beta, gather_cfg::cout, plot_hgcal_utils::dEdx, MillePedeFileConverter_cfg::e, CustomPhysics_cfi::gamma, cmsBatch::log, funct::m, PV3DBase< T, PVType, FrameType >::mag(), MaterialEffectsUpdator::mass(), AlCaHLTBitMon_ParallelJobs::p, particleFlowDisplacedVertex_cfi::ratio, mathSSE::sqrt(), MediumProperties::xi(), protons_cff::xi, and PV3DBase< T, PVType, FrameType >::z().

                                                                      {

 double theDeltaP=0., theDeltaCov=0;

  //
  // calculate absolute momentum and correction to path length from angle 
  // of incidence
  //
  double p = localP.mag();
  double xf = fabs(p/localP.z());
  // constants
  const double m = mass;            // use mass hypothesis from constructor

  const double emass = 0.511e-3;
  const double poti = 16.e-9 * 10.75; // = 16 eV * Z**0.9, for Si Z=14
  const double eplasma = 28.816e-9 * sqrt(2.33*0.498); // 28.816 eV * sqrt(rho*(Z/A)) for Si
  const double delta0 = 2*log(eplasma/poti) - 1.;

  // calculate general physics things
  double e     = sqrt(p*p + m*m);
  double beta  = p/e;
  double gamma = e/m;
  double eta2  = beta*gamma; eta2 *= eta2;
  //  double lnEta2 = log(eta2);
  double ratio = emass/m;
  double emax  = 2.*emass*eta2/(1. + 2.*ratio*gamma + ratio*ratio);
  // double delta = delta0 + lnEta2;

  // calculate the mean and sigma of energy loss
  // xi = d[g/cm2] * 0.307075MeV/(g/cm2) * Z/A * 1/2
  double xi = materialConstants.xi()*xf; xi /= (beta*beta);

//   double dEdx = xi*(log(2.*emass*eta2*emax/(poti*poti)) - 2.*(beta*beta));
  //double dEdx = xi*(log(2.*emass*emax/(poti*poti))+lnEta2 - 2.*(beta*beta) - delta);

  double dEdx = xi*(log(2.*emass*emax/(poti*poti)) - 2.*(beta*beta) - delta0);
  double dEdx2 = xi*emax*(1.-0.5*(beta*beta));
  double dP    = dEdx/beta;
  double sigp2 = dEdx2*e*e/(p*p*p*p*p*p);
  theDeltaP += -dP;
  theDeltaCov += sigp2;

  std::cout << "pion old " <<  theDeltaP << " " << theDeltaCov << std::endl;

}

void oldComputeElectrons	(	const LocalVector &	localP,
		const MediumProperties &	materialConstants,
		const PropagationDirection	propDir
	)

Definition at line 210 of file EnergyLossUpdator.cc.

References gather_cfg::cout, JetChargeProducer_cfi::exp, f, cmsBatch::log, PV3DBase< T, PVType, FrameType >::mag(), oppositeToMomentum, AlCaHLTBitMon_ParallelJobs::p, MediumProperties::radLen(), z, and PV3DBase< T, PVType, FrameType >::z().

                                                 {

  double theDeltaP=0., theDeltaCov=0;

  //
  // calculate absolute momentum and correction to path length from angle 
  // of incidence
  //
  double p = localP.mag();
  double normalisedPath = fabs(p/localP.z())*materialConstants.radLen();
  //
  // Energy loss and variance according to Bethe and Heitler, see also
  // Comp. Phys. Comm. 79 (1994) 157. 
  //
  double z = exp(-normalisedPath);
  double varz = exp(-normalisedPath*log(3.)/log(2.))- 
                z*z;
        // exp(-2*normalisedPath);

  if ( propDir==oppositeToMomentum ) {
    //
    // for backward propagation: delta(1/p) is linear in z=p_outside/p_inside
    // convert to obtain equivalent delta(p). Sign of deltaP is corrected
    // in method compute -> deltaP<0 at this place!!!
    //
    theDeltaP += -p*(1/z-1);
    theDeltaCov += varz/(p*p);
  }
  else {    
    //
    // for forward propagation: calculate in p (linear in 1/z=p_inside/p_outside),
    // then convert sig(p) to sig(1/p). 
    //
    theDeltaP += p*(z-1);
    //    double f = 1/p/z/z;
    // patch to ensure consistency between for- and backward propagation
    double f = 1./(p*z);
    theDeltaCov += f*f*varz;
  }

  std::cout << "electron old " <<  theDeltaP << " " << theDeltaCov << std::endl;


}