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

oldComputeBetheBloch()

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

Definition at line 143 of file EnergyLossUpdator.cc.

References HLT_2023v12_cff::beta, gather_cfg::cout, plot_hgcal_utils::dEdx, MillePedeFileConverter_cfg::e, HLT_2023v12_cff::eta2, CustomPhysics_cfi::gamma, dqm-mbProfile::log, visualization-live-secondInstance_cfg::m, PV3DBase< T, PVType, FrameType >::mag(), MaterialEffectsUpdator::mass(), AlCaHLTBitMon_ParallelJobs::p, particleFlowDisplacedVertex_cfi::ratio, mathSSE::sqrt(), protons_cff::xi, MediumProperties::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;
}

oldComputeElectrons()

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

Definition at line 189 of file EnergyLossUpdator.cc.

References gather_cfg::cout, JetChargeProducer_cfi::exp, f, dqm-mbProfile::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;
}