EnergyLossUpdator.cc

Go to the documentation of this file.

#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"



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

//
// Computation of contribution of energy loss to momentum and covariance
// matrix of local parameters based on Bethe-Bloch. For electrons 
// contribution of radiation acc. to Bethe & Heitler.
//
void EnergyLossUpdator::compute (const TrajectoryStateOnSurface& TSoS, 
                 const PropagationDirection propDir, Effect & effect) const
{
  //
  // Get surface
  //
  const Surface& surface = TSoS.surface();
  //
  //
  // Now get information on medium
  //
  if (surface.mediumProperties().isValid()) {
    //
    // Bethe-Bloch
    //
    if ( mass()>0.001 )
      computeBetheBloch(TSoS.localMomentum(),surface.mediumProperties(),effect);
    //
    // Special treatment for electrons (currently rather crude
    // distinction using mass)
    //
    else
      computeElectrons(TSoS.localMomentum(),surface.mediumProperties(),
               propDir,effect);
    if (propDir != alongMomentum) effect.deltaP *= -1.;
  }
}
//
// Computation of energy loss according to Bethe-Bloch
//
void
EnergyLossUpdator::computeBetheBloch (const LocalVector& localP,
                      const MediumProperties& materialConstants, Effect & effect) const {
  //
  // calculate absolute momentum and correction to path length from angle 
  // of incidence
  //

  typedef float Float;

  Float p2 = localP.mag2();
  Float xf = std::abs(std::sqrt(p2)/localP.z());
   
  // constants
  const Float m2 = mass()*mass();           // use mass hypothesis from constructor

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

  // calculate general physics things
  Float im2 = Float(1.)/m2;
  Float e2     = p2 + m2;
  Float e = std::sqrt(e2);
  Float beta2  = p2/e2;
  Float eta2  = p2*im2;
  Float ratio2 = (emass*emass)*im2;
  Float emax  = Float(2.)*emass*eta2/(Float(1.) + Float(2.)*emass*e*im2 + ratio2);
  
  Float xi = materialConstants.xi()*xf; xi /= beta2;
  
  Float dEdx = xi*(unsafe_logf<2>(Float(2.)*emass*emax/(poti*poti)) - Float(2.)*(beta2) - delta0);
  
  Float dEdx2 = xi*emax*(Float(1.)-Float(0.5)*beta2);
  Float dP    = dEdx/std::sqrt(beta2);
  Float sigp2 = dEdx2/(beta2*p2*p2);
  effect.deltaP += -dP;
  using namespace materialEffect;
  effect.deltaCov[elos] += sigp2;


  // std::cout << "pion new " <<  theDeltaP << " " << theDeltaCov(0,0) << std::endl;
  // oldComputeBetheBloch (localP, materialConstants, mass());

}
//
// Computation of energy loss for electrons
//
void 
EnergyLossUpdator::computeElectrons (const LocalVector& localP,
                     const MediumProperties& materialConstants,
                     const PropagationDirection propDir, Effect & effect) const {
  //
  // calculate absolute momentum and correction to path length from angle 
  // of incidence
  //
  float p2 = localP.mag2();
  float p = std::sqrt(p2);
  float normalisedPath = std::abs(p/localP.z())*materialConstants.radLen();
  //
  // Energy loss and variance according to Bethe and Heitler, see also
  // Comp. Phys. Comm. 79 (1994) 157. 
  //
  const float l3ol2 = std::log(3.)/std::log(2.);
  float z = unsafe_expf<3>(-normalisedPath);
  float varz = unsafe_expf<3>(-normalisedPath*l3ol2)- 
                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!!!
    //
    effect.deltaP += -p*(1.f/z-1.f);
    using namespace materialEffect;
    effect.deltaCov[elos]  += varz/p2;
  }
  else {    
    //
    // for forward propagation: calculate in p (linear in 1/z=p_inside/p_outside),
    // then convert sig(p) to sig(1/p). 
    //
    effect.deltaP += p*(z-1.f);
    //    float f = 1/p/z/z;
    // patch to ensure consistency between for- and backward propagation
    float f2 = 1.f/(p2*z*z);
    using namespace materialEffect;
    effect.deltaCov[elos] += f2*varz;
  }


  // std::cout << "electron new " <<  theDeltaP << " " << theDeltaCov(0,0) << std::endl;
  // oldComputeElectrons (localP, materialConstants, propDir);

}









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

 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) {

  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;


}