HiFJRhoProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    HiJetBackground/HiFJRhoProducer
// Class:      HiFJRhoProducer
//
//
// Original Author:  Marta Verweij
//         Created:  Thu, 16 Jul 2015 10:57:12 GMT
//
//

#include <memory>
#include <string>

#include "RecoHI/HiJetAlgos/plugins/HiFJRhoProducer.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"

#include "DataFormats/Common/interface/View.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/PatCandidates/interface/Jet.h"

//using namespace edm;
using namespace reco;
//using namespace pat;

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
HiFJRhoProducer::HiFJRhoProducer(const edm::ParameterSet& iConfig)
    : src_(iConfig.getParameter<edm::InputTag>("jetSource")),
      nExcl_(iConfig.getParameter<int>("nExcl")),
      etaMaxExcl_(iConfig.getParameter<double>("etaMaxExcl")),
      ptMinExcl_(iConfig.getParameter<double>("ptMinExcl")),
      nExcl2_(iConfig.getParameter<int>("nExcl2")),
      etaMaxExcl2_(iConfig.getParameter<double>("etaMaxExcl2")),
      ptMinExcl2_(iConfig.getParameter<double>("ptMinExcl2")),
      checkJetCand(true),
      usingPackedCand(false) {
  jetsToken_ = consumes<edm::View<reco::Jet>>(src_);

  //register your products
  produces<std::vector<double>>("mapEtaEdges");
  produces<std::vector<double>>("mapToRho");
  produces<std::vector<double>>("mapToRhoM");
  produces<std::vector<double>>("ptJets");
  produces<std::vector<double>>("areaJets");
  produces<std::vector<double>>("etaJets");
  etaRanges = iConfig.getParameter<std::vector<double>>("etaRanges");
}

HiFJRhoProducer::~HiFJRhoProducer() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}

// ------------ method called to produce the data  ------------
void HiFJRhoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  // Get the vector of jets
  edm::Handle<edm::View<reco::Jet>> jets;
  iEvent.getByToken(jetsToken_, jets);

  int neta = (int)etaRanges.size();
  auto mapEtaRangesOut = std::make_unique<std::vector<double>>(neta, -999.);

  for (int ieta = 0; ieta < neta; ieta++) {
    mapEtaRangesOut->at(ieta) = etaRanges[ieta];
  }
  auto mapToRhoOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
  auto mapToRhoMOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);

  int njets = jets->size();

  auto ptJetsOut = std::make_unique<std::vector<double>>(njets, 1e-6);
  auto areaJetsOut = std::make_unique<std::vector<double>>(njets, 1e-6);
  auto etaJetsOut = std::make_unique<std::vector<double>>(njets, 1e-6);

  double rhoVec[999];
  double rhomVec[999];
  double etaVec[999];

  // int neta = (int)mapEtaRangesOut->size();
  int nacc = 0;
  unsigned int njetsEx = 0;
  unsigned int njetsEx2 = 0;
  for (auto jet = jets->begin(); jet != jets->end(); ++jet) {
    if (njetsEx < nExcl_ && fabs(jet->eta()) < etaMaxExcl_ && jet->pt() > ptMinExcl_) {
      njetsEx++;
      continue;
    }
    if (njetsEx2 < nExcl2_ && fabs(jet->eta()) < etaMaxExcl2_ && fabs(jet->eta()) > etaMaxExcl_ &&
        jet->pt() > ptMinExcl2_) {
      njetsEx2++;
      continue;
    }
    float pt = jet->pt();
    float area = jet->jetArea();
    float eta = jet->eta();

    if (eta < mapEtaRangesOut->at(0) || eta > mapEtaRangesOut->at(neta - 1))
      continue;
    if (area > 0.) {
      rhoVec[nacc] = pt / area;
      rhomVec[nacc] = calcMd(&*jet) / area;
      etaVec[nacc] = eta;
      ptJetsOut->at(nacc) = pt;
      areaJetsOut->at(nacc) = area;
      etaJetsOut->at(nacc) = eta;
      ++nacc;
    }
  }

  ptJetsOut->resize(nacc);
  areaJetsOut->resize(nacc);
  etaJetsOut->resize(nacc);
  //calculate rho and rhom in eta ranges
  double radius = 0.2;  //distance kt clusters needs to be from edge
  for (int ieta = 0; ieta < (neta - 1); ieta++) {
    std::vector<double> rhoVecCur;
    std::vector<double> rhomVecCur;

    double etaMin = mapEtaRangesOut->at(ieta) + radius;
    double etaMax = mapEtaRangesOut->at(ieta + 1) - radius;

    int naccCur = 0;
    double rhoCurSum = 0.;
    double rhomCurSum = 0.;
    for (int i = 0; i < nacc; i++) {
      if (etaVec[i] >= etaMin && etaVec[i] < etaMax) {
        rhoVecCur.push_back(rhoVec[i]);
        rhomVecCur.push_back(rhomVec[i]);

        rhoCurSum += rhoVec[i];
        rhomCurSum += rhomVec[i];
        ++naccCur;
      }  //eta selection
    }    //accepted jet loop

    if (naccCur > 0) {
      double rhoCur = calcMedian(rhoVecCur);
      double rhomCur = calcMedian(rhomVecCur);
      mapToRhoOut->at(ieta) = rhoCur;
      mapToRhoMOut->at(ieta) = rhomCur;
    }
  }  //eta ranges

  iEvent.put(std::move(mapEtaRangesOut), "mapEtaEdges");
  iEvent.put(std::move(mapToRhoOut), "mapToRho");
  iEvent.put(std::move(mapToRhoMOut), "mapToRhoM");

  iEvent.put(std::move(ptJetsOut), "ptJets");
  iEvent.put(std::move(areaJetsOut), "areaJets");
  iEvent.put(std::move(etaJetsOut), "etaJets");

  return;
}

