diff --git a/Dumpers/BinaryDumpers/options/dump_geometry.py b/Dumpers/BinaryDumpers/options/dump_geometry.py new file mode 100644 index 0000000000000000000000000000000000000000..005f42edb8d7f27d1f9f698e5d63802bdc7125ca --- /dev/null +++ b/Dumpers/BinaryDumpers/options/dump_geometry.py @@ -0,0 +1,103 @@ +############################################################################### +# (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### +from GaudiConf import IOHelper +from GaudiPython.Bindings import AppMgr, gbl +from Configurables import LHCbApp, CondDB +from Configurables import GaudiSequencer +from Configurables import FTRawBankDecoder +from Configurables import (PrStoreUTHit, PrStoreFTHit) +from Configurables import ApplicationMgr +from Configurables import HistogramPersistencySvc +from Configurables import (AuditorSvc, SequencerTimerTool) +from Configurables import (DumpUTGeometry, DumpFTGeometry, DumpMuonTable, + DumpUTLookupTables, DumpMuonGeometry, + DumpCaloGeometry, DumpVPGeometry, DumpMagneticField, + DumpBeamline) +from Configurables import RootHistCnv__PersSvc +from Configurables import IODataManager +from Configurables import createODIN +from Configurables import TestMuonTable + +app = LHCbApp( + DataType="Upgrade", + EvtMax=50, + Simulation=True, + DDDBtag="dddb-20171122", + CondDBtag="sim-20180530-vc-md100") + +# Upgrade DBs +CondDB().Upgrade = True + +dec_seq = GaudiSequencer("DecodingSeq") +dec_seq.Members = [createODIN()] + +ecal_location = "/dd/Structure/LHCb/DownstreamRegion/Ecal" +hcal_location = "/dd/Structure/LHCb/DownstreamRegion/Hcal" + +ecal_geom = DumpCaloGeometry("DumpEcal", Location=ecal_location) +hcal_geom = DumpCaloGeometry("DumpHcal", Location=hcal_location) + +# Add the service that will dump the UT and FT geometry +ApplicationMgr().ExtSvc += [ + DumpMagneticField(), + DumpBeamline(), + DumpVPGeometry(), + DumpUTGeometry(), + DumpFTGeometry(), + DumpMuonGeometry(), + DumpMuonTable(), + DumpUTLookupTables(), ecal_geom, hcal_geom +] + +ApplicationMgr().TopAlg = [dec_seq] + +# Some extra stuff for timing table +ApplicationMgr().ExtSvc += ['ToolSvc', 'AuditorSvc'] +ApplicationMgr().AuditAlgorithms = True +AuditorSvc().Auditors += ['TimingAuditor'] +SequencerTimerTool().OutputLevel = 4 + +# Some extra stuff to save histograms +ApplicationMgr().HistogramPersistency = "ROOT" +RootHistCnv__PersSvc('RootHistCnv').ForceAlphaIds = True +HistogramPersistencySvc().OutputFile = "dump-histos.root" + +# No error messages when reading MDF +IODataManager().DisablePFNWarning = True + +IOHelper('MDF').inputFiles( + [ + # SciFi v5, minbias + 'PFN:/data/bfys/raaij/upgrade_scifi_v5_uncompressed/upgrade_mc_minbias_scifi_v5_%03d.mdf' + % i for i in range(5) + ], + clear=True) + +# GaudiPython +gaudi = AppMgr() +gaudi.initialize() + +TES = gaudi.evtSvc() +det = gaudi.detSvc() +ecal = det[ecal_location] +hcal = det[hcal_location] + +ecal_type = gbl.LHCb.RawBank.EcalPacked +hcal_type = gbl.LHCb.RawBank.HcalPacked + +while True: + gaudi.run(1) + if not TES['/Event']: + break + raw_event = TES['DAQ/RawEvent'] + ecal_banks = raw_event.banks(ecal_type) + hcal_banks = raw_event.banks(hcal_type) + break diff --git a/Rec/Allen/CMakeLists.txt b/Rec/Allen/CMakeLists.txt index 4d651bd4080ced818dc2ac7565596301edc5f666..4246396a5b94975fb93970f3a4a93de9e0a0dc51 100644 --- a/Rec/Allen/CMakeLists.txt +++ b/Rec/Allen/CMakeLists.txt @@ -57,6 +57,7 @@ include_directories(${CMAKE_SOURCE_DIR}/device/vertex_fit/common/include) include_directories(${CMAKE_SOURCE_DIR}/device/vertex_fit/vertex_fitter/include) include_directories(${CMAKE_SOURCE_DIR}/device/selections/line_types/include) include_directories(${CMAKE_SOURCE_DIR}/device/selections/Hlt1/include) +include_directories(${CMAKE_SOURCE_DIR}/device/calo/clustering/include) include_directories(${CMAKE_SOURCE_DIR}/device/calo/event_model/include) include_directories(${CMAKE_SOURCE_DIR}/backend/include) include_directories(${CMAKE_SOURCE_DIR}/test/contracts/include) diff --git a/Rec/Allen/src/AllenCaloToCaloClusters.cpp b/Rec/Allen/src/AllenCaloToCaloClusters.cpp new file mode 100644 index 0000000000000000000000000000000000000000..aa82b7ab92bc45cdd271f4d30b6918d62452f8ee --- /dev/null +++ b/Rec/Allen/src/AllenCaloToCaloClusters.cpp @@ -0,0 +1,189 @@ +/***************************************************************************** \ + * (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration * + * * + * This software is distributed under the terms of the GNU General Public * + * Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". * + * * + * In applying this licence, CERN does not waive the privileges and immunities * + * granted to it by virtue of its status as an Intergovernmental Organization * + * or submit itself to any jurisdiction. * +\*****************************************************************************/ + +/** + * Convert AllenCalo to CaloCluster v2 + * + * author Dorothea vom Bruch + * + */ + +#include "AllenCaloToCaloClusters.h" + +DECLARE_COMPONENT(AllenCaloToCaloClusters) + +AllenCaloToCaloClusters::AllenCaloToCaloClusters(const std::string& name, ISvcLocator* pSvcLocator) : + MultiTransformer( + name, + pSvcLocator, + // Inputs + {KeyValue {"AllenOutput", "Allen/Out/HostBuffers"}}, + // Outputs + {KeyValue {"AllenEcalClusters", "Allen/Calo/EcalCluster"}, + KeyValue {"AllenHcalClusters", "Allen/Calo/HcalCluster"}}) +{} + +std::tuple<LHCb::Event::Calo::Clusters, LHCb::Event::Calo::Clusters> AllenCaloToCaloClusters::operator()( + const HostBuffers& host_buffers) const +{ + // avoid long names + // using namespace LHCb::CaloDataFunctor; + LHCb::Event::Calo::Clusters EcalClusters; + LHCb::Event::Calo::Clusters HcalClusters; + // Make the clusters + const unsigned i_event = 0; + const unsigned number_of_events = 1; + + unsigned number_of_ecal_clusters = + host_buffers.host_ecal_cluster_offsets[number_of_events] - host_buffers.host_ecal_cluster_offsets[i_event]; + unsigned number_of_hcal_clusters = + host_buffers.host_hcal_cluster_offsets[number_of_events] - host_buffers.host_hcal_cluster_offsets[i_event]; + CaloCluster* ecal_clusters = (CaloCluster*) (host_buffers.host_ecal_clusters.data()); + CaloCluster* hcal_clusters = (CaloCluster*) (host_buffers.host_hcal_clusters.data()); + + if (msgLevel(MSG::DEBUG)) { + debug() << "Number of Ecal clusters to convert = " << number_of_ecal_clusters << endmsg; + debug() << "Number of Hcal clusters to convert = " << number_of_hcal_clusters << endmsg; + } + + EcalClusters.reserveForEntries(number_of_ecal_clusters); + HcalClusters.reserveForEntries(number_of_hcal_clusters); + + // Loop over Allen Ecal clusters and convert them + // Don't need to access them with offset since one event is processed at a time + int16_t iFirstEntry = 0; + for (unsigned i = 0; i < number_of_ecal_clusters; i++) { + const auto& cluster = ecal_clusters[i]; + + auto seedCellID = LHCb::Calo::DenseIndex::details::toCellID(cluster.center_id); + + if (msgLevel(MSG::DEBUG) && i < 50 && i % 5 == 0) { + if (LHCb::Calo::isValid(seedCellID)) + debug() << "ECAL CellID " << seedCellID << " corresponding to dense ID " << cluster.center_id << " is invalid!" + << endmsg; + debug() << " \t ECAL center_id = " << cluster.center_id << " cellID: " << seedCellID << ", e = " << cluster.e + << ", x = " << cluster.x << ", y = " << cluster.y; + for (unsigned j = 0; j < Calo::Constants::max_neighbours; ++j) + debug() << " " << cluster.digits[j]; + debug() << endmsg; + } + + // Add the all digits, marking the seed ones + EcalClusters.emplace_back( + seedCellID, + cluster.e, + 1.0, + LHCb::CaloDigitStatus::Status {LHCb::CaloDigitStatus::Mask::UseForEnergy, LHCb::CaloDigitStatus::Mask::SeedCell}); + for (unsigned j = 0; j < Calo::Constants::max_neighbours; ++j) { + auto cellID = LHCb::Calo::DenseIndex::details::toCellID(cluster.digits[j]); + EcalClusters.emplace_back( + cellID, + 0., + 1.0, + LHCb::CaloDigitStatus::Status {LHCb::CaloDigitStatus::Mask::UseForEnergy, + LHCb::CaloDigitStatus::Mask::OwnedCell}); + } + EcalClusters.emplace_back( + seedCellID, + LHCb::Event::Calo::Clusters::Type::Area3x3, + {iFirstEntry, Calo::Constants::max_neighbours + 1}, + cluster.e, + {cluster.x, cluster.y, 0.}); + + iFirstEntry += Calo::Constants::max_neighbours + 1; // seed digit+ associated digits making the cluster + } + + // Loop over Allen Hcal clusters and convert them + + iFirstEntry = 0; + for (unsigned i = 0; i < number_of_hcal_clusters; i++) { + const auto& cluster = hcal_clusters[i]; + + auto seedCellID = LHCb::Calo::DenseIndex::details::toCellID(cluster.center_id); + + if (msgLevel(MSG::DEBUG) && i < 50 && i % 5 == 0) { + if (LHCb::Calo::isValid(seedCellID)) + debug() << "HCAL CellID " << seedCellID << " corresponding to dense ID " << cluster.center_id << " is invalid!" + << endmsg; + debug() << " \t HCAL center_id = " << cluster.center_id << " cellID: " << seedCellID << ", e = " << cluster.e + << ", x = " << cluster.x << ", y = " << cluster.y; + for (unsigned j = 0; j < Calo::Constants::max_neighbours; ++j) + debug() << " " << cluster.digits[j]; + debug() << endmsg; + } + + // Add the all digits, marking the seed ones + HcalClusters.emplace_back( + seedCellID, + cluster.e, + LHCb::CaloDigitStatus::Status {LHCb::CaloDigitStatus::Mask::UseForEnergy, LHCb::CaloDigitStatus::Mask::SeedCell}); + for (unsigned j = 0; j < Calo::Constants::max_neighbours; ++j) { + auto cellID = LHCb::Calo::DenseIndex::details::toCellID(cluster.digits[j]); + HcalClusters.emplace_back( + cellID, + .0, + LHCb::CaloDigitStatus::Status {LHCb::CaloDigitStatus::Mask::UseForEnergy, + LHCb::CaloDigitStatus::Mask::OwnedCell}); + } + LHCb::Event::Calo::v2::Clusters::which_entries_t entries {iFirstEntry, Calo::Constants::max_neighbours + 1}; + Gaudi::XYZPoint position {cluster.x, cluster.y, 0.}; + + HcalClusters.emplace_back(seedCellID, LHCb::Event::Calo::Clusters::Type::Area3x3, entries, cluster.e, position); + + iFirstEntry += Calo::Constants::max_neighbours + 1; // seed digit+ associated digits making the cluster + } + if (msgLevel(MSG::DEBUG)) { + debug() << "Number of ecal seed clusters: " << EcalClusters.size() << endmsg; + uint i = 0; + for (const auto& Cluster : EcalClusters) { + auto cellID = Cluster.cellID(); + const double e = Cluster.energy(); + const double x = Cluster.position().x(); + const double y = Cluster.position().y(); + const double z = Cluster.position().z(); + // const double et = e * sqrt( x * x + y * y ) / sqrt( x * x + y * y + z * z ); + + if (i % 5 == 0) { + debug() << "Ecal cellID: " << cellID << " energy = " << e << ", x = " << x << ", y = " << y << ", z = " << z + << endmsg; + auto digits = Cluster.entries(); + for (const auto& digit : digits) + debug() << " cellID: " << digit.cellID() << " energy: " << digit.energy() + << " fraction: " << digit.fraction() << " Status: " << digit.status() << endmsg; + } + if (i > 50) break; + ++i; + } + debug() << "Number of hcal seed clusters: " << HcalClusters.size() << endmsg; + i = 0; + for (const auto& Cluster : HcalClusters) { + auto cellID = Cluster.cellID(); + const double e = Cluster.energy(); + const double x = Cluster.position().x(); + const double y = Cluster.position().y(); + const double z = Cluster.position().z(); + // const double et = e * sqrt( x * x + y * y ) / sqrt( x * x + y * y + z * z ); + + if (i % 5 == 0) { + debug() << "Hcal cellID: " << cellID << " energy = " << e << ", x = " << x << ", y = " << y << ", z = " << z + << endmsg; + auto digits = Cluster.entries(); + for (const auto& digit : digits) + debug() << " cellID: " << digit.cellID() << " energy: " << digit.energy() + << " fraction: " << digit.fraction() << " Status: " << digit.status() << endmsg; + } + if (i > 50) break; + ++i; + } + } + + return {EcalClusters, HcalClusters}; +} diff --git a/Rec/Allen/src/AllenCaloToCaloClusters.h b/Rec/Allen/src/AllenCaloToCaloClusters.h new file mode 100644 index 0000000000000000000000000000000000000000..5c43b659b05f4759216fdcbf904dca9cc8915b93 --- /dev/null +++ b/Rec/Allen/src/AllenCaloToCaloClusters.h @@ -0,0 +1,47 @@ +/***************************************************************************** \ + * (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration * + * * + * This software is distributed under the terms of the GNU General Public * + * Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". * + * * + * In applying this licence, CERN does not waive the privileges and immunities * + * granted to it by virtue of its status as an Intergovernmental Organization * + * or submit itself to any jurisdiction. * +\*****************************************************************************/ +#ifndef ALLENCALOTOCALOCLUSTERS_H +#define ALLENCALOTOCALOCLUSTERS_H + +// Gaudi +#include "GaudiAlg/Transformer.h" + +// LHCb +#include "Event/Track.h" + +// Allen +#include "HostBuffers.cuh" +#include "Logger.h" +#include "VeloConsolidated.cuh" +#include "CaloCluster.cuh" +#include "CaloConstants.cuh" +//#include "Event/CaloCluster.h" +#include "Event/CaloHypos_v2.h" +#include "Event/CaloPosition.h" +#include "Event/CaloClusters_v2.h" +#include "Kernel/CaloCellID.h" +#include "GaudiKernel/Point3DTypes.h" + +class AllenCaloToCaloClusters final + : public Gaudi::Functional::MultiTransformer<std::tuple<LHCb::Event::Calo::Clusters, LHCb::Event::Calo::Clusters>( + const HostBuffers&)> { +public: + /// Standard constructor + AllenCaloToCaloClusters(const std::string& name, ISvcLocator* pSvcLocator); + + /// Algorithm execution + std::tuple<LHCb::Event::Calo::Clusters, LHCb::Event::Calo::Clusters> operator()(const HostBuffers&) const override; + +private: + Gaudi::Property<float> m_EtCalo {this, "EtCalo", 400 * Gaudi::Units::MeV, "Default ET for Calo Clusters"}; +}; + +#endif diff --git a/configuration/constants/Calo.json b/configuration/constants/Calo.json new file mode 100644 index 0000000000000000000000000000000000000000..7119bb7d6c1673ae472e89be9977217e8bb946b0 --- /dev/null +++ b/configuration/constants/Calo.json @@ -0,0 +1,3 @@ +{ + "configured_lines": [] +} diff --git a/configuration/sequences/calo.py b/configuration/sequences/calo.py new file mode 100644 index 0000000000000000000000000000000000000000..1db5eb8c1d986df120cd023b9f46b3e495952d82 --- /dev/null +++ b/configuration/sequences/calo.py @@ -0,0 +1,12 @@ +############################################################################### +# (c) Copyright 2018-2020 CERN for the benefit of the LHCb Collaboration # +############################################################################### +from definitions.GECSequence import GECSequence +from definitions.CaloSequence import CaloSequence +from definitions.algorithms import compose_sequences + +gec_sequence = GECSequence() + +calo_sequence = CaloSequence(gec_sequence['initialize_lists']) + +compose_sequences(gec_sequence, calo_sequence).generate() diff --git a/configuration/sequences/definitions/CaloSequence.py b/configuration/sequences/definitions/CaloSequence.py new file mode 100644 index 0000000000000000000000000000000000000000..0a913417627820e0e1c1d4809557b5583da8ced8 --- /dev/null +++ b/configuration/sequences/definitions/CaloSequence.py @@ -0,0 +1,61 @@ +############################################################################### +# (c) Copyright 2018-2020 CERN for the benefit of the LHCb Collaboration # +############################################################################### +from definitions.algorithms import * + + +def CaloSequence(initialize_lists): + ecal_banks = data_provider_t(name="ecal_banks", bank_type="EcalPacked") + + calo_count_digits = calo_count_digits_t( + name="calo_count_digits", + host_number_of_events_t=initialize_lists.host_number_of_events_t(), + dev_number_of_events_t=initialize_lists.dev_number_of_events_t(), + dev_event_list_t=initialize_lists.dev_event_list_t()) + + prefix_sum_ecal_num_digits = host_prefix_sum_t( + name="prefix_sum_ecal_num_digits", + dev_input_buffer_t=calo_count_digits.dev_ecal_num_digits_t()) + + calo_decode = calo_decode_t( + name="calo_decode", + host_number_of_events_t=initialize_lists.host_number_of_events_t(), + host_ecal_number_of_digits_t=prefix_sum_ecal_num_digits. + host_total_sum_holder_t(), + dev_event_list_t=initialize_lists.dev_event_list_t(), + dev_ecal_raw_input_t=ecal_banks.dev_raw_banks_t(), + dev_ecal_raw_input_offsets_t=ecal_banks.dev_raw_offsets_t(), + dev_ecal_digits_offsets_t=prefix_sum_ecal_num_digits. + dev_output_buffer_t()) + + calo_seed_clusters = calo_seed_clusters_t( + name="calo_seed_clusters", + host_number_of_events_t=initialize_lists.host_number_of_events_t(), + dev_event_list_t=initialize_lists.dev_event_list_t(), + dev_ecal_digits_t=calo_decode.dev_ecal_digits_t(), + dev_ecal_digits_offsets_t=prefix_sum_ecal_num_digits. + dev_output_buffer_t()) + + prefix_sum_ecal_num_clusters = host_prefix_sum_t( + name="prefix_sum_ecal_num_clusters", + dev_input_buffer_t=calo_seed_clusters.dev_ecal_num_clusters_t()) + + calo_find_clusters = calo_find_clusters_t( + name="calo_find_clusters", + host_ecal_number_of_clusters_t=prefix_sum_ecal_num_clusters. + host_total_sum_holder_t(), + dev_event_list_t=initialize_lists.dev_event_list_t(), + dev_ecal_digits_t=calo_decode.dev_ecal_digits_t(), + dev_ecal_digits_offsets_t=prefix_sum_ecal_num_digits. + dev_output_buffer_t(), + dev_ecal_seed_clusters_t=calo_seed_clusters.dev_ecal_seed_clusters_t(), + dev_ecal_cluster_offsets_t=prefix_sum_ecal_num_clusters. + dev_output_buffer_t()) + + calo_sequence = Sequence(ecal_banks, calo_count_digits, + prefix_sum_ecal_num_digits, + calo_decode, + calo_seed_clusters, prefix_sum_ecal_num_clusters, + calo_find_clusters) + + return calo_sequence diff --git a/device/calo/CMakeLists.txt b/device/calo/CMakeLists.txt index a77630131e6c1008d85bda0376254c321ae581a0..32a93472a3226ec35de607871d7b9c526a3c0f65 100644 --- a/device/calo/CMakeLists.txt +++ b/device/calo/CMakeLists.txt @@ -9,9 +9,11 @@ # or submit itself to any jurisdiction. # ############################################################################### file(GLOB calo_decoding "decoding/src/*cu") +file(GLOB calo_clustering "clustering/src/*cu") include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES}) include_directories(${CMAKE_SOURCE_DIR}/device/event_model/calo/include) +include_directories(clustering/include) include_directories(decoding/include) include_directories(${CMAKE_SOURCE_DIR}/main/include) include_directories(${CMAKE_SOURCE_DIR}/backend/include) @@ -23,4 +25,5 @@ include_directories(${CPPGSL_INCLUDE_DIR}) allen_add_device_library(Calo STATIC ${calo_decoding} + ${calo_clustering} ) diff --git a/device/calo/clustering/include/CaloFindClusters.cuh b/device/calo/clustering/include/CaloFindClusters.cuh new file mode 100644 index 0000000000000000000000000000000000000000..1d63662ede59702677f1e034c4af4603d776c95f --- /dev/null +++ b/device/calo/clustering/include/CaloFindClusters.cuh @@ -0,0 +1,49 @@ +/*****************************************************************************\ +* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration * +\*****************************************************************************/ + +#pragma once + +#include "CaloGeometry.cuh" +#include "CaloDigit.cuh" +#include "CaloCluster.cuh" +#include "DeviceAlgorithm.cuh" + +namespace calo_find_clusters { + struct Parameters { + HOST_INPUT(host_ecal_number_of_clusters_t, unsigned) host_ecal_number_of_clusters; + HOST_INPUT(host_hcal_number_of_clusters_t, unsigned) host_hcal_number_of_clusters; + DEVICE_INPUT(dev_event_list_t, unsigned) dev_event_list; + DEVICE_INPUT(dev_ecal_digits_t, CaloDigit) dev_ecal_digits; + DEVICE_INPUT(dev_hcal_digits_t, CaloDigit) dev_hcal_digits; + DEVICE_INPUT(dev_ecal_digits_offsets_t, unsigned) dev_ecal_digits_offsets; + DEVICE_INPUT(dev_hcal_digits_offsets_t, unsigned) dev_hcal_digits_offsets; + DEVICE_INPUT(dev_ecal_seed_clusters_t, CaloSeedCluster) dev_ecal_seed_clusters; + DEVICE_INPUT(dev_hcal_seed_clusters_t, CaloSeedCluster) dev_hcal_seed_clusters; + DEVICE_INPUT(dev_ecal_cluster_offsets_t, unsigned) dev_ecal_cluster_offsets; + DEVICE_INPUT(dev_hcal_cluster_offsets_t, unsigned) dev_hcal_cluster_offsets; + DEVICE_OUTPUT(dev_ecal_clusters_t, CaloCluster) dev_ecal_clusters; + DEVICE_OUTPUT(dev_hcal_clusters_t, CaloCluster) dev_hcal_clusters; + PROPERTY(block_dim_x_t, "block_dim_x", "block dimension X", unsigned) block_dim; + }; + + // Global function + __global__ void + calo_find_clusters(Parameters parameters, const char* raw_ecal_geometry, const char* raw_hcal_geometry); + + // Algorithm + struct calo_find_clusters_t : public DeviceAlgorithm, Parameters { + void set_arguments_size(ArgumentReferences<Parameters>, const RuntimeOptions&, const Constants&, const HostBuffers&) + const; + + __host__ void operator()( + const ArgumentReferences<Parameters>& arguments, + const RuntimeOptions& runtime_options, + const Constants& constants, + HostBuffers&, + Allen::Context const&) const; + + private: + Property<block_dim_x_t> m_block_dim_x {this, 32}; + }; +} // namespace calo_find_clusters diff --git a/device/calo/clustering/include/CaloSeedClusters.cuh b/device/calo/clustering/include/CaloSeedClusters.cuh new file mode 100644 index 0000000000000000000000000000000000000000..cc49af82d0a2a1dd522455ce104793abe73ebc01 --- /dev/null +++ b/device/calo/clustering/include/CaloSeedClusters.cuh @@ -0,0 +1,58 @@ +/*****************************************************************************\ +* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration * +\*****************************************************************************/ + +#pragma once + +#include "CaloGeometry.cuh" +#include "CaloDigit.cuh" +#include "CaloCluster.cuh" +#include "DeviceAlgorithm.cuh" + +namespace calo_seed_clusters { + struct Parameters { + HOST_INPUT(host_number_of_events_t, unsigned) host_number_of_events; + DEVICE_INPUT(dev_event_list_t, unsigned) dev_event_list; + DEVICE_INPUT(dev_ecal_digits_t, CaloDigit) dev_ecal_digits; + DEVICE_INPUT(dev_hcal_digits_t, CaloDigit) dev_hcal_digits; + DEVICE_INPUT(dev_ecal_digits_offsets_t, unsigned) dev_ecal_digits_offsets; + DEVICE_INPUT(dev_hcal_digits_offsets_t, unsigned) dev_hcal_digits_offsets; + DEVICE_OUTPUT(dev_ecal_num_clusters_t, unsigned) dev_ecal_num_clusters; + DEVICE_OUTPUT(dev_hcal_num_clusters_t, unsigned) dev_hcal_num_clusters; + DEVICE_OUTPUT(dev_ecal_seed_clusters_t, CaloSeedCluster) dev_ecal_seed_clusters; + DEVICE_OUTPUT(dev_hcal_seed_clusters_t, CaloSeedCluster) dev_hcal_seed_clusters; + PROPERTY(block_dim_x_t, "block_dim_x", "block dimension X", unsigned) block_dim; + PROPERTY(ecal_min_adc_t, "ecal_min_adc", "ECal seed cluster minimum ADC", int16_t) ecal_min_adc; + PROPERTY(hcal_min_adc_t, "hcal_min_adc", "HCal seed cluster minimum ADC", int16_t) hcal_min_adc; + }; + + // Global function + __global__ void calo_seed_clusters( + Parameters parameters, + const char* raw_ecal_geometry, + const char* raw_hcal_geometry, + const int16_t ecal_min_adc, + const int16_t hcal_min_adc); + + // Algorithm + struct calo_seed_clusters_t : public DeviceAlgorithm, Parameters { + + void set_arguments_size( + ArgumentReferences<Parameters> arguments, + const RuntimeOptions& runtime_options, + const Constants&, + const HostBuffers&) const; + + void operator()( + const ArgumentReferences<Parameters>& arguments, + const RuntimeOptions& runtime_options, + const Constants& constants, + HostBuffers&, + Allen::Context const&) const; + + private: + Property<block_dim_x_t> m_block_dim_x {this, 32}; + Property<ecal_min_adc_t> m_ecal_min_adc {this, 20}; + Property<hcal_min_adc_t> m_hcal_min_adc {this, 20}; + }; +} // namespace calo_seed_clusters diff --git a/device/calo/clustering/src/CaloFindClusters.cu b/device/calo/clustering/src/CaloFindClusters.cu new file mode 100644 index 0000000000000000000000000000000000000000..122772215ba23d7913336638a8ccd66438b3bb7c --- /dev/null +++ b/device/calo/clustering/src/CaloFindClusters.cu @@ -0,0 +1,120 @@ +/*****************************************************************************\ +* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration * +* * +* This software is distributed under the terms of the Apache License * +* version 2 (Apache-2.0), copied verbatim in the file "LICENSE". * +* * +* In applying this licence, CERN does not waive the privileges and immunities * +* granted to it by virtue of its status as an Intergovernmental Organization * +* or submit itself to any jurisdiction. * +\*****************************************************************************/ +#include <CaloCluster.cuh> +#include <CaloFindClusters.cuh> + +__device__ void simple_clusters( + CaloDigit const* digits, + CaloSeedCluster const* seed_clusters, + CaloCluster* clusters, + unsigned const num_clusters, + const CaloGeometry& calo) +{ + for (unsigned c = threadIdx.x; c < num_clusters; c += blockDim.x) { + auto const& seed_cluster = seed_clusters[c]; + auto& cluster = clusters[c]; + cluster.center_id = seed_cluster.id; + cluster.e = calo.gain[seed_cluster.id] * (seed_cluster.adc - calo.pedestal); + cluster.x = seed_cluster.x; + cluster.y = seed_cluster.y; + + uint16_t const* neighbors = &(calo.neighbors[seed_cluster.id * Calo::Constants::max_neighbours]); + for (uint16_t n = 0; n < Calo::Constants::max_neighbours; n++) { + auto const n_id = neighbors[n]; + int16_t adc = digits[n_id].adc; + if (n_id != 0 && (adc != SHRT_MAX)) { + cluster.e += (adc - calo.pedestal) * calo.gain[n_id]; + cluster.digits[n] = n_id; + } + else { + cluster.digits[n] = 0; + } + } + + // Code to update position, double check + // for (uint16_t n = 0; n < Calo::Constants::max_neighbours; n++) { + // auto const n_id = neighbors[n]; + // if (n_id != 0 && ((auto const adc = digits[n_id].adc) != SHRT_MAX)) { + // float const adc_frac = float(adc) / float(cluster.e); + // cluster.x += adc_frac * (geometry.getX(n_id) - seed_cluster.x); + // cluster.y += adc_frac * (geometry.getY(n_id) - seed_cluster.y); + // } + // } + } +} + +__global__ void calo_find_clusters::calo_find_clusters( + calo_find_clusters::Parameters parameters, + const char* raw_ecal_geometry, + const char* raw_hcal_geometry) +{ + // Get proper geometry. + auto ecal_geometry = CaloGeometry(raw_ecal_geometry); + auto hcal_geometry = CaloGeometry(raw_hcal_geometry); + + unsigned const event_number = parameters.dev_event_list[blockIdx.x]; + + // Build simple 3x3 clusters from seed clusters + // Ecal + unsigned const ecal_digits_offset = parameters.dev_ecal_digits_offsets[event_number]; + unsigned const ecal_clusters_offset = parameters.dev_ecal_cluster_offsets[event_number]; + unsigned const ecal_num_clusters = parameters.dev_ecal_cluster_offsets[event_number + 1] - ecal_clusters_offset; + simple_clusters( + parameters.dev_ecal_digits + ecal_digits_offset, + parameters.dev_ecal_seed_clusters + Calo::Constants::ecal_max_index / 8 * event_number, + parameters.dev_ecal_clusters + ecal_clusters_offset, + ecal_num_clusters, + ecal_geometry); + + // Hcal + unsigned const hcal_digits_offset = parameters.dev_hcal_digits_offsets[event_number]; + unsigned const hcal_clusters_offset = parameters.dev_hcal_cluster_offsets[event_number]; + unsigned const hcal_num_clusters = parameters.dev_hcal_cluster_offsets[event_number + 1] - hcal_clusters_offset; + simple_clusters( + parameters.dev_hcal_digits + hcal_digits_offset, + parameters.dev_hcal_seed_clusters + Calo::Constants::hcal_max_index / 8 * event_number, + parameters.dev_hcal_clusters + hcal_clusters_offset, + hcal_num_clusters, + hcal_geometry); +} + +void calo_find_clusters::calo_find_clusters_t::set_arguments_size( + ArgumentReferences<Parameters> arguments, + const RuntimeOptions&, + const Constants&, + const HostBuffers&) const +{ + set_size<dev_ecal_clusters_t>(arguments, first<host_ecal_number_of_clusters_t>(arguments)); + set_size<dev_hcal_clusters_t>(arguments, first<host_hcal_number_of_clusters_t>(arguments)); +} + +__host__ void calo_find_clusters::calo_find_clusters_t::operator()( + const ArgumentReferences<Parameters>& arguments, + const RuntimeOptions& runtime_options, + const Constants& constants, + HostBuffers& host_buffers, + Allen::Context const& context) const +{ + // Find clusters. + initialize<dev_ecal_clusters_t>(arguments, SHRT_MAX, context); + initialize<dev_hcal_clusters_t>(arguments, SHRT_MAX, context); + + global_function(calo_find_clusters)( + dim3(size<dev_event_list_t>(arguments)), dim3(property<block_dim_x_t>().get()), context)( + arguments, constants.dev_ecal_geometry, constants.dev_hcal_geometry); + + if (runtime_options.do_check) { + safe_assign_to_host_buffer<dev_ecal_cluster_offsets_t>(host_buffers.host_ecal_cluster_offsets, arguments, context); + safe_assign_to_host_buffer<dev_hcal_cluster_offsets_t>(host_buffers.host_hcal_cluster_offsets, arguments, context); + safe_assign_to_host_buffer<dev_ecal_clusters_t>(host_buffers.host_ecal_clusters, arguments, context); + safe_assign_to_host_buffer<dev_hcal_clusters_t>(host_buffers.host_hcal_clusters, arguments, context); + } +} diff --git a/device/calo/clustering/src/CaloSeedClusters.cu b/device/calo/clustering/src/CaloSeedClusters.cu new file mode 100644 index 0000000000000000000000000000000000000000..d145a3df8f6c6124b51a2c023c2691b46c711713 --- /dev/null +++ b/device/calo/clustering/src/CaloSeedClusters.cu @@ -0,0 +1,107 @@ +/*****************************************************************************\ +* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration * +* * +* This software is distributed under the terms of the Apache License * +* version 2 (Apache-2.0), copied verbatim in the file "LICENSE". * +* * +* In applying this licence, CERN does not waive the privileges and immunities * +* granted to it by virtue of its status as an Intergovernmental Organization * +* or submit itself to any jurisdiction. * +\*****************************************************************************/ +#include <CaloSeedClusters.cuh> + +__device__ void seed_clusters( + CaloDigit const* digits, + unsigned const num_digits, + CaloSeedCluster* clusters, + unsigned* num_clusters, + const CaloGeometry& geometry, + const int16_t min_adc) +{ + // Loop over all CellIDs. + for (unsigned i = threadIdx.x; i < num_digits; i += blockDim.x) { + int16_t adc = digits[i].adc; + if (adc == SHRT_MAX || adc < min_adc) { + continue; + } + uint16_t* neighbors = &(geometry.neighbors[i * Calo::Constants::max_neighbours]); + bool is_max = true; + for (unsigned n = 0; n < Calo::Constants::max_neighbours; n++) { + auto const neighbor_adc = digits[neighbors[n]].adc; + is_max = is_max && (neighbor_adc == SHRT_MAX || adc > neighbor_adc); + } + if (is_max) { + auto const id = atomicAdd(num_clusters, 1); + clusters[id] = CaloSeedCluster(i, digits[i].adc, geometry.getX(i), geometry.getY(i)); + } + } +} + +__global__ void calo_seed_clusters::calo_seed_clusters( + calo_seed_clusters::Parameters parameters, + const char* raw_ecal_geometry, + const char* raw_hcal_geometry, + const int16_t ecal_min_adc, + const int16_t hcal_min_adc) +{ + unsigned const event_number = parameters.dev_event_list[blockIdx.x]; + + // Get geometry. + auto ecal_geometry = CaloGeometry(raw_ecal_geometry); + auto hcal_geometry = CaloGeometry(raw_hcal_geometry); + + // ECal + auto const ecal_digits_offset = parameters.dev_ecal_digits_offsets[event_number]; + seed_clusters( + ¶meters.dev_ecal_digits[ecal_digits_offset], + parameters.dev_ecal_digits_offsets[event_number + 1] - ecal_digits_offset, + ¶meters.dev_ecal_seed_clusters[Calo::Constants::ecal_max_index / 8 * event_number], + ¶meters.dev_ecal_num_clusters[event_number], + ecal_geometry, + ecal_min_adc); + + // HCal + auto const hcal_digits_offset = parameters.dev_hcal_digits_offsets[event_number]; + seed_clusters( + ¶meters.dev_hcal_digits[hcal_digits_offset], + parameters.dev_hcal_digits_offsets[event_number + 1] - hcal_digits_offset, + ¶meters.dev_hcal_seed_clusters[Calo::Constants::hcal_max_index / 8 * event_number], + ¶meters.dev_hcal_num_clusters[event_number], + hcal_geometry, + hcal_min_adc); +} + +void calo_seed_clusters::calo_seed_clusters_t::set_arguments_size( + ArgumentReferences<Parameters> arguments, + const RuntimeOptions&, + const Constants&, + const HostBuffers&) const +{ + auto const n_events = first<host_number_of_events_t>(arguments); + set_size<dev_ecal_num_clusters_t>(arguments, n_events); + set_size<dev_hcal_num_clusters_t>(arguments, n_events); + + // TODO: get this from the geometry too + set_size<dev_ecal_seed_clusters_t>(arguments, Calo::Constants::ecal_max_index / 8 * n_events); + set_size<dev_hcal_seed_clusters_t>(arguments, Calo::Constants::hcal_max_index / 8 * n_events); +} + +void calo_seed_clusters::calo_seed_clusters_t::operator()( + const ArgumentReferences<Parameters>& arguments, + const RuntimeOptions&, + const Constants& constants, + HostBuffers&, + Allen::Context const& context) const +{ + initialize<dev_ecal_num_clusters_t>(arguments, 0, context); + initialize<dev_hcal_num_clusters_t>(arguments, 0, context); + + // Find local maxima. + global_function(calo_seed_clusters)( + dim3(size<dev_event_list_t>(arguments)), dim3(property<block_dim_x_t>().get()), context)( + arguments, + constants.dev_ecal_geometry, + constants.dev_hcal_geometry, + property<ecal_min_adc_t>().get(), + property<hcal_min_adc_t>().get()); +} diff --git a/device/event_model/calo/include/CaloCluster.cuh b/device/event_model/calo/include/CaloCluster.cuh new file mode 100644 index 0000000000000000000000000000000000000000..d0bc0c09c644a547da0ab18eae7683b599621464 --- /dev/null +++ b/device/event_model/calo/include/CaloCluster.cuh @@ -0,0 +1,49 @@ +/*****************************************************************************\ +* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration * +* * +* This software is distributed under the terms of the Apache License * +* version 2 (Apache-2.0), copied verbatim in the file "LICENSE". * +* * +* In applying this licence, CERN does not waive the privileges and immunities * +* granted to it by virtue of its status as an Intergovernmental Organization * +* or submit itself to any jurisdiction. * +\*****************************************************************************/ +#pragma once + +#include <BackendCommon.h> +#include "CaloConstants.cuh" + +struct CaloDigitClusters { + uint16_t clustered_at_iteration = 0; + uint16_t clusters[Calo::Constants::digit_max_clusters]; + + __device__ __host__ CaloDigitClusters() + { + for (int i = 0; i < Calo::Constants::digit_max_clusters; i++) { + clusters[i] = 0; + } + } +}; + +struct CaloCluster { + float e = 0.f; + float x = 0.f; + float y = 0.f; + uint16_t center_id = 0; + uint16_t digits[Calo::Constants::max_neighbours] = {0, 0, 0, 0, 0, 0, 0, 0, 0}; + + __device__ __host__ CaloCluster(uint16_t index, float energy, float rX, float rY) : + e {energy}, x {rX}, y {rY}, center_id {index} + {} +}; + +struct CaloSeedCluster { + uint16_t id = 0; + int16_t adc = 0; // Is a double in the original algorithm, but is still integer here? + float x = 0.f; + float y = 0.f; + + __device__ __host__ CaloSeedCluster(uint16_t cellid, int16_t a, float rX, float rY) : + id {cellid}, adc {a}, x {rX}, y {rY} + {} +}; diff --git a/device/event_model/calo/include/CaloConstants.cuh b/device/event_model/calo/include/CaloConstants.cuh index f43e134c7654abb762721cbc96e0aca444de696a..b16d2bbcbcaba7c6436e3277f8f630c950eb6a5c 100644 --- a/device/event_model/calo/include/CaloConstants.cuh +++ b/device/event_model/calo/include/CaloConstants.cuh @@ -13,5 +13,16 @@ namespace Calo { namespace Constants { constexpr uint16_t card_channels = 32; + constexpr uint16_t ecal_max_index = 6016; + constexpr uint16_t hcal_max_index = 1488; + + // Max distance based on CellIDs is 64 steps away, so the iteration in which a cell is clustered can never be more + // than 64. + constexpr uint16_t unclustered = 65; + constexpr uint16_t ignore = unclustered + 1; + + constexpr uint16_t digit_max_clusters = 15; + constexpr uint16_t max_neighbours = 9; + } // namespace Constants } // namespace Calo diff --git a/main/src/Allen.cpp b/main/src/Allen.cpp index c7d589db99978b5c623913cea5efcda05329cde2..ad5290ce01b6b56033ce889a9b090622737aecd0 100644 --- a/main/src/Allen.cpp +++ b/main/src/Allen.cpp @@ -388,8 +388,13 @@ int allen( else { mep_layout = false; // The binary input provider expects the folders for the bank types as connections - std::vector<std::string> connections = { - folder_name_velopix_raw, folder_name_UT_raw, folder_name_SciFi_raw, folder_name_Muon_raw, folder_name_ODIN_raw}; + std::vector<std::string> connections = {folder_name_velopix_raw, + folder_name_UT_raw, + folder_name_SciFi_raw, + folder_name_Muon_raw, + folder_name_ecal_raw, + folder_name_hcal_raw, + folder_name_ODIN_raw}; input_provider = std::make_unique<BinaryProvider<BankTypes::VP, BankTypes::UT, BankTypes::FT, BankTypes::MUON, BankTypes::ODIN>>( number_of_slices, events_per_slice, n_events, std::move(connections), n_io_reps, file_list); diff --git a/stream/sequence/include/HostBuffers.cuh b/stream/sequence/include/HostBuffers.cuh index 7cb456619e49787111e157dd4361989f274c3d50..7d06b34757002383fa870857881bfae71b0c0f9d 100644 --- a/stream/sequence/include/HostBuffers.cuh +++ b/stream/sequence/include/HostBuffers.cuh @@ -10,6 +10,7 @@ #include <gsl/gsl> #include <CaloDigit.cuh> +#include <CaloCluster.cuh> #include <BackendCommon.h> @@ -109,6 +110,10 @@ struct HostBuffers { bool* host_match_upstream_muon; // Calo + gsl::span<unsigned> host_ecal_cluster_offsets = {}; + gsl::span<unsigned> host_hcal_cluster_offsets = {}; + gsl::span<CaloCluster> host_ecal_clusters = {}; + gsl::span<CaloCluster> host_hcal_clusters = {}; gsl::span<unsigned> host_ecal_digits_offsets = {}; gsl::span<unsigned> host_hcal_digits_offsets = {}; gsl::span<CaloDigit> host_ecal_digits = {};