#include <string>
#include <iostream>
#include <iomanip>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"

#include "JPhysics/JConstants.hh"
#include "JTools/JRange.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {
  
  using namespace JPP;
  
  static const double K =  0.307075;      // [MeV g^-1 cm^2]
  static const double Z = 10.0;           // atomic number [unit]
  static const double A = 18.0;           // atomic mass [g/mol]

  static const char* const electron = "electron";
  static const char* const positron = "positron";
  static const char* const muon     = "muon";
  static const char* const tau      = "tau";
}


/**
 * \file
 *
 * Program to determine the energy loss due to visible delta-rays.\n
 * Energy loss due to delta-rays and maximal kinetic energy (Tmax) are taken from reference
 * https://pdg.lbl.gov/2020/reviews/rpp2020-rev-passage-particles-matter.pdf
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;

  typedef JTOOLS::JRange<double> JRange_t;
  
  string   outputFile;
  JRange_t T_GeV;
  int      numberOfPoints;
  string   lepton;
  string   option;
  double   precision;
  int      debug;

  try {

    JParser<> zap("Program to determine the energy loss due to visible delta-rays.");

    zap['o'] = make_field(outputFile)     = "delta-rays.root";
    zap['T'] = make_field(T_GeV,          "kinetic energy range of electron [GeV]")  = JRange_t();
    zap['n'] = make_field(numberOfPoints, "points for integration")                  = 1000000;
    zap['L'] = make_field(lepton)         = muon, tau, positron, electron;
    zap['O'] = make_field(option)         = "";
    zap['e'] = make_field(precision)      = 1.0e-6;
    zap['d'] = make_field(debug)          = 1;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  using namespace JPP;

  double MASS_LEPTON;

  if      (lepton == muon)
    MASS_LEPTON = MASS_MUON;
  else if (lepton == tau)
    MASS_LEPTON = MASS_TAU;
  else if (lepton == positron)
    MASS_LEPTON = MASS_ELECTRON;
  else if (lepton == electron)
    MASS_LEPTON = MASS_ELECTRON;
  else
    FATAL("Invalid lepton " << lepton << endl);

  
  const double RATIO =  MASS_ELECTRON / MASS_LEPTON;

  
  // Minimal total   energy of electron for Cherenkov threshold [GeV]

  const double EMIN  =  MASS_ELECTRON * INDEX_OF_REFRACTION_WATER / sqrt(INDEX_OF_REFRACTION_WATER*INDEX_OF_REFRACTION_WATER - 1.0);

  // Minimal kinetic energy of electron for Cherenkov threshold [GeV]

  double Tmin  =  sqrt(EMIN*EMIN - MASS_ELECTRON*MASS_ELECTRON);
  
  NOTICE("Tmin [GeV] " << FIXED(8,6) << Tmin << ' ' << FIXED(8,6) << T_GeV.getLowerLimit() << endl);

  if (T_GeV.is_valid()) {
    if (Tmin < T_GeV.getLowerLimit()) {
      Tmin = T_GeV.getLowerLimit();
    }
  }

  
  TFile out(outputFile.c_str(), "recreate");

  TH1D h0("h0", NULL, 80, log10(MASS_LEPTON), +9.25);
  TH1D h1("h1", NULL, 80, log10(MASS_LEPTON), +9.25);
  
  for (int i = 1; i <= h0.GetNbinsX(); ++i) {

    const double x     = h0.GetBinCenter(i);
    const double E     = pow(10.0,x);         // particle energy [GeV]

    const double gamma = E / MASS_LEPTON;
    const double beta  = sqrt(1.0 - 1.0/(gamma*gamma));


    // Maximal kinetic energy of electron [GeV]
    
    double Tmax  = (lepton == electron ?
		    (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
		    (2.0*MASS_ELECTRON*beta*beta*gamma*gamma)  /  (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));

    if (T_GeV.is_valid()) {
      if (Tmax > T_GeV.getUpperLimit()) {
	Tmax = T_GeV.getUpperLimit();
      }
    }


    double weight = 0.0;
    double count  = 0.0;
      
    if (Tmax > Tmin) {
    
      const double a =  1.0;                  // Form factor formula
      const double b = -beta*beta/Tmax;       // 
      const double c =  0.5/(E*E);            // F(T)  =  a + b*T + c*T^2

      const double W =  0.5 * K * (Z/A) * (1.0/(beta*beta));
      
      {
	const double xmin = log(Tmin);        // T * 1/T^2 * dT = d log(T)
	const double xmax = log(Tmax);
	const double dx   = (xmax - xmin) / numberOfPoints;

	for (double x = xmin; x <= xmax; x += dx) {

	  const double T = exp(x);
	
	  const double F = a + T*(b + T*(c)); // Form factor formula
	  const double y = W * F * dx;

	  weight += y;
	}
      }
      {
	const double xmin = 1.0/Tmax;         // 1/T^2 * dT = d -1/T
	const double xmax = 1.0/Tmin;
	const double dx   = (xmax - xmin) / numberOfPoints;

	for (double x = xmin; x <= xmax; x += dx) {
	
	  const double T = 1.0/x;
	
	  const double F = a + T*(b + T*(c));   // Form factor formula
	  const double y = W * F * dx;

	  count += y;
	}
      }
    }
    
    DEBUG("E     [GeV]            " << SCIENTIFIC(8,2) << E      << endl);
    DEBUG("Tmin  [GeV]            " << SCIENTIFIC(8,2) << Tmin   << endl);
    DEBUG("Tmax  [GeV]            " << SCIENTIFIC(8,2) << Tmax   << endl);
    DEBUG("dE/dx [MeV g^-1 cm^2]  " << FIXED(5,4)      << weight << endl);
    DEBUG("dE/dx [g^-1 cm^2]      " << FIXED(5,2)      << count  << endl);

    h0.SetBinContent(i, weight);
    h1.SetBinContent(i, count);
  }


  TF1 f1("f1", "[0] + [1]*x + [2]*x*x + [3]*x*x*x");

  if (option.find('W') == string::npos) {
    option += "W";
  }
  
  f1.SetParameter(0,  h0.GetMinimum());
  f1.SetParameter(1, (h0.GetMaximum() - h0.GetMinimum()) / (h0.GetXaxis()->GetXmax() - h0.GetXaxis()->GetXmin()));
  f1.SetParameter(2,  0.0);
  f1.SetParameter(3,  0.0);

  f1.SetLineColor(kRed);
  
  // fit
  
  h0.Fit(&f1, option.c_str(), "same");


  cout << "\t// " << lepton << endl;
  cout << "\t// dE/dX = a + bx + cx^2 + dx^3; // [MeV g^-1 cm^2]; x = log10(E/GeV);" << endl;
  
  for (int i = 0; i != 4; ++i) {
    cout << "\tstatic const double " << (char) ('a' + i) << " = " << SCIENTIFIC(10,3) << f1.GetParameter(i) << ";" << endl;
  }

  {
    double Tmin  =  sqrt(EMIN*EMIN - MASS_ELECTRON*MASS_ELECTRON);

    double Emin  = 1 * MASS_LEPTON;
    double Emax  = 5 * MASS_LEPTON;

    for (double E = 0.5 * (Emin + Emax); ; E = 0.5 * (Emin + Emax)) {

      const double gamma = E / MASS_LEPTON;
      const double beta  = sqrt(1.0 - 1.0/(gamma*gamma));
      
      const double Tmax  = (lepton == electron ?
			    (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
			    (2.0*MASS_ELECTRON*beta*beta*gamma*gamma)  /  (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));

      if (fabs(Tmax - Tmin) < precision) {
	cout << "\tstatic const double Emin = " << FIXED(7,5) << E << "; // [GeV]" <<  endl;
	break;
      }
      
      if (Tmax > Tmin)
	Emax = E;
      else
	Emin = E;
    }
  }
  
  out.Write();
  out.Close();
}