LightFlavorRatios/swimming_correction/CutEfficiency_mjp.C

0001 #include "../util/binning.h"
0002 #include "../yield_and_ratios/LambdaModel.h"
0003 #include "../yield_and_ratios/KshortModel.h"
0004
0005 TH1F* n_truth_mass = new TH1F("n_truth_mass","num. truth mass",1000,1.05,1.25);
0006 TH1F* n_reco_mass = new TH1F("n_reco_mass","num. reco mass",1000,1.05,1.25);
0007 TH1F* d_truth_mass = new TH1F("d_truth_mass","denom. truth mass",1000,0.35,0.65);
0008 TH1F* d_reco_mass = new TH1F("d_reco_mass","denom. reco mass",1000,0.35,0.65);
0009
0010 // must have some permutation of the corect flavor of daughters
0011 std::string correct_truth_daughters_cut(std::vector<int> daughters, bool include_opposite)
0012 {
0013   std::vector<std::string> flavor_permutations;
0014   std::vector<std::string> opposite_flavor_permutations;
0015
0016   std::vector<int> indices(daughters.size());
0017   std::iota(indices.begin(),indices.end(),0);
0018
0019   do
0020   {
0021     std::string cut = "(";
0022     std::string opposite_parity_cut = "(";
0023     for(int i=0; i<indices.size(); i++)
0024     {
0025       cut += "track_"+std::to_string(i+1)+"_true_ID == "+std::to_string(daughters[indices[i]]);
0026       opposite_parity_cut += "track_"+std::to_string(i+1)+"_true_ID == "+std::to_string(-1*daughters[indices[i]]);
0027       if(i<indices.size()-1)
0028       {
0029         cut += " && ";
0030         opposite_parity_cut += " && ";
0031       }
0032     }
0033     cut += ")";
0034     opposite_parity_cut += ")";
0035     flavor_permutations.push_back(cut);
0036     opposite_flavor_permutations.push_back(opposite_parity_cut);
0037   } while(std::next_permutation(indices.begin(),indices.end()));
0038
0039   std::string total_cut = "";
0040   for(int i=0;i<flavor_permutations.size();i++)
0041   {
0042     total_cut += flavor_permutations[i];
0043     if(include_opposite)
0044     {
0045       total_cut += " || "+opposite_flavor_permutations[i];
0046     }
0047     if(i<flavor_permutations.size()-1)
0048     {
0049       total_cut += " || ";
0050     }
0051   }
0052   return total_cut;
0053 }
0054
0055 std::string correct_truth_parent_cut(int pdgid, int n_daughters, bool include_opposite)
0056 {
0057   std::vector<std::string> daughter_cuts;
0058
0059   for(int i=1;i<=n_daughters;i++)
0060   {
0061     std::string cut = "Sum$(";
0062     if(include_opposite)
0063     {
0064       cut += "abs(";
0065     }
0066     cut += "track_"+std::to_string(i)+"_true_track_history_PDG_ID";
0067     if(include_opposite)
0068     {
0069       cut += ")";
0070     }
0071     cut += " == "+std::to_string(pdgid)+") > 0";
0072     daughter_cuts.push_back(cut);
0073   }
0074
0075   std::string total_cut = "";
0076   for(int i=0;i<daughter_cuts.size();i++)
0077   {
0078     total_cut += daughter_cuts[i];
0079     if(i<daughter_cuts.size()-1)
0080     {
0081       total_cut += " && ";
0082     }
0083   }
0084   return total_cut;
0085 }
0086
0087 std::string correct_daughter_PDGID_assignment_cut(int n_daughters)
0088 {
0089   std::vector<std::string> daughter_cuts;
0090   for(int i=1;i<=n_daughters;i++)
0091   {
0092     std::string cut = "track_"+std::to_string(i)+"_PDG_ID == track_"+std::to_string(i)+"_true_ID";
0093     daughter_cuts.push_back(cut);
0094   }
0095
0096   std::string total_cut = "";
0097   for(int i=0;i<n_daughters;i++)
0098   {
0099     total_cut += daughter_cuts[i];
0100     if(i<n_daughters-1)
0101     {
0102       total_cut += " && ";
0103     }
0104   }
0105   return total_cut;
0106 }
0107
0108 // "primary parent" means that no other hadron/lepton is encountered in the post-hadronization PDGID history
0109 // in other words, the mother is the first hadron in the post-hadronization PDGID history
0110 std::pair<int,int> posthadronization_history_info(std::vector<int>* pdgid_history)
0111 {
0112   int n_entries = 0;
0113   int primary_parent = 0;
0114   for(int i=0;i<pdgid_history->size();i++)
0115   {
0116     const int pdgid = pdgid_history->at(i);
0117     // "nontrivial" here means "meson, baryon, or lepton"
0118     const std::string particle_class = TDatabasePDG::Instance()->GetParticle(pdgid)->ParticleClass();
0119     if(particle_class == "Meson" || particle_class == "Baryon" || particle_class == "Lepton")
0120     {
0121       n_entries++;
0122       primary_parent = pdgid;
0123     }
0124     else
0125     {
0126       break;
0127     }
0128   }
0129   return {n_entries,primary_parent};
0130 }
0131
0132 /*
0133 std::string primary_truth_parent_cut(const int pdgid, const int n_daughters, const bool include_opposite)
0134 {
0135   std::string cut;
0136   for(int i=1;i<=n_daughters;i++)
0137   {
0138     cut += "track_"+std::to_string(i)+"_true_track_history_PDGID[0] == "+std:to_string(pdgid)+" && ";
0139     cut += "[&](std::vector<int>
0140     if(include_opposite)
0141     {
0142       cut += " || has_primary_parent("+std::to_string(-1*pdgid)+",track_"+std::to_string(i)+"_true_track_history_PDGID)";
0143     }
0144     cut += ")";
0145     if(i<n_daughters)
0146     {
0147       cut += " && ";
0148     }
0149   }
0150   return cut;
0151 }
0152 */
0153
0154 std::string full_truth_cut(const std::string& name, const int pdgid, const std::vector<int> daughters, const bool include_opposite)
0155 {
0156   std::string correct_truth_daughters = correct_truth_daughters_cut(daughters,include_opposite);
0157   std::string correct_truth_parent = correct_truth_parent_cut(pdgid,daughters.size(),include_opposite);
0158   //std::string primary_truth_parent = primary_truth_parent_cut(pdgid,daughters.size(),include_opposite);
0159   std::string rapidity_cut = "fabs("+name+"_rapidity)<1.";
0160   std::string geoaccept_cut = "track_1_MVTX_nHits>0 && track_2_MVTX_nHits>0 && fabs(primary_vertex_z)<10.";
0161   std::string correct_daughter_pdgid = correct_daughter_PDGID_assignment_cut(daughters.size());
0162
0163   return "("+correct_truth_daughters+") && ("+correct_truth_parent+") && ("+geoaccept_cut+") && ("+correct_daughter_pdgid+") && ("+rapidity_cut+") && ("+BinInfo::fiducial_cuts(name,{BinInfo::final_pt_bins,BinInfo::final_eta_bins,BinInfo::final_phi_bins,BinInfo::final_rapidity_bins})+")";
0164 }
0165
0166 std::string full_reco_cut(const std::string& name, const int pdgid, const std::pair<float,float>& mass_window, const std::vector<int> daughters, const bool include_opposite, const std::map<std::string,HistogramInfo>& massbins_map)
0167 {
0168   std::string truth_cut = full_truth_cut(name,pdgid,daughters,include_opposite);
0169   const HistogramInfo massbins = massbins_map.at(name);
0170   std::string reco_cuts = massbins.cut_string;
0171   // for sideband-subtraction yield extraction, must also include the restriction of the yield to the signal window
0172   // std::string mass_window_cuts = name+"_mass>"+std::to_string(mass_window.first)+" && "+name+"_mass<"+std::to_string(mass_window.second);
0173   std::string correct_daughter_pdgid = correct_daughter_PDGID_assignment_cut(daughters.size());
0174   return truth_cut+" && "+reco_cuts+" && "+correct_daughter_pdgid;
0175 }
0176
0177 void post_draw_process(TTree* t, TH1F* h, const HistogramInfo& var, const std::string& name, TEntryList* elist, const int pdgid, const int n_daughters)
0178 {
0179   t->ResetBranchAddresses();
0180   t->SetBranchStatus("*",true);
0181   std::vector<std::vector<int>*> daughter_pdgid_histories;
0182   daughter_pdgid_histories.resize(n_daughters);
0183   float var_branch;
0184   float mass;
0185   for(int i=0;i<n_daughters;i++)
0186   {
0187     daughter_pdgid_histories[i] = nullptr;
0188     t->SetBranchAddress(("track_"+std::to_string(i+1)+"_true_track_history_PDG_ID").c_str(),&daughter_pdgid_histories[i]);
0189   }
0190   t->SetBranchAddress((name+"_"+var.name).c_str(),&var_branch);
0191
0192   t->SetBranchAddress((name+"_mass").c_str(),&mass);
0193
0194   t->SetEntryList(elist);
0195   for(size_t e=0;e<elist->GetN();e++)
0196   {
0197     //std::cout << elist->GetEntry(e) << std::endl;
0198     t->GetEntry(elist->GetEntry(e));
0199
0200     // fill reco mass histograms after truth cuts and after reco cuts
0201     if(std::string(h->GetName())=="Lambda0_truth_vspT") n_truth_mass->Fill(mass);
0202     if(std::string(h->GetName())=="Lambda0_reco_vspT") n_reco_mass->Fill(mass);
0203     if(std::string(h->GetName())=="K_S0_truth_vspT") d_truth_mass->Fill(mass);
0204     if(std::string(h->GetName())=="K_S0_reco_vspT") d_reco_mass->Fill(mass);
0205     // special case for Kshorts
0206     if(pdgid==310)
0207     {
0208       std::pair<int,int> daughter1_historyinfo = posthadronization_history_info(daughter_pdgid_histories[0]);
0209       std::pair<int,int> daughter2_historyinfo = posthadronization_history_info(daughter_pdgid_histories[1]);
0210       bool primary_Kshort = (daughter_pdgid_histories[0]->at(0) == 310 && daughter_pdgid_histories[1]->at(0) == 310 && daughter1_historyinfo.first == 1 && daughter2_historyinfo.first == 1);
0211       bool primary_K0 = daughter_pdgid_histories[0]->at(0) == 310 && daughter_pdgid_histories[1]->at(0) == 310 && daughter1_historyinfo.first == 2 && daughter2_historyinfo.first == 2 && ((daughter1_historyinfo.second == 311 && daughter2_historyinfo.second == 311) || (daughter1_historyinfo.second == -311 && daughter2_historyinfo.second == -311));
0212       if(primary_Kshort || primary_K0)
0213       {
0214         h->Fill(var_branch);
0215       }
0216     }
0217     else
0218     {
0219       bool all_daughters_primary_parents = std::all_of(daughter_pdgid_histories.begin(),daughter_pdgid_histories.end(),
0220              [&](std::vector<int>* history){return (history->at(0)==pdgid && posthadronization_history_info(history).first==1);});
0221       if(all_daughters_primary_parents)
0222       {
0223         h->Fill(var_branch);
0224       }
0225     }
0226   }
0227 }
0228
0229 void CutEfficiency_mjp(const std::string& numerator_name = "Lambda0", const int numerator_pdgid = 3122, const std::vector<int> numerator_daughters = {-211,2212},
0230                        const std::string& numerator_infile = "/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/closureTestSample/ppi_reco/outputKFParticle_ppi_reco_009501.root", const bool numerator_include_opposite = true,
0231                        const std::string& denominator_name = "K_S0", const int denominator_pdgid = 310, const std::vector<int> denominator_daughters = {211, -211},
0232                        const std::string& denominator_infile = "/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/closureTestSample/pipi_reco/outputKFParticle_pipi_reco_009501.root", const bool denominator_include_opposite = false,
0233                        const std::string& outfile = "test.root", const std::map<std::string,HistogramInfo>& massbins_map = BinInfo::mass_bins_MC)
0234 {
0235
0236   bool verbose = true;
0237
0238   TFile* f_numerator = TFile::Open(numerator_infile.c_str());
0239   TFile* f_denominator = TFile::Open(denominator_infile.c_str());
0240
0241   TTree* t_numerator = (TTree*)f_numerator->Get("DecayTree");
0242   TTree* t_denominator = (TTree*)f_denominator->Get("DecayTree");
0243
0244   std::vector<HistogramInfo> variables =
0245   {
0246     BinInfo::final_pt_bins,
0247     BinInfo::final_eta_bins,
0248     BinInfo::final_phi_bins,
0249     BinInfo::final_rapidity_bins
0250   };
0251
0252   if(verbose)
0253   {
0254     std::string numerator_correct_truth_daughters = correct_truth_daughters_cut(numerator_daughters,numerator_include_opposite);
0255     std::string denominator_correct_truth_daughters = correct_truth_daughters_cut(denominator_daughters,denominator_include_opposite);
0256
0257     std::cout << "correct truth daughter cuts:" << std::endl
0258      << "numerator:" << std::endl
0259      << numerator_correct_truth_daughters << std::endl
0260      << "denominator:" << std::endl
0261      << denominator_correct_truth_daughters << std::endl;
0262
0263     std::string numerator_correct_truth_parent = correct_truth_parent_cut(numerator_pdgid,numerator_daughters.size(),numerator_include_opposite);
0264     std::string denominator_correct_truth_parent = correct_truth_parent_cut(denominator_pdgid,denominator_daughters.size(),denominator_include_opposite);
0265
0266     std::cout << "correct truth parent cuts:" << std::endl
0267      << "numerator:" << std::endl
0268      << numerator_correct_truth_parent << std::endl
0269      << "denominator:" << std::endl
0270      << denominator_correct_truth_parent << std::endl;
0271 /*
0272     std::string numerator_primary_truth_parent = primary_truth_parent_cut(numerator_pdgid,numerator_daughters.size(),numerator_include_opposite);
0273     std::string denominator_primary_truth_parent = primary_truth_parent_cut(denominator_pdgid,denominator_daughters.size(),denominator_include_opposite);
0274
0275     std::cout << "primary truth parent cuts:" << std::endl
0276      << "numerator:" << std::endl
0277      << numerator_primary_truth_parent << std::endl
0278      << "denominator:" << std::endl
0279      << denominator_primary_truth_parent << std::endl;
0280 */
0281   }
0282
0283   LambdaModel lambdamodel(massbins_map.at("Lambda0"));
0284   KshortModel kshortmodel(massbins_map.at("K_S0"));
0285
0286   std::string numerator_truth_cut = full_truth_cut(numerator_name,numerator_pdgid,numerator_daughters,numerator_include_opposite);
0287   std::string numerator_reco_cut = full_reco_cut(numerator_name,numerator_pdgid,{lambdamodel.left_sideband.second,lambdamodel.right_sideband.first},numerator_daughters,numerator_include_opposite,massbins_map);
0288
0289   std::string denominator_truth_cut = full_truth_cut(denominator_name,denominator_pdgid,denominator_daughters,denominator_include_opposite);
0290   std::string denominator_reco_cut = full_reco_cut(denominator_name,denominator_pdgid,{kshortmodel.left_sideband.second,kshortmodel.right_sideband.first},denominator_daughters,denominator_include_opposite,massbins_map);
0291
0292   if(verbose)
0293   {
0294     std::cout << "full cuts:" << std::endl
0295      << numerator_name << " truth:" << std::endl
0296      << numerator_truth_cut << std::endl
0297      << numerator_name << " reco:" << std::endl
0298      << numerator_reco_cut << std::endl
0299      << denominator_name << " truth:" << std::endl
0300      << denominator_truth_cut << std::endl
0301      << denominator_name << " reco:" << std::endl
0302      << denominator_reco_cut << std::endl;
0303   }
0304
0305   TFile* fout = new TFile(outfile.c_str(),"RECREATE");
0306 /*
0307   TEntryList* numerator_truth_elist = new TEntryList("numerator_truth_elist","");
0308   TEntryList* numerator_reco_elist = new TEntryList("numerator_reco_elist","");
0309   TEntryList* denominator_truth_elist = new TEntryList("denominator_truth_elist","");
0310   TEntryList* denominator_reco_elist = new TEntryList("denominator_reco_elist","");
0311 */
0312   t_numerator->Draw(">>numerator_truth_elist",numerator_truth_cut.c_str(),"entrylist goff");
0313   t_numerator->Draw(">>numerator_reco_elist",numerator_reco_cut.c_str(),"entrylist goff");
0314   t_denominator->Draw(">>denominator_truth_elist",denominator_truth_cut.c_str(),"entrylist goff");
0315   t_denominator->Draw(">>denominator_reco_elist",denominator_reco_cut.c_str(),"entrylist goff");
0316
0317   TEntryList* numerator_truth_elist = (TEntryList*)gDirectory->Get("numerator_truth_elist");
0318   TEntryList* numerator_reco_elist = (TEntryList*)gDirectory->Get("numerator_reco_elist");
0319   TEntryList* denominator_truth_elist = (TEntryList*)gDirectory->Get("denominator_truth_elist");
0320   TEntryList* denominator_reco_elist = (TEntryList*)gDirectory->Get("denominator_reco_elist");
0321
0322   std::cout << "pre-primary selection: " << numerator_truth_elist->GetN() << " " << numerator_reco_elist->GetN() << " " << denominator_truth_elist->GetN() << " " << denominator_reco_elist->GetN() << std::endl;
0323
0324   for(HistogramInfo& var : variables)
0325   {
0326     TH1F* numerator_truth_h = makeHistogram((numerator_name+"_truth").c_str(),(numerator_name+" truth candidates").c_str(),var);
0327     TH1F* numerator_reco_h = makeHistogram((numerator_name+"_reco").c_str(),(numerator_name+" reco candidates").c_str(),var);
0328     TH1F* denominator_truth_h = makeHistogram((denominator_name+"_truth").c_str(),(denominator_name+" truth candidates").c_str(),var);
0329     TH1F* denominator_reco_h = makeHistogram((denominator_name+"_reco").c_str(),(denominator_name+" reco candidates").c_str(),var);
0330
0331     numerator_truth_h->Sumw2();
0332     numerator_reco_h->Sumw2();
0333     denominator_truth_h->Sumw2();
0334     denominator_reco_h->Sumw2();
0335
0336     post_draw_process(t_numerator,numerator_truth_h,var,numerator_name,numerator_truth_elist,numerator_pdgid,numerator_daughters.size());
0337     post_draw_process(t_numerator,numerator_reco_h,var,numerator_name,numerator_reco_elist,numerator_pdgid,numerator_daughters.size());
0338     post_draw_process(t_denominator,denominator_truth_h,var,denominator_name,denominator_truth_elist,denominator_pdgid,denominator_daughters.size());
0339     post_draw_process(t_denominator,denominator_reco_h,var,denominator_name,denominator_reco_elist,denominator_pdgid,denominator_daughters.size());
0340
0341     if(verbose)
0342     {
0343       std::cout << var.name << ":" << std::endl;
0344       std::cout << "numerator: passed " << numerator_reco_h->GetEntries() << " / " << numerator_truth_h->GetEntries() << std::endl;
0345       std::cout << "denominator: passed " << denominator_reco_h->GetEntries() << " / " << denominator_truth_h->GetEntries() << std::endl;
0346     }
0347
0348     numerator_truth_h->Write();
0349     numerator_reco_h->Write();
0350     denominator_truth_h->Write();
0351     denominator_reco_h->Write();
0352   }
0353
0354   n_truth_mass->Write();
0355   n_reco_mass->Write();
0356   d_truth_mass->Write();
0357   d_reco_mass->Write();
0358 }