//____________________________________________________________________________ /*! \class PropaPair \brief A class to compute muon energy loss: pair production \author Carla Distefano LNS-INFN, Catania \created December 9, 2015 \cpright Copyright (c) 2015-2019, The KM3NeT Collaboration For the full text of the license see $GSEAGEN/LICENSE */ //____________________________________________________________________________ #include #include #include #include #include "PropaUtils.h" #include "PropaPair.h" using namespace std; using namespace propamuon; double PairXSecDiff(double Energy, double Z, double v, double p){ // return diff xsec in cm^2 double vmin = 4.*kElectronMass/Energy; double vmax = 1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1./3.); if (vvmax) return 0; if (Energy<0) return 0; double pmax = (1.-6.*kMuonMass2/(Energy*Energy*(1.-v)))*TMath::Sqrt(1.-4.*kElectronMass/(Energy*v)); if(p>pmax) return 0; if(p<-pmax) return 0.; double p2=p*p; double beta=v*v/2./(1.-v); double Xi=kMuonMass2/kElectronMass2*v*v/4.*(1.-p2)/(1.-v); double R=189.; double Ye1=5.-p2+4.*beta*(1.+p2); double Ye2=2.*(1.+3.*beta)*TMath::Log(3.+1./Xi)-p2-2.*beta*(2.-p2); double Ye=Ye1/Ye2; double Ym1=4.+p2+3.*beta*(1.+p2); double Ym2=(1.+p2)*(3./2.+2.*beta)*TMath::Log(3.+Xi)+1.-3./2.*p2; double Ym=Ym1/Ym2; double LE1a=R*TMath::Power(Z,-1./3.)*TMath::Sqrt((1.+Xi)*(1.+Ye)); double LE1b=1.+(2.*kElectronMass*kSqrte*R*TMath::Power(Z,-1./3.)*(1+Xi)*(1.+Ye))/(Energy*v*(1-p2)); double LE1= TMath::Log((LE1a)/(LE1b)); // double LE2=-0.5*TMath::Log(1.+9./4.*kElectronMass2/kMuonMass2*TMath::Power(Z,2./3.)*(1.+Xi)*(1.+Ye)); double LE2=-0.5*TMath::Log(1.+9./4./kMuonElectronMass2*TMath::Power(Z,2./3.)*(1.+Xi)*(1.+Ye)); double LE=LE1+LE2; double LMa=2./3.*kMuonMass/kElectronMass*R*TMath::Power(Z,-2./3.); // double LMa=2./3.*kMuonElectronMass*R*TMath::Power(Z,-2./3.); double LMb=1.+(2.*kElectronMass*kSqrte*R*TMath::Power(Z,-1./3.)*(1+Xi)*(1.+Ym))/(Energy*v*(1-p2)); double LM=TMath::Log(LMa/LMb);; double PhiE1=((2.+p2)*(1.+beta)+Xi*(3.+p2))*TMath::Log(1.+1./Xi); double PhiE2=(1.-p2-beta)/(1.+Xi); double PhiE3=-(3.+p2); double PhiE=(PhiE1+PhiE2+PhiE3)*LE; double PhiM1=((1+p2)*(1.+3/2.*beta)-1./Xi*(1.+2.*beta)*(1.-p2))*TMath::Log(1.+Xi); ; double PhiM2=Xi*(1.-p2-beta)/(1.+Xi); double PhiM3=(1.+2.*beta)*(1.-p2); double PhiM=(PhiM1+PhiM2+PhiM3)*LM; double d2s_dvdp=kAem4*2./(3.*kPi)*Z*(Z+1)*kLe2*(1-v)/v*(PhiE+PhiM/kMuonElectronMass2); return d2s_dvdp; } //____________________________________________________________________________ PropaPair::PropaPair() { } //____________________________________________________________________________ PropaPair::~PropaPair(){ } //____________________________________________________________________________ double PropaPair::BCoeff (double Energy, double Z, double A){ double Vmin = 4.*kElectronMass/Energy; double Vmax = 1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1./3.); return this->BCoeff(Energy, Z, A, Vmin, Vmax); } //____________________________________________________________________________ double PropaPair::BCoeff (double Energy, double Z, double A, double Vmin, double Vmax){ Vmin=TMath::Max(Vmin,4.*kElectronMass/Energy); Vmax=TMath::Min(Vmax,1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1./3.)); ROOT::Math::IBaseFunctionOneDim * integrand = new PropaPairBCoeff(Energy,Z); ROOT::Math::IntegrationOneDim::Type ig_type = ROOT::Math::IntegrationOneDim::kADAPTIVE; double abstol = 1; double reltol = 1E-4; int nmaxeval = 100000; ROOT::Math::Integrator ig(*integrand,ig_type,abstol,reltol,nmaxeval); double b = kNA/A*ig.Integral(Vmin, Vmax); // in 1/(g/cm^2) delete integrand; return b; } //____________________________________________________________________________ double PropaPair::dE_dx(double Energy, double Z, double A) { return this->BCoeff(Energy,Z,A)*Energy; } //____________________________________________________________________________ double PropaPair::dE_dx(double Energy, double Z, double A, double Vmin, double Vmax) { Vmin=TMath::Max(Vmin,4.*kElectronMass/Energy); Vmax=TMath::Min(Vmax,1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1./3.)); return this->BCoeff(Energy,Z,A,Vmin,Vmax)*Energy; } //____________________________________________________________________________ double PropaPair::BCoeff (double Energy, map MediaComposition){ double Vmin = 4.*kElectronMass/Energy; double b=0.; map::iterator iter = MediaComposition.begin(); for ( ; iter != MediaComposition.end(); ++iter) { int TgtPdgCode=iter->first; double Percentage=iter->second; double Z=propamuon::PdgToZ(TgtPdgCode); double A=propamuon::PdgToA(TgtPdgCode); double Vmax = 1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.); b += Percentage* this->BCoeff(Energy,Z,A,Vmin,Vmax); } return b; } //____________________________________________________________________________ double PropaPair::BCoeff (double Energy, map MediaComposition, double Vmin, double Vmax){ Vmin=TMath::Max(Vmin,4.*kElectronMass/Energy); double b=0.; map::iterator iter = MediaComposition.begin(); for ( ; iter != MediaComposition.end(); ++iter) { int TgtPdgCode=iter->first; double Percentage=iter->second; double Z=propamuon::PdgToZ(TgtPdgCode); double A=propamuon::PdgToA(TgtPdgCode); Vmax=TMath::Min(Vmax,1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.)); b += Percentage* this->BCoeff(Energy,Z,A,Vmin,Vmax); } return b; } //____________________________________________________________________________ double PropaPair::dE_dx(double Energy, map MediaComposition) { return this->BCoeff(Energy,MediaComposition)*Energy; } //____________________________________________________________________________ double PropaPair::dE_dx(double Energy, map MediaComposition, double Vmin, double Vmax) { return this->BCoeff(Energy,MediaComposition,Vmin,Vmax)*Energy; } //____________________________________________________________________________ double PropaPair::SigmaTot (double Energy, double Z){ double Vmin = 4.*kElectronMass/Energy; double Vmax = 1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.); return this->SigmaTot(Energy, Z, Vmin, Vmax); } //____________________________________________________________________________ double PropaPair::SigmaTot(double Energy, double Z, double Vmin, double Vmax){ Vmin=TMath::Max(Vmin,4.*kElectronMass/Energy); Vmax=TMath::Min(Vmax,1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1./3.)); double Pmin=-1.; double Pmax=1.; double Xmin[2] = { Vmin, Pmin }; double Xmax[2] = { Vmax, Pmax }; ROOT::Math::IBaseFunctionMultiDim * integrand = new PropaPairXSec2Dim(Energy,Z); ROOT::Math::IntegrationMultiDim::Type ig_type = ROOT::Math::IntegrationMultiDim::kADAPTIVE; double abstol = 1; double reltol = 1E-4; int nmaxeval = 500000; ROOT::Math::IntegratorMultiDim ig(*integrand,ig_type,abstol,reltol,nmaxeval); double sig = ig.Integral(Xmin, Xmax); // in cm^2 delete integrand; return sig; } //____________________________________________________________________________ double PropaPair::SigmaTot(double Energy, map MediaComposition) { double Vmin = 4.*kElectronMass/Energy; double sigma=0.; map::iterator iter = MediaComposition.begin(); for ( ; iter != MediaComposition.end(); ++iter) { int TgtPdgCode=iter->first; double Percentage=iter->second; double Z=propamuon::PdgToZ(TgtPdgCode); double Vmax = 1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.); sigma += Percentage* this->SigmaTot(Energy,Z,Vmin,Vmax); } return sigma; } //____________________________________________________________________________ double PropaPair::SigmaTot(double Energy, map MediaComposition, double Vmin, double Vmax) { Vmin=TMath::Max(Vmin,4.*kElectronMass/Energy); double sigma=0.; map::iterator iter = MediaComposition.begin(); for ( ; iter != MediaComposition.end(); ++iter) { int TgtPdgCode=iter->first; double Percentage=iter->second; double Z=propamuon::PdgToZ(TgtPdgCode); Vmax=TMath::Min(Vmax,1.-0.75*kSqrte*(kMuonMass/Energy)*TMath::Power(Z,1/3.)); sigma += Percentage*this->SigmaTot(Energy,Z,Vmin,Vmax); } return sigma; } //____________________________________________________________________________ double PropaPair::Lambda(double Energy, double Z, double A) { double Vmin = 4.*kElectronMass/Energy; double Vmax = 1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.); double lambda=this->Lambda(Energy,Z,A,Vmin,Vmax); return lambda; } //____________________________________________________________________________ double PropaPair::Lambda(double Energy, double Z, double A, double Vmin, double Vmax) { Vmin=TMath::Max(Vmin,4.*kElectronMass/Energy); Vmax=TMath::Min(Vmax,1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.)); double lambda=this->SigmaTot(Energy,Z,Vmin,Vmax); lambda*=kNA/A; lambda=1./lambda; return lambda; } //____________________________________________________________________________ double PropaPair::Lambda(double Energy, map MediaComposition) { double Vmin = 4.*kElectronMass/Energy; double lambda=0.; map::iterator iter = MediaComposition.begin(); for ( ; iter != MediaComposition.end(); ++iter) { int TgtPdgCode=iter->first; double Percentage=iter->second; double Z=propamuon::PdgToZ(TgtPdgCode); double A=propamuon::PdgToA(TgtPdgCode); double Vmax = 1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.); lambda += Percentage* this->SigmaTot(Energy,Z,Vmin,Vmax)/A; // lambda += Percentage* this->SigmaTot(Energy,Z,Vmin,Vmax); } lambda*=kNA; lambda=1./lambda; return lambda; } //____________________________________________________________________________ double PropaPair::Lambda(double Energy, map MediaComposition, double Vmin, double Vmax) { Vmin=TMath::Max(Vmin,4.*kElectronMass/Energy); double lambda=0.; map::iterator iter = MediaComposition.begin(); for ( ; iter != MediaComposition.end(); ++iter) { int TgtPdgCode=iter->first; double Percentage=iter->second; double Z=propamuon::PdgToZ(TgtPdgCode); double A=propamuon::PdgToA(TgtPdgCode); Vmax=TMath::Min(Vmax,1.-0.75*kSqrte*(kMuonMass/Energy)*TMath::Power(Z,1/3.)); lambda += Percentage*this->SigmaTot(Energy,Z,Vmin,Vmax)/A; } lambda*=kNA; lambda=1./lambda; return lambda; } //____________________________________________________________________________ PropaPairBCoeff::PropaPairBCoeff(double Energy, double Z) : ROOT::Math::IBaseFunctionOneDim() { fEnergy = Energy; fZ = Z; } //____________________________________________________________________________ PropaPairBCoeff::~PropaPairBCoeff() { } //____________________________________________________________________________ unsigned int PropaPairBCoeff::NDim(void) const { return 1; } //____________________________________________________________________________ double PropaPairBCoeff::DoEval(const double xin) const { // Returns v*(ds/dv) double v=xin; double pmax_v = (1.-6.*kMuonMass2/(fEnergy*fEnergy*(1.-v)))*TMath::Sqrt(1.-4.*kElectronMass/(fEnergy*v)); double pmin_v = 0.; ROOT::Math::IBaseFunctionOneDim * integrand = new PropaPairXSecV(fEnergy,fZ,v); ROOT::Math::IntegrationOneDim::Type ig_type = ROOT::Math::IntegrationOneDim::kADAPTIVE; double abstol = 1; double reltol = 1E-4; int nmaxeval = 100000; ROOT::Math::Integrator ig(*integrand,ig_type,abstol,reltol,nmaxeval); double vds_dv = 2.*v*ig.Integral(pmin_v, pmax_v); // in cm2 delete integrand; return vds_dv; } //____________________________________________________________________________ ROOT::Math::IBaseFunctionOneDim * PropaPairBCoeff::Clone(void) const { return new PropaPairBCoeff(fEnergy, fZ); } //____________________________________________________________________________ PropaPairXSec::PropaPairXSec(double Energy, double Z) : ROOT::Math::IBaseFunctionOneDim() { fEnergy = Energy; fZ = Z; } //____________________________________________________________________________ PropaPairXSec::~PropaPairXSec() { } //____________________________________________________________________________ unsigned int PropaPairXSec::NDim(void) const { return 1; } //____________________________________________________________________________ double PropaPairXSec::DoEval(const double xin) const { // Returns (ds/dv) double v=xin; double pmax_v = (1.-6.*kMuonMass2/(fEnergy*fEnergy*(1.-v)))*TMath::Sqrt(1.-4.*kElectronMass/(fEnergy*v)); double pmin_v = 0.; ROOT::Math::IBaseFunctionOneDim * integrand = new PropaPairXSecV(fEnergy,fZ,v); ROOT::Math::IntegrationOneDim::Type ig_type = ROOT::Math::IntegrationOneDim::kADAPTIVE; double abstol = 1; double reltol = 1E-4; int nmaxeval = 100000; ROOT::Math::Integrator ig(*integrand,ig_type,abstol,reltol,nmaxeval); double ds_dv = 2.*ig.Integral(pmin_v, pmax_v); // in cm2 delete integrand; return ds_dv; } //____________________________________________________________________________ ROOT::Math::IBaseFunctionOneDim * PropaPairXSec::Clone(void) const { return new PropaPairXSec(fEnergy, fZ); } //____________________________________________________________________________ PropaPairXSecV::PropaPairXSecV(double Energy, double Z, double v) : ROOT::Math::IBaseFunctionOneDim() { fEnergy = Energy; fZ = Z; fV = v; } //____________________________________________________________________________ PropaPairXSecV::~PropaPairXSecV() { } //____________________________________________________________________________ unsigned int PropaPairXSecV::NDim(void) const { return 1; } //____________________________________________________________________________ double PropaPairXSecV::DoEval(double xin) const { double ds_dv=PairXSecDiff(fEnergy,fZ,fV, xin); return ds_dv; } //____________________________________________________________________________ ROOT::Math::IBaseFunctionOneDim * PropaPairXSecV::Clone(void) const { return new PropaPairXSecV(fEnergy, fZ, fV); } //____________________________________________________________________________ PropaPairXSec2Dim::PropaPairXSec2Dim(double Energy, double Z) : ROOT::Math::IBaseFunctionMultiDim() { fEnergy = Energy; fZ = Z; } //____________________________________________________________________________ PropaPairXSec2Dim::~PropaPairXSec2Dim() { } //____________________________________________________________________________ unsigned int PropaPairXSec2Dim::NDim(void) const { return 2; } //____________________________________________________________________________ double PropaPairXSec2Dim::DoEval(const double * xin) const { // Returns (ds/dv) double v=xin[0]; double p=xin[1]; double d2s_dvdp=PairXSecDiff(fEnergy, fZ, v, p); return d2s_dvdp; } //____________________________________________________________________________ ROOT::Math::IBaseFunctionMultiDim * PropaPairXSec2Dim::Clone(void) const { return new PropaPairXSec2Dim(fEnergy, fZ); }