IsolationAlgos/src/IsoDepositVetoFactory.cc

0001 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositVetoFactory.h"
0002
0003 #include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
0004 #include "PhysicsTools/IsolationAlgos/interface/EventDependentAbsVetos.h"
0005 #include <regex>
0006
0007 // ---------- FIRST DEFINE NEW VETOS ------------
0008 namespace reco {
0009   namespace isodeposit {
0010
0011     class SwitchingEcalVeto : public AbsVeto {
0012     public:
0013       // creates SwitchingEcalVeto from another AbsVeto (which becomes owned by this veto)
0014       SwitchingEcalVeto(AbsVeto *veto, bool isBarrel) : veto_(veto), barrel_(isBarrel) {}
0015       bool veto(double eta, double phi, float value) const override {
0016         return (fabs(eta) < 1.479) == (barrel_) ? veto_->veto(eta, phi, value) : false;
0017       }
0018       void centerOn(double eta, double phi) override { veto_->centerOn(eta, phi); }
0019
0020     private:
0021       std::unique_ptr<AbsVeto> veto_;
0022       bool barrel_;
0023     };
0024
0025     class NumCrystalVeto : public AbsVeto {
0026     public:
0027       NumCrystalVeto(Direction dir, double iR) : vetoDir_(dir), iR_(iR) {}
0028       bool veto(double eta, double phi, float value) const override {
0029         if (fabs(vetoDir_.eta()) < 1.479) {
0030           return (vetoDir_.deltaR(Direction(eta, phi)) < 0.0174 * iR_);
0031         } else {
0032           return (vetoDir_.deltaR(Direction(eta, phi)) < 0.00864 * fabs(sinh(eta)) * iR_);
0033         }
0034       }
0035       void centerOn(double eta, double phi) override { vetoDir_ = Direction(eta, phi); }
0036
0037     private:
0038       Direction vetoDir_;
0039       float iR_;
0040     };
0041
0042     class NumCrystalEtaPhiVeto : public AbsVeto {
0043     public:
0044       NumCrystalEtaPhiVeto(const math::XYZVectorD &dir, double iEta, double iPhi)
0045           : vetoDir_(dir.eta(), dir.phi()), iEta_(iEta), iPhi_(iPhi) {}
0046       NumCrystalEtaPhiVeto(Direction dir, double iEta, double iPhi)
0047           : vetoDir_(dir.eta(), dir.phi()), iEta_(iEta), iPhi_(iPhi) {}
0048       bool veto(double eta, double phi, float value) const override {
0049         double dPhi = phi - vetoDir_.phi();
0050         double dEta = eta - vetoDir_.eta();
0051         while (dPhi < -M_PI)
0052           dPhi += 2 * M_PI;
0053         while (dPhi >= M_PI)
0054           dPhi -= 2 * M_PI;
0055         if (fabs(vetoDir_.eta()) < 1.479) {
0056           return ((fabs(dEta) < 0.0174 * iEta_) && (fabs(dPhi) < 0.0174 * iPhi_));
0057         } else {
0058           return ((fabs(dEta) < 0.00864 * fabs(sinh(eta)) * iEta_) && (fabs(dPhi) < 0.00864 * fabs(sinh(eta)) * iPhi_));
0059         }
0060       }
0061       void centerOn(double eta, double phi) override { vetoDir_ = Direction(eta, phi); }
0062
0063     private:
0064       Direction vetoDir_;
0065       double iEta_, iPhi_;
0066     };
0067
0068   }  // namespace isodeposit
0069 }  // namespace reco
0070
0071 // ---------- THEN THE ACTUAL FACTORY CODE ------------
0072 reco::isodeposit::AbsVeto *IsoDepositVetoFactory::make(const char *string, edm::ConsumesCollector &iC) {
0073   reco::isodeposit::EventDependentAbsVeto *evdep = nullptr;
0074   std::unique_ptr<reco::isodeposit::AbsVeto> ret(make(string, evdep, iC));
0075   if (evdep != nullptr) {
0076     throw cms::Exception("Configuration") << "The resulting AbsVeto depends on the edm::Event.\n"
0077                                           << "Please use the two-arguments IsoDepositVetoFactory::make.\n";
0078   }
0079   return ret.release();
0080 }
0081
0082 reco::isodeposit::AbsVeto *IsoDepositVetoFactory::make(const char *string,
0083                                                        reco::isodeposit::EventDependentAbsVeto *&evdep,
0084                                                        edm::ConsumesCollector &iC) {
0085   using namespace reco::isodeposit;
0086   static const std::regex ecalSwitch("^Ecal(Barrel|Endcaps):(.*)"), threshold("Threshold\\((\\d+\\.\\d+)\\)"),
0087       thresholdtransverse("ThresholdFromTransverse\\((\\d+\\.\\d+)\\)"),
0088       absthreshold("AbsThreshold\\((\\d+\\.\\d+)\\)"),
0089       absthresholdtransverse("AbsThresholdFromTransverse\\((\\d+\\.\\d+)\\)"), cone("ConeVeto\\((\\d+\\.\\d+)\\)"),
0090       angleCone("AngleCone\\((\\d+\\.\\d+)\\)"), angleVeto("AngleVeto\\((\\d+\\.\\d+)\\)"),
0091       rectangularEtaPhiVeto(
0092           "RectangularEtaPhiVeto\\(([+-]?\\d+\\.\\d+),([+-]?\\d+\\.\\d+),([+-]?\\d+\\.\\d+),([+-]?\\d+\\.\\d+)\\)"),
0093       numCrystal("NumCrystalVeto\\((\\d+\\.\\d+)\\)"),
0094       numCrystalEtaPhi("NumCrystalEtaPhiVeto\\((\\d+\\.\\d+),(\\d+\\.\\d+)\\)"),
0095       otherCandidatesDR("OtherCandidatesByDR\\((\\w+:?\\w*:?\\w*),\\s*(\\d+\\.?|\\d*\\.\\d*)\\)"),
0096       otherJetConstituentsDR(
0097           "OtherJetConstituentsDeltaRVeto\\((\\w+:?\\w*:?\\w*),\\s*(\\d+\\.?|\\d*\\.\\d*),\\s*(\\w+:?\\w*:?\\w*),\\s*("
0098           "\\d+\\.?|\\d*\\.\\d*)\\)"),
0099       otherCand("^(.*?):(.*)"), number("^(\\d+\\.?|\\d*\\.\\d*)$");
0100   std::cmatch match;
0101
0102   evdep = nullptr;  // by default it does not depend on this
0103   if (regex_match(string, match, ecalSwitch)) {
0104     return new SwitchingEcalVeto(make(match[2].first, iC), (match[1] == "Barrel"));
0105   } else if (regex_match(string, match, threshold)) {
0106     return new ThresholdVeto(atof(match[1].first));
0107   } else if (regex_match(string, match, thresholdtransverse)) {
0108     return new ThresholdVetoFromTransverse(atof(((std::string)match[1]).c_str()));
0109   } else if (regex_match(string, match, absthreshold)) {
0110     return new AbsThresholdVeto(atof(match[1].first));
0111   } else if (regex_match(string, match, absthresholdtransverse)) {
0112     return new AbsThresholdVetoFromTransverse(atof(((std::string)match[1]).c_str()));
0113   } else if (regex_match(string, match, cone)) {
0114     return new ConeVeto(Direction(), atof(match[1].first));
0115   } else if (regex_match(string, match, number)) {
0116     return new ConeVeto(Direction(), atof(match[1].first));
0117   } else if (regex_match(string, match, angleCone)) {
0118     return new AngleCone(Direction(), atof(match[1].first));
0119   } else if (regex_match(string, match, angleVeto)) {
0120     return new AngleConeVeto(Direction(), atof(match[1].first));
0121   } else if (regex_match(string, match, rectangularEtaPhiVeto)) {
0122     return new RectangularEtaPhiVeto(
0123         Direction(), atof(match[1].first), atof(match[2].first), atof(match[3].first), atof(match[4].first));
0124   } else if (regex_match(string, match, numCrystal)) {
0125     return new NumCrystalVeto(Direction(), atof(match[1].first));
0126   } else if (regex_match(string, match, numCrystalEtaPhi)) {
0127     return new NumCrystalEtaPhiVeto(Direction(), atof(match[1].first), atof(match[2].first));
0128   } else if (regex_match(string, match, otherCandidatesDR)) {
0129     OtherCandidatesDeltaRVeto *ret = new OtherCandidatesDeltaRVeto(edm::InputTag(match[1]), atof(match[2].first), iC);
0130     evdep = ret;
0131     return ret;
0132   } else if (regex_match(string, match, otherJetConstituentsDR)) {
0133     OtherJetConstituentsDeltaRVeto *ret = new OtherJetConstituentsDeltaRVeto(
0134         Direction(), edm::InputTag(match[1]), atof(match[2].first), edm::InputTag(match[3]), atof(match[4].first), iC);
0135     evdep = ret;
0136     return ret;
0137   } else if (regex_match(string, match, otherCand)) {
0138     OtherCandVeto *ret = new OtherCandVeto(edm::InputTag(match[1]), make(match[2].first, iC), iC);
0139     evdep = ret;
0140     return ret;
0141   } else {
0142     throw cms::Exception("Not Implemented") << "Veto " << string << " not implemented yet...";
0143   }
0144 }