Common/interface/TrackingParticleSelector.h

0001 #ifndef SimTracker_Common_TrackingParticleSelector_h
0002 #define SimTracker_Common_TrackingParticleSelector_h
0003 /* \class TrackingParticleSelector
0004  *
0005  * \author Giuseppe Cerati, INFN
0006  *
0007  *  $Date: 2013/05/14 15:46:46 $
0008  *  $Revision: 1.5.4.2 $
0009  *
0010  */
0011 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0012 #include "DataFormats/Math/interface/PtEtaPhiMass.h"
0013 #include "DataFormats/Math/interface/deltaPhi.h"
0014 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0015
0016 class TrackingParticleSelector {
0017 public:
0018   TrackingParticleSelector() {}
0019   TrackingParticleSelector(double ptMin,
0020                            double ptMax,
0021                            double minRapidity,
0022                            double maxRapidity,
0023                            double tip,
0024                            double lip,
0025                            int minHit,
0026                            bool signalOnly,
0027                            bool intimeOnly,
0028                            bool chargedOnly,
0029                            bool stableOnly,
0030                            const std::vector<int> &pdgId = std::vector<int>(),
0031                            bool invertRapidityCut = false,
0032                            double minPhi = -3.2,
0033                            double maxPhi = 3.2)
0034       : ptMin2_(ptMin * ptMin),
0035         ptMax2_(ptMax * ptMax),
0036         minRapidity_(minRapidity),
0037         maxRapidity_(maxRapidity),
0038         meanPhi_((minPhi + maxPhi) / 2.),
0039         rangePhi_((maxPhi - minPhi) / 2.),
0040         tip2_(tip * tip),
0041         lip_(lip),
0042         minHit_(minHit),
0043         signalOnly_(signalOnly),
0044         intimeOnly_(intimeOnly),
0045         chargedOnly_(chargedOnly),
0046         stableOnly_(stableOnly),
0047         pdgId_(pdgId),
0048         invertRapidityCut_(invertRapidityCut) {
0049     if (minPhi >= maxPhi) {
0050       throw cms::Exception("Configuration")
0051           << "TrackingParticleSelector: minPhi (" << minPhi << ") must be smaller than maxPhi (" << maxPhi
0052           << "). The range is constructed from minPhi to maxPhi around their "
0053              "average.";
0054     }
0055     if (minPhi >= M_PI) {
0056       throw cms::Exception("Configuration") << "TrackingParticleSelector: minPhi (" << minPhi
0057                                             << ") must be smaller than PI. The range is constructed from minPhi "
0058                                                "to maxPhi around their average.";
0059     }
0060     if (maxPhi <= -M_PI) {
0061       throw cms::Exception("Configuration") << "TrackingParticleSelector: maxPhi (" << maxPhi
0062                                             << ") must be larger than -PI. The range is constructed from minPhi "
0063                                                "to maxPhi around their average.";
0064     }
0065   }
0066
0067   bool isCharged(const TrackingParticle *tp) const { return (tp->charge() == 0 ? false : true); }
0068
0069   bool isInTime(const TrackingParticle *tp) const { return (tp->eventId().bunchCrossing() == 0); }
0070
0071   bool isSignal(const TrackingParticle *tp) const {
0072     return (tp->eventId().bunchCrossing() == 0 && tp->eventId().event() == 0);
0073   }
0074
0075   bool isStable(const TrackingParticle *tp) const {
0076     for (TrackingParticle::genp_iterator j = tp->genParticle_begin(); j != tp->genParticle_end(); ++j) {
0077       if (j->get() == nullptr || j->get()->status() != 1) {
0078         return false;
0079       }
0080     }
0081     // test for remaining unstabled due to lack of genparticle pointer
0082     int pdgid = std::abs(tp->pdgId());
0083     if (tp->status() == -99 && (pdgid != 11 && pdgid != 13 && pdgid != 211 && pdgid != 321 && pdgid != 2212 &&
0084                                 pdgid != 3112 && pdgid != 3222 && pdgid != 3312 && pdgid != 3334)) {
0085       return false;
0086     }
0087     return true;
0088   }
0089
0090   /// Operator() performs the selection: e.g. if (tPSelector(tp)) {...}
0091   /// https://stackoverflow.com/questions/14466620/c-template-specialization-calling-methods-on-types-that-could-be-pointers-or/14466705
0092   bool operator()(const TrackingParticle &tp) const { return select(&tp); }
0093   bool operator()(const TrackingParticle *tp) const { return select(tp); }
0094
0095   bool select(const TrackingParticle *tp) const {
0096     // signal only means no PU particles
0097     if (signalOnly_ && !isSignal(tp))
0098       return false;
0099     // intime only means no OOT PU particles
0100     if (intimeOnly_ && !isInTime(tp))
0101       return false;
0102
0103     // select only if charge!=0
0104     if (chargedOnly_ && !isCharged(tp))
0105       return false;
0106
0107     // select for particle type
0108     if (!selectParticleType(tp)) {
0109       return false;
0110     }
0111
0112     // select only stable particles
0113     if (stableOnly_ && !isStable(tp)) {
0114       return false;
0115     }
0116
0117     return selectKinematics(tp);
0118   }
0119
0120   bool selectKinematics(const TrackingParticle *tp) const {
0121     auto etaOk = [&](const TrackingParticle *p) -> bool {
0122       float eta = etaFromXYZ(p->px(), p->py(), p->pz());
0123       if (!invertRapidityCut_)
0124         return (eta >= minRapidity_) && (eta <= maxRapidity_);
0125       else
0126         return (eta < minRapidity_ || eta > maxRapidity_);
0127     };
0128     auto phiOk = [&](const TrackingParticle *p) {
0129       float dphi = deltaPhi(atan2f(p->py(), p->px()), meanPhi_);
0130       return dphi >= -rangePhi_ && dphi <= rangePhi_;
0131     };
0132     auto ptOk = [&](const TrackingParticle *p) {
0133       double pt2 = tp->p4().perp2();
0134       return pt2 >= ptMin2_ && pt2 <= ptMax2_;
0135     };
0136     return (tp->numberOfTrackerLayers() >= minHit_ && ptOk(tp) && etaOk(tp) && phiOk(tp) &&
0137             std::abs(tp->vertex().z()) <= lip_ &&  // vertex last to avoid to load it if not striclty
0138                                                    // necessary...
0139             tp->vertex().perp2() <= tip2_);
0140   }
0141
0142   bool selectParticleType(const TrackingParticle *tp) const {
0143     auto pdgid = tp->pdgId();
0144     if (!pdgId_.empty()) {
0145       for (auto id : pdgId_) {
0146         if (id == pdgid) {
0147           return true;
0148         }
0149       }
0150     } else {
0151       return true;
0152     }
0153     return false;
0154   }
0155
0156 private:
0157   double ptMin2_;
0158   double ptMax2_;
0159   float minRapidity_;
0160   float maxRapidity_;
0161   float meanPhi_;
0162   float rangePhi_;
0163   double tip2_;
0164   double lip_;
0165   int minHit_;
0166   bool signalOnly_;
0167   bool intimeOnly_;
0168   bool chargedOnly_;
0169   bool stableOnly_;
0170   std::vector<int> pdgId_;
0171   bool invertRapidityCut_;
0172 };
0173
0174 #include "CommonTools/UtilAlgos/interface/ParameterAdapter.h"
0175 #include "FWCore/Framework/interface/ConsumesCollector.h"
0176
0177 namespace reco {
0178   namespace modules {
0179
0180     template <>
0181     struct ParameterAdapter<TrackingParticleSelector> {
0182       static TrackingParticleSelector make(const edm::ParameterSet &cfg, edm::ConsumesCollector &iC) {
0183         return make(cfg);
0184       }
0185
0186       static TrackingParticleSelector make(const edm::ParameterSet &cfg) {
0187         return TrackingParticleSelector(cfg.getParameter<double>("ptMin"),
0188                                         cfg.getParameter<double>("ptMax"),
0189                                         cfg.getParameter<double>("minRapidity"),
0190                                         cfg.getParameter<double>("maxRapidity"),
0191                                         cfg.getParameter<double>("tip"),
0192                                         cfg.getParameter<double>("lip"),
0193                                         cfg.getParameter<int>("minHit"),
0194                                         cfg.getParameter<bool>("signalOnly"),
0195                                         cfg.getParameter<bool>("intimeOnly"),
0196                                         cfg.getParameter<bool>("chargedOnly"),
0197                                         cfg.getParameter<bool>("stableOnly"),
0198                                         cfg.getParameter<std::vector<int>>("pdgId"),
0199                                         cfg.getParameter<bool>("invertRapidityCut"),
0200                                         cfg.getParameter<double>("minPhi"),
0201                                         cfg.getParameter<double>("maxPhi"));
0202       }
0203     };
0204
0205   }  // namespace modules
0206 }  // namespace reco
0207
0208 #endif