File indexing completed on 2021-10-04 02:55:41
0001 #include "PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.h"
0002 #include "FWCore/Utilities/interface/isFinite.h"
0003 #include "FWCore/Utilities/interface/Likely.h"
0004
0005 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "MagneticField/Engine/interface/MagneticField.h"
0008
0009
0010
0011 namespace {
0012 inline double sqr(double a) { return a * a; }
0013 }
0014
0015 PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::PhotonConversionTrajectorySeedProducerFromSingleLegAlgo(
0016 const edm::ParameterSet& conf, edm::ConsumesCollector&& iC)
0017 : theHitsGenerator(new CombinedHitPairGeneratorForPhotonConversion(
0018 conf.getParameter<edm::ParameterSet>("OrderedHitsFactoryPSet"), iC)),
0019 theSeedCreator(new SeedForPhotonConversion1Leg(conf.getParameter<edm::ParameterSet>("SeedCreatorPSet"), iC)),
0020 theRegionProducer(
0021 new GlobalTrackingRegionProducerFromBeamSpot(conf.getParameter<edm::ParameterSet>("RegionFactoryPSet"), iC)),
0022 theClusterCheck(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet"), iC),
0023 theSilentOnClusterCheck(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet")
0024 .getUntrackedParameter<bool>("silentClusterCheck", false)),
0025 _vtxMinDoF(conf.getParameter<double>("vtxMinDoF")),
0026 _maxDZSigmas(conf.getParameter<double>("maxDZSigmas")),
0027 _maxNumSelVtx(conf.getParameter<uint32_t>("maxNumSelVtx")),
0028 _applyTkVtxConstraint(conf.getParameter<bool>("applyTkVtxConstraint")),
0029 _countSeedTracks(0),
0030 _primaryVtxInputTag(conf.getParameter<edm::InputTag>("primaryVerticesTag")),
0031 _beamSpotInputTag(conf.getParameter<edm::InputTag>("beamSpotInputTag")) {
0032 token_vertex = iC.consumes<reco::VertexCollection>(_primaryVtxInputTag);
0033 token_bs = iC.consumes<reco::BeamSpot>(_beamSpotInputTag);
0034 token_refitter = iC.consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("TrackRefitter"));
0035 token_magField = iC.esConsumes();
0036 }
0037
0038 PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::~PhotonConversionTrajectorySeedProducerFromSingleLegAlgo() {}
0039
0040 void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::find(const edm::Event& event,
0041 const edm::EventSetup& setup,
0042 TrajectorySeedCollection& output) {
0043 myEsetup = &setup;
0044 myEvent = &event;
0045
0046 assert(seedCollection == nullptr);
0047
0048 seedCollection = &output;
0049
0050 size_t clustsOrZero = theClusterCheck.tooManyClusters(event);
0051 if (clustsOrZero) {
0052 if (!theSilentOnClusterCheck)
0053 edm::LogError("TooManyClusters") << "Found too many clusters (" << clustsOrZero << "), bailing out.\n";
0054 seedCollection = nullptr;
0055 return;
0056 }
0057
0058 const auto& magField = setup.getData(token_magField);
0059 if (UNLIKELY(magField.inTesla(GlobalPoint(0., 0., 0.)).z() < 0.01)) {
0060 seedCollection = nullptr;
0061 return;
0062 }
0063
0064 _IdealHelixParameters.setMagnField(&magField);
0065
0066 event.getByToken(token_vertex, vertexHandle);
0067 if (!vertexHandle.isValid() || vertexHandle->empty()) {
0068 edm::LogError("PhotonConversionFinderFromTracks")
0069 << "Error! Can't get the product primary Vertex Collection " << _primaryVtxInputTag << "\n";
0070 seedCollection = nullptr;
0071 return;
0072 }
0073
0074 event.getByToken(token_bs, recoBeamSpotHandle);
0075
0076 regions = theRegionProducer->regions(event, setup);
0077
0078
0079 loopOnTracks();
0080
0081 #ifdef debugTSPFSLA
0082 std::stringstream ss;
0083 ss.str("");
0084 ss << "\n++++++++++++++++++\n";
0085 ss << "seed collection size " << seedCollection->size();
0086 for (auto const& tjS : *seedCollection) {
0087 po.print(ss, tjS);
0088 }
0089 edm::LogInfo("debugTrajSeedFromSingleLeg") << ss.str();
0090
0091 #endif
0092
0093
0094 theHitsGenerator->clearCache();
0095
0096 seedCollection = nullptr;
0097 }
0098
0099 void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::loopOnTracks() {
0100
0101 myEvent->getByToken(token_refitter, trackCollectionH);
0102
0103 if (trackCollectionH.isValid() == 0) {
0104 edm::LogError("MissingInput")
0105 << " could not find track collecion in PhotonConversionTrajectorySeedProducerFromSingleLegAlgo";
0106 return;
0107 }
0108
0109 size_t idx = 0, sel = 0;
0110 _countSeedTracks = 0;
0111
0112 #ifdef debugTSPFSLA
0113 ss.str("");
0114 #endif
0115
0116 for (reco::TrackCollection::const_iterator tr = trackCollectionH->begin(); tr != trackCollectionH->end();
0117 tr++, idx++) {
0118 if (rejectTrack(*tr))
0119 continue;
0120 std::vector<reco::Vertex> selectedPriVtxCompatibleWithTrack;
0121 if (!_applyTkVtxConstraint) {
0122 selectedPriVtxCompatibleWithTrack.push_back(*(vertexHandle->begin()));
0123 } else {
0124 if (!selectPriVtxCompatibleWithTrack(*tr, selectedPriVtxCompatibleWithTrack))
0125 continue;
0126 }
0127
0128 sel++;
0129 loopOnPriVtx(*tr, selectedPriVtxCompatibleWithTrack);
0130 }
0131 #ifdef debugTSPFSLA
0132 edm::LogInfo("debugTrajSeedFromSingleLeg") << ss.str();
0133 edm::LogInfo("debugTrajSeedFromSingleLeg") << "Inspected " << sel << " tracks over " << idx
0134 << " tracks. \t # tracks providing at least one seed " << _countSeedTracks;
0135 #endif
0136 }
0137
0138 bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::selectPriVtxCompatibleWithTrack(
0139 const reco::Track& tk, std::vector<reco::Vertex>& selectedPriVtxCompatibleWithTrack) {
0140 std::vector<std::pair<double, short> > idx;
0141 short count = -1;
0142
0143 double cosPhi = tk.px() / tk.pt();
0144 double sinPhi = tk.py() / tk.pt();
0145 double sphi2 = tk.covariance(2, 2);
0146 double stheta2 = tk.covariance(1, 1);
0147
0148 for (const reco::Vertex& vtx : *vertexHandle) {
0149 count++;
0150 if (vtx.ndof() <= _vtxMinDoF)
0151 continue;
0152
0153 double _dz = tk.dz(vtx.position());
0154 double _dzError = tk.dzError();
0155
0156 double cotTheta = tk.pz() / tk.pt();
0157 double dx = vtx.position().x();
0158 double dy = vtx.position().y();
0159 double sx2 = vtx.covariance(0, 0);
0160 double sy2 = vtx.covariance(1, 1);
0161
0162 double sxy2 = sqr(cosPhi * cotTheta) * sx2 + sqr(sinPhi * cotTheta) * sy2 +
0163 sqr(cotTheta * (-dx * sinPhi + dy * cosPhi)) * sphi2 +
0164 sqr((1 + cotTheta * cotTheta) * (dx * cosPhi + dy * sinPhi)) * stheta2;
0165
0166 _dzError = sqrt(
0167 _dzError * _dzError + vtx.covariance(2, 2) +
0168 sxy2);
0169
0170 #ifdef debugTSPFSLA
0171 ss << " primary vtx " << vtx.position() << " \tk vz " << tk.vz() << " vx " << tk.vx() << " vy " << tk.vy()
0172 << " pz/pt " << tk.pz() / tk.pt() << " \t dz " << _dz << " \t " << _dzError << " sxy2 " << sxy2
0173 << " \t dz/dzErr " << _dz / _dzError << std::endl;
0174 #endif
0175
0176 if (fabs(_dz) / _dzError > _maxDZSigmas)
0177 continue;
0178
0179 idx.push_back(std::pair<double, short>(fabs(_dz), count));
0180 }
0181 if (idx.empty()) {
0182 #ifdef debugTSPFSLA
0183 ss << "no vertex selected " << std::endl;
0184 #endif
0185 return false;
0186 }
0187
0188 std::stable_sort(
0189 idx.begin(), idx.end(), [](std::pair<double, short> a, std::pair<double, short> b) { return a.first < b.first; });
0190
0191 for (size_t i = 0; i < _maxNumSelVtx && i < idx.size(); ++i) {
0192 selectedPriVtxCompatibleWithTrack.push_back((*vertexHandle)[idx[i].second]);
0193 #ifdef debugTSPFSLA
0194 ss << "selected vtx dz " << idx[0].first << " position" << idx[0].second << std::endl;
0195 #endif
0196 }
0197
0198 return true;
0199 }
0200
0201 void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::loopOnPriVtx(
0202 const reco::Track& tk, const std::vector<reco::Vertex>& selectedPriVtxCompatibleWithTrack) {
0203 bool foundAtLeastASeedCand = false;
0204 for (auto const& vtx : selectedPriVtxCompatibleWithTrack) {
0205 math::XYZPoint primaryVertexPoint = math::XYZPoint(vtx.position());
0206
0207 for (IR ir = regions.begin(), irEnd = regions.end(); ir < irEnd; ++ir) {
0208 const TrackingRegion& region = **ir;
0209
0210 #ifdef debugTSPFSLA
0211 ss << "[PrintRegion] " << region.print() << std::endl;
0212 #endif
0213
0214
0215
0216
0217
0218 if (inspectTrack(&tk, region, primaryVertexPoint) and !foundAtLeastASeedCand) {
0219 foundAtLeastASeedCand = true;
0220 _countSeedTracks++;
0221 }
0222 }
0223 }
0224 }
0225
0226 bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::rejectTrack(const reco::Track& track) {
0227 math::XYZVector beamSpot;
0228 if (recoBeamSpotHandle.isValid()) {
0229 beamSpot = math::XYZVector(recoBeamSpotHandle->position());
0230
0231 _IdealHelixParameters.setData(&track, beamSpot);
0232 if (_IdealHelixParameters.GetTangentPoint().r() == 0) {
0233
0234 return true;
0235 }
0236
0237 float rMin = 2.;
0238 if (_IdealHelixParameters.GetTangentPoint().rho() < rMin) {
0239
0240
0241
0242
0243 return true;
0244 }
0245 }
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280 return false;
0281 }
0282
0283 bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::inspectTrack(const reco::Track* track,
0284 const TrackingRegion& region,
0285 math::XYZPoint& primaryVertexPoint) {
0286 _IdealHelixParameters.setData(track, primaryVertexPoint);
0287
0288 if (edm::isNotFinite(_IdealHelixParameters.GetTangentPoint().r()) ||
0289 (_IdealHelixParameters.GetTangentPoint().r() == 0)) {
0290
0291 return false;
0292 }
0293
0294 float rMin = 3.;
0295 if (_IdealHelixParameters.GetTangentPoint().rho() < rMin) {
0296
0297
0298
0299
0300 return false;
0301 }
0302
0303 float ptmin = 0.5;
0304 float originRBound = 3;
0305 float originZBound = 3.;
0306
0307 GlobalPoint originPos;
0308 originPos = GlobalPoint(_IdealHelixParameters.GetTangentPoint().x(),
0309 _IdealHelixParameters.GetTangentPoint().y(),
0310 _IdealHelixParameters.GetTangentPoint().z());
0311 float cotTheta;
0312 if (std::abs(_IdealHelixParameters.GetMomentumAtTangentPoint().rho()) > 1.e-4f) {
0313 cotTheta =
0314 _IdealHelixParameters.GetMomentumAtTangentPoint().z() / _IdealHelixParameters.GetMomentumAtTangentPoint().rho();
0315 } else {
0316 if (_IdealHelixParameters.GetMomentumAtTangentPoint().z() > 0)
0317 cotTheta = 99999.f;
0318 else
0319 cotTheta = -99999.f;
0320 }
0321 GlobalVector originBounds(originRBound, originRBound, originZBound);
0322
0323 GlobalPoint pvtxPoint(primaryVertexPoint.x(), primaryVertexPoint.y(), primaryVertexPoint.z());
0324
0325 ConversionRegion convRegion(originPos, pvtxPoint, cotTheta, track->thetaError(), -1 * track->charge());
0326
0327 #ifdef debugTSPFSLA
0328 ss << "\nConversion Point " << originPos << " " << originPos.perp() << "\n";
0329 #endif
0330
0331 const OrderedSeedingHits& hitss = theHitsGenerator->run(convRegion, region, *myEvent, *myEsetup);
0332
0333 unsigned int nHitss = hitss.size();
0334
0335 if (nHitss == 0)
0336 return false;
0337
0338 #ifdef debugTSPFSLA
0339 ss << "\n nHitss " << nHitss << "\n";
0340 #endif
0341
0342 if (seedCollection->empty())
0343 seedCollection->reserve(
0344 nHitss);
0345
0346 for (unsigned int iHits = 0; iHits < nHitss; ++iHits) {
0347 #ifdef debugTSPFSLA
0348 ss << "\n iHits " << iHits << "\n";
0349 #endif
0350 const SeedingHitSet& hits = hitss[iHits];
0351 theSeedCreator->trajectorySeed(
0352 *seedCollection, hits, originPos, originBounds, ptmin, *myEsetup, convRegion.cotTheta(), ss);
0353 }
0354 return true;
0355 }