cascadeAnalysis/rawYieldFit/XiKS0_SingleBin_Calculator.C

0001 #include <TFile.h>
0002 #include <TTree.h>
0003 #include <TH1F.h>
0004 #include <TH2F.h>
0005 #include <TCanvas.h>
0006 #include <TStyle.h>
0007 #include <TLegend.h>
0008 #include <TPaveText.h>
0009 #include <TLine.h>
0010 #include <TGraphAsymmErrors.h>
0011 #include <RooRealVar.h>
0012 #include <RooGaussian.h>
0013 #include <RooChebychev.h>
0014 #include <RooPolynomial.h>
0015 #include <RooAddPdf.h>
0016 #include <RooDataSet.h>
0017 #include <RooFitResult.h>
0018 #include <RooPlot.h>
0019 #include <RooArgList.h>
0020 #include <TMath.h>
0021 #include <iostream>
0022 #include <fstream>
0023 #include <vector>
0024 #include <string>
0025 #include <cmath>
0026
0027 using namespace RooFit;
0028
0029 std::string getDate()
0030 {
0031     std::time_t t = std::time(0);   // get time now
0032     std::tm* now = std::localtime(&t);
0033
0034     std::stringstream date;
0035     date << (now->tm_mon + 1) << '/'
0036          <<  now->tm_mday << '/'
0037          << (now->tm_year + 1900);
0038     return date.str();
0039 }
0040
0041 template <typename T>
0042 string to_string_with_precision(const T a_value, const int n = 2)
0043 {
0044     ostringstream out;
0045     out.precision(n);
0046     out << fixed << a_value;
0047     return out.str();
0048 }
0049
0050 std::pair<double,double> GetKS0Yield(TTree* inTree, bool savePlots, std::string binType, int binNumber)
0051 {
0052   RooRealVar mass("K_S0_mass", "mass", 0.45, 0.55);
0053   RooDataSet dataSet("K_S0_mass", "data", mass, Import(*inTree));
0054
0055   // DCB Signal
0056   /*
0057   RooRealVar mean_KS0("KS0_mean", "KS0_mean", 0.494, 0.48, 0.50);
0058   RooRealVar sigma_KS0("sigma_KS0", "sigma_KS0", 0.0032, 1e-3, 0.015);
0059   RooRealVar aLow("aL", "aLow", 1.6135); aLow.setConstant(true);
0060   RooRealVar nLow("nL", "nLow", 4.9563); nLow.setConstant(true);
0061   RooRealVar aHigh("aH", "aHigh", 1.3797); aHigh.setConstant(true);
0062   RooRealVar nHigh("nH", "nHigh", 2.9996); nHigh.setConstant(true);
0063   RooRealVar frac("fDCBcore", "fDCBcore", 0.80260); frac.setConstant(true);
0064   RooDoubleCB signal_KS0("signal_KS0", "signal_KS0", mass, mean_KS0, sigma_KS0, aLow, nLow, aHigh, nHigh, frac);
0065   RooRealVar fSig_KS0("fSig_KS0", "fSig_KS0", 0.5*dataSet.numEntries(), 0., 1.2*dataSet.numEntries());
0066   */
0067   // Triple Gaussian Signal
0068
0069   RooRealVar  mean1_KS0("KS0_mean1", "KS0_mean1", 0.49, 0.48, 0.50);
0070   RooRealVar  sigma1_KS0("sigma1_KS0", "sigma1_KS0", 0.01, 0.0005, 0.5);
0071   RooGaussian gaus1("gaus1", "gaus1", mass, mean1_KS0, sigma1_KS0);
0072   RooRealVar  mean2_KS0("KS0_mean2", "KS0_mean2", 0.49, 0.48, 0.50);
0073   RooRealVar  sigma2_KS0("sigma2_KS0", "sigma2_KS0", 0.01, 0.0005, 0.5);
0074   RooGaussian gaus2("gaus2", "gaus2", mass, mean2_KS0, sigma2_KS0);
0075   RooRealVar  mean3_KS0("KS0_mean3", "KS0_mean3", 0.49, 0.48, 0.50);
0076   RooRealVar  sigma3_KS0("sigma3_KS0", "sigma3_KS0", 0.01, 0.0005, 0.5);
0077   RooGaussian gaus3("gaus3", "gaus3", mass, mean3_KS0, sigma3_KS0);
0078   RooRealVar f1("f1","f1",0.6,0.0,1.0);
0079   RooRealVar f2("f2","f2",0.3,0.0,1.0);
0080   RooAddPdf signal_KS0("signal_KS0","signal_KS0",RooArgList(gaus1,gaus2,gaus3),RooArgList(f1,f2));
0081   RooRealVar  fSig_KS0("fSig_KS0", "fSig_KS0", 0.9*dataSet.numEntries(), 0.5*dataSet.numEntries(), 1.2*dataSet.numEntries());
0082
0083
0084   // Linear Background
0085   RooRealVar p1("p1","coeff #1", 0., -100., 100.);
0086   RooPolynomial background("background","background", mass, RooArgList(p1));
0087
0088   // Quadratic Chebyshev Polynomial Background
0089   //RooRealVar c1("c1", "c1", 0.0, -1.0, 1.0);
0090   //RooRealVar c2("c2", "c2", 0.0, -1.0, 1.0);
0091   //RooChebychev background("background","background", mass, RooArgList(c1, c2));
0092
0093   RooRealVar fBkg("fBkg", "fBkg", 0.5*dataSet.numEntries(), 0., dataSet.numEntries());
0094
0095   RooArgList fitModelList(signal_KS0, background), fitFracList(fSig_KS0, fBkg);
0096   RooAddPdf model("model", "model", fitModelList, fitFracList);
0097   RooFitResult* r  = model.fitTo(dataSet, Save(true), Extended(true), PrintLevel(-1));
0098   r->Print();
0099
0100   double fitYield = fSig_KS0.getValV();
0101   double fitYieldErr = fSig_KS0.getError();
0102
0103   if (savePlots)
0104   {
0105     RooPlot* frame = mass.frame(Title(""));       //creating the frame
0106
0107     RooBinning bins(0.45, 0.55);
0108     int nBins = 40;
0109     bins.addUniform(nBins, 0.45, 0.55);
0110
0111     dataSet.plotOn(frame, Binning(bins), XErrorSize(0), DataError(RooAbsData::SumW2));
0112     model.plotOn(frame, Components(background), LineColor(kGray), DrawOption("F"), FillColor(kGray));
0113     model.plotOn(frame, Components(RooArgSet(signal_KS0, background)), LineColor(kAzure+8), DrawOption("F"), FillColor(kAzure+8), MoveToBack());
0114     model.plotOn(frame, LineColor(kBlack));        //plotting the fit onto the frame
0115     dataSet.plotOn(frame, DrawOption("PE1"), Binning(bins),XErrorSize(0), DataError(RooAbsData::SumW2));
0116     RooHist* pull = frame->pullHist();
0117     RooPlot* frame2 = mass.frame(Title(""));
0118     frame2->addPlotable(pull,"PE1");
0119
0120     TCanvas* c = new TCanvas("massFitCanvas", "massFitCanvas",800, 800);
0121
0122     TPad mainPad("mainPad", "mainPad", 0., 0.3, 1., 1.);
0123     mainPad.SetBottomMargin(0.0);
0124
0125     mainPad.Draw();
0126
0127     TPad pullPad("pullPad", "pullPad", 0., 0.0, 1., 0.3);
0128     pullPad.SetBottomMargin(0.5);
0129     pullPad.SetTopMargin(0);
0130
0131     pullPad.Draw();
0132
0133     mainPad.cd();
0134
0135     frame->SetMarkerStyle(kCircle);
0136     frame->SetMarkerSize(0.02);
0137     frame->SetLineWidth(1);
0138     frame->GetXaxis()->SetNdivisions(505);
0139     frame->GetXaxis()->SetTitle("m(#pi^{+}#pi^{-}) [GeV/c^{2}]");
0140     frame->GetXaxis()->SetTitleOffset(0.9);
0141     frame->GetXaxis()->SetTitleFont(42);
0142     frame->GetXaxis()->SetLabelFont(42);
0143     frame->GetYaxis()->SetTitleFont(42);
0144     frame->GetYaxis()->SetLabelFont(42);
0145     double maxVal = frame->GetMaximum();
0146     frame->GetYaxis()->SetRangeUser(1.0, maxVal*1.5);
0147     float binWidth = 1000.*(0.55 - 0.45) / nBins;
0148     string yAxisTitle = "Candidates / (" + to_string_with_precision(binWidth, 0) + " MeV)";
0149     frame->GetYaxis()->SetTitle(yAxisTitle.c_str());
0150     frame->Draw();
0151
0152     //string fitWidthVal_KS0 = to_string_with_precision(1e3 * sigma_KS0.getValV(), 1);
0153     //string fitWidthErr_KS0 = to_string_with_precision(1e3 * sigma_KS0.getError(), 1);
0154     //string fitWidth_KS0 = "#sigma(K_{S}^{0}) = " + fitWidthVal_KS0 + "#pm " + fitWidthErr_KS0 + " MeV";
0155
0156     TPaveText *pt;
0157     pt = new TPaveText(0.52,0.59,0.90,0.90, "NDC"); //For 4 GeV
0158     pt->SetFillColor(0);
0159     pt->SetFillStyle(0);
0160     pt->SetTextFont(42);
0161     TText *pt_LaTex = pt->AddText("#it{#bf{sPHENIX}} Internal");
0162     pt_LaTex = pt->AddText("#it{p}+#it{p} #sqrt{s} = 200 GeV");
0163     //pt_LaTex = pt->AddText(fitWidth_KS0.c_str());
0164     pt->SetBorderSize(0);
0165     pt->DrawClone();
0166
0167     TPaveText *ptDate;
0168     ptDate = new TPaveText(0.67,0.92,1.05,1.00, "NDC");
0169     ptDate->SetFillColor(0);
0170     ptDate->SetFillStyle(0);
0171     ptDate->SetTextFont(42);
0172     std::string compilation_date = getDate();
0173     TText *pt_LaTexDate = ptDate->AddText(compilation_date.c_str());
0174     ptDate->SetBorderSize(0);
0175     ptDate->Draw();
0176     gPad->Modified();
0177
0178     gPad->SetTopMargin(0.08);
0179
0180     frame->Print();
0181     TLegend *legend = new TLegend(0.19,0.70,0.52,0.90);
0182     legend->AddEntry(frame->findObject("h_K_S0_mass"),"Data","PE2");
0183     legend->AddEntry(frame->findObject("model_Norm[K_S0_mass]"),"Fit","L");
0184     legend->AddEntry(frame->findObject("model_Norm[K_S0_mass]_Comp[signal_KS0,background]"),"K_{S}^{0}#rightarrow#pi^{+}#pi^{-}","f");
0185     legend->AddEntry(frame->findObject("model_Norm[K_S0_mass]_Comp[background]"),"Comb. Bkg.","f");
0186     legend->SetFillColor(0);
0187     legend->SetFillStyle(0);
0188     legend->SetBorderSize(0);
0189     legend->SetTextSize(0.05);
0190     legend->Draw();
0191
0192     gPad->Modified();
0193     pullPad.cd();
0194
0195     frame2->SetMarkerStyle(kCircle);
0196     frame2->SetMarkerSize(0.02);
0197     frame2->SetTitle("");
0198     frame2->GetXaxis()->SetNdivisions(505);
0199     frame2->GetXaxis()->SetTitle("m(#pi^{+}#pi^{-}) [GeV/c^{2}]");
0200     frame2->GetXaxis()->SetTitleOffset(0.9);
0201     frame2->GetXaxis()->SetTitleFont(42);
0202     frame2->GetXaxis()->SetTitleSize(0.12);
0203     frame2->GetXaxis()->SetLabelFont(42);
0204     frame2->GetXaxis()->SetLabelSize(0.12);
0205     frame2->GetYaxis()->SetTitle("Pull");
0206     frame2->GetYaxis()->SetTitleOffset(0.55);
0207     frame2->GetYaxis()->SetTitleFont(42);
0208     frame2->GetYaxis()->SetTitleSize(0.12);
0209     frame2->GetYaxis()->SetLabelFont(42);
0210     frame2->GetYaxis()->SetLabelSize(0.12);
0211     frame2->GetYaxis()->SetRangeUser(-6, 6);
0212     frame2->GetYaxis()->SetNdivisions(5);
0213     TF1 *plusThreeLine = new TF1("plusThreeLine",  "pol1", 0.45, 0.55);
0214     plusThreeLine->SetParameters(3, 0);
0215     plusThreeLine->SetLineColor(1);
0216     plusThreeLine->SetLineStyle(2);
0217     plusThreeLine->SetLineWidth(2);
0218     TF1 *zeroLine = new TF1("zeroLine",  "pol1",0.45, 0.55);
0219     zeroLine->SetParameters(0, 0);
0220     zeroLine->SetLineColor(1);
0221     zeroLine->SetLineStyle(2);
0222     zeroLine->SetLineWidth(2);
0223
0224     TF1 *minusThreeLine = new TF1("minusThreeLine", "pol1", 0.45, 0.55);
0225     minusThreeLine->SetParameters(-3, 0);
0226     minusThreeLine->SetLineColor(1);
0227     minusThreeLine->SetLineStyle(2);
0228     minusThreeLine->SetLineWidth(2);
0229     frame2->Draw();
0230     plusThreeLine->Draw("same");
0231     zeroLine->Draw("same");
0232     minusThreeLine->Draw("same");
0233
0234     string outputFile = "./yield_fits/KS0_yield_" + binType + "_bin_" + std::to_string(binNumber) + ".png";
0235     c->SaveAs(outputFile.c_str());
0236   }
0237
0238   return {fitYield, fitYieldErr};
0239 }
0240
0241 std::pair<double,double> GetXiYield(TTree* inTree, bool savePlots, std::string binType, int binNumber)
0242 {
0243   //RooRealVar mass("Ximinus_mass", "mass", 1.3, 1.34);
0244   //RooDataSet dataSet("Ximinus_mass", "data", mass, Import(*inTree));
0245
0246   // DCB Signal
0247   /*
0248   RooRealVar mean_Xi("Xi_mean", "Xi_mean", 1.32, 1.31, 1.33);
0249   RooRealVar sigma_Xi("sigma_Xi", "sigma_Xi", 0.002, 1e-3, 0.005);
0250   RooRealVar aLow("aL", "aLow", 1.882); aLow.setConstant(true);
0251   RooRealVar nLow("nL", "nLow", 1.5159); nLow.setConstant(true);
0252   RooRealVar aHigh("aH", "aHigh", 1.4578); aHigh.setConstant(true);
0253   RooRealVar nHigh("nH", "nHigh", 2.5526); nHigh.setConstant(true);
0254   RooRealVar frac("fDCBcore", "fDCBcore", 0.53753); frac.setConstant(true);
0255   RooDoubleCB signal_Xi("signal_Xi", "signal_Xi", mass, mean_Xi, sigma_Xi, aLow, nLow, aHigh, nHigh, frac);
0256   RooRealVar fSig_Xi("fSig_Xi", "fSig_Xi", 0.5*dataSet.numEntries(), 0., 1.2*dataSet.numEntries());
0257   */
0258
0259   // Triple Gaussian signal
0260   /*
0261   RooRealVar  mean1_Xi("Xi_mean1", "Xi_mean1", 1.32, 1.315, 1.325);
0262   RooRealVar  sigma1_Xi("sigma1_Xi", "sigma1_Xi", 0.001, 0.0005, 0.002);
0263   RooGaussian gaus1("gaus1", "gaus1", mass, mean1_Xi, sigma1_Xi);
0264   //RooRealVar  mean2_Xi("Xi_mean2", "Xi_mean2", 1.32, 1.315, 1.325);
0265   RooRealVar  sigma2_Xi("sigma2_Xi", "sigma2_Xi", 0.002, 0.001, 0.003);
0266   RooGaussian gaus2("gaus2", "gaus2", mass, mean1_Xi, sigma2_Xi);
0267   //RooRealVar  mean3_Xi("Xi_mean3", "Xi_mean3", 1.32, 1.315, 1.325);
0268   //RooRealVar  sigma3_Xi("sigma3_Xi", "sigma3_Xi", 0.007, 0.005, 0.01);
0269   //RooGaussian gaus3("gaus3", "gaus3", mass, mean1_Xi, sigma3_Xi);
0270   RooRealVar f1("f1","f1",0.6,0.0,1.0);
0271   RooRealVar f2("f2","f2",0.3,0.0,1.0);
0272   //RooAddPdf signal_Xi("signal_Xi","signal_Xi",RooArgList(gaus1,gaus2,gaus3),RooArgList(f1,f2));
0273   RooAddPdf signal_Xi("signal_Xi","signal_Xi",RooArgList(gaus1,gaus2),RooArgList(f1,f2));
0274   RooRealVar  fSig_Xi("fSig_Xi", "fSig_Xi", 0.9*dataSet.numEntries(), 0.5*dataSet.numEntries(), 1.2*dataSet.numEntries());
0275   */
0276
0277   /*
0278   // Linear Background
0279   RooRealVar p1("p1","coeff #1", -10, -1000., 0.);
0280   RooPolynomial background("background","background", mass, RooArgList(p1));
0281
0282   // Quadratic Chebyshev Polynomial Background
0283   //RooRealVar c1("c1", "c1", 0.0, -1.0, 1.0);
0284   //RooRealVar c2("c2", "c2", 0.0, -1.0, 1.0);
0285   //RooChebychev background("background","background", mass, RooArgList(c1, c2));
0286
0287   RooRealVar fBkg("fBkg", "fBkg", 0.5*dataSet.numEntries(), 0., dataSet.numEntries());
0288
0289   RooArgList fitModelList(signal_Xi, background), fitFracList(fSig_Xi, fBkg);
0290   RooAddPdf model("model", "model", fitModelList, fitFracList);
0291   RooFitResult* r  = model.fitTo(dataSet, Save(true), Extended(true), PrintLevel(-1));
0292   r->Print();
0293
0294   double fitYield = fSig_Xi.getValV();
0295   double fitYieldErr = fSig_Xi.getError();
0296   */
0297
0298   TH1F* h_Xi_mass = new TH1F("h_Xi_mass", ";m(#Lambda^{0}#pi^{-}) [GeV/c^{2}];Candidates / (1 MeV)", 40, 1.3, 1.34);
0299   inTree->Draw("Ximinus_mass>>h_Xi_mass", "", "goff");
0300
0301   TF1* fitFunc = new TF1("gausPlusLinear","[0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x",1.3, 1.34);
0302   fitFunc->SetLineColor(kRed);
0303   fitFunc->SetParName(0, "Const"); fitFunc->SetParameter(0, 100); fitFunc->SetParLimits(0, 0, 1e4);
0304   fitFunc->SetParName(1, "Mean"); fitFunc->SetParameter(1, 1.32); fitFunc->SetParLimits(1, 1.31, 1.33);
0305   fitFunc->SetParName(2, "Width"); fitFunc->SetParameter(2, 0.002); fitFunc->SetParLimits(2, 0, 0.01);
0306   fitFunc->SetParName(3, "Intercept"); fitFunc->SetParameter(3, 100.); fitFunc->SetParLimits(3, 0., 1e4);
0307   fitFunc->SetParName(4, "Slope");     fitFunc->SetParameter(4, -10.); fitFunc->SetParLimits(4, -1e4, 0.);
0308
0309   double binWidth = h_Xi_mass->GetBinWidth(1);
0310
0311   TCanvas *c1  = new TCanvas("myCanvas", "myCanvas",800,800);
0312
0313   h_Xi_mass->Draw("PE1");
0314   h_Xi_mass->Fit(fitFunc, "ER");
0315
0316   double fitMean = fitFunc->GetParameter(1);
0317   double fitMeanErr = fitFunc->GetParError(1);
0318   double fitWidth = fitFunc->GetParameter(2);
0319   double fitWidthErr = fitFunc->GetParError(2);
0320
0321   std::cout << "Mean : " << fitMean << "+/-" << fitMeanErr << std::endl;
0322   std::cout << "Width : " << fitWidth << "+/-" << fitWidthErr << std::endl;
0323
0324   double A     = fitFunc->GetParameter(0);
0325   double sigma = fitFunc->GetParameter(2);
0326   double Aerr  = fitFunc->GetParError(0);
0327   double sigerr = fitFunc->GetParError(2);
0328
0329   double fitYield    = (A * sigma * sqrt(2 * TMath::Pi()))/binWidth;
0330   double fitYieldErr = sqrt(2 * TMath::Pi()) / binWidth
0331                   * sqrt(pow(sigma * Aerr, 2) + pow(A * sigerr, 2));
0332
0333   string outputFile = "./yield_fits/Xi_yield_" + binType + "_bin_" + std::to_string(binNumber) + ".png";
0334   c1->SaveAs(outputFile.c_str());
0335
0336   if (savePlots)
0337   {
0338     /*
0339     RooPlot* frame = mass.frame(Title(""));       //creating the frame
0340
0341     RooBinning bins(1.3, 1.34);
0342     int nBins = 40;
0343     bins.addUniform(nBins, 1.3, 1.34);
0344
0345     dataSet.plotOn(frame, Binning(bins), XErrorSize(0), DataError(RooAbsData::SumW2));
0346     model.plotOn(frame, Components(background), LineColor(kGray), DrawOption("F"), FillColor(kGray));
0347     model.plotOn(frame, Components(RooArgSet(signal_Xi, background)), LineColor(kAzure+8), DrawOption("F"), FillColor(kAzure+8), MoveToBack());
0348     model.plotOn(frame, LineColor(kBlack));        //plotting the fit onto the frame
0349     dataSet.plotOn(frame, DrawOption("PE1"), Binning(bins),XErrorSize(0), DataError(RooAbsData::SumW2));
0350     RooHist* pull = frame->pullHist();
0351     RooPlot* frame2 = mass.frame(Title(""));
0352     frame2->addPlotable(pull,"PE1");
0353
0354     TCanvas* c = new TCanvas("massFitCanvas", "massFitCanvas",800, 800);
0355
0356     TPad mainPad("mainPad", "mainPad", 0., 0.3, 1., 1.);
0357     mainPad.SetBottomMargin(0.0);
0358
0359     mainPad.Draw();
0360
0361     TPad pullPad("pullPad", "pullPad", 0., 0.0, 1., 0.3);
0362     pullPad.SetBottomMargin(0.5);
0363     pullPad.SetTopMargin(0);
0364
0365     pullPad.Draw();
0366
0367     mainPad.cd();
0368
0369     frame->SetMarkerStyle(kCircle);
0370     frame->SetMarkerSize(0.02);
0371     frame->SetLineWidth(1);
0372     frame->GetXaxis()->SetNdivisions(505);
0373     frame->GetXaxis()->SetTitle("m(#Lambda^{0}#pi^{-}) [GeV/c^{2}]");
0374     frame->GetXaxis()->SetTitleOffset(0.9);
0375     frame->GetXaxis()->SetTitleFont(42);
0376     frame->GetXaxis()->SetLabelFont(42);
0377     frame->GetYaxis()->SetTitleFont(42);
0378     frame->GetYaxis()->SetLabelFont(42);
0379     double maxVal = frame->GetMaximum();
0380     frame->GetYaxis()->SetRangeUser(1.0, maxVal*1.5);
0381     float binWidth = 1000.*(1.34 - 1.30) / nBins;
0382     string yAxisTitle = "Candidates / (" + to_string_with_precision(binWidth, 0) + " MeV)";
0383     frame->GetYaxis()->SetTitle(yAxisTitle.c_str());
0384     frame->Draw();
0385
0386     //string fitWidthVal_Xi = to_string_with_precision(1e3 * sigma_Xi.getValV(), 1);
0387     //string fitWidthErr_Xi = to_string_with_precision(1e3 * sigma_Xi.getError(), 1);
0388     //string fitWidth_Xi = "#sigma(#Xi^{-}) = " + fitWidthVal_Xi + "#pm " + fitWidthErr_Xi + " MeV";
0389
0390     TPaveText *pt;
0391     pt = new TPaveText(0.52,0.59,0.90,0.90, "NDC"); //For 4 GeV
0392     pt->SetFillColor(0);
0393     pt->SetFillStyle(0);
0394     pt->SetTextFont(42);
0395     TText *pt_LaTex = pt->AddText("#it{#bf{sPHENIX}} Internal");
0396     pt_LaTex = pt->AddText("#it{p}+#it{p} #sqrt{s} = 200 GeV");
0397     //pt_LaTex = pt->AddText(fitWidth_Xi.c_str());
0398     pt->SetBorderSize(0);
0399     pt->DrawClone();
0400
0401     TPaveText *ptDate;
0402     ptDate = new TPaveText(0.67,0.92,1.05,1.00, "NDC");
0403     ptDate->SetFillColor(0);
0404     ptDate->SetFillStyle(0);
0405     ptDate->SetTextFont(42);
0406     std::string compilation_date = getDate();
0407     TText *pt_LaTexDate = ptDate->AddText(compilation_date.c_str());
0408     ptDate->SetBorderSize(0);
0409     ptDate->Draw();
0410     gPad->Modified();
0411
0412     gPad->SetTopMargin(0.08);
0413
0414     frame->Print();
0415     TLegend *legend = new TLegend(0.19,0.70,0.52,0.90);
0416     legend->AddEntry(frame->findObject("h_Ximinus_mass"),"Data","PE2");
0417     legend->AddEntry(frame->findObject("model_Norm[Ximinus_mass]"),"Fit","L");
0418     legend->AddEntry(frame->findObject("model_Norm[Ximinus_mass]_Comp[signal_Xi,background]"),"#Xi^{-}#rightarrow#Lambda^{0}#pi^{-}","f");
0419     legend->AddEntry(frame->findObject("model_Norm[Ximinus_mass]_Comp[background]"),"Comb. Bkg.","f");
0420     legend->SetFillColor(0);
0421     legend->SetFillStyle(0);
0422     legend->SetBorderSize(0);
0423     legend->SetTextSize(0.05);
0424     legend->Draw();
0425
0426     gPad->Modified();
0427     pullPad.cd();
0428
0429     frame2->SetMarkerStyle(kCircle);
0430     frame2->SetMarkerSize(0.02);
0431     frame2->SetTitle("");
0432     frame2->GetXaxis()->SetNdivisions(505);
0433     frame2->GetXaxis()->SetTitle("m(#Lambda^{0}#pi^{-}) [GeV/c^{2}]");
0434     frame2->GetXaxis()->SetTitleOffset(0.9);
0435     frame2->GetXaxis()->SetTitleFont(42);
0436     frame2->GetXaxis()->SetTitleSize(0.12);
0437     frame2->GetXaxis()->SetLabelFont(42);
0438     frame2->GetXaxis()->SetLabelSize(0.12);
0439     frame2->GetYaxis()->SetTitle("Pull");
0440     frame2->GetYaxis()->SetTitleOffset(0.55);
0441     frame2->GetYaxis()->SetTitleFont(42);
0442     frame2->GetYaxis()->SetTitleSize(0.12);
0443     frame2->GetYaxis()->SetLabelFont(42);
0444     frame2->GetYaxis()->SetLabelSize(0.12);
0445     frame2->GetYaxis()->SetRangeUser(-6, 6);
0446     frame2->GetYaxis()->SetNdivisions(5);
0447     TF1 *plusThreeLine = new TF1("plusThreeLine",  "pol1", 1.3, 1.34);
0448     plusThreeLine->SetParameters(3, 0);
0449     plusThreeLine->SetLineColor(1);
0450     plusThreeLine->SetLineStyle(2);
0451     plusThreeLine->SetLineWidth(2);
0452     TF1 *zeroLine = new TF1("zeroLine",  "pol1",1.3, 1.34);
0453     zeroLine->SetParameters(0, 0);
0454     zeroLine->SetLineColor(1);
0455     zeroLine->SetLineStyle(2);
0456     zeroLine->SetLineWidth(2);
0457
0458     TF1 *minusThreeLine = new TF1("minusThreeLine", "pol1", 1.3, 1.34);
0459     minusThreeLine->SetParameters(-3, 0);
0460     minusThreeLine->SetLineColor(1);
0461     minusThreeLine->SetLineStyle(2);
0462     minusThreeLine->SetLineWidth(2);
0463     frame2->Draw();
0464     plusThreeLine->Draw("same");
0465     zeroLine->Draw("same");
0466     minusThreeLine->Draw("same");
0467
0468     string outputFile = "./yield_fits/Xi_yield_" + binType + "_bin_" + std::to_string(binNumber) + ".png";
0469     c->SaveAs(outputFile.c_str());
0470     */
0471   }
0472
0473   return {fitYield, fitYieldErr};
0474 }
0475
0476 void XiKS0_SingleBin_Calculator(std::string binType = "pT", int binNumber = 0, bool saveYieldPlots = false)
0477 {
0478   TFile* KS0DataFile = new TFile("/sphenix/user/rosstom/TrackingAnalysis/Physics_Val_TF/data/analysisMacros/KShort_fullDataset_finalCuts_0p2pTCut_rapidity1p0Cut_BCOcut_.root");
0479   TTree* KS0DataTree = (TTree*)KS0DataFile->Get("DecayTree");
0480
0481   TFile* XiDataFile = new TFile("/sphenix/user/rosstom/TrackingAnalysis/Physics_Val_TF/data/analysisMacros/Ximinus_fullDataset_finalCuts_0p2pTCut_rapidity1p0Cut_BCOcut_.root");
0482   TTree* XiDataTree = (TTree*)XiDataFile->Get("DecayTree");
0483
0484   const int nBins_pT = 4;
0485   double pTbins[nBins_pT + 1] = {0.75, 1.07, 1.25, 1.49, 2.0};
0486   const int nBins_phi = 6;
0487   //double phibins[nBins_phi + 1] = {-M_PI,-2.14,-1.37,-0.61,0.35,2.21,M_PI};
0488   double phibins[nBins_phi + 1] = {-M_PI,-2.09439,-1.04719,0.0,1.04721,2.09441,M_PI};
0489   const int nBins_eta = 6;
0490   //double etabins[nBins_eta + 1] = {-1.1,-0.82,-0.56,-0.20,0.30,0.72,1.1};
0491   double etabins[nBins_eta + 1] = {-1.1,-0.733333,-0.366666,0.0,0.366668,0.733335,1.1};
0492   const int nBins_rap = 6;
0493   //double rapbins[nBins_rap + 1] = {-1.0,-0.64,-0.41,-0.15,0.21,0.55,1.0};
0494   double rapbins[nBins_rap + 1] = {-1.0,-0.66,-0.33,0.0,0.33,0.66,1.0};
0495
0496   double KS02pipi_BR_PDG = 0.6920;
0497   double KS02pipi_BR_uncert_PDG = 0.0005;
0498   double Xi2Lambdapi_BR_PDG = 0.99887;
0499   double Xi2Lambdapi_BR_uncert_PDG = 0.00035;
0500   double Lambda2ppi_BR_PDG = 0.641;
0501   double Lambda2ppi_BR_uncert_PDG = 0.005;
0502
0503   std::string geomAccept_dir = "/sphenix/user/rosstom/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots_evaluator/";
0504   std::string cutEff_dir = "/sphenix/user/rosstom/TrackingAnalysis/Physics_Val_TF/sim/pTComp/";
0505   std::string trackAccept_dir = "/sphenix/user/rosstom/analysis/LightFlavorRatios/bachelor_pi_efficiency/analysis/";
0506
0507   std::string tempFileString;
0508
0509   double ratioValue;
0510   double ratioUncertainty;
0511   if (binType == "pT")
0512   {
0513     if (binNumber >= nBins_pT)
0514     {
0515       std::cout << "ERROR: Bin value out of range for this bin type" << std::endl;
0516       return;
0517     }
0518     std::string binCut_KS0;
0519     std::string binCut_Xi;
0520     if (binNumber == nBins_pT - 1)
0521     {
0522       binCut_KS0 = Form("K_S0_pT >= %f && K_S0_pT <= %f", pTbins[binNumber],pTbins[binNumber+1]);
0523       binCut_Xi = Form("Ximinus_pT >= %f && Ximinus_pT <= %f && Lambda0_track_1_charge == track_3_charge", pTbins[binNumber],pTbins[binNumber+1]);
0524     }
0525     else
0526     {
0527       binCut_KS0 = Form("K_S0_pT >= %f && K_S0_pT < %f", pTbins[binNumber],pTbins[binNumber+1]);
0528       binCut_Xi = Form("Ximinus_pT >= %f && Ximinus_pT < %f && Lambda0_track_1_charge == track_3_charge", pTbins[binNumber],pTbins[binNumber+1]);
0529     }
0530
0531     tempFileString = "temp_cut_trees_trash_" + binType + "_" + std::to_string(binNumber) + ".root";
0532     TFile* tempFile = new TFile(tempFileString.c_str(), "RECREATE");
0533     TTree* KS0DataTree_cut = KS0DataTree->CopyTree(binCut_KS0.c_str());
0534     TTree* XiDataTree_cut = XiDataTree->CopyTree(binCut_Xi.c_str());
0535
0536     std::cout << "Extracting KS0 yield" << std::endl;
0537     std::pair<double,double> KS0_yield_uncert = GetKS0Yield(KS0DataTree_cut, saveYieldPlots, binType, binNumber);
0538     std::cout << "KS0 yield = " << KS0_yield_uncert.first << " +/- " << KS0_yield_uncert.second << std::endl;
0539     std::cout << "Extracting Xi yield" << std::endl;
0540     std::pair<double,double> Xi_yield_uncert = GetXiYield(XiDataTree_cut, saveYieldPlots, binType, binNumber);
0541     std::cout << "Xi yield = " << Xi_yield_uncert.first << " +/- " << Xi_yield_uncert.second << std::endl;
0542
0543     tempFile->Close();
0544
0545     // raw ratio value
0546     ratioValue = Xi_yield_uncert.first/(2*KS0_yield_uncert.first);
0547     ratioUncertainty = std::sqrt(std::pow((Xi_yield_uncert.second/(2*KS0_yield_uncert.first)),2) +
0548                                  std::pow(((KS0_yield_uncert.second*Xi_yield_uncert.first)/(2*std::pow(KS0_yield_uncert.first,2))),2));
0549     std::cout << "Raw ratio value : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0550
0551     std::string outfile_name = "Xi_KS0_ratio_TriGauss_Linear_" + binType + "_bin_" + std::to_string(binNumber)  + ".txt";
0552     std::ofstream outFile(outfile_name.c_str());
0553     outFile << KS0_yield_uncert.first << "\n";
0554     outFile << KS0_yield_uncert.second << "\n";
0555     outFile << Xi_yield_uncert.first << "\n";
0556     outFile << Xi_yield_uncert.second << "\n";
0557     outFile.close();
0558
0559     //correct for geometric acceptance and efficiency
0560     std::string geoAccFile = geomAccept_dir + "Ximinus_to_KS0_geometric_acceptance_ratio_onlyPrimaries_pT.root";
0561     TFile* f_geoAccFile = TFile::Open(geoAccFile.c_str(), "READ");
0562     TCanvas* c_geoAccFile = (TCanvas*)f_geoAccFile->Get("myCanvas");
0563     TH1F* h_geoAccFile = (TH1F*)c_geoAccFile->GetPrimitive("Ximinus_inGeo_pT");
0564     double binContent = h_geoAccFile->GetBinContent(binNumber+1);
0565     double binError   = h_geoAccFile->GetBinError(binNumber+1);
0566     ratioValue = ratioValue/binContent;
0567     ratioUncertainty = ratioUncertainty/binContent;
0568     std::cout << "Ratio with geometric acceptance correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0569     f_geoAccFile->Close();
0570
0571     //correct for cut efficiency
0572     std::string cutEffFile = cutEff_dir + "XiKS0ratio_pT_cutEff.root";
0573     TFile* f_cutEffFile = TFile::Open(cutEffFile.c_str(), "READ");
0574     TCanvas* c_cutEffFile = (TCanvas*)f_cutEffFile->Get("myCanvas");
0575     TGraphAsymmErrors* g_cutEffFile = (TGraphAsymmErrors*)c_cutEffFile->GetPrimitive("Graph");
0576     double binCenter;
0577     g_cutEffFile->GetPoint(binNumber, binCenter, binContent);
0578     binError   = g_cutEffFile->GetErrorYhigh(binNumber);
0579     ratioValue = ratioValue/binContent;
0580     ratioUncertainty = ratioUncertainty/binContent;
0581     std::cout << "Ratio with cut efficiency correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0582     f_cutEffFile->Close();
0583
0584     //correct for tracking acceptance and efficiency
0585     std::string trackAccFile = trackAccept_dir + "trackingAcceptanceEfficiencyCorrections_withSystematics.root";
0586     TFile* f_trackAccFile = TFile::Open(trackAccFile.c_str(), "READ");
0587     TH1F* h_trackAccFile = (TH1F*)f_trackAccFile->Get("h_pT_corrections");
0588     binContent = h_trackAccFile->GetBinContent(binNumber+1);
0589     binError   = h_trackAccFile->GetBinError(binNumber+1);
0590     ratioValue = ratioValue/binContent;
0591     ratioUncertainty = ratioUncertainty/binContent;
0592     std::cout << "Ratio with tracking acceptance and efficiency correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0593     f_trackAccFile->Close();
0594
0595     //correct for branching ratios
0596     ratioValue = ratioValue*(KS02pipi_BR_PDG/(Xi2Lambdapi_BR_PDG*Lambda2ppi_BR_PDG));
0597     ratioUncertainty = ratioUncertainty*(KS02pipi_BR_PDG/(Xi2Lambdapi_BR_PDG*Lambda2ppi_BR_PDG));
0598     std::cout << "Ratio with branching ratio correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0599   }
0600   else if (binType == "eta")
0601   {
0602     if (binNumber >= nBins_eta)
0603     {
0604       std::cout << "ERROR: Bin value out of range for this bin type" << std::endl;
0605       return;
0606     }
0607     std::string binCut_KS0;
0608     std::string binCut_Xi;
0609     if (binNumber == nBins_eta - 1)
0610     {
0611       binCut_KS0 = Form("K_S0_pseudorapidity >= %f && K_S0_pseudorapidity <= %f", etabins[binNumber],etabins[binNumber+1]);
0612       binCut_Xi = Form("Ximinus_pseudorapidity >= %f && Ximinus_pseudorapidity <= %f && Lambda0_track_1_charge == track_3_charge", etabins[binNumber],etabins[binNumber+1]);
0613     }
0614     else
0615     {
0616       binCut_KS0 = Form("K_S0_pseudorapidity >= %f && K_S0_pseudorapidity < %f", etabins[binNumber],etabins[binNumber+1]);
0617       binCut_Xi = Form("Ximinus_pseudorapidity >= %f && Ximinus_pseudorapidity < %f && Lambda0_track_1_charge == track_3_charge", etabins[binNumber],etabins[binNumber+1]);
0618     }
0619
0620     tempFileString = "temp_cut_trees_trash_" + binType + "_" + std::to_string(binNumber) + ".root";
0621     TFile* tempFile = new TFile(tempFileString.c_str(), "RECREATE");
0622     TTree* KS0DataTree_cut = KS0DataTree->CopyTree(binCut_KS0.c_str());
0623     TTree* XiDataTree_cut = XiDataTree->CopyTree(binCut_Xi.c_str());
0624
0625     std::cout << "Extracting KS0 yield" << std::endl;
0626     std::pair<double,double> KS0_yield_uncert = GetKS0Yield(KS0DataTree_cut, saveYieldPlots, binType, binNumber);
0627     std::cout << "KS0 yield = " << KS0_yield_uncert.first << " +/- " << KS0_yield_uncert.second << std::endl;
0628     std::cout << "Extracting Xi yield" << std::endl;
0629     std::pair<double,double> Xi_yield_uncert = GetXiYield(XiDataTree_cut, saveYieldPlots, binType, binNumber);
0630     std::cout << "Xi yield = " << Xi_yield_uncert.first << " +/- " << Xi_yield_uncert.second << std::endl;
0631
0632     tempFile->Close();
0633
0634     // raw ratio value
0635     ratioValue = Xi_yield_uncert.first/(2*KS0_yield_uncert.first);
0636     ratioUncertainty = std::sqrt(std::pow((Xi_yield_uncert.second/(2*KS0_yield_uncert.first)),2) +
0637                                  std::pow(((KS0_yield_uncert.second*Xi_yield_uncert.first)/(2*std::pow(KS0_yield_uncert.first,2))),2));
0638     std::cout << "Raw ratio value : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0639
0640     std::string outfile_name = "Xi_KS0_ratio_TriGauss_Linear_" + binType + "_bin_" + std::to_string(binNumber)  + ".txt";
0641     std::ofstream outFile(outfile_name.c_str());
0642     outFile << KS0_yield_uncert.first << "\n";
0643     outFile << KS0_yield_uncert.second << "\n";
0644     outFile << Xi_yield_uncert.first << "\n";
0645     outFile << Xi_yield_uncert.second << "\n";
0646     outFile.close();
0647
0648     //correct for geometric acceptance and efficiency
0649     std::string geoAccFile = geomAccept_dir + "Ximinus_to_KS0_geometric_acceptance_ratio_onlyPrimaries_eta.root";
0650     TFile* f_geoAccFile = TFile::Open(geoAccFile.c_str(), "READ");
0651     TCanvas* c_geoAccFile = (TCanvas*)f_geoAccFile->Get("myCanvas");
0652     TH1F* h_geoAccFile = (TH1F*)c_geoAccFile->GetPrimitive("Ximinus_inGeo_#eta");
0653     double binContent = h_geoAccFile->GetBinContent(binNumber+1);
0654     double binError   = h_geoAccFile->GetBinError(binNumber+1);
0655     ratioValue = ratioValue/binContent;
0656     ratioUncertainty = ratioUncertainty/binContent;
0657     std::cout << "Ratio with geometric acceptance correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0658     f_geoAccFile->Close();
0659
0660     //correct for cut efficiency
0661     std::string cutEffFile = cutEff_dir + "XiKS0ratio_eta_cutEff.root";
0662     TFile* f_cutEffFile = TFile::Open(cutEffFile.c_str(), "READ");
0663     TCanvas* c_cutEffFile = (TCanvas*)f_cutEffFile->Get("myCanvas");
0664     TGraphAsymmErrors* g_cutEffFile = (TGraphAsymmErrors*)c_cutEffFile->GetPrimitive("Graph");
0665     double binCenter;
0666     g_cutEffFile->GetPoint(binNumber, binCenter, binContent);
0667     binError   = g_cutEffFile->GetErrorYhigh(binNumber);
0668     ratioValue = ratioValue/binContent;
0669     ratioUncertainty = ratioUncertainty/binContent;
0670     std::cout << "Ratio with cut efficiency correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0671     f_cutEffFile->Close();
0672
0673     //correct for tracking acceptance and efficiency
0674     std::string trackAccFile = trackAccept_dir + "trackingAcceptanceEfficiencyCorrections_withSystematics.root";
0675     TFile* f_trackAccFile = TFile::Open(trackAccFile.c_str(), "READ");
0676     TH1F* h_trackAccFile = (TH1F*)f_trackAccFile->Get("h_eta_corrections");
0677     binContent = h_trackAccFile->GetBinContent(binNumber+1);
0678     binError   = h_trackAccFile->GetBinError(binNumber+1);
0679     ratioValue = ratioValue/binContent;
0680     ratioUncertainty = ratioUncertainty/binContent;
0681     std::cout << "Ratio with tracking acceptance and efficiency correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0682     f_trackAccFile->Close();
0683
0684     //correct for branching ratios
0685     ratioValue = ratioValue*(KS02pipi_BR_PDG/(Xi2Lambdapi_BR_PDG*Lambda2ppi_BR_PDG));
0686     ratioUncertainty = ratioUncertainty*(KS02pipi_BR_PDG/(Xi2Lambdapi_BR_PDG*Lambda2ppi_BR_PDG));
0687     std::cout << "Ratio with branching ratio correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0688   }
0689   else if (binType == "phi")
0690   {
0691     if (binNumber >= nBins_phi)
0692     {
0693       std::cout << "ERROR: Bin value out of range for this bin type" << std::endl;
0694       return;
0695     }
0696     std::string binCut_KS0;
0697     std::string binCut_Xi;
0698     if (binNumber == nBins_phi - 1)
0699     {
0700       binCut_KS0 = Form("K_S0_phi >= %f && K_S0_phi <= %f", phibins[binNumber],phibins[binNumber+1]);
0701       binCut_Xi = Form("Ximinus_phi >= %f && Ximinus_phi <= %f && Lambda0_track_1_charge == track_3_charge", phibins[binNumber],phibins[binNumber+1]);
0702     }
0703     else
0704     {
0705       binCut_KS0 = Form("K_S0_phi >= %f && K_S0_phi < %f", phibins[binNumber],phibins[binNumber+1]);
0706       binCut_Xi = Form("Ximinus_phi >= %f && Ximinus_phi < %f && Lambda0_track_1_charge == track_3_charge", phibins[binNumber],phibins[binNumber+1]);
0707     }
0708
0709     tempFileString = "temp_cut_trees_trash_" + binType + "_" + std::to_string(binNumber) + ".root";
0710     TFile* tempFile = new TFile(tempFileString.c_str(), "RECREATE");
0711     TTree* KS0DataTree_cut = KS0DataTree->CopyTree(binCut_KS0.c_str());
0712     TTree* XiDataTree_cut = XiDataTree->CopyTree(binCut_Xi.c_str());
0713
0714     std::cout << "Extracting KS0 yield" << std::endl;
0715     std::pair<double,double> KS0_yield_uncert = GetKS0Yield(KS0DataTree_cut, saveYieldPlots, binType, binNumber);
0716     std::cout << "KS0 yield = " << KS0_yield_uncert.first << " +/- " << KS0_yield_uncert.second << std::endl;
0717     std::cout << "Extracting Xi yield" << std::endl;
0718     std::pair<double,double> Xi_yield_uncert = GetXiYield(XiDataTree_cut, saveYieldPlots, binType, binNumber);
0719     std::cout << "Xi yield = " << Xi_yield_uncert.first << " +/- " << Xi_yield_uncert.second << std::endl;
0720
0721     tempFile->Close();
0722
0723     // raw ratio value
0724     ratioValue = Xi_yield_uncert.first/(2*KS0_yield_uncert.first);
0725     ratioUncertainty = std::sqrt(std::pow((Xi_yield_uncert.second/(2*KS0_yield_uncert.first)),2) +
0726                                  std::pow(((KS0_yield_uncert.second*Xi_yield_uncert.first)/(2*std::pow(KS0_yield_uncert.first,2))),2));
0727     std::cout << "Raw ratio value : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0728
0729     std::string outfile_name = "Xi_KS0_ratio_TriGauss_Linear_" + binType + "_bin_" + std::to_string(binNumber)  + ".txt";
0730     std::ofstream outFile(outfile_name.c_str());
0731     outFile << KS0_yield_uncert.first << "\n";
0732     outFile << KS0_yield_uncert.second << "\n";
0733     outFile << Xi_yield_uncert.first << "\n";
0734     outFile << Xi_yield_uncert.second << "\n";
0735     outFile.close();
0736
0737     //correct for geometric acceptance and efficiency
0738     std::string geoAccFile = geomAccept_dir + "Ximinus_to_KS0_geometric_acceptance_ratio_onlyPrimaries_phi.root";
0739     TFile* f_geoAccFile = TFile::Open(geoAccFile.c_str(), "READ");
0740     TCanvas* c_geoAccFile = (TCanvas*)f_geoAccFile->Get("myCanvas");
0741     TH1F* h_geoAccFile = (TH1F*)c_geoAccFile->GetPrimitive("Ximinus_inGeo_#phi");
0742     double binContent = h_geoAccFile->GetBinContent(binNumber+1);
0743     double binError   = h_geoAccFile->GetBinError(binNumber+1);
0744     ratioValue = ratioValue/binContent;
0745     ratioUncertainty = ratioUncertainty/binContent;
0746     std::cout << "Ratio with geometric acceptance correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0747     f_geoAccFile->Close();
0748
0749     //correct for cut efficiency
0750     std::string cutEffFile = cutEff_dir + "XiKS0ratio_phi_cutEff.root";
0751     TFile* f_cutEffFile = TFile::Open(cutEffFile.c_str(), "READ");
0752     TCanvas* c_cutEffFile = (TCanvas*)f_cutEffFile->Get("myCanvas");
0753     TGraphAsymmErrors* g_cutEffFile = (TGraphAsymmErrors*)c_cutEffFile->GetPrimitive("Graph");
0754     double binCenter;
0755     g_cutEffFile->GetPoint(binNumber, binCenter, binContent);
0756     binError   = g_cutEffFile->GetErrorYhigh(binNumber);
0757     ratioValue = ratioValue/binContent;
0758     ratioUncertainty = ratioUncertainty/binContent;
0759     std::cout << "Ratio with cut efficiency correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0760     f_cutEffFile->Close();
0761
0762     //correct for tracking acceptance and efficiency
0763     std::string trackAccFile = trackAccept_dir + "trackingAcceptanceEfficiencyCorrections_withSystematics.root";
0764     TFile* f_trackAccFile = TFile::Open(trackAccFile.c_str(), "READ");
0765     TH1F* h_trackAccFile = (TH1F*)f_trackAccFile->Get("h_phi_corrections");
0766     binContent = h_trackAccFile->GetBinContent(binNumber+1);
0767     binError   = h_trackAccFile->GetBinError(binNumber+1);
0768     ratioValue = ratioValue/binContent;
0769     ratioUncertainty = ratioUncertainty/binContent;
0770     std::cout << "Ratio with tracking acceptance and efficiency correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0771     f_trackAccFile->Close();
0772
0773     //correct for branching ratios
0774     ratioValue = ratioValue*(KS02pipi_BR_PDG/(Xi2Lambdapi_BR_PDG*Lambda2ppi_BR_PDG));
0775     ratioUncertainty = ratioUncertainty*(KS02pipi_BR_PDG/(Xi2Lambdapi_BR_PDG*Lambda2ppi_BR_PDG));
0776     std::cout << "Ratio with branching ratio correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0777   }
0778   else if (binType == "rap")
0779   {
0780     if (binNumber >= nBins_rap)
0781     {
0782       std::cout << "ERROR: Bin value out of range for this bin type" << std::endl;
0783       return;
0784     }
0785     std::string binCut_KS0;
0786     std::string binCut_Xi;
0787     if (binNumber == nBins_rap - 1)
0788     {
0789       binCut_KS0 = Form("K_S0_rapidity >= %f && K_S0_rapidity <= %f", rapbins[binNumber],rapbins[binNumber+1]);
0790       binCut_Xi = Form("Ximinus_rapidity >= %f && Ximinus_rapidity <= %f && Lambda0_track_1_charge == track_3_charge", rapbins[binNumber],rapbins[binNumber+1]);
0791     }
0792     else
0793     {
0794       binCut_KS0 = Form("K_S0_rapidity >= %f && K_S0_rapidity < %f", rapbins[binNumber],rapbins[binNumber+1]);
0795       binCut_Xi = Form("Ximinus_rapidity >= %f && Ximinus_rapidity < %f && Lambda0_track_1_charge == track_3_charge", rapbins[binNumber],rapbins[binNumber+1]);
0796     }
0797
0798     tempFileString = "temp_cut_trees_trash_" + binType + "_" + std::to_string(binNumber) + ".root";
0799     TFile* tempFile = new TFile(tempFileString.c_str(), "RECREATE");
0800     TTree* KS0DataTree_cut = KS0DataTree->CopyTree(binCut_KS0.c_str());
0801     TTree* XiDataTree_cut = XiDataTree->CopyTree(binCut_Xi.c_str());
0802
0803     std::cout << "Extracting KS0 yield" << std::endl;
0804     std::pair<double,double> KS0_yield_uncert = GetKS0Yield(KS0DataTree_cut, saveYieldPlots, binType, binNumber);
0805     std::cout << "KS0 yield = " << KS0_yield_uncert.first << " +/- " << KS0_yield_uncert.second << std::endl;
0806     std::cout << "Extracting Xi yield" << std::endl;
0807     std::pair<double,double> Xi_yield_uncert = GetXiYield(XiDataTree_cut, saveYieldPlots, binType, binNumber);
0808     std::cout << "Xi yield = " << Xi_yield_uncert.first << " +/- " << Xi_yield_uncert.second << std::endl;
0809
0810     tempFile->Close();
0811
0812     // raw ratio value
0813     ratioValue = Xi_yield_uncert.first/(2*KS0_yield_uncert.first);
0814     ratioUncertainty = std::sqrt(std::pow((Xi_yield_uncert.second/(2*KS0_yield_uncert.first)),2) +
0815                                  std::pow(((KS0_yield_uncert.second*Xi_yield_uncert.first)/(2*std::pow(KS0_yield_uncert.first,2))),2));
0816     std::cout << "Raw ratio value : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0817
0818     std::string outfile_name = "Xi_KS0_ratio_TriGauss_Linear_" + binType + "_bin_" + std::to_string(binNumber)  + ".txt";
0819     std::ofstream outFile(outfile_name.c_str());
0820     outFile << KS0_yield_uncert.first << "\n";
0821     outFile << KS0_yield_uncert.second << "\n";
0822     outFile << Xi_yield_uncert.first << "\n";
0823     outFile << Xi_yield_uncert.second << "\n";
0824     outFile.close();
0825
0826     //correct for geometric acceptance and efficiency
0827     std::string geoAccFile = geomAccept_dir + "Ximinus_to_KS0_geometric_acceptance_ratio_onlyPrimaries_rap.root";
0828     TFile* f_geoAccFile = TFile::Open(geoAccFile.c_str(), "READ");
0829     TCanvas* c_geoAccFile = (TCanvas*)f_geoAccFile->Get("myCanvas");
0830     TH1F* h_geoAccFile = (TH1F*)c_geoAccFile->GetPrimitive("Ximinus_inGeo_y");
0831     double binContent = h_geoAccFile->GetBinContent(binNumber+1);
0832     double binError   = h_geoAccFile->GetBinError(binNumber+1);
0833     ratioValue = ratioValue/binContent;
0834     ratioUncertainty = ratioUncertainty/binContent;
0835     std::cout << "Ratio with geometric acceptance correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0836     f_geoAccFile->Close();
0837
0838     //correct for cut efficiency
0839     std::string cutEffFile = cutEff_dir + "XiKS0ratio_y_cutEff.root";
0840     TFile* f_cutEffFile = TFile::Open(cutEffFile.c_str(), "READ");
0841     TCanvas* c_cutEffFile = (TCanvas*)f_cutEffFile->Get("myCanvas");
0842     TGraphAsymmErrors* g_cutEffFile = (TGraphAsymmErrors*)c_cutEffFile->GetPrimitive("Graph");
0843     double binCenter;
0844     g_cutEffFile->GetPoint(binNumber, binCenter, binContent);
0845     binError   = g_cutEffFile->GetErrorYhigh(binNumber);
0846     ratioValue = ratioValue/binContent;
0847     ratioUncertainty = ratioUncertainty/binContent;
0848     std::cout << "Ratio with cut efficiency correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0849     f_cutEffFile->Close();
0850
0851     //correct for tracking acceptance and efficiency
0852     std::string trackAccFile = trackAccept_dir + "trackingAcceptanceEfficiencyCorrections_withSystematics.root";
0853     TFile* f_trackAccFile = TFile::Open(trackAccFile.c_str(), "READ");
0854     TH1F* h_trackAccFile = (TH1F*)f_trackAccFile->Get("h_rap_corrections");
0855     binContent = h_trackAccFile->GetBinContent(binNumber+1);
0856     binError   = h_trackAccFile->GetBinError(binNumber+1);
0857     ratioValue = ratioValue/binContent;
0858     ratioUncertainty = ratioUncertainty/binContent;
0859     std::cout << "Ratio with tracking acceptance and efficiency correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0860     f_trackAccFile->Close();
0861
0862     //correct for branching ratios
0863     ratioValue = ratioValue*(KS02pipi_BR_PDG/(Xi2Lambdapi_BR_PDG*Lambda2ppi_BR_PDG));
0864     ratioUncertainty = ratioUncertainty*(KS02pipi_BR_PDG/(Xi2Lambdapi_BR_PDG*Lambda2ppi_BR_PDG));
0865     std::cout << "Ratio with branching ratio correction : " << ratioValue << "+/-" << ratioUncertainty << std::endl;
0866   }
0867   //std::string outfile_name = "Xi_KS0_ratio_TriGauss_Cheb_" + binType + "_bin_" + std::to_string(binNumber)  + ".txt";
0868   //std::ofstream outFile(outfile_name.c_str());
0869   //outFile << ratioValue << "\n";
0870   //outFile << ratioUncertainty << "\n";
0871   //outFile.close();
0872
0873   gSystem->Unlink(tempFileString.c_str());
0874 }