//____________________________________________________________________________
/*!

\class   PropaNucl

\brief   A class to compute muon energy loss: photo-nuclear 

\author  Carla Distefano <distefano_c@lns.infn.it>
         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 <TMath.h>
#include <Math/Integrator.h>

#include "PropaUtils.h"
#include "PropaNucl.h"

using namespace std;
using namespace propamuon;

double NuclXSecDiff(double Energy, double A, double v){
  // return diff xsec in cm^2

  if (v<0) return 0;
  if (v>1) return 0;
  if (Energy<0) return 0;

  double t      = kMuonMass2 *v*v/(1-v);
  double k      = 1.-2./v+2./(v*v);
  double A13    = TMath::Power(A,1./3.);
  double M1_2   = 0.54; // m1^2 in photonuclear diff. xsec formula (in GeV^2)
  double M2_2   = 1.80; // m2^2 in photonuclear diff. xsec formula (in GeV^2)
  double M1_2_t = M1_2/t;
  double M2_2_t = M2_2/t;
  double mmu2_t = kMuonMass2/t;
  double d      = M1_2/(t+M1_2);

  // Calculate the cross section (in ub) for photonuclear interaction
  double Ep   = v*Energy; // photon energy (GeV)
  double loge = TMath::Log(0.0213*Ep); // factor 0.0213 has units of GeV^-1
  double sig  = 114.3+1.647*loge*loge; // in ub

  // Calculate the factor G (dimensionless) appearing in the differential
  // photonuclear interaction cross section
  double x   = 0.00282*A13*sig;  // factor 0.00282 has units of ub^-1
  double x2  = x*x;
  double x3  = x2*x;
  double G   = 3*(0.5*x2-1.+(1.+x)*TMath::Exp(-x))/x3;

  // Calculate the differential cross section ds/dv for muon nuclear
  // interaction based on the Bezrukov-Bugaev formula.
  double bbA = 0.5*(kAem/kPi) * A * sig * v;
  double bbB = 0.75*G     * ( k*TMath::Log(1.+M1_2_t) - k*d - 2.*mmu2_t );
  double bbC = 0.25       * ( k*TMath::Log(1.+M2_2_t) - 2.*mmu2_t );
  double bbD = 0.5*mmu2_t * ( 0.75*G*d + 0.25*M2_2_t*TMath::Log(1.+1./M2_2_t) );

  double ds_dv  = bbA*(bbB+bbC+bbD)*1E-30;    //units::ub/units::cm2; // in cm2

  return ds_dv;
}
//____________________________________________________________________________
PropaNucl::PropaNucl() { }
//____________________________________________________________________________
PropaNucl::~PropaNucl(){ }
//____________________________________________________________________________
double PropaNucl::BCoeff (double Energy, double Z, double A){
  double Vmin = 0.;
  double Vmax = 1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.);
  return this->BCoeff(Energy, Z, A, Vmin, Vmax);
}
//____________________________________________________________________________
double PropaNucl::BCoeff (double Energy, double Z, double A, double Vmin, double Vmax){

  Vmin=TMath::Max(Vmin,0.);
  Vmax=TMath::Min(Vmax,1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.));

  ROOT::Math::IBaseFunctionOneDim * integrand = new PropaNuclBCoeff(Energy,A);
  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);	// b in 1/(g/cm^2)
  delete integrand;  
  return b;
}
//____________________________________________________________________________
double PropaNucl::dE_dx(double Energy, double Z, double A) {
  return this->BCoeff(Energy,Z,A)*Energy;
}


//____________________________________________________________________________
double PropaNucl::dE_dx(double Energy, double Z, double A, double Vmin, double Vmax) {

  Vmin=TMath::Max(Vmin,0.);
  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 PropaNucl::BCoeff (double Energy, map<int,double> MediaComposition){
  double Vmin = 0.;
  double b=0.;
  map<int,double>::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 PropaNucl::BCoeff (double Energy, map<int,double> MediaComposition, double Vmin, double Vmax){

  Vmin=TMath::Max(Vmin,0.);
  double b=0.;
  map<int,double>::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 PropaNucl::dE_dx(double Energy, map<int,double> MediaComposition) {
  return this->BCoeff(Energy,MediaComposition)*Energy;

}
//____________________________________________________________________________
double PropaNucl::dE_dx(double Energy, map<int,double> MediaComposition, double Vmin, double Vmax) {

  return this->BCoeff(Energy,MediaComposition,Vmin,Vmax)*Energy;

}
//____________________________________________________________________________
double PropaNucl::SigmaTot (double Energy, double Z, double A){
  double Vmin = 0.;
  double Vmax = 1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.);
  return this->SigmaTot(Energy, Z, A, Vmin, Vmax);
}
//____________________________________________________________________________
double PropaNucl::SigmaTot (double Energy, double Z, double A, double Vmin, double Vmax){

  Vmin=TMath::Max(Vmin,0.);
  Vmax=TMath::Min(Vmax,1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.));

  ROOT::Math::IBaseFunctionOneDim * integrand = new PropaNuclXSec(Energy,A);
  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 sig = ig.Integral(Vmin, Vmax);  // in cm2
  delete integrand;  
  return sig;
}


//____________________________________________________________________________
double PropaNucl::SigmaTot (double Energy, map<int,double> MediaComposition){
  double Vmin = 0.;

  double sigma=0.;  

  map<int,double>::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.);
    sigma += Percentage* this->SigmaTot(Energy, Z, A, Vmin, Vmax);
  }
  
  return sigma;
}

//____________________________________________________________________________
double PropaNucl::SigmaTot (double Energy, map<int,double> MediaComposition, double Vmin, double Vmax){

  Vmin=TMath::Max(Vmin,0.);
  double sigma=0.;  

  map<int,double>::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.));
    sigma += Percentage* this->SigmaTot(Energy, Z, A, Vmin, Vmax);
  }
  return sigma;
}




//____________________________________________________________________________
double PropaNucl::Lambda(double Energy, double Z, double A) {
  double Vmin = 0.;
  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 PropaNucl::Lambda(double Energy, double Z, double A, double Vmin, double Vmax) {

  Vmin=TMath::Max(Vmin,0.);
  Vmax=TMath::Min(Vmax,1. - 0.75*kSqrte* (kMuonMass/Energy) * TMath::Power(Z,1/3.));

  double lambda=this->SigmaTot(Energy,Z,A,Vmin,Vmax);
  lambda*=kNA/A;
  lambda=1./lambda;  
  return lambda;
}
//____________________________________________________________________________
double PropaNucl::Lambda(double Energy, map<int,double> MediaComposition) {
  double Vmin = 0.;
  double lambda=0.;
  map<int,double>::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,A,Vmin,Vmax)/A;
  }
  lambda*=kNA;
  lambda=1./lambda;  
  return lambda;
}
//____________________________________________________________________________
double PropaNucl::Lambda(double Energy, map<int,double> MediaComposition, double Vmin, double Vmax) {

  Vmin=TMath::Max(Vmin,0.);


  double lambda=0.;
  map<int,double>::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,A,Vmin,Vmax)/A;
//    lambda += Percentage* this->SigmaTot(Energy,Z,A,Vmin,Vmax);    
//    cout<<Z<<" "<<A<<" "<<Percentage<<" "<<this->SigmaTot(Energy,Z,A,Vmin,Vmax)<<" "<<lambda<<" "<<kNA<<endl;
  }
  lambda*=kNA;
  lambda=1./lambda;  
  return lambda;
}
//____________________________________________________________________________
PropaNuclBCoeff::PropaNuclBCoeff(double Energy, double A) : ROOT::Math::IBaseFunctionOneDim() {
  fEnergy = Energy;
  fA = A;
}
//____________________________________________________________________________
PropaNuclBCoeff::~PropaNuclBCoeff() { }
//____________________________________________________________________________
unsigned int PropaNuclBCoeff::NDim(void) const {
  return 1;
}
//____________________________________________________________________________
double PropaNuclBCoeff::DoEval(double xin) const {
// Calculate the Bethe-Heitler cross section ds/dv for muon bremsstrahlung
// Returns v*(ds/dv)
  double vds_dv = xin*NuclXSecDiff(fEnergy,fA,xin);
  return vds_dv; 
}
//____________________________________________________________________________
ROOT::Math::IBaseFunctionOneDim * PropaNuclBCoeff::Clone(void) const { 
  return new PropaNuclBCoeff(fEnergy, fA);
}
//____________________________________________________________________________
PropaNuclXSec::PropaNuclXSec(double Energy, double A) : ROOT::Math::IBaseFunctionOneDim() {
  fEnergy = Energy;
  fA = A;
}
//____________________________________________________________________________
PropaNuclXSec::~PropaNuclXSec() { }
//____________________________________________________________________________
unsigned int PropaNuclXSec::NDim(void) const {
  return 1;
}
//____________________________________________________________________________
double PropaNuclXSec::DoEval(double xin) const {
// Calculate the Bethe-Heitler cross section ds/dv for muon bremsstrahlung
// Returns (ds/dv)
  double ds_dv=NuclXSecDiff(fEnergy,fA,xin);
  return ds_dv; 
}
//____________________________________________________________________________
ROOT::Math::IBaseFunctionOneDim * PropaNuclXSec::Clone(void) const { 
  return new PropaNuclXSec(fEnergy, fA);
}
//____________________________________________________________________________