RecHitsSortedInPhi.h

Go to the documentation of this file.

#ifndef RecHitsSortedInPhi_H
#define RecHitsSortedInPhi_H

#include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"

#include <vector>
#include<array>

#include<cassert>

class RecHitsSortedInPhi {
public:

  typedef BaseTrackerRecHit const * Hit;

  // A RecHit extension that caches the phi angle for fast access
  class HitWithPhi {
  public:
    HitWithPhi( const Hit & hit) : theHit(hit), thePhi(hit->globalPosition().barePhi()) {}
    HitWithPhi( const Hit & hit,float phi) : theHit(hit), thePhi(phi) {}
    HitWithPhi( float phi) : theHit(nullptr), thePhi(phi) {}
    float phi() const {return thePhi;}
    Hit const & hit() const { return theHit;}
  private:
    Hit   theHit;
    float thePhi;
  };

  struct HitLessPhi {
    bool operator()( const HitWithPhi& a, const HitWithPhi& b) { return a.phi() < b.phi(); }
  };
  typedef std::vector<HitWithPhi>::const_iterator      HitIter;
  typedef std::pair<HitIter,HitIter>            Range;

  using DoubleRange = std::array<int,4>;
  
  RecHitsSortedInPhi(const std::vector<Hit>& hits, GlobalPoint const & origin, DetLayer const * il);

  bool empty() const { return theHits.empty(); }
  std::size_t size() const { return theHits.size();}


  // Returns the hits in the phi range (phi in radians).
  //  The phi interval ( phiMin, phiMax) defined as the signed path along the 
  //  trigonometric circle from the point at phiMin to the point at phiMax
  //  must be positive and smaller than pi.
  //  At least one of phiMin, phiMax must be in (-pi,pi) range.
  //  Examples of correct intervals: (-3,-2), (-4,-3), (3.1,3.2), (3,-3).
  //  Examples of WRONG intervals: (-5,-4),(3,2), (4,3), (3.2,3.1), (-3,3), (4,5).
  //  Example of use: myHits = recHitsSortedInPhi( phi-deltaPhi, phi+deltaPhi);
  //
  std::vector<Hit> hits( float phiMin, float phiMax) const;

  // Same as above but the result is allocated by the caller and passed by reference.
  //  The caller is responsible for clearing of the container "result".
  //  This interface is not nice and not safe, but is much faster, since the
  //  dominant CPU time of the "nice" method hits(phimin,phimax) is spent in
  //  memory allocation of the result!
  //
  void hits( float phiMin, float phiMax, std::vector<Hit>& result) const;

  // some above, just double range of indeces..
  DoubleRange doubleRange(float phiMin, float phiMax) const;

  // Fast access to the hits in the phi interval (phi in radians).
  //  The arguments must satisfy -pi <= phiMin < phiMax <= pi
  //  No check is made for this.
  //
  Range unsafeRange( float phiMin, float phiMax) const;

  std::vector<Hit> hits() const {
    std::vector<Hit> result; result.reserve(theHits.size());
    for (HitIter i=theHits.begin(); i!=theHits.end(); i++) result.push_back(i->hit());
    return result;
  }


  Range all() const {
    return Range(theHits.begin(), theHits.end());
  }

public:
  float       phi(int i) const { return theHits[i].phi();}
  float       gv(int i) const { return isBarrel ? z[i] : gp(i).perp();}  // global v
  float       rv(int i) const { return isBarrel ? u[i] : v[i];}  // dispaced r
  GlobalPoint gp(int i) const { return GlobalPoint(x[i],y[i],z[i]);}

public:

  GlobalPoint theOrigin;

  std::vector<HitWithPhi> theHits;

  DetLayer const * layer;
  bool isBarrel;

  std::vector<float> x;
  std::vector<float> y;
  std::vector<float> z;
  std::vector<float> drphi;

  // barrel: u=r, v=z, forward the opposite...
  std::vector<float> u;
  std::vector<float> v;
  std::vector<float> du;
  std::vector<float> dv;
  std::vector<float> lphi;

  static void copyResult( const Range& range, std::vector<Hit>& result) {
    result.reserve(result.size()+(range.second-range.first));
    for (HitIter i = range.first; i != range.second; i++) result.push_back( i->hit());
  }

};



/*
 *   a collection of hit pairs issued by a doublet search
 * replace HitPairs as a communication mean between doublet and triplet search algos
 *
 */
class HitDoublets {
public:
  enum layer { inner=0, outer=1};

  using HitLayer = RecHitsSortedInPhi;
  using Hit=RecHitsSortedInPhi::Hit;
  using ADoublet = std::pair<int,int>;
  
  HitDoublets(  RecHitsSortedInPhi const & in,
        RecHitsSortedInPhi const & out) :
    layers{{&in,&out}}{}
  
  HitDoublets(HitDoublets && rh) : layers(std::move(rh.layers)), indeces(std::move(rh.indeces)){}
  
  void reserve(std::size_t s) { indeces.reserve(s);}
  std::size_t size() const { return indeces.size();}
  bool empty() const { return indeces.empty();}
  void clear() { indeces.clear();}
  void shrink_to_fit() {
    indeces.shrink_to_fit();
  }
  
  void add (int il, int ol) {
    indeces.emplace_back(il,ol);
  }

  int index(int i, layer l) const { return l==inner ? innerHitId(i) : outerHitId(i);}
  DetLayer const * detLayer(layer l) const { return layers[l]->layer; }
  HitLayer const & innerLayer() const { return *layers[inner];}
  HitLayer const & outerLayer() const { return *layers[outer];}
  int innerHitId(int i) const {return indeces[i].first;}
  int outerHitId(int i) const {return indeces[i].second;}
  Hit const & hit(int i, layer l) const { return layers[l]->theHits[index(i,l)].hit();}
  float       phi(int i, layer l) const { return layers[l]->phi(index(i,l));}
  float       rv(int i, layer l) const { return layers[l]->rv(index(i,l));}
  float       r(int i, layer l) const { float xp = x(i,l); float yp = y(i,l);  return std::sqrt (xp*xp + yp*yp);}
  float        z(int i, layer l) const { return layers[l]->z[index(i,l)];}
  float        x(int i, layer l) const { return layers[l]->x[index(i,l)];}
  float        y(int i, layer l) const { return layers[l]->y[index(i,l)];}
  GlobalPoint gp(int i, layer l) const { return GlobalPoint(x(i,l),y(i,l),z(i,l));}

private:

  std::array<RecHitsSortedInPhi const *,2> layers;


  std::vector<ADoublet> indeces; // naturally sorted by outerId

};

#endif