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