BezrukovBugaevModel.cxx

Go to the documentation of this file.

#include <cassert>

#include "BezrukovBugaevModel.h"
#include "Integrator.h"
#include "Constants.h"

ClassImp(BezrukovBugaevModel)

//_________________________________________________________________________________________
BezrukovBugaevModel::BezrukovBugaevModel() :
ProcessModel()
{

}
//_________________________________________________________________________________________
BezrukovBugaevModel::BezrukovBugaevModel(const Material & material) :
ProcessModel(material, Process::eNuclearInteraction)
{

}
//_________________________________________________________________________________________
BezrukovBugaevModel::~BezrukovBugaevModel()
{

}
//_________________________________________________________________________________________
ValidityRange_t BezrukovBugaevModel::ValidityRange(void) const
{
  ValidityRange_t validity_range;

  validity_range.Emin = 0.75*sqrt(MuELoss::e)*(MuELoss::Mm)*pow(fMaterial.Z(), 1./3.);
  validity_range.Emax = 10000000.; // MeV

  return validity_range;
}
//_________________________________________________________________________________________
double BezrukovBugaevModel::dE_dx(double E) const
{
// Calculate the muon -dE/dx due to muon nuclear interaction in MeV / (gr / cm^2 )
// Input  : E, muon energy in MeV

  return bnuc(E) * E;
}
//_________________________________________________________________________________________
double BezrukovBugaevModel::bnuc(double E) const
{
// Calculate the b factor for muon nuclear interaction ( -dE/dx = b*E ) in cm^2 / gr
// Input  : E, muon energy in MeV

  const int kNV = 600;

  double Vmin = 0;
  double Vmax = V_max(E);
  double dV   = (Vmax - Vmin) / kNV;

  Integrator * integrator = Integrator::Instance();

  double * vds = new double[kNV];

  for(int iv = 0; iv < kNV; iv++)  vds[iv] = (Vmin+(iv+1)*dV) * ds_dv(Vmin+(iv+1)*dV, E);

  double integral = integrator->Simpson(vds, kNV, dV);

  delete [] vds;

  // 1e-7 a dimensional factor:
  // v*ds/dv (in microbarns) --> v*ds/dv*10^-30 cm^2, N -> N*10^23

  return 1e-7*(MuELoss::Na/fMaterial.A())*integral;
}
//_________________________________________________________________________________________
double BezrukovBugaevModel::V_max(double E) const
{
// Calculate the maximum fraction of energy carried to the photon.
// Input  : E, muon energy in MeV
//
  return ( 1 - (3/4.)*sqrt(MuELoss::e)*(MuELoss::Mm/E)* pow(fMaterial.Z(),1/3.) );
}
//_________________________________________________________________________________________
double BezrukovBugaevModel::ds_dv(double v, double E) const
{
// Calculate the differential cross section ds/dv (in microbarn) for muon nuclear
// interaction based on the Bezrukov, Bugaev formula.
// Input  : E, muon energy in MeV
//          v, fraction of energy transferred to the photon
//
  assert(E >= 0);
  assert(v >  0);

  double sig = PhotonuclXSec(v,E);
  double g   = 3.*G(v,E)/4.;
  double v2  = v*v;
  double t   = MuELoss::Mm_2 *v2/(1-v);
  double k   = 1. - 2./v + 2./v2;

  assert(t >  0);

  double m1_2_t = MuELoss::m1_2 / t;
  double m2_2_t = MuELoss::m2_2 / t;
  double mm_2_t = MuELoss::Mm_2 / t;
  double d      = MuELoss::m1_2 / ( t + MuELoss::m1_2 );

  double A = MuELoss::a * fMaterial.A() * sig * v / (2 * MuELoss::pi);

  double B = g          * ( k*log(1.+m1_2_t) - k*d - 2.*mm_2_t );
  double C = 0.25       * ( k*log(1.+m2_2_t) - 2.*mm_2_t );
  double D = 0.5*mm_2_t * ( g*d + 0.25*m2_2_t*log(1.+1./m2_2_t) );

  return A*(B+C+D); // microbarns
}
//___________________________________________________________________________________________
double BezrukovBugaevModel::G(double v, double E) const
{
// Calculate the factor G (dimensionless) appearing in the differential photonuclear
// interaction xsec
// Input  : E, muon energy in MeV
//          v, fraction of energy transferred to the photon
//
  assert(E >= 0);
  assert(v >  0);
  double A13 = pow(fMaterial.A(),1/3.);
  double x   = 0.00282 * A13 * PhotonuclXSec(v,E);  // factor 0.00282 has units of ub^-1
  double x2  = x*x;
  double x3  = x2*x;

  return 3. * ( x2/2. - 1. + (1.+x)*exp(-x) ) / x3;
}
//___________________________________________________________________________________________
double BezrukovBugaevModel::PhotonuclXSec(double v, double E) const
{
// Calculate the cross section (in microbarns) for photonuclear interaction
// Input  : E, muon energy in MeV
//          v, fraction of energy transferred to the photon
//
  assert(E >= 0);
  assert(v >  0);

  double Ep   = v*E; // photon energy
  double loge = log(0.0213*1e-3*Ep); // factor 0.0213 has units of GeV^-1
  
  return 114.3 + 1.647 * loge*loge;
}
//___________________________________________________________________________________________

---