PYTHIA8/pythia8175/examples/JetMatching.h

   1 // JetMatching.h is a part of the PYTHIA event generator.
   2 // Copyright (C) 2013 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 // 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 // This file provides the classes to perform MLM matching of
  11 // Alpgen or MadGraph 5 input.
  12 // Example usage is shown in main32.cc, and further details
  13 // can be found in the 'Jet Matching Style' manual page.
  14
  15 #ifndef Pythia8_JetMatching_H
  16 #define Pythia8_JetMatching_H
  17
  18 // Includes
  19 #include "Pythia.h"
  20 #include "GeneratorInput.h"
  21 using namespace Pythia8;
  22
  23 //==========================================================================
  24
  25 // Declaration of main JetMatching class to perform MLM matching.
  26 // Note that it is defined with virtual inheritance, so that it can
  27 // be combined with other UserHooks classes, see e.g. main33.cc.
  28
  29 class JetMatching : virtual public UserHooks {
  30
  31 public:
  32
  33   // Constructor and destructor
  34  JetMatching() : cellJet(NULL), slowJet(NULL) {}
  35   ~JetMatching() {
  36     if (cellJet) delete cellJet;
  37     if (slowJet) delete slowJet;
  38   }
  39
  40   // Initialisation
  41   virtual bool initAfterBeams() = 0;
  42
  43   // Process level vetos
  44   bool canVetoProcessLevel() { return doMerge; }
  45   bool doVetoProcessLevel(Event& process) {
  46     eventProcessOrig = process;
  47     return false;
  48   }
  49
  50   // Parton level vetos (before beam remnants and resonance decays)
  51   bool canVetoPartonLevelEarly() { return doMerge; }
  52   bool doVetoPartonLevelEarly(const Event& event);
  53
  54 protected:
  55
  56   // Constants to be changed for debug printout or extra checks.
  57   static const bool MATCHINGDEBUG, MATCHINGCHECK;
  58
  59   // Different steps of the matching algorithm.
  60   virtual void sortIncomingProcess(const Event &)=0;
  61   virtual void jetAlgorithmInput(const Event &, int)=0;
  62   virtual void runJetAlgorithm()=0;
  63   virtual bool matchPartonsToJets(int)=0;
  64   virtual int  matchPartonsToJetsLight()=0;
  65   virtual int  matchPartonsToJetsHeavy()=0;
  66
  67   enum vetoStatus { NONE, LESS_JETS, MORE_JETS, HARD_JET, UNMATCHED_PARTON };
  68   enum partonTypes { ID_CHARM=4, ID_BOT=5, ID_TOP=6, ID_LEPMIN=11,
  69     ID_LEPMAX=16, ID_GLUON=21, ID_PHOTON=22 };
  70
  71   // Master switch for merging
  72   bool   doMerge;
  73
  74   // Maximum and current number of jets
  75   int    nJetMax, nJet;
  76
  77   // Jet algorithm parameters
  78   int    jetAlgorithm;
  79   double eTjetMin, coneRadius, etaJetMax, etaJetMaxAlgo;
  80
  81   // Internal jet algorithms
  82   CellJet* cellJet;
  83   SlowJet* slowJet;
  84
  85   // SlowJet specific
  86   int    slowJetPower;
  87
  88   // Event records to store original incoming process, final-state of the
  89   // incoming process and what will be passed to the jet algorithm.
  90   // Not completely necessary to store all steps, but makes tracking the
  91   // steps of the algorithm a lot easier.
  92   Event eventProcessOrig, eventProcess, workEventJet;
  93
  94   // Sort final-state of incoming process into light/heavy jets and 'other'
  95   vector<int> typeIdx[3];
  96   set<int>    typeSet[3];
  97
  98   // Momenta output of jet algorithm (to provide same output regardless of
  99   // the selected jet algorithm)
 100   vector<Vec4> jetMomenta;
 101
 102   // CellJet specific
 103   int    nEta, nPhi;
 104   double eTseed, eTthreshold;
 105
 106   // Merging procedure parameters
 107   int    jetAllow, jetMatch, exclusiveMode;
 108   double coneMatchLight, coneRadiusHeavy, coneMatchHeavy;
 109   bool   exclusive;
 110
 111   // Store the minimum eT/pT of matched light jets
 112   double eTpTlightMin;
 113
 114 };
 115
 116 //==========================================================================
 117
 118 // Declaration of main UserHooks class to perform Alpgen matching.
 119
 120 class JetMatchingAlpgen : virtual public JetMatching {
 121
 122 public:
 123
 124   // Constructor and destructor
 125   JetMatchingAlpgen() { }
 126   ~JetMatchingAlpgen() { }
 127
 128   // Initialisation
 129   bool initAfterBeams();
 130
 131 private:
 132
 133   // Different steps of the matching algorithm.
 134   void sortIncomingProcess(const Event &);
 135   void jetAlgorithmInput(const Event &, int);
 136   void runJetAlgorithm();
 137   bool matchPartonsToJets(int);
 138   int  matchPartonsToJetsLight();
 139   int  matchPartonsToJetsHeavy();
 140
 141   // Sorting utility
 142   void sortTypeIdx(vector < int > &vecIn);
 143
 144   // Constants
 145   static const double GHOSTENERGY, ZEROTHRESHOLD;
 146
 147 };
 148
 149 //==========================================================================
 150
 151 // Declaration of main UserHooks class to perform Madgraph matching.
 152
 153 class JetMatchingMadgraph : virtual public JetMatching {
 154
 155 public:
 156
 157   // Constructor and destructor
 158   JetMatchingMadgraph() { }
 159   ~JetMatchingMadgraph() { }
 160
 161   // Initialisation
 162   bool initAfterBeams();
 163
 164 protected:
 165
 166   // Different steps of the matching algorithm.
 167   void sortIncomingProcess(const Event &);
 168   void jetAlgorithmInput(const Event &, int);
 169   void runJetAlgorithm();
 170   bool matchPartonsToJets(int);
 171   int  matchPartonsToJetsLight();
 172   int  matchPartonsToJetsHeavy();
 173
 174   // Variables.
 175   int    nQmatch;
 176   double qCut, qCutSq, clFact;
 177
 178 };
 179
 180 //==========================================================================
 181
 182 // Main implementation of JetMatching class.
 183 // This may be split out to a separate C++ file if desired,
 184 // but currently included here for ease of use.
 185
 186 //--------------------------------------------------------------------------
 187
 188 // Constants to be changed for debug printout or extra checks.
 189 const bool JetMatching::MATCHINGDEBUG = false;
 190 const bool JetMatching::MATCHINGCHECK = false;
 191
 192 //--------------------------------------------------------------------------
 193
 194 // Early parton level veto (before beam remnants and resonance showers)
 195
 196 bool JetMatching::doVetoPartonLevelEarly(const Event& event) {
 197
 198   // 1) Sort the original incoming process. After this step is performed,
 199   //    the following assignments have been made:
 200   //    eventProcessOrig - the original incoming process
 201   //    eventProcess     - the final-state of the incoming process with
 202   //                       resonance decays removed (and resonances
 203   //                       themselves now with positive status code)
 204   //    typeIdx[0/1/2]   - Indices into 'eventProcess' of
 205   //                       light jets/heavy jets/other
 206   //    typeSet[0/1/2]   - Indices into 'event' of light jets/heavy jets/other
 207   //    workEvent        - partons from the hardest subsystem + ISR + FSR only
 208   sortIncomingProcess(event);
 209
 210   // Debug printout.
 211   if (MATCHINGDEBUG) {
 212     // Begin
 213     cout << endl << "-------- Begin Madgraph Debug --------" << endl;
 214     // Original incoming process
 215     cout << endl << "Original incoming process:";
 216     eventProcessOrig.list();
 217     // Final-state of original incoming process
 218     cout << endl << "Final-state incoming process:";
 219     eventProcess.list();
 220     // List categories of sorted particles
 221     for (size_t i = 0; i < typeIdx[0].size(); i++)
 222       cout << ((i == 0) ? "Light jets: " : ", ")   << setw(3) << typeIdx[0][i];
 223     if( typeIdx[0].size()== 0 )
 224       cout << "Light jets: None";
 225
 226     for (size_t i = 0; i < typeIdx[1].size(); i++)
 227       cout << ((i == 0) ? "\nHeavy jets: " : ", ") << setw(3) << typeIdx[1][i];
 228     for (size_t i = 0; i < typeIdx[2].size(); i++)
 229       cout << ((i == 0) ? "\nOther:      " : ", ") << setw(3) << typeIdx[2][i];
 230     // Full event at this stage
 231     cout << endl << endl << "Event:";
 232     event.list();
 233     // Work event (partons from hardest subsystem + ISR + FSR)
 234     cout << endl << "Work event:";
 235     workEvent.list();
 236   }
 237
 238   // 2) Light/heavy jets: iType = 0 (light jets), 1 (heavy jets)
 239   int iTypeEnd = (typeIdx[1].empty()) ? 1 : 2;
 240   for (int iType = 0; iType < iTypeEnd; iType++) {
 241
 242     // 2a) Find particles which will be passed from the jet algorithm.
 243     //     Input from 'workEvent' and output in 'workEventJet'.
 244     jetAlgorithmInput(event, iType);
 245
 246     // Debug printout.
 247     if (MATCHINGDEBUG) {
 248       // Jet algorithm event
 249       cout << endl << "Jet algorithm event (iType = " << iType << "):";
 250       workEventJet.list();
 251     }
 252
 253     // 2b) Run jet algorithm on 'workEventJet'.
 254     //     Output is stored in jetMomenta.
 255     runJetAlgorithm();
 256
 257     // 2c) Match partons to jets and decide if veto is necessary
 258     if (matchPartonsToJets(iType) == true) {
 259       // Debug printout.
 260       if (MATCHINGDEBUG) {
 261         cout << endl << "Event vetoed" << endl
 262              << "----------  End MLM Debug  ----------" << endl;
 263       }
 264       return true;
 265     }
 266   }
 267
 268   // Debug printout.
 269   if (MATCHINGDEBUG) {
 270     cout << endl << "Event accepted" << endl
 271          << "----------  End MLM Debug  ----------" << endl;
 272   }
 273
 274   // If we reached here, then no veto
 275   return false;
 276
 277 }
 278
 279 //==========================================================================
 280
 281 // Main implementation of Alpgen UserHooks class.
 282 // This may be split out to a separate C++ file if desired,
 283 // but currently included here for ease of use.
 284
 285 //--------------------------------------------------------------------------
 286
 287 // Constants: could be changed here if desired, but normally should not.
 288 // These are of technical nature, as described for each.
 289
 290 // The energy of ghost particles. For technical reasons, this cannot be
 291 // set arbitrarily low, see 'Particle::TINY' in 'Event.cc' for details.
 292 const double JetMatchingAlpgen::GHOSTENERGY   = 1e-15;
 293
 294 // A zero threshold value for double comparisons.
 295 const double JetMatchingAlpgen::ZEROTHRESHOLD = 1e-10;
 296
 297 //--------------------------------------------------------------------------
 298
 299 // Function to sort typeIdx vectors into descending eT/pT order.
 300 // Uses a selection sort, as number of partons generally small
 301 // and so efficiency not a worry.
 302
 303 void JetMatchingAlpgen::sortTypeIdx(vector < int > &vecIn) {
 304   for (size_t i = 0; i < vecIn.size(); i++) {
 305     size_t jMax = i;
 306     double vMax = (jetAlgorithm == 1) ?
 307       eventProcess[vecIn[i]].eT() :
 308       eventProcess[vecIn[i]].pT();
 309     for (size_t j = i + 1; j < vecIn.size(); j++) {
 310       double vNow = (jetAlgorithm == 1) ?
 311         eventProcess[vecIn[j]].eT() :
 312         eventProcess[vecIn[j]].pT();
 313       if (vNow > vMax) {
 314         vMax = vNow;
 315         jMax = j;
 316       }
 317     }
 318     if (jMax != i) swap(vecIn[i], vecIn[jMax]);
 319   }
 320 }
 321
 322 //--------------------------------------------------------------------------
 323
 324 // Initialisation routine automatically called from Pythia::init().
 325 // Setup all parts needed for the merging.
 326
 327 bool JetMatchingAlpgen::initAfterBeams() {
 328
 329   // Read in parameters
 330   doMerge         = settingsPtr->flag("JetMatching:merge");
 331   jetAlgorithm    = settingsPtr->mode("JetMatching:jetAlgorithm");
 332   nJet            = settingsPtr->mode("JetMatching:nJet");
 333   nJetMax         = settingsPtr->mode("JetMatching:nJetMax");
 334   eTjetMin        = settingsPtr->parm("JetMatching:eTjetMin");
 335   coneRadius      = settingsPtr->parm("JetMatching:coneRadius");
 336   etaJetMax       = settingsPtr->parm("JetMatching:etaJetMax");
 337
 338   // Use etaJetMax + coneRadius in input to jet algorithms
 339   etaJetMaxAlgo   = etaJetMax + coneRadius;
 340
 341   // CellJet specific
 342   nEta            = settingsPtr->mode("JetMatching:nEta");
 343   nPhi            = settingsPtr->mode("JetMatching:nPhi");
 344   eTseed          = settingsPtr->parm("JetMatching:eTseed");
 345   eTthreshold     = settingsPtr->parm("JetMatching:eTthreshold");
 346
 347   // SlowJet specific
 348   slowJetPower    = settingsPtr->mode("JetMatching:slowJetPower");
 349   coneMatchLight  = settingsPtr->parm("JetMatching:coneMatchLight");
 350   coneRadiusHeavy = settingsPtr->parm("JetMatching:coneRadiusHeavy");
 351   if (coneRadiusHeavy < 0.) coneRadiusHeavy = coneRadius;
 352   coneMatchHeavy  = settingsPtr->parm("JetMatching:coneMatchHeavy");
 353
 354   // Matching procedure
 355   jetAllow        = settingsPtr->mode("JetMatching:jetAllow");
 356   jetMatch        = settingsPtr->mode("JetMatching:jetMatch");
 357   exclusiveMode   = settingsPtr->mode("JetMatching:exclusive");
 358
 359   // If not merging, then done
 360   if (!doMerge) return true;
 361
 362   // Exclusive mode; if set to 2, then set based on nJet/nJetMax
 363   if (exclusiveMode == 2) {
 364
 365     // No nJet or nJetMax, so default to exclusive mode
 366     if (nJet < 0 || nJetMax < 0) {
 367       infoPtr->errorMsg("Warning in JetMatchingAlpgen:init: "
 368           "missing jet multiplicity information; running in exclusive mode");
 369       exclusive = true;
 370
 371     // Inclusive if nJet == nJetMax, exclusive otherwise
 372     } else {
 373       exclusive = (nJet == nJetMax) ? false : true;
 374     }
 375
 376   // Otherwise, just set as given
 377   } else {
 378     exclusive = (exclusiveMode == 0) ? false : true;
 379   }
 380
 381   // Initialise chosen jet algorithm. CellJet.
 382   if (jetAlgorithm == 1) {
 383
 384     // Extra options for CellJet. nSel = 1 means that all final-state
 385     // particles are taken and we retain control of what to select.
 386     // smear/resolution/upperCut are not used and are set to default values.
 387     int    nSel = 2, smear = 0;
 388     double resolution = 0.5, upperCut = 2.;
 389     cellJet = new CellJet(etaJetMaxAlgo, nEta, nPhi, nSel,
 390                           smear, resolution, upperCut, eTthreshold);
 391
 392   // SlowJet
 393   } else if (jetAlgorithm == 2) {
 394     slowJet = new SlowJet(slowJetPower, coneRadius, eTjetMin, etaJetMaxAlgo);
 395   }
 396
 397   // Check the jetMatch parameter; option 2 only works with SlowJet
 398   if (jetAlgorithm == 1 && jetMatch == 2) {
 399     infoPtr->errorMsg("Warning in JetMatchingAlpgen:init: "
 400         "jetMatch = 2 only valid with SlowJet algorithm. "
 401         "Reverting to jetMatch = 1");
 402     jetMatch = 1;
 403   }
 404
 405   // Setup local event records
 406   eventProcessOrig.init("(eventProcessOrig)", particleDataPtr);
 407   eventProcess.init("(eventProcess)", particleDataPtr);
 408   workEventJet.init("(workEventJet)", particleDataPtr);
 409
 410   // Print information
 411   string jetStr  = (jetAlgorithm ==  1) ? "CellJet" :
 412                    (slowJetPower == -1) ? "anti-kT" :
 413                    (slowJetPower ==  0) ? "C/A"     :
 414                    (slowJetPower ==  1) ? "kT"      : "unknown";
 415   string modeStr = (exclusive)         ? "exclusive" : "inclusive";
 416   stringstream nJetStr, nJetMaxStr;
 417   if (nJet >= 0)    nJetStr    << nJet;    else nJetStr    << "unknown";
 418   if (nJetMax >= 0) nJetMaxStr << nJetMax; else nJetMaxStr << "unknown";
 419   cout << endl
 420        << " *-------  MLM matching parameters  -------*" << endl
 421        << " |  nJet                |  " << setw(14)
 422        << nJetStr.str() << "  |" << endl
 423        << " |  nJetMax             |  " << setw(14)
 424        << nJetMaxStr.str() << "  |" << endl
 425        << " |  Jet algorithm       |  " << setw(14)
 426        << jetStr << "  |" << endl
 427        << " |  eTjetMin            |  " << setw(14)
 428        << eTjetMin << "  |" << endl
 429        << " |  coneRadius          |  " << setw(14)
 430        << coneRadius << "  |" << endl
 431        << " |  etaJetMax           |  " << setw(14)
 432        << etaJetMax << "  |" << endl
 433        << " |  jetAllow            |  " << setw(14)
 434        << jetAllow << "  |" << endl
 435        << " |  jetMatch            |  " << setw(14)
 436        << jetMatch << "  |" << endl
 437        << " |  coneMatchLight      |  " << setw(14)
 438        << coneMatchLight << "  |" << endl
 439        << " |  coneRadiusHeavy     |  " << setw(14)
 440        << coneRadiusHeavy << "  |" << endl
 441        << " |  coneMatchHeavy      |  " << setw(14)
 442        << coneMatchHeavy << "  |" << endl
 443        << " |  Mode                |  " << setw(14)
 444        << modeStr << "  |" << endl
 445        << " *-----------------------------------------*" << endl;
 446
 447   return true;
 448 }
 449
 450 //--------------------------------------------------------------------------
 451
 452 // Step (1): sort the incoming particles
 453
 454 void JetMatchingAlpgen::sortIncomingProcess(const Event &event) {
 455
 456   // Remove resonance decays from original process and keep only final
 457   // state. Resonances will have positive status code after this step.
 458   omitResonanceDecays(eventProcessOrig, true);
 459   eventProcess = workEvent;
 460
 461   // Sort original process final state into light/heavy jets and 'other'.
 462   // Criteria:
 463   //   1 <= ID <= 5 and massless, or ID == 21 --> light jet (typeIdx[0])
 464   //   4 <= ID <= 6 and massive               --> heavy jet (typeIdx[1])
 465   //   All else                               --> other     (typeIdx[2])
 466   // Note that 'typeIdx' stores indices into 'eventProcess' (after resonance
 467   // decays are omitted), while 'typeSet' stores indices into the original
 468   // process record, 'eventProcessOrig', but these indices are also valid
 469   // in 'event'.
 470   for (int i = 0; i < 3; i++) {
 471     typeIdx[i].clear();
 472     typeSet[i].clear();
 473   }
 474   for (int i = 0; i < eventProcess.size(); i++) {
 475     // Ignore nonfinal and default to 'other'
 476     if (!eventProcess[i].isFinal()) continue;
 477     int idx = 2;
 478
 479     // Light jets
 480     if (eventProcess[i].id() == ID_GLUON
 481       || (eventProcess[i].idAbs() <= ID_BOT
 482       && abs(eventProcess[i].m()) < ZEROTHRESHOLD)) idx = 0;
 483
 484     // Heavy jets
 485     else if (eventProcess[i].idAbs() >= ID_CHARM
 486       && eventProcess[i].idAbs() <= ID_TOP) idx = 1;
 487
 488     // Store
 489     typeIdx[idx].push_back(i);
 490     typeSet[idx].insert(eventProcess[i].daughter1());
 491   }
 492
 493   // Extract partons from hardest subsystem + ISR + FSR only into
 494   // workEvent. Note no resonance showers or MPIs.
 495   subEvent(event);
 496 }
 497
 498 //--------------------------------------------------------------------------
 499
 500 // Step (2a): pick which particles to pass to the jet algorithm
 501
 502 void JetMatchingAlpgen::jetAlgorithmInput(const Event &event, int iType) {
 503
 504   // Take input from 'workEvent' and put output in 'workEventJet'
 505   workEventJet = workEvent;
 506
 507   // Loop over particles and decide what to pass to the jet algorithm
 508   for (int i = 0; i < workEventJet.size(); ++i) {
 509     if (!workEventJet[i].isFinal()) continue;
 510
 511     // jetAllow option to disallow certain particle types
 512     if (jetAllow == 1) {
 513
 514       // Original AG+Py6 algorithm explicitly excludes tops,
 515       // leptons and photons.
 516       int id = workEventJet[i].idAbs();
 517       if ( (id >= ID_LEPMIN && id <= ID_LEPMAX) || id == ID_TOP
 518         || id == ID_PHOTON) {
 519         workEventJet[i].statusNeg();
 520         continue;
 521       }
 522     }
 523
 524     // Get the index of this particle in original event
 525     int idx = workEventJet[i].daughter1();
 526
 527     // Start with particle idx, and afterwards track mothers
 528     while (true) {
 529
 530       // Light jets
 531       if (iType == 0) {
 532
 533         // Do not include if originates from heavy jet or 'other'
 534         if (typeSet[1].find(idx) != typeSet[1].end() ||
 535             typeSet[2].find(idx) != typeSet[2].end()) {
 536           workEventJet[i].statusNeg();
 537           break;
 538         }
 539
 540         // Made it to start of event record so done
 541         if (idx == 0) break;
 542         // Otherwise next mother and continue
 543         idx = event[idx].mother1();
 544
 545       // Heavy jets
 546       } else if (iType == 1) {
 547
 548         // Only include if originates from heavy jet
 549         if (typeSet[1].find(idx) != typeSet[1].end()) break;
 550
 551         // Made it to start of event record with no heavy jet mother,
 552         // so DO NOT include particle
 553         if (idx == 0) {
 554           workEventJet[i].statusNeg();
 555           break;
 556         }
 557
 558         // Otherwise next mother and continue
 559         idx = event[idx].mother1();
 560
 561       } // if (iType)
 562     } // while (true)
 563   } // for (i)
 564
 565   // For jetMatch = 2, insert ghost particles corresponding to
 566   // each hard parton in the original process
 567   if (jetMatch == 2) {
 568     for (int i = 0; i < int(typeIdx[iType].size()); i++) {
 569       // Get y/phi of the parton
 570       Vec4   pIn = eventProcess[typeIdx[iType][i]].p();
 571       double y   = pIn.rap();
 572       double phi = pIn.phi();
 573
 574       // Create a ghost particle and add to the workEventJet
 575       double e   = GHOSTENERGY;
 576       double e2y = exp(2. * y);
 577       double pz  = e * (e2y - 1.) / (e2y + 1.);
 578       double pt  = sqrt(e*e - pz*pz);
 579       double px  = pt * cos(phi);
 580       double py  = pt * sin(phi);
 581       workEventJet.append( ID_GLUON, 99, 0, 0, 0, 0, 0, 0, px, py, pz, e);
 582
 583       // Extra check on reconstructed y/phi values. If many warnings
 584       // of this type, GHOSTENERGY may be set too low.
 585       if (MATCHINGCHECK) {
 586       int lastIdx = workEventJet.size() - 1;
 587       if (abs(y   - workEventJet[lastIdx].y())   > ZEROTHRESHOLD ||
 588           abs(phi - workEventJet[lastIdx].phi()) > ZEROTHRESHOLD)
 589         infoPtr->errorMsg("Warning in JetMatchingAlpgen:jetAlgorithmInput: "
 590             "ghost particle y/phi mismatch");
 591       }
 592
 593     } // for (i)
 594   } // if (jetMatch == 2)
 595 }
 596
 597 //--------------------------------------------------------------------------
 598
 599 // Step (2b): run jet algorithm and provide common output
 600
 601 void JetMatchingAlpgen::runJetAlgorithm() {
 602
 603   // Run the jet clustering algorithm
 604   if (jetAlgorithm == 1)
 605     cellJet->analyze(workEventJet, eTjetMin, coneRadius, eTseed);
 606   else
 607     slowJet->analyze(workEventJet);
 608
 609   // Extract four-momenta of jets with |eta| < etaJetMax and
 610   // put into jetMomenta. Note that this is done backwards as
 611   // jets are removed with SlowJet.
 612   jetMomenta.clear();
 613   int iJet = (jetAlgorithm == 1) ? cellJet->size() - 1:
 614                                    slowJet->sizeJet() - 1;
 615   for (int i = iJet; i > -1; i--) {
 616     Vec4 jetMom = (jetAlgorithm == 1) ? cellJet->pMassive(i) :
 617                                         slowJet->p(i);
 618     double eta = jetMom.eta();
 619
 620     if (abs(eta) > etaJetMax) {
 621       if (jetAlgorithm == 2) slowJet->removeJet(i);
 622       continue;
 623     }
 624     jetMomenta.push_back(jetMom);
 625   }
 626
 627   // Reverse jetMomenta to restore eT/pT ordering
 628   reverse(jetMomenta.begin(), jetMomenta.end());
 629 }
 630
 631 //--------------------------------------------------------------------------
 632
 633 // Step (2c): veto decision (returning true vetoes the event)
 634
 635 bool JetMatchingAlpgen::matchPartonsToJets(int iType) {
 636
 637   // Use two different routines for light/heavy jets as
 638   // different veto conditions and for clarity
 639   if (iType == 0) return (matchPartonsToJetsLight() > 0);
 640   else            return (matchPartonsToJetsHeavy() > 0);
 641 }
 642
 643 //--------------------------------------------------------------------------
 644
 645 // Step(2c): light jets
 646 // Return codes are given indicating the reason for a veto.
 647 // Although not currently used, they are a useful debugging tool:
 648 //   0 = no veto
 649 //   1 = veto as number of jets less than number of partons
 650 //   2 = veto as exclusive mode and number of jets greater than
 651 //       number of partons
 652 //   3 = veto as inclusive mode and there would be an extra jet
 653 //       that is harder than any matched soft jet
 654 //   4 = veto as there is a parton which does not match a jet
 655
 656 int JetMatchingAlpgen::matchPartonsToJetsLight() {
 657
 658   // Always veto if number of jets is less than original number of jets
 659   if (jetMomenta.size() < typeIdx[0].size()) return LESS_JETS;
 660   // Veto if in exclusive mode and number of jets bigger than original
 661   if (exclusive && jetMomenta.size() > typeIdx[0].size()) return MORE_JETS;
 662
 663   // Sort partons by eT/pT
 664   sortTypeIdx(typeIdx[0]);
 665
 666   // Number of hard partons
 667   int nParton = typeIdx[0].size();
 668
 669   // Keep track of which jets have been assigned a hard parton
 670   vector < bool > jetAssigned;
 671   jetAssigned.assign(jetMomenta.size(), false);
 672
 673   // Jet matching procedure: (1) deltaR between partons and jets
 674   if (jetMatch == 1) {
 675
 676     // Loop over light hard partons and get 4-momentum
 677     for (int i = 0; i < nParton; i++) {
 678       Vec4 p1 = eventProcess[typeIdx[0][i]].p();
 679
 680       // Track which jet has the minimal dR measure with this parton
 681       int    jMin  = -1;
 682       double dRmin = 0.;
 683
 684       // Loop over all jets (skipping those already assigned).
 685       for (int j = 0; j < int(jetMomenta.size()); j++) {
 686         if (jetAssigned[j]) continue;
 687
 688         // DeltaR between parton/jet and store if minimum
 689         double dR = (jetAlgorithm == 1) ?
 690           REtaPhi(p1, jetMomenta[j]) : RRapPhi(p1, jetMomenta[j]);
 691         if (jMin < 0 || dR < dRmin) {
 692           dRmin = dR;
 693           jMin  = j;
 694         }
 695       } // for (j)
 696
 697       // Check for jet-parton match
 698       if (jMin >= 0 && dRmin < coneRadius * coneMatchLight) {
 699
 700         // If the matched jet is not one of the nParton hardest jets,
 701         // the extra left over jet would be harder than some of the
 702         // matched jets. This is disallowed, so veto.
 703         if (jMin >= nParton) return HARD_JET;
 704
 705         // Mark jet as assigned.
 706         jetAssigned[jMin] = true;
 707
 708       // If no match, then event will be vetoed in all cases
 709       } else return UNMATCHED_PARTON;
 710
 711     } // for (i)
 712
 713   // Jet matching procedure: (2) ghost particles in SlowJet
 714   } else {
 715
 716     // Loop over added 'ghost' particles and find if assigned to a jet
 717     for (int i = workEventJet.size() - nParton;
 718         i < workEventJet.size(); i++) {
 719       int jMin = slowJet->jetAssignment(i);
 720
 721       // Veto if:
 722       //  1) not one of nParton hardest jets
 723       //  2) not assigned to a jet
 724       //  3) jet has already been assigned
 725       if (jMin >= nParton)               return HARD_JET;
 726       if (jMin < 0 || jetAssigned[jMin]) return UNMATCHED_PARTON;
 727
 728       // Mark jet as assigned
 729       jetAssigned[jMin] = true;
 730
 731     } // for (i)
 732   } // if (jetMatch)
 733
 734   // Minimal eT/pT (CellJet/SlowJet) of matched light jets. Needed
 735   // later for heavy jet vetos in inclusive mode.
 736   if (nParton > 0)
 737     eTpTlightMin = (jetAlgorithm == 1) ? jetMomenta[nParton - 1].eT()
 738                                        : jetMomenta[nParton - 1].pT();
 739   else
 740     eTpTlightMin = -1.;
 741
 742   // No veto
 743   return NONE;
 744 }
 745
 746 //--------------------------------------------------------------------------
 747
 748 // Step(2c): heavy jets
 749 // Return codes are given indicating the reason for a veto.
 750 // Although not currently used, they are a useful debugging tool:
 751 //   0 = no veto as there are no extra jets present
 752 //   1 = veto as in exclusive mode and extra jets present
 753 //   2 = veto as in inclusive mode and extra jets were harder
 754 //       than any matched light jet
 755
 756 int JetMatchingAlpgen::matchPartonsToJetsHeavy() {
 757
 758   // If there are no extra jets, then accept
 759   if (jetMomenta.empty()) return NONE;
 760
 761   // Number of hard partons
 762   int nParton = typeIdx[1].size();
 763
 764   // Remove jets that are close to heavy quarks
 765   set < int > removeJets;
 766
 767   // Jet matching procedure: (1) deltaR between partons and jets
 768   if (jetMatch == 1) {
 769
 770     // Loop over heavy hard partons and get 4-momentum
 771     for (int i = 0; i < nParton; i++) {
 772       Vec4 p1 = eventProcess[typeIdx[1][i]].p();
 773
 774       // Loop over all jets, find dR and mark for removal if match
 775       for (int j = 0; j < int(jetMomenta.size()); j++) {
 776         double dR = (jetAlgorithm == 1) ?
 777             REtaPhi(p1, jetMomenta[j]) : RRapPhi(p1, jetMomenta[j]);
 778         if (dR < coneRadiusHeavy * coneMatchHeavy)
 779           removeJets.insert(j);
 780
 781       } // for (j)
 782     } // for (i)
 783
 784   // Jet matching procedure: (2) ghost particles in SlowJet
 785   } else {
 786
 787     // Loop over added 'ghost' particles and if assigned to a jet
 788     // then mark this jet for removal
 789     for (int i = workEventJet.size() - nParton;
 790         i < workEventJet.size(); i++) {
 791       int jMin = slowJet->jetAssignment(i);
 792       if (jMin >= 0) removeJets.insert(jMin);
 793     }
 794
 795   }
 796
 797   // Remove jets (backwards order to not disturb indices)
 798   for (set < int >::reverse_iterator it  = removeJets.rbegin();
 799                                      it != removeJets.rend(); it++)
 800     jetMomenta.erase(jetMomenta.begin() + *it);
 801
 802   // Handle case if there are still extra jets
 803   if (!jetMomenta.empty()) {
 804
 805     // Exclusive mode, so immediate veto
 806     if (exclusive) return MORE_JETS;
 807
 808     // Inclusive mode; extra jets must be softer than any matched light jet
 809     else if (eTpTlightMin >= 0.)
 810       for (size_t j = 0; j < jetMomenta.size(); j++) {
 811         // CellJet uses eT, SlowJet uses pT
 812         if ( (jetAlgorithm == 1 && jetMomenta[j].eT() > eTpTlightMin) ||
 813              (jetAlgorithm == 2 && jetMomenta[j].pT() > eTpTlightMin) )
 814           return HARD_JET;
 815       }
 816
 817   } // if (!jetMomenta.empty())
 818
 819   // No extra jets were present so no veto
 820   return NONE;
 821 }
 822
 823 //==========================================================================
 824
 825 // Main implementation of Madgraph UserHooks class.
 826 // This may be split out to a separate C++ file if desired,
 827 // but currently included here for ease of use.
 828
 829 //--------------------------------------------------------------------------
 830
 831 // Initialisation routine automatically called from Pythia::init().
 832 // Setup all parts needed for the merging.
 833
 834 bool JetMatchingMadgraph::initAfterBeams() {
 835
 836   // Read in Madgraph specific configuration variables
 837   bool setMad    = settingsPtr->flag("JetMatching:setMad");
 838
 839   // If Madgraph parameters are present, then parse in MadgraphPar object
 840   MadgraphPar par(infoPtr);
 841   string parStr = infoPtr->header("MGRunCard");
 842   if (!parStr.empty()) {
 843     par.parse(parStr);
 844     par.printParams();
 845   }
 846
 847   // Set Madgraph merging parameters from the file if requested
 848   if (setMad) {
 849     if ( par.haveParam("xqcut")    && par.haveParam("maxjetflavor")
 850       && par.haveParam("alpsfact") && par.haveParam("ickkw") ) {
 851       settingsPtr->flag("JetMatching:merge", par.getParam("ickkw"));
 852       settingsPtr->parm("JetMatching:qCut", par.getParam("xqcut"));
 853       settingsPtr->mode("JetMatching:nQmatch",
 854         par.getParamAsInt("maxjetflavor"));
 855       settingsPtr->parm("JetMatching:clFact",
 856         clFact = par.getParam("alpsfact"));
 857       if (par.getParamAsInt("ickkw") == 0)
 858         infoPtr->errorMsg("Error in JetMatchingMadgraph:init: "
 859           "Madgraph file parameters are not set for merging");
 860
 861     // Warn if setMad requested, but one or more parameters not present
 862     } else {
 863        infoPtr->errorMsg("Warning in JetMatchingMadgraph:init: "
 864           "Madgraph merging parameters not found");
 865        if (!par.haveParam("xqcut")) infoPtr->errorMsg("Warning in "
 866           "JetMatchingMadgraph:init: No xqcut");
 867        if (!par.haveParam("ickkw")) infoPtr->errorMsg("Warning in "
 868           "JetMatchingMadgraph:init: No ickkw");
 869        if (!par.haveParam("maxjetflavor")) infoPtr->errorMsg("Warning in "
 870           "JetMatchingMadgraph:init: No maxjetflavor");
 871        if (!par.haveParam("alpsfact")) infoPtr->errorMsg("Warning in "
 872           "JetMatchingMadgraph:init: No alpsfact");
 873     }
 874   }
 875
 876   // Read in Madgraph merging parameters
 877   doMerge      = settingsPtr->flag("JetMatching:merge");
 878   qCut         = settingsPtr->parm("JetMatching:qCut");
 879   nQmatch      = settingsPtr->mode("JetMatching:nQmatch");
 880   clFact       = settingsPtr->parm("JetMatching:clFact");
 881
 882   // Read in jet algorithm parameters
 883   jetAlgorithm   = settingsPtr->mode("JetMatching:jetAlgorithm");
 884   nJetMax        = settingsPtr->mode("JetMatching:nJetMax");
 885   eTjetMin       = settingsPtr->parm("JetMatching:eTjetMin");
 886   coneRadius     = settingsPtr->parm("JetMatching:coneRadius");
 887   etaJetMax      = settingsPtr->parm("JetMatching:etaJetMax");
 888   slowJetPower   = settingsPtr->mode("JetMatching:slowJetPower");
 889
 890   // Matching procedure
 891   jetAllow       = settingsPtr->mode("JetMatching:jetAllow");
 892   exclusiveMode  = settingsPtr->mode("JetMatching:exclusive");
 893   qCutSq         = pow(qCut,2);
 894   etaJetMaxAlgo  = etaJetMax;
 895
 896   // If not merging, then done
 897   if (!doMerge) return true;
 898
 899   // Exclusive mode; if set to 2, then set based on nJet/nJetMax
 900   if (exclusiveMode == 2) {
 901
 902     // No nJet or nJetMax, so default to exclusive mode
 903     if (nJetMax < 0) {
 904       infoPtr->errorMsg("Warning in JetMatchingMadgraph:init: "
 905         "missing jet multiplicity information; running in exclusive mode");
 906       exclusiveMode = 1;
 907     }
 908   }
 909
 910   // Initialise chosen jet algorithm.
 911   // Currently, this only supports the kT-algorithm in SlowJet.
 912   // Use the QCD distance measure by default.
 913   jetAlgorithm = 2;
 914   slowJetPower = 1;
 915   slowJet = new SlowJet(slowJetPower, coneRadius, eTjetMin,
 916     etaJetMaxAlgo, 2, 2, NULL, false);
 917
 918   // Setup local event records
 919   eventProcessOrig.init("(eventProcessOrig)", particleDataPtr);
 920   eventProcess.init("(eventProcess)", particleDataPtr);
 921   workEventJet.init("(workEventJet)", particleDataPtr);
 922
 923   // Print information
 924   string jetStr  = (jetAlgorithm ==  1) ? "CellJet" :
 925                    (slowJetPower == -1) ? "anti-kT" :
 926                    (slowJetPower ==  0) ? "C/A"     :
 927                    (slowJetPower ==  1) ? "kT"      : "unknown";
 928   string modeStr = (exclusiveMode)         ? "exclusive" : "inclusive";
 929   cout << endl
 930        << " *-----  Madgraph matching parameters  -----*" << endl
 931        << " |  qCut                |  " << setw(14)
 932        << qCut << "  |" << endl
 933        << " |  nQmatch             |  " << setw(14)
 934        << nQmatch << "  |" << endl
 935        << " |  clFact              |  " << setw(14)
 936        << clFact << "  |" << endl
 937        << " |  Jet algorithm       |  " << setw(14)
 938        << jetStr << "  |" << endl
 939        << " |  eTjetMin            |  " << setw(14)
 940        << eTjetMin << "  |" << endl
 941        << " |  etaJetMax           |  " << setw(14)
 942        << etaJetMax << "  |" << endl
 943        << " |  jetAllow            |  " << setw(14)
 944        << jetAllow << "  |" << endl
 945        << " |  Mode                |  " << setw(14)
 946        << modeStr << "  |" << endl
 947        << " *-----------------------------------------*" << endl;
 948
 949   return true;
 950 }
 951
 952 //--------------------------------------------------------------------------
 953
 954 // Step (1): sort the incoming particles
 955
 956 void JetMatchingMadgraph::sortIncomingProcess(const Event &event) {
 957
 958   // Remove resonance decays from original process and keep only final
 959   // state. Resonances will have positive status code after this step.
 960   omitResonanceDecays(eventProcessOrig, true);
 961   eventProcess = workEvent;
 962
 963   // Sort original process final state into light/heavy jets and 'other'.
 964   // Criteria:
 965   //   1 <= ID <= nQmatch, or ID == 21         --> light jet (typeIdx[0])
 966   //   nQMatch < ID                            --> heavy jet (typeIdx[1])
 967   //   All else                                --> other     (typeIdx[2])
 968   // Note that 'typeIdx' stores indices into 'eventProcess' (after resonance
 969   // decays are omitted), while 'typeSet' stores indices into the original
 970   // process record, 'eventProcessOrig', but these indices are also valid
 971   // in 'event'.
 972   for (int i = 0; i < 3; i++) {
 973     typeIdx[i].clear();
 974     typeSet[i].clear();
 975   }
 976   for (int i = 0; i < eventProcess.size(); i++) {
 977     // Ignore non-final state and default to 'other'
 978     if (!eventProcess[i].isFinal()) continue;
 979     int idx = 2;
 980
 981     // Light jets: all gluons and quarks with id less than or equal to nQmatch
 982     if (eventProcess[i].id() == ID_GLUON
 983       || (eventProcess[i].idAbs() <= nQmatch) ) idx = 0;
 984
 985     // Heavy jets:  all quarks with id greater than nQmatch
 986     else if (eventProcess[i].idAbs() > nQmatch
 987       && eventProcess[i].idAbs() <= ID_TOP) idx = 1;
 988
 989     // Store
 990     typeIdx[idx].push_back(i);
 991     typeSet[idx].insert(eventProcess[i].daughter1());
 992   }
 993
 994   // Exclusive mode; if set to 2, then set based on nJet/nJetMax
 995   if (exclusiveMode == 2) {
 996
 997     // Inclusive if nJet == nJetMax, exclusive otherwise
 998     int nParton = typeIdx[0].size();
 999     exclusive = (nParton == nJetMax) ? false : true;
1000
1001   // Otherwise, just set as given
1002   } else {
1003     exclusive = (exclusiveMode == 0) ? false : true;
1004   }
1005
1006   // Extract partons from hardest subsystem + ISR + FSR only into
1007   // workEvent. Note no resonance showers or MPIs.
1008   subEvent(event);
1009 }
1010
1011 //--------------------------------------------------------------------------
1012
1013 // Step (2a): pick which particles to pass to the jet algorithm
1014
1015 void JetMatchingMadgraph::jetAlgorithmInput(const Event &event, int iType) {
1016
1017   // Take input from 'workEvent' and put output in 'workEventJet'
1018   workEventJet = workEvent;
1019
1020   // Loop over particles and decide what to pass to the jet algorithm
1021   for (int i = 0; i < workEventJet.size(); ++i) {
1022     if (!workEventJet[i].isFinal()) continue;
1023
1024     // jetAllow option to disallow certain particle types
1025     if (jetAllow == 1) {
1026
1027       // Original AG+Py6 algorithm explicitly excludes tops,
1028       // leptons and photons.
1029       int id = workEventJet[i].idAbs();
1030       if ((id >= ID_LEPMIN && id <= ID_LEPMAX) || id == ID_TOP
1031       || id == ID_PHOTON) {
1032         workEventJet[i].statusNeg();
1033         continue;
1034       }
1035     }
1036
1037     // Get the index of this particle in original event
1038     int idx = workEventJet[i].daughter1();
1039
1040     // Start with particle idx, and afterwards track mothers
1041     while (true) {
1042
1043       // Light jets
1044       if (iType == 0) {
1045
1046         // Do not include if originates from heavy jet or 'other'
1047         if (typeSet[1].find(idx) != typeSet[1].end() ||
1048             typeSet[2].find(idx) != typeSet[2].end()) {
1049           workEventJet[i].statusNeg();
1050           break;
1051         }
1052
1053         // Made it to start of event record so done
1054         if (idx == 0) break;
1055         // Otherwise next mother and continue
1056         idx = event[idx].mother1();
1057
1058       // Heavy jets
1059       } else if (iType == 1) {
1060
1061         // Only include if originates from heavy jet
1062         if (typeSet[1].find(idx) != typeSet[1].end()) break;
1063
1064         // Made it to start of event record with no heavy jet mother,
1065         // so DO NOT include particle
1066         if (idx == 0) {
1067           workEventJet[i].statusNeg();
1068           break;
1069         }
1070
1071         // Otherwise next mother and continue
1072         idx = event[idx].mother1();
1073
1074       } // if (iType)
1075     } // while (true)
1076   } // for (i)
1077 }
1078
1079 //--------------------------------------------------------------------------
1080
1081 // Step (2b): run jet algorithm and provide common output
1082 // This does nothing, because the jet algorithm is run several times
1083 //  in the matching algorithm.
1084
1085 void JetMatchingMadgraph::runJetAlgorithm() {; }
1086
1087 //--------------------------------------------------------------------------
1088
1089 // Step (2c): veto decision (returning true vetoes the event)
1090
1091 bool JetMatchingMadgraph::matchPartonsToJets(int iType) {
1092
1093   // Use two different routines for light/heavy jets as
1094   // different veto conditions and for clarity
1095   if (iType == 0) return (matchPartonsToJetsLight() > 0);
1096   else            return (matchPartonsToJetsHeavy() > 0);
1097 }
1098
1099 //--------------------------------------------------------------------------
1100
1101 // Step(2c): light jets
1102 // Return codes are given indicating the reason for a veto.
1103 // Although not currently used, they are a useful debugging tool:
1104 //   0 = no veto
1105 //   1 = veto as number of jets less than number of partons
1106 //   2 = veto as exclusive mode and number of jets greater than
1107 //       number of partons
1108 //   3 = veto as inclusive mode and there would be an extra jet
1109 //       that is harder than any matched soft jet
1110 //   4 = veto as there is a parton which does not match a jet
1111
1112 int JetMatchingMadgraph::matchPartonsToJetsLight() {
1113
1114   // Count the number of hard partons
1115   int nParton = typeIdx[0].size();
1116
1117   // Initialize SlowJet with current working event
1118   if (!slowJet->setup(workEventJet) ) {
1119     infoPtr->errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
1120       "Light: the SlowJet algorithm failed on setup");
1121     return NONE;
1122   }
1123   double localQcutSq = qCutSq;
1124   double dOld = 0.0;
1125   // Cluster in steps to find all hadronic jets at the scale qCut
1126   while ( slowJet->sizeAll() - slowJet->sizeJet() > 0 ) {
1127     if( slowJet->dNext() > localQcutSq ) break;
1128     dOld = slowJet->dNext();
1129     slowJet->doStep();
1130   }
1131   int nJets = slowJet->sizeJet();
1132   int nClus = slowJet->sizeAll();
1133
1134   // Debug printout.
1135   if (MATCHINGDEBUG) slowJet->list(true);
1136
1137   // Count of the number of hadronic jets in SlowJet accounting
1138   int nCLjets = nClus - nJets;
1139   // Compare number of hadronic jets to number of partons
1140   // Veto event if too few hadronic jets
1141   if ( nCLjets < nParton ) return LESS_JETS;
1142   // In exclusive mode, do not allow more hadronic jets than partons
1143   if ( exclusive ) {
1144     if ( nCLjets > nParton ) return MORE_JETS;
1145
1146   // In inclusive mode, there can be more hadronic jets than partons,
1147   //  provided that all partons have a matching hadronic jet
1148   } else {
1149     if (!slowJet->setup(workEventJet) ) {
1150       infoPtr->errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
1151         "Light: the SlowJet algorithm failed on setup");
1152       return NONE;
1153     }
1154     // Cluster into hadronic jets until there are the same number as partons
1155     while ( slowJet->sizeAll() - slowJet->sizeJet() > nParton )
1156       slowJet->doStep();
1157     // Sort partons in pT.  Update local qCut value.
1158     //  Hadronic jets are already sorted in pT.
1159     localQcutSq = dOld;
1160     if ( clFact >= 0. && nParton > 0 ) {
1161        vector<double> partonPt;
1162        for (int i = 0; i < nParton; ++i)
1163          partonPt.push_back( eventProcess[typeIdx[0][i]].pT2() );
1164        sort( partonPt.begin(), partonPt.end());
1165        localQcutSq = max( qCutSq, partonPt[0]);
1166     }
1167     nJets = slowJet->sizeJet();
1168     nClus = slowJet->sizeAll();
1169   }
1170   // Update scale if clustering factor is non-zero
1171   if ( clFact != 0. ) localQcutSq *= pow2(clFact);
1172
1173   Event tempEvent;
1174   tempEvent.init( "(tempEvent)", particleDataPtr);
1175   int nPass = 0;
1176   double pTminEstimate = -1.;
1177   // Construct a master copy of the event containing only the
1178   // hardest nParton hadronic clusters. While constructing the event,
1179   // the parton type (ID_GLUON) and status (98,99) are arbitrary.
1180   for (int i = nJets; i < nClus; ++i) {
1181     tempEvent.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0, slowJet->p(i).px(),
1182       slowJet->p(i).py(), slowJet->p(i).pz(), slowJet->p(i).e() );
1183     ++nPass;
1184     pTminEstimate = max( pTminEstimate, slowJet->pT(i));
1185     if(nPass == nParton) break;
1186   }
1187   size_t tempSize = tempEvent.size();
1188   // This keeps track of which hadronic jets are matched to parton
1189   vector<bool> jetAssigned;
1190   jetAssigned.assign( tempSize, false);
1191
1192   // Take the list of unmatched hadronic jets and append a parton, one at
1193   // a time. The parton will be clustered with the "closest" hadronic jet
1194   // or the beam if the distance measure is below the cut. When a hadronic
1195   // jet is matched to a parton, it is removed from the list of unmatched
1196   // hadronic jets. This process continues until all hadronic jets are
1197   // matched to partons or it is not possible to make a match.
1198   int iNow = 0;
1199   while ( iNow < nParton ) {
1200     Event tempEventJet;
1201     tempEventJet.init("(tempEventJet)", particleDataPtr);
1202     for (size_t i = 0; i < tempSize; ++i) {
1203       if (jetAssigned[i]) continue;
1204       Vec4 pIn = tempEvent[i].p();
1205       // Append unmatched hadronic jets
1206       tempEventJet.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0,
1207         pIn.px(), pIn.py(), pIn.pz(), pIn.e() );
1208     }
1209     Vec4 pIn = eventProcess[typeIdx[0][iNow]].p();
1210     // Append the current parton
1211     tempEventJet.append( ID_GLUON, 99, 0, 0, 0, 0, 0, 0,
1212       pIn.px(), pIn.py(), pIn.pz(), pIn.e() );
1213     if ( !slowJet->setup(tempEventJet) ) {
1214       infoPtr->errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
1215         "Light: the SlowJet algorithm failed on setup");
1216     } return NONE;
1217     // These are the conditions for the hadronic jet to match the parton
1218     //  at the local qCut scale
1219     if ( slowJet->iNext() == tempEventJet.size() - 1
1220       && slowJet->jNext() > -1 && slowJet->dNext() < localQcutSq ) {
1221       int iKnt = -1;
1222       for (size_t i = 0; i != tempSize; ++i) {
1223         if (jetAssigned[i]) continue;
1224         ++iKnt;
1225         // Identify the hadronic jet that matches the parton
1226         if (iKnt == slowJet->jNext() ) jetAssigned[i] = true;
1227       }
1228     } else { return UNMATCHED_PARTON; }
1229     ++iNow;
1230   }
1231
1232   // Minimal eT/pT (CellJet/SlowJet) of matched light jets.
1233   // Needed later for heavy jet vetos in inclusive mode.
1234   // This information is not used currently.
1235   if (nParton > 0 && pTminEstimate > 0) eTpTlightMin = pTminEstimate;
1236   else eTpTlightMin = -1.;
1237
1238   // No veto
1239   return NONE;
1240 }
1241
1242 //--------------------------------------------------------------------------
1243
1244 // Step(2c): heavy jets
1245 // Return codes are given indicating the reason for a veto.
1246 // Although not currently used, they are a useful debugging tool:
1247 //   0 = no veto as there are no extra jets present
1248 //   1 = veto as in exclusive mode and extra jets present
1249 //   2 = veto as in inclusive mode and extra jets were harder
1250 //       than any matched light jet
1251
1252 int JetMatchingMadgraph::matchPartonsToJetsHeavy() {
1253
1254   // Currently, heavy jets are unmatched
1255   // If there are no extra jets, then accept
1256   if (jetMomenta.empty()) return NONE;
1257
1258   // No extra jets were present so no veto
1259   return NONE;
1260 }
1261
1262 //==========================================================================
1263
1264 #endif // end Pythia8_JetMatching_H