#include <JetMatchingMGFastJet.h>

Inheritance diagram for gen::JetMatchingMGFastJet:

Public Member Functions
const std::vector< int > *	getPartonList ()
	JetMatchingMGFastJet (const edm::ParameterSet &params)
	~JetMatchingMGFastJet ()
Protected Member Functions
void	beforeHadronisation (const lhef::LHEEvent *)
void	beforeHadronisationExec ()
virtual void	init (const lhef::LHERunInfo *runInfo)
bool	initAfterBeams ()
int	match (const lhef::LHEEvent partonLevel, const std::vector< fastjet::PseudoJet > jetInput)
Private Types
enum	partonTypes { ID_TOP = 6, ID_GLUON = 21, ID_PHOTON = 22 }
enum	vetoStatus { NONE, LESS_JETS, MORE_JETS, HARD_JET, UNMATCHED_PARTON }
Private Attributes
double	clFact
double	coneRadius
bool	doMerge
double	etaJetMax
std::vector< fastjet::PseudoJet >	fClusJets
int	fDJROutFlag
std::ofstream	fDJROutput
bool	fExcLocal
bool	fIsInit
fastjet::JetDefinition *	fJetFinder
std::vector< fastjet::PseudoJet >	fPtSortedJets
int	jetAlgoPower
int	nJetMax
int	nJetMin
int	nQmatch
double	qCut
double	qCutSq
std::vector< int >	typeIdx [3]

Detailed Description

Definition at line 36 of file JetMatchingMGFastJet.h.

Member Enumeration Documentation

enum gen::JetMatchingMGFastJet::partonTypes [private]

Enumerator:

ID_TOP
ID_GLUON
ID_PHOTON

Definition at line 61 of file JetMatchingMGFastJet.h.

{ ID_TOP=6, ID_GLUON=21, ID_PHOTON=22 };

enum gen::JetMatchingMGFastJet::vetoStatus [private]

Enumerator:

NONE
LESS_JETS
MORE_JETS
HARD_JET
UNMATCHED_PARTON

Definition at line 60 of file JetMatchingMGFastJet.h.

{ NONE, LESS_JETS, MORE_JETS, HARD_JET, UNMATCHED_PARTON };

Constructor & Destructor Documentation

gen::JetMatchingMGFastJet::JetMatchingMGFastJet ( const edm::ParameterSet & params ) [inline]

Definition at line 41 of file JetMatchingMGFastJet.h.

References fDJROutFlag, and edm::ParameterSet::getParameter().

                                                          : JetMatchingMadgraph(params), 
                                                              fJetFinder(0),  
                                                              fIsInit(false) 
                                                              { fDJROutFlag = params.getParameter<int>("outTree_flag"); }

gen::JetMatchingMGFastJet::~JetMatchingMGFastJet ( ) [inline]

Definition at line 45 of file JetMatchingMGFastJet.h.

References fJetFinder.

{ if (fJetFinder) delete fJetFinder; }

Member Function Documentation

void gen::JetMatchingMGFastJet::beforeHadronisation ( const lhef::LHEEvent * lhee ) [protected, virtual]

Reimplemented from gen::JetMatchingMadgraph.

Definition at line 76 of file JetMatchingMGFastJet.cc.

References Exception, gen::JetMatchingMadgraph::exclusive, fExcLocal, lhef::LHEEvent::getHEPEUP(), i, ID_GLUON, ID_TOP, lhef::HEPEUP::IDUP, customizeTrackingMonitorSeedNumber::idx, lhef::HEPEUP::ISTUP, lhef::HEPEUP::MOTHUP, nJetMax, nQmatch, lhef::HEPEUP::NUP, gen::JetMatchingMadgraph::runInitialized, gen::JetMatchingMadgraph::soup, and typeIdx.

{

   if (!runInitialized)
                throw cms::Exception("Generator|PartonShowerVeto")
                        << "Run not initialized in JetMatchingMadgraph"
                        << std::endl;
   
   for (int i = 0; i < 3; i++) 
   {
    typeIdx[i].clear();
   }

   // 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])
   //
   const lhef::HEPEUP& hepeup = *lhee->getHEPEUP(); 
   int idx = 2; 
   for ( int i=0; i < hepeup.NUP; i++ )
   {
      if ( hepeup.ISTUP[i] < 0 ) continue;
      if ( hepeup.MOTHUP[i].first != 1 && hepeup.MOTHUP[i].second !=2 ) continue; // this way we skip entries that come
                                                                                  // from resonance decays;
                                                                                  // we only take those that descent
                                                                                  // directly from "incoming partons"
      idx = 2;
      if ( hepeup.IDUP[i] == ID_GLUON || (fabs(hepeup.IDUP[i]) <= nQmatch) ) // light jet
         // light jet
         idx = 0;
      else if ( fabs(hepeup.IDUP[i]) > nQmatch && fabs(hepeup.IDUP[i]) <= ID_TOP) // heavy jet
         idx = 1;
      // Store
      typeIdx[idx].push_back(i); 
   } 
   
   // NOTE: In principle, I should use exclusive, inclusive, or soup !!!   
   // should be like this:
   if ( soup )
   {
      int NPartons = typeIdx[0].size(); 
      fExcLocal = ( NPartons < nJetMax );
   }
   else
      fExcLocal = exclusive;
        
   return;

}

