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