Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-05-20 02:50:33 

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
[ 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!