--- /dev/null
+// MLMhooks.h is a part of the PYTHIA event generator.
+// Copyright (C) 2012 Torbjorn Sjostrand.
+// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
+// Please respect the MCnet Guidelines, see GUIDELINES for details.
+
+// Author: Richard Corke (richard.corke@thep.lu.se)
+// This file provides the MLMhooks class to perform MLM merging.
+// Example usage is shown in main32.cc, and further details
+// can be found in the 'Alpgen and MLM Merging' manual page.
+
+#ifndef MLMHOOKS_H
+#define MLMHOOKS_H
+
+// Includes
+#include "Pythia.h"
+using namespace Pythia8;
+
+//==========================================================================
+
+// Preprocessor settings
+
+// Debug flag to give verbose information at all stages of the merging.
+//#define MLMDEBUG
+// Debug flag to enable some extra checks in the code
+//#define MLMCHECK
+
+//==========================================================================
+
+// Declaration of main MLMhooks class to perform MLM matching.
+// Note that it is defined with virtual inheritance, so that it can
+// be combined with other UserHooks classes, see e.g. main32.cc.
+
+class MLMhooks : virtual public UserHooks {
+
+public:
+
+ // Constructor and destructor
+ MLMhooks() : cellJet(NULL), slowJet(NULL) {}
+ ~MLMhooks() {
+ if (cellJet) delete cellJet;
+ if (slowJet) delete slowJet;
+ }
+
+ // Initialisation
+ virtual bool initAfterBeams();
+
+ // Process level vetos
+ virtual bool canVetoProcessLevel();
+ virtual bool doVetoProcessLevel(Event &);
+
+ // Parton level vetos (before beam remnants and resonance decays)
+ virtual bool canVetoPartonLevelEarly();
+ virtual bool doVetoPartonLevelEarly(const Event &);
+
+private:
+
+ // Different steps of the MLM matching algorithm
+ void sortIncomingProcess(const Event &);
+ void jetAlgorithmInput(const Event &, int);
+ void runJetAlgorithm();
+ bool matchPartonsToJets(int);
+ int matchPartonsToJetsLight();
+ int matchPartonsToJetsHeavy();
+
+ // DeltaR between two 4-vectors (eta and y variants)
+ inline double Vec4eta(const Vec4 &pIn) {
+ return -log(tan(pIn.theta() / 2.));
+ }
+ inline double Vec4y(const Vec4 &pIn) {
+ return 0.5 * log((pIn.e() + pIn.pz()) / (pIn.e() - pIn.pz()));
+ }
+ inline double deltaReta(const Vec4 &p1, const Vec4 &p2) {
+ double dEta = abs(Vec4eta(p1) - Vec4eta(p2));
+ double dPhi = abs(p1.phi() - p2.phi());
+ if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
+ return sqrt(dEta*dEta + dPhi*dPhi);
+ }
+ inline double deltaRy(const Vec4 &p1, const Vec4 &p2) {
+ double dy = abs(Vec4y(p1) - Vec4y(p2));
+ double dPhi = abs(p1.phi() - p2.phi());
+ if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
+ return sqrt(dy*dy + dPhi*dPhi);
+ }
+
+ // Function to sort typeIdx vectors into descending eT/pT order.
+ // Uses a selection sort, as number of partons generally small
+ // and so efficiency not a worry.
+ void sortTypeIdx(vector < int > &vecIn) {
+ for (size_t i = 0; i < vecIn.size(); i++) {
+ size_t jMax = i;
+ double vMax = (jetAlgorithm == 1) ?
+ eventProcess[vecIn[i]].eT() :
+ eventProcess[vecIn[i]].pT();
+ for (size_t j = i + 1; j < vecIn.size(); j++) {
+ double vNow = (jetAlgorithm == 1) ?
+ eventProcess[vecIn[j]].eT() :
+ eventProcess[vecIn[j]].pT();
+ if (vNow > vMax) {
+ vMax = vNow;
+ jMax = j;
+ }
+ }
+ if (jMax != i) swap(vecIn[i], vecIn[jMax]);
+ }
+ }
+
+ // Master switch for merging
+ bool doMerge;
+
+ // Maximum and current number of jets
+ int nJetMax, nJet;
+
+ // Jet algorithm parameters
+ int jetAlgorithm;
+ double eTjetMin, coneRadius, etaJetMax, etaJetMaxAlgo;
+
+ // CellJet specific
+ int nEta, nPhi;
+ double eTseed, eTthreshold;
+
+ // SlowJet specific
+ int slowJetPower;
+
+ // Merging procedure parameters
+ int jetAllow, jetMatch, exclusiveMode;
+ double coneMatchLight, coneRadiusHeavy, coneMatchHeavy;
+ bool exclusive;
+
+ // Event records to store original incoming process, final-state of the
+ // incoming process and what will be passed to the jet algorithm.
+ // Not completely necessary to store all steps, but makes tracking the
+ // steps of the algorithm a lot easier.
+ Event eventProcessOrig, eventProcess, workEventJet;
+
+ // Internal jet algorithms
+ CellJet *cellJet;
+ SlowJet *slowJet;
+
+ // Sort final-state of incoming process into light/heavy jets and 'other'
+ vector < int > typeIdx[3];
+ set < int > typeSet[3];
+
+ // Momenta output of jet algorithm (to provide same output regardless of
+ // the selected jet algorithm)
+ vector < Vec4 > jetMomenta;
+
+ // Store the minimum eT/pT of matched light jets
+ double eTpTlightMin;
+
+ // Constants
+ static const double GHOSTENERGY, ZEROTHRESHOLD;
+};
+
+//--------------------------------------------------------------------------
+
+// Main implementation of MLMhooks class. This may be split out to a
+// separate C++ file if desired, but currently included here for ease
+// of use.
+
+// Constants: could be changed here if desired, but normally should not.
+// These are of technical nature, as described for each.
+
+// The energy of ghost particles. For technical reasons, this cannot be
+// set arbitrarily low, see 'Particle::TINY' in 'Event.cc' for details.
+const double MLMhooks::GHOSTENERGY = 1e-15;
+
+// A zero threshold value for double comparisons.
+const double MLMhooks::ZEROTHRESHOLD = 1e-10;
+
+//--------------------------------------------------------------------------
+
+// Initialisation routine automatically called from Pythia::init().
+// Setup all parts needed for the merging.
+
+bool MLMhooks::initAfterBeams() {
+
+ // Read in parameters
+ doMerge = settingsPtr->flag("MLM:merge");
+ nJet = settingsPtr->mode("MLM:nJet");
+ nJetMax = settingsPtr->mode("MLM:nJetMax");
+ jetAlgorithm = settingsPtr->mode("MLM:jetAlgorithm");
+ eTjetMin = settingsPtr->parm("MLM:eTjetMin");
+ coneRadius = settingsPtr->parm("MLM:coneRadius");
+ etaJetMax = settingsPtr->parm("MLM:etaJetMax");
+ // Use etaJetMax + coneRadius in input to jet algorithms
+ etaJetMaxAlgo = etaJetMax + coneRadius;
+ // CellJet specific
+ nEta = settingsPtr->mode("MLM:nEta");
+ nPhi = settingsPtr->mode("MLM:nPhi");
+ eTseed = settingsPtr->parm("MLM:eTseed");
+ eTthreshold = settingsPtr->parm("MLM:eTthreshold");
+ // SlowJet specific
+ slowJetPower = settingsPtr->mode("MLM:slowJetPower");
+ // Matching procedure
+ jetAllow = settingsPtr->mode("MLM:jetAllow");
+ jetMatch = settingsPtr->mode("MLM:jetMatch");
+ coneMatchLight = settingsPtr->parm("MLM:coneMatchLight");
+ coneRadiusHeavy = settingsPtr->parm("MLM:coneRadiusHeavy");
+ if (coneRadiusHeavy < 0.) coneRadiusHeavy = coneRadius;
+ coneMatchHeavy = settingsPtr->parm("MLM:coneMatchHeavy");
+ exclusiveMode = settingsPtr->mode("MLM:exclusive");
+
+ // If not merging, then done
+ if (!doMerge) return true;
+
+ // Exclusive mode; if set to 2, then set based on nJet/nJetMax
+ if (exclusiveMode == 2) {
+
+ // No nJet or nJetMax, so default to exclusive mode
+ if (nJet < 0 || nJetMax < 0) {
+ infoPtr->errorMsg("Warning in MLMhooks:init: "
+ "missing jet multiplicity information; running in exclusive mode");
+ exclusive = true;
+
+ // Inclusive if nJet == nJetMax, exclusive otherwise
+ } else {
+ exclusive = (nJet == nJetMax) ? false : true;
+ }
+
+ // Otherwise, just set as given
+ } else {
+ exclusive = (exclusiveMode == 0) ? false : true;
+ }
+
+ // Initialise chosen jet algorithm. CellJet.
+ if (jetAlgorithm == 1) {
+
+ // Extra options for CellJet. nSel = 1 means that all final-state
+ // particles are taken and we retain control of what to select.
+ // smear/resolution/upperCut are not used and are set to default values.
+ int nSel = 2, smear = 0;
+ double resolution = 0.5, upperCut = 2.;
+ cellJet = new CellJet(etaJetMaxAlgo, nEta, nPhi, nSel,
+ smear, resolution, upperCut, eTthreshold);
+
+ // SlowJet
+ } else if (jetAlgorithm == 2) {
+ slowJet = new SlowJet(slowJetPower, coneRadius, eTjetMin, etaJetMaxAlgo);
+ }
+
+ // Check the jetMatch parameter; option 2 only works with SlowJet
+ if (jetAlgorithm == 1 && jetMatch == 2) {
+ infoPtr->errorMsg("Warning in MLMhooks:init: "
+ "jetMatch = 2 only valid with SlowJet algorithm. "
+ "Reverting to jetMatch = 1.");
+ jetMatch = 1;
+ }
+
+ // Setup local event records
+ eventProcessOrig.init("(eventProcessOrig)", particleDataPtr);
+ eventProcess.init("(eventProcess)", particleDataPtr);
+ workEventJet.init("(workEventJet)", particleDataPtr);
+
+ // Print information
+ string jetStr = (jetAlgorithm == 1) ? "CellJet" :
+ (slowJetPower == -1) ? "anti-kT" :
+ (slowJetPower == 0) ? "C/A" :
+ (slowJetPower == 1) ? "kT" : "unknown";
+ string modeStr = (exclusive) ? "exclusive" : "inclusive";
+ stringstream nJetStr, nJetMaxStr;
+ if (nJet >= 0) nJetStr << nJet; else nJetStr << "unknown";
+ if (nJetMax >= 0) nJetMaxStr << nJetMax; else nJetMaxStr << "unknown";
+ cout << endl
+ << " *------- MLM matching parameters -------*" << endl
+ << " | nJet | " << setw(14)
+ << nJetStr.str() << " |" << endl
+ << " | nJetMax | " << setw(14)
+ << nJetMaxStr.str() << " |" << endl
+ << " | Jet algorithm | " << setw(14)
+ << jetStr << " |" << endl
+ << " | eTjetMin | " << setw(14)
+ << eTjetMin << " |" << endl
+ << " | coneRadius | " << setw(14)
+ << coneRadius << " |" << endl
+ << " | etaJetMax | " << setw(14)
+ << etaJetMax << " |" << endl
+ << " | jetAllow | " << setw(14)
+ << jetAllow << " |" << endl
+ << " | jetMatch | " << setw(14)
+ << jetMatch << " |" << endl
+ << " | coneMatchLight | " << setw(14)
+ << coneMatchLight << " |" << endl
+ << " | coneRadiusHeavy | " << setw(14)
+ << coneRadiusHeavy << " |" << endl
+ << " | coneMatchHeavy | " << setw(14)
+ << coneMatchHeavy << " |" << endl
+ << " | Mode | " << setw(14)
+ << modeStr << " |" << endl
+ << " *-----------------------------------------*" << endl;
+
+ return true;
+}
+
+//--------------------------------------------------------------------------
+
+// Process level veto. Stores incoming event for later.
+
+bool MLMhooks::canVetoProcessLevel() { return doMerge; }
+
+bool MLMhooks::doVetoProcessLevel(Event& process) {
+
+ // Copy incoming process
+ eventProcessOrig = process;
+ return false;
+}
+
+//--------------------------------------------------------------------------
+
+// Early parton level veto (before beam remnants and resonance showers)
+
+bool MLMhooks::canVetoPartonLevelEarly() { return doMerge; }
+
+bool MLMhooks::doVetoPartonLevelEarly(const Event& event) {
+
+ // 1) Sort the original incoming process. After this step is performed,
+ // the following assignments have been made:
+ // eventProcessOrig - the original incoming process
+ // eventProcess - the final-state of the incoming process with
+ // resonance decays removed (and resonances
+ // themselves now with positive status code)
+ // typeIdx[0/1/2] - indices into 'eventProcess' of
+ // light jets/heavy jets/other
+ // typeSet[0/1/2] - indices into 'event' of
+ // light jets/heavy jets/other
+ // workEvent - partons from the hardest subsystem
+ // + ISR + FSR only
+ sortIncomingProcess(event);
+
+ // Debug
+#ifdef MLMDEBUG
+ // Begin
+ cout << endl << "---------- Begin MLM Debug ----------" << endl;
+
+ // Original incoming process
+ cout << endl << "Original incoming process:";
+ eventProcessOrig.list();
+
+ // Final-state of original incoming process
+ cout << endl << "Final-state incoming process:";
+ eventProcess.list();
+
+ // List categories of sorted particles
+ for (size_t i = 0; i < typeIdx[0].size(); i++)
+ cout << ((i == 0) ? "Light jets: " : ", ")
+ << setw(3) << typeIdx[0][i];
+ for (size_t i = 0; i < typeIdx[1].size(); i++)
+ cout << ((i == 0) ? "\nHeavy jets: " : ", ")
+ << setw(3) << typeIdx[1][i];
+ for (size_t i = 0; i < typeIdx[2].size(); i++)
+ cout << ((i == 0) ? "\nOther: " : ", ")
+ << setw(3) << typeIdx[2][i];
+
+ // Full event at this stage
+ cout << endl << endl << "Event:";
+ event.list();
+
+ // Work event (partons from hardest subsystem + ISR + FSR)
+ cout << endl << "Work event:";
+ workEvent.list();
+#endif
+
+ // 2) Light/heavy jets: iType = 0 (light jets), 1 (heavy jets)
+ int iTypeEnd = (typeIdx[1].empty()) ? 1 : 2;
+ for (int iType = 0; iType < iTypeEnd; iType++) {
+
+ // 2a) Find particles which will be passed from the jet algorithm.
+ // Input from 'workEvent' and output in 'workEventJet'.
+ jetAlgorithmInput(event, iType);
+
+ // Debug
+#ifdef MLMDEBUG
+ // Jet algorithm event
+ cout << endl << "Jet algorithm event (iType = " << iType << "):";
+ workEventJet.list();
+#endif
+
+ // 2b) Run jet algorithm on 'workEventJet'.
+ // Output is stored in jetMomenta.
+ runJetAlgorithm();
+
+ // 2c) Match partons to jets and decide if veto is necessary
+ if (matchPartonsToJets(iType) == true) {
+ // Debug
+#ifdef MLMDEBUG
+ cout << endl << "Event vetoed" << endl
+ << "---------- End MLM Debug ----------" << endl;
+#endif
+ return true;
+ }
+ }
+
+ // Debug
+#ifdef MLMDEBUG
+ cout << endl << "Event accepted" << endl
+ << "---------- End MLM Debug ----------" << endl;
+#endif
+
+ // If we reached here, then no veto
+ return false;
+}
+
+//--------------------------------------------------------------------------
+
+// Step (1): sort the incoming particles
+
+void MLMhooks::sortIncomingProcess(const Event &event) {
+
+ // Remove resonance decays from original process and keep only final
+ // state. Resonances will have positive status code after this step.
+ omitResonanceDecays(eventProcessOrig, true);
+ eventProcess = workEvent;
+
+ // Sort original process final state into light/heavy jets and 'other'.
+ // Criteria:
+ // 1 <= ID <= 5 and massless, or ID == 21 --> light jet (typeIdx[0])
+ // 4 <= ID <= 6 and massive --> heavy jet (typeIdx[1])
+ // All else --> other (typeIdx[2])
+ // Note that 'typeIdx' stores indices into 'eventProcess' (after resonance
+ // decays are omitted), while 'typeSet' stores indices into the original
+ // process record, 'eventProcessOrig', but these indices are also valid
+ // in 'event'.
+ for (int i = 0; i < 3; i++) {
+ typeIdx[i].clear();
+ typeSet[i].clear();
+ }
+ for (int i = 0; i < eventProcess.size(); i++) {
+ // Ignore nonfinal and default to 'other'
+ if (!eventProcess[i].isFinal()) continue;
+ int idx = 2;
+
+ // Light jets
+ if (eventProcess[i].id() == 21 || (eventProcess[i].idAbs() <= 5 &&
+ abs(eventProcess[i].m()) < ZEROTHRESHOLD))
+ idx = 0;
+
+ // Heavy jets
+ else if (eventProcess[i].idAbs() >= 4 && eventProcess[i].idAbs() <= 6)
+ idx = 1;
+
+ // Store
+ typeIdx[idx].push_back(i);
+ typeSet[idx].insert(eventProcess[i].daughter1());
+ }
+
+ // Extract partons from hardest subsystem + ISR + FSR only into
+ // workEvent. Note no resonance showers or MPIs.
+ subEvent(event);
+}
+
+//--------------------------------------------------------------------------
+
+// Step (2a): pick which particles to pass to the jet algorithm
+
+void MLMhooks::jetAlgorithmInput(const Event &event, int iType) {
+
+ // Take input from 'workEvent' and put output in 'workEventJet'
+ workEventJet = workEvent;
+
+ // Loop over particles and decide what to pass to the jet algorithm
+ for (int i = 0; i < workEventJet.size(); ++i) {
+ if (!workEventJet[i].isFinal()) continue;
+
+ // jetAllow option to disallow certain particle types
+ if (jetAllow == 1) {
+
+ // Original AG+Py6 algorithm explicitly excludes tops,
+ // leptons and photons.
+ int id = workEventJet[i].idAbs();
+ if ((id >= 11 && id <= 16) || id == 6 || id == 22) {
+ workEventJet[i].statusNeg();
+ continue;
+ }
+ }
+
+ // Get the index of this particle in original event
+ int idx = workEventJet[i].daughter1();
+
+ // Start with particle idx, and afterwards track mothers
+ while (true) {
+
+ // Light jets
+ if (iType == 0) {
+
+ // Do not include if originates from heavy jet or 'other'
+ if (typeSet[1].find(idx) != typeSet[1].end() ||
+ typeSet[2].find(idx) != typeSet[2].end()) {
+ workEventJet[i].statusNeg();
+ break;
+ }
+
+ // Made it to start of event record so done
+ if (idx == 0) break;
+ // Otherwise next mother and continue
+ idx = event[idx].mother1();
+
+ // Heavy jets
+ } else if (iType == 1) {
+
+ // Only include if originates from heavy jet
+ if (typeSet[1].find(idx) != typeSet[1].end()) break;
+
+ // Made it to start of event record with no heavy jet mother,
+ // so DO NOT include particle
+ if (idx == 0) {
+ workEventJet[i].statusNeg();
+ break;
+ }
+
+ // Otherwise next mother and continue
+ idx = event[idx].mother1();
+
+ } // if (iType)
+ } // while (true)
+ } // for (i)
+
+ // For jetMatch = 2, insert ghost particles corresponding to
+ // each hard parton in the original process
+ if (jetMatch == 2) {
+ for (int i = 0; i < int(typeIdx[iType].size()); i++) {
+ // Get y/phi of the parton
+ Vec4 pIn = eventProcess[typeIdx[iType][i]].p();
+ double y = Vec4y(pIn);
+ double phi = pIn.phi();
+
+ // Create a ghost particle and add to the workEventJet
+ double e = GHOSTENERGY;
+ double e2y = exp(2. * y);
+ double pz = e * (e2y - 1.) / (e2y + 1.);
+ double pt = sqrt(e*e - pz*pz);
+ double px = pt * cos(phi);
+ double py = pt * sin(phi);
+ workEventJet.append(Particle(21, 99, 0, 0, 0, 0, 0, 0,
+ px, py, pz, e, 0., 0, 9.));
+
+ // Extra check on reconstructed y/phi values. If many warnings
+ // of this type, GHOSTENERGY may be set too low.
+#ifdef MLMCHECK
+ int lastIdx = workEventJet.size() - 1;
+ if (abs(y - workEventJet[lastIdx].y()) > ZEROTHRESHOLD ||
+ abs(phi - workEventJet[lastIdx].phi()) > ZEROTHRESHOLD)
+ infoPtr->errorMsg("Warning in MLMhooks:jetAlgorithmInput: "
+ "ghost particle y/phi mismatch");
+#endif
+
+ } // for (i)
+ } // if (jetMatch == 2)
+}
+
+//--------------------------------------------------------------------------
+
+// Step (2b): run jet algorithm and provide common output
+
+void MLMhooks::runJetAlgorithm() {
+
+ // Run the jet clustering algorithm
+ if (jetAlgorithm == 1)
+ cellJet->analyze(workEventJet, eTjetMin, coneRadius, eTseed);
+ else
+ slowJet->analyze(workEventJet);
+
+ // Extract four-momenta of jets with |eta| < etaJetMax and
+ // put into jetMomenta. Note that this is done backwards as
+ // jets are removed with SlowJet.
+ jetMomenta.clear();
+ int iJet = (jetAlgorithm == 1) ? cellJet->size() - 1:
+ slowJet->sizeJet() - 1;
+ for (int i = iJet; i > -1; i--) {
+ Vec4 jetMom = (jetAlgorithm == 1) ? cellJet->pMassive(i) :
+ slowJet->p(i);
+ double eta = Vec4eta(jetMom);
+
+ if (abs(eta) > etaJetMax) {
+ if (jetAlgorithm == 2) slowJet->removeJet(i);
+ continue;
+ }
+ jetMomenta.push_back(jetMom);
+ }
+
+ // Reverse jetMomenta to restore eT/pT ordering
+ reverse(jetMomenta.begin(), jetMomenta.end());
+}
+
+//--------------------------------------------------------------------------
+
+// Step (2c): veto decision (returning true vetoes the event)
+
+bool MLMhooks::matchPartonsToJets(int iType) {
+
+ // Use two different routines for light/heavy jets as
+ // different veto conditions and for clarity
+ if (iType == 0) return (matchPartonsToJetsLight() > 0);
+ else return (matchPartonsToJetsHeavy() > 0);
+}
+
+//--------------------------------------------------------------------------
+
+// Step(2c): light jets
+// Return codes are given indicating the reason for a veto.
+// Although not currently used, they are a useful debugging tool:
+// 0 = no veto
+// 1 = veto as number of jets less than number of partons
+// 2 = veto as exclusive mode and number of jets greater than
+// number of partons
+// 3 = veto as inclusive mode and there would be an extra jet
+// that is harder than any matched soft jet
+// 4 = veto as there is a parton which does not match a jet
+
+int MLMhooks::matchPartonsToJetsLight() {
+
+ // Always veto if number of jets is less than original number of jets
+ if (jetMomenta.size() < typeIdx[0].size()) return 1;
+ // Veto if in exclusive mode and number of jets bigger than original
+ if (exclusive && jetMomenta.size() > typeIdx[0].size()) return 2;
+
+ // Sort partons by eT/pT
+ sortTypeIdx(typeIdx[0]);
+
+ // Number of hard partons
+ int nParton = typeIdx[0].size();
+
+ // Keep track of which jets have been assigned a hard parton
+ vector < bool > jetAssigned;
+ jetAssigned.assign(jetMomenta.size(), false);
+
+ // Jet matching procedure: (1) deltaR between partons and jets
+ if (jetMatch == 1) {
+
+ // Loop over light hard partons and get 4-momentum
+ for (int i = 0; i < nParton; i++) {
+ Vec4 p1 = eventProcess[typeIdx[0][i]].p();
+
+ // Track which jet has the minimal dR measure with this parton
+ int jMin = -1;
+ double dRmin = 0.;
+
+ // Loop over all jets (skipping those already assigned).
+ for (int j = 0; j < int(jetMomenta.size()); j++) {
+ if (jetAssigned[j]) continue;
+
+ // DeltaR between parton/jet and store if minimum
+ double dR = (jetAlgorithm == 1) ?
+ deltaReta(p1, jetMomenta[j]) : deltaRy(p1, jetMomenta[j]);
+ if (jMin < 0 || dR < dRmin) {
+ dRmin = dR;
+ jMin = j;
+ }
+ } // for (j)
+
+ // Check for jet-parton match
+ if (jMin >= 0 && dRmin < coneRadius * coneMatchLight) {
+
+ // If the matched jet is not one of the nParton hardest jets,
+ // the extra left over jet would be harder than some of the
+ // matched jets. This is disallowed, so veto.
+ if (jMin >= nParton) return 3;
+
+ // Mark jet as assigned.
+ jetAssigned[jMin] = true;
+
+ // If no match, then event will be vetoed in all cases
+ } else return 4;
+
+ } // for (i)
+
+ // Jet matching procedure: (2) ghost particles in SlowJet
+ } else {
+
+ // Loop over added 'ghost' particles and find if assigned to a jet
+ for (int i = workEventJet.size() - nParton;
+ i < workEventJet.size(); i++) {
+ int jMin = slowJet->jetAssignment(i);
+
+ // Veto if:
+ // 1) not one of nParton hardest jets
+ // 2) not assigned to a jet
+ // 3) jet has already been assigned
+ if (jMin >= nParton) return 3;
+ if (jMin < 0 || jetAssigned[jMin]) return 4;
+
+ // Mark jet as assigned
+ jetAssigned[jMin] = true;
+
+ } // for (i)
+ } // if (jetMatch)
+
+ // Minimal eT/pT (CellJet/SlowJet) of matched light jets. Needed
+ // later for heavy jet vetos in inclusive mode.
+ if (nParton > 0)
+ eTpTlightMin = (jetAlgorithm == 1) ? jetMomenta[nParton - 1].eT()
+ : jetMomenta[nParton - 1].pT();
+ else
+ eTpTlightMin = -1.;
+
+ // No veto
+ return 0;
+}
+
+//--------------------------------------------------------------------------
+
+// Step(2c): heavy jets
+// Return codes are given indicating the reason for a veto.
+// Although not currently used, they are a useful debugging tool:
+// 0 = no veto as there are no extra jets present
+// 1 = veto as in exclusive mode and extra jets present
+// 2 = veto as in inclusive mode and extra jets were harder
+// than any matched light jet
+
+int MLMhooks::matchPartonsToJetsHeavy() {
+
+ // If there are no extra jets, then accept
+ if (jetMomenta.empty()) return 0;
+
+ // Number of hard partons
+ int nParton = typeIdx[1].size();
+
+ // Remove jets that are close to heavy quarks
+ set < int > removeJets;
+
+ // Jet matching procedure: (1) deltaR between partons and jets
+ if (jetMatch == 1) {
+
+ // Loop over heavy hard partons and get 4-momentum
+ for (int i = 0; i < nParton; i++) {
+ Vec4 p1 = eventProcess[typeIdx[1][i]].p();
+
+ // Loop over all jets, find dR and mark for removal if match
+ for (int j = 0; j < int(jetMomenta.size()); j++) {
+ double dR = (jetAlgorithm == 1) ?
+ deltaReta(p1, jetMomenta[j]) : deltaRy(p1, jetMomenta[j]);
+ if (dR < coneRadiusHeavy * coneMatchHeavy)
+ removeJets.insert(j);
+
+ } // for (j)
+ } // for (i)
+
+ // Jet matching procedure: (2) ghost particles in SlowJet
+ } else {
+
+ // Loop over added 'ghost' particles and if assigned to a jet
+ // then mark this jet for removal
+ for (int i = workEventJet.size() - nParton;
+ i < workEventJet.size(); i++) {
+ int jMin = slowJet->jetAssignment(i);
+ if (jMin >= 0) removeJets.insert(jMin);
+ }
+
+ }
+
+ // Remove jets (backwards order to not disturb indices)
+ for (set < int >::reverse_iterator it = removeJets.rbegin();
+ it != removeJets.rend(); it++)
+ jetMomenta.erase(jetMomenta.begin() + *it);
+
+ // Handle case if there are still extra jets
+ if (!jetMomenta.empty()) {
+
+ // Exclusive mode, so immediate veto
+ if (exclusive) return 1;
+
+ // Inclusive mode; extra jets must be softer than any matched light jet
+ else if (eTpTlightMin >= 0.)
+ for (size_t j = 0; j < jetMomenta.size(); j++) {
+ // CellJet uses eT, SlowJet uses pT
+ if ( (jetAlgorithm == 1 && jetMomenta[j].eT() > eTpTlightMin) ||
+ (jetAlgorithm == 2 && jetMomenta[j].pT() > eTpTlightMin) )
+ return 2;
+ }
+
+ } // if (!jetMomenta.empty())
+
+ // No extra jets were present so no veto
+ return 0;
+}
+
+//==========================================================================
+
+#endif // MLMHOOKS_H
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff