NjettinessAdder.cc

Go to the documentation of this file.

#include "RecoJets/JetProducers/interface/NjettinessAdder.h"

#include "FWCore/Framework/interface/MakerMacros.h"

NjettinessAdder::NjettinessAdder(const edm::ParameterSet& iConfig)
    : src_(iConfig.getParameter<edm::InputTag>("src")),
      src_token_(consumes<edm::View<reco::Jet>>(src_)),
      Njets_(iConfig.getParameter<std::vector<unsigned>>("Njets")),
      measureDefinition_(iConfig.getParameter<unsigned>("measureDefinition")),
      beta_(iConfig.getParameter<double>("beta")),
      R0_(iConfig.getParameter<double>("R0")),
      Rcutoff_(iConfig.getParameter<double>("Rcutoff")),
      axesDefinition_(iConfig.getParameter<unsigned>("axesDefinition")),
      nPass_(iConfig.getParameter<int>("nPass")),
      akAxesR0_(iConfig.getParameter<double>("akAxesR0")) {
  for (std::vector<unsigned>::const_iterator n = Njets_.begin(); n != Njets_.end(); ++n) {
    std::ostringstream tauN_str;
    tauN_str << "tau" << *n;

    produces<edm::ValueMap<float>>(tauN_str.str());
  }

  // Get the measure definition
  fastjet::contrib::NormalizedMeasure normalizedMeasure(beta_, R0_);
  fastjet::contrib::UnnormalizedMeasure unnormalizedMeasure(beta_);
  fastjet::contrib::OriginalGeometricMeasure geometricMeasure(beta_);  // changed in 1.020
  fastjet::contrib::NormalizedCutoffMeasure normalizedCutoffMeasure(beta_, R0_, Rcutoff_);
  fastjet::contrib::UnnormalizedCutoffMeasure unnormalizedCutoffMeasure(beta_, Rcutoff_);
  //fastjet::contrib::GeometricCutoffMeasure     geometricCutoffMeasure   (beta_,Rcutoff_); // removed in 1.020

  fastjet::contrib::MeasureDefinition const* measureDef = nullptr;
  switch (measureDefinition_) {
    case UnnormalizedMeasure:
      measureDef = &unnormalizedMeasure;
      break;
    case OriginalGeometricMeasure:
      measureDef = &geometricMeasure;
      break;  // changed in 1.020
    case NormalizedCutoffMeasure:
      measureDef = &normalizedCutoffMeasure;
      break;
    case UnnormalizedCutoffMeasure:
      measureDef = &unnormalizedCutoffMeasure;
      break;
    //case GeometricCutoffMeasure : measureDef = &geometricCutoffMeasure; break; // removed in 1.020
    case NormalizedMeasure:
    default:
      measureDef = &normalizedMeasure;
      break;
  }

  // Get the axes definition
  fastjet::contrib::KT_Axes kt_axes;
  fastjet::contrib::CA_Axes ca_axes;
  fastjet::contrib::AntiKT_Axes antikt_axes(akAxesR0_);
  fastjet::contrib::WTA_KT_Axes wta_kt_axes;
  fastjet::contrib::WTA_CA_Axes wta_ca_axes;
  fastjet::contrib::OnePass_KT_Axes onepass_kt_axes;
  fastjet::contrib::OnePass_CA_Axes onepass_ca_axes;
  fastjet::contrib::OnePass_AntiKT_Axes onepass_antikt_axes(akAxesR0_);
  fastjet::contrib::OnePass_WTA_KT_Axes onepass_wta_kt_axes;
  fastjet::contrib::OnePass_WTA_CA_Axes onepass_wta_ca_axes;
  fastjet::contrib::MultiPass_Axes multipass_axes(nPass_);

  fastjet::contrib::AxesDefinition const* axesDef = nullptr;
  switch (axesDefinition_) {
    case KT_Axes:
    default:
      axesDef = &kt_axes;
      break;
    case CA_Axes:
      axesDef = &ca_axes;
      break;
    case AntiKT_Axes:
      axesDef = &antikt_axes;
      break;
    case WTA_KT_Axes:
      axesDef = &wta_kt_axes;
      break;
    case WTA_CA_Axes:
      axesDef = &wta_ca_axes;
      break;
    case OnePass_KT_Axes:
      axesDef = &onepass_kt_axes;
      break;
    case OnePass_CA_Axes:
      axesDef = &onepass_ca_axes;
      break;
    case OnePass_AntiKT_Axes:
      axesDef = &onepass_antikt_axes;
      break;
    case OnePass_WTA_KT_Axes:
      axesDef = &onepass_wta_kt_axes;
      break;
    case OnePass_WTA_CA_Axes:
      axesDef = &onepass_wta_ca_axes;
      break;
    case MultiPass_Axes:
      axesDef = &multipass_axes;
      break;
  };

  routine_ = std::unique_ptr<fastjet::contrib::Njettiness>(new fastjet::contrib::Njettiness(*axesDef, *measureDef));
}

void NjettinessAdder::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  // read input collection
  edm::Handle<edm::View<reco::Jet>> jets;
  iEvent.getByToken(src_token_, jets);

  for (std::vector<unsigned>::const_iterator n = Njets_.begin(); n != Njets_.end(); ++n) {
    std::ostringstream tauN_str;
    tauN_str << "tau" << *n;

    // prepare room for output
    std::vector<float> tauN;
    tauN.reserve(jets->size());

    for (typename edm::View<reco::Jet>::const_iterator jetIt = jets->begin(); jetIt != jets->end(); ++jetIt) {
      edm::Ptr<reco::Jet> jetPtr = jets->ptrAt(jetIt - jets->begin());

      float t = getTau(*n, jetPtr);

      tauN.push_back(t);
    }

    auto outT = std::make_unique<edm::ValueMap<float>>();
    edm::ValueMap<float>::Filler fillerT(*outT);
    fillerT.insert(jets, tauN.begin(), tauN.end());
    fillerT.fill();

    iEvent.put(std::move(outT), tauN_str.str());
  }
}

float NjettinessAdder::getTau(unsigned num, const edm::Ptr<reco::Jet>& object) const {
  std::vector<fastjet::PseudoJet> FJparticles;
  for (unsigned k = 0; k < object->numberOfDaughters(); ++k) {
    const reco::CandidatePtr& dp = object->daughterPtr(k);
    if (dp.isNonnull() && dp.isAvailable())
      FJparticles.push_back(fastjet::PseudoJet(dp->px(), dp->py(), dp->pz(), dp->energy()));
    else
      edm::LogWarning("MissingJetConstituent") << "Jet constituent required for N-subjettiness computation is missing!";
  }

  return routine_->getTau(num, FJparticles);
}

DEFINE_FWK_MODULE(NjettinessAdder);