// ------------ method called once each stream before processing any runs, lumis or events  ------------
void HiFJRhoProducer::beginStream(edm::StreamID) {}

// ------------ method called once each stream after processing all runs, lumis and events  ------------
void HiFJRhoProducer::endStream() {}

double HiFJRhoProducer::calcMd(const reco::Jet* jet) {
  //
  //get md as defined in http://arxiv.org/pdf/1211.2811.pdf
  //

  //Loop over the jet constituents
  double sum = 0.;
  for (auto daughter : jet->getJetConstituentsQuick()) {
    if (isPackedCandidate(daughter)) {  //packed candidate situation
      auto part = static_cast<const pat::PackedCandidate*>(daughter);
      sum += sqrt(part->mass() * part->mass() + part->pt() * part->pt()) - part->pt();
    } else {
      auto part = static_cast<const reco::PFCandidate*>(daughter);
      sum += sqrt(part->mass() * part->mass() + part->pt() * part->pt()) - part->pt();
    }
  }

  return sum;
}

bool HiFJRhoProducer::isPackedCandidate(const reco::Candidate* candidate) {
  if (checkJetCand) {
    if (typeid(pat::PackedCandidate) == typeid(*candidate))
      usingPackedCand = true;
    else if (typeid(reco::PFCandidate) == typeid(*candidate))
      usingPackedCand = false;
    else
      throw cms::Exception("WrongJetCollection", "Jet constituents are not particle flow candidates");
    checkJetCand = false;
  }
  return usingPackedCand;
}

void HiFJRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("jetSource", edm::InputTag("kt4PFJets"));
  desc.add<int>("nExcl", 2);
  desc.add<double>("etaMaxExcl", 2.);
  desc.add<double>("ptMinExcl", 20.);
  desc.add<int>("nExcl2", 2);
  desc.add<double>("etaMaxExcl2", 2.);
  desc.add<double>("ptMinExcl2", 20.);
  desc.add<std::vector<double>>("etaRanges", {});
  descriptions.add("hiFJRhoProducer", desc);
}

//--------- method to calculate median ------------------
double HiFJRhoProducer::calcMedian(std::vector<double>& v) {
  //post-condition: After returning, the elements in v may be reordered and the resulting order is implementation defined.
  //works for even and odd collections
  if (v.empty()) {
    return 0.0;
  }
  auto n = v.size() / 2;
  std::nth_element(v.begin(), v.begin() + n, v.end());
  auto med = v[n];
  if (!(v.size() & 1)) {  //If the set size is even
    auto max_it = std::max_element(v.begin(), v.begin() + n);
    med = (*max_it + med) / 2.0;
  }
  return med;
}

//define this as a plug-in
DEFINE_FWK_MODULE(HiFJRhoProducer);