Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:18:56

0001 #ifndef RecoLocalCalo_HcalRecProducers_src_KernelHelpers_h
0002 #define RecoLocalCalo_HcalRecProducers_src_KernelHelpers_h
0003 
0004 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalConstants.h"
0005 
0006 #include "DeclsForKernels.h"
0007 
0008 namespace hcal {
0009   namespace reconstruction {
0010 
0011     // this is from HcalTimeSlew.
0012     // HcalTimeSlew are values that come in from ESProducer that takes them
0013     // from a python config. see DeclsForKernels for more explanation
0014     __forceinline__ __device__ float compute_time_slew_delay(float const fC,
0015                                                              float const tzero,
0016                                                              float const slope,
0017                                                              float const tmax) {
0018       auto const rawDelay = tzero + slope * std::log(fC);
0019       return rawDelay < 0 ? 0 : (rawDelay > tmax ? tmax : rawDelay);
0020     }
0021 
0022     // HcalQIEShapes are hardcoded in HcalQIEData.cc basically
0023     // + some logic to generate 128 and 256 value arrays...
0024     __constant__ float const qie8shape[129] = {
0025         -1,   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,    10,   11,   12,   13,   14,   16,
0026         18,   20,   22,   24,   26,   28,   31,   34,   37,   40,   44,   48,   52,   57,   62,   57,   62,
0027         67,   72,   77,   82,   87,   92,   97,   102,  107,  112,  117,  122,  127,  132,  142,  152,  162,
0028         172,  182,  192,  202,  217,  232,  247,  262,  282,  302,  322,  347,  372,  347,  372,  397,  422,
0029         447,  472,  497,  522,  547,  572,  597,  622,  647,  672,  697,  722,  772,  822,  872,  922,  972,
0030         1022, 1072, 1147, 1222, 1297, 1372, 1472, 1572, 1672, 1797, 1922, 1797, 1922, 2047, 2172, 2297, 2422,
0031         2547, 2672, 2797, 2922, 3047, 3172, 3297, 3422, 3547, 3672, 3922, 4172, 4422, 4672, 4922, 5172, 5422,
0032         5797, 6172, 6547, 6922, 7422, 7922, 8422, 9047, 9672, 10297};
0033 
0034     __constant__ float const qie11shape[257] = {
0035         -0.5,    0.5,     1.5,     2.5,     3.5,     4.5,     5.5,     6.5,     7.5,     8.5,     9.5,     10.5,
0036         11.5,    12.5,    13.5,    14.5,    15.5,    17.5,    19.5,    21.5,    23.5,    25.5,    27.5,    29.5,
0037         31.5,    33.5,    35.5,    37.5,    39.5,    41.5,    43.5,    45.5,    47.5,    49.5,    51.5,    53.5,
0038         55.5,    59.5,    63.5,    67.5,    71.5,    75.5,    79.5,    83.5,    87.5,    91.5,    95.5,    99.5,
0039         103.5,   107.5,   111.5,   115.5,   119.5,   123.5,   127.5,   131.5,   135.5,   139.5,   147.5,   155.5,
0040         163.5,   171.5,   179.5,   187.5,   171.5,   179.5,   187.5,   195.5,   203.5,   211.5,   219.5,   227.5,
0041         235.5,   243.5,   251.5,   259.5,   267.5,   275.5,   283.5,   291.5,   299.5,   315.5,   331.5,   347.5,
0042         363.5,   379.5,   395.5,   411.5,   427.5,   443.5,   459.5,   475.5,   491.5,   507.5,   523.5,   539.5,
0043         555.5,   571.5,   587.5,   603.5,   619.5,   651.5,   683.5,   715.5,   747.5,   779.5,   811.5,   843.5,
0044         875.5,   907.5,   939.5,   971.5,   1003.5,  1035.5,  1067.5,  1099.5,  1131.5,  1163.5,  1195.5,  1227.5,
0045         1259.5,  1291.5,  1355.5,  1419.5,  1483.5,  1547.5,  1611.5,  1675.5,  1547.5,  1611.5,  1675.5,  1739.5,
0046         1803.5,  1867.5,  1931.5,  1995.5,  2059.5,  2123.5,  2187.5,  2251.5,  2315.5,  2379.5,  2443.5,  2507.5,
0047         2571.5,  2699.5,  2827.5,  2955.5,  3083.5,  3211.5,  3339.5,  3467.5,  3595.5,  3723.5,  3851.5,  3979.5,
0048         4107.5,  4235.5,  4363.5,  4491.5,  4619.5,  4747.5,  4875.5,  5003.5,  5131.5,  5387.5,  5643.5,  5899.5,
0049         6155.5,  6411.5,  6667.5,  6923.5,  7179.5,  7435.5,  7691.5,  7947.5,  8203.5,  8459.5,  8715.5,  8971.5,
0050         9227.5,  9483.5,  9739.5,  9995.5,  10251.5, 10507.5, 11019.5, 11531.5, 12043.5, 12555.5, 13067.5, 13579.5,
0051         12555.5, 13067.5, 13579.5, 14091.5, 14603.5, 15115.5, 15627.5, 16139.5, 16651.5, 17163.5, 17675.5, 18187.5,
0052         18699.5, 19211.5, 19723.5, 20235.5, 20747.5, 21771.5, 22795.5, 23819.5, 24843.5, 25867.5, 26891.5, 27915.5,
0053         28939.5, 29963.5, 30987.5, 32011.5, 33035.5, 34059.5, 35083.5, 36107.5, 37131.5, 38155.5, 39179.5, 40203.5,
0054         41227.5, 43275.5, 45323.5, 47371.5, 49419.5, 51467.5, 53515.5, 55563.5, 57611.5, 59659.5, 61707.5, 63755.5,
0055         65803.5, 67851.5, 69899.5, 71947.5, 73995.5, 76043.5, 78091.5, 80139.5, 82187.5, 84235.5, 88331.5, 92427.5,
0056         96523.5, 100620,  104716,  108812,  112908};
0057 
0058     // Conditions are transferred once per IOV
0059     // Access is performed based on the det id which is converted to a linear index
0060     // 2 funcs below are taken from HcalTopology (reimplemented here).
0061     // Inputs are constants that are also taken from HcalTopology
0062     // but passed to the kernel as arguments using the HclaTopology itself
0063     constexpr int32_t IPHI_MAX = 72;
0064 
0065     __forceinline__ __device__ uint32_t did2linearIndexHB(
0066         uint32_t const didraw, int const maxDepthHB, int const firstHBRing, int const lastHBRing, int const nEtaHB) {
0067       HcalDetId did{didraw};
0068       uint32_t const value = (did.depth() - 1) + maxDepthHB * (did.iphi() - 1);
0069       return did.ieta() > 0 ? value + maxDepthHB * hcal::reconstruction::IPHI_MAX * (did.ieta() - firstHBRing)
0070                             : value + maxDepthHB * hcal::reconstruction::IPHI_MAX * (did.ieta() + lastHBRing + nEtaHB);
0071     }
0072 
0073     __forceinline__ __device__ uint32_t did2linearIndexHE(uint32_t const didraw,
0074                                                           int const maxDepthHE,
0075                                                           int const maxPhiHE,
0076                                                           int const firstHERing,
0077                                                           int const lastHERing,
0078                                                           int const nEtaHE) {
0079       HcalDetId did{didraw};
0080       uint32_t const value = (did.depth() - 1) + maxDepthHE * (did.iphi() - 1);
0081       return did.ieta() > 0 ? value + maxDepthHE * maxPhiHE * (did.ieta() - firstHERing)
0082                             : value + maxDepthHE * maxPhiHE * (did.ieta() + lastHERing + nEtaHE);
0083     }
0084 
0085     __forceinline__ __device__ uint32_t get_qiecoder_index(uint32_t const capid, uint32_t const range) {
0086       return capid * 4 + range;
0087     }
0088 
0089     __forceinline__ __device__ float compute_reco_correction_factor(float const par1,
0090                                                                     float const par2,
0091                                                                     float const par3,
0092                                                                     float const x) {
0093       return par3 * x * x + par2 * x + par1;
0094     }
0095 
0096     // compute the charge using the adc, qie type and the appropriate qie shape array
0097     __forceinline__ __device__ float compute_coder_charge(
0098         int const qieType, uint8_t const adc, uint8_t const capid, float const* qieOffsets, float const* qieSlopes) {
0099       auto const range = qieType == 0 ? (adc >> 5) & 0x3 : (adc >> 6) & 0x3;
0100       auto const* qieShapeToUse = qieType == 0 ? qie8shape : qie11shape;
0101       auto const nbins = qieType == 0 ? 32 : 64;
0102       auto const center = adc % nbins == nbins - 1 ? 0.5 * (3 * qieShapeToUse[adc] - qieShapeToUse[adc - 1])
0103                                                    : 0.5 * (qieShapeToUse[adc] + qieShapeToUse[adc + 1]);
0104       auto const index = get_qiecoder_index(capid, range);
0105       return (center - qieOffsets[index]) / qieSlopes[index];
0106     }
0107 
0108     // this is from
0109     //  https://github.com/cms-sw/cmssw/blob/master/RecoLocalCalo/HcalRecProducers/src/HBHEPhase1Reconstructor.cc#L140
0110 
0111     __forceinline__ __device__ float compute_diff_charge_gain(int const qieType,
0112                                                               uint8_t adc,
0113                                                               uint8_t const capid,
0114                                                               float const* qieOffsets,
0115                                                               float const* qieSlopes,
0116                                                               bool const isqie11) {
0117       constexpr uint32_t mantissaMaskQIE8 = 0x1fu;
0118       constexpr uint32_t mantissaMaskQIE11 = 0x3f;
0119       auto const mantissaMask = isqie11 ? mantissaMaskQIE11 : mantissaMaskQIE8;
0120       auto const q = compute_coder_charge(qieType, adc, capid, qieOffsets, qieSlopes);
0121       auto const mantissa = adc & mantissaMask;
0122 
0123       if (mantissa == 0u || mantissa == mantissaMask - 1u)
0124         return compute_coder_charge(qieType, adc + 1u, capid, qieOffsets, qieSlopes) - q;
0125       else if (mantissa == 1u || mantissa == mantissaMask)
0126         return q - compute_coder_charge(qieType, adc - 1u, capid, qieOffsets, qieSlopes);
0127       else {
0128         auto const qup = compute_coder_charge(qieType, adc + 1u, capid, qieOffsets, qieSlopes);
0129         auto const qdown = compute_coder_charge(qieType, adc - 1u, capid, qieOffsets, qieSlopes);
0130         auto const upgain = qup - q;
0131         auto const downgain = q - qdown;
0132         auto const averagegain = (qup - qdown) / 2.f;
0133         if (std::abs(upgain - downgain) < 0.01f * averagegain)
0134           return averagegain;
0135         else {
0136           auto const q2up = compute_coder_charge(qieType, adc + 2u, capid, qieOffsets, qieSlopes);
0137           auto const q2down = compute_coder_charge(qieType, adc - 2u, capid, qieOffsets, qieSlopes);
0138           auto const upgain2 = q2up - qup;
0139           auto const downgain2 = qdown - q2down;
0140           if (std::abs(upgain2 - upgain) < std::abs(downgain2 - downgain))
0141             return upgain;
0142           else
0143             return downgain;
0144         }
0145       }
0146     }
0147 
0148     // TODO: remove what's not needed
0149     // originally from from RecoLocalCalo/HcalRecAlgos/src/PulseShapeFunctor.cc
0150     __forceinline__ __device__ float compute_pulse_shape_value(float const pulse_time,
0151                                                                int const sample,
0152                                                                int const shift,
0153                                                                float const* acc25nsVec,
0154                                                                float const* diff25nsItvlVec,
0155                                                                float const* accVarLenIdxMinusOneVec,
0156                                                                float const* diffVarItvlIdxMinusOneVec,
0157                                                                float const* accVarLenIdxZeroVec,
0158                                                                float const* diffVarItvlIdxZeroVec) {
0159       // constants
0160       constexpr float slew = 0.f;
0161       constexpr auto ns_per_bx = hcal::constants::nsPerBX;
0162 
0163       // FIXME: clean up all the rounding... this is coming from original cpu version
0164       float const i_start_float = -hcal::constants::iniTimeShift - pulse_time - slew > 0.f
0165                                       ? 0.f
0166                                       : std::abs(-hcal::constants::iniTimeShift - pulse_time - slew) + 1.f;
0167       int i_start = static_cast<int>(i_start_float);
0168       float offset_start = static_cast<float>(i_start) - hcal::constants::iniTimeShift - pulse_time - slew;
0169       // FIXME: do we need a check for nan???
0170 #ifdef HCAL_MAHI_GPUDEBUG
0171       if (shift == 0)
0172         printf("i_start_float = %f i_start = %d offset_start = %f\n", i_start_float, i_start, offset_start);
0173 #endif
0174 
0175       // boundary
0176       if (offset_start == 1.0f) {
0177         offset_start = 0.f;
0178         i_start -= 1;
0179       }
0180 
0181 #ifdef HCAL_MAHI_GPUDEBUG
0182       if (shift == 0)
0183         printf("i_start_float = %f i_start = %d offset_start = %f\n", i_start_float, i_start, offset_start);
0184 #endif
0185 
0186       int const bin_start = static_cast<int>(offset_start);
0187       auto const bin_start_up = static_cast<float>(bin_start) + 0.5f;
0188       int const bin_0_start = offset_start < bin_start_up ? bin_start - 1 : bin_start;
0189       int const its_start = i_start / ns_per_bx;
0190       int const distTo25ns_start = hcal::constants::nsPerBX - 1 - i_start % ns_per_bx;
0191       auto const factor = offset_start - static_cast<float>(bin_0_start) - 0.5;
0192 
0193 #ifdef HCAL_MAHI_GPUDEBUG
0194       if (shift == 0) {
0195         printf("bin_start = %d bin_0_start = %d its_start = %d distTo25ns_start = %d factor = %f\n",
0196                bin_start,
0197                bin_0_start,
0198                its_start,
0199                distTo25ns_start,
0200                factor);
0201       }
0202 #endif
0203 
0204       auto const sample_over10ts = sample + shift;
0205       float value = 0.0f;
0206       if (sample_over10ts == its_start) {
0207         value = bin_0_start == -1
0208                     ? accVarLenIdxMinusOneVec[distTo25ns_start] + factor * diffVarItvlIdxMinusOneVec[distTo25ns_start]
0209                     : accVarLenIdxZeroVec[distTo25ns_start] + factor * diffVarItvlIdxZeroVec[distTo25ns_start];
0210       } else if (sample_over10ts > its_start) {
0211         int const bin_idx = distTo25ns_start + 1 + (sample_over10ts - its_start - 1) * ns_per_bx + bin_0_start;
0212         value = acc25nsVec[bin_idx] + factor * diff25nsItvlVec[bin_idx];
0213       }
0214       return value;
0215     }
0216 
0217   }  // namespace reconstruction
0218 }  // namespace hcal
0219 
0220 #endif  // RecoLocalCalo_HcalRecProducers_src_KernelHelpers_h