JetDeltaRValueMapProducer

Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
/* \class JetDeltaRValueMapProducer
 *
 * Associates jets using delta-R matching, and writes out
 * a valuemap of single float variables based on a StringObjectFunction. 
 * This is used for miniAOD. 
 *
 * \author: Sal Rappoccio
 *
 *
 * for more details about the cut syntax, see the documentation
 * page below:
 *
 *   https://twiki.cern.ch/twiki/bin/view/CMS/SWGuidePhysicsCutParser
 *
 *
 */

#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/JetReco/interface/Jet.h"
#include "DataFormats/JetReco/interface/PFJet.h"
#include "DataFormats/JetReco/interface/CaloJet.h"

#include "FWCore/Framework/interface/MakerMacros.h"
#include "CommonTools/Utils/interface/StringObjectFunction.h"
#include "FWCore/Framework/interface/Event.h"
#include "DataFormats/Common/interface/ValueMap.h"
#include "DataFormats/Math/interface/deltaR.h"

template <class T, class C = T, typename TN = float>
class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> {
public:
  typedef TN value_type;
  typedef edm::ValueMap<value_type> JetValueMap;
  typedef std::vector<value_type> ValueVector;

  JetDeltaRValueMapProducer(edm::ParameterSet const& params)
      : srcToken_(consumes<typename edm::View<T> >(params.getParameter<edm::InputTag>("src"))),
        matchedToken_(consumes<typename edm::View<C> >(params.getParameter<edm::InputTag>("matched"))),
        distMax_(params.getParameter<double>("distMax")),
        value_(params.existsAs<std::string>("value") ? params.getParameter<std::string>("value") : ""),
        values_(params.existsAs<std::vector<std::string> >("values")
                    ? params.getParameter<std::vector<std::string> >("values")
                    : std::vector<std::string>()),
        valueLabels_(params.existsAs<std::vector<std::string> >("valueLabels")
                         ? params.getParameter<std::vector<std::string> >("valueLabels")
                         : std::vector<std::string>()),
        lazyParser_(params.existsAs<bool>("lazyParser") ? params.getParameter<bool>("lazyParser") : false),
        multiValue_(false) {
    if (!value_.empty()) {
      if (value_ != "index") {
        evaluationMap_.insert(std::make_pair(
            value_, std::unique_ptr<StringObjectFunction<C> >(new StringObjectFunction<C>(value_, lazyParser_))));
      }
      produces<JetValueMap>();
    }

    if (!valueLabels_.empty() || !values_.empty()) {
      if (valueLabels_.size() == values_.size()) {
        multiValue_ = true;
        for (size_t i = 0; i < valueLabels_.size(); ++i) {
          evaluationMap_.insert(std::make_pair(
              valueLabels_[i],
              std::unique_ptr<StringObjectFunction<C> >(new StringObjectFunction<C>(values_[i], lazyParser_))));
          produces<JetValueMap>(valueLabels_[i]);
        }
      } else
        edm::LogWarning("ValueLabelMismatch")
            << "The number of value labels does not match the number of values. Values will not be evaluated.";
    }
  }

  ~JetDeltaRValueMapProducer() override {}

private:
  void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override {
    edm::Handle<typename edm::View<T> > h_jets1;
    iEvent.getByToken(srcToken_, h_jets1);
    edm::Handle<typename edm::View<C> > h_jets2;
    iEvent.getByToken(matchedToken_, h_jets2);

    ValueVector values(h_jets1->size(), -99999);

    std::map<std::string, std::vector<float> > valuesMap;
    if (multiValue_) {
      for (size_t i = 0; i < valueLabels_.size(); ++i)
        valuesMap.insert(std::make_pair(valueLabels_[i], std::vector<float>(h_jets1->size(), -99999)));
    }
    std::vector<bool> jets1_locks(h_jets1->size(), false);

    for (typename edm::View<C>::const_iterator ibegin = h_jets2->begin(), iend = h_jets2->end(), ijet = ibegin;
         ijet != iend;
         ++ijet) {
      float matched_dR2 = 1e9;
      int matched_index = -1;
      int iindex = ijet - ibegin;

      for (typename edm::View<T>::const_iterator jbegin = h_jets1->begin(), jend = h_jets1->end(), jjet = jbegin;
           jjet != jend;
           ++jjet) {
        int jindex = jjet - jbegin;

        if (jets1_locks.at(jindex))
          continue;  // skip jets that have already been matched

        float temp_dR2 = reco::deltaR2(ijet->eta(), ijet->phi(), jjet->eta(), jjet->phi());
        if (temp_dR2 < matched_dR2) {
          matched_dR2 = temp_dR2;
          matched_index = jindex;
        }
      }  // end loop over src jets

      if (matched_index >= 0) {
        if (matched_dR2 > distMax_ * distMax_)
          edm::LogInfo("MatchedJetsFarApart") << "Matched jets separated by dR greater than distMax=" << distMax_;
        else {
          jets1_locks.at(matched_index) = true;
          if (!value_.empty()) {
            if (value_ == "index") {
              values.at(matched_index) = iindex;
            } else {
              values.at(matched_index) = (*(evaluationMap_.at(value_)))(*ijet);
            }
          }
          if (multiValue_) {
            for (size_t i = 0; i < valueLabels_.size(); ++i)
              valuesMap.at(valueLabels_[i]).at(matched_index) = (*(evaluationMap_.at(valueLabels_[i])))(*ijet);
          }
        }
      }
    }  // end loop over matched jets

    if (!value_.empty()) {
      std::unique_ptr<JetValueMap> jetValueMap(new JetValueMap());

      typename JetValueMap::Filler filler(*jetValueMap);
      filler.insert(h_jets1, values.begin(), values.end());
      filler.fill();

      // put in Event
      iEvent.put(std::move(jetValueMap));
    }
    if (multiValue_) {
      for (size_t i = 0; i < valueLabels_.size(); ++i) {
        std::unique_ptr<JetValueMap> jetValueMap(new JetValueMap());

        typename JetValueMap::Filler filler(*jetValueMap);
        filler.insert(h_jets1, valuesMap.at(valueLabels_[i]).begin(), valuesMap.at(valueLabels_[i]).end());
        filler.fill();

        // put in Event
        iEvent.put(std::move(jetValueMap), valueLabels_[i]);
      }
    }
  }

  const edm::EDGetTokenT<typename edm::View<T> > srcToken_;
  const edm::EDGetTokenT<typename edm::View<C> > matchedToken_;
  const double distMax_;
  const std::string value_;
  const std::vector<std::string> values_;
  const std::vector<std::string> valueLabels_;
  const bool lazyParser_;
  bool multiValue_;
  std::map<std::string, std::unique_ptr<const StringObjectFunction<C> > > evaluationMap_;
};

typedef JetDeltaRValueMapProducer<reco::Jet> RecoJetDeltaRValueMapProducer;
typedef JetDeltaRValueMapProducer<reco::Jet, pat::Jet> RecoJetToPatJetDeltaRValueMapProducer;
typedef JetDeltaRValueMapProducer<pat::Jet> PatJetDeltaRValueMapProducer;
typedef JetDeltaRValueMapProducer<reco::Jet, reco::GenJet, int> RecoJetToGenJetDeltaRValueMapProducer;

DEFINE_FWK_MODULE(RecoJetDeltaRValueMapProducer);
DEFINE_FWK_MODULE(RecoJetToPatJetDeltaRValueMapProducer);
DEFINE_FWK_MODULE(PatJetDeltaRValueMapProducer);
DEFINE_FWK_MODULE(RecoJetToGenJetDeltaRValueMapProducer);