File indexing completed on 2025-02-05 23:51:09
0001 #ifndef RecoSelectors_CosmicTrackingParticleSelector_h
0002 #define RecoSelectors_CosmicTrackingParticleSelector_h
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/ESHandle.h"
0016 #include "FWCore/Framework/interface/ConsumesCollector.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0019 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0020 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0021 #include "MagneticField/Engine/interface/MagneticField.h"
0022 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0023
0024 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0025 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0026 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0027 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0028
0029 #include "DataFormats/GeometrySurface/interface/Surface.h"
0030 #include "DataFormats/GeometrySurface/interface/GloballyPositioned.h"
0031 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0032 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0033 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0034
0035 #include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
0036 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0037
0038 class TrajectoryStateClosestToBeamLineBuilder;
0039
0040 class CosmicTrackingParticleSelector {
0041 public:
0042 typedef TrackingParticleCollection collection;
0043 typedef std::vector<const TrackingParticle*> container;
0044 typedef container::const_iterator const_iterator;
0045
0046 CosmicTrackingParticleSelector() {}
0047
0048 CosmicTrackingParticleSelector(double ptMin,
0049 double minRapidity,
0050 double maxRapidity,
0051 double tip,
0052 double lip,
0053 int minHit,
0054 bool chargedOnly,
0055 const std::vector<int>& pdgId = std::vector<int>())
0056 : ptMin_(ptMin),
0057 minRapidity_(minRapidity),
0058 maxRapidity_(maxRapidity),
0059 tip_(tip),
0060 lip_(lip),
0061 minHit_(minHit),
0062 chargedOnly_(chargedOnly),
0063 pdgId_(pdgId) {}
0064
0065 CosmicTrackingParticleSelector(const edm::ParameterSet& cfg, edm::ConsumesCollector&& iC)
0066 : ptMin_(cfg.getParameter<double>("ptMin")),
0067 minRapidity_(cfg.getParameter<double>("minRapidity")),
0068 maxRapidity_(cfg.getParameter<double>("maxRapidity")),
0069 tip_(cfg.getParameter<double>("tip")),
0070 lip_(cfg.getParameter<double>("lip")),
0071 minHit_(cfg.getParameter<int>("minHit")),
0072 chargedOnly_(cfg.getParameter<bool>("chargedOnly")),
0073 pdgId_(cfg.getParameter<std::vector<int> >("pdgId")),
0074 beamSpotToken_(iC.consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"))),
0075 globalTrackingGeomToken_(iC.esConsumes()),
0076 theMFToken_(iC.esConsumes()) {}
0077
0078 static void fillPSetDescription(edm::ParameterSetDescription& desc) {
0079 desc.add<double>("ptMin", 0.9);
0080 desc.add<double>("minRapidity", -2.4);
0081 desc.add<double>("maxRapidity", 2.4);
0082 desc.add<double>("tip", 100.0);
0083 desc.add<double>("lip", 100.0);
0084 desc.add<int>("minHit", 0);
0085 desc.add<bool>("chargedOnly", true);
0086 desc.add<std::vector<int> >("pdgId", {});
0087 }
0088
0089 void select(const edm::Handle<collection>& c, const edm::Event& event, const edm::EventSetup& setup) {
0090 selected_.clear();
0091 edm::Handle<reco::BeamSpot> beamSpot;
0092 event.getByToken(beamSpotToken_, beamSpot);
0093 for (TrackingParticleCollection::const_iterator itp = c->begin(); itp != c->end(); ++itp)
0094 if (operator()(TrackingParticleRef(c, itp - c->begin()), beamSpot.product(), event, setup)) {
0095 selected_.push_back(&*itp);
0096 }
0097 }
0098
0099 const_iterator begin() const { return selected_.begin(); }
0100 const_iterator end() const { return selected_.end(); }
0101
0102 void initEvent(edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssocToSet) const {
0103 simHitsTPAssoc = simHitsTPAssocToSet;
0104 }
0105
0106
0107 bool operator()(const TrackingParticleRef tpr,
0108 const reco::BeamSpot* bs,
0109 const edm::Event& iEvent,
0110 const edm::EventSetup& iSetup) const {
0111 if (chargedOnly_ && tpr->charge() == 0)
0112 return false;
0113
0114
0115
0116
0117
0118
0119
0120 edm::ESHandle<GlobalTrackingGeometry> theGeometry = iSetup.getHandle(globalTrackingGeomToken_);
0121
0122 GlobalVector finalGV(0, 0, 0);
0123 GlobalPoint finalGP(0, 0, 0);
0124 GlobalVector momentum(0, 0, 0);
0125 GlobalPoint vertex(0, 0, 0);
0126 double radius(9999);
0127 bool found(false);
0128
0129 int ii = 0;
0130 DetId::Detector det;
0131 int subdet;
0132
0133 edm::LogVerbatim("CosmicTrackingParticleSelector")
0134 << "TOT Number of PSimHits = " << tpr->numberOfHits()
0135 << ", Number of Tracker PSimHits = " << tpr->numberOfTrackerHits() << "\n";
0136
0137 if (simHitsTPAssoc.isValid() == 0) {
0138 edm::LogError("CosmicTrackingParticleSelector") << "Invalid handle!";
0139 return false;
0140 }
0141 std::pair<TrackingParticleRef, TrackPSimHitRef> clusterTPpairWithDummyTP(
0142 tpr, TrackPSimHitRef());
0143
0144 auto range = std::equal_range(simHitsTPAssoc->begin(),
0145 simHitsTPAssoc->end(),
0146 clusterTPpairWithDummyTP,
0147 SimHitTPAssociationProducer::simHitTPAssociationListGreater);
0148 for (auto ip = range.first; ip != range.second; ++ip) {
0149 TrackPSimHitRef it = ip->second;
0150 ++ii;
0151 const GeomDet* tmpDet = theGeometry->idToDet(DetId(it->detUnitId()));
0152 if (!tmpDet) {
0153 edm::LogVerbatim("CosmicTrackingParticleSelector")
0154 << "***WARNING: PSimHit " << ii << ", no GeomDet for: " << it->detUnitId() << ". Skipping it.";
0155 continue;
0156 } else {
0157 det = DetId(it->detUnitId()).det();
0158 subdet = DetId(it->detUnitId()).subdetId();
0159 }
0160
0161 LocalVector lv = it->momentumAtEntry();
0162 Local3DPoint lp = it->localPosition();
0163 GlobalVector gv = tmpDet->surface().toGlobal(lv);
0164 GlobalPoint gp = tmpDet->surface().toGlobal(lp);
0165 edm::LogVerbatim("CosmicTrackingParticleSelector")
0166 << "PSimHit " << ii << ", Detector = " << det << ", subdet = " << subdet << "\t Radius = " << gp.perp()
0167 << ", z = " << gp.z() << "\t pt = " << gv.perp() << ", pz = " << gv.z();
0168 edm::LogVerbatim("CosmicTrackingParticleSelector")
0169 << "\t trackId = " << it->trackId() << ", particleType = " << it->particleType()
0170 << ", processType = " << it->processType();
0171
0172
0173 if (it->processType() != 0)
0174 continue;
0175
0176 if (gp.perp() < radius) {
0177 found = true;
0178 radius = gp.perp();
0179 finalGV = gv;
0180 finalGP = gp;
0181 }
0182 }
0183 edm::LogVerbatim("CosmicTrackingParticleSelector")
0184 << "\n"
0185 << "FINAL State at InnerMost Hit: Radius = " << finalGP.perp() << ", z = " << finalGP.z()
0186 << ", pt = " << finalGV.perp() << ", pz = " << finalGV.z();
0187
0188 if (!found)
0189 return false;
0190 else {
0191 FreeTrajectoryState ftsAtProduction(finalGP, finalGV, TrackCharge(tpr->charge()), &iSetup.getData(theMFToken_));
0192 TSCBLBuilderNoMaterial tscblBuilder;
0193
0194 TrajectoryStateClosestToBeamLine tsAtClosestApproach = tscblBuilder(ftsAtProduction, *bs);
0195 if (!tsAtClosestApproach.isValid()) {
0196 edm::LogVerbatim("CosmicTrackingParticleSelector")
0197 << "*** WARNING in CosmicTrackingParticleSelector: tsAtClosestApproach is not valid."
0198 << "\n";
0199 return false;
0200 } else {
0201 momentum = tsAtClosestApproach.trackStateAtPCA().momentum();
0202 vertex = tsAtClosestApproach.trackStateAtPCA().position();
0203
0204 edm::LogVerbatim("CosmicTrackingParticleSelector")
0205 << "FINAL State extrapolated at PCA: Radius = " << vertex.perp() << ", z = " << vertex.z()
0206 << ", pt = " << momentum.perp() << ", pz = " << momentum.z() << "\n";
0207
0208 return (tpr->numberOfTrackerLayers() >= minHit_ && sqrt(momentum.perp2()) >= ptMin_ &&
0209 momentum.eta() >= minRapidity_ && momentum.eta() <= maxRapidity_ && sqrt(vertex.perp2()) <= tip_ &&
0210 fabs(vertex.z()) <= lip_);
0211 }
0212 }
0213 }
0214
0215 size_t size() const { return selected_.size(); }
0216
0217 private:
0218 double ptMin_;
0219 double minRapidity_;
0220 double maxRapidity_;
0221 double tip_;
0222 double lip_;
0223 int minHit_;
0224 bool chargedOnly_;
0225 std::vector<int> pdgId_;
0226 container selected_;
0227 edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0228 edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> globalTrackingGeomToken_;
0229 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theMFToken_;
0230
0231 mutable edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
0232 };
0233
0234 #endif