JetMatchingMadgraph Class Reference

#include <JetMatchingPy8Internal.h>

Inheritance diagram for JetMatchingMadgraph:

Public Member Functions
bool	canVetoStep ()
bool	doVetoStep (int, int, int, const Pythia8::Event &)
vector< double >	GetDJR ()
bool	initAfterBeams ()
	JetMatchingMadgraph ()
vector< int >	nMEPartons ()
int	numberVetoStep ()
	~JetMatchingMadgraph ()
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 ()

Public Attributes
Pythia8::SlowJet *	slowJetDJR

Protected Member Functions
void	ClearDJR ()
void	ClearnME ()
bool	doShowerKtVeto (double pTfirst)
void	jetAlgorithmInput (const Pythia8::Event &, int)
bool	matchPartonsToJets (int)
int	matchPartonsToJetsHeavy ()
int	matchPartonsToJetsLight ()
void	runJetAlgorithm ()
void	SetDJR (const Pythia8::Event &event)
void	SetnME ()
void	sortIncomingProcess (const Pythia8::Event &)

Protected Attributes
double	clFact
vector< double >	DJR
bool	doFxFx
vector< int >	nME
int	nPartonsNow
int	nQmatch
vector< int >	origTypeIdx [3]
double	qCut
double	qCutME
double	qCutMESq
double	qCutSq
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

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 }
Static Protected Attributes inherited from JetMatching
static const bool	MATCHINGCHECK = false
static const bool	MATCHINGDEBUG = false

Detailed Description

Definition at line 165 of file JetMatchingPy8Internal.h.

Constructor & Destructor Documentation

JetMatchingMadgraph::JetMatchingMadgraph ( )

inline

Definition at line 170 of file JetMatchingPy8Internal.h.

{ }

JetMatchingMadgraph::~JetMatchingMadgraph ( )

inline

Definition at line 171 of file JetMatchingPy8Internal.h.

{ }

Member Function Documentation

bool JetMatchingMadgraph::canVetoStep ( )

inline

Definition at line 178 of file JetMatchingPy8Internal.h.

References JetMatching::doShowerKt.

{ return doShowerKt; }

void JetMatchingMadgraph::ClearDJR ( )

inlineprotected

Definition at line 199 of file JetMatchingPy8Internal.h.

References DJR.

Referenced by SetDJR().

{ DJR.resize(0);}

void JetMatchingMadgraph::ClearnME ( )

inlineprotected

Definition at line 200 of file JetMatchingPy8Internal.h.

References nME.

Referenced by SetnME().

{ nME.resize(0);}

bool JetMatchingMadgraph::doShowerKtVeto ( double  pTfirst )

protected

Definition at line 855 of file JetMatchingPy8Internal.cc.

References i, bookConverter::min, and funct::pow().

                                                       {

  // Only check veto in the shower-kT scheme.
  if ( !doShowerKt ) return false;

  // Reset veto code
  bool doVeto = false;

  // Find the (kinematical) pT of the softest (light) parton in the hard
  // process.   
  int nParton = typeIdx[0].size();
  double pTminME=1e10;  
  for ( int i = 0; i < nParton; ++i)
    pTminME = min(pTminME,eventProcess[typeIdx[0][i]].pT());

  // Veto if the softest hard process parton is below Qcut.
  if ( nParton > 0 && pow(pTminME,2) < qCutSq ) doVeto = true;

  // For non-highest multiplicity, veto if the hardest emission is harder
  // than Qcut.
  if ( exclusive && pow(pTfirst,2) > qCutSq ) {
    doVeto = true;
  // For highest multiplicity sample, veto if the hardest emission is harder
  // than the hard process parton.
  } else if ( !exclusive && nParton > 0 && pTfirst > pTminME ) {
    doVeto = true;      
  }

  // Return veto
  return doVeto;

}

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

Definition at line 799 of file JetMatchingPy8Internal.cc.

