Double_t fWave(Double_t *x,Double_t *parPtr) {
  Double_t fitVal,xVal;

  xVal = *x;
  fitVal = 0.0;

  double a0     = parPtr[0];

  double cosVal = xVal;
  double sinVal = sqrt(1-cosVal*cosVal);
  double cosHalf = sqrt((1+cosVal)/2);
  double sinHalf = sqrt((1-cosVal)/2);

  double a10Mag  = parPtr[1];
  double phase10 = parPtr[2];
  double a10re   = a10Mag*cos(phase10);
  double a10im   = a10Mag*sin(phase10);
  double f10 = 0.5*(1+cosVal)*cosHalf;

  double a11Mag  = parPtr[3];
  double phase11 = parPtr[4];
  double a11re   = a11Mag*cos(phase11);
  double a11im   = a11Mag*sin(phase11);
  double f11 = -sqrt(3)*0.5*(1+cosVal)*sinHalf;

  double a12Mag  = parPtr[5];
  double phase12 = parPtr[6];
  double a12re   = a12Mag*cos(phase12);
  double a12im   = a12Mag*sin(phase12);
  double f12 = sqrt(3)*0.5*(1-cosVal)*sinHalf;

  double a13Mag  = parPtr[7];
  double phase13 = parPtr[8];
  double a13re   = a13Mag*cos(phase13);
  double a13im   = a13Mag*sin(phase13);
  double f13 = -0.5*(1-cosVal)*sinHalf;

  double a14Mag  = parPtr[9];
  double phase14 = parPtr[10];
  double a14re   = a14Mag*cos(phase14);
  double a14im   = a14Mag*sin(phase14);
  double f14 = 0.5*(3*cosVal-1)*cosHalf;

  double a15Mag  = parPtr[11];
  double phase15 = parPtr[12];
  double a15re   = a15Mag*cos(phase15);
  double a15im   = a15Mag*sin(phase15);
  double f15 = -0.5*(3*cosVal+1)*sinHalf;
  

  double fitValRe = a0 + a10re*f10 + a11re*f11 +
    a12re*f12 + a13re*f13 + a14re*f14 + a15re*f15;
  
  double fitValIm = a10im*f10 + a11im*f11 +
    a12im*f12 + a13im*f13 + a14im*f14 + a15im*f15;

  fitVal = pow(fitValRe,2) + pow(fitValIm,2);

}
Double_t fTest(Double_t *x,Double_t *parPtr) {
  Double_t fitVal,xVal;

  xVal = *x;
  fitVal = 0.0;

  double gMag   = parPtr[0];
  double center = parPtr[1];
  double sigma  = parPtr[2];

  double gFun = TMath::Gaus(xVal,center,sigma)*gMag;

  double bkg = parPtr[3] + parPtr[4]*xVal + parPtr[5]*pow(xVal,2);

  double sCenter1 = parPtr[6];
  double sWidth1 = parPtr[7];
  double sFun1 = 1.0/(exp((-xVal + sCenter1)*sWidth1) + 1);

  //double sCenter2 = parPtr[8];
  //double sWidth2 = parPtr[9];
  //double sFun2 = 1.0/(exp((xVal - sCenter2)*sWidth2) + 1);

  //double gMag2   = parPtr[10];
  //double center2 = parPtr[11];
  //double width2  = parPtr[12];
  //double gFun2 = TMath::Voigt(xVal-center2,sigma,width2)*gMag2;
  //double gMag3   = parPtr[13];
  //double center3 = parPtr[14];
  //double width3  = parPtr[15];
  //double gFun3 = TMath::Gaus(xVal,center3,width3)*gMag3;
  //fitVal = fabs(gFun) + sFun1*bkg*sFun2 + fabs(gFun2) + fabs(gFun3);

  fitVal = fabs(gFun) + sFun1*bkg;

  return fitVal;

}



