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);
214 jetsToken_(consumes<edm::
View<
reco::Jet> >( iConfig.getParameter<edm::InputTag>(
"jets")) ),
215 groomedJetsToken_(mayConsume<edm::
View<
reco::Jet> >( iConfig.exists(
"groomedJets") ? iConfig.getParameter<edm::InputTag>(
"groomedJets") : edm::InputTag() )),
216 subjetsToken_(mayConsume<edm::
View<
reco::Jet> >( iConfig.exists(
"subjets") ? iConfig.getParameter<edm::InputTag>(
"subjets") : edm::InputTag() )),
220 leptonsToken_(mayConsume<
reco::
GenParticleRefVector>( iConfig.exists(
"leptons") ? iConfig.getParameter<edm::InputTag>(
"leptons") : edm::InputTag() )),
221 jetAlgorithm_(iConfig.getParameter<std::
string>(
"jetAlgorithm")),
222 rParam_(iConfig.getParameter<double>(
"rParam")),
224 ghostRescaling_(iConfig.exists(
"ghostRescaling") ? iConfig.getParameter<double>(
"ghostRescaling") : 1
e-18),
225 relPtTolerance_(iConfig.exists(
"relPtTolerance") ? iConfig.getParameter<double>(
"relPtTolerance") : 1
e-03),
226 hadronFlavourHasPriority_(iConfig.getParameter<bool>(
"hadronFlavourHasPriority")),
227 useSubjets_(iConfig.exists(
"groomedJets") && iConfig.exists(
"subjets")),
228 useLeptons_(iConfig.exists(
"leptons"))
232 produces<reco::JetFlavourInfoMatchingCollection>();
234 produces<reco::JetFlavourInfoMatchingCollection>(
"SubJets");
244 throw cms::Exception(
"InvalidJetAlgorithm") <<
"Jet clustering algorithm is invalid: " <<
jetAlgorithm_ <<
", use CambridgeAachen | Kt | AntiKt" << std::endl;
290 std::auto_ptr<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfos;
295 std::vector<fastjet::PseudoJet> fjInputs;
299 std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
300 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
301 for( m = constituents.begin(); m != constituents.end(); ++
m )
304 if(constit->pt() == 0)
306 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
309 fjInputs.push_back(fastjet::PseudoJet(constit->px(),constit->py(),constit->pz(),constit->energy()));
325 std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(
fjClusterSeq_->inclusive_jets(
jetPtMin_) );
327 if( inclusiveJets.size() < jets->size() )
328 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.";
331 std::vector<int> reclusteredIndices;
335 std::vector<int> groomedIndices;
338 if( groomedJets->size() > jets->size() )
339 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.";
345 std::vector<std::vector<int> > subjetIndices;
348 matchSubjets(groomedIndices,groomedJets,subjets,subjetIndices);
352 for(
size_t i=0;
i<jets->size(); ++
i)
360 if( reclusteredIndices.at(
i) < 0 )
363 (*jetFlavourInfos)[jets->refAt(
i)] =
reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
368 if( ( fabs( inclusiveJets.at(reclusteredIndices.at(
i)).
pt() - jets->at(
i).pt() ) / jets->at(
i).pt() ) >
relPtTolerance_ )
370 if( jets->at(
i).pt() < 10. )
371 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"
372 <<
"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"
373 <<
"Since the mismatch is at low Pt (Pt<10 GeV), it is ignored and only a warning is issued.\n"
374 <<
"\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.";
376 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"
377 <<
"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"
378 <<
"\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.";
382 std::vector<fastjet::PseudoJet> constituents = fastjet::sorted_by_pt( inclusiveJets.at(reclusteredIndices.at(
i)).constituents() );
385 for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
387 if( !it->has_user_info() )
continue;
407 int hadronFlavour = 0;
408 int partonFlavour = 0;
411 setFlavours(clusteredbHadrons, clusteredcHadrons, clusteredPartons, hadronFlavour, partonFlavour);
414 (*jetFlavourInfos)[jets->refAt(
i)] =
reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, hadronFlavour, partonFlavour);
420 if( subjetIndices.at(
i).size()==0 )
continue;
429 assignToSubjets(clusteredbHadrons, subjets, subjetIndices.at(
i), assignedbHadrons);
431 assignToSubjets(clusteredcHadrons, subjets, subjetIndices.at(
i), assignedcHadrons);
433 assignToSubjets(clusteredPartons, subjets, subjetIndices.at(
i), assignedPartons);
436 assignToSubjets(clusteredLeptons, subjets, subjetIndices.at(
i), assignedLeptons);
439 for(
size_t sj=0; sj<subjetIndices.at(
i).size(); ++sj)
441 int subjetHadronFlavour = 0;
442 int subjetPartonFlavour = 0;
445 setFlavours(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
448 (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(
i).at(sj))] =
reco::JetFlavourInfo(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), assignedLeptons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
454 iEvent.
put( jetFlavourInfos );
457 iEvent.
put( subjetFlavourInfos,
"SubJets" );
463 const double ghostRescaling,
464 const bool isHadron,
const bool isbHadron,
const bool isParton,
const bool isLepton,
465 std::vector<fastjet::PseudoJet>& constituents)
472 edm::LogInfo(
"NullTransverseMomentum") <<
"dropping input ghost candidate with pt=0";
475 fastjet::PseudoJet
p((*it)->px(),(*it)->py(),(*it)->pz(),(*it)->energy());
477 p.set_user_info(
new GhostInfo(isHadron, isbHadron, isParton, isLepton, *it));
478 constituents.push_back(p);
485 const std::vector<fastjet::PseudoJet>& reclusteredJets,
486 std::vector<int>& matchedIndices)
488 std::vector<bool> matchedLocks(reclusteredJets.size(),
false);
490 for(
size_t j=0;
j<
jets->size(); ++
j)
492 double matchedDR2 = 1e9;
495 for(
size_t rj=0; rj<reclusteredJets.size(); ++rj)
497 if( matchedLocks.at(rj) )
continue;
499 double tempDR2 =
reco::deltaR2(
jets->at(
j).rapidity(),
jets->at(
j).phi(), reclusteredJets.at(rj).rapidity(), reclusteredJets.at(rj).phi_std() );
500 if( tempDR2 < matchedDR2 )
502 matchedDR2 = tempDR2;
511 edm::LogError(
"JetMatchingFailed") <<
"Matched reclustered jet " << matchedIdx <<
" and original jet " <<
j <<
" are separated by dR=" <<
sqrt(matchedDR2) <<
" which is greater than the jet size R=" << rParam_ <<
".\n"
512 <<
"This is not expected so please check that the jet algorithm and jet size match those used for the original jet collection.";
515 matchedLocks.at(matchedIdx) =
true;
518 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.";
520 matchedIndices.push_back(matchedIdx);
528 std::vector<int>& matchedIndices)
530 std::vector<bool> jetLocks(
jets->size(),
false);
531 std::vector<int> jetIndices;
533 for(
size_t gj=0; gj<groomedJets->size(); ++gj)
535 double matchedDR2 = 1e9;
538 if( groomedJets->at(gj).pt()>0. )
540 for(
size_t j=0;
j<
jets->size(); ++
j)
542 if( jetLocks.at(
j) )
continue;
544 double tempDR2 =
reco::deltaR2(
jets->at(
j).rapidity(),
jets->at(
j).phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi() );
545 if( tempDR2 < matchedDR2 )
547 matchedDR2 = tempDR2;
557 edm::LogWarning(
"MatchedJetsFarApart") <<
"Matched groomed jet " << gj <<
" and original jet " << matchedIdx <<
" are separated by dR=" <<
sqrt(matchedDR2) <<
" which is greater than the jet size R=" << rParam_ <<
".\n"
558 <<
"This is not expected so the matching of these two jets has been discarded. Please check that the two jet collections belong to each other.";
562 jetLocks.at(matchedIdx) =
true;
564 jetIndices.push_back(matchedIdx);
567 for(
size_t j=0;
j<
jets->size(); ++
j)
569 std::vector<int>::iterator matchedIndex =
std::find( jetIndices.begin(), jetIndices.end(),
j );
571 matchedIndices.push_back( matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(),matchedIndex) : -1 );
580 std::vector<std::vector<int> >& matchedIndices)
582 for(
size_t g=0;
g<groomedIndices.size(); ++
g)
584 std::vector<int> subjetIndices;
586 if( groomedIndices.at(
g)>=0 )
588 for(
size_t s=0;
s<groomedJets->at(groomedIndices.at(
g)).numberOfDaughters(); ++
s)
592 for(
size_t sj=0; sj<subjets->size(); ++sj)
596 subjetIndices.push_back(sj);
602 if( subjetIndices.size() == 0 )
603 edm::LogError(
"SubjetMatchingFailed") <<
"Matching subjets to original jets failed. Please check that the groomed jet and subjet collections belong to each other.";
605 matchedIndices.push_back(subjetIndices);
608 matchedIndices.push_back(subjetIndices);
628 if( hardestParton.
isNull() )
629 hardestParton = (*it);
631 if( hardestLightParton.
isNull() )
634 hardestLightParton = (*it);
637 if( flavourParton.
isNull() && (
abs( (*it)->pdgId() ) == 4 ) )
638 flavourParton = (*it);
640 if(
abs( (*it)->pdgId() ) == 5 )
642 if( flavourParton.
isNull() )
643 flavourParton = (*it);
644 else if(
abs( flavourParton->pdgId() ) != 5 )
645 flavourParton = (*it);
650 if( clusteredbHadrons.
size()>0 )
652 else if( clusteredcHadrons.
size()>0 && clusteredbHadrons.
size()==0 )
655 if( flavourParton.
isNull() )
657 if( hardestParton.
isNonnull() ) partonFlavour = hardestParton->pdgId();
660 partonFlavour = flavourParton->pdgId();
665 if( hadronFlavour==0 && (
abs(partonFlavour)==4 ||
abs(partonFlavour)==5) )
666 partonFlavour = ( hardestLightParton.
isNonnull() ? hardestLightParton->pdgId() : 0 );
667 else if( hadronFlavour!=0 &&
abs(partonFlavour)!=hadronFlavour )
668 partonFlavour = hadronFlavour;
676 const std::vector<int>& subjetIndices,
677 std::vector<reco::GenParticleRefVector>& assignedParticles)
682 std::vector<double> dR2toSubjets;
684 for(
size_t sj=0; sj<subjetIndices.size(); ++sj)
685 dR2toSubjets.push_back(
reco::deltaR2( (*it)->rapidity(), (*it)->phi(), subjets->at(subjetIndices.at(sj)).rapidity(), subjets->at(subjetIndices.at(sj)).
phi() ) );
688 int closestSubjetIdx = std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
690 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)
double deltaR2(const T1 &t1, const T2 &t2)
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)
const double relPtTolerance_
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_