{

gStyle->SetPalette(1);

// I change the Npbins to 1100 to produce the 1d spectrum and to 110 to produce
// the 2d spectrum
Int_t Npbins = 1100;
Int_t Ncbins = 100;

TH2D* hMpcos;
hMpcos = new TH2D("hMpcos","",Npbins,0.0,55.,Ncbins,-1.0,1.0);
hMpcos->SetStats(kFALSE);
hMpcos->GetZaxis()->SetTitle("Events");
hMpcos->GetZaxis()->SetTitleSize(0.07);
hMpcos->GetZaxis()->SetTitleOffset(0.9);
hMpcos->GetZaxis()->SetNdivisions(505,kTRUE);
hMpcos->GetZaxis()->CenterTitle();
hMpcos->GetZaxis()->SetLabelSize(0.07);
hMpcos->GetXaxis()->SetTitle("Momentum (MeV/c)");
hMpcos->GetXaxis()->SetTitleSize(0.07);
hMpcos->GetXaxis()->SetTitleOffset(1.0);
hMpcos->GetXaxis()->SetNdivisions(505,kTRUE);
hMpcos->GetXaxis()->CenterTitle();
hMpcos->GetXaxis()->SetLabelSize(0.07);
hMpcos->GetYaxis()->SetTitle("cos #theta");
hMpcos->GetYaxis()->SetTitleSize(0.07);
hMpcos->GetYaxis()->SetTitleOffset(1.5);
hMpcos->GetYaxis()->SetNdivisions(505,kTRUE);
hMpcos->GetYaxis()->CenterTitle();
hMpcos->GetYaxis()->SetLabelSize(0.07);

TH1D* hMp;
hMp = new TH1D("hMp","",Npbins,0.0,55.);
hMp->SetStats(kFALSE);
hMp->GetYaxis()->SetTitle("Events");
hMp->GetYaxis()->SetTitleSize(0.1);
hMp->GetYaxis()->SetTitleOffset(0.15);
hMp->GetYaxis()->SetNdivisions(500,kTRUE);
hMp->GetYaxis()->CenterTitle();
hMp->GetXaxis()->SetTitle("Momentum (MeV/c)");
hMp->GetXaxis()->SetTitleSize(0.1);
hMp->GetXaxis()->SetTitleOffset(0.8);
hMp->GetXaxis()->SetLabelSize(0.1);
hMp->GetXaxis()->CenterTitle();
hMp->SetLineWidth(3);

TH1D* hMdiff;
hMdiff = new TH1D("hMdiff","",Npbins,0.0,55.);
hMdiff->SetStats(kFALSE);
hMdiff->GetYaxis()->SetTitle("Events");
hMdiff->GetYaxis()->SetTitleSize(0.1);
hMdiff->GetYaxis()->SetTitleOffset(0.15);
hMdiff->GetYaxis()->SetNdivisions(500,kTRUE);
hMdiff->GetYaxis()->CenterTitle();
hMdiff->GetXaxis()->SetTitle("Momentum (MeV/c)");
hMdiff->GetXaxis()->SetTitleSize(0.1);
hMdiff->GetXaxis()->SetTitleOffset(0.8);
hMdiff->GetXaxis()->SetLabelSize(0.1);
hMdiff->GetXaxis()->CenterTitle();
hMdiff->SetLineWidth(3);

TH2D* hframe;
hframe = new TH2D("hframe","",Npbins,0.0,55.,11,-0.2,0.6);
hframe->SetStats(kFALSE);
hframe->GetYaxis()->SetTitle("Events");
hframe->GetYaxis()->SetTitleSize(0.1);
hframe->GetYaxis()->SetTitleOffset(0.3);
hframe->GetYaxis()->SetNdivisions(100,kTRUE);
hframe->GetYaxis()->CenterTitle();
hframe->GetXaxis()->SetTitle("Momentum (MeV/c)");
hframe->GetXaxis()->SetTitleSize(0.1);
hframe->GetXaxis()->SetTitleOffset(0.8);
hframe->GetXaxis()->SetLabelSize(0.1);
hframe->GetXaxis()->CenterTitle();
TPaveText *plZero;
plZero = new TPaveText(-2.15,-0.1,-0.15,0.1);
plZero->SetFillStyle(0);
plZero->SetBorderSize(0);
plZero->SetTextSize(0.1);
TText *lZero = plZero->AddText("0");

TH1D* hMrho;
hMrho = new TH1D("hMrho","",Npbins,0.0,55.);
hMrho->SetStats(kFALSE);
hMrho->SetLineWidth(3);
hMrho->SetLineColor(2);
hMrho->SetLineStyle(2);

TH1D* hMeta;
hMeta = new TH1D("hMeta","",Npbins,0.0,55.);
hMeta->SetStats(kFALSE);
hMeta->SetLineWidth(4);
hMeta->SetLineColor(4);
hMeta->SetLineStyle(3);

TH2D* hfrdelta;
hfrdelta = new TH2D("hfrdelta","",Npbins,0.0,55.,10,-0.3,0.8);
hfrdelta->SetStats(kFALSE);
hfrdelta->GetYaxis()->SetTitle("Events");
hfrdelta->GetYaxis()->SetTitleSize(0.1);
hfrdelta->GetYaxis()->SetTitleOffset(0.3);
hfrdelta->GetYaxis()->SetNdivisions(100,kTRUE);
hfrdelta->GetYaxis()->CenterTitle();
hfrdelta->GetXaxis()->SetTitle("Momentum (MeV/c)");
hfrdelta->GetXaxis()->SetTitleSize(0.1);
hfrdelta->GetXaxis()->SetTitleOffset(0.8);
hfrdelta->GetXaxis()->SetLabelSize(0.1);
hfrdelta->GetXaxis()->CenterTitle();
TPaveText *plZero;
plZero = new TPaveText(-2.15,-0.1,-0.15,0.1);
plZero->SetFillStyle(0);
plZero->SetBorderSize(0);
plZero->SetTextSize(0.1);
TText *lZero = plZero->AddText("0");

TH1D* hMdelta;
hMdelta = new TH1D("hMdelta","",Npbins,0.0,55.);
hMdelta->SetStats(kFALSE);
hMdelta->SetLineWidth(3);
hMdelta->SetLineColor(2);
hMdelta->SetLineStyle(2);

Double_t pbinwidth = 55.0/Npbins;
Double_t cbinwidth = 2.0/Ncbins;

Double_t p;
Double_t x;
Double_t c;
Double_t Gp;
Double_t Gpc;
Double_t dxrho;
Double_t dxeta;
Double_t xnot = 9.97e-3;
Double_t Gdiff;
Double_t dxdelta;

for(Int_t pbin=1;pbin<=Npbins;pbin++){
  // Calculate bin content for momentum spectrum
  p = 0.0 + (pbin-0.5)*pbinwidth;
  x = p/55.0;
  Gp = 1.5*pow(x,2) - pow(x,3);
  hMp->SetBinContent(pbin,Gp);
  dxrho = 300.*0.0026*((8.0/3.0)*pow(x,3) - 2.0*pow(x,2));
  hMrho->SetBinContent(pbin,dxrho);
  dxeta = 300.*0.013*(3*xnot*(x - pow(x,2)));
  hMeta->SetBinContent(pbin,dxeta);
  Gdiff = pow(x,3) - 0.5*pow(x,2);
  hMdiff->SetBinContent(pbin,Gdiff);
  dxdelta = 300*0.0038*((8.0/3.0)*pow(x,3) - 2.0*pow(x,2));
  hMdelta->SetBinContent(pbin,dxdelta);
//   for(Int_t cbin=1;cbin<=Ncbins;cbin++){
//     // Calculate bin content for cos(theta) vs. momentum spectrum
//     c = -1.0 + (cbin-0.5)*cbinwidth;
//     Gpc = 1000*(Gp + c*(pow(x,3) - 0.5*pow(x,2)));
//     hMpcos->SetBinContent(pbin,cbin,Gpc);
//   }
}

// TCanvas *cp;
// cp = new TCanvas("cp","1d Michel distribution",600,300);
// cp->SetFillColor(0);
// cp->SetBorderSize(0.0);
// cp->SetBorderMode(0.0);
// cp->SetLeftMargin(0.05);
// cp->SetRightMargin(0.01);
// cp->SetTopMargin(0.01);
// cp->SetBottomMargin(0.2);
// hMp->Draw();

// TCanvas *cpcos;
// cpcos = new TCanvas("cpcos","2d Michel distribution",700,400);
// cpcos->SetFillColor(0);
// cpcos->SetBorderSize(0.0);
// cpcos->SetBorderMode(0.0);
// cpcos->SetLeftMargin(0.12);
// cpcos->SetRightMargin(0.17);
// cpcos->SetTopMargin(0.01);
// cpcos->SetBottomMargin(0.1);
// hMpcos->Draw("surf2z");

Double_t xZero[2] = {0.0001,54.9999};
Double_t yZero[2] = {0.0,0.0};

TGraph *gZero;
gZero = new TGraph(2,xZero,yZero);
gZero->SetLineColor(16);

// TLegend *leg;
// leg = new TLegend(0.1,0.6,0.4,1.0);
// leg->AddEntry(hMp,"Michel spectrum","L");
// leg->AddEntry(hMrho,"300X #Delta#rho effect","L");
// leg->AddEntry(hMeta,"300X #Delta#eta effect","L");
// leg->SetFillColor(0);
// leg->SetTextAlign(22);
// leg->SetFillStyle(0);
// leg->SetBorderSize(0);
// leg->SetTextSize(0.1);
// leg->SetTextFont(62);

// TCanvas *crhoeta;
// crhoeta = new TCanvas("crhoeta","Effect of rho and eta",600,300);
// crhoeta->SetFillColor(0);
// crhoeta->SetBorderSize(0.0);
// crhoeta->SetBorderMode(0.0);
// crhoeta->SetLeftMargin(0.08);
// crhoeta->SetRightMargin(0.01);
// crhoeta->SetTopMargin(0.01);
// crhoeta->SetBottomMargin(0.2);
// hframe->Draw();
// gZero->Draw("C");
// hMp->Draw("same");
// hMrho->Draw("same");
// hMeta->Draw("same");
// plZero->Draw();
// leg->Draw();

TLegend *legdelta;
legdelta = new TLegend(0.1,0.6,0.4,1.0);
legdelta->AddEntry(hMp,"Forward decays - backward decays","L");
legdelta->AddEntry(hMrho,"300X #Delta#delta effect","L");
legdelta->SetFillColor(0);
legdelta->SetTextAlign(22);
legdelta->SetFillStyle(0);
legdelta->SetBorderSize(0);
legdelta->SetTextSize(0.1);
legdelta->SetTextFont(62);

TCanvas *cdelta;
cdelta = new TCanvas("cdelta","Effect of delta",600,300);
cdelta->SetFillColor(0);
cdelta->SetBorderSize(0.0);
cdelta->SetBorderMode(0.0);
cdelta->SetLeftMargin(0.08);
cdelta->SetRightMargin(0.01);
cdelta->SetTopMargin(0.01);
cdelta->SetBottomMargin(0.2);
hfrdelta->Draw();
gZero->Draw("C");
hMdiff->Draw("same");
hMdelta->Draw("same");
plZero->Draw();
legdelta->Draw();
}