Inheritance diagram for JetPartonMatcher:

Public Member Functions
	JetPartonMatcher (const edm::ParameterSet &)
	~JetPartonMatcher ()
Private Member Functions
int	fillAlgoritDefinition (const Jet &)
int	fillPhysicsDefinition (const Jet &)
virtual void	produce (edm::Event &, const edm::EventSetup &)
Private Attributes
double	coneSizeToAssociate
bool	doPriority
edm::InputTag	m_jetsSrc
edm::InputTag	m_ParticleSrc
Handle< GenParticleRefVector >	particles
bool	physDefinition
vector< int >	priorityList
int	theHardest
int	theHeaviest
int	theNearest2
int	theNearest3

Detailed Description

Definition at line 94 of file JetPartonMatcher.cc.

Constructor & Destructor Documentation

JetPartonMatcher::JetPartonMatcher ( const edm::ParameterSet & iConfig )

Definition at line 122 of file JetPartonMatcher.cc.

References coneSizeToAssociate, doPriority, edm::ParameterSet::exists(), edm::ParameterSet::getParameter(), m_jetsSrc, m_ParticleSrc, and priorityList.

                                                                   :
  theHeaviest(0),
  theNearest2(0),
  theNearest3(0),
  theHardest(0)
{
    produces<JetMatchedPartonsCollection>();
    m_jetsSrc           = iConfig.getParameter<edm::InputTag>("jets");
    m_ParticleSrc       = iConfig.getParameter<edm::InputTag>("partons");
    coneSizeToAssociate = iConfig.getParameter<double>("coneSizeToAssociate");
    if ( iConfig.exists("doPriority") ) {
      doPriority = iConfig.getParameter<bool>("doPriority");
      priorityList = iConfig.getParameter<vector<int> >("priorityList");
    } else {
      doPriority = false;
      priorityList.clear();
    }
}

JetPartonMatcher::~JetPartonMatcher ( )

Definition at line 143 of file JetPartonMatcher.cc.

{
}

Member Function Documentation

int JetPartonMatcher::fillAlgoritDefinition ( const Jet & theJet ) [private]

Definition at line 215 of file JetPartonMatcher.cc.

References abs, coneSizeToAssociate, doPriority, spr::find(), m, reco::LeafCandidate::p4(), reco::Candidate::p4(), particles, reco::Candidate::pdgId(), priorityList, reco::Candidate::pt(), reco::Candidate::status(), theHardest, theHeaviest, and theNearest2.

Referenced by produce().

                                                               {

  int tempParticle = -1;
  int tempPartonHighestPt = -1;
  int tempNearest = -1;
  float maxPt = 0;
  float minDr = 1000;
  bool foundPriority = false;

  // Loop over the particles in question until we find a priority
  // "hit", or if we find none in the priority list (or we don't want
  // to consider priority), then loop through all particles and fill
  // standard definition. 

  // 
  // Matching:
  //
  // 1) First try to match by hand. The "priority list" is given
  //    by the user. The algorithm finds any particles in that 
  //    "priority list" that are within the specified cone size.
  //    If it finds one, it counts the object as associating to that
  //    particle.
  //    NOTE! The objects in the priority list are given in order
  //    of priority. So if the user specifies:
  //    6, 24, 21,
  //    then it will first look for top quarks, then W bosons, then gluons.
  // 2) If no priority items are found, do the default "standard"
  //    matching. 
  for( size_t m = 0; m != particles->size() && !foundPriority; ++ m ) {
    const Candidate & aParton = *(particles->at(m).get());

    // "Priority" behavoir:
    // Associate to the first particle found in the priority list, regardless
    // of delta R. 
    if ( doPriority ) {
      int ipdgId = aParton.pdgId();
      vector<int>::const_iterator ipriority = find( priorityList.begin(), priorityList.end(), abs(ipdgId) );
      // Check to see if this particle is in our priority list
      if ( ipriority != priorityList.end() ) {
        // This particle is on our priority list. If it matches,
        // we break, since we take in order of priority, not deltaR
        double dist = DeltaR( theJet.p4(), aParton.p4() );
        if( dist <= coneSizeToAssociate ) {
          tempParticle = m;
          foundPriority = true;
        }
      }
    }
    // Here we do not want to do priority matching. Ensure the "foundPriority" swtich
    // is turned off. 
    else {
      foundPriority = false;
    }

    // Default behavior: 
    // Look for partons before the color string to associate. 
    // Do this if we don't want to do priority matching, or if
    // we didn't find a priority object. 
    
    if( !foundPriority && aParton.status() != 3 ) {
      double dist = DeltaR( theJet.p4(), aParton.p4() );
      if( dist <= coneSizeToAssociate ) {
        if( dist < minDr ) {
          minDr = dist;
          tempNearest = m;
        }
        if( tempParticle == -1 && ( abs( aParton.pdgId() ) == 4 )  ) tempParticle = m;
        if(                         abs( aParton.pdgId() ) == 5    ) tempParticle = m;
        if( aParton.pt() > maxPt ) {
          maxPt = aParton.pt();
          tempPartonHighestPt = m;
        }
      }
    }
  }

  if ( foundPriority ) {
    theHeaviest = tempParticle; // The priority-matched particle
    theHardest  = -1;  //  set the rest to -1
    theNearest2 = -1;  // "                  "
  } else {
    theHeaviest = tempParticle;
    theHardest  = tempPartonHighestPt;
    theNearest2 = tempNearest;
    if ( tempParticle == -1 ) tempParticle = tempPartonHighestPt;
  }
  return tempParticle;
}

