TkHitPairs/src/HitPairGeneratorFromLayerPair.cc

0001 #include "RecoTracker/TkHitPairs/interface/HitPairGeneratorFromLayerPair.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0005 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0006 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0007
0008 #include "RecoTracker/TkTrackingRegions/interface/HitRZCompatibility.h"
0009 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
0010 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegionBase.h"
0011 #include "RecoTracker/TkHitPairs/interface/OrderedHitPairs.h"
0012 #include "RecoTracker/TkHitPairs/src/InnerDeltaPhi.h"
0013
0014 #include "FWCore/Framework/interface/Event.h"
0015
0016 using namespace GeomDetEnumerators;
0017 using namespace std;
0018
0019 typedef PixelRecoRange<float> Range;
0020
0021 namespace {
0022   template <class T>
0023   inline T sqr(T t) {
0024     return t * t;
0025   }
0026 }  // namespace
0027
0028 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0029 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0030 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0031 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0032 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisationMaker.h"
0033 #include "RecoTracker/Record/interface/TrackerMultipleScatteringRecord.h"
0034
0035 HitPairGeneratorFromLayerPair::HitPairGeneratorFromLayerPair(
0036     edm::ConsumesCollector iC, unsigned int inner, unsigned int outer, LayerCacheType* layerCache, unsigned int max)
0037     : theLayerCache(layerCache),
0038       theFieldToken(iC.esConsumes()),
0039       theMSMakerToken(iC.esConsumes()),
0040       theOuterLayer(outer),
0041       theInnerLayer(inner),
0042       theMaxElement(max) {}
0043
0044 HitPairGeneratorFromLayerPair::~HitPairGeneratorFromLayerPair() {}
0045
0046 // devirtualizer
0047 #include <tuple>
0048 namespace {
0049
0050   template <typename Algo>
0051   struct Kernel {
0052     using Base = HitRZCompatibility;
0053     void set(Base const* a) {
0054       assert(a->algo() == Algo::me);
0055       checkRZ = reinterpret_cast<Algo const*>(a);
0056     }
0057
0058     void operator()(int b, int e, const RecHitsSortedInPhi& innerHitsMap, bool* ok) const {
0059       constexpr float nSigmaRZ = 3.46410161514f;  // std::sqrt(12.f);
0060       for (int i = b; i != e; ++i) {
0061         Range allowed = checkRZ->range(innerHitsMap.u[i]);
0062         float vErr = nSigmaRZ * innerHitsMap.dv[i];
0063         Range hitRZ(innerHitsMap.v[i] - vErr, innerHitsMap.v[i] + vErr);
0064         Range crossRange = allowed.intersection(hitRZ);
0065         ok[i - b] = !crossRange.empty();
0066       }
0067     }
0068     Algo const* checkRZ;
0069   };
0070
0071   template <typename... Args>
0072   using Kernels = std::tuple<Kernel<Args>...>;
0073
0074 }  // namespace
0075
0076 bool HitPairGeneratorFromLayerPair::hitPairs(const TrackingRegion& region,
0077                                              OrderedHitPairs& result,
0078                                              const edm::Event& iEvent,
0079                                              const edm::EventSetup& iSetup,
0080                                              Layers layers) {
0081   auto const& ds = doublets(region, iEvent, iSetup, layers);
0082   if (not ds) {
0083     return false;
0084   }
0085   for (std::size_t i = 0; i != ds->size(); ++i) {
0086     result.push_back(OrderedHitPair(ds->hit(i, HitDoublets::inner), ds->hit(i, HitDoublets::outer)));
0087   }
0088   if (theMaxElement != 0 && result.size() >= theMaxElement) {
0089     result.clear();
0090     edm::LogError("TooManyPairs") << "number of pairs exceed maximum, no pairs produced";
0091     return false;
0092   }
0093   return true;
0094 }
0095
0096 std::optional<HitDoublets> HitPairGeneratorFromLayerPair::doublets(const TrackingRegion& region,
0097                                                                    const edm::Event& iEvent,
0098                                                                    const edm::EventSetup& iSetup,
0099                                                                    const Layer& innerLayer,
0100                                                                    const Layer& outerLayer,
0101                                                                    LayerCacheType& layerCache) {
0102   const RecHitsSortedInPhi& innerHitsMap = layerCache(innerLayer, region);
0103   if (innerHitsMap.empty())
0104     return HitDoublets(innerHitsMap, innerHitsMap);
0105
0106   const RecHitsSortedInPhi& outerHitsMap = layerCache(outerLayer, region);
0107   if (outerHitsMap.empty())
0108     return HitDoublets(innerHitsMap, outerHitsMap);
0109   const auto& field = iSetup.getData(theFieldToken);
0110   const auto& msmaker = iSetup.getData(theMSMakerToken);
0111   HitDoublets result(innerHitsMap, outerHitsMap);
0112   result.reserve(std::max(innerHitsMap.size(), outerHitsMap.size()));
0113   bool succeeded = doublets(region,
0114                             *innerLayer.detLayer(),
0115                             *outerLayer.detLayer(),
0116                             innerHitsMap,
0117                             outerHitsMap,
0118                             field,
0119                             msmaker,
0120                             theMaxElement,
0121                             result);
0122   if (succeeded) {
0123     return result;
0124   } else {
0125     return std::nullopt;
0126   }
0127 }
0128
0129 bool HitPairGeneratorFromLayerPair::doublets(const TrackingRegion& region,
0130                                              const DetLayer& innerHitDetLayer,
0131                                              const DetLayer& outerHitDetLayer,
0132                                              const RecHitsSortedInPhi& innerHitsMap,
0133                                              const RecHitsSortedInPhi& outerHitsMap,
0134                                              const MagneticField& field,
0135                                              const MultipleScatteringParametrisationMaker& msmaker,
0136                                              const unsigned int theMaxElement,
0137                                              HitDoublets& result) {
0138   //  HitDoublets result(innerHitsMap,outerHitsMap); result.reserve(std::max(innerHitsMap.size(),outerHitsMap.size()));
0139   typedef RecHitsSortedInPhi::Hit Hit;
0140   InnerDeltaPhi deltaPhi(outerHitDetLayer, innerHitDetLayer, region, field, msmaker);
0141
0142   // std::cout << "layers " << theInnerLayer.detLayer()->seqNum()  << " " << outerLayer.detLayer()->seqNum() << std::endl;
0143
0144   // constexpr float nSigmaRZ = std::sqrt(12.f);
0145   constexpr float nSigmaPhi = 3.f;
0146   for (int io = 0; io != int(outerHitsMap.theHits.size()); ++io) {
0147     if (!deltaPhi.prefilter(outerHitsMap.x[io], outerHitsMap.y[io]))
0148       continue;
0149     Hit const& ohit = outerHitsMap.theHits[io].hit();
0150     PixelRecoRange<float> phiRange =
0151         deltaPhi(outerHitsMap.x[io], outerHitsMap.y[io], outerHitsMap.z[io], nSigmaPhi * outerHitsMap.drphi[io]);
0152
0153     if (phiRange.empty())
0154       continue;
0155
0156     std::unique_ptr<const HitRZCompatibility> checkRZ =
0157         region.checkRZ(&innerHitDetLayer,
0158                        ohit,
0159                        &outerHitDetLayer,
0160                        outerHitsMap.rv(io),
0161                        outerHitsMap.z[io],
0162                        outerHitsMap.isBarrel ? outerHitsMap.du[io] : outerHitsMap.dv[io],
0163                        outerHitsMap.isBarrel ? outerHitsMap.dv[io] : outerHitsMap.du[io]);
0164     if (!checkRZ)
0165       continue;
0166
0167     Kernels<HitZCheck, HitRCheck, HitEtaCheck> kernels;
0168
0169     auto innerRange = innerHitsMap.doubleRange(phiRange.min(), phiRange.max());
0170     LogDebug("HitPairGeneratorFromLayerPair")
0171         << "preparing for combination of: " << innerRange[1] - innerRange[0] + innerRange[3] - innerRange[2]
0172         << " inner and: " << outerHitsMap.theHits.size() << " outter";
0173     for (int j = 0; j < 3; j += 2) {
0174       auto b = innerRange[j];
0175       auto e = innerRange[j + 1];
0176       if (e == b)
0177         continue;
0178       bool ok[e - b];
0179       switch (checkRZ->algo()) {
0180         case (HitRZCompatibility::zAlgo):
0181           std::get<0>(kernels).set(checkRZ.get());
0182           std::get<0>(kernels)(b, e, innerHitsMap, ok);
0183           break;
0184         case (HitRZCompatibility::rAlgo):
0185           std::get<1>(kernels).set(checkRZ.get());
0186           std::get<1>(kernels)(b, e, innerHitsMap, ok);
0187           break;
0188         case (HitRZCompatibility::etaAlgo):
0189           std::get<2>(kernels).set(checkRZ.get());
0190           std::get<2>(kernels)(b, e, innerHitsMap, ok);
0191           break;
0192       }
0193       for (int i = 0; i != e - b; ++i) {
0194         if (!ok[i])
0195           continue;
0196         if (theMaxElement != 0 && result.size() >= theMaxElement) {
0197           result.clear();
0198           edm::LogError("TooManyPairs") << "number of pairs exceed maximum, no pairs produced";
0199           return false;
0200         }
0201         result.add(b + i, io);
0202       }
0203     }
0204   }
0205   LogDebug("HitPairGeneratorFromLayerPair") << " total number of pairs provided back: " << result.size();
0206   result.shrink_to_fit();
0207   return true;
0208 }