#include <JetMatchingPy8Internal.h>

Inheritance diagram for JetMatchingAlpgen:

Public Member Functions
bool	initAfterBeams ()

	JetMatchingAlpgen ()

	~JetMatchingAlpgen ()

Public Member Functions inherited from JetMatching
bool	canVetoPartonLevelEarly ()

bool	canVetoProcessLevel ()

bool	canVetoStep ()

bool	doVetoPartonLevelEarly (const Pythia8::Event &event)

bool	doVetoProcessLevel (Pythia8::Event &process)

bool	doVetoStep (int, int, int, const Pythia8::Event &)

	JetMatching ()

int	numberVetoStep ()

	~JetMatching ()

Private Member Functions
void	jetAlgorithmInput (const Pythia8::Event &, int)

bool	matchPartonsToJets (int)

int	matchPartonsToJetsHeavy ()

int	matchPartonsToJetsLight ()

void	runJetAlgorithm ()

void	sortIncomingProcess (const Pythia8::Event &)

void	sortTypeIdx (vector< int > &vecIn)

Static Private Attributes
static const double	GHOSTENERGY = 1e-15

static const double	ZEROTHRESHOLD = 1e-10

Additional Inherited Members
Protected Types inherited from JetMatching
enum	partonTypes { ID_CHARM =4, ID_BOT =5, ID_TOP =6, ID_LEPMIN =11, ID_LEPMAX =16, ID_GLUON =21, ID_PHOTON =22 }

enum	vetoStatus { NONE, LESS_JETS, MORE_JETS, HARD_JET, UNMATCHED_PARTON }

Protected Attributes inherited from JetMatching
Pythia8::CellJet *	cellJet

double	coneMatchHeavy

double	coneMatchLight

double	coneRadius

double	coneRadiusHeavy

bool	doMerge

bool	doShowerKt

double	etaJetMax

double	etaJetMaxAlgo

double	eTjetMin

double	eTpTlightMin

double	eTseed

double	eTthreshold

Pythia8::Event	eventProcess

Pythia8::Event	eventProcessOrig

bool	exclusive

int	exclusiveMode

int	jetAlgorithm

int	jetAllow

int	jetMatch

vector< Pythia8::Vec4 >	jetMomenta

int	nEta

int	nJet

int	nJetMax

int	nPhi

Pythia8::SlowJet *	slowJet

Pythia8::SlowJet *	slowJetHard

int	slowJetPower

vector< int >	typeIdx [3]

set< int >	typeSet [3]

Pythia8::Event	workEventJet

Static Protected Attributes inherited from JetMatching
static const bool	MATCHINGCHECK = false

static const bool	MATCHINGDEBUG = false

Detailed Description

Definition at line 132 of file JetMatchingPy8Internal.h.

Constructor & Destructor Documentation

JetMatchingAlpgen::JetMatchingAlpgen ( )

inline

Definition at line 137 of file JetMatchingPy8Internal.h.

137 { }

JetMatchingAlpgen::~JetMatchingAlpgen ( )

inline

Definition at line 138 of file JetMatchingPy8Internal.h.

138 { }

Member Function Documentation

bool JetMatchingAlpgen::initAfterBeams ( )

virtual

Implements JetMatching.

Definition at line 155 of file JetMatchingPy8Internal.cc.

References gather_cfg::cout, and dtDQMClient_cfg::resolution.

Referenced by JetMatchingAlpgenInputAlpgen::initAfterBeams().

                                        {
 
   // Read in parameters
   doMerge         = settingsPtr->flag("JetMatching:merge");
   jetAlgorithm    = settingsPtr->mode("JetMatching:jetAlgorithm");
   nJet            = settingsPtr->mode("JetMatching:nJet");
   nJetMax         = settingsPtr->mode("JetMatching:nJetMax");
   eTjetMin        = settingsPtr->parm("JetMatching:eTjetMin");
   coneRadius      = settingsPtr->parm("JetMatching:coneRadius");
   etaJetMax       = settingsPtr->parm("JetMatching:etaJetMax");
 
   // Use etaJetMax + coneRadius in input to jet algorithms
   etaJetMaxAlgo   = etaJetMax + coneRadius;
 
   // CellJet specific
   nEta            = settingsPtr->mode("JetMatching:nEta");
   nPhi            = settingsPtr->mode("JetMatching:nPhi");
   eTseed          = settingsPtr->parm("JetMatching:eTseed");
   eTthreshold     = settingsPtr->parm("JetMatching:eTthreshold");
 
   // SlowJet specific
   slowJetPower    = settingsPtr->mode("JetMatching:slowJetPower");
   coneMatchLight  = settingsPtr->parm("JetMatching:coneMatchLight");
   coneRadiusHeavy = settingsPtr->parm("JetMatching:coneRadiusHeavy");
   if (coneRadiusHeavy < 0.) coneRadiusHeavy = coneRadius;
   coneMatchHeavy  = settingsPtr->parm("JetMatching:coneMatchHeavy");
 
   // Matching procedure
   jetAllow        = settingsPtr->mode("JetMatching:jetAllow");
   jetMatch        = settingsPtr->mode("JetMatching:jetMatch");
   exclusiveMode   = settingsPtr->mode("JetMatching: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 JetMatchingAlpgen: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 JetMatchingAlpgen: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;
 }

void JetMatchingAlpgen::jetAlgorithmInput	(	const Pythia8::Event &	,
		int
	)

privatevirtual

Implements JetMatching.

Definition at line 330 of file JetMatchingPy8Internal.cc.

References funct::abs(), funct::cos(), alignCSCRings::e, end, create_public_lumi_plots::exp, spr::find(), i, customizeTrackingMonitorSeedNumber::idx, phi, RecoTauCleanerPlugins::pt, funct::sin(), findQualityFiles::size, mathSSE::sqrt(), and detailsBasic3DVector::y.

                                                                        {
 
   // 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 >= ID_LEPMIN && id <= ID_LEPMAX) || id == ID_TOP
         || id == ID_PHOTON) {
         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   = pIn.rap();
       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( ID_GLUON, 99, 0, 0, 0, 0, 0, 0, px, py, pz, e);
 
       // Extra check on reconstructed y/phi values. If many warnings
       // of this type, GHOSTENERGY may be set too low.
       if (MATCHINGCHECK) {
       int lastIdx = workEventJet.size() - 1;
       if (abs(y   - workEventJet[lastIdx].y())   > ZEROTHRESHOLD ||
           abs(phi - workEventJet[lastIdx].phi()) > ZEROTHRESHOLD)
         infoPtr->errorMsg("Warning in JetMatchingAlpgen:jetAlgorithmInput: "
             "ghost particle y/phi mismatch");
       }
 
     } // for (i)
   } // if (jetMatch == 2)
 }

