Double_t d1half(Int_t m1,Int_t m2,Double_t theta){
  Double_t d1halfVal = -100;
  if (m1 ==  1 && m2 ==  1) d1halfVal =  cos(theta/2);
  if (m1 == -1 && m2 ==  1) d1halfVal = -sin(theta/2);
  if (m1 ==  1 && m2 == -1) d1halfVal =  sin(theta/2);
  if (m1 == -1 && m2 == -1) d1halfVal =  cos(theta/2);
  return d1halfVal;
}
Double_t d3half(Int_t m1,Int_t m2,Double_t theta){
  Double_t d3halfVal = -100;
  if (m1 ==  1 && m2 ==  1) d3halfVal =  (1.0/2)*cos(theta/2)*(3*cos(theta) - 1);
  if (m1 == -1 && m2 ==  1) d3halfVal = -(1.0/2)*sin(theta/2)*(3*cos(theta) + 1);
  if (m1 ==  1 && m2 == -1) d3halfVal =  (1.0/2)*sin(theta/2)*(3*cos(theta) + 1);
  if (m1 == -1 && m2 == -1) d3halfVal =  (1.0/2)*cos(theta/2)*(3*cos(theta) - 1);

  if (m1 ==  3 && m2 ==  1) d3halfVal =  sqrt(3.0)*sin(theta/2)*pow(cos(theta/2),2);
  if (m1 == -3 && m2 ==  1) d3halfVal =  sqrt(3.0)*pow(sin(theta/2),2)*cos(theta/2);
  if (m1 ==  3 && m2 == -1) d3halfVal =  sqrt(3.0)*pow(sin(theta/2),2)*cos(theta/2);
  if (m1 == -3 && m2 == -1) d3halfVal = -sqrt(3.0)*sin(theta/2)*pow(cos(theta/2),2);

  if (m1 ==  1 && m2 ==  3) d3halfVal = -sqrt(3.0)*sin(theta/2)*pow(cos(theta/2),2);
  if (m1 == -1 && m2 ==  3) d3halfVal =  sqrt(3.0)*pow(sin(theta/2),2)*cos(theta/2);
  if (m1 ==  1 && m2 == -3) d3halfVal =  sqrt(3.0)*pow(sin(theta/2),2)*cos(theta/2);
  if (m1 == -1 && m2 == -3) d3halfVal =  sqrt(3.0)*sin(theta/2)*pow(cos(theta/2),2);

  if (m1 ==  3 && m2 ==  3) d3halfVal =  pow(cos(theta/2),3);
  if (m1 == -3 && m2 ==  3) d3halfVal = -pow(sin(theta/2),3);
  if (m1 ==  3 && m2 == -3) d3halfVal =  pow(sin(theta/2),3);
  if (m1 == -3 && m2 == -3) d3halfVal =  pow(cos(theta/2),3);
  return d3halfVal;
}
Double_t d5half(Int_t m1,Int_t m2,Double_t theta){
  Double_t d5halfVal = -100;

  //if (m1 ==  1 && m2 ==  1) d5halfVal =  cos(theta/2)*(2.5*pow(cos(theta) - 1,2) + 4*(cos(theta - 1) + 1);
  //if (m1 == -1 && m2 ==  1) d5halfVal = -sin(theta/2)*(2.5*pow(cos(theta) - 1,2) + 6*(cos(theta - 1) + 3);
  //if (m1 ==  1 && m2 == -1) d5halfVal =  sin(theta/2)*(2.5*pow(cos(theta) - 1,2) + 6*(cos(theta - 1) + 3);
  //if (m1 == -1 && m2 == -1) d5halfVal =  cos(theta/2)*(2.5*pow(cos(theta) - 1,2) + 4*(cos(theta - 1) + 1);

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

  xVal = *x;
  fitVal = 0.0;

  double thetaVal = acos(xVal);
  double a0     = parPtr[0];
  /////////////////////////
  //    j=1/2 and l=0    //
  /////////////////////////
  //j=1/2, m = 1/2, l = 0, s = 1/2 (lambda = 1/2 + lambda = -1/2)
  double a1101Mag  = parPtr[1];
  double phase1101 = parPtr[2];
  double d1101Val  = d1half(1,1,thetaVal)+d1half(1,-1,thetaVal);
  double f1101re   = d1101Val*a1101Mag*cos(phase1101);
  double f1101im   = d1101Val*a1101Mag*sin(phase1101);

  //j=1/2, m = -1/2, l = 0, s = 1/2 (lambda = 1/2 + lambda = -1/2)
  double a1n101Mag  = parPtr[3];
  double phase1n101 = parPtr[4];
  double d1n101Val  = d1half(-1,1,thetaVal)+d1half(-1,-1,thetaVal);
  double f1n101re   = d1n101Val*a1n101Mag*cos(phase1n101);
  double f1n101im   = d1n101Val*a1n101Mag*sin(phase1n101);
  /////////////////////////
  //    j=1/2 and l=1    //
  /////////////////////////
  //j=1/2, m = 1/2, l = 0, s = 1/2 (lambda = 1/2 + lambda = -1/2)
  double a1111Mag  = parPtr[5];
  double phase1111 = parPtr[6];
  double d1111Val  = -sqrt(1.0/3)*d1half(1,1,thetaVal)+sqrt(1.0/3)*d1half(1,-1,thetaVal);
  double f1111re   = d1111Val*a1111Mag*cos(phase1111);
  double f1111im   = d1111Val*a1111Mag*sin(phase1111);

  //j=1/2, m = -1/2, l = 0, s = 1/2 (lambda = 1/2 + lambda = -1/2)
  double a1n111Mag  = parPtr[7];
  double phase1n111 = parPtr[8];
  double d1n111Val  = -sqrt(1.0/3)*d1half(-1,1,thetaVal)+sqrt(1.0/3)*d1half(-1,-1,thetaVal);
  double f1n111re   = d1n111Val*a1n111Mag*cos(phase1n111);
  double f1n111im   = d1n111Val*a1n111Mag*sin(phase1n111);
  /////////////////////////
  //    j=3/2 and l=1    //
  /////////////////////////
  //j=3/2, m = 3/2, l = 1, s = 1/2 (lambda = 1/2 + lambda = -1/2)
  double a3311Mag  = parPtr[9];
  double phase3311 = parPtr[10];
  double d3311Val  = -sqrt(1.0/3)*d3half(3,1,thetaVal)+sqrt(1.0/3)*d3half(3,-1,thetaVal);
  double f3311re   = d3311Val*a3311Mag*cos(phase3311);
  double f3311im   = d3311Val*a3311Mag*sin(phase3311);
  //j=3/2, m = -3/2, l = 1, s = 1/2 (lambda = 1/2 + lambda = -1/2)
  double a3n311Mag  = parPtr[11];
  double phase3n311 = parPtr[12];
  double d3n311Val  = -sqrt(1.0/3)*d3half(-3,1,thetaVal)+sqrt(1.0/3)*d3half(-3,-1,thetaVal);
  double f3n311re   = d3n311Val*a3n311Mag*cos(phase3n311);
  double f3n311im   = d3n311Val*a3n311Mag*sin(phase3n311);
  //j=3/2, m = 1/2, l = 1, s = 1/2 (lambda = 1/2 + lambda = -1/2)
  double a3111Mag  = parPtr[13];
  double phase3111 = parPtr[14];
  double d3111Val  = -sqrt(1.0/3)*d3half(1,1,thetaVal)+sqrt(1.0/3)*d3half(1,-1,thetaVal);
  double f3111re   = d3111Val*a3111Mag*cos(phase3111);
  double f3111im   = d3111Val*a3111Mag*sin(phase3111);
  //j=3/2, m = -1/2, l = 1, s = 1/2 (lambda = 1/2 + lambda = -1/2)
  double a3n111Mag  = parPtr[15];
  double phase3n111 = parPtr[16];
  double d3n111Val  = -sqrt(1.0/3)*d3half(1,1,thetaVal)+sqrt(1.0/3)*d3half(1,-1,thetaVal);
  double f3n111re   = d3n111Val*a3n111Mag*cos(phase3n111);
  double f3n111im   = d3n111Val*a3n111Mag*sin(phase3n111);


  double fre = f1101re + f1n101re + f1111re + f1n111re;
  double fim = f1101im + f1n101im + f1111im + f1n111im;
  fre += f3311re + f3n311re;
  fim += f3311im + f3n311im;
  fre += f3111re + f3n111re;
  fim += f3111im + f3n111im;

  //f1111re = 0.2;
  //fre = f1111re;
  //fim = 0;
  //cout<<"f3311re = "<<f3311re<<endl;
  //cout<<"a3311Mag = "<<a3311Mag<<endl;
  //cout<<"d3half(1,1,thetaVal) = "<<d3half(1,1,thetaVal)<<endl;
  //cout<<"d3half(1,-1,thetaVal) = "<<d3half(1,-1,thetaVal)<<endl;
  
  fitVal = pow(a0,2) + pow(fre,2) + pow(fim,2);

  return fitVal;
}


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 cosThetaGJ(int iShow,int sBin){
  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("hCosThetaGJzVsCasStar_MC_Kin_4K",hMassXY);
  TH2D *hMassXYMC; hFileMC->GetObject("hCosThetaGJzVsCasStar_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");
  TF1 *f3 = new TF1("f3",fWave,-1,1,17);

  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); 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);
  hMassX6->Rebin(2); hB6->Rebin(2); hS6->Rebin(2);
  hMassX7->Rebin(2); hB7->Rebin(2); hS7->Rebin(2);
  hMassX8->Rebin(2); hB8->Rebin(2); hS8->Rebin(2);
  hMassX9->Rebin(2); hB9->Rebin(2); hS9->Rebin(2);
  hMassX10->Rebin(2); hB10->Rebin(2); hS10->Rebin(2);

  //hMassX1->Rebin(2);

  double bFrac = -1.0;

  double fitLow = 1.45;
  double fitHigh = 2.2;
  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;
  double fitLow6 = fitLow+0.04;
  double fitHigh6 = fitHigh;
  double fitLow7 = fitLow;
  double fitHigh7 = fitHigh;
  double fitLow8 = fitLow;
  double fitHigh8 = fitHigh;
  double fitLow9 = fitLow;
  double fitHigh9 = fitHigh;
  double fitLow10 = fitLow;
  double fitHigh10 = fitHigh;
  //f2->SetParameter(0,20);
  f2->FixParameter(1,1.53);
  //f2->SetParLimits(1,1.52,1.54);
  f2->FixParameter(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,bFrac);

  int binLow = sBin;
  int binHigh = binLow;
  if (sBin < 0){ 
    binLow = 6;
    binHigh = 9;
  }
  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;
  } 


  hMassX6->Fit("f2","IB","R",fitLow6,fitHigh6);
  hMassX6->GetXaxis()->SetRangeUser(fitLow6,fitHigh6);
  f2->SetParameter(0,0.0);
  hB6->Eval(f2);
  hS6->Add(hMassX6,hB6,1.0,bFrac);
  //binLow = 6;
  //binHigh = 8;
  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,bFrac);
  //binLow = 6;
  //binHigh = 8;
  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,bFrac);
  //binLow = 6;
  //binHigh = 8;
  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,bFrac);
  //binLow = 6;
  //binHigh = 8;
  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;
  //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,bFrac);
  //binLow = 7;
  //binHigh = 9;
  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");
  //f3->FixParameter(1,0);
  f3->FixParameter(0,0);
  // j=1/2, m=+1/2, l=0, s=1/2
  f3->SetParameter(1,1);
  f3->SetParameter(2,0);
  f3->SetParLimits(2,-3.14,3.14);
  // j=1/2, m=-1/2, l=0, s=1/2
  f3->SetParameter(3,1);
  f3->SetParameter(4,0);
  f3->SetParLimits(4,-3.14,3.14);
  // j=1/2, m=+1/2, l=1, s=1/2
  f3->SetParameter(5,1);
  f3->SetParameter(6,0);
  f3->SetParLimits(6,-3.14,3.14);
  // j=1/2, m=-1/2, l=1, s=1/2
  f3->SetParameter(7,1);
  f3->SetParameter(8,0);
  f3->SetParLimits(8,-3.14,3.14);
  // j=3/2, m=+3/2, l=1, s=1/2
  f3->SetParameter(9,1);
  f3->SetParameter(10,0);
  f3->SetParLimits(10,-3.14,3.14);
  // j=3/2, m=-3/2, l=1, s=1/2
  f3->SetParameter(11,1);
  f3->SetParameter(12,0);
  f3->SetParLimits(12,-3.14,3.14);
  // j=3/2, m=+1/2, l=1, s=1/2
  f3->SetParameter(13,1);
  f3->SetParameter(14,0);
  f3->SetParLimits(14,-3.14,3.14);
  // j=3/2, m=-1/2, l=1, s=1/2
  f3->SetParameter(15,1);
  f3->SetParameter(16,0);
  f3->SetParLimits(16,-3.14,3.14);

  f3->FixParameter(2,0);//Set one phase to zero

  //f3->FixParameter(1,0);
  //f3->FixParameter(3,0);
  //f3->FixParameter(5,0);
  //f3->FixParameter(7,0);
  //f3->FixParameter(9,0);
  //f3->FixParameter(11,0);
  //f3->FixParameter(11,0);
  //f3->FixParameter(11,0);



  hDiv->Fit("f3","IB");
  double aVal = f3->GetParameter(0);
  double aErr = f3->GetParError(0);
  double bVal = f3->GetParameter(2);
  double baFrac = 2*bVal/(3*aVal);
  //cout<<"baFrac = "<<baFrac<<endl;

  cout<<" J    M   L   S   AmpSquared"<<endl;
  cout<<"1/2 +1/2  0  1/2 "<<pow(f3->GetParameter(1),2)<<endl;
  cout<<"1/2 -1/2  0  1/2 "<<pow(f3->GetParameter(3),2)<<endl;
  cout<<"1/2 +1/2  1  1/2 "<<pow(f3->GetParameter(5),2)<<endl;
  cout<<"1/2 -1/2  1  1/2 "<<pow(f3->GetParameter(7),2)<<endl;
  cout<<"3/2 +3/2  1  1/2 "<<pow(f3->GetParameter(9),2)<<endl;
  cout<<"3/2 -3/2  1  1/2 "<<pow(f3->GetParameter(11),2)<<endl;
  cout<<"3/2 +1/2  1  1/2 "<<pow(f3->GetParameter(13),2)<<endl;
  cout<<"3/2 -1/2  1  1/2 "<<pow(f3->GetParameter(15),2)<<endl;

  //hDiv->Fit("pol4","IB");
  hDiv->SetMarkerStyle(21);
  hDiv->Draw("e1");
  f3->SetLineColor(2);
  //f3->Draw("same");

  //f3->SetParameter(1,0);
  //f3->SetParameter(3,0);
  //f3->SetParameter(5,0);
  //f3->SetParameter(7,0);
  //f3->SetParameter(9,0);
  //f3->SetParameter(11,0);
  //f3->SetParameter(13,0);
  //f3->SetParameter(15,0);
  //f3->SetLineColor(2);
  //f3->Draw("same");

  double intErr;
  double intVal = hDiv->IntegralAndError(1,10,intErr);

  if (bVal < 0){
    aVal = aVal - fabs(bVal);
  }
  
  double aFrac = aVal*10/intVal;
  double aFracErr = aErr*10/intVal;
  double binC = hMassX5->GetBinCenter(sBin);
  
  //cout<<"aFrac = "<<aFrac<<" +/- "<<aFracErr<<endl;
  cout<<"binCenter = "<<binC<<endl;
}