LightFlavorRatios/yield_and_ratios/Lambda_Kshort_ratio_MC.C

0001 #include "../bco_correction/V0DuplicateReader_mod.h"
0002
0003 #include "../util/RooFit_import_TTree.h"
0004
0005 #include "../corrections/EfficiencyCorrection.h"
0006 #include "../corrections/TrivialLambdaFeedDownCorrection.h"
0007 #include "../corrections/LambdaFeedDownCorrection.h"
0008 #include "../corrections/GeoAcceptanceCorrection.h"
0009 #include "../corrections/TrivialEfficiencyCorrection.h"
0010 #include "../corrections/CutEfficiencyCorrection.h"
0011
0012 #include "ResonanceRatio.h"
0013 //#include "calculate_ratios.C"
0014 #include "LambdaModel.h"
0015 #include "KshortModel.h"
0016
0017 void Lambda_Kshort_ratio_MC()
0018 {
0019   TFile* lambda_file = TFile::Open("/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/merged_lambda_MC.root");
0020   TFile* Ks_file = TFile::Open("/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/merged_Kshort_MC.root");
0021
0022   //TFile* lambda_file = TFile::Open("/sphenix/tg/tg01/hf/mjpeters/lambdaKshortMB/lambdaKshort_20260422_DetroitMB_CR_2_mode_pTref_1p4/ppi_reco/merged_lambda.root");
0023   //TFile* Ks_file = TFile::Open("/sphenix/tg/tg01/hf/mjpeters/lambdaKshortMB/lambdaKshort_20260422_DetroitMB_CR_2_mode_pTref_1p4/pipi_reco/merged_kshort.root");
0024
0025   //TFile* lambda_file = TFile::Open("/gpfs/mnt/gpfs02/sphenix/user/cdean/software/analysis/LightFlavorRatios/geometric_acceptance/simulation/outputKFParticle_Lambda2ppi_reco_Usman_patch.root");
0026   //TFile* Ks_file = TFile::Open("/gpfs/mnt/gpfs02/sphenix/user/cdean/software/analysis/LightFlavorRatios/geometric_acceptance/simulation/outputKFParticle_Kshort2pipi_reco_Usman_patch.root");
0027
0028   //TFile* Ks_file = TFile::Open("/sphenix/tg/tg01/hf/cdean/LF_analysis/data_nTuples/output_Kshort_run3pp_looseCuts_20260608.root");
0029   //TFile* lambda_file = TFile::Open("/sphenix/tg/tg01/hf/cdean/LF_analysis/data_nTuples/output_Lambda0_run3pp_looseCuts_20260608.root");
0030
0031   //TFile* Ks_file = TFile::Open("/sphenix/tg/tg01/hf/aopatton/SVLooseJun4/6RunsCombinedKShortSVLoose.root");
0032   //TFile* lambda_file = TFile::Open("/sphenix/tg/tg01/hf/aopatton/SVLooseJun4/6RunsCombinedLambdaSVLoose.root");
0033
0034   //TFile* Ks_file = TFile::Open("/sphenix/tg/tg01/hf/mjpeters/LightFlavorResults/KShort6RunCombined.root");
0035   //TFile* lambda_file = TFile::Open("/sphenix/tg/tg01/hf/mjpeters/LightFlavorResults/Lambda6RunCombined.root");
0036
0037   //TFile* Ks_file = TFile::Open("/sphenix/tg/tg01/hf/mjpeters/LightFlavorResults/Kshort_3runs.root");
0038   //TFile* lambda_file = TFile::Open("/sphenix/tg/tg01/hf/mjpeters/LightFlavorResults/Lambda_3runs.root");
0039
0040   //TTree* Ks_tree = (TTree*)Ks_file->Get("DecayTree");
0041   //TTree* lambda_tree = (TTree*)lambda_file->Get("DecayTree");
0042
0043   TH1F* integrated_lambda_mass = (TH1F*)lambda_file->Get("Lambda0_mass");
0044   TH1F* integrated_kshort_mass = (TH1F*)Ks_file->Get("K_S0_mass");
0045
0046   std::vector<HistogramInfo> diff_variables =
0047   {
0048     BinInfo::final_pt_bins,
0049     BinInfo::final_eta_bins,
0050     BinInfo::final_rapidity_bins,
0051     BinInfo::final_phi_bins,
0052   };
0053
0054   std::map<std::string,HistogramInfo> massbins_map = BinInfo::mass_bins_MC;
0055
0056   std::vector<DifferentialContainer> diff_lambda_data;
0057   std::vector<DifferentialContainer> diff_ks_data;
0058
0059   for(HistogramInfo& hinfo : diff_variables)
0060   {
0061     diff_lambda_data.push_back(DifferentialContainer(lambda_file,"Lambda0",massbins_map,hinfo));
0062     diff_ks_data.push_back(DifferentialContainer(Ks_file,"K_S0",massbins_map,hinfo));
0063   }
0064 /*
0065   RooArgList Ks_args;
0066   RooArgList lambda_args;
0067
0068   RooRealVar m_ks("K_S0_mass","K_S0_mass",0.4,0.6);
0069   RooRealVar m_lambda("Lambda0_mass","Lambda0_mass",1.08,1.15);
0070
0071   Ks_args.add(m_ks);
0072   lambda_args.add(m_lambda);
0073
0074   std::vector<RooRealVar> Ks_diffvars;
0075   std::vector<RooRealVar> lambda_diffvars;
0076
0077   std::vector<RooRealVar> Ks_cutvars;
0078   std::vector<RooRealVar> lambda_cutvars;
0079
0080   std::vector<RooRealVar> Ks_cutvars_int;
0081   std::vector<RooRealVar> lambda_cutvars_int;
0082
0083   for(HistogramInfo& hinfo : diff_variables)
0084   {
0085     std::string Ks_branchname = "K_S0_"+hinfo.name;
0086     std::string lambda_branchname = "Lambda0_"+hinfo.name;
0087     std::cout << Ks_branchname << " " << lambda_branchname << std::endl;
0088     Ks_diffvars.push_back(make_var(Ks_branchname,Ks_branchname,Ks_tree));
0089     lambda_diffvars.push_back(make_var(lambda_branchname,lambda_branchname,lambda_tree));
0090   }
0091
0092   for(int i=0; i<Ks_diffvars.size(); i++)
0093   {
0094     Ks_args.add(Ks_diffvars[i]);
0095     lambda_args.add(lambda_diffvars[i]);
0096   }
0097 */
0098   HistogramInfo Ks_massbins = massbins_map.at("K_S0");
0099   HistogramInfo Lambda_massbins = massbins_map.at("Lambda0");
0100 /*
0101   for(const std::string& cutvar : Ks_massbins.get_cutvars(Ks_tree))
0102   {
0103     if(isIntBranch(Ks_tree->GetBranch(cutvar.c_str())))
0104     {
0105       std::cout << "branch " << cutvar << " is of integral type" << std::endl;
0106       Ks_cutvars_int.push_back(make_var(cutvar,cutvar,Ks_tree));
0107     }
0108     else
0109     {
0110       Ks_cutvars.push_back(make_var(cutvar,cutvar,Ks_tree));
0111     }
0112   }
0113
0114   for(const std::string& cutvar : Lambda_massbins.get_cutvars(lambda_tree))
0115   {
0116     if(isIntBranch(lambda_tree->GetBranch(cutvar.c_str())))
0117     {
0118       std::cout << "branch " << cutvar << " is of integral type" << std::endl;
0119       lambda_cutvars_int.push_back(make_var(cutvar,cutvar,lambda_tree));
0120     }
0121     else
0122     {
0123       lambda_cutvars.push_back(make_var(cutvar,cutvar,lambda_tree));
0124     }
0125   }
0126
0127   for(int i=0; i<Ks_cutvars.size(); i++)
0128   {
0129     Ks_args.add(Ks_cutvars[i]);
0130   }
0131
0132   for(int i=0; i<Ks_cutvars_int.size(); i++)
0133   {
0134     Ks_args.add(Ks_cutvars_int[i]);
0135   }
0136
0137   for(int i=0; i<lambda_cutvars.size(); i++)
0138   {
0139     lambda_args.add(lambda_cutvars[i]);
0140   }
0141
0142   for(int i=0; i<lambda_cutvars_int.size(); i++)
0143   {
0144     lambda_args.add(lambda_cutvars_int[i]);
0145   }
0146
0147   std::string Ks_cuts = Ks_massbins.cut_string;
0148   std::string Lambda_cuts = Lambda_massbins.cut_string;
0149
0150   Ks_args.Print();
0151   lambda_args.Print();
0152
0153   RooDataSet* Ks_ds = new RooDataSet("K_S0","K_S0",Ks_args,RooFit::Import(*Ks_tree));
0154   RooDataSet* lambda_ds = new RooDataSet("Lambda0","Lambda0",lambda_args,RooFit::Import(*lambda_tree));
0155
0156   V0DuplicateReader ks_reader(Ks_tree, V0DuplicateReader::ParticleType::K0s);
0157   V0DuplicateReader lambda_reader(lambda_tree, V0DuplicateReader::ParticleType::Lambda);
0158
0159   ks_reader.enableDeltaBCOCut(0, 350);
0160   lambda_reader.enableDeltaBCOCut(0, 350);
0161
0162   for (Long64_t i = 0; i < ks_reader.entries(); ++i)
0163   {
0164     if(i % 10000 == 0) std::cout << "processing BCO for Kshorts entry " << i << " / " << ks_reader.entries() << std::endl;
0165     ks_reader.loadEntry(i);
0166
0167     if (!ks_reader.passesDeltaBCOCut()) continue;
0168     if (!ks_reader.isCurrentEntryUnique()) continue;
0169
0170     m_ks.setVal(ks_reader.get<float>("K_S0_mass"));
0171
0172     for(size_t idiff = 0; idiff < diff_variables.size(); idiff++)
0173     {
0174       Ks_diffvars[idiff].setVal(ks_reader.get<float>(Ks_diffvars[idiff].GetName()));
0175     }
0176
0177     for(size_t icut = 0; icut < Ks_cutvars.size(); icut++)
0178     {
0179       Ks_cutvars[icut].setVal(ks_reader.get<float>(Ks_cutvars[icut].GetName()));
0180     }
0181
0182     for(size_t icut_int = 0; icut_int < Ks_cutvars_int.size(); icut_int++)
0183     {
0184       Ks_cutvars_int[icut_int].setVal(ks_reader.get<int>(Ks_cutvars_int[icut_int].GetName()));
0185     }
0186
0187     Ks_ds->add(Ks_args);
0188   }
0189
0190   for (Long64_t i = 0; i < lambda_reader.entries(); ++i)
0191   {
0192     if(i % 10000 == 0) std::cout << "processing BCO for lambda entry " << i << " / " << lambda_reader.entries() << std::endl;
0193     lambda_reader.loadEntry(i);
0194
0195     if (!lambda_reader.passesDeltaBCOCut()) continue;
0196     if (!lambda_reader.isCurrentEntryUnique()) continue;
0197
0198     m_lambda.setVal(lambda_reader.get<float>("Lambda0_mass"));
0199
0200     for(size_t idiff = 0; idiff < diff_variables.size(); idiff++)
0201     {
0202       lambda_diffvars[idiff].setVal(lambda_reader.get<float>(lambda_diffvars[idiff].GetName()));
0203     }
0204
0205     for(size_t icut = 0; icut < lambda_cutvars.size(); icut++)
0206     {
0207       lambda_cutvars[icut].setVal(lambda_reader.get<float>(lambda_cutvars[icut].GetName()));
0208     }
0209
0210     for(size_t icut_int = 0; icut_int < lambda_cutvars_int.size(); icut_int++)
0211     {
0212       lambda_cutvars_int[icut_int].setVal(lambda_reader.get<int>(lambda_cutvars_int[icut_int].GetName()));
0213     }
0214
0215     lambda_ds->add(lambda_args);
0216   }
0217
0218   RooDataSet* Ks_ds_withcuts = (RooDataSet*)Ks_ds->reduce(Ks_args,Ks_cuts.c_str());
0219   RooDataSet* lambda_ds_withcuts = (RooDataSet*)lambda_ds->reduce(lambda_args,Lambda_cuts.c_str());
0220
0221   std::cout << "Ks_cuts " << Ks_cuts << std::endl;
0222   std::cout << "Lamdba_cuts " << Lambda_cuts << std::endl;
0223 */
0224
0225   KshortModel kshort_model(Ks_massbins);
0226   LambdaModel lambda_model(Lambda_massbins);
0227
0228   std::string fd_filename = "/sphenix/tg/tg01/hf/hjheng/HF-analysis/simulation/Pythia_ppMinBias/cascade_feeddown/Cascade_feeddown_fraction.root";
0229   std::vector<std::vector<std::shared_ptr<CorrectionHistogram1D>>> corrections(diff_variables.size());
0230   // pT
0231   corrections[0].push_back(std::make_shared<TrivialLambdaFeedDownCorrection>(fd_filename,"h_feeddown_frac_xi_all"));
0232   corrections[0].push_back(std::make_shared<TrivialEfficiencyCorrection>(""));
0233 //  corrections[0].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/u/cdean/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots/Lambda0_to_KS0_geometric_acceptance_ratio_pT.root","Lambda0_inGeo_pT"));
0234 //  corrections[0].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/tg/tg01/hf/gregoryottino/lightFlavorPpg16/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots_systemtics/Lambda0_to_KS0_geometric_acceptance_ratio_pT.root","Lambda0_inGeo_pT"));
0235 //  corrections[0].push_back(std::make_shared<CutEfficiencyCorrection>("../swimming_correction/LamdbaKsCutEfficiency_200MeV_hists.root","hEffRatio_pT"));
0236   corrections[0].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/geometricAcceptanceCorrection/corrections/geo_acceptance_inclusive.root","Lambda0_over_K_S0_geo_acceptance_correction_vspT"));
0237   corrections[0].push_back(std::make_shared<CutEfficiencyCorrection>("/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/cutEfficiencyCorrection/cut_efficiency_correction.root","Lambda0_over_K_S0_cuteff_correction_vspT"));
0238   // eta
0239   corrections[1].push_back(std::make_shared<TrivialLambdaFeedDownCorrection>(fd_filename,"h_feeddown_frac_xi_eta_all"));
0240   corrections[1].push_back(std::make_shared<TrivialEfficiencyCorrection>(""));
0241 //  corrections[1].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/u/cdean/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots/Lambda0_to_KS0_geometric_acceptance_ratio_eta.root","Lambda0_inGeo_#eta"));
0242 //  corrections[1].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/tg/tg01/hf/gregoryottino/lightFlavorPpg16/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots_systemtics/Lambda0_to_KS0_geometric_acceptance_ratio_eta.root","Lambda0_inGeo_#eta"));
0243 //  corrections[1].push_back(std::make_shared<CutEfficiencyCorrection>("../swimming_correction/LamdbaKsCutEfficiency_200MeV_hists.root","hEffRatio_eta"));
0244   corrections[1].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/geometricAcceptanceCorrection/corrections/geo_acceptance_inclusive.root","Lambda0_over_K_S0_geo_acceptance_correction_vspseudorapidity"));
0245   corrections[1].push_back(std::make_shared<CutEfficiencyCorrection>("/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/cutEfficiencyCorrection/cut_efficiency_correction.root","Lambda0_over_K_S0_cuteff_correction_vspseudorapidity"));
0246
0247
0248   // rapidity
0249   corrections[2].push_back(std::make_shared<TrivialLambdaFeedDownCorrection>(fd_filename,"h_feeddown_frac_xi_rapidity_all"));
0250   corrections[2].push_back(std::make_shared<TrivialEfficiencyCorrection>(""));
0251 //  corrections[2].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/u/cdean/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots/Lambda0_to_KS0_geometric_acceptance_ratio_rap.root","Lambda0_inGeo_y"));
0252 //  corrections[2].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/tg/tg01/hf/gregoryottino/lightFlavorPpg16/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots_systemtics/Lambda0_to_KS0_geometric_acceptance_ratio_rap.root","Lambda0_inGeo_y"));
0253 //  corrections[2].push_back(std::make_shared<CutEfficiencyCorrection>("../swimming_correction/LamdbaKsCutEfficiency_200MeV_hists.root","hEffRatio_y"));
0254   corrections[2].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/geometricAcceptanceCorrection/corrections/geo_acceptance_inclusive.root","Lambda0_over_K_S0_geo_acceptance_correction_vsrapidity"));
0255   corrections[2].push_back(std::make_shared<CutEfficiencyCorrection>("/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/cutEfficiencyCorrection/cut_efficiency_correction.root","Lambda0_over_K_S0_cuteff_correction_vsrapidity"));
0256
0257
0258   // phi
0259   corrections[3].push_back(std::make_shared<TrivialLambdaFeedDownCorrection>(fd_filename,"h_feeddown_frac_xi_phi_all"));
0260   corrections[3].push_back(std::make_shared<TrivialEfficiencyCorrection>(""));
0261 //  corrections[3].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/u/cdean/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots/Lambda0_to_KS0_geometric_acceptance_ratio_phi.root","Lambda0_inGeo_#phi"));
0262 //  corrections[3].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/tg/tg01/hf/gregoryottino/lightFlavorPpg16/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots_systemtics/Lambda0_to_KS0_geometric_acceptance_ratio_phi.root","Lambda0_inGeo_#phi"));
0263 //  corrections[3].push_back(std::make_shared<CutEfficiencyCorrection>("../swimming_correction/LamdbaKsCutEfficiency_200MeV_hists.root","hEffRatio_phi"));
0264   corrections[3].push_back(std::make_shared<GeoAcceptanceCorrection>("/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/geometricAcceptanceCorrection/corrections/geo_acceptance_inclusive.root","Lambda0_over_K_S0_geo_acceptance_correction_vsphi"));
0265   corrections[3].push_back(std::make_shared<CutEfficiencyCorrection>("/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/cutEfficiencyCorrection/cut_efficiency_correction.root","Lambda0_over_K_S0_cuteff_correction_vsphi"));
0266
0267   TFile* fout = new TFile("fits_MC.root","RECREATE");
0268
0269   ResonanceRatio analyzer(lambda_model,kshort_model,massbins_map,
0270                           fout,"lambdaKsratio","(#Lambda^{0}+#bar{#Lambda^{0}})/2K_{S}^{0} ratio",1./2.,false,
0271                           diff_variables,corrections);
0272
0273   analyzer.calculate_ratios_binned(integrated_lambda_mass,diff_lambda_data,integrated_kshort_mass,diff_ks_data);
0274 }