CandAlgos/interface/GenParticleCustomSelector.h

0001 #ifndef RecoSelectors_GenParticleCustomSelector_h
0002 #define RecoSelectors_GenParticleCustomSelector_h
0003 /* \class GenParticleCustomSelector

0004  *

0005  * \author Giuseppe Cerati, UCSD

0006  *

0007  */
0008
0009 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0011
0012 class GenParticleCustomSelector {
0013 public:
0014   GenParticleCustomSelector() {}
0015   GenParticleCustomSelector(double ptMin,
0016                             double minRapidity,
0017                             double maxRapidity,
0018                             double tip,
0019                             double lip,
0020                             bool chargedOnly,
0021                             int status,
0022                             const std::vector<int>& pdgId = std::vector<int>(),
0023                             bool invertRapidityCut = false,
0024                             double minPhi = -3.2,
0025                             double maxPhi = 3.2)
0026       : ptMin_(ptMin),
0027         minRapidity_(minRapidity),
0028         maxRapidity_(maxRapidity),
0029         meanPhi_((minPhi + maxPhi) / 2.),
0030         rangePhi_((maxPhi - minPhi) / 2.),
0031         tip_(tip),
0032         lip_(lip),
0033         chargedOnly_(chargedOnly),
0034         status_(status),
0035         pdgId_(pdgId),
0036         invertRapidityCut_(invertRapidityCut) {
0037     if (minPhi >= maxPhi) {
0038       throw cms::Exception("Configuration")
0039           << "GenParticleCustomSelector: minPhi (" << minPhi << ") must be smaller than maxPhi (" << maxPhi
0040           << "). The range is constructed from minPhi to maxPhi around their "
0041              "average.";
0042     }
0043     if (minPhi >= M_PI) {
0044       throw cms::Exception("Configuration") << "GenParticleCustomSelector: minPhi (" << minPhi
0045                                             << ") must be smaller than PI. The range is constructed from minPhi "
0046                                                "to maxPhi around their average.";
0047     }
0048     if (maxPhi <= -M_PI) {
0049       throw cms::Exception("Configuration") << "GenParticleCustomSelector: maxPhi (" << maxPhi
0050                                             << ") must be larger than -PI. The range is constructed from minPhi "
0051                                                "to maxPhi around their average.";
0052     }
0053   }
0054
0055   /// Operator() performs the selection: e.g. if (tPSelector(tp)) {...}

0056   bool operator()(const reco::GenParticle& tp) const {
0057     if (chargedOnly_ && tp.charge() == 0)
0058       return false;  //select only if charge!=0

0059     bool testId = false;
0060     unsigned int idSize = pdgId_.size();
0061     if (idSize == 0)
0062       testId = true;
0063     else
0064       for (unsigned int it = 0; it != idSize; ++it) {
0065         if (tp.pdgId() == pdgId_[it])
0066           testId = true;
0067       }
0068
0069     auto etaOk = [&](const reco::GenParticle& p) -> bool {
0070       float eta = p.eta();
0071       if (!invertRapidityCut_)
0072         return (eta >= minRapidity_) && (eta <= maxRapidity_);
0073       else
0074         return (eta < minRapidity_ || eta > maxRapidity_);
0075     };
0076     auto phiOk = [&](const reco::GenParticle& p) {
0077       float dphi = deltaPhi(atan2f(p.py(), p.px()), meanPhi_);
0078       return dphi >= -rangePhi_ && dphi <= rangePhi_;
0079     };
0080     auto ptOk = [&](const reco::GenParticle& p) {
0081       double pt = p.pt();
0082       return pt >= ptMin_;
0083     };
0084
0085     return (ptOk(tp) && etaOk(tp) && phiOk(tp) && sqrt(tp.vertex().perp2()) <= tip_ && fabs(tp.vertex().z()) <= lip_ &&
0086             tp.status() == status_ && testId);
0087   }
0088
0089 private:
0090   double ptMin_;
0091   double minRapidity_;
0092   double maxRapidity_;
0093   float meanPhi_;
0094   float rangePhi_;
0095   double tip_;
0096   double lip_;
0097   bool chargedOnly_;
0098   int status_;
0099   std::vector<int> pdgId_;
0100   bool invertRapidityCut_;
0101 };
0102
0103 #include "FWCore/Framework/interface/ConsumesCollector.h"
0104 #include "CommonTools/UtilAlgos/interface/ParameterAdapter.h"
0105
0106 namespace reco {
0107   namespace modules {
0108
0109     template <>
0110     struct ParameterAdapter<GenParticleCustomSelector> {
0111       static GenParticleCustomSelector make(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC) {
0112         return make(cfg);
0113       }
0114
0115       static GenParticleCustomSelector make(const edm::ParameterSet& cfg) {
0116         return GenParticleCustomSelector(cfg.getParameter<double>("ptMin"),
0117                                          cfg.getParameter<double>("minRapidity"),
0118                                          cfg.getParameter<double>("maxRapidity"),
0119                                          cfg.getParameter<double>("tip"),
0120                                          cfg.getParameter<double>("lip"),
0121                                          cfg.getParameter<bool>("chargedOnly"),
0122                                          cfg.getParameter<int>("status"),
0123                                          cfg.getParameter<std::vector<int> >("pdgId"),
0124                                          cfg.getParameter<bool>("invertRapidityCut"),
0125                                          cfg.getParameter<double>("minPhi"),
0126                                          cfg.getParameter<double>("maxPhi"));
0127       }
0128
0129       static void fillPSetDescription(edm::ParameterSetDescription& desc) {
0130         desc.add<double>("ptMin", 0.9);
0131         desc.add<double>("minRapidity", -2.4);
0132         desc.add<double>("maxRapidity", 2.4);
0133         desc.add<double>("tip", 3.5);
0134         desc.add<double>("lip", 30.0);
0135         desc.add<bool>("chargedOnly", true);
0136         desc.add<int>("status", 1);
0137         desc.add<std::vector<int> >("pdgId", {});
0138         desc.add<bool>("invertRapidityCut", false);
0139         desc.add<double>("minPhi", -3.2);
0140         desc.add<double>("maxPhi", 3.2);
0141       }
0142     };
0143
0144   }  // namespace modules

0145 }  // namespace reco

0146
0147 #endif