References i, and findQualityFiles::size.

                       {

  // Do not perform any veto if not in the Shower-kT scheme.
  if ( !doShowerKt ) return false;

  // Default to no veto in case the hard input matrix element already has too
  // many partons (same as in Pythia6).
  if ( int(typeIdx[0].size()) > nJetMax ) return false;

  // Do nothing for emissions after the first one.
  if ( nISR + nFSR > 1 ) return false;

  // Do nothing in resonance decay showers.
  if (iPos == 5) return false;

  // Clear the event of MPI systems and resonace decay products. Store trimmed
  // event in workEvent.
  sortIncomingProcess(event);

  // Get (kinematical) pT of first emission
  double pTfirst = 0.;
  for (int i = 0; i < workEvent.size(); i++){
    // Since this event only contains one emission, this test is enough
    // to isolate this emission.
    if ( workEvent[i].isFinal()
      && (workEvent[i].statusAbs()==43 || workEvent[i].statusAbs()==51)) {
      // Only check partons originating from QCD splittings.
      int iPos = 1;
      bool QCDemission = true;
      while ( workEvent[iPos].statusAbs() > 23 ) {
        if ( workEvent[iPos].id() == 22 || workEvent[iPos].id() == 23
          || workEvent[iPos].idAbs() == 24){
          QCDemission = false;
          break;
        }
        iPos = workEvent[iPos].mother1();
      }
      // Get kinematical pT.
      if (QCDemission) {
        pTfirst = workEvent[i].pT();
        break;
      }
    }
  }

  // Check veto.
  if ( doShowerKtVeto(pTfirst) ) return true;

  // No veto if come this far.
  return false;

}

vector<double> JetMatchingMadgraph::GetDJR ( )

inline

Definition at line 185 of file JetMatchingPy8Internal.h.

References DJR.

{ return DJR;}

bool JetMatchingMadgraph::initAfterBeams ( )

virtual

Implements JetMatching.

Definition at line 662 of file JetMatchingPy8Internal.cc.

References gather_cfg::cout, MadgraphPar::getParam(), MadgraphPar::getParamAsInt(), MadgraphPar::haveParam(), NULL, MadgraphPar::parse(), funct::pow(), and MadgraphPar::printParams().