void gen::JetMatchingMGFastJet::beforeHadronisationExec ( ) [inline, protected, virtual]

Reimplemented from gen::JetMatchingMadgraph.

Definition at line 54 of file JetMatchingMGFastJet.h.

{ return; }

const std::vector<int>* gen::JetMatchingMGFastJet::getPartonList ( ) [inline, virtual]

Reimplemented from gen::JetMatching.

Definition at line 47 of file JetMatchingMGFastJet.h.

References typeIdx.

{ return typeIdx; }

virtual void gen::JetMatchingMGFastJet::init ( const lhef::LHERunInfo * runInfo ) [inline, protected, virtual]

Reimplemented from gen::JetMatchingMadgraph.

Definition at line 50 of file JetMatchingMGFastJet.h.

References fIsInit, and initAfterBeams().

                                                         { if (fIsInit) return; JetMatchingMadgraph::init(runInfo); 
                                                             initAfterBeams(); fIsInit=true; return; }

bool gen::JetMatchingMGFastJet::initAfterBeams ( ) [protected, virtual]

Reimplemented from gen::JetMatching.

Definition at line 32 of file JetMatchingMGFastJet.cc.

References clFact, coneRadius, doMerge, gen::MEMAIN::etaclmax, etaJetMax, fClusJets, fExcLocal, fIsInit, fJetFinder, fPtSortedJets, gen::UPPRIV::ickkw, jetAlgoPower, gen::MEMAIN::maxjets, gen::memain_, gen::MEMAIN::minjets, nJetMax, nJetMin, nQmatch, gen::MEMAIN::nqmatch, funct::pow(), qCut, gen::MEMAIN::qcut, qCutSq, and gen::uppriv_.

Referenced by init().

{
  
  if ( fIsInit ) return true;
  
  //
  doMerge = uppriv_.ickkw; 
  // doMerge = true;
  qCut    = memain_.qcut; // 
  nQmatch        = memain_.nqmatch;
  clFact         = 1.; // Default value 
                      // NOTE: ME2pythia seems to default to 1.5 - need to check !!!
                      // in general, needs to read key ALPSFACT from LHE file - fix CMSSW code !!!
  nJetMin        = memain_.minjets;
  nJetMax        = memain_.maxjets;

  etaJetMax      = memain_.etaclmax; 

  coneRadius     = 1.0; 
  jetAlgoPower   = 1; //   this is the kT algorithm !!! 

  // Matching procedure
  //
  qCutSq         = pow(qCut,2);  
  // this should be something like memaev_.iexc
  fExcLocal = true;


   // If not merging, then done (?????)
   //
   // if (!doMerge) return true;
  
   // Initialise chosen jet algorithm. 
   //
   fJetFinder = new fastjet::JetDefinition( fastjet::kt_algorithm, coneRadius );
   fClusJets.clear();
   fPtSortedJets.clear();  
   
   fIsInit = true;
   
   return true;

}

int gen::JetMatchingMGFastJet::match	(	const lhef::LHEEvent *	partonLevel,
		const std::vector< fastjet::PseudoJet > *	jetInput
	)		`[protected, virtual]`

Reimplemented from gen::JetMatchingMadgraph.

Definition at line 128 of file JetMatchingMGFastJet.cc.

References clFact, funct::false, fClusJets, fDJROutFlag, fDJROutput, fExcLocal, fJetFinder, fPtSortedJets, lhef::LHEEvent::getHEPEUP(), i, customizeTrackingMonitorSeedNumber::idx, LESS_JETS, max(), MORE_JETS, NONE, convertSQLitetoXML_cfg::output, funct::pow(), lhef::HEPEUP::PUP, qCut, qCutSq, mathSSE::sqrt(), typeIdx, and UNMATCHED_PARTON.

