#include <KokoulinPetrukhinModel.h>

Inheritance diagram for genie::mueloss::gsl::KokoulinPetrukhinIntegrand:

Collaboration diagram for genie::mueloss::gsl::KokoulinPetrukhinIntegrand:

Public Member Functions
	KokoulinPetrukhinIntegrand (double E, double Z)

	~KokoulinPetrukhinIntegrand ()

unsigned int	NDim (void) const

double	DoEval (const double *xin) const

ROOT::Math::IBaseFunctionMultiDim *	Clone (void) const

Private Attributes
double	fE

double	fZ

Detailed Description

Definition at line 80 of file KokoulinPetrukhinModel.h.

Constructor & Destructor Documentation

gsl::KokoulinPetrukhinIntegrand::KokoulinPetrukhinIntegrand	(	double	E,
		double	Z
	)

Definition at line 93 of file KokoulinPetrukhinModel.cxx.

References fE, and fZ.

                                                                             :
 ROOT::Math::IBaseFunctionMultiDim()
 {
   fE = E;
   fZ = Z;
 }

gsl::KokoulinPetrukhinIntegrand::~KokoulinPetrukhinIntegrand ( )

Definition at line 100 of file KokoulinPetrukhinModel.cxx.

 {
 
 }

Member Function Documentation

ROOT::Math::IBaseFunctionMultiDim * gsl::KokoulinPetrukhinIntegrand::Clone ( void ) const

Definition at line 181 of file KokoulinPetrukhinModel.cxx.

 {
   return new gsl::KokoulinPetrukhinIntegrand(fE,fZ);
 }

double gsl::KokoulinPetrukhinIntegrand::DoEval ( const double * xin ) const

Definition at line 110 of file KokoulinPetrukhinModel.cxx.

References genie::units::b, genie::constants::kAem, genie::constants::kElectronMass, genie::constants::kElectronMass2, genie::constants::kLe, genie::constants::kMuonMass, genie::constants::kMuonMass2, genie::constants::kPi, and genie::constants::kSqrtNapierConst.

 {
 // Returns v*(d^2s/dvdp)
 
   double v  = xin[0]; // v, the fraction of energy transfered to the photon
   double p  = xin[1]; //
 
   if (! (v >0)) return 0;
   if (   v >1)  return 0;
   if (! (fE>0)) return 0;
 
   double pmax_v = (1. - 6.*kMuonMass2 / (fE*fE*(1.-v)) ) *
                                  TMath::Sqrt(1.-4.*kElectronMass/(fE*v));
   if(p>pmax_v) return 0;
 
   double v2      = TMath::Power(v,2.);
   double p2      = TMath::Power(p,2.);
   double R       = 189.;
   double a4      = TMath::Power(kAem,4.);
   double pi      = kPi;
   double ZLe2    = TMath::Power(fZ*kLe,2.);
   double Zm13    = TMath::Power(fZ,-1./3.);
   double Zm23    = TMath::Power(fZ,-2./3.);
   double Z13     = TMath::Power(fZ,1./3.);
   double me      = kElectronMass;
   double me2     = kElectronMass2;
   double mmu     = kMuonMass;
   double mmu2    = kMuonMass2;
   double memu2   = me2/mmu2;
   double memu    = me/mmu;
   double mume    = mmu/me;
 
   double b    = 0.5*v2/(1.-v);
   double xi   = (1.-p2) * TMath::Power(0.5*v*mume, 2.) / (1.-v);
 
   // Approximate the FIm factor (dimensionless) for the Kokoulin Petrukhin
   // pair production cross section
 
   double FImA = ( (1.+p2)*(1.+3.*b/2.) - (1.+2.*b)*(1.-p2)/xi ) * TMath::Log(1.+xi);
   double FImB = xi*(1.-p2-b) / (1.+xi);
   double FImC = (1.+2.*b) * (1.-p2);
   double Ym   = (4. + p2 + 3.*b*(1.+p2)) /
                       ((1.+p2)*(1.5+2.*b)*TMath::Log(3.+xi) + 1. - 1.5*p2);
   double LmA  = (2./3.) * mume * R * Zm23;
   double LmB  = 1. + (2.*me*R * kSqrtNapierConst * Zm13 * (1+xi) * (1+Ym)) / (fE*v*(1-p2) );
   double Lm   = TMath::Log(LmA/LmB);
   double FIm  = (FImA+FImB+FImC)*Lm;
 
   // Approximate the FIe factor (dimensionless) for the Kokoulin Petrukhin
   // pair production cross section
 
   double FIeA = ( (2.+p2) * (1.+b) + xi*(3.+p2) ) * log(1.+1./xi);
   double FIeB = (1.-p2-b) / (1.+xi);
   double FIeC = 3. + p2;
   double Ye   = (5. - p2 + 4*b*(1+p2)) /
                     (2.*(1.+3.*b)*TMath::Log(3.+1./xi) - p2 - 2.*b*(2.-p2));
   double x_Y  = (1+xi)*(1+Ye);
   double LeA  = R*Zm13*TMath::Sqrt(x_Y);
   double LeB  = 1. + (2.*me*R*kSqrtNapierConst*Zm13*x_Y) / (fE*v*(1-p2));
   double LeC  = 1.5 * memu * Z13;
   double LeD  = 1 + TMath::Power(LeC,2.)*x_Y;
   double Le   = TMath::Log(LeA/LeB) - 0.5*TMath::Log(LeD);
   double FIe  = (FIeA+FIeB-FIeC)*Le;
 
   double d2s_dvdp = a4 * (2./(3.*pi)) * ZLe2 * ((1.-v)/v) * (FIe + FIm*memu2);
 
   double vd2s_dvdp = v*d2s_dvdp;
   return vd2s_dvdp;
 }