//*************************************************************************************
//    START fitting MACROS 
//**************************************************************************************

  TString one="1";                       /* do two rapidity ranges here*/
  TString two="2";

tfitter4(Int_t sample,TString raprange){ 

gStyle->SetNdivisions(505,"Y");//as per Horner's suggestion
gStyle->SetNdivisions(505,"X");
gStyle->SetPadBorderMode(0);
gStyle->SetPadBorderSize(0);
gStyle->SetPadColor(kWhite);
gStyle->SetPadTickX(1);
gStyle->SetPadTickY(1);
gROOT->SetStyle("Plain");

//define variables here
  TString tplot = "tp"+raprange;
  cout << "Tplot = "<<tplot << endl;
  TString MCtplot = "tMC"+raprange;
  cout << "MCTplot = "<<MCtplot << endl;
  TString tbgplot = "tbg"+raprange;
  TString tfitseps  = "tfits"+raprange+".eps";
  //TString tfitsgif  = "tfits"+raprange+".gif";
  TString thistseps = "thists"+raprange+".eps";
  TString thistsgif = "thists"+raprange+".eps";

  TString tlabel = "XnXn t for 0 < |y| < 0.5";
  if (sample ==1 && raprange == "2") tlabel = "XnXn t for 0.5 < |y| < 1.0";
  if (sample ==2 && raprange == "1") tlabel = "Topo t for 0.05 < |y| < 0.5";
  if (sample ==2 && raprange == "2") tlabel = "Topo t for 0.5 < |y| < 1.0";
  cout << "TLabel = "<<tlabel<<endl;

if (sample ==1) {
     TFile* DataFile =new TFile("mbhist.root","READ");
     TFile* IntFile =new TFile("xnxnint70hist.root","READ");
     TFile* NointFile =new TFile("xnxnnoint70hist.root","READ");
     TFile* IntFile2 =new TFile("xnxnint70hist2.root","READ");
     TFile* NointFile2 =new TFile("xnxnnoint70hist2.root","READ");
     cout << "XnXn files are open.  Starting fitting histograms"<<endl;
}

if (sample ==2) {
     TFile* DataFile =new TFile("topohist.root","READ");
     TFile* IntFile =new TFile("intzhist.root","READ");
     TFile* NointFile =new TFile("nointzhist.root","READ");
     TFile* IntFile2 =new TFile("intzhist2.root","READ");
     TFile* NointFile2 =new TFile("nointzhist2.root","READ");
     cout << "Topology files are open.  Starting fitting histograms"<<endl;
}

gStyle->SetStatW(0.2);
gStyle->SetStatH(0.2);

gStyle->SetOptStat(1111);
gStyle->SetOptFit(1111111);

//gStyle->SetOptStat(0);
//gStyle->SetOptFit(0);

// skip wide t plots

TCanvas* c228 = new TCanvas("tcomparisons",tlabel,0,0,600,400);
c228->cd(1);

//gPad->SetLogy(1);
/*gStyle->SetOptStat(0);
gStyle->SetOptFit(0);*/

gStyle->SetStatW(0.25);
gStyle->SetStatH(0.25);

//  dN/dt
cout << tplot << "tplot" << endl;
TH1F * thist = DataFile->Get(tplot);
thist->GetYaxis()->SetTitleOffset(0.85);
thist->GetYaxis()->SetTitleSize(0.06);
thist->GetXaxis()->SetTitleSize(0.055);
thist->GetXaxis()->SetTitleOffset(0.8);
thist->GetXaxis()->SetTitle("t ~ p_{T}^{2} (GeV^{2})");
thist->GetYaxis()->SetTitle("dN/dt (GeV^{-2})");
thist->Sumw2();
thist->SetMarkerColor(2);
thist->SetLineColor(2);

//  set axis range

//thist->SetMaximum(90.);
thist->Draw("EP");

TH1F * tbghist = DataFile->Get(tbgplot);
tbghist->DrawCopy("same");

// integral for normalization - use full spectrum

Double_t tintegral = thist->Integral();

TH1F * tinthist = IntFile->Get(tplot);
tinthist->SetLineColor(4);
//  Normalize total event count

tinthist->Sumw2();
Double_t tintintegral = tinthist->Integral();
tinthist->Scale(tintegral/tintintegral);

tinthist->DrawCopy("histsame");

TH1F * tnointhist = NointFile->Get(tplot);
tnointhist->SetLineColor(3);
tnointhist->Sumw2();
Double_t tnointintegral = tnointhist->Integral();
tnointhist->Scale((tintegral/tnointintegral));
tnointhist->DrawCopy("histSAME");

//  *******************************************************************
//Start Fitting - take ratio int/noint and use in fit

 tfits:

TString tlabel2=tlabel+" Fits";

TCanvas* c233 = new TCanvas("Ratio fits",tlabel2,0,0,900,800);
c233->Divide(2,3);
c233->cd(2);

// Compare STARlight Int and Noint; take ratio; fit ratio to polynomial

TH1F *  tSTARint = IntFile2->Get(MCtplot);
TH1F *  tSTARnoint = NointFile2->Get(MCtplot);

// Normalize so int and noint have the same area under the curve
// This is NOT the same as requiring that STARrat average to 1

Double_t Intint = tSTARint->Integral();
Double_t Nointint = tSTARnoint->Integral();
tSTARint->Scale(Nointint/Intint);

TH1F * STARrat = new TH1F("STARrat","Int/Noint STARlight",50,0.,.01);
//STARrat->Sumw2();
tSTARint->Sumw2();
tSTARnoint->Sumw2();

STARrat->Divide(tSTARint,tSTARnoint,1.,1.);
STARrat->GetYaxis()->SetTitleSize(0.06);
STARrat->GetYaxis()->SetTitleOffset(0.9);
STARrat->GetXaxis()->SetTitleSize(0.055);
STARrat->GetXaxis()->SetTitleOffset(0.8);
STARrat->GetXaxis()->SetTitle("t ~ p_{T}^{2} (GeV^{2})");
STARrat->GetYaxis()->SetTitle("Ratio");
STARrat->SetMarkerColor(2);
STARrat->SetLineColor(2);
STARrat->SetMaximum(1.5);

// fit STARrat and extract the parameters

//Double_t par[5];

//STARfunc = new TF1("STARfunc","pol6",0.,.01);
//gROOT->LoadMacro("sfun.C");
TF1 *STARfunc = new TF1("STARfunc","1++(1/(x+0.012))++(1/((x+0.012)*(x+0.012)))++(1/((x+0.012)*(x+0.012)*(x+0.012)))++(1/((x+0.012)*(x+0.012)*(x+0.012)*(x+0.012)))",0.0,0.01);
cout<<" here "<<endl;
STARfunc->SetLineColor(2);
STARfunc->SetLineWidth(1);
//STARfunc->FixParameter(0,0.012);
STARrat->Fit("STARfunc");

cout << "STARrat parameters" << endl;
STARfunc->Print();

Double_t STARpar[9];
Double_t STARerr[9];

cout << "STARlight Ratio = ";
  for (Int_t i=0;i<= 5; i++) {
    STARpar[i] = STARfunc->GetParameter(i);
    STARerr[i] = STARfunc->GetParError(i);
    cout << STARpar[i] << "+/-"<<STARerr[i]<< " x^"<<i<<"+";
  }
 cout<<endl;

STARrat->DrawCopy();

// Compare STARlight + GEANT, bfc Int and Noint; take ratio; fit ratio to polynomial

c233->cd(1);

TH1F *  tGEANTint = IntFile->Get(tplot);
TH1F *  tGEANTnoint = NointFile->Get(tplot);

//These also need to be normalized

Double_t intGint=tGEANTint->Integral();
Double_t intGnoint=tGEANTnoint->Integral();
tGEANTint->Scale(intGnoint/intGint);

TH1F * GEANTrat = new TH1F("GEANTrat","Int/Noint- GEANT",50,0.,.01);
//GEANTrat->Sumw2();
tGEANTint->Sumw2();
tGEANTnoint->Sumw2();

GEANTrat->Divide(tGEANTint,tGEANTnoint,1.,1.);
Double_t GEANTint = GEANTrat->Integral();

GEANTrat->Scale(GEANTint/GEANTint);
GEANTrat->SetMarkerColor(4);
GEANTrat->SetLineColor(4);
GEANTrat->SetMaximum(1.5);
GEANTrat->GetYaxis()->SetTitle("Ratio");
GEANTrat->GetYaxis()->SetTitleSize(0.06);
GEANTrat->GetYaxis()->SetTitleOffset(0.9);
GEANTrat->GetXaxis()->SetTitleSize(0.055);
GEANTrat->GetXaxis()->SetTitleOffset(0.8);
GEANTrat->GetXaxis()->SetTitle("t ~ p_{T}^{2} (GeV^{2})");

// fit GEANT ratio
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
THE GEANT RATIO IS FIT BELOW
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

TF1 *Gfunc = new TF1("Gfunc","1++(1/(x+0.012))++(1/((x+0.012)*(x+0.012)))++(1/((x+0.012)*(x+0.012)*(x+0.012)))++(1/((x+0.012)*(x+0.012)*(x+0.012)*(x+0.012)))",0.0,0.01);
Gfunc->SetLineColor(4);
Gfunc->SetLineWidth(1);
GEANTrat->Fit("Gfunc");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PARAMETERS ARE RECORDED FROM THE RATIO FIT TO BE APPLIED TO THE OVERALL FIT WHICH INCLUDES R(t)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Double_t Gpar[9];
Double_t Gerr[9];

cout << "GEANT Ratio = ";
  for (Int_t i=0;i<= 5; i++) {
    Gpar[i] = Gfunc->GetParameter(i);
    Gerr[i] = Gfunc->GetParError(i);
    cout << Gpar[i] << "+/-"<<Gerr[i]<< " x^"<<i<<"+";
  }
cout<<endl;

GEANTrat->Draw("same");

c233->cd(3);

// take the ratio of ratios

//TH1F * ratrat = new TH1F("ratrat","Ratio of ratios",50,0.,.01);
//ratrat->Divide(STARrat,GEANTrat,1.,1.);
//ratrat->Draw();

//  Calculate efficiencies, with and without interference; with goes first

TH1F * tintgeant = IntFile->Get(tplot);
 tintgeant->Sumw2();
 TH1F * tintMC = IntFile2->Get(MCtplot);
 tintMC->Sumw2();

TH1F * inteff = new TH1F("Inteff"," Efficiency",50.,0.,.01);
 inteff->Divide(tintgeant,tintMC,1.0,1.0);
 inteff->SetLineColor(4);
 // inteff->Sumw2();
 inteff->GetYaxis()->SetTitleSize(0.06);
 inteff->GetYaxis()->SetTitleOffset(0.9);
 inteff->GetXaxis()->SetTitleSize(0.055);
 inteff->GetXaxis()->SetTitleOffset(0.8);
 inteff->GetXaxis()->SetTitle("t ~ p_{T}^{2} (GeV^{2})");
 inteff->DrawCopy("hist");

// now for no interference

TH1F * tnointgeant = NointFile->Get(tplot);
TH1F * tnointMC = NointFile2->Get(MCtplot);

 tnointgeant->Sumw2();
 tnointMC->Sumw2();

TH1F * nointeff = new TH1F("Nointeff"," Efficiency",50.,0.,.01);
nointeff->Divide(tnointgeant,tnointMC,1.0,1.0);
nointeff->SetLineColor(3);
//nointeff->Sumw2();
nointeff->DrawCopy("samehist");

TLegend * legend = new TLegend(0.15,0.80,0.65,0.95);
legend->SetFillColor(0);
legend->AddEntry(nointeff,"nointefficiency");
legend->AddEntry(inteff,"intefficiency");
legend->Draw();

c233->cd(4);

// now for fitting.  First raw data + corrected Monte Carlo  (+ backgrounds)

// define function with polynomial to 6th order, then fix parameters

//Gfix = new TF1("Gfix","[0]*exp(-[1]*x)*(1+[2]*([3]-1+[4]*x+[5]*x*x+[6]*x^3+[7]*x^4+[8]*x^5+[9]*x^6))");

TF1 *Gfix = new TF1("Gfix","[0]*exp(-[1]*x)*(1+[2]*([3] -1 + [4]/(x+0.012) + [5]/((x+0.012)*(x+0.012)) + [6]/((x+0.012)*(x+0.012)*(x+0.012)) + [7]/((x+0.012)*(x+0.012)*(x+0.012)*(x+0.012))))");

// TF1 *Gfix = new TF1("Gfix","[0]*exp(-[1]*x)*(1+[2]*([4] -1 + [5]/(x+[3]) + [6]/((x+[3])*(x+[3])) + [7]/((x+[3])*(x+[3])*(x+[3])) + [8]/((x+[3])*(x+[3])*(x+0.012)*(x+[3]))))");

Gfix->FixParameter(3,Gpar[0]);
Gfix->FixParameter(4,Gpar[1]);
Gfix->FixParameter(5,Gpar[2]);
Gfix->FixParameter(6,Gpar[3]);
Gfix->FixParameter(7,Gpar[4]);
//Gfix->FixParameter(8,Gpar[5]);
//Gfix->FixParameter(9,Gpar[6]);

/*Gfix->SetParameter(1,250.);
Gfix->SetParameter(2,1.);*/

Gfix->SetParameter(0,7000);
Gfix->SetParameter(1,300.);
Gfix->SetParameter(2,1.);


Gfix->SetLineWidth(1);
Gfix->SetLineColor(2);

TH1F * Datahist = DataFile->Get(tplot);
Datahist->SetTitle("t (uncorrected)");
Datahist->Sumw2();
Int_t Nevts = Datahist->Integral();

Datahist->Fit("Gfix");
Datahist->Fit("Gfix");
Datahist->SetMarkerStyle(22);
Datahist->SetMarkerColor(2);
Datahist->SetMarkerSize(0.8);
Datahist->DrawCopy();

tinthist->DrawCopy("histsame");
tnointhist->DrawCopy("histsame");

// background

TH1F * bgds = DataFile->Get(tbgplot);
bgds->DrawCopy("same");
Int_t nbgds = bgds->Integral();

cout << "Number of events = "<<Nevts<<".  Number of Background = " <<nbgds<<endl;

c233->cd(5);

// first create an efficiency corrected t
// Datahist->Sumw2();
 
TH1F * tcorr = new TH1F("tcorr","t (Eff. Corr.)",50.,0.,.01);
tcorr->SetLineColor(2);
tcorr->SetMarkerStyle(22);
tcorr->SetMarkerColor(2);
tcorr->SetMarkerSize(0.8);
tcorr->Divide(Datahist,inteff,1.0.,1.0);
tcorr->GetYaxis()->SetTitleSize(0.06);
tcorr->GetYaxis()->SetTitleOffset(0.85);
tcorr->GetXaxis()->SetTitleSize(0.055);
tcorr->GetXaxis()->SetTitleOffset(0.8);
tcorr->GetXaxis()->SetTitle("t ~ p_{T}^{2} (GeV^{2})");
tcorr->GetYaxis()->SetTitle("dN/dt (GeV^{-2})");

// now fit the corrected data + raw STARlight, for comparison

//SLfix = new TF1("SLfix","[0]*exp(-[1]*x)*(1+[2]*([3]-1+[4]*x+[5]*x*x+[6]*x^3+[7]*x^4+[8]*x^5+[9]*x^6))");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
THE OVERALL FIT TO THE EFFICIENCY CORRECTED DATA IS SHOWN BELOW
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


TF1 *SLfix = new TF1("SLfix","[0]*exp(-[1]*x)*(1+[2]*([3] -1 + [4]/(x+0.012) + [5]/((x+0.012)*(x+0.012)) + [6]/((x+0.012)*(x+0.012)*(x+0.012)) + [7]/((x+0.012)*(x+0.012)*(x+0.012)*(x+0.012))))");

 //TF1 *SLfix = new TF1("SLfix","[0]*exp(-[1]*x)*(1+[2]*([4] -1 + [5]/(x+[3]) + [6]/((x+[3])*(x+[3])) + [7]/((x+[3])*(x+[3])*(x+[3])) + [8]/((x+[3])*(x+[3])*(x+[3])*(x+[3]))))");

SLfix->FixParameter(3,STARpar[0]);
SLfix->FixParameter(4,STARpar[1]);
SLfix->FixParameter(5,STARpar[2]);
SLfix->FixParameter(6,STARpar[3]);
SLfix->FixParameter(7,STARpar[4]);
//SLfix->FixParameter(8,STARpar[5]);
//SLfix->FixParameter(9,STARpar[6]);

cout << "STARrat parameters" << endl;
STARfunc->Print();

//new TCanvas();

cout << "SLfix parameters" << endl;
SLfix->Print();
SLfix->SetParameter(0,7000);
SLfix->SetParameter(1,300.);
SLfix->SetParameter(2,1.);

SLfix->SetLineWidth(1);
SLfix->SetLineColor(2);

tcorr->DrawCopy();
//Call fit twice, otherwise it fails sometimes with just initial parameters
tcorr->Fit("SLfix");
tcorr->Fit("SLfix");
//tSTARint->Fit("SLfix");

//tcorr->SetStats(0);
tcorr->SetTitle("Total Fit");
tcorr->GetYaxis()->SetTitleSize(0.06);
tcorr->GetYaxis()->SetTitleOffset(0.85);
tcorr->GetXaxis()->SetTitleSize(0.055);
tcorr->GetXaxis()->SetTitleOffset(0.8);
tcorr->GetYaxis()->SetTitle("dN/dt (GeV^{-2})");
tcorr->GetXaxis()->SetTitle("t ~ p_{T}^{2} (GeV^{2})");

/*get Covariance Matrix elements
cout << "covariance Matrix here " << endl;
Double_t matrix[3][3];
gMinuit->mnemat(&matrix[0][0],3);
for (Int_t i = 0; i < 3; i++)
        printf("%5.2f, %5.2f, %5.2f\n",matrix[i][0],matrix[i][1],matrix[i][2]);

TVirtualFitter *fitter = TVirtualFitter::GetFitter();
TMatrixDSym cov;
cov.Use(fitter->GetNumberTotalParameters(),fitter->GetCovarianceMatrix());
TVectorD var = TMatrixDDiag(cov);
printf("%5.2f, %5.2f, %5.2f %5.2f\n",var[0],var[1],var[2],var[3]);
*/

tSTARint->SetLineColor(4);
tSTARint->DrawCopy("histsame");
tSTARnoint->SetLineColor(3);
tSTARnoint->DrawCopy("histsame");

 nextstep:
//cout << "Writing files" << tfitseps << " and "<< tfitsgif <<endl;
cout << "Writing files" << tfitseps << endl;

 TString tlabel11=tlabel+"fit of uncorrected t_{MC} distribution";



 // using R(t) from detector uncorrected distributions 
TCanvas* c11 = new TCanvas("tcom11", tlabel11,0,0,600,400);
c11->cd(1);

tSTARint->Fit("SLfix");
tSTARint->SetLineColor(4);
tSTARint->DrawCopy("histsame");
tSTARnoint->SetLineColor(3);
tSTARnoint->DrawCopy("histsame");



// now let's get the R(t) from uncorrected distributions 

 if(raprange == "1")
   {
     TH1F * t_dst_i = IntFile->Get("tp1");
     TH1F * t_dst_n = NointFile->Get("tp1");
     t_dst_i->Sumw2();
     t_dst_n->Sumw2();
   }


 if(raprange == "2")
   {
     TH1F * t_dst_i = IntFile->Get("tp2");
     TH1F * t_dst_n = NointFile->Get("tp2");
     t_dst_i->Sumw2();
     t_dst_n->Sumw2();
   }


TH1F * tratio = new TH1F("tratio1","t ratio for R(t)",50,0.,.01);


tratio->Divide(t_dst_i,t_dst_n,1.,1.);

tratio->GetYaxis()->SetTitleSize(0.06);
tratio->GetYaxis()->SetTitleOffset(0.85);
tratio->GetXaxis()->SetTitleSize(0.055);
tratio->GetXaxis()->SetTitleOffset(0.8);
tratio->GetYaxis()->SetTitle("dN/dt (GeV^{-2})");
tratio->GetXaxis()->SetTitle("t ~ p_{T}^{2} (GeV^{2})");
tratio->SetMarkerColor(2);
tratio->SetLineColor(2);
tratio->SetMaximum(1.5);

TF1 *tratiofunc = new TF1("tratiofunc","1++(1/(x+0.012))++(1/((x+0.012)*(x+0.012)))++(1/((x+0.012)*(x+0.012)*(x+0.012)))++(1/((x+0.012)*(x+0.012)*(x+0.012)*(x+0.012)))",0.0,0.01);
TCanvas* c12 = new TCanvas("tcom12", tlabel11,0,0,450,850);
 c12->Divide(1,4);
c12->cd(1);
tratiofunc->SetLineColor(2);
tratiofunc->SetLineWidth(1);
tratio->Fit("tratiofunc");


Double_t tratiopar[9];
Double_t tratioerr[9];


  for (Int_t i=0;i<= 5; i++) {
    tratiopar[i] = tratiofunc->GetParameter(i);
    tratioerr[i] = tratiofunc->GetParError(i);
     }
 

TF1 *SLfixratio = new TF1("SLfixratio","[0]*exp(-[1]*x)*(1+[2]*([3] -1 + [4]/(x+0.012) + [5]/((x+0.012)*(x+0.012)) + [6]/((x+0.012)*(x+0.012)*(x+0.012)) + [7]/((x+0.012)*(x+0.012)*(x+0.012)*(x+0.012))))");


SLfixratio->FixParameter(3,STARpar[0]);
SLfixratio->FixParameter(4,STARpar[1]);
SLfixratio->FixParameter(5,STARpar[2]);
SLfixratio->FixParameter(6,STARpar[3]);
SLfixratio->FixParameter(7,STARpar[4]);

SLfixratio->SetParameter(0,7000);
SLfixratio->SetParameter(1,300.);
SLfixratio->SetParameter(2,1.);

SLfixratio->SetLineWidth(1);
SLfixratio->SetLineColor(2);
 c12->cd(2);

tcorr->DrawCopy();
//Call fit twice, otherwise it fails sometimes with just initial parameters
tcorr->Fit("SLfixratio");
tcorr->Fit("SLfixratio");

 c12->cd(3);
tSTARint->GetYaxis()->SetTitleSize(0.06);
tSTARint->GetYaxis()->SetTitleOffset(0.85);
tSTARint->GetXaxis()->SetTitleSize(0.055);
tSTARint->GetXaxis()->SetTitleOffset(0.8);
tSTARint->GetYaxis()->SetTitle("dN/dt (GeV^{-2})");
tSTARint->GetXaxis()->SetTitle("t ~ p_{T}^{2} (GeV^{2})");
tSTARint->DrawCopy();
tSTARint->Fit("SLfixratio");
tSTARint->Fit("SLfixratio");

 c12->cd(4);
t_dst_i->DrawCopy();
t_dst_i->Fit("SLfix");
t_dst_i->Fit("SLfix");

}