Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-01 22:41:02

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 void 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   for (std::size_t i = 0; i != ds.size(); ++i) {
0083     result.push_back(OrderedHitPair(ds.hit(i, HitDoublets::inner), ds.hit(i, HitDoublets::outer)));
0084   }
0085   if (theMaxElement != 0 && result.size() >= theMaxElement) {
0086     result.clear();
0087     edm::LogError("TooManyPairs") << "number of pairs exceed maximum, no pairs produced";
0088   }
0089 }
0090 
0091 HitDoublets HitPairGeneratorFromLayerPair::doublets(const TrackingRegion& region,
0092                                                     const edm::Event& iEvent,
0093                                                     const edm::EventSetup& iSetup,
0094                                                     const Layer& innerLayer,
0095                                                     const Layer& outerLayer,
0096                                                     LayerCacheType& layerCache) {
0097   const RecHitsSortedInPhi& innerHitsMap = layerCache(innerLayer, region);
0098   if (innerHitsMap.empty())
0099     return HitDoublets(innerHitsMap, innerHitsMap);
0100 
0101   const RecHitsSortedInPhi& outerHitsMap = layerCache(outerLayer, region);
0102   if (outerHitsMap.empty())
0103     return HitDoublets(innerHitsMap, outerHitsMap);
0104   const auto& field = iSetup.getData(theFieldToken);
0105   const auto& msmaker = iSetup.getData(theMSMakerToken);
0106   HitDoublets result(innerHitsMap, outerHitsMap);
0107   result.reserve(std::max(innerHitsMap.size(), outerHitsMap.size()));
0108   doublets(region,
0109            *innerLayer.detLayer(),
0110            *outerLayer.detLayer(),
0111            innerHitsMap,
0112            outerHitsMap,
0113            field,
0114            msmaker,
0115            theMaxElement,
0116            result);
0117 
0118   return result;
0119 }
0120 
0121 void HitPairGeneratorFromLayerPair::doublets(const TrackingRegion& region,
0122                                              const DetLayer& innerHitDetLayer,
0123                                              const DetLayer& outerHitDetLayer,
0124                                              const RecHitsSortedInPhi& innerHitsMap,
0125                                              const RecHitsSortedInPhi& outerHitsMap,
0126                                              const MagneticField& field,
0127                                              const MultipleScatteringParametrisationMaker& msmaker,
0128                                              const unsigned int theMaxElement,
0129                                              HitDoublets& result) {
0130   //  HitDoublets result(innerHitsMap,outerHitsMap); result.reserve(std::max(innerHitsMap.size(),outerHitsMap.size()));
0131   typedef RecHitsSortedInPhi::Hit Hit;
0132   InnerDeltaPhi deltaPhi(outerHitDetLayer, innerHitDetLayer, region, field, msmaker);
0133 
0134   // std::cout << "layers " << theInnerLayer.detLayer()->seqNum()  << " " << outerLayer.detLayer()->seqNum() << std::endl;
0135 
0136   // constexpr float nSigmaRZ = std::sqrt(12.f);
0137   constexpr float nSigmaPhi = 3.f;
0138   for (int io = 0; io != int(outerHitsMap.theHits.size()); ++io) {
0139     if (!deltaPhi.prefilter(outerHitsMap.x[io], outerHitsMap.y[io]))
0140       continue;
0141     Hit const& ohit = outerHitsMap.theHits[io].hit();
0142     PixelRecoRange<float> phiRange =
0143         deltaPhi(outerHitsMap.x[io], outerHitsMap.y[io], outerHitsMap.z[io], nSigmaPhi * outerHitsMap.drphi[io]);
0144 
0145     if (phiRange.empty())
0146       continue;
0147 
0148     std::unique_ptr<const HitRZCompatibility> checkRZ =
0149         region.checkRZ(&innerHitDetLayer,
0150                        ohit,
0151                        &outerHitDetLayer,
0152                        outerHitsMap.rv(io),
0153                        outerHitsMap.z[io],
0154                        outerHitsMap.isBarrel ? outerHitsMap.du[io] : outerHitsMap.dv[io],
0155                        outerHitsMap.isBarrel ? outerHitsMap.dv[io] : outerHitsMap.du[io]);
0156     if (!checkRZ)
0157       continue;
0158 
0159     Kernels<HitZCheck, HitRCheck, HitEtaCheck> kernels;
0160 
0161     auto innerRange = innerHitsMap.doubleRange(phiRange.min(), phiRange.max());
0162     LogDebug("HitPairGeneratorFromLayerPair")
0163         << "preparing for combination of: " << innerRange[1] - innerRange[0] + innerRange[3] - innerRange[2]
0164         << " inner and: " << outerHitsMap.theHits.size() << " outter";
0165     for (int j = 0; j < 3; j += 2) {
0166       auto b = innerRange[j];
0167       auto e = innerRange[j + 1];
0168       if (e == b)
0169         continue;
0170       bool ok[e - b];
0171       switch (checkRZ->algo()) {
0172         case (HitRZCompatibility::zAlgo):
0173           std::get<0>(kernels).set(checkRZ.get());
0174           std::get<0>(kernels)(b, e, innerHitsMap, ok);
0175           break;
0176         case (HitRZCompatibility::rAlgo):
0177           std::get<1>(kernels).set(checkRZ.get());
0178           std::get<1>(kernels)(b, e, innerHitsMap, ok);
0179           break;
0180         case (HitRZCompatibility::etaAlgo):
0181           std::get<2>(kernels).set(checkRZ.get());
0182           std::get<2>(kernels)(b, e, innerHitsMap, ok);
0183           break;
0184       }
0185       for (int i = 0; i != e - b; ++i) {
0186         if (!ok[i])
0187           continue;
0188         if (theMaxElement != 0 && result.size() >= theMaxElement) {
0189           result.clear();
0190           edm::LogError("TooManyPairs") << "number of pairs exceed maximum, no pairs produced";
0191           return;
0192         }
0193         result.add(b + i, io);
0194       }
0195     }
0196   }
0197   LogDebug("HitPairGeneratorFromLayerPair") << " total number of pairs provided back: " << result.size();
0198   result.shrink_to_fit();
0199 }