GenJetClosestMatchSelectorDefinition.h

Go to the documentation of this file.

#ifndef PhysicsTools_PFCandProducer_GenJetClosestMatchSelectorDefinition
#define PhysicsTools_PFCandProducer_GenJetClosestMatchSelectorDefinition
 
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
 
#include "DataFormats/JetReco/interface/GenJet.h"
 
#include "DataFormats/Math/interface/deltaR.h"
 
#include <iostream>
 
struct GenJetClosestMatchSelectorDefinition {
  typedef reco::GenJetCollection collection;
  typedef edm::Handle<collection> HandleToCollection;
  typedef std::vector<reco::GenJet *> container;
  typedef container::const_iterator const_iterator;
 
  GenJetClosestMatchSelectorDefinition(const edm::ParameterSet &cfg, edm::ConsumesCollector &&iC) {
    matchTo_ = iC.consumes<edm::View<reco::Candidate>>(cfg.getParameter<edm::InputTag>("MatchTo"));
  }
 
  const_iterator begin() const { return selected_.begin(); }
 
  const_iterator end() const { return selected_.end(); }
 
  void select(const HandleToCollection &hc, const edm::Event &e, const edm::EventSetup &s) {
    selected_.clear();
 
    edm::Handle<edm::View<reco::Candidate>> matchCandidates;
    e.getByToken(matchTo_, matchCandidates);
 
    unsigned key = 0;
 
    //    std::cout<<"number of candidates
    //    "<<matchCandidates->size()<<std::endl;
 
    typedef edm::View<reco::Candidate>::const_iterator IC;
    for (IC ic = matchCandidates->begin(); ic != matchCandidates->end(); ++ic) {
      double eta2 = ic->eta();
      double phi2 = ic->phi();
 
      //      std::cout<<"cand "<<eta2<<" "<<phi2<<std::endl;
 
      // look for the closest gen jet
      double deltaR2Min = 9999;
      collection::const_iterator closest = hc->end();
      for (collection::const_iterator genjet = hc->begin(); genjet != hc->end(); ++genjet, ++key) {
        reco::GenJetRef genJetRef(hc, key);
 
        // is it matched?
 
        double eta1 = genjet->eta();
        double phi1 = genjet->phi();
 
        double deltaR2 = reco::deltaR2(eta1, phi1, eta2, phi2);
 
        // std::cout<<"  genjet "<<eta1<<" "<<phi1<<" "<<deltaR2<<std::endl;
 
        // cut should be a parameter
        if (deltaR2 < deltaR2Min) {
          deltaR2Min = deltaR2;
          closest = genjet;
        }
      }
 
      if (deltaR2Min < 0.01) {
        // std::cout<<deltaR2Min<<std::endl;
        selected_.push_back(new reco::GenJet(*closest));
      }
    }  // end collection iteration
 
    // std::cout<<selected_.size()<<std::endl;
  }  // end select()
 
  size_t size() const { return selected_.size(); }
 
private:
  container selected_;
  edm::EDGetTokenT<edm::View<reco::Candidate>> matchTo_;
};
 
#endif