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_;
196  const double relPtTolerance_;
198  const bool useSubjets_;
199  const bool useLeptons_;
200 
203 };
204 
205 //
206 // static data member definitions
207 //
208 
209 //
210 // constructors and destructor
211 //
213 
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() )),
217  bHadronsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("bHadrons") )),
218  cHadronsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("cHadrons") )),
219  partonsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("partons") )),
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")),
223  jetPtMin_(0.), // hardcoded to 0. since we simply want to recluster all input jets which already had some PtMin applied
224  ghostRescaling_(iConfig.exists("ghostRescaling") ? iConfig.getParameter<double>("ghostRescaling") : 1e-18),
225  relPtTolerance_(iConfig.exists("relPtTolerance") ? iConfig.getParameter<double>("relPtTolerance") : 1e-03), // 0.1% relative difference in Pt should be sufficient to detect possible misconfigurations
226  hadronFlavourHasPriority_(iConfig.getParameter<bool>("hadronFlavourHasPriority")),
227  useSubjets_(iConfig.exists("groomedJets") && iConfig.exists("subjets")),
228  useLeptons_(iConfig.exists("leptons"))
229 
230 {
231  // register your products
232  produces<reco::JetFlavourInfoMatchingCollection>();
233  if( useSubjets_ )
234  produces<reco::JetFlavourInfoMatchingCollection>("SubJets");
235 
236  // set jet algorithm
237  if (jetAlgorithm_=="Kt")
238  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::kt_algorithm, rParam_) );
239  else if (jetAlgorithm_=="CambridgeAachen")
240  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::cambridge_algorithm, rParam_) );
241  else if (jetAlgorithm_=="AntiKt")
242  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::antikt_algorithm, rParam_) );
243  else
244  throw cms::Exception("InvalidJetAlgorithm") << "Jet clustering algorithm is invalid: " << jetAlgorithm_ << ", use CambridgeAachen | Kt | AntiKt" << std::endl;
245 }
246 
247 
249 {
250 
251  // do anything here that needs to be done at desctruction time
252  // (e.g. close files, deallocate resources etc.)
253 
254 }
255 
256 
257 //
258 // member functions
259 //
260 
261 // ------------ method called to produce the data ------------
262 void
264 {
266  iEvent.getByToken(jetsToken_, jets);
267 
268  edm::Handle<edm::View<reco::Jet> > groomedJets;
270  if( useSubjets_ )
271  {
272  iEvent.getByToken(groomedJetsToken_, groomedJets);
273  iEvent.getByToken(subjetsToken_, subjets);
274  }
275 
277  iEvent.getByToken(bHadronsToken_, bHadrons);
278 
280  iEvent.getByToken(cHadronsToken_, cHadrons);
281 
283  iEvent.getByToken(partonsToken_, partons);
284 
286  if( useLeptons_ )
287  iEvent.getByToken(leptonsToken_, leptons);
288 
289  std::auto_ptr<reco::JetFlavourInfoMatchingCollection> jetFlavourInfos( new reco::JetFlavourInfoMatchingCollection(reco::JetRefBaseProd(jets)) );
290  std::auto_ptr<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfos;
291  if( useSubjets_ )
292  subjetFlavourInfos = std::auto_ptr<reco::JetFlavourInfoMatchingCollection>( new reco::JetFlavourInfoMatchingCollection(reco::JetRefBaseProd(subjets)) );
293 
294  // vector of constituents for reclustering jets and "ghosts"
295  std::vector<fastjet::PseudoJet> fjInputs;
296  // loop over all input jets and collect all their constituents
297  for(edm::View<reco::Jet>::const_iterator it = jets->begin(); it != jets->end(); ++it)
298  {
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 )
302  {
303  reco::CandidatePtr constit = *m;
304  if(constit->pt() == 0)
305  {
306  edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
307  continue;
308  }
309  fjInputs.push_back(fastjet::PseudoJet(constit->px(),constit->py(),constit->pz(),constit->energy()));
310  }
311  }
312  // insert "ghost" b hadrons in the vector of constituents
313  insertGhosts(bHadrons, ghostRescaling_, true, true, false, false, fjInputs);
314  // insert "ghost" c hadrons in the vector of constituents
315  insertGhosts(cHadrons, ghostRescaling_, true, false, false, false, fjInputs);
316  // insert "ghost" partons in the vector of constituents
317  insertGhosts(partons, ghostRescaling_, false, false, true, false, fjInputs);
318  // if used, insert "ghost" leptons in the vector of constituents
319  if( useLeptons_ )
320  insertGhosts(leptons, ghostRescaling_, false, false, false, true, fjInputs);
321 
322  // define jet clustering sequence
323  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs, *fjJetDefinition_ ) );
324  // recluster jet constituents and inserted "ghosts"
325  std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt( fjClusterSeq_->inclusive_jets(jetPtMin_) );
326 
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.";
329 
330  // match reclustered and original jets
331  std::vector<int> reclusteredIndices;
332  matchReclusteredJets(jets,inclusiveJets,reclusteredIndices);
333 
334  // match groomed and original jets
335  std::vector<int> groomedIndices;
336  if( useSubjets_ )
337  {
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.";
340 
341  matchGroomedJets(jets,groomedJets,groomedIndices);
342  }
343 
344  // match subjets and original jets
345  std::vector<std::vector<int> > subjetIndices;
346  if( useSubjets_ )
347  {
348  matchSubjets(groomedIndices,groomedJets,subjets,subjetIndices);
349  }
350 
351  // determine jet flavour
352  for(size_t i=0; i<jets->size(); ++i)
353  {
354  reco::GenParticleRefVector clusteredbHadrons;
355  reco::GenParticleRefVector clusteredcHadrons;
356  reco::GenParticleRefVector clusteredPartons;
357  reco::GenParticleRefVector clusteredLeptons;
358 
359  // if matching reclustered to original jets failed
360  if( reclusteredIndices.at(i) < 0 )
361  {
362  // set an empty JetFlavourInfo for this jet
363  (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
364  }
365  else
366  {
367  // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original jets should in principle stay the same
368  if( ( fabs( inclusiveJets.at(reclusteredIndices.at(i)).pt() - jets->at(i).pt() ) / jets->at(i).pt() ) > relPtTolerance_ )
369  {
370  if( jets->at(i).pt() < 10. ) // special handling for low-Pt jets (Pt<10 GeV)
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.";
375  else
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.";
379  }
380 
381  // get jet constituents (sorted by Pt)
382  std::vector<fastjet::PseudoJet> constituents = fastjet::sorted_by_pt( inclusiveJets.at(reclusteredIndices.at(i)).constituents() );
383 
384  // loop over jet constituents and try to find "ghosts"
385  for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
386  {
387  if( !it->has_user_info() ) continue; // skip if not a "ghost"
388 
389  // "ghost" hadron
390  if( it->user_info<GhostInfo>().isHadron() )
391  {
392  // "ghost" b hadron
393  if( it->user_info<GhostInfo>().isbHadron() )
394  clusteredbHadrons.push_back(it->user_info<GhostInfo>().particleRef());
395  // "ghost" c hadron
396  else
397  clusteredcHadrons.push_back(it->user_info<GhostInfo>().particleRef());
398  }
399  // "ghost" parton
400  else if( it->user_info<GhostInfo>().isParton() )
401  clusteredPartons.push_back(it->user_info<GhostInfo>().particleRef());
402  // "ghost" lepton
403  else if( it->user_info<GhostInfo>().isLepton() )
404  clusteredLeptons.push_back(it->user_info<GhostInfo>().particleRef());
405  }
406 
407  int hadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
408  int partonFlavour = 0; // default parton flavour set to 0 (= undefined)
409 
410  // set hadron- and parton-based flavours
411  setFlavours(clusteredbHadrons, clusteredcHadrons, clusteredPartons, hadronFlavour, partonFlavour);
412 
413  // set the JetFlavourInfo for this jet
414  (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, hadronFlavour, partonFlavour);
415  }
416 
417  // if subjets are used, determine their flavour
418  if( useSubjets_ )
419  {
420  if( subjetIndices.at(i).size()==0 ) continue; // continue if the original jet does not have subjets assigned
421 
422  // define vectors of GenParticleRefVectors for hadrons and partons assigned to different subjets
423  std::vector<reco::GenParticleRefVector> assignedbHadrons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
424  std::vector<reco::GenParticleRefVector> assignedcHadrons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
425  std::vector<reco::GenParticleRefVector> assignedPartons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
426  std::vector<reco::GenParticleRefVector> assignedLeptons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
427 
428  // loop over clustered b hadrons and assign them to different subjets based on smallest dR
429  assignToSubjets(clusteredbHadrons, subjets, subjetIndices.at(i), assignedbHadrons);
430  // loop over clustered c hadrons and assign them to different subjets based on smallest dR
431  assignToSubjets(clusteredcHadrons, subjets, subjetIndices.at(i), assignedcHadrons);
432  // loop over clustered partons and assign them to different subjets based on smallest dR
433  assignToSubjets(clusteredPartons, subjets, subjetIndices.at(i), assignedPartons);
434  // if used, loop over clustered leptons and assign them to different subjets based on smallest dR
435  if( useLeptons_ )
436  assignToSubjets(clusteredLeptons, subjets, subjetIndices.at(i), assignedLeptons);
437 
438  // loop over subjets and determine their flavour
439  for(size_t sj=0; sj<subjetIndices.at(i).size(); ++sj)
440  {
441  int subjetHadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
442  int subjetPartonFlavour = 0; // default parton flavour set to 0 (= undefined)
443 
444  // set hadron- and parton-based flavours
445  setFlavours(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
446 
447  // set the JetFlavourInfo for this subjet
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);
449  }
450  }
451  }
452 
453  // put jet flavour infos in the event
454  iEvent.put( jetFlavourInfos );
455  // put subjet flavour infos in the event
456  if( useSubjets_ )
457  iEvent.put( subjetFlavourInfos, "SubJets" );
458 }
459 
460 // ------------ method that inserts "ghost" particles in the vector of jet constituents ------------
461 void
463  const double ghostRescaling,
464  const bool isHadron, const bool isbHadron, const bool isParton, const bool isLepton,
465  std::vector<fastjet::PseudoJet>& constituents)
466 {
467  // insert "ghost" particles in the vector of jet constituents
468  for(reco::GenParticleRefVector::const_iterator it = particles->begin(); it != particles->end(); ++it)
469  {
470  if((*it)->pt() == 0)
471  {
472  edm::LogInfo("NullTransverseMomentum") << "dropping input ghost candidate with pt=0";
473  continue;
474  }
475  fastjet::PseudoJet p((*it)->px(),(*it)->py(),(*it)->pz(),(*it)->energy());
476  p*=ghostRescaling; // rescale particle momentum
477  p.set_user_info(new GhostInfo(isHadron, isbHadron, isParton, isLepton, *it));
478  constituents.push_back(p);
479  }
480 }
481 
482 // ------------ method that matches reclustered and original jets based on minimum dR ------------
483 void
485  const std::vector<fastjet::PseudoJet>& reclusteredJets,
486  std::vector<int>& matchedIndices)
487 {
488  std::vector<bool> matchedLocks(reclusteredJets.size(),false);
489 
490  for(size_t j=0; j<jets->size(); ++j)
491  {
492  double matchedDR2 = 1e9;
493  int matchedIdx = -1;
494 
495  for(size_t rj=0; rj<reclusteredJets.size(); ++rj)
496  {
497  if( matchedLocks.at(rj) ) continue; // skip jets that have already been matched
498 
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 )
501  {
502  matchedDR2 = tempDR2;
503  matchedIdx = rj;
504  }
505  }
506 
507  if( matchedIdx>=0 )
508  {
509  if ( matchedDR2 > rParam_*rParam_ )
510  {
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.";
513  }
514  else
515  matchedLocks.at(matchedIdx) = true;
516  }
517  else
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.";
519 
520  matchedIndices.push_back(matchedIdx);
521  }
522 }
523 
524 // ------------ method that matches groomed and original jets based on minimum dR ------------
525 void
527  const edm::Handle<edm::View<reco::Jet> >& groomedJets,
528  std::vector<int>& matchedIndices)
529 {
530  std::vector<bool> jetLocks(jets->size(),false);
531  std::vector<int> jetIndices;
532 
533  for(size_t gj=0; gj<groomedJets->size(); ++gj)
534  {
535  double matchedDR2 = 1e9;
536  int matchedIdx = -1;
537 
538  if( groomedJets->at(gj).pt()>0. ) // skip pathological cases of groomed jets with Pt=0
539  {
540  for(size_t j=0; j<jets->size(); ++j)
541  {
542  if( jetLocks.at(j) ) continue; // skip jets that have already been matched
543 
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 )
546  {
547  matchedDR2 = tempDR2;
548  matchedIdx = j;
549  }
550  }
551  }
552 
553  if( matchedIdx>=0 )
554  {
555  if ( matchedDR2 > rParam_*rParam_ )
556  {
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.";
559  matchedIdx = -1;
560  }
561  else
562  jetLocks.at(matchedIdx) = true;
563  }
564  jetIndices.push_back(matchedIdx);
565  }
566 
567  for(size_t j=0; j<jets->size(); ++j)
568  {
569  std::vector<int>::iterator matchedIndex = std::find( jetIndices.begin(), jetIndices.end(), j );
570 
571  matchedIndices.push_back( matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(),matchedIndex) : -1 );
572  }
573 }
574 
575 // ------------ method that matches subjets and original jets ------------
576 void
577 JetFlavourClustering::matchSubjets(const std::vector<int>& groomedIndices,
578  const edm::Handle<edm::View<reco::Jet> >& groomedJets,
579  const edm::Handle<edm::View<reco::Jet> >& subjets,
580  std::vector<std::vector<int> >& matchedIndices)
581 {
582  for(size_t g=0; g<groomedIndices.size(); ++g)
583  {
584  std::vector<int> subjetIndices;
585 
586  if( groomedIndices.at(g)>=0 )
587  {
588  for(size_t s=0; s<groomedJets->at(groomedIndices.at(g)).numberOfDaughters(); ++s)
589  {
590  const edm::Ptr<reco::Candidate> & subjet = groomedJets->at(groomedIndices.at(g)).daughterPtr(s);
591 
592  for(size_t sj=0; sj<subjets->size(); ++sj)
593  {
594  if( subjet == edm::Ptr<reco::Candidate>(subjets->ptrAt(sj)) )
595  {
596  subjetIndices.push_back(sj);
597  break;
598  }
599  }
600  }
601 
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.";
604 
605  matchedIndices.push_back(subjetIndices);
606  }
607  else
608  matchedIndices.push_back(subjetIndices);
609  }
610 }
611 
612 // ------------ method that sets hadron- and parton-based flavours ------------
613 void
615  const reco::GenParticleRefVector& clusteredcHadrons,
616  const reco::GenParticleRefVector& clusteredPartons,
617  int& hadronFlavour,
618  int& partonFlavour)
619 {
620  reco::GenParticleRef hardestParton;
621  reco::GenParticleRef hardestLightParton;
622  reco::GenParticleRef flavourParton;
623 
624  // loop over clustered partons (already sorted by Pt)
625  for(reco::GenParticleRefVector::const_iterator it = clusteredPartons.begin(); it != clusteredPartons.end(); ++it)
626  {
627  // hardest parton
628  if( hardestParton.isNull() )
629  hardestParton = (*it);
630  // hardest light-flavour parton
631  if( hardestLightParton.isNull() )
632  {
633  if( CandMCTagUtils::isLightParton( *(*it) ) )
634  hardestLightParton = (*it);
635  }
636  // c flavour
637  if( flavourParton.isNull() && ( abs( (*it)->pdgId() ) == 4 ) )
638  flavourParton = (*it);
639  // b flavour gets priority
640  if( abs( (*it)->pdgId() ) == 5 )
641  {
642  if( flavourParton.isNull() )
643  flavourParton = (*it);
644  else if( abs( flavourParton->pdgId() ) != 5 )
645  flavourParton = (*it);
646  }
647  }
648 
649  // set hadron-based flavour
650  if( clusteredbHadrons.size()>0 )
651  hadronFlavour = 5;
652  else if( clusteredcHadrons.size()>0 && clusteredbHadrons.size()==0 )
653  hadronFlavour = 4;
654  // set parton-based flavour
655  if( flavourParton.isNull() )
656  {
657  if( hardestParton.isNonnull() ) partonFlavour = hardestParton->pdgId();
658  }
659  else
660  partonFlavour = flavourParton->pdgId();
661 
662  // if enabled, check for conflicts between hadron- and parton-based flavours and give priority to the hadron-based flavour
664  {
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;
669  }
670 }
671 
672 // ------------ method that assigns clustered particles to subjets ------------
673 void
675  const edm::Handle<edm::View<reco::Jet> >& subjets,
676  const std::vector<int>& subjetIndices,
677  std::vector<reco::GenParticleRefVector>& assignedParticles)
678 {
679  // loop over clustered particles and assign them to different subjets based on smallest dR
680  for(reco::GenParticleRefVector::const_iterator it = clusteredParticles.begin(); it != clusteredParticles.end(); ++it)
681  {
682  std::vector<double> dR2toSubjets;
683 
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() ) );
686 
687  // find the closest subjet
688  int closestSubjetIdx = std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
689 
690  assignedParticles.at(closestSubjetIdx).push_back( *it );
691  }
692 }
693 
694 // ------------ method called once each job just before starting event loop ------------
695 void
697 {
698 }
699 
700 // ------------ method called once each job just after ending the event loop ------------
701 void
703 }
704 
705 // ------------ method called when starting to processes a run ------------
706 void
708 {
709 }
710 
711 // ------------ method called when ending the processing of a run ------------
712 void
714 {
715 }
716 
717 // ------------ method called when starting to processes a luminosity block ------------
718 void
720 {
721 }
722 
723 // ------------ method called when ending the processing of a luminosity block ------------
724 void
726 {
727 }
728 
729 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
730 void
732  //The following says we do not know what parameters are allowed so do no validation
733  // Please change this to state exactly what you do use, even if it is no parameters
735  desc.setUnknown();
736  descriptions.addDefault(desc);
737 }
738 
739 //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:230
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
T sqrt(T t)
Definition: SSEVec.h:48
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
double deltaR2(const T1 &t1, const T2 &t2)
Definition: deltaR.h:36
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
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