Referenced by JetMatchingMadgraphInputAlpgen::initAfterBeams().

                                         {

  // Read in Madgraph specific configuration variables
  bool setMad    = settingsPtr->flag("JetMatching:setMad");

  // If Madgraph parameters are present, then parse in MadgraphPar object
  MadgraphPar par(infoPtr);
  string parStr = infoPtr->header("MGRunCard");
  if (!parStr.empty()) {
    par.parse(parStr);
    par.printParams();
  }
 
  // Set Madgraph merging parameters from the file if requested
  if (setMad) {
    if ( par.haveParam("xqcut")    && par.haveParam("maxjetflavor")
      && par.haveParam("alpsfact") && par.haveParam("ickkw") ) {
      settingsPtr->flag("JetMatching:merge", par.getParam("ickkw"));
      settingsPtr->parm("JetMatching:qCut", par.getParam("xqcut"));
      settingsPtr->mode("JetMatching:nQmatch",
        par.getParamAsInt("maxjetflavor"));
      settingsPtr->parm("JetMatching:clFact",
        clFact = par.getParam("alpsfact"));
      if (par.getParamAsInt("ickkw") == 0)
        infoPtr->errorMsg("Error in JetMatchingMadgraph:init: "
          "Madgraph file parameters are not set for merging");

    // Warn if setMad requested, but one or more parameters not present
    } else {
       infoPtr->errorMsg("Warning in JetMatchingMadgraph:init: "
          "Madgraph merging parameters not found");
       if (!par.haveParam("xqcut")) infoPtr->errorMsg("Warning in "
          "JetMatchingMadgraph:init: No xqcut");
       if (!par.haveParam("ickkw")) infoPtr->errorMsg("Warning in "
          "JetMatchingMadgraph:init: No ickkw");
       if (!par.haveParam("maxjetflavor")) infoPtr->errorMsg("Warning in "
          "JetMatchingMadgraph:init: No maxjetflavor");
       if (!par.haveParam("alpsfact")) infoPtr->errorMsg("Warning in "
          "JetMatchingMadgraph:init: No alpsfact");
    }
  }

  // Read in FxFx matching parameters
  doFxFx       = settingsPtr->flag("JetMatching:doFxFx");
  nPartonsNow  = settingsPtr->mode("JetMatching:nPartonsNow");
  qCutME       = settingsPtr->parm("JetMatching:qCutME");
  qCutMESq     = pow(qCutME,2);

  // Read in Madgraph merging parameters
  doMerge      = settingsPtr->flag("JetMatching:merge");
  doShowerKt   = settingsPtr->flag("JetMatching:doShowerKt");
  qCut         = settingsPtr->parm("JetMatching:qCut");
  nQmatch      = settingsPtr->mode("JetMatching:nQmatch");
  clFact       = settingsPtr->parm("JetMatching:clFact");

  // Read in jet algorithm parameters
  jetAlgorithm   = settingsPtr->mode("JetMatching:jetAlgorithm");
  nJetMax        = settingsPtr->mode("JetMatching:nJetMax");
  eTjetMin       = settingsPtr->parm("JetMatching:eTjetMin");
  coneRadius     = settingsPtr->parm("JetMatching:coneRadius");
  etaJetMax      = settingsPtr->parm("JetMatching:etaJetMax");
  slowJetPower   = settingsPtr->mode("JetMatching:slowJetPower");

  // Matching procedure
  jetAllow       = settingsPtr->mode("JetMatching:jetAllow");
  exclusiveMode  = settingsPtr->mode("JetMatching:exclusive");
  qCutSq         = pow(qCut,2);
  etaJetMaxAlgo  = etaJetMax;

  // 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 (nJetMax < 0) {
      infoPtr->errorMsg("Warning in JetMatchingMadgraph:init: "
        "missing jet multiplicity information; running in exclusive mode");
      exclusiveMode = 1;
    }
  }

  // Initialise chosen jet algorithm.
  // Currently, this only supports the kT-algorithm in SlowJet.
  // Use the QCD distance measure by default.
  jetAlgorithm = 2;
  slowJetPower = 1;
  slowJet = new SlowJet(slowJetPower, coneRadius, eTjetMin,
    etaJetMaxAlgo, 2, 2, NULL, false);

  // For FxFx, also initialise jet algorithm to define matrix element jets.
  // Currently, this only supports the kT-algorithm in SlowJet.
  // Use the QCD distance measure by default.
  slowJetHard = new SlowJet(slowJetPower, coneRadius, qCutME,
    etaJetMaxAlgo, 2, 2, NULL, false);

  // To access the DJR's
  slowJetDJR = new SlowJet(slowJetPower, coneRadius, 0.0/*qCutME*/,
    etaJetMaxAlgo, 2, 2, NULL, false);

  // 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 = (exclusiveMode)         ? "exclusive" : "inclusive";
  cout << endl
       << " *-----  Madgraph matching parameters  -----*" << endl
       << " |  qCut                |  " << setw(14)
       << qCut << "  |" << endl
       << " |  nQmatch             |  " << setw(14)
       << nQmatch << "  |" << endl
       << " |  clFact              |  " << setw(14)
       << clFact << "  |" << endl
       << " |  Jet algorithm       |  " << setw(14)
       << jetStr << "  |" << endl
       << " |  eTjetMin            |  " << setw(14)
       << eTjetMin << "  |" << endl
       << " |  etaJetMax           |  " << setw(14)
       << etaJetMax << "  |" << endl
       << " |  jetAllow            |  " << setw(14)
       << jetAllow << "  |" << endl
       << " |  Mode                |  " << setw(14)
       << modeStr << "  |" << endl
       << " *-----------------------------------------*" << endl;

  return true;
}

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

protectedvirtual

Implements JetMatching.

Definition at line 1038 of file JetMatchingPy8Internal.cc.

