KokoulinPetrukhinModel.cxx

Go to the documentation of this file.

#include <cassert>

#include "KokoulinPetrukhinModel.h"
#include "Constants.h"
#include "Integrator.h"

ClassImp(KokoulinPetrukhinModel)

//_________________________________________________________________________________________
KokoulinPetrukhinModel::KokoulinPetrukhinModel() :
ProcessModel()
{

}
//_________________________________________________________________________________________
KokoulinPetrukhinModel::KokoulinPetrukhinModel(const Material & material) :
ProcessModel(material, Process::ePairProduction)
{

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

}
//_________________________________________________________________________________________
ValidityRange_t KokoulinPetrukhinModel::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 KokoulinPetrukhinModel::dE_dx(double E) const
{
// Calculates the muon -dE/dx in ( MeV cm^2 / gr ) due to e+e- pair production
// Input  : E, muon energy in MeV

  return bpair(E) * E;
}
//_________________________________________________________________________________________
double KokoulinPetrukhinModel::bpair(double E) const
{
// Purpose: get the b factor for e+e- pair production from muons ( -dE/dx = b*E )
// Input  : E, muon energy in MeV
// Output : b factor in cm^2 / gr

  const int kNV  = 300;
  const int kNP  = 300;

  assert(E>=0);

  Integrator * integrator = Integrator::Instance();

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

  double * vds_dv = new double[kNV];

  for(int iv = 0; iv < kNV; iv++) {

     double v = Vmin + (iv+1) * dV;

     double Pmin = 0;
     double Pmax = P_max(v, E);

     if(Pmax < 0) return 0; // put this here for now so as not to hit an assert() below... 
                            // Eventually I need to solve Pmax(E) > 0 for E, and combine it
                            // with the E limit from Vmax in the validity range of the model

     double dP   = (Pmax-Pmin)/kNP;

     double * d2s = new double[kNP];

     for(int ip = 0; ip < kNP; ip++)  d2s[ip] = d2s_dvdp(v, (Pmin + (ip+1) * dP) ,E);

     vds_dv[iv] = v * integrator->Simpson(d2s, kNP, dP);

     delete [] d2s;

  } // v

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

  delete [] vds_dv;

  // factor 10 is for dimensions:
  //    10 :   v*ds/dv(integral) in 10^22 cm2 --  N->N*10^(-23) -> b in 0.1 cm2/gr

  return 10* ( 2*MuELoss::Na / fMaterial.A() ) * integral;
}
//___________________________________________________________________________________________
double KokoulinPetrukhinModel::d2s_dvdp(double v, double p, double E) const
{
// Calculate the differential cross section d^2s/dpdv (in 10^22 cm^2) for e+e- pair production
// by muons based on the Kokoulin,Petrukhin parameterization
// Input  : v, the fraction of the energy transfered to the e+e- pair
//          p, the asymmetry parameter of the e+e- pair p = (E(+) - E(-)) / (E(+) + E(-))
//          E, muon energy in MeV
//
  assert(v>=0);
  assert(p>=0);
  assert(E>=0);

  double a4      = MuELoss::a_4;
  double pi      = MuELoss::pi;
  double ZLe2    = pow(fMaterial.Z()*MuELoss::Le,2.);
  double me_mm_2 = MuELoss::Me_Mm_2;
  double fie     = FIe(v,p,E);
  double fim     = FIm(v,p,E);

  return a4 * (2./(3.*pi)) * ZLe2 * ( (1.-v)/v ) * ( fie + fim * me_mm_2 );
}
//___________________________________________________________________________________________
double KokoulinPetrukhinModel::FIe(double v, double p, double E) const
{
// Approximate the FIe factor (dimensionless) from pair production cross section
// Input  : v, the fraction of the energy transfered to the e+e- pair
//          p, the asymmetry parameter of the e+e- pair p = (E(+) - E(-)) / (E(+) + E(-))
//          E, muon energy in MeV
//
  double b  = v*v / (2*(1-v));
  double p2 = p*p;
  double xi = (1.-p2) * pow( 0.5 *v *MuELoss::Mm_Me , 2.) / (1-v);
  double A  = ( (2.+p2) * (1.+b) + xi*(3.+p2) ) * log(1.+1./xi);
  double B  = (1.-p2-b) / (1.+xi);
  double C  = 3. + p2;

  return (A+B-C)*Le(v,p,E);
}
//___________________________________________________________________________________________
double KokoulinPetrukhinModel::FIm(double v, double p, double E) const
{
// Approximates the factor FIm factor (dimensionless) from pair production cross section
// Input  : v, the fraction of the energy transfered to the e+e- pair
//          p, the asymmetry parameter of the e+e- pair p = (E(+) - E(-)) / (E(+) + E(-))
//          E, muon energy in MeV
//
  double b  = v*v / (2*(1-v));
  double p2 = p*p;
  double xi = (1.-p2) * pow( 0.5 *v *MuELoss::Mm_Me , 2.) / (1-v);
  double A  = ( (1.+p2)*(1.+3.*b/2.) - (1.+2.*b)*(1.-p2)/xi ) * log(1.+xi);
  double B  = xi * (1. - p2 - b) / (1.+xi);
  double C  = (1.+2.*b) * (1.-p2);

  return (A+B+C)*Lm(v,p,E);
}
//___________________________________________________________________________________________
double KokoulinPetrukhinModel::Le(double v, double p, double E) const
{
// Approximates the Le factor (dimensionless) appearing in the factor FIe at pair production xsec
// Input  : v, the fraction of the energy transfered to the e+e- pair
//          p, the asymmetry parameter of the e+e- pair p = (E(+) - E(-)) / (E(+) + E(-))
//          E, muon energy in MeV
//
  double p2   = p*p;
  double xi   = (1.-p2) * pow( 0.5 *v *MuELoss::Mm_Me , 2.) / (1-v);
  double Y    = Ye(v,p);
  double Z_13 = pow(fMaterial.Z(),-1./3.);
  double Z13  = pow(fMaterial.Z(), 1./3.);
  double MeR  = MuELoss::Me * MuELoss::R;
  double x_Y  = (1+xi)*(1+Y);

  double A    = MuELoss::R * Z_13 * sqrt( x_Y );
  double B    = 1. + (2. * MeR * MuELoss::sqrt_e * Z_13 * x_Y) / ( E*v*(1-p2) );
  double C    = 1.5 * MuELoss::Me_Mm * Z13;
  double D    = 1 + C*C*x_Y;

  return log(A/B) - 0.5 * log(D);
}
//___________________________________________________________________________________________
double KokoulinPetrukhinModel::Lm(double v, double p, double E) const
{
// Approximates the Lm factor (dimensioness) appearing in the factor FIm at pair production xsec
// Input  : v, the fraction of the energy transfered to the e+e- pair
//          p, the asymmetry parameter of the e+e- pair p = (E(+) - E(-)) / (E(+) + E(-))
//          E, muon energy in MeV

  double p2   = p*p;
  double xi   = (1.-p2) * pow( 0.5 *v *MuELoss::Mm_Me , 2.) / (1-v);
  double Y    = Ym(v,p);
  double Z_13 = pow(fMaterial.Z(),-1/3.);
  double Z_23 = pow(fMaterial.Z(),-2/3.);
  double MeR  = MuELoss::Me * MuELoss::R;
  double A    = (2./3.) * MuELoss::Mm_Me * MuELoss::R * Z_23;
  double B    = 1. + ( 2. * MeR * MuELoss::sqrt_e * Z_13 * (1+xi) * (1+Y) ) / ( E*v*(1-p2) );

  return log(A/B);
}
//___________________________________________________________________________________________
double KokoulinPetrukhinModel::Ye(double v, double p) const
{
// Approximates the Ye factor (dimensionless) in the factor Le at pair production xsec
// Input  : v, the fraction of the energy transfered to the e+e- pair
//          p, the asymmetry parameter of the e+e- pair p = (E(+) - E(-)) / (E(+) + E(-))
//
  double b  = v*v / (2*(1-v));
  double p2 = p*p;
  double xi   = (1.-p2) * pow( 0.5 *v *MuELoss::Mm_Me , 2.) / (1-v);

  return ( 5. - p2 + 4*b*(1+p2) ) / ( 2.*(1.+3.*b)*log(3.+1./xi) - p2 - 2.*b*(2.-p2) );
}
//___________________________________________________________________________________________
double KokoulinPetrukhinModel::Ym(double v, double p) const
{
// Approximates the Ym factor (dimensionless) in the factor Lm at pair production xsec
// Input  : v, the fraction of the energy transfered to the e+e- pair
//          p, the asymmetry parameter of the e+e- pair p = (E(+) - E(-)) / (E(+) + E(-))
//
  double b  = v*v / (2*(1-v));
  double p2 = p*p;
  double xi   = (1.-p2) * pow( 0.5 *v *MuELoss::Mm_Me , 2.) / (1-v);

  return ( 4. + p2 + 3.*b*(1.+p2) )  /  ( (1.+p2)*(1.5+2.*b)*log(3.+xi) + 1. - 1.5*p2 );
}
//___________________________________________________________________________________________
double KokoulinPetrukhinModel::V_min(double E) const
{
// Minimum value of parameter v = E(+)+(E-)/E

  return 4.*MuELoss::Me/E;
}
//___________________________________________________________________________________________
double KokoulinPetrukhinModel::V_max(double E) const
{
// Maximum value of parameter v = E(+)+(E-)/E

  return 1.- 3.* MuELoss::sqrt_e * MuELoss::Mm * pow(fMaterial.Z(), 1./3.) / (4.*E);
}
//___________________________________________________________________________________________
double KokoulinPetrukhinModel::P_max(double v, double E) const
{
  return (1. - 6.*MuELoss::Mm_2 / (E*E*(1.-v)) ) * sqrt(1.-4.*MuELoss::Me/(E*v));
}
//___________________________________________________________________________________________

---