Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / examples / MLMhooks.h
1 // MLMhooks.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2012 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6 // Author: Richard Corke (richard.corke@thep.lu.se)
7 // This file provides the MLMhooks class to perform MLM merging.
8 // Example usage is shown in main32.cc, and further details
9 // can be found in the 'Alpgen and MLM Merging' manual page.
10
11 #ifndef MLMHOOKS_H
12 #define MLMHOOKS_H
13
14 // Includes
15 #include "Pythia.h"
16 using namespace Pythia8;
17
18 //==========================================================================
19
20 // Preprocessor settings
21
22 // Debug flag to give verbose information at all stages of the merging.
23 //#define MLMDEBUG
24 // Debug flag to enable some extra checks in the code
25 //#define MLMCHECK
26
27 //==========================================================================
28
29 // Declaration of main MLMhooks class to perform MLM matching.
30 // Note that it is defined with virtual inheritance, so that it can
31 // be combined with other UserHooks classes, see e.g. main32.cc.
32
33 class MLMhooks : virtual public UserHooks {
34
35 public:
36
37   // Constructor and destructor
38   MLMhooks() : cellJet(NULL), slowJet(NULL) {}
39   ~MLMhooks() {
40     if (cellJet) delete cellJet;
41     if (slowJet) delete slowJet; 
42   }
43
44   // Initialisation
45   virtual bool initAfterBeams();
46
47   // Process level vetos
48   virtual bool canVetoProcessLevel();
49   virtual bool doVetoProcessLevel(Event &);
50
51   // Parton level vetos (before beam remnants and resonance decays)
52   virtual bool canVetoPartonLevelEarly();
53   virtual bool doVetoPartonLevelEarly(const Event &);
54
55 private:
56
57   // Different steps of the MLM matching algorithm
58   void sortIncomingProcess(const Event &);
59   void jetAlgorithmInput(const Event &, int);
60   void runJetAlgorithm();
61   bool matchPartonsToJets(int);
62   int  matchPartonsToJetsLight();
63   int  matchPartonsToJetsHeavy();
64
65   // DeltaR between two 4-vectors (eta and y variants)
66   inline double Vec4eta(const Vec4 &pIn) {
67     return -log(tan(pIn.theta() / 2.));
68   }
69   inline double Vec4y(const Vec4 &pIn) {
70     return 0.5 * log((pIn.e() + pIn.pz()) / (pIn.e() - pIn.pz()));
71   }
72   inline double deltaReta(const Vec4 &p1, const Vec4 &p2) {
73     double dEta = abs(Vec4eta(p1) - Vec4eta(p2));
74     double dPhi = abs(p1.phi() - p2.phi());
75     if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
76     return sqrt(dEta*dEta + dPhi*dPhi);
77   }
78   inline double deltaRy(const Vec4 &p1, const Vec4 &p2) {
79     double dy   = abs(Vec4y(p1) - Vec4y(p2));
80     double dPhi = abs(p1.phi() - p2.phi());
81     if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
82     return sqrt(dy*dy + dPhi*dPhi);
83   }
84
85   // Function to sort typeIdx vectors into descending eT/pT order.
86   // Uses a selection sort, as number of partons generally small
87   // and so efficiency not a worry.
88   void sortTypeIdx(vector < int > &vecIn) {
89     for (size_t i = 0; i < vecIn.size(); i++) {
90       size_t jMax = i;
91       double vMax = (jetAlgorithm == 1) ?
92                     eventProcess[vecIn[i]].eT() :
93                     eventProcess[vecIn[i]].pT();
94       for (size_t j = i + 1; j < vecIn.size(); j++) {
95         double vNow = (jetAlgorithm == 1) ?
96                       eventProcess[vecIn[j]].eT() :
97                       eventProcess[vecIn[j]].pT();
98         if (vNow > vMax) {
99           vMax = vNow;
100           jMax = j;
101         }
102       }
103       if (jMax != i) swap(vecIn[i], vecIn[jMax]);
104     }
105   }
106
107   // Master switch for merging
108   bool   doMerge;
109
110   // Maximum and current number of jets
111   int    nJetMax, nJet;
112
113   // Jet algorithm parameters
114   int    jetAlgorithm;
115   double eTjetMin, coneRadius, etaJetMax, etaJetMaxAlgo;
116
117   // CellJet specific
118   int    nEta, nPhi;
119   double eTseed, eTthreshold;
120
121   // SlowJet specific
122   int    slowJetPower;
123
124   // Merging procedure parameters
125   int    jetAllow, jetMatch, exclusiveMode;
126   double coneMatchLight, coneRadiusHeavy, coneMatchHeavy;
127   bool   exclusive;
128
129   // Event records to store original incoming process, final-state of the
130   // incoming process and what will be passed to the jet algorithm.
131   // Not completely necessary to store all steps, but makes tracking the
132   // steps of the algorithm a lot easier.
133   Event eventProcessOrig, eventProcess, workEventJet;
134
135   // Internal jet algorithms
136   CellJet *cellJet;
137   SlowJet *slowJet;
138
139   // Sort final-state of incoming process into light/heavy jets and 'other'
140   vector < int > typeIdx[3];
141   set    < int > typeSet[3];
142
143   // Momenta output of jet algorithm (to provide same output regardless of
144   // the selected jet algorithm)
145   vector < Vec4 > jetMomenta;
146
147   // Store the minimum eT/pT of matched light jets
148   double eTpTlightMin;
149
150   // Constants
151   static const double GHOSTENERGY, ZEROTHRESHOLD;
152 };
153
154 //--------------------------------------------------------------------------
155
156 // Main implementation of MLMhooks class. This may be split out to a
157 // separate C++ file if desired, but currently included here for ease
158 // of use.
159
160 // Constants: could be changed here if desired, but normally should not.
161 // These are of technical nature, as described for each.
162
163 // The energy of ghost particles. For technical reasons, this cannot be
164 // set arbitrarily low, see 'Particle::TINY' in 'Event.cc' for details.
165 const double MLMhooks::GHOSTENERGY   = 1e-15;
166
167 // A zero threshold value for double comparisons.
168 const double MLMhooks::ZEROTHRESHOLD = 1e-10;
169
170 //--------------------------------------------------------------------------
171
172 // Initialisation routine automatically called from Pythia::init().
173 // Setup all parts needed for the merging.
174
175 bool MLMhooks::initAfterBeams() {
176
177   // Read in parameters
178   doMerge         = settingsPtr->flag("MLM:merge");
179   nJet            = settingsPtr->mode("MLM:nJet");
180   nJetMax         = settingsPtr->mode("MLM:nJetMax");
181   jetAlgorithm    = settingsPtr->mode("MLM:jetAlgorithm");
182   eTjetMin        = settingsPtr->parm("MLM:eTjetMin");
183   coneRadius      = settingsPtr->parm("MLM:coneRadius");
184   etaJetMax       = settingsPtr->parm("MLM:etaJetMax");
185   // Use etaJetMax + coneRadius in input to jet algorithms
186   etaJetMaxAlgo   = etaJetMax + coneRadius;
187   // CellJet specific
188   nEta            = settingsPtr->mode("MLM:nEta");
189   nPhi            = settingsPtr->mode("MLM:nPhi");
190   eTseed          = settingsPtr->parm("MLM:eTseed");
191   eTthreshold     = settingsPtr->parm("MLM:eTthreshold");
192   // SlowJet specific
193   slowJetPower    = settingsPtr->mode("MLM:slowJetPower");
194   // Matching procedure
195   jetAllow        = settingsPtr->mode("MLM:jetAllow");
196   jetMatch        = settingsPtr->mode("MLM:jetMatch");
197   coneMatchLight  = settingsPtr->parm("MLM:coneMatchLight");
198   coneRadiusHeavy = settingsPtr->parm("MLM:coneRadiusHeavy");
199   if (coneRadiusHeavy < 0.) coneRadiusHeavy = coneRadius;
200   coneMatchHeavy  = settingsPtr->parm("MLM:coneMatchHeavy");
201   exclusiveMode   = settingsPtr->mode("MLM:exclusive");
202
203   // If not merging, then done
204   if (!doMerge) return true;
205
206   // Exclusive mode; if set to 2, then set based on nJet/nJetMax
207   if (exclusiveMode == 2) {
208
209     // No nJet or nJetMax, so default to exclusive mode
210     if (nJet < 0 || nJetMax < 0) {
211       infoPtr->errorMsg("Warning in MLMhooks:init: "
212           "missing jet multiplicity information; running in exclusive mode");
213       exclusive = true;
214
215     // Inclusive if nJet == nJetMax, exclusive otherwise
216     } else {
217       exclusive = (nJet == nJetMax) ? false : true;
218     }
219
220   // Otherwise, just set as given
221   } else {
222     exclusive = (exclusiveMode == 0) ? false : true;
223   }
224
225   // Initialise chosen jet algorithm. CellJet.
226   if (jetAlgorithm == 1) {
227
228     // Extra options for CellJet. nSel = 1 means that all final-state
229     // particles are taken and we retain control of what to select.
230     // smear/resolution/upperCut are not used and are set to default values.
231     int    nSel = 2, smear = 0;
232     double resolution = 0.5, upperCut = 2.;
233     cellJet = new CellJet(etaJetMaxAlgo, nEta, nPhi, nSel,
234                           smear, resolution, upperCut, eTthreshold);
235
236   // SlowJet
237   } else if (jetAlgorithm == 2) {
238     slowJet = new SlowJet(slowJetPower, coneRadius, eTjetMin, etaJetMaxAlgo);
239   }
240
241   // Check the jetMatch parameter; option 2 only works with SlowJet
242   if (jetAlgorithm == 1 && jetMatch == 2) {
243     infoPtr->errorMsg("Warning in MLMhooks:init: "
244         "jetMatch = 2 only valid with SlowJet algorithm. "
245         "Reverting to jetMatch = 1.");
246     jetMatch = 1;
247   }
248
249   // Setup local event records
250   eventProcessOrig.init("(eventProcessOrig)", particleDataPtr);
251   eventProcess.init("(eventProcess)", particleDataPtr);
252   workEventJet.init("(workEventJet)", particleDataPtr);
253
254   // Print information
255   string jetStr  = (jetAlgorithm ==  1) ? "CellJet" :
256                    (slowJetPower == -1) ? "anti-kT" :
257                    (slowJetPower ==  0) ? "C/A"     :
258                    (slowJetPower ==  1) ? "kT"      : "unknown";
259   string modeStr = (exclusive)         ? "exclusive" : "inclusive";
260   stringstream nJetStr, nJetMaxStr;
261   if (nJet >= 0)    nJetStr    << nJet;    else nJetStr    << "unknown";
262   if (nJetMax >= 0) nJetMaxStr << nJetMax; else nJetMaxStr << "unknown";
263   cout << endl
264        << " *-------  MLM matching parameters  -------*" << endl
265        << " |  nJet                |  " << setw(14)
266        << nJetStr.str() << "  |" << endl
267        << " |  nJetMax             |  " << setw(14)
268        << nJetMaxStr.str() << "  |" << endl
269        << " |  Jet algorithm       |  " << setw(14)
270        << jetStr << "  |" << endl
271        << " |  eTjetMin            |  " << setw(14)
272        << eTjetMin << "  |" << endl
273        << " |  coneRadius          |  " << setw(14)
274        << coneRadius << "  |" << endl
275        << " |  etaJetMax           |  " << setw(14)
276        << etaJetMax << "  |" << endl
277        << " |  jetAllow            |  " << setw(14)
278        << jetAllow << "  |" << endl
279        << " |  jetMatch            |  " << setw(14)
280        << jetMatch << "  |" << endl
281        << " |  coneMatchLight      |  " << setw(14)
282        << coneMatchLight << "  |" << endl
283        << " |  coneRadiusHeavy     |  " << setw(14)
284        << coneRadiusHeavy << "  |" << endl
285        << " |  coneMatchHeavy      |  " << setw(14)
286        << coneMatchHeavy << "  |" << endl
287        << " |  Mode                |  " << setw(14)
288        << modeStr << "  |" << endl
289        << " *-----------------------------------------*" << endl;
290
291   return true;
292 }
293
294 //--------------------------------------------------------------------------
295
296 // Process level veto. Stores incoming event for later.
297
298 bool MLMhooks::canVetoProcessLevel() { return doMerge; }
299
300 bool MLMhooks::doVetoProcessLevel(Event& process) { 
301
302   // Copy incoming process
303   eventProcessOrig = process;
304   return false;
305 }
306
307 //--------------------------------------------------------------------------
308
309 // Early parton level veto (before beam remnants and resonance showers)
310
311 bool MLMhooks::canVetoPartonLevelEarly() { return doMerge; }
312
313 bool MLMhooks::doVetoPartonLevelEarly(const Event& event) {
314
315   // 1) Sort the original incoming process. After this step is performed,
316   //    the following assignments have been made:
317   //      eventProcessOrig - the original incoming process
318   //      eventProcess     - the final-state of the incoming process with
319   //                         resonance decays removed (and resonances
320   //                         themselves now with positive status code)
321   //      typeIdx[0/1/2]   - indices into 'eventProcess' of
322   //                         light jets/heavy jets/other
323   //      typeSet[0/1/2]   - indices into 'event' of
324   //                         light jets/heavy jets/other
325   //      workEvent        - partons from the hardest subsystem
326   //                         + ISR + FSR only
327   sortIncomingProcess(event);
328
329   // Debug
330 #ifdef MLMDEBUG
331   // Begin
332   cout << endl << "---------- Begin MLM Debug ----------" << endl;
333
334   // Original incoming process
335   cout << endl << "Original incoming process:";
336   eventProcessOrig.list();
337
338   // Final-state of original incoming process
339   cout << endl << "Final-state incoming process:";
340   eventProcess.list();
341
342   // List categories of sorted particles
343   for (size_t i = 0; i < typeIdx[0].size(); i++) 
344     cout << ((i == 0) ? "Light jets: " : ", ")
345          << setw(3) << typeIdx[0][i];
346   for (size_t i = 0; i < typeIdx[1].size(); i++) 
347     cout << ((i == 0) ? "\nHeavy jets: " : ", ")
348          << setw(3) << typeIdx[1][i];
349   for (size_t i = 0; i < typeIdx[2].size(); i++) 
350     cout << ((i == 0) ? "\nOther:      " : ", ")
351          << setw(3) << typeIdx[2][i];
352
353   // Full event at this stage
354   cout << endl << endl << "Event:";
355   event.list();
356
357   // Work event (partons from hardest subsystem + ISR + FSR)
358   cout << endl << "Work event:";
359   workEvent.list();
360 #endif
361
362   // 2) Light/heavy jets: iType = 0 (light jets), 1 (heavy jets)
363   int iTypeEnd = (typeIdx[1].empty()) ? 1 : 2;
364   for (int iType = 0; iType < iTypeEnd; iType++) {
365
366     // 2a) Find particles which will be passed from the jet algorithm.
367     //     Input from 'workEvent' and output in 'workEventJet'.
368     jetAlgorithmInput(event, iType);
369
370     // Debug
371 #ifdef MLMDEBUG
372     // Jet algorithm event
373     cout << endl << "Jet algorithm event (iType = " << iType << "):";
374     workEventJet.list();
375 #endif
376
377     // 2b) Run jet algorithm on 'workEventJet'.
378     //     Output is stored in jetMomenta.
379     runJetAlgorithm();
380
381     // 2c) Match partons to jets and decide if veto is necessary
382     if (matchPartonsToJets(iType) == true) {
383       // Debug
384 #ifdef MLMDEBUG
385       cout << endl << "Event vetoed" << endl
386            << "----------  End MLM Debug  ----------" << endl;
387 #endif
388       return true;
389     }
390   }
391
392   // Debug
393 #ifdef MLMDEBUG
394   cout << endl << "Event accepted" << endl
395        << "----------  End MLM Debug  ----------" << endl;
396 #endif
397
398   // If we reached here, then no veto
399   return false;
400 }
401
402 //--------------------------------------------------------------------------
403
404 // Step (1): sort the incoming particles
405
406 void MLMhooks::sortIncomingProcess(const Event &event) {
407
408   // Remove resonance decays from original process and keep only final
409   // state. Resonances will have positive status code after this step.
410   omitResonanceDecays(eventProcessOrig, true);
411   eventProcess = workEvent;
412
413   // Sort original process final state into light/heavy jets and 'other'.
414   // Criteria:
415   //   1 <= ID <= 5 and massless, or ID == 21 --> light jet (typeIdx[0])
416   //   4 <= ID <= 6 and massive               --> heavy jet (typeIdx[1])
417   //   All else                               --> other     (typeIdx[2])
418   // Note that 'typeIdx' stores indices into 'eventProcess' (after resonance
419   // decays are omitted), while 'typeSet' stores indices into the original
420   // process record, 'eventProcessOrig', but these indices are also valid
421   // in 'event'.
422   for (int i = 0; i < 3; i++) {
423     typeIdx[i].clear();
424     typeSet[i].clear();
425   }
426   for (int i = 0; i < eventProcess.size(); i++) {
427     // Ignore nonfinal and default to 'other'
428     if (!eventProcess[i].isFinal()) continue;
429     int idx = 2;
430
431     // Light jets
432     if (eventProcess[i].id() == 21 || (eventProcess[i].idAbs() <= 5 &&
433         abs(eventProcess[i].m()) < ZEROTHRESHOLD))
434       idx = 0;
435
436     // Heavy jets
437     else if (eventProcess[i].idAbs() >= 4 && eventProcess[i].idAbs() <= 6)
438       idx = 1;
439
440     // Store
441     typeIdx[idx].push_back(i);
442     typeSet[idx].insert(eventProcess[i].daughter1());
443   }
444
445   // Extract partons from hardest subsystem + ISR + FSR only into
446   // workEvent. Note no resonance showers or MPIs.
447   subEvent(event);
448 }
449
450 //--------------------------------------------------------------------------
451
452 // Step (2a): pick which particles to pass to the jet algorithm
453
454 void MLMhooks::jetAlgorithmInput(const Event &event, int iType) {
455
456   // Take input from 'workEvent' and put output in 'workEventJet'
457   workEventJet = workEvent;
458
459   // Loop over particles and decide what to pass to the jet algorithm
460   for (int i = 0; i < workEventJet.size(); ++i) {
461     if (!workEventJet[i].isFinal()) continue;
462
463     // jetAllow option to disallow certain particle types
464     if (jetAllow == 1) {
465
466       // Original AG+Py6 algorithm explicitly excludes tops,
467       // leptons and photons.
468       int id = workEventJet[i].idAbs();
469       if ((id >= 11 && id <= 16) || id == 6 || id == 22) {
470         workEventJet[i].statusNeg();
471         continue;
472       }
473     }
474
475     // Get the index of this particle in original event
476     int idx = workEventJet[i].daughter1();
477
478     // Start with particle idx, and afterwards track mothers
479     while (true) {
480
481       // Light jets
482       if (iType == 0) {
483
484         // Do not include if originates from heavy jet or 'other'
485         if (typeSet[1].find(idx) != typeSet[1].end() ||
486             typeSet[2].find(idx) != typeSet[2].end()) {
487           workEventJet[i].statusNeg();
488           break;
489         }
490
491         // Made it to start of event record so done
492         if (idx == 0) break;
493         // Otherwise next mother and continue
494         idx = event[idx].mother1();
495
496       // Heavy jets
497       } else if (iType == 1) {
498
499         // Only include if originates from heavy jet
500         if (typeSet[1].find(idx) != typeSet[1].end()) break;
501
502         // Made it to start of event record with no heavy jet mother,
503         // so DO NOT include particle
504         if (idx == 0) {
505           workEventJet[i].statusNeg();
506           break;
507         }
508
509         // Otherwise next mother and continue
510         idx = event[idx].mother1();
511
512       } // if (iType)
513     } // while (true)
514   } // for (i)
515
516   // For jetMatch = 2, insert ghost particles corresponding to
517   // each hard parton in the original process
518   if (jetMatch == 2) {
519     for (int i = 0; i < int(typeIdx[iType].size()); i++) {
520       // Get y/phi of the parton
521       Vec4   pIn = eventProcess[typeIdx[iType][i]].p();
522       double y   = Vec4y(pIn);
523       double phi = pIn.phi();
524
525       // Create a ghost particle and add to the workEventJet
526       double e   = GHOSTENERGY;
527       double e2y = exp(2. * y);
528       double pz  = e * (e2y - 1.) / (e2y + 1.);
529       double pt  = sqrt(e*e - pz*pz);
530       double px  = pt * cos(phi);
531       double py  = pt * sin(phi);
532       workEventJet.append(Particle(21, 99, 0, 0, 0, 0, 0, 0,
533                                 px, py, pz, e, 0., 0, 9.));
534
535       // Extra check on reconstructed y/phi values. If many warnings
536       // of this type, GHOSTENERGY may be set too low.
537 #ifdef MLMCHECK
538       int lastIdx = workEventJet.size() - 1;
539       if (abs(y   - workEventJet[lastIdx].y())   > ZEROTHRESHOLD ||
540           abs(phi - workEventJet[lastIdx].phi()) > ZEROTHRESHOLD)
541         infoPtr->errorMsg("Warning in MLMhooks:jetAlgorithmInput: "
542             "ghost particle y/phi mismatch");
543 #endif
544
545     } // for (i)
546   } // if (jetMatch == 2)
547 }
548
549 //--------------------------------------------------------------------------
550
551 // Step (2b): run jet algorithm and provide common output
552
553 void MLMhooks::runJetAlgorithm() {
554
555   // Run the jet clustering algorithm
556   if (jetAlgorithm == 1)
557     cellJet->analyze(workEventJet, eTjetMin, coneRadius, eTseed);
558   else
559     slowJet->analyze(workEventJet);
560
561   // Extract four-momenta of jets with |eta| < etaJetMax and
562   // put into jetMomenta. Note that this is done backwards as
563   // jets are removed with SlowJet.
564   jetMomenta.clear();
565   int iJet = (jetAlgorithm == 1) ? cellJet->size() - 1:
566                                    slowJet->sizeJet() - 1;
567   for (int i = iJet; i > -1; i--) {
568     Vec4 jetMom = (jetAlgorithm == 1) ? cellJet->pMassive(i) :
569                                         slowJet->p(i);
570     double eta = Vec4eta(jetMom);
571
572     if (abs(eta) > etaJetMax) {
573       if (jetAlgorithm == 2) slowJet->removeJet(i);
574       continue;
575     }
576     jetMomenta.push_back(jetMom);
577   }
578
579   // Reverse jetMomenta to restore eT/pT ordering
580   reverse(jetMomenta.begin(), jetMomenta.end());
581 }
582
583 //--------------------------------------------------------------------------
584
585 // Step (2c): veto decision (returning true vetoes the event)
586
587 bool MLMhooks::matchPartonsToJets(int iType) {
588
589   // Use two different routines for light/heavy jets as
590   // different veto conditions and for clarity
591   if (iType == 0) return (matchPartonsToJetsLight() > 0);
592   else            return (matchPartonsToJetsHeavy() > 0);
593 }
594
595 //--------------------------------------------------------------------------
596
597 // Step(2c): light jets
598 // Return codes are given indicating the reason for a veto.
599 // Although not currently used, they are a useful debugging tool:
600 //   0 = no veto
601 //   1 = veto as number of jets less than number of partons
602 //   2 = veto as exclusive mode and number of jets greater than
603 //       number of partons
604 //   3 = veto as inclusive mode and there would be an extra jet
605 //       that is harder than any matched soft jet
606 //   4 = veto as there is a parton which does not match a jet
607
608 int MLMhooks::matchPartonsToJetsLight() {
609
610   // Always veto if number of jets is less than original number of jets
611   if (jetMomenta.size() < typeIdx[0].size()) return 1;
612   // Veto if in exclusive mode and number of jets bigger than original
613   if (exclusive && jetMomenta.size() > typeIdx[0].size()) return 2;
614
615   // Sort partons by eT/pT
616   sortTypeIdx(typeIdx[0]);
617
618   // Number of hard partons
619   int nParton = typeIdx[0].size();
620
621   // Keep track of which jets have been assigned a hard parton
622   vector < bool > jetAssigned;
623   jetAssigned.assign(jetMomenta.size(), false);
624
625   // Jet matching procedure: (1) deltaR between partons and jets
626   if (jetMatch == 1) {
627
628     // Loop over light hard partons and get 4-momentum
629     for (int i = 0; i < nParton; i++) {
630       Vec4 p1 = eventProcess[typeIdx[0][i]].p();
631
632       // Track which jet has the minimal dR measure with this parton
633       int    jMin  = -1;
634       double dRmin = 0.;
635
636       // Loop over all jets (skipping those already assigned).
637       for (int j = 0; j < int(jetMomenta.size()); j++) {
638         if (jetAssigned[j]) continue;
639
640         // DeltaR between parton/jet and store if minimum
641         double dR = (jetAlgorithm == 1) ?
642             deltaReta(p1, jetMomenta[j]) : deltaRy(p1, jetMomenta[j]);
643         if (jMin < 0 || dR < dRmin) {
644           dRmin = dR;
645           jMin  = j;
646         }
647       } // for (j)
648
649       // Check for jet-parton match
650       if (jMin >= 0 && dRmin < coneRadius * coneMatchLight) {
651
652         // If the matched jet is not one of the nParton hardest jets,
653         // the extra left over jet would be harder than some of the
654         // matched jets. This is disallowed, so veto.
655         if (jMin >= nParton) return 3;
656
657         // Mark jet as assigned.
658         jetAssigned[jMin] = true;
659
660       // If no match, then event will be vetoed in all cases
661       } else return 4;
662
663     } // for (i)
664
665   // Jet matching procedure: (2) ghost particles in SlowJet
666   } else {
667
668     // Loop over added 'ghost' particles and find if assigned to a jet
669     for (int i = workEventJet.size() - nParton;
670         i < workEventJet.size(); i++) {
671       int jMin = slowJet->jetAssignment(i);
672
673       // Veto if:
674       //  1) not one of nParton hardest jets
675       //  2) not assigned to a jet
676       //  3) jet has already been assigned
677       if (jMin >= nParton)               return 3;
678       if (jMin < 0 || jetAssigned[jMin]) return 4;
679
680       // Mark jet as assigned
681       jetAssigned[jMin] = true;
682
683     } // for (i)
684   } // if (jetMatch)
685
686   // Minimal eT/pT (CellJet/SlowJet) of matched light jets. Needed
687   // later for heavy jet vetos in inclusive mode.
688   if (nParton > 0)
689     eTpTlightMin = (jetAlgorithm == 1) ? jetMomenta[nParton - 1].eT()
690                                        : jetMomenta[nParton - 1].pT();
691   else
692     eTpTlightMin = -1.;
693
694   // No veto
695   return 0;
696 }
697
698 //--------------------------------------------------------------------------
699
700 // Step(2c): heavy jets
701 // Return codes are given indicating the reason for a veto.
702 // Although not currently used, they are a useful debugging tool:
703 //   0 = no veto as there are no extra jets present
704 //   1 = veto as in exclusive mode and extra jets present
705 //   2 = veto as in inclusive mode and extra jets were harder
706 //       than any matched light jet
707
708 int MLMhooks::matchPartonsToJetsHeavy() {
709
710   // If there are no extra jets, then accept
711   if (jetMomenta.empty()) return 0;
712
713   // Number of hard partons
714   int nParton = typeIdx[1].size();
715
716   // Remove jets that are close to heavy quarks
717   set < int > removeJets;
718
719   // Jet matching procedure: (1) deltaR between partons and jets
720   if (jetMatch == 1) {
721
722     // Loop over heavy hard partons and get 4-momentum
723     for (int i = 0; i < nParton; i++) {
724       Vec4 p1 = eventProcess[typeIdx[1][i]].p();
725
726       // Loop over all jets, find dR and mark for removal if match
727       for (int j = 0; j < int(jetMomenta.size()); j++) {
728         double dR = (jetAlgorithm == 1) ?
729             deltaReta(p1, jetMomenta[j]) : deltaRy(p1, jetMomenta[j]);
730         if (dR < coneRadiusHeavy * coneMatchHeavy)
731           removeJets.insert(j);
732
733       } // for (j)
734     } // for (i)
735
736   // Jet matching procedure: (2) ghost particles in SlowJet
737   } else {
738
739     // Loop over added 'ghost' particles and if assigned to a jet
740     // then mark this jet for removal
741     for (int i = workEventJet.size() - nParton;
742         i < workEventJet.size(); i++) {
743       int jMin = slowJet->jetAssignment(i);
744       if (jMin >= 0) removeJets.insert(jMin);
745     }
746       
747   }
748
749   // Remove jets (backwards order to not disturb indices)
750   for (set < int >::reverse_iterator it  = removeJets.rbegin();
751                                      it != removeJets.rend(); it++)
752     jetMomenta.erase(jetMomenta.begin() + *it);
753
754   // Handle case if there are still extra jets
755   if (!jetMomenta.empty()) {
756
757     // Exclusive mode, so immediate veto
758     if (exclusive) return 1;
759
760     // Inclusive mode; extra jets must be softer than any matched light jet
761     else if (eTpTlightMin >= 0.)
762       for (size_t j = 0; j < jetMomenta.size(); j++) {
763         // CellJet uses eT, SlowJet uses pT
764         if ( (jetAlgorithm == 1 && jetMomenta[j].eT() > eTpTlightMin) ||
765              (jetAlgorithm == 2 && jetMomenta[j].pT() > eTpTlightMin) )
766           return 2;
767       }
768
769   } // if (!jetMomenta.empty())
770
771   // No extra jets were present so no veto
772   return 0;
773 }
774
775 //==========================================================================
776
777 #endif // MLMHOOKS_H