References end, spr::find(), i, and customizeTrackingMonitorSeedNumber::idx.

                                                                         {

  // 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 || (id > nQmatch && id!=21)) {
        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)
}

bool JetMatchingMadgraph::matchPartonsToJets ( int  iType )

protectedvirtual

Implements JetMatching.

Definition at line 1114 of file JetMatchingPy8Internal.cc.

                                                      {

  // Use two different routines for light/heavy jets as
  // different veto conditions and for clarity
  if (iType == 0) {
        SetDJR(workEventJet);
        SetnME();
        return (matchPartonsToJetsLight() > 0);
  }     
  else            return (matchPartonsToJetsHeavy() > 0);
}

int JetMatchingMadgraph::matchPartonsToJetsHeavy ( )

protectedvirtual

Implements JetMatching.

Definition at line 1376 of file JetMatchingPy8Internal.cc.

References NONE.

                                                 {

  // Currently, heavy jets are unmatched
  // If there are no extra jets, then accept
  if (jetMomenta.empty()) return NONE;

  // No extra jets were present so no veto
  return NONE;
}

int JetMatchingMadgraph::matchPartonsToJetsLight ( )

protectedvirtual

Implements JetMatching.

Definition at line 1139 of file JetMatchingPy8Internal.cc.

References alignCSCRings::e, i, max(), NONE, findQualityFiles::size, and python.multivaluedict::sort().

                                                 {

  // Count the number of hard partons
  int nParton = typeIdx[0].size();
  if((int)origTypeIdx[0].size()>nJetMax)return 1;

  // Initialize SlowJet with current working event
  if (!slowJet->setup(workEventJet) ) {
    infoPtr->errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
      "Light: the SlowJet algorithm failed on setup");
    return NONE;
  }
  double localQcutSq = qCutSq;
  double dOld = 0.0;
  // Cluster in steps to find all hadronic jets at the scale qCut
  while ( slowJet->sizeAll() - slowJet->sizeJet() > 0 ) {
    if( slowJet->dNext() > localQcutSq ) break;
    dOld = slowJet->dNext();
    slowJet->doStep();
  }
  int nJets = slowJet->sizeJet();
  int nClus = slowJet->sizeAll();

  // Debug printout.
  if (MATCHINGDEBUG) slowJet->list(true);

  // Count of the number of hadronic jets in SlowJet accounting
  int nCLjets = nClus - nJets;
  // Get number of partons. Different for MLM and FxFx schemes.
  int nRequested = (doFxFx) ? nPartonsNow : nParton;

  // Veto event if too few hadronic jets
  if ( nCLjets < nRequested ) return LESS_JETS;

  // In exclusive mode, do not allow more hadronic jets than partons
  if ( exclusive && !doFxFx ) {
    if ( nCLjets > nRequested ) return MORE_JETS;
  } else {

    // For FxFx, in the non-highest multipicity, all jets need to matched to
    // partons. For nCLjets > nRequested, this is not possible. Hence, we can
    // veto here already.
    if ( doFxFx && nRequested < nJetMax && nCLjets > nRequested )
      return MORE_JETS;

    // Now continue in inclusive mode.
    // In inclusive mode, there can be more hadronic jets than partons,
    // provided that all partons are properly matched to hadronic jets.
    // Start by setting up the jet algorithm.
    if (!slowJet->setup(workEventJet) ) {
      infoPtr->errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
        "Light: the SlowJet algorithm failed on setup");
      return NONE;
    }

    // For FxFx, continue clustering as long as the jet separation is above
    // qCut.
    if (doFxFx) {
      while ( slowJet->sizeAll() - slowJet->sizeJet() > 0 ) {
        if( slowJet->dNext() > localQcutSq ) break;
        slowJet->doStep();
      }
    // For MLM, cluster into hadronic jets until there are the same number as
    // partons.
    } else {
      while ( slowJet->sizeAll() - slowJet->sizeJet() > nParton )
        slowJet->doStep();
    }

    // Sort partons in pT.  Update local qCut value.
    //  Hadronic jets are already sorted in pT.
    localQcutSq = dOld;
    if ( clFact >= 0. && nParton > 0 ) {
       vector<double> partonPt;
       for (int i = 0; i < nParton; ++i)
         partonPt.push_back( eventProcess[typeIdx[0][i]].pT2() );
       sort( partonPt.begin(), partonPt.end());
       localQcutSq = max( qCutSq, partonPt[0]);
    }
    nJets = slowJet->sizeJet();
    nClus = slowJet->sizeAll();
  }
  // Update scale if clustering factor is non-zero
  if ( clFact != 0. ) localQcutSq *= pow2(clFact);

  Event tempEvent;
  tempEvent.init( "(tempEvent)", particleDataPtr);
  int nPass = 0;
  double pTminEstimate = -1.;
  // Construct a master copy of the event containing only the
  // hardest nParton hadronic clusters. While constructing the event,
  // the parton type (ID_GLUON) and status (98,99) are arbitrary.
  for (int i = nJets; i < nClus; ++i) {
    tempEvent.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0, slowJet->p(i).px(),
      slowJet->p(i).py(), slowJet->p(i).pz(), slowJet->p(i).e() );
    ++nPass;
    pTminEstimate = max( pTminEstimate, slowJet->pT(i));
    if(nPass == nRequested) break;
  }

  int tempSize = tempEvent.size();
  // This keeps track of which hadronic jets are matched to parton
  vector<bool> jetAssigned;
  jetAssigned.assign( tempSize, false);

  // This keeps track of which partons are matched to which hadronic
  // jets.
  vector< vector<bool> > partonMatchesJet;
  for (int i=0; i < nParton; ++i )
    partonMatchesJet.push_back( vector<bool>(tempEvent.size(),false) );

  // Begin matching.
  // Do jet matching for FxFx.
  // Make sure that the nPartonsNow hardest hadronic jets are matched to any
  // of the nPartonsNow (+1) partons. This matching is done by attaching a jet
  // from the list of unmatched hadronic jets, and appending a jet from the
  // list of partonic jets, one at a time. The partonic jet will be clustered
  // with the hadronic jet or the beam if the distance measure is below the
  // cut. The hadronic jet is matched once this happens. Otherwise, another
  // partonic jet is tried. When a hadronic jet is matched to a partonic jet,
  // it is removed from the list of unmatched hadronic jets. This process
  // continues until the nPartonsNow hardest hadronic jets are matched to
  // partonic jets, or it is not possible to make a match for a hadronic jet.
  int iNow = 0;
  while ( doFxFx && iNow < tempSize ) {
 
    // Check if this shower jet matches any partonic jet.
    Event tempEventJet;
    tempEventJet.init("(tempEventJet)", particleDataPtr);
    for (int i=0; i < nParton; ++i ) {

      //for (int j=0; j < tempSize; ++j )
      //  if ( partonMatchesJet[i][j]) continue;

      // Attach a single hadronic jet.
      tempEventJet.clear();
      tempEventJet.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0,
        tempEvent[iNow].px(), tempEvent[iNow].py(),
        tempEvent[iNow].pz(), tempEvent[iNow].e() );
      // Attach the current parton.
      Vec4 pIn = eventProcess[typeIdx[0][i]].p();
      tempEventJet.append( ID_GLUON, 99, 0, 0, 0, 0, 0, 0,
        pIn.px(), pIn.py(), pIn.pz(), pIn.e() );

      // Setup jet algorithm.
      if ( !slowJet->setup(tempEventJet) ) {
        infoPtr->errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
          "Light: the SlowJet algorithm failed on setup");
        return NONE;
      }

      // These are the conditions for the hadronic jet to match the parton
      //  at the local qCut scale
      if ( slowJet->iNext() == tempEventJet.size() - 1
        && slowJet->jNext() > -1 && slowJet->dNext() < localQcutSq ) {
        jetAssigned[iNow] = true;
        partonMatchesJet[i][iNow] = true;
      }

    } // End loop over hard partons.

    // Veto if the jet could not be assigned to any parton.
    if ( !jetAssigned[iNow] ) return UNMATCHED_PARTON;

    // Continue;
    ++iNow;
  }

  // Do jet matching for MLM.
  // Take the list of unmatched hadronic jets and append a parton, one at
  // a time. The parton will be clustered with the "closest" hadronic jet
  // or the beam if the distance measure is below the cut. When a hadronic
  // jet is matched to a parton, it is removed from the list of unmatched
  // hadronic jets. This process continues until all hadronic jets are
  // matched to partons or it is not possible to make a match.
  iNow = 0;
  while (!doFxFx && iNow < nParton ) {
    Event tempEventJet;
    tempEventJet.init("(tempEventJet)", particleDataPtr);
    for (int i = 0; i < tempSize; ++i) {
      if (jetAssigned[i]) continue;
      Vec4 pIn = tempEvent[i].p();
      // Append unmatched hadronic jets
      tempEventJet.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0,
        pIn.px(), pIn.py(), pIn.pz(), pIn.e() );
    }

    Vec4 pIn = eventProcess[typeIdx[0][iNow]].p();
    // Append the current parton
    tempEventJet.append( ID_GLUON, 99, 0, 0, 0, 0, 0, 0,
      pIn.px(), pIn.py(), pIn.pz(), pIn.e() );
    if ( !slowJet->setup(tempEventJet) ) {
      infoPtr->errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
        "Light: the SlowJet algorithm failed on setup");
      return NONE;
    }
    // These are the conditions for the hadronic jet to match the parton
    //  at the local qCut scale
    if ( slowJet->iNext() == tempEventJet.size() - 1
      && slowJet->jNext() > -1 && slowJet->dNext() < localQcutSq ) {
      int iKnt = -1;
      for (int i = 0; i != tempSize; ++i) {
        if (jetAssigned[i]) continue;
        ++iKnt;
        // Identify the hadronic jet that matches the parton
        if (iKnt == slowJet->jNext() ) jetAssigned[i] = true;
      }
    } else {
      return UNMATCHED_PARTON;
    }
    ++iNow;
  }

  // Minimal eT/pT (CellJet/SlowJet) of matched light jets.
  // Needed later for heavy jet vetos in inclusive mode.
  // This information is not used currently.
  if (nParton > 0 && pTminEstimate > 0) eTpTlightMin = pTminEstimate;
  else eTpTlightMin = -1.;

  //SetDJR(workEventJet);
  //SetnME();
  // No veto
  //std::cout<<"No VETO"<<std::endl;
  return NONE;
}

