TkSeedGenerator/plugins/MultiHitGeneratorFromChi2.cc

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 }  // namespace
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")),  //extra window in ThirdHitRZPrediction range
0054       extraHitRPhitolerance(cfg.getParameter<double>(
0055           "extraHitRPhitolerance")),  //extra window in ThirdHitPredictionFromCircle range (divide by R to get phi)
0056       extraZKDBox(
0057           cfg.getParameter<double>("extraZKDBox")),  //extra windown in Z when building the KDTree box (used in barrel)
0058       extraRKDBox(
0059           cfg.getParameter<double>("extraRKDBox")),  //extra windown in R when building the KDTree box (used in endcap)
0060       extraPhiKDBox(cfg.getParameter<double>("extraPhiKDBox")),  //extra windown in Phi when building the KDTree box
0061       fnSigmaRZ(cfg.getParameter<double>(
0062           "fnSigmaRZ")),  //this multiplies the max hit error on the layer when building the KDTree box
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   //if (detIdsToDebug.size()<3) //fixme
0079 #else
0080   detIdsToDebug.push_back(0);
0081   detIdsToDebug.push_back(0);
0082   detIdsToDebug.push_back(0);
0083 #endif
0084   // 2014/02/11 mia:
0085   // we should get rid of the boolean parameter useSimpleMF,
0086   // and use only a string magneticField [instead of SimpleMagneticField]
0087   // or better an edm::ESInputTag (at the moment HLT does not handle ESInputTag)
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   // fixed phi filtering
0113   desc.add<bool>("useFixedPreFiltering", false);
0114   desc.add<double>("phiPreFiltering", 0.3);
0115
0116   // box properties
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   // refit&filter hits
0125   desc.add<bool>("refitHits", true);
0126   desc.add<std::string>("ClusterShapeHitFilterName", "ClusterShapeHitFilter");
0127   desc.add<std::string>("TTRHBuilder", "WithTrackAngle");
0128
0129   // chi2 cuts
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   // debugging
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);  // more than enough (never actually used)
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 }  // namespace
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     //  LogDebug("MultiHitGeneratorFromChi2") << "empy pairs";
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   //gc: these are all the layers compatible with the layer pairs (as defined in the config file)
0229
0230   //gc: initialize a KDTree per each 3rd layer
0231   std::vector<KDTreeNodeInfo<RecHitsSortedInPhi::HitIter, 2> > layerTree;  // re-used throughout
0232   std::vector<RecHitsSortedInPhi::HitIter> foundNodes;                     // re-used thoughout
0233   foundNodes.reserve(100);
0234   declareDynArray(KDTreeLinkerAlgo<RecHitsSortedInPhi::HitIter>, nThirdLayers, hitTree);
0235   declareDynArray(LayerRZPredictions, nThirdLayers, mapPred);
0236   float rzError[nThirdLayers];  //save maximum errors
0237
0238   const float maxDelphi = region.ptMin() < 0.3f ? float(M_PI) / 4.f : float(M_PI) / 8.f;  // FIXME move to config??
0239   const float maxphi = M_PI + maxDelphi, minphi = -maxphi;  // increase to cater for any range
0240   const float safePhi = M_PI - maxDelphi;                   // sideband
0241
0242   //gc: loop over each layer
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     //gc: now we take all hits in the layer and fill the KDTree
0252     auto const& layer3 = *thirdHitMap[il];  // Get iterators
0253     layerTree.clear();
0254     float minz = 999999.0f, maxz = -minz;  // Initialise to extreme values in case no hits
0255     float maxErr = 0.0f;
0256     if (!layer3.empty()) {
0257       minz = layer3.v[0];
0258       maxz = minz;  //In case there's only one hit on the layer
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         //use (phi,r) for endcaps rather than (phi,z)
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));  // save it
0282         // populate side-bands
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);  // declare our bounds
0292     //add fudge factors in case only one hit and also for floating-point inaccuracy
0293     hitTree[il].build(layerTree, phiZ);  // make KDtree
0294     rzError[il] = maxErr;                //save error
0295   }
0296   //gc: now we have initialized the KDTrees and we are out of the layer loop
0297
0298   //gc: this sets the minPt of the triplet
0299   auto curv = PixelRecoUtilities::curvature(1. / region.ptMin(), *bfield);
0300
0301   LogTrace("MultiHitGeneratorFromChi2") << "doublet size=" << doublets.size() << std::endl;
0302
0303   //fixme add pixels
0304   auto filterHit = [&](TrackingRecHit const* hit, GlobalVector const& dir) -> bool {
0305     auto hh = reinterpret_cast<BaseTrackerRecHit const*>(hit);
0306     if (  //hh->geographicalId().subdetId() > 2
0307         hh->geographicalId().subdetId() == SiStripDetId::TIB || hh->geographicalId().subdetId() == SiStripDetId::TID
0308         //|| hh->geographicalId().subdetId()==SiStripDetId::TOB
0309         //|| hh->geographicalId().subdetId()==SiStripDetId::TEC
0310     ) {
0311       // carefull " for matched and projected local of tsos != local for individual clusters...
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;  //FIXME??
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   //gc: now we loop over all pairs
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       // refit hits
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       // not refit clone anyhow
0386       hit0.reset((BaseTrackerRecHit*)hit0->clone());
0387       hit1.reset((BaseTrackerRecHit*)hit1->clone());
0388     }
0389
0390     assert(hit0.isOwn());
0391     assert(hit1.isOwn());
0392
0393     //gc: create the RZ line for the pair
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     //gc: this is the curvature of the two hits assuming the region
0400     Range pairCurvature = predictionRPhi.curvature(region.originRBound());
0401     //gc: intersect not only returns a bool but may change pairCurvature to intersection with curv
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     //gc: loop over all third layers compatible with the pair
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;  // Don't bother if no hits
0424       }
0425
0426       cacheHitPointer bestL2;
0427       float chi2FromThisLayer = std::numeric_limits<float>::max();
0428
0429       const DetLayer* layer = thirdLayerDetLayer[il];
0430       // bool barrelLayer = layer->location() == GeomDetEnumerators::barrel;
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       //gc: this takes the z at R-thick/2 and R+thick/2 according to
0442       //    the line from the two points and the adds the extra tolerance
0443       Range rzRange = predRZ.line();
0444
0445       if (rzRange.first >= rzRange.second) {
0446         IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << "rzRange empty";
0447         continue;
0448       }
0449       //gc: check that rzRange is compatible with detector bounds
0450       //    note that intersect may change rzRange to intersection with bounds
0451       if (!intersect(rzRange, predRZ.line.detSize())) {  // theDetSize = Range(-maxZ, maxZ);
0452         IfLogTrace(debugPair, "MultiHitGeneratorFromChi2") << "rzRange and detector do not intersect";
0453         continue;
0454       }
0455       Range radius = barrelLayer ? predRZ.line.detRange() : rzRange;
0456
0457       //gc: define the phi range of the hits
0458       Range phiRange;
0459       if (useFixedPreFiltering) {
0460         //gc: in this case it takes as range the phi of the outer
0461         //    hit +/- the phiPreFiltering value from cfg
0462         float phi0 = oriHit0->globalPosition().phi();
0463         phiRange = Range(phi0 - dphi, phi0 + dphi);
0464       } else {
0465         //gc: predictionRPhi uses the cosine rule to find the phi of the 3rd point at radius, assuming the pairCurvature range [-c,+c]
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       //gc: this is the place where hits in the compatible region are put in the foundNodes
0490       using Hit = RecHitsSortedInPhi::Hit;
0491       foundNodes.clear();  // Now recover hits in bounding box...
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       //gc: now we loop over the hits in the box for this layer
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) {  //fixme
0537
0538           //fitting all 3 hits takes too much time... do it quickly only for 3rd hit
0539           GlobalVector initMomentum(gp2 - gp1);
0540           initMomentum /= initMomentum.perp();  //set pT=1
0541
0542           //fixme add pixels
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           // not refit clone anyhow
0552           hit2.reset((BaseTrackerRecHit*)oriHit2->clone());
0553         }
0554
0555         //gc: add the chi2 cut
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         // should fix nan
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           // apparently this takes too much time...
0592           //      if (chi2_cuts.size()>1) {
0593           //        int ncuts = chi2_cuts.size();
0594           //        if ( pt<=pt_interv[0] && chi2 > chi2_cuts[0] ) continue;
0595           //        bool pass = true;
0596           //        for (int icut=1; icut<ncuts-1; icut++){
0597           //          if ( pt>pt_interv[icut-1] && pt<=pt_interv[icut] && chi2 > chi2_cuts[icut] ) pass=false;
0598           //        }
0599           //        if (!pass) continue;
0600           //        if ( pt>pt_interv[ncuts-2] && chi2 > chi2_cuts[ncuts-1] ) continue;
0601           //        if (hit0->rawId()==debug_Id0 && hit1->rawId()==debug_Id1 && hit2->rawId()==debug_Id2) {
0602           //          LogTrace("MultiHitGeneratorFromChi2") << "triplet passed chi2 vs pt cut" << std::endl;
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         //result.push_back(SeedingHitSet(hit0, hit1, hit2));
0614         /* no refit so keep only hit2
0615     assert(tripletFromThisLayer.empty());
0616     assert(hit0.isOwn()); assert(hit1.isOwn());assert(hit2.isOwn());
0617     tripletFromThisLayer.emplace_back(std::move(hit0));
0618     tripletFromThisLayer.emplace_back(std::move(hit1));
0619     tripletFromThisLayer.emplace_back(std::move(hit2));
0620     assert(hit0.isEmpty()); assert(hit1.isEmpty());assert(hit2.isEmpty());
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       }  //loop over hits in KDTree
0631
0632       if (usePair)
0633         break;
0634       else {
0635         //if there is one triplet in more than one layer, try picking the one with best chi2
0636         if (chi2FromThisLayer < minChi2) {
0637           bestH2 = std::move(bestL2);
0638           minChi2 = chi2FromThisLayer;
0639         }
0640         /*
0641     else {
0642       if (!bestH2 && foundTripletsFromPair>0)
0643         LogTrace("MultiHitGeneratorFromChi2") << "what?? " <<  minChi2 << ' '  << chi2FromThisLayer;
0644     }
0645     */
0646       }
0647
0648     }  //loop over layers
0649
0650     if (foundTripletsFromPair == 0)
0651       continue;
0652
0653     //push back only (max) once per pair
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     // LogTrace("MultiHitGeneratorFromChi2") << (usePair ? "pair " : "triplet ") << minChi2 <<' ' << refittedHitStorage.size();
0671
0672   }  //loop over pairs
0673   LogTrace("MultiHitGeneratorFromChi2") << "triplet size=" << result.size();
0674   // std::cout << "MultiHitGeneratorFromChi2 " << "triplet size=" << result.size() << std::endl;
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   //these need to be sorted in R
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   //if (isDebug) { cout << "vec20.eta=" << vec20.eta() << endl; }
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   //check sign according to scalar product
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   //now set z component
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   //get charge from vectorial product
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   // hit1.reset((SeedingHitSet::RecHitPointer)(cloner(*hit1,state1)));
0729
0730   TrajectoryStateOnSurface(GlobalTrajectoryParameters(gp2, p2, q, &ufield), *hit2->surface()).swap(state2);
0731   // hit2.reset((SeedingHitSet::RecHitPointer)(cloner(*hit2,state2)));
0732 }
0733
0734 /*
0735 void MultiHitGeneratorFromChi2::refit3Hits(HitOwnPtr & hit0,
0736                        HitOwnPtr & hit1,
0737                        HitOwnPtr & hit2,
0738                        TrajectoryStateOnSurface& state0,
0739                        TrajectoryStateOnSurface& state1,
0740                        TrajectoryStateOnSurface& state2,
0741                        float nomField, bool isDebug) {
0742
0743   //these need to be sorted in R
0744   GlobalPoint gp0 = hit0->globalPosition();
0745   GlobalPoint gp1 = hit1->globalPosition();
0746   GlobalPoint gp2 = hit2->globalPosition();
0747
0748   IfLogTrace(isDebug, "MultiHitGeneratorFromChi2") << "positions before refitting: " << hit0->globalPosition() << " " << hit1->globalPosition() << " " << hit2->globalPosition();
0749   }
0750
0751   FastCircle theCircle(gp2,gp1,gp0);
0752   GlobalPoint cc(theCircle.x0(),theCircle.y0(),0);
0753   float tesla0 = 0.1*nomField;
0754   float rho = theCircle.rho();
0755   float cm2GeV = 0.01 * 0.3*tesla0;
0756   float pt = cm2GeV * rho;
0757
0758   GlobalVector vec20 = gp2-gp0;
0759   //IfLogTrace(isDebug, "MultiHitGeneratorFromChi2") << "vec20.eta=" << vec20.eta();
0760
0761   GlobalVector p0( gp0.y()-cc.y(), -gp0.x()+cc.x(), 0. );
0762   p0 = p0*pt/p0.perp();
0763   GlobalVector p1( gp1.y()-cc.y(), -gp1.x()+cc.x(), 0. );
0764   p1 = p1*pt/p1.perp();
0765   GlobalVector p2( gp2.y()-cc.y(), -gp2.x()+cc.x(), 0. );
0766   p2 = p2*pt/p2.perp();
0767
0768   //check sign according to scalar product
0769   if ( (p0.x()*(gp1.x()-gp0.x())+p0.y()*(gp1.y()-gp0.y()) ) < 0 ) {
0770     p0*=-1.;
0771     p1*=-1.;
0772     p2*=-1.;
0773   }
0774
0775   //now set z component
0776   p0 = GlobalVector(p0.x(),p0.y(),p0.perp()/tan(vec20.theta()));
0777   p1 = GlobalVector(p1.x(),p1.y(),p1.perp()/tan(vec20.theta()));
0778   p2 = GlobalVector(p2.x(),p2.y(),p2.perp()/tan(vec20.theta()));
0779
0780   //get charge from vectorial product
0781   TrackCharge q = 1;
0782   if ((gp1-cc).x()*p1.y() - (gp1-cc).y()*p1.x() > 0) q =-q;
0783
0784   GlobalTrajectoryParameters kine0 = GlobalTrajectoryParameters(gp0, p0, q, &*bfield);
0785   state0 = TrajectoryStateOnSurface(kine0,*hit0->surface());
0786   hit0 = hit0->clone(state0);
0787
0788   GlobalTrajectoryParameters kine1 = GlobalTrajectoryParameters(gp1, p1, q, &*bfield);
0789   state1 = TrajectoryStateOnSurface(kine1,*hit1->surface());
0790   hit1 = hit1->clone(state1);
0791
0792   GlobalTrajectoryParameters kine2 = GlobalTrajectoryParameters(gp2, p2, q, &*bfield);
0793   state2 = TrajectoryStateOnSurface(kine2,*hit2->surface());
0794   hit2 = hit2->clone(state2);
0795
0796   IfLogTrace(isDebug, "MultiHitGeneratorFromChi2") << "charge=" << q << std::endl
0797                                                    << "state0 pt=" << state0.globalMomentum().perp() << " eta=" << state0.globalMomentum().eta()  << " phi=" << state0.globalMomentum().phi() << std::endl
0798                                                    << "state1 pt=" << state1.globalMomentum().perp() << " eta=" << state1.globalMomentum().eta()  << " phi=" << state1.globalMomentum().phi() << std::endl
0799                                                    << "state2 pt=" << state2.globalMomentum().perp() << " eta=" << state2.globalMomentum().eta()  << " phi=" << state2.globalMomentum().phi() << std::endl
0800                                                    << "positions after refitting: " << hit0->globalPosition() << " " << hit1->globalPosition() << " " << hit2->globalPosition();
0801 }
0802 */