CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
JetFlavourClustering.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: JetFlavourClustering
4 // Class: JetFlavourClustering
5 //
70 //
71 // Original Author: Dinko Ferencek
72 // Created: Wed Nov 6 00:49:55 CET 2013
73 //
74 //
75 
76 
77 // system include files
78 #include <memory>
79 
80 // user include files
83 
86 
88 
97 
98 #include "fastjet/JetDefinition.hh"
99 #include "fastjet/ClusterSequence.hh"
100 #include "fastjet/Selector.hh"
101 #include "fastjet/PseudoJet.hh"
102 
103 //
104 // constants, enums and typedefs
105 //
106 typedef boost::shared_ptr<fastjet::ClusterSequence> ClusterSequencePtr;
107 typedef boost::shared_ptr<fastjet::JetDefinition> JetDefPtr;
108 
109 //
110 // class declaration
111 //
112 class GhostInfo : public fastjet::PseudoJet::UserInfoBase{
113  public:
114  GhostInfo(const bool & isHadron,
115  const bool & isbHadron,
116  const bool & isParton,
117  const bool & isLepton,
119  m_particleRef(particleRef)
120  {
121  m_type = 0;
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 );
126  }
127 
128  const bool isHadron() const { return ( m_type & (1 << 0) ); }
129  const bool isbHadron() const { return ( m_type & (1 << 1) ); }
130  const bool isParton() const { return ( m_type & (1 << 2) ); }
131  const bool isLepton() const { return ( m_type & (1 << 3) ); }
132  const reco::GenParticleRef & particleRef() const { return m_particleRef; }
133 
134  protected:
136  int m_type;
137 };
138 
140  public:
141  explicit JetFlavourClustering(const edm::ParameterSet&);
143 
144  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
145 
146  private:
147  virtual void beginJob() ;
148  virtual void produce(edm::Event&, const edm::EventSetup&);
149  virtual void endJob() ;
150 
151  virtual void beginRun(edm::Run&, edm::EventSetup const&);
152  virtual void endRun(edm::Run&, edm::EventSetup const&);
155 
157  const double ghostRescaling,
158  const bool isHadron, const bool isbHadron, const bool isParton, const bool isLepton,
159  std::vector<fastjet::PseudoJet>& constituents);
160 
162  const std::vector<fastjet::PseudoJet>& matchedJets,
163  std::vector<int>& matchedIndices);
165  const edm::Handle<edm::View<reco::Jet> >& matchedJets,
166  std::vector<int>& matchedIndices);
167  void matchSubjets(const std::vector<int>& groomedIndices,
168  const edm::Handle<edm::View<reco::Jet> >& groomedJets,
169  const edm::Handle<edm::View<reco::Jet> >& subjets,
170  std::vector<std::vector<int> >& matchedIndices);
171 
172  void setFlavours(const reco::GenParticleRefVector& clusteredbHadrons,
173  const reco::GenParticleRefVector& clusteredcHadrons,
174  const reco::GenParticleRefVector& clusteredPartons,
175  int& hadronFlavour,
176  int& partonFlavour);
177 
178  void assignToSubjets(const reco::GenParticleRefVector& clusteredParticles,
179  const edm::Handle<edm::View<reco::Jet> >& subjets,
180  const std::vector<int>& subjetIndices,
181  std::vector<reco::GenParticleRefVector>& assignedParticles);
182 
183  // ----------member data ---------------------------
184  const edm::EDGetTokenT<edm::View<reco::Jet> > jetsToken_; // Input jet collection
185  const edm::EDGetTokenT<edm::View<reco::Jet> > groomedJetsToken_; // Input groomed jet collection
186  const edm::EDGetTokenT<edm::View<reco::Jet> > subjetsToken_; // Input subjet collection
191 
193  const double rParam_;
194  const double jetPtMin_;
195  const double ghostRescaling_;
197  const bool useSubjets_;
198  const bool useLeptons_;
199 
202 };
203 
204 //
205 // static data member definitions
206 //
207 
208 //
209 // constructors and destructor
210 //
212 
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() )),
216  bHadronsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("bHadrons") )),
217  cHadronsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("cHadrons") )),
218  partonsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("partons") )),
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")),
222  jetPtMin_(0.), // hardcoded to 0. since we simply want to recluster all input jets which already had some PtMin applied
223  ghostRescaling_(iConfig.exists("ghostRescaling") ? iConfig.getParameter<double>("ghostRescaling") : 1e-18),
224  hadronFlavourHasPriority_(iConfig.getParameter<bool>("hadronFlavourHasPriority")),
225  useSubjets_(iConfig.exists("groomedJets") && iConfig.exists("subjets")),
226  useLeptons_(iConfig.exists("leptons"))
227 
228 {
229  // register your products
230  produces<reco::JetFlavourInfoMatchingCollection>();
231  if( useSubjets_ )
232  produces<reco::JetFlavourInfoMatchingCollection>("SubJets");
233 
234  // set jet algorithm
235  if (jetAlgorithm_=="Kt")
236  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::kt_algorithm, rParam_) );
237  else if (jetAlgorithm_=="CambridgeAachen")
238  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::cambridge_algorithm, rParam_) );
239  else if (jetAlgorithm_=="AntiKt")
240  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::antikt_algorithm, rParam_) );
241  else
242  throw cms::Exception("InvalidJetAlgorithm") << "Jet clustering algorithm is invalid: " << jetAlgorithm_ << ", use CambridgeAachen | Kt | AntiKt" << std::endl;
243 }
244 
245 
247 {
248 
249  // do anything here that needs to be done at desctruction time
250  // (e.g. close files, deallocate resources etc.)
251 
252 }
253 
254 
255 //
256 // member functions
257 //
258 
259 // ------------ method called to produce the data ------------
260 void
262 {
264  iEvent.getByToken(jetsToken_, jets);
265 
266  edm::Handle<edm::View<reco::Jet> > groomedJets;
268  if( useSubjets_ )
269  {
270  iEvent.getByToken(groomedJetsToken_, groomedJets);
271  iEvent.getByToken(subjetsToken_, subjets);
272  }
273 
275  iEvent.getByToken(bHadronsToken_, bHadrons);
276 
278  iEvent.getByToken(cHadronsToken_, cHadrons);
279 
281  iEvent.getByToken(partonsToken_, partons);
282 
284  if( useLeptons_ )
285  iEvent.getByToken(leptonsToken_, leptons);
286 
287  std::auto_ptr<reco::JetFlavourInfoMatchingCollection> jetFlavourInfos( new reco::JetFlavourInfoMatchingCollection(reco::JetRefBaseProd(jets)) );
288  std::auto_ptr<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfos;
289  if( useSubjets_ )
290  subjetFlavourInfos = std::auto_ptr<reco::JetFlavourInfoMatchingCollection>( new reco::JetFlavourInfoMatchingCollection(reco::JetRefBaseProd(subjets)) );
291 
292  // vector of constituents for reclustering jets and "ghosts"
293  std::vector<fastjet::PseudoJet> fjInputs;
294  // loop over all input jets and collect all their constituents
295  for(edm::View<reco::Jet>::const_iterator it = jets->begin(); it != jets->end(); ++it)
296  {
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 )
300  {
301  reco::CandidatePtr constit = *m;
302  if(constit->pt() == 0)
303  {
304  edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
305  continue;
306  }
307  fjInputs.push_back(fastjet::PseudoJet(constit->px(),constit->py(),constit->pz(),constit->energy()));
308  }
309  }
310  // insert "ghost" b hadrons in the vector of constituents
311  insertGhosts(bHadrons, ghostRescaling_, true, true, false, false, fjInputs);
312  // insert "ghost" c hadrons in the vector of constituents
313  insertGhosts(cHadrons, ghostRescaling_, true, false, false, false, fjInputs);
314  // insert "ghost" partons in the vector of constituents
315  insertGhosts(partons, ghostRescaling_, false, false, true, false, fjInputs);
316  // if used, insert "ghost" leptons in the vector of constituents
317  if( useLeptons_ )
318  insertGhosts(leptons, ghostRescaling_, false, false, false, true, fjInputs);
319 
320  // define jet clustering sequence
321  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs, *fjJetDefinition_ ) );
322  // recluster jet constituents and inserted "ghosts"
323  std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt( fjClusterSeq_->inclusive_jets(jetPtMin_) );
324 
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.";
327 
328  // match reclustered and original jets
329  std::vector<int> reclusteredIndices;
330  matchReclusteredJets(jets,inclusiveJets,reclusteredIndices);
331 
332  // match groomed and original jets
333  std::vector<int> groomedIndices;
334  if( useSubjets_ )
335  {
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.";
338 
339  matchGroomedJets(jets,groomedJets,groomedIndices);
340  }
341 
342  // match subjets and original jets
343  std::vector<std::vector<int> > subjetIndices;
344  if( useSubjets_ )
345  {
346  matchSubjets(groomedIndices,groomedJets,subjets,subjetIndices);
347  }
348 
349  // determine jet flavour
350  for(size_t i=0; i<jets->size(); ++i)
351  {
352  // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original jets should in principle stay the same
353  if( ( fabs( inclusiveJets.at(reclusteredIndices.at(i)).pt() - jets->at(i).pt() ) / jets->at(i).pt() ) > 1e-3 ) // 0.1% difference in Pt should be sufficient to detect possible misconfigurations
354  {
355  if( jets->at(i).pt() < 10. ) // special handling for low-Pt jets (Pt<10 GeV)
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.";
360  else
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.";
364  }
365 
366  reco::GenParticleRefVector clusteredbHadrons;
367  reco::GenParticleRefVector clusteredcHadrons;
368  reco::GenParticleRefVector clusteredPartons;
369  reco::GenParticleRefVector clusteredLeptons;
370 
371  // get jet constituents (sorted by Pt)
372  std::vector<fastjet::PseudoJet> constituents = fastjet::sorted_by_pt( inclusiveJets.at(reclusteredIndices.at(i)).constituents() );
373 
374  // loop over jet constituents and try to find "ghosts"
375  for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
376  {
377  if( !it->has_user_info() ) continue; // skip if not a "ghost"
378 
379  // "ghost" hadron
380  if( it->user_info<GhostInfo>().isHadron() )
381  {
382  // "ghost" b hadron
383  if( it->user_info<GhostInfo>().isbHadron() )
384  clusteredbHadrons.push_back(it->user_info<GhostInfo>().particleRef());
385  // "ghost" c hadron
386  else
387  clusteredcHadrons.push_back(it->user_info<GhostInfo>().particleRef());
388  }
389  // "ghost" parton
390  else if( it->user_info<GhostInfo>().isParton() )
391  clusteredPartons.push_back(it->user_info<GhostInfo>().particleRef());
392  // "ghost" lepton
393  else if( it->user_info<GhostInfo>().isLepton() )
394  clusteredLeptons.push_back(it->user_info<GhostInfo>().particleRef());
395  }
396 
397  int hadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
398  int partonFlavour = 0; // default parton flavour set to 0 (= undefined)
399 
400  // set hadron- and parton-based flavours
401  setFlavours(clusteredbHadrons, clusteredcHadrons, clusteredPartons, hadronFlavour, partonFlavour);
402 
403  // set the JetFlavourInfo for this jet
404  (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, hadronFlavour, partonFlavour);
405 
406  // if subjets are used, determine their flavour
407  if( useSubjets_ )
408  {
409  if( subjetIndices.at(i).size()==0 ) continue; // continue if the original jet does not have subjets assigned
410 
411  // define vectors of GenParticleRefVectors for hadrons and partons assigned to different subjets
412  std::vector<reco::GenParticleRefVector> assignedbHadrons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
413  std::vector<reco::GenParticleRefVector> assignedcHadrons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
414  std::vector<reco::GenParticleRefVector> assignedPartons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
415  std::vector<reco::GenParticleRefVector> assignedLeptons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
416 
417  // loop over clustered b hadrons and assign them to different subjets based on smallest dR
418  assignToSubjets(clusteredbHadrons, subjets, subjetIndices.at(i), assignedbHadrons);
419  // loop over clustered c hadrons and assign them to different subjets based on smallest dR
420  assignToSubjets(clusteredcHadrons, subjets, subjetIndices.at(i), assignedcHadrons);
421  // loop over clustered partons and assign them to different subjets based on smallest dR
422  assignToSubjets(clusteredPartons, subjets, subjetIndices.at(i), assignedPartons);
423  // if used, loop over clustered leptons and assign them to different subjets based on smallest dR
424  if( useLeptons_ )
425  assignToSubjets(clusteredLeptons, subjets, subjetIndices.at(i), assignedLeptons);
426 
427  // loop over subjets and determine their flavour
428  for(size_t sj=0; sj<subjetIndices.at(i).size(); ++sj)
429  {
430  int subjetHadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
431  int subjetPartonFlavour = 0; // default parton flavour set to 0 (= undefined)
432 
433  // set hadron- and parton-based flavours
434  setFlavours(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
435 
436  // set the JetFlavourInfo for this subjet
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);
438  }
439  }
440  }
441 
442  // put jet flavour infos in the event
443  iEvent.put( jetFlavourInfos );
444  // put subjet flavour infos in the event
445  if( useSubjets_ )
446  iEvent.put( subjetFlavourInfos, "SubJets" );
447 }
448 
449 // ------------ method that inserts "ghost" particles in the vector of jet constituents ------------
450 void
452  const double ghostRescaling,
453  const bool isHadron, const bool isbHadron, const bool isParton, const bool isLepton,
454  std::vector<fastjet::PseudoJet>& constituents)
455 {
456  // insert "ghost" particles in the vector of jet constituents
457  for(reco::GenParticleRefVector::const_iterator it = particles->begin(); it != particles->end(); ++it)
458  {
459  fastjet::PseudoJet p((*it)->px(),(*it)->py(),(*it)->pz(),(*it)->energy());
460  p*=ghostRescaling; // rescale particle momentum
461  p.set_user_info(new GhostInfo(isHadron, isbHadron, isParton, isLepton, *it));
462  constituents.push_back(p);
463  }
464 }
465 
466 // ------------ method that matches reclustered and original jets based on minimum dR ------------
467 void
469  const std::vector<fastjet::PseudoJet>& reclusteredJets,
470  std::vector<int>& matchedIndices)
471 {
472  std::vector<bool> matchedLocks(reclusteredJets.size(),false);
473 
474  for(size_t j=0; j<jets->size(); ++j)
475  {
476  double matchedDR = 1e9;
477  int matchedIdx = -1;
478 
479  for(size_t rj=0; rj<reclusteredJets.size(); ++rj)
480  {
481  if( matchedLocks.at(rj) ) continue; // skip jets that have already been matched
482 
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 )
485  {
486  matchedDR = tempDR;
487  matchedIdx = rj;
488  }
489  }
490 
491  if( matchedIdx>=0 ) matchedLocks.at(matchedIdx) = true;
492  matchedIndices.push_back(matchedIdx);
493  }
494 
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.";
497 }
498 
499 // ------------ method that matches groomed and original jets based on minimum dR ------------
500 void
502  const edm::Handle<edm::View<reco::Jet> >& groomedJets,
503  std::vector<int>& matchedIndices)
504 {
505  std::vector<bool> jetLocks(jets->size(),false);
506  std::vector<int> jetIndices;
507 
508  for(size_t gj=0; gj<groomedJets->size(); ++gj)
509  {
510  double matchedDR = 1e9;
511  int matchedIdx = -1;
512 
513  for(size_t j=0; j<jets->size(); ++j)
514  {
515  if( jetLocks.at(j) ) continue; // skip jets that have already been matched
516 
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 )
519  {
520  matchedDR = tempDR;
521  matchedIdx = j;
522  }
523  }
524 
525  if( matchedIdx>=0 ) jetLocks.at(matchedIdx) = true;
526  jetIndices.push_back(matchedIdx);
527  }
528 
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.";
531 
532  for(size_t j=0; j<jets->size(); ++j)
533  {
534  std::vector<int>::iterator matchedIndex = std::find( jetIndices.begin(), jetIndices.end(), j );
535 
536  matchedIndices.push_back( matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(),matchedIndex) : -1 );
537  }
538 }
539 
540 // ------------ method that matches subjets and original jets ------------
541 void
542 JetFlavourClustering::matchSubjets(const std::vector<int>& groomedIndices,
543  const edm::Handle<edm::View<reco::Jet> >& groomedJets,
544  const edm::Handle<edm::View<reco::Jet> >& subjets,
545  std::vector<std::vector<int> >& matchedIndices)
546 {
547  for(size_t g=0; g<groomedIndices.size(); ++g)
548  {
549  std::vector<int> subjetIndices;
550 
551  if( groomedIndices.at(g)>=0 )
552  {
553  for(size_t s=0; s<groomedJets->at(groomedIndices.at(g)).numberOfDaughters(); ++s)
554  {
555  const edm::Ptr<reco::Candidate> & subjet = groomedJets->at(groomedIndices.at(g)).daughterPtr(s);
556 
557  for(size_t sj=0; sj<subjets->size(); ++sj)
558  {
559  if( subjet == edm::Ptr<reco::Candidate>(subjets->ptrAt(sj)) )
560  {
561  subjetIndices.push_back(sj);
562  break;
563  }
564  }
565  }
566 
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.";
569 
570  matchedIndices.push_back(subjetIndices);
571  }
572  else
573  matchedIndices.push_back(subjetIndices);
574  }
575 }
576 
577 // ------------ method that sets hadron- and parton-based flavours ------------
578 void
580  const reco::GenParticleRefVector& clusteredcHadrons,
581  const reco::GenParticleRefVector& clusteredPartons,
582  int& hadronFlavour,
583  int& partonFlavour)
584 {
585  reco::GenParticleRef hardestParton;
586  reco::GenParticleRef hardestLightParton;
587  reco::GenParticleRef flavourParton;
588 
589  // loop over clustered partons (already sorted by Pt)
590  for(reco::GenParticleRefVector::const_iterator it = clusteredPartons.begin(); it != clusteredPartons.end(); ++it)
591  {
592  // hardest parton
593  if( hardestParton.isNull() )
594  hardestParton = (*it);
595  // hardest light-flavour parton
596  if( hardestLightParton.isNull() )
597  {
598  if( CandMCTagUtils::isLightParton( *(*it) ) )
599  hardestLightParton = (*it);
600  }
601  // c flavour
602  if( flavourParton.isNull() && ( abs( (*it)->pdgId() ) == 4 ) )
603  flavourParton = (*it);
604  // b flavour gets priority
605  if( abs( (*it)->pdgId() ) == 5 )
606  {
607  if( flavourParton.isNull() )
608  flavourParton = (*it);
609  else if( abs( flavourParton->pdgId() ) != 5 )
610  flavourParton = (*it);
611  }
612  }
613 
614  // set hadron-based flavour
615  if( clusteredbHadrons.size()>0 )
616  hadronFlavour = 5;
617  else if( clusteredcHadrons.size()>0 && clusteredbHadrons.size()==0 )
618  hadronFlavour = 4;
619  // set parton-based flavour
620  if( flavourParton.isNull() )
621  {
622  if( hardestParton.isNonnull() ) partonFlavour = hardestParton->pdgId();
623  }
624  else
625  partonFlavour = flavourParton->pdgId();
626 
627  // if enabled, check for conflicts between hadron- and parton-based flavours and give priority to the hadron-based flavour
629  {
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;
634  }
635 }
636 
637 // ------------ method that assigns clustered particles to subjets ------------
638 void
640  const edm::Handle<edm::View<reco::Jet> >& subjets,
641  const std::vector<int>& subjetIndices,
642  std::vector<reco::GenParticleRefVector>& assignedParticles)
643 {
644  // loop over clustered particles and assign them to different subjets based on smallest dR
645  for(reco::GenParticleRefVector::const_iterator it = clusteredParticles.begin(); it != clusteredParticles.end(); ++it)
646  {
647  std::vector<double> dR2toSubjets;
648 
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() ) );
651 
652  // find the closest subjet
653  int closestSubjetIdx = std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
654 
655  assignedParticles.at(closestSubjetIdx).push_back( *it );
656  }
657 }
658 
659 // ------------ method called once each job just before starting event loop ------------
660 void
662 {
663 }
664 
665 // ------------ method called once each job just after ending the event loop ------------
666 void
668 }
669 
670 // ------------ method called when starting to processes a run ------------
671 void
673 {
674 }
675 
676 // ------------ method called when ending the processing of a run ------------
677 void
679 {
680 }
681 
682 // ------------ method called when starting to processes a luminosity block ------------
683 void
685 {
686 }
687 
688 // ------------ method called when ending the processing of a luminosity block ------------
689 void
691 {
692 }
693 
694 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
695 void
697  //The following says we do not know what parameters are allowed so do no validation
698  // Please change this to state exactly what you do use, even if it is no parameters
700  desc.setUnknown();
701  descriptions.addDefault(desc);
702 }
703 
704 //define this as a plug-in
const reco::GenParticleRef & particleRef() const
bool isLightParton(const reco::Candidate &c)
Definition: CandMCTag.cc:63
boost::shared_ptr< fastjet::JetDefinition > JetDefPtr
int i
Definition: DBlmapReader.cc:9
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:81
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.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
boost::shared_ptr< fastjet::ClusterSequence > ClusterSequencePtr
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
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)
Definition: pdgIdUtils.h:19
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.
Definition: RefVector.h:249
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
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
Definition: Activities.doc:4
const bool isParton() const
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:244
const bool isLepton() const
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const edm::EDGetTokenT< reco::GenParticleRefVector > cHadronsToken_
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
int iEvent
Definition: GenABIO.cc:243
bool isNull() const
Checks for null.
Definition: Ref.h:247
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.
Definition: Event.h:116
vector< PseudoJet > jets
Class storing the jet flavour information.
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:30
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)
Definition: CandMCTag.cc:48
T1 deltaR2(T1 eta1, T2 phi1, T3 eta2, T4 phi2)
Definition: deltaR.h:58
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&lt;C, T&gt; to the RefVector.
Definition: RefVector.h:64
size_type size() const
Size of the RefVector.
Definition: RefVector.h:89
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_
Definition: Run.h:41
Definition: DDAxes.h:10