Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:28:41 

0001 // author: Florent Lacroix (UIC)
0002 // date: 07/14/2009
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   // parameters for the cut:
0033   // sqrt(DeltaMEX**2+DeltaMEY**2)>DeltaMEXsigma*sigma
0034   // with sigma=sigma_a+sigma_b*sqrt(SET)+sigma_c*SET
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   // std::cout << "FL: beginJob" << std::endl;
0099 }
0100 
0101 bool PFMETFilter::filter(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0102   // std::cout << "FL: filter" << std::endl;
0103   // std::cout << "FL: Mins = " << min_ << std::endl;
0104 
0105   if (!checkInput())
0106     return true;  // no filtering !
0107 
0108   bool skip = false;
0109 
0110   for (unsigned int varc = 0; varc < collections_.size(); ++varc) {
0111     // std::cout << "FL: var[" << varc << "] = " << collections_[varc] <<
0112     // std::endl; std::cout << "FL: var[0] = " << collections_[0] << std::endl;
0113 
0114     // if the collection is collection1-collection2:
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       // std::cout << "collection1 = " << collection1 << std::endl;
0120       std::string collection2;
0121       collection2.assign(collections_[varc], minuspos + 1, collections_[varc].size());
0122       // std::cout << "collection2 = " << collection2 << std::endl;
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       // std::cout << "FL: var10.pt = " << var10->et() << std::endl;
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       // std::cout << "FL: coll_var1[" << variables_[varc] << "] = " <<
0146       // coll_var1 << std::endl;
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       // std::cout << "FL: var20.pt = " << var20->et() << std::endl;
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       // std::cout << "FL: coll_var2[" << variables_[varc] << "] = " <<
0170       // coll_var2 << std::endl; std::cout << "FL: max_[varc] = " << max_[varc]
0171       // << std::endl; std::cout << "FL: min_[varc] = " << min_[varc] <<
0172       // std::endl;
0173 
0174       // Delta computation
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       // cuts
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       // std::cout << "skip = " << skip << std::endl;
0210     } else {
0211       // get the variable:
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       // std::cout << "FL: var00.pt = " << var00->et() << std::endl;
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         // std::cout << "FL: SETc = " << SETc << std::endl;
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       // std::cout << "FL: coll_var[" << variables_[varc] << "] = " << coll_var
0280       // << std::endl; std::cout << "FL: max_[varc] = " << max_[varc] <<
0281       // std::endl; std::cout << "FL: min_[varc] = " << min_[varc] << std::endl;
0282 
0283       // cuts
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       // std::cout << "skip = " << skip << std::endl;
0294     }
0295   }
0296   // std::cout << "final skip = " << skip << std::endl;
0297   return !skip;
0298 }
0299 
0300 void PFMETFilter::endJob() {
0301   // std::cout << "FL: endJob" << std::endl;
0302 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!