#include "TCanvas.h"
#include "TStyle.h"
#include "TH1.h"
#include "TGaxis.h"
#include "TRandom.h"

     Double_t expo (Double_t *x, Double_t *par)
     {return par[0]*exp(-x[0]/par[1])+par[2]*exp(-x[0]/par[3])+par[4];}


    void fitsum(Int_t maxchan, Double_t tlow, Double_t tmax)
    {
      gROOT->Reset(); 
   TCanvas *c1 = new TCanvas("c1","Fit esponenziale", 200,10,1000,800);
   c1->Divide(0,2);

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

   gStyle->SetOptFit(1112);


   string fileNameString = "../camac/dati_050606.txt";
   char* fileName = fileNameString.c_str();
   cout<<" *************** reading from file "<<fileName<<endl;
   ifstream inputFile;
   inputFile.open(fileName);

   Int_t ncanale = 2;
   Int_t nP = 2;
   //Float_t formato1 = 0;
   //Float_t formato2 = 0;
   //Float_t bin1 = 0;
   //Float_t bin2 = 0;
   Float_t x[2500];
   Float_t dx[2500];
   Float_t y[2500];
   Float_t dy[2500]; 
   Float_t x2[2500];
   Float_t dx2[2500];
   Float_t y2[2500];
   Float_t dy2[2500]; 
   Float_t conteggio;
   
   //inputFile >> formato1 >> formato2 >> bin1 >> bin2;
   
   //cout<<"numero canali ADC = " << formato1 <<" num.colonne =  " <<formato2 <<endl;
   //cout <<" counts bin1 =  "<< bin1 <<" counts bin2 = "<< bin2 <<endl;
    

    while (1) {
      //inputFile >> formato1 >>formato2 >> bin1 >>bin2;
     inputFile >> conteggio;
     x[ncanale]= ncanale;
     x2[ncanale] = ncanale;
     dx [ncanale]= 0.;
     dx2 [ncanale]= 0.;
     y[ncanale]= conteggio;
     y2 [ncanale] = conteggio;
     dy [ncanale]= sqrt(conteggio);
     dy2 [ncanale]= sqrt(conteggio);
     if (ncanale>maxchan) break;
       
       
     // cout <<" ncanale = " << ncanale  <<";  counts = " << conteggio <<endl;
        ncanale++;        
   }

        
    Int_t nP = 2;
    while (1)
      {
    if  (nP<=150)
       x2[nP]= nP/138.20;

    if (nP>150 && nP <=720)
      x2[nP]=nP/142.6397;

    if (nP>720 && nP <=1435)
      x2[nP]=nP/143.3775;

    if (nP>1435)
     x2[nP]=nP/137.6391;
        
    if (nP>maxchan) break;
	nP++;
      }
 
   
   TGraphErrors *gr1 = new TGraphErrors(ncanale,x,y,dx,dy);
  
   //gr2 = new TGraphErrors(ncanale,x,y,dx,dy)
   //string MyTitle = "ADC counts";
   gr1->SetTitle("Counts vs ADC channel");
   gr1->SetMarkerColor(8);
   gr1->SetMarkerStyle(22);
   gr1->GetXaxis()->SetTitle(" ADC channel ");
   gr1->GetYaxis()->SetTitle("counts");
   //gr->Draw("ALP");
   c1->cd(1);
   gr1->Draw("AP");
  
   // costruisco un secondo grafico che avrà tempi in x e la corretta calibrazione
   //TF1 *f1 = new TF1 ("f1", "expo", minchan, maxchan);
   //gr1->Fit("f1","R");
   //f1->SetLineColor(kRed);
   //TF1 *f1 = new TF1 ("f1", "[0]*exp(-(x/[1]))", minchan,maxchan);   
   //gr1->Fit("f1");
   
   /* Double_t chi2 = f1->GetChisquare(); */
/*    Double_t p0 = f1-> GetParameter(0); */
/*    Double_t p1   = f1->GetParameter(1); */
   
/*    Double_t ep0   = f1->GetParError(0); */
/*    Double_t ep1   = f1->GetParError(1); */
   
/*    cout <<" Chi2  = "<<chi2<<endl; */
/*    cout <<" Costante = "<<p0<<" +/- "<< ep0<<endl; */
/*    cout <<" Tempo di decadimento "<<p1<<" +/- "<<ep1<<endl; */


   TGraphErrors *gr2 = new TGraphErrors(nP,x2,y2,dx2,dy2);
   gr2->SetTitle("Counts vs time");
   gr2->SetMarkerColor(9);
   gr2->SetMarkerStyle(23);
   gr2->GetXaxis()->SetTitle(" time [#mus] ");
   gr2->GetYaxis()->SetTitle("counts");
   c1->WaitPrimitive();
   c1->cd(2);
   gr2->Draw("AP");

   TF1 *func = new TF1("func",expo,tlow,tmax,5);
   func->SetParNames ("Constant","Decay_time","constant2", "decay_time2", "background");
   func->SetParameter(0,-4);
   func->SetParameter(1,2.2);
   func->SetParameter(2,30);
   func->SetParameter(3,1);
   func->SetParameter(4,2);
   //func->SetParameter(2,7)
   func->SetLineColor(2);
   gr2->Fit("func","R"); 

   /* TF1 *f2 = new TF1 ("f2", "expo", tlow, tmax); */
/*    f2->SetLineColor(2); */
/*    f2->SetLineWidth(4); */
/*    gr2->Fit("f2","R"); */

/*    Double_t chi2 = f2->GetChisquare();  */
/*    Double_t  p0 = f2-> GeParameter(0); */
/*    Double_t p1   = f2->GetParameter(1); */
/*    Double_t tau = -1./p1; */
   
/*    Double_t ep0   = f2->GetParError(0); */
/*    Double_t ep1   = f2->GetParError(1); */
/*    Double_t sigma = ep1/(p1*p1); */
/*    cout <<" Chi2  = "<<chi2<<endl; */
/*    cout <<" Costante = "<<p0<<" +/- "<< ep0<<endl; */
/*    cout <<" Tempo di decadimento "<<tau<<" +/- "<<sigma<<endl; */

   c1->Update();
   }