HGCalValidation/interface/CaloParticleSelector.h

0001 #ifndef Validation_HGCalValidation_CaloParticleSelector_h
0002 #define Validation_HGCalValidation_CaloParticleSelector_h
0003
0004 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0005 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0006 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0007 #include "DataFormats/Math/interface/PtEtaPhiMass.h"
0008 #include "DataFormats/Math/interface/deltaPhi.h"
0009
0010 class CaloParticleSelector {
0011 public:
0012   CaloParticleSelector() {}
0013   CaloParticleSelector(double ptMin,
0014                        double ptMax,
0015                        double minRapidity,
0016                        double maxRapidity,
0017                        double lip,
0018                        double tip,
0019                        int minHit,
0020                        unsigned int maxSimClusters,
0021                        bool signalOnly,
0022                        bool intimeOnly,
0023                        bool chargedOnly,
0024                        bool stableOnly,
0025                        bool notConvertedOnly,
0026                        const std::vector<int>& pdgId = std::vector<int>(),
0027                        double minPhi = -3.2,
0028                        double maxPhi = 3.2)
0029       : ptMin2_(ptMin * ptMin),
0030         ptMax2_(ptMax * ptMax),
0031         minRapidity_(minRapidity),
0032         maxRapidity_(maxRapidity),
0033         lip_(lip),
0034         tip2_(tip * tip),
0035         meanPhi_((minPhi + maxPhi) / 2.),
0036         rangePhi_((maxPhi - minPhi) / 2.),
0037         minHit_(minHit),
0038         maxSimClusters_(maxSimClusters),
0039         signalOnly_(signalOnly),
0040         intimeOnly_(intimeOnly),
0041         chargedOnly_(chargedOnly),
0042         stableOnly_(stableOnly),
0043         notConvertedOnly_(notConvertedOnly),
0044         pdgId_(pdgId) {
0045     if (minPhi >= maxPhi) {
0046       throw cms::Exception("Configuration")
0047           << "CaloParticleSelector: minPhi (" << minPhi << ") must be smaller than maxPhi (" << maxPhi
0048           << "). The range is constructed from minPhi to maxPhi around their average.";
0049     }
0050     if (minPhi >= M_PI) {
0051       throw cms::Exception("Configuration")
0052           << "CaloParticleSelector: minPhi (" << minPhi
0053           << ") must be smaller than PI. The range is constructed from minPhi to maxPhi around their average.";
0054     }
0055     if (maxPhi <= -M_PI) {
0056       throw cms::Exception("Configuration")
0057           << "CaloParticleSelector: maxPhi (" << maxPhi
0058           << ") must be larger than -PI. The range is constructed from minPhi to maxPhi around their average.";
0059     }
0060   }
0061
0062   // Operator() performs the selection: e.g. if (cPSelector(cp)) {...}
0063   // For the moment there shouldn't be any SimTracks from different crossings in the CaloParticle.
0064   bool operator()(const CaloParticle& cp, std::vector<SimVertex> const& simVertices) const {
0065     // signal only means no PU particles
0066     if (signalOnly_ && !(cp.eventId().bunchCrossing() == 0 && cp.eventId().event() == 0))
0067       return false;
0068     // intime only means no OOT PU particles
0069     if (intimeOnly_ && !(cp.eventId().bunchCrossing() == 0))
0070       return false;
0071
0072     if (cp.simClusters().size() > maxSimClusters_)
0073       return false;
0074
0075     auto pdgid = cp.pdgId();
0076     if (!pdgId_.empty()) {
0077       bool testId = false;
0078       for (auto id : pdgId_) {
0079         if (id == pdgid) {
0080           testId = true;
0081           break;
0082         }
0083       }
0084       if (!testId)
0085         return false;
0086     }
0087
0088     if (chargedOnly_ && cp.charge() == 0)
0089       return false;  //select only if charge!=0
0090
0091     // select only stable particles
0092     if (stableOnly_) {
0093       for (CaloParticle::genp_iterator j = cp.genParticle_begin(); j != cp.genParticle_end(); ++j) {
0094         if (j->get() == nullptr || j->get()->status() != 1) {
0095           return false;
0096         }
0097       }
0098
0099       // test for remaining unstabled due to lack of genparticle pointer
0100       std::vector<int> pdgids{11, 13, 211, 321, 2212, 3112, 3222, 3312, 3334};
0101       if (cp.status() == -99 && (!std::binary_search(pdgids.begin(), pdgids.end(), std::abs(pdgid)))) {
0102         return false;
0103       }
0104     }
0105
0106     // select only particles which did not convert/decay before the calorimeter
0107     // in case of electrons, bremsstrahlung is usually the cause, thus this selection is skipped
0108     if (std::abs(pdgid) != 11) {
0109       if (notConvertedOnly_) {
0110         if (cp.g4Tracks()[0].getPositionAtBoundary() == math::XYZTLorentzVectorF(0, 0, 0, 0)) {
0111           return false;
0112         }
0113         if (cp.g4Tracks()[0].getMomentumAtBoundary() == math::XYZTLorentzVectorF(0, 0, 0, 0)) {
0114           return false;
0115         }
0116       }
0117     }
0118
0119     auto etaOk = [&](const CaloParticle& p) -> bool {
0120       float eta = etaFromXYZ(p.px(), p.py(), p.pz());
0121       return (eta >= minRapidity_) & (eta <= maxRapidity_);
0122     };
0123     auto phiOk = [&](const CaloParticle& p) {
0124       float dphi = deltaPhi(atan2f(p.py(), p.px()), meanPhi_);
0125       return dphi >= -rangePhi_ && dphi <= rangePhi_;
0126     };
0127     auto ptOk = [&](const CaloParticle& p) {
0128       double pt2 = cp.p4().perp2();
0129       return pt2 >= ptMin2_ && pt2 <= ptMax2_;
0130     };
0131
0132     return (ptOk(cp) && etaOk(cp) && phiOk(cp));
0133   }
0134
0135 private:
0136   double ptMin2_;
0137   double ptMax2_;
0138   float minRapidity_;
0139   float maxRapidity_;
0140   double lip_;
0141   double tip2_;
0142   float meanPhi_;
0143   float rangePhi_;
0144   int minHit_;
0145   unsigned int maxSimClusters_;
0146   bool signalOnly_;
0147   bool intimeOnly_;
0148   bool chargedOnly_;
0149   bool stableOnly_;
0150   bool notConvertedOnly_;
0151   std::vector<int> pdgId_;
0152 };
0153
0154 #include "FWCore/Framework/interface/ConsumesCollector.h"
0155 #include "CommonTools/UtilAlgos/interface/ParameterAdapter.h"
0156
0157 namespace reco {
0158   namespace modules {
0159
0160     template <>
0161     struct ParameterAdapter<CaloParticleSelector> {
0162       static CaloParticleSelector make(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC) { return make(cfg); }
0163
0164       static CaloParticleSelector make(const edm::ParameterSet& cfg) {
0165         return CaloParticleSelector(cfg.getParameter<double>("ptMinCP"),
0166                                     cfg.getParameter<double>("ptMaxCP"),
0167                                     cfg.getParameter<double>("minRapidityCP"),
0168                                     cfg.getParameter<double>("maxRapidityCP"),
0169                                     cfg.getParameter<double>("lip"),
0170                                     cfg.getParameter<double>("tip"),
0171                                     cfg.getParameter<int>("minHitCP"),
0172                                     cfg.getParameter<int>("maxSimClustersCP"),
0173                                     cfg.getParameter<bool>("signalOnlyCP"),
0174                                     cfg.getParameter<bool>("intimeOnlyCP"),
0175                                     cfg.getParameter<bool>("chargedOnlyCP"),
0176                                     cfg.getParameter<bool>("stableOnlyCP"),
0177                                     cfg.getParameter<bool>("notConvertedOnlyCP"),
0178                                     cfg.getParameter<std::vector<int> >("pdgIdCP"),
0179                                     cfg.getParameter<double>("minPhiCP"),
0180                                     cfg.getParameter<double>("maxPhiCP"));
0181       }
0182     };
0183
0184   }  // namespace modules
0185 }  // namespace reco
0186
0187 #endif