Type1MET/src/METCorrectionAlgorithm.cc

0001 #include "JetMETCorrections/Type1MET/interface/METCorrectionAlgorithm.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004
0005 #include <TString.h>
0006
0007 #include <string>
0008
0009 METCorrectionAlgorithm::METCorrectionAlgorithm(const edm::ParameterSet& cfg,
0010                                                edm::ConsumesCollector&& iConsumesCollector) {
0011   applyType1Corrections_ = cfg.getParameter<bool>("applyType1Corrections");
0012   if (applyType1Corrections_) {
0013     vInputTag srcType1Corrections = cfg.getParameter<vInputTag>("srcType1Corrections");
0014     for (vInputTag::const_iterator inputTag = srcType1Corrections.begin(); inputTag != srcType1Corrections.end();
0015          ++inputTag) {
0016       type1Tokens_.push_back(iConsumesCollector.consumes<CorrMETData>(*inputTag));
0017     }
0018   }
0019
0020   applyType2Corrections_ = cfg.getParameter<bool>("applyType2Corrections");
0021   if (applyType2Corrections_) {
0022     vInputTag srcUnclEnergySums = cfg.getParameter<vInputTag>("srcUnclEnergySums");
0023
0024     if (cfg.exists("type2Binning")) {
0025       typedef std::vector<edm::ParameterSet> vParameterSet;
0026       vParameterSet cfgType2Binning = cfg.getParameter<vParameterSet>("type2Binning");
0027       for (vParameterSet::const_iterator cfgType2BinningEntry = cfgType2Binning.begin();
0028            cfgType2BinningEntry != cfgType2Binning.end();
0029            ++cfgType2BinningEntry) {
0030         type2Binning_.push_back(
0031             new type2BinningEntryType(*cfgType2BinningEntry, srcUnclEnergySums, iConsumesCollector));
0032       }
0033     } else {
0034       std::string type2CorrFormula = cfg.getParameter<std::string>("type2CorrFormula");
0035       edm::ParameterSet type2CorrParameter = cfg.getParameter<edm::ParameterSet>("type2CorrParameter");
0036       type2Binning_.push_back(
0037           new type2BinningEntryType(type2CorrFormula, type2CorrParameter, srcUnclEnergySums, iConsumesCollector));
0038     }
0039   }
0040
0041   applyType0Corrections_ =
0042       cfg.exists("applyType0Corrections") ? cfg.getParameter<bool>("applyType0Corrections") : false;
0043   if (applyType0Corrections_) {
0044     vInputTag srcCHSSums = cfg.getParameter<vInputTag>("srcCHSSums");
0045     for (vInputTag::const_iterator inputTag = srcCHSSums.begin(); inputTag != srcCHSSums.end(); ++inputTag) {
0046       chsSumTokens_.push_back(iConsumesCollector.consumes<CorrMETData>(*inputTag));
0047     }
0048
0049     type0Rsoft_ = cfg.getParameter<double>("type0Rsoft");
0050     type0Cuncl_ = 1.0;
0051     if (applyType2Corrections_) {
0052       if (cfg.exists("type2Binning"))
0053         throw cms::Exception("Invalid Arg")
0054             << "Currently, applyType0Corrections and type2Binning cannot be used together!";
0055       std::string type2CorrFormula = cfg.getParameter<std::string>("type2CorrFormula");
0056       if (!(type2CorrFormula == "A"))
0057         throw cms::Exception("Invalid Arg") << "type2CorrFormula must be \"A\" if applyType0Corrections!";
0058       edm::ParameterSet type2CorrParameter = cfg.getParameter<edm::ParameterSet>("type2CorrParameter");
0059       type0Cuncl_ = type2CorrParameter.getParameter<double>("A");
0060     }
0061   }
0062 }
0063
0064 METCorrectionAlgorithm::~METCorrectionAlgorithm() {
0065   for (std::vector<type2BinningEntryType*>::const_iterator it = type2Binning_.begin(); it != type2Binning_.end();
0066        ++it) {
0067     delete (*it);
0068   }
0069 }
0070
0071 CorrMETData METCorrectionAlgorithm::compMETCorrection(edm::Event& evt) {
0072   CorrMETData metCorr;
0073   metCorr.mex = 0.;
0074   metCorr.mey = 0.;
0075   metCorr.sumet = 0.;
0076
0077   if (applyType0Corrections_) {
0078     //--- sum all Type 0 MET correction terms
0079     edm::Handle<CorrMETData> chsSum;
0080     for (std::vector<edm::EDGetTokenT<CorrMETData> >::const_iterator corrToken = chsSumTokens_.begin();
0081          corrToken != chsSumTokens_.end();
0082          ++corrToken) {
0083       evt.getByToken(*corrToken, chsSum);
0084
0085       metCorr.mex += type0Cuncl_ * (1 - type0Rsoft_) * chsSum->mex;
0086       metCorr.mey += type0Cuncl_ * (1 - type0Rsoft_) * chsSum->mey;
0087       metCorr.sumet += type0Cuncl_ * (1 - type0Rsoft_) * chsSum->sumet;
0088     }
0089   }
0090
0091   if (applyType1Corrections_) {
0092     //--- sum all Type 1 MET correction terms
0093     edm::Handle<CorrMETData> type1Correction;
0094     for (std::vector<edm::EDGetTokenT<CorrMETData> >::const_iterator corrToken = type1Tokens_.begin();
0095          corrToken != type1Tokens_.end();
0096          ++corrToken) {
0097       evt.getByToken(*corrToken, type1Correction);
0098
0099       metCorr.mex += type1Correction->mex;
0100       metCorr.mey += type1Correction->mey;
0101       metCorr.sumet += type1Correction->sumet;
0102     }
0103   }
0104
0105   if (applyType2Corrections_) {
0106     //--- compute momentum sum of all "unclustered energy" in the event
0107     //
0108     //    NOTE: calibration factors/formulas for Type 2 MET correction may depend on eta
0109     //         (like the jet energy correction factors do)
0110     //
0111
0112     for (std::vector<type2BinningEntryType*>::const_iterator type2BinningEntry = type2Binning_.begin();
0113          type2BinningEntry != type2Binning_.end();
0114          ++type2BinningEntry) {
0115       CorrMETData unclEnergySum;
0116       edm::Handle<CorrMETData> unclEnergySummand;
0117       for (std::vector<edm::EDGetTokenT<CorrMETData> >::const_iterator corrToken =
0118                (*type2BinningEntry)->corrTokens_.begin();
0119            corrToken != (*type2BinningEntry)->corrTokens_.end();
0120            ++corrToken) {
0121         evt.getByToken(*corrToken, unclEnergySummand);
0122
0123         unclEnergySum.mex += unclEnergySummand->mex;
0124         unclEnergySum.mey += unclEnergySummand->mey;
0125         unclEnergySum.sumet += unclEnergySummand->sumet;
0126       }
0127
0128       //--- calibrate "unclustered energy"
0129       double unclEnergySumPt = sqrt(unclEnergySum.mex * unclEnergySum.mex + unclEnergySum.mey * unclEnergySum.mey);
0130       double unclEnergyScaleFactor = (*type2BinningEntry)->binCorrFormula_->Eval(unclEnergySumPt);
0131
0132       //--- MET balances momentum of reconstructed particles,
0133       //    hence correction to "unclustered energy" and corresponding Type 2 MET correction are of opposite sign
0134       metCorr.mex -= (unclEnergyScaleFactor - 1.) * unclEnergySum.mex;
0135       metCorr.mey -= (unclEnergyScaleFactor - 1.) * unclEnergySum.mey;
0136       metCorr.sumet += (unclEnergyScaleFactor - 1.) * unclEnergySum.sumet;
0137     }
0138   }
0139
0140   return metCorr;
0141 }