{
         
   // Number of hard partons
   //
   int NPartons = typeIdx[0].size();
   
   fClusJets.clear();
   fPtSortedJets.clear();
   
   int ClusSeqNJets = 0;

   fastjet::ClusterSequence ClusSequence( *jetInput, *fJetFinder );
   
   if ( fExcLocal )
   {
      fClusJets = ClusSequence.exclusive_jets( qCutSq );
   }
   else
   {
      fClusJets = ClusSequence.inclusive_jets( qCut ); 
   }
   
   ClusSeqNJets = fClusJets.size(); 
   
   
   if ( ClusSeqNJets < NPartons ) return LESS_JETS;
   
   double localQcutSq = qCutSq;
   
   if ( fExcLocal ) // exclusive
   {
      if( ClusSeqNJets > NPartons ) return MORE_JETS;
   }
   else // inclusive 
   {
      fPtSortedJets = fastjet::sorted_by_pt( *jetInput );
      localQcutSq = std::max( qCutSq, fPtSortedJets[0].pt2() );
      fClusJets = ClusSequence.exclusive_jets( NPartons ); // override
      ClusSeqNJets = NPartons;
   }
   
  if( clFact != 0 ) localQcutSq *= pow(clFact,2);
  
  std::vector<fastjet::PseudoJet> MatchingInput;
  
  std::vector<bool> jetAssigned;
  jetAssigned.assign( fClusJets.size(), false );

  int counter = 0;
  
  const lhef::HEPEUP& hepeup = *partonLevel->getHEPEUP();
    
  while ( counter < NPartons )
  {

     MatchingInput.clear();

     for ( int i=0; i<ClusSeqNJets; i++ )
     {
        if ( jetAssigned[i] ) continue;
        if ( i == NPartons )  break;
        //
        // this looks "awkward" but this way we do NOT pass in cluster_hist_index
        // which would actually indicate position in "history" of the "master" ClusSeq
        //
        fastjet::PseudoJet exjet = fClusJets[i];
        MatchingInput.push_back( fastjet::PseudoJet( exjet.px(), exjet.py(), exjet.pz(), exjet.e() ) );
        MatchingInput.back().set_user_index(i);
     }

     int idx = typeIdx[0][counter];
     MatchingInput.push_back( fastjet::PseudoJet( hepeup.PUP[idx][0],
                                                  hepeup.PUP[idx][1],
                                                  hepeup.PUP[idx][2],
                                                  hepeup.PUP[idx][3]) );

     //
     // in principle, one can use ClusterSequence::n_particles() 
     //
     int NBefore = MatchingInput.size();

     // Create new clustering object - which includes the 1st clustering run !!!
     // NOTE-1: it better be a new object for each try, or history will be "too crowded".
     // NOTE-2: when created, the it ALWAYS makes at least 1 clustering step, so in most
     //         cases at least 1 new jet object will be added; although in some cases 
     //         no jet is added, if the system doesn't cluster to anything (for example,
     //         input jet(s) and parton(s) are going opposite directions)
     //
     fastjet::ClusterSequence ClusSeq( MatchingInput, *fJetFinder );
     
     const std::vector<fastjet::PseudoJet>&        output = ClusSeq.jets();
     int NClusJets = output.size() - NBefore;
     
     // 
     // JVY - I think this is the right idea:
     // at least 1 (one) new jet needs to be added
     // however, need to double check details and refine, especially for inclusive mode
     //
     // 
     if ( NClusJets < 1 )
     {
        return UNMATCHED_PARTON;
     }
     //
     // very unlikely case but let's do it just to be safe
     //
     if ( NClusJets >= NBefore )
     {
        return MORE_JETS;
     }
      
     // Now browse history and see how close the clustering distance
     //
     // NOTE: Remember, there maybe more than one new jet in the list (for example,
     //       for process=2,3,...);
     //       in this case we take the ** first ** one that satisfies the distance/cut,
     //       which is ** typically ** also the best one 
     //
     bool matched = false;
     const std::vector<fastjet::ClusterSequence::history_element>&  history = ClusSeq.history();
              
     // for ( unsigned int i=nBefore; i<history.size(); i++ )
     for ( unsigned int i=NBefore; i<output.size(); i++ )
     {

        int hidx = output[i].cluster_sequence_history_index();
        double dNext = history[hidx].dij;
        if ( dNext < localQcutSq )
        {
           //
           // the way we form input, parent1 is always jet, and parent2 can be any,
           // but if it's a good match/cluster, it should point at the parton
           //
           int parent1 = history[hidx].parent1;
           int parent2 = history[hidx].parent2;
           if ( parent1 < 0 || parent2 < 0 ) break; // bad cluster, no match
           //
           // pull up jet's "global" index 
           //
           int pidx = MatchingInput[parent1].user_index();
           jetAssigned[pidx] = true; 
           matched = true;
           break;
        }
     }
     if ( !matched )
     {
        return UNMATCHED_PARTON;
     }
     
     counter++;   
  }

   // Now save some info for DJR analysis (if requested).
   // This mimics what is done in ME2pythia.f
   // Basically, NJets and matching scale for these 4 cases:
   // 1->0, 2->1, 3->2, and 4->3
   //
   if ( fDJROutFlag > 0 )
   {   
      std::vector<double> MergingScale;
      MergingScale.clear();
      for ( int nj=0; nj<4; nj++ )
      {
         double dmscale2 = ClusSequence.exclusive_dmerge( nj );
         double dmscale = sqrt( dmscale2 );
         MergingScale.push_back( dmscale );
      }
      fDJROutput.open( "events.tree", std::ios_base::app );  
      double dNJets = (double)NPartons;
      fDJROutput << " " << dNJets << " " << MergingScale[0] << " " 
                                         << MergingScale[1] << " "
                                         << MergingScale[2] << " "
                                         << MergingScale[3] << std::endl;  
      fDJROutput.close();
   }
 
  return NONE;

}