void cosThetaH(int iShow){
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);


  char  hId[120];
  TFile *hFileMC = new TFile("./allBrandonMC.root");
  TFile *hFile = new TFile("./allSp18.root");


  TH2D *hMassXY; hFile->GetObject("hCosThetaGJyVsCasStar_MC_Kin_4K",hMassXY);
  TH2D *hMassXYMC; hFileMC->GetObject("hCosThetaGJyVsCasStar_MC_Kin_4K",hMassXYMC);

  TH1D *hAcc = (TH1D*) hMassXYMC->ProjectionY("hAcc",1,40,"");
  
  TH1D *hMassY = (TH1D*) hMassXY->ProjectionY("hY",13,15,"");

  TH1D *hMassX1 = (TH1D*) hMassXY->ProjectionX("hX1",1,2,"");
  TH1D *hMassX2 = (TH1D*) hMassXY->ProjectionX("hX2",3,4,"");
  TH1D *hMassX3 = (TH1D*) hMassXY->ProjectionX("hX3",5,6,"");
  TH1D *hMassX4 = (TH1D*) hMassXY->ProjectionX("hX4",7,8,"");
  TH1D *hMassX5 = (TH1D*) hMassXY->ProjectionX("hX5",9,10,"");
  TH1D *hMassX6 = (TH1D*) hMassXY->ProjectionX("hX6",11,12,"");
  TH1D *hMassX7 = (TH1D*) hMassXY->ProjectionX("hX7",13,14,"");
  TH1D *hMassX8 = (TH1D*) hMassXY->ProjectionX("hX8",15,16,"");
  TH1D *hMassX9 = (TH1D*) hMassXY->ProjectionX("hX9",17,18,"");
  TH1D *hMassX10 = (TH1D*) hMassXY->ProjectionX("hX10",19,20,"");

  TF1 *f2 = new TF1("f2",fTest,0,100,8);
  TF1 *f1 = new TF1("f1","gaus(0)+pol2(3)");
  TF1 *f0 = new TF1("f0","gaus");

  TH1D *hDiv = new TH1D("hDiv","",10,-1.0,1.0);
  TH1D *hY = new TH1D("hY","",10,-1.0,1.0);

  TH1D *hB1 = new TH1D("hB1","",140,1.4,2.8);
  TH1D *hB2 = new TH1D("hB2","",140,1.4,2.8);
  TH1D *hB3 = new TH1D("hB3","",140,1.4,2.8);
  TH1D *hB4 = new TH1D("hB4","",140,1.4,2.8);
  TH1D *hB5 = new TH1D("hB5","",140,1.4,2.8);
  TH1D *hB6 = new TH1D("hB6","",140,1.4,2.8);
  TH1D *hB7 = new TH1D("hB7","",140,1.4,2.8);
  TH1D *hB8 = new TH1D("hB8","",140,1.4,2.8);
  TH1D *hB9 = new TH1D("hB9","",140,1.4,2.8);
  TH1D *hB10 = new TH1D("hB10","",140,1.4,2.8);

  hB1->SetLineWidth(2);
  hB1->SetLineColor(2);
  hB2->SetLineWidth(2);
  hB2->SetLineColor(2);
  hB3->SetLineWidth(2);
  hB3->SetLineColor(2);
  hB4->SetLineWidth(2);
  hB4->SetLineColor(2);
  hB5->SetLineWidth(2);
  hB5->SetLineColor(2);
  hB6->SetLineWidth(2);
  hB6->SetLineColor(2);
  hB7->SetLineWidth(2);
  hB7->SetLineColor(2);
  hB8->SetLineWidth(2);
  hB8->SetLineColor(2);
  hB9->SetLineWidth(2);
  hB9->SetLineColor(2);
  hB10->SetLineWidth(2);
  hB10->SetLineColor(2);


  TH1D *hS1 = new TH1D("hS1","",140,1.4,2.8);
  TH1D *hS2 = new TH1D("hS2","",140,1.4,2.8);
  TH1D *hS3 = new TH1D("hS3","",140,1.4,2.8);
  TH1D *hS4 = new TH1D("hS4","",140,1.4,2.8);
  TH1D *hS5 = new TH1D("hS5","",140,1.4,2.8);
  TH1D *hS6 = new TH1D("hS6","",140,1.4,2.8);
  TH1D *hS7 = new TH1D("hS7","",140,1.4,2.8);
  TH1D *hS8 = new TH1D("hS8","",140,1.4,2.8);
  TH1D *hS9 = new TH1D("hS9","",140,1.4,2.8);
  TH1D *hS10 = new TH1D("hS10","",140,1.4,2.8);

  hS1->SetLineColor(3);
  hS1->SetFillStyle(3001);
  hS1->SetFillColor(3);

  hS2->SetLineColor(3);
  hS2->SetFillStyle(3001);
  hS2->SetFillColor(3);

  hS3->SetLineColor(3);
  hS3->SetFillStyle(3001);
  hS3->SetFillColor(3);

  hS4->SetLineColor(3);
  hS4->SetFillStyle(3001);
  hS4->SetFillColor(3);

  hS5->SetLineColor(3);
  hS5->SetFillStyle(3001);
  hS5->SetFillColor(3);

  hS6->SetLineColor(3);
  hS6->SetFillStyle(3001);
  hS6->SetFillColor(3);

  hS7->SetLineColor(3);
  hS7->SetFillStyle(3001);
  hS7->SetFillColor(3);

  hS8->SetLineColor(3);
  hS8->SetFillStyle(3001);
  hS8->SetFillColor(3);

  hS9->SetLineColor(3);
  hS9->SetFillStyle(3001);
  hS9->SetFillColor(3);

  hS10->SetLineColor(3);
  hS10->SetFillStyle(3001);
  hS10->SetFillColor(3);

  //hMassX1->Rebin(2);



  double fitLow = 1.45;
  double fitHigh = 1.65;
  double fitLow1 = fitLow;
  double fitHigh1 = fitHigh;
  double fitLow2 = fitLow;
  double fitHigh2 = fitHigh;
  double fitLow3 = fitLow;
  double fitHigh3 = fitHigh;
  double fitLow4 = fitLow+0.051;
  double fitHigh4 = fitHigh;
  double fitLow5 = fitLow;
  double fitHigh5 = fitHigh;
  double fitLow6 = fitLow+0.031;
  double fitHigh6 = fitHigh-0.007;
  double fitLow7 = fitLow;
  double fitHigh7 = fitHigh+0.01;
  double fitLow8 = fitLow;
  double fitHigh8 = fitHigh;
  double fitLow9 = fitLow;
  double fitHigh9 = fitHigh-0.05;
  double fitLow10 = fitLow;
  double fitHigh10 = fitHigh;
  //f2->SetParameter(0,20);
  f2->SetParameter(1,1.53);
  f2->SetParLimits(1,1.52,1.54);
  f2->SetParameter(2,0.01);
  f2->SetParLimits(2,0.008,0.01);
  f2->SetParameter(6,1.45);
  f2->SetParameter(7,200);

  hMassX1->Fit("f2","IB","R",fitLow1,fitHigh1);
  hMassX1->GetXaxis()->SetRangeUser(fitLow1,fitHigh1);
  f2->SetParameter(0,0.0);
  hB1->Eval(f2);
  hS1->Add(hMassX1,hB1,1.0,-1.0);
  int binLow = 11;
  int binHigh = 16;
  hS1->GetXaxis()->SetRange(binLow,binHigh);
  double yErr1;
  double yVal1 = hS1->IntegralAndError(binLow,binHigh,yErr1);

  if (iShow == 1){
    hMassX1->Draw("e");
    hB1->Draw("csame");
    hS1->Draw("same");
    cout<<"y1 = "<<yVal1<<" +/- "<<yErr1<<endl;
    return;
  }   


  hMassX2->Fit("f2","IB","R",fitLow2,fitHigh2);
  hMassX2->GetXaxis()->SetRangeUser(fitLow2,fitHigh2);
  f2->SetParameter(0,0.0);
  hB2->Eval(f2);
  hS2->Add(hMassX2,hB2,1.0,-1.0);
  binLow = 12;
  binHigh = 16;
  hS2->GetXaxis()->SetRange(binLow,binHigh);
  double yErr2;
  double yVal2 = hS2->IntegralAndError(binLow,binHigh,yErr2);

  if (iShow == 2){
    hMassX2->Draw("e");
    hB2->Draw("csame");
    hS2->Draw("same");
    cout<<"y2 = "<<yVal2<<" +/- "<<yErr2<<endl;
    return;
  } 


  hMassX3->Fit("f2","IB","R",fitLow3,fitHigh3);
  hMassX3->GetXaxis()->SetRangeUser(fitLow3,fitHigh3);
  f2->SetParameter(0,0.0);
  hB3->Eval(f2);
  hS3->Add(hMassX3,hB3,1.0,-1.0);
  binLow = 11;
  binHigh = 14;
  hS3->GetXaxis()->SetRange(binLow,binHigh);
  double yErr3;
  double yVal3 = hS3->IntegralAndError(binLow,binHigh,yErr3);

  if (iShow == 3){
    hMassX3->Draw("e");
    hB3->Draw("csame");
    hS3->Draw("same");
    cout<<"y3 = "<<yVal3<<" +/- "<<yErr3<<endl;
    return;
  } 
  hMassX4->Fit("f2","IB","R",fitLow4,fitHigh4);
  hMassX4->GetXaxis()->SetRangeUser(fitLow4,fitHigh4);
  f2->SetParameter(0,0.0);
  hB4->Eval(f2);
  hS4->Add(hMassX4,hB4,1.0,-1.0);
  binLow = 13;
  binHigh = 16;
  hS4->GetXaxis()->SetRange(binLow,binHigh);
  double yErr4;
  double yVal4 = hS4->IntegralAndError(binLow,binHigh,yErr4);

  if (iShow == 4){
    hMassX4->Draw("e");
    hB4->Draw("csame");
    hS4->Draw("same");
    cout<<"y4 = "<<yVal4<<" +/- "<<yErr4<<endl;
    return;
  } 

  fitHigh = 1.65;
  f2->SetParameter(1,1.53);
  f2->SetParLimits(1,1.52,1.54);
  f2->SetParameter(2,0.001);
  f2->SetParLimits(2,0.008,0.01);

  hMassX5->Fit("f2","IB","R",fitLow5,fitHigh5);
  hMassX5->GetXaxis()->SetRangeUser(fitLow5,fitHigh5);
  f2->SetParameter(0,0.0);
  hB5->Eval(f2);
  hS5->Add(hMassX5,hB5,1.0,-1.0);
  binLow = 12;
  binHigh = 15;
  hS5->GetXaxis()->SetRange(binLow,binHigh);
  double yErr5;
  double yVal5 = hS5->IntegralAndError(binLow,binHigh,yErr5);

  if (iShow == 5){
    hMassX5->Draw("e");
    hB5->Draw("csame");
    hS5->Draw("same");
    cout<<"y5 = "<<yVal5<<" +/- "<<yErr5<<endl;
    return;
  } 

  fitLow = 1.48;
  fitHigh = 1.59;
  hMassX6->Fit("f2","IB","R",fitLow6,fitHigh6);
  fitHigh = 1.58;

  hMassX6->GetXaxis()->SetRangeUser(fitLow6,fitHigh6);
  f2->SetParameter(0,0.0);
  hB6->Eval(f2);
  hS6->Add(hMassX6,hB6,1.0,-1.0);
  binLow = 12;
  binHigh = 17;
  hS6->GetXaxis()->SetRange(binLow,binHigh);
  double yErr6;
  double yVal6 = hS6->IntegralAndError(binLow,binHigh,yErr6);

  if (iShow == 6){
    hMassX6->Draw("e");
    hB6->Draw("csame");
    hS6->Draw("same");
    cout<<"y6 = "<<yVal6<<" +/- "<<yErr6<<endl;
    return;
  } 

  hMassX7->Fit("f2","IB","R",fitLow7,fitHigh7);
  hMassX7->GetXaxis()->SetRangeUser(fitLow7,fitHigh7);
  f2->SetParameter(0,0.0);
  hB7->Eval(f2);
  hS7->Add(hMassX7,hB7,1.0,-1.0);
  binLow = 12;
  binHigh = 16;
  hS7->GetXaxis()->SetRange(binLow,binHigh);
  double yErr7;
  double yVal7 = hS7->IntegralAndError(binLow,binHigh,yErr7);

  if (iShow == 7){
    hMassX7->Draw("e");
    hB7->Draw("csame");
    hS7->Draw("same");
    cout<<"y7 = "<<yVal7<<" +/- "<<yErr7<<endl;
    return;
  } 

  f2->SetParLimits(1,1.52,1.58);
  //fitHigh = 1.65;
  hMassX8->Fit("f2","IB","R",fitLow8,fitHigh8);
  hMassX8->GetXaxis()->SetRangeUser(fitLow8,fitHigh8);
  f2->SetParameter(0,0.0);
  hB8->Eval(f2);
  hS8->Add(hMassX8,hB8,1.0,-1.0);
  binLow = 12;
  binHigh = 15;
  hS8->GetXaxis()->SetRange(binLow,binHigh);
  double yErr8;
  double yVal8 = hS8->IntegralAndError(binLow,binHigh,yErr8);

  if (iShow == 8){
    hMassX8->Draw("e");
    hB8->Draw("csame");
    hS8->Draw("same");
    cout<<"y8 = "<<yVal8<<" +/- "<<yErr8<<endl;
    return;
  } 

  hMassX9->Fit("f2","IB","R",fitLow9,fitHigh9);
  hMassX9->GetXaxis()->SetRangeUser(fitLow9,fitHigh9);
  f2->SetParameter(0,0.0);
  hB9->Eval(f2);
  hS9->Add(hMassX9,hB9,1.0,-1.0);
  binLow = 12;
  binHigh = 16;
  hS9->GetXaxis()->SetRange(binLow,binHigh);
  double yErr9;
  double yVal9 = hS9->IntegralAndError(binLow,binHigh,yErr9);

  if (iShow == 9){
    hMassX9->Draw("e");
    hB9->Draw("csame");
    hS9->Draw("same");
    cout<<"y9 = "<<yVal9<<" +/- "<<yErr9<<endl;
    return;
  } 

  //fitLow = 1.47;
  //fitHigh = 1.75;
  //hMassX10->Rebin(2);
  //hB10->Rebin(2);
  //hS10->Rebin(2);
  //f2->FixParameter(0,3);
  //f2->FixParameter(1,1.53);
  //f2->FixParameter(2,0.02);
  hMassX10->Fit("f2","IB","R",fitLow10,fitHigh10);
  hMassX10->GetXaxis()->SetRangeUser(fitLow10,fitHigh10);
  f2->SetParameter(0,0.0);
  hB10->Eval(f2);
  hS10->Add(hMassX10,hB10,1.0,-1.0);
  binLow = 12;
  binHigh = 17;
  hS10->GetXaxis()->SetRange(binLow,binHigh);
  double yErr10;
  double yVal10 = hS10->IntegralAndError(binLow,binHigh,yErr10);

  if (iShow == 10){
    hMassX10->Draw("e");
    hB10->Draw("csame");
    hS10->Draw("same");
    cout<<"y10 = "<<yVal10<<" +/- "<<yErr10<<endl;
    return;
  } 
  hY->SetBinContent(1,yVal1);
  hY->SetBinError(1,yErr1);

  hY->SetBinContent(2,yVal2);
  hY->SetBinError(2,yErr2);

  hY->SetBinContent(3,yVal3);
  hY->SetBinError(3,yErr3);

  hY->SetBinContent(4,yVal4);
  hY->SetBinError(4,yErr4);

  hY->SetBinContent(5,yVal5);
  hY->SetBinError(5,yErr5);

  hY->SetBinContent(6,yVal6);
  hY->SetBinError(6,yErr6);

  hY->SetBinContent(7,yVal7);
  hY->SetBinError(7,yErr7);

  hY->SetBinContent(8,yVal8);
  hY->SetBinError(8,yErr8);

  hY->SetBinContent(9,yVal9);
  hY->SetBinError(9,yErr9);

  hY->SetBinContent(10,yVal10);
  hY->SetBinError(10,yErr10);

  int rebinFac = 2;
  hAcc->Rebin(rebinFac);
  //hMassY->Rebin(rebinFac);
  //hDiv->Rebin(rebinFac);

  hDiv->Sumw2();

  //hDiv->Divide(hMassY,hAcc);
  hDiv->Divide(hY,hAcc);

  hDiv->SetMinimum(0);
  
  hDiv->Fit("pol2");
  hDiv->SetMarkerStyle(21);
  hDiv->Draw("e1");



}