#include "TH1.h"
#include "TF1.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "Riostream.h"
// fattori di calibrazione sono 66.20+-0.022 channels/microsecondo e
//								3.151 +- 0.2739 channels.
// then t=(nchan-3.151)/66.20

// calibration function, convert adc channel to time
//Double_t t_calibrate(Double_t chan) {return (chan-3.151)/66.20;}  // calib sabrina
//Double_t t_calibrate(Double_t chan) {return (chan+4.351)/66.88;}  //calib Stefania (>microsec)
//Double_t t_calibrate(Double_t chan) {return (chan-2.745)/65.06;}  //calib Stefania (<1microsec) better
//Double_t t_calibrate(Double_t chan) {return (chan+0.1836)/69.49;}  //calib Stefania (<1microsec 14-05-2005) 
Double_t t_calibrate(Double_t chan) {return (chan-1.944)/65.79;}  //calib Stefania (all; Dt moved on vert.scale)
//Double_t t_calibrate(Double_t chan)
//{
//    if (chan>250.) return (chan+4.351)/66.88;
//    else return (chan-2.745)/65.06;
//}


// flat background function
Double_t background(Double_t *x, Double_t *par) {return par[0];}

// exponential function
Double_t exponential(Double_t *x, Double_t *par) {return (par[0]*par[1])*(exp(-x[0]/par[1]));}

// exponential + background function
Double_t expobckg(Double_t *x, Double_t *par) {return background(x,par) + exponential(x,&par[1]);}

// Sum of background and two exponential functions
Double_t fitFunction(Double_t *x, Double_t *par) {
  return background(x,par) + exponential(x,&par[1])+exponential(x,&par[3]);
}
// Sum  two exponential functions
Double_t twoexpo(Double_t *x, Double_t *par) {return exponential(x,&par[0])+exponential(x,&par[2]);}

