{
//This provides formulae for ttbar xs versus top mass

//3 mass points
Int_t n=3;


Double_t xRun1[n] = { 160., 175., 190.};
Double_t yRun1[n] = { 7.78, 4.75, 2.97};

TF1 *func= new TF1("func","[0]*exp((175.0-x)/[1])",160.0, 190.0);
func->SetParameters(0,6.70);
func->SetParameters(1,32.0);
func->SetLineColor(6);


//Numbers from Kidonakis and Vogt
Double_t x[n] = { 170., 175., 180.};
Double_t y1[n] = { 7.91, 6.76, 5.80}; //central
Double_t y2[n] = { 6.79, 5.82, 5.00}; //min MLM
Double_t y3[n] = { 8.69, 7.41, 6.34}; //max MLM

//Store fit function results
Double_t parRun1[2];
Double_t par1[2];
Double_t par2[2];
Double_t par3[2];

Double_t ex[n], ey[n];
for(int i=0; i<n;i++){
  ex[i]= 0.0001;
  ey[i]= 0.01;
}

TGraphErrors *grRun1 = new TGraphErrors(n,xRun1,yRun1,ex,ey);
TGraphErrors *gr1 = new TGraphErrors(n,x,y1,ex,ey);
TGraphErrors *gr2 = new TGraphErrors(n,x,y2,ex,ey);
TGraphErrors *gr3 = new TGraphErrors(n,x,y3,ex,ey);

//Fit each set of points and store params
grRun1->Fit("func"); 
func->GetParameters(parRun1);
gr1->Fit("func"); 
func->GetParameters(par1);
gr2->Fit("func"); 
func->GetParameters(par2);
gr3->Fit("func"); 
func->GetParameters(par3);


//Draw each set of points and fitted function
nslide(2,3);
p1->cd(1);
grRun1->SetMarkerStyle(21);
grRun1->SetTitle("Run 1");
grRun1->GetXaxis()->SetTitle("Top mass");
grRun1->GetYaxis()->SetTitle("Top pair xs");
grRun1->Draw("ALP");

p1->cd(2);
gr1->SetMarkerStyle(21);
gr1->SetTitle("Run 2 central");
gr1->GetXaxis()->SetTitle("Top mass");
gr1->GetYaxis()->SetTitle("Top pair xs");
gr1->Draw("ALP");

p1->cd(3);
gr2->SetMarkerStyle(21);
gr2->SetTitle("Min");
gr2->GetXaxis()->SetTitle("Top mass");
gr2->GetYaxis()->SetTitle("Top pair xs");
gr2->Draw("ALP");

p1->cd(4);
gr3->SetMarkerStyle(21);
gr3->SetTitle("Max");
gr3->GetXaxis()->SetTitle("Top mass");
gr3->GetYaxis()->SetTitle("Top pair xs");
gr3->Draw("ALP");

//Some printout
 
std::cout << "Run I Central formula is "<< parRun1[0] << " *exp((175.0-mtop)/" << parRun1[1] << ")"<< std::endl;
std::cout << "Run II Min     formula is "<< par2[0] << " *exp((175.0-mtop)/" << par2[1] << ")"<< std::endl;
std::cout << "Run II Central formula is "<< par1[0] << " *exp((175.0-mtop)/" << par1[1] << ")"<< std::endl;
std::cout << "Run II Max     formula is "<< par3[0] << " *exp((175.0-mtop)/" << par3[1] << ")"<< std::endl;

double mtop=175;
double xs = 0.0;
double xs2 = 0.0;
double xs3 = 0.0;


std::cout << "mtop  xs: min      central   max " << std::endl;  
for(int i=0;i<9;i++){
  mtop = 160.0 + 5.0*i;
  xs = par1[0]*exp((175.0-mtop)/par1[1]);
  xs2 = par2[0]*exp((175.0-mtop)/par2[1]);
  xs3 = par3[0]*exp((175.0-mtop)/par3[1]);
  std::cout << mtop << "       " <<  xs2 << "  " << xs << "  " << xs3 << std::endl;  
}


//Now make functions corresponding to RUN II min, central and max xs
p1->cd(5);
//Vogt
TF1 *xstheory= new TF1("xstheory","6.77*exp((175.0-x)/32.19)",160.0, 190.0);
//Mangano
//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(4);

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(3);
xstheorymin->SetLineStyle(2);

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(3);
xstheorymax->SetLineStyle(2);

//Draw them
xstheorymax->Draw();
xstheory->Draw("same");
xstheorymin->Draw("same");
//sanity check - draw input points
gr1->Draw("same");
gr2->Draw("same");
gr3->Draw("same");


}