Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:26:22

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     // apply the lorentz offset correction
0238     float xoff = 0.5f * float(detParams.nRows) * comParams.thePitchX;
0239     float yoff = 0.5f * float(detParams.nCols) * comParams.thePitchY;
0240 
0241     //correction for bigpixels for phase1
0242     xoff = xoff + TrackerTraits::bigPixXCorrection * comParams.thePitchX;
0243     yoff = yoff + TrackerTraits::bigPixYCorrection * comParams.thePitchY;
0244 
0245     // apply the lorentz offset correction
0246     auto xPos = detParams.shiftX + (comParams.thePitchX * 0.5f * float(mx)) - xoff;
0247     auto yPos = detParams.shiftY + (comParams.thePitchY * 0.5f * float(my)) - yoff;
0248 
0249     float cotalpha = 0, cotbeta = 0;
0250 
0251     computeAnglesFromDet(detParams, xPos, yPos, cotalpha, cotbeta);
0252 
0253     auto thickness = detParams.isBarrel ? comParams.theThicknessB : comParams.theThicknessE;
0254 
0255     auto xcorr = correction(cp.maxRow[ic] - cp.minRow[ic],
0256                             cp.q_f_X[ic],
0257                             cp.q_l_X[ic],
0258                             llxl,
0259                             urxl,
0260                             detParams.chargeWidthX,  // lorentz shift in cm
0261                             thickness,
0262                             cotalpha,
0263                             comParams.thePitchX,
0264                             TrackerTraits::isBigPixX(cp.minRow[ic]),
0265                             TrackerTraits::isBigPixX(cp.maxRow[ic]));
0266 
0267     auto ycorr = correction(cp.maxCol[ic] - cp.minCol[ic],
0268                             cp.q_f_Y[ic],
0269                             cp.q_l_Y[ic],
0270                             llyl,
0271                             uryl,
0272                             detParams.chargeWidthY,  // lorentz shift in cm
0273                             thickness,
0274                             cotbeta,
0275                             comParams.thePitchY,
0276                             TrackerTraits::isBigPixY(cp.minCol[ic]),
0277                             TrackerTraits::isBigPixY(cp.maxCol[ic]));
0278 
0279     cp.xpos[ic] = xPos + xcorr;
0280     cp.ypos[ic] = yPos + ycorr;
0281   }
0282 
0283   template <typename TrackerTraits>
0284   constexpr inline void errorFromSize(CommonParams const& __restrict__ comParams,
0285                                       DetParams const& __restrict__ detParams,
0286                                       ClusParams& cp,
0287                                       uint32_t ic) {
0288     // Edge cluster errors
0289     cp.xerr[ic] = 0.0050;
0290     cp.yerr[ic] = 0.0085;
0291 
0292     // FIXME these are errors form Run1
0293     float xerr_barrel_l1_def = TrackerTraits::xerr_barrel_l1_def;
0294     float yerr_barrel_l1_def = TrackerTraits::yerr_barrel_l1_def;
0295     float xerr_barrel_ln_def = TrackerTraits::xerr_barrel_ln_def;
0296     float yerr_barrel_ln_def = TrackerTraits::yerr_barrel_ln_def;
0297     float xerr_endcap_def = TrackerTraits::xerr_endcap_def;
0298     float yerr_endcap_def = TrackerTraits::yerr_endcap_def;
0299 
0300     constexpr float xerr_barrel_l1[] = {0.00115, 0.00120, 0.00088};  //TODO MOVE THESE SOMEWHERE ELSE
0301     constexpr float yerr_barrel_l1[] = {
0302         0.00375, 0.00230, 0.00250, 0.00250, 0.00230, 0.00230, 0.00210, 0.00210, 0.00240};
0303     constexpr float xerr_barrel_ln[] = {0.00115, 0.00120, 0.00088};
0304     constexpr float yerr_barrel_ln[] = {
0305         0.00375, 0.00230, 0.00250, 0.00250, 0.00230, 0.00230, 0.00210, 0.00210, 0.00240};
0306     constexpr float xerr_endcap[] = {0.0020, 0.0020};
0307     constexpr float yerr_endcap[] = {0.00210};
0308 
0309     auto sx = cp.maxRow[ic] - cp.minRow[ic];
0310     auto sy = cp.maxCol[ic] - cp.minCol[ic];
0311 
0312     // is edgy ?
0313     bool isEdgeX = cp.xsize[ic] < 1;
0314     bool isEdgeY = cp.ysize[ic] < 1;
0315 
0316     // is one and big?
0317     bool isBig1X = ((0 == sx) && TrackerTraits::isBigPixX(cp.minRow[ic]));
0318     bool isBig1Y = ((0 == sy) && TrackerTraits::isBigPixY(cp.minCol[ic]));
0319 
0320     if (!isEdgeX && !isBig1X) {
0321       if (not detParams.isBarrel) {
0322         cp.xerr[ic] = sx < std::size(xerr_endcap) ? xerr_endcap[sx] : xerr_endcap_def;
0323       } else if (detParams.layer == 1) {
0324         cp.xerr[ic] = sx < std::size(xerr_barrel_l1) ? xerr_barrel_l1[sx] : xerr_barrel_l1_def;
0325       } else {
0326         cp.xerr[ic] = sx < std::size(xerr_barrel_ln) ? xerr_barrel_ln[sx] : xerr_barrel_ln_def;
0327       }
0328     }
0329 
0330     if (!isEdgeY && !isBig1Y) {
0331       if (not detParams.isBarrel) {
0332         cp.yerr[ic] = sy < std::size(yerr_endcap) ? yerr_endcap[sy] : yerr_endcap_def;
0333       } else if (detParams.layer == 1) {
0334         cp.yerr[ic] = sy < std::size(yerr_barrel_l1) ? yerr_barrel_l1[sy] : yerr_barrel_l1_def;
0335       } else {
0336         cp.yerr[ic] = sy < std::size(yerr_barrel_ln) ? yerr_barrel_ln[sy] : yerr_barrel_ln_def;
0337       }
0338     }
0339   }
0340 
0341   template <typename TrackerTraits>
0342   constexpr inline void errorFromDB(CommonParams const& __restrict__ comParams,
0343                                     DetParams const& __restrict__ detParams,
0344                                     ClusParams& cp,
0345                                     uint32_t ic) {
0346     // Edge cluster errors
0347     cp.xerr[ic] = 0.0050f;
0348     cp.yerr[ic] = 0.0085f;
0349 
0350     auto sx = cp.maxRow[ic] - cp.minRow[ic];
0351     auto sy = cp.maxCol[ic] - cp.minCol[ic];
0352 
0353     // is edgy ?  (size is set negative: see above)
0354     bool isEdgeX = cp.xsize[ic] < 1;
0355     bool isEdgeY = cp.ysize[ic] < 1;
0356     // is one and big?
0357     bool isOneX = (0 == sx);
0358     bool isOneY = (0 == sy);
0359     bool isBigX = TrackerTraits::isBigPixX(cp.minRow[ic]);
0360     bool isBigY = TrackerTraits::isBigPixY(cp.minCol[ic]);
0361 
0362     auto ch = cp.charge[ic];
0363     auto bin = 0;
0364     for (; bin < kGenErrorQBins - 1; ++bin)
0365       // find first bin which minimum charge exceeds cluster charge
0366       if (ch < detParams.minCh[bin + 1])
0367         break;
0368 
0369     // in detParams qBins are reversed bin0 -> smallest charge, bin4-> largest charge
0370     // whereas in CondFormats/SiPixelTransient/src/SiPixelGenError.cc it is the opposite
0371     // so we reverse the bin here -> kGenErrorQBins - 1 - bin
0372     cp.status[ic].qBin = kGenErrorQBins - 1 - bin;
0373     cp.status[ic].isOneX = isOneX;
0374     cp.status[ic].isBigX = (isOneX & isBigX) | isEdgeX;
0375     cp.status[ic].isOneY = isOneY;
0376     cp.status[ic].isBigY = (isOneY & isBigY) | isEdgeY;
0377 
0378     auto xoff = -float(TrackerTraits::xOffset) * comParams.thePitchX;
0379     int low_value = 0;
0380     int high_value = kNumErrorBins - 1;
0381     int bin_value = float(kNumErrorBins) * (cp.xpos[ic] + xoff) / (2 * xoff);
0382     // return estimated bin value truncated to [0, 15]
0383     int jx = std::clamp(bin_value, low_value, high_value);
0384 
0385     auto toCM = [](uint8_t x) { return float(x) * 1.e-4f; };
0386 
0387     if (not isEdgeX) {
0388       cp.xerr[ic] = isOneX ? toCM(isBigX ? detParams.sx2 : detParams.sigmax1[jx])
0389                            : detParams.xfact[bin] * toCM(detParams.sigmax[jx]);
0390     }
0391 
0392     auto ey = cp.ysize[ic] > 8 ? detParams.sigmay[std::min(cp.ysize[ic] - 9, 15)] : detParams.sy1;
0393     if (not isEdgeY) {
0394       cp.yerr[ic] = isOneY ? toCM(isBigY ? detParams.sy2 : detParams.sy1) : detParams.yfact[bin] * toCM(ey);
0395     }
0396   }
0397 
0398   //for Phase2 -> fallback to error from size
0399   template <>
0400   constexpr inline void errorFromDB<pixelTopology::Phase2>(CommonParams const& __restrict__ comParams,
0401                                                            DetParams const& __restrict__ detParams,
0402                                                            ClusParams& cp,
0403                                                            uint32_t ic) {
0404     errorFromSize<pixelTopology::Phase2>(comParams, detParams, cp, ic);
0405   }
0406 
0407   template <typename TrackerTopology>
0408   struct ParamsOnDeviceT {
0409     using LayerGeometry = LayerGeometryT<TrackerTopology>;
0410     using AverageGeometry = pixelTopology::AverageGeometryT<TrackerTopology>;
0411 
0412     CommonParams m_commonParams;
0413     // Will contain an array of DetParams instances
0414     DetParams m_detParams[TrackerTopology::numberOfModules];
0415     LayerGeometry m_layerGeometry;
0416     AverageGeometry m_averageGeometry;
0417 
0418     constexpr CommonParams const& __restrict__ commonParams() const { return m_commonParams; }
0419     constexpr DetParams const& __restrict__ detParams(int i) const { return m_detParams[i]; }
0420     constexpr LayerGeometry const& __restrict__ layerGeometry() const { return m_layerGeometry; }
0421     constexpr AverageGeometry const& __restrict__ averageGeometry() const { return m_averageGeometry; }
0422 
0423     CommonParams& commonParams() { return m_commonParams; }
0424     DetParams& detParams(int i) { return m_detParams[i]; }
0425     LayerGeometry& layerGeometry() { return m_layerGeometry; }
0426     AverageGeometry& averageGeometry() { return m_averageGeometry; }
0427 
0428     constexpr uint8_t layer(uint16_t id) const { return m_layerGeometry.layer[id / TrackerTopology::maxModuleStride]; };
0429   };
0430 
0431 }  // namespace pixelCPEforDevice
0432 
0433 #endif  // RecoLocalTracker_SiPixelRecHits_interface_pixelCPEforDevice_h