{
gROOT->Reset();

int nrap = 3;

float rap[3] = {0.0,-1.45,-1.95};
float raperr_plot[3] = {0.35,0.25,0.25};
float raperr[3] = {3*0};

//RUN2
///////float sigrap[3] = {61.8,45.2,30.5};  // wrong oficial numbers
float sigrap[3] = {61.8,49.6,28.5};
// statistical
//////float sigraperr[3] = {9.8,9.9,6.4};   // wrong official numbers  
float sigraperr[3] = {9.8,9.9,6.3};   
// sys1 error for electron point is 15.3% - final
// sys1 error for muon points is 13% & 5% folded in quadrature = 13.9% - final
//////float sigrapsys1[3] = {9.4,6.28,4.24};  // wrong official numbers  
float sigrapsys1[3] = {9.4,6.89,3.96};
// sys2 error for all is from abs luminosity (9.6%) and JPsi bias of minbias 
// trigger (take 3% from pi0 trigger measurement). It is final at 10.1%. 
//////float sigrapsys2[3] = {6.24,4.57,3.08};   // wrong official numbers  
float sigrapsys2[3] = {6.24,5.0,2.88};   

// run3pp
float run3rap[6] = {0.0,-1.45,-1.95,1.4,1.8,2.2};
float run3raperr_plot[6] = {0.35,0.25,0.25,0.2,0.2,0.2};  
float run3raperr[6] = {6*0};
float run3sigrap[6] = {61.8,35.90,22.72,42.66,28.91,16.49}; // 
//float run3sigrap[6] = {61.8,48.51,36.35,42.01,62.66,15.52}; // south : 0.21 north :0.52 from nicki
//float run3sigrap[5] = {61.8,57.87,39.34,84.33,40.21};
// statistical
float run3sigraperr[6] = {9.8,5.14,3.31,4.81,2.51,3.06};
//float run3sigraperr[6] = {9.8,6.95,5.30,4.73,5.43,2.88};
//float run3sigraperr[5] = {9.8,4.29,3.03,5.60,4.12};
// sys1 error for electron point is 15.3% - final
// sys1 error for muon points is 13% & 5% folded in quadrature = 13.9% - final
float run3sigrapsys1[6] = {0}; // unknown
// sys2 error for all is from abs luminosity (9.6%) and JPsi bias of minbias
// trigger (take 3% from pi0 trigger measurement). It is final at 10.1%.
float run3sigrapsys2[6] = {0}; // unknown




gStyle->SetOptStat(0);
gStyle->SetOptTitle(0);
gStyle->SetPadBottomMargin(0.17);
gStyle->SetPadLeftMargin(0.18);

float lblsize = 0.07;
float txtsize1 = 0.06;
float txtsize2 = 0.08;
float txtsize3 = 0.04;

// This is the one we fit, with zero error on rapidity

grap = new TGraphAsymmErrors(nrap,rap,sigrap,raperr,raperr,sigraperr,sigraperr);

// This is the one we plot

grap_plot = new TGraphAsymmErrors(nrap,rap,sigrap,raperr_plot,raperr_plot,sigraperr,sigraperr);

grap_plotrun3 = new TGraphAsymmErrors(6,run3rap,run3sigrap,run3raperr_plot,run3raperr_plot,run3sigraperr,run3sigraperr);

// For plotting systematic errors
float sysplot_up[3];
float sysplot_dn[3];
for(int i=0;i<3;i++)
{
  sysplot_up[i]=sigrap[i]+sigrapsys1[i];
  sysplot_dn[i]=sigrap[i]-sigrapsys1[i];

  cout << i << " " << sysplot_up[i] << " " << sysplot_dn[i] << endl;
}

cc = new TCanvas("cc", "", 800, 500);
gPad->SetLogy(0);

TH1F *h = new TH1F("dummy","dummy",100,-3.0,3.0);
//TH1F *h = new TH1F("dummy","dummy",100,-3.7,3.7);
h->SetMaximum(75);
//h->SetMaximum(120);
h->SetMinimum(0);
h->SetNdivisions(7,"X");
h->SetNdivisions(5,"Y");
h->SetLabelSize(lblsize,"X");
h->SetLabelSize(lblsize,"Y");
h->Draw();

grap_plot->SetLineColor(4);
grap_plot->SetLineWidth(2);
grap_plot->SetMarkerColor(4);
grap_plot->SetMarkerStyle(24);

grap_plot->Draw("P");

grap_plotrun3->SetLineColor(2);
grap_plotrun3->SetLineWidth(2);
grap_plotrun3->SetMarkerColor(2);
grap_plotrun3->SetMarkerStyle(20);
grap_plotrun3->SetMarkerSize(1.1);

grap_plotrun3->Draw("P");


tex = new TLatex(-3.9,25.0,"B_{ll} #frac{d#sigma}{dy} (nb)");
tex->SetTextFont(132);
tex->SetTextSize(txtsize2);
tex->SetTextAngle(90.0);
tex->SetLineWidth(2);
tex->Draw();

tex = new TLatex(1.0,-15,"rapidity");
tex->SetTextFont(132);
tex->SetTextSize(txtsize2);
tex->SetLineWidth(2);
tex->Draw();


//////////////////////////////////////////////////////////////////////
// Some different pythia rapidity plots
/////////////////////////////////////////////////////////////////////

float xmin = -3.75;  // Pythia goes to zero by 3.5
float xmax = 3.75;
int nbin = 120;

if(1)
{
  // From Hiroki
  
  // Use this Pythia distribution normalized to the data to get the 
  // total cross section 
  
  TFile *fin = new TFile("/phenix/data25/satohiro/pythia/Dimuons_nochi_noselect_GRV94HO_10000.root");
  
  //TFile *fin = new TFile("/phenix/data25/satohiro/pythia/Dimuons_nochi_noselect_noptcut_kt1.5_10000.root");
  //TFile *fin = new TFile("/phenix/data25/satohiro/pythia/Dimuons_nochi_noselect_noptcut_kt04_10000.root");  
  // Here are some alternative Pythia distributions
  //TFile *fin = new TFile("/phenix/data25/satohiro/pythia/Dimuons_nochi_noselect_CTEQ4M_10000.root");
  //TFile *fin = new TFile("/phenix/data25/satohiro/pythia/Dimuons_nochi_noselect_noptcut_kt082.root");
  //TFile *fin = new TFile("/phenix/data25/satohiro/pythia/Dimuons_nochi_noselect_noptcut_kt1.root");
  
  
  TH1F* h101 = new TH1F("h101","dimuon mass",nbin,xmin,xmax);
  TH1F* h101p = new TH1F("h101p","dimuon mass",nbin,xmin,xmax);
  TH1F* h101n = new TH1F("h101n","dimuon mass",nbin,xmin,xmax);
  
  float norm1 = 3585*0.06; //GRV94LO
  float norm2 = 3547*0.06; // CTEQ4M
  
  TCut cut2 = "mass>3.09 && mass<3.10";
  
  cout << "Project" << endl;
  
  Dimuons->Project("h101p","abs(y)","mass>3.09 && mass<3.10 && (pid1*pid2<0)");
  Dimuons->Project("h101n","-abs(y)","mass>3.09 && mass<3.10 && (pid1*pid2<0)");
  
  h101->Add(h101p,h101n,0.5,0.5);
  float normfactor = (norm1/((xmax-xmin)*nbin));
  normfactor = 1.09*normfactor;
  h101->Scale(normfactor);
  h101->SetLineWidth(3);
  //h101->Draw();
  //h101->Draw("same");
  
  // Fit the pythia result with a smooth curve
  
  TF1 *f1 = new TF1("f1","[0]+[1]*x**2+[2]*x**4+[3]*x**6",-3.5,3.5);
  h101->Fit("f1","RN");
  f1->SetLineColor(2);
  //f1->Draw("same");
  
  // Normalize the pythia curve to the data

  double p0 = f1->GetParameter(0);
  double p1 = f1->GetParameter(1);
  double p2 = f1->GetParameter(2);
  double p3 = f1->GetParameter(3);
  //double p4 = f1->GetParameter(4);
  TF1 *f2 = new TF1("f2","[0]*([1]+[2]*x**2+[3]*x**4+[4]*x**6)",-3.5,3.5);
  f2->SetParameter(0,1.0);
  f2->SetParameter(1,p0);
  f2->SetParameter(2,p1);
  f2->SetParameter(3,p2);
  f2->SetParameter(4,p3);
  // fix the parameter values from fit to Pythia
  f2->SetParLimits(1,p0,p0);
  f2->SetParLimits(2,p1,p1);
  f2->SetParLimits(3,p2,p2);
  f2->SetParLimits(4,p3,p3);
  grap->Fit("f2","RN");
  f2->SetLineStyle(4);  // dot-dashed
  f2->SetLineWidth(3);  // thick
  f2->SetLineColor(2);  // red
  f2->Draw("same");

  // Get the total cross section from the normalized Pythia curve

  double norm = f2->GetParameter(0);
  double norm_err = f2->GetParError(0);

  double sigtot = f2->Integral(-3.5,3.5);
  double sigtot_err = sigtot * norm_err/norm;

  double bll = 0.0599;

  cout << "Total cross section (microbarns) = " << sigtot/(bll*1000) 
       << " stat error = " << sigtot_err/(bll*1000) << endl;

  tex = new TLatex(-0.8,9.5,"Scaled PYTHIA (GRV94HO)");
  //tex = new TLatex(-0.8,9.5,"Scaled PYTHIA (noptcut_kt1.5)");
  //tex = new TLatex(-0.8,9.5,"Scaled PYTHIA (noptcut_kt04)");
  tex->SetTextFont(131);
  tex->SetTextSize(0.05);
  //tex->Draw();

  a1 = new TArrow();
  a1.SetArrowSize(0.01);
  a1.SetLineWidth(2);
  //a1.DrawArrow(2.1,10.5,2.75,10.5,0.02);

}


// Plot Systematic errors

float ln = 0.10;
float ln2 = 2.0;
float lw = 1;
float lc = 1;

for (int i = 0;i<3;i++)
{
  TLine *line1 = new TLine(rap[i]-ln,sysplot_up[i],rap[i]+ln,sysplot_up[i]);
  line1->SetLineColor(lc);
  line1->SetLineWidth(lw);
  line1->Draw();

  TLine *line3 = new TLine(rap[i]-ln,sysplot_up[i],rap[i]-ln,sysplot_up[i]-ln2);
  line3->SetLineColor(lc);
  line3->SetLineWidth(lw);
  line3->Draw();

  TLine *line5 = new TLine(rap[i]+ln,sysplot_up[i],rap[i]+ln,sysplot_up[i]-ln2);
  line5->SetLineColor(lc);
  line5->SetLineWidth(lw);
  line5->Draw();


  TLine *line2 = new TLine(rap[i]-ln,sysplot_dn[i],rap[i]+ln,sysplot_dn[i]);
  line2->SetLineColor(lc);
  line2->SetLineWidth(lw);
  line2->Draw();

  TLine *line4 = new TLine(rap[i]-ln,sysplot_dn[i],rap[i]-ln,sysplot_dn[i]+ln2);
  line4->SetLineColor(lc);
  line4->SetLineWidth(lw);
  line4->Draw();

  TLine *line6 = new TLine(rap[i]+ln,sysplot_dn[i],rap[i]+ln,sysplot_dn[i]+ln2);
  line6->SetLineColor(lc);
  line6->SetLineWidth(lw);
  line6->Draw();

}

cc->Update();

////////////////////////////////////////////////////////////////////////////
// Estimate the sys1 error on the total cross section
// The sys2 error is just the absolute normalization, so it is easy
////////////////////////////////////////////////////////////////////////////

// Make two new graphs with the points at the limits of the sys1 errors

// First up

gsys_up = new TGraphAsymmErrors(nrap,rap,sysplot_up,raperr,raperr,sigraperr,sigraperr);

// Normalize the pythia curve to the data

p0 = f1->GetParameter(0);
p1 = f1->GetParameter(1);
p2 = f1->GetParameter(2);
p3 = f1->GetParameter(3);
TF1 *f11 = new TF1("f11","[0]*([1]+[2]*x**2+[3]*x**4+[4]*x**6)",-3.5,3.5);
f11->SetParameter(0,1.0);
f11->SetParameter(1,p0);
f11->SetParameter(2,p1);
f11->SetParameter(3,p2);
f11->SetParameter(4,p3);
// fix the parameter values from fit to Pythia
f11->SetParLimits(1,p0,p0);
f11->SetParLimits(2,p1,p1);
f11->SetParLimits(3,p2,p2);
f11->SetParLimits(4,p3,p3);
gsys_up->Fit("f11","RN");
f11->SetLineColor(2);
//f11->Draw("same");

// Get the total cross section from the normalized Pythia curve

double norm_up = f11->GetParameter(0);
double norm_up_err = f11->GetParError(0);

double sigtot_up = f11->Integral(-3.5,3.5);
double sigtot_up_err = sigtot * norm_err/norm;

cout << "Total sys1_up cross section (microbarns) = " << sigtot_up/(bll*1000) 
  << " stat error = " << sigtot_up_err/(bll*1000) << endl;

// Now down

gsys_up = new TGraphAsymmErrors(nrap,rap,sysplot_dn,raperr,raperr,sigraperr,sigraperr);

// Normalize the pythia curve to the data

p0 = f1->GetParameter(0);
p1 = f1->GetParameter(1);
p2 = f1->GetParameter(2);
TF1 *f12 = new TF1("f12","[0]*([1]+[2]*x**2+[3]*x**4+[4]*x**6)",-3.5,3.5);
f12->SetParameter(0,1.0);
f12->SetParameter(1,p0);
f12->SetParameter(2,p1);
f12->SetParameter(3,p2);
f12->SetParameter(4,p3);
// fix the parameter values from fit to Pythia
f12->SetParLimits(1,p0,p0);
f12->SetParLimits(2,p1,p1);
f12->SetParLimits(3,p2,p2);
f12->SetParLimits(4,p3,p3);
gsys_up->Fit("f12","RN");
f12->SetLineColor(2);
//f12->Draw("same");

// Get the total cross section from the normalized Pythia curve
// 3.5 is the rapidity at which Pythia goes to zero, and is also the 
// limit of the validity of the fit to Pythia

double norm_dn = f12->GetParameter(0);
double norm_dn_err = f12->GetParError(0);

double sigtot_dn = f12->Integral(-3.5,3.5);
double sigtot_dn_err = sigtot * norm_err/norm;

cout << "Total sys1_dn cross section (microbarns) = " << sigtot_dn/(bll*1000) 
  << " stat error = " << sigtot_dn_err/(bll*1000) << endl;

double sigtot_sys1 = (sigtot_up-sigtot_dn)/2;

cout << " sigtot_sys1 = " << sigtot_sys1/(bll*1000) << endl;

/////////////////////////////////////////////////////////////////////////
// Output a summary of the cross section results
/////////////////////////////////////////////////////////////////////////

cout << endl << "   ******** Summary of cross section results: ******" << endl;

cout << "Rapidity         B*dsigma/dy (nb)       " << endl;

for (int i=0;i<3;i++)
{
  cout << rap[i] << "          " <<  sigrap[i] << " +/- " << sigraperr[i] 
       << " (stat) +/- " << sigrapsys1[i] << " (sys) +/- " << sigrapsys2[i] 
       << " abs " << endl;
}

// The absolute luminosity and JPsi trigger bias error 
// final value. It is the same for all points

float sigabs = sigtot * sigrapsys2[0]/sigrap[0];

cout << "Total B * sigma (nb):" << endl;

cout << sigtot << " +/- " << sigtot_err << " (stat) +/- " << sigtot_sys1 
<< " (sys) +/- " << sigabs << " (abs) " << endl;

cout << "Total sigma (microbarns):" << endl;

cout << sigtot/(bll*1000) << " +/- " << sigtot_err/(bll*1000) 
<< " (stat) +/- " << sigtot_sys1/(bll*1000) << " (sys) +/- " 
<< sigabs/(bll*1000) << " (abs) " << endl;


/////////////////////////////////////////////////////////////////////////
// Output a summary of the dN/dy results
/////////////////////////////////////////////////////////////////////////

// All we need is (sigma_inelastic * eff_minbias) = 21.8 mb = 21.8 * 10**6 nb
// to convert sigma to N

// Wait, this just gives yield per minbias trigger - we need eff_minbias
// to get yield/inelastic collision

// Get eff_minbias from luminosity note !!

float sigbbc = 21.8 * 1000000;

cout << endl << "   ******** Summary of dN/dy results: ******" << endl;

cout << "Rapidity         B*dN/dy        " << endl;

for (int i=0;i<3;i++)
{
  cout << rap[i] << "          " <<  sigrap[i]/sigbbc << " +/- " 
       << sigraperr[i]/sigbbc 
       << " (stat) +/- " << sigrapsys1[i]/sigbbc << " (sys) +/- " 
       << sigrapsys2[i]/sigbbc 
       << " abs " << endl;
}

float ntot = sigtot/sigbbc;
float ntot_err = sigtot_err/sigbbc;

// Systematic error on the total cross section
float nsys = ntot * sigtot_sys1/sigtot;

// The absolute luminosity and JPsi trigger bias error - final value
// same for all points
float nabs = ntot * sigrapsys2[0]/sigrap[0];

cout << "Total B * JPsi yield/event:" << endl;

cout << ntot << " +/- " << ntot_err << " (stat) +/- " << nsys 
<< " (sys) +/- " << nabs << " (abs) " << endl;

cout << "Total JPsi yield/event:" << endl;

cout << ntot/bll << " +/- " << ntot_err/bll << " (stat) +/- " 
<< nsys/bll << " (sys) +/- " << nabs/bll << " (abs) " << endl;

if(0)
{
  // Ramona Vogt rapidity distributions from CEM

  // start now with (MRST HO, m_c=1.2, mu = 2m)

  float rvrap[36];
  float rvdsdy[36];
  float rvdsdy_err[36];
  float rvrap_err[36]={36*0};

  ifstream rvin;
  rvin.open("RVogt_cem_rap_dist.txt");

  for(int i=0;i<36;i++)
    {
      rvin >> rvrap[i] >> rvdsdy[i] >> rvdsdy_err[i];
      cout << " i " << " rap " << rvrap[i] << " ds/dy " << rvdsdy[i] << endl;
    }

  graprv1 = new TGraphErrors(36,rvrap,rvdsdy,rvrap_err,rvdsdy_err);
  //graprv1->Draw("same");

  // Fit with a polynomial

  TF1 *f31 = new TF1("f31","[0]+[1]*x**2+[2]*x**4+[3]*x**6",xmin,xmax);

  graprv1->Fit("f31","RN");
  //f31->Draw("same");

  // Normalize the curve to the data

  p0 = f31->GetParameter(0);
  p1 = f31->GetParameter(1);
  p2 = f31->GetParameter(2);
  p3 = f31->GetParameter(3);
  TF1 *f32 = new TF1("f32","[0]*([1]+[2]*x**2+[3]*x**4+[4]*x**6)",xmin,xmax);
  f32->SetParameter(0,1.0);
  f32->SetParameter(1,p0);
  f32->SetParameter(2,p1);
  f32->SetParameter(3,p2);
  f32->SetParameter(4,p3);
  // fix the parameter values from fit to Pythia. leaving only the scale to vary
  f32->SetParLimits(1,p0,p0);
  f32->SetParLimits(2,p1,p1);
  f32->SetParLimits(3,p2,p2);
  f32->SetParLimits(4,p3,p3);
  grap->Fit("f32","RN");
  f32->SetLineColor(1);   // black
  f32->SetLineStyle(4);
//  f32->Draw("same");  //dj

  texrv1 = new TLatex(-1.25,25,"CEM (MRST HO)");
  texrv1->SetTextFont(132);
  texrv1->SetTextSize(0.05);
  texrv1->Draw();

  texrv1a = new TLatex(-1.25,20,"Scaled");
  texrv1a->SetTextFont(132);
  texrv1a->SetTextSize(0.05);
  texrv1a->Draw();

  arv1 = new TArrow();
  arv1.SetArrowSize(0.01);
  arv1.SetLineWidth(2);
  arv1.DrawArrow(0.0,30,0.6,52,0.02);

  // Get the total cross section from the normalized curve

  norm = f32->GetParameter(0);
  norm_err = f32->GetParError(0);

  sigtot = f32->Integral(-3.5,3.5);
  cout << " RV1 cross section is " << sigtot/(bll*1000) << endl;

}

if(0)
{
  // Now do (MRST HO, m=1.4, mu = m)
  
  for(int i=0;i<36;i++)
    {
      rvin >> rvrap[i] >> rvdsdy[i] >> rvdsdy_err[i];
      cout << " i " << " rap " << rvrap[i] << " ds/dy " << rvdsdy[i] << endl;
    }

  graprv2 = new TGraphErrors(36,rvrap,rvdsdy,rvrap_err,rvdsdy_err);
  //graprv2->Draw("same");

  // Fit with a polynomial

  TF1 *f41 = new TF1("f41","[0]+[1]*x**2+[2]*x**4+[3]*x**6",xmin,xmax);

  graprv2->Fit("f41","RN");
  //f41->Draw("same");

  // Normalize the curve to the data

  p0 = f41->GetParameter(0);
  p1 = f41->GetParameter(1);
  p2 = f41->GetParameter(2);
  p3 = f41->GetParameter(3);
  TF1 *f42 = new TF1("f42","[0]*([1]+[2]*x**2+[3]*x**4+[4]*x**6)",xmin,xmax);
  f42->SetParameter(0,1.0);
  f42->SetParameter(1,p0);
  f42->SetParameter(2,p1);
  f42->SetParameter(3,p2);
  f42->SetParameter(4,p3);
  // fix the parameter values from fit to shape. leaving only the scale to vary
  f42->SetParLimits(1,p0,p0);
  f42->SetParLimits(2,p1,p1);
  f42->SetParLimits(3,p2,p2);
  f42->SetParLimits(4,p3,p3);
  grap->Fit("f42","RN");
  f42->SetLineColor(3);  // green
  f42->Draw("same");

  texg = new TLatex(-1.25,30,"MRST HO m=1.4 mu=m");
  texg->SetTextFont(132);
  texg->SetTextSize(txtsize3);
  texg->SetLineWidth(2);
  texg->Draw();

  // Get the total cross section from the normalized curve

  norm = f42->GetParameter(0);
  norm_err = f42->GetParError(0);

  sigtot = f42->Integral(-3.5,3.5);
  cout << " RV2 cross section is " << sigtot/(bll*1000) << endl;

  // Now do (CTEQ 5M, m=1.2, mu = 2m)

  for(int i=0;i<36;i++)
    {
      rvin >> rvrap[i] >> rvdsdy[i] >> rvdsdy_err[i];
      cout << " i " << " rap " << rvrap[i] << " ds/dy " << rvdsdy[i] << endl;
    }

  graprv3 = new TGraphErrors(36,rvrap,rvdsdy,rvrap_err,rvdsdy_err);
  //graprv3->Draw("same");

  // Fit with a polynomial

  TF1 *f51 = new TF1("f51","[0]+[1]*x**2+[2]*x**4+[3]*x**6",xmin,xmax);

  graprv3->Fit("f51","RN");
  //f51->Draw("same");

  // Normalize the curve to the data

  p0 = f51->GetParameter(0);
  p1 = f51->GetParameter(1);
  p2 = f51->GetParameter(2);
  p3 = f51->GetParameter(3);
  TF1 *f52 = new TF1("f52","[0]*([1]+[2]*x**2+[3]*x**4+[4]*x**6)",xmin,xmax);
  f52->SetParameter(0,1.0);
  f52->SetParameter(1,p0);
  f52->SetParameter(2,p1);
  f52->SetParameter(3,p2);
  f52->SetParameter(4,p3);
  // fix the parameter values from fit to shape. leaving only the scale to vary
  f52->SetParLimits(1,p0,p0);
  f52->SetParLimits(2,p1,p1);
  f52->SetParLimits(3,p2,p2);
  f52->SetParLimits(4,p3,p3);
  grap->Fit("f52","RN");
  f52->SetLineColor(7);  // cyan
  f52->Draw("same");

  texg = new TLatex(-1.25,25,"CTEQ 5M, m=1.2, mu = 2m");
  texg->SetTextFont(132);
  texg->SetTextSize(txtsize3);
  texg->SetLineWidth(2);
  texg->Draw();

  // Get the total cross section from the normalized curve

  norm = f52->GetParameter(0);
  norm_err = f52->GetParError(0);

  sigtot = f52->Integral(-3.5,3.5);
  cout << " RV3 cross section is " << sigtot/(bll*1000) << endl;

// Now do (GRV98HO, m=1.3, mu = m)

  for(int i=0;i<36;i++)
    {
      rvin >> rvrap[i] >> rvdsdy[i] >> rvdsdy_err[i];
      cout << " i " << " rap " << rvrap[i] << " ds/dy " << rvdsdy[i] << endl;
    }

  graprv4 = new TGraphErrors(36,rvrap,rvdsdy,rvrap_err,rvdsdy_err);
  //graprv4->Draw("same");

  // Fit with a polynomial

  TF1 *f61 = new TF1("f61","[0]+[1]*x**2+[2]*x**4+[3]*x**6",xmin,xmax);

  graprv4->Fit("f61","RN");
  //f61->Draw("same");

  // Normalize the curve to the data

  p0 = f61->GetParameter(0);
  p1 = f61->GetParameter(1);
  p2 = f61->GetParameter(2);
  p3 = f61->GetParameter(3);
  TF1 *f62 = new TF1("f62","[0]*([1]+[2]*x**2+[3]*x**4+[4]*x**6)",xmin,xmax);
  f62->SetParameter(0,1.0);
  f62->SetParameter(1,p0);
  f62->SetParameter(2,p1);
  f62->SetParameter(3,p2);
  f62->SetParameter(4,p3);
  // fix the parameter values from fit to shape. leaving only the scale to vary
  f62->SetParLimits(1,p0,p0);
  f62->SetParLimits(2,p1,p1);
  f62->SetParLimits(3,p2,p2);
  f62->SetParLimits(4,p3,p3);
  grap->Fit("f62","RN");
  f62->SetLineColor(8);  // dark green
  f62->Draw("same");

  texg = new TLatex(-1.75,15,"GRV98HO, m=1.3, mu = m");
  texg->SetTextFont(132);
  texg->SetTextSize(txtsize3);
  texg->SetLineWidth(2);
  texg->Draw();

  // Get the total cross section from the normalized curve

  norm = f62->GetParameter(0);
  norm_err = f62->GetParError(0);

  sigtot = f62->Integral(-3.75,3.75);

  cout << " RV4 cross section is " << sigtot/(bll*1000) << endl;

}

if(1)
{
  // Read in Hiroki's latest COM with GRV98HO rapidity distributions 

  ifstream rapdist;
  rapdist.open("com_rapidity_GRV98HO_scale_Mjpsi_Mc_1.48.dat");

  float rapin, dsdy;
  float rapx[100],rapy[100];
  int i = 0;
  int nrap = 100;

  for(int irap=0;irap<nrap/2;irap++)
    {
      rapdist >> rapin >> dsdy;

      cout << " irap " << irap << " rapin " << rapin 
	   << " dsdy " << dsdy << endl;

      rapx[nrap/2+irap] = rapin;
      rapy[nrap/2+irap] = dsdy;
      rapx[nrap/2-irap-1] = -rapin;
      rapy[nrap/2-irap-1] = dsdy;
    }

  grap1 = new TGraph(nrap,rapx,rapy);

  // Fit with a polynomial

  TF1 *f71 = new TF1("f71","[0]+[1]*x**2+[2]*x**4+[3]*x**6",xmin,xmax);

  grap1->Fit("f71","RN");

  // Normalize the curve to the data

  p0 = f71->GetParameter(0);
  p1 = f71->GetParameter(1);
  p2 = f71->GetParameter(2);
  p3 = f71->GetParameter(3);
  TF1 *f72c = new TF1("f72c","[0]*([1]+[2]*x**2+[3]*x**4+[4]*x**6)",xmin,xmax);
  f72c->SetParameter(0,1.0);
  f72c->SetParameter(1,p0);
  f72c->SetParameter(2,p1);
  f72c->SetParameter(3,p2);
  f72c->SetParameter(4,p3);
  // fix the parameter values from fit to shape. leaving only the scale to vary
  f72c->SetParLimits(1,p0,p0);
  f72c->SetParLimits(2,p1,p1);
  f72c->SetParLimits(3,p2,p2);
  f72c->SetParLimits(4,p3,p3);
  grap->Fit("f72c","RN");
  f72c->SetLineColor(1);  // black
  f72c->SetLineStyle(1);  // solid
  f72c->SetLineWidth(3);  // ordinary

  f72c->Draw("same");

  // Get the cross section from the integral of the normalized function

  float sigest = f72c->Integral(xmin,xmax)/(bll*1000);
  cout << " Cross section estimate is: " << sigest << endl;

  tex = new TLatex(-1.4,21.0,"Scaled COM (GRV98HO)");
  //tex = new TLatex(-1.4,21.0,"COM (GRV98HO)");
  tex->SetTextFont(131);
  tex->SetTextSize(0.05);
  //tex->Draw();

  a1 = new TArrow();
  a1.SetArrowSize(0.01);
  a1.SetLineWidth(2);
  //a1.DrawArrow(1.30,22,2.55,25.0,0.02);

  //////////////////////////////////////////////////////////////////////////
  // The following was a one time thing to see what the absolutely 
  // normalized COM prediction would look like - it is not used in the figure
  //////////////////////////////////////////////////////////////////////////

  // The cross section for this set in fig. 4b is 4.632 microbarns.
  // renormalize the fit so that the integral matches

  float rawint = f71->Integral(xmin,xmax)/(bll*1000);

  float renormfact = 4.632/rawint;

  cout << "GRVHO rawint = " << rawint << " from " << xmin << " to " << xmax 
       << " renormfact = " << renormfact << endl;

  TF1 *f73c= new TF1("f73c","1018.49*([0]+[1]*x**2+[2]*x**4+[3]*x**6)",xmin,xmax);

  f73c->SetParameter(0,p0);
  f73c->SetParameter(1,p1);
  f73c->SetParameter(2,p2);
  f73c->SetParameter(3,p3);

  f73c->SetLineColor(1);  // black
  f73c->SetLineStyle(1);  // solid
  f73c->SetLineWidth(3);  // ordinary

  //f73c->Draw("same");
}

if(1)
{
  // Read in Hiroki's latest COM with MRST2001 rapidity distributions 

  ifstream rapdist;
  //rapdist.open("com_rapidity_MRST2001Set1_scale_Mjpsi_Mc_1.48.dat");
  rapdist.open("com_rapidity_MRST2001Set1_scale_2.3_Mc_1.55.dat");

  float rapin, dsdy;
  float rapx[100],rapy[100];
  int i = 0;
  int nrap = 100;

  for(int irap=0;irap<nrap/2;irap++)
    {
      rapdist >> rapin >> dsdy;

      rapx[nrap/2+irap] = rapin;
      rapy[nrap/2+irap] = dsdy;
      rapx[nrap/2-irap-1] = -rapin;
      rapy[nrap/2-irap-1] = dsdy;
    }

  grap1 = new TGraph(nrap,rapx,rapy);

  // Fit with a polynomial

  TF1 *f71 = new TF1("f71","[0]+[1]*x**2+[2]*x**4+[3]*x**6",xmin,xmax);

  grap1->Fit("f71","RN");

  // Normalize the curve to the data

  p0 = f71->GetParameter(0);
  p1 = f71->GetParameter(1);
  p2 = f71->GetParameter(2);
  p3 = f71->GetParameter(3);
  TF1 *f72a = new TF1("f72a","[0]*([1]+[2]*x**2+[3]*x**4+[4]*x**6)",xmin,xmax);
  f72a->SetParameter(0,1.0);
  f72a->SetParameter(1,p0);
  f72a->SetParameter(2,p1);
  f72a->SetParameter(3,p2);
  f72a->SetParameter(4,p3);
  // fix the parameter values from fit to shape. leaving only the scale to vary
  f72a->SetParLimits(1,p0,p0);
  f72a->SetParLimits(2,p1,p1);
  f72a->SetParLimits(3,p2,p2);
  f72a->SetParLimits(4,p3,p3);
  grap->Fit("f72a","RN");
  f72a->SetLineColor(4);  // blue
  f72a->SetLineStyle(2);  // dashed
  f72a->SetLineWidth(3);  // ordinary
  f72a->Draw("same"); // 

  // Get the cross section from the integral of the normalized function

  float sigest = f72a->Integral(xmin,xmax)/(bll*1000);
  cout << " Cross section estimate is: " << sigest << endl;

  tex = new TLatex(-1.2,32.0,"Scaled COM (MRST2001HO)");
  //tex = new TLatex(-1.2,32.0,"COM (MRST2001)");
  tex->SetTextFont(131);
  tex->SetTextSize(0.05);
  //tex->Draw();

  a1 = new TArrow();
  a1.SetArrowSize(0.01);
  a1.SetLineWidth(2);
  //a1.DrawArrow(0.0,36.0,0.6,53.5,0.02);

  //////////////////////////////////////////////////////////////////////////
  // The following was a one time thing to see what the absolutely 
  // normalized COM prediction would look like - it is not used in the figure
  //////////////////////////////////////////////////////////////////////////

  // The cross section for this set in fig. 4b is 3.579 microbarns.
  // renormalize the fit so that the integral matches

  ////float rawint = f71->Integral(xmin,xmax)/(bll*1000);

  ////float renormfact = 3.579/rawint;

  ////cout << "MRST2001 rawint = " << rawint << " from " << xmin << " to " << xmax 
  ////     << " renormfact = " << renormfact << endl;

  TF1 *f73a = new TF1("f73a","1019.92*([0]+[1]*x**2+[2]*x**4+[3]*x**6)",xmin,xmax);

  f73a->SetParameter(0,p0);
  f73a->SetParameter(1,p1);
  f73a->SetParameter(2,p2);
  f73a->SetParameter(3,p3);

  f73a->SetLineColor(4);  // blue
  f73a->SetLineStyle(2);  // dashed
  f73a->SetLineWidth(3);

  //f73a->Draw("same");
}


leg = new TLegend( 0.35, 0.25, 0.75, 0.45);
leg->SetBorderSize(0);
leg->SetFillColor(0);

leg->AddEntry( f2, "Pythia (GRV94NLO)", "l");
//leg->AddEntry( f72b, "COM (GRV94NLO)", "l");
leg->AddEntry( f72c, "COM (GRV98NLO)", "l");
leg->AddEntry( f72a, "COM (MRST2001NLO)", "l");

leg->Draw();


}