vector<int> JetMatchingMadgraph::nMEPartons ( )

inline

Definition at line 186 of file JetMatchingPy8Internal.h.

References nME.

{return nME;} 

int JetMatchingMadgraph::numberVetoStep ( )

inline

Definition at line 177 of file JetMatchingPy8Internal.h.

{return 1;}

void JetMatchingMadgraph::runJetAlgorithm ( )

protectedvirtual

Implements JetMatching.

Definition at line 1108 of file JetMatchingPy8Internal.cc.

{; }

void JetMatchingMadgraph::SetDJR ( const Pythia8::Event &  event )

inlineprotected

Definition at line 202 of file JetMatchingPy8Internal.h.

References ClearDJR(), DJR, i, query::result, slowJetDJR, and mathSSE::sqrt().

                                          {

   // Clear members.
   ClearDJR();

   vector<double> result;

    // Initialize SlowJetDJR jet algorithm with event
    if (!slowJetDJR->setup(event) ) {
      infoPtr->errorMsg("Warning in JetMatchingMadgraph:iGetDJR"
        ": the SlowJet algorithm failed on setup");
      return;
    }

    // Cluster in steps to find all hadronic jets
    while ( slowJetDJR->sizeAll() - slowJetDJR->sizeJet() > 0 ) {
      // Save the next clustering scale.
      result.push_back(slowJetDJR->dNext());
      // Perform step.
      slowJetDJR->doStep();
    }

    // Save clustering scales in reserve order.
    for (int i=int(result.size())-1; i > 0; --i){
//  std::cout<<"Saving DJR "<<i<<" "<<log10(sqrt(result[i]))<<std::endl;
      DJR.push_back(log10(sqrt(result[i])));
    }
  }