void mufit(int limit=1000,int fact=4, float lowl=2., float highl=14.0){

  gROOT->Reset();
  gStyle->SetOptStat(1);
  gStyle->SetOptFit(1);

  gBenchmark->Start("canvas");
  

// create and set canvas style
   TCanvas *c1 = new TCanvas("c1","Fitting Example",10,10,1000,800);
//   c1->SetFillColor(33);
   c1->SetFillColor(17);
//   c1->SetFrameFillColor(41);
   c1->SetFrameFillColor(40);
   c1->SetGrid();
   c1->Divide(2,2);
   c1->cd(1);
   c1->SetLogy(1);
   
   
   TH1F *histo = new TH1F("histo","Laboratorio Alte Energie",limit,0.0,limit);
   histo->SetMarkerStyle(21);
   histo->SetMarkerSize(0.8);
   TH1F *resid = new TH1F("resid","Laboratorio Alte Energie",limit,0.0,limit);
   resid->SetMarkerStyle(21);
   resid->SetMarkerSize(0.8);
   float halfbin=fact*(t_calibrate(11)-t_calibrate(10))/2.;
   float halfbin_ch=fact/2.;

   printf("\n If rebinning by %d halfbin IN TIME  is %8f microseconds \n ",fact,halfbin);
   printf("\n If rebinning by %d halfbin IN ADCch is %8f  \n ",fact,halfbin_ch);
   //edo   TH1F *time = new TH1F("time","Time calibrated",limit/fact,t_calibrate(0.0)-halfbin,t_calibrate(limit)-halfbin);
   float xmin = t_calibrate(0.0-halfbin_ch);
   float xmax = t_calibrate(limit-1-halfbin_ch);
   int nbin = (int)limit/fact;
   TH1F *time = new TH1F("time","ADC counts", nbin, xmin, xmax);
   time->SetXTitle("Time / micros");
   time->SetYTitle("yields / ch");
   
   printf("Histo. min, max, Nbin are  %8f %8f %d \n ",xmin, xmax, nbin);
   resid->SetMarkerStyle(21);
   resid->SetMarkerSize(0.8);

  
// here open the data file and fills the histogram
   Float_t Entries;
   Int_t nlines = 0;
   ifstream in;
   //in.open("edo.dat");
   //   in.open("dati_veto_16052005_linux.text"); // fittando da 2 a 14 1 exp bin 1/1 calib stefania Dt on vert scale OK 
   in.open("dati_veto_16052005_linux.text"); // fittando da 0.5(0.2) a 14 2 exp bin 1/1 calib stefania Dt on vert scale OK
   //   in.open("stefania.dat"); // fittando da 0.5microsec fino a 14 2exp + bin 1/1 + mixed calib stef. OK 
   //   in.open("stefania.dat"); // fittando da 2microsec fino a 14 1exp + bin 1/1 (or 2/1, 4/1)OK 
   //   in.open("stefania.dat"); // fittando da 0.5microsec fino a 14 2exp + bin 1/1 (calib stef. <1us)OK 
   //   in.open("sabrina.dat"); // fittando da 2microsec fino a 14 1exp + bin 1/1 (or 2/1, 4/1)OK 
   while (1) {
      in >> Entries;
      if (!in.good()) break;
      if (nlines < 10) printf("Reading from input file \n line = ch = %d Entries=%8f, \n",nlines, Entries);
          histo->Fill(nlines,Entries);
	  time->Fill(t_calibrate(nlines),Entries);
	  if (nlines < 10) printf("Calibrated time, Entries=%8f, %8f \n",t_calibrate(nlines),Entries);
	  nlines++;
   }
   printf(" found %d points\n",nlines);
   in.close();
   Stat_t chmax = time->GetBinContent(2);
   time->SetMaximum(chmax*(1.3));
   

//   c1->Divide(2,2);

//   c1->cd(1);
//   c1->SetLogy(1);
// create a background function and fit the long range 
//   TF1 *expoFcn = new TF1("expoFcn",background,600,limit,1);
   TF1 *expoFcn = new TF1("expoFcn",background,10.0,15.0,1);
   expoFcn->SetNpx(500);
   expoFcn->SetLineWidth(2);
   expoFcn->SetLineColor(kRed);
//   expoFcn->SetParameter(0,1.5); // noise level
     
   time->Fit("expoFcn","lr");
//   TH1F *h1=time;
//   h1->Draw();
//   time->Draw("e1");   

   time->Draw("e1");
   //c1->Modified();
   //c1->Update();
   printf(" First fit done background\n");
   c1->WaitPrimitive();
   c1->cd(2);
   
   
//   c1->cd(2);
// create a exponential +backg fixed function and fit the long range 
//   TF1 *expobckg = new TF1("expobckg",expobckg,0,600,3);
   TF1 *expobckg = new TF1("expobckg",expobckg,lowl,highl,3);
   expobckg->SetNpx(500);
   expobckg->SetLineWidth(2);
   expobckg->SetLineColor(kBlue);
   expobckg->SetParameter(0,expoFcn->GetParameter(0)); // bckg level from previous fit
   expobckg->SetParLimits(0,expoFcn->GetParameter(0),expoFcn->GetParameter(0)); // fixes bckg level
   expobckg->SetParameter(1,20); // integral
   expobckg->SetParameter(2,2.2); // tau
     
   printf("now Fitting (1exp - bck fixed).....\n");
   time->Fit("expobckg","lr");
   //time->Draw("e1 same");   
   time->Draw("e1");   
   printf(" Second fit done - single exponential\n");
   //   c1->Modified();
   //   c1->Update();        
   c1->WaitPrimitive();
   c1->cd(3);
   
/*   
      for (int i=1;i!=limit;i++)
	     resid->Fill(i,(histo->GetBinContent(i)-expobckg->Eval((float)i)));
   c1->SetLogy(0);
    resid->Draw();
	c1->Modified();
   c1->Update(); 
	 c1->WaitPrimitive();
      c1->SetLogy(1);
*/   
//   c1->cd(3);
// create a 2 exponential +backg function and fit the whole range 
//   TF1 *fitfun = new TF1("fitfun",fitFunction,5,600,5);
   
   TF1 *fitfun = new TF1("fitfun",fitFunction,lowl,highl,5);
   fitfun->SetNpx(500);
   fitfun->SetLineWidth(2);
   fitfun->SetLineColor(kGreen);
   fitfun->SetLineStyle(1);
   fitfun->SetParameter(0,expobckg->GetParameter(0)); // bckg level
   fitfun->SetParLimits(0,expobckg->GetParameter(0),expobckg->GetParameter(0)); // fixes bckg level
   fitfun->SetParameter(1,expobckg->GetParameter(1)); // integral
   fitfun->SetParameter(2,expobckg->GetParameter(2)); // tau
   fitfun->SetParLimits(2,expobckg->GetParameter(2),expobckg->GetParameter(2)); // fixes previously fitted exp.
   fitfun->SetParameter(3,expobckg->GetParameter(1)); // integral
   fitfun->SetParameter(4,0.55); // tau
   printf("now Fitting  (2exp, bck and first-tau fixed).....\n");
   time->Fit("fitfun","lr");
   time->Draw("e1");   
   //   c1->Modified();
   //   c1->Update();        
   c1->WaitPrimitive();
   c1->cd(4);
   
   //
   TF1 *fit2exp = new TF1("fit2exp",fitFunction,lowl,highl,5);
   fit2exp->SetNpx(500);
   fit2exp->SetLineWidth(3);
   fit2exp->SetLineColor(kGreen);
   fit2exp->SetLineStyle(2);
   fit2exp->SetParameter(0,expobckg->GetParameter(0)); // bckg level
   fit2exp->SetParLimits(0,expobckg->GetParameter(0),expobckg->GetParameter(0)); // fixes bckg level
   fit2exp->SetParameter(1,expobckg->GetParameter(1)); // integral
   fit2exp->SetParameter(2,expobckg->GetParameter(2)); // tau
   fit2exp->SetParameter(3,expobckg->GetParameter(1)); // integral
   fit2exp->SetParameter(4,0.55); // tau
   printf("now Fitting (2exp., only bkg fixed).....\n");
   time->Fit("fit2exp","lr");
   time->Draw("e1");
   c1->Modified();
   c1->Update();        

  gBenchmark->Show("canvas");
/*
//   c1->Modified();
   c1->Update(); 
   c1->WaitPrimitive();
   c1->cd(4);        
   // create a 2 exponential function and fit the whole range 
   TF1 *twoexpo = new TF1("twoexpo",twoexpo,0,10,4);
   twoexpo->SetNpx(500);
   twoexpo->SetLineWidth(1);
   twoexpo->SetLineColor(kYellow);
   twoexpo->SetParameter(0,fitfun->GetParameter(1)/2); //integral
   twoexpo->SetParameter(1,fitfun->GetParameter(2)); // tau
   twoexpo->SetParameter(2,110.0); // integral
   twoexpo->SetParameter(3,2.2); // tau
   printf("now Fitting.....\n");
   time->Fit("twoexpo","lr");
//   time->Draw("e1");  
*/
}