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