LightFlavorRatios/swimming_correction/CutEfficiency_systematics.C

0001 #include "CutEfficiency_mjp.C"
0002 #include "PVSVUncertainty.h"
0003
0004 // cuts producing the minimum and maximum cut efficiency allowed by the PV and SV uncertainty
0005 std::pair<CutSettings,CutSettings> uncertainty_modified_cuts(const CutSettings& cuts, const PVSVUncertainty& pvsv)
0006 {
0007   CutSettings new_cuts_min = cuts;
0008   CutSettings new_cuts_max = cuts;
0009
0010   double mean_PV_sigma = pvsv.PV_sigma->GetMean();
0011   double mean_PV_sigma_xy = pvsv.PV_sigma_xy->GetMean();
0012   double mean_SV_sigma = pvsv.SV_sigma->GetMean();
0013   double mean_track_1_PV_DCA_sigma = pvsv.track_1_PV_DCA_sigma->GetMean();
0014   double mean_track_2_PV_DCA_sigma = pvsv.track_2_PV_DCA_sigma->GetMean();
0015   double mean_track_1_PV_DCA_xy_sigma = pvsv.track_1_PV_DCA_xy_sigma->GetMean();
0016   double mean_track_2_PV_DCA_xy_sigma = pvsv.track_2_PV_DCA_xy_sigma->GetMean();
0017
0018   new_cuts_min.track_1_PV_min_DCA_xy = cuts.track_1_PV_min_DCA_xy + mean_PV_sigma_xy + mean_track_1_PV_DCA_xy_sigma;
0019   new_cuts_min.track_2_PV_min_DCA_xy = cuts.track_2_PV_min_DCA_xy + mean_PV_sigma_xy + mean_track_2_PV_DCA_xy_sigma;
0020   new_cuts_min.track_1_track_2_max_DCA = cuts.track_1_track_2_max_DCA - mean_track_1_PV_DCA_sigma - mean_track_2_PV_DCA_sigma;
0021   new_cuts_min.track_1_track_2_max_DCA_xy = cuts.track_1_track_2_max_DCA_xy - mean_track_1_PV_DCA_xy_sigma - mean_track_2_PV_DCA_xy_sigma;
0022   // DIRA cut is more complicated than just the changing of a limit, so we put it in the misc_cutstring
0023   new_cuts_min.min_DIRA = -FLT_MAX;
0024   std::string min_min_DIRA_cut = "cos(acos("+cuts.mother_name+"_DIRA) + atan("+std::to_string(mean_PV_sigma)+"/"+cuts.mother_name+"_decayLength) + atan("+std::to_string(mean_SV_sigma)+"/"+cuts.mother_name+"_decayLength)) >= "+std::to_string(cuts.min_DIRA);
0025   if(!cuts.misc_cutstring.empty()) new_cuts_min.misc_cutstring += " && ";
0026   new_cuts_min.misc_cutstring += min_min_DIRA_cut;
0027
0028   new_cuts_max.track_1_PV_min_DCA_xy = cuts.track_1_PV_min_DCA_xy - mean_PV_sigma_xy - mean_track_1_PV_DCA_xy_sigma;
0029   new_cuts_max.track_2_PV_min_DCA_xy = cuts.track_2_PV_min_DCA_xy - mean_PV_sigma_xy - mean_track_1_PV_DCA_xy_sigma;
0030   new_cuts_max.track_1_track_2_max_DCA = cuts.track_1_track_2_max_DCA + mean_track_1_PV_DCA_sigma + mean_track_2_PV_DCA_sigma;
0031   new_cuts_max.track_1_track_2_max_DCA_xy = cuts.track_1_track_2_max_DCA_xy + mean_track_1_PV_DCA_sigma + mean_track_2_PV_DCA_sigma;
0032   new_cuts_max.min_DIRA = -FLT_MAX;
0033   std::string max_min_DIRA_cut = "cos(acos("+cuts.mother_name+"_DIRA) - atan("+std::to_string(mean_PV_sigma)+"/"+cuts.mother_name+"_decayLength) - atan("+std::to_string(mean_SV_sigma)+"/"+cuts.mother_name+"_decayLength)) >= "+std::to_string(cuts.min_DIRA);
0034   if(!cuts.misc_cutstring.empty()) new_cuts_max.misc_cutstring += " && ";
0035   new_cuts_max.misc_cutstring += max_min_DIRA_cut;
0036
0037   return {new_cuts_min, new_cuts_max};
0038 }
0039
0040 void CutEfficiency_systematics(const std::string& merged_infile = "/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/cutEfficiencyCorrection/merged_cut_efficiency.root",
0041                                const std::string& numerator_name = "Lambda0",
0042                                const int numerator_pdgid = 3122,
0043                                const std::vector<int> numerator_daughters = {-211,2212},
0044                                const std::string& numerator_infile = "/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/closureTestSample/ppi_reco/outputKFParticle_ppi_reco_000001.root",
0045                                const bool numerator_include_opposite = true,
0046                                const CutSettings& numerator_cuts = StandardCuts::MC_Lambda0_cuts,
0047                                const HistogramInfo& numerator_massbins = BinInfo::Lambda0_MC_mass_bins,
0048                                const std::string& denominator_name = "K_S0",
0049                                const int denominator_pdgid = 310,
0050                                const std::vector<int> denominator_daughters = {211, -211},
0051                                const std::string& denominator_infile = "/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/closureTestSample/pipi_reco/outputKFParticle_pipi_reco_000001.root",
0052                                const bool denominator_include_opposite = false,
0053                                const CutSettings& denominator_cuts = StandardCuts::MC_K_S0_cuts,
0054                                const HistogramInfo& denominator_massbins = BinInfo::K_S0_MC_mass_bins,
0055                                const std::string& min_outfile = "test_min.root",
0056                                const std::string& max_outfile = "test_max.root")
0057 {
0058   TFile* infile = TFile::Open(merged_infile.c_str());
0059
0060   if(!infile)
0061   {
0062     std::cout << "ERROR: no nominal cut efficiency merged file detected; cut efficiency systematics determination requires that the nominal cut effficiency calculation already exists!" << std::endl;
0063     return;
0064   }
0065
0066   PVSVUncertainty numerator_pvsv(infile,numerator_name);
0067   PVSVUncertainty denominator_pvsv(infile,denominator_name);
0068
0069   std::pair<CutSettings,CutSettings> numerator_modified_cuts = uncertainty_modified_cuts(numerator_cuts,numerator_pvsv);
0070   std::pair<CutSettings,CutSettings> denominator_modified_cuts = uncertainty_modified_cuts(denominator_cuts,denominator_pvsv);
0071
0072   CutEfficiency_mjp(numerator_name,numerator_pdgid,numerator_daughters,numerator_infile,numerator_include_opposite,
0073                     numerator_modified_cuts.first,numerator_massbins,
0074                     denominator_name,denominator_pdgid,denominator_daughters,denominator_infile,denominator_include_opposite,
0075                     denominator_modified_cuts.first,denominator_massbins,
0076                     min_outfile);
0077
0078   CutEfficiency_mjp(numerator_name,numerator_pdgid,numerator_daughters,numerator_infile,numerator_include_opposite,
0079                     numerator_modified_cuts.second,numerator_massbins,
0080                     denominator_name,denominator_pdgid,denominator_daughters,denominator_infile,denominator_include_opposite,
0081                     denominator_modified_cuts.second,denominator_massbins,
0082                     max_outfile);
0083 }