cascadeAnalysis/finalPlots/finalRatioPlots.C

0001 #include <fstream>
0002 #include <iostream>
0003 #include "stdio.h"
0004
0005 #include <TChain.h>
0006 #include <TFile.h>
0007 #include <TGraphAsymmErrors.h>
0008 #include <TGraphErrors.h>
0009 #include <TLatex.h>
0010 #include <TLegend.h>
0011 #include <TLine.h>
0012 #include <TROOT.h>
0013 #include <TMath.h>
0014 #include <TString.h>
0015 #include <TTree.h>
0016 #include <TVectorD.h>
0017 #include <TVirtualFitter.h>
0018 #include <cassert>
0019 #include <cmath>
0020 #include <vector>
0021
0022 void finalRatioPlots()
0023 {
0024   TFile* rawYieldFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/systematics/rawDataYield.root", "READ");
0025
0026   const std::vector<std::string> kinVars  = {"pT", "phi", "eta", "rap"};
0027   std::map<std::string, std::string> xAxisLabels  = {{"pT","Mother p_{T} [GeV/c]"}, {"phi","Mother #phi [rad]"}, {"eta","Mother #eta"}, {"rap","Mother Rapidity"}};
0028   std::map<std::string, TH1F*> hMap_raw;
0029   for (const auto& kin : kinVars)
0030   {
0031     for (const auto& par : {"KS0", "Xi"})
0032     {
0033       std::string name    = "h_" + kin + "_base_" + par;
0034       TH1F* h             = (TH1F*)rawYieldFile->Get(name.c_str());
0035       h->SetDirectory(0);
0036       hMap_raw[kin + "_" + par] = h;
0037     }
0038   }
0039   rawYieldFile->Close();
0040
0041   TFile* geoAccEffFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots_evaluator/geometricAcceptanceEfficiencyCorrections.root","READ");
0042   std::map<std::string, TH1F*> hMap_geoAccEff;
0043   for (const auto& kin : kinVars)
0044   {
0045     for (const auto& par : {"KS0", "Xi"})
0046     {
0047       std::string name    = "h_geoAccEff_" + kin + "_" + par;
0048       TH1F* h             = (TH1F*)geoAccEffFile->Get(name.c_str());
0049       h->SetDirectory(0);
0050       hMap_geoAccEff[kin + "_" + par] = h;
0051     }
0052   }
0053   geoAccEffFile->Close();
0054
0055   TFile* cutEffFile = TFile::Open("/sphenix/user/rosstom/TrackingAnalysis/Physics_Val_TF/sim/pTComp/cutEfficiencyCorrections.root","READ");
0056   std::map<std::string, TH1F*> hMap_cutEff;
0057   for (const auto& kin : kinVars)
0058   {
0059     for (const auto& par : {"KS0", "Xi"})
0060     {
0061       std::string name    = "h_cutEff_" + kin + "_" + par;
0062       TH1F* h             = (TH1F*)cutEffFile->Get(name.c_str());
0063       h->SetDirectory(0);
0064       hMap_cutEff[kin + "_" + par] = h;
0065     }
0066   }
0067   cutEffFile->Close();
0068
0069   TFile* trackAccFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/bachelor_pi_efficiency/analysis/trackingAcceptanceEfficiencyCorrections_withSystematics.root","READ");
0070   std::map<std::string, TH1F*> hMap_trackAcc;
0071   for (const auto& kin : kinVars)
0072   {
0073     for (const auto& par : {"KS0", "Xi"})
0074     {
0075       std::string name    = "h_trackAcc_" + kin + "_" + par;
0076       TH1F* h             = (TH1F*)trackAccFile->Get(name.c_str());
0077       h->SetDirectory(0);
0078       hMap_trackAcc[kin + "_" + par] = h;
0079     }
0080   }
0081   trackAccFile->Close();
0082
0083   double KS02pipi_BR_PDG = 0.6920;
0084   double KS02pipi_BR_uncert_PDG = 0.0005;
0085   double Xi2Lambdapi_BR_PDG = 0.99887;
0086   double Xi2Lambdapi_BR_uncert_PDG = 0.00035;
0087   double Lambda2ppi_BR_PDG = 0.641;
0088   double Lambda2ppi_BR_uncert_PDG = 0.005;
0089
0090   TFile* geoAccEffSysFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/systematics/geoAccEffSystematicUncertainties.root","READ");
0091   std::map<std::string, TH1F*> hMap_geoAccEffSys;
0092   for (const auto& kin : kinVars)
0093   {
0094     std::string name    = "h_ratio_base_" + kin;
0095     TH1F* h1            = (TH1F*)geoAccEffSysFile->Get(name.c_str());
0096     name                = "h_syst_" + kin;
0097     TH1F* h2            = (TH1F*)geoAccEffSysFile->Get(name.c_str());
0098     h1->SetDirectory(0);
0099     h2->SetDirectory(0);
0100
0101     TH1F* h_syst = (TH1F*)h2->Clone(("h_geoAccEffSys_" + kin).c_str());
0102     h_syst->SetDirectory(0);
0103     h_syst->Divide(h_syst, h1, 1., 1.);
0104
0105     hMap_geoAccEffSys[kin] = h_syst;
0106   }
0107   geoAccEffSysFile->Close();
0108
0109   TFile* cutEffSysFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/systematics/cutEffSystematicUncertainties.root","READ");
0110   std::map<std::string, TH1F*> hMap_cutEffSys;
0111   for (const auto& kin : kinVars)
0112   {
0113     std::string name    = "h_ratio_base_" + kin;
0114     TH1F* h1            = (TH1F*)cutEffSysFile->Get(name.c_str());
0115     name                = "h_syst_" + kin;
0116     TH1F* h2            = (TH1F*)cutEffSysFile->Get(name.c_str());
0117     h1->SetDirectory(0);
0118     h2->SetDirectory(0);
0119
0120     TH1F* h_syst = (TH1F*)h2->Clone(("h_cutEffSys_" + kin).c_str());
0121     h_syst->SetDirectory(0);
0122     h_syst->Divide(h_syst, h1, 1., 1.);
0123
0124     hMap_cutEffSys[kin] = h_syst;
0125   }
0126   cutEffSysFile->Close();
0127
0128   TFile* trackAccSysFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/systematics/trackAccSystematicUncertainties.root","READ");
0129   std::map<std::string, TH1F*> hMap_trackAccSys;
0130   for (const auto& kin : kinVars)
0131   {
0132     std::string name    = "h_ratio_base_" + kin;
0133     TH1F* h1            = (TH1F*)trackAccSysFile->Get(name.c_str());
0134     name                = "h_syst_" + kin;
0135     TH1F* h2            = (TH1F*)trackAccSysFile->Get(name.c_str());
0136     h1->SetDirectory(0);
0137     h2->SetDirectory(0);
0138
0139     TH1F* h_syst = (TH1F*)h2->Clone(("h_trackAccSys_" + kin).c_str());
0140     h_syst->SetDirectory(0);
0141     h_syst->Divide(h_syst, h1, 1., 1.);
0142
0143     hMap_trackAccSys[kin] = h_syst;
0144   }
0145   trackAccSysFile->Close();
0146
0147   TFile* cutVarSysFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/systematics/cutVarianceSystematicUncertainties.root","READ");
0148   std::map<std::string, TH1F*> hMap_cutVarSys;
0149   for (const auto& kin : kinVars)
0150   {
0151     std::string name    = "h_ratio_base_" + kin;
0152     TH1F* h1            = (TH1F*)cutVarSysFile->Get(name.c_str());
0153     name                = "h_syst_" + kin;
0154     TH1F* h2            = (TH1F*)cutVarSysFile->Get(name.c_str());
0155     h1->SetDirectory(0);
0156     h2->SetDirectory(0);
0157
0158     TH1F* h_syst = (TH1F*)h2->Clone(("h_cutVarSys_" + kin).c_str());
0159     h_syst->SetDirectory(0);
0160     h_syst->Divide(h_syst, h1, 1., 1.);
0161
0162     hMap_cutVarSys[kin] = h_syst;
0163   }
0164   cutVarSysFile->Close();
0165
0166   TFile* fitFuncSysFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/systematics/fitFunctionSystematicUncertainties.root","READ");
0167   std::map<std::string, TH1F*> hMap_fitFuncSys;
0168   for (const auto& kin : kinVars)
0169   {
0170     std::string name    = "h_ratio_base_" + kin;
0171     TH1F* h1            = (TH1F*)fitFuncSysFile->Get(name.c_str());
0172     name                = "h_syst_" + kin;
0173     TH1F* h2            = (TH1F*)fitFuncSysFile->Get(name.c_str());
0174     h1->SetDirectory(0);
0175     h2->SetDirectory(0);
0176
0177     TH1F* h_syst = (TH1F*)h2->Clone(("h_fitFuncSys_" + kin).c_str());
0178     h_syst->SetDirectory(0);
0179     h_syst->Divide(h_syst, h1, 1., 1.);
0180
0181     hMap_fitFuncSys[kin] = h_syst;
0182   }
0183   fitFuncSysFile->Close();
0184
0185   for (const auto& kin : kinVars)
0186   {
0187     TH1F* h_KS0 = (TH1F*)hMap_raw[kin + "_KS0"]->Clone(("h_KS0" + kin).c_str());
0188     TH1F* h_Xi = (TH1F*)hMap_raw[kin + "_Xi"]->Clone(("h_Xi" + kin).c_str());
0189     h_KS0->SetDirectory(0);
0190     h_Xi->SetDirectory(0);
0191
0192     int nBins = h_KS0->GetNbinsX();
0193
0194     for (int b = 1; b <= nBins; b++)
0195     {
0196       double KS0_yield    = h_KS0->GetBinContent(b);
0197       double KS0_statErr  = h_KS0->GetBinError(b);
0198       double KS0_geoAcc   = hMap_geoAccEff[kin + "_KS0"]->GetBinContent(b);
0199       double KS0_cutEff   = hMap_cutEff[kin + "_KS0"]->GetBinContent(b);
0200       double KS0_trackAcc = hMap_trackAcc[kin + "_KS0"]->GetBinContent(b);
0201
0202       double KS0_totalCorr  = KS0_geoAcc * KS0_cutEff * KS0_trackAcc * KS02pipi_BR_PDG;
0203       double KS0_corrYield  = KS0_yield / KS0_totalCorr;
0204       double KS0_corrErr    = KS0_statErr / KS0_totalCorr;
0205       h_KS0->SetBinContent(b, KS0_corrYield);
0206       h_KS0->SetBinError(b, KS0_corrErr);
0207
0208       double Xi_yield    = h_Xi->GetBinContent(b);
0209       double Xi_statErr  = h_Xi->GetBinError(b);
0210       double Xi_geoAcc   = hMap_geoAccEff[kin + "_Xi"]->GetBinContent(b);
0211       double Xi_cutEff   = hMap_cutEff[kin + "_Xi"]->GetBinContent(b);
0212       double Xi_trackAcc = hMap_trackAcc[kin + "_Xi"]->GetBinContent(b);
0213
0214       double Xi_totalCorr  = Xi_geoAcc * Xi_cutEff * Xi_trackAcc * Xi2Lambdapi_BR_PDG * Lambda2ppi_BR_PDG;
0215       double Xi_corrYield  = Xi_yield / Xi_totalCorr;
0216       double Xi_corrErr    = Xi_statErr / Xi_totalCorr;
0217       h_Xi->SetBinContent(b, Xi_corrYield);
0218       h_Xi->SetBinError(b, Xi_corrErr);
0219     }
0220
0221     std::string ratioName  = "h_ratio_" + kin;
0222     TH1F* h_ratio = (TH1F*)h_Xi->Clone(ratioName.c_str());
0223     h_ratio->SetDirectory(0);
0224
0225     for (int b = 1; b <= nBins; b++)
0226     {
0227       double Xi_val    = h_Xi->GetBinContent(b);
0228       double Xi_err    = h_Xi->GetBinError(b);
0229       double KS0_val   = h_KS0->GetBinContent(b);
0230       double KS0_err   = h_KS0->GetBinError(b);
0231
0232       double ratio    = Xi_val / (2.0 * KS0_val);
0233
0234       double relErrXi  = Xi_err  / Xi_val;
0235       double relErrKS0 = KS0_err / KS0_val;
0236       double ratioErr  = ratio * sqrt(relErrXi * relErrXi + relErrKS0 * relErrKS0);
0237
0238       h_ratio->SetBinContent(b, ratio);
0239       h_ratio->SetBinError(b, ratioErr);
0240     }
0241
0242     TGraphErrors* g_stat = new TGraphErrors(nBins);
0243     TGraphErrors* g_syst = new TGraphErrors(nBins);
0244     for (int b = 1; b <= nBins; b++)
0245     {
0246       double x        = h_ratio->GetBinCenter(b);
0247       double binWidth = h_ratio->GetBinWidth(b);
0248       double val      = h_ratio->GetBinContent(b);
0249       double statErr  = h_ratio->GetBinError(b);
0250
0251       double geoAccEffSyst = hMap_geoAccEffSys[kin]->GetBinContent(b) * val;
0252       double cutEffSyst    = hMap_cutEffSys[kin]->GetBinContent(b) * val;
0253       double trackAccSyst  = hMap_trackAccSys[kin]->GetBinContent(b) * val;
0254       double cutVarSyst    = hMap_cutVarSys[kin]->GetBinContent(b) * val;
0255       double fitFuncSyst   = hMap_fitFuncSys[kin]->GetBinContent(b) * val;
0256
0257       double totalSyst = sqrt(geoAccEffSyst * geoAccEffSyst +
0258                               cutEffSyst * cutEffSyst    +
0259                               trackAccSyst * trackAccSyst  +
0260                               cutVarSyst * cutVarSyst    +
0261                               fitFuncSyst * fitFuncSyst);
0262
0263       g_stat->SetPoint(b - 1, x, val);
0264       g_stat->SetPointError(b - 1, 0, statErr);
0265
0266       g_syst->SetPoint(b - 1, x, val);
0267       g_syst->SetPointError(b - 1, binWidth / 2.0, totalSyst);
0268     }
0269
0270     g_stat->SetMarkerStyle(20);
0271     g_stat->SetMarkerSize(0.8);
0272     g_stat->SetMarkerColor(kBlack);
0273     g_stat->SetLineColor(kBlack);
0274     g_stat->SetLineWidth(2);
0275
0276     g_syst->SetMarkerStyle(0);
0277     g_syst->SetFillColorAlpha(kBlue, 0.35);
0278     g_syst->SetFillStyle(1001);
0279     g_syst->SetLineColor(kBlue);
0280     g_syst->SetLineWidth(1);
0281
0282     TCanvas* c = new TCanvas(("c_" + kin).c_str(), ("c_" + kin).c_str(), 800, 600);
0283
0284     g_syst->Draw("A E2");
0285     g_stat->Draw("P E same");
0286
0287     std::string axisLabel = xAxisLabels[kin];
0288     g_syst->GetXaxis()->SetTitle(axisLabel.c_str());
0289     g_syst->GetYaxis()->SetTitle("(#Xi^{-} + #bar{#Xi}^{+}) / (2K_{S}^{0})");
0290     g_syst->SetMinimum(0.0);
0291     g_syst->SetMaximum(0.2);
0292
0293     TLegend* leg = new TLegend(0.2, 0.78, 0.6, 0.9);
0294     leg->SetBorderSize(0);
0295     leg->AddEntry(g_stat, "#it{#bf{sPHENIX}} p+p #sqrt{s} = 200 GeV", "lpe");
0296     //leg->AddEntry(g_syst, "Systematic uncertainty",  "f");
0297
0298     if (kin == "pT")
0299     {
0300       const double m0_Xi  = 1.3217;
0301       const double m0_KS0 = 0.4976;
0302
0303       // STAR Fit Values
0304       double C_Xim  = 0.237066; double eC_Xim  = 0.0131962;
0305       double dNdy_Xim = 0.00280512; double edNdy_Xim = 0.000345498;
0306       double n_Xim  = 14.2103; double en_Xim  = 1.95642;
0307
0308       double C_Xip  = 0.235707; double eC_Xip  = 0.0108784;
0309       double dNdy_Xip = 0.00276249; double edNdy_Xip = 0.000266289;
0310       double n_Xip  = 13.9646; double en_Xip  = 1.71798;
0311
0312       double C_KS0  = 0.196422; double eC_KS0  = 0.00376875;
0313       double dNdy_KS0 = 0.131203; double edNdy_KS0 = 0.00698827;
0314       double n_KS0  = 12.2662; double en_KS0  = 0.90061;
0315
0316       // ALICE Fit Values
0317       double C_Xi_LHC  = 0.566197; double eC_Xi_LHC  = 0.00458052;
0318       double dNdy_Xi_LHC = 0.183844; double edNdy_Xi_LHC = 0.00159675;
0319       double n_Xi_LHC  = 10.5974; double en_Xi_LHC  = 0.319702;
0320
0321       double C_KS0_LHC = 0.424953; double eC_KS0_LHC = 0.00987472;
0322       double dNdy_KS0_LHC = 1.73686; double edNdy_KS0_LHC = 0.0161898;
0323       double n_KS0_LHC = 8.15947; double en_KS0_LHC = 0.466203;
0324
0325       // --- Levy evaluation lambda ---
0326       auto levy = [](double pT, double m0, double C, double dNdy, double n) {
0327     return dNdy * (((n-1)*(n-2))/(2*M_PI*n*C*(n*C + m0*(n-2))))
0328          * pow(1 + (std::sqrt(pT*pT + m0*m0) - m0)/(n*C), -n);
0329       };
0330
0331       auto levyErr = [](double pT,   double m0,
0332               double dNdy, double edNdy,
0333               double C,    double eC,
0334               double n,    double en) -> double
0335       {
0336     double mT      = sqrt(pT * pT + m0 * m0);
0337     double u       = (mT - m0) / (n * C);          // (mT-m0)/(nC)
0338     double base    = 1.0 + u;                       // 1 + (mT-m0)/(nC)
0339     double numer   = n * C + m0 * (n - 2);          // nC + m0(n-2)
0340
0341     // Normalization prefactor A and power term P
0342     double A = (n - 1) * (n - 2) / (2 * M_PI * n * C * numer);
0343     double P = pow(base, -n);
0344     double f = dNdy * A * P;
0345
0346     // d/d(dNdy) -- simplest, just f/dNdy
0347     double dfdNdy = f / dNdy;
0348
0349     // d/dC
0350     // Three contributions: from A's 1/C, from A's 1/numer * n, from P's base
0351     double dfdC = f * (
0352       - 1.0 / C
0353       - n * n / numer                               // d/dC of 1/(nC+m0(n-2)) * numer factor
0354       + n * (mT - m0) / (n * C * C * base)          // d/dC of (1 + (mT-m0)/(nC))^-n
0355     );
0356
0357         // d/dn
0358     // Contributions from A's (n-1), (n-2), n, numer, and from P
0359     double dPdn = P * (-(mT - m0) / (n * C * base) - log(base));  // d/dn of base^(-n)
0360     double dAdn = A * (
0361       1.0 / (n - 1)
0362       + 1.0 / (n - 2)
0363       - 1.0 / n
0364       - (C + m0) / numer                            // d/dn of nC + m0(n-2) = C + m0
0365     );
0366     double dfdn = dNdy * (dAdn * P + A * dPdn);
0367
0368     return sqrt(pow(dfdNdy * edNdy, 2) +
0369             pow(dfdC   * eC,    2) +
0370             pow(dfdn   * en,    2));
0371       };
0372
0373       // --- Build ratio and uncertainty band ---
0374       const int    nPts  = 500;
0375       const double pTmin = 0.5;
0376       const double pTmax = 4.0;
0377       const double step  = (pTmax - pTmin) / nPts;
0378
0379       TGraph*      g_ratio    = new TGraph(nPts);
0380       TGraphErrors* g_band    = new TGraphErrors(nPts);
0381       TGraph*      g_ratio_LHC    = new TGraph(nPts);
0382       TGraphErrors* g_band_LHC    = new TGraphErrors(nPts);
0383
0384       for (int i = 0; i < nPts; i++)
0385       {
0386     double pT = pTmin + i * step;
0387
0388     double fXim = levy(pT, m0_Xi,  C_Xim,  dNdy_Xim,  n_Xim);
0389     double fXip = levy(pT, m0_Xi,  C_Xip,  dNdy_Xip,  n_Xip);
0390     double fKS0 = levy(pT, m0_KS0, C_KS0,  dNdy_KS0,  n_KS0);
0391
0392         double fXi_LHC = levy(pT, m0_Xi,  C_Xi_LHC,  dNdy_Xi_LHC,  n_Xi_LHC);
0393         double fKS0_LHC = levy(pT, m0_KS0,  C_KS0_LHC,  dNdy_KS0_LHC,  n_KS0_LHC);
0394
0395     double numerator   = fXim + fXip;
0396     double denominator = 2.0 * fKS0;
0397
0398         double numerator_LHC   = fXi_LHC;
0399         double denominator_LHC = 2.0 * fKS0_LHC;
0400
0401     if (denominator <= 0) continue;
0402     double ratio = numerator / denominator;
0403     double ratio_LHC = numerator_LHC / denominator_LHC;
0404
0405         // --- Error propagation on ratio R = (Xim + Xip) / (2*KS0) ---
0406     // dR/d(Xim) =  1 / denom
0407     // dR/d(Xip) =  1 / denom
0408     // dR/d(KS0) = -numerator / denom^2 * 2 = -R / KS0_val
0409     double eXim = levyErr(pT, m0_Xi,  C_Xim, eC_Xim, dNdy_Xim, edNdy_Xim, n_Xim, en_Xim);
0410     double eXip = levyErr(pT, m0_Xi,  C_Xip, eC_Xip, dNdy_Xip, edNdy_Xip, n_Xip, en_Xip);
0411     double eKS0 = levyErr(pT, m0_KS0, C_KS0, eC_KS0, dNdy_KS0, edNdy_KS0, n_KS0, en_KS0);
0412     double eXi_LHC = levyErr(pT, m0_Xi, C_Xi_LHC, eC_Xi_LHC, dNdy_Xi_LHC, edNdy_Xi_LHC, n_Xi_LHC, en_Xi_LHC);
0413     double eKS0_LHC = levyErr(pT, m0_KS0, C_KS0_LHC, eC_KS0_LHC, dNdy_KS0_LHC, edNdy_KS0_LHC, n_KS0_LHC, en_KS0_LHC);
0414
0415     double ratioErr = sqrt(
0416       pow(eXim / denominator,       2) +
0417       pow(eXip / denominator,       2) +
0418       pow(ratio * eKS0 / fKS0,      2)
0419     );
0420
0421         double ratioErr_LHC = sqrt(
0422       pow(eXi_LHC / denominator_LHC,       2) +
0423       pow(ratio_LHC * eKS0_LHC / fKS0_LHC,      2)
0424     );
0425
0426     g_ratio->SetPoint(i,      pT, ratio);
0427     g_band->SetPoint(i,       pT, ratio);
0428     g_band->SetPointError(i,  0,  ratioErr);
0429     g_ratio_LHC->SetPoint(i,      pT, ratio_LHC);
0430     g_band_LHC->SetPoint(i,       pT, ratio_LHC);
0431     g_band_LHC->SetPointError(i,  0,  ratioErr_LHC);
0432       }
0433
0434       // Draw band first so curve sits on top
0435       g_band->SetFillColorAlpha(kRed, 0.3);
0436       g_band->SetFillStyle(1001);
0437       g_band->SetLineColorAlpha(kRed, 0.0); // hide band outline
0438       g_band->Draw("E3 same");                 // E3 = filled uncertainty band
0439
0440       g_ratio->SetLineColor(kRed);
0441       g_ratio->SetLineWidth(2);
0442       g_ratio->Draw("L same");
0443
0444       leg->AddEntry(g_ratio, "STAR PRC 75 (2007) 064901", "l");
0445
0446       g_band_LHC->SetFillColorAlpha(kGreen, 0.3);
0447       g_band_LHC->SetFillStyle(1001);
0448       g_band_LHC->SetLineColorAlpha(kGreen, 0.0); // hide band outline
0449       g_band_LHC->Draw("E3 same");                 // E3 = filled uncertainty band
0450
0451       g_ratio_LHC->SetLineColor(kGreen);
0452       g_ratio_LHC->SetLineWidth(2);
0453       g_ratio_LHC->Draw("L same");
0454
0455       leg->AddEntry(g_ratio_LHC, "ALICE EPJC 81 (2021) 256", "l");
0456
0457       for (int q = 0; q < 4; ++q)
0458       {
0459         double x, y;
0460         g_stat->GetPoint(q, x, y);
0461         double stat_err = g_stat->GetErrorY(q);
0462         double syst_err = g_syst->GetErrorY(q);
0463         double fXim = levy(x, m0_Xi,  C_Xim,  dNdy_Xim,  n_Xim);
0464         double fXip = levy(x, m0_Xi,  C_Xip,  dNdy_Xip,  n_Xip);
0465         double fKS0 = levy(x, m0_KS0, C_KS0,  dNdy_KS0,  n_KS0);
0466         double fXi_LHC = levy(x, m0_Xi,  C_Xi_LHC,  dNdy_Xi_LHC,  n_Xi_LHC);
0467         double fKS0_LHC = levy(x, m0_KS0,  C_KS0_LHC,  dNdy_KS0_LHC,  n_KS0_LHC);
0468         double numerator   = fXim + fXip;
0469         double denominator = 2.0 * fKS0;
0470         double numerator_LHC   = fXi_LHC;
0471         double denominator_LHC = 2.0 * fKS0_LHC;
0472         double ratio = numerator / denominator;
0473         double ratio_LHC = numerator_LHC / denominator_LHC;
0474         double eXim = levyErr(x, m0_Xi,  C_Xim, eC_Xim, dNdy_Xim, edNdy_Xim, n_Xim, en_Xim);
0475         double eXip = levyErr(x, m0_Xi,  C_Xip, eC_Xip, dNdy_Xip, edNdy_Xip, n_Xip, en_Xip);
0476         double eKS0 = levyErr(x, m0_KS0, C_KS0, eC_KS0, dNdy_KS0, edNdy_KS0, n_KS0, en_KS0);
0477         double eXi_LHC = levyErr(x, m0_Xi, C_Xi_LHC, eC_Xi_LHC, dNdy_Xi_LHC, edNdy_Xi_LHC, n_Xi_LHC, en_Xi_LHC);
0478         double eKS0_LHC = levyErr(x, m0_KS0, C_KS0_LHC, eC_KS0_LHC, dNdy_KS0_LHC, edNdy_KS0_LHC, n_KS0_LHC, en_KS0_LHC);
0479         double ratioErr = sqrt(
0480           pow(eXim / denominator,       2) +
0481           pow(eXip / denominator,       2) +
0482           pow(ratio * eKS0 / fKS0,      2)
0483         );
0484         double ratioErr_LHC = sqrt(
0485           pow(eXi_LHC / denominator_LHC,       2) +
0486           pow(ratio_LHC * eKS0_LHC / fKS0_LHC,      2)
0487         );
0488         double sigma = std::sqrt(pow(stat_err,2) + pow(syst_err,2) + pow(ratioErr,2));
0489         double pull = std::abs(ratio - y)/sigma;
0490         double sigma_LHC = std::sqrt(pow(stat_err,2) + pow(syst_err,2) + pow(ratioErr_LHC,2));
0491         double pull_LHC = std::abs(ratio_LHC - y)/sigma_LHC;
0492         std::cout << "pT: " << x << ", pull_STAR = " << pull << ", pull_LHC = " << pull_LHC << std::endl;
0493       }
0494     }
0495     leg->Draw();
0496
0497     TPaveText *pt;
0498     pt = new TPaveText(0.75,0.8,0.95,1., "NDC");
0499     pt->SetFillColor(0);
0500     pt->SetFillStyle(0);
0501     pt->SetTextFont(42);
0502     string label = "#it{#bf{sPHENIX}} Internal";
0503     TText *pt_LaTex = pt->AddText(label.c_str());
0504     pt->SetBorderSize(0);
0505     pt->Draw();
0506     gPad->Modified();
0507
0508     //c->SaveAs(("finalRatio_" + kin + ".png").c_str());
0509   }
0510 }