bool JetMatchingAlpgen::matchPartonsToJets ( int iType )

privatevirtual

Implements JetMatching.

Definition at line 463 of file JetMatchingPy8Internal.cc.

                                                     {
 
   // 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);
 }

int JetMatchingAlpgen::matchPartonsToJetsHeavy ( )

privatevirtual

Implements JetMatching.

Definition at line 584 of file JetMatchingPy8Internal.cc.

References PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi::dR, i, j, NONE, and p1.

                                                {
 
   // If there are no extra jets, then accept
   if (jetMomenta.empty()) return NONE;
 
   // 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) ?
             REtaPhi(p1, jetMomenta[j]) : RRapPhi(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 MORE_JETS;
 
     // 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 HARD_JET;
       }
 
   } // if (!jetMomenta.empty())
 
   // No extra jets were present so no veto
   return NONE;
 }

int JetMatchingAlpgen::matchPartonsToJetsLight ( )

privatevirtual

Implements JetMatching.

Definition at line 484 of file JetMatchingPy8Internal.cc.

References PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi::dR, i, j, NONE, and p1.

                                                {
 
   // Always veto if number of jets is less than original number of jets
   if (jetMomenta.size() < typeIdx[0].size()) return LESS_JETS;
   // Veto if in exclusive mode and number of jets bigger than original
   if (exclusive && jetMomenta.size() > typeIdx[0].size()) return MORE_JETS;
 
   // 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) 
           ? REtaPhi(p1, jetMomenta[j]) : RRapPhi(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 HARD_JET;
 
         // Mark jet as assigned.
         jetAssigned[jMin] = true;
 
       // If no match, then event will be vetoed in all cases
       } else return UNMATCHED_PARTON;
 
     } // 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 HARD_JET;
       if (jMin < 0 || jetAssigned[jMin]) return UNMATCHED_PARTON;
 
       // 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 NONE;
 }

void JetMatchingAlpgen::runJetAlgorithm ( )

privatevirtual

Implements JetMatching.

Definition at line 429 of file JetMatchingPy8Internal.cc.

References funct::abs(), eta(), and i.

                                         {
 
   // 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 = jetMom.eta();
 
     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());
 }

void JetMatchingAlpgen::sortIncomingProcess ( const Pythia8::Event & )

privatevirtual

Implements JetMatching.

Definition at line 282 of file JetMatchingPy8Internal.cc.

References funct::abs(), i, customizeTrackingMonitorSeedNumber::idx, and m.

                                                               {
 
   // 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() == ID_GLUON
       || (eventProcess[i].idAbs() <= ID_BOT
       && abs(eventProcess[i].m()) < ZEROTHRESHOLD)) idx = 0;
 
     // Heavy jets
     else if (eventProcess[i].idAbs() >= ID_CHARM
       && eventProcess[i].idAbs() <= ID_TOP) 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);
 }

void JetMatchingAlpgen::sortTypeIdx ( vector< int > & vecIn )

private

Definition at line 132 of file JetMatchingPy8Internal.cc.

References i, j, and swap().

                                                          {
   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]);
   }
 }

Member Data Documentation

const double JetMatchingAlpgen::GHOSTENERGY = 1e-15

staticprivate

Definition at line 157 of file JetMatchingPy8Internal.h.

const double JetMatchingAlpgen::ZEROTHRESHOLD = 1e-10

staticprivate

Definition at line 157 of file JetMatchingPy8Internal.h.

Public Member Functions

Private Member Functions

Static Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation