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 //
72 //
73 // Original Author: Dinko Ferencek
74 // Created: Wed Nov 6 00:49:55 CET 2013
75 //
76 //
77 
78 
79 // system include files
80 #include <memory>
81 
82 // user include files
85 
88 
90 
99 
100 #include "fastjet/JetDefinition.hh"
101 #include "fastjet/ClusterSequence.hh"
102 #include "fastjet/Selector.hh"
103 #include "fastjet/PseudoJet.hh"
104 
105 //
106 // constants, enums and typedefs
107 //
108 typedef boost::shared_ptr<fastjet::ClusterSequence> ClusterSequencePtr;
109 typedef boost::shared_ptr<fastjet::JetDefinition> JetDefPtr;
110 
111 //
112 // class declaration
113 //
114 class GhostInfo : public fastjet::PseudoJet::UserInfoBase{
115  public:
116  GhostInfo(const bool & isHadron,
117  const bool & isbHadron,
118  const bool & isParton,
119  const bool & isLepton,
121  m_particleRef(particleRef)
122  {
123  m_type = 0;
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 );
128  }
129 
130  const bool isHadron() const { return ( m_type & (1 << 0) ); }
131  const bool isbHadron() const { return ( m_type & (1 << 1) ); }
132  const bool isParton() const { return ( m_type & (1 << 2) ); }
133  const bool isLepton() const { return ( m_type & (1 << 3) ); }
134  const reco::GenParticleRef & particleRef() const { return m_particleRef; }
135 
136  protected:
138  int m_type;
139 };
140 
142  public:
143  explicit JetFlavourClustering(const edm::ParameterSet&);
145 
146  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
147 
148  private:
149  virtual void produce(edm::Event&, const edm::EventSetup&);
150 
152  const double ghostRescaling,
153  const bool isHadron, const bool isbHadron, const bool isParton, const bool isLepton,
154  std::vector<fastjet::PseudoJet>& constituents);
155 
157  const std::vector<fastjet::PseudoJet>& matchedJets,
158  std::vector<int>& matchedIndices);
160  const edm::Handle<edm::View<reco::Jet> >& matchedJets,
161  std::vector<int>& matchedIndices);
162  void matchSubjets(const std::vector<int>& groomedIndices,
163  const edm::Handle<edm::View<reco::Jet> >& groomedJets,
164  const edm::Handle<edm::View<reco::Jet> >& subjets,
165  std::vector<std::vector<int> >& matchedIndices);
166 
167  void setFlavours(const reco::GenParticleRefVector& clusteredbHadrons,
168  const reco::GenParticleRefVector& clusteredcHadrons,
169  const reco::GenParticleRefVector& clusteredPartons,
170  int& hadronFlavour,
171  int& partonFlavour);
172 
173  void assignToSubjets(const reco::GenParticleRefVector& clusteredParticles,
174  const edm::Handle<edm::View<reco::Jet> >& subjets,
175  const std::vector<int>& subjetIndices,
176  std::vector<reco::GenParticleRefVector>& assignedParticles);
177 
178  // ----------member data ---------------------------
179  const edm::EDGetTokenT<edm::View<reco::Jet> > jetsToken_; // Input jet collection
180  const edm::EDGetTokenT<edm::View<reco::Jet> > groomedJetsToken_; // Input groomed jet collection
181  const edm::EDGetTokenT<edm::View<reco::Jet> > subjetsToken_; // Input subjet collection
186 
188  const double rParam_;
189  const double jetPtMin_;
190  const double ghostRescaling_;
191  const double relPtTolerance_;
193  const bool useSubjets_;
194  const bool useLeptons_;
195 
198 };
199 
200 //
201 // static data member definitions
202 //
203 
204 //
205 // constructors and destructor
206 //
208 
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() )),
212  bHadronsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("bHadrons") )),
213  cHadronsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("cHadrons") )),
214  partonsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("partons") )),
215  leptonsToken_(mayConsume<reco::GenParticleRefVector>( iConfig.exists("leptons") ? iConfig.getParameter<edm::InputTag>("leptons") : edm::InputTag() )),
216  jetAlgorithm_(iConfig.getParameter<std::string>("jetAlgorithm")),
217  rParam_(iConfig.getParameter<double>("rParam")),
218  jetPtMin_(0.), // hardcoded to 0. since we simply want to recluster all input jets which already had some PtMin applied
219  ghostRescaling_(iConfig.exists("ghostRescaling") ? iConfig.getParameter<double>("ghostRescaling") : 1e-18),
220  relPtTolerance_(iConfig.exists("relPtTolerance") ? iConfig.getParameter<double>("relPtTolerance") : 1e-03), // 0.1% relative difference in Pt should be sufficient to detect possible misconfigurations
221  hadronFlavourHasPriority_(iConfig.getParameter<bool>("hadronFlavourHasPriority")),
222  useSubjets_(iConfig.exists("groomedJets") && iConfig.exists("subjets")),
223  useLeptons_(iConfig.exists("leptons"))
224 
225 {
226  // register your products
227  produces<reco::JetFlavourInfoMatchingCollection>();
228  if( useSubjets_ )
229  produces<reco::JetFlavourInfoMatchingCollection>("SubJets");
230 
231  // set jet algorithm
232  if (jetAlgorithm_=="Kt")
233  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::kt_algorithm, rParam_) );
234  else if (jetAlgorithm_=="CambridgeAachen")
235  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::cambridge_algorithm, rParam_) );
236  else if (jetAlgorithm_=="AntiKt")
237  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::antikt_algorithm, rParam_) );
238  else
239  throw cms::Exception("InvalidJetAlgorithm") << "Jet clustering algorithm is invalid: " << jetAlgorithm_ << ", use CambridgeAachen | Kt | AntiKt" << std::endl;
240 }
241 
242 
244 {
245 
246  // do anything here that needs to be done at desctruction time
247  // (e.g. close files, deallocate resources etc.)
248 
249 }
250 
251 
252 //
253 // member functions
254 //
255 
256 // ------------ method called to produce the data ------------
257 void
259 {
261  iEvent.getByToken(jetsToken_, jets);
262 
263  edm::Handle<edm::View<reco::Jet> > groomedJets;
265  if( useSubjets_ )
266  {
267  iEvent.getByToken(groomedJetsToken_, groomedJets);
268  iEvent.getByToken(subjetsToken_, subjets);
269  }
270 
272  iEvent.getByToken(bHadronsToken_, bHadrons);
273 
275  iEvent.getByToken(cHadronsToken_, cHadrons);
276 
278  iEvent.getByToken(partonsToken_, partons);
279 
281  if( useLeptons_ )
282  iEvent.getByToken(leptonsToken_, leptons);
283 
284  std::auto_ptr<reco::JetFlavourInfoMatchingCollection> jetFlavourInfos( new reco::JetFlavourInfoMatchingCollection(reco::JetRefBaseProd(jets)) );
285  std::auto_ptr<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfos;
286  if( useSubjets_ )
287  subjetFlavourInfos = std::auto_ptr<reco::JetFlavourInfoMatchingCollection>( new reco::JetFlavourInfoMatchingCollection(reco::JetRefBaseProd(subjets)) );
288 
289  // vector of constituents for reclustering jets and "ghosts"
290  std::vector<fastjet::PseudoJet> fjInputs;
291  // loop over all input jets and collect all their constituents
292  for(edm::View<reco::Jet>::const_iterator it = jets->begin(); it != jets->end(); ++it)
293  {
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 )
297  {
298  reco::CandidatePtr constit = *m;
299  if(constit->pt() == 0)
300  {
301  edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
302  continue;
303  }
304  fjInputs.push_back(fastjet::PseudoJet(constit->px(),constit->py(),constit->pz(),constit->energy()));
305  }
306  }
307  // insert "ghost" b hadrons in the vector of constituents
308  insertGhosts(bHadrons, ghostRescaling_, true, true, false, false, fjInputs);
309  // insert "ghost" c hadrons in the vector of constituents
310  insertGhosts(cHadrons, ghostRescaling_, true, false, false, false, fjInputs);
311  // insert "ghost" partons in the vector of constituents
312  insertGhosts(partons, ghostRescaling_, false, false, true, false, fjInputs);
313  // if used, insert "ghost" leptons in the vector of constituents
314  if( useLeptons_ )
315  insertGhosts(leptons, ghostRescaling_, false, false, false, true, fjInputs);
316 
317  // define jet clustering sequence
318  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs, *fjJetDefinition_ ) );
319  // recluster jet constituents and inserted "ghosts"
320  std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt( fjClusterSeq_->inclusive_jets(jetPtMin_) );
321 
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.";
324 
325  // match reclustered and original jets
326  std::vector<int> reclusteredIndices;
327  matchReclusteredJets(jets,inclusiveJets,reclusteredIndices);
328 
329  // match groomed and original jets
330  std::vector<int> groomedIndices;
331  if( useSubjets_ )
332  {
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.";
335 
336  matchGroomedJets(jets,groomedJets,groomedIndices);
337  }
338 
339  // match subjets and original jets
340  std::vector<std::vector<int> > subjetIndices;
341  if( useSubjets_ )
342  {
343  matchSubjets(groomedIndices,groomedJets,subjets,subjetIndices);
344  }
345 
346  // determine jet flavour
347  for(size_t i=0; i<jets->size(); ++i)
348  {
349  reco::GenParticleRefVector clusteredbHadrons;
350  reco::GenParticleRefVector clusteredcHadrons;
351  reco::GenParticleRefVector clusteredPartons;
352  reco::GenParticleRefVector clusteredLeptons;
353 
354  // if matching reclustered to original jets failed
355  if( reclusteredIndices.at(i) < 0 )
356  {
357  // set an empty JetFlavourInfo for this jet
358  (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
359  }
360  else if( jets->at(i).pt() == 0 )
361  {
362  edm::LogWarning("NullTransverseMomentum") << "The original jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
363 
364  // set an empty JetFlavourInfo for this jet
365  (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
366 
367  // if subjets are used
368  if( useSubjets_ && subjetIndices.at(i).size()>0 )
369  {
370  // loop over subjets
371  for(size_t sj=0; sj<subjetIndices.at(i).size(); ++sj)
372  {
373  // set an empty JetFlavourInfo for this subjet
374  (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(i).at(sj))] = reco::JetFlavourInfo(reco::GenParticleRefVector(), reco::GenParticleRefVector(), reco::GenParticleRefVector(), reco::GenParticleRefVector(), 0, 0);
375  }
376  }
377  }
378  else
379  {
380  // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original jets should in principle stay the same
381  if( ( std::abs( inclusiveJets.at(reclusteredIndices.at(i)).pt() - jets->at(i).pt() ) / jets->at(i).pt() ) > relPtTolerance_ )
382  {
383  if( jets->at(i).pt() < 10. ) // special handling for low-Pt jets (Pt<10 GeV)
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.";
388  else
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.";
392  }
393 
394  // get jet constituents (sorted by Pt)
395  std::vector<fastjet::PseudoJet> constituents = fastjet::sorted_by_pt( inclusiveJets.at(reclusteredIndices.at(i)).constituents() );
396 
397  // loop over jet constituents and try to find "ghosts"
398  for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
399  {
400  if( !it->has_user_info() ) continue; // skip if not a "ghost"
401 
402  // "ghost" hadron
403  if( it->user_info<GhostInfo>().isHadron() )
404  {
405  // "ghost" b hadron
406  if( it->user_info<GhostInfo>().isbHadron() )
407  clusteredbHadrons.push_back(it->user_info<GhostInfo>().particleRef());
408  // "ghost" c hadron
409  else
410  clusteredcHadrons.push_back(it->user_info<GhostInfo>().particleRef());
411  }
412  // "ghost" parton
413  else if( it->user_info<GhostInfo>().isParton() )
414  clusteredPartons.push_back(it->user_info<GhostInfo>().particleRef());
415  // "ghost" lepton
416  else if( it->user_info<GhostInfo>().isLepton() )
417  clusteredLeptons.push_back(it->user_info<GhostInfo>().particleRef());
418  }
419 
420  int hadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
421  int partonFlavour = 0; // default parton flavour set to 0 (= undefined)
422 
423  // set hadron- and parton-based flavours
424  setFlavours(clusteredbHadrons, clusteredcHadrons, clusteredPartons, hadronFlavour, partonFlavour);
425 
426  // set the JetFlavourInfo for this jet
427  (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, hadronFlavour, partonFlavour);
428  }
429 
430  // if subjets are used, determine their flavour
431  if( useSubjets_ )
432  {
433  if( subjetIndices.at(i).size()==0 ) continue; // continue if the original jet does not have subjets assigned
434 
435  // define vectors of GenParticleRefVectors for hadrons and partons assigned to different subjets
436  std::vector<reco::GenParticleRefVector> assignedbHadrons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
437  std::vector<reco::GenParticleRefVector> assignedcHadrons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
438  std::vector<reco::GenParticleRefVector> assignedPartons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
439  std::vector<reco::GenParticleRefVector> assignedLeptons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
440 
441  // loop over clustered b hadrons and assign them to different subjets based on smallest dR
442  assignToSubjets(clusteredbHadrons, subjets, subjetIndices.at(i), assignedbHadrons);
443  // loop over clustered c hadrons and assign them to different subjets based on smallest dR
444  assignToSubjets(clusteredcHadrons, subjets, subjetIndices.at(i), assignedcHadrons);
445  // loop over clustered partons and assign them to different subjets based on smallest dR
446  assignToSubjets(clusteredPartons, subjets, subjetIndices.at(i), assignedPartons);
447  // if used, loop over clustered leptons and assign them to different subjets based on smallest dR
448  if( useLeptons_ )
449  assignToSubjets(clusteredLeptons, subjets, subjetIndices.at(i), assignedLeptons);
450 
451  // loop over subjets and determine their flavour
452  for(size_t sj=0; sj<subjetIndices.at(i).size(); ++sj)
453  {
454  int subjetHadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
455  int subjetPartonFlavour = 0; // default parton flavour set to 0 (= undefined)
456 
457  // set hadron- and parton-based flavours
458  setFlavours(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
459 
460  // set the JetFlavourInfo for this subjet
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);
462  }
463  }
464  }
465 
466  // put jet flavour infos in the event
467  iEvent.put( jetFlavourInfos );
468  // put subjet flavour infos in the event
469  if( useSubjets_ )
470  iEvent.put( subjetFlavourInfos, "SubJets" );
471 }
472 
473 // ------------ method that inserts "ghost" particles in the vector of jet constituents ------------
474 void
476  const double ghostRescaling,
477  const bool isHadron, const bool isbHadron, const bool isParton, const bool isLepton,
478  std::vector<fastjet::PseudoJet>& constituents)
479 {
480  // insert "ghost" particles in the vector of jet constituents
481  for(reco::GenParticleRefVector::const_iterator it = particles->begin(); it != particles->end(); ++it)
482  {
483  if((*it)->pt() == 0)
484  {
485  edm::LogInfo("NullTransverseMomentum") << "dropping input ghost candidate with pt=0";
486  continue;
487  }
488  fastjet::PseudoJet p((*it)->px(),(*it)->py(),(*it)->pz(),(*it)->energy());
489  p*=ghostRescaling; // rescale particle momentum
490  p.set_user_info(new GhostInfo(isHadron, isbHadron, isParton, isLepton, *it));
491  constituents.push_back(p);
492  }
493 }
494 
495 // ------------ method that matches reclustered and original jets based on minimum dR ------------
496 void
498  const std::vector<fastjet::PseudoJet>& reclusteredJets,
499  std::vector<int>& matchedIndices)
500 {
501  std::vector<bool> matchedLocks(reclusteredJets.size(),false);
502 
503  for(size_t j=0; j<jets->size(); ++j)
504  {
505  double matchedDR2 = 1e9;
506  int matchedIdx = -1;
507 
508  for(size_t rj=0; rj<reclusteredJets.size(); ++rj)
509  {
510  if( matchedLocks.at(rj) ) continue; // skip jets that have already been matched
511 
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 )
514  {
515  matchedDR2 = tempDR2;
516  matchedIdx = rj;
517  }
518  }
519 
520  if( matchedIdx>=0 )
521  {
522  if ( matchedDR2 > rParam_*rParam_ )
523  {
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.";
526  }
527  else
528  matchedLocks.at(matchedIdx) = true;
529  }
530  else
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.";
532 
533  matchedIndices.push_back(matchedIdx);
534  }
535 }
536 
537 // ------------ method that matches groomed and original jets based on minimum dR ------------
538 void
540  const edm::Handle<edm::View<reco::Jet> >& groomedJets,
541  std::vector<int>& matchedIndices)
542 {
543  std::vector<bool> jetLocks(jets->size(),false);
544  std::vector<int> jetIndices;
545 
546  for(size_t gj=0; gj<groomedJets->size(); ++gj)
547  {
548  double matchedDR2 = 1e9;
549  int matchedIdx = -1;
550 
551  if( groomedJets->at(gj).pt()>0. ) // skip pathological cases of groomed jets with Pt=0
552  {
553  for(size_t j=0; j<jets->size(); ++j)
554  {
555  if( jetLocks.at(j) ) continue; // skip jets that have already been matched
556 
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 )
559  {
560  matchedDR2 = tempDR2;
561  matchedIdx = j;
562  }
563  }
564  }
565 
566  if( matchedIdx>=0 )
567  {
568  if ( matchedDR2 > rParam_*rParam_ )
569  {
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.";
572  matchedIdx = -1;
573  }
574  else
575  jetLocks.at(matchedIdx) = true;
576  }
577  jetIndices.push_back(matchedIdx);
578  }
579 
580  for(size_t j=0; j<jets->size(); ++j)
581  {
582  std::vector<int>::iterator matchedIndex = std::find( jetIndices.begin(), jetIndices.end(), j );
583 
584  matchedIndices.push_back( matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(),matchedIndex) : -1 );
585  }
586 }
587 
588 // ------------ method that matches subjets and original jets ------------
589 void
590 JetFlavourClustering::matchSubjets(const std::vector<int>& groomedIndices,
591  const edm::Handle<edm::View<reco::Jet> >& groomedJets,
592  const edm::Handle<edm::View<reco::Jet> >& subjets,
593  std::vector<std::vector<int> >& matchedIndices)
594 {
595  for(size_t g=0; g<groomedIndices.size(); ++g)
596  {
597  std::vector<int> subjetIndices;
598 
599  if( groomedIndices.at(g)>=0 )
600  {
601  for(size_t s=0; s<groomedJets->at(groomedIndices.at(g)).numberOfDaughters(); ++s)
602  {
603  const edm::Ptr<reco::Candidate> & subjet = groomedJets->at(groomedIndices.at(g)).daughterPtr(s);
604 
605  for(size_t sj=0; sj<subjets->size(); ++sj)
606  {
607  if( subjet == edm::Ptr<reco::Candidate>(subjets->ptrAt(sj)) )
608  {
609  subjetIndices.push_back(sj);
610  break;
611  }
612  }
613  }
614 
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.";
617 
618  matchedIndices.push_back(subjetIndices);
619  }
620  else
621  matchedIndices.push_back(subjetIndices);
622  }
623 }
624 
625 // ------------ method that sets hadron- and parton-based flavours ------------
626 void
628  const reco::GenParticleRefVector& clusteredcHadrons,
629  const reco::GenParticleRefVector& clusteredPartons,
630  int& hadronFlavour,
631  int& partonFlavour)
632 {
633  reco::GenParticleRef hardestParton;
634  reco::GenParticleRef hardestLightParton;
635  reco::GenParticleRef flavourParton;
636 
637  // loop over clustered partons (already sorted by Pt)
638  for(reco::GenParticleRefVector::const_iterator it = clusteredPartons.begin(); it != clusteredPartons.end(); ++it)
639  {
640  // hardest parton
641  if( hardestParton.isNull() )
642  hardestParton = (*it);
643  // hardest light-flavour parton
644  if( hardestLightParton.isNull() )
645  {
646  if( CandMCTagUtils::isLightParton( *(*it) ) )
647  hardestLightParton = (*it);
648  }
649  // c flavour
650  if( flavourParton.isNull() && ( std::abs( (*it)->pdgId() ) == 4 ) )
651  flavourParton = (*it);
652  // b flavour gets priority
653  if( std::abs( (*it)->pdgId() ) == 5 )
654  {
655  if( flavourParton.isNull() )
656  flavourParton = (*it);
657  else if( std::abs( flavourParton->pdgId() ) != 5 )
658  flavourParton = (*it);
659  }
660  }
661 
662  // set hadron-based flavour
663  if( clusteredbHadrons.size()>0 )
664  hadronFlavour = 5;
665  else if( clusteredcHadrons.size()>0 && clusteredbHadrons.size()==0 )
666  hadronFlavour = 4;
667  // set parton-based flavour
668  if( flavourParton.isNull() )
669  {
670  if( hardestParton.isNonnull() ) partonFlavour = hardestParton->pdgId();
671  }
672  else
673  partonFlavour = flavourParton->pdgId();
674 
675  // if enabled, check for conflicts between hadron- and parton-based flavours and give priority to the hadron-based flavour
677  {
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;
682  }
683 }
684 
685 // ------------ method that assigns clustered particles to subjets ------------
686 void
688  const edm::Handle<edm::View<reco::Jet> >& subjets,
689  const std::vector<int>& subjetIndices,
690  std::vector<reco::GenParticleRefVector>& assignedParticles)
691 {
692  // loop over clustered particles and assign them to different subjets based on smallest dR
693  for(reco::GenParticleRefVector::const_iterator it = clusteredParticles.begin(); it != clusteredParticles.end(); ++it)
694  {
695  std::vector<double> dR2toSubjets;
696 
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() ) );
699 
700  // find the closest subjet
701  int closestSubjetIdx = std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
702 
703  assignedParticles.at(closestSubjetIdx).push_back( *it );
704  }
705 }
706 
707 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
708 void
710  //The following says we do not know what parameters are allowed so do no validation
711  // Please change this to state exactly what you do use, even if it is no parameters
713  desc.setUnknown();
714  descriptions.addDefault(desc);
715 }
716 
717 //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
virtual void produce(edm::Event &, const edm::EventSetup &)
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:252
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:457
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:255
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:250
const bool isLepton() const
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const edm::EDGetTokenT< reco::GenParticleRefVector > cHadronsToken_
int iEvent
Definition: GenABIO.cc:230
void addDefault(ParameterSetDescription const &psetDescription)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:115
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
JetFlavourClustering(const edm::ParameterSet &)
bool isNull() const
Checks for null.
Definition: Ref.h:249
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
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:85
const bool isbHadron() 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:62
size_type size() const
Size of the RefVector.
Definition: RefVector.h:99
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_