#include "TMath.h"
#include "TFile.h"
#include "TH1.h"
#include "THStack.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include "TF1.h"

//
// CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers
// D. Heck, J. Knapp, J.N. Capdevielle, G. Schatz, T. Thouw
// Forschungszentrum Karlsruhe Report FZKA 6019 (1998)
// https://web.ikp.kit.edu/corsika/physics_description/corsika_phys.pdf
// "7.2.2 Nishimura-Kamata-Greisen option" -> pages 51 to 54
// "Lateral electron distribution" -> pages 53 and 54

//
 char datafile1[200]= "G:/cors/gamma/gammatot/cors_plot-gammatot.root";
 char datafile2[200]= "G:/cors/prot/prottot/cors_plot-prottot.root";
 char dataoutfit[200]="G:/cors/gamma/gammaprot_distr_elelat_fit.png";


double NKG_rho_e_e(const double *x, const double *p) {
  double r = x[0];
  double N_e = p[0];
  double s = p[1];
  double r_mol = p[2];
  double a=p[3];
  //double r_mol=100;
  double rho_e_e = 1.0;
#if 0 // 0 or 1 
  rho_e /= TMath::TwoPi();
  rho_e *=
    TMath::Gamma(4.5 - s) /
    TMath::Gamma(s) /
    TMath::Gamma(4.5 - 2.0 * s);
  r_mol *= 0.78 - 0.21 * s; // r_mol * s_m
#endif // 0 or 1 
  r /= r_mol;
  //rho_e_e *= TMath::Power(r, s - 2.0) * TMath::Power(1.0 + r, s - 4.5);
   rho_e_e *= TMath::Power(r,  -a) * TMath::Power(1.0 + r, -s-a);
  rho_e_e *= N_e / r_mol / r_mol;
  return rho_e_e;
}

void distrelelat()
{
   TFile *f1 = TFile::Open(datafile1);
   if ((!f1) || f1->IsZombie()) {
      // if we cannot open the file, print an error message and return immediatly
      printf("Error: cannot open cors_plot.root!\n");
      delete f1;
      return;
   }
   f1->ls();
   TFile *f2 = TFile::Open(datafile2);
   if ((!f2) || f2->IsZombie()) {
      // if we cannot open the file, print an error message and return immediatly
      printf("Error: cannot open cors_plot.root!\n");
      delete f2;
      return;
   }
   f2->ls();
  
   TH1F *helelat1; f1->GetObject("helelat", helelat1);
   if (!(helelat)) {
      printf("Error getting an histogram 1 from the file!\n");
      delete f1; delete f2;  
      return;
   }
   TH1F *helelat2; f2->GetObject("helelat", helelat2);
    if (!(helelat)) {
      printf("Error getting  histogram 2 from the file!\n");
      delete f1; delete f2;  
      return;
   }
  
   helelat1->SetLineColor(kBlue);
   helelat1->SetName("EAS ind. #gamma");
   helelat2->SetName("EAS ind. #font[12]{p}");
   helelat2->SetLineColor(kGreen);
   
   THStack *hs = new THStack("hs", "Distribuzione laterale degli e^{#pm} negli EAS indotti da #gamma e #font[12]{p};Distanza dal core (m);Densita' (e^{#pm}/m^{2}_{})");
   hs->Add(helelat1);
   hs->Add(helelat2);
   
   gStyle->SetOptStat("nemr"); // "nemr" or "neMR"
   gStyle->SetOptFit(112);
   
   float offx=1.0;
   float offy=1.2;
   float n=3;
   float marg=0.1;
   TCanvas *c = new TCanvas("c5","hists with different scales",1280,1000);
   c->SetRightMargin(marg);
   c->SetLogx();
   c->SetLogy();
   hs->Draw("NOSTACK");
   c->Update();
   hs->GetXaxis()->SetTitleOffset(offx);
   hs->GetYaxis()->SetTitleOffset(offy);
   helelat1->Draw("sames");
   helelat2->Draw("sames");
  
   c->Update();
    TLegend* leg = new TLegend(0.7, 0.7, .6, .80);
   leg->SetNColumns(1);
   leg->AddEntry(helelat1, "EAS indotti da #gamma", "l");
   leg->AddEntry(helelat2, "EAS indotti da #font[12]{p}", "l");
   
   leg->Draw();
   // now retrieve each stats box and reposition them
   TPaveStats *stats1 = (TPaveStats*)helelat1->GetListOfFunctions()->FindObject("stats");
   TPaveStats *stats2 = (TPaveStats*)helelat2->GetListOfFunctions()->FindObject("stats");
 
   stats1->SetTextColor(kBlue);
   stats2->SetTextColor(kGreen);
   
   stats1->SetX1NDC(0.80); stats1->SetX2NDC(0.98);
   stats1->SetY1NDC(0.77); stats1->SetY2NDC(0.92);
   stats2->SetX1NDC(0.80); stats2->SetX2NDC(0.98);
   stats2->SetY1NDC(0.60); stats2->SetY2NDC(0.75);
  
   c->Update();

   c->Print(dataoutfit);
   Double_t r_mole = 89.;
   Double_t s = 1.5;
   Double_t a=5;
   
  
   TF1 *fitdistrlatele = new TF1("fitdistrlatele", NKG_rho_e_e, 0, 1e4, 3);
   fitdistrlatele->SetParNames("c*N_{e^{#pm}}", "s", "r_{M}");
   fitdistrlatele->SetParameters(1, s, r_mole,a);
   fitdistrlatele->FixParameter(2, r_mole); // fix "r_mol"
 
   const char *fit_options = ""; // e.g. "" or "L" or "W" or "WW"
   fitdistrlatele->SetParameters(1, s, r_mole,a);
   fitdistrlatele->SetParameter(0, helelat->GetBinContent(1)/fitdistrlatele->Eval(helelat->GetBinCenter(1)));
   //fitdistrlatele->SetLineStyle(8);
   helelat1->Fit(fitdistrlatele, fit_options);
   helelat2->Fit(fitdistrlatele, fit_options);
 
   c->Update();
   c->Print(dataoutfit);
   delete leg;
   delete hs;
   delete c;
   delete f1; delete f2; // automatically deletes all histograms
   delete fitdistrlatele;
   exit();
}