#include <string>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TF1.h>
#include <TPad.h>


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 fitPlateau_showRate_root6(string Id="XX"){ 
   gROOT->Reset(); 

   if (Id=="XX")
     {
       gSystem->Exec("ls -lart hv_plateau*.dat");
       //       std::cout<<"Check *.dat with \nls -lart hv_plateau*.dat"<<std::endl;
       std::cout<<"Please run with one of the above input files\n Example: \n.x fitPlateau_showRate(\"15atLS5_doppie_LS1_LS3\")"<<std::endl;
       return;
     }
   // create first graph

   std::string fileNameString = "hv_plateau_scint_"+Id+".dat";
   //char* fileName = fileNameString.c_str();
   ifstream inputFile; //abilita la lettura da file inputFile è un comando per caricare il file

   //inputFile.open(fileName);
   inputFile.open(fileNameString.c_str());
   if (!inputFile.good())
     {
       cout<<"errore nell'apertura del file "<<endl;
       return;
     } 

   Int_t   nP = 0;
   Float_t x[100];
   Float_t dx[100];
   Float_t eff_v[100];
   Float_t deff_v[100];
   Float_t single_v[100];
   Float_t dsingle_v[100];
   while (1) {
       Float_t vth = 0; Float_t doppie = 0.; Float_t triple = 0; Float_t singole = 0.; 
       Float_t time=0.;
       inputFile >> vth >> doppie >> triple >> singole >> time;
       if (!inputFile.good()) break;
       Float_t eff = triple/doppie;
       Float_t deff = sqrt(eff*(1.-eff)/doppie);
       Float_t rate = singole/time;
       //cout<<"reading from file "<<fileName<<" HV "
       //    <<vth<<" t/q/s  "<<triple<<" "<<quadruple<<" "<<singole<<endl;
       cout<<" HV "<<vth<<" eff "<<eff<<"+/-"<<deff<<" singole "<<singole<<" rate di singola "<< rate <<"Hz" <<endl;
       // qui si riempiono vettori di eff, err-su-eff,
       // rate di singola e poi li si plotta come i vettori in basso
       x[nP]      = vth;
       //dx[nP]     = 0.001; // trascuriamo l'rrore sulla tensione di discriminazione
       dx[nP]     = 0.; // quando l'errore e` piccolo ma non 0, il fit talvolta non converge perche'
                        // una derivata alta della funzione che moltiplica sigma_x al denominatore del chi^2
                        // rende il chi^2 piccolo anche in presenza di residui grandi !!! 
       eff_v[nP]  = eff;
       deff_v[nP] = deff;
       single_v[nP] = rate*0.001;
       dsingle_v[nP] = 0.001*sqrt(singole)/time; 
       nP++;
       cout<<"nP = "<<nP<<std::endl;
   }
   inputFile.close();
   TCanvas* c1 = new TCanvas("c1","SCINTILLATORS",200,10,700,500);
   TPad *pad = new TPad("pad","",0,0,1,1);

   //   pad->SetFillColor(42);
   pad->SetGrid();
   pad->Draw();
   pad->cd();

   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)



   // draw a frame to define the range
   
   float hvmax = x[nP-1]+0.1;
   float effmax = eff_v[nP-1]+0.05;
   TH1F *hr = c1->DrawFrame(1.1,0.,hvmax,effmax);
   string MyTitle = "scintillatore "+Id+"";
   hr->SetTitle(MyTitle.c_str());
   hr->SetXTitle("HV / kV");
   hr->SetYTitle("Efficienza");
   hr->GetXaxis()->SetLabelOffset(0.015);
   hr->GetYaxis()->SetLabelOffset(0.01);
   hr->GetYaxis()->SetTitleOffset(1.3);
   hr->GetXaxis()->SetTitleOffset(1.3);
   //   pad->GetFrame()->SetFillColor(21);
   pad->GetFrame()->SetBorderSize(12);

   //   Style_t tFont=hr->GetXaxis()->GetTitleFont();
   Style_t tFont=hr->GetXaxis()->GetTitleFont();
   Float_t tSize=hr->GetXaxis()->GetTitleSize();
 
     



   TGraphErrors* gr     = new TGraphErrors(nP,x,eff_v,dx,deff_v);//grafico con barre di errore
   TGraphErrors* grRate = new TGraphErrors(nP,x,single_v,dx,dsingle_v);//grafico con barre di errore

   gr->SetMarkerColor(kBlue);
   gr->SetMarkerSize(0.6);
   gr->SetMarkerStyle(20);
   gr->Draw("P");

   // definiamo una nostra funzione entro il range definito da lowerLimit e upperLimit
   Double_t lowerLimit = 0.725;
   Double_t upperLimit = 1.95;
   // la funzione myfunction e` denominate "f1" - ha 4 parametri - e` definita entro i limiti fissati
   TF1 *f1 = new TF1("f1",myfunction,lowerLimit,upperLimit,3); //3 è il numero dei parametri e tf1 è un oggetto usato per i fit
   // fissiamo i valori (iniziali) dei parametri 
   f1->SetParameter(0,eff_v[nP-1]);//setta il parametro 0 a 1
   f1->SetParameter(1,1.6);
   //f1->SetParLimits(1,1.1,1.3);//setta i limiti del parametro 1 fra 1.1 1.3
   f1->SetParLimits(1,1.2,1.7);//l'intervallo è stato usato solo per lo scintillatore 7
   f1->SetParameter(2,0.2);
   std::cout<<"Fit starting values are: "<< eff_v[nP-1]<<" "<<" 1.4  0.2"<<std::endl;
   //f1->SetParLimits(2,0.05,0.195);
   f1->SetLineColor(7);
   //f1->Draw("L");
   //fitta entro il range (opzione "R") specificato nella definizione della funzione
   gr->Fit("f1","R"); //fitta l'istogramma f1 e con l'opzione R rispetta i parametri stabiliti in f1

   // Set stat options
   gStyle->SetStatY(0.9);                
   // Set y-position (fraction of pad size)
   gStyle->SetStatX(0.45);                
   // Set x-position (fraction of pad size)
   gStyle->SetStatW(0.17);                
   // Set width of stat-box (fraction of pad size)
   gStyle->SetStatH(0.12);                
   // Set height of stat-box (fraction of pad size)
   cout<<" ************************************************* Fit done "<<endl;


      // create second graph


   Int_t n2 = 10;
   Float_t x2[]  = {-0.28, 0.005, 0.19, 0.29, 0.45, 0.56,0.65,0.80,0.90,1.01};
   Float_t y2[]  = {0.82,3.86,7,9,10,10.55,9.64,7.26,5.42,2};
   Float_t ex2[] = {.04,.12,.08,.06,.05,.04,.07,.06,.08,.04};
   Float_t ey2[] = {.6,.8,.7,.4,.3,.3,.4,.5,.6,.7};


   for (Int_t i=0;i<n2;i++) {
      y2[i] *= 100;
      ey2[i] *= 100;
   } 
   //create a transparent pad drawn on top of the main pad    c1->cd();
   TPad *overlay = new TPad("overlay","",0,0,1,1);
   //overlay->SetLogy();

   overlay->SetFillStyle(4000);
   overlay->SetFillColor(0);
   overlay->SetFrameFillStyle(4000);
   overlay->Draw();
   overlay->cd();



   grRate->SetMarkerColor(kRed);
   grRate->SetLineColor(6);
   grRate->SetLineWidth(2);
   grRate->SetLineStyle(2);
   grRate->SetMarkerStyle(23);
   grRate->SetMarkerSize(1.);
   //   grRate->SetName("gr2");
   Double_t xmin = pad->GetUxmin();
   Double_t ymin = 0;

   Double_t xmax = pad->GetUxmax(); 
   Double_t ymax = single_v[nP-1]*1.1;//500000.; 
   TH1F *hframe = overlay->DrawFrame(xmin,ymin,xmax,ymax); 
   //hframe->GetYaxis()->SetTickLength->(0.); 
   hframe->GetXaxis()->SetLabelOffset(99);    
   hframe->GetYaxis()->SetLabelOffset(99);    
   grRate->Draw("LPE");
   //Draw an axis on the right side 
   TGaxis *axis = new TGaxis(xmax,ymin,xmax, ymax,ymin,ymax,510,"+L");

   axis->SetTitleFont (tFont);
   axis->SetTitleSize (tSize);
   axis->SetTitleColor(kRed);
   axis->SetLineColor(kRed);
   axis->SetLabelColor(kRed);
   axis->SetLabelFont(tFont);
   axis->SetLabelSize(tSize);
   axis->SetTitle("Rate in singola / kHz");
   axis->SetLabelOffset(0.01);
   axis->SetTitleOffset(1.3);
   axis->Draw();
   c1->SaveAs(std::string(fileNameString+".png").c_str());
}