TkHitPairs/interface/RecHitsSortedInPhi.h

0001 #ifndef RecHitsSortedInPhi_H
0002 #define RecHitsSortedInPhi_H
0003
0004 #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"
0005 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0006
0007 #include <vector>
0008 #include <array>
0009
0010 #include <cassert>
0011
0012 /** A RecHit container sorted in phi.
0013  *  Provides fast access for hits in a given phi window
0014  *  using binary search.
0015  */
0016
0017 class RecHitsSortedInPhi {
0018 public:
0019   typedef BaseTrackerRecHit const* Hit;
0020
0021   // A RecHit extension that caches the phi angle for fast access
0022   class HitWithPhi {
0023   public:
0024     HitWithPhi(const Hit& hit) : theHit(hit), thePhi(hit->globalPosition().barePhi()) {}
0025     HitWithPhi(const Hit& hit, float phi) : theHit(hit), thePhi(phi) {}
0026     HitWithPhi(float phi) : theHit(nullptr), thePhi(phi) {}
0027     float phi() const { return thePhi; }
0028     Hit const& hit() const { return theHit; }
0029
0030   private:
0031     Hit theHit;
0032     float thePhi;
0033   };
0034
0035   struct HitLessPhi {
0036     bool operator()(const HitWithPhi& a, const HitWithPhi& b) { return a.phi() < b.phi(); }
0037   };
0038   typedef std::vector<HitWithPhi>::const_iterator HitIter;
0039   typedef std::pair<HitIter, HitIter> Range;
0040
0041   using DoubleRange = std::array<int, 4>;
0042
0043   RecHitsSortedInPhi(const std::vector<Hit>& hits, GlobalPoint const& origin, DetLayer const* il);
0044
0045   bool empty() const { return theHits.empty(); }
0046   std::size_t size() const { return theHits.size(); }
0047
0048   // Returns the hits in the phi range (phi in radians).
0049   //  The phi interval ( phiMin, phiMax) defined as the signed path along the
0050   //  trigonometric circle from the point at phiMin to the point at phiMax
0051   //  must be positive and smaller than pi.
0052   //  At least one of phiMin, phiMax must be in (-pi,pi) range.
0053   //  Examples of correct intervals: (-3,-2), (-4,-3), (3.1,3.2), (3,-3).
0054   //  Examples of WRONG intervals: (-5,-4),(3,2), (4,3), (3.2,3.1), (-3,3), (4,5).
0055   //  Example of use: myHits = recHitsSortedInPhi( phi-deltaPhi, phi+deltaPhi);
0056   //
0057   std::vector<Hit> hits(float phiMin, float phiMax) const;
0058
0059   // Same as above but the result is allocated by the caller and passed by reference.
0060   //  The caller is responsible for clearing of the container "result".
0061   //  This interface is not nice and not safe, but is much faster, since the
0062   //  dominant CPU time of the "nice" method hits(phimin,phimax) is spent in
0063   //  memory allocation of the result!
0064   //
0065   void hits(float phiMin, float phiMax, std::vector<Hit>& result) const;
0066
0067   // some above, just double range of indeces..
0068   DoubleRange doubleRange(float phiMin, float phiMax) const;
0069
0070   // Fast access to the hits in the phi interval (phi in radians).
0071   //  The arguments must satisfy -pi <= phiMin < phiMax <= pi
0072   //  No check is made for this.
0073   //
0074   Range unsafeRange(float phiMin, float phiMax) const;
0075
0076   std::vector<Hit> hits() const {
0077     std::vector<Hit> result;
0078     result.reserve(theHits.size());
0079     for (HitIter i = theHits.begin(); i != theHits.end(); i++)
0080       result.push_back(i->hit());
0081     return result;
0082   }
0083
0084   Range all() const { return Range(theHits.begin(), theHits.end()); }
0085
0086 public:
0087   float phi(int i) const { return theHits[i].phi(); }
0088   float gv(int i) const { return isBarrel ? z[i] : gp(i).perp(); }  // global v
0089   float rv(int i) const { return isBarrel ? u[i] : v[i]; }          // dispaced r
0090   GlobalPoint gp(int i) const { return GlobalPoint(x[i], y[i], z[i]); }
0091
0092 public:
0093   GlobalPoint theOrigin;
0094
0095   std::vector<HitWithPhi> theHits;
0096
0097   DetLayer const* layer;
0098   bool isBarrel;
0099
0100   std::vector<float> x;
0101   std::vector<float> y;
0102   std::vector<float> z;
0103   std::vector<float> drphi;
0104
0105   // barrel: u=r, v=z, forward the opposite...
0106   std::vector<float> u;
0107   std::vector<float> v;
0108   std::vector<float> du;
0109   std::vector<float> dv;
0110   std::vector<float> lphi;
0111
0112   static void copyResult(const Range& range, std::vector<Hit>& result) {
0113     result.reserve(result.size() + (range.second - range.first));
0114     for (HitIter i = range.first; i != range.second; i++)
0115       result.push_back(i->hit());
0116   }
0117 };
0118
0119 /*
0120  *   a collection of hit pairs issued by a doublet search
0121  * replace HitPairs as a communication mean between doublet and triplet search algos
0122  *
0123  */
0124 class HitDoublets {
0125 public:
0126   enum layer { inner = 0, outer = 1 };
0127
0128   using HitLayer = RecHitsSortedInPhi;
0129   using Hit = RecHitsSortedInPhi::Hit;
0130   using ADoublet = std::pair<int, int>;
0131
0132   HitDoublets(RecHitsSortedInPhi const& in, RecHitsSortedInPhi const& out) : layers{{&in, &out}} {}
0133
0134   HitDoublets(HitDoublets&& rh) : layers(rh.layers), indeces(std::move(rh.indeces)) {}
0135
0136   void reserve(std::size_t s) { indeces.reserve(s); }
0137   std::size_t size() const { return indeces.size(); }
0138   bool empty() const { return indeces.empty(); }
0139   void clear() { indeces.clear(); }
0140   void shrink_to_fit() { indeces.shrink_to_fit(); }
0141
0142   void add(int il, int ol) { indeces.emplace_back(il, ol); }
0143
0144   int index(int i, layer l) const { return l == inner ? innerHitId(i) : outerHitId(i); }
0145   DetLayer const* detLayer(layer l) const { return layers[l]->layer; }
0146   HitLayer const& innerLayer() const { return *layers[inner]; }
0147   HitLayer const& outerLayer() const { return *layers[outer]; }
0148   int innerHitId(int i) const { return indeces[i].first; }
0149   int outerHitId(int i) const { return indeces[i].second; }
0150   Hit const& hit(int i, layer l) const { return layers[l]->theHits[index(i, l)].hit(); }
0151   float phi(int i, layer l) const { return layers[l]->phi(index(i, l)); }
0152   float rv(int i, layer l) const { return layers[l]->rv(index(i, l)); }
0153   float r(int i, layer l) const {
0154     float xp = x(i, l);
0155     float yp = y(i, l);
0156     return std::sqrt(xp * xp + yp * yp);
0157   }
0158   float z(int i, layer l) const { return layers[l]->z[index(i, l)]; }
0159   float x(int i, layer l) const { return layers[l]->x[index(i, l)]; }
0160   float y(int i, layer l) const { return layers[l]->y[index(i, l)]; }
0161   GlobalPoint gp(int i, layer l) const { return GlobalPoint(x(i, l), y(i, l), z(i, l)); }
0162
0163 private:
0164   std::array<RecHitsSortedInPhi const*, 2> layers;
0165
0166   std::vector<ADoublet> indeces;  // naturally sorted by outerId
0167 };
0168
0169 #endif