PixelLowPtUtilities/interface/ClusterShapeHitFilter.h

0001 #ifndef RecoTracker_PixelLowPtUtilities_ClusterShapeHitFilter_h
0002 #define RecoTracker_PixelLowPtUtilities_ClusterShapeHitFilter_h
0003
0004 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilter.h"
0005
0006 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0007 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0008
0009 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0010 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0011 #include "DataFormats/SiPixelCluster/interface/SiPixelClusterShapeCache.h"
0012
0013 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0014
0015 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterData.h"
0016
0017 #include <utility>
0018 #include <unordered_map>
0019 #include <cstring>
0020
0021 /*****************************************************************************/
0022 class PixelKeys {
0023 public:
0024   PixelKeys(int part, int dx, int dy) : key((part == 0) ? barrelPacking(dx, dy) : endcapPacking(dx, dy)) {}
0025
0026   operator unsigned int() const { return key; }
0027
0028   static unsigned char endcapPacking(int dx, int dy) {
0029     if (dx < 0 || dy < 0)
0030       return N;
0031     if (dx > 10 || dy > 4)
0032       return N;
0033     return N_barrel + dx * offset_endcap_dy + dy;  // max 11*5 = 55
0034   }
0035
0036   static unsigned char barrelPacking(int dx, int dy) {
0037     if (dx < 0 || dy < 0)
0038       return N;
0039     if (dx > 10 || dy > 15)
0040       return N;
0041     if (dx < 8)
0042       return dx * 16 + dy;  // max 8*16=128
0043     if (dy > 2)
0044       return N;
0045     return 128 + (dx - 8) * 3 + dy;  // max = 128+9 = 137
0046   }
0047
0048   bool isValid() const { return key < N; }
0049
0050   static const int offset_endcap_dy = 5;
0051   static const int offset_endcap_dx = 10;
0052   static const int N_endcap = 55;
0053   static const int N_barrel = 137;
0054   static const int N = N_barrel + N_endcap;
0055
0056   bool operator<(const PixelKeys& right) const { return key < right.key; }
0057
0058 private:
0059   unsigned char key;
0060 };
0061
0062 class StripKeys {
0063 public:
0064   static const int N = 40;
0065
0066   StripKeys(int width) : key(width > 0 ? width - 1 : N) {}
0067
0068   operator unsigned int() const { return key; }
0069
0070   bool isValid() const { return key < N; }
0071
0072   bool operator<(const StripKeys& right) const { return key < right.key; }
0073
0074 private:
0075   unsigned char key;  // max 40;
0076 };
0077
0078 struct PixelLimits {
0079   PixelLimits() {
0080     // init to make sure inside is true;
0081     auto limit = data[0];
0082     limit[0][0] = -10e10;
0083     limit[0][1] = 10e10;
0084     limit[1][0] = -10e10;
0085     limit[1][1] = 10e10;
0086     limit = data[1];
0087     limit[0][0] = -10e10;
0088     limit[0][1] = 10e10;
0089     limit[1][0] = -10e10;
0090     limit[1][1] = 10e10;
0091   }
0092
0093   float data[2][2][2];
0094
0095   bool isInside(const std::pair<float, float>& pred) const {
0096     auto limit = data[0];
0097     bool one = (pred.first > limit[0][0]) && (pred.first < limit[0][1]) && (pred.second > limit[1][0]) &&
0098                (pred.second < limit[1][1]);
0099
0100     limit = data[1];
0101     bool two = (pred.first > limit[0][0]) && (pred.first < limit[0][1]) && (pred.second > limit[1][0]) &&
0102                (pred.second < limit[1][1]);
0103
0104     return one || two;
0105   }
0106 };
0107
0108 struct StripLimits {
0109   StripLimits() {
0110     data[0][0] = -10e10;
0111     data[0][1] = 10e10;
0112     data[1][0] = -10e10;
0113     data[1][1] = 10e10;
0114   }
0115
0116   float data[2][2];
0117
0118   bool isInside(float pred) const {
0119     float const* limit = data[0];
0120     bool one = pred > limit[0] && pred < limit[1];
0121     limit = data[1];
0122     bool two = pred > limit[0] && pred < limit[1];
0123
0124     return one || two;
0125   }
0126 };
0127
0128 /*****************************************************************************/
0129 namespace edm {
0130   class EventSetup;
0131 }
0132
0133 class TrackerGeometry;
0134 class TrackerTopology;
0135 class MagneticField;
0136 class SiPixelLorentzAngle;
0137 class SiStripLorentzAngle;
0138 class PixelGeomDetUnit;
0139 class StripGeomDetUnit;
0140 class StripTopology;
0141
0142 // Function for testing ClusterShapeHitFilter
0143 namespace test {
0144   namespace ClusterShapeHitFilterTest {
0145     int test();
0146   }
0147 }  // namespace test
0148
0149 class ClusterShapeHitFilter {
0150   // For tests
0151   friend int test::ClusterShapeHitFilterTest::test();
0152
0153 public:
0154   struct PixelData {
0155     const PixelGeomDetUnit* det;
0156     unsigned short part;
0157     unsigned short layer;
0158     std::pair<float, float> drift;
0159     std::pair<float, float> cotangent;
0160   };
0161
0162   struct StripData {
0163     const StripGeomDetUnit* det;
0164     StripTopology const* topology;
0165     float drift;
0166     float thickness;
0167     int nstrips;
0168   };
0169
0170   typedef TrajectoryFilter::Record Record;
0171
0172   ClusterShapeHitFilter(const TrackerGeometry* theTracker_,
0173                         const TrackerTopology* theTkTopol_,
0174                         const MagneticField* theMagneticField_,
0175                         const SiPixelLorentzAngle* theSiPixelLorentzAngle_,
0176                         const SiStripLorentzAngle* theSiStripLorentzAngle_,
0177                         const std::string& pixelShapeFile_,
0178                         const std::string& pixelShapeFileL1_);
0179
0180   ~ClusterShapeHitFilter();
0181
0182   void setShapeCuts(bool cutOnPixelShape, bool cutOnStripShape) {
0183     cutOnPixelShape_ = cutOnPixelShape;
0184     cutOnStripShape_ = cutOnStripShape;
0185   }
0186
0187   void setChargeCuts(bool cutOnPixelCharge, float minGoodPixelCharge, bool cutOnStripCharge, float minGoodStripCharge) {
0188     cutOnPixelCharge_ = cutOnPixelCharge;
0189     minGoodPixelCharge_ = minGoodPixelCharge;
0190     cutOnStripCharge_ = cutOnStripCharge;
0191     minGoodStripCharge_ = minGoodStripCharge;
0192   }
0193
0194   bool getSizes(const SiPixelRecHit& recHit,
0195                 const LocalVector& ldir,
0196                 const SiPixelClusterShapeCache& clusterShapeCache,
0197                 int& part,
0198                 ClusterData::ArrayType& meas,
0199                 std::pair<float, float>& predr,
0200                 PixelData const* pd = nullptr) const;
0201   bool isCompatible(const SiPixelRecHit& recHit,
0202                     const LocalVector& ldir,
0203                     const SiPixelClusterShapeCache& clusterShapeCache,
0204                     PixelData const* pd = nullptr) const;
0205   bool isCompatible(const SiPixelRecHit& recHit,
0206                     const GlobalVector& gdir,
0207                     const SiPixelClusterShapeCache& clusterShapeCache,
0208                     PixelData const* pd = nullptr) const;
0209
0210   bool getSizes(DetId detId,
0211                 const SiStripCluster& cluster,
0212                 const LocalPoint& lpos,
0213                 const LocalVector& ldir,
0214                 int& meas,
0215                 float& pred) const;
0216   bool getSizes(
0217       const SiStripRecHit2D& recHit, const LocalPoint& lpos, const LocalVector& ldir, int& meas, float& pred) const {
0218     return getSizes(recHit.geographicalId(), recHit.stripCluster(), lpos, ldir, meas, pred);
0219   }
0220   bool isCompatible(DetId detId, const SiStripCluster& cluster, const LocalPoint& lpos, const LocalVector& ldir) const;
0221   bool isCompatible(DetId detId, const SiStripCluster& cluster, const LocalVector& ldir) const {
0222     return isCompatible(detId, cluster, LocalPoint(0, 0, 0), ldir);
0223   }
0224
0225   bool isCompatible(DetId detId, const SiStripCluster& cluster, const GlobalPoint& gpos, const GlobalVector& gdir) const;
0226   bool isCompatible(DetId detId, const SiStripCluster& cluster, const GlobalVector& gdir) const;
0227
0228   bool isCompatible(const SiStripRecHit2D& recHit, const LocalPoint& lpos, const LocalVector& ldir) const {
0229     return isCompatible(recHit.geographicalId(), recHit.stripCluster(), lpos, ldir);
0230   }
0231   bool isCompatible(const SiStripRecHit2D& recHit, const LocalVector& ldir) const {
0232     return isCompatible(recHit.geographicalId(), recHit.stripCluster(), ldir);
0233   }
0234   bool isCompatible(const SiStripRecHit2D& recHit, const GlobalPoint& gpos, const GlobalVector& gdir) const {
0235     return isCompatible(recHit.geographicalId(), recHit.stripCluster(), gpos, gdir);
0236   }
0237   bool isCompatible(const SiStripRecHit2D& recHit, const GlobalVector& gdir) const {
0238     return isCompatible(recHit.geographicalId(), recHit.stripCluster(), gdir);
0239   }
0240
0241 private:
0242   // for testing purposes only
0243   ClusterShapeHitFilter(std::string const& f1, std::string const& f2) {
0244     loadPixelLimits(f1, pixelLimits);
0245     loadPixelLimits(f2, pixelLimitsL1);
0246     loadStripLimits();
0247   }
0248
0249   const PixelData& getpd(const SiPixelRecHit& recHit, PixelData const* pd = nullptr) const {
0250     if (pd)
0251       return *pd;
0252     // Get detector
0253     DetId id = recHit.geographicalId();
0254     auto p = pixelData.find(id);
0255     return (*p).second;
0256   }
0257
0258   const StripData& getsd(DetId id) const { return stripData.find(id)->second; }
0259
0260   void loadPixelLimits(std::string const& file, PixelLimits* plim);
0261   void loadStripLimits();
0262   void fillPixelData();
0263   void fillStripData();
0264
0265   std::pair<float, float> getCotangent(const PixelGeomDetUnit* pixelDet) const;
0266   float getCotangent(const ClusterShapeHitFilter::StripData& sd, const LocalPoint& p) const;
0267
0268   std::pair<float, float> getDrift(const PixelGeomDetUnit* pixelDet) const;
0269   float getDrift(const StripGeomDetUnit* stripDet) const;
0270
0271   bool isNormalOriented(const GeomDetUnit* geomDet) const;
0272
0273   const TrackerGeometry* theTracker;
0274   const TrackerTopology* theTkTopol;
0275   const MagneticField* theMagneticField;
0276
0277   const SiPixelLorentzAngle* theSiPixelLorentzAngle;
0278   const SiStripLorentzAngle* theSiStripLorentzAngle;
0279
0280   std::unordered_map<unsigned int, PixelData> pixelData;
0281   std::unordered_map<unsigned int, StripData> stripData;
0282
0283   PixelLimits pixelLimits[PixelKeys::N + 1];    // [2][2][2]
0284   PixelLimits pixelLimitsL1[PixelKeys::N + 1];  // for BPIX1
0285
0286   StripLimits stripLimits[StripKeys::N + 1];  // [2][2]
0287
0288   float theAngle[6];
0289   bool cutOnPixelCharge_, cutOnStripCharge_;
0290   float minGoodPixelCharge_, minGoodStripCharge_;
0291   bool cutOnPixelShape_, cutOnStripShape_;
0292   bool checkClusterCharge(DetId detId, const SiStripCluster& cluster, const LocalVector& ldir) const;
0293   bool checkClusterCharge(DetId detId, const SiPixelCluster& cluster, const LocalVector& ldir) const;
0294 };
0295
0296 #endif