//
// Theoretical curves showing ttbar xs dependence vs top mass from:
//
// 1) M. Cacciari et al  JHEP 0404:068 (2004)
//    see also PLB 378 (1996) p329 for the functional relation xs(TopMass).
//
// 2) N. Kidonakis, R. Vogt Phys. Rev. D 68, 114014 (2003).


{
gROOT->Reset();


Int_t MAXRESULTS = 30;
Float_t analysisIndex[MAXRESULTS];
Float_t ZeroArray[MAXRESULTS];
Float_t SmallArray[MAXRESULTS];

TString analysisDescription[MAXRESULTS];
Float_t mass[MAXRESULTS];
Float_t xsecValue[MAXRESULTS];
Float_t statErrorPos[MAXRESULTS];
Float_t statErrorNeg[MAXRESULTS];
Float_t systErrorPos[MAXRESULTS];
Float_t systErrorNeg[MAXRESULTS];
Float_t statSystErrorPos[MAXRESULTS];
Float_t statSystErrorNeg[MAXRESULTS];



for (Int_t j = 0; j < MAXRESULTS; j++) {
  analysisIndex[j] = j + 1;
  ZeroArray[j] = 0.;
  SmallArray[j] = 0.1;
}

Int_t i = -1;

//Results here

//200pb

/*
i++;
analysisDescription[i] = "";
mass[i] = 170.0;
xsecValue[i] = 6.24;
statErrorPos[i] = 0.90;
statErrorNeg[i] = 0.90;
statSystErrorPos[i] = 1.20;
statSystErrorNeg[i] = 1.20;

i++;
analysisDescription[i] = "";
mass[i] = 175.0;
xsecValue[i] = 6.02;
statErrorPos[i] = 0.90;
statErrorNeg[i] = 0.90;
statSystErrorPos[i] = 1.18;
statSystErrorNeg[i] = 1.18;

i++;
analysisDescription[i] = "";
mass[i] = 180.0;
xsecValue[i] = 5.82;
statErrorPos[i] = 0.90;
statErrorNeg[i] = 0.90;
statSystErrorPos[i] = 1.16;
statSystErrorNeg[i] = 1.16;

*/
/*
//350pb-1 results
i++;
analysisDescription[i] = "";
mass[i] = 170.0;
xsecValue[i] = 7.44;
statErrorPos[i] = 0.60;
statErrorNeg[i] = 0.60;
statSystErrorPos[i] = 1.08;
statSystErrorNeg[i] = 1.08;

i++;
analysisDescription[i] = "";
mass[i] = 175.0;
xsecValue[i] = 7.14;
statErrorPos[i] = 0.60;
statErrorNeg[i] = 0.60;
statSystErrorPos[i] = 1.05;
statSystErrorNeg[i] = 1.05;


i++;
analysisDescription[i] = "";
mass[i] = 180.0;
xsecValue[i] = 6.70;
statErrorPos[i] = 0.60;
statErrorNeg[i] = 0.60;
statSystErrorPos[i] = 1.00;
statSystErrorNeg[i] = 1.00;
*/


//750pb-1 results
i++;
analysisDescription[i] = "";
mass[i] = 170.0;
xsecValue[i] = 7.68;
statErrorPos[i] = 0.47;
statErrorNeg[i] = 0.47;
statSystErrorPos[i] = 0.88;
statSystErrorNeg[i] = 0.88;

i++;
analysisDescription[i] = "";
mass[i] = 175.0;
xsecValue[i] = 7.32;
statErrorPos[i] = 0.47;
statErrorNeg[i] = 0.47;
statSystErrorPos[i] = 0.85;
statSystErrorNeg[i] = 0.85;


i++;
analysisDescription[i] = "";
mass[i] = 180.0;
xsecValue[i] = 6.95;
statErrorPos[i] = 0.47;
statErrorNeg[i] = 0.47;
statSystErrorPos[i] = 0.82;
statSystErrorNeg[i] = 0.82;

//////////////////////////////////




//Fit mass dependence of combined results
grStatSyst = new TGraphAsymmErrors(i+1, mass, xsecValue, 
				   //ZeroArray, ZeroArray, statSystErrorNeg, statSystErrorPos);
			       ZeroArray, ZeroArray, SmallArray, SmallArray);

grStatSyst->SetTitle("Results (stat. +syst. errors)");

TF1 *func= new TF1("func","[0]+(175.0-x)*[1]",160.0, 190.0);
func->SetParameters(0,6.70);
func->SetParameters(1,0.0);
func->SetLineColor(6);
grStatSyst->Fit("func");
Double_t par1[2];
func->GetParameters(par1);

 TF1 *xscombo= new TF1("xscombo","[0]+(175.0-x)*[1]",169.3, 174.7); //Top Mass range
xscombo->SetParameter(0,par1[0]);
xscombo->SetParameter(1,par1[1]);

 TF1 *xscombostat= new TF1("xscombostat","[0]+(175.0-x)*[1]",170.4, 173.6); //Top Mass range
xscombostat->SetParameter(0,par1[0]);
xscombostat->SetParameter(1,par1[1]);


grErr = new TGraphAsymmErrors(i+1, mass, statSystErrorPos, 
				   //ZeroArray, ZeroArray, statSystErrorNeg, statSystErrorPos);
			       ZeroArray, ZeroArray, SmallArray, SmallArray);


grErr->Fit("func");
Double_t par2[2];
func->GetParameters(par2);

 TF1 *xserr= new TF1("xserr","[0]+(175.0-x)*[1]",169.3, 174.7); //Top Mass range
xserr->SetParameter(0,par2[0]);
xserr->SetParameter(1,par2[1]);

//Draw it!
gROOT->SetStyle("Plain");
gStyle->SetOptStat(0);
TCanvas *myCanvas = new TCanvas("myCanvas","Dependence of Combined result on top mass", 700, 500);
myCanvas->SetTicks();
myCanvas->SetBottomMargin(0.15);

TH2F *nullTH2F = new TH2F("nullTH2F", "", 20, 160., 180., 20, 0., 12.);
gStyle->SetTitleFont(22);
nullTH2F->GetXaxis()->SetTitle("Top Quark Mass (GeV/c^{2})");
nullTH2F->GetYaxis()->SetTitle("#sigma(p#bar{p} #rightarrow t#bar{t}) (pb)");
//nullTH2F->GetXaxis()->CenterTitle();
nullTH2F->Draw();


//All the theory curves - Kidonakis first

// These two lines should go through the following points:
// mt=175 syst low:  (1800., 4.51), (1960., 5.82) => y=0.0082x - 10.25
// mt=175 syst high: (1800., 5.71), (1960., 7.41) => y=0.0106x - 13.37
Float_t theoryMass[7] = {
 //1.5000E+02,
 //1.5500E+02,
 1.6000E+02,
 1.6500E+02,
 1.7000E+02,
 1.7500E+02,
 1.8000E+02,
 1.8500E+02,
 1.9000E+02//,
 //1.9500E+02,
 //200.0
 };

Float_t theoryNLO[7] = {
//1.4863E+01,
//1.2508E+01,
1.0570E+01,
8.9653E+00,
7.6306E+00,
6.5151E+00,
5.5786E+00,
4.7893E+00,
4.1216E+00,
//3.5548E+00,
//3.0721E+00
};

Float_t theoryNLO_PIM[7] = {
//1.4130E+01,
//1.1965E+01,
1.0166E+01,
8.6638E+00,
7.4058E+00,
6.3473E+00,
5.4534E+00,
4.6958E+00,
4.0523E+00,
//3.5037E+00,
//3.0348E+00
};

Float_t theoryNLO_1PI[7] = {
//1.6455E+01,
//1.3831E+01,
1.1674E+01,
9.8922E+00,
8.4114E+00,
7.1754E+00,
6.1392E+00,
5.2669E+00,
4.5297E+00,
//3.9046E+00,
//3.3727E+00
};

Float_t theoryNLO_AVG[7] = {
//1.5292E+01,
//1.2898E+01,
1.0920E+01,
9.2780E+00,
7.9086E+00,
6.7614E+00,
5.7963E+00,
4.9814E+00,
4.2910E+00,
//3.7042E+00,
//3.2038E+00
};

TGraph *theoryPlotNLO = new TGraph(7,theoryMass,theoryNLO);
theoryPlotNLO->SetLineStyle(1);
//theoryPlotNLO->Draw();

TGraph *theoryPlotNLO_PIM = new TGraph(7,theoryMass,theoryNLO_PIM);
theoryPlotNLO_PIM->SetLineStyle(4);
theoryPlotNLO_PIM->SetLineColor(1);
theoryPlotNLO_PIM->SetLineWidth(2.0);
theoryPlotNLO_PIM->Draw();

TGraph *theoryPlotNLO_1PI = new TGraph(7,theoryMass,theoryNLO_1PI);
theoryPlotNLO_1PI->SetLineStyle(3);
theoryPlotNLO_1PI->SetLineColor(1);
theoryPlotNLO_1PI->SetLineWidth(2.0);
theoryPlotNLO_1PI->Draw();

TGraph *theoryPlotNLO_AVG = new TGraph(7,theoryMass,theoryNLO_AVG);
theoryPlotNLO_AVG->SetLineStyle(1);
//theoryPlotNLO_AVG->Draw();

TF1 *xstheoryvogt= new TF1("xstheoryvogt","6.77*exp((175.0-x)/32.19)",160.0, 190.0);

//
//All the theory curves - Cacciari next
//
//

TF1 *xstheory= new TF1("xstheory","6.70*exp((175.0-x)/32.35)",160.0, 190.0);
//TF1 *xstheory= new TF1("xstheory","[0]*exp((175.0-x)/[1])",160.0, 190.0);
//xstheory->SetParameter(0,par1[0]);
//xstheory->SetParameter(1,par1[1]);
xstheory->SetLineColor(1);
xstheory->SetLineWidth(2.0);

TF1 *xstheorymin= new TF1("xstheorymin","5.82*exp((175.0-x)/32.65)",160.0, 190.0);
//TF1 *xstheorymin= new TF1("xstheorymin","[0]*exp((175.0-x)/[1])",160.0, 190.0);
//xstheorymin->SetParameter(0,par2[0]);
//xstheorymin->SetParameter(1,par2[1]);
xstheorymin->SetLineColor(1);
xstheorymin->SetLineStyle(2);
xstheorymin->SetLineWidth(2.0);

TF1 *xstheorymax= new TF1("xstheorymax","7.42*exp((175.0-x)/31.68)",160.0, 190.0);
//TF1 *xstheorymax= new TF1("xstheorymax","[0]*exp((175.0-x)/[1])",160.0, 190.0);
//xstheorymax->SetParameter(0,par3[0]);
//xstheorymax->SetParameter(1,par3[1]);
xstheorymax->SetLineColor(1);
xstheorymax->SetLineStyle(2);
xstheorymax->SetLineWidth(2.0);

//Draw these ...
xstheorymax->Draw("same");
xstheory->Draw("same");
xstheorymin->Draw("same");  


TLegend *tleg = new TLegend(0.35,0.2,0.88,0.38);
tleg->SetFillColor(0);
tleg->SetBorderSize(1);
tleg->AddEntry(xstheory, "Cacciari et al. JHEP 0404:068 (2004)", "l");
tleg->AddEntry(xstheorymax, "Cacciari et al. #pm uncertainty", "l");
tleg->AddEntry(theoryPlotNLO_PIM, "Kidonakis,Vogt PIM PRD 68 114014 (2003)", "l");
tleg->AddEntry(theoryPlotNLO_1PI, "Kidonakis,Vogt 1PI", "l");
tleg->Draw();

//////////////////////////////////////////////////////////////////////

//Draw result for world average mass, with mass dependence for +- 1 std dev of world average

Float_t massa[MAXRESULTS];
Float_t xseca[MAXRESULTS];
Float_t posa[MAXRESULTS];
Float_t nega[MAXRESULTS];
Float_t stata[MAXRESULTS];

 massa[0] = 172.0;
 xseca[0] = xscombo->Eval(massa[0]);
 posa[0]  = xserr->Eval(massa[0]);
 nega[0] = posa[0];
 stata[0] = statErrorPos[0];
 grComboStatSyst = new TGraphAsymmErrors(1, massa, xseca, ZeroArray, ZeroArray, posa, nega);
 grComboStat = new TGraphAsymmErrors(1, massa, xseca, ZeroArray, ZeroArray, stata, stata);


xscombo->SetLineColor(4);
xscombo->SetLineWidth(3);
xscombo->Draw("same");

xscombostat->SetLineColor(2);
xscombostat->SetLineWidth(3);
xscombostat->Draw("same");


grComboStatSyst->SetLineColor(4);
grComboStatSyst->SetLineWidth(3);
grComboStatSyst->SetMarkerSize(1.6);
grComboStatSyst->SetMarkerColor(4);
grComboStatSyst->SetMarkerStyle(20);
grComboStatSyst->Draw("P");

grComboStat->SetLineColor(2);
grComboStat->SetLineWidth(3);
grComboStat->SetMarkerSize(1.6);
grComboStat->SetMarkerColor(2);
grComboStat->SetMarkerStyle(20);
grComboStat->Draw("P");






TLatex *legendLabel = new TLatex();
legendLabel->SetTextSize(0.055);
legendLabel->SetTextColor(2);
legendLabel->DrawLatex(168.0, 10.5, "CDF II Preliminary 760 pb^{-1}");

 std::cout << " xs = " << xseca[0] << " +- " << posa[0] << " pb " << std::endl;

 // TString analysisResult = "";
 //char cv[20];
 //analysisResult += "Combined #sigma(p#bar{p} #rightarrow t#bar{t})=" ;
 //sprintf(cv,"%4.1f",xseca[0]);
 //analysisResult += cv;
 //analysisResult +=  "#pm";
 //sprintf(cv,"%4.1f",posa[0]);
 //analysisResult += cv;
 //analysisResult +=  " pb at m_{top} = 172.7 #pm 2.9 GeV/c^{2}";
 //legendLabel->DrawLatex(165.0, 2.0, analysisResult);




 myCanvas->Update();
/*Sanity check for fitted line

grStatSyst->SetLineColor(4);
grStatSyst->SetLineWidth(2);
grStatSyst->SetMarkerSize(1.4);
grStatSyst->SetMarkerColor(4);
grStatSyst_Ht->SetLineColor(4);
grStatSyst_Ht->SetLineWidth(2);
grStatSyst_Ht->SetMarkerSize(1.4);
grStatSyst_Ht->SetMarkerColor(4);
grStatSyst_Ht->SetMarkerStyle(22);
grStatSyst->Draw("P");
grStatSyst_Ht->Draw("P");
*/

}