Back to home page

sPhenix code displayed by LXR

 
 

    


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   //const std::vector<std::string> variants  = {"", "_max"};
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             //totalCorr  = geoAccEff * cutEff * trackAcc * KS02pipi_BR_PDG;
0220             totalCorr  = cutEff * KS02pipi_BR_PDG;
0221           }
0222           else
0223           {
0224             //totalCorr  = geoAccEff * cutEff * trackAcc * Xi2Lambdapi_BR_PDG * Lambda2ppi_BR_PDG;
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         //h->Divide(corrMap_geoAccEff.at(par + "_" + dif));
0234         //h->Divide(corrMap_cutEff.at(par + "_" + dif + var));
0235         //h->Divide(corrMap_trackAcc.at(par + "_" + dif));
0236         //if (par == "KS0")
0237         //{
0238         //  h->Scale(1.0/(KS02pipi_BR_PDG));
0239         //}
0240         //else
0241         //{
0242         //  h->Scale((1.0/(Xi2Lambdapi_BR_PDG*Lambda2ppi_BR_PDG)));
0243         //}
0244         //h->SetName((histName + var).c_str());
0245         h->SetName((histName).c_str());
0246         outputFile->cd();
0247         h->Write();
0248       }
0249     }
0250   }
0251 }