gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
Py8EventConversions.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
18 
19 #pragma once
20 
22 
23 //#define COLLIDERBIT_DEBUG
24 
25 
26 namespace Gambit
27 {
28 
29  namespace ColliderBit
30  {
31 
32  using namespace EventConversion;
33 
36  template<typename EventT>
37  void convertParticleEvent(const EventT& pevt, HEPUtils::Event& result, double antiktR, double jet_pt_min)
38  {
39  result.clear();
40 
41  std::vector<FJNS::PseudoJet> bhadrons; //< for input to FastJet b-tagging
42  std::vector<HEPUtils::Particle> bpartons, cpartons, tauCandidates;
43  HEPUtils::P4 pout; //< Sum of momenta outside acceptance
44  std::vector<FJNS::PseudoJet> jetparticles;
45 
46  for (decltype(pevt.size()) i = 0; i < pevt.size(); i++)
47  {
48  const auto &p = pevt[i];
49  const int pid = get_unified_pid(p);
50  const int apid = abs(pid);
51  const HEPUtils::P4 p4 = get_unified_momentum(p);
52 
53  //b, c and tau idenitification:
54 
55  // Find last b-hadrons in b decay chains as the best proxy for b-tagging
57  if (apid == 5)
58  {
59  bool isGoodB = true;
60  std::vector<int> childIDs;
61  get_unified_child_ids(p, pevt, childIDs);
62  for (int childID : childIDs)
63  {
64  if (abs(childID) == 5) isGoodB = false;
65  }
66  if (isGoodB) bpartons.push_back(HEPUtils::Particle(p4,pid));
67  }
68 
69  // Find last c-hadrons in decay chains as the best proxy for c-tagging
71  if (apid == 4)
72  {
73  bool isGoodC = true;
74  std::vector<int> childIDs;
75  get_unified_child_ids(p, pevt, childIDs);
76  for (int childID : childIDs)
77  {
78  if (abs(childID) == 4) isGoodC = false;
79  }
80  if (isGoodC) cpartons.push_back(HEPUtils::Particle(p4,pid));
81  }
82 
83  // Find tau candidates
84  if (apid == MCUtils::PID::TAU)
85  {
86  bool isGoodTau=true;
87  std::vector<int> childIDs;
88  get_unified_child_ids(p, pevt, childIDs);
89  int abschildID;
90  for (int childID : childIDs)
91  {
92  // Veto leptonic taus
94  abschildID = abs(childID);
95  if (abschildID == MCUtils::PID::ELECTRON || abschildID == MCUtils::PID::MUON ||
96  abschildID == MCUtils::PID::WPLUSBOSON || abschildID == MCUtils::PID::TAU)
97  {
98  isGoodTau = false;
99  }
100  }
101  if (isGoodTau) tauCandidates.push_back(HEPUtils::Particle(p4, pid));
102  }
103 
104  //We only want final state particles:
105  if (!get_unified_isFinal(p)) continue;
106 
107  //Check there's no partons.
108  if (pid == 21 || abs(pid) <= 6)
109  {
110  std::ostringstream sid;
111  bool gotmother = false;
112  //HepMC seems to have no equivalent of the .mother1, .mother2 call, so the HepMC3 mother function will just
113  //return 0, and gotmother will always be false - which means it won't try to print non-existent event info.
114  if (get_unified_mother1(p) != 0 ){gotmother = true; sid << get_unified_mother1_pid(p, pevt);}
115  if (get_unified_mother2(p) != 0 ){gotmother = true; sid << get_unified_mother2_pid(p, pevt);}
116  if (gotmother) sid << " -> ";
117  sid << pid;
118  ColliderBit_error().forced_throw(LOCAL_INFO, "Found final-state parton " + sid.str() + " in particle-level event converter: "
119  "reconfigure your generator to include hadronization, or Gambit to use the partonic event converter.");
120  }
121 
122  // Add particle outside ATLAS/CMS acceptance to MET and then ignore said particle.
124  if (std::abs(get_unified_eta(p)) > 5.0)
125  {
126  pout += p4;
127  continue;
128  }
129 
130  // Promptness: for leptons and photons we're only interested if they don't come from hadron/tau decays
131  const bool prompt = !get_unified_fromHadron(p, pevt, i);
132  const bool visible = MCUtils::PID::isStrongInteracting(pid) || MCUtils::PID::isEMInteracting(pid);
133 
134  // Add prompt and invisible particles as individual particles
135  if (prompt || !visible)
136  {
137  HEPUtils::Particle* gp = new HEPUtils::Particle(p4, pid);
138  gp->set_prompt();
139  result.add_particle(gp);
140  }
141 
142  // All particles other than invisibles and muons are jet constituents
143  if (visible && apid != MCUtils::PID::MUON)
144  {
145  jetparticles.push_back(get_unified_pseudojet(p));
146  }
147  }
148 
151  const FJNS::JetDefinition jet_def(FJNS::antikt_algorithm, antiktR);
152  FJNS::ClusterSequence cseq(jetparticles, jet_def);
153  std::vector<FJNS::PseudoJet> pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min));
154 
158  for (auto& pj : pjets)
159  {
160  HEPUtils::P4 jetMom = HEPUtils::mk_p4(pj);
162  bool isB = false;
163  for (HEPUtils::Particle& pb : bpartons)
164  {
165  if (jetMom.deltaR_eta(pb.mom()) < 0.4)
166  {
167  isB = true;
168  break;
169  }
170  }
171 
172  bool isC = false;
173  for (HEPUtils::Particle& pc : cpartons)
174  {
175  if (jetMom.deltaR_eta(pc.mom()) < 0.4)
176  {
177  isC = true;
178  break;
179  }
180  }
181 
182  bool isTau = false;
183  for (HEPUtils::Particle& ptau : tauCandidates)
184  {
185  if (jetMom.deltaR_eta(ptau.mom()) < 0.5)
186  {
187  isTau = true;
188  break;
189  }
190  }
191 
192  // Add to the event (use jet momentum for tau)
193  if (isTau)
194  {
195  HEPUtils::Particle* gp = new HEPUtils::Particle(HEPUtils::mk_p4(pj), MCUtils::PID::TAU);
196  gp->set_prompt();
197  result.add_particle(gp);
198  }
199  result.add_jet(new HEPUtils::Jet(HEPUtils::mk_p4(pj), isB, isC));
200  }
201 
203  //
204  // From balance of all visible momenta (requires isolation)
205  // const std::vector<Particle*> visibles = result.visible_particles();
206  // HEPUtils::P4 pvis;
207  // for (size_t i = 0; i < visibles.size(); ++i)
208  // {
209  // pvis += visibles[i]->mom();
210  // }
211  // for (size_t i = 0; i < result.jets.size(); ++i)
212  // {
213  // pvis += result.jets[i]->mom();
214  // }
215  // set_missingmom(-pvis);
216  //
217  // From sum of invisibles, including those out of range
218  for (size_t i = 0; i < result.invisible_particles().size(); ++i)
219  {
220  pout += result.invisible_particles()[i]->mom();
221  }
222  result.set_missingmom(pout);
223 
224  #ifdef COLLIDERBIT_DEBUG
225  // Print event summary
226  cout << " MET = " << result.met() << " GeV" << endl;
227  cout << " #e = " << result.electrons().size() << endl;
228  cout << " #mu = " << result.muons().size() << endl;
229  cout << " #tau = " << result.taus().size() << endl;
230  cout << " #jet = " << result.jets().size() << endl;
231  cout << " #pho = " << result.photons().size() << endl;
232  cout << endl;
233  #endif
234  }
235 
236 
238  template<typename EventT>
239  void convertPartonEvent(const EventT& pevt, HEPUtils::Event& result, double antiktR, double jet_pt_min)
240  {
241  result.clear();
242 
243  std::vector<HEPUtils::Particle> tauCandidates;
244 
245  // Make a first pass of non-final particles to gather taus
246  for (int i = 0; i < pevt.size(); ++i)
247  {
248  const auto& p = pevt[i];
249 
250  // Find last tau in prompt tau replica chains as a proxy for tau-tagging
251  if (p.idAbs() == MCUtils::PID::TAU) {
252  std::vector<int> tauDaughterList = p.daughterList();
253  HEPUtils::P4 tmpMomentum;
254  bool isGoodTau=true;
255 
256  for (size_t daughter = 0; daughter < tauDaughterList.size(); daughter++)
257  {
258  const auto& pDaughter = pevt[tauDaughterList[daughter]];
259  int daughterID = pDaughter.idAbs();
260  if (daughterID == MCUtils::PID::ELECTRON || daughterID == MCUtils::PID::MUON ||
261  daughterID == MCUtils::PID::WPLUSBOSON || daughterID == MCUtils::PID::TAU)
262  isGoodTau = false;
263  if (daughterID != MCUtils::PID::TAU) tmpMomentum += mk_p4(pDaughter.p());
264  }
265 
266  if (isGoodTau) {
267  tauCandidates.push_back(HEPUtils::Particle(mk_p4(p.p()), p.id()));
268  }
269  }
270  }
271 
272  std::vector<FJNS::PseudoJet> jetparticles; //< Pseudojets for input to FastJet
273  HEPUtils::P4 pout; //< Sum of momenta outside acceptance
274 
275  // Make a single pass over the event to gather final leptons, partons, and photons
276  for (int i = 0; i < pevt.size(); ++i)
277  {
278  const auto& p = pevt[i];
279 
280  // We only use "final" partons, i.e. those with no children. So Py8 must have hadronization disabled
281  if (!p.isFinal()) continue;
282 
283  // Only consider partons within ATLAS/CMS acceptance
285  if (std::abs(p.eta()) > 5.0)
286  {
287  pout += mk_p4(p.p());
288  continue;
289  }
290 
291  // Find electrons/muons/taus/photons to be treated as prompt (+ invisibles)
294  const bool prompt = isFinalPhoton(i, pevt) || (isFinalLepton(i, pevt)); // && std::abs(p.id()) != MCUtils::PID::TAU);
295  const bool visible = MCUtils::PID::isStrongInteracting(p.id()) || MCUtils::PID::isEMInteracting(p.id());
296  if (prompt || !visible)
297  {
298  HEPUtils::Particle* gp = new HEPUtils::Particle(mk_p4(p.p()), p.id());
299  gp->set_prompt();
300  result.add_particle(gp);
301  }
302 
303  // Everything other than invisibles and muons, including taus & partons are jet constituents
305  // if (visible && (isFinalParton(i, pevt) || isFinalTau(i, pevt))) {
306  if (visible && p.idAbs() != MCUtils::PID::MUON)
307  {
308  FJNS::PseudoJet pj = mk_pseudojet(p.p());
309  //pj.set_user_index(std::abs(p.id()));
310  jetparticles.push_back(pj);
311  }
312 
313  }
314 
317  const FJNS::JetDefinition jet_def(FJNS::antikt_algorithm, antiktR);
318  FJNS::ClusterSequence cseq(jetparticles, jet_def);
319  std::vector<FJNS::PseudoJet> pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min));
320 
321  // Add to the event, with b-tagging info"
322  for (const FJNS::PseudoJet& pj : pjets)
323  {
324  // Do jet b-tagging, etc. by looking for b quark constituents (i.e. user index = |parton ID| = 5)
326  const bool isB = HEPUtils::any(pj.constituents(),
327  [](const FJNS::PseudoJet& c){ return c.user_index() == MCUtils::PID::BQUARK; });
328  const bool isC = HEPUtils::any(pj.constituents(),
329  [](const FJNS::PseudoJet& c){ return c.user_index() == MCUtils::PID::CQUARK; });
330  result.add_jet(new HEPUtils::Jet(HEPUtils::mk_p4(pj), isB, isC));
331 
332  bool isTau=false;
333  for (auto& ptau : tauCandidates)
334  {
335  HEPUtils::P4 jetMom = HEPUtils::mk_p4(pj);
336  if (jetMom.deltaR_eta(ptau.mom()) < 0.5)
337  {
338  isTau=true;
339  break;
340  }
341  }
342  // Add to the event (use jet momentum for tau)
343  if (isTau)
344  {
345  HEPUtils::Particle* gp = new HEPUtils::Particle(HEPUtils::mk_p4(pj), MCUtils::PID::TAU);
346  gp->set_prompt();
347  result.add_particle(gp);
348  }
349  }
350 
352  //
353  // From balance of all visible momenta (requires isolation)
354  // const std::vector<Particle*> visibles = result.visible_particles();
355  // HEPUtils::P4 pvis;
356  // for (size_t i = 0; i < visibles.size(); ++i) {
357  // pvis += visibles[i]->mom();
358  // }
359  // for (size_t i = 0; i < result.jets.size(); ++i) {
360  // pvis += result.jets[i]->mom();
361  // }
362  // set_missingmom(-pvis);
363  //
364  // From sum of invisibles, including those out of range
365  for (const HEPUtils::Particle* p : result.invisible_particles()) pout += p->mom();
366  result.set_missingmom(pout);
367  }
368 
369  }
370 
371 }
void get_unified_child_ids(ParticleP &p, EventT &pevt, std::vector< int > &unified_child_id_results)
#define LOCAL_INFO
Definition: local_info.hpp:34
FJNS::PseudoJet mk_pseudojet(const Vec4T &p)
Definition: Py8Utils.hpp:39
HEPUtils::P4 mk_p4(const Vec4T &p)
Definition: Py8Utils.hpp:45
void convertPartonEvent(const EventT &pevt, HEPUtils::Event &result, double antiktR, double jet_pt_min)
Convert a partonic (no hadrons) EventT into an unsmeared HEPUtils::Event.
int get_unified_mother2_pid(ParticleP &p, EventT &pevt)
HEPUtils::P4 get_unified_momentum(ParticleP p)
void convertParticleEvent(const EventT &pevt, HEPUtils::Event &result, double antiktR, double jet_pt_min)
Convert a hadron-level EventT into an unsmeared HEPUtils::Event.
Helper functions for converting between different event types.
FJNS::PseudoJet get_unified_pseudojet(ParticleP p)
int get_unified_mother1_pid(ParticleP &p, EventT &pevt)
bool isFinalPhoton(int n, const EventT &evt)
Definition: Py8Utils.hpp:172
bool get_unified_fromHadron(ParticleP &, const EventT &pevt, int i)
TODO: see if we can use this one:
Definition: Analysis.hpp:33
bool isFinalLepton(int n, const EventT &evt)
Definition: Py8Utils.hpp:185