Type1MET/plugins/MultShiftMETcorrInputProducer.cc

0001 #include "JetMETCorrections/Type1MET/plugins/MultShiftMETcorrInputProducer.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004
0005 #include "DataFormats/METReco/interface/CorrMETData.h"
0006 #include "DataFormats/VertexReco/interface/Vertex.h"
0007 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0008 #include "DataFormats/Common/interface/View.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include <TString.h>
0012
0013 #include <memory>
0014
0015 int MultShiftMETcorrInputProducer::translateTypeToAbsPdgId(reco::PFCandidate::ParticleType type) {
0016   switch (type) {
0017     case reco::PFCandidate::ParticleType::h:
0018       return 211;  // pi+
0019     case reco::PFCandidate::ParticleType::e:
0020       return 11;
0021     case reco::PFCandidate::ParticleType::mu:
0022       return 13;
0023     case reco::PFCandidate::ParticleType::gamma:
0024       return 22;
0025     case reco::PFCandidate::ParticleType::h0:
0026       return 130;  // K_L0
0027     case reco::PFCandidate::ParticleType::h_HF:
0028       return 1;  // dummy pdg code
0029     case reco::PFCandidate::ParticleType::egamma_HF:
0030       return 2;  // dummy pdg code
0031     case reco::PFCandidate::ParticleType::X:
0032     default:
0033       return 0;
0034   }
0035 }
0036
0037 MultShiftMETcorrInputProducer::MultShiftMETcorrInputProducer(const edm::ParameterSet& cfg)
0038     : moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
0039   pflow_ = consumes<edm::View<reco::Candidate>>(cfg.getParameter<edm::InputTag>("srcPFlow"));
0040   vertices_ = consumes<edm::View<reco::Vertex>>(cfg.getParameter<edm::InputTag>("vertexCollection"));
0041
0042   edm::InputTag srcWeights = cfg.getParameter<edm::InputTag>("srcWeights");
0043   if (!srcWeights.label().empty())
0044     weightsToken_ = consumes<edm::ValueMap<float>>(srcWeights);
0045
0046   cfgCorrParameters_ = cfg.getParameter<std::vector<edm::ParameterSet>>("parameters");
0047   etaMin_.clear();
0048   etaMax_.clear();
0049   type_.clear();
0050   varType_.clear();
0051
0052   produces<CorrMETData>();
0053
0054   for (std::vector<edm::ParameterSet>::const_iterator v = cfgCorrParameters_.begin(); v != cfgCorrParameters_.end();
0055        v++) {
0056     TString corrPxFormula = v->getParameter<std::string>("fx");
0057     TString corrPyFormula = v->getParameter<std::string>("fy");
0058     std::vector<double> corrPxParams = v->getParameter<std::vector<double>>("px");
0059     std::vector<double> corrPyParams = v->getParameter<std::vector<double>>("py");
0060
0061     formula_x_.push_back(std::make_unique<TF1>(
0062         std::string(moduleLabel_).append("_").append(v->getParameter<std::string>("name")).append("_corrPx").c_str(),
0063         v->getParameter<std::string>("fx").c_str()));
0064     formula_y_.push_back(std::make_unique<TF1>(
0065         std::string(moduleLabel_).append("_").append(v->getParameter<std::string>("name")).append("_corrPy").c_str(),
0066         v->getParameter<std::string>("fy").c_str()));
0067
0068     for (unsigned i = 0; i < corrPxParams.size(); i++)
0069       formula_x_.back()->SetParameter(i, corrPxParams[i]);
0070     for (unsigned i = 0; i < corrPyParams.size(); i++)
0071       formula_y_.back()->SetParameter(i, corrPyParams[i]);
0072
0073     counts_.push_back(0);
0074     sumPt_.push_back(0.);
0075     etaMin_.push_back(v->getParameter<double>("etaMin"));
0076     etaMax_.push_back(v->getParameter<double>("etaMax"));
0077     type_.push_back(v->getParameter<int>("type"));
0078     varType_.push_back(v->getParameter<int>("varType"));
0079   }
0080 }
0081
0082 MultShiftMETcorrInputProducer::~MultShiftMETcorrInputProducer() {}
0083
0084 void MultShiftMETcorrInputProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0085   //get primary vertices
0086   edm::Handle<edm::View<reco::Vertex>> hpv;
0087   evt.getByToken(vertices_, hpv);
0088   if (!hpv.isValid()) {
0089     edm::LogError("MultShiftMETcorrInputProducer::produce") << "could not find vertex collection ";
0090   }
0091   uint ngoodVertices = 0;
0092   for (unsigned i = 0; i < hpv->size(); i++) {
0093     if ((*hpv)[i].ndof() > 4 && (fabs((*hpv)[i].z()) <= 24.) && (fabs((*hpv)[i].position().rho()) <= 2.0))
0094       ngoodVertices += 1;
0095   }
0096   uint nVertices = hpv->size();
0097
0098   for (unsigned i = 0; i < counts_.size(); i++)
0099     counts_[i] = 0;
0100   for (unsigned i = 0; i < sumPt_.size(); i++)
0101     sumPt_[i] = 0.;
0102
0103   edm::Handle<edm::View<reco::Candidate>> particleFlow;
0104   evt.getByToken(pflow_, particleFlow);
0105
0106   edm::Handle<edm::ValueMap<float>> weights;
0107   if (!weightsToken_.isUninitialized())
0108     evt.getByToken(weightsToken_, weights);
0109   if (std::find(varType_.begin(), varType_.end(), 0) != varType_.end() ||
0110       std::find(varType_.begin(), varType_.end(), 2) != varType_.end()) {
0111     for (unsigned i = 0; i < particleFlow->size(); ++i) {
0112       const reco::Candidate& c = particleFlow->at(i);
0113       for (unsigned j = 0; j < type_.size(); j++) {
0114         if (abs(c.pdgId()) == translateTypeToAbsPdgId(reco::PFCandidate::ParticleType(type_[j]))) {
0115           if ((c.eta() > etaMin_[j]) and (c.eta() < etaMax_[j])) {
0116             float weight = (!weightsToken_.isUninitialized()) ? (*weights)[particleFlow->ptrAt(i)] : 1.0;
0117             counts_[j] += (weight > 0);
0118             sumPt_[j] += c.pt() * weight;
0119             continue;
0120           }
0121         }
0122       }
0123     }
0124   }
0125
0126   //MM: loop over all constituent types and sum each correction
0127   std::unique_ptr<CorrMETData> metCorr(new CorrMETData());
0128
0129   double corx = 0.;
0130   double cory = 0.;
0131
0132   for (std::vector<edm::ParameterSet>::const_iterator v = cfgCorrParameters_.begin(); v != cfgCorrParameters_.end();
0133        v++) {
0134     unsigned j = v - cfgCorrParameters_.begin();
0135
0136     double val(0.);
0137     if (varType_[j] == 0) {
0138       val = counts_[j];
0139     } else if (varType_[j] == 1) {
0140       val = ngoodVertices;
0141     } else if (varType_[j] == 2) {
0142       val = sumPt_[j];
0143     } else if (varType_[j] == 3) {
0144       val = nVertices;
0145     }
0146
0147     corx -= formula_x_[j]->Eval(val);
0148     cory -= formula_y_[j]->Eval(val);
0149
0150   }  //end loop over corrections
0151
0152   metCorr->mex = corx;
0153   metCorr->mey = cory;
0154   evt.put(std::move(metCorr), "");
0155 }
0156
0157 #include "FWCore/Framework/interface/MakerMacros.h"
0158
0159 DEFINE_FWK_MODULE(MultShiftMETcorrInputProducer);