gpuPixelRecHits Namespace Reference

Functions
	__attribute__ ((always_inline)) void getHits(pixelCPEforGPU

Function Documentation

__attribute__()

gpuPixelRecHits::__attribute__ ( (always_inline)  )

inline

Definition at line 18 of file gpuPixelRecHits.h.

References cms::cudacompat::__syncthreads(), SiPixelDigisCUDASOAView::adc(), cms::cuda::assert(), cms::cudacompat::atomicAdd(), cms::cudacompat::atomicMax(), cms::cudacompat::atomicMin(), cms::cudacompat::blockDim, cms::cudacompat::blockIdx, cms::cuda::bs, haddnano::cl, SiPixelDigisCUDASOAView::clus(), bsc_activity_cfg::clusters, pixelCPEforGPU::errorFromDB(), pixelCPEforGPU::errorFromSize(), first, h, hfClusterShapes_cfi::hits, mps_fire::i, gpuClustering::invalidModuleId, HLT_2022v15_cff::isPhase2, dqmdumpme::k, SiStripPI::max, pixelCPEforGPU::MaxHitsInIter, hlt_dqm_clientPB-live_cfg::me, SiStripPI::min, SiPixelDigisCUDASOAView::moduleInd(), phase1PixelTopology::numberOfLaddersInBarrel, phase2PixelTopology::numberOfLaddersInBarrel, gpuClustering::numElements, phits, pixelCPEforGPU::position(), mathSSE::sqrt(), cms::cudacompat::threadIdx, x, SiPixelDigisCUDASOAView::xx(), and SiPixelDigisCUDASOAView::yy().

                                                          {
    // FIXME
    // the compiler seems NOT to optimize loads from views (even in a simple test case)
    // The whole gimnastic here of copying or not is a pure heuristic exercise that seems to produce the fastest code with the above signature
    // not using views (passing a gazzilion of array pointers) seems to produce the fastest code (but it is harder to mantain)
    assert(phits);
    assert(cpeParams);
    auto& hits = *phits;

    auto const& clusters = *pclusters;
    auto isPhase2 = cpeParams->commonParams().isPhase2;
    // copy average geometry corrected by beamspot . FIXME (move it somewhere else???)
    if (0 == blockIdx.x) {
      auto& agc = hits.averageGeometry();
      auto const& ag = cpeParams->averageGeometry();
      auto nLadders =
          isPhase2 ? phase2PixelTopology::numberOfLaddersInBarrel : phase1PixelTopology::numberOfLaddersInBarrel;

      for (int il = threadIdx.x, nl = nLadders; il < nl; il += blockDim.x) {
        agc.ladderZ[il] = ag.ladderZ[il] - bs->z;
        agc.ladderX[il] = ag.ladderX[il] - bs->x;
        agc.ladderY[il] = ag.ladderY[il] - bs->y;
        agc.ladderR[il] = sqrt(agc.ladderX[il] * agc.ladderX[il] + agc.ladderY[il] * agc.ladderY[il]);
        agc.ladderMinZ[il] = ag.ladderMinZ[il] - bs->z;
        agc.ladderMaxZ[il] = ag.ladderMaxZ[il] - bs->z;
      }

      if (0 == threadIdx.x) {
        agc.endCapZ[0] = ag.endCapZ[0] - bs->z;
        agc.endCapZ[1] = ag.endCapZ[1] - bs->z;
        //         printf("endcapZ %f %f\n",agc.endCapZ[0],agc.endCapZ[1]);
      }
    }

    // to be moved in common namespace...
    using gpuClustering::invalidModuleId;
    constexpr int32_t MaxHitsInIter = pixelCPEforGPU::MaxHitsInIter;

    using ClusParams = pixelCPEforGPU::ClusParams;

    // as usual one block per module
    __shared__ ClusParams clusParams;

    auto me = clusters.moduleId(blockIdx.x);
    int nclus = clusters.clusInModule(me);

    if (0 == nclus)
      return;
#ifdef GPU_DEBUG
    if (threadIdx.x == 0) {
      auto k = clusters.moduleStart(1 + blockIdx.x);
      while (digis.moduleInd(k) == invalidModuleId)
        ++k;
      assert(digis.moduleInd(k) == me);
    }
#endif
#ifdef GPU_DEBUG
    if (me % 100 == 1)
      if (threadIdx.x == 0)
        printf("hitbuilder: %d clusters in module %d. will write at %d\n", nclus, me, clusters.clusModuleStart(me));
#endif
    for (int startClus = 0, endClus = nclus; startClus < endClus; startClus += MaxHitsInIter) {
      int nClusInIter = std::min(MaxHitsInIter, endClus - startClus);
      int lastClus = startClus + nClusInIter;
      assert(nClusInIter <= nclus);
      assert(nClusInIter > 0);
      assert(lastClus <= nclus);

      assert(nclus > MaxHitsInIter || (0 == startClus && nClusInIter == nclus && lastClus == nclus));

      // init
      for (int ic = threadIdx.x; ic < nClusInIter; ic += blockDim.x) {
        clusParams.minRow[ic] = std::numeric_limits<uint32_t>::max();
        clusParams.maxRow[ic] = 0;
        clusParams.minCol[ic] = std::numeric_limits<uint32_t>::max();
        clusParams.maxCol[ic] = 0;
        clusParams.charge[ic] = 0;
        clusParams.q_f_X[ic] = 0;
        clusParams.q_l_X[ic] = 0;
        clusParams.q_f_Y[ic] = 0;
        clusParams.q_l_Y[ic] = 0;
      }

      __syncthreads();

      // one thread per "digi"
      auto first = clusters.moduleStart(1 + blockIdx.x) + threadIdx.x;
      for (int i = first; i < numElements; i += blockDim.x) {
        auto id = digis.moduleInd(i);
        if (id == invalidModuleId)
          continue;  // not valid
        if (id != me)
          break;  // end of module
        auto cl = digis.clus(i);
        if (cl < startClus || cl >= lastClus)
          continue;
        cl -= startClus;
        assert(cl >= 0);
        assert(cl < MaxHitsInIter);
        auto x = digis.xx(i);
        auto y = digis.yy(i);
        atomicMin(&clusParams.minRow[cl], x);
        atomicMax(&clusParams.maxRow[cl], x);
        atomicMin(&clusParams.minCol[cl], y);
        atomicMax(&clusParams.maxCol[cl], y);
      }

      __syncthreads();

      auto pixmx = cpeParams->detParams(me).pixmx;
      for (int i = first; i < numElements; i += blockDim.x) {
        auto id = digis.moduleInd(i);
        if (id == invalidModuleId)
          continue;  // not valid
        if (id != me)
          break;  // end of module
        auto cl = digis.clus(i);
        if (cl < startClus || cl >= lastClus)
          continue;
        cl -= startClus;
        assert(cl >= 0);
        assert(cl < MaxHitsInIter);
        auto x = digis.xx(i);
        auto y = digis.yy(i);
        auto ch = digis.adc(i);
        atomicAdd(&clusParams.charge[cl], ch);
        ch = std::min(ch, pixmx);
        if (clusParams.minRow[cl] == x)
          atomicAdd(&clusParams.q_f_X[cl], ch);
        if (clusParams.maxRow[cl] == x)
          atomicAdd(&clusParams.q_l_X[cl], ch);
        if (clusParams.minCol[cl] == y)
          atomicAdd(&clusParams.q_f_Y[cl], ch);
        if (clusParams.maxCol[cl] == y)
          atomicAdd(&clusParams.q_l_Y[cl], ch);
      }

      __syncthreads();

      // next one cluster per thread...

      first = clusters.clusModuleStart(me) + startClus;
      for (int ic = threadIdx.x; ic < nClusInIter; ic += blockDim.x) {
        auto h = first + ic;  // output index in global memory

        assert(h < hits.nHits());
        assert(h < clusters.clusModuleStart(me + 1));

        pixelCPEforGPU::position(cpeParams->commonParams(), cpeParams->detParams(me), clusParams, ic);
        if (!isPhase2)
          pixelCPEforGPU::errorFromDB(cpeParams->commonParams(), cpeParams->detParams(me), clusParams, ic);
        else
          pixelCPEforGPU::errorFromSize(cpeParams->commonParams(), cpeParams->detParams(me), clusParams, ic);

        // store it
        hits.setChargeAndStatus(h, clusParams.charge[ic], clusParams.status[ic]);
        hits.detectorIndex(h) = me;

        float xl, yl;
        hits.xLocal(h) = xl = clusParams.xpos[ic];
        hits.yLocal(h) = yl = clusParams.ypos[ic];

        hits.clusterSizeX(h) = clusParams.xsize[ic];
        hits.clusterSizeY(h) = clusParams.ysize[ic];

        hits.xerrLocal(h) = clusParams.xerr[ic] * clusParams.xerr[ic] + cpeParams->detParams(me).apeXX;
        hits.yerrLocal(h) = clusParams.yerr[ic] * clusParams.yerr[ic] + cpeParams->detParams(me).apeYY;

        // keep it local for computations
        float xg, yg, zg;
        // to global and compute phi...
        cpeParams->detParams(me).frame.toGlobal(xl, yl, xg, yg, zg);
        // here correct for the beamspot...
        xg -= bs->x;
        yg -= bs->y;
        zg -= bs->z;

        hits.xGlobal(h) = xg;
        hits.yGlobal(h) = yg;
        hits.zGlobal(h) = zg;

        hits.rGlobal(h) = std::sqrt(xg * xg + yg * yg);
        hits.iphi(h) = unsafe_atan2s<7>(yg, xg);
      }
      __syncthreads();
    }  // end loop on batches
  }