#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;
}
Double_t myfunction1(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.+ exp( (par[1]-x[0])/par[2] ) );
    return plat;
}
 
void plateau_singole(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);

   Float_t framemax = 1.1;
   if (Id == "8" || Id == "9") framemax = 0.7;
   TH1F *hr = c1->DrawFrame(.95,0.,1.65,framemax);
   hr->GetXaxis()->SetTitle("HV / kV");
   hr->GetYaxis()->SetTitle("efficiency");

   string fileNameString = "../data/hv_scintillatore"+Id+".txt";
   //string fileNameString = "../dati/eff_vs_tensione/hvscan_scint"+Id+"_th18.dat";
   //string fileNameString = "../dati/hvscan_scint"+Id+"_th18.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];
   Float_t ssingle_v[100];
   Float_t sdsingle_v[100];
   Float_t maxrate=0;
   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)/10.; // per 100 s
       if (maxrate < singole+dsingle_v[nP]) maxrate = singole +dsingle_v[nP];
       nP++;
   }
   if (nP == 0) 
   {
       cout<<" ************************** File "<<fileName<<" not readable ************ "<<endl;
       return;
   }
   
       
   gr = new TGraphErrors(nP,x,eff_v,dx,deff_v);
   string MyTitle = "scintillatore n. "+Id+" - discrimination thr. 18 mV";
   gr->SetTitle(MyTitle.c_str());
   hr->SetTitle(MyTitle.c_str());
   gr->SetMarkerColor(4);
   gr->SetMarkerStyle(21);
   gr->GetXaxis()->SetTitle("HV / kV");
   gr->GetYaxis()->SetTitle("efficiency");
   //gr->Draw("ALP");
   gr->Draw("P");
   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);
   TF1 *f1 = new TF1("f1",myfunction1,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();
   

   Float_t rightmax = 1.1*maxrate;
   if (Id == "5") rightmax = 1200.;
   Float_t scalefact = framemax/rightmax;
   cout<<" framemax "<<framemax<<" rightmax "<<rightmax<<endl;
   for(int i=0; i<nP; i++)
   {
       ssingle_v[i] = single_v[i]*scalefact;
       sdsingle_v[i]=dsingle_v[i]*scalefact;
       cout<<" i "<<i<<" single "<<single_v[i]<<" scaled single "<<ssingle_v[i]<<endl;
   }
   grate = new TGraphErrors(nP,x,ssingle_v,dx,sdsingle_v);
   grate->SetTitle(MyTitle.c_str());
   grate->SetMarkerColor(3);
   grate->SetMarkerStyle(23);
   grate->GetYaxis()->SetTitle("Single Rate /Hz");
   //   grate->Scale(scalefact);
   grate->SetLineColor(kGreen);
   grate->Draw("PL");

   //draw an axis on the right side
   TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
         gPad->GetUxmax(), gPad->GetUymax(),0,rightmax,510,"+L");
   //   axis->SetLineColor(kRed);
   //   axis->SetLabelColor(kRed);
   axis->SetTitle("Single Rate /Hz");
   axis->Draw();
   
   cout<<" ************************************************* Plot drawn "<<endl;

}