//
// 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];

TString analysisDescription_Ht[MAXRESULTS];
Float_t mass_Ht[MAXRESULTS];
Float_t xsecValue_Ht[MAXRESULTS];
Float_t statErrorPos_Ht[MAXRESULTS];
Float_t statErrorNeg_Ht[MAXRESULTS];
Float_t systErrorPos_Ht[MAXRESULTS];
Float_t systErrorNeg_Ht[MAXRESULTS];
Float_t statSystErrorPos_Ht[MAXRESULTS];
Float_t statSystErrorNeg_Ht[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

i++;
analysisDescription[i] = "ANN";
mass[i] = 160.0;
xsecValue[i] = 7.93;
statErrorPos[i] = 1.33;
statErrorNeg[i] = 1.33;
systErrorPos[i] = 1.52;
systErrorNeg[i] = 1.52;

i++;
analysisDescription[i] = "";
mass[i] = 165.0;
xsecValue[i] = 7.48;
statErrorPos[i] = 1.27;
statErrorNeg[i] = 1.27;
systErrorPos[i] = 1.52;
systErrorNeg[i] = 1.52;

i++;
analysisDescription[i] = "";
mass[i] = 170.0;
xsecValue[i] = 6.98;
statErrorPos[i] = 1.19;
statErrorNeg[i] = 1.19;
systErrorPos[i] = 1.52;
systErrorNeg[i] = 1.52;

i++;
analysisDescription[i] = "";
mass[i] = 175.0;
xsecValue[i] = 6.64;
statErrorPos[i] = 1.13;
statErrorNeg[i] = 1.13;
systErrorPos[i] = 1.52;
systErrorNeg[i] = 1.52;

i++;
analysisDescription[i] = "";
mass[i] = 180.0;
xsecValue[i] = 6.25;
statErrorPos[i] = 1.07;
statErrorNeg[i] = 1.07;
systErrorPos[i] = 1.52;
systErrorNeg[i] = 1.52;

i++;
analysisDescription[i] = "";
mass[i] = 185.0;
xsecValue[i] = 5.86;
statErrorPos[i] = 1.03;
statErrorNeg[i] = 1.03;
systErrorPos[i] = 1.52;
systErrorNeg[i] = 1.52;

i++;
analysisDescription[i] = "";
mass[i] = 190.0;
xsecValue[i] = 5.67;
statErrorPos[i] = 1.0;
statErrorNeg[i] = 1.0;
systErrorPos[i] = 1.52;
systErrorNeg[i] = 1.52;


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

Int_t l = -1;

l++;
analysisDescription_Ht[l] = "Ht";
mass_Ht[l] = 160.0;
xsecValue_Ht[l] = 5.17;
statErrorPos_Ht[l] = 2.09;
statErrorNeg_Ht[l] = 2.09;
systErrorPos_Ht[l] = 1.81;
systErrorNeg_Ht[l] = 1.81;

l++;
analysisDescription_Ht[l] = "";
mass_Ht[l] = 170.0;
xsecValue_Ht[l] = 4.85;
statErrorPos_Ht[l] = 1.72;
statErrorNeg_Ht[l] = 1.72;
systErrorPos_Ht[l] = 1.81;
systErrorNeg_Ht[l] = 1.81;


l++;
analysisDescription_Ht[l] = "";
mass_Ht[l] = 165.0;
xsecValue_Ht[l] = 5.11;
statErrorPos_Ht[l] = 1.89;
statErrorNeg_Ht[l] = 1.89;
systErrorPos_Ht[l] = 1.81;
systErrorNeg_Ht[l] = 1.81;

l++;
analysisDescription_Ht[l] = "";
mass_Ht[l] = 175.0;
xsecValue_Ht[l] = 4.75;
statErrorPos_Ht[l] = 1.62;
statErrorNeg_Ht[l] = 1.62;
systErrorPos_Ht[l] = 1.81;
systErrorNeg_Ht[l] = 1.81;

l++;
analysisDescription_Ht[l] = "";
mass_Ht[l] = 180.0;
xsecValue_Ht[l] = 4.5;
statErrorPos_Ht[l] = 1.42;
statErrorNeg_Ht[l] = 1.42;
systErrorPos_Ht[l] = 1.81;
systErrorNeg_Ht[l] = 1.81;

l++;
analysisDescription_Ht[l] = "";
mass_Ht[l] = 185.0;
xsecValue_Ht[l] = 4.38;
statErrorPos_Ht[l] = 1.31;
statErrorNeg_Ht[l] = 1.31;
systErrorPos_Ht[l] = 1.81;
systErrorNeg_Ht[l] = 1.81;

l++;
analysisDescription_Ht[l] = "";
mass_Ht[l] = 190.0;
xsecValue_Ht[l] = 4.18;
statErrorPos_Ht[l] = 1.22;
statErrorNeg_Ht[l] = 1.22;
systErrorPos_Ht[l] = 1.81;
systErrorNeg_Ht[l] = 1.81;


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

for (Int_t j = 0; j < MAXRESULTS; j++) {
  statSystErrorNeg[j] = sqrt(statErrorNeg[j]**2 + systErrorNeg[j]**2);
  statSystErrorPos[j] = sqrt(statErrorPos[j]**2 + systErrorPos[j]**2);
  statSystErrorNeg_Ht[j] = sqrt(statErrorNeg_Ht[j]**2 + systErrorNeg_Ht[j]**2);
  statSystErrorPos_Ht[j] = sqrt(statErrorPos_Ht[j]**2 + systErrorPos_Ht[j]**2);
}



//Fit mass dependence of ANN 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 *xsann= new TF1("xsann","[0]+(175.0-x)*[1]",169.4, 186.6);
xsann->SetParameter(0,par1[0]);
xsann->SetParameter(1,par1[1]);

//Fit mass dependence of HT results

grStatSyst_Ht = new TGraphAsymmErrors(l+1, mass_Ht, xsecValue_Ht, 
			     //ZeroArray, ZeroArray, statSystErrorNeg_Ht, statSystErrorPos_Ht);
			       ZeroArray, ZeroArray, SmallArray, SmallArray);
grStatSyst_Ht->SetTitle("Results (stat.+syst. errors)");

grStatSyst_Ht->Fit("func"); 
Double_t par2[2];
func->GetParameters(par2);
TF1 *xs_Ht= new TF1("xs_Ht","[0]+(175.0-x)*[1]",169.4, 186.6);
xs_Ht->SetParameter(0,par2[0]);
xs_Ht->SetParameter(1,par2[1]);



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

TH2F *nullTH2F = new TH2F("nullTH2F", "", 20, 160., 190., 20, 0., 13.);
gStyle->SetTitleFont(22);
nullTH2F->GetXaxis()->SetTitle("Top 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();

//
//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.7,0.88,0.88);
tleg->SetFillColor(0);
tleg->SetBorderSize(1);
tleg->AddEntry(xstheory, "Cacciari et al. JHEP 0404:068 (2004)", "l");
tleg->AddEntry(xstheorymax, "Cacciari et al. #pm 1#sigma", "l");
tleg->AddEntry(theoryPlotNLO_PIM, "Kidonakis,Vogt PIM PRD 68 114014 (2003)", "l");
tleg->AddEntry(theoryPlotNLO_1PI, "Kidonakis,Vogt 1PI", "l");
tleg->Draw();

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

//Draw 175 results, with mass dependence for +- 2 std dev of world average
//Offset them a little, so that errors don't overlap

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

 massa[0] = mass[3]+0.25;
 xseca[0] = xsecValue[3];
 nega[0]  = statSystErrorNeg[3];
 posa[0]  = statSystErrorPos[3];
 grAnn = new TGraphAsymmErrors(1, massa, xseca, ZeroArray, ZeroArray, nega, posa);
 nega[0]  = statErrorNeg[3];
 posa[0]  = statErrorPos[3];
 grAnnStat = new TGraphAsymmErrors(1, massa, xseca, ZeroArray, ZeroArray, nega, posa);


 massa[0] = mass_Ht[3]-0.25;
 xseca[0] = xsecValue_Ht[3];
 nega[0]  = statSystErrorNeg_Ht[3];
 posa[0]  = statSystErrorPos_Ht[3];
 gr_Ht = new TGraphAsymmErrors(1, massa, xseca, ZeroArray, ZeroArray, nega, posa);
 nega[0]  = statErrorNeg_Ht[3];
 posa[0]  = statErrorPos_Ht[3];
 gr_HtStat = new TGraphAsymmErrors(1, massa, xseca, ZeroArray, ZeroArray, nega, posa);


grAnn->SetLineColor(4);
grAnn->SetLineWidth(3);
grAnn->SetMarkerSize(1.6);
grAnn->SetMarkerColor(2);
grAnn->SetMarkerStyle(20);
grAnnStat->SetLineColor(2);
grAnnStat->SetLineWidth(3);
grAnnStat->SetMarkerSize(1.6);
grAnnStat->SetMarkerColor(2);
grAnnStat->SetMarkerStyle(20);
grAnn->Draw("P");
xsann->SetLineColor(4);
xsann->SetLineWidth(3);
xsann->Draw("same");
grAnnStat->Draw("P");

gr_Ht->SetLineColor(4);
gr_Ht->SetLineWidth(3);
gr_Ht->SetMarkerSize(1.6);
gr_Ht->SetMarkerColor(2);
gr_Ht->SetMarkerStyle(22);
gr_HtStat->SetLineColor(2);
gr_HtStat->SetLineWidth(3);
gr_HtStat->SetMarkerSize(1.6);
gr_HtStat->SetMarkerColor(2);
gr_HtStat->SetMarkerStyle(22);
gr_Ht->Draw("P");
xs_Ht->SetLineColor(4);
xs_Ht->SetLineWidth(3);
xs_Ht->Draw("same");
gr_HtStat->Draw("P");

TLegend *tleg = new TLegend(0.15,0.20,0.27,0.30);
tleg->SetFillColor(0);
tleg->SetBorderSize(1);
tleg->AddEntry(grAnn, "ANN fit", "p");
tleg->AddEntry(gr_Ht, "Ht fit", "p");
tleg->Draw();

/*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");
*/

}