//____________________________________________________________________________
/*!

\class   PropaBrem

\brief   A class to compute muon energy loss: bremsstrahlung 

\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 "PropaBrem.h"

using namespace std;
using namespace propamuon;

double B(double Z){
  double B;

  if(Z==1)B=202.4;
  else if(Z==2)B=151.9;
  else if(Z==3)B=159.9;
  else if(Z==4)B=172.3;
  else if(Z==5)B=177.9;
  else if(Z==6)B=178.3;
  else if(Z==7)B=176.6;
  else if(Z==8)B=173.4;  
  else if(Z==9)B=170.0;
  else if(Z==10)B=165.8;
  else if(Z==11)B=165.8;
  else if(Z==12)B=167.1;
  else if(Z==13)B=169.1;
  else if(Z==14)B=170.8;
  else if(Z==15)B=172.2;
  else if(Z==16)B=173.4;
  else if(Z==17)B=174.3;
  else if(Z==18)B=174.8;
  else if(Z==19)B=175.1;
  else if(Z==20)B=175.6;
  else if(Z==21)B=176.2;
  else if(Z==22)B=176.8;
  else if(Z==26)B=175.8;
  else if(Z==29)B=173.1;
  else if(Z==32)B=173.0;
  else if(Z==35)B=173.5;
  else if(Z==42)B=175.9;
  else if(Z==50)B=177.4;
  else if(Z==53)B=178.6;
  else if(Z==74)B=177.6;
  else if(Z==82)B=178.0;
  else if(Z==92)B=179.8;
  else B=182.7;

  return B;
}


double BremXSecDiff(double Energy, double Z, double A, double v, int Model){
  // return diff xsec in cm^2
  if (v<0) return 0;
  if (v>1) return 0;
  if (Energy<0) return 0;
  
  double ds_dv;
// Models:
// 1   Petrukhin and Shestakov
// 2   Kelner, Kokoulin and Petrukhin
  
  if(Model==1){   
    // Calculate the minimum monentum transfer to the nucleus (in GeV)
    double delta = (kMuonMass2/Energy) * 0.5*v/(1.-v);
    // Calculate the fi(delta) according to the Petrukhin/Shestakov formula 
    double a  = ((Z<10) ? 189.*kMuonElectronMass*TMath::Power(Z,-1./3.) : 189.*kMuonElectronMass*(2./3.)*TMath::Power(Z,-2./3.));
    double b  = 1.+(189./kElectronMass)*TMath::Power(Z,-1./3.)*kSqrte*delta;
    double fi = TMath::Log(a/b);
    // Calculate the Bethe-Heitler cross section ds/dv for muon bremsstrahlung in cm^2  
    ds_dv  = kAem3*(4*Z*(Z+1))*kLe2/kMuonElectronMass2*(4/3.-4*v/3.+v*v)/v*fi;
    
//    ds_dv  = kAem/v*TMath::Power((2.*Z*kElectronMass/kMuonMass*kRe),2)*(4./3.-4./3*v+v*v)*fi;
    
  }
  else if(Model==2){
    // Calculate the minimum monentum transfer to the nucleus (in GeV)
    
    
    double delta = (kMuonMass2/Energy) * 0.5*v/(1.-v);    
    double a=B(Z)*kMuonMass*TMath::Power(Z,-1./3.)/kElectronMass;
    double b=1.+delta*kSqrte*B(Z)*TMath::Power(Z,-1./3.)/kElectronMass;

    double Dn=1.54*TMath::Power(A,0.27);    
    double Delta_n= TMath::Log(Dn/(1.+delta*(Dn*kSqrte-2.)/kMuonMass));
    
    
    double fi=TMath::Log(a/b)-Delta_n;

    ds_dv  = kAem3*(4*Z*(Z+1))*kLe2/kMuonElectronMass2*(4/3.-4*v/3.+v*v)/v*fi;  
   
//    ds_dv  = kAem/v*TMath::Power((2.*Z*kElectronMass/kMuonMass*kRe),2)*(4./3.-4./3*v+v*v)*fi;
    
  }
  
  
  return ds_dv;
}
//____________________________________________________________________________
PropaBrem::PropaBrem() { 
  fModel=1;
}
//____________________________________________________________________________
PropaBrem::PropaBrem(int Model) { 
  fModel=Model;
}
//____________________________________________________________________________
PropaBrem::~PropaBrem(){ }


//____________________________________________________________________________
double PropaBrem::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 PropaBrem::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 PropaBremBCoeff(Energy,Z,A,fModel);
  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 PropaBrem::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 PropaBrem::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 PropaBrem::dE_dx(double Energy, double Z, double A) {
  return this->BCoeff(Energy,Z,A)*Energy;
}
//____________________________________________________________________________
double PropaBrem::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 PropaBrem::dE_dx(double Energy, map<int,double> MediaComposition) {
  return this->BCoeff(Energy,MediaComposition)*Energy;
}
//____________________________________________________________________________
double PropaBrem::dE_dx(double Energy, map<int,double> MediaComposition, double Vmin, double Vmax) {
  return this->BCoeff(Energy,MediaComposition,Vmin,Vmax)*Energy;
}
//____________________________________________________________________________
double PropaBrem::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 PropaBrem::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 PropaBrem::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 PropaBrem::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 PropaBremXSec(Energy,Z,A,fModel);
  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 PropaBrem::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 PropaBrem::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 PropaBrem::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 PropaBrem::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<<" "<<Percentage*this->SigmaTot(Energy,Z,A,Vmin,Vmax)/A<<" "<<lambda<<endl;
  }
  lambda*=kNA;
  lambda=1./lambda;  
  return lambda;
}
//____________________________________________________________________________
PropaBremBCoeff::PropaBremBCoeff(double Energy, double Z, double A, int Model) : ROOT::Math::IBaseFunctionOneDim() {
  fModel = Model;
  fEnergy = Energy;
  fZ = Z;
  fA = A;
}
//____________________________________________________________________________
PropaBremBCoeff::~PropaBremBCoeff() { }
//____________________________________________________________________________
unsigned int PropaBremBCoeff::NDim(void) const {
  return 1;
}
//____________________________________________________________________________
double PropaBremBCoeff::DoEval(double xin) const {
  double vds_dv = xin*BremXSecDiff(fEnergy,fZ,fA,xin,fModel);
  return vds_dv; 
}
//____________________________________________________________________________
ROOT::Math::IBaseFunctionOneDim * PropaBremBCoeff::Clone(void) const { 
  return new PropaBremBCoeff(fEnergy, fZ, fA, fModel);
}
//____________________________________________________________________________
PropaBremXSec::PropaBremXSec(double Energy, double Z, double A, int Model) : ROOT::Math::IBaseFunctionOneDim() {
  fModel = Model;
  fEnergy = Energy;
  fZ = Z;
  fA = A;
}
//____________________________________________________________________________
PropaBremXSec::~PropaBremXSec() { }
//____________________________________________________________________________
unsigned int PropaBremXSec::NDim(void) const {
  return 1;
}
//____________________________________________________________________________
double PropaBremXSec::DoEval(double xin) const {
  double ds_dv=BremXSecDiff(fEnergy,fZ,fA,xin,fModel);
  return ds_dv; 
}
//____________________________________________________________________________
ROOT::Math::IBaseFunctionOneDim * PropaBremXSec::Clone(void) const { 
  return new PropaBremXSec(fEnergy, fZ, fA, fModel);
}
//____________________________________________________________________________