File indexing completed on 2023-03-17 11:28:41
0001
0002
0003
0004 #include "DataFormats/Candidate/interface/Candidate.h"
0005 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "DataFormats/METReco/interface/MET.h"
0008 #include "FWCore/Framework/interface/global/EDFilter.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/Utilities/interface/InputTag.h"
0014
0015 class PFMETFilter : public edm::global::EDFilter<> {
0016 public:
0017 explicit PFMETFilter(const edm::ParameterSet &);
0018 ~PFMETFilter() override;
0019
0020 bool filter(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0021 void beginJob() override;
0022 void endJob() override;
0023 bool checkInput() const;
0024
0025 private:
0026 std::vector<std::string> collections_;
0027 std::vector<std::string> variables_;
0028 std::vector<double> min_;
0029 std::vector<double> max_;
0030 std::vector<int> doMin_;
0031 std::vector<int> doMax_;
0032
0033
0034
0035 std::string TrueMET_;
0036 double DeltaMEXsigma_;
0037 double sigma_a_;
0038 double sigma_b_;
0039 double sigma_c_;
0040 bool verbose_;
0041 };
0042
0043 #include "FWCore/Framework/interface/MakerMacros.h"
0044 DEFINE_FWK_MODULE(PFMETFilter);
0045
0046 PFMETFilter::PFMETFilter(const edm::ParameterSet &iConfig) {
0047 collections_ = iConfig.getParameter<std::vector<std::string>>("Collections");
0048 variables_ = iConfig.getParameter<std::vector<std::string>>("Variables");
0049 min_ = iConfig.getParameter<std::vector<double>>("Mins");
0050 max_ = iConfig.getParameter<std::vector<double>>("Maxs");
0051 doMin_ = iConfig.getParameter<std::vector<int>>("DoMin");
0052 doMax_ = iConfig.getParameter<std::vector<int>>("DoMax");
0053 TrueMET_ = iConfig.getParameter<std::string>("TrueMET");
0054 DeltaMEXsigma_ = iConfig.getParameter<double>("DeltaMEXsigma");
0055 sigma_a_ = iConfig.getParameter<double>("sigma_a");
0056 sigma_b_ = iConfig.getParameter<double>("sigma_b");
0057 sigma_c_ = iConfig.getParameter<double>("sigma_c");
0058 verbose_ = iConfig.getParameter<bool>("verbose");
0059 }
0060
0061 PFMETFilter::~PFMETFilter() {}
0062
0063 bool PFMETFilter::checkInput() const {
0064 if (collections_.size() != min_.size()) {
0065 std::cout << "Error: in PFMETFilter: collections_.size()!=min_.size()" << std::endl;
0066 std::cout << "collections_.size() = " << collections_.size() << std::endl;
0067 std::cout << "min_.size() = " << min_.size() << std::endl;
0068 return false;
0069 }
0070 if (collections_.size() != max_.size()) {
0071 std::cout << "Error: in PFMETFilter: collections_.size()!=max_.size()" << std::endl;
0072 std::cout << "collections_.size() = " << collections_.size() << std::endl;
0073 std::cout << "max_.size() = " << max_.size() << std::endl;
0074 return false;
0075 }
0076 if (collections_.size() != doMin_.size()) {
0077 std::cout << "Error: in PFMETFilter: collections_.size()!=min_.size()" << std::endl;
0078 std::cout << "collections_.size() = " << collections_.size() << std::endl;
0079 std::cout << "doMin_.size() = " << doMin_.size() << std::endl;
0080 return false;
0081 }
0082 if (collections_.size() != doMax_.size()) {
0083 std::cout << "Error: in PFMETFilter: collections_.size()!=min_.size()" << std::endl;
0084 std::cout << "collections_.size() = " << collections_.size() << std::endl;
0085 std::cout << "doMax_.size() = " << doMax_.size() << std::endl;
0086 return false;
0087 }
0088 if (collections_.size() != variables_.size()) {
0089 std::cout << "Error: in PFMETFilter: collections_.size()!=variables_.size()" << std::endl;
0090 std::cout << "collections_.size() = " << collections_.size() << std::endl;
0091 std::cout << "variables_.size() = " << variables_.size() << std::endl;
0092 return false;
0093 }
0094 return true;
0095 }
0096
0097 void PFMETFilter::beginJob() {
0098
0099 }
0100
0101 bool PFMETFilter::filter(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0102
0103
0104
0105 if (!checkInput())
0106 return true;
0107
0108 bool skip = false;
0109
0110 for (unsigned int varc = 0; varc < collections_.size(); ++varc) {
0111
0112
0113
0114
0115 const unsigned int minuspos = collections_[varc].find('-');
0116 if (minuspos < collections_[varc].size()) {
0117 std::string collection1;
0118 collection1.assign(collections_[varc], 0, minuspos);
0119
0120 std::string collection2;
0121 collection2.assign(collections_[varc], minuspos + 1, collections_[varc].size());
0122
0123
0124 const edm::View<reco::Candidate> *var1;
0125 edm::Handle<edm::View<reco::Candidate>> var1_hnd;
0126 bool isVar1 = iEvent.getByLabel(collection1, var1_hnd);
0127 if (!isVar1) {
0128 std::cout << "Warning : no " << collection1 << " in input !" << std::endl;
0129 return false;
0130 }
0131 var1 = var1_hnd.product();
0132 const reco::Candidate *var10 = &(*var1)[0];
0133
0134 double coll_var1;
0135 if (variables_[varc] == "et")
0136 coll_var1 = var10->et();
0137 else if (variables_[varc] == "phi")
0138 coll_var1 = var10->phi();
0139 else if (variables_[varc] == "eta")
0140 coll_var1 = var10->eta();
0141 else {
0142 std::cout << "Error: PFMETFilter: variable unknown: " << variables_[varc] << std::endl;
0143 return true;
0144 }
0145
0146
0147
0148 const edm::View<reco::Candidate> *var2;
0149 edm::Handle<edm::View<reco::Candidate>> var2_hnd;
0150 bool isVar2 = iEvent.getByLabel(collection2, var2_hnd);
0151 if (!isVar2) {
0152 std::cout << "Warning : no " << collection2 << " in input !" << std::endl;
0153 return false;
0154 }
0155 var2 = var2_hnd.product();
0156 const reco::Candidate *var20 = &(*var2)[0];
0157
0158 double coll_var2;
0159 if (variables_[varc] == "et")
0160 coll_var2 = var20->et();
0161 else if (variables_[varc] == "phi")
0162 coll_var2 = var20->phi();
0163 else if (variables_[varc] == "eta")
0164 coll_var2 = var20->eta();
0165 else {
0166 std::cout << "Error: PFMETFilter: variable unknown: " << variables_[varc] << std::endl;
0167 return true;
0168 }
0169
0170
0171
0172
0173
0174
0175 double delta = coll_var1 - coll_var2;
0176 if (variables_[varc] == "phi") {
0177 if (coll_var1 > M_PI)
0178 coll_var1 -= ceil((coll_var1 - M_PI) / (2 * M_PI)) * 2 * M_PI;
0179 if (coll_var1 <= -M_PI)
0180 coll_var1 += ceil((coll_var1 + M_PI) / (-2. * M_PI)) * 2. * M_PI;
0181 if (coll_var2 > M_PI)
0182 coll_var2 -= ceil((coll_var2 - M_PI) / (2 * M_PI)) * 2 * M_PI;
0183 if (coll_var2 <= -M_PI)
0184 coll_var2 += ceil((coll_var2 + M_PI) / (-2. * M_PI)) * 2 * M_PI;
0185
0186 double deltaphi = -999.0;
0187 if (fabs(coll_var1 - coll_var2) < M_PI) {
0188 deltaphi = (coll_var1 - coll_var2);
0189 } else {
0190 if ((coll_var1 - coll_var2) > 0.0) {
0191 deltaphi = (2 * M_PI - fabs(coll_var1 - coll_var2));
0192 } else {
0193 deltaphi = -(2 * M_PI - fabs(coll_var1 - coll_var2));
0194 }
0195 }
0196 delta = deltaphi;
0197 }
0198
0199
0200 if (doMin_[varc] && doMax_[varc] && max_[varc] < min_[varc]) {
0201 if (delta > max_[varc] && delta < min_[varc])
0202 skip = true;
0203 } else {
0204 if (doMin_[varc] && delta < min_[varc])
0205 skip = true;
0206 if (doMax_[varc] && delta > max_[varc])
0207 skip = true;
0208 }
0209
0210 } else {
0211
0212 const edm::View<reco::Candidate> *var0;
0213 edm::Handle<edm::View<reco::Candidate>> var0_hnd;
0214 bool isVar0 = iEvent.getByLabel(collections_[varc], var0_hnd);
0215 if (!isVar0) {
0216 std::cout << "Warning : no " << collections_[varc] << " in input !" << std::endl;
0217 return false;
0218 }
0219 var0 = var0_hnd.product();
0220 const reco::Candidate *var00 = &(*var0)[0];
0221
0222 double coll_var;
0223 if (variables_[varc] == "et")
0224 coll_var = var00->et();
0225 else if (variables_[varc] == "phi")
0226 coll_var = var00->phi();
0227 else if (variables_[varc] == "eta")
0228 coll_var = var00->eta();
0229 else if (variables_[varc] == "DeltaMEXcut") {
0230 const edm::View<reco::Candidate> *truevar0;
0231 edm::Handle<edm::View<reco::Candidate>> truevar0_hnd;
0232 bool istrueVar0 = iEvent.getByLabel(TrueMET_, truevar0_hnd);
0233 if (!istrueVar0) {
0234 std::cout << "Warning : no " << TrueMET_ << " in input !" << std::endl;
0235 return false;
0236 }
0237 truevar0 = truevar0_hnd.product();
0238 const reco::Candidate *truevar00 = &(*truevar0)[0];
0239
0240 const double DeltaMEX = var00->px() - truevar00->px();
0241 const double DeltaMEY = var00->py() - truevar00->py();
0242 const double cutvalc = sqrt(DeltaMEX * DeltaMEX + DeltaMEY * DeltaMEY);
0243 const reco::MET *met = static_cast<const reco::MET *>(truevar00);
0244 const double SETc = met->sumEt();
0245
0246 const double sigmac = sigma_a_ + sigma_b_ * sqrt(SETc) + sigma_c_ * SETc;
0247 if (cutvalc > DeltaMEXsigma_ * sigmac) {
0248 if (verbose_) {
0249 std::cout << "DeltaMET = " << var00->et() - truevar00->et() << std::endl;
0250 std::cout << "trueSET = " << SETc << std::endl;
0251 std::cout << "pfMET = " << var00->et() << std::endl;
0252 std::cout << "trueMET = " << truevar00->et() << std::endl;
0253 std::cout << "DeltaMEX = " << DeltaMEX << std::endl;
0254 std::cout << "DeltaMEY = " << DeltaMEY << std::endl;
0255 std::cout << "cutvalc = " << cutvalc << std::endl;
0256 std::cout << "sigmac = " << sigmac << std::endl;
0257 std::cout << "cutvalc/sigmac = " << cutvalc / sigmac << std::endl;
0258 }
0259 return true;
0260 } else {
0261 if (verbose_ && (var00->et() - truevar00->et()) > 300.0) {
0262 std::cout << "EVENT NOT KEPT:" << std::endl;
0263 std::cout << "DeltaMET = " << var00->et() - truevar00->et() << std::endl;
0264 std::cout << "SETc = " << SETc << std::endl;
0265 std::cout << "pfMET = " << var00->et() << std::endl;
0266 std::cout << "trueMET = " << truevar00->et() << std::endl;
0267 std::cout << "DeltaMEX = " << DeltaMEX << std::endl;
0268 std::cout << "DeltaMEY = " << DeltaMEY << std::endl;
0269 std::cout << "cutvalc = " << cutvalc << std::endl;
0270 std::cout << "sigmac = " << sigmac << std::endl;
0271 std::cout << "cutvalc/sigmac = " << cutvalc / sigmac << std::endl;
0272 }
0273 return false;
0274 }
0275 } else {
0276 std::cout << "Error: PFMETFilter: variable unknown: " << variables_[varc] << std::endl;
0277 return true;
0278 }
0279
0280
0281
0282
0283
0284 if (doMin_[varc] && doMax_[varc] && max_[varc] < min_[varc]) {
0285 if (coll_var > max_[varc] && coll_var < min_[varc])
0286 skip = true;
0287 } else {
0288 if (doMin_[varc] && coll_var < min_[varc])
0289 skip = true;
0290 if (doMax_[varc] && coll_var > max_[varc])
0291 skip = true;
0292 }
0293
0294 }
0295 }
0296
0297 return !skip;
0298 }
0299
0300 void PFMETFilter::endJob() {
0301
0302 }