int JetPartonMatcher::fillPhysicsDefinition ( const Jet & theJet ) [private]

Definition at line 326 of file JetPartonMatcher.cc.

References abs, coneSizeToAssociate, doPriority, spr::find(), reco::flavour(), m, reco::Candidate::numberOfDaughters(), reco::LeafCandidate::p4(), reco::Candidate::p4(), particles, reco::Candidate::pdgId(), benchmark_cfg::pdgId, priorityList, reco::Candidate::status(), theHardest, theHeaviest, theNearest2, and theNearest3.

Referenced by produce().

                                                               {

  float TheBiggerConeSize = 0.7; // In HepMC it's 0.3 --> it's a mistake: value has to be 0.7
  int tempParticle = -1;
  int nInTheCone = 0;
  int tempNearest = -1;
  float minDr = 1000;
  bool foundPriority = false;

  vector<const reco::Candidate *> theContaminations;
  theContaminations.clear();

  // Loop over the particles in question until we find a priority
  // "hit", or if we find none in the priority list (or we don't want
  // to consider priority), then loop through all particles and fill
  // standard definition. 

  // 
  // Matching:
  //
  // 1) First try to match by hand. The "priority list" is given
  //    by the user. The algorithm finds any particles in that 
  //    "priority list" that are within the specified cone size.
  //    If it finds one, it counts the object as associating to that
  //    particle.
  //    NOTE! The objects in the priority list are given in order
  //    of priority. So if the user specifies:
  //    6, 24, 21,
  //    then it will first look for top quarks, then W bosons, then gluons.
  // 2) If no priority items are found, do the default "standard"
  //    matching. 
  for( size_t m = 0; m != particles->size() && !foundPriority; ++ m ) {

    const Candidate & aParticle = *(particles->at(m).get());

    // "Priority" behavoir:
    // Associate to the first particle found in the priority list, regardless
    // of delta R. 
    if ( doPriority ) {
      int ipdgId = aParticle.pdgId();
      vector<int>::const_iterator ipriority = find( priorityList.begin(), priorityList.end(), abs(ipdgId) );
      // Check to see if this particle is in our priority list
      if ( ipriority != priorityList.end() ) {
        // This particle is on our priority list. If it matches,
        // we break, since we take in order of priority, not deltaR
        double dist = DeltaR( theJet.p4(), aParticle.p4() );
        if( dist <= coneSizeToAssociate ) {
          tempParticle = m;
          foundPriority = true;
        }
      }
    }
    // Here we do not want to do priority matching. Ensure the "foundPriority" swtich
    // is turned off. 
    else {
      foundPriority = false;
    }

    // Default behavior: 
    if ( !foundPriority ) {

      // skipping all particle but udscbg (is this correct/enough?!?!)
      bool isAParton = false;
      int flavour = abs(aParticle.pdgId());
      if(flavour == 1 || 
         flavour == 2 ||
         flavour == 3 ||
         flavour == 4 ||
         flavour == 5 ||
         flavour == 21 ) isAParton = true;
      if(!isAParton) continue;
      double dist = DeltaR( theJet.p4(), aParticle.p4() );
      if( aParticle.status() == 3 && dist < minDr ) {
        minDr = dist;
        tempNearest = m;
      }
      if( aParticle.status() == 3 && dist <= coneSizeToAssociate ) {
        //cout << "particle in small cone=" << aParticle.pdgId() << endl;
        tempParticle = m;
        nInTheCone++;
      }
      // Look for heavy partons in TheBiggerConeSize now
      if( aParticle.numberOfDaughters() > 0 && aParticle.status() != 3 ) {
        if( flavour ==  1 ||
            flavour ==  2 ||
            flavour ==  3 ||
            flavour == 21 ) continue;
        if( dist < TheBiggerConeSize ) theContaminations.push_back( &aParticle );
      }
    }
  }


  // Here's the default behavior for assignment if there is no priority. 
  if ( !foundPriority ) {
    theNearest3 = tempNearest;

    if(nInTheCone != 1) return -1; // rejected --> only one initialParton requested
    if(theContaminations.size() == 0 ) return tempParticle; //no contamination
    int initialPartonFlavour = abs( (particles->at(tempParticle).get()) ->pdgId() );

    vector<const Candidate *>::const_iterator itCont = theContaminations.begin();
    for( ; itCont != theContaminations.end(); itCont++ ) {
      int contaminatingFlavour = abs( (*itCont)->pdgId() );
      if( (*itCont)->numberOfMothers()>0 && (*itCont)->mother(0) == particles->at(tempParticle).get() ) continue; // mother is the initialParton --> OK
      if( initialPartonFlavour == 4 ) {  
        if( contaminatingFlavour == 4 ) continue; // keep association --> the initialParton is a c --> the contaminated parton is a c
        tempParticle = -1; // all the other cases reject!
        return tempParticle;
      }
    } 
  }
  // If there is priority, then just set the heaviest to priority, the rest are -1. 
  else {
    theHeaviest = tempParticle; // Set the heaviest to tempParticle
    theNearest2 = -1; //  Set the rest to -1
    theNearest3 = -1; // "                  "
    theHardest = -1;  // "                  "
  }

  return tempParticle;   
}

