Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-31 22:27:10 

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