#include <string>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TGraph.h>
#include <TH2F.h>
#include <TF1.h>
#include <TPad.h>
#include <iostream>


Double_t myfunction(Double_t *x, Double_t *par)
{
  //     par(0) -> asymptotic efficiency
  //     par(1) -> HV of 50% of asymptotic efficiency
  //     par(2) -> HV of 90% - HV of 10% of asymptotic efficiency
  Double_t plat=par[0]/( 1.+ pow(81.,  (par[1]-x[0])/par[2] ) );
  //Double_t plat=par[0]/( 1.+ 81.**( (par[1]-x[0])/par[2] ) );
  return plat;
}


void StoppingPowerIntegrateIdEdx(Float_t Estart = 2.3, Float_t Estep=-1.)
{
  gROOT->Reset();
  Bool_t conf2d  = true;
  Bool_t m_debug = false;

  std::string fileNameD = "NIST_stoppingPower_e_Diamond_finegrid.dat";
  std::string fileNameB = "NIST_stoppingPower_e_Brass_finegrid.dat";

  ifstream inputFileD; //abilita la lettura da file inputFile
  ifstream inputFileB; //abilita la lettura da file inputFile

  inputFileD.open(fileNameD.c_str());
  if (!inputFileD.good())
    {
      cout<<"Error when opening file named  "<<fileNameD<<endl;
      return;
    }
  inputFileB.open(fileNameB.c_str());
  if (!inputFileB.good())
    {
      cout<<"Error when opening file named  "<<fileNameB<<endl;
      return;
    }


  Float_t DiamondDensity = 3.5;  // g/cm3
  Float_t BrassDensity  = 8.5;  // g/cm3



  Int_t   nPD = 0;
  Int_t   nEpointsD = 0;
  Float_t EnergyD    [1000];
  Float_t rad_dedxD  [1000];
  Float_t col_dedxD  [1000];
  Float_t tot_dedxD  [1000];
  Int_t   nPB = 0;
  Int_t   nEpointsB = 0;
  Float_t EnergyB    [1000];
  Float_t rad_dedxB  [1000];
  Float_t col_dedxB  [1000];
  Float_t tot_dedxB  [1000];


  std::string line;
  Int_t nrow =0;
  std::cout<<std::endl<<"**************************************************"<<std::endl;
  std::cout<<"Reading from file "<<fileNameD<<std::endl;
  std::cout<<"**************************************************"<<std::endl<<std::endl;
  while ( (!inputFileD.eof()) && nrow<8)
    {
      getline (inputFileD,line);
      std::cout<<line<<std::endl;
      ++nrow;
    }
  while (1) {
    Float_t e = 0;
    Float_t cdedx = 0.; Float_t rdedx = 0; Float_t tdedx = 0.;
    Float_t range = 0.; Float_t ry = 0.;   Float_t dc = 0.;
    inputFileD >> e >> cdedx >> rdedx >> tdedx;
    if (!inputFileD.good()) break;
    EnergyD[nPD]=e;
    rad_dedxD[nPD]=rdedx*DiamondDensity;
    col_dedxD[nPD]=cdedx*DiamondDensity;
    tot_dedxD[nPD]=tdedx*DiamondDensity;
    //       density_c[nP]=dc;
    cout<<"nPD = "<<nPD<<" energy = "<<EnergyD[nPD]<<" t_dE/dx "<<tot_dedxD[nPD]<<std::endl;
    ++nPD;
  }
  nEpointsD = nPD;
  inputFileD.close();
  nrow =0;
  std::cout<<std::endl<<"**************************************************"<<std::endl;
  std::cout<<"Reading from file "<<fileNameB<<std::endl;
  std::cout<<"**************************************************"<<std::endl<<std::endl;
  while ( (!inputFileB.eof()) && nrow<8)
    {
      getline (inputFileB,line);
      std::cout<<line<<std::endl;
      ++nrow;
    }
  while (1) {
    Float_t e = 0;
    Float_t cdedx = 0.; Float_t rdedx = 0; Float_t tdedx = 0.;
    Float_t range = 0.; Float_t ry = 0.;   Float_t dc = 0.;
    inputFileB >> e >> cdedx >> rdedx >> tdedx;
    if (!inputFileB.good()) break;
    EnergyB[nPB]=e;
    rad_dedxB[nPB]=rdedx*BrassDensity;
    col_dedxB[nPB]=cdedx*BrassDensity;
    tot_dedxB[nPB]=tdedx*BrassDensity;
    //       density_c[nP]=dc;
    cout<<"nPB = "<<nPB<<" energy = "<<EnergyB[nPB]<<" t_dE/dx "<<tot_dedxB[nPB]<<std::endl;
    ++nPB;
  }
  nEpointsB = nPB;
  inputFileB.close();
  //   return;
  TCanvas* c1 = new TCanvas("stopping_power_Diamond","Stopping Power in Diamond",200,10,700,500);
//  TCanvas* c2 = new TCanvas("stopping_power_Brass",  "Stopping Power in Brass",  200,10,700,500);

  string MyTitle;
  TGraph* gr_tdedxD     = new TGraph(nEpointsD, EnergyD, tot_dedxD);
  TGraph* gr_cdedxD     = new TGraph(nEpointsD, EnergyD, col_dedxD);
  MyTitle = "Stopping power in Diamond (3.5g/cm3) and Brass60-40 (8.5g/cm3)";
  gr_tdedxD->SetTitle(MyTitle.c_str());
  gr_tdedxD->GetXaxis()->SetTitle("Energy / MeV");
  gr_tdedxD->GetYaxis()->SetTitle("Stopping power in Diamond [Mev cm2 / g]");
  gr_tdedxD->GetXaxis()->SetLabelOffset(0.015);
  gr_tdedxD->GetYaxis()->SetLabelOffset(0.01);
  gr_tdedxD->GetYaxis()->SetTitleOffset(1.3);
  gr_tdedxD->GetXaxis()->SetTitleOffset(1.3);

  TGraph* gr_tdedxB     = new TGraph(nEpointsB, EnergyB, tot_dedxB);
  //MyTitle = fileNameB;
  gr_tdedxB->SetTitle(MyTitle.c_str());
  gr_tdedxB->GetXaxis()->SetTitle("Energy / MeV");
  gr_tdedxB->GetYaxis()->SetTitle("Stopping power in Diamond [Mev cm2 / g]");
  gr_tdedxB->GetXaxis()->SetLabelOffset(0.015);
  gr_tdedxB->GetYaxis()->SetLabelOffset(0.01);
  gr_tdedxB->GetYaxis()->SetTitleOffset(1.3);
  gr_tdedxB->GetXaxis()->SetTitleOffset(1.3);


  c1->SetLogy();
  c1->SetLogx();
  c1->cd();
  //c1->cd(1);
  TH2F* yy = new TH2F("a","a",100,0.05,100.,100.,4.,40.);
  yy->SetTitle(MyTitle.c_str());
  yy->GetXaxis()->SetTitle("Energy / MeV");
  yy->GetYaxis()->SetTitle("Stopping power [Mev / cm]");
  yy->GetXaxis()->SetLabelOffset(0.015);
  yy->GetYaxis()->SetLabelOffset(0.01);
  yy->GetYaxis()->SetTitleOffset(1.3);
  yy->GetXaxis()->SetTitleOffset(1.3);
  yy->Draw();
  gr_cdedxD->Draw("SAME");
  gr_tdedxD->Draw("SAME");
  c1->SaveAs(std::string(fileNameD+"tdedx.png").c_str());
  /*
    c2->SetLogy();
    c2->SetLogx();
    c2->cd();
    //c1->cd(2);
    yy->Draw();
  */
  gr_tdedxB->Draw("SAME");
  //c2->SaveAs(std::string(fileNameB+"range.png").c_str());
  c1->SaveAs(std::string("StoppingPower_Diamond_Brass.png").c_str());


  Float_t BrassThickness = 0.05; //cm
  Float_t xStep = 0.0001; // 1. micros



  TH2F* tt;
  if (BrassThickness>0.0001) tt = new TH2F("box","box",100,0.,0.2,100,0.5,2.);
  else tt = new TH2F("box","box",100,0.,0.2,100,0.5,0.8);
  gStyle->SetOptStat(0);
  TCanvas* c3 = new TCanvas("Elostvsdepth","E lost vs depth",200,10,700,500);
  //c3->SetLogy();

  tt->SetTitle("Energy loss/micron vs depth [KeV]");
  tt->GetXaxis()->SetTitle("depth / cm");
  tt->GetYaxis()->SetTitle("E lost in 1 um [KeV]");
  tt->GetXaxis()->SetLabelOffset(0.015);
  tt->GetYaxis()->SetLabelOffset(0.01);
  tt->GetYaxis()->SetTitleOffset(1.3);
  tt->GetXaxis()->SetTitleOffset(1.3);
  tt->Draw();



  Bool_t firstTimeInCalibD = true;
  Bool_t firstTimeInBrass = true;
  Bool_t firstTimeInTrgD = true;
  Bool_t first = true;
  Float_t ElostinDUT = 0.;
  Float_t ElostinCD = 0.;
  Float_t ElostinBrass = 0.;
  Float_t ElostinTrgD = 0.;
  Float_t finalEnergy = 0.;
  Float_t DiamondThickness = 0.05;

  while(Estart>0.7) {
    std::cout<<"---------------------------------------------"<<std::endl;
    Int_t nsamples = 0;
    Float_t ElossInSample[10000];
    Float_t DepthOfSample[10000];
    Float_t Ene = Estart;
    Double_t Elost = 0;
    Double_t depth = 0.;
    Float_t totThickness = 3.*DiamondThickness+BrassThickness;
    if (conf2d) totThickness = 2.*DiamondThickness+BrassThickness;
    while (Ene >0.){
      if (!conf2d) {
	if (depth < DiamondThickness) {// in DUT
	  Elost = xStep * gr_tdedxD->Eval(Ene);
	  if (m_debug) std::cout<<" In DUT         Diamond ";
	}
	else if (depth < 2.*DiamondThickness)
	  {// in Single Crystal Diamond
	    if (firstTimeInCalibD)
	      {
		firstTimeInCalibD = false;
		ElostinDUT = Estart - Ene;
	      }
	    Elost = xStep * gr_tdedxD->Eval(Ene);
	    if (m_debug) std::cout<<" In Calibration Diamond ";
	  }
	else if ( depth < (2.*DiamondThickness+BrassThickness) )
	  {// in Brass
	    if (firstTimeInBrass)
	      {
		firstTimeInBrass = false;
		ElostinCD = - Ene + (Estart -ElostinDUT) ;
	      }
	    Elost = xStep * gr_tdedxB->Eval(Ene);
	    if (m_debug) std::cout<<" In energy filter Brass ";
	  }
	else {// in Triggering Diamond
	  if (first)
	    {
	      first = false;
	      ElostinBrass = - Ene + (Estart - ElostinDUT - ElostinCD) ;
	      std::cout<<" Elost in Brass = "<<ElostinBrass<<"MeV;  Brass thickness = "<<BrassThickness<<" cm"<<std::endl;
	    }
	  Elost = xStep * gr_tdedxD->Eval(Ene);
	  if (m_debug) std::cout<<" In Triggering  Diamond ";
	}
	depth = depth + xStep;
	if (depth > totThickness) {
          finalEnergy = Ene;
          std::cout<<" elecron survives the materials of the setup; final energy = "<<Ene<<" starting Energy = "<<Estart<<" MeV"<<std::endl;
          break;
        }
	Ene = Ene - Elost;
	ElossInSample[nsamples]=Elost/(xStep*10.);
	DepthOfSample[nsamples]=depth;
	++nsamples;


	if (m_debug) std::cout<<" initial_E = "<<Estart<<" depth so far "<<depth<<" Elost"<<Elost<<" Energy "<<Ene<<" MeV"<<std::endl;
      }
      else if(conf2d)
	{
	  ElostinCD = 0.;
	  if (depth < DiamondThickness) {// in DUT
	    Elost = xStep * gr_tdedxD->Eval(Ene);
	    if (m_debug) std::cout<<" In DUT         Diamond ";
	  }
	  else if ( depth < (DiamondThickness+BrassThickness) )
	    {// in Brass
	      if (firstTimeInBrass)
		{
		  firstTimeInBrass = false;
		  ElostinDUT = - Ene + Estart ;
		  std::cout<<"2d: First Time in Brass - Estart, Ene => ElostInDUT "<<Estart<<" "<<Ene<<" "<<ElostinDUT<<std::endl;
		}
	      Elost = xStep * gr_tdedxB->Eval(Ene);
	      if (m_debug) std::cout<<" In energy filter Brass ";
	    }
	  else
	    {// in Triggering Diamond
	      if (first)
		{
		  first = false;
		  ElostinBrass = - Ene + (Estart - ElostinDUT) ;
		  std::cout<<"AAAAAAAAAAAAAAAAAAA Elost in Brass = "<<ElostinBrass<<"MeV;  Brass thickness = "<<BrassThickness<<" cm"<<std::endl;
		  std::cout<<"2d: First Time in Trg Diamond - Estart, Ene, ElostInDUT => ElostInBrass "<<Estart<<" "<<Ene<<" "<<ElostinDUT<<" "<<ElostinBrass<<std::endl;
		}
	      Elost = xStep * gr_tdedxD->Eval(Ene);
	      if (m_debug) std::cout<<" In Triggering  Diamond ";
	    }
	  depth = depth + xStep;
	  if (depth > totThickness) {
	    finalEnergy = Ene;
	    std::cout<<" elecron survives the materials of the setup; final energy = "<<Ene<<" starting Energy = "<<Estart<<" MeV"<<std::endl;
	    break;
	  }
	  Ene = Ene - Elost;
	  ElossInSample[nsamples]=Elost/(xStep*10.);
	  DepthOfSample[nsamples]=depth;
	  ++nsamples;


	  if (m_debug) std::cout<<" initial_E = "<<Estart<<" depth so far "<<depth<<" Elost"<<Elost<<" Energy "<<Ene<<" MeV"<<std::endl;
	}
    }
    ElostinTrgD = - Ene + (Estart - ElostinDUT - ElostinCD - ElostinBrass);
    std::cout<<" Estart "<<    Estart <<" MeV\n Elost in DUT "<<ElostinDUT<<" MeV\n in Calib Diamond "<<ElostinCD<<" MeV\n in Brass "<<ElostinBrass
	     <<" MeV Brass thickness = "<<BrassThickness<<" cm\n in Triggering Diamond "<<ElostinTrgD<<" MeV"<<std::endl;
    std::cout<<" ElostinDUT+ElostinCD+ElostinCD+ElostinTrgD = "<<ElostinDUT+ElostinCD+ElostinCD+ElostinTrgD<<" Summing  FinalEnergy = "<<ElostinDUT+ElostinCD+ElostinCD+ElostinTrgD+finalEnergy
	     <<std::endl;
    std::cout
      <<"\n Estart                                                                "<< Estart
      <<" MeV \n Estart - ElostinDUT                                                   "<< Estart - ElostinDUT
      <<" MeV \n Estart - ElostinDUT - ElostinCD                                       "<< Estart - ElostinDUT - ElostinCD
      <<" MeV \n Estart - ElostinDUT - ElostinCD - ElostinBrass                        "<< Estart - ElostinDUT - ElostinCD -ElostinBrass
      <<" MeV \n Estart - ElostinDUT - ElostinCD - ElostinBrass - ElostinTrgD =        "<<Estart - ElostinDUT - ElostinCD - ElostinBrass - ElostinTrgD<<" MeV"<<std::endl;
    if (!conf2d) std::cout<<" Effective track length in triggering Diamond = "<<depth - 2.*DiamondThickness - BrassThickness<<" cm"<<std::endl;
    else std::cout<<" Effective track length in triggering Diamond = "<<depth - DiamondThickness - BrassThickness<<" cm"<<std::endl;


    std::string MyTitleBragg = "E loss vs depth";
    TGraph* gr_bragg     = new TGraph(nsamples, DepthOfSample, ElossInSample);
    MyTitle = MyTitleBragg;
    gr_bragg->SetTitle(MyTitle.c_str());
    gr_bragg->GetXaxis()->SetTitle("depth / cm");
    gr_bragg->GetYaxis()->SetTitle("E lost in 1 um [KeV]");
    gr_bragg->GetXaxis()->SetLabelOffset(0.015);
    gr_bragg->GetYaxis()->SetLabelOffset(0.01);
    gr_bragg->GetYaxis()->SetTitleOffset(1.3);
    gr_bragg->GetXaxis()->SetTitleOffset(1.3);

    gr_bragg->Draw("SAME");
    std::cout<<"At E = "<<Estart<<" Eloss spread in DUT = "<<(gr_bragg->Eval(0.05)-gr_bragg->Eval(0.))/gr_bragg->Eval(0.)<<std::endl;
    if (Estep<0.) break;
    Estart = Estart - Estep;
    std::cout<<" Redefine Estart = "<<Estart<<std::endl;
  }
  c3->SaveAs("Elost_vs_depth.png");

  return;

  //gStyle->SetOptFit(1112); // probabilità , chiquadro, errore, nome e valore parametri (0=niente,1=stampa,per l'ultimo 2=stampa tutti i parametri anche quelli di valore fissato)
}