File indexing completed on 2024-10-08 23:10:02
0001 #include "RecoTracker/TkSeedGenerator/plugins/MultiHitGeneratorFromChi2.h"
0002 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0003
0004 #include "RecoTracker/PixelSeeding/interface/ThirdHitPredictionFromCircle.h"
0005 #include "RecoTracker/PixelSeeding/interface/ThirdHitRZPrediction.h"
0006 #include "FWCore/Utilities/interface/ESInputTag.h"
0007
0008 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0009 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0010 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0011 #include "RecoTracker/TkHitPairs/interface/RecHitsSortedInPhi.h"
0012
0013 #include "CommonTools/RecoAlgos/interface/KDTreeLinkerAlgo.h"
0014
0015 #include "RecoTracker/PixelTrackFitting/interface/RZLine.h"
0016 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0017 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0018 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0019 #include "RecoTracker/TkHitPairs/interface/HitPairGeneratorFromLayerPair.h"
0020
0021 #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"
0022 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0023 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0024 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0025 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0026
0027 #include "DataFormats/Math/interface/normalizedPhi.h"
0028
0029 #include "FWCore/Utilities/interface/isFinite.h"
0030 #include "CommonTools/Utils/interface/DynArray.h"
0031
0032 #include <algorithm>
0033 #include <iostream>
0034 #include <vector>
0035 #include <cmath>
0036 #include <map>
0037 #include <limits>
0038
0039 using namespace std;
0040
0041 typedef PixelRecoRange<float> Range;
0042
0043 namespace {
0044 struct LayerRZPredictions {
0045 ThirdHitRZPrediction<SimpleLineRZ> line;
0046 };
0047 }
0048
0049 MultiHitGeneratorFromChi2::MultiHitGeneratorFromChi2(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0050 : MultiHitGeneratorFromPairAndLayers(cfg),
0051 useFixedPreFiltering(cfg.getParameter<bool>("useFixedPreFiltering")),
0052 extraHitRZtolerance(
0053 cfg.getParameter<double>("extraHitRZtolerance")),
0054 extraHitRPhitolerance(cfg.getParameter<double>(
0055 "extraHitRPhitolerance")),
0056 extraZKDBox(
0057 cfg.getParameter<double>("extraZKDBox")),
0058 extraRKDBox(
0059 cfg.getParameter<double>("extraRKDBox")),
0060 extraPhiKDBox(cfg.getParameter<double>("extraPhiKDBox")),
0061 fnSigmaRZ(cfg.getParameter<double>(
0062 "fnSigmaRZ")),
0063 chi2VsPtCut(cfg.getParameter<bool>("chi2VsPtCut")),
0064 maxChi2(cfg.getParameter<double>("maxChi2")),
0065 refitHits(cfg.getParameter<bool>("refitHits")),
0066 filterName_(cfg.getParameter<std::string>("ClusterShapeHitFilterName")),
0067 builderName_(cfg.getParameter<std::string>("TTRHBuilder")),
0068 useSimpleMF_(false),
0069 mfName_("") {
0070 if (useFixedPreFiltering)
0071 dphi = cfg.getParameter<double>("phiPreFiltering");
0072 if (chi2VsPtCut) {
0073 pt_interv = cfg.getParameter<std::vector<double> >("pt_interv");
0074 chi2_cuts = cfg.getParameter<std::vector<double> >("chi2_cuts");
0075 }
0076 #ifdef EDM_ML_DEBUG
0077 detIdsToDebug = cfg.getParameter<std::vector<int> >("detIdsToDebug");
0078
0079 #else
0080 detIdsToDebug.push_back(0);
0081 detIdsToDebug.push_back(0);
0082 detIdsToDebug.push_back(0);
0083 #endif
0084
0085
0086
0087
0088 if (cfg.exists("SimpleMagneticField")) {
0089 useSimpleMF_ = true;
0090 mfName_ = cfg.getParameter<std::string>("SimpleMagneticField");
0091 }
0092 filter = nullptr;
0093 bfield = nullptr;
0094 nomField = -1.;
0095
0096 if (useSimpleMF_) {
0097 magneticFieldESToken_ = iC.esConsumes(edm::ESInputTag("", mfName_));
0098 } else {
0099 magneticFieldESToken_ = iC.esConsumes();
0100 }
0101 if (refitHits) {
0102 clusterShapeHitFilterESToken_ = iC.esConsumes(edm::ESInputTag("", filterName_));
0103 transientTrackingRecHitBuilderESToken_ = iC.esConsumes(edm::ESInputTag("", builderName_));
0104 }
0105 }
0106
0107 MultiHitGeneratorFromChi2::~MultiHitGeneratorFromChi2() {}
0108
0109 void MultiHitGeneratorFromChi2::fillDescriptions(edm::ParameterSetDescription& desc) {
0110 MultiHitGeneratorFromPairAndLayers::fillDescriptions(desc);
0111
0112
0113 desc.add<bool>("useFixedPreFiltering", false);
0114 desc.add<double>("phiPreFiltering", 0.3);
0115
0116
0117 desc.add<double>("extraHitRPhitolerance", 0);
0118 desc.add<double>("extraHitRZtolerance", 0);
0119 desc.add<double>("extraZKDBox", 0.2);
0120 desc.add<double>("extraRKDBox", 0.2);
0121 desc.add<double>("extraPhiKDBox", 0.005);
0122 desc.add<double>("fnSigmaRZ", 2.0);
0123
0124
0125 desc.add<bool>("refitHits", true);
0126 desc.add<std::string>("ClusterShapeHitFilterName", "ClusterShapeHitFilter");
0127 desc.add<std::string>("TTRHBuilder", "WithTrackAngle");
0128
0129
0130 desc.add<double>("maxChi2", 5.0);
0131 desc.add<bool>("chi2VsPtCut", true);
0132 desc.add<std::vector<double> >("pt_interv", std::vector<double>{{0.4, 0.7, 1.0, 2.0}});
0133 desc.add<std::vector<double> >("chi2_cuts", std::vector<double>{{3.0, 4.0, 5.0, 5.0}});
0134
0135
0136 desc.add<std::vector<int> >("detIdsToDebug", std::vector<int>{{0, 0, 0}});
0137 }
0138
0139 void MultiHitGeneratorFromChi2::initES(const edm::EventSetup& es) {
0140 bfield = &es.getData(magneticFieldESToken_);
0141 nomField = bfield->nominalValue();
0142 ufield.set(nomField);
0143
0144 if (refitHits) {
0145 filter = &es.getData(clusterShapeHitFilterESToken_);
0146 auto const& builderRef = es.getData(transientTrackingRecHitBuilderESToken_);
0147 builder = (TkTransientTrackingRecHitBuilder const*)(&builderRef);
0148 cloner = (*builder).cloner();
0149 }
0150 }
0151
0152 namespace {
0153 inline bool intersect(Range& range, const Range& second) {
0154 if (range.first > second.max() || range.second < second.min())
0155 return false;
0156 if (range.first < second.min())
0157 range.first = second.min();
0158 if (range.second > second.max())
0159 range.second = second.max();
0160 return range.first < range.second;
0161 }
0162 }
0163
0164 void MultiHitGeneratorFromChi2::hitSets(const TrackingRegion& region,
0165 OrderedMultiHits& result,
0166 const edm::Event& ev,
0167 const edm::EventSetup& es,
0168 SeedingLayerSetsHits::SeedingLayerSet pairLayers,
0169 std::vector<SeedingLayerSetsHits::SeedingLayer> thirdLayers) {
0170 LogDebug("MultiHitGeneratorFromChi2") << "pair: " << thePairGenerator->innerLayer(pairLayers).name() << "+"
0171 << thePairGenerator->outerLayer(pairLayers).name()
0172 << " 3rd lay size: " << thirdLayers.size();
0173
0174 auto const& doublets = thePairGenerator->doublets(region, ev, es, pairLayers);
0175 LogTrace("MultiHitGeneratorFromChi2") << "";
0176 if (not doublets or doublets->empty()) {
0177
0178 return;
0179 }
0180
0181 assert(theLayerCache);
0182 hitSets(region, result, *doublets, thirdLayers, *theLayerCache, cache);
0183 }
0184
0185 void MultiHitGeneratorFromChi2::hitSets(const TrackingRegion& region,
0186 OrderedMultiHits& result,
0187 const HitDoublets& doublets,
0188 const std::vector<SeedingLayerSetsHits::SeedingLayer>& thirdLayers,
0189 LayerCacheType& layerCache,
0190 cacheHits& refittedHitStorage) {
0191 int size = thirdLayers.size();
0192 const RecHitsSortedInPhi* thirdHitMap[size];
0193 vector<const DetLayer*> thirdLayerDetLayer(size, nullptr);
0194 for (int il = 0; il < size; ++il) {
0195 thirdHitMap[il] = &layerCache(thirdLayers[il], region);
0196
0197 thirdLayerDetLayer[il] = thirdLayers[il].detLayer();
0198 }
0199 hitSets(region, result, doublets, thirdHitMap, thirdLayerDetLayer, size, refittedHitStorage);
0200 }
0201
0202 void MultiHitGeneratorFromChi2::hitTriplets(const TrackingRegion& region,
0203 OrderedMultiHits& result,
0204 const HitDoublets& doublets,
0205 const RecHitsSortedInPhi** thirdHitMap,
0206 const std::vector<const DetLayer*>& thirdLayerDetLayer,
0207 const int nThirdLayers) {
0208 hitSets(region, result, doublets, thirdHitMap, thirdLayerDetLayer, nThirdLayers, cache);
0209 }
0210
0211 void MultiHitGeneratorFromChi2::hitSets(const TrackingRegion& region,
0212 OrderedMultiHits& result,
0213 const HitDoublets& doublets,
0214 const RecHitsSortedInPhi** thirdHitMap,
0215 const std::vector<const DetLayer*>& thirdLayerDetLayer,
0216 const int nThirdLayers,
0217 cacheHits& refittedHitStorage) {
0218 #ifdef EDM_ML_DEBUG
0219 unsigned int debug_Id0 = detIdsToDebug[0];
0220 unsigned int debug_Id1 = detIdsToDebug[1];
0221 unsigned int debug_Id2 = detIdsToDebug[2];
0222 #endif
0223
0224 std::array<bool, 3> bl;
0225 bl[0] = doublets.innerLayer().isBarrel;
0226 bl[1] = doublets.outerLayer().isBarrel;
0227
0228
0229
0230
0231 std::vector<KDTreeNodeInfo<RecHitsSortedInPhi::HitIter, 2> > layerTree;
0232 std::vector<RecHitsSortedInPhi::HitIter> foundNodes;
0233 foundNodes.reserve(100);
0234 declareDynArray(KDTreeLinkerAlgo<RecHitsSortedInPhi::HitIter>, nThirdLayers, hitTree);
0235 declareDynArray(LayerRZPredictions, nThirdLayers, mapPred);
0236 std::vector<float> rzError(nThirdLayers);
0237
0238 const float maxDelphi = region.ptMin() < 0.3f ? float(M_PI) / 4.f : float(M_PI) / 8.f;
0239 const float maxphi = M_PI + maxDelphi, minphi = -maxphi;
0240 const float safePhi = M_PI - maxDelphi;
0241
0242
0243 for (int il = 0; il < nThirdLayers; il++) {
0244 LogTrace("MultiHitGeneratorFromChi2")
0245 << "considering third layer: with hits: " << thirdHitMap[il]->all().second - thirdHitMap[il]->all().first;
0246 const DetLayer* layer = thirdLayerDetLayer[il];
0247 LayerRZPredictions& predRZ = mapPred[il];
0248 predRZ.line.initLayer(layer);
0249 predRZ.line.initTolerance(extraHitRZtolerance);
0250
0251
0252 auto const& layer3 = *thirdHitMap[il];
0253 layerTree.clear();
0254 float minz = 999999.0f, maxz = -minz;
0255 float maxErr = 0.0f;
0256 if (!layer3.empty()) {
0257 minz = layer3.v[0];
0258 maxz = minz;
0259 for (auto i = 0U; i < layer3.size(); ++i) {
0260 auto hi = layer3.theHits.begin() + i;
0261 auto angle = layer3.phi(i);
0262 auto myz = layer3.v[i];
0263 #ifdef EDM_ML_DEBUG
0264 IfLogTrace(hi->hit()->rawId() == debug_Id2, "MultiHitGeneratorFromChi2")
0265 << "filling KDTree with hit in id=" << debug_Id2 << " with pos: " << hi->hit()->globalPosition()
0266 << " phi=" << hi->hit()->globalPosition().phi() << " z=" << hi->hit()->globalPosition().z()
0267 << " r=" << hi->hit()->globalPosition().perp();
0268 #endif
0269
0270 if (myz < minz) {
0271 minz = myz;
0272 } else {
0273 if (myz > maxz) {
0274 maxz = myz;
0275 }
0276 }
0277 auto myerr = layer3.dv[i];
0278 if (myerr > maxErr) {
0279 maxErr = myerr;
0280 }
0281 layerTree.push_back(KDTreeNodeInfo<RecHitsSortedInPhi::HitIter, 2>(hi, angle, myz));
0282
0283 if (angle > safePhi)
0284 layerTree.push_back(
0285 KDTreeNodeInfo<RecHitsSortedInPhi::HitIter, 2>(hi, float(angle - Geom::twoPi()), float(myz)));
0286 else if (angle < -safePhi)
0287 layerTree.push_back(
0288 KDTreeNodeInfo<RecHitsSortedInPhi::HitIter, 2>(hi, float(angle + Geom::twoPi()), float(myz)));
0289 }
0290 }
0291 KDTreeBox phiZ(minphi, maxphi, minz - 0.01f, maxz + 0.01f);
0292
0293 hitTree[il].build(layerTree, phiZ);
0294 rzError[il] = maxErr;
0295 }
0296
0297
0298
0299 auto curv = PixelRecoUtilities::curvature(1. / region.ptMin(), *bfield);
0300
0301 LogTrace("MultiHitGeneratorFromChi2") << "doublet size=" << doublets.size() << std::endl;
0302
0303
0304 auto filterHit = [&](TrackingRecHit const* hit, GlobalVector const& dir) -> bool {
0305 auto hh = reinterpret_cast<BaseTrackerRecHit const*>(hit);
0306 if (
0307 hh->geographicalId().subdetId() == SiStripDetId::TIB || hh->geographicalId().subdetId() == SiStripDetId::TID
0308
0309
0310 ) {
0311
0312 if (hh->isMatched()) {
0313 const SiStripMatchedRecHit2D* matchedHit = reinterpret_cast<const SiStripMatchedRecHit2D*>(hh);
0314 if (filter->isCompatible(DetId(matchedHit->monoId()), matchedHit->monoCluster(), dir) == 0 ||
0315 filter->isCompatible(DetId(matchedHit->stereoId()), matchedHit->stereoCluster(), dir) == 0)
0316 return false;
0317 } else if (hh->isProjected()) {
0318 const ProjectedSiStripRecHit2D* precHit = reinterpret_cast<const ProjectedSiStripRecHit2D*>(hh);
0319 if (filter->isCompatible(precHit->originalHit(), dir) == 0)
0320 return false;
0321 } else if (hh->isSingle()) {
0322 const SiStripRecHit2D* recHit = reinterpret_cast<const SiStripRecHit2D*>(hh);
0323 if (filter->isCompatible(*recHit, dir) == 0)
0324 return false;
0325 }
0326 }
0327 return true;
0328 };
0329
0330
0331 for (std::size_t ip = 0; ip != doublets.size(); ip++) {
0332 int foundTripletsFromPair = 0;
0333 bool usePair = false;
0334 cacheHitPointer bestH2;
0335 float minChi2 = std::numeric_limits<float>::max();
0336
0337 SeedingHitSet::ConstRecHitPointer oriHit0 = doublets.hit(ip, HitDoublets::inner);
0338 SeedingHitSet::ConstRecHitPointer oriHit1 = doublets.hit(ip, HitDoublets::outer);
0339
0340 HitOwnPtr hit0(*oriHit0);
0341 HitOwnPtr hit1(*oriHit1);
0342 GlobalPoint gp0 = doublets.gp(ip, HitDoublets::inner);
0343 GlobalPoint gp1 = doublets.gp(ip, HitDoublets::outer);
0344
0345 #ifdef EDM_ML_DEBUG
0346 bool debugPair = oriHit0->rawId() == debug_Id0 && oriHit1->rawId() == debug_Id1;
0347 #endif
0348 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << endl
0349 << endl
0350 << "found new pair with ids " << oriHit0->rawId() << " "
0351 << oriHit1->rawId() << " with pos: " << gp0 << " " << gp1;
0352
0353 if (refitHits) {
0354 TrajectoryStateOnSurface tsos0, tsos1;
0355 assert(!hit0.isOwn());
0356 assert(!hit1.isOwn());
0357 #ifdef EDM_ML_DEBUG
0358 refit2Hits(hit0, hit1, tsos0, tsos1, region, nomField, debugPair);
0359 #else
0360 refit2Hits(hit0, hit1, tsos0, tsos1, region, nomField, false);
0361 #endif
0362
0363 bool passFilterHit0 = filterHit(hit0->hit(), tsos0.globalMomentum());
0364 IfLogTrace(debugPair && !passFilterHit0, "MultiHitGeneratorFromChi2") << "hit0 did not pass cluster shape filter";
0365 if (!passFilterHit0)
0366 continue;
0367 bool passFilterHit1 = filterHit(hit1->hit(), tsos1.globalMomentum());
0368 IfLogTrace(debugPair && !passFilterHit1, "MultiHitGeneratorFromChi2") << "hit1 did not pass cluster shape filter";
0369 if (!passFilterHit1)
0370 continue;
0371
0372 hit0.reset((SeedingHitSet::RecHitPointer)(cloner(*hit0, tsos0)));
0373 hit1.reset((SeedingHitSet::RecHitPointer)(cloner(*hit1, tsos1)));
0374
0375 #ifdef EDM_ML_DEBUG
0376 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2")
0377 << "charge=" << tsos0.charge() << std::endl
0378 << "state1 pt=" << tsos0.globalMomentum().perp() << " eta=" << tsos0.globalMomentum().eta()
0379 << " phi=" << tsos0.globalMomentum().phi() << std::endl
0380 << "state2 pt=" << tsos1.globalMomentum().perp() << " eta=" << tsos1.globalMomentum().eta()
0381 << " phi=" << tsos1.globalMomentum().phi() << std::endl
0382 << "positions after refitting: " << hit0->globalPosition() << " " << hit1->globalPosition();
0383 #endif
0384 } else {
0385
0386 hit0.reset((BaseTrackerRecHit*)hit0->clone());
0387 hit1.reset((BaseTrackerRecHit*)hit1->clone());
0388 }
0389
0390 assert(hit0.isOwn());
0391 assert(hit1.isOwn());
0392
0393
0394 SimpleLineRZ line(PixelRecoPointRZ(gp0.perp(), gp0.z()), PixelRecoPointRZ(gp1.perp(), gp1.z()));
0395 ThirdHitPredictionFromCircle predictionRPhi(gp0, gp1, extraHitRPhitolerance);
0396
0397 auto toPos = std::signbit(gp1.z() - gp0.z());
0398
0399
0400 Range pairCurvature = predictionRPhi.curvature(region.originRBound());
0401
0402 if (!intersect(pairCurvature, Range(-curv, curv))) {
0403 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2")
0404 << "curvature cut: curv=" << curv << " gc=(" << pairCurvature.first << ", " << pairCurvature.second << ")";
0405 continue;
0406 }
0407
0408 std::array<GlobalPoint, 3> gp;
0409 std::array<GlobalError, 3> ge;
0410 gp[0] = hit0->globalPosition();
0411 ge[0] = hit0->globalPositionError();
0412 gp[1] = hit1->globalPosition();
0413 ge[1] = hit1->globalPositionError();
0414
0415
0416 for (int il = 0; (il < nThirdLayers) & (!usePair); il++) {
0417 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2")
0418 << "cosider layer:"
0419 << " for this pair. Location: " << thirdLayerDetLayer[il]->location();
0420
0421 if (hitTree[il].empty()) {
0422 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << "empty hitTree";
0423 continue;
0424 }
0425
0426 cacheHitPointer bestL2;
0427 float chi2FromThisLayer = std::numeric_limits<float>::max();
0428
0429 const DetLayer* layer = thirdLayerDetLayer[il];
0430
0431 auto const& layer3 = *thirdHitMap[il];
0432 bool barrelLayer = layer3.isBarrel;
0433 bl[2] = layer3.isBarrel;
0434
0435 if ((!barrelLayer) & (toPos != std::signbit(layer->position().z())))
0436 continue;
0437
0438 LayerRZPredictions& predRZ = mapPred[il];
0439 predRZ.line.initPropagator(&line);
0440
0441
0442
0443 Range rzRange = predRZ.line();
0444
0445 if (rzRange.first >= rzRange.second) {
0446 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << "rzRange empty";
0447 continue;
0448 }
0449
0450
0451 if (!intersect(rzRange, predRZ.line.detSize())) {
0452 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << "rzRange and detector do not intersect";
0453 continue;
0454 }
0455 Range radius = barrelLayer ? predRZ.line.detRange() : rzRange;
0456
0457
0458 Range phiRange;
0459 if (useFixedPreFiltering) {
0460
0461
0462 float phi0 = oriHit0->globalPosition().phi();
0463 phiRange = Range(phi0 - dphi, phi0 + dphi);
0464 } else {
0465
0466 if (pairCurvature.first < 0. && pairCurvature.second < 0.) {
0467 radius.swap();
0468 } else if (pairCurvature.first >= 0. && pairCurvature.second >= 0.) {
0469 ;
0470 } else {
0471 radius.first = radius.second;
0472 }
0473 auto phi12 = predictionRPhi.phi(pairCurvature.first, radius.first);
0474 auto phi22 = predictionRPhi.phi(pairCurvature.second, radius.second);
0475 phi12 = normalizedPhi(phi12);
0476 phi22 = proxim(phi22, phi12);
0477 phiRange = Range(phi12, phi22);
0478 phiRange.sort();
0479 }
0480
0481 float prmin = phiRange.min(), prmax = phiRange.max();
0482
0483 if (prmax - prmin > maxDelphi) {
0484 auto prm = phiRange.mean();
0485 prmin = prm - 0.5f * maxDelphi;
0486 prmax = prm + 0.5f * maxDelphi;
0487 }
0488
0489
0490 using Hit = RecHitsSortedInPhi::Hit;
0491 foundNodes.clear();
0492
0493 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << "defining kd tree box";
0494
0495 if (barrelLayer) {
0496 KDTreeBox phiZ(prmin - extraPhiKDBox,
0497 prmax + extraPhiKDBox,
0498 float(rzRange.min() - fnSigmaRZ * rzError[il] - extraZKDBox),
0499 float(rzRange.max() + fnSigmaRZ * rzError[il] + extraZKDBox));
0500 hitTree[il].search(phiZ, foundNodes);
0501
0502 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2")
0503 << "kd tree box bounds, phi: " << prmin - extraPhiKDBox << "," << prmax + extraPhiKDBox
0504 << " z: " << rzRange.min() - fnSigmaRZ * rzError[il] - extraZKDBox << ","
0505 << rzRange.max() + fnSigmaRZ * rzError[il] + extraZKDBox << " rzRange: " << rzRange.min() << ","
0506 << rzRange.max();
0507
0508 } else {
0509 KDTreeBox phiR(prmin - extraPhiKDBox,
0510 prmax + extraPhiKDBox,
0511 float(rzRange.min() - fnSigmaRZ * rzError[il] - extraRKDBox),
0512 float(rzRange.max() + fnSigmaRZ * rzError[il] + extraRKDBox));
0513 hitTree[il].search(phiR, foundNodes);
0514
0515 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2")
0516 << "kd tree box bounds, phi: " << prmin - extraPhiKDBox << "," << prmax + extraPhiKDBox
0517 << " r: " << rzRange.min() - fnSigmaRZ * rzError[il] - extraRKDBox << ","
0518 << rzRange.max() + fnSigmaRZ * rzError[il] + extraRKDBox << " rzRange: " << rzRange.min() << ","
0519 << rzRange.max();
0520 }
0521
0522 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << "kd tree box size: " << foundNodes.size();
0523
0524
0525 for (std::vector<RecHitsSortedInPhi::HitIter>::iterator ih = foundNodes.begin();
0526 ih != foundNodes.end() && !usePair;
0527 ++ih) {
0528 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << endl << "triplet candidate";
0529
0530 const RecHitsSortedInPhi::HitIter KDdata = *ih;
0531
0532 auto oriHit2 = KDdata->hit();
0533 auto kk = KDdata - layer3.theHits.begin();
0534 cacheHitPointer hit2;
0535 auto gp2 = layer3.gp(kk);
0536 if (refitHits) {
0537
0538
0539 GlobalVector initMomentum(gp2 - gp1);
0540 initMomentum /= initMomentum.perp();
0541
0542
0543 bool passFilterHit2 = filterHit(oriHit2->hit(), initMomentum);
0544 if (!passFilterHit2)
0545 continue;
0546 TrajectoryStateOnSurface state(GlobalTrajectoryParameters(gp2, initMomentum, 1, &ufield),
0547 *oriHit2->surface());
0548 hit2.reset((SeedingHitSet::RecHitPointer)(cloner(*oriHit2, state)));
0549
0550 } else {
0551
0552 hit2.reset((BaseTrackerRecHit*)oriHit2->clone());
0553 }
0554
0555
0556 gp[2] = hit2->globalPosition();
0557 ge[2] = hit2->globalPositionError();
0558 RZLine rzLine(gp, ge, bl);
0559 float chi2 = rzLine.chi2();
0560
0561 #ifdef EDM_ML_DEBUG
0562 bool debugTriplet = debugPair && hit2->rawId() == debug_Id2;
0563 #endif
0564 IfLogTrace(debugTriplet, "MultiHitGeneratorFromChi2")
0565 << endl
0566 << "triplet candidate in debug id" << std::endl
0567 << "hit in id=" << hit2->rawId() << " (from KDTree) with pos: " << KDdata->hit()->globalPosition()
0568 << " refitted: " << hit2->globalPosition() << " chi2: " << chi2;
0569
0570 if ((chi2 > maxChi2) | edm::isNotFinite(chi2))
0571 continue;
0572
0573 if (chi2VsPtCut) {
0574 FastCircle theCircle(gp[2], gp[1], gp[0]);
0575 float tesla0 = 0.1f * nomField;
0576 float rho = theCircle.rho();
0577 float cm2GeV = 0.01f * 0.3f * tesla0;
0578 float pt = cm2GeV * rho;
0579 IfLogTrace(debugTriplet, "MultiHitGeneratorFromChi2") << "triplet pT=" << pt;
0580 if (pt < region.ptMin())
0581 continue;
0582
0583 if (chi2_cuts.size() == 4) {
0584 if ((pt > pt_interv[0] && pt <= pt_interv[1] && chi2 > chi2_cuts[0]) ||
0585 (pt > pt_interv[1] && pt <= pt_interv[2] && chi2 > chi2_cuts[1]) ||
0586 (pt > pt_interv[2] && pt <= pt_interv[3] && chi2 > chi2_cuts[2]) ||
0587 (pt > pt_interv[3] && chi2 > chi2_cuts[3]))
0588 continue;
0589 }
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599
0600
0601
0602
0603
0604
0605 }
0606
0607 if (theMaxElement != 0 && result.size() >= theMaxElement) {
0608 result.clear();
0609 edm::LogError("TooManyTriplets") << " number of triples exceed maximum. no triplets produced.";
0610 return;
0611 }
0612 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << "triplet made";
0613
0614
0615
0616
0617
0618
0619
0620
0621
0622 bestL2 = std::move(hit2);
0623 chi2FromThisLayer = chi2;
0624 foundTripletsFromPair++;
0625 if (foundTripletsFromPair >= 2) {
0626 usePair = true;
0627 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << "using pair";
0628 break;
0629 }
0630 }
0631
0632 if (usePair)
0633 break;
0634 else {
0635
0636 if (chi2FromThisLayer < minChi2) {
0637 bestH2 = std::move(bestL2);
0638 minChi2 = chi2FromThisLayer;
0639 }
0640
0641
0642
0643
0644
0645
0646 }
0647
0648 }
0649
0650 if (foundTripletsFromPair == 0)
0651 continue;
0652
0653
0654 IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << "Done seed #" << result.size();
0655 if (usePair)
0656 result.push_back(SeedingHitSet(oriHit0, oriHit1));
0657 else {
0658 assert(1 == foundTripletsFromPair);
0659 assert(bestH2);
0660 result.emplace_back(&*hit0, &*hit1, &*bestH2);
0661 assert(hit0.isOwn());
0662 assert(hit1.isOwn());
0663 refittedHitStorage.emplace_back(const_cast<BaseTrackerRecHit*>(hit0.release()));
0664 refittedHitStorage.emplace_back(const_cast<BaseTrackerRecHit*>(hit1.release()));
0665 refittedHitStorage.emplace_back(std::move(bestH2));
0666 assert(hit0.empty());
0667 assert(hit1.empty());
0668 assert(!bestH2);
0669 }
0670
0671
0672 }
0673 LogTrace("MultiHitGeneratorFromChi2") << "triplet size=" << result.size();
0674
0675 }
0676
0677 void MultiHitGeneratorFromChi2::refit2Hits(HitOwnPtr& hit1,
0678 HitOwnPtr& hit2,
0679 TrajectoryStateOnSurface& state1,
0680 TrajectoryStateOnSurface& state2,
0681 const TrackingRegion& region,
0682 float nomField,
0683 bool isDebug) {
0684
0685 const GlobalPoint& gp0 = region.origin();
0686 GlobalPoint gp1 = hit1->globalPosition();
0687 GlobalPoint gp2 = hit2->globalPosition();
0688
0689 IfLogTrace(isDebug, "MultiHitGeneratorFromChi2")
0690 << "positions before refitting: " << hit1->globalPosition() << " " << hit2->globalPosition();
0691
0692 FastCircle theCircle(gp2, gp1, gp0);
0693 GlobalPoint cc(theCircle.x0(), theCircle.y0(), 0);
0694 float tesla0 = 0.1f * nomField;
0695 float rho = theCircle.rho();
0696 float cm2GeV = 0.01f * 0.3f * tesla0;
0697 float pt = cm2GeV * rho;
0698
0699 GlobalVector vec20 = gp2 - gp0;
0700
0701
0702 GlobalVector p0(gp0.y() - cc.y(), -gp0.x() + cc.x(), 0.);
0703 p0 = p0 * pt / p0.perp();
0704 GlobalVector p1(gp1.y() - cc.y(), -gp1.x() + cc.x(), 0.);
0705 p1 = p1 * pt / p1.perp();
0706 GlobalVector p2(gp2.y() - cc.y(), -gp2.x() + cc.x(), 0.);
0707 p2 = p2 * pt / p2.perp();
0708
0709
0710 if ((p0.x() * (gp1.x() - gp0.x()) + p0.y() * (gp1.y() - gp0.y())) < 0) {
0711 p0 *= -1.;
0712 p1 *= -1.;
0713 p2 *= -1.;
0714 }
0715
0716
0717 auto zv = vec20.z() / vec20.perp();
0718 p0 = GlobalVector(p0.x(), p0.y(), p0.perp() * zv);
0719 p1 = GlobalVector(p1.x(), p1.y(), p1.perp() * zv);
0720 p2 = GlobalVector(p2.x(), p2.y(), p2.perp() * zv);
0721
0722
0723 TrackCharge q = 1;
0724 if ((gp1 - cc).x() * p1.y() - (gp1 - cc).y() * p1.x() > 0)
0725 q = -q;
0726
0727 TrajectoryStateOnSurface(GlobalTrajectoryParameters(gp1, p1, q, &ufield), *hit1->surface()).swap(state1);
0728
0729
0730 TrajectoryStateOnSurface(GlobalTrajectoryParameters(gp2, p2, q, &ufield), *hit2->surface()).swap(state2);
0731
0732 }
0733
0734
0735
0736
0737
0738
0739
0740
0741
0742
0743
0744
0745
0746
0747
0748
0749
0750
0751
0752
0753
0754
0755
0756
0757
0758
0759
0760
0761
0762
0763
0764
0765
0766
0767
0768
0769
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782
0783
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802