void JetMatchingMadgraph::SetnME ( )

inlineprotected

Definition at line 230 of file JetMatchingPy8Internal.h.

References ClearnME(), nME, origTypeIdx, query::result, findQualityFiles::size, and JetMatching::typeIdx.

                {
    ClearnME();
    vector<int> result;
        nME.push_back(origTypeIdx[0].size());
    nME.push_back(typeIdx[0].size());
    //std::cout<<"Number of partons "<<origTypeIdx[0].size()<<" "<<typeIdx[0].size()<<std::endl;
  } 

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

protectedvirtual

Implements JetMatching.

Definition at line 892 of file JetMatchingPy8Internal.cc.

References reco::flavour(), i, customizeTrackingMonitorSeedNumber::idx, j, contentValuesFiles::parts, pileupReCalc_HLTpaths::scale, and mathSSE::sqrt().

                                                                {

  // Remove resonance decays from original process and keep only final
  // state. Resonances will have positive status code after this step.
  omitResonanceDecays(eventProcessOrig, true);

  ClearDJR();
  ClearnME();

  // For FxFx, pre-cluster partons in the event into jets.
  if (doFxFx) {

    // Get final state partons
    eventProcess.clear();
    workEventJet.clear();
    for( int i=0; i < workEvent.size(); ++i) {
      // Original AG+Py6 algorithm explicitly excludes tops,
      // leptons and photons.
      int id = workEvent[i].idAbs();
      if ((id >= ID_LEPMIN && id <= ID_LEPMAX) || id == ID_TOP
        || id == ID_PHOTON || id == 23 || id == 24 || id == 25) {
        eventProcess.append(workEvent[i]);
      } else {
        workEventJet.append(workEvent[i]);
      }
    }

    // Initialize SlowJetHard jet algorithm with current working event
    if (!slowJetHard->setup(workEventJet) ) {
      infoPtr->errorMsg("Warning in JetMatchingMadgraph:sortIncomingProcess"
        ": the SlowJet algorithm failed on setup");
      return;
    }

    // Get matrix element cut scale.
    double localQcutSq = qCutMESq;
    // Cluster in steps to find all hadronic jets at the scale qCutME
    while ( slowJetHard->sizeAll() - slowJetHard->sizeJet() > 0 ) {
      // Done if next step is above qCut
      if( slowJetHard->dNext() > localQcutSq ) break;
      // Done if we're at or below the number of partons in the Born state.
      if( slowJetHard->sizeAll()-slowJetHard->sizeJet() <= nPartonsNow) break;
      slowJetHard->doStep();
    }

    // Construct a master copy of the event containing only the
    // hardest nPartonsNow hadronic clusters. While constructing the event,
    // the parton type (ID_GLUON) and status (98,99) are arbitrary.
    int nJets = slowJetHard->sizeJet();
    int nClus = slowJetHard->sizeAll();
    int nNow = 0;
    for (int i = nJets; i < nClus; ++i) {
      vector<int> parts;
      if (i < nClus-nJets) parts = slowJetHard->clusConstituents(i);
      else parts = slowJetHard->constituents(nClus-nJets-i);
      int flavour = ID_GLUON;
      for(int j=0; j < int(parts.size()); ++j)
        if (workEventJet[parts[j]].id() == ID_BOT)
          flavour = ID_BOT;
      eventProcess.append( flavour, 98,
        workEventJet[parts.back()].mother1(),
        workEventJet[parts.back()].mother2(),
        workEventJet[parts.back()].daughter1(),
        workEventJet[parts.back()].daughter2(),
        0, 0, slowJetHard->p(i).px(), slowJetHard->p(i).py(),
        slowJetHard->p(i).pz(), slowJetHard->p(i).e() );
      nNow++;
    }

    // Done. Clean-up
    workEventJet.clear();

  // For MLM matching, simply take hard process state from workEvent,
  // without any preclustering.
  } else {
    eventProcess = workEvent;
  }

  // Sort original process final state into light/heavy jets and 'other'.
  // Criteria:
  //   1 <= ID <= nQmatch, or ID == 21         --> light jet (typeIdx[0])
  //   nQMatch < ID                            --> 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();
    origTypeIdx[i].clear();
  }
  for (int i = 0; i < eventProcess.size(); i++) {
    // Ignore non-final state and default to 'other'
    if (!eventProcess[i].isFinal()) continue;
    int idx = 2;
    int orig_idx = 2;

    // Light jets: all gluons and quarks with id less than or equal to nQmatch
    if (eventProcess[i].id() == ID_GLUON
      || (eventProcess[i].idAbs() <= nQmatch) ) {
      orig_idx = 0;
      // Crucial point: MG puts the scale of a non-QCD particle to eCM. For
      // such particles, we should keep the default "2"
      if ( eventProcess[i].scale() < 1.999*sqrt(infoPtr->eA()*infoPtr->eB()) )
        idx = 0;
    }

    // Heavy jets:  all quarks with id greater than nQmatch
    else if (eventProcess[i].idAbs() > nQmatch
      && eventProcess[i].idAbs() <= ID_TOP) {
      idx = 1;
      orig_idx = 1;

    } else {
      idx = 2;
      orig_idx = 2;
    }

    // Store
    typeIdx[idx].push_back(i);
    typeSet[idx].insert(eventProcess[i].daughter1());
    origTypeIdx[orig_idx].push_back(i);
  }

  // Exclusive mode; if set to 2, then set based on nJet/nJetMax
  if (exclusiveMode == 2) {

    // Inclusive if nJet == nJetMax, exclusive otherwise
    int nParton = origTypeIdx[0].size();
    exclusive = (nParton == nJetMax) ? false : true;

  // Otherwise, just set as given
  } else {
    exclusive = (exclusiveMode == 0) ? false : true;
  }

  // Extract partons from hardest subsystem + ISR + FSR only into
  // workEvent. Note no resonance showers or MPIs.
  subEvent(event);
}

