PYTHIA  8.312
JetMatching.h
1 // JetMatching.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Authors: Richard Corke (implementation of MLM matching as
7 // in Alpgen for Alpgen input)
8 // and Stephen Mrenna (implementation of MLM-style matching as
9 // in Madgraph for Alpgen or Madgraph 5 input.)
10 // and Simon de Visscher, Stefan Prestel (implementation of shower-kT
11 // MLM-style matching and flavour treatment for Madgraph input)
12 // and Stefan Prestel (FxFx NLO jet matching with aMC@NLO.)
13 // This file provides the classes to perform MLM matching of
14 // Alpgen or MadGraph 5 input.
15 // Example usage is shown in main32.cc, and further details
16 // can be found in the 'Jet Matching Style' manual page.
17 
18 #ifndef Pythia8_JetMatching_H
19 #define Pythia8_JetMatching_H
20 
21 // Includes
22 #include "Pythia8/Pythia.h"
23 #include "Pythia8Plugins/GeneratorInput.h"
24 
25 namespace Pythia8 {
26 
27 //==========================================================================
28 
29 class HJSlowJet: public SlowJet {
30 
31  public:
32 
33  HJSlowJet(int powerIn, double Rin, double pTjetMinIn = 0.,
34  double etaMaxIn = 25., int selectIn = 1, int massSetIn = 2,
35  SlowJetHook* sjHookPtrIn = 0, bool useFJcoreIn = false,
36  bool useStandardRin = true) :
37  SlowJet(powerIn, Rin, pTjetMinIn, etaMaxIn, selectIn, massSetIn,
38  sjHookPtrIn, useFJcoreIn, useStandardRin) {}
39 
40  // Note: The functions below have been made public to ease the generation
41  // of Python bindings.
42  //protected:
43 
44  void findNext();
45 
46 };
47 
48 //--------------------------------------------------------------------------
49 
50 // Find next cluster pair to join.
51 
53 
54  // Find smallest of diB, dij.
55  if (clSize > 0) {
56  iMin = 0;
57  jMin = -1;
58  dMin = 1.0/TINY;
59  // Remove the possibility of choosing a beam clustering
60  for (int i = 1; i < clSize; ++i) {
61  for (int j = 0; j < i; ++j) {
62  if (dij[i*(i-1)/2 + j] < dMin) {
63  iMin = i;
64  jMin = j;
65  dMin = dij[i*(i-1)/2 + j];
66  }
67  }
68  }
69 
70  // If no clusters left then instead default values.
71  } else {
72  iMin = -1;
73  jMin = -1;
74  dMin = 0.;
75  }
76 
77 }
78 
79 //==========================================================================
80 
81 // Declaration of main JetMatching class to perform MLM matching.
82 // Note that it is defined with virtual inheritance, so that it can
83 // be combined with other UserHooks classes, see e.g. main33.cc.
84 
85 class JetMatching : virtual public UserHooks {
86 
87 public:
88 
89  // Constructor and destructor
90  JetMatching() : cellJet(nullptr), slowJet(nullptr), slowJetHard(nullptr),
91  hjSlowJet(nullptr) {}
93  if (cellJet) delete cellJet;
94  if (slowJet) delete slowJet;
95  if (slowJetHard) delete slowJetHard;
96  if (hjSlowJet) delete hjSlowJet;
97  // Print error statistics before exiting. Printing code
98  // basically copied from Info class.
99  // Header.
100  cout << "\n *------- JetMatching Error and Warning Messages Statistics"
101  << " -----------------------------------------------------* \n"
102  << " | "
103  << " | \n"
104  << " | times message "
105  << " | \n"
106  << " | "
107  << " | \n";
108 
109  // Loop over all messages
110  map<string, int>::iterator messageEntry = messages.begin();
111  if (messageEntry == messages.end())
112  cout << " | 0 no errors or warnings to report "
113  << " | \n";
114  while (messageEntry != messages.end()) {
115  // Message printout.
116  string temp = messageEntry->first;
117  int len = temp.length();
118  temp.insert( len, max(0, 102 - len), ' ');
119  cout << " | " << setw(6) << messageEntry->second << " "
120  << temp << " | \n";
121  ++messageEntry;
122  }
123 
124  // Done.
125  cout << " | "
126  << " | \n"
127  << " *------- End JetMatching Error and Warning Messages "
128  << "Statistics -------------------------------------------------* "
129  << endl;
130  }
131 
132  // Initialisation
133  virtual bool initAfterBeams() = 0;
134 
135  // Process level vetos
136  bool canVetoProcessLevel() { return doMerge; }
137  bool doVetoProcessLevel(Event& process) {
138  eventProcessOrig = process;
139  return false;
140  }
141 
142  // Parton level vetos (before beam remnants and resonance decays)
143  bool canVetoPartonLevelEarly() { return doMerge; }
144  bool doVetoPartonLevelEarly(const Event& event);
145 
146  // Shower step vetoes (after the first emission, for Shower-kT scheme)
147  int numberVetoStep() {return 1;}
148  bool canVetoStep() { return false; }
149  bool doVetoStep(int, int, int, const Event& ) { return false; }
150 
151  // Note: The functions below have been made public to ease the generation
152  // of Python bindings.
153  //protected:
154 
155  // Different steps of the matching algorithm.
156  virtual void sortIncomingProcess(const Event &)=0;
157  virtual void jetAlgorithmInput(const Event &, int)=0;
158  virtual void runJetAlgorithm()=0;
159  virtual bool matchPartonsToJets(int)=0;
160  virtual int matchPartonsToJetsLight()=0;
161  virtual int matchPartonsToJetsHeavy()=0;
162 
163  // Print a message the first few times. Insert in database.
164  void errorMsg(string messageIn) {
165  // Recover number of times message occured. Also inserts new string.
166  int times = messages[messageIn];
167  ++messages[messageIn];
168  // Print message the first few times.
169  if (times < TIMESTOPRINT) cout << " PYTHIA " << messageIn << endl;
170  }
171 
172 protected:
173 
174  // Constants to be changed for debug printout or extra checks.
175  static const bool MATCHINGDEBUG, MATCHINGCHECK;
176 
177  enum vetoStatus { NONE, LESS_JETS, MORE_JETS, HARD_JET,
178  UNMATCHED_PARTON, INCLUSIVE_VETO };
179  enum partonTypes { ID_CHARM=4, ID_BOT=5, ID_TOP=6, ID_LEPMIN=11,
180  ID_LEPMAX=16, ID_GLUON=21, ID_PHOTON=22 };
181 
182  // Master switch for merging
183  bool doMerge;
184  // Switch for merging in the shower-kT scheme. Needed here because
185  // the scheme uses different UserHooks functionality.
187 
188  // Maximum and current number of jets
189  int nJetMax, nJet;
190 
191  // Jet algorithm parameters
193  double eTjetMin, coneRadius, etaJetMax, etaJetMaxAlgo;
194 
195  // Internal jet algorithms
197  SlowJet* slowJet;
198  SlowJet* slowJetHard;
199  HJSlowJet* hjSlowJet;
200 
201  // SlowJet specific
203 
204  // Event records to store original incoming process, final-state of the
205  // incoming process and what will be passed to the jet algorithm.
206  // Not completely necessary to store all steps, but makes tracking the
207  // steps of the algorithm a lot easier.
208  Event eventProcessOrig, eventProcess, workEventJet;
209 
210  // Sort final-state of incoming process into light/heavy jets and 'other'
211  vector<int> typeIdx[3];
212  set<int> typeSet[3];
213 
214  // Momenta output of jet algorithm (to provide same output regardless of
215  // the selected jet algorithm)
216  vector<Vec4> jetMomenta;
217 
218  // CellJet specific
219  int nEta, nPhi;
220  double eTseed, eTthreshold;
221 
222  // Merging procedure parameters
223  int jetAllow, jetMatch, exclusiveMode;
224  double coneMatchLight, coneRadiusHeavy, coneMatchHeavy;
225  bool exclusive;
226 
227  // Store the minimum eT/pT of matched light jets
228  double eTpTlightMin;
229 
230  // Map for all error messages.
231  map<string, int> messages;
232  // Number of times the same error message is repeated, unless overridden.
233  static const int TIMESTOPRINT = 1;
234 
235 };
236 
237 //==========================================================================
238 
239 // Declaration of main UserHooks class to perform Alpgen matching.
240 
241 class JetMatchingAlpgen : virtual public JetMatching {
242 
243 public:
244 
245  // Constructor and destructor
247  ~JetMatchingAlpgen() { }
248 
249  // Initialisation
250  bool initAfterBeams();
251 
252 private:
253 
254  // Different steps of the matching algorithm.
255  void sortIncomingProcess(const Event &);
256  void jetAlgorithmInput(const Event &, int);
257  void runJetAlgorithm();
258  bool matchPartonsToJets(int);
259  int matchPartonsToJetsLight();
260  int matchPartonsToJetsHeavy();
261 
262  // Sorting utility
263  void sortTypeIdx(vector < int > &vecIn);
264 
265  // Constants
266  static const double GHOSTENERGY, ZEROTHRESHOLD;
267 
268 };
269 
270 //==========================================================================
271 
272 // Declaration of main UserHooks class to perform Madgraph matching.
273 
274 class JetMatchingMadgraph : virtual public JetMatching {
275 
276 public:
277 
278  // Constructor and destructor
279  JetMatchingMadgraph() : slowJetDJR(nullptr) { }
280  ~JetMatchingMadgraph() { if (slowJetDJR) delete slowJetDJR; }
281 
282  // Initialisation
283  bool initAfterBeams();
284 
285  // Process level vetos
286  bool canVetoProcessLevel() { return doMerge; }
287  bool doVetoProcessLevel(Event& process);
288 
289  // Shower step vetoes (after the first emission, for Shower-kT scheme)
290  int numberVetoStep() {return 1;}
291  bool canVetoStep() { return doShowerKt; }
292  bool doVetoStep(int, int, int, const Event& );
293 
294  // Jet algorithm to access the jet separations in the cleaned event
295  // after showering.
297  // Functions to return the jet clustering scales and number of ME partons.
298  // These are useful to investigate the matching systematics.
299  vector<double> getDJR() { return DJR;}
300  pair<int,int> nMEpartons() { return nMEpartonsSave;}
301 
302  // For systematic variations of the jet matching parameters, it is helpful
303  // to decouple the jet matching veto from the internal book-keeping. The
304  // veto can then be applied in hindsight by an expert user. The functions
305  // below return all the necessary information to do this.
306  Event getWorkEventJet() { return workEventJetSave; }
307  Event getProcessSubset() { return processSubsetSave; }
308  bool getExclusive() { return exclusive; }
309  double getPTfirst() { return pTfirstSave; }
310 
311  // Note: The functions below have been made public to ease the generation
312  // of Python bindings.
313  //protected:
314 
315  // Different steps of the matching algorithm.
316  void sortIncomingProcess(const Event &);
317  void jetAlgorithmInput(const Event &, int);
318  void runJetAlgorithm();
319  bool matchPartonsToJets(int);
320  int matchPartonsToJetsLight();
321  int matchPartonsToJetsHeavy();
322  int matchPartonsToJetsOther();
323  bool doShowerKtVeto(double pTfirst);
324 
325  // Functions to clear and set the jet clustering scales.
326  void clearDJR() { DJR.resize(0);}
327  void setDJR( const Event& event);
328  // Functions to clear and set the jet clustering scales.
329  void clear_nMEpartons() { nMEpartonsSave.first = nMEpartonsSave.second =-1;}
330  void set_nMEpartons( const int nOrig, const int nMatch) {
331  clear_nMEpartons();
332  nMEpartonsSave.first = nOrig;
333  nMEpartonsSave.second = nMatch;
334  };
335 
336  // Function to get the current number of partons in the Born state, as
337  // read from LHE.
338  int npNLO();
339 
340 private:
341 
342  // Stored values of all inputs necessary to perform the jet matching, as
343  // needed when the veto is applied externally.
344  Event processSubsetSave;
345  Event workEventJetSave;
346  double pTfirstSave;
347 
348  // Save if code should apply the veto, or simply store the things necessary
349  // to perform the veto externally.
350  bool performVeto;
351 
352  // Variables.
353  vector<int> origTypeIdx[3];
354  int nQmatch;
355  double qCut, qCutSq, clFact;
356  bool doFxFx;
357  int nPartonsNow;
358  double qCutME, qCutMESq;
359 
360  // Vector to store the jet clustering scales.
361  vector<double> DJR;
362  // Pair of integers giving the number of ME partons read from LHEF and used
363  // in the matching (can be different if some partons should not be matched)
364  pair<int,int> nMEpartonsSave;
365 
366 };
367 
368 //==========================================================================
369 
370 // Main implementation of JetMatching class.
371 // This may be split out to a separate C++ file if desired,
372 // but currently included here for ease of use.
373 
374 //--------------------------------------------------------------------------
375 
376 // Constants to be changed for debug printout or extra checks.
377 const bool JetMatching::MATCHINGDEBUG = false;
378 const bool JetMatching::MATCHINGCHECK = false;
379 
380 //--------------------------------------------------------------------------
381 
382 // Early parton level veto (before beam remnants and resonance showers)
383 
384 inline bool JetMatching::doVetoPartonLevelEarly(const Event& event) {
385 
386  // 1) Sort the original incoming process. After this step is performed,
387  // the following assignments have been made:
388  // eventProcessOrig - the original incoming process
389  // eventProcess - the final-state of the incoming process with
390  // resonance decays removed (and resonances
391  // themselves now with positive status code)
392  // typeIdx[0/1/2] - Indices into 'eventProcess' of
393  // light jets/heavy jets/other
394  // typeSet[0/1/2] - Indices into 'event' of light jets/heavy jets/other
395  // workEvent - partons from the hardest subsystem + ISR + FSR only
396  sortIncomingProcess(event);
397 
398  // For the shower-kT scheme, do not perform any veto here, as any vetoing
399  // will already have taken place in doVetoStep.
400  if ( doShowerKt ) return false;
401 
402  // Debug printout.
403  if (MATCHINGDEBUG) {
404  // Begin
405  cout << endl << "-------- Begin Madgraph Debug --------" << endl;
406  // Original incoming process
407  cout << endl << "Original incoming process:";
408  eventProcessOrig.list();
409  // Final-state of original incoming process
410  cout << endl << "Final-state incoming process:";
411  eventProcess.list();
412  // List categories of sorted particles
413  for (size_t i = 0; i < typeIdx[0].size(); i++)
414  cout << ((i == 0) ? "Light jets: " : ", ") << setw(3) << typeIdx[0][i];
415  if( typeIdx[0].size()== 0 )
416  cout << "Light jets: None";
417 
418  for (size_t i = 0; i < typeIdx[1].size(); i++)
419  cout << ((i == 0) ? "\nHeavy jets: " : ", ") << setw(3) << typeIdx[1][i];
420  for (size_t i = 0; i < typeIdx[2].size(); i++)
421  cout << ((i == 0) ? "\nOther: " : ", ") << setw(3) << typeIdx[2][i];
422  // Full event at this stage
423  cout << endl << endl << "Event:";
424  event.list();
425  // Work event (partons from hardest subsystem + ISR + FSR)
426  cout << endl << "Work event:";
427  workEvent.list();
428  }
429 
430  // 2) Light/heavy jets: iType = 0 (light jets), 1 (heavy jets)
431  int iTypeEnd = (typeIdx[2].empty()) ? 2 : 3;
432  for (int iType = 0; iType < iTypeEnd; iType++) {
433 
434  // 2a) Find particles which will be passed from the jet algorithm.
435  // Input from 'workEvent' and output in 'workEventJet'.
436  jetAlgorithmInput(event, iType);
437 
438  // Debug printout.
439  if (MATCHINGDEBUG) {
440  // Jet algorithm event
441  cout << endl << "Jet algorithm event (iType = " << iType << "):";
442  workEventJet.list();
443  }
444 
445  // 2b) Run jet algorithm on 'workEventJet'.
446  // Output is stored in jetMomenta.
447  runJetAlgorithm();
448 
449  // 2c) Match partons to jets and decide if veto is necessary
450  if (matchPartonsToJets(iType) == true) {
451  // Debug printout.
452  if (MATCHINGDEBUG) {
453  cout << endl << "Event vetoed" << endl
454  << "---------- End MLM Debug ----------" << endl;
455  }
456  return true;
457  }
458  }
459 
460  // Debug printout.
461  if (MATCHINGDEBUG) {
462  cout << endl << "Event accepted" << endl
463  << "---------- End MLM Debug ----------" << endl;
464  }
465 
466  // If we reached here, then no veto
467  return false;
468 
469 }
470 
471 //==========================================================================
472 
473 // Main implementation of Alpgen UserHooks class.
474 // This may be split out to a separate C++ file if desired,
475 // but currently included here for ease of use.
476 
477 //--------------------------------------------------------------------------
478 
479 // Constants: could be changed here if desired, but normally should not.
480 // These are of technical nature, as described for each.
481 
482 // The energy of ghost particles. For technical reasons, this cannot be
483 // set arbitrarily low, see 'Particle::TINY' in 'Event.cc' for details.
484 const double JetMatchingAlpgen::GHOSTENERGY = 1e-15;
485 
486 // A zero threshold value for double comparisons.
487 const double JetMatchingAlpgen::ZEROTHRESHOLD = 1e-10;
488 
489 //--------------------------------------------------------------------------
490 
491 // Function to sort typeIdx vectors into descending eT/pT order.
492 // Uses a selection sort, as number of partons generally small
493 // and so efficiency not a worry.
494 
495 inline void JetMatchingAlpgen::sortTypeIdx(vector < int > &vecIn) {
496  for (size_t i = 0; i < vecIn.size(); i++) {
497  size_t jMax = i;
498  double vMax = (jetAlgorithm == 1) ?
499  eventProcess[vecIn[i]].eT() :
500  eventProcess[vecIn[i]].pT();
501  for (size_t j = i + 1; j < vecIn.size(); j++) {
502  double vNow = (jetAlgorithm == 1)
503  ? eventProcess[vecIn[j]].eT() : eventProcess[vecIn[j]].pT();
504  if (vNow > vMax) {
505  vMax = vNow;
506  jMax = j;
507  }
508  }
509  if (jMax != i) swap(vecIn[i], vecIn[jMax]);
510  }
511 }
512 
513 //--------------------------------------------------------------------------
514 
515 // Initialisation routine automatically called from Pythia::init().
516 // Setup all parts needed for the merging.
517 
519 
520  // Read in parameters
521  doMerge = settingsPtr->flag("JetMatching:merge");
522  jetAlgorithm = settingsPtr->mode("JetMatching:jetAlgorithm");
523  nJet = settingsPtr->mode("JetMatching:nJet");
524  nJetMax = settingsPtr->mode("JetMatching:nJetMax");
525  eTjetMin = settingsPtr->parm("JetMatching:eTjetMin");
526  coneRadius = settingsPtr->parm("JetMatching:coneRadius");
527  etaJetMax = settingsPtr->parm("JetMatching:etaJetMax");
528  doShowerKt = settingsPtr->flag("JetMatching:doShowerKt");
529 
530  // Use etaJetMax + coneRadius in input to jet algorithms
531  etaJetMaxAlgo = etaJetMax + coneRadius;
532 
533  // CellJet specific
534  nEta = settingsPtr->mode("JetMatching:nEta");
535  nPhi = settingsPtr->mode("JetMatching:nPhi");
536  eTseed = settingsPtr->parm("JetMatching:eTseed");
537  eTthreshold = settingsPtr->parm("JetMatching:eTthreshold");
538 
539  // SlowJet specific
540  slowJetPower = settingsPtr->mode("JetMatching:slowJetPower");
541  coneMatchLight = settingsPtr->parm("JetMatching:coneMatchLight");
542  coneRadiusHeavy = settingsPtr->parm("JetMatching:coneRadiusHeavy");
543  if (coneRadiusHeavy < 0.) coneRadiusHeavy = coneRadius;
544  coneMatchHeavy = settingsPtr->parm("JetMatching:coneMatchHeavy");
545 
546  // Matching procedure
547  jetAllow = settingsPtr->mode("JetMatching:jetAllow");
548  jetMatch = settingsPtr->mode("JetMatching:jetMatch");
549  exclusiveMode = settingsPtr->mode("JetMatching:exclusive");
550 
551  // If not merging, then done
552  if (!doMerge) return true;
553 
554  // Exclusive mode; if set to 2, then set based on nJet/nJetMax
555  if (exclusiveMode == 2) {
556 
557  // No nJet or nJetMax, so default to exclusive mode
558  if (nJet < 0 || nJetMax < 0) {
559  errorMsg("Warning in JetMatchingAlpgen:init: "
560  "missing jet multiplicity information; running in exclusive mode");
561  exclusive = true;
562 
563  // Inclusive if nJet == nJetMax, exclusive otherwise
564  } else {
565  exclusive = (nJet == nJetMax) ? false : true;
566  }
567 
568  // Otherwise, just set as given
569  } else {
570  exclusive = (exclusiveMode == 0) ? false : true;
571  }
572 
573  // Initialise chosen jet algorithm. CellJet.
574  if (jetAlgorithm == 1) {
575 
576  // Extra options for CellJet. nSel = 1 means that all final-state
577  // particles are taken and we retain control of what to select.
578  // smear/resolution/upperCut are not used and are set to default values.
579  int nSel = 2, smear = 0;
580  double resolution = 0.5, upperCut = 2.;
581  cellJet = new CellJet(etaJetMaxAlgo, nEta, nPhi, nSel,
582  smear, resolution, upperCut, eTthreshold);
583 
584  // SlowJet
585  } else if (jetAlgorithm == 2) {
586  slowJet = new SlowJet(slowJetPower, coneRadius, eTjetMin, etaJetMaxAlgo);
587  }
588 
589  // Check the jetMatch parameter; option 2 only works with SlowJet
590  if (jetAlgorithm == 1 && jetMatch == 2) {
591  errorMsg("Warning in JetMatchingAlpgen:init: "
592  "jetMatch = 2 only valid with SlowJet algorithm. "
593  "Reverting to jetMatch = 1");
594  jetMatch = 1;
595  }
596 
597  // Setup local event records
598  eventProcessOrig.init("(eventProcessOrig)", particleDataPtr);
599  eventProcess.init("(eventProcess)", particleDataPtr);
600  workEventJet.init("(workEventJet)", particleDataPtr);
601 
602  // Print information
603  string jetStr = (jetAlgorithm == 1) ? "CellJet" :
604  (slowJetPower == -1) ? "anti-kT" :
605  (slowJetPower == 0) ? "C/A" :
606  (slowJetPower == 1) ? "kT" : "unknown";
607  string modeStr = (exclusive) ? "exclusive" : "inclusive";
608  stringstream nJetStr, nJetMaxStr;
609  if (nJet >= 0) nJetStr << nJet; else nJetStr << "unknown";
610  if (nJetMax >= 0) nJetMaxStr << nJetMax; else nJetMaxStr << "unknown";
611  cout << endl
612  << " *------- MLM matching parameters -------*" << endl
613  << " | nJet | " << setw(14)
614  << nJetStr.str() << " |" << endl
615  << " | nJetMax | " << setw(14)
616  << nJetMaxStr.str() << " |" << endl
617  << " | Jet algorithm | " << setw(14)
618  << jetStr << " |" << endl
619  << " | eTjetMin | " << setw(14)
620  << eTjetMin << " |" << endl
621  << " | coneRadius | " << setw(14)
622  << coneRadius << " |" << endl
623  << " | etaJetMax | " << setw(14)
624  << etaJetMax << " |" << endl
625  << " | jetAllow | " << setw(14)
626  << jetAllow << " |" << endl
627  << " | jetMatch | " << setw(14)
628  << jetMatch << " |" << endl
629  << " | coneMatchLight | " << setw(14)
630  << coneMatchLight << " |" << endl
631  << " | coneRadiusHeavy | " << setw(14)
632  << coneRadiusHeavy << " |" << endl
633  << " | coneMatchHeavy | " << setw(14)
634  << coneMatchHeavy << " |" << endl
635  << " | Mode | " << setw(14)
636  << modeStr << " |" << endl
637  << " *-----------------------------------------*" << endl;
638 
639  return true;
640 }
641 
642 //--------------------------------------------------------------------------
643 
644 // Step (1): sort the incoming particles
645 
646 inline void JetMatchingAlpgen::sortIncomingProcess(const Event &event) {
647 
648  // Remove resonance decays from original process and keep only final
649  // state. Resonances will have positive status code after this step.
650  omitResonanceDecays(eventProcessOrig, true);
651  eventProcess = workEvent;
652 
653  // Sort original process final state into light/heavy jets and 'other'.
654  // Criteria:
655  // 1 <= ID <= 5 and massless, or ID == 21 --> light jet (typeIdx[0])
656  // 4 <= ID <= 6 and massive --> heavy jet (typeIdx[1])
657  // All else --> other (typeIdx[2])
658  // Note that 'typeIdx' stores indices into 'eventProcess' (after resonance
659  // decays are omitted), while 'typeSet' stores indices into the original
660  // process record, 'eventProcessOrig', but these indices are also valid
661  // in 'event'.
662  for (int i = 0; i < 3; i++) {
663  typeIdx[i].clear();
664  typeSet[i].clear();
665  }
666  for (int i = 0; i < eventProcess.size(); i++) {
667  // Ignore nonfinal and default to 'other'
668  if (!eventProcess[i].isFinal()) continue;
669  int idx = 2;
670 
671  // Light jets
672  if (eventProcess[i].id() == ID_GLUON
673  || (eventProcess[i].idAbs() <= ID_BOT
674  && abs(eventProcess[i].m()) < ZEROTHRESHOLD)) idx = 0;
675 
676  // Heavy jets
677  else if (eventProcess[i].idAbs() >= ID_CHARM
678  && eventProcess[i].idAbs() <= ID_TOP) idx = 1;
679 
680  // Store
681  typeIdx[idx].push_back(i);
682  typeSet[idx].insert(eventProcess[i].daughter1());
683  }
684 
685  // Extract partons from hardest subsystem + ISR + FSR only into
686  // workEvent. Note no resonance showers or MPIs.
687  subEvent(event);
688 }
689 
690 //--------------------------------------------------------------------------
691 
692 // Step (2a): pick which particles to pass to the jet algorithm
693 
694 inline void JetMatchingAlpgen::jetAlgorithmInput(const Event &event,
695  int iType) {
696 
697  // Take input from 'workEvent' and put output in 'workEventJet'
698  workEventJet = workEvent;
699 
700  // Loop over particles and decide what to pass to the jet algorithm
701  for (int i = 0; i < workEventJet.size(); ++i) {
702  if (!workEventJet[i].isFinal()) continue;
703 
704  // jetAllow option to disallow certain particle types
705  if (jetAllow == 1) {
706 
707  // Original AG+Py6 algorithm explicitly excludes tops,
708  // leptons and photons.
709  int id = workEventJet[i].idAbs();
710  if ( (id >= ID_LEPMIN && id <= ID_LEPMAX) || id == ID_TOP
711  || id == ID_PHOTON) {
712  workEventJet[i].statusNeg();
713  continue;
714  }
715  }
716 
717  // Get the index of this particle in original event
718  int idx = workEventJet[i].daughter1();
719 
720  // Start with particle idx, and afterwards track mothers
721  while (true) {
722 
723  // Light jets
724  if (iType == 0) {
725 
726  // Do not include if originates from heavy jet or 'other'
727  if (typeSet[1].find(idx) != typeSet[1].end() ||
728  typeSet[2].find(idx) != typeSet[2].end()) {
729  workEventJet[i].statusNeg();
730  break;
731  }
732 
733  // Made it to start of event record so done
734  if (idx == 0) break;
735  // Otherwise next mother and continue
736  idx = event[idx].mother1();
737 
738  // Heavy jets
739  } else if (iType == 1) {
740 
741  // Only include if originates from heavy jet
742  if (typeSet[1].find(idx) != typeSet[1].end()) break;
743 
744  // Made it to start of event record with no heavy jet mother,
745  // so DO NOT include particle
746  if (idx == 0) {
747  workEventJet[i].statusNeg();
748  break;
749  }
750 
751  // Otherwise next mother and continue
752  idx = event[idx].mother1();
753 
754  // Other jets
755  } else if (iType == 2) {
756 
757  // Only include if originates from other jet
758  if (typeSet[2].find(idx) != typeSet[2].end()) break;
759 
760  // Made it to start of event record with no heavy jet mother,
761  // so DO NOT include particle
762  if (idx == 0) {
763  workEventJet[i].statusNeg();
764  break;
765  }
766 
767  // Otherwise next mother and continue
768  idx = event[idx].mother1();
769 
770  } // if (iType)
771  } // while (true)
772  } // for (i)
773 
774  // For jetMatch = 2, insert ghost particles corresponding to
775  // each hard parton in the original process
776  if (jetMatch == 2) {
777  for (int i = 0; i < int(typeIdx[iType].size()); i++) {
778  // Get y/phi of the parton
779  Vec4 pIn = eventProcess[typeIdx[iType][i]].p();
780  double y = pIn.rap();
781  double phi = pIn.phi();
782 
783  // Create a ghost particle and add to the workEventJet
784  double e = GHOSTENERGY;
785  double e2y = exp(2. * y);
786  double pz = e * (e2y - 1.) / (e2y + 1.);
787  double pt = sqrt(e*e - pz*pz);
788  double px = pt * cos(phi);
789  double py = pt * sin(phi);
790  workEventJet.append( ID_GLUON, 99, 0, 0, 0, 0, 0, 0, px, py, pz, e);
791 
792  // Extra check on reconstructed y/phi values. If many warnings
793  // of this type, GHOSTENERGY may be set too low.
794  if (MATCHINGCHECK) {
795  int lastIdx = workEventJet.size() - 1;
796  if (abs(y - workEventJet[lastIdx].y()) > ZEROTHRESHOLD ||
797  abs(phi - workEventJet[lastIdx].phi()) > ZEROTHRESHOLD)
798  errorMsg("Warning in JetMatchingAlpgen:jetAlgorithmInput: "
799  "ghost particle y/phi mismatch");
800  }
801 
802  } // for (i)
803  } // if (jetMatch == 2)
804 }
805 
806 //--------------------------------------------------------------------------
807 
808 // Step (2b): run jet algorithm and provide common output
809 
810 inline void JetMatchingAlpgen::runJetAlgorithm() {
811 
812  // Run the jet clustering algorithm
813  if (jetAlgorithm == 1)
814  cellJet->analyze(workEventJet, eTjetMin, coneRadius, eTseed);
815  else
816  slowJet->analyze(workEventJet);
817 
818  // Extract four-momenta of jets with |eta| < etaJetMax and
819  // put into jetMomenta. Note that this is done backwards as
820  // jets are removed with SlowJet.
821  jetMomenta.clear();
822  int iJet = (jetAlgorithm == 1) ? cellJet->size() - 1:
823  slowJet->sizeJet() - 1;
824  for (int i = iJet; i > -1; i--) {
825  Vec4 jetMom = (jetAlgorithm == 1) ? cellJet->pMassive(i) :
826  slowJet->p(i);
827  double eta = jetMom.eta();
828 
829  if (abs(eta) > etaJetMax) {
830  if (jetAlgorithm == 2) slowJet->removeJet(i);
831  continue;
832  }
833  jetMomenta.push_back(jetMom);
834  }
835 
836  // Reverse jetMomenta to restore eT/pT ordering
837  reverse(jetMomenta.begin(), jetMomenta.end());
838 }
839 
840 //--------------------------------------------------------------------------
841 
842 // Step (2c): veto decision (returning true vetoes the event)
843 
844 inline bool JetMatchingAlpgen::matchPartonsToJets(int iType) {
845 
846  // Use two different routines for light/heavy jets as
847  // different veto conditions and for clarity
848  if (iType == 0) return (matchPartonsToJetsLight() > 0);
849  else if (iType == 1) return (matchPartonsToJetsHeavy() > 0);
850  else if (iType == 2) return false;
851  return true;
852 }
853 
854 //--------------------------------------------------------------------------
855 
856 // Step(2c): light jets
857 // Return codes are given indicating the reason for a veto.
858 // Although not currently used, they are a useful debugging tool:
859 // 0 = no veto
860 // 1 = veto as number of jets less than number of partons
861 // 2 = veto as exclusive mode and number of jets greater than
862 // number of partons
863 // 3 = veto as inclusive mode and there would be an extra jet
864 // that is harder than any matched soft jet
865 // 4 = veto as there is a parton which does not match a jet
866 
867 inline int JetMatchingAlpgen::matchPartonsToJetsLight() {
868 
869  // Always veto if number of jets is less than original number of jets
870  if (jetMomenta.size() < typeIdx[0].size()) return LESS_JETS;
871  // Veto if in exclusive mode and number of jets bigger than original
872  if (exclusive && jetMomenta.size() > typeIdx[0].size()) return MORE_JETS;
873 
874  // Sort partons by eT/pT
875  sortTypeIdx(typeIdx[0]);
876 
877  // Number of hard partons
878  int nParton = typeIdx[0].size();
879 
880  // Keep track of which jets have been assigned a hard parton
881  vector < bool > jetAssigned;
882  jetAssigned.assign(jetMomenta.size(), false);
883 
884  // Jet matching procedure: (1) deltaR between partons and jets
885  if (jetMatch == 1) {
886 
887  // Loop over light hard partons and get 4-momentum
888  for (int i = 0; i < nParton; i++) {
889  Vec4 p1 = eventProcess[typeIdx[0][i]].p();
890 
891  // Track which jet has the minimal dR measure with this parton
892  int jMin = -1;
893  double dRmin = 0.;
894 
895  // Loop over all jets (skipping those already assigned).
896  for (int j = 0; j < int(jetMomenta.size()); j++) {
897  if (jetAssigned[j]) continue;
898 
899  // DeltaR between parton/jet and store if minimum
900  double dR = (jetAlgorithm == 1)
901  ? REtaPhi(p1, jetMomenta[j]) : RRapPhi(p1, jetMomenta[j]);
902  if (jMin < 0 || dR < dRmin) {
903  dRmin = dR;
904  jMin = j;
905  }
906  } // for (j)
907 
908  // Check for jet-parton match
909  if (jMin >= 0 && dRmin < coneRadius * coneMatchLight) {
910 
911  // If the matched jet is not one of the nParton hardest jets,
912  // the extra left over jet would be harder than some of the
913  // matched jets. This is disallowed, so veto.
914  if (jMin >= nParton) return HARD_JET;
915 
916  // Mark jet as assigned.
917  jetAssigned[jMin] = true;
918 
919  // If no match, then event will be vetoed in all cases
920  } else return UNMATCHED_PARTON;
921 
922  } // for (i)
923 
924  // Jet matching procedure: (2) ghost particles in SlowJet
925  } else {
926 
927  // Loop over added 'ghost' particles and find if assigned to a jet
928  for (int i = workEventJet.size() - nParton;
929  i < workEventJet.size(); i++) {
930  int jMin = slowJet->jetAssignment(i);
931 
932  // Veto if:
933  // 1) not one of nParton hardest jets
934  // 2) not assigned to a jet
935  // 3) jet has already been assigned
936  if (jMin >= nParton) return HARD_JET;
937  if (jMin < 0 || jetAssigned[jMin]) return UNMATCHED_PARTON;
938 
939  // Mark jet as assigned
940  jetAssigned[jMin] = true;
941 
942  } // for (i)
943  } // if (jetMatch)
944 
945  // Minimal eT/pT (CellJet/SlowJet) of matched light jets. Needed
946  // later for heavy jet vetos in inclusive mode.
947  if (nParton > 0)
948  eTpTlightMin = (jetAlgorithm == 1) ? jetMomenta[nParton - 1].eT()
949  : jetMomenta[nParton - 1].pT();
950  else
951  eTpTlightMin = -1.;
952 
953  // No veto
954  return NONE;
955 }
956 
957 //--------------------------------------------------------------------------
958 
959 // Step(2c): heavy jets
960 // Return codes are given indicating the reason for a veto.
961 // Although not currently used, they are a useful debugging tool:
962 // 0 = no veto as there are no extra jets present
963 // 1 = veto as in exclusive mode and extra jets present
964 // 2 = veto as in inclusive mode and extra jets were harder
965 // than any matched light jet
966 
967 inline int JetMatchingAlpgen::matchPartonsToJetsHeavy() {
968 
969  // If there are no extra jets, then accept
970  if (jetMomenta.empty()) return NONE;
971 
972  // Number of hard partons
973  int nParton = typeIdx[1].size();
974 
975  // Remove jets that are close to heavy quarks
976  set < int > removeJets;
977 
978  // Jet matching procedure: (1) deltaR between partons and jets
979  if (jetMatch == 1) {
980 
981  // Loop over heavy hard partons and get 4-momentum
982  for (int i = 0; i < nParton; i++) {
983  Vec4 p1 = eventProcess[typeIdx[1][i]].p();
984 
985  // Loop over all jets, find dR and mark for removal if match
986  for (int j = 0; j < int(jetMomenta.size()); j++) {
987  double dR = (jetAlgorithm == 1) ?
988  REtaPhi(p1, jetMomenta[j]) : RRapPhi(p1, jetMomenta[j]);
989  if (dR < coneRadiusHeavy * coneMatchHeavy)
990  removeJets.insert(j);
991 
992  } // for (j)
993  } // for (i)
994 
995  // Jet matching procedure: (2) ghost particles in SlowJet
996  } else {
997 
998  // Loop over added 'ghost' particles and if assigned to a jet
999  // then mark this jet for removal
1000  for (int i = workEventJet.size() - nParton;
1001  i < workEventJet.size(); i++) {
1002  int jMin = slowJet->jetAssignment(i);
1003  if (jMin >= 0) removeJets.insert(jMin);
1004  }
1005 
1006  }
1007 
1008  // Remove jets (backwards order to not disturb indices)
1009  for (set < int >::reverse_iterator it = removeJets.rbegin();
1010  it != removeJets.rend(); it++)
1011  jetMomenta.erase(jetMomenta.begin() + *it);
1012 
1013  // Handle case if there are still extra jets
1014  if (!jetMomenta.empty()) {
1015 
1016  // Exclusive mode, so immediate veto
1017  if (exclusive) return MORE_JETS;
1018 
1019  // Inclusive mode; extra jets must be softer than any matched light jet
1020  else if (eTpTlightMin >= 0.)
1021  for (size_t j = 0; j < jetMomenta.size(); j++) {
1022  // CellJet uses eT, SlowJet uses pT
1023  if ( (jetAlgorithm == 1 && jetMomenta[j].eT() > eTpTlightMin) ||
1024  (jetAlgorithm == 2 && jetMomenta[j].pT() > eTpTlightMin) )
1025  return HARD_JET;
1026  }
1027 
1028  } // if (!jetMomenta.empty())
1029 
1030  // No extra jets were present so no veto
1031  return NONE;
1032 }
1033 
1034 //==========================================================================
1035 
1036 // Main implementation of Madgraph UserHooks class.
1037 // This may be split out to a separate C++ file if desired,
1038 // but currently included here for ease of use.
1039 
1040 //--------------------------------------------------------------------------
1041 
1042 // Initialisation routine automatically called from Pythia::init().
1043 // Setup all parts needed for the merging.
1044 
1046 
1047  // Initialise values for stored jet matching veto inputs.
1048  pTfirstSave = -1.;
1049  processSubsetSave.init("(eventProcess)", particleDataPtr);
1050  workEventJetSave.init("(workEventJet)", particleDataPtr);
1051 
1052  // Read in Madgraph specific configuration variables
1053  bool setMad = settingsPtr->flag("JetMatching:setMad");
1054 
1055  // If Madgraph parameters are present, then parse in MadgraphPar object
1056  MadgraphPar par;
1057  string parStr = infoPtr->header("MGRunCard");
1058  if (!parStr.empty()) {
1059  par.parse(parStr);
1060  par.printParams();
1061  }
1062 
1063  // Set Madgraph merging parameters from the file if requested
1064  if (setMad) {
1065  if ( par.haveParam("xqcut") && par.haveParam("maxjetflavor")
1066  && par.haveParam("alpsfact") && par.haveParam("ickkw") ) {
1067  settingsPtr->flag("JetMatching:merge", par.getParam("ickkw"));
1068  settingsPtr->parm("JetMatching:qCut", par.getParam("xqcut"));
1069  settingsPtr->mode("JetMatching:nQmatch",
1070  par.getParamAsInt("maxjetflavor"));
1071  settingsPtr->parm("JetMatching:clFact",
1072  clFact = par.getParam("alpsfact"));
1073  if (par.getParamAsInt("ickkw") == 0)
1074  errorMsg("Error in JetMatchingMadgraph:init: "
1075  "Madgraph file parameters are not set for merging");
1076 
1077  // Warn if setMad requested, but one or more parameters not present
1078  } else {
1079  errorMsg("Warning in JetMatchingMadgraph:init: "
1080  "Madgraph merging parameters not found");
1081  if (!par.haveParam("xqcut")) errorMsg("Warning in "
1082  "JetMatchingMadgraph:init: No xqcut");
1083  if (!par.haveParam("ickkw")) errorMsg("Warning in "
1084  "JetMatchingMadgraph:init: No ickkw");
1085  if (!par.haveParam("maxjetflavor")) errorMsg("Warning in "
1086  "JetMatchingMadgraph:init: No maxjetflavor");
1087  if (!par.haveParam("alpsfact")) errorMsg("Warning in "
1088  "JetMatchingMadgraph:init: No alpsfact");
1089  }
1090  }
1091 
1092  // Read in FxFx matching parameters
1093  doFxFx = settingsPtr->flag("JetMatching:doFxFx");
1094  nPartonsNow = settingsPtr->mode("JetMatching:nPartonsNow");
1095  qCutME = settingsPtr->parm("JetMatching:qCutME");
1096  qCutMESq = pow(qCutME,2);
1097 
1098  // Read in Madgraph merging parameters
1099  doMerge = settingsPtr->flag("JetMatching:merge");
1100  doShowerKt = settingsPtr->flag("JetMatching:doShowerKt");
1101  qCut = settingsPtr->parm("JetMatching:qCut");
1102  nQmatch = settingsPtr->mode("JetMatching:nQmatch");
1103  clFact = settingsPtr->parm("JetMatching:clFact");
1104 
1105  // Read in jet algorithm parameters
1106  jetAlgorithm = settingsPtr->mode("JetMatching:jetAlgorithm");
1107  nJetMax = settingsPtr->mode("JetMatching:nJetMax");
1108  eTjetMin = settingsPtr->parm("JetMatching:eTjetMin");
1109  coneRadius = settingsPtr->parm("JetMatching:coneRadius");
1110  etaJetMax = settingsPtr->parm("JetMatching:etaJetMax");
1111  slowJetPower = settingsPtr->mode("JetMatching:slowJetPower");
1112 
1113  // Matching procedure
1114  jetAllow = settingsPtr->mode("JetMatching:jetAllow");
1115  exclusiveMode = settingsPtr->mode("JetMatching:exclusive");
1116  qCutSq = pow(qCut,2);
1117  etaJetMaxAlgo = etaJetMax;
1118 
1119  // Read if veto should be performed internally.
1120  performVeto = settingsPtr->flag("JetMatching:doVeto");
1121 
1122  // If not merging, then done
1123  if (!doMerge) return true;
1124 
1125  // Exclusive mode; if set to 2, then set based on nJet/nJetMax
1126  if (exclusiveMode == 2) {
1127 
1128  // No nJet or nJetMax, so default to exclusive mode
1129  if (nJetMax < 0) {
1130  errorMsg("Warning in JetMatchingMadgraph:init: "
1131  "missing jet multiplicity information; running in exclusive mode");
1132  exclusiveMode = 1;
1133  }
1134  }
1135 
1136  // Initialise chosen jet algorithm.
1137  // Currently, this only supports the kT-algorithm in SlowJet.
1138  // Use the QCD distance measure by default.
1139  jetAlgorithm = 2;
1140  slowJetPower = 1;
1141  slowJet = new SlowJet(slowJetPower, coneRadius, eTjetMin,
1142  etaJetMaxAlgo, 2, 2, nullptr, false);
1143 
1144  // For FxFx, also initialise jet algorithm to define matrix element jets.
1145  // Currently, this only supports the kT-algorithm in SlowJet.
1146  // Use the QCD distance measure by default.
1147  slowJetHard = new SlowJet(slowJetPower, coneRadius, qCutME,
1148  etaJetMaxAlgo, 2, 2, nullptr, false);
1149 
1150  // To access the DJR's
1151  slowJetDJR = new SlowJet(slowJetPower, coneRadius, qCutME,
1152  etaJetMaxAlgo, 2, 2, nullptr, false);
1153 
1154  // A special version of SlowJet to handle heavy and other partons
1155  hjSlowJet = new HJSlowJet(slowJetPower, coneRadius, 0.0,
1156  100.0, 1, 2, nullptr, false, true);
1157 
1158  // Setup local event records
1159  eventProcessOrig.init("(eventProcessOrig)", particleDataPtr);
1160  eventProcess.init("(eventProcess)", particleDataPtr);
1161  workEventJet.init("(workEventJet)", particleDataPtr);
1162 
1163  // Print information
1164  string jetStr = (jetAlgorithm == 1) ? "CellJet" :
1165  (slowJetPower == -1) ? "anti-kT" :
1166  (slowJetPower == 0) ? "C/A" :
1167  (slowJetPower == 1) ? "kT" : "unknown";
1168  string modeStr = (exclusiveMode) ? "exclusive" : "inclusive";
1169  cout << endl
1170  << " *----- Madgraph matching parameters -----*" << endl
1171  << " | qCut | " << setw(14)
1172  << qCut << " |" << endl
1173  << " | nQmatch | " << setw(14)
1174  << nQmatch << " |" << endl
1175  << " | clFact | " << setw(14)
1176  << clFact << " |" << endl
1177  << " | Jet algorithm | " << setw(14)
1178  << jetStr << " |" << endl
1179  << " | eTjetMin | " << setw(14)
1180  << eTjetMin << " |" << endl
1181  << " | etaJetMax | " << setw(14)
1182  << etaJetMax << " |" << endl
1183  << " | jetAllow | " << setw(14)
1184  << jetAllow << " |" << endl
1185  << " | Mode | " << setw(14)
1186  << modeStr << " |" << endl
1187  << " *-----------------------------------------*" << endl;
1188 
1189  return true;
1190 }
1191 
1192 //--------------------------------------------------------------------------
1193 
1194 // Process level vetos
1195 
1197 
1198  eventProcessOrig = process;
1199 
1200  // Setup for veto if hard ME has too many partons.
1201  // This is done to achieve consistency with the Pythia6 implementation.
1202 
1203  // Clear the event of MPI systems and resonace decay products. Store trimmed
1204  // event in workEvent.
1205  sortIncomingProcess(process);
1206 
1207  // Veto in case the hard input matrix element already has too many partons.
1208  if ( !doFxFx && int(typeIdx[0].size()) > nJetMax )
1209  return true;
1210  if ( doFxFx && npNLO() < nJetMax && int(typeIdx[0].size()) > nJetMax )
1211  return true;
1212 
1213  // Done
1214  return false;
1215 
1216 }
1217 
1218 //--------------------------------------------------------------------------
1219 
1220 inline bool JetMatchingMadgraph::doVetoStep(int iPos, int nISR, int nFSR,
1221  const Event& event) {
1222 
1223  // Do not perform any veto if not in the Shower-kT scheme.
1224  if ( !doShowerKt ) return false;
1225 
1226  // Do nothing for emissions after the first one.
1227  if ( nISR + nFSR > 1 ) return false;
1228 
1229  // Do nothing in resonance decay showers.
1230  if (iPos == 5) return false;
1231 
1232  // Clear the event of MPI systems and resonace decay products. Store trimmed
1233  // event in workEvent.
1234  sortIncomingProcess(event);
1235 
1236  // Get (kinematical) pT of first emission
1237  double pTfirst = 0.;
1238 
1239  // Get weak bosons, for later checks if the emission is a "QCD emission".
1240  vector<int> weakBosons;
1241  for (int i = 0; i < event.size(); i++) {
1242  if ( event[i].id() == 22
1243  && event[i].id() == 23
1244  && event[i].idAbs() == 24)
1245  weakBosons.push_back(i);
1246  }
1247 
1248  for (int i = workEvent.size()-1; i > 0; --i) {
1249  if ( workEvent[i].isFinal() && workEvent[i].colType() != 0
1250  && (workEvent[i].statusAbs() == 43 || workEvent[i].statusAbs() == 51)) {
1251  // Check if any of the EW bosons are ancestors of this parton. This
1252  // should never happen for the first non-resonance shower emission.
1253  // Check just to be sure.
1254  bool QCDemission = true;
1255  // Get position of this parton in the actual event (workEvent does
1256  // not contain right mother-daughter relations). Stored in daughters.
1257  int iPosOld = workEvent[i].daughter1();
1258  for (int j = 0; i < int(weakBosons.size()); ++i)
1259  if ( event[iPosOld].isAncestor(j)) {
1260  QCDemission = false;
1261  break;
1262  }
1263  // Done for a QCD emission.
1264  if (QCDemission){
1265  pTfirst = workEvent[i].pT();
1266  break;
1267  }
1268  }
1269  }
1270 
1271  // Store things that are necessary to perform the shower-kT veto externally.
1272  pTfirstSave = pTfirst;
1273  // Done if only inputs for an external vetoing procedure should be stored.
1274  if (!performVeto) return false;
1275 
1276  // Check veto.
1277  if ( doShowerKtVeto(pTfirst) ) return true;
1278 
1279  // No veto if come this far.
1280  return false;
1281 
1282 }
1283 
1284 //--------------------------------------------------------------------------
1285 
1286 inline bool JetMatchingMadgraph::doShowerKtVeto(double pTfirst) {
1287 
1288  // Only check veto in the shower-kT scheme.
1289  if ( !doShowerKt ) return false;
1290 
1291  // Reset veto code
1292  bool doVeto = false;
1293 
1294  // Find the (kinematical) pT of the softest (light) parton in the hard
1295  // process.
1296  int nParton = typeIdx[0].size();
1297  double pTminME=1e10;
1298  for ( int i = 0; i < nParton; ++i)
1299  pTminME = min(pTminME,eventProcess[typeIdx[0][i]].pT());
1300 
1301  // Veto if the softest hard process parton is below Qcut.
1302  if ( nParton > 0 && pow(pTminME,2) < qCutSq ) doVeto = true;
1303 
1304  // For non-highest multiplicity, veto if the hardest emission is harder
1305  // than Qcut.
1306  if ( exclusive && pow(pTfirst,2) > qCutSq ) {
1307  doVeto = true;
1308  // For highest multiplicity sample, veto if the hardest emission is harder
1309  // than the hard process parton.
1310  } else if ( !exclusive && nParton > 0 && pTfirst > pTminME ) {
1311  doVeto = true;
1312  }
1313 
1314  // Return veto
1315  return doVeto;
1316 
1317 }
1318 
1319 //--------------------------------------------------------------------------
1320 
1321 // Function to set the jet clustering scales (to be used as output)
1322 
1323 inline void JetMatchingMadgraph::setDJR( const Event& event) {
1324 
1325  // Clear members.
1326  clearDJR();
1327  vector<double> result;
1328 
1329  // Initialize SlowJetDJR jet algorithm with event
1330  if (!slowJetDJR->setup(event) ) {
1331  errorMsg("Warning in JetMatchingMadgraph:setDJR"
1332  ": the SlowJet algorithm failed on setup");
1333  return;
1334  }
1335 
1336  // Cluster in steps to find all hadronic jets
1337  while ( slowJetDJR->sizeAll() - slowJetDJR->sizeJet() > 0 ) {
1338  // Save the next clustering scale.
1339  result.push_back(sqrt(slowJetDJR->dNext()));
1340  // Perform step.
1341  slowJetDJR->doStep();
1342  }
1343 
1344  // Save clustering scales in reserve order.
1345  for (int i=int(result.size())-1; i >= 0; --i)
1346  DJR.push_back(result[i]);
1347 
1348 }
1349 
1350 //--------------------------------------------------------------------------
1351 
1352 // Function to get the current number of partons in the Born state, as
1353 // read from LHE.
1354 
1356  string npIn = infoPtr->getEventAttribute("npNLO",true);
1357  int np = (npIn != "") ? atoi((char*)npIn.c_str()) : -1;
1358  if ( np < 0 ) { ; }
1359  else return np;
1360  return nPartonsNow;
1361 }
1362 
1363 //--------------------------------------------------------------------------
1364 
1365 // Step (1): sort the incoming particles
1366 
1368 
1369  // Remove resonance decays from original process and keep only final
1370  // state. Resonances will have positive status code after this step.
1371  omitResonanceDecays(eventProcessOrig, true);
1372  clearDJR();
1373  clear_nMEpartons();
1374 
1375  // Note: Compared to older versions (cf. arXiv:2108.07826):
1376  // No more preclustering for FxFx.
1377 
1378  // Note: Compared to older versions (cf. arXiv:2108.07826):
1379  // Set the same eventProcess for both MLM and FxFx. This was done separately
1380  // for FxFx before. Add the type 2 selection also for FxFx
1381 
1382  // For jet matching, simply take hard process state from workEvent.
1383  eventProcess = workEvent;
1384 
1385  // Sort original process final state into light/heavy jets and 'other'.
1386  // Criteria:
1387  // 1 <= ID <= nQmatch, or ID == 21 --> light jet (typeIdx[0])
1388  // nQMatch < ID --> heavy jet (typeIdx[1])
1389  // All else that is colored --> other (typeIdx[2])
1390  // Note that 'typeIdx' stores indices into 'eventProcess' (after resonance
1391  // decays are omitted), while 'typeSet' stores indices into the original
1392  // process record, 'eventProcessOrig', but these indices are also valid
1393  // in 'event'.
1394  for (int i = 0; i < 3; i++) {
1395  typeIdx[i].clear();
1396  typeSet[i].clear();
1397  origTypeIdx[i].clear();
1398  }
1399  for (int i = 0; i < eventProcess.size(); i++) {
1400  // Ignore non-final state and default to 'other'
1401  if (!eventProcess[i].isFinal()) continue;
1402  int idx = -1;
1403  int orig_idx = -1;
1404 
1405  // Light jets: all gluons and quarks with id less than or equal to nQmatch
1406  if (eventProcess[i].isGluon()
1407  || (eventProcess[i].idAbs() <= nQmatch) ) {
1408  orig_idx = 0;
1409  if (doFxFx) {
1410  // Crucial point FxFx: MG5 puts the scale of a not-to-be-matched
1411  // quark 1 MeV lower than scalup. For such particles, we should
1412  // keep the default "2"
1413  idx = ( trunc(1000.*eventProcess[i].scale())
1414  == trunc(1000.*infoPtr->scalup()) ) ? 0 : 2;
1415 
1416  } else {
1417  // Crucial point: MG puts the scale of a non-QCD particle to eCM. For
1418  // such particles, we should keep the default "2"
1419  idx = ( eventProcess[i].scale() < 1.999 * sqrt(infoPtr->eA()
1420  * infoPtr->eB()) ) ? 0 : 2;
1421  }
1422  }
1423 
1424  // Heavy jets: all quarks with id greater than nQmatch
1425  else if (eventProcess[i].idAbs() > nQmatch
1426  && eventProcess[i].idAbs() <= ID_TOP) {
1427  idx = 1;
1428  orig_idx = 1;
1429  // Update to include non-SM colored particles
1430  } else if (eventProcess[i].colType() != 0
1431  && eventProcess[i].idAbs() > ID_TOP) {
1432  idx = 1;
1433  orig_idx = 1;
1434  }
1435  if( idx < 0 ) continue;
1436  // Store
1437  typeIdx[idx].push_back(i);
1438  typeSet[idx].insert(eventProcess[i].daughter1());
1439  origTypeIdx[orig_idx].push_back(i);
1440  }
1441 
1442  // Exclusive mode; if set to 2, then set based on nJet/nJetMax
1443  if (exclusiveMode == 2) {
1444 
1445  // Inclusive if nJet == nJetMax, exclusive otherwise
1446  int nParton = origTypeIdx[0].size();
1447  exclusive = (nParton == nJetMax) ? false : true;
1448 
1449  // Otherwise, just set as given
1450  } else {
1451  exclusive = (exclusiveMode == 0) ? false : true;
1452  }
1453 
1454  // Extract partons from hardest subsystem + ISR + FSR only into
1455  // workEvent. Note no resonance showers or MPIs.
1456  subEvent(event);
1457 
1458  // Store things that are necessary to perform the kT-MLM veto externally.
1459  int nParton = typeIdx[0].size();
1460  processSubsetSave.clear();
1461  for ( int i = 0; i < nParton; ++i)
1462  processSubsetSave.append( eventProcess[typeIdx[0][i]] );
1463 
1464 }
1465 
1466 //--------------------------------------------------------------------------
1467 
1468 // Step (2a): pick which particles to pass to the jet algorithm
1469 
1471  int iType) {
1472 
1473  // Take input from 'workEvent' and put output in 'workEventJet'
1474  workEventJet = workEvent;
1475 
1476  // Loop over particles and decide what to pass to the jet algorithm
1477  for (int i = 0; i < workEventJet.size(); ++i) {
1478  if (!workEventJet[i].isFinal()) continue;
1479 
1480  // jetAllow option to disallow certain particle types
1481  if (jetAllow == 1) {
1482  // Remove all non-QCD partons from veto list
1483  if( workEventJet[i].colType() == 0 ) {
1484  workEventJet[i].statusNeg();
1485  continue;
1486  }
1487  }
1488 
1489  // Get the index of this particle in original event
1490  int idx = workEventJet[i].daughter1();
1491 
1492  // Start with particle idx, and afterwards track mothers
1493  while (true) {
1494 
1495  // Light jets
1496  if (iType == 0) {
1497 
1498  // Do not include if originates from heavy jet or 'other'
1499  if (typeSet[1].find(idx) != typeSet[1].end() ||
1500  typeSet[2].find(idx) != typeSet[2].end()) {
1501  workEventJet[i].statusNeg();
1502  break;
1503  }
1504 
1505  // Made it to start of event record so done
1506  if (idx == 0) break;
1507  // Otherwise next mother and continue
1508  idx = event[idx].mother1();
1509 
1510  // Heavy jets
1511  } else if (iType == 1) {
1512 
1513  // Only include if originates from heavy jet
1514  if (typeSet[1].find(idx) != typeSet[1].end()) break;
1515 
1516  // Made it to start of event record with no heavy jet mother,
1517  // so DO NOT include particle
1518  if (idx == 0) {
1519  workEventJet[i].statusNeg();
1520  break;
1521  }
1522 
1523  // Otherwise next mother and continue
1524  idx = event[idx].mother1();
1525 
1526  // Other jets
1527  } else if (iType == 2) {
1528 
1529  // Only include if originates from other jet
1530  if (typeSet[2].find(idx) != typeSet[2].end()) break;
1531 
1532  // Made it to start of event record with no heavy jet mother,
1533  // so DO NOT include particle
1534  if (idx == 0) {
1535  workEventJet[i].statusNeg();
1536  break;
1537  }
1538 
1539  // Otherwise next mother and continue
1540  idx = event[idx].mother1();
1541 
1542  } // if (iType)
1543  } // while (true)
1544  } // for (i)
1545 }
1546 
1547 //--------------------------------------------------------------------------
1548 
1549 // Step (2b): run jet algorithm and provide common output
1550 // This does nothing, because the jet algorithm is run several times
1551 // in the matching algorithm.
1552 
1554 
1555 //--------------------------------------------------------------------------
1556 
1557 // Step (2c): veto decision (returning true vetoes the event)
1558 
1560 
1561  // Use different routines for light/heavy/other jets as
1562  // different veto conditions and for clarity
1563  if (iType == 0) {
1564  // Record the jet separations here, also if matchPartonsToJetsLight
1565  // returns preemptively.
1566  setDJR(workEventJet);
1567  set_nMEpartons(origTypeIdx[0].size(), typeIdx[0].size());
1568  // Perform jet matching.
1569  return (matchPartonsToJetsLight() > 0);
1570  } else if (iType == 1) {
1571  return (matchPartonsToJetsHeavy() > 0);
1572  } else {
1573  return (matchPartonsToJetsOther() > 0);
1574  }
1575 
1576 }
1577 
1578 //--------------------------------------------------------------------------
1579 
1580 // Step(2c): light jets
1581 // Return codes are given indicating the reason for a veto.
1582 // Although not currently used, they are a useful debugging tool:
1583 // 0 = no veto
1584 // 1 = veto as number of jets less than number of partons
1585 // 2 = veto as exclusive mode and number of jets greater than
1586 // number of partons
1587 // 3 = veto as inclusive mode and there would be an extra jet
1588 // that is harder than any matched soft jet
1589 // 4 = veto as there is a parton which does not match a jet
1590 
1592 
1593  // Store things that are necessary to perform the kT-MLM veto externally.
1594  workEventJetSave = workEventJet;
1595  // Done if only inputs for an external vetoing procedure should be stored.
1596  if (!performVeto) return false;
1597 
1598  // Count the number of hard partons
1599  int nParton = typeIdx[0].size();
1600 
1601  // Initialize SlowJet with current working event
1602  if (!slowJet->setup(workEventJet) ) {
1603  errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
1604  "Light: the SlowJet algorithm failed on setup");
1605  return NONE;
1606  }
1607  double localQcutSq = qCutSq;
1608  double dOld = 0.0;
1609  // Cluster in steps to find all hadronic jets at the scale qCut
1610  while ( slowJet->sizeAll() - slowJet->sizeJet() > 0 ) {
1611  if( slowJet->dNext() > localQcutSq ) break;
1612  dOld = slowJet->dNext();
1613  slowJet->doStep();
1614  }
1615  int nJets = slowJet->sizeJet();
1616  int nClus = slowJet->sizeAll();
1617 
1618  // Debug printout.
1619  if (MATCHINGDEBUG) slowJet->list(true);
1620 
1621  // Count of the number of hadronic jets in SlowJet accounting
1622  int nCLjets = nClus - nJets;
1623  // Get number of partons. Different for MLM and FxFx schemes.
1624  // Note: Difference compared to old versions of FxFx (cf. arXiv:2108.07826):
1625  // For FxFx, change nRequested subtracting the typeIdx[2] partons
1626  // Exclude the highest multiplicity sample in the case of real emissions,
1627  // exclude all typeIdx[2] particles (npNLO=multiplicity,
1628  // nJetMax=njmax(shower_card), typeIdx[2]="Weak" jets)
1629  int nRequested = ( doFxFx
1630  && !(npNLO()==nJetMax
1631  && npNLO()==((int)typeIdx[2].size()-1) ))
1632  ? npNLO()-typeIdx[2].size()
1633  : nParton;
1634 
1635  // Note: Difference compared to old versions of FxFx (cf. arXiv:2108.07826):
1636  // For FxFx veto the real emissions that have only typeIdx=2 partons
1637  // Exclude the highest multiplicity sample
1638  if ( doFxFx && npNLO()<nJetMax && typeIdx[2].size()>0
1639  && npNLO()==((int)typeIdx[2].size()-1)) {
1640  return MORE_JETS;
1641  }
1642 
1643  // Veto event if too few hadronic jets
1644  if ( nCLjets < nRequested ) return LESS_JETS;
1645 
1646  // In exclusive mode, do not allow more hadronic jets than partons
1647  if ( exclusive && !doFxFx ) {
1648  if ( nCLjets > nRequested ) return MORE_JETS;
1649  } else {
1650 
1651  // For FxFx, in the non-highest multipicity, all jets need to matched to
1652  // partons. For nCLjets > nRequested, this is not possible. Hence, we can
1653  // veto here already.
1654  // Note: Difference wrt. old versions of FxFx (cf. arXiv:2108.07826):
1655  // Change the nRequested to npNLO() in the first condition. This way we
1656  // correctly select the non-highest multipicity regardless the nRequested
1657  if ( doFxFx && npNLO() < nJetMax && nCLjets > nRequested )
1658  return MORE_JETS;
1659 
1660  // Now continue in inclusive mode.
1661  // In inclusive mode, there can be more hadronic jets than partons,
1662  // provided that all partons are properly matched to hadronic jets.
1663  // Start by setting up the jet algorithm.
1664  if (!slowJet->setup(workEventJet) ) {
1665  errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
1666  "Light: the SlowJet algorithm failed on setup");
1667  return NONE;
1668  }
1669 
1670  // For FxFx, continue clustering as long as the jet separation is above
1671  // qCut.
1672  if (doFxFx) {
1673  while ( slowJet->sizeAll() - slowJet->sizeJet() > 0 ) {
1674  if( slowJet->dNext() > localQcutSq ) break;
1675  slowJet->doStep();
1676  }
1677  // For MLM, cluster into hadronic jets until there are the same number as
1678  // partons.
1679  } else {
1680  while ( slowJet->sizeAll() - slowJet->sizeJet() > nParton )
1681  slowJet->doStep();
1682  }
1683 
1684  // Sort partons in pT. Update local qCut value.
1685  // Hadronic jets are already sorted in pT.
1686  localQcutSq = dOld;
1687  if ( clFact >= 0. && nParton > 0 ) {
1688  vector<double> partonPt;
1689  for (int i = 0; i < nParton; ++i)
1690  partonPt.push_back( eventProcess[typeIdx[0][i]].pT2() );
1691  sort( partonPt.begin(), partonPt.end());
1692  localQcutSq = max( qCutSq, partonPt[0]);
1693  }
1694  nJets = slowJet->sizeJet();
1695  nClus = slowJet->sizeAll();
1696  }
1697  // Update scale if clustering factor is non-zero
1698  if ( clFact != 0. ) localQcutSq *= pow2(clFact);
1699 
1700  Event tempEvent;
1701  tempEvent.init( "(tempEvent)", particleDataPtr);
1702  int nPass = 0;
1703  double pTminEstimate = -1.;
1704  // Construct a master copy of the event containing only the
1705  // hardest nParton hadronic clusters. While constructing the event,
1706  // the parton type (ID_GLUON) and status (98,99) are arbitrary.
1707  for (int i = nJets; i < nClus; ++i) {
1708  tempEvent.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0, slowJet->p(i).px(),
1709  slowJet->p(i).py(), slowJet->p(i).pz(), slowJet->p(i).e() );
1710  ++nPass;
1711  pTminEstimate = max( pTminEstimate, slowJet->pT(i));
1712  if(nPass == nRequested) break;
1713  }
1714 
1715  int tempSize = tempEvent.size();
1716  // This keeps track of which hadronic jets are matched to parton
1717  vector<bool> jetAssigned;
1718  jetAssigned.assign( tempSize, false);
1719 
1720  // This keeps track of which partons are matched to which hadronic
1721  // jets.
1722  vector< vector<bool> > partonMatchesJet;
1723  for (int i=0; i < nParton; ++i )
1724  partonMatchesJet.push_back( vector<bool>(tempEvent.size(),false) );
1725 
1726  // Begin matching.
1727  // Do jet matching for FxFx.
1728  // Make sure that the nPartonsNow hardest hadronic jets are matched to any
1729  // of the nPartonsNow (+1) partons. This matching is done by attaching a jet
1730  // from the list of unmatched hadronic jets, and appending a jet from the
1731  // list of partonic jets, one at a time. The partonic jet will be clustered
1732  // with the hadronic jet or the beam if the distance measure is below the
1733  // cut. The hadronic jet is matched once this happens. Otherwise, another
1734  // partonic jet is tried. When a hadronic jet is matched to a partonic jet,
1735  // it is removed from the list of unmatched hadronic jets. This process
1736  // continues until the nPartonsNow hardest hadronic jets are matched to
1737  // partonic jets, or it is not possible to make a match for a hadronic jet.
1738  int iNow = 0;
1739  int nMatched = 0;
1740  while ( doFxFx && iNow < tempSize ) {
1741 
1742  // Check if this shower jet matches any partonic jet.
1743  Event tempEventJet;
1744  tempEventJet.init("(tempEventJet)", particleDataPtr);
1745  for (int i=0; i < nParton; ++i ) {
1746 
1747  // Only assign a parton once.
1748  //for (int j=0; j < tempSize; ++j )
1749  // if ( partonMatchesJet[i][j]) continue;
1750 
1751  // Attach a single hadronic jet.
1752  tempEventJet.clear();
1753  tempEventJet.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0,
1754  tempEvent[iNow].px(), tempEvent[iNow].py(),
1755  tempEvent[iNow].pz(), tempEvent[iNow].e() );
1756  // Attach the current parton.
1757  Vec4 pIn = eventProcess[typeIdx[0][i]].p();
1758  tempEventJet.append( ID_GLUON, 99, 0, 0, 0, 0, 0, 0,
1759  pIn.px(), pIn.py(), pIn.pz(), pIn.e() );
1760 
1761  // Setup jet algorithm.
1762  if ( !slowJet->setup(tempEventJet) ) {
1763  errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
1764  "Light: the SlowJet algorithm failed on setup");
1765  return NONE;
1766  }
1767 
1768  // These are the conditions for the hadronic jet to match the parton
1769  // at the local qCut scale
1770  if ( slowJet->iNext() == tempEventJet.size() - 1
1771  && slowJet->jNext() > -1 && slowJet->dNext() < localQcutSq ) {
1772  jetAssigned[iNow] = true;
1773  partonMatchesJet[i][iNow] = true;
1774  }
1775 
1776  } // End loop over hard partons.
1777 
1778  // Veto if the jet could not be assigned to any parton.
1779  if ( jetAssigned[iNow] ) nMatched++;
1780 
1781  // Continue;
1782  ++iNow;
1783  }
1784 
1785  // Jet matching veto for FxFx
1786  if (doFxFx) {
1787  // Note: Difference wrt. old versions of FxFx (cf. arXiv:2108.07826):
1788  // Change the nRequested to npNLO() in the first condition (cf. above).
1789  // This way we correctly select the non-highest multipicity regardless
1790  // the nRequested
1791  if ( npNLO() < nJetMax && nMatched != nRequested )
1792  return UNMATCHED_PARTON;
1793  if ( npNLO() == nJetMax && nMatched < nRequested )
1794  return UNMATCHED_PARTON;
1795  }
1796 
1797  // Do jet matching for MLM.
1798  // Take the list of unmatched hadronic jets and append a parton, one at
1799  // a time. The parton will be clustered with the "closest" hadronic jet
1800  // or the beam if the distance measure is below the cut. When a hadronic
1801  // jet is matched to a parton, it is removed from the list of unmatched
1802  // hadronic jets. This process continues until all hadronic jets are
1803  // matched to partons or it is not possible to make a match.
1804  iNow = 0;
1805  while (!doFxFx && iNow < nParton ) {
1806  Event tempEventJet;
1807  tempEventJet.init("(tempEventJet)", particleDataPtr);
1808  for (int i = 0; i < tempSize; ++i) {
1809  if (jetAssigned[i]) continue;
1810  Vec4 pIn = tempEvent[i].p();
1811  // Append unmatched hadronic jets
1812  tempEventJet.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0,
1813  pIn.px(), pIn.py(), pIn.pz(), pIn.e() );
1814  }
1815 
1816  Vec4 pIn = eventProcess[typeIdx[0][iNow]].p();
1817  // Append the current parton
1818  tempEventJet.append( ID_GLUON, 99, 0, 0, 0, 0, 0, 0,
1819  pIn.px(), pIn.py(), pIn.pz(), pIn.e() );
1820  if ( !slowJet->setup(tempEventJet) ) {
1821  errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
1822  "Light: the SlowJet algorithm failed on setup");
1823  return NONE;
1824  }
1825  // These are the conditions for the hadronic jet to match the parton
1826  // at the local qCut scale
1827  if ( slowJet->iNext() == tempEventJet.size() - 1
1828  && slowJet->jNext() > -1 && slowJet->dNext() < localQcutSq ) {
1829  int iKnt = -1;
1830  for (int i = 0; i != tempSize; ++i) {
1831  if (jetAssigned[i]) continue;
1832  ++iKnt;
1833  // Identify the hadronic jet that matches the parton
1834  if (iKnt == slowJet->jNext() ) jetAssigned[i] = true;
1835  }
1836  } else {
1837  return UNMATCHED_PARTON;
1838  }
1839  ++iNow;
1840  }
1841 
1842  // Minimal eT/pT (CellJet/SlowJet) of matched light jets.
1843  // Needed later for heavy jet vetos in inclusive mode.
1844  // This information is not used currently.
1845  if (nParton > 0 && pTminEstimate > 0) eTpTlightMin = pTminEstimate;
1846  else eTpTlightMin = -1.;
1847 
1848  // Record the jet separations.
1849  setDJR(workEventJet);
1850 
1851  // No veto
1852  return NONE;
1853 }
1854 
1855 //--------------------------------------------------------------------------
1856 
1857 // Step(2c): heavy jets
1858 // Return codes are given indicating the reason for a veto.
1859 // Although not currently used, they are a useful debugging tool:
1860 // 0 = no veto as there are no extra jets present
1861 // 1 = veto as in exclusive mode and extra jets present
1862 // 2 = veto as in inclusive mode and extra jets were harder
1863 // than any matched light jet
1864 
1866 
1867  // Currently, heavy jets are unmatched
1868  // If there are no extra jets, then accept
1869  // jetMomenta is NEVER used by MadGraph and is always empty.
1870  // This check does nothing.
1871  // Rather, if there is any heavy flavor that is harder than
1872  // what is present at the LHE level, then the event should
1873  // be vetoed.
1874 
1875  // if (jetMomenta.empty()) return NONE;
1876  // Count the number of hard partons
1877  int nParton = typeIdx[1].size();
1878 
1879  Event tempEventJet(workEventJet);
1880 
1881  double scaleF(1.0);
1882  // Rescale the heavy partons that are from the hard process to
1883  // have pT=collider energy. Soft/collinear gluons will cluster
1884  // onto them, leaving a remnant of hard emissions.
1885  for( int i=0; i<nParton; ++i) {
1886  scaleF = eventProcessOrig[0].e()/workEventJet[typeIdx[1][i]].pT();
1887  tempEventJet[typeIdx[1][i]].rescale5(scaleF);
1888  }
1889 
1890  if (!hjSlowJet->setup(tempEventJet) ) {
1891  errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
1892  "Heavy: the SlowJet algorithm failed on setup");
1893  return NONE;
1894  }
1895 
1896 
1897  while ( hjSlowJet->sizeAll() - hjSlowJet->sizeJet() > 0 ) {
1898  if( hjSlowJet->dNext() > qCutSq ) break;
1899  hjSlowJet->doStep();
1900  }
1901 
1902  int nCLjets(0);
1903  // Count the number of clusters with pT>qCut. This includes the
1904  // original hard partons plus any hard emissions.
1905  for(int idx=0 ; idx< hjSlowJet->sizeAll(); ++idx) {
1906  if( hjSlowJet->pT(idx) > sqrt(qCutSq) ) nCLjets++;
1907  }
1908 
1909  // Debug printout.
1910  if (MATCHINGDEBUG) hjSlowJet->list(true);
1911 
1912  // Count of the number of hadronic jets in SlowJet accounting
1913  // int nCLjets = nClus - nJets;
1914  // Get number of partons. Different for MLM and FxFx schemes.
1915  int nRequested = nParton;
1916 
1917  // Veto event if too few hadronic jets
1918  if ( nCLjets < nRequested ) {
1919  if (MATCHINGDEBUG) cout << "veto : hvy LESS_JETS " << endl;
1920  if (MATCHINGDEBUG) cout << "nCLjets = " << nCLjets << "; nRequest = "
1921  << nRequested << endl;
1922  return LESS_JETS;
1923  }
1924 
1925  // In exclusive mode, do not allow more hadronic jets than partons
1926  if ( exclusive ) {
1927  if ( nCLjets > nRequested ) {
1928  if (MATCHINGDEBUG) cout << "veto : excl hvy MORE_JETS " << endl;
1929  return MORE_JETS;
1930  }
1931  }
1932 
1933  // No extra jets were present so no veto
1934  return NONE;
1935 }
1936 
1937 //--------------------------------------------------------------------------
1938 
1939 // Step(2c): other jets
1940 // Return codes are given indicating the reason for a veto.
1941 // Although not currently used, they are a useful debugging tool:
1942 // 0 = no veto as there are no extra jets present
1943 // 1 = veto as in exclusive mode and extra jets present
1944 // 2 = veto as in inclusive mode and extra jets were harder
1945 // than any matched light jet
1946 
1948 
1949  // Currently, heavy jets are unmatched
1950  // If there are no extra jets, then accept
1951  // jetMomenta is NEVER used by MadGraph and is always empty.
1952  // This check does nothing.
1953  // Rather, if there is any heavy flavor that is harder than
1954  // what is present at the LHE level, then the event should
1955  // be vetoed.
1956 
1957  // if (jetMomenta.empty()) return NONE;
1958  // Count the number of hard partons
1959  int nParton = typeIdx[2].size();
1960 
1961  Event tempEventJet(workEventJet);
1962 
1963  double scaleF(1.0);
1964  // Rescale the heavy partons that are from the hard process to
1965  // have pT=collider energy. Soft/collinear gluons will cluster
1966  // onto them, leaving a remnant of hard emissions.
1967  for( int i=0; i<nParton; ++i) {
1968  scaleF = eventProcessOrig[0].e()/workEventJet[typeIdx[2][i]].pT();
1969  tempEventJet[typeIdx[2][i]].rescale5(scaleF);
1970  }
1971 
1972  if (!hjSlowJet->setup(tempEventJet) ) {
1973  errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
1974  "Heavy: the SlowJet algorithm failed on setup");
1975  return NONE;
1976  }
1977 
1978 
1979  while ( hjSlowJet->sizeAll() - hjSlowJet->sizeJet() > 0 ) {
1980  if( hjSlowJet->dNext() > qCutSq ) break;
1981  hjSlowJet->doStep();
1982  }
1983 
1984  int nCLjets(0);
1985  // Count the number of clusters with pT>qCut. This includes the
1986  // original hard partons plus any hard emissions.
1987  for(int idx=0 ; idx< hjSlowJet->sizeAll(); ++idx) {
1988  if( hjSlowJet->pT(idx) > sqrt(qCutSq) ) nCLjets++;
1989  }
1990 
1991  // Debug printout.
1992  if (MATCHINGDEBUG) hjSlowJet->list(true);
1993 
1994  // Count of the number of hadronic jets in SlowJet accounting
1995  // int nCLjets = nClus - nJets;
1996  // Get number of partons. Different for MLM and FxFx schemes.
1997  int nRequested = nParton;
1998 
1999  // Veto event if too few hadronic jets
2000  if ( nCLjets < nRequested ) {
2001  if (MATCHINGDEBUG) cout << "veto : other LESS_JETS " << endl;
2002  if (MATCHINGDEBUG) cout << "nCLjets = " << nCLjets << "; nRequest = "
2003  << nRequested << endl;
2004  return LESS_JETS;
2005  }
2006 
2007  // In exclusive mode, do not allow more hadronic jets than partons
2008  if ( exclusive ) {
2009  if ( nCLjets > nRequested ) {
2010  if (MATCHINGDEBUG) cout << "veto : excl other MORE_JETS" << endl;
2011  return MORE_JETS;
2012  }
2013  }
2014 
2015  // No extra jets were present so no veto
2016  return NONE;
2017 }
2018 
2019 //==========================================================================
2020 
2021 } // end namespace Pythia8
2022 
2023 #endif // end Pythia8_JetMatching_H
bool initAfterBeams()
Initialisation.
Definition: JetMatching.h:518
void sortIncomingProcess(const Event &)
Different steps of the matching algorithm.
Definition: JetMatching.h:1367
Definition: Analysis.h:307
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:182
bool doVetoProcessLevel(Event &process)
Definition: JetMatching.h:137
JetMatchingAlpgen()
Constructor and destructor.
Definition: JetMatching.h:246
int nJetMax
Maximum and current number of jets.
Definition: JetMatching.h:189
static const bool MATCHINGDEBUG
Constants to be changed for debug printout or extra checks.
Definition: JetMatching.h:175
void runJetAlgorithm()
Definition: JetMatching.h:1553
int matchPartonsToJetsLight()
Definition: JetMatching.h:1591
The Event class holds all info on the generated event.
Definition: Event.h:453
vector< double > getDJR()
Definition: JetMatching.h:299
void jetAlgorithmInput(const Event &, int)
Step (2a): pick which particles to pass to the jet algorithm.
Definition: JetMatching.h:1470
Event getWorkEventJet()
Definition: JetMatching.h:306
Declaration of main UserHooks class to perform Alpgen matching.
Definition: JetMatching.h:241
double RRapPhi(const Vec4 &v1, const Vec4 &v2)
R is distance in cylindrical (y/eta, phi) coordinates.
Definition: Basics.cc:753
void errorMsg(string messageIn)
Print a message the first few times. Insert in database.
Definition: JetMatching.h:164
bool doVetoStep(int, int, int, const Event &)
Definition: JetMatching.h:1220
void clear()
Clear event record.
Definition: Event.h:473
void findNext()
Find next cluster pair to join.
Definition: JetMatching.h:52
bool doVetoStep(int, int, int, const Event &)
Definition: JetMatching.h:149
void setDJR(const Event &event)
Function to set the jet clustering scales (to be used as output)
Definition: JetMatching.h:1323
Event eventProcessOrig
Definition: JetMatching.h:208
bool canVetoPartonLevelEarly()
Parton level vetos (before beam remnants and resonance decays)
Definition: JetMatching.h:143
int slowJetPower
SlowJet specific.
Definition: JetMatching.h:202
void clear_nMEpartons()
Functions to clear and set the jet clustering scales.
Definition: JetMatching.h:329
map< string, int > messages
Map for all error messages.
Definition: JetMatching.h:231
~JetMatching()
Definition: JetMatching.h:92
SlowJet(int powerIn, double Rin, double pTjetMinIn=0., double etaMaxIn=25., int selectIn=2, int massSetIn=2, SlowJetHook *sjHookPtrIn=0, bool useFJcoreIn=true, bool useStandardRin=true)
Constructor.
Definition: Analysis.h:427
int append(Particle entryIn)
Put a new particle at the end of the event record; return index.
Definition: Event.h:507
int numberVetoStep()
Shower step vetoes (after the first emission, for Shower-kT scheme)
Definition: JetMatching.h:147
double REtaPhi(const Vec4 &v1, const Vec4 &v2)
Distance in cylindrical (eta, phi) coordinates.
Definition: Basics.cc:764
static const double TINY
Small number to avoid division by zero.
Definition: Analysis.h:524
int nEta
CellJet specific.
Definition: JetMatching.h:219
vector< Vec4 > jetMomenta
Definition: JetMatching.h:216
JetMatching()
Constructor and destructor.
Definition: JetMatching.h:90
Definition: Analysis.h:422
bool parse(const string paramStr)
Parse an incoming Madgraph parameter file string.
Definition: GeneratorInput.h:1139
Declaration of main UserHooks class to perform Madgraph matching.
Definition: JetMatching.h:274
CellJet * cellJet
Internal jet algorithms.
Definition: JetMatching.h:196
static const int TIMESTOPRINT
Constants: could only be changed in the code itself.
Definition: Analysis.h:523
Definition: JetMatching.h:29
int jetAllow
Merging procedure parameters.
Definition: JetMatching.h:223
double eTpTlightMin
Store the minimum eT/pT of matched light jets.
Definition: JetMatching.h:228
Definition: JetMatching.h:85
void printParams()
Print parameters read from the &#39;.par&#39; file.
Definition: GeneratorInput.h:1190
int size() const
Event record size.
Definition: Event.h:504
int matchPartonsToJetsHeavy()
Definition: JetMatching.h:1865
bool doVetoPartonLevelEarly(const Event &event)
Early parton level veto (before beam remnants and resonance showers)
Definition: JetMatching.h:384
int jetAlgorithm
Jet algorithm parameters.
Definition: JetMatching.h:192
bool canVetoStep()
Definition: JetMatching.h:291
bool canVetoStep()
Definition: JetMatching.h:148
bool canVetoProcessLevel()
Process level vetos.
Definition: JetMatching.h:136
bool doShowerKtVeto(double pTfirst)
Definition: JetMatching.h:1286
SlowJet * slowJetDJR
Definition: JetMatching.h:296
int matchPartonsToJetsOther()
Definition: JetMatching.h:1947
double phi(const Vec4 &v1, const Vec4 &v2)
phi is azimuthal angle between v1 and v2 around z axis.
Definition: Basics.cc:693
int numberVetoStep()
Shower step vetoes (after the first emission, for Shower-kT scheme)
Definition: JetMatching.h:290
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
double m(const Vec4 &v1)
Invariant mass and its square.
Definition: Basics.cc:595
UserHooks is base class for user access to program execution.
Definition: UserHooks.h:32
bool doShowerKt
Definition: JetMatching.h:186
JetMatchingMadgraph()
Constructor and destructor.
Definition: JetMatching.h:279
bool haveParam(const string &paramIn)
Check if a parameter exists.
Definition: GeneratorInput.h:170
Definition: GeneratorInput.h:156
bool doVetoProcessLevel(Event &process)
Process level vetos.
Definition: JetMatching.h:1196
double getParam(const string &paramIn)
Definition: GeneratorInput.h:175
void init(string headerIn="", ParticleData *particleDataPtrIn=0, int startColTagIn=100)
Initialize header for event listing, particle data table, and colour.
Definition: Event.h:467
bool matchPartonsToJets(int)
Step (2c): veto decision (returning true vetoes the event)
Definition: JetMatching.h:1559
Definition: Analysis.h:371
bool canVetoProcessLevel()
Process level vetos.
Definition: JetMatching.h:286
Definition: Basics.h:32
bool doMerge
Master switch for merging.
Definition: JetMatching.h:183
bool initAfterBeams()
Initialisation.
Definition: JetMatching.h:1045
void clearDJR()
Functions to clear and set the jet clustering scales.
Definition: JetMatching.h:326
int npNLO()
Definition: JetMatching.h:1355