Type1MET/plugins/Type2CorrectionProducer.cc

0001 // -*- C++ -*-
0002
0003 //____________________________________________________________________________||
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/stream/EDProducer.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Framework/interface/ConsumesCollector.h"
0010
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012
0013 #include "DataFormats/METReco/interface/CorrMETData.h"
0014
0015 #include <TFormula.h>
0016
0017 #include <iostream>
0018 #include <vector>
0019
0020 //____________________________________________________________________________||
0021 class Type2CorrectionProducer : public edm::stream::EDProducer<> {
0022 public:
0023   explicit Type2CorrectionProducer(const edm::ParameterSet&);
0024   ~Type2CorrectionProducer() override {}
0025
0026 private:
0027   typedef std::vector<edm::InputTag> vInputTag;
0028
0029   typedef std::vector<std::string> vstring;
0030   struct type2BinningEntryType {
0031     type2BinningEntryType(const std::string& binCorrformula,
0032                           const edm::ParameterSet& binCorrParameter,
0033                           const vInputTag& srcUnclEnergySums,
0034                           edm::ConsumesCollector&& iConsumesCollector)
0035         : binLabel_(""), binCorrFormula_(nullptr) {
0036       for (vInputTag::const_iterator inputTag = srcUnclEnergySums.begin(); inputTag != srcUnclEnergySums.end();
0037            ++inputTag) {
0038         corrTokens_.push_back(iConsumesCollector.consumes<CorrMETData>(*inputTag));
0039       }
0040
0041       initialize(binCorrformula, binCorrParameter);
0042     }
0043     type2BinningEntryType(const edm::ParameterSet& cfg,
0044                           const vInputTag& srcUnclEnergySums,
0045                           edm::ConsumesCollector&& iConsumesCollector)
0046         : binLabel_(cfg.getParameter<std::string>("binLabel")), binCorrFormula_(nullptr) {
0047       for (vInputTag::const_iterator srcUnclEnergySum = srcUnclEnergySums.begin();
0048            srcUnclEnergySum != srcUnclEnergySums.end();
0049            ++srcUnclEnergySum) {
0050         std::string instanceLabel = srcUnclEnergySum->instance();
0051         if (!instanceLabel.empty() && !binLabel_.empty())
0052           instanceLabel.append("#");
0053         instanceLabel.append(binLabel_);
0054         edm::InputTag inputTag(srcUnclEnergySum->label(), instanceLabel);
0055         corrTokens_.push_back(iConsumesCollector.consumes<CorrMETData>(inputTag));
0056       }
0057
0058       std::string binCorrFormula = cfg.getParameter<std::string>("binCorrFormula");
0059
0060       edm::ParameterSet binCorrParameter = cfg.getParameter<edm::ParameterSet>("binCorrParameter");
0061
0062       initialize(binCorrFormula, binCorrParameter);
0063     }
0064     void initialize(const std::string& binCorrFormula, const edm::ParameterSet& binCorrParameter) {
0065       TString formula = binCorrFormula;
0066
0067       vstring parNames = binCorrParameter.getParameterNamesForType<double>();
0068       int numParameter = parNames.size();
0069       binCorrParameter_.resize(numParameter);
0070       for (int parIndex = 0; parIndex < numParameter; ++parIndex) {
0071         const std::string& parName = parNames[parIndex];
0072
0073         double parValue = binCorrParameter.getParameter<double>(parName);
0074         binCorrParameter_[parIndex] = parValue;
0075
0076         TString parName_internal = Form("[%i]", parIndex);
0077         formula = formula.ReplaceAll(parName.data(), parName_internal);
0078       }
0079
0080       std::string binCorrFormulaName = std::string("binCorrFormula").append("_").append(binLabel_);
0081       binCorrFormula_ = new TFormula(binCorrFormulaName.data(), formula);
0082
0083       //--- check that syntax of formula string is valid
0084       //   (i.e. that TFormula "compiled" without errors)
0085       if (!(binCorrFormula_->GetNdim() <= 1 && binCorrFormula_->GetNpar() == numParameter))
0086         throw cms::Exception("METCorrectionAlgorithm2")
0087             << "Formula for Type 2 correction has invalid syntax = " << formula << " !!\n";
0088
0089       for (int parIndex = 0; parIndex < numParameter; ++parIndex) {
0090         binCorrFormula_->SetParameter(parIndex, binCorrParameter_[parIndex]);
0091         binCorrFormula_->SetParName(parIndex, parNames[parIndex].data());
0092       }
0093     }
0094     ~type2BinningEntryType() { delete binCorrFormula_; }
0095     std::string binLabel_;
0096     std::vector<edm::EDGetTokenT<CorrMETData> > corrTokens_;
0097     TFormula* binCorrFormula_;
0098     std::vector<double> binCorrParameter_;
0099   };
0100
0101   std::vector<type2BinningEntryType*> type2Binning_;
0102
0103   void produce(edm::Event&, const edm::EventSetup&) override;
0104 };
0105
0106 //____________________________________________________________________________||
0107 Type2CorrectionProducer::Type2CorrectionProducer(const edm::ParameterSet& cfg) {
0108   vInputTag srcUnclEnergySums = cfg.getParameter<vInputTag>("srcUnclEnergySums");
0109
0110   if (cfg.exists("type2Binning")) {
0111     typedef std::vector<edm::ParameterSet> vParameterSet;
0112     vParameterSet cfgType2Binning = cfg.getParameter<vParameterSet>("type2Binning");
0113     for (vParameterSet::const_iterator cfgType2BinningEntry = cfgType2Binning.begin();
0114          cfgType2BinningEntry != cfgType2Binning.end();
0115          ++cfgType2BinningEntry) {
0116       type2Binning_.push_back(new type2BinningEntryType(*cfgType2BinningEntry, srcUnclEnergySums, consumesCollector()));
0117     }
0118   } else {
0119     std::string type2CorrFormula = cfg.getParameter<std::string>("type2CorrFormula");
0120     edm::ParameterSet type2CorrParameter = cfg.getParameter<edm::ParameterSet>("type2CorrParameter");
0121     type2Binning_.push_back(
0122         new type2BinningEntryType(type2CorrFormula, type2CorrParameter, srcUnclEnergySums, consumesCollector()));
0123   }
0124   produces<CorrMETData>("");
0125 }
0126
0127 //____________________________________________________________________________||
0128 void Type2CorrectionProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0129   CorrMETData product;
0130
0131   for (std::vector<type2BinningEntryType*>::const_iterator type2BinningEntry = type2Binning_.begin();
0132        type2BinningEntry != type2Binning_.end();
0133        ++type2BinningEntry) {
0134     CorrMETData unclEnergySum;
0135
0136     edm::Handle<CorrMETData> unclEnergySummand;
0137     for (std::vector<edm::EDGetTokenT<CorrMETData> >::const_iterator corrToken =
0138              (*type2BinningEntry)->corrTokens_.begin();
0139          corrToken != (*type2BinningEntry)->corrTokens_.end();
0140          ++corrToken) {
0141       evt.getByToken(*corrToken, unclEnergySummand);
0142       unclEnergySum += (*unclEnergySummand);
0143     }
0144
0145     double unclEnergySumPt = sqrt(unclEnergySum.mex * unclEnergySum.mex + unclEnergySum.mey * unclEnergySum.mey);
0146     double unclEnergyScaleFactor = (*type2BinningEntry)->binCorrFormula_->Eval(unclEnergySumPt);
0147
0148     unclEnergySum.mex = -unclEnergySum.mex;
0149     unclEnergySum.mey = -unclEnergySum.mey;
0150
0151     product += (unclEnergyScaleFactor - 1.) * unclEnergySum;
0152   }
0153
0154   std::unique_ptr<CorrMETData> pprod(new CorrMETData(product));
0155   evt.put(std::move(pprod), "");
0156 }
0157
0158 //____________________________________________________________________________||
0159
0160 DEFINE_FWK_MODULE(Type2CorrectionProducer);