File indexing completed on 2026-07-16 08:11:23
0001 void CalculateCorrectedYield_VaryCorrection(std::string inputFile)
0002 {
0003 enum CorrIdx
0004 {
0005 kPt = 0,
0006 kPhi = 1,
0007 kEta = 2,
0008 kRap = 3,
0009 kPt_min = 4,
0010 kPhi_min = 5,
0011 kEta_min = 6,
0012 kRap_min = 7,
0013 kPt_max = 8,
0014 kPhi_max = 9,
0015 kEta_max = 10,
0016 kRap_max = 11,
0017 kNCorrTypes
0018 };
0019
0020 std::vector<TH1F*> h_geoAccEff_KS0(kNCorrTypes, nullptr);
0021 std::vector<TH1F*> h_geoAccEff_Xi(kNCorrTypes, nullptr);
0022 std::vector<TH1F*> h_cutEff_KS0(kNCorrTypes, nullptr);
0023 std::vector<TH1F*> h_cutEff_Xi(kNCorrTypes, nullptr);
0024 std::vector<TH1F*> h_trackAcc_KS0(kNCorrTypes, nullptr);
0025 std::vector<TH1F*> h_trackAcc_Xi(kNCorrTypes, nullptr);
0026
0027 const std::vector<std::string> corrNames = {"pT", "phi", "eta", "rap"};
0028
0029 TFile* correctionFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots_evaluator/geometricAcceptanceEfficiencyCorrections.root", "READ");
0030 for (int i = 0; i < 4; i++)
0031 {
0032 TH1F* h_KS0 = (TH1F*)correctionFile->Get(("h_geoAccEff_" + corrNames[i] + "_KS0").c_str());
0033 h_geoAccEff_KS0[i] = (TH1F*)h_KS0->Clone(("h_geoAccEff_" + corrNames[i] + "_KS0_copy").c_str());
0034 h_geoAccEff_KS0[i]->SetDirectory(0);
0035
0036 TH1F* h_Xi = (TH1F*)correctionFile->Get(("h_geoAccEff_" + corrNames[i] + "_Xi").c_str());
0037 h_geoAccEff_Xi[i] = (TH1F*)h_Xi->Clone(("h_geoAccEff_" + corrNames[i] + "_Xi_copy").c_str());
0038 h_geoAccEff_Xi[i]->SetDirectory(0);
0039 }
0040 correctionFile->Close();
0041
0042 correctionFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots_evaluator/geometricAcceptanceEfficiencyCorrections_pTref0p3.root", "READ");
0043 for (int i = 0; i < 4; i++)
0044 {
0045 TH1F* h_KS0 = (TH1F*)correctionFile->Get(("h_geoAccEff_" + corrNames[i] + "_KS0").c_str());
0046 h_geoAccEff_KS0[i+4] = (TH1F*)h_KS0->Clone(("h_geoAccEff_" + corrNames[i] + "_KS0_min").c_str());
0047 h_geoAccEff_KS0[i+4]->SetDirectory(0);
0048
0049 TH1F* h_Xi = (TH1F*)correctionFile->Get(("h_geoAccEff_" + corrNames[i] + "_Xi").c_str());
0050 h_geoAccEff_Xi[i+4] = (TH1F*)h_Xi->Clone(("h_geoAccEff_" + corrNames[i] + "_Xi_min").c_str());
0051 h_geoAccEff_Xi[i+4]->SetDirectory(0);
0052
0053 h_geoAccEff_KS0[i+8] = (TH1F*)h_KS0->Clone(("h_geoAccEff_" + corrNames[i] + "_KS0_max").c_str());
0054 h_geoAccEff_KS0[i+8]->SetDirectory(0);
0055
0056 h_geoAccEff_Xi[i+8] = (TH1F*)h_Xi->Clone(("h_geoAccEff_" + corrNames[i] + "_Xi_max").c_str());
0057 h_geoAccEff_Xi[i+8]->SetDirectory(0);
0058 }
0059 correctionFile->Close();
0060
0061 correctionFile = TFile::Open("/sphenix/user/rosstom/TrackingAnalysis/Physics_Val_TF/sim/pTComp/cutEfficiencyCorrections.root", "READ");
0062 for (int i = 0; i < 4; i++)
0063 {
0064 TH1F* h_KS0 = (TH1F*)correctionFile->Get(("h_cutEff_" + corrNames[i] + "_KS0").c_str());
0065 h_cutEff_KS0[i] = (TH1F*)h_KS0->Clone(("h_cutEff_" + corrNames[i] + "_KS0_copy").c_str());
0066 h_cutEff_KS0[i]->SetDirectory(0);
0067
0068 TH1F* h_Xi = (TH1F*)correctionFile->Get(("h_cutEff_" + corrNames[i] + "_Xi").c_str());
0069 h_cutEff_Xi[i] = (TH1F*)h_Xi->Clone(("h_cutEff_" + corrNames[i] + "_Xi_copy").c_str());
0070 h_cutEff_Xi[i]->SetDirectory(0);
0071
0072 h_KS0 = (TH1F*)correctionFile->Get(("h_cutEff_" + corrNames[i] + "_KS0_max").c_str());
0073 h_cutEff_KS0[i+8] = (TH1F*)h_KS0->Clone(("h_cutEff_" + corrNames[i] + "_KS0_max_copy").c_str());
0074 h_cutEff_KS0[i+8]->SetDirectory(0);
0075
0076 h_Xi = (TH1F*)correctionFile->Get(("h_cutEff_" + corrNames[i] + "_Xi_max").c_str());
0077 h_cutEff_Xi[i+8] = (TH1F*)h_Xi->Clone(("h_cutEff_" + corrNames[i] + "_Xi_max_copy").c_str());
0078 h_cutEff_Xi[i+8]->SetDirectory(0);
0079 }
0080 correctionFile->Close();
0081
0082 correctionFile = TFile::Open("/sphenix/user/rosstom/analysis/LightFlavorRatios/bachelor_pi_efficiency/analysis/trackingAcceptanceEfficiencyCorrections_withSystematics.root", "READ");
0083 std::string tailMod = "";
0084 for (int i = 0; i < kNCorrTypes; i++)
0085 {
0086 if (i == 4)
0087 {
0088 tailMod = "_high";
0089 }
0090 else if (i == 8)
0091 {
0092 tailMod = "_low";
0093 }
0094 TH1F* h_KS0 = (TH1F*)correctionFile->Get(("h_trackAcc_" + corrNames[i%4] + "_KS0" + tailMod).c_str());
0095 h_trackAcc_KS0[i] = (TH1F*)h_KS0->Clone(("h_trackAcc_" + corrNames[i%4] + "_KS0" + tailMod + "_copy").c_str());
0096 h_trackAcc_KS0[i]->SetDirectory(0);
0097
0098 TH1F* h_Xi = (TH1F*)correctionFile->Get(("h_trackAcc_" + corrNames[i%4] + "_Xi" + tailMod).c_str());
0099 h_trackAcc_Xi[i] = (TH1F*)h_Xi->Clone(("h_trackAcc_" + corrNames[i%4] + "_Xi" + tailMod + "_copy").c_str());
0100 h_trackAcc_Xi[i]->SetDirectory(0);
0101 }
0102 correctionFile->Close();
0103
0104 double KS02pipi_BR_PDG = 0.6920;
0105 double KS02pipi_BR_uncert_PDG = 0.0005;
0106 double Xi2Lambdapi_BR_PDG = 0.99887;
0107 double Xi2Lambdapi_BR_uncert_PDG = 0.00035;
0108 double Lambda2ppi_BR_PDG = 0.641;
0109 double Lambda2ppi_BR_uncert_PDG = 0.005;
0110
0111 TFile* inputRawYield = new TFile(inputFile.c_str());
0112 const std::vector<std::string> variants = {"", "_max", "_min"};
0113
0114 const std::vector<std::string> particles = {"KS0", "Xi"};
0115 const std::vector<std::string> diffNames = {"pT", "phi", "eta", "rap"};
0116
0117 std::map<std::string, TH1F*> corrMap_geoAccEff =
0118 {
0119 {"KS0_pT", h_geoAccEff_KS0[kPt]},
0120 {"KS0_phi", h_geoAccEff_KS0[kPhi]},
0121 {"KS0_eta", h_geoAccEff_KS0[kEta]},
0122 {"KS0_rap", h_geoAccEff_KS0[kRap]},
0123 {"Xi_pT", h_geoAccEff_Xi[kPt]},
0124 {"Xi_phi", h_geoAccEff_Xi[kPhi]},
0125 {"Xi_eta", h_geoAccEff_Xi[kEta]},
0126 {"Xi_rap", h_geoAccEff_Xi[kRap]},
0127 {"KS0_pT_min", h_geoAccEff_KS0[kPt_min]},
0128 {"KS0_phi_min", h_geoAccEff_KS0[kPhi_min]},
0129 {"KS0_eta_min", h_geoAccEff_KS0[kEta_min]},
0130 {"KS0_rap_min", h_geoAccEff_KS0[kRap_min]},
0131 {"Xi_pT_min", h_geoAccEff_Xi[kPt_min]},
0132 {"Xi_phi_min", h_geoAccEff_Xi[kPhi_min]},
0133 {"Xi_eta_min", h_geoAccEff_Xi[kEta_min]},
0134 {"Xi_rap_min", h_geoAccEff_Xi[kRap_min]},
0135 {"KS0_pT_max", h_geoAccEff_KS0[kPt_max]},
0136 {"KS0_phi_max", h_geoAccEff_KS0[kPhi_max]},
0137 {"KS0_eta_max", h_geoAccEff_KS0[kEta_max]},
0138 {"KS0_rap_max", h_geoAccEff_KS0[kRap_max]},
0139 {"Xi_pT_max", h_geoAccEff_Xi[kPt_max]},
0140 {"Xi_phi_max", h_geoAccEff_Xi[kPhi_max]},
0141 {"Xi_eta_max", h_geoAccEff_Xi[kEta_max]},
0142 {"Xi_rap_max", h_geoAccEff_Xi[kRap_max]},
0143 };
0144
0145 std::map<std::string, TH1F*> corrMap_cutEff =
0146 {
0147 {"KS0_pT", h_cutEff_KS0[kPt]},
0148 {"KS0_phi", h_cutEff_KS0[kPhi]},
0149 {"KS0_eta", h_cutEff_KS0[kEta]},
0150 {"KS0_rap", h_cutEff_KS0[kRap]},
0151 {"Xi_pT", h_cutEff_Xi[kPt]},
0152 {"Xi_phi", h_cutEff_Xi[kPhi]},
0153 {"Xi_eta", h_cutEff_Xi[kEta]},
0154 {"Xi_rap", h_cutEff_Xi[kRap]},
0155 {"KS0_pT_max", h_cutEff_KS0[kPt_max]},
0156 {"KS0_phi_max", h_cutEff_KS0[kPhi_max]},
0157 {"KS0_eta_max", h_cutEff_KS0[kEta_max]},
0158 {"KS0_rap_max", h_cutEff_KS0[kRap_max]},
0159 {"Xi_pT_max", h_cutEff_Xi[kPt_max]},
0160 {"Xi_phi_max", h_cutEff_Xi[kPhi_max]},
0161 {"Xi_eta_max", h_cutEff_Xi[kEta_max]},
0162 {"Xi_rap_max", h_cutEff_Xi[kRap_max]},
0163 };
0164
0165 std::map<std::string, TH1F*> corrMap_trackAcc =
0166 {
0167 {"KS0_pT", h_trackAcc_KS0[kPt]},
0168 {"KS0_phi", h_trackAcc_KS0[kPhi]},
0169 {"KS0_eta", h_trackAcc_KS0[kEta]},
0170 {"KS0_rap", h_trackAcc_KS0[kRap]},
0171 {"Xi_pT", h_trackAcc_Xi[kPt]},
0172 {"Xi_phi", h_trackAcc_Xi[kPhi]},
0173 {"Xi_eta", h_trackAcc_Xi[kEta]},
0174 {"Xi_rap", h_trackAcc_Xi[kRap]},
0175 {"KS0_pT_min", h_trackAcc_KS0[kPt_min]},
0176 {"KS0_phi_min", h_trackAcc_KS0[kPhi_min]},
0177 {"KS0_eta_min", h_trackAcc_KS0[kEta_min]},
0178 {"KS0_rap_min", h_trackAcc_KS0[kRap_min]},
0179 {"Xi_pT_min", h_trackAcc_Xi[kPt_min]},
0180 {"Xi_phi_min", h_trackAcc_Xi[kPhi_min]},
0181 {"Xi_eta_min", h_trackAcc_Xi[kEta_min]},
0182 {"Xi_rap_min", h_trackAcc_Xi[kRap_min]},
0183 {"KS0_pT_max", h_trackAcc_KS0[kPt_max]},
0184 {"KS0_phi_max", h_trackAcc_KS0[kPhi_max]},
0185 {"KS0_eta_max", h_trackAcc_KS0[kEta_max]},
0186 {"KS0_rap_max", h_trackAcc_KS0[kRap_max]},
0187 {"Xi_pT_max", h_trackAcc_Xi[kPt_max]},
0188 {"Xi_phi_max", h_trackAcc_Xi[kPhi_max]},
0189 {"Xi_eta_max", h_trackAcc_Xi[kEta_max]},
0190 {"Xi_rap_max", h_trackAcc_Xi[kRap_max]},
0191 };
0192
0193 TFile* outputFile = TFile::Open("cutVarianceSystematics.root", "RECREATE");
0194
0195 for (const auto& var: variants)
0196 {
0197 for (const auto& par: particles)
0198 {
0199 for (const auto& dif: diffNames)
0200 {
0201 std::string histName = "h_" + dif + "_base_" + par + var;
0202 std::cout << "Correcting histogram : " << histName << std::endl;
0203 TH1F* h = (TH1F*)inputRawYield->Get(histName.c_str());
0204 h->SetDirectory(0);
0205 int nBins = h->GetNbinsX();
0206
0207 for (int b = 1; b <= nBins; b++)
0208 {
0209 double yield = h->GetBinContent(b);
0210 double statErr = h->GetBinError(b);
0211
0212 double geoAccEff = corrMap_geoAccEff.at(par + "_" + dif)->GetBinContent(b);
0213 double cutEff = corrMap_cutEff.at(par + "_" + dif)->GetBinContent(b);
0214 double trackAcc = corrMap_trackAcc.at(par + "_" + dif)->GetBinContent(b);
0215
0216 double totalCorr;
0217 if (par == "KS0")
0218 {
0219
0220 totalCorr = cutEff * KS02pipi_BR_PDG;
0221 }
0222 else
0223 {
0224
0225 totalCorr = cutEff * Xi2Lambdapi_BR_PDG * Lambda2ppi_BR_PDG;
0226 }
0227 double corrYield = yield / totalCorr;
0228 double corrErr = statErr / totalCorr;
0229 h->SetBinContent(b, corrYield);
0230 h->SetBinError(b, corrErr);
0231 }
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245 h->SetName((histName).c_str());
0246 outputFile->cd();
0247 h->Write();
0248 }
0249 }
0250 }
0251 }