File indexing completed on 2025-02-05 23:51:43
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
0020
0021
0022 constexpr int kGenErrorQBins = 5;
0023
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
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];
0055 int16_t ysize[N];
0056
0057 Status status[N];
0058 };
0059
0060
0061 struct CommonParams {
0062 float theThicknessB;
0063 float theThicknessE;
0064
0065 uint16_t maxModuleStride;
0066 uint8_t numberOfLaddersInBarrel;
0067 };
0068
0069 struct DetParams {
0070 bool isBarrel;
0071 bool isPosZ;
0072 uint16_t layer;
0073 uint16_t index;
0074 uint32_t rawId;
0075
0076 float shiftX;
0077 float shiftY;
0078 float thePitchX;
0079 float thePitchY;
0080 float chargeWidthX;
0081 float chargeWidthY;
0082 uint16_t pixmx;
0083
0084 uint16_t nRowsRoc;
0085 uint16_t nColsRoc;
0086 uint16_t nRows;
0087 uint16_t nCols;
0088
0089 uint32_t numPixsInModule;
0090
0091 float x0, y0, z0;
0092
0093 float apeXX, apeYY;
0094 uint8_t sx2, sy1, sy2;
0095 uint8_t sigmax[kNumErrorBins], sigmax1[kNumErrorBins],
0096 sigmay[kNumErrorBins];
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
0116 auto gvx = x - detParams.x0;
0117 auto gvy = y - detParams.y0;
0118 auto gvz = -1.f / detParams.z0;
0119
0120
0121 cotalpha = gvx * gvz;
0122 cotbeta = gvy * gvz;
0123 }
0124
0125 constexpr inline float correction(int sizeM1,
0126 int q_f,
0127 int q_l,
0128 uint16_t upper_edge_first_pix,
0129 uint16_t lower_edge_last_pix,
0130 float lorentz_shift,
0131 float theThickness,
0132 float cot_angle,
0133 float pitch,
0134 bool first_is_big,
0135 bool last_is_big)
0136 {
0137 if (0 == sizeM1)
0138 return 0;
0139
0140 float w_eff = 0;
0141 bool simple = true;
0142 if (1 == sizeM1) {
0143
0144
0145
0146 auto w_inner = pitch * float(lower_edge_last_pix - upper_edge_first_pix);
0147
0148
0149 auto w_pred = theThickness * cot_angle
0150 - lorentz_shift;
0151
0152 w_eff = std::abs(w_pred) - w_inner;
0153
0154
0155
0156
0157
0158
0159
0160 simple = (w_eff < 0.0f) | (w_eff > pitch);
0161 }
0162
0163 if (simple) {
0164
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;
0171 }
0172
0173
0174 float qdiff = q_l - q_f;
0175 float qsum = q_l + q_f;
0176
0177
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
0191 uint16_t llx = cp.minRow[ic] + 1;
0192 uint16_t lly = cp.minCol[ic] + 1;
0193
0194
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);
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
0238
0239
0240
0241 float xPos = std::fmaf(0.5f * ((float)mx - (float)detParams.nRows) - TrackerTraits::bigPixXCorrection,
0242 detParams.thePitchX,
0243 detParams.shiftX);
0244 float yPos = std::fmaf(0.5f * ((float)my - (float)detParams.nCols) - TrackerTraits::bigPixYCorrection,
0245 detParams.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,
0260 thickness,
0261 cotalpha,
0262 detParams.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,
0272 thickness,
0273 cotbeta,
0274 detParams.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
0288 cp.xerr[ic] = 0.0050;
0289 cp.yerr[ic] = 0.0085;
0290
0291
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};
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
0312 bool isEdgeX = cp.xsize[ic] < 1;
0313 bool isEdgeY = cp.ysize[ic] < 1;
0314
0315
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
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
0353 bool isEdgeX = cp.xsize[ic] < 1;
0354 bool isEdgeY = cp.ysize[ic] < 1;
0355
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
0365 if (ch < detParams.minCh[bin + 1])
0366 break;
0367
0368
0369
0370
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) * detParams.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
0382
0383
0384
0385 int tmp_max = std::max<int>(bin_value, low_value);
0386 int jx = std::min<int>(tmp_max, high_value);
0387
0388 auto toCM = [](uint8_t x) { return float(x) * 1.e-4f; };
0389
0390 if (not isEdgeX) {
0391 cp.xerr[ic] = isOneX ? toCM(isBigX ? detParams.sx2 : detParams.sigmax1[jx])
0392 : detParams.xfact[bin] * toCM(detParams.sigmax[jx]);
0393 }
0394
0395 auto ey = cp.ysize[ic] > 8 ? detParams.sigmay[std::min(cp.ysize[ic] - 9, 15)] : detParams.sy1;
0396 if (not isEdgeY) {
0397 cp.yerr[ic] = isOneY ? toCM(isBigY ? detParams.sy2 : detParams.sy1) : detParams.yfact[bin] * toCM(ey);
0398 }
0399 }
0400
0401
0402 template <>
0403 constexpr inline void errorFromDB<pixelTopology::Phase2>(CommonParams const& __restrict__ comParams,
0404 DetParams const& __restrict__ detParams,
0405 ClusParams& cp,
0406 uint32_t ic) {
0407 errorFromSize<pixelTopology::Phase2>(comParams, detParams, cp, ic);
0408 }
0409
0410 template <typename TrackerTopology>
0411 struct ParamsOnDeviceT {
0412 using LayerGeometry = LayerGeometryT<TrackerTopology>;
0413 using AverageGeometry = pixelTopology::AverageGeometryT<TrackerTopology>;
0414
0415 CommonParams m_commonParams;
0416
0417 DetParams m_detParams[TrackerTopology::numberOfModules];
0418 LayerGeometry m_layerGeometry;
0419 AverageGeometry m_averageGeometry;
0420
0421 constexpr CommonParams const& __restrict__ commonParams() const { return m_commonParams; }
0422 constexpr DetParams const& __restrict__ detParams(int i) const { return m_detParams[i]; }
0423 constexpr LayerGeometry const& __restrict__ layerGeometry() const { return m_layerGeometry; }
0424 constexpr AverageGeometry const& __restrict__ averageGeometry() const { return m_averageGeometry; }
0425
0426 CommonParams& commonParams() { return m_commonParams; }
0427 DetParams& detParams(int i) { return m_detParams[i]; }
0428 LayerGeometry& layerGeometry() { return m_layerGeometry; }
0429 AverageGeometry& averageGeometry() { return m_averageGeometry; }
0430
0431 constexpr uint8_t layer(uint16_t id) const { return m_layerGeometry.layer[id / TrackerTopology::maxModuleStride]; };
0432 };
0433
0434 }
0435
0436 #endif