Type1MET/plugins/MultShiftMETcorrDBInputProducer.cc

0001 #include "JetMETCorrections/Type1MET/plugins/MultShiftMETcorrDBInputProducer.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "CondFormats/JetMETObjects/interface/Utilities.h"
0006
0007 #include "DataFormats/METReco/interface/CorrMETData.h"
0008 #include "DataFormats/VertexReco/interface/Vertex.h"
0009 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0010 #include "DataFormats/Common/interface/View.h"
0011
0012 #include <TString.h>
0013
0014 #include <memory>
0015
0016 int MultShiftMETcorrDBInputProducer::translateTypeToAbsPdgId(reco::PFCandidate::ParticleType type) {
0017   switch (type) {
0018     case reco::PFCandidate::ParticleType::h:
0019       return 211;  // pi+
0020     case reco::PFCandidate::ParticleType::e:
0021       return 11;
0022     case reco::PFCandidate::ParticleType::mu:
0023       return 13;
0024     case reco::PFCandidate::ParticleType::gamma:
0025       return 22;
0026     case reco::PFCandidate::ParticleType::h0:
0027       return 130;  // K_L0
0028     case reco::PFCandidate::ParticleType::h_HF:
0029       return 1;  // dummy pdg code
0030     case reco::PFCandidate::ParticleType::egamma_HF:
0031       return 2;  // dummy pdg code
0032     case reco::PFCandidate::ParticleType::X:
0033     default:
0034       return 0;
0035   }
0036 }
0037
0038 MultShiftMETcorrDBInputProducer::MultShiftMETcorrDBInputProducer(const edm::ParameterSet& cfg)
0039     : moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
0040   mPayloadName = cfg.getParameter<std::string>("payloadName");
0041   mSampleType = (cfg.exists("sampleType")) ? cfg.getParameter<std::string>("sampleType") : "MC";
0042   mIsData = cfg.getParameter<bool>("isData");
0043
0044   pflow_ = consumes<edm::View<reco::Candidate>>(cfg.getParameter<edm::InputTag>("srcPFlow"));
0045   vertices_ = consumes<edm::View<reco::Vertex>>(cfg.getParameter<edm::InputTag>("vertexCollection"));
0046   mMEtXYcorParaColl_ =
0047       esConsumes<MEtXYcorrectParametersCollection, MEtXYcorrectRecord>(edm::ESInputTag("", mPayloadName));
0048
0049   edm::InputTag srcWeights = cfg.getParameter<edm::InputTag>("srcWeights");
0050   if (!srcWeights.label().empty())
0051     weightsToken_ = consumes<edm::ValueMap<float>>(srcWeights);
0052
0053   etaMin_.clear();
0054   etaMax_.clear();
0055
0056   produces<CorrMETData>();
0057 }
0058
0059 MultShiftMETcorrDBInputProducer::~MultShiftMETcorrDBInputProducer() {}
0060
0061 void MultShiftMETcorrDBInputProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0062   // Get para.s from DB
0063   edm::ESHandle<MEtXYcorrectParametersCollection> MEtXYcorParaColl = es.getHandle(mMEtXYcorParaColl_);
0064
0065   // get the sections from Collection (pair of section and METCorr.Par class)
0066   std::vector<MEtXYcorrectParametersCollection::key_type> keys;
0067   // save level to keys for each METParameter in METParameter collection
0068   MEtXYcorParaColl->validKeys(keys);
0069
0070   //get primary vertices
0071   edm::Handle<edm::View<reco::Vertex>> hpv;
0072   evt.getByToken(vertices_, hpv);
0073   if (!hpv.isValid()) {
0074     edm::LogError("MultShiftMETcorrDBInputProducer::produce") << "could not find vertex collection ";
0075   }
0076   std::vector<reco::Vertex> goodVertices;
0077   for (unsigned i = 0; i < hpv->size(); i++) {
0078     if ((*hpv)[i].ndof() > 4 && (fabs((*hpv)[i].z()) <= 24.) && (fabs((*hpv)[i].position().rho()) <= 2.0))
0079       goodVertices.push_back((*hpv)[i]);
0080   }
0081   int ngoodVertices = goodVertices.size();
0082
0083   edm::Handle<edm::View<reco::Candidate>> particleFlow;
0084   evt.getByToken(pflow_, particleFlow);
0085
0086   edm::Handle<edm::ValueMap<float>> weights;
0087   if (!weightsToken_.isUninitialized())
0088     evt.getByToken(weightsToken_, weights);
0089
0090   //loop over all constituent types and sum each correction
0091   //std::unique_ptr<CorrMETData> metCorr(new CorrMETData());
0092   std::unique_ptr<CorrMETData> metCorr(new CorrMETData());
0093
0094   double corx = 0.;
0095   double cory = 0.;
0096
0097   // check DB
0098   for (std::vector<MEtXYcorrectParametersCollection::key_type>::const_iterator ikey = keys.begin(); ikey != keys.end();
0099        ++ikey) {
0100     if (mIsData) {
0101       if (!MEtXYcorParaColl->isShiftData(*ikey))
0102         throw cms::Exception("MultShiftMETcorrDBInputProducer::produce")
0103             << "DB is not for Data. Set proper option: \"corrPfMetXYMultDB.isData\" !!\n";
0104     } else {
0105       if (MEtXYcorParaColl->isShiftData(*ikey))
0106         throw cms::Exception("MultShiftMETcorrDBInputProducer::produce")
0107             << "DB is for Data. Set proper option: \"corrPfMetXYMultDB.isData\" !!\n";
0108     }
0109   }
0110
0111   for (std::vector<MEtXYcorrectParametersCollection::key_type>::const_iterator ikey = keys.begin(); ikey != keys.end();
0112        ++ikey) {
0113     if (!mIsData) {
0114       if (mSampleType == "MC") {
0115         if (!MEtXYcorParaColl->isShiftMC(*ikey))
0116           continue;
0117       } else if (mSampleType == "DY") {
0118         if (!MEtXYcorParaColl->isShiftDY(*ikey))
0119           continue;
0120       } else if (mSampleType == "TTJets") {
0121         if (!MEtXYcorParaColl->isShiftTTJets(*ikey))
0122           continue;
0123       } else if (mSampleType == "WJets") {
0124         if (!MEtXYcorParaColl->isShiftWJets(*ikey))
0125           continue;
0126       } else
0127         throw cms::Exception("MultShiftMETcorrDBInputProducer::produce")
0128             << "SampleType: " << mSampleType << " is not reserved !!!\n";
0129     }
0130
0131     std::string sectionName = MEtXYcorParaColl->findLabel(*ikey);
0132     MEtXYcorrectParameters const& MEtXYcorParams = (*MEtXYcorParaColl)[*ikey];
0133
0134     counts_ = 0;
0135     sumPt_ = 0;
0136
0137     for (unsigned i = 0; i < particleFlow->size(); ++i) {
0138       const reco::Candidate& c = particleFlow->at(i);
0139       if (abs(c.pdgId()) ==
0140           translateTypeToAbsPdgId(reco::PFCandidate::ParticleType(MEtXYcorParams.definitions().PtclType()))) {
0141         if ((c.eta() > MEtXYcorParams.record(0).xMin(0)) and (c.eta() < MEtXYcorParams.record(0).xMax(0))) {
0142           float weight = (!weightsToken_.isUninitialized()) ? (*weights)[particleFlow->ptrAt(i)] : 1.0;
0143           counts_ += (weight > 0);
0144           sumPt_ += c.pt() * weight;
0145           continue;
0146         }
0147       }
0148     }
0149     double val(0.);
0150     unsigned parVar = MEtXYcorParams.definitions().parVar(0);
0151
0152     if (parVar == 0) {
0153       val = counts_;
0154
0155     } else if (parVar == 1) {
0156       val = ngoodVertices;
0157
0158     } else if (parVar == 2) {
0159       val = sumPt_;
0160
0161     } else {
0162       throw cms::Exception("MultShiftMETcorrDBInputProducer::produce")
0163           << "parVar: " << parVar << " is not reserved !!!\n";
0164     }
0165
0166     formula_x_ = std::make_unique<TF1>("corrPx", MEtXYcorParams.definitions().formula().c_str());
0167     formula_y_ = std::make_unique<TF1>("corrPy", MEtXYcorParams.definitions().formula().c_str());
0168
0169     for (unsigned i(0); i < MEtXYcorParams.record(0).nParameters(); i++) {
0170       formula_x_->SetParameter(i, MEtXYcorParams.record(0).parameter(i));
0171     }
0172     for (unsigned i(0); i < MEtXYcorParams.record(1).nParameters(); i++) {
0173       formula_y_->SetParameter(i, MEtXYcorParams.record(1).parameter(i));
0174     }
0175
0176     corx -= formula_x_->Eval(val);
0177     cory -= formula_y_->Eval(val);
0178
0179   }  //end loop over corrections
0180
0181   metCorr->mex = corx;
0182   metCorr->mey = cory;
0183   evt.put(std::move(metCorr), "");
0184 }
0185
0186 #include "FWCore/Framework/interface/MakerMacros.h"
0187
0188 DEFINE_FWK_MODULE(MultShiftMETcorrDBInputProducer);