98 #include "fastjet/JetDefinition.hh"
99 #include "fastjet/ClusterSequence.hh"
100 #include "fastjet/Selector.hh"
101 #include "fastjet/PseudoJet.hh"
107 typedef boost::shared_ptr<fastjet::JetDefinition>
JetDefPtr;
112 class GhostInfo :
public fastjet::PseudoJet::UserInfoBase{
122 if( isHadron )
m_type |= ( 1 << 0 );
123 if( isbHadron )
m_type |= ( 1 << 1 );
124 if( isParton )
m_type |= ( 1 << 2 );
125 if( isLepton )
m_type |= ( 1 << 3 );
157 const double ghostRescaling,
158 const bool isHadron,
const bool isbHadron,
const bool isParton,
const bool isLepton,
159 std::vector<fastjet::PseudoJet>& constituents);
162 const std::vector<fastjet::PseudoJet>& matchedJets,
163 std::vector<int>& matchedIndices);
166 std::vector<int>& matchedIndices);
167 void matchSubjets(
const std::vector<int>& groomedIndices,
170 std::vector<std::vector<int> >& matchedIndices);
180 const std::vector<int>& subjetIndices,
181 std::vector<reco::GenParticleRefVector>& assignedParticles);
213 jetsToken_(consumes<edm::
View<
reco::Jet> >( iConfig.getParameter<edm::InputTag>(
"jets")) ),
214 groomedJetsToken_(mayConsume<edm::
View<
reco::Jet> >( iConfig.exists(
"groomedJets") ? iConfig.getParameter<edm::InputTag>(
"groomedJets") : edm::InputTag() )),
215 subjetsToken_(mayConsume<edm::
View<
reco::Jet> >( iConfig.exists(
"subjets") ? iConfig.getParameter<edm::InputTag>(
"subjets") : edm::InputTag() )),
219 leptonsToken_(mayConsume<
reco::
GenParticleRefVector>( iConfig.exists(
"leptons") ? iConfig.getParameter<edm::InputTag>(
"leptons") : edm::InputTag() )),
220 jetAlgorithm_(iConfig.getParameter<std::
string>(
"jetAlgorithm")),
221 rParam_(iConfig.getParameter<double>(
"rParam")),
223 ghostRescaling_(iConfig.exists(
"ghostRescaling") ? iConfig.getParameter<double>(
"ghostRescaling") : 1
e-18),
224 hadronFlavourHasPriority_(iConfig.getParameter<bool>(
"hadronFlavourHasPriority")),
225 useSubjets_(iConfig.exists(
"groomedJets") && iConfig.exists(
"subjets")),
226 useLeptons_(iConfig.exists(
"leptons"))
230 produces<reco::JetFlavourInfoMatchingCollection>();
232 produces<reco::JetFlavourInfoMatchingCollection>(
"SubJets");
242 throw cms::Exception(
"InvalidJetAlgorithm") <<
"Jet clustering algorithm is invalid: " <<
jetAlgorithm_ <<
", use CambridgeAachen | Kt | AntiKt" << std::endl;
288 std::auto_ptr<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfos;
293 std::vector<fastjet::PseudoJet> fjInputs;
297 std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
298 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
299 for( m = constituents.begin(); m != constituents.end(); ++
m )
302 if(constit->pt() == 0)
304 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
307 fjInputs.push_back(fastjet::PseudoJet(constit->px(),constit->py(),constit->pz(),constit->energy()));
323 std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(
fjClusterSeq_->inclusive_jets(
jetPtMin_) );
325 if( inclusiveJets.size() < jets->size() )
326 edm::LogError(
"TooFewReclusteredJets") <<
"There are fewer reclustered (" << inclusiveJets.size() <<
") than original jets (" << jets->size() <<
"). Please check that the jet algorithm and jet size match those used for the original jet collection.";
329 std::vector<int> reclusteredIndices;
333 std::vector<int> groomedIndices;
336 if( groomedJets->size() > jets->size() )
337 edm::LogError(
"TooManyGroomedJets") <<
"There are more groomed (" << groomedJets->size() <<
") than original jets (" << jets->size() <<
"). Please check that the two jet collections belong to each other.";
343 std::vector<std::vector<int> > subjetIndices;
346 matchSubjets(groomedIndices,groomedJets,subjets,subjetIndices);
350 for(
size_t i=0;
i<jets->size(); ++
i)
353 if( ( fabs( inclusiveJets.at(reclusteredIndices.at(
i)).
pt() - jets->at(
i).pt() ) / jets->at(
i).pt() ) > 1
e-3 )
355 if( jets->at(
i).pt() < 10. )
356 edm::LogWarning(
"JetPtMismatchAtLowPt") <<
"The reclustered and original jet " <<
i <<
" have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << jets->at(
i).pt() <<
" GeV, respectively).\n"
357 <<
"Please check that the jet algorithm and jet size match those used for the original jet collection and also make sure the original jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n"
358 <<
"Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n"
359 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
361 edm::LogError(
"JetPtMismatch") <<
"The reclustered and original jet " <<
i <<
" have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << jets->at(
i).pt() <<
" GeV, respectively).\n"
362 <<
"Please check that the jet algorithm and jet size match those used for the original jet collection and also make sure the original jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n"
363 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
372 std::vector<fastjet::PseudoJet> constituents = fastjet::sorted_by_pt( inclusiveJets.at(reclusteredIndices.at(
i)).constituents() );
375 for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
377 if( !it->has_user_info() )
continue;
397 int hadronFlavour = 0;
398 int partonFlavour = 0;
401 setFlavours(clusteredbHadrons, clusteredcHadrons, clusteredPartons, hadronFlavour, partonFlavour);
404 (*jetFlavourInfos)[jets->refAt(
i)] =
reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, hadronFlavour, partonFlavour);
409 if( subjetIndices.at(
i).size()==0 )
continue;
418 assignToSubjets(clusteredbHadrons, subjets, subjetIndices.at(
i), assignedbHadrons);
420 assignToSubjets(clusteredcHadrons, subjets, subjetIndices.at(
i), assignedcHadrons);
422 assignToSubjets(clusteredPartons, subjets, subjetIndices.at(
i), assignedPartons);
425 assignToSubjets(clusteredLeptons, subjets, subjetIndices.at(
i), assignedLeptons);
428 for(
size_t sj=0; sj<subjetIndices.at(
i).size(); ++sj)
430 int subjetHadronFlavour = 0;
431 int subjetPartonFlavour = 0;
434 setFlavours(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
437 (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(
i).at(sj))] =
reco::JetFlavourInfo(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), assignedLeptons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
443 iEvent.
put( jetFlavourInfos );
446 iEvent.
put( subjetFlavourInfos,
"SubJets" );
452 const double ghostRescaling,
453 const bool isHadron,
const bool isbHadron,
const bool isParton,
const bool isLepton,
454 std::vector<fastjet::PseudoJet>& constituents)
459 fastjet::PseudoJet
p((*it)->px(),(*it)->py(),(*it)->pz(),(*it)->energy());
461 p.set_user_info(
new GhostInfo(isHadron, isbHadron, isParton, isLepton, *it));
462 constituents.push_back(p);
469 const std::vector<fastjet::PseudoJet>& reclusteredJets,
470 std::vector<int>& matchedIndices)
472 std::vector<bool> matchedLocks(reclusteredJets.size(),
false);
474 for(
size_t j=0;
j<
jets->size(); ++
j)
476 double matchedDR = 1e9;
479 for(
size_t rj=0; rj<reclusteredJets.size(); ++rj)
481 if( matchedLocks.at(rj) )
continue;
483 double tempDR =
reco::deltaR(
jets->at(
j).rapidity(),
jets->at(
j).phi(), reclusteredJets.at(rj).rapidity(), reclusteredJets.at(rj).phi_std() );
484 if( tempDR < matchedDR )
491 if( matchedIdx>=0 ) matchedLocks.at(matchedIdx) =
true;
492 matchedIndices.push_back(matchedIdx);
495 if(
std::find( matchedIndices.begin(), matchedIndices.end(), -1 ) != matchedIndices.end() )
496 edm::LogError(
"JetMatchingFailed") <<
"Matching reclustered to original jets failed. Please check that the jet algorithm and jet size match those used for the original jet collection.";
503 std::vector<int>& matchedIndices)
505 std::vector<bool> jetLocks(
jets->size(),
false);
506 std::vector<int> jetIndices;
508 for(
size_t gj=0; gj<groomedJets->size(); ++gj)
510 double matchedDR = 1e9;
513 for(
size_t j=0;
j<
jets->size(); ++
j)
515 if( jetLocks.at(
j) )
continue;
517 double tempDR =
reco::deltaR(
jets->at(
j).rapidity(),
jets->at(
j).phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi() );
518 if( tempDR < matchedDR )
525 if( matchedIdx>=0 ) jetLocks.at(matchedIdx) =
true;
526 jetIndices.push_back(matchedIdx);
529 if(
std::find( jetIndices.begin(), jetIndices.end(), -1 ) != jetIndices.end() )
530 edm::LogError(
"JetMatchingFailed") <<
"Matching groomed to original jets failed. Please check that the two jet collections belong to each other.";
532 for(
size_t j=0;
j<
jets->size(); ++
j)
534 std::vector<int>::iterator matchedIndex =
std::find( jetIndices.begin(), jetIndices.end(),
j );
536 matchedIndices.push_back( matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(),matchedIndex) : -1 );
545 std::vector<std::vector<int> >& matchedIndices)
547 for(
size_t g=0;
g<groomedIndices.size(); ++
g)
549 std::vector<int> subjetIndices;
551 if( groomedIndices.at(
g)>=0 )
553 for(
size_t s=0;
s<groomedJets->at(groomedIndices.at(
g)).numberOfDaughters(); ++
s)
557 for(
size_t sj=0; sj<subjets->size(); ++sj)
561 subjetIndices.push_back(sj);
567 if( subjetIndices.size() == 0 )
568 edm::LogError(
"SubjetMatchingFailed") <<
"Matching subjets to original jets failed. Please check that the groomed jet and subjet collections belong to each other.";
570 matchedIndices.push_back(subjetIndices);
573 matchedIndices.push_back(subjetIndices);
593 if( hardestParton.
isNull() )
594 hardestParton = (*it);
596 if( hardestLightParton.
isNull() )
599 hardestLightParton = (*it);
602 if( flavourParton.
isNull() && (
abs( (*it)->pdgId() ) == 4 ) )
603 flavourParton = (*it);
605 if(
abs( (*it)->pdgId() ) == 5 )
607 if( flavourParton.
isNull() )
608 flavourParton = (*it);
609 else if(
abs( flavourParton->pdgId() ) != 5 )
610 flavourParton = (*it);
615 if( clusteredbHadrons.
size()>0 )
617 else if( clusteredcHadrons.
size()>0 && clusteredbHadrons.
size()==0 )
620 if( flavourParton.
isNull() )
622 if( hardestParton.
isNonnull() ) partonFlavour = hardestParton->pdgId();
625 partonFlavour = flavourParton->pdgId();
630 if( hadronFlavour==0 && (
abs(partonFlavour)==4 ||
abs(partonFlavour)==5) )
631 partonFlavour = ( hardestLightParton.
isNonnull() ? hardestLightParton->pdgId() : 0 );
632 else if( hadronFlavour!=0 &&
abs(partonFlavour)!=hadronFlavour )
633 partonFlavour = hadronFlavour;
641 const std::vector<int>& subjetIndices,
642 std::vector<reco::GenParticleRefVector>& assignedParticles)
647 std::vector<double> dR2toSubjets;
649 for(
size_t sj=0; sj<subjetIndices.size(); ++sj)
650 dR2toSubjets.push_back(
reco::deltaR2( (*it)->rapidity(), (*it)->phi(), subjets->at(subjetIndices.at(sj)).rapidity(), subjets->at(subjetIndices.at(sj)).
phi() ) );
653 int closestSubjetIdx = std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
655 assignedParticles.at(closestSubjetIdx).push_back( *it );
const reco::GenParticleRef & particleRef() const
bool isLightParton(const reco::Candidate &c)
boost::shared_ptr< fastjet::JetDefinition > JetDefPtr
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
virtual void produce(edm::Event &, const edm::EventSetup &)
virtual void endLuminosityBlock(edm::LuminosityBlock &, edm::EventSetup const &)
void assignToSubjets(const reco::GenParticleRefVector &clusteredParticles, const edm::Handle< edm::View< reco::Jet > > &subjets, const std::vector< int > &subjetIndices, std::vector< reco::GenParticleRefVector > &assignedParticles)
Clusters hadrons, partons, and jet contituents to determine the jet flavour.
const bool hadronFlavourHasPriority_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
boost::shared_ptr< fastjet::ClusterSequence > ClusterSequencePtr
#define DEFINE_FWK_MODULE(type)
void insertGhosts(const edm::Handle< reco::GenParticleRefVector > &particles, const double ghostRescaling, const bool isHadron, const bool isbHadron, const bool isParton, const bool isLepton, std::vector< fastjet::PseudoJet > &constituents)
GhostInfo(const bool &isHadron, const bool &isbHadron, const bool &isParton, const bool &isLepton, const reco::GenParticleRef &particleRef)
const reco::GenParticleRef m_particleRef
void setFlavours(const reco::GenParticleRefVector &clusteredbHadrons, const reco::GenParticleRefVector &clusteredcHadrons, const reco::GenParticleRefVector &clusteredPartons, int &hadronFlavour, int &partonFlavour)
const edm::EDGetTokenT< reco::GenParticleRefVector > bHadronsToken_
bool isLepton(const Candidate &part)
void matchGroomedJets(const edm::Handle< edm::View< reco::Jet > > &jets, const edm::Handle< edm::View< reco::Jet > > &matchedJets, std::vector< int > &matchedIndices)
const_iterator end() const
Termination of iteration.
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
const double ghostRescaling_
const edm::EDGetTokenT< reco::GenParticleRefVector > partonsToken_
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
const bool isParton() const
const_iterator begin() const
Initialize an iterator over the RefVector.
JetDefPtr fjJetDefinition_
const bool isLepton() const
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const edm::EDGetTokenT< reco::GenParticleRefVector > cHadronsToken_
bool isNonnull() const
Checks for non-null.
bool isNull() const
Checks for null.
void addDefault(ParameterSetDescription const &psetDescription)
virtual void beginRun(edm::Run &, edm::EventSetup const &)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Class storing the jet flavour information.
Abs< T >::type abs(const T &t)
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
virtual void endRun(edm::Run &, edm::EventSetup const &)
JetFlavourClustering(const edm::ParameterSet &)
void matchReclusteredJets(const edm::Handle< edm::View< reco::Jet > > &jets, const std::vector< fastjet::PseudoJet > &matchedJets, std::vector< int > &matchedIndices)
const std::string jetAlgorithm_
const edm::EDGetTokenT< edm::View< reco::Jet > > groomedJetsToken_
edm::RefVector< GenParticleCollection > GenParticleRefVector
vector of reference to GenParticle in the same collection
const edm::EDGetTokenT< edm::View< reco::Jet > > jetsToken_
bool isParton(const reco::Candidate &c)
T1 deltaR2(T1 eta1, T2 phi1, T3 eta2, T4 phi2)
const bool isbHadron() const
virtual void beginLuminosityBlock(edm::LuminosityBlock &, edm::EventSetup const &)
const bool isHadron() const
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
size_type size() const
Size of the RefVector.
const edm::EDGetTokenT< reco::GenParticleRefVector > leptonsToken_
void matchSubjets(const std::vector< int > &groomedIndices, const edm::Handle< edm::View< reco::Jet > > &groomedJets, const edm::Handle< edm::View< reco::Jet > > &subjets, std::vector< std::vector< int > > &matchedIndices)
ClusterSequencePtr fjClusterSeq_
const edm::EDGetTokenT< edm::View< reco::Jet > > subjetsToken_