From 164e86b8a48d722b4986890a35f24f7d9789e321 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Fri, 24 May 2024 09:25:59 +0200 Subject: [PATCH] First heterogeneous CLUE workflow. --- .../python/heterogeneousCLUE_cff.py | 5 + DataFormats/HGCalReco/BuildFile.xml | 4 + .../HGCalReco/interface/HGCalSoACells.h | 27 ++ .../interface/HGCalSoACellsHostCollection.h | 10 + .../HGCalReco/interface/HGCalSoACellsOut.h | 20 ++ .../HGCalSoACellsOutHostCollection.h | 10 + .../HGCalReco/interface/HGCalSoAClusters.h | 20 ++ .../HGCalSoAClustersHostCollection.h | 9 + .../interface/HGCalSoAClustersService.h | 19 ++ .../HGCalSoAClustersServiceHostCollection.h | 9 + .../alpaka/HGCalSoACellsDeviceCollection.h | 15 + .../alpaka/HGCalSoACellsOutDeviceCollection.h | 15 + .../alpaka/HGCalSoAClustersDeviceCollection.h | 14 + .../HGCalSoAClustersServiceDeviceCollection.h | 14 + .../HGCalReco/src/alpaka/classes_cuda.h | 7 + .../HGCalReco/src/alpaka/classes_cuda_def.xml | 13 + .../HGCalReco/src/alpaka/classes_rocm.h | 7 + .../HGCalReco/src/alpaka/classes_rocm_def.xml | 13 + DataFormats/HGCalReco/src/classes.h | 6 + DataFormats/HGCalReco/src/classes_def.xml | 54 ++++ .../scripts/portableHostCollectionHints | 22 +- .../hgcalLayerClustersFromSoAProducer_cfi.py | 14 + .../hgcalSoALayerClustersProducer_cfi.py | 13 + ...gcalSoARecHitsLayerClustersProducer_cfi.py | 10 + .../modules/hgcalSoARecHitsProducer_cfi.py | 19 ++ .../sequences/hgcalLocalRecoSequence_cfi.py | 30 +- .../python/HLT_75e33_timing_cff.py | 10 + .../interface/DumpClustersDetails.h | 186 +++++++++++ .../interface/HGCalLayerTiles.h | 18 +- .../interface/HGCalTilesConstants.h | 6 + .../HGCalRecProducers/plugins/BuildFile.xml | 71 +++-- .../plugins/HGCalCLUEAlgo.cc | 32 +- .../HGCalRecProducers/plugins/HGCalCLUEAlgo.h | 7 + .../HGCalLayerClusterHeterogeneousDumper.cc | 41 +++ ...HGCalLayerClusterHeterogeneousSoADumper.cc | 51 ++++ .../plugins/HGCalLayerClusterProducer.cc | 26 +- .../HGCalLayerClustersFromSoAProducer.cc | 165 ++++++++++ .../plugins/MergeClusterProducer.cc | 18 +- .../HGCalLayerClustersAlgoWrapper.dev.cc | 38 +++ .../alpaka/HGCalLayerClustersAlgoWrapper.h | 24 ++ .../HGCalLayerClustersSoAAlgoWrapper.dev.cc | 288 ++++++++++++++++++ .../alpaka/HGCalLayerClustersSoAAlgoWrapper.h | 30 ++ .../alpaka/HGCalSoALayerClustersProducer.cc | 104 +++++++ .../HGCalSoARecHitsLayerClustersProducer.cc | 69 +++++ .../plugins/alpaka/HGCalSoARecHitsProducer.cc | 221 ++++++++++++++ 45 files changed, 1750 insertions(+), 54 deletions(-) create mode 100644 Configuration/ProcessModifiers/python/heterogeneousCLUE_cff.py create mode 100644 DataFormats/HGCalReco/interface/HGCalSoACells.h create mode 100644 DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h create mode 100644 DataFormats/HGCalReco/interface/HGCalSoACellsOut.h create mode 100644 DataFormats/HGCalReco/interface/HGCalSoACellsOutHostCollection.h create mode 100644 DataFormats/HGCalReco/interface/HGCalSoAClusters.h create mode 100644 DataFormats/HGCalReco/interface/HGCalSoAClustersHostCollection.h create mode 100644 DataFormats/HGCalReco/interface/HGCalSoAClustersService.h create mode 100644 DataFormats/HGCalReco/interface/HGCalSoAClustersServiceHostCollection.h create mode 100644 DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsDeviceCollection.h create mode 100644 DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsOutDeviceCollection.h create mode 100644 DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersDeviceCollection.h create mode 100644 DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersServiceDeviceCollection.h create mode 100644 DataFormats/HGCalReco/src/alpaka/classes_cuda.h create mode 100644 DataFormats/HGCalReco/src/alpaka/classes_cuda_def.xml create mode 100644 DataFormats/HGCalReco/src/alpaka/classes_rocm.h create mode 100644 DataFormats/HGCalReco/src/alpaka/classes_rocm_def.xml mode change 100755 => 100644 DataFormats/Portable/scripts/portableHostCollectionHints create mode 100644 HLTrigger/Configuration/python/HLT_75e33/modules/hgcalLayerClustersFromSoAProducer_cfi.py create mode 100644 HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoALayerClustersProducer_cfi.py create mode 100644 HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoARecHitsLayerClustersProducer_cfi.py create mode 100644 HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoARecHitsProducer_cfi.py create mode 100644 RecoLocalCalo/HGCalRecProducers/interface/DumpClustersDetails.h create mode 100644 RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterHeterogeneousDumper.cc create mode 100644 RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterHeterogeneousSoADumper.cc create mode 100644 RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClustersFromSoAProducer.cc create mode 100644 RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.dev.cc create mode 100644 RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.h create mode 100644 RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.dev.cc create mode 100644 RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.h create mode 100644 RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoALayerClustersProducer.cc create mode 100644 RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsLayerClustersProducer.cc create mode 100644 RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsProducer.cc diff --git a/Configuration/ProcessModifiers/python/heterogeneousCLUE_cff.py b/Configuration/ProcessModifiers/python/heterogeneousCLUE_cff.py new file mode 100644 index 0000000000000..42e51034fb7b2 --- /dev/null +++ b/Configuration/ProcessModifiers/python/heterogeneousCLUE_cff.py @@ -0,0 +1,5 @@ +import FWCore.ParameterSet.Config as cms + +# This modifier runs the heterogeneous version of CLUE on the EE part of HGCAL +# in the HLT Phase2 Simplified menu. All seeded HLT Paths will keep on using the regular CPU version. +heterogeneousCLUE = cms.Modifier() diff --git a/DataFormats/HGCalReco/BuildFile.xml b/DataFormats/HGCalReco/BuildFile.xml index ea59e3506573a..a13c487df9bc3 100644 --- a/DataFormats/HGCalReco/BuildFile.xml +++ b/DataFormats/HGCalReco/BuildFile.xml @@ -5,6 +5,10 @@ + + + + diff --git a/DataFormats/HGCalReco/interface/HGCalSoACells.h b/DataFormats/HGCalReco/interface/HGCalSoACells.h new file mode 100644 index 0000000000000..57ee3e6549e1f --- /dev/null +++ b/DataFormats/HGCalReco/interface/HGCalSoACells.h @@ -0,0 +1,27 @@ +#ifndef DataFormats_HGCalReco_interface_HGCalSoACells_h +#define DataFormats_HGCalReco_interface_HGCalSoACells_h + +#include +#include + +#include "DataFormats/SoATemplate/interface/SoACommon.h" +#include "DataFormats/SoATemplate/interface/SoALayout.h" +#include "DataFormats/SoATemplate/interface/SoAView.h" + +// SoA layout with dim1, dim2, weight, sigmaNoise, recHitsIndex layer and cellsCount fields +GENERATE_SOA_LAYOUT(HGCalSoACellsLayout, + // columns: one value per element + SOA_COLUMN(float, dim1), + SOA_COLUMN(float, dim2), + SOA_COLUMN(float, dim3), + SOA_COLUMN(int, layer), + SOA_COLUMN(float, weight), + SOA_COLUMN(float, sigmaNoise), + SOA_COLUMN(unsigned int, recHitIndex), + SOA_COLUMN(uint32_t, detid), + SOA_COLUMN(float, time), + SOA_COLUMN(float, timeError)) + +using HGCalSoACells = HGCalSoACellsLayout<>; + +#endif // DataFormats_PortableTestObjects_interface_TestSoA_h diff --git a/DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h b/DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h new file mode 100644 index 0000000000000..e579847724862 --- /dev/null +++ b/DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h @@ -0,0 +1,10 @@ +#ifndef DataFormats_HGCalReco_interface_HGCalSoACellsHostCollection_h +#define DataFormats_HGCalReco_interface_HGCalSoACellsHostCollection_h + +#include "DataFormats/Portable/interface/PortableHostCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACells.h" + +// SoA with x, y, z, id fields in host memory +using HGCalSoACellsHostCollection = PortableHostCollection; + +#endif // DataFormats_HGCalReco_interface_HGCalSoACellsHostCollection_h diff --git a/DataFormats/HGCalReco/interface/HGCalSoACellsOut.h b/DataFormats/HGCalReco/interface/HGCalSoACellsOut.h new file mode 100644 index 0000000000000..cfefed2d844b0 --- /dev/null +++ b/DataFormats/HGCalReco/interface/HGCalSoACellsOut.h @@ -0,0 +1,20 @@ +#ifndef DataFormats_HGCalReco_interface_HGCalSoACellsOut_h +#define DataFormats_HGCalReco_interface_HGCalSoACellsOut_h + +#include "DataFormats/SoATemplate/interface/SoACommon.h" +#include "DataFormats/SoATemplate/interface/SoALayout.h" +#include "DataFormats/SoATemplate/interface/SoAView.h" + +// SoA layout with delta, rho, weight, nearestHigher, clusterIndex, layer, isSeed, and cellsCount fields +GENERATE_SOA_LAYOUT(HGCalSoACellsOutLayout, + // columns: one value per element + SOA_COLUMN(float, delta), + SOA_COLUMN(float, rho), + SOA_COLUMN(unsigned int, nearestHigher), + SOA_COLUMN(int, clusterIndex), + SOA_COLUMN(uint8_t, isSeed), + SOA_SCALAR(unsigned int, numberOfClustersScalar)) + +using HGCalSoACellsOut = HGCalSoACellsOutLayout<>; + +#endif diff --git a/DataFormats/HGCalReco/interface/HGCalSoACellsOutHostCollection.h b/DataFormats/HGCalReco/interface/HGCalSoACellsOutHostCollection.h new file mode 100644 index 0000000000000..1dc6014fc170f --- /dev/null +++ b/DataFormats/HGCalReco/interface/HGCalSoACellsOutHostCollection.h @@ -0,0 +1,10 @@ +#ifndef DataFormats_HGCalReco_interface_HGCalSoACellsOutHostCollection_h +#define DataFormats_HGCalReco_interface_HGCalSoACellsOutHostCollection_h + +#include "DataFormats/Portable/interface/PortableHostCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACellsOut.h" + +// SoA with delta, rho, weight, nearestHigher, clusterIndex, layer, isSeed, and cellsCount fields in host memory +using HGCalSoACellsOutHostCollection = PortableHostCollection; + +#endif // DataFormats_HGCalReco_interface_HGCalSoACellsOutHostCollection_h diff --git a/DataFormats/HGCalReco/interface/HGCalSoAClusters.h b/DataFormats/HGCalReco/interface/HGCalSoAClusters.h new file mode 100644 index 0000000000000..65466663dea85 --- /dev/null +++ b/DataFormats/HGCalReco/interface/HGCalSoAClusters.h @@ -0,0 +1,20 @@ +#ifndef DataFormats_HGCalReco_interface_HGCalSoAClusters_h +#define DataFormats_HGCalReco_interface_HGCalSoAClusters_h + +#include "DataFormats/SoATemplate/interface/SoACommon.h" +#include "DataFormats/SoATemplate/interface/SoALayout.h" +#include "DataFormats/SoATemplate/interface/SoAView.h" + +GENERATE_SOA_LAYOUT(HGCalSoAClustersLayout, + // columns: one value per element + SOA_COLUMN(float, x), + SOA_COLUMN(float, y), + SOA_COLUMN(float, z), + SOA_COLUMN(float, energy), + SOA_COLUMN(int, cells), // number of hits in the cluster + SOA_COLUMN(int, seed) // This is the index of the seed of each cluster inside the RecHit SoA +) + +using HGCalSoAClusters = HGCalSoAClustersLayout<>; + +#endif diff --git a/DataFormats/HGCalReco/interface/HGCalSoAClustersHostCollection.h b/DataFormats/HGCalReco/interface/HGCalSoAClustersHostCollection.h new file mode 100644 index 0000000000000..c88e021832c40 --- /dev/null +++ b/DataFormats/HGCalReco/interface/HGCalSoAClustersHostCollection.h @@ -0,0 +1,9 @@ +#ifndef DataFormats_HGCalReco_interface_HGCalSoAClustersHostCollection_h +#define DataFormats_HGCalReco_interface_HGCalSoAClustersHostCollection_h + +#include "DataFormats/Portable/interface/PortableHostCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoAClusters.h" + +using HGCalSoAClustersHostCollection = PortableHostCollection; + +#endif // DataFormats_HGCalReco_interface_HGCalSoAClustersHostCollection_h diff --git a/DataFormats/HGCalReco/interface/HGCalSoAClustersService.h b/DataFormats/HGCalReco/interface/HGCalSoAClustersService.h new file mode 100644 index 0000000000000..c4dce7799956b --- /dev/null +++ b/DataFormats/HGCalReco/interface/HGCalSoAClustersService.h @@ -0,0 +1,19 @@ +#ifndef DataFormats_HGCalReco_interface_HGCalSoAClustersService_h +#define DataFormats_HGCalReco_interface_HGCalSoAClustersService_h + +#include "DataFormats/SoATemplate/interface/SoACommon.h" +#include "DataFormats/SoATemplate/interface/SoALayout.h" +#include "DataFormats/SoATemplate/interface/SoAView.h" + +GENERATE_SOA_LAYOUT( + HGCalSoAClustersServiceLayout, + // columns: one value per element + SOA_COLUMN(float, total_weight), + SOA_COLUMN(float, total_weight_log), + SOA_COLUMN(float, maxEnergyValue), + SOA_COLUMN(int, maxEnergyIndex) // Index in the RecHitSoA of the rechit with highest energy in each cluster +) + +using HGCalSoAClustersService = HGCalSoAClustersServiceLayout<>; + +#endif diff --git a/DataFormats/HGCalReco/interface/HGCalSoAClustersServiceHostCollection.h b/DataFormats/HGCalReco/interface/HGCalSoAClustersServiceHostCollection.h new file mode 100644 index 0000000000000..bd5338b32e9fc --- /dev/null +++ b/DataFormats/HGCalReco/interface/HGCalSoAClustersServiceHostCollection.h @@ -0,0 +1,9 @@ +#ifndef DataFormats_HGCalReco_interface_HGCalSoAClustersServiceHostCollection_h +#define DataFormats_HGCalReco_interface_HGCalSoAClustersServiceHostCollection_h + +#include "DataFormats/Portable/interface/PortableHostCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoAClustersService.h" + +using HGCalSoAClustersServiceHostCollection = PortableHostCollection; + +#endif // DataFormats_HGCalReco_interface_HGCalSoAClustersServiceHostCollection_h diff --git a/DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsDeviceCollection.h b/DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsDeviceCollection.h new file mode 100644 index 0000000000000..2c20503dc7fcb --- /dev/null +++ b/DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsDeviceCollection.h @@ -0,0 +1,15 @@ +#ifndef DataFormats_PortableTestObjects_interface_alpaka_HGCalSoACellsDeviceCollection_h +#define DataFormats_PortableTestObjects_interface_alpaka_HGCalSoACellsDeviceCollection_h + +#include "DataFormats/Portable/interface/alpaka/PortableCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACells.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + // SoA with x, y, z, id fields in device global memory + using HGCalSoACellsDeviceCollection = PortableCollection; + +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#endif // DataFormats_PortableTestObjects_interface_alpaka_HGCalSoACellsDeviceCollection_h diff --git a/DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsOutDeviceCollection.h b/DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsOutDeviceCollection.h new file mode 100644 index 0000000000000..1122ddd404962 --- /dev/null +++ b/DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsOutDeviceCollection.h @@ -0,0 +1,15 @@ +#ifndef DataFormats_PortableTestObjects_interface_alpaka_HGCalSoACellsOutDeviceCollection_h +#define DataFormats_PortableTestObjects_interface_alpaka_HGCalSoACellsOutDeviceCollection_h + +#include "DataFormats/Portable/interface/alpaka/PortableCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACellsOut.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + // SoA with delta, rho, weight, nearestHigher, clusterIndex, layer, isSeed, and cellsCount fields in device global memory + using HGCalSoACellsOutDeviceCollection = PortableCollection; + +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#endif // DataFormats_PortableTestObjects_interface_alpaka_HGCalSoACellsOutDeviceCollection_h diff --git a/DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersDeviceCollection.h b/DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersDeviceCollection.h new file mode 100644 index 0000000000000..72834691e5f2e --- /dev/null +++ b/DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersDeviceCollection.h @@ -0,0 +1,14 @@ +#ifndef DataFormats_PortableTestObjects_interface_alpaka_HGCalSoAClustersDeviceCollection_h +#define DataFormats_PortableTestObjects_interface_alpaka_HGCalSoAClustersDeviceCollection_h + +#include "DataFormats/Portable/interface/alpaka/PortableCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoAClusters.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + using HGCalSoAClustersDeviceCollection = PortableCollection; + +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#endif // DataFormats_PortableTestObjects_interface_alpaka_HGCalSoAClustersDeviceCollection_h diff --git a/DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersServiceDeviceCollection.h b/DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersServiceDeviceCollection.h new file mode 100644 index 0000000000000..12e2031427df0 --- /dev/null +++ b/DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersServiceDeviceCollection.h @@ -0,0 +1,14 @@ +#ifndef DataFormats_PortableTestObjects_interface_alpaka_HGCalSoAClustersServiceDeviceCollection_h +#define DataFormats_PortableTestObjects_interface_alpaka_HGCalSoAClustersServiceDeviceCollection_h + +#include "DataFormats/Portable/interface/alpaka/PortableCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoAClustersService.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + using HGCalSoAClustersServiceDeviceCollection = PortableCollection; + +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#endif // DataFormats_PortableTestObjects_interface_alpaka_HGCalSoAClustersServiceDeviceCollection_h diff --git a/DataFormats/HGCalReco/src/alpaka/classes_cuda.h b/DataFormats/HGCalReco/src/alpaka/classes_cuda.h new file mode 100644 index 0000000000000..e191a760f6aa0 --- /dev/null +++ b/DataFormats/HGCalReco/src/alpaka/classes_cuda.h @@ -0,0 +1,7 @@ +#include "DataFormats/Common/interface/DeviceProduct.h" +#include "DataFormats/Common/interface/Wrapper.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACells.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACellsOut.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsOutDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersDeviceCollection.h" diff --git a/DataFormats/HGCalReco/src/alpaka/classes_cuda_def.xml b/DataFormats/HGCalReco/src/alpaka/classes_cuda_def.xml new file mode 100644 index 0000000000000..f186fe5755633 --- /dev/null +++ b/DataFormats/HGCalReco/src/alpaka/classes_cuda_def.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/DataFormats/HGCalReco/src/alpaka/classes_rocm.h b/DataFormats/HGCalReco/src/alpaka/classes_rocm.h new file mode 100644 index 0000000000000..e191a760f6aa0 --- /dev/null +++ b/DataFormats/HGCalReco/src/alpaka/classes_rocm.h @@ -0,0 +1,7 @@ +#include "DataFormats/Common/interface/DeviceProduct.h" +#include "DataFormats/Common/interface/Wrapper.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACells.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACellsOut.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsOutDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersDeviceCollection.h" diff --git a/DataFormats/HGCalReco/src/alpaka/classes_rocm_def.xml b/DataFormats/HGCalReco/src/alpaka/classes_rocm_def.xml new file mode 100644 index 0000000000000..b5bc5ed6ee095 --- /dev/null +++ b/DataFormats/HGCalReco/src/alpaka/classes_rocm_def.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/DataFormats/HGCalReco/src/classes.h b/DataFormats/HGCalReco/src/classes.h index d871bfb485a71..12468224ce5a7 100644 --- a/DataFormats/HGCalReco/src/classes.h +++ b/DataFormats/HGCalReco/src/classes.h @@ -5,3 +5,9 @@ #include "DataFormats/HGCalReco/interface/TICLSeedingRegion.h" #include "DataFormats/HGCalReco/interface/TICLCandidate.h" #include "DataFormats/Common/interface/Wrapper.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACells.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACellsOut.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACellsOutHostCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoAClusters.h" +#include "DataFormats/HGCalReco/interface/HGCalSoAClustersHostCollection.h" diff --git a/DataFormats/HGCalReco/src/classes_def.xml b/DataFormats/HGCalReco/src/classes_def.xml index 18c2cf05b934d..37539d3e9b8ea 100644 --- a/DataFormats/HGCalReco/src/classes_def.xml +++ b/DataFormats/HGCalReco/src/classes_def.xml @@ -55,4 +55,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/DataFormats/Portable/scripts/portableHostCollectionHints b/DataFormats/Portable/scripts/portableHostCollectionHints old mode 100755 new mode 100644 index d92e9cb2f132e..6e10ffc4cc5d1 --- a/DataFormats/Portable/scripts/portableHostCollectionHints +++ b/DataFormats/Portable/scripts/portableHostCollectionHints @@ -29,5 +29,25 @@ if len(layouts) > 1: print("") print(" ") print(" "% collectionName) +print() +print(" ") +print(" ") +print(" ") +print(" \n") print(" \" splitLevel=\"0\"/>"% collectionName) -print("") +print("") \ No newline at end of file diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalLayerClustersFromSoAProducer_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalLayerClustersFromSoAProducer_cfi.py new file mode 100644 index 0000000000000..73fc7bad201ff --- /dev/null +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalLayerClustersFromSoAProducer_cfi.py @@ -0,0 +1,14 @@ +import FWCore.ParameterSet.Config as cms +from ..psets.hgcal_reco_constants_cfi import HGCAL_reco_constants as HGCAL_reco_constants + +hgCalLayerClustersFromSoAProducer = cms.EDProducer("HGCalLayerClustersFromSoAProducer", + detector = cms.string('EE'), + hgcalRecHitsLayerClustersSoA = cms.InputTag("hgCalSoARecHitsLayerClustersProducer"), + hgcalRecHitsSoA = cms.InputTag("hgCalSoARecHitsProducer"), + mightGet = cms.optional.untracked.vstring, + nHitsTime = cms.uint32(3), + srcDevice = cms.InputTag("hgCalSoALayerClustersProducer"), + timeClname = cms.string('timeLayerCluster') +) + + diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoALayerClustersProducer_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoALayerClustersProducer_cfi.py new file mode 100644 index 0000000000000..573967558fda3 --- /dev/null +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoALayerClustersProducer_cfi.py @@ -0,0 +1,13 @@ +import FWCore.ParameterSet.Config as cms +from ..psets.hgcal_reco_constants_cfi import HGCAL_reco_constants as HGCAL_reco_constants + +hgCalSoALayerClustersProducer = cms.EDProducer("HGCalSoALayerClustersProducer@alpaka", + alpaka = cms.untracked.PSet( + backend = cms.untracked.string('cuda_async') + ), + hgcalRecHitsLayerClustersSoA = cms.InputTag("hgCalSoARecHitsLayerClustersProducer"), + hgcalRecHitsSoA = cms.InputTag("hgCalSoARecHitsProducer"), + mightGet = cms.optional.untracked.vstring, + positionDeltaRho2 = cms.double(1.69), + thresholdW0 = cms.double(2.9) +) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoARecHitsLayerClustersProducer_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoARecHitsLayerClustersProducer_cfi.py new file mode 100644 index 0000000000000..ec9a90ff64b75 --- /dev/null +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoARecHitsLayerClustersProducer_cfi.py @@ -0,0 +1,10 @@ +import FWCore.ParameterSet.Config as cms +from ..psets.hgcal_reco_constants_cfi import HGCAL_reco_constants as HGCAL_reco_constants + +hgCalSoARecHitsLayerClustersProducer = cms.EDProducer("HGCalSoARecHitsLayerClustersProducer@alpaka", + alpaka = cms.untracked.PSet( + backend = cms.untracked.string('cuda_async') + ), + hgcalRecHitsSoA = cms.InputTag("hgCalSoARecHitsProducer"), + mightGet = cms.optional.untracked.vstring +) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoARecHitsProducer_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoARecHitsProducer_cfi.py new file mode 100644 index 0000000000000..2f56f4db85d80 --- /dev/null +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/hgcalSoARecHitsProducer_cfi.py @@ -0,0 +1,19 @@ +import FWCore.ParameterSet.Config as cms +from ..psets.hgcal_reco_constants_cfi import HGCAL_reco_constants as HGCAL_reco_constants + +hgCalSoARecHitsProducer = cms.EDProducer("HGCalSoARecHitsProducer@alpaka", + alpaka = cms.untracked.PSet( + backend = cms.untracked.string('cuda_async') + ), + dEdXweights = HGCAL_reco_constants.dEdXweights, + detector = cms.string('EE'), + ecut = cms.double(3), + fcPerEle = HGCAL_reco_constants.fcPerEle, + fcPerMip = HGCAL_reco_constants.fcPerMip, + maxNumberOfThickIndices = HGCAL_reco_constants.maxNumberOfThickIndices, + mightGet = cms.optional.untracked.vstring, + noises = HGCAL_reco_constants.noises, + recHits = cms.InputTag("HGCalRecHit","HGCEERecHits"), + thicknessCorrection = HGCAL_reco_constants.thicknessCorrection, +) + diff --git a/HLTrigger/Configuration/python/HLT_75e33/sequences/hgcalLocalRecoSequence_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/sequences/hgcalLocalRecoSequence_cfi.py index 40cff277bb2a0..6fad347d2c262 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/sequences/hgcalLocalRecoSequence_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/sequences/hgcalLocalRecoSequence_cfi.py @@ -6,5 +6,33 @@ from ..modules.hgcalMergeLayerClusters_cfi import * from ..modules.HGCalRecHit_cfi import * from ..modules.HGCalUncalibRecHit_cfi import * +# Heterogeneous HGCAL EE layer clusters +from ..modules.hgcalSoARecHitsProducer_cfi import * +from ..modules.hgcalSoARecHitsLayerClustersProducer_cfi import * +from ..modules.hgcalSoALayerClustersProducer_cfi import * +from ..modules.hgcalLayerClustersFromSoAProducer_cfi import * -hgcalLocalRecoSequence = cms.Sequence(HGCalUncalibRecHit+HGCalRecHit+hgcalLayerClustersEE+hgcalLayerClustersHSci+hgcalLayerClustersHSi+hgcalMergeLayerClusters) +hgcalLocalRecoSequence = cms.Sequence( + HGCalUncalibRecHit+ + HGCalRecHit+ + hgcalLayerClustersEE+ + hgcalLayerClustersHSci+ + hgcalLayerClustersHSi+ + hgcalMergeLayerClusters) + +_hgcalLocalRecoSequence_heterogeneous = cms.Sequence( + HGCalUncalibRecHit+ + HGCalRecHit+ + hgCalSoARecHitsProducer+ + hgCalSoARecHitsLayerClustersProducer+ + hgCalSoALayerClustersProducer+ + hgCalLayerClustersFromSoAProducer+ + hgcalLayerClustersHSci+ + hgcalLayerClustersHSi+ + hgcalMergeLayerClusters) + +from Configuration.ProcessModifiers.heterogeneousCLUE_cff import heterogeneousCLUE +heterogeneousCLUE.toReplaceWith(hgcalLocalRecoSequence, _hgcalLocalRecoSequence_heterogeneous) +heterogeneousCLUE.toModify(hgcalMergeLayerClusters, + layerClustersEE = cms.InputTag("hgCalLayerClustersFromSoAProducer"), + time_layerclustersEE = cms.InputTag("hgCalLayerClustersFromSoAProducer", "timeLayerCluster")) diff --git a/HLTrigger/Configuration/python/HLT_75e33_timing_cff.py b/HLTrigger/Configuration/python/HLT_75e33_timing_cff.py index 52b6f8ab90b19..6c6cbd387ad36 100644 --- a/HLTrigger/Configuration/python/HLT_75e33_timing_cff.py +++ b/HLTrigger/Configuration/python/HLT_75e33_timing_cff.py @@ -2,6 +2,16 @@ fragment = cms.ProcessFragment("HLT") +### ACCELERATORS, in case they are needed +## Load explicitly +# One ProcessAccelerator for each accelerator technology +fragment.load("Configuration.StandardSequences.Accelerators_cff") + +# And one ProcessAccelerator for Alpaka +# (eventually to be absorbed to Accelerators_cff) +fragment.load("HeterogeneousCore.AlpakaCore.ProcessAcceleratorAlpaka_cfi") + + ### Non HLT-specific event-setups fragment.load("CalibMuon/CSCCalibration/CSCChannelMapper_cfi") fragment.load("CalibMuon/CSCCalibration/CSCIndexer_cfi") diff --git a/RecoLocalCalo/HGCalRecProducers/interface/DumpClustersDetails.h b/RecoLocalCalo/HGCalRecProducers/interface/DumpClustersDetails.h new file mode 100644 index 0000000000000..6996ef80ee231 --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/interface/DumpClustersDetails.h @@ -0,0 +1,186 @@ +#ifndef RecoLocalCalo_HGCalRecProducers_DumpClustersDetails_h +#define RecoLocalCalo_HGCalRecProducers_DumpClustersDetails_h + +#include +#include +#include +#include +#include +#include +#include + +namespace hgcalUtils { + + static bool sortByDetId(const std::pair& pair1, const std::pair& pair2) { + return pair1.first < pair2.first; + } + + class DumpClustersDetails { + public: + DumpClustersDetails(){}; + + template + void dumpInfos(const T& clusters, bool dumpCellsDetId = false) const { + // Get the process ID + pid_t pid = getpid(); + + // Create the filename using the PID + std::ostringstream filename; + // Seed the random number generator + srand(time(0)); + // Generate a random number between 100 and 999 + int random_number = 100 + rand() % 900; + filename << "CLUSTERS_" << pid << "_" << random_number << ".txt"; + // Open the file + std::ofstream outfile(filename.str()); + int count = 0; + for (auto& i : clusters) { + outfile << fmt::format( + "Seed: {}, Idx: {}, energy: {:.{}f}, x: {:.{}f}, y: {:.{}f}, z: {:.{}f}, eta: {:.{}f}, phi: {:.{}f}", + i.seed().rawId(), + count++, + (float)i.energy(), + std::numeric_limits::max_digits10, + i.x(), + std::numeric_limits::max_digits10, + i.y(), + std::numeric_limits::max_digits10, + i.z(), + std::numeric_limits::max_digits10, + i.eta(), + std::numeric_limits::max_digits10, + i.phi(), + std::numeric_limits::max_digits10); + if (dumpCellsDetId) { + auto sorted = i.hitsAndFractions(); // copy... + std::stable_sort(std::begin(sorted), std::end(sorted), sortByDetId); + for (auto const& c : sorted) { + outfile << fmt::format(" ({}, {:.{}f})", c.first, c.second, std::numeric_limits::max_digits10); + } // loop on hits and fractions + } else { + outfile << fmt::format(" ({} cells)", i.hitsAndFractions().size()); + } + outfile << std::endl; + } // loop on clusters + outfile.close(); + } + }; + + class DumpClustersSoADetails { + public: + DumpClustersSoADetails(){}; + + template + void dumpInfos(const T& clustersSoA) const { + // Get the process ID + pid_t pid = getpid(); + + // Create the filename using the PID + std::ostringstream filename; + filename << "CLUSTERS_UTILS_SOA_" << pid << ".txt"; + // Open the file + std::ofstream outfile(filename.str()); + for (int i = 0; i < clustersSoA->metadata().size(); ++i) { + auto clusterSoAV = clustersSoA.view()[i]; + outfile << fmt::format("Idx: {}, delta: {:.{}f}, rho: {:.{}f}, nearest: {}, clsIdx: {}, isSeed: {}", + i, + clusterSoAV.delta(), + std::numeric_limits::max_digits10, + clusterSoAV.rho(), + std::numeric_limits::max_digits10, + clusterSoAV.nearestHigher(), + clusterSoAV.clusterIndex(), + clusterSoAV.isSeed()) + << std::endl; + } + outfile.close(); + } + }; + + class DumpCellsSoADetails { + public: + DumpCellsSoADetails(){}; + + template + void dumpInfos(const T& cells) const { + // Get the process ID + pid_t pid = getpid(); + + // Create the filename using the PID + std::ostringstream filename; + filename << "RECHITS_SOA_" << pid << ".txt"; + // Open the file + std::ofstream outfile(filename.str()); + for (int i = 0; i < cells->metadata().size(); ++i) { + auto cellSoAV = cells.view()[i]; + outfile + << fmt::format( + "Idx Cell: {}, x: {:.{}f}, y: {:.{}f}, layer: {}, weight: {:.{}f}, sigmaNoise: {:.{}f}, detid: {}", + i, + (cellSoAV.dim1()), + std::numeric_limits::max_digits10, + (cellSoAV.dim2()), + std::numeric_limits::max_digits10, + cellSoAV.layer(), + (cellSoAV.weight()), + std::numeric_limits::max_digits10, + (cellSoAV.sigmaNoise()), + std::numeric_limits::max_digits10, + cellSoAV.detid()) + << std::endl; + } + outfile.close(); + } + }; + + class DumpLegacySoADetails { + public: + DumpLegacySoADetails(){}; + + template + void dumpInfos(T& cells) const { + // Get the process ID + pid_t pid = getpid(); + + // Create the filename using the PID + std::ostringstream filename; + filename << "RECHITS_LEGACY_" << pid << ".txt"; + // Open the file + std::ofstream outfile(filename.str()); + for (unsigned int l = 0; l < cells.size(); l++) { + for (unsigned int i = 0; i < cells.at(l).dim1.size(); ++i) { + outfile << fmt::format( + "Idx Cell: {}, x: {:.{}f}, y: {:.{}f}, layer: {}, weight: {:.{}f}, sigmaNoise: {:.{}f}, " + "delta: {:.{}f}, rho: {:.{}f}, nearestHigher: {}, clsIdx: {}, isSeed: {}, detid: {}", + i, + (cells.at(l).dim1.at(i)), + std::numeric_limits::max_digits10, + (cells.at(l).dim2.at(i)), + std::numeric_limits::max_digits10, + l, + (cells.at(l).weight.at(i)), + std::numeric_limits::max_digits10, + (cells.at(l).sigmaNoise.at(i)), + std::numeric_limits::max_digits10, + cells.at(l).delta.at(i), + std::numeric_limits::max_digits10, + cells.at(l).rho.at(i), + std::numeric_limits::max_digits10, + cells.at(l).nearestHigher.at(i), + cells.at(l).clusterIndex.at(i), + cells.at(l).isSeed.at(i), + cells.at(l).detid.at(i)) + << std::endl; + } + } + outfile.close(); + } + }; + + using DumpClusters = DumpClustersDetails; + using DumpClustersSoA = DumpClustersSoADetails; + using DumpCellsSoA = DumpCellsSoADetails; + using DumpLegacySoA = DumpLegacySoADetails; +} // namespace hgcalUtils + +#endif diff --git a/RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h b/RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h index d06b452b6b1d1..a45e33eefca17 100644 --- a/RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h +++ b/RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h @@ -20,11 +20,11 @@ class HGCalLayerTilesT { public: typedef T type; /** - * @brief fill the tile - * + * @brief fill the tile + * * @param[in] dim1 represents x or eta * @param[in] dim2 represents y or phils - * + * */ void fill(const std::vector& dim1, const std::vector& dim2) { auto cellsSize = dim1.size(); @@ -33,9 +33,9 @@ class HGCalLayerTilesT { tiles_[idx].push_back(i); } } - /** + /** * @brief compute bin for dim1 (x or eta) - * + * * @param[in] dim for bining * @return computed bin */ @@ -48,9 +48,9 @@ class HGCalLayerTilesT { return dimBin; } - /** + /** * @brief compute bin for dim2 (y or phi) - * + * * @param[in] dim for bining * @return computed bin */ @@ -78,8 +78,10 @@ class HGCalLayerTilesT { float d2 = dim2Cell1 - dim2Cell2; if constexpr (std::is_same_v) { d2 = reco::deltaPhi(dim2Cell1, dim2Cell2); + d2 *= d2; + return std::fmaf(d1, d1, d2); } - return (d1 * d1 + d2 * d2); + return std::fmaf(d1, d1, d2 * d2); } int getGlobalBin(float dim1, float dim2) const { return getDim1Bin(dim1) + getDim2Bin(dim2) * T::nColumns; } diff --git a/RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h b/RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h index 29a65d52a4d2a..7ec24cb72ec01 100644 --- a/RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h +++ b/RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h @@ -19,6 +19,9 @@ struct HGCalSiliconTilesConstants { static constexpr int nColumns = reco::ceil((maxDim1 - minDim1) / tileSize); static constexpr int nRows = reco::ceil((maxDim2 - minDim2) / tileSize); static constexpr int nTiles = nColumns * nRows; + static constexpr float invDim1BinSize = nColumns / (maxDim1 - minDim1); + static constexpr float invDim2BinSize = nRows / (maxDim2 - minDim2); + static constexpr int maxTileDepth = 64; // For accelerators. }; struct HGCalScintillatorTilesConstants { @@ -30,6 +33,9 @@ struct HGCalScintillatorTilesConstants { static constexpr int nColumns = reco::ceil((maxDim1 - minDim1) / tileSize); static constexpr int nRows = reco::ceil(2. * M_PI / tileSize); static constexpr int nTiles = nColumns * nRows; + static constexpr float invDim1BinSize = nColumns / (maxDim1 - minDim1); + static constexpr float invDim2BinSize = nRows / (maxDim2 - minDim2); + static constexpr int maxTileDepth = 32; // For accelerators. }; #endif diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/BuildFile.xml b/RecoLocalCalo/HGCalRecProducers/plugins/BuildFile.xml index cfddea6d9ae20..d8672094f0d9c 100644 --- a/RecoLocalCalo/HGCalRecProducers/plugins/BuildFile.xml +++ b/RecoLocalCalo/HGCalRecProducers/plugins/BuildFile.xml @@ -1,25 +1,52 @@ - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.cc b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.cc index 5c62c7ac4391b..cdbebfc023488 100644 --- a/RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.cc +++ b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.cc @@ -14,6 +14,12 @@ #include #include "DataFormats/DetId/interface/DetId.h" +#define DEBUG_CLUSTERS_ALPAKA 0 + +#if DEBUG_CLUSTERS_ALPAKA +#include "RecoLocalCalo/HGCalRecProducers/interface/DumpClustersDetails.h" +#endif + using namespace hgcal_clustering; template @@ -127,6 +133,10 @@ void HGCalCLUEAlgoT::makeClusters() { numberOfClustersPerLayer_[i] = findAndAssignClusters(i, delta); }); }); +#if DEBUG_CLUSTERS_ALPAKA + hgcalUtils::DumpLegacySoA dumperLegacySoA; + dumperLegacySoA.dumpInfos(cells_); +#endif } template @@ -215,7 +225,8 @@ std::vector HGCalCLUEAlgoT::getClusters(bool) { x *= inv_tot_weight; y *= inv_tot_weight; } else { - x = y = 0.f; + x = cellsOnLayer.dim1[maxEnergyCellIndex]; + y = cellsOnLayer.dim2[maxEnergyCellIndex]; } math::XYZPoint position = math::XYZPoint(x, y, z); @@ -348,6 +359,7 @@ void HGCalCLUEAlgoT::calculateDistanceToHigher(const T& lt, const u float maxDelta = std::numeric_limits::max(); float i_delta = maxDelta; int i_nearestHigher = -1; + float rho_max = 0.f; auto range = outlierDeltaFactor_ * delta; std::array search_box = lt.searchBox(cellsOnLayer.dim1[i] - range, cellsOnLayer.dim1[i] + range, @@ -366,14 +378,22 @@ void HGCalCLUEAlgoT::calculateDistanceToHigher(const T& lt, const u // loop over all hits in this bin for (unsigned int j = 0; j < binSize; j++) { unsigned int otherId = lt[binId][j]; - float dist = distance(lt, i, otherId, layerId); + float dist = distance2(lt, i, otherId, layerId); bool foundHigher = (cellsOnLayer.rho[otherId] > cellsOnLayer.rho[i]) || (cellsOnLayer.rho[otherId] == cellsOnLayer.rho[i] && cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]); - if (foundHigher && dist <= i_delta) { - // update i_delta + if (foundHigher && dist < i_delta) { + rho_max = cellsOnLayer.rho[otherId]; + i_delta = dist; + i_nearestHigher = otherId; + } else if (foundHigher && dist == i_delta && cellsOnLayer.rho[otherId] > rho_max) { + rho_max = cellsOnLayer.rho[otherId]; + i_delta = dist; + i_nearestHigher = otherId; + } else if (foundHigher && dist == i_delta && cellsOnLayer.rho[otherId] == rho_max && + cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]) { + rho_max = cellsOnLayer.rho[otherId]; i_delta = dist; - // update i_nearestHigher i_nearestHigher = otherId; } } @@ -381,7 +401,7 @@ void HGCalCLUEAlgoT::calculateDistanceToHigher(const T& lt, const u } bool foundNearestHigherInSearchBox = (i_delta != maxDelta); if (foundNearestHigherInSearchBox) { - cellsOnLayer.delta[i] = i_delta; + cellsOnLayer.delta[i] = std::sqrt(i_delta); cellsOnLayer.nearestHigher[i] = i_nearestHigher; } else { // otherwise delta is guaranteed to be larger outlierDeltaFactor_*delta_c diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.h b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.h index b940ec18786bc..cd099413ae3dc 100644 --- a/RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.h +++ b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.h @@ -198,6 +198,13 @@ class HGCalCLUEAlgoT : public HGCalClusteringAlgoBase { std::vector numberOfClustersPerLayer_; + inline float distance2(const TILE& lt, int cell1, int cell2, int layerId) const { // 2-d distance on the layer (x-y) + return (lt.distance2(cells_[layerId].dim1[cell1], + cells_[layerId].dim2[cell1], + cells_[layerId].dim1[cell2], + cells_[layerId].dim2[cell2])); + } + inline float distance(const TILE& lt, int cell1, int cell2, int layerId) const { // 2-d distance on the layer (x-y) return std::sqrt(lt.distance2(cells_[layerId].dim1[cell1], cells_[layerId].dim2[cell1], diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterHeterogeneousDumper.cc b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterHeterogeneousDumper.cc new file mode 100644 index 0000000000000..153ad8591166b --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterHeterogeneousDumper.cc @@ -0,0 +1,41 @@ +#include "DataFormats/HGCalReco/interface/HGCalSoACellsOutHostCollection.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDAnalyzer.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Utilities/interface/Exception.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include + +class HGCalLayerClusterHeterogeneousDumper : public edm::global::EDAnalyzer<> { +public: + HGCalLayerClusterHeterogeneousDumper(edm::ParameterSet const& iConfig) + : deviceToken_{consumes(iConfig.getParameter("srcDevice"))} {} + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("srcDevice", edm::InputTag("hgCalSoARecHitsLayerClustersProducer")); + descriptions.addWithDefaultLabel(desc); + } + + void analyze(edm::StreamID iStream, edm::Event const& iEvent, edm::EventSetup const& iSetup) const override { + auto const& deviceData = iEvent.get(deviceToken_); + + auto const deviceView = deviceData.view(); + std::cout << fmt::format("view.numberOfClusters() = {}", deviceView.numberOfClustersScalar()) << std::endl; + for (int i = 0; i < deviceData->metadata().size(); ++i) { + std::cout << fmt::format("view[{}].clusterIndex() = {}", i, deviceView.clusterIndex(i)) << std::endl; + } + } + +private: + edm::EDGetTokenT const deviceToken_; +}; + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(HGCalLayerClusterHeterogeneousDumper); diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterHeterogeneousSoADumper.cc b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterHeterogeneousSoADumper.cc new file mode 100644 index 0000000000000..0429c70dd10a2 --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterHeterogeneousSoADumper.cc @@ -0,0 +1,51 @@ +#include "DataFormats/HGCalReco/interface/HGCalSoAClustersHostCollection.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDAnalyzer.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Utilities/interface/Exception.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include + +class HGCalLayerClusterHeterogeneousSoADumper : public edm::global::EDAnalyzer<> { +public: + HGCalLayerClusterHeterogeneousSoADumper(edm::ParameterSet const& iConfig) + : deviceToken_{consumes(iConfig.getParameter("srcDevice"))} {} + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("srcDevice", edm::InputTag("hgCalSoALayerClustersProducer")); + descriptions.addWithDefaultLabel(desc); + } + + void analyze(edm::StreamID iStream, edm::Event const& iEvent, edm::EventSetup const& iSetup) const override { + auto const& deviceData = iEvent.get(deviceToken_); + + auto const deviceView = deviceData.view(); + std::cout << fmt::format("hgcalSoALayerClustersProducer size = {}", deviceView.metadata().size()) << std::endl; + for (int i = 0; i < deviceData->metadata().size(); ++i) { + std::cout << fmt::format("CLUSTERS_SOA {}, energy = {:.{}f}, x = {:.{}f}, y = {:.{}f}, z= {:.{}f}", + i, + deviceView.energy(i), + std::numeric_limits::max_digits10, + deviceView.x(i), + std::numeric_limits::max_digits10, + deviceView.y(i), + std::numeric_limits::max_digits10, + deviceView.z(i), + std::numeric_limits::max_digits10) + << std::endl; + } + } + +private: + edm::EDGetTokenT const deviceToken_; +}; + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(HGCalLayerClusterHeterogeneousSoADumper); diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterProducer.cc b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterProducer.cc index 742e388ac1589..b21fb07b7c904 100644 --- a/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterProducer.cc +++ b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterProducer.cc @@ -2,6 +2,8 @@ // Date: 03/2023 // @file create layer clusters +#define DEBUG_CLUSTERS_ALPAKA 0 + // user include files #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/stream/EDProducer.h" @@ -32,27 +34,31 @@ #include "FWCore/Framework/interface/ConsumesCollector.h" +#if DEBUG_CLUSTERS_ALPAKA +#include "RecoLocalCalo/HGCalRecProducers/interface/DumpClustersDetails.h" +#endif + class HGCalLayerClusterProducer : public edm::stream::EDProducer<> { public: /** * @brief Constructor with parameter settings - which can be changed in hgcalLayerCluster_cff.py. - * Constructor will set all variables by input param ps. + * Constructor will set all variables by input param ps. * algoID variables will be set accordingly to the detector type. - * + * * @param[in] ps parametr set to set variables */ HGCalLayerClusterProducer(const edm::ParameterSet&); ~HGCalLayerClusterProducer() override {} /** * @brief Method fill description which will be used in pyhton file. - * + * * @param[out] description to be fill */ static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); /** * @brief Method run the algoritm to get clusters. - * + * * @param[in, out] evt from get info and put result * @param[in] es to get event setup info */ @@ -83,7 +89,7 @@ class HGCalLayerClusterProducer : public edm::stream::EDProducer<> { /** * @brief Counts position for all points in the cluster - * + * * @param[in] hitmap hitmap to find correct RecHit * @param[in] hitsAndFraction all hits in the cluster * @return counted position @@ -93,7 +99,7 @@ class HGCalLayerClusterProducer : public edm::stream::EDProducer<> { /** * @brief Counts time for all points in the cluster - * + * * @param[in] hitmap hitmap to find correct RecHit only for silicon (not for BH-HSci) * @param[in] hitsAndFraction all hits in the cluster * @return counted time @@ -195,7 +201,7 @@ math::XYZPoint HGCalLayerClusterProducer::calculatePosition( float inv_tot_weight = 1.f / total_weight; return math::XYZPoint(x * inv_tot_weight, y * inv_tot_weight, positionMaxEnergy.z()); } else { - return math::XYZPoint(0.f, 0.f, 0.f); + return {positionMaxEnergy.x(), positionMaxEnergy.y(), positionMaxEnergy.z()}; } } @@ -263,6 +269,12 @@ void HGCalLayerClusterProducer::produce(edm::Event& evt, const edm::EventSetup& } } +#if DEBUG_CLUSTERS_ALPAKA + hgcalUtils::DumpClusters dumper; + + dumper.dumpInfos(*clusters, true); +#endif + auto clusterHandle = evt.put(std::move(clusters)); if (detector_ == "HFNose") { diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClustersFromSoAProducer.cc b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClustersFromSoAProducer.cc new file mode 100644 index 0000000000000..ddeb4c16a4893 --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClustersFromSoAProducer.cc @@ -0,0 +1,165 @@ +#include "FWCore/Framework/interface/stream/EDProducer.h" + +#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h" + +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Utilities/interface/EDGetToken.h" + +#include "DataFormats/Common/interface/ValueMap.h" +#include "DataFormats/Math/interface/Point3D.h" +#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h" +#include "DataFormats/HGCalReco/interface/HGCalSoAClustersHostCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACellsOutHostCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h" + +#include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h" + +#define DEBUG_CLUSTERS_ALPAKA 0 + +#if DEBUG_CLUSTERS_ALPAKA +#include "RecoLocalCalo/HGCalRecProducers/interface/DumpClustersDetails.h" +#endif + +class HGCalLayerClustersFromSoAProducer : public edm::stream::EDProducer<> { +public: + HGCalLayerClustersFromSoAProducer(edm::ParameterSet const& config) + : getTokenSoAClusters_(consumes(config.getParameter("srcDevice"))), + getTokenSoACells_(consumes(config.getParameter("hgcalRecHitsSoA"))), + getTokenSoACellsOut_(consumes(config.getParameter("hgcalRecHitsLayerClustersSoA"))), + detector_(config.getParameter("detector")), + hitsTime_(config.getParameter("nHitsTime")), + timeClname_(config.getParameter("timeClname")) { + if (detector_ == "HFNose") { + algoId_ = reco::CaloCluster::hfnose; + } else if (detector_ == "EE") { + algoId_ = reco::CaloCluster::hgcal_em; + } else { //for FH or BH + algoId_ = reco::CaloCluster::hgcal_had; + } + + produces>("InitialLayerClustersMask"); + produces>(); + produces>>(timeClname_); + } + + ~HGCalLayerClustersFromSoAProducer() override = default; + + void produce(edm::Event& iEvent, edm::EventSetup const& iSetup) override { + auto const& deviceData = iEvent.get(getTokenSoAClusters_); + + auto const& deviceSoACellsOut = iEvent.get(getTokenSoACellsOut_); + auto const soaCellsOut_v = deviceSoACellsOut.view(); + + auto const& deviceSoACells = iEvent.get(getTokenSoACells_); + auto const soaCells_v = deviceSoACells.view(); + + auto const deviceView = deviceData.view(); + + std::unique_ptr> clusters(new std::vector); + clusters->reserve(deviceData->metadata().size()); + + for (int i = 0; i < deviceData->metadata().size(); ++i) { + std::vector> thisCluster; + thisCluster.reserve(deviceView.cells(i)); + clusters->emplace_back(deviceView.energy(i), + math::XYZPoint(deviceView.x(i), deviceView.y(i), deviceView.z(i)), + reco::CaloID::DET_HGCAL_ENDCAP, + std::move(thisCluster), + algoId_); + clusters->back().setSeed(deviceView.seed(i)); + } + + // Populate hits and fractions required to compute the cluster's time. + // This procedure is complex and involves two SoAs: the original RecHits + // SoA and the clustering algorithm's output SoA. Both SoAs have the same + // cardinality, and crucially, the output SoA includes the cluster index. + + // Create a vector of locations, where each location holds a + // vector of 16 floats. These vectors are used to compute the time for + // each cluster. The size of 16 is a heuristic to minimize memory + // reallocation and the associated overhead when the vector's capacity + // needs to be extended. + std::vector> times(clusters->size(), std::vector(16, 0.0f)); + std::vector> timeErrors(clusters->size(), std::vector(16, 0.0f)); + for (int32_t i = 0; i < soaCellsOut_v.metadata().size(); ++i) { + if (soaCellsOut_v[i].clusterIndex() == -1) { + continue; + } + assert(soaCellsOut_v[i].clusterIndex() < (int)clusters->size()); + (*clusters)[soaCellsOut_v[i].clusterIndex()].addHitAndFraction(soaCells_v[i].detid(), 1.f); + if (soaCells_v[i].timeError() < 0.f) { + continue; + } + times[soaCellsOut_v[i].clusterIndex()].push_back(soaCells_v[i].time()); + timeErrors[soaCellsOut_v[i].clusterIndex()].push_back(1.f / + (soaCells_v[i].timeError() * soaCells_v[i].timeError())); + } + + // Finally, compute and assign the time to each cluster. + std::vector> cluster_times; + hgcalsimclustertime::ComputeClusterTime timeEstimator; + for (unsigned i = 0; i < clusters->size(); ++i) { + if (detector_ != "BH") { + cluster_times.push_back(timeEstimator.fixSizeHighestDensity(times[i], timeErrors[i], hitsTime_)); + } else { + cluster_times.push_back(std::pair(-99.f, -1.f)); + } + } + +#if DEBUG_CLUSTERS_ALPAKA + // hgcalUtils::DumpCellsSoA dumperCellsSoA; + // dumperCellsSoA.dumpInfos(cells); + + hgcalUtils::DumpClusters dumper; + dumper.dumpInfos(*clusters, true); + +// hgcalUtils::DumpClustersSoA dumperSoA; +// dumperSoA.dumpInfos(clustersSoA); +#endif + + auto clusterHandle = iEvent.put(std::move(clusters)); + + auto timeCl = std::make_unique>>(); + edm::ValueMap>::Filler filler(*timeCl); + filler.insert(clusterHandle, cluster_times.begin(), cluster_times.end()); + filler.fill(); + iEvent.put(std::move(timeCl), timeClname_); + + // The layerClusterMask for the HGCAL detector is created at a later + // stage, when the layer clusters from the different components of HGCAL + // are merged together into a unique collection. For the case of HFNose, + // since there is no further merging step needed, we create the + // layerClustersMask directly here. + if (detector_ == "HFNose") { + std::unique_ptr> layerClustersMask(new std::vector); + layerClustersMask->resize(clusterHandle->size(), 1.0); + iEvent.put(std::move(layerClustersMask), "InitialLayerClustersMask"); + } + } + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("srcDevice", edm::InputTag("hgCalSoALayerClustersProducer")); + desc.add("hgcalRecHitsLayerClustersSoA", edm::InputTag("TO BE DEFINED")); + desc.add("hgcalRecHitsSoA", edm::InputTag("TO BE DEFINED")); + desc.add("nHitsTime", 3); + desc.add("timeClname", "timeLayerCluster"); + desc.add("detector", "EE"); + descriptions.addWithDefaultLabel(desc); + } + +private: + edm::EDGetTokenT const getTokenSoAClusters_; + edm::EDGetTokenT const getTokenSoACells_; + edm::EDGetTokenT const getTokenSoACellsOut_; + std::string detector_; + unsigned int hitsTime_; + std::string timeClname_; + reco::CaloCluster::AlgoId algoId_; +}; +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(HGCalLayerClustersFromSoAProducer); diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/MergeClusterProducer.cc b/RecoLocalCalo/HGCalRecProducers/plugins/MergeClusterProducer.cc index 093ca0200cf02..29111af73285e 100644 --- a/RecoLocalCalo/HGCalRecProducers/plugins/MergeClusterProducer.cc +++ b/RecoLocalCalo/HGCalRecProducers/plugins/MergeClusterProducer.cc @@ -16,22 +16,22 @@ class MergeClusterProducer : public edm::stream::EDProducer<> { public: /** * @brief Constructor with parameter settings - which can be changed in ...todo. - * Constructor will set all variables by input param ps. - * + * Constructor will set all variables by input param ps. + * * @param[in] ps parametr set to set variables */ MergeClusterProducer(const edm::ParameterSet &); ~MergeClusterProducer() override {} /** * @brief Method fill description which will be used in pyhton file. - * + * * @param[out] description to be fill */ static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); /** * @brief Method will merge the producers and put them back to event - * + * * @param[in, out] evt from get info and put result * @param[in] es to get event setup info */ @@ -49,7 +49,7 @@ class MergeClusterProducer : public edm::stream::EDProducer<> { /** * @brief method merge three vectors of reco::CaloCluster to one - * + * * @param[out] merge the vector into which others vectors will be merge * @param[in] EE vector for Electromagnetic silicon * @param[in] HSi vector for Hardon silicon @@ -62,7 +62,7 @@ class MergeClusterProducer : public edm::stream::EDProducer<> { /** * @brief copy all values from vm to to - * + * * @param[in] vm Value map with values * @param[out] to vector to will be copy value map */ @@ -74,7 +74,7 @@ class MergeClusterProducer : public edm::stream::EDProducer<> { } /** * @brief Merge value map of time for all parts of detector together to vector times - * + * * @param[in] evt Event to get time value maps * @param[in] size of all 3 value maps * @param[out] times vector of merged time vectors @@ -93,9 +93,9 @@ class MergeClusterProducer : public edm::stream::EDProducer<> { } /** * @brief get info form event and then call merge - * + * * it is used for merge and clusters and time - * + * * @param[in] evt Event * @param[in] EE_token token for Electromagnetic silicon * @param[in] HSi_token token for Hardon silicon diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.dev.cc b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.dev.cc new file mode 100644 index 0000000000000..6239b1f2be437 --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.dev.cc @@ -0,0 +1,38 @@ +// Check that ALPAKA_HOST_ONLY is not defined during device compilation: +#ifdef ALPAKA_HOST_ONLY +#error ALPAKA_HOST_ONLY defined in device compilation +#endif + +#include "RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h" + +#include "HGCalLayerClustersAlgoWrapper.h" + +#include "CLUEAlgoAlpaka.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + using namespace cms::alpakatools; + + void HGCalLayerClustersAlgoWrapper::run(Queue& queue, + const unsigned int size, + const HGCalSoACellsDeviceCollection::ConstView inputs, + HGCalSoACellsOutDeviceCollection::View outputs) const { + CLUEAlgoAlpaka algoStandalone( + queue, 1.3f, 9.f, 2.f, false); + + algoStandalone.makeClustersCMSSW(size, + inputs.dim1(), + inputs.dim2(), + inputs.layer(), + inputs.weight(), + inputs.sigmaNoise(), + inputs.detid(), + outputs.rho(), + outputs.delta(), + outputs.nearestHigher(), + outputs.clusterIndex(), + outputs.isSeed(), + &outputs.numberOfClustersScalar()); + } + +} // namespace ALPAKA_ACCELERATOR_NAMESPACE diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.h b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.h new file mode 100644 index 0000000000000..50064d74aa45f --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.h @@ -0,0 +1,24 @@ +#ifndef RecoLocalCalo_HGCalRecProducers_plugins_alpaka_HGCalLayerClustersAlgoWrapper_h +#define RecoLocalCalo_HGCalRecProducers_plugins_alpaka_HGCalLayerClustersAlgoWrapper_h + +#include + +#include "DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsOutDeviceCollection.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "HeterogeneousCore/AlpakaInterface/interface/traits.h" +#include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + class HGCalLayerClustersAlgoWrapper { + public: + void run(Queue& queue, + const unsigned int size, + const HGCalSoACellsDeviceCollection::ConstView inputs, + HGCalSoACellsOutDeviceCollection::View outputs) const; + }; +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#endif // RecoLocalCalo_HGCalRecProducers_plugins_alpaka_HGCalLayerClustersAlgoWrapper_h diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.dev.cc b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.dev.cc new file mode 100644 index 0000000000000..bf8fc22477a0c --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.dev.cc @@ -0,0 +1,288 @@ +// Check that ALPAKA_HOST_ONLY is not defined during device compilation: +#ifdef ALPAKA_HOST_ONLY +#error ALPAKA_HOST_ONLY defined in device compilation +#endif + +#include "RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h" + +#include "HGCalLayerClustersSoAAlgoWrapper.h" + +#include "CLUEAlgoAlpaka.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + using namespace cms::alpakatools; + + // Set energy and number of hits in each clusters + class HGCalLayerClustersSoAAlgoKernelEnergy { + public: + template >> + ALPAKA_FN_ACC void operator()(TAcc const& acc, + const unsigned int numer_of_clusters, + const HGCalSoACellsDeviceCollection::ConstView input_rechits_soa, + const HGCalSoACellsOutDeviceCollection::ConstView input_clusters_soa, + HGCalSoAClustersDeviceCollection::View outputs) const { + // make a strided loop over the kernel grid, covering up to "size" elements + for (int32_t i : uniform_elements(acc, input_rechits_soa.metadata().size())) { + // Skip unassigned rechits + if (input_clusters_soa[i].clusterIndex() == -1) { + continue; + } + auto clIdx = input_clusters_soa[i].clusterIndex(); + alpaka::atomicAdd(acc, &outputs[clIdx].energy(), input_rechits_soa[i].weight()); + alpaka::atomicAdd(acc, &outputs[clIdx].cells(), 1); + } + } + }; + + class HGCalLayerClustersSoAAlgoKernelPosition { + public: + template >> + ALPAKA_FN_ACC void operator()(TAcc const& acc, + const unsigned int numer_of_clusters, + float thresholdW0, + float positionDeltaRho2, + const HGCalSoACellsDeviceCollection::ConstView input_rechits_soa, + const HGCalSoACellsOutDeviceCollection::ConstView input_clusters_soa, + HGCalSoAClustersDeviceCollection::View outputs) const { + // global index of the thread within the grid + const int32_t thread_idx = alpaka::getIdx(acc)[0u]; + + outputs[thread_idx].x() = 0.f; + outputs[thread_idx].y() = 0.f; + outputs[thread_idx].z() = 0.f; + float maxEnergyValue = 0.f; + float total_weight = 0.f; + unsigned int maxEnergyIndex = 0; + + for (auto hit = 0; hit < input_rechits_soa.metadata().size(); ++hit) { + if (input_clusters_soa[hit].clusterIndex() != thread_idx) { + continue; + } + total_weight += input_rechits_soa[hit].weight(); + if (input_rechits_soa[hit].weight() > maxEnergyValue) { + maxEnergyValue = input_rechits_soa[hit].weight(); + maxEnergyIndex = hit; + } + } + + float total_weight_log = 0.f; + float reference_x = input_rechits_soa[maxEnergyIndex].dim1(); + float reference_y = input_rechits_soa[maxEnergyIndex].dim2(); + float reference_z = input_rechits_soa[maxEnergyIndex].dim3(); + for (auto hit = 0; hit < input_rechits_soa.metadata().size(); ++hit) { + if (input_clusters_soa[hit].clusterIndex() != thread_idx) { + continue; + } + //for silicon only just use 1+6 cells = 1.3cm for all thicknesses + const float d1 = input_rechits_soa[hit].dim1() - reference_x; + const float d2 = input_rechits_soa[hit].dim2() - reference_y; + if ((d1 * d1 + d2 * d2) > positionDeltaRho2) { + continue; + } + float Wi = std::max(thresholdW0 + std::log(input_rechits_soa[hit].weight() / total_weight), 0.f); + outputs[thread_idx].x() += input_rechits_soa[hit].dim1() * Wi; + outputs[thread_idx].y() += input_rechits_soa[hit].dim2() * Wi; + total_weight_log += Wi; + } + total_weight = total_weight_log; + float inv_tot_weight = 1.f / total_weight; + outputs[thread_idx].x() *= inv_tot_weight; + outputs[thread_idx].y() *= inv_tot_weight; + outputs[thread_idx].z() = reference_z; + } + }; + + // Kernel to find the max for every cluster + class HGCalLayerClustersSoAAlgoKernelPositionByHits { + public: + template >> + ALPAKA_FN_ACC void operator()(TAcc const& acc, + const unsigned int numer_of_clusters, + float thresholdW0, + float positionDeltaRho2, + const HGCalSoACellsDeviceCollection::ConstView input_rechits_soa, + const HGCalSoACellsOutDeviceCollection::ConstView input_clusters_soa, + HGCalSoAClustersDeviceCollection::View outputs, + HGCalSoAClustersServiceDeviceCollection::View outputs_service) const { + // make a strided loop over the kernel grid, covering up to "size" elements + for (int32_t hit_index : uniform_elements(acc, input_rechits_soa.metadata().size())) { + const int cluster_index = input_clusters_soa[hit_index].clusterIndex(); + + // Bail out if you are not part of any cluster + if (cluster_index == -1) { + continue; + } + + alpaka::atomicAdd(acc, &outputs_service[cluster_index].total_weight(), input_rechits_soa[hit_index].weight()); + // Read the current seed index, and the associated energy. + int clusterSeed = outputs_service[cluster_index].maxEnergyIndex(); + float clusterEnergy = (clusterSeed == -1) ? 0. : input_rechits_soa[clusterSeed].weight(); + + while (input_rechits_soa[hit_index].weight() > clusterEnergy) { + // If output_service[cluster_index].maxEnergyIndex() did not change, + // store the new value and exit the loop. Otherwise return the value + // that has been updated, and decide again if the maximum needs to be + // updated. + int seed = alpaka::atomicCas(acc, &outputs_service[cluster_index].maxEnergyIndex(), clusterSeed, hit_index); + if (seed == hit_index) { + // atomicCas has stored the new value. + break; + } else { + // Update the seed index and re-read the associated energy. + clusterSeed = seed; + clusterEnergy = input_rechits_soa[clusterSeed].weight(); + } + } // CAS + } // uniform_elements + } // operator() + }; + + // Real Kernel position + class HGCalLayerClustersSoAAlgoKernelPositionByHits2 { + public: + template >> + ALPAKA_FN_ACC void operator()(TAcc const& acc, + const unsigned int numer_of_clusters, + float thresholdW0, + float positionDeltaRho2, + const HGCalSoACellsDeviceCollection::ConstView input_rechits_soa, + const HGCalSoACellsOutDeviceCollection::ConstView input_clusters_soa, + HGCalSoAClustersDeviceCollection::View outputs, + HGCalSoAClustersServiceDeviceCollection::View outputs_service) const { + // make a strided loop over the kernel grid, covering up to "size" elements + for (int32_t hit_index : uniform_elements(acc, input_rechits_soa.metadata().size())) { + const int cluster_index = input_clusters_soa[hit_index].clusterIndex(); + const int max_energy_index = outputs_service[cluster_index].maxEnergyIndex(); + + // Bail out if you are not part of any cluster + if (cluster_index == -1) { + continue; + } + + //for silicon only just use 1+6 cells = 1.3cm for all thicknesses + const float d1 = input_rechits_soa[hit_index].dim1() - input_rechits_soa[max_energy_index].dim1(); + const float d2 = input_rechits_soa[hit_index].dim2() - input_rechits_soa[max_energy_index].dim2(); + if ((d1 * d1 + d2 * d2) > positionDeltaRho2) { + continue; + } + float Wi = std::max(thresholdW0 + std::log(input_rechits_soa[hit_index].weight() / + outputs_service[cluster_index].total_weight()), + 0.f); + alpaka::atomicAdd(acc, &outputs[cluster_index].x(), input_rechits_soa[hit_index].dim1() * Wi); + alpaka::atomicAdd(acc, &outputs[cluster_index].y(), input_rechits_soa[hit_index].dim2() * Wi); + alpaka::atomicAdd(acc, &outputs_service[cluster_index].total_weight_log(), Wi); + } // uniform_elements + } // operator() + }; + + // Besides the final position, add also the DetId of the seed of each cluster + class HGCalLayerClustersSoAAlgoKernelPositionByHits3 { + public: + template >> + ALPAKA_FN_ACC void operator()(TAcc const& acc, + const unsigned int numer_of_clusters, + float thresholdW0, + float positionDeltaRho2, + const HGCalSoACellsDeviceCollection::ConstView input_rechits_soa, + const HGCalSoACellsOutDeviceCollection::ConstView input_clusters_soa, + HGCalSoAClustersDeviceCollection::View outputs, + HGCalSoAClustersServiceDeviceCollection::View outputs_service) const { + // make a strided loop over the kernel grid, covering up to "size" elements + for (int32_t cluster_index : uniform_elements(acc, outputs.metadata().size())) { + const int max_energy_index = outputs_service[cluster_index].maxEnergyIndex(); + + if (outputs_service[cluster_index].total_weight_log() > 0.f) { + float inv_tot_weight = 1.f / outputs_service[cluster_index].total_weight_log(); + outputs[cluster_index].x() *= inv_tot_weight; + outputs[cluster_index].y() *= inv_tot_weight; + } else { + outputs[cluster_index].x() = input_rechits_soa[max_energy_index].dim1(); + outputs[cluster_index].y() = input_rechits_soa[max_energy_index].dim2(); + } + outputs[cluster_index].z() = input_rechits_soa[max_energy_index].dim3(); + outputs[cluster_index].seed() = input_rechits_soa[max_energy_index].detid(); + } // uniform_elements + } // operator() + }; + + void HGCalLayerClustersSoAAlgoWrapper::run(Queue& queue, + const unsigned int size, + float thresholdW0, + float positionDeltaRho2, + const HGCalSoACellsDeviceCollection::ConstView input_rechits_soa, + const HGCalSoACellsOutDeviceCollection::ConstView input_clusters_soa, + HGCalSoAClustersDeviceCollection::View outputs, + HGCalSoAClustersServiceDeviceCollection::View outputs_service) const { + auto x = cms::alpakatools::make_device_view(alpaka::getDev(queue), outputs.x(), size); + alpaka::memset(queue, x, 0x0); + auto y = cms::alpakatools::make_device_view(alpaka::getDev(queue), outputs.y(), size); + alpaka::memset(queue, y, 0x0); + auto z = cms::alpakatools::make_device_view(alpaka::getDev(queue), outputs.z(), size); + alpaka::memset(queue, z, 0x0); + auto seed = cms::alpakatools::make_device_view(alpaka::getDev(queue), outputs.seed(), size); + alpaka::memset(queue, seed, 0x0); + auto energy = cms::alpakatools::make_device_view(alpaka::getDev(queue), outputs.energy(), size); + alpaka::memset(queue, energy, 0x0); + auto cells = cms::alpakatools::make_device_view(alpaka::getDev(queue), outputs.cells(), size); + alpaka::memset(queue, cells, 0x0); + auto total_weight = + cms::alpakatools::make_device_view(alpaka::getDev(queue), outputs_service.total_weight(), size); + alpaka::memset(queue, total_weight, 0x0); + auto total_weight_log = + cms::alpakatools::make_device_view(alpaka::getDev(queue), outputs_service.total_weight_log(), size); + alpaka::memset(queue, total_weight_log, 0x0); + auto maxEnergyValue = + cms::alpakatools::make_device_view(alpaka::getDev(queue), outputs_service.maxEnergyValue(), size); + alpaka::memset(queue, maxEnergyValue, 0x0); + auto maxEnergyIndex = + cms::alpakatools::make_device_view(alpaka::getDev(queue), outputs_service.maxEnergyIndex(), size); + alpaka::memset(queue, maxEnergyIndex, 0xff); + + // use 64 items per group (this value is arbitrary, but it's a reasonable starting point) + uint32_t items = 64; + + // use as many groups as needed to cover the whole problem + uint32_t groups = divide_up_by(input_rechits_soa.metadata().size(), items); + + // map items to + // - threads with a single element per thread on a GPU backend + // - elements within a single thread on a CPU backend + auto workDiv = make_workdiv(groups, items); + + alpaka::exec( + queue, workDiv, HGCalLayerClustersSoAAlgoKernelEnergy{}, size, input_rechits_soa, input_clusters_soa, outputs); + alpaka::exec(queue, + workDiv, + HGCalLayerClustersSoAAlgoKernelPositionByHits{}, + size, + thresholdW0, + positionDeltaRho2, + input_rechits_soa, + input_clusters_soa, + outputs, + outputs_service); + alpaka::exec(queue, + workDiv, + HGCalLayerClustersSoAAlgoKernelPositionByHits2{}, + size, + thresholdW0, + positionDeltaRho2, + input_rechits_soa, + input_clusters_soa, + outputs, + outputs_service); + uint32_t group_clusters = divide_up_by(size, items); + auto workDivClusters = make_workdiv(group_clusters, items); + alpaka::exec(queue, + workDivClusters, + HGCalLayerClustersSoAAlgoKernelPositionByHits3{}, + size, + thresholdW0, + positionDeltaRho2, + input_rechits_soa, + input_clusters_soa, + outputs, + outputs_service); + } +} // namespace ALPAKA_ACCELERATOR_NAMESPACE diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.h b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.h new file mode 100644 index 0000000000000..4b3606a668ce7 --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.h @@ -0,0 +1,30 @@ +#ifndef RecoLocalCalo_HGCalRecProducers_plugins_alpaka_HGCalLayerClustersSoAAlgoWrapper_h +#define RecoLocalCalo_HGCalRecProducers_plugins_alpaka_HGCalLayerClustersSoAAlgoWrapper_h + +#include + +#include "DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsOutDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersServiceDeviceCollection.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "HeterogeneousCore/AlpakaInterface/interface/traits.h" +#include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + class HGCalLayerClustersSoAAlgoWrapper { + public: + void run(Queue& queue, + const unsigned int numer_of_clusters, + float thresholdW0, + float positionDeltaRho2, + const HGCalSoACellsDeviceCollection::ConstView input_rechits_soa, + const HGCalSoACellsOutDeviceCollection::ConstView input_clusters_soa, + HGCalSoAClustersDeviceCollection::View outputs, + HGCalSoAClustersServiceDeviceCollection::View outputs_service) const; + }; +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#endif // RecoLocalCalo_HGCalRecProducers_plugins_alpaka_HGCalLayerClustersSoAAlgoWrapper_h diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoALayerClustersProducer.cc b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoALayerClustersProducer.cc new file mode 100644 index 0000000000000..195d8dbeef608 --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoALayerClustersProducer.cc @@ -0,0 +1,104 @@ +#include "DataFormats/PortableTestObjects/interface/alpaka/TestDeviceCollection.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/stream/SynchronizingEDProducer.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/EDPutToken.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESGetToken.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "HeterogeneousCore/AlpakaTest/interface/AlpakaESTestRecords.h" +#include "HeterogeneousCore/AlpakaTest/interface/alpaka/AlpakaESTestData.h" +#include "RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h" +#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h" +#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h" + +#include "DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoAClustersDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsOutDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/HGCalSoAClusters.h" +#include "DataFormats/HGCalReco/interface/HGCalSoAClustersService.h" +#include "HGCalLayerClustersSoAAlgoWrapper.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + class HGCalSoALayerClustersProducer : public stream::SynchronizingEDProducer<> { + public: + HGCalSoALayerClustersProducer(edm::ParameterSet const& config) + : getTokenDeviceRecHits_{consumes(config.getParameter("hgcalRecHitsSoA"))}, + getTokenDeviceClusters_{consumes(config.getParameter("hgcalRecHitsLayerClustersSoA"))}, + deviceTokenSoAClusters_{produces()}, + thresholdW0_(config.getParameter("thresholdW0")), + positionDeltaRho2_(config.getParameter("positionDeltaRho2")) {} + + ~HGCalSoALayerClustersProducer() override = default; + + void acquire(device::Event const& iEvent, device::EventSetup const& iSetup) override { + // Get LayerClusters almost-SoA on device: this has still the same + // cardinality as the RecHitsSoA, but has all the required information + // to assemble the clusters, i.e., it has the cluster index assigned to + // each rechit. + auto const& deviceInputClusters = iEvent.get(getTokenDeviceClusters_); + auto const inputClusters_v = deviceInputClusters.view(); + // + // Allocate output SoA for the clusters, one entry for each cluster + auto device_numclusters = cms::alpakatools::make_device_view( + alpaka::getDev(iEvent.queue()), inputClusters_v.numberOfClustersScalar()); + auto host_numclusters = cms::alpakatools::make_host_view(num_clusters_); + alpaka::memcpy(iEvent.queue(), host_numclusters, device_numclusters); + } + + void produce(device::Event& iEvent, device::EventSetup const& iSetup) override { + // Get RecHitsSoA on the device + auto const& deviceInputRecHits = iEvent.get(getTokenDeviceRecHits_); + auto const inputRechits_v = deviceInputRecHits.view(); + + // Get LayerClusters almost-SoA on device: this has still the same + // cardinality as the RecHitsSoA, but has all the required information + // to assemble the clusters, i.e., it has the cluster index assigned to + // each rechit. + auto const& deviceInputClusters = iEvent.get(getTokenDeviceClusters_); + auto const inputClusters_v = deviceInputClusters.view(); + + ALPAKA_ACCELERATOR_NAMESPACE::PortableCollection output(num_clusters_, iEvent.queue()); + auto output_v = output.view(); + // Allocate service SoA cluster + ALPAKA_ACCELERATOR_NAMESPACE::PortableCollection outputService(num_clusters_, + iEvent.queue()); + auto output_service_v = outputService.view(); + + algo_.run(iEvent.queue(), + num_clusters_, + thresholdW0_, + positionDeltaRho2_, + inputRechits_v, + inputClusters_v, + output_v, + output_service_v); + iEvent.emplace(deviceTokenSoAClusters_, std::move(output)); + } + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("hgcalRecHitsLayerClustersSoA", edm::InputTag("TO BE DEFINED")); + desc.add("hgcalRecHitsSoA", edm::InputTag("TO BE DEFINED")); + desc.add("thresholdW0", 2.9); + desc.add("positionDeltaRho2", 1.69); + descriptions.addWithDefaultLabel(desc); + } + + private: + device::EDGetToken> const getTokenDeviceRecHits_; + device::EDGetToken> const getTokenDeviceClusters_; + device::EDPutToken> const deviceTokenSoAClusters_; + HGCalLayerClustersSoAAlgoWrapper algo_; + unsigned int num_clusters_; + float thresholdW0_; + float positionDeltaRho2_; + }; + +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/MakerMacros.h" +DEFINE_FWK_ALPAKA_MODULE(HGCalSoALayerClustersProducer); diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsLayerClustersProducer.cc b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsLayerClustersProducer.cc new file mode 100644 index 0000000000000..e7082d1f264b9 --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsLayerClustersProducer.cc @@ -0,0 +1,69 @@ +#include "DataFormats/PortableTestObjects/interface/alpaka/TestDeviceCollection.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/stream/EDProducer.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/EDPutToken.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESGetToken.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "HeterogeneousCore/AlpakaTest/interface/AlpakaESTestRecords.h" +#include "HeterogeneousCore/AlpakaTest/interface/alpaka/AlpakaESTestData.h" +#include "RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h" +#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h" +#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h" + +#include "DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsDeviceCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsOutDeviceCollection.h" +#include "HGCalLayerClustersAlgoWrapper.h" + +// Processes the input RecHit SoA collection and generates an output SoA +// containing all the necessary information to build the clusters. +// Specifically, this producer does not create the clusters in any format. +// Instead, it fills a SoA (HGCalSoACellsOut) with the same size as the input +// RecHit SoA. This output SoA includes all the data needed to assemble the +// clusters and assigns a clusterId to each cell that belongs to a cluster. +// Consequently, this producer must be used by another downstream producer to +// either build traditional clusters or to create a SoA representing the +// clusters, complete with all required information (e.g., energy, position). +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + class HGCalSoARecHitsLayerClustersProducer : public stream::EDProducer<> { + public: + HGCalSoARecHitsLayerClustersProducer(edm::ParameterSet const& config) + : getTokenDevice_{consumes(config.getParameter("hgcalRecHitsSoA"))}, deviceToken_{produces()} {} + + ~HGCalSoARecHitsLayerClustersProducer() override = default; + + void produce(device::Event& iEvent, device::EventSetup const& iSetup) override { + auto const& deviceInput = iEvent.get(getTokenDevice_); + //std::cout << "Size of device collection: " << deviceInput->metadata().size() << std::endl; + auto const input_v = deviceInput.view(); + // Allocate output SoA + ALPAKA_ACCELERATOR_NAMESPACE::PortableCollection output(deviceInput->metadata().size(), + iEvent.queue()); + auto output_v = output.view(); + + algo_.run(iEvent.queue(), deviceInput->metadata().size(), input_v, output_v); + iEvent.emplace(deviceToken_, std::move(output)); + } + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("hgcalRecHitsSoA", edm::InputTag("TO BE DEFINED")); + descriptions.addWithDefaultLabel(desc); + } + + private: + // use device::EDGetToken to read from device memory space + device::EDGetToken> const getTokenDevice_; + device::EDPutToken> const deviceToken_; + HGCalLayerClustersAlgoWrapper algo_; + }; + +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/MakerMacros.h" +DEFINE_FWK_ALPAKA_MODULE(HGCalSoARecHitsLayerClustersProducer); diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsProducer.cc b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsProducer.cc new file mode 100644 index 0000000000000..6566e1f919d8b --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsProducer.cc @@ -0,0 +1,221 @@ +#include "DataFormats/PortableTestObjects/interface/alpaka/TestDeviceCollection.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/stream/EDProducer.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/EDPutToken.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESGetToken.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "HeterogeneousCore/AlpakaTest/interface/AlpakaESTestRecords.h" +#include "HeterogeneousCore/AlpakaTest/interface/alpaka/AlpakaESTestData.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" +#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h" +#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h" + +#include "DataFormats/HGCalReco/interface/HGCalSoACellsHostCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoACellsDeviceCollection.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + class HGCalSoARecHitsProducer : public stream::EDProducer<> { + public: + HGCalSoARecHitsProducer(edm::ParameterSet const& config) + : caloGeomToken_(consumesCollector().esConsumes()), + detector_(config.getParameter("detector")), + deviceToken_{produces()}, + initialized_(false) { + hits_token_ = consumes(config.getParameter("recHits")); + isNose_ = false; + if (detector_ == "HFNose") { + isNose_ = true; + } + maxNumberOfThickIndices_ = config.getParameter("maxNumberOfThickIndices"); + fcPerEle_ = config.getParameter("fcPerEle"); + fcPerMip_ = config.getParameter>("fcPerMip"); + nonAgedNoises_ = config.getParameter>("noises"); + ecut_ = config.getParameter("ecut"); + thicknessCorrection_ = config.getParameter>("thicknessCorrection"); + dEdXweights_ = config.getParameter>("dEdXweights"); + } + + ~HGCalSoARecHitsProducer() override = default; + + void produce(device::Event& iEvent, device::EventSetup const& iSetup) override { + edm::Handle hits_h; + + edm::ESHandle geom = iSetup.getHandle(caloGeomToken_); + rhtools_.setGeometry(*geom); + maxlayer_ = rhtools_.lastLayer(isNose_); + + hits_h = iEvent.getHandle(hits_token_); + auto const& hits = *(hits_h.product()); + computeThreshold(); + + // Count effective hits above threshold + uint32_t index = 0; + for (unsigned int i = 0; i < hits.size(); ++i) { + const HGCRecHit& hgrh = hits[i]; + DetId detid = hgrh.detid(); + unsigned int layerOnSide = (rhtools_.getLayerWithOffset(detid) - 1); + + // set sigmaNoise default value 1 to use kappa value directly in case of + // sensor-independent thresholds + int thickness_index = rhtools_.getSiThickIndex(detid); + if (thickness_index == -1) { + thickness_index = maxNumberOfThickIndices_; + } + double storedThreshold = thresholds_[layerOnSide][thickness_index]; + if (hgrh.energy() < storedThreshold) + continue; // this sets the ZS threshold at ecut times the sigma noise + index++; + } + + // Allocate Host SoA will contain one entry for each RecHit above threshold + HGCalSoACellsHostCollection cells(index, iEvent.queue()); // TO BE VERIFIED + auto cellsView = cells.view(); + + // loop over all hits and create the Hexel structure, skip energies below ecut + // for each layer and wafer calculate the thresholds (sigmaNoise and energy) + // once + index = 0; + for (unsigned int i = 0; i < hits.size(); ++i) { + const HGCRecHit& hgrh = hits[i]; + DetId detid = hgrh.detid(); + unsigned int layerOnSide = (rhtools_.getLayerWithOffset(detid) - 1); + + // set sigmaNoise default value 1 to use kappa value directly in case of + // sensor-independent thresholds + float sigmaNoise = 1.f; + int thickness_index = rhtools_.getSiThickIndex(detid); + if (thickness_index == -1) { + thickness_index = maxNumberOfThickIndices_; + } + double storedThreshold = thresholds_[layerOnSide][thickness_index]; + if (detid.det() == DetId::HGCalHSi || detid.subdetId() == HGCHEF) { + storedThreshold = thresholds_.at(layerOnSide).at(thickness_index + deltasi_index_regemfac_); + } + sigmaNoise = v_sigmaNoise_.at(layerOnSide).at(thickness_index); + + if (hgrh.energy() < storedThreshold) + continue; // this sets the ZS threshold at ecut times the sigma noise + // for the sensor + + const GlobalPoint position(rhtools_.getPosition(detid)); + int offset = ((rhtools_.zside(detid) + 1) >> 1) * maxlayer_; + int layer = layerOnSide + offset; + auto entryInSoA = cellsView[index]; + if (detector_ == "BH") { + entryInSoA.dim1() = position.eta(); + entryInSoA.dim2() = position.phi(); + } // else, isSilicon == true and eta phi values will not be used + else { + entryInSoA.dim1() = position.x(); + entryInSoA.dim2() = position.y(); + } + entryInSoA.dim3() = position.z(); + entryInSoA.weight() = hgrh.energy(); + entryInSoA.sigmaNoise() = sigmaNoise; + entryInSoA.layer() = layer; + entryInSoA.recHitIndex() = i; + entryInSoA.detid() = detid.rawId(); + entryInSoA.time() = hgrh.time(); + entryInSoA.timeError() = hgrh.timeError(); + index++; + } +#if 0 + std::cout << "Size: " << cells->metadata().size() << " count cells: " << index + << " i.e. " << cells->metadata().size() << std::endl; +#endif + + if constexpr (!std::is_same_v) { + // Trigger copy async to GPU + //std::cout << "GPU" << std::endl; + HGCalSoACellsDeviceCollection deviceProduct{cells->metadata().size(), iEvent.queue()}; // QUEUE TO BE VERIFIED + alpaka::memcpy(iEvent.queue(), deviceProduct.buffer(), cells.const_buffer()); + iEvent.emplace(deviceToken_, std::move(deviceProduct)); + } else { + //std::cout << "CPU" << std::endl; + iEvent.emplace(deviceToken_, std::move(cells)); + } + } + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("detector", "EE")->setComment("options EE, FH, BH, HFNose; other value defaults to EE"); + desc.add("recHits", edm::InputTag("HGCalRecHit", "HGCEERecHits")); + desc.add("maxNumberOfThickIndices", 6); + desc.add("fcPerEle", 0.00016020506); + desc.add>("fcPerMip"); + desc.add>("thicknessCorrection"); + desc.add>("noises"); + desc.add>("dEdXweights"); + desc.add("ecut", 3.); + descriptions.addWithDefaultLabel(desc); + } + + private: + hgcal::RecHitTools rhtools_; + edm::ESGetToken caloGeomToken_; + bool isNose_; + std::string detector_; + edm::EDGetTokenT hits_token_; + device::EDPutToken> const deviceToken_; + + bool initialized_; + + unsigned maxNumberOfThickIndices_; + std::vector> thresholds_; + std::vector> v_sigmaNoise_; + unsigned int maxlayer_; + double sciThicknessCorrection_; + std::vector fcPerMip_; + double fcPerEle_; + std::vector nonAgedNoises_; + double ecut_; + std::vector dEdXweights_; + std::vector thicknessCorrection_; + int deltasi_index_regemfac_; + + void computeThreshold() { + // To support the TDR geometry and also the post-TDR one (v9 onwards), we + // need to change the logic of the vectors containing signal to noise and + // thresholds. The first 3 indices will keep on addressing the different + // thicknesses of the Silicon detectors in CE_E , the next 3 indices will address + // the thicknesses of the Silicon detectors in CE_H, while the last one, number 6 (the + // seventh) will address the Scintillators. This change will support both + // geometries at the same time. + + if (initialized_) + return; // only need to calculate thresholds once + + initialized_ = true; + + std::vector dummy; + + dummy.resize(maxNumberOfThickIndices_ + !isNose_, 0); // +1 to accomodate for the Scintillators + thresholds_.resize(maxlayer_, dummy); + v_sigmaNoise_.resize(maxlayer_, dummy); + + for (unsigned ilayer = 1; ilayer <= maxlayer_; ++ilayer) { + for (unsigned ithick = 0; ithick < maxNumberOfThickIndices_; ++ithick) { + float sigmaNoise = 0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] / + (fcPerMip_[ithick] * thicknessCorrection_[ithick]); + thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_; + v_sigmaNoise_[ilayer - 1][ithick] = sigmaNoise; +#if 0 + std::cout << "ilayer: " << ilayer << " nonAgedNoises: " << nonAgedNoises_[ithick] + << " fcPerEle: " << fcPerEle_ << " fcPerMip: " << fcPerMip_[ithick] + << " noiseMip: " << fcPerEle_ * nonAgedNoises_[ithick] / fcPerMip_[ithick] + << " sigmaNoise: " << sigmaNoise << "\n"; +#endif + } + } + } + }; + +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/MakerMacros.h" +DEFINE_FWK_ALPAKA_MODULE(HGCalSoARecHitsProducer);