//____________________________________________________________________________
/*!

\class   PropaIon

\brief   A class to compute muon energy loss: ionization 

\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 <cassert>
#include <cstdlib>
#include <string>
#include <vector>
#include <cstdio>
#include <iostream>
#include <stdio.h>

#include <TMath.h>

#include "PropaUtils.h"
#include "PropaIon.h"

using namespace std;
using namespace propamuon;

//____________________________________________________________________________
PropaIon::PropaIon() { }
//____________________________________________________________________________
PropaIon::~PropaIon(){ }
//____________________________________________________________________________
double PropaIon::ACoeff (double Energy, double Z, double A){
  
  double Energy2 = TMath::Power(Energy,2);
  double beta    = TMath::Sqrt(Energy2-kMuonMass2)/Energy;
  double beta2   = TMath::Power(beta,2);
  double gamma   = Energy/kMuonMass;
  double gamma2  = TMath::Power(gamma,2);

  double p2      = Energy2-kMuonMass2;
  double EMaxT   = 2.*kElectronMass*p2 / (kElectronMass2 + kMuonMass2 + 2.*kElectronMass*Energy);
  double EMaxT2  = TMath::Power(EMaxT,2);

  SterCoeff coeff = this->CorrFactors(Z);
  double I2      = TMath::Power(coeff.I,2); // in GeV^2
  double X  =  TMath::Log10(beta*gamma);

  double delta = 0;
  if(coeff.X0<X && X<coeff.X1) delta = 4.6052*X + coeff.a*TMath::Power(coeff.X1-X,coeff.m) + coeff.C;
  if(X>coeff.X1)	       delta = 4.6052*X + coeff.C;

  double acoeff = kAem2*(2*kPi*kNA*kLe2)*Z/A*(kElectronMass/beta2)*(TMath::Log(2.*kElectronMass*beta2*gamma2*EMaxT/I2)-2.*beta2+0.25*(EMaxT2/Energy2)-delta);
  return acoeff;
}
//____________________________________________________________________________
double PropaIon::dE_dx(double Energy, double Z, double A) {
  return this->ACoeff(Energy,Z,A);
}
//____________________________________________________________________________
double PropaIon::ACoeff (double Energy, map<int,double> MediaComposition){
  double a=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);
    a += Percentage*this->ACoeff(Energy,Z,A);      
  }
  return a;
}
//____________________________________________________________________________
double PropaIon::dE_dx(double Energy, map<int,double> MediaComposition) {
  return this->ACoeff(Energy,MediaComposition);
}
//___________________________________________________________________________
SterCoeff PropaIon::CorrFactors(int Z){
//  returns the Sternheimer coeff.
// The ionization potential convereted in GeV before exiting

  SterCoeff coeff;

  switch(Z){
    case  1:
      coeff.name = "H";
      coeff.I    = 21.8; // in eV
      coeff.X0   = 0.4400; coeff.X1 = 1.8856; coeff.a = 0.1348; coeff.m = 5.6249; coeff.C = -3.0977; 
      break;      
    case  6:
      coeff.name = "C";
      coeff.I    = 78.0; // in eV
      coeff.X0   = -0.0090; coeff.X1 = 2.4817; coeff.a = 0.2076; coeff.m = 2.9532; coeff.C = -2.8926;
      break;      
    case  7:
      coeff.name = "N";
      coeff.I    = 82.0; // in eV
      coeff.X0   = 0.3039; coeff.X1 = 2.0000; coeff.a = 0.5329; coeff.m = 3.0000; coeff.C = -3.9996;
      break;      
    case  8:
      coeff.name = "O";
      coeff.I    = 95.0; // in eV
      coeff.X0   = 0.2868; coeff.X1 = 2.0000; coeff.a = 0.5223; coeff.m = 3.0000; coeff.C = -3.9471;
      break;      
    case  11:
      coeff.name = "Na";
      coeff.I    = 149.0; // in eV
      coeff.X0   = 0.2880; coeff.X1 = 3.1962; coeff.a = 0.0777; coeff.m = 3.6452; coeff.C = -5.0526;
      break;      
    case  12:
      coeff.name = "Mg";
      coeff.I    = 156.0; // in eV
      coeff.X0   = 0.1499; coeff.X1 = 3.0668; coeff.a = 0.0816; coeff.m = 3.6166; coeff.C = -4.5297;
      break;      
    case  13:
      coeff.name = "Al";
      coeff.I    = 166.0; // in eV
      coeff.X0   = 0.1708; coeff.X1 = 3.0127; coeff.a = 0.0802; coeff.m = 3.6345; coeff.C = -4.2395;
      break;      
    case  14:
      coeff.name = "C";
      coeff.I    = 173.0; // in eV
      coeff.X0   = 0.2015; coeff.X1 = 2.8716; coeff.a = 0.1492; coeff.m = 3.2546; coeff.C = -4.4355;
      break;      
    case  16:
      coeff.name = "S";
      coeff.I    = 180.0; // in eV
      coeff.X0   = 0.1580; coeff.X1 = 2.7159; coeff.a = 0.3399; coeff.m = 2.6456; coeff.C = -4.6659;
      break;      
    case  17:
      coeff.name = "Cl";
      coeff.I    = 174.0; // in eV
      coeff.X0   = 0.2000; coeff.X1 = 3.0000; coeff.a = 0.1802; coeff.m = 3.0000; coeff.C = -4.8776;
      break;      
    case  19:
      coeff.name = "K";
      coeff.I    = 190.0; // in eV
      coeff.X0   = 0.3851; coeff.X1 = 3.1724; coeff.a = 0.1983; coeff.m = 2.9233; coeff.C = -5.6423;
      break;      
    case  20:
      coeff.name = "Ca";
      coeff.I    = 191.0; // in eV
      coeff.X0   = 0.3228; coeff.X1 = 3.1191; coeff.a = 0.1564; coeff.m = 3.0745; coeff.C = -5.0396;
      break;      
    case  22:
      coeff.name = "Ti";
      coeff.I    = 233.0; // in eV
      coeff.X0   = 0.0957; coeff.X1 = 3.0386; coeff.a = 0.1566; coeff.m = 3.0302; coeff.C = -4.4450;
      break;      
    case  26:
      coeff.name = "Fe";
      coeff.I    = 286.0; // in eV
      coeff.X0   = -0.0012; coeff.X1 = 3.1531; coeff.a = 0.1468; coeff.m = 2.9632; coeff.C = -4.2911;
      break;      
    case  28:
      coeff.name = "Ni";
      coeff.I    = 311.0; // in eV
      coeff.X0   = -0.0566; coeff.X1 = 3.1851; coeff.a = 0.1650; coeff.m = 2.8430; coeff.C = -4.3115;
      break;      
    case  35:
      coeff.name = "Br";
      coeff.I    = 357.0; // in eV
      coeff.X0   = 0.3669; coeff.X1 = 3.0000; coeff.a = 0.2211; coeff.m = 3.0000; coeff.C = -5.7268;
      break;      
    default:
      cout<<"PropaIon: No Sternheimer coeff found for element with Z="<<Z<<" ...exiting.";
      exit(1);
  }  

  coeff.I/=1.E9; // in GeV

  return coeff;
}