43 #include "filval/filval.hpp" 44 #include "filval_root/filval_root.hpp" 50 #define W_MASS 80.385 // GeV/c^2 51 #define Z_MASS 91.188 // GeV/c^2 52 #define T_MASS 172.44 // GeV/c^2 58 void enable_branches(MiniTreeDataSet& mt){
59 mt.fChain->SetBranchStatus(
"*",
false);
61 mt.track_branch<
int>(
"nLepGood");
62 mt.track_branch_vec<
int >(
"nLepGood",
"LepGood_pdgId");
63 mt.track_branch_vec<
float>(
"nLepGood",
"LepGood_pt");
64 mt.track_branch_vec<
float>(
"nLepGood",
"LepGood_eta");
65 mt.track_branch_vec<
float>(
"nLepGood",
"LepGood_phi");
66 mt.track_branch_vec<
float>(
"nLepGood",
"LepGood_mass");
67 mt.track_branch_vec<
int >(
"nLepGood",
"LepGood_charge");
68 mt.track_branch_vec<
int >(
"nLepGood",
"LepGood_mcMatchId");
69 mt.track_branch_vec<
int >(
"nLepGood",
"LepGood_mcMatchPdgId");
70 mt.track_branch_vec<
int >(
"nLepGood",
"LepGood_mcMatchAny");
71 mt.track_branch_vec<
int >(
"nLepGood",
"LepGood_mcMatchTau");
72 mt.track_branch_vec<
int >(
"nLepGood",
"LepGood_mcPt");
74 mt.track_branch<
int>(
"nJet");
75 mt.track_branch_vec<
float>(
"nJet",
"Jet_pt");
76 mt.track_branch_vec<
float>(
"nJet",
"Jet_eta");
77 mt.track_branch_vec<
float>(
"nJet",
"Jet_phi");
78 mt.track_branch_vec<
float>(
"nJet",
"Jet_mass");
79 mt.track_branch_vec<
float>(
"nJet",
"Jet_btagCMVA");
80 mt.track_branch_vec<
int >(
"nJet",
"Jet_mcMatchFlav");
81 mt.track_branch_vec<
int >(
"nJet",
"Jet_mcMatchId");
82 mt.track_branch_vec<
int >(
"nJet",
"Jet_mcFlavour");
85 mt.track_branch<
int>(
"nGenPart");
86 mt.track_branch_vec<
int >(
"nGenPart",
"GenPart_pdgId");
87 mt.track_branch_vec<
int >(
"nGenPart",
"GenPart_motherIndex");
88 mt.track_branch_vec<
int >(
"nGenPart",
"GenPart_motherId");
89 mt.track_branch_vec<
float>(
"nGenPart",
"GenPart_pt");
90 mt.track_branch_vec<
float>(
"nGenPart",
"GenPart_eta");
91 mt.track_branch_vec<
float>(
"nGenPart",
"GenPart_phi");
92 mt.track_branch_vec<
float>(
"nGenPart",
"GenPart_mass");
93 mt.track_branch_vec<
int >(
"nGenPart",
"GenPart_status");
95 mt.track_branch<
int>(
"nBJetLoose40");
96 mt.track_branch<
int>(
"nBJetMedium40");
97 mt.track_branch<
int>(
"nBJetTight40");
100 mt.track_branch<
int>(
"nVert");
102 mt.track_branch<
int >(
"run" );
103 mt.track_branch<
int >(
"lumi");
104 mt.track_branch<
int >(
"evt" );
105 mt.track_branch<
float>(
"xsec");
114 Jet(
const TLorentzVector& v,
int idx,
int pdgid,
float b_cmva)
115 :v(v),idx(idx),pdgid(pdgid),b_cmva(b_cmva) { }
117 static Jet reco(
const TLorentzVector& v,
int idx,
float b_cmva){
118 return Jet(v, idx, 0, b_cmva);
121 static Jet mc(
const TLorentzVector& v,
int idx,
int pdgid){
122 return Jet(v, idx, pdgid, 0);
126 void declare_values(MiniTreeDataSet& mt){
129 lorentz_vectors(
"LepGood_pt",
"LepGood_eta",
"LepGood_phi",
"LepGood_mass",
"LepGood_4v");
130 lorentz_vectors(
"GenPart_pt",
"GenPart_eta",
"GenPart_phi",
"GenPart_mass",
"GenPart_4v");
131 lorentz_vectors(
"Jet_pt",
"Jet_eta",
"Jet_phi",
"Jet_mass",
"Jet_4v" );
133 energies(
"GenPart_4v",
"GenPart_energy");
137 auto& b_mva_filter = GenFunction::register_function<bool(Jet)>(
"b_mva_filter",
138 FUNC(([cut=0.0](
const Jet& j){
139 return j.b_cmva > cut;
141 auto& b_pdgid_filter = GenFunction::register_function<bool(Jet)>(
"b_pdgid_filter",
142 FUNC(([](
const Jet& j){
143 return j.pdgid == 5 || j.pdgid==-5;
145 auto& w_mass_filter = GenFunction::register_function<bool(Jet, Jet)>(
"w_mass_filter",
146 FUNC(([win_l=W_MASS-10, win_h=W_MASS+10](
const Jet& j1,
const Jet& j2){
147 float inv_mass = (j1.v + j2.v).M();
148 return inv_mass > win_l && inv_mass < win_h;
150 auto& dup_filter = GenFunction::register_function<bool(std::tuple<Jet,Jet>,Jet)>(
"dup_filter",
151 FUNC(([](
const std::tuple<Jet,Jet>& w,
const Jet& b){
153 int j1 = std::get<0>(w).idx;
154 int j2 = std::get<1>(w).idx;
155 return (j0 != j1) && (j0 != j2) && (j1 != j2);
157 auto& qg_id_filter = GenFunction::register_function<bool(Jet, Jet)>(
"qg_id_filter",
158 FUNC(([](
const Jet& j1,
const Jet& j2){
160 int id1 = abs(j1.pdgid);
161 int id2 = abs(j2.pdgid);
162 return ((id1 >=1 && id1 <= 5) || id1 == 21) &&
163 ((id2 >=1 && id2 <= 5) || id2 == 21);
167 auto jets = apply(GenFunction::register_function<std::vector<Jet>(std::vector<TLorentzVector>,std::vector<float>)>(
"build_reco_jets",
168 FUNC(([](
const std::vector<TLorentzVector>& vs,
const std::vector<float>& b_cmvas){
169 std::vector<Jet> jets;
170 for(
int i=0; i<vs.size(); i++){
171 jets.push_back(Jet::reco(vs[i],i, b_cmvas[i]));
174 }))), fv::tuple(lookup<std::vector<TLorentzVector>>(
"Jet_4v"), lookup<std::vector<float>>(
"Jet_btagCMVA")),
"reco_jets");
177 auto b_jets = filter(b_mva_filter, jets,
"reco_b_jets");
178 auto w_dijets = tup_filter<Jet,Jet>(w_mass_filter, combinations<Jet,2>(jets,
"reco_dijets"));
180 auto top_cands = cart_product<std::tuple<Jet,Jet>, Jet>(w_dijets, b_jets);
182 top_cands = tup_filter(dup_filter, top_cands);
184 auto& t_mass = GenFunction::register_function<float(std::tuple<Jet,Jet>,Jet)>(
"t_mass",
185 FUNC(([](
const std::tuple<Jet,Jet>& w,
const Jet& b){
186 return (std::get<0>(w).v+std::get<1>(w).v+b.v).M();
189 fv::map(t_mass, top_cands,
"reco_top_mass");
192 jets = apply(GenFunction::register_function<std::vector<Jet>(std::vector<TLorentzVector>,std::vector<int>)>(
"build_mcjets",
193 FUNC(([](
const std::vector<TLorentzVector>& vs,
const std::vector<int>& pdgid){
194 std::vector<Jet> jets;
195 for(
int i=0; i<vs.size(); i++){
196 jets.push_back(Jet::mc(vs[i],i, pdgid[i]));
199 }))), fv::tuple(lookup<std::vector<TLorentzVector>>(
"GenPart_4v"), lookup<std::vector<int>>(
"GenPart_pdgId")),
"mcjets");
201 b_jets = filter(b_pdgid_filter, jets);
203 w_dijets = tup_filter(qg_id_filter, combinations<Jet,2>(jets));
204 w_dijets = tup_filter(w_mass_filter, w_dijets);
206 top_cands = cart_product<std::tuple<Jet,Jet>, Jet>(w_dijets, b_jets);
208 top_cands = tup_filter(dup_filter, top_cands);
211 fv::map(t_mass, top_cands,
"mc_top_mass");
215 auto& inv_mass2 = GenFunction::register_function<float(Jet, Jet)>(
"inv_mass2",
216 FUNC(([] (
const Jet& j1,
const Jet& j2){
217 TLorentzVector sum = j1.v + j2.v;
218 return (
float)sum.M();
220 fv::map(inv_mass2, lookup<std::vector<std::tuple<Jet,Jet>>>(
"reco_dijets"),
"dijet_inv_mass");
224 count<float>(GenFunction::register_function<bool(float)>(
"bJet_Selection",
227 }))),
"Jet_btagCMVA",
"b_jet_count");
229 auto &is_electron = GenFunction::register_function<bool(int)>(
"is_electron",
230 FUNC(([](
int pdgId) {
231 return abs(pdgId) == 11;
233 auto &is_muon = GenFunction::register_function<bool(int)>(
"is_muon",
234 FUNC(([](
int pdgId) {
235 return abs(pdgId) == 13;
237 auto &is_lepton = GenFunction::register_function<bool(int)>(
"is_lepton",
238 FUNC(([ie=&is_electron, im=&is_muon](
int pdgId) {
239 return (*ie)(pdgId) || (*im)(pdgId);
242 count<int>(is_electron,
"GenPart_pdgId",
"genEle_count");
243 count<int>(is_muon,
"GenPart_pdgId",
"genMu_count");
244 count<int>(is_lepton,
"GenPart_pdgId",
"genLep_count");
246 count<int>(is_electron,
"LepGood_pdgId",
"recEle_count");
247 count<int>(is_muon,
"LepGood_pdgId",
"recMu_count");
248 count<int>(is_lepton,
"LepGood_pdgId",
"recLep_count");
251 fv::pair<int, int>(
"genEle_count",
"recEle_count",
"genEle_count_v_recEle_count");
252 fv::pair<int, int>(
"genMu_count",
"recMu_count",
"genMu_count_v_recMu_count");
253 fv::pair<int, int>(
"genLep_count",
"recLep_count",
"genLep_count_v_recLep_count");
267 obs_filter(
"trilepton", FUNC(([nLepGood=lookup<int>(
"nLepGood")]()
269 return dynamic_cast<Value<int>*
>(nLepGood)->get_value() == 3;
272 obs_filter(
"os-dilepton", FUNC(([LepGood_charge=lookup<vector<int>>(
"LepGood_charge")]()
274 auto& charges = LepGood_charge->get_value();
275 return charges.size()==2 && (charges[0] != charges[1]);
278 obs_filter(
"ss-dilepton", FUNC(([LepGood_charge=lookup<vector<int>>(
"LepGood_charge")]()
280 auto& charges = LepGood_charge->get_value();
281 return charges.size()==2 && (charges[0] == charges[1]);
286 void declare_containers(MiniTreeDataSet& mt){
288 mt.register_container<ContainerTH1<int>>(
"lepton_count",
"Lepton Multiplicity", lookup<int>(
"nLepGood"), 8, 0, 8);
290 mt.register_container<ContainerTH1Many<float>>(
"Jet_pt_dist",
"Jet P_T",
291 lookup<vector<float>>(
"Jet_pt"), 50, 0, 500,
293 mt.register_container<ContainerTH1Many<float>>(
"Jet_eta_dist",
"Jet Eta",
294 lookup<vector<float>>(
"Jet_eta"), 50, -3, 3,
296 mt.register_container<ContainerTH1Many<float>>(
"Jet_phi_dist",
"Jet Phi",
297 lookup<vector<float>>(
"Jet_phi"), 20, -PI, PI,
299 mt.register_container<ContainerTH1Many<float>>(
"Jet_mass_dist",
"Jet Mass",
300 lookup<vector<float>>(
"Jet_mass"), 50, 0, 200,
303 mt.register_container<ContainerTH1Many<float>>(
"dijet_inv_mass",
"Di-Jet Inv. Mass - All",
304 lookup<vector<float>>(
"dijet_inv_mass"), 100, 0, 500,
306 mt.register_container<ContainerTH1Many<float>>(
"dijet_inv_mass_osdilepton",
"Di-Jet Inv. Mass - OS Dilepton",
307 lookup<vector<float>>(
"dijet_inv_mass"), 100, 0, 500,
308 "Di-Jet Mass")->add_filter(lookup_obs_filter(
"os-dilepton"));
309 mt.register_container<ContainerTH1Many<float>>(
"dijet_inv_mass_ssdilepton",
"Di-Jet Inv. Mass - SS Dilepton",
310 lookup<vector<float>>(
"dijet_inv_mass"), 100, 0, 500,
311 "Di-Jet Mass")->add_filter(lookup_obs_filter(
"ss-dilepton"));
312 mt.register_container<ContainerTH1Many<float>>(
"dijet_inv_mass_trilepton",
"Di-Jet Inv. Mass - Trilepton",
313 lookup<vector<float>>(
"dijet_inv_mass"), 100, 0, 500,
314 "Di-Jet Mass")->add_filter(lookup_obs_filter(
"trilepton"));
316 mt.register_container<ContainerTH1Many<float>>(
"reco_top_mass",
"Reconstructed Top mass",
317 lookup<vector<float>>(
"reco_top_mass"), 100, 0, 500,
318 "Tri-jet invarient mass");
323 mt.register_container<ContainerTH1Many<float>>(
"mc_top_mass",
"Reconstructed MC Top mass",
324 lookup<vector<float>>(
"mc_top_mass"), 100, 0, 500,
325 "Tri-jet invarient mass");
328 mt.register_container<ContainerTH2<int>>(
"nLepvsnJet",
"Number of Leptons vs Number of Jets",
329 fv::pair<int, int>(
"nLepGood",
"nJet"),
331 "Number of Leptons",
"Number of Jets");
334 mt.register_container<ContainerTH2<int>>(
"genEle_count_v_recEle_count",
"Number of Generated Electrons v. Number of Reconstructed Electrons",
335 lookup<std::pair<int,int>>(
"genEle_count_v_recEle_count"),
336 10, 0, 10, 10, 0, 10,
337 "Generated Electrons",
"Reconstructed Electrons");
339 mt.register_container<ContainerTH2<int>>(
"genMu_count_v_recMu_count",
"Number of Generated Muons v. Number of Reconstructed Muons",
340 lookup<std::pair<int,int>>(
"genMu_count_v_recMu_count"),
341 10, 0, 10, 10, 0, 10,
342 "Generated Muons",
"Reconstructed Muons");
344 mt.register_container<ContainerTH2<int>>(
"genLep_count_v_recLep_count",
"Number of Generated Leptons v. Number of Reconstructed Leptons (e & mu only)",
345 lookup<std::pair<int,int>>(
"genLep_count_v_recLep_count"),
346 10, 0, 10, 10, 0, 10,
347 "Generated Leptons",
"Reconstructed Leptons");
349 mt.register_container<ContainerTH1<int>>(
"b_jet_count",
"B-Jet Multiplicity", lookup<int>(
"b_jet_count"), 10, 0, 10);
352 mt.register_container<ContainerTH1<int>>(
"jet_count_os_dilepton",
"Jet Multiplicity - OS Dilepton Events",
353 lookup<int>(
"nJet"), 14, 0, 14)->add_filter(lookup_obs_filter(
"os-dilepton"));
354 mt.register_container<ContainerTH1<int>>(
"jet_count_ss_dilepton",
"Jet Multiplicity - SS Dilepton Events",
355 lookup<int>(
"nJet"), 14, 0, 14)->add_filter(lookup_obs_filter(
"ss-dilepton"));
356 mt.register_container<ContainerTH1<int>>(
"jet_count_trilepton",
"Jet Multiplicity - Trilepton Events",
357 lookup<int>(
"nJet"), 14, 0, 14)->add_filter(lookup_obs_filter(
"trilepton"));
359 mt.register_container<ContainerTH1<int>>(
"nVert",
"Number of Primary Vertices", lookup<int>(
"nVert"), 50, 0, 50);
361 mt.register_container<
CounterMany<int>>(
"GenPart_pdgId_counter", lookup<vector<int>>(
"GenPart_pdgId"));
364 mt.register_container<Vector<std::vector< int >>>(
"GenPart_pdgId", lookup<std::vector< int >>(
"GenPart_pdgId"));
365 mt.register_container<Vector<std::vector< int >>>(
"GenPart_motherIndex", lookup<std::vector< int >>(
"GenPart_motherIndex"));
366 mt.register_container<Vector<std::vector< int >>>(
"GenPart_motherId", lookup<std::vector< int >>(
"GenPart_motherId"));
367 mt.register_container<Vector<std::vector<float>>>(
"GenPart_pt", lookup<std::vector<float>>(
"GenPart_pt"));
368 mt.register_container<Vector<std::vector<float>>>(
"GenPart_eta", lookup<std::vector<float>>(
"GenPart_eta"));
369 mt.register_container<Vector<std::vector<float>>>(
"GenPart_phi", lookup<std::vector<float>>(
"GenPart_phi"));
370 mt.register_container<Vector<std::vector<float>>>(
"GenPart_mass", lookup<std::vector<float>>(
"GenPart_mass"));
371 mt.register_container<Vector<std::vector<float>>>(
"GenPart_energy", lookup<std::vector<float>>(
"GenPart_energy"));
372 mt.register_container<Vector<std::vector< int >>>(
"GenPart_status", lookup<std::vector< int >>(
"GenPart_status"));
374 mt.register_container<Vector<vector< int >>>(
"LepGood_mcMatchId", lookup<vector< int >>(
"LepGood_mcMatchId"));
375 mt.register_container<Vector<vector< int >>>(
"LepGood_mcMatchPdgId", lookup<vector< int >>(
"LepGood_mcMatchPdgId"));
377 mt.register_container<Vector< int >>(
"run", lookup< int >(
"run") );
378 mt.register_container<Vector< int >>(
"lumi", lookup< int >(
"lumi"));
379 mt.register_container<Vector< int >>(
"evt", lookup< int >(
"evt") );
380 mt.register_container<Vector<float>>(
"xsec", lookup<float>(
"xsec"));
382 mt.register_container<Vector<std::vector< int >>>(
"Jet_mcMatchFlav", lookup<std::vector< int >>(
"Jet_mcMatchFlav"));
383 mt.register_container<Vector<std::vector< int >>>(
"Jet_mcMatchId", lookup<std::vector< int >>(
"Jet_mcMatchId"));
384 mt.register_container<Vector<std::vector< int >>>(
"Jet_mcFlavour", lookup<std::vector< int >>(
"Jet_mcFlavour"));
386 mt.register_container<Vector<std::vector<float>>>(
"Jet_pt", lookup<std::vector<float>>(
"Jet_pt"));
387 mt.register_container<Vector<std::vector<float>>>(
"Jet_eta", lookup<std::vector<float>>(
"Jet_eta"));
388 mt.register_container<Vector<std::vector<float>>>(
"Jet_phi", lookup<std::vector<float>>(
"Jet_phi"));
389 mt.register_container<Vector<std::vector<float>>>(
"Jet_mass", lookup<std::vector<float>>(
"Jet_mass"));
394 void run_analysis(
const std::string& input_filename,
bool silent){
395 gSystem->Load(
"libfilval.so");
396 auto replace_suffix = [](
const std::string& input,
const std::string& new_suffix){
397 return input.substr(0, input.find_last_of(
".")) + new_suffix;
399 string log_filename = replace_suffix(input_filename,
"_result.log");
400 fv::util::Log::init_logger(log_filename, fv::util::LogPriority::kLogDebug);
402 string output_filename = replace_suffix(input_filename,
"_result.root");
403 MiniTreeDataSet mt(output_filename, input_filename);
407 declare_containers(mt);
413 int main(
int argc,
char * argv[])
415 fv::util::ArgParser args(argc, argv);
416 if(!args.cmdOptionExists(
"-f")) {
417 cout <<
"Usage: ./main (-s) -f input_minitree.root" << endl;
420 bool silent = args.cmdOptionExists(
"-s");
421 string input_filename = args.getCmdOption(
"-f");
422 run_analysis(input_filename, silent);
Same as Counter but accepts multiple values per fill.
The namespace containing all filval classes and functions.