100 #include "fastjet/JetDefinition.hh"
101 #include "fastjet/ClusterSequence.hh"
102 #include "fastjet/Selector.hh"
103 #include "fastjet/PseudoJet.hh"
109 typedef boost::shared_ptr<fastjet::JetDefinition>
JetDefPtr;
114 class GhostInfo :
public fastjet::PseudoJet::UserInfoBase{
124 if( isHadron )
m_type |= ( 1 << 0 );
125 if( isbHadron )
m_type |= ( 1 << 1 );
126 if( isParton )
m_type |= ( 1 << 2 );
127 if( isLepton )
m_type |= ( 1 << 3 );
152 const double ghostRescaling,
153 const bool isHadron,
const bool isbHadron,
const bool isParton,
const bool isLepton,
154 std::vector<fastjet::PseudoJet>& constituents);
157 const std::vector<fastjet::PseudoJet>& matchedJets,
158 std::vector<int>& matchedIndices);
161 std::vector<int>& matchedIndices);
162 void matchSubjets(
const std::vector<int>& groomedIndices,
165 std::vector<std::vector<int> >& matchedIndices);
175 const std::vector<int>& subjetIndices,
176 std::vector<reco::GenParticleRefVector>& assignedParticles);
209 jetsToken_(consumes<edm::
View<
reco::
Jet> >( iConfig.getParameter<edm::
InputTag>(
"jets")) ),
210 groomedJetsToken_(mayConsume<edm::
View<
reco::
Jet> >( iConfig.exists(
"groomedJets") ? iConfig.getParameter<edm::
InputTag>(
"groomedJets") : edm::
InputTag() )),
211 subjetsToken_(mayConsume<edm::
View<
reco::
Jet> >( iConfig.exists(
"subjets") ? iConfig.getParameter<edm::
InputTag>(
"subjets") : edm::
InputTag() )),
216 jetAlgorithm_(iConfig.getParameter<std::
string>(
"jetAlgorithm")),
217 rParam_(iConfig.getParameter<double>(
"rParam")),
219 ghostRescaling_(iConfig.exists(
"ghostRescaling") ? iConfig.getParameter<double>(
"ghostRescaling") : 1
e-18),
220 relPtTolerance_(iConfig.exists(
"relPtTolerance") ? iConfig.getParameter<double>(
"relPtTolerance") : 1
e-03),
221 hadronFlavourHasPriority_(iConfig.getParameter<bool>(
"hadronFlavourHasPriority")),
222 useSubjets_(iConfig.exists(
"groomedJets") && iConfig.exists(
"subjets")),
223 useLeptons_(iConfig.exists(
"leptons"))
227 produces<reco::JetFlavourInfoMatchingCollection>();
229 produces<reco::JetFlavourInfoMatchingCollection>(
"SubJets");
239 throw cms::Exception(
"InvalidJetAlgorithm") <<
"Jet clustering algorithm is invalid: " <<
jetAlgorithm_ <<
", use CambridgeAachen | Kt | AntiKt" << std::endl;
285 std::auto_ptr<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfos;
290 std::vector<fastjet::PseudoJet> fjInputs;
294 std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
295 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
296 for( m = constituents.begin(); m != constituents.end(); ++
m )
299 if(constit->pt() == 0)
301 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
304 fjInputs.push_back(fastjet::PseudoJet(constit->px(),constit->py(),constit->pz(),constit->energy()));
320 std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(
fjClusterSeq_->inclusive_jets(
jetPtMin_) );
322 if( inclusiveJets.size() < jets->size() )
323 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.";
326 std::vector<int> reclusteredIndices;
330 std::vector<int> groomedIndices;
333 if( groomedJets->size() > jets->size() )
334 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.";
340 std::vector<std::vector<int> > subjetIndices;
343 matchSubjets(groomedIndices,groomedJets,subjets,subjetIndices);
347 for(
size_t i=0;
i<jets->size(); ++
i)
355 if( reclusteredIndices.at(
i) < 0 )
358 (*jetFlavourInfos)[jets->refAt(
i)] =
reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
360 else if( jets->at(
i).pt() == 0 )
362 edm::LogWarning(
"NullTransverseMomentum") <<
"The original jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
365 (*jetFlavourInfos)[jets->refAt(
i)] =
reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
371 for(
size_t sj=0; sj<subjetIndices.at(
i).size(); ++sj)
383 if( jets->at(
i).pt() < 10. )
384 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"
385 <<
"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"
386 <<
"Since the mismatch is at low Pt (Pt<10 GeV), it is ignored and only a warning is issued.\n"
387 <<
"\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.";
389 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"
390 <<
"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"
391 <<
"\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.";
395 std::vector<fastjet::PseudoJet> constituents = fastjet::sorted_by_pt( inclusiveJets.at(reclusteredIndices.at(
i)).constituents() );
398 for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
400 if( !it->has_user_info() )
continue;
420 int hadronFlavour = 0;
421 int partonFlavour = 0;
424 setFlavours(clusteredbHadrons, clusteredcHadrons, clusteredPartons, hadronFlavour, partonFlavour);
427 (*jetFlavourInfos)[jets->refAt(
i)] =
reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, hadronFlavour, partonFlavour);
433 if( subjetIndices.at(
i).size()==0 )
continue;
442 assignToSubjets(clusteredbHadrons, subjets, subjetIndices.at(
i), assignedbHadrons);
444 assignToSubjets(clusteredcHadrons, subjets, subjetIndices.at(
i), assignedcHadrons);
446 assignToSubjets(clusteredPartons, subjets, subjetIndices.at(
i), assignedPartons);
449 assignToSubjets(clusteredLeptons, subjets, subjetIndices.at(
i), assignedLeptons);
452 for(
size_t sj=0; sj<subjetIndices.at(
i).size(); ++sj)
454 int subjetHadronFlavour = 0;
455 int subjetPartonFlavour = 0;
458 setFlavours(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
461 (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(
i).at(sj))] =
reco::JetFlavourInfo(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), assignedLeptons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
467 iEvent.
put( jetFlavourInfos );
470 iEvent.
put( subjetFlavourInfos,
"SubJets" );
476 const double ghostRescaling,
477 const bool isHadron,
const bool isbHadron,
const bool isParton,
const bool isLepton,
478 std::vector<fastjet::PseudoJet>& constituents)
485 edm::LogInfo(
"NullTransverseMomentum") <<
"dropping input ghost candidate with pt=0";
488 fastjet::PseudoJet
p((*it)->px(),(*it)->py(),(*it)->pz(),(*it)->energy());
490 p.set_user_info(
new GhostInfo(isHadron, isbHadron, isParton, isLepton, *it));
491 constituents.push_back(p);
498 const std::vector<fastjet::PseudoJet>& reclusteredJets,
499 std::vector<int>& matchedIndices)
501 std::vector<bool> matchedLocks(reclusteredJets.size(),
false);
503 for(
size_t j=0;
j<
jets->size(); ++
j)
505 double matchedDR2 = 1e9;
508 for(
size_t rj=0; rj<reclusteredJets.size(); ++rj)
510 if( matchedLocks.at(rj) )
continue;
512 double tempDR2 =
reco::deltaR2(
jets->at(
j).rapidity(),
jets->at(
j).phi(), reclusteredJets.at(rj).rapidity(), reclusteredJets.at(rj).phi_std() );
513 if( tempDR2 < matchedDR2 )
515 matchedDR2 = tempDR2;
524 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"
525 <<
"This is not expected so please check that the jet algorithm and jet size match those used for the original jet collection.";
528 matchedLocks.at(matchedIdx) =
true;
531 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.";
533 matchedIndices.push_back(matchedIdx);
541 std::vector<int>& matchedIndices)
543 std::vector<bool> jetLocks(
jets->size(),
false);
544 std::vector<int> jetIndices;
546 for(
size_t gj=0; gj<groomedJets->size(); ++gj)
548 double matchedDR2 = 1e9;
551 if( groomedJets->at(gj).pt()>0. )
553 for(
size_t j=0;
j<
jets->size(); ++
j)
555 if( jetLocks.at(
j) )
continue;
557 double tempDR2 =
reco::deltaR2(
jets->at(
j).rapidity(),
jets->at(
j).phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi() );
558 if( tempDR2 < matchedDR2 )
560 matchedDR2 = tempDR2;
570 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"
571 <<
"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.";
575 jetLocks.at(matchedIdx) =
true;
577 jetIndices.push_back(matchedIdx);
580 for(
size_t j=0;
j<
jets->size(); ++
j)
582 std::vector<int>::iterator matchedIndex =
std::find( jetIndices.begin(), jetIndices.end(),
j );
584 matchedIndices.push_back( matchedIndex != jetIndices.end() ?
std::distance(jetIndices.begin(),matchedIndex) : -1 );
593 std::vector<std::vector<int> >& matchedIndices)
595 for(
size_t g=0;
g<groomedIndices.size(); ++
g)
597 std::vector<int> subjetIndices;
599 if( groomedIndices.at(
g)>=0 )
601 for(
size_t s=0;
s<groomedJets->at(groomedIndices.at(
g)).numberOfDaughters(); ++
s)
605 for(
size_t sj=0; sj<subjets->size(); ++sj)
609 subjetIndices.push_back(sj);
615 if( subjetIndices.size() == 0 )
616 edm::LogError(
"SubjetMatchingFailed") <<
"Matching subjets to original jets failed. Please check that the groomed jet and subjet collections belong to each other.";
618 matchedIndices.push_back(subjetIndices);
621 matchedIndices.push_back(subjetIndices);
641 if( hardestParton.
isNull() )
642 hardestParton = (*it);
644 if( hardestLightParton.
isNull() )
647 hardestLightParton = (*it);
650 if( flavourParton.
isNull() && (
std::abs( (*it)->pdgId() ) == 4 ) )
651 flavourParton = (*it);
653 if(
std::abs( (*it)->pdgId() ) == 5 )
655 if( flavourParton.
isNull() )
656 flavourParton = (*it);
657 else if(
std::abs( flavourParton->pdgId() ) != 5 )
658 flavourParton = (*it);
663 if( clusteredbHadrons.
size()>0 )
665 else if( clusteredcHadrons.
size()>0 && clusteredbHadrons.
size()==0 )
668 if( flavourParton.
isNull() )
670 if( hardestParton.
isNonnull() ) partonFlavour = hardestParton->pdgId();
673 partonFlavour = flavourParton->pdgId();
678 if( hadronFlavour==0 && (
std::abs(partonFlavour)==4 ||
std::abs(partonFlavour)==5) )
679 partonFlavour = ( hardestLightParton.
isNonnull() ? hardestLightParton->pdgId() : 0 );
680 else if( hadronFlavour!=0 &&
std::abs(partonFlavour)!=hadronFlavour )
681 partonFlavour = hadronFlavour;
689 const std::vector<int>& subjetIndices,
690 std::vector<reco::GenParticleRefVector>& assignedParticles)
695 std::vector<double> dR2toSubjets;
697 for(
size_t sj=0; sj<subjetIndices.size(); ++sj)
698 dR2toSubjets.push_back(
reco::deltaR2( (*it)->rapidity(), (*it)->phi(), subjets->at(subjetIndices.at(sj)).rapidity(), subjets->at(subjetIndices.at(sj)).
phi() ) );
701 int closestSubjetIdx =
std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
703 assignedParticles.at(closestSubjetIdx).push_back( *it );
const reco::GenParticleRef & particleRef() const
bool isLightParton(const reco::Candidate &c)
boost::shared_ptr< fastjet::JetDefinition > JetDefPtr
virtual void produce(edm::Event &, const edm::EventSetup &)
bool isNonnull() const
Checks for non-null.
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_
void addDefault(ParameterSetDescription const &psetDescription)
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)
JetFlavourClustering(const edm::ParameterSet &)
bool isNull() const
Checks for null.
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_
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
const bool isbHadron() 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_