Member Data Documentation

double JetMatchingMadgraph::clFact

protected

Definition at line 241 of file JetMatchingPy8Internal.h.

vector<double> JetMatchingMadgraph::DJR

protected

Definition at line 247 of file JetMatchingPy8Internal.h.

Referenced by ClearDJR(), GetDJR(), and SetDJR().

bool JetMatchingMadgraph::doFxFx

protected

Definition at line 242 of file JetMatchingPy8Internal.h.

vector<int> JetMatchingMadgraph::nME

protected

Definition at line 248 of file JetMatchingPy8Internal.h.

Referenced by ClearnME(), nMEPartons(), and SetnME().

int JetMatchingMadgraph::nPartonsNow

protected

Definition at line 243 of file JetMatchingPy8Internal.h.

int JetMatchingMadgraph::nQmatch

protected

Definition at line 240 of file JetMatchingPy8Internal.h.

vector<int> JetMatchingMadgraph::origTypeIdx[3]

protected

Definition at line 239 of file JetMatchingPy8Internal.h.

Referenced by SetnME().

double JetMatchingMadgraph::qCut

protected

Definition at line 241 of file JetMatchingPy8Internal.h.

double JetMatchingMadgraph::qCutME

protected

Definition at line 244 of file JetMatchingPy8Internal.h.

double JetMatchingMadgraph::qCutMESq

protected

Definition at line 244 of file JetMatchingPy8Internal.h.

double JetMatchingMadgraph::qCutSq

protected

Definition at line 241 of file JetMatchingPy8Internal.h.

Pythia8::SlowJet* JetMatchingMadgraph::slowJetDJR

Definition at line 183 of file JetMatchingPy8Internal.h.

Referenced by SetDJR().