SiPixelRecHits/interface/pixelCPEforDevice.h

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