PixelLowPtUtilities/src/ClusterShape.cc

0001 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterShape.h"
0002 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterData.h"
0003
0004 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0005 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0006
0007 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0008
0009 #include <vector>
0010 #include <fstream>
0011
0012 using namespace std;
0013
0014 /*****************************************************************************/
0015 ClusterShape::ClusterShape() {}
0016
0017 /*****************************************************************************/
0018 ClusterShape::~ClusterShape() {}
0019
0020 /*****************************************************************************/
0021 int ClusterShape::getDirection(int low, int hig, int olow, int ohig) {
0022   if (hig == ohig && low == olow)
0023     return 0;
0024   if (hig >= ohig && low >= olow)
0025     return 1;
0026   if (hig <= ohig && low <= olow)
0027     return -1;
0028
0029   return -2;
0030 }
0031
0032 /*****************************************************************************/
0033 bool ClusterShape::processColumn(pair<int, int> pos, bool inTheLoop) {
0034   if (x[1] > -1) {  // Process previous column
0035     if (low < y[0] || x[1] == x[0])
0036       y[0] = low;
0037     if (hig > y[1] || x[1] == x[0])
0038       y[1] = hig;
0039
0040     if (x[1] > x[0]) {  // Determine direction
0041       int dir = getDirection(low, hig, olow, ohig);
0042
0043       // no direction
0044       if (dir == -2)
0045         return false;
0046
0047       if (x[1] > x[0] + 1) {  // Check if direction changes
0048         if (odir * dir == -1) {
0049           odir = -2;
0050           return false;
0051         }
0052       }
0053
0054       if (x[1] <= x[0] + 1 || odir == 0)
0055         odir = dir;
0056     }
0057
0058     olow = low;
0059     ohig = hig;
0060   } else {  // Very first column, initialize
0061     x[0] = pos.first;
0062   }
0063
0064   // Open new column
0065   if (inTheLoop) {
0066     x[1] = pos.first;
0067     low = pos.second;
0068     hig = pos.second;
0069   }
0070
0071   return (true);
0072 }
0073
0074 /*****************************************************************************/
0075 void ClusterShape::determineShape(const PixelGeomDetUnit& pixelDet, const SiPixelRecHit& recHit, ClusterData& data) {
0076   determineShape(pixelDet, *(recHit.cluster()), data);
0077 }
0078
0079 void ClusterShape::determineShape(const PixelGeomDetUnit& pixelDet, const SiPixelCluster& cluster, ClusterData& data) {
0080   // Topology
0081   const PixelTopology* theTopology = (&(pixelDet.specificTopology()));
0082
0083   // Initialize
0084   data.isStraight = true;
0085   data.isComplete = true;
0086
0087   x[0] = -1;
0088   x[1] = -1;
0089   y[0] = -1;
0090   y[1] = -1;
0091   olow = -2;
0092   ohig = -2;
0093   odir = 0;
0094   low = 0;
0095   hig = 0;
0096
0097   pair<int, int> pos;
0098
0099   // Get sorted pixels
0100   size_t npixels = cluster.pixelADC().size();
0101   pixels_.reserve(npixels);
0102   for (size_t i = 0; i < npixels; ++i) {
0103     pixels_.push_back(cluster.pixel(i));
0104   }
0105   sort(pixels_.begin(), pixels_.end(), [](auto const& a, auto const& b) {
0106     // slightly faster by avoiding branches
0107     return (a.x < b.x) | ((a.x == b.x) & (a.y < b.y));
0108   });
0109
0110   // Look at all the pixels
0111   for (const auto& pixel : pixels_) {
0112     // Position
0113     pos.first = (int)pixel.x;
0114     pos.second = (int)pixel.y;
0115
0116     // Check if at the edge or big
0117     if (theTopology->isItEdgePixelInX(pos.first) || theTopology->isItEdgePixelInY(pos.second)) {
0118       data.isComplete = false;
0119     }  // break; }
0120
0121     // Check if straight
0122     if (pos.first > x[1]) {  // column ready
0123       if (processColumn(pos, true) == false) {
0124         data.isStraight = false;
0125       }  // break; }
0126     } else {                     // increasing column
0127       if (pos.second > hig + 1)  // at least a pixel is missing
0128       {
0129         data.isStraight = false;
0130       }  // break; }
0131
0132       hig = pos.second;
0133     }
0134   }
0135   pixels_.clear();
0136
0137   // Check if straight, process last column
0138   if (processColumn(pos, false) == false)
0139     data.isStraight = false;
0140
0141   // Treat clusters with big pixel(s) inside
0142   const int minPixelRow = cluster.minPixelRow();
0143   const int maxPixelRow = cluster.maxPixelRow();
0144   for (int ix = minPixelRow + 1; ix < maxPixelRow; ix++)
0145     x[1] += theTopology->isItBigPixelInX(ix);
0146
0147   const int minPixelCol = cluster.minPixelCol();
0148   const int maxPixelCol = cluster.maxPixelCol();
0149   for (int iy = minPixelCol + 1; iy < maxPixelCol; iy++)
0150     y[1] += theTopology->isItBigPixelInY(iy);
0151
0152   // Treat clusters with bix pixel(s) outside, FIXME FIXME
0153   unsigned int px = 0;
0154   px += theTopology->isItBigPixelInX(minPixelRow);
0155   px += theTopology->isItBigPixelInX(maxPixelRow);
0156
0157   unsigned int py = 0;
0158   py += theTopology->isItBigPixelInY(minPixelCol);
0159   py += theTopology->isItBigPixelInY(maxPixelCol);
0160
0161   data.hasBigPixelsOnlyInside = (px <= 0 && py <= 0);
0162
0163   //if( (px > 0 || py > 0) && odir == 0)
0164   if (!data.hasBigPixelsOnlyInside && odir == 0) {
0165     // if outside and don't know the direction FIXME?
0166     data.isComplete = false;
0167   }
0168   // else
0169   {  // FIXME do it
0170     assert((px + 1) * (py + 1) <= data.size.capacity());
0171     const int pre_dx = x[1] - x[0];
0172     const int pre_dy = y[1] - y[0];
0173     for (unsigned int ax = 0; ax <= px; ax++)
0174       for (unsigned int ay = 0; ay <= py; ay++) {
0175         int dx = pre_dx + ax;
0176         int dy = pre_dy + ay;
0177         if (odir != 0)
0178           dy *= odir;
0179
0180         pair<int, int> s(dx, dy);
0181         data.size.push_back_unchecked(s);
0182       }
0183   }
0184 }