PixelSeeding/src/CACell.h

0001 #ifndef RecoTracker_PixelSeeding_src_CACell_h
0002 #define RecoTracker_PixelSeeding_src_CACell_h
0003
0004 #include <array>
0005 #include <cmath>
0006
0007 #include "DataFormats/Math/interface/deltaPhi.h"
0008 #include "RecoTracker/TkHitPairs/interface/RecHitsSortedInPhi.h"
0009 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
0010 #include "TrackingTools/TransientTrackingRecHit/interface/SeedingLayerSetsHits.h"
0011 #include "RecoTracker/PixelSeeding/interface/CACut.h"
0012
0013 class CACellStatus {
0014 public:
0015   unsigned char getCAState() const { return theCAState; }
0016
0017   // if there is at least one left neighbor with the same state (friend), the state has to be increased by 1.
0018   void updateState() { theCAState += hasSameStateNeighbors; }
0019
0020   bool isRootCell(const unsigned int minimumCAState) const { return (theCAState >= minimumCAState); }
0021
0022 public:
0023   unsigned char theCAState = 0;
0024   unsigned char hasSameStateNeighbors = 0;
0025 };
0026
0027 class CACell {
0028 public:
0029   using Hit = RecHitsSortedInPhi::Hit;
0030   using CAntuple = std::vector<unsigned int>;
0031   using CAntuplet = std::vector<unsigned int>;
0032   using CAColl = std::vector<CACell>;
0033   using CAStatusColl = std::vector<CACellStatus>;
0034
0035   CACell(const HitDoublets* doublets, int doubletId, const int innerHitId, const int outerHitId)
0036       : theDoublets(doublets),
0037         theDoubletId(doubletId),
0038         theInnerR(doublets->rv(doubletId, HitDoublets::inner)),
0039         theInnerZ(doublets->z(doubletId, HitDoublets::inner)) {}
0040
0041   Hit const& getInnerHit() const { return theDoublets->hit(theDoubletId, HitDoublets::inner); }
0042
0043   Hit const& getOuterHit() const { return theDoublets->hit(theDoubletId, HitDoublets::outer); }
0044
0045   int getInnerLayer() const { return theDoublets->detLayer(HitDoublets::inner)->seqNum(); }
0046
0047   int getOuterLayer() const { return theDoublets->detLayer(HitDoublets::outer)->seqNum(); }
0048
0049   float getInnerX() const { return theDoublets->x(theDoubletId, HitDoublets::inner); }
0050
0051   float getOuterX() const { return theDoublets->x(theDoubletId, HitDoublets::outer); }
0052
0053   float getInnerY() const { return theDoublets->y(theDoubletId, HitDoublets::inner); }
0054
0055   float getOuterY() const { return theDoublets->y(theDoubletId, HitDoublets::outer); }
0056
0057   float getInnerZ() const { return theInnerZ; }
0058
0059   float getOuterZ() const { return theDoublets->z(theDoubletId, HitDoublets::outer); }
0060
0061   float getInnerR() const { return theInnerR; }
0062
0063   float getOuterR() const { return theDoublets->rv(theDoubletId, HitDoublets::outer); }
0064
0065   float getInnerPhi() const { return theDoublets->phi(theDoubletId, HitDoublets::inner); }
0066
0067   float getOuterPhi() const { return theDoublets->phi(theDoubletId, HitDoublets::outer); }
0068
0069   void evolve(unsigned int me, CAStatusColl& allStatus) {
0070     allStatus[me].hasSameStateNeighbors = 0;
0071     auto mystate = allStatus[me].theCAState;
0072
0073     for (auto oc : theOuterNeighbors) {
0074       if (allStatus[oc].getCAState() == mystate) {
0075         allStatus[me].hasSameStateNeighbors = 1;
0076
0077         break;
0078       }
0079     }
0080   }
0081
0082   void checkAlignmentAndAct(CAColl& allCells,
0083                             CAntuple& innerCells,
0084                             const float ptmin,
0085                             const float region_origin_x,
0086                             const float region_origin_y,
0087                             const float region_origin_radius,
0088                             const CACut::CAValuesByInnerLayerIds& thetaCutByInnerLayer,
0089                             const CACut::CAValuesByInnerLayerIds& phiCutByInnerLayer,
0090                             const float hardPtCut,
0091                             std::vector<CACell::CAntuplet>* foundTriplets) {
0092     int ncells = innerCells.size();
0093     int constexpr VSIZE = 16;
0094     int ok[VSIZE];
0095     float r1[VSIZE];
0096     float z1[VSIZE];
0097     float thetaCut[VSIZE];
0098     float phiCut[VSIZE];
0099     auto ro = getOuterR();
0100     auto zo = getOuterZ();
0101     unsigned int cellId = this - &allCells.front();
0102     auto loop = [&](int i, int vs) {
0103       for (int j = 0; j < vs; ++j) {
0104         auto koc = innerCells[i + j];
0105         auto& oc = allCells[koc];
0106         r1[j] = oc.getInnerR();
0107         z1[j] = oc.getInnerZ();
0108         thetaCut[j] = thetaCutByInnerLayer.at(oc.getInnerLayer());
0109         phiCut[j] = phiCutByInnerLayer.at(oc.getInnerLayer());
0110       }
0111       // this vectorize!
0112       for (int j = 0; j < vs; ++j)
0113         ok[j] = areAlignedRZ(r1[j], z1[j], ro, zo, ptmin, thetaCut[j]);
0114       for (int j = 0; j < vs; ++j) {
0115         auto koc = innerCells[i + j];
0116         auto& oc = allCells[koc];
0117         if (ok[j] && haveSimilarCurvature(
0118                          oc, ptmin, region_origin_x, region_origin_y, region_origin_radius, phiCut[j], hardPtCut)) {
0119           if (foundTriplets)
0120             foundTriplets->emplace_back(CACell::CAntuplet{koc, cellId});
0121           else {
0122             oc.tagAsOuterNeighbor(cellId);
0123           }
0124         }
0125       }
0126     };
0127     auto lim = VSIZE * (ncells / VSIZE);
0128     for (int i = 0; i < lim; i += VSIZE)
0129       loop(i, VSIZE);
0130     loop(lim, ncells - lim);
0131   }
0132
0133   void checkAlignmentAndTag(CAColl& allCells,
0134                             CAntuple& innerCells,
0135                             const float ptmin,
0136                             const float region_origin_x,
0137                             const float region_origin_y,
0138                             const float region_origin_radius,
0139                             const CACut::CAValuesByInnerLayerIds& thetaCutByInnerLayer,
0140                             const CACut::CAValuesByInnerLayerIds& phiCutByInnerLayer,
0141                             const float hardPtCut) {
0142     checkAlignmentAndAct(allCells,
0143                          innerCells,
0144                          ptmin,
0145                          region_origin_x,
0146                          region_origin_y,
0147                          region_origin_radius,
0148                          thetaCutByInnerLayer,
0149                          phiCutByInnerLayer,
0150                          hardPtCut,
0151                          nullptr);
0152   }
0153   void checkAlignmentAndPushTriplet(CAColl& allCells,
0154                                     CAntuple& innerCells,
0155                                     std::vector<CACell::CAntuplet>& foundTriplets,
0156                                     const float ptmin,
0157                                     const float region_origin_x,
0158                                     const float region_origin_y,
0159                                     const float region_origin_radius,
0160                                     const CACut::CAValuesByInnerLayerIds& thetaCutByInnerLayer,
0161                                     const CACut::CAValuesByInnerLayerIds& phiCutByInnerLayer,
0162                                     const float hardPtCut) {
0163     checkAlignmentAndAct(allCells,
0164                          innerCells,
0165                          ptmin,
0166                          region_origin_x,
0167                          region_origin_y,
0168                          region_origin_radius,
0169                          thetaCutByInnerLayer,
0170                          phiCutByInnerLayer,
0171                          hardPtCut,
0172                          &foundTriplets);
0173   }
0174
0175   int areAlignedRZ(float r1, float z1, float ro, float zo, const float ptmin, const float thetaCut) const {
0176     float radius_diff = std::abs(r1 - ro);
0177     float distance_13_squared = radius_diff * radius_diff + (z1 - zo) * (z1 - zo);
0178
0179     float pMin = ptmin * std::sqrt(distance_13_squared);  //this needs to be divided by radius_diff later
0180
0181     float tan_12_13_half_mul_distance_13_squared =
0182         fabs(z1 * (getInnerR() - ro) + getInnerZ() * (ro - r1) + zo * (r1 - getInnerR()));
0183     return tan_12_13_half_mul_distance_13_squared * pMin <= thetaCut * distance_13_squared * radius_diff;
0184   }
0185
0186   void tagAsOuterNeighbor(unsigned int otherCell) { theOuterNeighbors.push_back(otherCell); }
0187
0188   bool haveSimilarCurvature(const CACell& otherCell,
0189                             const float ptmin,
0190                             const float region_origin_x,
0191                             const float region_origin_y,
0192                             const float region_origin_radius,
0193                             const float phiCut,
0194                             const float hardPtCut) const {
0195     auto x1 = otherCell.getInnerX();
0196     auto y1 = otherCell.getInnerY();
0197
0198     auto x2 = getInnerX();
0199     auto y2 = getInnerY();
0200
0201     auto x3 = getOuterX();
0202     auto y3 = getOuterY();
0203
0204     float distance_13_squared = (x1 - x3) * (x1 - x3) + (y1 - y3) * (y1 - y3);
0205     float tan_12_13_half_mul_distance_13_squared = std::abs(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
0206     // high pt : just straight
0207     if (tan_12_13_half_mul_distance_13_squared * ptmin <= 1.0e-4f * distance_13_squared) {
0208       float distance_3_beamspot_squared =
0209           (x3 - region_origin_x) * (x3 - region_origin_x) + (y3 - region_origin_y) * (y3 - region_origin_y);
0210
0211       float dot_bs3_13 = ((x1 - x3) * (region_origin_x - x3) + (y1 - y3) * (region_origin_y - y3));
0212       float proj_bs3_on_13_squared = dot_bs3_13 * dot_bs3_13 / distance_13_squared;
0213
0214       float distance_13_beamspot_squared = distance_3_beamspot_squared - proj_bs3_on_13_squared;
0215
0216       return distance_13_beamspot_squared < (region_origin_radius + phiCut) * (region_origin_radius + phiCut);
0217     }
0218
0219     //87 cm/GeV = 1/(3.8T * 0.3)
0220
0221     //take less than radius given by the hardPtCut and reject everything below
0222     float minRadius = hardPtCut * 87.f;  // FIXME move out and use real MagField
0223
0224     auto det = (x1 - x2) * (y2 - y3) - (x2 - x3) * (y1 - y2);
0225
0226     auto offset = x2 * x2 + y2 * y2;
0227
0228     auto bc = (x1 * x1 + y1 * y1 - offset) * 0.5f;
0229
0230     auto cd = (offset - x3 * x3 - y3 * y3) * 0.5f;
0231
0232     auto idet = 1.f / det;
0233
0234     auto x_center = (bc * (y2 - y3) - cd * (y1 - y2)) * idet;
0235     auto y_center = (cd * (x1 - x2) - bc * (x2 - x3)) * idet;
0236
0237     auto radius = std::sqrt((x2 - x_center) * (x2 - x_center) + (y2 - y_center) * (y2 - y_center));
0238
0239     if (radius < minRadius)
0240       return false;  // hard cut on pt
0241
0242     auto centers_distance_squared = (x_center - region_origin_x) * (x_center - region_origin_x) +
0243                                     (y_center - region_origin_y) * (y_center - region_origin_y);
0244     auto region_origin_radius_plus_tolerance = region_origin_radius + phiCut;
0245     auto minimumOfIntersectionRange =
0246         (radius - region_origin_radius_plus_tolerance) * (radius - region_origin_radius_plus_tolerance);
0247
0248     if (centers_distance_squared >= minimumOfIntersectionRange) {
0249       auto maximumOfIntersectionRange =
0250           (radius + region_origin_radius_plus_tolerance) * (radius + region_origin_radius_plus_tolerance);
0251       return centers_distance_squared <= maximumOfIntersectionRange;
0252     }
0253
0254     return false;
0255   }
0256
0257   // trying to free the track building process from hardcoded layers, leaving the visit of the graph
0258   // based on the neighborhood connections between cells.
0259
0260   void findNtuplets(CAColl& allCells,
0261                     std::vector<CAntuplet>& foundNtuplets,
0262                     CAntuplet& tmpNtuplet,
0263                     const unsigned int minHitsPerNtuplet) const {
0264     // the building process for a track ends if:
0265     // it has no outer neighbor
0266     // it has no compatible neighbor
0267     // the ntuplets is then saved if the number of hits it contains is greater than a threshold
0268
0269     if (tmpNtuplet.size() == minHitsPerNtuplet - 1) {
0270       foundNtuplets.push_back(tmpNtuplet);
0271     } else {
0272       unsigned int numberOfOuterNeighbors = theOuterNeighbors.size();
0273       for (unsigned int i = 0; i < numberOfOuterNeighbors; ++i) {
0274         tmpNtuplet.push_back((theOuterNeighbors[i]));
0275         allCells[theOuterNeighbors[i]].findNtuplets(allCells, foundNtuplets, tmpNtuplet, minHitsPerNtuplet);
0276         tmpNtuplet.pop_back();
0277       }
0278     }
0279   }
0280
0281 private:
0282   CAntuple theOuterNeighbors;
0283
0284   const HitDoublets* theDoublets;
0285   const int theDoubletId;
0286
0287   const float theInnerR;
0288   const float theInnerZ;
0289 };
0290
0291 #endif  // RecoTracker_PixelSeeding_src_CACell_h