void JetPartonMatcher::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iEs
	)		`[private, virtual]`

Implements edm::EDProducer.

Definition at line 149 of file JetPartonMatcher.cc.

References reco::LeafCandidate::eta(), fillAlgoritDefinition(), fillPhysicsDefinition(), edm::Event::getByLabel(), j, m, m_jetsSrc, m_ParticleSrc, particles, reco::LeafCandidate::pdgId(), reco::LeafCandidate::phi(), reco::LeafCandidate::pt(), edm::Event::put(), reco::LeafCandidate::status(), theHeaviest, theNearest2, and theNearest3.

{
  edm::Handle <edm::View <reco::Jet> > jets_h;
  iEvent.getByLabel(m_jetsSrc,     jets_h    );
  iEvent.getByLabel(m_ParticleSrc, particles );

  edm::LogVerbatim("JetPartonMatcher") << "=== Partons size:" << particles->size();

  for( size_t m = 0; m != particles->size(); ++ m ) {
    const GenParticle & aParton = *(particles->at(m).get());
    edm::LogVerbatim("JetPartonMatcher") <<  aParton.status() << " " <<
                                             aParton.pdgId()  << " " <<
                                             aParton.pt()     << " " << 
                                             aParton.eta()    << " " <<
                                             aParton.phi()    << endl;
  }

  auto_ptr<reco::JetMatchedPartonsCollection> jetMatchedPartons( new JetMatchedPartonsCollection(reco::JetRefBaseProd(jets_h)));
  
  for (size_t j = 0; j < jets_h->size(); j++) {

    const int theMappedPartonAlg = fillAlgoritDefinition( (*jets_h)[j] );
    const int theMappedPartonPhy = fillPhysicsDefinition( (*jets_h)[j] );

    GenParticleRef pHV;
    GenParticleRef pN2;
    GenParticleRef pN3;
    GenParticleRef pPH;
    GenParticleRef pAL;

    if(theHeaviest>=0)        pHV = particles->at( theHeaviest        );
    if(theNearest2>=0)        pN2 = particles->at( theNearest2        );
    if(theNearest3>=0)        pN3 = particles->at( theNearest3        );
    if(theMappedPartonPhy>=0) pPH = particles->at( theMappedPartonPhy );
    if(theMappedPartonAlg>=0) pAL = particles->at( theMappedPartonAlg );

    (*jetMatchedPartons)[jets_h->refAt(j)]=MatchedPartons(pHV,pN2,pN3,pPH,pAL);
  }
  
  iEvent.put(  jetMatchedPartons );

}