SiPixelRecHits/interface/pixelCPEforGPU.h

0001 #ifndef RecoLocalTracker_SiPixelRecHits_pixelCPEforGPU_h
0002 #define RecoLocalTracker_SiPixelRecHits_pixelCPEforGPU_h
0003
0004 #include <cassert>
0005 #include <cmath>
0006 #include <cstdint>
0007 #include <iterator>
0008
0009 #include "CUDADataFormats/SiPixelCluster/interface/gpuClusteringConstants.h"
0010 #include "DataFormats/GeometrySurface/interface/SOARotation.h"
0011 #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h"
0012 #include "HeterogeneousCore/CUDAUtilities/interface/cudaCompat.h"
0013 #include "CUDADataFormats/TrackingRecHit/interface/SiPixelHitStatus.h"
0014
0015 namespace CPEFastParametrisation {
0016   // From https://cmssdt.cern.ch/dxr/CMSSW/source/CondFormats/SiPixelTransient/src/SiPixelGenError.cc#485-486
0017   // qbin: int (0-4) describing the charge of the cluster
0018   // [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin]
0019   constexpr int kGenErrorQBins = 5;
0020   // arbitrary number of bins for sampling errors
0021   constexpr int kNumErrorBins = 16;
0022 }  // namespace CPEFastParametrisation
0023
0024 namespace pixelCPEforGPU {
0025
0026   using Status = SiPixelHitStatus;
0027   using Frame = SOAFrame<float>;
0028   using Rotation = SOARotation<float>;
0029
0030   // all modules are identical!
0031   struct CommonParams {
0032     float theThicknessB;
0033     float theThicknessE;
0034     float thePitchX;
0035     float thePitchY;
0036
0037     uint16_t maxModuleStride;
0038     uint8_t numberOfLaddersInBarrel;
0039   };
0040
0041   struct DetParams {
0042     bool isBarrel;
0043     bool isPosZ;
0044     uint16_t layer;
0045     uint16_t index;
0046     uint32_t rawId;
0047
0048     float shiftX;
0049     float shiftY;
0050     float chargeWidthX;
0051     float chargeWidthY;
0052     uint16_t pixmx;  // max pix charge
0053
0054     uint16_t nRowsRoc;  //we don't need 2^16 columns, is worth to use 15 + 1 for sign
0055     uint16_t nColsRoc;
0056     uint16_t nRows;
0057     uint16_t nCols;
0058
0059     uint32_t numPixsInModule;
0060
0061     float x0, y0, z0;  // the vertex in the local coord of the detector
0062
0063     float apeXX, apeYY;  // ape^2
0064     uint8_t sx2, sy1, sy2;
0065     uint8_t sigmax[CPEFastParametrisation::kNumErrorBins], sigmax1[CPEFastParametrisation::kNumErrorBins],
0066         sigmay[CPEFastParametrisation::kNumErrorBins];  // in micron
0067     float xfact[CPEFastParametrisation::kGenErrorQBins], yfact[CPEFastParametrisation::kGenErrorQBins];
0068     int minCh[CPEFastParametrisation::kGenErrorQBins];
0069
0070     Frame frame;
0071   };
0072
0073   template <typename TrackerTopology>
0074   struct LayerGeometryT {
0075     uint32_t layerStart[TrackerTopology::numberOfLayers + 1];
0076     uint8_t layer[pixelTopology::layerIndexSize<TrackerTopology>];
0077     uint16_t maxModuleStride;
0078   };
0079
0080   // using LayerGeometry = LayerGeometryT<pixelTopology::Phase1>;
0081   // using LayerGeometryPhase2 = LayerGeometryT<pixelTopology::Phase2>;
0082
0083   template <typename TrackerTopology>
0084   struct ParamsOnGPUT {
0085     using LayerGeometry = LayerGeometryT<TrackerTopology>;
0086     using AverageGeometry = pixelTopology::AverageGeometryT<TrackerTopology>;
0087
0088     CommonParams const* m_commonParams;
0089     DetParams const* m_detParams;
0090     LayerGeometry const* m_layerGeometry;
0091     AverageGeometry const* m_averageGeometry;
0092
0093     constexpr CommonParams const& __restrict__ commonParams() const {
0094       CommonParams const* __restrict__ l = m_commonParams;
0095       return *l;
0096     }
0097     constexpr DetParams const& __restrict__ detParams(int i) const {
0098       DetParams const* __restrict__ l = m_detParams;
0099       return l[i];
0100     }
0101     constexpr LayerGeometry const& __restrict__ layerGeometry() const { return *m_layerGeometry; }
0102     constexpr AverageGeometry const& __restrict__ averageGeometry() const { return *m_averageGeometry; }
0103
0104     __device__ uint8_t layer(uint16_t id) const {
0105       return __ldg(m_layerGeometry->layer + id / m_layerGeometry->maxModuleStride);
0106     };
0107   };
0108
0109   // SOA (on device)
0110   template <uint32_t N>
0111   struct ClusParamsT {
0112     uint32_t minRow[N];
0113     uint32_t maxRow[N];
0114     uint32_t minCol[N];
0115     uint32_t maxCol[N];
0116
0117     int32_t q_f_X[N];
0118     int32_t q_l_X[N];
0119     int32_t q_f_Y[N];
0120     int32_t q_l_Y[N];
0121
0122     int32_t charge[N];
0123
0124     float xpos[N];
0125     float ypos[N];
0126
0127     float xerr[N];
0128     float yerr[N];
0129
0130     int16_t xsize[N];  // (*8) clipped at 127 if negative is edge....
0131     int16_t ysize[N];
0132
0133     Status status[N];
0134   };
0135
0136   constexpr int32_t MaxHitsInIter = gpuClustering::maxHitsInIter();
0137   using ClusParams = ClusParamsT<MaxHitsInIter>;
0138
0139   constexpr inline void computeAnglesFromDet(
0140       DetParams const& __restrict__ detParams, float const x, float const y, float& cotalpha, float& cotbeta) {
0141     // x,y local position on det
0142     auto gvx = x - detParams.x0;
0143     auto gvy = y - detParams.y0;
0144     auto gvz = -1.f / detParams.z0;
0145     // normalization not required as only ratio used...
0146     // calculate angles
0147     cotalpha = gvx * gvz;
0148     cotbeta = gvy * gvz;
0149   }
0150
0151   constexpr inline float correction(int sizeM1,
0152                                     int q_f,                        //!< Charge in the first pixel.
0153                                     int q_l,                        //!< Charge in the last pixel.
0154                                     uint16_t upper_edge_first_pix,  //!< As the name says.
0155                                     uint16_t lower_edge_last_pix,   //!< As the name says.
0156                                     float lorentz_shift,            //!< L-shift at half thickness
0157                                     float theThickness,             //detector thickness
0158                                     float cot_angle,                //!< cot of alpha_ or beta_
0159                                     float pitch,                    //!< thePitchX or thePitchY
0160                                     bool first_is_big,              //!< true if the first is big
0161                                     bool last_is_big)               //!< true if the last is big
0162   {
0163     if (0 == sizeM1)  // size 1
0164       return 0;
0165
0166     float w_eff = 0;
0167     bool simple = true;
0168     if (1 == sizeM1) {  // size 2
0169       //--- Width of the clusters minus the edge (first and last) pixels.
0170       //--- In the note, they are denoted x_F and x_L (and y_F and y_L)
0171       // assert(lower_edge_last_pix >= upper_edge_first_pix);
0172       auto w_inner = pitch * float(lower_edge_last_pix - upper_edge_first_pix);  // in cm
0173
0174       //--- Predicted charge width from geometry
0175       auto w_pred = theThickness * cot_angle  // geometric correction (in cm)
0176                     - lorentz_shift;          // (in cm) &&& check fpix!
0177
0178       w_eff = std::abs(w_pred) - w_inner;
0179
0180       //--- If the observed charge width is inconsistent with the expectations
0181       //--- based on the track, do *not* use w_pred-w_inner.  Instead, replace
0182       //--- it with an *average* effective charge width, which is the average
0183       //--- length of the edge pixels.
0184
0185       // this can produce "large" regressions for very small numeric differences
0186       simple = (w_eff < 0.0f) | (w_eff > pitch);
0187     }
0188
0189     if (simple) {
0190       //--- Total length of the two edge pixels (first+last)
0191       float sum_of_edge = 2.0f;
0192       if (first_is_big)
0193         sum_of_edge += 1.0f;
0194       if (last_is_big)
0195         sum_of_edge += 1.0f;
0196       w_eff = pitch * 0.5f * sum_of_edge;  // ave. length of edge pixels (first+last) (cm)
0197     }
0198
0199     //--- Finally, compute the position in this projection
0200     float qdiff = q_l - q_f;
0201     float qsum = q_l + q_f;
0202
0203     //--- Temporary fix for clusters with both first and last pixel with charge = 0
0204     if (qsum == 0)
0205       qsum = 1.0f;
0206
0207     return 0.5f * (qdiff / qsum) * w_eff;
0208   }
0209
0210   template <typename TrackerTraits>
0211   constexpr inline void position(CommonParams const& __restrict__ comParams,
0212                                  DetParams const& __restrict__ detParams,
0213                                  ClusParams& cp,
0214                                  uint32_t ic) {
0215     constexpr int maxSize = TrackerTraits::maxSizeCluster;
0216     //--- Upper Right corner of Lower Left pixel -- in measurement frame
0217     uint16_t llx = cp.minRow[ic] + 1;
0218     uint16_t lly = cp.minCol[ic] + 1;
0219
0220     //--- Lower Left corner of Upper Right pixel -- in measurement frame
0221     uint16_t urx = cp.maxRow[ic];
0222     uint16_t ury = cp.maxCol[ic];
0223
0224     uint16_t llxl = llx, llyl = lly, urxl = urx, uryl = ury;
0225
0226     llxl = TrackerTraits::localX(llx);
0227     llyl = TrackerTraits::localY(lly);
0228     urxl = TrackerTraits::localX(urx);
0229     uryl = TrackerTraits::localY(ury);
0230
0231     auto mx = llxl + urxl;
0232     auto my = llyl + uryl;
0233
0234     int xsize = int(urxl) + 2 - int(llxl);
0235     int ysize = int(uryl) + 2 - int(llyl);
0236     assert(xsize >= 0);  // 0 if bixpix...
0237     assert(ysize >= 0);
0238
0239     if (TrackerTraits::isBigPixX(cp.minRow[ic]))
0240       ++xsize;
0241     if (TrackerTraits::isBigPixX(cp.maxRow[ic]))
0242       ++xsize;
0243     if (TrackerTraits::isBigPixY(cp.minCol[ic]))
0244       ++ysize;
0245     if (TrackerTraits::isBigPixY(cp.maxCol[ic]))
0246       ++ysize;
0247
0248     int unbalanceX = 8.f * std::abs(float(cp.q_f_X[ic] - cp.q_l_X[ic])) / float(cp.q_f_X[ic] + cp.q_l_X[ic]);
0249     int unbalanceY = 8.f * std::abs(float(cp.q_f_Y[ic] - cp.q_l_Y[ic])) / float(cp.q_f_Y[ic] + cp.q_l_Y[ic]);
0250
0251     xsize = 8 * xsize - unbalanceX;
0252     ysize = 8 * ysize - unbalanceY;
0253
0254     cp.xsize[ic] = std::min(xsize, maxSize);
0255     cp.ysize[ic] = std::min(ysize, maxSize);
0256
0257     if (cp.minRow[ic] == 0 || cp.maxRow[ic] == uint32_t(detParams.nRows - 1))
0258       cp.xsize[ic] = -cp.xsize[ic];
0259
0260     if (cp.minCol[ic] == 0 || cp.maxCol[ic] == uint32_t(detParams.nCols - 1))
0261       cp.ysize[ic] = -cp.ysize[ic];
0262
0263     // apply the lorentz offset correction
0264     float xoff = 0.5f * float(detParams.nRows) * comParams.thePitchX;
0265     float yoff = 0.5f * float(detParams.nCols) * comParams.thePitchY;
0266
0267     //correction for bigpixels for phase1
0268     xoff = xoff + TrackerTraits::bigPixXCorrection * comParams.thePitchX;
0269     yoff = yoff + TrackerTraits::bigPixYCorrection * comParams.thePitchY;
0270
0271     // apply the lorentz offset correction
0272     auto xPos = detParams.shiftX + (comParams.thePitchX * 0.5f * float(mx)) - xoff;
0273     auto yPos = detParams.shiftY + (comParams.thePitchY * 0.5f * float(my)) - yoff;
0274
0275     float cotalpha = 0, cotbeta = 0;
0276
0277     computeAnglesFromDet(detParams, xPos, yPos, cotalpha, cotbeta);
0278
0279     auto thickness = detParams.isBarrel ? comParams.theThicknessB : comParams.theThicknessE;
0280
0281     auto xcorr = correction(cp.maxRow[ic] - cp.minRow[ic],
0282                             cp.q_f_X[ic],
0283                             cp.q_l_X[ic],
0284                             llxl,
0285                             urxl,
0286                             detParams.chargeWidthX,  // lorentz shift in cm
0287                             thickness,
0288                             cotalpha,
0289                             comParams.thePitchX,
0290                             TrackerTraits::isBigPixX(cp.minRow[ic]),
0291                             TrackerTraits::isBigPixX(cp.maxRow[ic]));
0292
0293     auto ycorr = correction(cp.maxCol[ic] - cp.minCol[ic],
0294                             cp.q_f_Y[ic],
0295                             cp.q_l_Y[ic],
0296                             llyl,
0297                             uryl,
0298                             detParams.chargeWidthY,  // lorentz shift in cm
0299                             thickness,
0300                             cotbeta,
0301                             comParams.thePitchY,
0302                             TrackerTraits::isBigPixY(cp.minCol[ic]),
0303                             TrackerTraits::isBigPixY(cp.maxCol[ic]));
0304
0305     cp.xpos[ic] = xPos + xcorr;
0306     cp.ypos[ic] = yPos + ycorr;
0307   }
0308
0309   template <typename TrackerTraits>
0310   constexpr inline void errorFromSize(CommonParams const& __restrict__ comParams,
0311                                       DetParams const& __restrict__ detParams,
0312                                       ClusParams& cp,
0313                                       uint32_t ic) {
0314     // Edge cluster errors
0315     cp.xerr[ic] = 0.0050;
0316     cp.yerr[ic] = 0.0085;
0317
0318     // FIXME these are errors form Run1
0319     float xerr_barrel_l1_def = TrackerTraits::xerr_barrel_l1_def;
0320     float yerr_barrel_l1_def = TrackerTraits::yerr_barrel_l1_def;
0321     float xerr_barrel_ln_def = TrackerTraits::xerr_barrel_ln_def;
0322     float yerr_barrel_ln_def = TrackerTraits::yerr_barrel_ln_def;
0323     float xerr_endcap_def = TrackerTraits::xerr_endcap_def;
0324     float yerr_endcap_def = TrackerTraits::yerr_endcap_def;
0325
0326     constexpr float xerr_barrel_l1[] = {0.00115, 0.00120, 0.00088};  //TODO MOVE THESE SOMEWHERE ELSE
0327     constexpr float yerr_barrel_l1[] = {
0328         0.00375, 0.00230, 0.00250, 0.00250, 0.00230, 0.00230, 0.00210, 0.00210, 0.00240};
0329     constexpr float xerr_barrel_ln[] = {0.00115, 0.00120, 0.00088};
0330     constexpr float yerr_barrel_ln[] = {
0331         0.00375, 0.00230, 0.00250, 0.00250, 0.00230, 0.00230, 0.00210, 0.00210, 0.00240};
0332     constexpr float xerr_endcap[] = {0.0020, 0.0020};
0333     constexpr float yerr_endcap[] = {0.00210};
0334
0335     auto sx = cp.maxRow[ic] - cp.minRow[ic];
0336     auto sy = cp.maxCol[ic] - cp.minCol[ic];
0337
0338     // is edgy ?
0339     bool isEdgeX = cp.xsize[ic] < 1;
0340     bool isEdgeY = cp.ysize[ic] < 1;
0341
0342     // is one and big?
0343     bool isBig1X = ((0 == sx) && TrackerTraits::isBigPixX(cp.minRow[ic]));
0344     bool isBig1Y = ((0 == sy) && TrackerTraits::isBigPixY(cp.minCol[ic]));
0345
0346     if (!isEdgeX && !isBig1X) {
0347       if (not detParams.isBarrel) {
0348         cp.xerr[ic] = sx < std::size(xerr_endcap) ? xerr_endcap[sx] : xerr_endcap_def;
0349       } else if (detParams.layer == 1) {
0350         cp.xerr[ic] = sx < std::size(xerr_barrel_l1) ? xerr_barrel_l1[sx] : xerr_barrel_l1_def;
0351       } else {
0352         cp.xerr[ic] = sx < std::size(xerr_barrel_ln) ? xerr_barrel_ln[sx] : xerr_barrel_ln_def;
0353       }
0354     }
0355
0356     if (!isEdgeY && !isBig1Y) {
0357       if (not detParams.isBarrel) {
0358         cp.yerr[ic] = sy < std::size(yerr_endcap) ? yerr_endcap[sy] : yerr_endcap_def;
0359       } else if (detParams.layer == 1) {
0360         cp.yerr[ic] = sy < std::size(yerr_barrel_l1) ? yerr_barrel_l1[sy] : yerr_barrel_l1_def;
0361       } else {
0362         cp.yerr[ic] = sy < std::size(yerr_barrel_ln) ? yerr_barrel_ln[sy] : yerr_barrel_ln_def;
0363       }
0364     }
0365   }
0366
0367   template <typename TrackerTraits>
0368   constexpr inline void errorFromDB(CommonParams const& __restrict__ comParams,
0369                                     DetParams const& __restrict__ detParams,
0370                                     ClusParams& cp,
0371                                     uint32_t ic) {
0372     // Edge cluster errors
0373     cp.xerr[ic] = 0.0050f;
0374     cp.yerr[ic] = 0.0085f;
0375
0376     auto sx = cp.maxRow[ic] - cp.minRow[ic];
0377     auto sy = cp.maxCol[ic] - cp.minCol[ic];
0378
0379     // is edgy ?  (size is set negative: see above)
0380     bool isEdgeX = cp.xsize[ic] < 1;
0381     bool isEdgeY = cp.ysize[ic] < 1;
0382     // is one and big?
0383     bool isOneX = (0 == sx);
0384     bool isOneY = (0 == sy);
0385     bool isBigX = TrackerTraits::isBigPixX(cp.minRow[ic]);
0386     bool isBigY = TrackerTraits::isBigPixY(cp.minCol[ic]);
0387
0388     auto ch = cp.charge[ic];
0389     auto bin = 0;
0390     for (; bin < CPEFastParametrisation::kGenErrorQBins - 1; ++bin)
0391       // find first bin which minimum charge exceeds cluster charge
0392       if (ch < detParams.minCh[bin + 1])
0393         break;
0394
0395     // in detParams qBins are reversed bin0 -> smallest charge, bin4-> largest charge
0396     // whereas in CondFormats/SiPixelTransient/src/SiPixelGenError.cc it is the opposite
0397     // so we reverse the bin here -> kGenErrorQBins - 1 - bin
0398     cp.status[ic].qBin = CPEFastParametrisation::kGenErrorQBins - 1 - bin;
0399     cp.status[ic].isOneX = isOneX;
0400     cp.status[ic].isBigX = (isOneX & isBigX) | isEdgeX;
0401     cp.status[ic].isOneY = isOneY;
0402     cp.status[ic].isBigY = (isOneY & isBigY) | isEdgeY;
0403
0404     auto xoff = -float(TrackerTraits::xOffset) * comParams.thePitchX;
0405     int low_value = 0;
0406     int high_value = CPEFastParametrisation::kNumErrorBins - 1;
0407     int bin_value = float(CPEFastParametrisation::kNumErrorBins) * (cp.xpos[ic] + xoff) / (2 * xoff);
0408     // return estimated bin value truncated to [0, 15]
0409     int jx = std::clamp(bin_value, low_value, high_value);
0410
0411     auto toCM = [](uint8_t x) { return float(x) * 1.e-4f; };
0412
0413     if (not isEdgeX) {
0414       cp.xerr[ic] = isOneX ? toCM(isBigX ? detParams.sx2 : detParams.sigmax1[jx])
0415                            : detParams.xfact[bin] * toCM(detParams.sigmax[jx]);
0416     }
0417
0418     auto ey = cp.ysize[ic] > 8 ? detParams.sigmay[std::min(cp.ysize[ic] - 9, 15)] : detParams.sy1;
0419     if (not isEdgeY) {
0420       cp.yerr[ic] = isOneY ? toCM(isBigY ? detParams.sy2 : detParams.sy1) : detParams.yfact[bin] * toCM(ey);
0421     }
0422   }
0423
0424   //for Phase2 -> fallback to error from size
0425   template <>
0426   constexpr inline void errorFromDB<pixelTopology::Phase2>(CommonParams const& __restrict__ comParams,
0427                                                            DetParams const& __restrict__ detParams,
0428                                                            ClusParams& cp,
0429                                                            uint32_t ic) {
0430     errorFromSize<pixelTopology::Phase2>(comParams, detParams, cp, ic);
0431   }
0432
0433 }  // namespace pixelCPEforGPU
0434
0435 #endif  // RecoLocalTracker_SiPixelRecHits_pixelCPEforGPU_h