Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #ifndef RecoLocalTracker_SiPixelClusterizer_plugins_alpaka_ClusterChargeCut_h
0002 #define RecoLocalTracker_SiPixelClusterizer_plugins_alpaka_ClusterChargeCut_h
0003 
0004 #include <cstdint>
0005 #include <cstdio>
0006 
0007 #include <alpaka/alpaka.hpp>
0008 
0009 #include "DataFormats/SiPixelClusterSoA/interface/ClusteringConstants.h"
0010 #include "DataFormats/SiPixelClusterSoA/interface/SiPixelClustersSoA.h"
0011 #include "DataFormats/SiPixelDigiSoA/interface/SiPixelDigisSoA.h"
0012 #include "HeterogeneousCore/AlpakaInterface/interface/prefixScan.h"
0013 #include "RecoLocalTracker/SiPixelClusterizer/interface/SiPixelClusterThresholds.h"
0014 
0015 //#define GPU_DEBUG
0016 
0017 namespace pixelClustering {
0018 
0019   template <typename TrackerTraits>
0020   struct ClusterChargeCut {
0021     template <typename TAcc>
0022     ALPAKA_FN_ACC void operator()(TAcc const& acc,
0023                                   SiPixelDigisSoAView digi_view,
0024                                   SiPixelClustersSoAView clus_view,
0025                                   // charge cut on cluster in electrons (for layer 1 and for other layers)
0026                                   SiPixelClusterThresholds clusterThresholds,
0027                                   const uint32_t numElements) const {
0028       constexpr int32_t maxNumClustersPerModules = TrackerTraits::maxNumClustersPerModules;
0029 
0030       auto& charge = alpaka::declareSharedVar<int32_t[maxNumClustersPerModules], __COUNTER__>(acc);
0031       auto& ok = alpaka::declareSharedVar<uint8_t[maxNumClustersPerModules], __COUNTER__>(acc);
0032       auto& newclusId = alpaka::declareSharedVar<uint16_t[maxNumClustersPerModules], __COUNTER__>(acc);
0033 
0034       constexpr int startBPIX2 = TrackerTraits::layerStart[1];
0035 
0036       ALPAKA_ASSERT_ACC(TrackerTraits::numberOfModules < maxNumModules);
0037       ALPAKA_ASSERT_ACC(startBPIX2 < TrackerTraits::numberOfModules);
0038 
0039       auto endModule = clus_view[0].moduleStart();
0040       for (auto module : cms::alpakatools::independent_groups(acc, endModule)) {
0041         auto firstPixel = clus_view[1 + module].moduleStart();
0042         auto thisModuleId = digi_view[firstPixel].moduleId();
0043         while (thisModuleId == invalidModuleId and firstPixel < numElements) {
0044           // skip invalid or duplicate pixels
0045           ++firstPixel;
0046           thisModuleId = digi_view[firstPixel].moduleId();
0047         }
0048         if (firstPixel >= numElements) {
0049           // reached the end of the input while skipping the invalid pixels, nothing left to do
0050           break;
0051         }
0052         if (thisModuleId != clus_view[module].moduleId()) {
0053           // reached the end of the module while skipping the invalid pixels, skip this module
0054           continue;
0055         }
0056         ALPAKA_ASSERT_ACC(thisModuleId < TrackerTraits::numberOfModules);
0057 
0058         uint32_t nclus = clus_view[thisModuleId].clusInModule();
0059         if (nclus == 0)
0060           return;
0061 
0062         if (cms::alpakatools::once_per_block(acc) && nclus > maxNumClustersPerModules)
0063           printf("Warning: too many clusters in module %u in block %u: %u > %d\n",
0064                  thisModuleId,
0065                  module,
0066                  nclus,
0067                  maxNumClustersPerModules);
0068 
0069         if (nclus > maxNumClustersPerModules) {
0070           // remove excess  FIXME find a way to cut charge first....
0071           for (auto i : cms::alpakatools::independent_group_elements(acc, firstPixel, numElements)) {
0072             if (digi_view[i].moduleId() == invalidModuleId)
0073               continue;  // not valid
0074             if (digi_view[i].moduleId() != thisModuleId)
0075               break;  // end of module
0076             if (digi_view[i].clus() >= maxNumClustersPerModules) {
0077               digi_view[i].moduleId() = invalidModuleId;
0078               digi_view[i].clus() = invalidModuleId;
0079             }
0080           }
0081           nclus = maxNumClustersPerModules;
0082         }
0083 
0084 #ifdef GPU_DEBUG
0085         if (thisModuleId % 100 == 1)
0086           if (cms::alpakatools::once_per_block(acc))
0087             printf("start cluster charge cut for module %d in block %d\n", thisModuleId, module);
0088 #endif
0089 
0090         ALPAKA_ASSERT_ACC(nclus <= maxNumClustersPerModules);
0091         for (auto i : cms::alpakatools::independent_group_elements(acc, nclus)) {
0092           charge[i] = 0;
0093         }
0094         alpaka::syncBlockThreads(acc);
0095 
0096         for (auto i : cms::alpakatools::independent_group_elements(acc, firstPixel, numElements)) {
0097           if (digi_view[i].moduleId() == invalidModuleId)
0098             continue;  // not valid
0099           if (digi_view[i].moduleId() != thisModuleId)
0100             break;  // end of module
0101           alpaka::atomicAdd(acc,
0102                             &charge[digi_view[i].clus()],
0103                             static_cast<int32_t>(digi_view[i].adc()),
0104                             alpaka::hierarchy::Threads{});
0105         }
0106         alpaka::syncBlockThreads(acc);
0107 
0108         auto chargeCut = clusterThresholds.getThresholdForLayerOnCondition(thisModuleId < startBPIX2);
0109 
0110         bool good = true;
0111         for (auto i : cms::alpakatools::independent_group_elements(acc, nclus)) {
0112           newclusId[i] = ok[i] = (charge[i] >= chargeCut) ? 1 : 0;
0113           if (0 == ok[i])
0114             good = false;
0115         }
0116 
0117         // if all clusters are above threshold, do nothing
0118         if (alpaka::syncBlockThreadsPredicate<alpaka::BlockAnd>(acc, good))
0119           continue;
0120 
0121         // renumber
0122         auto& ws = alpaka::declareSharedVar<uint16_t[32], __COUNTER__>(acc);
0123         // FIXME this value should come from cms::alpakatools::blockPrefixScan itself
0124         constexpr uint32_t maxThreads = 1024;
0125         auto minClust = std::min(nclus, maxThreads);
0126 
0127         cms::alpakatools::blockPrefixScan(acc, newclusId, minClust, ws);
0128         if constexpr (maxNumClustersPerModules > maxThreads)  //only if needed
0129         {
0130           for (uint32_t offset = maxThreads; offset < nclus; offset += maxThreads) {
0131             cms::alpakatools::blockPrefixScan(acc, newclusId + offset, nclus - offset, ws);
0132             for (uint32_t i : cms::alpakatools::independent_group_elements(acc, offset, nclus)) {
0133               uint32_t prevBlockEnd = (i / maxThreads) * maxThreads - 1;
0134               newclusId[i] += newclusId[prevBlockEnd];
0135             }
0136             alpaka::syncBlockThreads(acc);
0137           }
0138         }
0139         ALPAKA_ASSERT_ACC(nclus >= newclusId[nclus - 1]);
0140 
0141         clus_view[thisModuleId].clusInModule() = newclusId[nclus - 1];
0142 
0143         // reassign id
0144         for (auto i : cms::alpakatools::independent_group_elements(acc, firstPixel, numElements)) {
0145           if (digi_view[i].moduleId() == invalidModuleId)
0146             continue;  // not valid
0147           if (digi_view[i].moduleId() != thisModuleId)
0148             break;  // end of module
0149           if (0 == ok[digi_view[i].clus()])
0150             digi_view[i].moduleId() = digi_view[i].clus() = invalidModuleId;
0151           else
0152             digi_view[i].clus() = newclusId[digi_view[i].clus()] - 1;
0153         }
0154 
0155         // done
0156         alpaka::syncBlockThreads(acc);
0157       }
0158     }
0159   };
0160 
0161 }  // namespace pixelClustering
0162 
0163 #endif  // RecoLocalTracker_SiPixelClusterizer_plugins_alpaka_ClusterChargeCut_h