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