#include <string>
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.+ 81.**( (par[1]-x[0])/par[2] ) );
    return plat;
}
 
void fitplateau(string Id="1")
{
    cout<<"\n ****************** processing scintillator "<<Id<<" *****"<<endl;
   //
   // To see the output of this macro, click  here. 
   //
   gROOT->Reset();
   //c1 = new TCanvas("c1","Efficiency plateau",200,10,700,500);
   c1 = new TCanvas("c1","Efficiency plateau");


   gPad->SetFillColor(0);
   gPad->SetGrid();
   gPad->GetFrame()->SetFillColor(21);
   gPad->GetFrame()->SetBorderSize(12);

   gStyle->SetOptFit(1112);


   string fileNameString = "../data/hv_scintillatore"+Id+".txt";
   char* fileName = fileNameString.c_str();
   
   cout<<" *************************************** reading from file "<<fileName<<endl;
   ifstream inputFile;
   inputFile.open(fileName);
   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];
   // vth (in kV), triple, doppie, rate (Hz)
   Float_t Dt = 100.; // 100 s di tempo di misura 
   while (1) {
       Float_t  vth = 0; Float_t triple = 0.; Float_t quadruple = 0; Float_t singole = 0.;
       inputFile >> vth >> triple >> quadruple >> singole;
       if (!inputFile.good()) break;
       Float_t eff = quadruple/triple;
       Float_t deff = sqrt(eff*(1.-eff)/triple);
       //cout<<"reading from file "<<fileName<<" HV "
       //    <<vth<<" t/q/s  "<<triple<<" "<<quadruple<<" "<<singole<<endl;
       cout<<" HV "
           <<vth<<" eff "<<eff<<"+/-"<<deff<<" single rate "<<singole<<" 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] = singole;
       dsingle_v[nP] = sqrt(singole)/sqrt(Dt); 
       nP++;
   }
   
   gr = new TGraphErrors(nP,x,eff_v,dx,deff_v);
   string MyTitle = "scintillatore n. "+Id+" - discrimination thr. 18 mV";
   gr->SetTitle(MyTitle.c_str());
   gr->SetMarkerColor(4);
   gr->SetMarkerStyle(21);
   gr->GetXaxis()->SetTitle("HV / kV");
   gr->GetYaxis()->SetTitle("efficiency");
   //gr->Draw("ALP");
   gr->Draw("AP");
   cout<<" ************************************************* Plot drawn "<<endl;


   // definiamo una nostra funzione entro il range definito da lowerLimit e upperLimit
   Double_t lowerLimit = 1.0;
   Double_t upperLimit = 1.6;
   // la funzione myfunction e` denominate "f1" - ha 4 parametri - e` definita entro i limiti fissati
   TF1 *f1 = new TF1("f1",myfunction,lowerLimit,upperLimit,3);
   // fissiamo i valori (iniziali) dei parametri 
   f1->SetParameter(0,1.);
   f1->SetParameter(1,1.246);
   f1->SetParLimits(1,1.1,1.3);
   f1->SetParameter(2,0.1875);
   //f1->SetParLimits(2,0.05,0.195);
   f1->SetLineColor(kRed);
   //f1->Draw("L");
   //fitta entro il range (opzione "R") specificato nella definizione della funzione
   gr->Fit("f1","R"); 
   cout<<" ************************************************* Fit done "<<endl;

   c1->Update();
}