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 cosThetaGJ22(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_4",hMassXY);
  TH2D *hMassXYMC; hFileMC->GetObject("hCosThetaGJyVsCasStar_MC_Kin_4",hMassXYMC);

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

  TH1D *hMassX1 = (TH1D*) hMassXY->ProjectionX("hX1",1,4,"");
  TH1D *hMassX2 = (TH1D*) hMassXY->ProjectionX("hX2",5,8,"");
  TH1D *hMassX3 = (TH1D*) hMassXY->ProjectionX("hX3",9,12,"");
  TH1D *hMassX4 = (TH1D*) hMassXY->ProjectionX("hX4",13,16,"");
  TH1D *hMassX5 = (TH1D*) hMassXY->ProjectionX("hX5",17,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","",5,-1.0,1.0);
  TH1D *hY = new TH1D("hY","",5,-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);

  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);


  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);

  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);


  //hMassX1->Rebin(2); hB1->Rebin(2); hS1->Rebin(2);
  //hMassX2->Rebin(2); hB2->Rebin(2); hS2->Rebin(2);
  //hMassX3->Rebin(2); hB3->Rebin(2); hS3->Rebin(2);
  //hMassX4->Rebin(2); hB4->Rebin(2); hS4->Rebin(2);
  //hMassX5->Rebin(2); hB5->Rebin(2); hS5->Rebin(2);

  //hMassX1->Rebin(2);



  double fitLow = 1.4;
  double fitHigh = 1.75;
  double fitLow1 = fitLow;
  double fitHigh1 = fitHigh;
  double fitLow2 = fitLow;
  double fitHigh2 = fitHigh;
  double fitLow3 = fitLow;
  double fitHigh3 = fitHigh;
  double fitLow4 = fitLow;
  double fitHigh4 = fitHigh;
  double fitLow5 = fitLow;
  double fitHigh5 = fitHigh;

  //f2->SetParameter(0,20);
  f2->FixParameter(1,1.53);
  f2->SetParLimits(1,1.5,1.55);
  f2->SetParameter(2,0.01);
  f2->SetParLimits(2,0.008,0.013);
  f2->SetParameter(6,1.45);
  f2->SetParameter(7,200);

  double bFrac = -0;

  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,bFrac);
  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,bFrac);
  //binLow = 6;
  //binHigh = 8;
  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,bFrac);
  //binLow = 6;
  //binHigh = 8;
  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,bFrac);
  //binLow = 7;
  //binHigh = 9;
  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,bFrac);
  //binLow = 6;
  //binHigh = 8;
  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;
  } 


  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);

  int rebinFac = 4;
  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");



}