GFEnergyLossBetheBloch.cxx

Go to the documentation of this file.

/* Copyright 2008-2009, Technische Universitaet Muenchen,
   Authors: Christian Hoeppner & Sebastian Neubert

   This file is part of GENFIT.

   GENFIT is free software: you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as published
   by the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   GENFIT is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public License
   along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "larreco/Genfit/GFEnergyLossBetheBloch.h"
#include "larreco/Genfit/GFException.h"
#include "math.h"

genf::GFEnergyLossBetheBloch::~GFEnergyLossBetheBloch(){
}

double genf::GFEnergyLossBetheBloch::energyLoss(const double& step,
                                          const double& mom,
                                          const int&    pdg,              
                                          const double& matDensity,
                                          const double& matZ,
                                          const double& matA,
                                          const double& /* radiationLength */,
                                          const double& meanExcitationEnergy,
                                          const bool&   doNoise,
                                                TMatrixT<Double_t>* noise,
                                          const TMatrixT<Double_t>* /* jacobian */,
                                          const TVector3* /* directionBefore */,
                                          const TVector3* /* directionAfter */){

  static const double me = getParticleMass(11); // electron mass (GeV)  0.000519...

  double charge, mass;
  getParticleParameters(pdg, charge, mass);
  
  const double beta = mom/sqrt(mass*mass+mom*mom);
  double dedx = 0.307075*matZ/matA*matDensity/(beta*beta)*charge*charge;

  //for numerical stability
  double gammaSquare = 1.-beta*beta;
  if(gammaSquare>1.E-10) gammaSquare = 1./gammaSquare;
  else gammaSquare = 1.E10;
  double gamma = sqrt(gammaSquare);

  double massRatio = me/mass;
  double argument = gammaSquare*beta*beta*me*1.E3*2./((1.E-6*meanExcitationEnergy) * sqrt(1+2*sqrt(gammaSquare)*massRatio + massRatio*massRatio));
  if (argument <= exp(beta*beta)) 
    dedx = 0.;
  else{
    dedx *= (log(argument)-beta*beta); // Bethe-Bloch [MeV/cm]
    dedx *= 1.E-3;  // in GeV/cm, hence 1.e-3
    if (dedx <= 0.)
      throw GFException("GFEnergyLossBetheBloch::energyLoss(): non-positive dE/dx", __LINE__, __FILE__).setFatal();
  }
  
  double DE = step * dedx; //always positive
  double momLoss = sqrt(mom*mom+2.*sqrt(mom*mom+mass*mass)*DE+DE*DE) - mom; //always positive
  
  //in vacuum it can numerically happen that momLoss becomes a small negative number. A cut-off at 0.01 eV for momentum loss seems reasonable
  if(fabs(momLoss)<1.E-11)momLoss=1.E-11;

    
  if (doNoise) {
    if (!noise)
      throw GFException("GFEnergyLossBetheBloch::energyLoss(): no noise matrix specified!", __LINE__, __FILE__).setFatal();
    
    // ENERGY LOSS FLUCTUATIONS; calculate sigma^2(E); 
    double sigma2E = 0.;
    double zeta  = 153.4E3 * charge*charge/(beta*beta) * matZ/matA * matDensity * step; // eV
    double Emax  = 2.E9*me*beta*beta*gammaSquare / (1. + 2.*gamma*me/mass + (me/mass)*(me/mass) ); // eV
    double kappa = zeta/Emax;
    
    if (kappa > 0.01) { // Vavilov-Gaussian regime
      sigma2E += zeta*Emax*(1.-beta*beta/2.);  // eV^2
    }
    else { // Urban/Landau approximation
      double alpha = 0.996;
      double sigmaalpha = 15.76;
      // calculate number of collisions Nc
      double I = 16. * pow(matZ, 0.9); // eV
      double f2 = 0.;
      if (matZ > 2.) f2 = 2./matZ;
      double f1 = 1. - f2;
      double e2 = 10.*matZ*matZ; // eV
      double e1 = pow( (I/pow(e2,f2)), 1./f1);  // eV
      
      double mbbgg2 = 2.E9*mass*beta*beta*gammaSquare; // eV
      double Sigma1 = dedx*1.0E9 * f1/e1 * (log(mbbgg2 / e1) - beta*beta) / (log(mbbgg2 / I) - beta*beta) * 0.6; // 1/cm
      double Sigma2 = dedx*1.0E9 * f2/e2 * (log(mbbgg2 / e2) - beta*beta) / (log(mbbgg2 / I) - beta*beta) * 0.6; // 1/cm
      double Sigma3 = dedx*1.0E9 * Emax / ( I*(Emax+I)*log((Emax+I)/I) ) * 0.4; // 1/cm
        
      double Nc = (Sigma1 + Sigma2 + Sigma3)*step;
      
      if (Nc > 50.) { // truncated Landau distribution 
        // calculate sigmaalpha  (see GEANT3 manual W5013)
        double RLAMED = -0.422784 - beta*beta - log(zeta/Emax);
        double RLAMAX =  0.60715 + 1.1934*RLAMED +(0.67794 + 0.052382*RLAMED)*exp(0.94753+0.74442*RLAMED);
        // from lambda max to sigmaalpha=sigma (empirical polynomial)
        if(RLAMAX <= 1010.) {
           sigmaalpha =  1.975560  
                        +9.898841e-02 *RLAMAX  
                        -2.828670e-04 *RLAMAX*RLAMAX
                        +5.345406e-07 *pow(RLAMAX,3.)   
                        -4.942035e-10 *pow(RLAMAX,4.) 
                        +1.729807e-13 *pow(RLAMAX,5.); 
        }
        else { sigmaalpha = 1.871887E+01 + 1.296254E-02 *RLAMAX; }
        // alpha=54.6  corresponds to a 0.9996 maximum cut
        if(sigmaalpha > 54.6) sigmaalpha=54.6;
        sigma2E += sigmaalpha*sigmaalpha * zeta*zeta;  // eV^2
      }
      else { // Urban model
        double Ealpha  = I / (1.-(alpha*Emax/(Emax+I)));   // eV
        double meanE32 = I*(Emax+I)/Emax * (Ealpha-I);     // eV^2
        sigma2E += step * (Sigma1*e1*e1 + Sigma2*e2*e2 + Sigma3*meanE32); // eV^2
      }
    }
        
    sigma2E*=1.E-18; // eV -> GeV
    
    // update noise matrix
    (*noise)[6][6] += (mom*mom+mass*mass)/pow(mom,6.)*sigma2E; 
  }
  
  return momLoss;
}

//ClassImp(GFEnergyLossBetheBloch)