diff --git a/Dumpers/BinaryDumpers/CMakeLists.txt b/Dumpers/BinaryDumpers/CMakeLists.txt index f9cfc428b3475c9d614f246e542e7de9892795e2..12f161691d033c0c5290d5d64f9f1e8ec536697e 100644 --- a/Dumpers/BinaryDumpers/CMakeLists.txt +++ b/Dumpers/BinaryDumpers/CMakeLists.txt @@ -34,7 +34,8 @@ gaudi_add_library(BinaryDumpers LHCb::RecEvent LHCb::RichFutureDAQLib LHCb::RichFutureKernel - ) + LHCb::RichDetectorsLib + LHCb::RichDetLib) gaudi_add_module(BinaryDumpersModule SOURCES @@ -59,6 +60,7 @@ gaudi_add_module(BinaryDumpersModule src/DumpVPGeometry.cpp src/DumpRichCableMapping.cpp src/DumpRichPDMDBMapping.cpp + src/DumpRichGeometry.cpp src/PVDumper.cpp src/ProvideConstants.cpp src/TestMuonTable.cpp @@ -79,7 +81,8 @@ gaudi_add_module(BinaryDumpersModule LHCb::CaloDetLib LHCb::FTDetLib LHCb::MuonDAQLib - Rec::PrKernel + LHCb::RichDetLib + Rec::PrKernel yaml-cpp::yaml-cpp BinaryDumpers) diff --git a/Dumpers/BinaryDumpers/include/Dumpers/Identifiers.h b/Dumpers/BinaryDumpers/include/Dumpers/Identifiers.h index 352bca5456dd06a72622eff1201a33013eeaa46b..767a1a82a2e6717035dd9932dc7a314bb92fa3bc 100644 --- a/Dumpers/BinaryDumpers/include/Dumpers/Identifiers.h +++ b/Dumpers/BinaryDumpers/include/Dumpers/Identifiers.h @@ -122,5 +122,19 @@ namespace Allen { inline static std::string const id = "RichCableMapping"; }; + /** @class Rich1 + * Identifier for the RICH 1 objects for Allen + */ + struct Rich1Geometry : Identifier { + inline static std::string const id = "Rich1Geometry"; + }; + + /** @class Rich2 + * Identifier for the RICH 2 objects for Allen + */ + struct Rich2Geometry : Identifier { + inline static std::string const id = "Rich2Geometry"; + }; + } // namespace NonEventData } // namespace Allen diff --git a/Dumpers/BinaryDumpers/src/DumpRichCableMapping.cpp b/Dumpers/BinaryDumpers/src/DumpRichCableMapping.cpp index 9bd7c74627bb143d7ac3c8a8b634b71acc982069..f952d2f8bb7219b86147b7f8ca6227c667e763e5 100644 --- a/Dumpers/BinaryDumpers/src/DumpRichCableMapping.cpp +++ b/Dumpers/BinaryDumpers/src/DumpRichCableMapping.cpp @@ -33,13 +33,61 @@ // Dumper #include "Dumper.h" #include +#include "RichTel40CableMapping.cuh" namespace { struct RichCableMapping { RichCableMapping() = default; RichCableMapping(std::vector& data, const Rich::Future::DAQ::Tel40CableMapping& tel40Maps) { - Rich::Future::DAQ::Allen::Tel40CableMapping allenTel40Maps {tel40Maps}; + Allen::Rich::Decoding::Tel40CableMapping allenTel40Maps; + + allenTel40Maps.m_isInitialised = tel40Maps.isInitialised(); + allenTel40Maps.m_mappingVer = tel40Maps.version(); + + if (tel40Maps.isInitialised()) { + for (unsigned int i = 0; i < tel40Maps.tel40ModuleData().size(); ++i) { + for (unsigned int j = 0; j < tel40Maps.tel40ModuleData()[i].size(); ++j) { + for (unsigned int k = 0; k < tel40Maps.tel40ModuleData()[i][j].size(); ++k) { + auto& allenData = allenTel40Maps.m_tel40ModuleData[i][j][k]; + const auto& data = tel40Maps.tel40ModuleData()[i][j][k]; + + allenData.smartID = LHCb::RichSmartID::KeyType {data.smartID}; + allenData.moduleNum = data.moduleNum.data(); + allenData.sourceID = data.sourceID.data(); + allenData.connector = data.connector.data(); + allenData.pdmdbNum = data.pdmdbNum.data(); + allenData.linkNum = data.linkNum.data(); + allenData.isHType = data.isHType; + allenData.isActive = data.isActive; + } + } + } + + for (unsigned int i = 0; i < tel40Maps.tel40ConnData().size(); ++i) { + for (unsigned int j = 0; j < tel40Maps.tel40ConnData()[i].size(); ++j) { + for (unsigned int k = 0; k < tel40Maps.tel40ConnData()[i][j].size(); ++k) { + const auto& link_data = tel40Maps.tel40ConnData()[i][j][k]; + allenTel40Maps.m_tel40ConnMeta[i][j][k] = Allen::Rich::Decoding::Tel40CableMapping::Tel40MetaData { + static_cast(link_data.nActiveLinks), link_data.hasInactiveLinks}; + for (unsigned int l = 0; l < link_data.size(); ++l) { + auto& allenData = allenTel40Maps.m_tel40ConnData[i][j][k][l]; + const auto& data = link_data[l]; + + allenData.smartID = LHCb::RichSmartID::KeyType {data.smartID}; + allenData.moduleNum = data.moduleNum.data(); + allenData.sourceID = data.sourceID.data(); + allenData.connector = data.connector.data(); + allenData.pdmdbNum = data.pdmdbNum.data(); + allenData.linkNum = data.linkNum.data(); + allenData.isHType = data.isHType; + allenData.isActive = data.isActive; + } + } + } + } + } + DumpUtils::Writer output {}; output.write(allenTel40Maps); data = output.buffer(); diff --git a/Dumpers/BinaryDumpers/src/DumpRichGeometry.cpp b/Dumpers/BinaryDumpers/src/DumpRichGeometry.cpp new file mode 100644 index 0000000000000000000000000000000000000000..5cf5626a433539c8621424362db832ed49c716bc --- /dev/null +++ b/Dumpers/BinaryDumpers/src/DumpRichGeometry.cpp @@ -0,0 +1,144 @@ +/*****************************************************************************\ + * (c) Copyright 2000-2019 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. * + \*****************************************************************************/ + +// Gaudi Array properties ( must be first ...) +#include "GaudiKernel/StdArrayAsProperty.h" + +// Rich Kernel +#include "RichFutureKernel/RichAlgBase.h" +#include "Kernel/RichDetectorType.h" + +// Gaudi Functional +#include "LHCbAlgs/Transformer.h" + +// Rich Utils +#include "RichFutureUtils/RichDecodedData.h" +#include "RichUtils/RichException.h" +#include "RichUtils/RichHashMap.h" +#include "RichUtils/RichMap.h" +#include "RichUtils/RichSmartIDSorter.h" + +// LHCb +#include +#include + +// Dumper +#include "Dumper.h" +#include +#include "Rich.cuh" + +namespace { + template + struct RichGeometry_t { + RichGeometry_t() = default; + RichGeometry_t(std::vector& data, const RichT& rich) + { + Allen::Rich::RichDetector allenRich; + allenRich.m_type = RichID; + for (size_t k = 0; k < rich.pdPanels().size(); k++) { + auto& panel = rich.pdPanels()[k]; + auto& allenPanel = allenRich.m_panels[k]; + + allenPanel.m_modNumOffset = panel.modNumOffset(); + allenPanel.m_panelID = {panel.rich(), panel.side(), 0}; + allenPanel.m_rich = (Allen::Rich::Detector::Type) panel.rich(); + allenPanel.m_side = (Allen::Rich::Detector::Side) panel.side(); + panel.globalToPDPanel().GetComponents(allenPanel.m_gloToPDPanelM.begin()); + for (size_t i = 0; i < panel.pdModules().size(); i++) { + for (size_t j = 0; j < panel.pdModules()[i].size(); j++) { + if (panel.pdModules()[i][j] != nullptr) { + auto& allenPD = allenPanel.m_PDs[i][j]; + auto& pd = panel.pdModules()[i][j]; + + allenPD.m_pdSmartID = pd->pdSmartID(); + allenPD.m_effPixelArea = pd->effectivePixelArea(); + allenPD.m_numPixels = pd->effectiveNumActivePixels(); + allenPD.m_isHType = pd->isHType(); + allenPD.m_localZcoord = Rich::Detector::scalar(pd->localZCoord()); + allenPD.m_numPixColFrac = Rich::Detector::scalar(pd->getNumPixColFrac()); + allenPD.m_numPixRowFrac = Rich::Detector::scalar(pd->getNumPixRowFrac()); + allenPD.m_effectivePixelXSize = Rich::Detector::scalar(pd->getEffectivePixelXSize()); + allenPD.m_effectivePixelYSize = Rich::Detector::scalar(pd->getEffectivePixelYSize()); + pd->localToGlobal().GetComponents(allenPD.m_locToGloM.begin(), allenPD.m_locToGloM.end()); + pd->centrePointPanel().GetCoordinates(allenPD.m_zeroInPanelFrame.begin()); + } + else { + // Create null PD where HLT2 has nullptr + allenPanel.m_PDs[i][j].setIsNull(true); + } + } + } + } + + DumpUtils::Writer output {}; + output.write(allenRich); + data = output.buffer(); + } + }; + + using Rich1Geometry_t = RichGeometry_t; + using Rich2Geometry_t = RichGeometry_t; +} // namespace + +/** + * @brief Dump RICH detector object information. + */ +class DumpRichGeometry final : public Allen::Dumpers::Dumper< + void(Rich1Geometry_t const&, Rich2Geometry_t const&), + LHCb::DetDesc::usesConditions> { +public: + DumpRichGeometry(const std::string& name, ISvcLocator* svcLoc); + void operator()(const Rich1Geometry_t& Rich1Geo, const Rich2Geometry_t& Rich2Geo) const override; + StatusCode initialize() override; + +private: + std::vector m_Rich1Data; + std::vector m_Rich2Data; +}; + +DECLARE_COMPONENT(DumpRichGeometry) + +DumpRichGeometry::DumpRichGeometry(const std::string& name, ISvcLocator* svcLoc) : + Dumper( + name, + svcLoc, + {KeyValue {"Rich1GeometryLocation", location(name, "rich1geometry")}, + KeyValue {"Rich2GeometryLocation", location(name, "rich2Geometry")}}) +{} + +StatusCode DumpRichGeometry::initialize() +{ + return Dumper::initialize().andThen([&, this] { + register_producer(Allen::NonEventData::Rich1Geometry::id, "rich_1_geometry", m_Rich1Data); + Rich::Detector::Rich1::addConditionDerivation(this); + addConditionDerivation( + {Rich::Detector::Rich1::DefaultConditionKey}, + inputLocation(), + [&](const Rich::Detector::Rich1& det) { + Rich1Geometry_t Rich1Geo {m_Rich1Data, det}; + dump(); + return Rich1Geo; + }); + + register_producer(Allen::NonEventData::Rich2Geometry::id, "rich_2_geometry", m_Rich2Data); + Rich::Detector::Rich2::addConditionDerivation(this); + addConditionDerivation( + {Rich::Detector::Rich2::DefaultConditionKey}, + inputLocation(), + [&](const Rich::Detector::Rich2& det) { + Rich2Geometry_t Rich2Geo {m_Rich2Data, det}; + dump(); + return Rich2Geo; + }); + }); +} + +void DumpRichGeometry::operator()(const Rich1Geometry_t&, const Rich2Geometry_t&) const {} diff --git a/Dumpers/BinaryDumpers/src/DumpRichPDMDBMapping.cpp b/Dumpers/BinaryDumpers/src/DumpRichPDMDBMapping.cpp index af2dc548377c1cca5189c65c868fbf3bbbb08743..463c36d4d7f4f6ce69d6cf11a729123516020118 100644 --- a/Dumpers/BinaryDumpers/src/DumpRichPDMDBMapping.cpp +++ b/Dumpers/BinaryDumpers/src/DumpRichPDMDBMapping.cpp @@ -33,15 +33,24 @@ // Dumper #include "Dumper.h" #include +#include "RichPDMDBDecodeMapping.cuh" namespace { struct RichPDMDBMapping { RichPDMDBMapping() = default; - RichPDMDBMapping(std::vector& data, const Rich::Future::DAQ::PDMDBDecodeMapping& tel40Maps) + RichPDMDBMapping(std::vector& data, const Rich::Future::DAQ::PDMDBDecodeMapping& mapping) { - Rich::Future::DAQ::Allen::PDMDBDecodeMapping allenTel40Maps {tel40Maps}; + Allen::Rich::Decoding::PDMDBDecodeMapping allenPDMDBMapping; + + allenPDMDBMapping.m_pdmDataR = + std::bit_cast(mapping.pdmDataR()); + allenPDMDBMapping.m_pdmDataH = + std::bit_cast(mapping.pdmDataH()); + allenPDMDBMapping.m_isInitialised = mapping.isInitialised(); + allenPDMDBMapping.m_mappingVer = mapping.version(); + DumpUtils::Writer output {}; - output.write(allenTel40Maps); + output.write(allenPDMDBMapping); data = output.buffer(); } }; diff --git a/Rec/Allen/CMakeLists.txt b/Rec/Allen/CMakeLists.txt index 53c981529a7d1816e5deb8ff8d4883a6e5a46417..1abe385c8be3189ff492a5738fc76c9f46252cc8 100644 --- a/Rec/Allen/CMakeLists.txt +++ b/Rec/Allen/CMakeLists.txt @@ -39,7 +39,8 @@ gaudi_add_module(AllenWrapper src/CompareRecAllenFTClusters.cpp src/CompareRecAllenVPHits.cpp src/TestVeloClusters.cpp - src/TestAllenRichPixels.cpp + src/CompareRecAllenRichSmartIDs.cpp + src/CompareRecAllenRichPixels.cpp src/IOMonitor.cpp LINK AllenLib diff --git a/Rec/Allen/python/Allen/config.py b/Rec/Allen/python/Allen/config.py index 572b7ea16ea34f169d32b6c5e3e0b872efc5b427..e09a9f92fde53f1ddd674ef10ecb62e397d63eb5 100755 --- a/Rec/Allen/python/Allen/config.py +++ b/Rec/Allen/python/Allen/config.py @@ -21,7 +21,7 @@ from PyConf.Algorithms import ( DumpMagneticFieldPolarity, DumpVPGeometry, DumpCrossingAngles, DumpFTGeometry, DumpUTGeometry, DumpUTLookupTables, DumpMuonGeometry, DumpMuonTable, AllenODINProducer, DumpRichPDMDBMapping, - DumpRichCableMapping) + DumpRichCableMapping, DumpRichGeometry) from DDDB.CheckDD4Hep import UseDD4Hep from PyConf.reading import get_generator_BeamParameters from GaudiConf.LbExec import Options as DefaultOptions @@ -122,7 +122,9 @@ def setup_allen_non_event_data_service(allen_event_loop=False, 'Rich': [(DumpRichPDMDBMapping, 'DeviceRichPDMDBMapping', {}, 'rich_pdmdbmaps'), (DumpRichCableMapping, 'DeviceRichCableMapping', {}, - 'rich_tel40maps')] + 'rich_tel40maps'), + (DumpRichGeometry, 'DeviceRichGeometry', {}, + 'rich_geometry')] } else: converter_types = { @@ -147,7 +149,9 @@ def setup_allen_non_event_data_service(allen_event_loop=False, 'Rich': [(DumpRichPDMDBMapping, 'DeviceRichPDMDBMapping', {}, 'rich_pdmdbmaps'), (DumpRichCableMapping, 'DeviceRichCableMapping', {}, - 'rich_tel40maps')] + 'rich_tel40maps'), + (DumpRichGeometry, 'DeviceRichGeometry', {}, + 'rich_geometry')] } detector_names = { diff --git a/Rec/Allen/src/CompareRecAllenRichPixels.cpp b/Rec/Allen/src/CompareRecAllenRichPixels.cpp new file mode 100644 index 0000000000000000000000000000000000000000..675fddde9245cf57916210c120d91cce321dc606 --- /dev/null +++ b/Rec/Allen/src/CompareRecAllenRichPixels.cpp @@ -0,0 +1,260 @@ +/*****************************************************************************\ * (c) Copyright 2018-2020 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. * + \*****************************************************************************/ + +// Gaudi +#include "GaudiAlg/Consumer.h" +#include "Gaudi/Accumulators.h" + +// Rec +#include "RichFutureRecEvent/RichRecSIMDPixels.h" + +// Allen +#include "RichMakePixels.cuh" +#include "RichPixel.cuh" + +// std +#include +#include + +using AllenRichPixel = Allen::Rich::PixelReco::Pixel; + +enum ReturnState { IS_NULL, NOT_EXISTS, EXISTS }; + +void printPixelAttributes( + const std::string& label, + float gx, + float gy, + float gz, + float lx, + float ly, + float lz, + uint32_t smartIDKey, + float effArea, + int rich, + int side, + float timeWindow, + bool isInnerRegion) +{ + std::cout << std::fixed << std::setprecision(std::numeric_limits::max_digits10); + std::cout << label << ": "; + std::cout << "GP=(" << gx << "," << gy << "," << gz << "), "; + std::cout << "LP=(" << lx << "," << ly << "," << lz << "), "; + std::cout << "SID=" << smartIDKey << ", "; + std::cout << "EA=" << effArea << ", "; + std::cout << "R=" << rich << ", "; + std::cout << "S=" << side << ", "; + std::cout << "TW=" << timeWindow << ", "; + std::cout << "IR=" << isInnerRegion << "\n"; +} + +// This test verifies that all valid HLT2 pixels exist in Allen, and all valid Allen pixels exist in HLT2 +// Having the same number of m_allen_found_in_hlt2, and m_hlt2_found_in_allen means success. +class CompareRecAllenRichPixels final : public Gaudi::Functional::Consumer&, + const std::vector&, + const Rich::Future::Rec::SIMDPixelSummaries&)> { + +public: + /// Standard constructor + CompareRecAllenRichPixels(const std::string& name, ISvcLocator* pSvcLocator); + + /// Algorithm execution + void operator()( + const std::vector&, + const std::vector&, + const Rich::Future::Rec::SIMDPixelSummaries&) const override; + + /// Compare the attributes of an Allen Pixel and a Rec Pixel + bool matchPixels(const AllenRichPixel& allenPixel, const Rich::Future::Rec::SIMDPixel& recPixelSummary, size_t i) + const + { + auto equal = [](float a, float b, float tol = 1e-3f) { return std::abs(a - b) < tol; }; + + return ( + equal(allenPixel.gloPos().x, recPixelSummary.gloPos().X()[i]) && + equal(allenPixel.gloPos().y, recPixelSummary.gloPos().Y()[i]) && + equal(allenPixel.gloPos().z, recPixelSummary.gloPos().Z()[i]) && + + equal(allenPixel.locPos().x, recPixelSummary.locPos().X()[i]) && + equal(allenPixel.locPos().y, recPixelSummary.locPos().Y()[i]) && + equal(allenPixel.locPos().z, recPixelSummary.locPos().Z()[i]) && + + allenPixel.smartID().key() == recPixelSummary.smartID()[i].key() && + equal(allenPixel.effArea(), recPixelSummary.effArea()[i]) && + + static_cast(allenPixel.rich()) == static_cast(recPixelSummary.rich()) && + static_cast(allenPixel.side()) == static_cast(recPixelSummary.side()) && + + equal(allenPixel.timeWindow(), recPixelSummary.timeWindow()[i]) && + allenPixel.isInnerRegion() == recPixelSummary.isInnerRegion()[i]); + } + +private: + mutable Gaudi::Accumulators::Counter<> m_allen_in_rec {this, "Allen Pixels found in HLT2"}; + mutable Gaudi::Accumulators::Counter<> m_allen_not_in_rec {this, "Allen Pixels not found in HLT2"}; + mutable Gaudi::Accumulators::Counter<> m_rec_in_allen {this, "HLT2 Pixels found in Allen"}; + mutable Gaudi::Accumulators::Counter<> m_rec_not_in_allen {this, "HLT2 Pixels not found in Allen"}; + mutable Gaudi::Accumulators::Counter<> m_allen_null {this, "Null Allen Pixels"}; + mutable Gaudi::Accumulators::Counter<> m_rec_null {this, "Null HLT2 Pixels"}; + mutable Gaudi::Accumulators::Counter<> m_allen_reviewed {this, "Allen Pixels reviewed"}; + mutable Gaudi::Accumulators::Counter<> m_rec_reviewed {this, "HLT2 Pixels reviewed"}; +}; + +DECLARE_COMPONENT(CompareRecAllenRichPixels) + +CompareRecAllenRichPixels::CompareRecAllenRichPixels(const std::string& name, ISvcLocator* pSvcLocator) : + Consumer( + name, + pSvcLocator, + {KeyValue {"rich1_pixels", ""}, KeyValue {"rich2_pixels", ""}, KeyValue {"SIMDPixelSummaries", ""}}) +{} + +// When reading this code, keep in mind that recPixelSummaries contain multiple pixels, while allenRichPixels contain +// individual pixels +void CompareRecAllenRichPixels::operator()( + const std::vector& allenRich1Pixels, + const std::vector& allenRich2Pixels, + const Rich::Future::Rec::SIMDPixelSummaries& recPixelSummaries) const +{ + // Concatenate Allen Pixel vectors + std::vector allenPixels; + allenPixels.reserve(allenRich1Pixels.size() + allenRich2Pixels.size()); + allenPixels.insert(allenPixels.end(), allenRich1Pixels.begin(), allenRich1Pixels.end()); + allenPixels.insert(allenPixels.end(), allenRich2Pixels.begin(), allenRich2Pixels.end()); + + // functor to check if allen pixels exist in HLT2 + auto allenPixelExistsInRec = [&](const AllenRichPixel& allenPixel) -> ReturnState { + ++m_allen_reviewed; + // Check if pixel is valid + if (allenPixel.isNull()) { + ++m_allen_null; + return ReturnState::IS_NULL; + } + + // iterate over all pixel summaries + for (const auto& recPixelSummary : recPixelSummaries) { + // arbitralilly use any of the vectors in a pixel summary to get the pixel count and iterate that many times. + for (size_t i = 0; i < recPixelSummary.gloPos().X().size(); i++) { + // ensure HLT2 pixel validity + if (recPixelSummary.validMask()[i]) { + // check for match + if (matchPixels(allenPixel, recPixelSummary, i)) { + ++m_allen_in_rec; + return ReturnState::EXISTS; + } + } // invalid HLT2 pix + } // didn't find pix match + } // covered all HLT2 pixels + ++m_allen_not_in_rec; + error() << "Allen pixel " << allenPixel.smartID().key() << " not found in HLT2" << endmsg; + return ReturnState::NOT_EXISTS; + }; + + // functor to check if HLT2 pixels exist in Allen + auto recPixelExistsInAllen = [&](const Rich::Future::Rec::SIMDPixel& recPixelSummary) -> std::vector { + std::vector summaryStates; + bool found_current_pixel = true; + + // arbitralilly use any of the vector in a pixel summary to get the pixel count and iterate that many times. + for (size_t i = 0; i < recPixelSummary.gloPos().X().size(); i++) { + ++m_rec_reviewed; + if (found_current_pixel) { + found_current_pixel = false; + // ensure HLT2 pixel validity + if (recPixelSummary.validMask()[i]) { + // iterate over Allen pixels + for (const auto& allenPixel : allenPixels) { + // check allen pix validity + if (!allenPixel.isNull()) { + // check for match + if (matchPixels(allenPixel, recPixelSummary, i)) { + ++m_rec_in_allen; + found_current_pixel = true; + summaryStates.push_back(ReturnState::EXISTS); + continue; + } // found a match + } // ignore null Allen pix + } // covered all Allen pixels + } + else { + ++m_rec_null; + found_current_pixel = true; // assume correctness on invalid pix to ignore it. + + summaryStates.push_back(ReturnState::IS_NULL); + } // invalid HLT2 pix + } + else { // didn't find pix match + ++m_rec_not_in_allen; + summaryStates.push_back(ReturnState::NOT_EXISTS); + error() << "HLT2 pixel " << recPixelSummary.smartID()[i].key() << " not found in Allen" << endmsg; + } + } + return summaryStates; + }; + + // call allen in hlt2 functor + for (auto& allenPixel : allenPixels) { + ReturnState state = allenPixelExistsInRec(allenPixel); + if (state == ReturnState::NOT_EXISTS) { + printPixelAttributes( + "Allen", + allenPixel.gloPos().x, + allenPixel.gloPos().y, + allenPixel.gloPos().z, + allenPixel.locPos().x, + allenPixel.locPos().y, + allenPixel.locPos().z, + allenPixel.smartID().key(), + allenPixel.effArea(), + static_cast(allenPixel.rich()), + static_cast(allenPixel.side()), + allenPixel.timeWindow(), + allenPixel.isInnerRegion()); + } + } + + // call hlt2 in allen functor + + for (auto& recPixelSummary : recPixelSummaries) { + std::vector summaryStates = recPixelExistsInAllen(recPixelSummary); + for (size_t i = 0; i < summaryStates.size(); ++i) { + if (summaryStates[i] == ReturnState::NOT_EXISTS) { + printPixelAttributes( + "Rec", + recPixelSummary.gloPos().X()[i], + recPixelSummary.gloPos().Y()[i], + recPixelSummary.gloPos().Z()[i], + recPixelSummary.locPos().X()[i], + recPixelSummary.locPos().Y()[i], + recPixelSummary.locPos().Z()[i], + recPixelSummary.smartID()[i].key(), + recPixelSummary.effArea()[i], + static_cast(recPixelSummary.rich()), + static_cast(recPixelSummary.side()), + recPixelSummary.timeWindow()[i], + recPixelSummary.isInnerRegion()[i]); + } + } + } + + // verify that all Allen pixels were accounted for + if ((m_allen_in_rec.value() + m_allen_null.value()) != m_allen_reviewed.value()) { + error() << "Found " << m_allen_in_rec.value() << " Allen pixels in HLT2, and " << m_allen_null.value() + << " Allen null pixels totalling " << m_allen_null.value() + m_allen_in_rec.value() << ". Expected " + << m_allen_reviewed.value() << endmsg; + } + // verify that all HLT2 pixels were accounted for + if ((m_rec_in_allen.value() + m_rec_null.value()) != m_rec_reviewed.value()) { + error() << "Found " << m_rec_in_allen.value() << " HLT2 pixels in Allen, and " << m_rec_null.value() + << " HLT2 null pixels totalling " << m_rec_null.value() + m_rec_in_allen.value() << ". Expected " + << m_rec_reviewed.value() << endmsg; + } +} diff --git a/Rec/Allen/src/TestAllenRichPixels.cpp b/Rec/Allen/src/CompareRecAllenRichSmartIDs.cpp similarity index 75% rename from Rec/Allen/src/TestAllenRichPixels.cpp rename to Rec/Allen/src/CompareRecAllenRichSmartIDs.cpp index 98d6f414c23e199bc5a61ec54f85e2c901c24f41..e30e7839b21f29b98ec6e5ce01203ea57f0372f5 100644 --- a/Rec/Allen/src/TestAllenRichPixels.cpp +++ b/Rec/Allen/src/CompareRecAllenRichSmartIDs.cpp @@ -19,25 +19,27 @@ #include "RichDecoding.cuh" #include "RichPDMDBDecodeMapping.cuh" -class TestAllenRichPixels final : public Gaudi::Functional::Consumer&, - const std::vector&, - const Rich::Future::DAQ::DecodedData&)> { +using AllenRichSmartID = Allen::Rich::Decoding::SmartID; + +class CompareRecAllenRichSmartIDs final : public Gaudi::Functional::Consumer&, + const std::vector&, + const Rich::Future::DAQ::DecodedData&)> { public: /// Standard constructor - TestAllenRichPixels(const std::string& name, ISvcLocator* pSvcLocator); + CompareRecAllenRichSmartIDs(const std::string& name, ISvcLocator* pSvcLocator); /// Algorithm execution void operator()( - const std::vector&, - const std::vector&, + const std::vector&, + const std::vector&, const Rich::Future::DAQ::DecodedData&) const override; }; -DECLARE_COMPONENT(TestAllenRichPixels) +DECLARE_COMPONENT(CompareRecAllenRichSmartIDs) -TestAllenRichPixels::TestAllenRichPixels(const std::string& name, ISvcLocator* pSvcLocator) : +CompareRecAllenRichSmartIDs::CompareRecAllenRichSmartIDs(const std::string& name, ISvcLocator* pSvcLocator) : Consumer( name, pSvcLocator, @@ -46,9 +48,9 @@ TestAllenRichPixels::TestAllenRichPixels(const std::string& name, ISvcLocator* p KeyValue {"RichDecodedData", Rich::Future::DAQ::DecodedDataLocation::Default}}) {} -void TestAllenRichPixels::operator()( - const std::vector& allen_rich1_smart_ids, - const std::vector& allen_rich2_smart_ids, +void CompareRecAllenRichSmartIDs::operator()( + const std::vector& allen_rich1_smart_ids, + const std::vector& allen_rich2_smart_ids, const Rich::Future::DAQ::DecodedData& rec_rich_pixels) const { auto allen_rich1_ids = allen_rich1_smart_ids; @@ -56,15 +58,15 @@ void TestAllenRichPixels::operator()( auto allen_rich2_ids = allen_rich2_smart_ids; std::sort(allen_rich2_ids.begin(), allen_rich2_ids.end(), [](auto a, auto b) { return a.key() < b.key(); }); - std::vector rec_rich1_ids; - std::vector rec_rich2_ids; + std::vector rec_rich1_ids; + std::vector rec_rich2_ids; auto r1i = Rich::Rich1, r2i = Rich::Rich2; for (auto& [rD, rec_ids] : std::array {std::tie(r1i, rec_rich1_ids), std::tie(r2i, rec_rich2_ids)}) { for (const auto& pD : rec_rich_pixels[rD]) { for (const auto& mD : pD) { for (const auto& pd : mD) { std::transform(pd.smartIDs().begin(), pd.smartIDs().end(), std::back_inserter(rec_ids), [](auto id) { - return Allen::RichSmartID {id.key()}; + return AllenRichSmartID {id.key()}; }); } } @@ -78,7 +80,7 @@ void TestAllenRichPixels::operator()( if (rec_ids.size() != allen_ids.size()) { error() << "Allen and Rec " << rich << " Smart ID containers are not the same size" << endmsg; } - std::vector only_allen, only_rec; + std::vector only_allen, only_rec; std::set_difference( allen_ids.begin(), allen_ids.end(), diff --git a/backend/include/CPUBackend.h b/backend/include/CPUBackend.h index 2a3fd6369a6af6dbc95e284d425876f44eb7c083..e814912c28d7f1bbe8002afd429bf7b028a37dbd 100644 --- a/backend/include/CPUBackend.h +++ b/backend/include/CPUBackend.h @@ -37,7 +37,7 @@ using std::signbit; #define __constant__ #define __syncthreads() #define __threadfence() -#define __syncwarp() +#define __syncwarp(...) #define __launch_bounds__(_i) #define __popc __builtin_popcount #define __popcll __builtin_popcountll @@ -65,6 +65,37 @@ inline uint32_t __brev(uint32_t x) return ((x & 0xAAAAAAAA) >> 1) | ((x & 0x55555555) << 1); } +// Find the position of the n-th set to 1 bit in a 32-bit integer +inline unsigned __fns(unsigned mask, unsigned base, int offset) +{ + // https://docs.nvidia.com/cuda/parallel-thread-execution/index.html#integer-arithmetic-instructions-fns + unsigned d = 0xffffffff; + if (offset == 0) { + if (mask & (1 << base)) { + d = base; + } + } + else { + int pos = base; + int count = std::abs(offset) - 1; + int inc = (offset > 0) ? 1 : -1; + + while ((pos >= 0) && (pos < 32)) { + if (mask & (1 << pos)) { + if (count == 0) { + d = pos; + break; + } + else { + count--; + } + } + pos += inc; + } + } + return d; +} + unsigned inline __ballot_sync(unsigned mask, int predicate) { return predicate & mask; } // Support for dynamic shared memory buffers diff --git a/configuration/python/AllenConf/HLT1.py b/configuration/python/AllenConf/HLT1.py index 92f56b1177d4cd07c713d20b9fe24ce6f3adfd20..e11b6823f3bdb349b32ddefa48ba947607ab562a 100644 --- a/configuration/python/AllenConf/HLT1.py +++ b/configuration/python/AllenConf/HLT1.py @@ -1582,11 +1582,13 @@ def setup_hlt1_node(enablePhysics=True, if with_rich: hlt1_node = CompositeNode( "AllenWithRich", [ - hlt1_node, - reconstructed_objects["decoded_rich"]["dev_smart_ids"].producer + hlt1_node, reconstructed_objects["rich1_pixels"] + ["dev_rich_pixels"].producer, + reconstructed_objects["rich2_pixels"]["dev_rich_pixels"]. + producer ], NodeLogic.NONLAZY_AND, - force_order=True) + force_order=False) if enableRateValidator: hlt1_node = CompositeNode( diff --git a/configuration/python/AllenConf/hlt1_reconstruction.py b/configuration/python/AllenConf/hlt1_reconstruction.py index 4facbc4411ad9650bbc343f55d2e51772f75b1de..6e0960b8ea724d0c4c1d533428f667770b3d27bb 100644 --- a/configuration/python/AllenConf/hlt1_reconstruction.py +++ b/configuration/python/AllenConf/hlt1_reconstruction.py @@ -32,6 +32,8 @@ from AllenConf.persistency import make_gather_selections, make_sel_report_writer from AllenConf.filters import make_gec from AllenConf.best_track_creator import best_track_creator from AllenConf.enum_types import TrackingType +from AllenConf.secondary_vertex_reconstruction import make_kalman_long +from AllenConf.rich_reconstruction import make_pixels def hlt1_reconstruction(algorithm_name='', @@ -497,8 +499,13 @@ def hlt1_reconstruction(algorithm_name='', }) if with_rich: - from AllenConf.rich_reconstruction import decode_rich - output.update({"decoded_rich": decode_rich()}) + rich1_pixels = make_pixels(rich=1) + rich2_pixels = make_pixels(rich=2) + + output.update({ + "rich1_pixels": rich1_pixels, + "rich2_pixels": rich2_pixels + }) if with_AC_split: velo_tracks_A_side, velo_tracks_C_side = make_velo_tracks_ACsplit( diff --git a/configuration/python/AllenConf/rich_reconstruction.py b/configuration/python/AllenConf/rich_reconstruction.py index 3452cd084b4fb1db3c65ecd14e38fa6f13a0d441..182f02f1dc2a1e36380ff127da9306e08468b1fd 100644 --- a/configuration/python/AllenConf/rich_reconstruction.py +++ b/configuration/python/AllenConf/rich_reconstruction.py @@ -8,14 +8,20 @@ # granted to it by virtue of its status as an Intergovernmental Organization # # or submit itself to any jurisdiction. # ############################################################################### -from AllenCore.algorithms import (data_provider_t, rich_decoding_t) +from AllenCore.algorithms import (data_provider_t, rich_decoding_t, + rich_make_pixels_t) from AllenConf.utils import initialize_number_of_events from AllenCore.generator import make_algorithm +RICH_1 = 1 +RICH_2 = 2 -def decode_rich(rich=1): - if rich not in (1, 2): - raise ValueError("rich must be 1 or 2") +VALID_RICHS = [RICH_1, RICH_2] + + +def decode_rich(rich=RICH_1): + if rich not in VALID_RICHS: + raise ValueError(f"rich must be one of {VALID_RICHS}") number_of_events = initialize_number_of_events() @@ -30,10 +36,35 @@ def decode_rich(rich=1): dev_rich_raw_input_t=rich_banks.dev_raw_banks_t, dev_rich_raw_input_offsets_t=rich_banks.dev_raw_offsets_t, dev_rich_raw_input_sizes_t=rich_banks.dev_raw_sizes_t, - dev_rich_raw_input_types_t=rich_banks.dev_raw_types_t, - block_dim_x=128) + dev_rich_raw_input_types_t=rich_banks.dev_raw_types_t) return { - "dev_smart_ids": rich_decoding.dev_smart_ids_t, - "dev_rich_hit_offsets": rich_decoding.dev_rich_hit_offsets_t + "dev_smart_ids": + rich_decoding.dev_smart_ids_t, + "dev_rich_hit_offsets": + rich_decoding.dev_rich_hit_offsets_t, + "host_rich_total_number_of_hits": + rich_decoding.host_rich_total_number_of_hits_t } + + +def make_pixels(decoded_rich=None, rich=RICH_1): + if rich not in VALID_RICHS: + raise ValueError(f"rich must be one of {VALID_RICHS}") + + if decoded_rich is None: + decoded_rich = decode_rich(rich) + + number_of_hits = decoded_rich["host_rich_total_number_of_hits"] + smart_ids = decoded_rich["dev_smart_ids"] + hit_offsets = decoded_rich["dev_rich_hit_offsets"] + + rich_pixels = make_algorithm( + rich_make_pixels_t, + name=f"rich{rich}_make_pixels", + host_rich_total_number_of_hits_t=number_of_hits, + dev_smart_ids_t=smart_ids, + dev_rich_hit_offsets_t=hit_offsets, + current_rich=rich) + + return {"dev_rich_pixels": rich_pixels.dev_rich_pixels_t} diff --git a/configuration/python/AllenSequences/hlt1_pp_forward_then_matching_and_downstream_with_parkf_rich.py b/configuration/python/AllenSequences/hlt1_pp_forward_then_matching_and_downstream_with_parkf_rich.py new file mode 100644 index 0000000000000000000000000000000000000000..8f4d3e4b053fca6f85059328dcfc422fb9741b89 --- /dev/null +++ b/configuration/python/AllenSequences/hlt1_pp_forward_then_matching_and_downstream_with_parkf_rich.py @@ -0,0 +1,20 @@ +############################################################################### +# (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. # +############################################################################### +from AllenConf.HLT1 import setup_hlt1_node +from AllenCore.generator import generate +from AllenConf.enum_types import TrackingType + +hlt1_node = setup_hlt1_node( + tracking_type=TrackingType.FORWARD_THEN_MATCHING, + enableDownstream=True, + with_fullKF=True, + with_rich=True) +generate(hlt1_node) diff --git a/configuration/python/AllenSequences/rich.py b/configuration/python/AllenSequences/rich.py index 3592c2f89f80ffee849bfe2dc8a8dd268be24833..638ea281ca5e752d1d9bfd4c29c82e761df98307 100644 --- a/configuration/python/AllenSequences/rich.py +++ b/configuration/python/AllenSequences/rich.py @@ -8,13 +8,39 @@ # granted to it by virtue of its status as an Intergovernmental Organization # # or submit itself to any jurisdiction. # ############################################################################### -from AllenConf.rich_reconstruction import decode_rich +from AllenConf.rich_reconstruction import decode_rich, make_pixels from AllenCore.generator import generate from PyConf.control_flow import NodeLogic, CompositeNode -rich_decoding = CompositeNode( - "RichDecoding", [decode_rich()["dev_smart_ids"].producer], - NodeLogic.NONLAZY_AND, +decoded_rich1 = decode_rich(rich=1) +decoded_rich2 = decode_rich(rich=2) + +rich1_decoding = CompositeNode( + "Rich1Decoding", [decoded_rich1["dev_smart_ids"].producer], + NodeLogic.NONLAZY_OR, + force_order=True) + +rich2_decoding = CompositeNode( + "Rich2Decoding", [decoded_rich2["dev_smart_ids"].producer], + NodeLogic.NONLAZY_OR, + force_order=True) + +rich1_make_pixels = CompositeNode( + "Rich1MakePixels", + [make_pixels(decoded_rich1, rich=1)["dev_rich_pixels"].producer], + NodeLogic.NONLAZY_OR, + force_order=True) + +rich2_make_pixels = CompositeNode( + "Rich2MakePixels", + [make_pixels(decoded_rich2, rich=2)["dev_rich_pixels"].producer], + NodeLogic.NONLAZY_OR, + force_order=True) + +rich_node = CompositeNode( + "RichNode", + [rich1_decoding, rich1_make_pixels, rich2_decoding, rich2_make_pixels], + NodeLogic.NONLAZY_OR, force_order=True) -generate(rich_decoding) +generate(rich_node) diff --git a/device/event_model/CMakeLists.txt b/device/event_model/CMakeLists.txt index 93d7e315a51c87f67ae18933e51f0de32bf4b633..49be1f19fa712c99fd2c40a629fcbd4336fa6121 100644 --- a/device/event_model/CMakeLists.txt +++ b/device/event_model/CMakeLists.txt @@ -27,6 +27,7 @@ target_include_directories(EventModel INTERFACE $ $ $ + $ $) target_link_libraries(EventModel INTERFACE AllenCommon) diff --git a/device/event_model/rich/include/Rich.cuh b/device/event_model/rich/include/Rich.cuh new file mode 100644 index 0000000000000000000000000000000000000000..c8e47ff6f1df0022a87b7e551381695a0032d56c --- /dev/null +++ b/device/event_model/rich/include/Rich.cuh @@ -0,0 +1,30 @@ +/*****************************************************************************\ + * (c) Copyright 2023 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 "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. * + \*****************************************************************************/ +#pragma once +#include "RichDefinitions.cuh" +#include "RichPhotonDetectorPanel.cuh" + +namespace Allen::Rich { + struct RichDetector { + __host__ __device__ inline auto rich() const { return m_type; } + + __host__ __device__ inline const auto& pdPanels() const { return m_panels; } + + /// Access PD Panel for a given side + __host__ __device__ inline const auto& pdPanel(const Detector::Side side) const noexcept + { + return pdPanels()[side]; + } + + std::array m_panels {}; + Detector::Type m_type {Detector::Type::InvalidDetector}; + }; +} // namespace Allen::Rich diff --git a/device/event_model/rich/include/RichDefinitions.cuh b/device/event_model/rich/include/RichDefinitions.cuh new file mode 100644 index 0000000000000000000000000000000000000000..e0edb99a7fb7303e0293ac5ae10b0e40575b0993 --- /dev/null +++ b/device/event_model/rich/include/RichDefinitions.cuh @@ -0,0 +1,129 @@ +/*****************************************************************************\ + * (c) Copyright 2023 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 "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. * + \*****************************************************************************/ +#pragma once + +#include +#include +#include +#include "BackendCommon.h" + +namespace Allen::Rich { + using FP = float; + using Point = float3; + using DataType = std::uint32_t; + using Transform3D = std::array; + using KeyType = std::uint64_t; + using BitPackType = KeyType; + using ADCTimeType = std::uint16_t; + + //// Constants related to pixel reconstruction, value definition based on those in HLT2's + /// Rec/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDSummaryPixels.cpp + + /// Enabled 4D reconstruction + [[maybe_unused]] __device__ constexpr bool m_enable4D[2] = {false, false}; + + /// Average expected hit time for signal in each RICH in ns + [[maybe_unused]] __device__ constexpr float m_avHitTime[2] = {13.03, 52.94}; + + /// Course (pixel) Time window for each RICH in ns + [[maybe_unused]] __device__ constexpr float m_timeWindow[2] = {3.0, 3.0}; + + /// Enable the override of inner and out regions + [[maybe_unused]] __device__ constexpr bool m_overrideRegions[2] = {false, false}; + + /// Size in X defining the inner pixels for each RICH + [[maybe_unused]] __device__ constexpr double m_innerPixX[2] = {250.0, 99999.9}; + + /// Size in Y defining the inner pixels for each RICH + [[maybe_unused]] __device__ constexpr double m_innerPixY[2] = {300.0, 300.0}; + + /// Time resolution for inner regions in ns + [[maybe_unused]] __device__ constexpr float m_innerTimeWindow[2] = {0.15, 0.15}; + + /// Time resolution for outer regions in ns + [[maybe_unused]] __device__ constexpr float m_outerTimeWindow[2] = {0.3, 0.3}; + + /// Matrix-point transform3D, adapted from ROOT::Math + __host__ __device__ inline Point transform3DTimesPoint(Transform3D fM, Point point) + { + return Point {fM[0] * point.x + fM[1] * point.y + fM[2] * point.z + fM[3], + fM[4] * point.x + fM[5] * point.y + fM[6] * point.z + fM[7], + fM[8] * point.x + fM[9] * point.y + fM[10] * point.z + fM[11]}; + } +} // namespace Allen::Rich + +namespace Allen::Rich::Detector { + enum Type : std::int8_t { + InvalidDetector = -1, //< Unspecified Detector + Rich1 = 0, //< RICH1 detector + Rich2 = 1, //< RICH2 detector + Rich = 1 //< Single RICH detector + }; + + // Detector side enum + enum Side : std::int8_t { + InvalidSide = -1, //< Invalid side + // RICH1 + top = 0, //< Upper panel in RICH1 + bottom = 1, //< Lower panel in RICH1 + // RICH2 + left = 0, //< Left panel in RICH2 + right = 1, //< Right panel in RICH2 + aside = 0, //< A-Side panel in RICH2 + cside = 1, //< C-Side panel in RICH2 + // Generic + firstSide = 0, //< Upper panel in RICH1 or Left panel in RICH2 + secondSide = 1 //< Lower panel in RICH1 or Right panel in RICH2 + }; +} // namespace Allen::Rich::Detector + +namespace Allen::Rich::Decoding { + // Helper class for RichDecoding + class PackedFrameSizes final { + public: + // Packed type + using IntType = std::uint8_t; + + private: + // Bits for each Size + static const IntType Bits0 = 4; + static const IntType Bits1 = 4; + // shifts + static const IntType Shift0 = 0; + static const IntType Shift1 = Shift0 + Bits0; + // masks + static const IntType Mask0 = (IntType)((1 << Bits0) - 1) << Shift0; + static const IntType Mask1 = (IntType)((1 << Bits1) - 1) << Shift1; + // max values + static const IntType Max0 = (1 << Bits0) - 1; + static const IntType Max1 = (1 << Bits1) - 1; + + public: + // Contructor from a single word + __host__ __device__ explicit PackedFrameSizes(const IntType d) : m_data(d) {} + + // Get the overall data + __host__ __device__ inline IntType data() const noexcept { return m_data; } + + // Get first size word + __host__ __device__ inline IntType size0() const noexcept { return ((data() & Mask0) >> Shift0); } + + // Get second size word + __host__ __device__ inline IntType size1() const noexcept { return ((data() & Mask1) >> Shift1); } + + // Get the total size + __host__ __device__ inline auto totalSize() const noexcept { return size0() + size1(); } + + private: + // The data word + IntType m_data {0}; + }; +} // namespace Allen::Rich::Decoding diff --git a/device/rich/decoding/include/RichPDMDBDecodeMapping.cuh b/device/event_model/rich/include/RichPDMDBDecodeMapping.cuh similarity index 80% rename from device/rich/decoding/include/RichPDMDBDecodeMapping.cuh rename to device/event_model/rich/include/RichPDMDBDecodeMapping.cuh index 79ca73274d7487684448e632304a1988e16fd82c..0be24a92285957d85eea0dbd2adeece100ead521 100644 --- a/device/rich/decoding/include/RichPDMDBDecodeMapping.cuh +++ b/device/event_model/rich/include/RichPDMDBDecodeMapping.cuh @@ -21,16 +21,14 @@ #include #include -namespace Rich::Future::DAQ::Allen { +namespace Allen::Rich::Decoding { /// Helper class for RICH PDMDB readout mapping - class PDMDBDecodeMapping final { - public: + struct PDMDBDecodeMapping final { // data types /// The data for each anode - class BitData final { - public: + struct BitData final { /// The EC number (0-3) int8_t ec; /// The PMT number in EC @@ -38,7 +36,6 @@ namespace Rich::Future::DAQ::Allen { /// The Anode index (0-63) int8_t anode; - public: /// Default constructor BitData() = default; /// Constructor from values @@ -51,7 +48,8 @@ namespace Rich::Future::DAQ::Allen { {} }; - private: + using FrameBitmask = std::array; // first half (6bytes) + second half (5bytes) + // defines /// Max Number of frames per PDMDB @@ -78,7 +76,6 @@ namespace Rich::Future::DAQ::Allen { /// R-Type Module data for each RICH using RTypeRichData = std::array; - private: // methods /// Get the PDMDB data for given RICH, PDMDB and frame @@ -100,7 +97,6 @@ namespace Rich::Future::DAQ::Allen { } } - public: /// mapping version __host__ __device__ inline auto version() const { return m_mappingVer; } @@ -110,11 +106,22 @@ namespace Rich::Future::DAQ::Allen { /// Get PDMDB data for given Tel40 data __host__ __device__ inline const auto& getFrameData(const Tel40CableMapping::Tel40LinkData& cData) const { - const auto rich = cData.smartID.getData(::Allen::RichSmartID::ShiftRich, ::Allen::RichSmartID::MaskRich); + const auto rich = cData.smartID.getData(SmartID::ShiftRich, SmartID::MaskRich); return getFrameData(rich, cData.pdmdbNum, cData.linkNum, cData.isHType); } - private: + __host__ __device__ inline const auto& getFrameValidMask(const Tel40CableMapping::Tel40LinkData& cData) const + { + const auto rich = cData.smartID.getData(SmartID::ShiftRich, SmartID::MaskRich); + if (!cData.isHType) { + // R type PMT + return m_pdmMaskR[rich][cData.pdmdbNum][cData.linkNum]; + } + else { + return m_pdmMaskH[cData.pdmdbNum][cData.linkNum]; + } + } + // data /// R type data @@ -128,6 +135,9 @@ namespace Rich::Future::DAQ::Allen { /// Mapping version int m_mappingVer {-1}; + + std::array, PDMDBPerModule>, 2> m_pdmMaskR {}; + std::array, PDMDBPerModule> m_pdmMaskH {}; }; -} // namespace Rich::Future::DAQ::Allen +} // namespace Allen::Rich::Decoding diff --git a/device/event_model/rich/include/RichPhotonDetector.cuh b/device/event_model/rich/include/RichPhotonDetector.cuh new file mode 100644 index 0000000000000000000000000000000000000000..d04b158a9f3226f3928ecf519cb6464640636b61 --- /dev/null +++ b/device/event_model/rich/include/RichPhotonDetector.cuh @@ -0,0 +1,87 @@ +/*****************************************************************************\ + * (c) Copyright 2023 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 "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. * + \*****************************************************************************/ + +#pragma once + +#include +#include +#include + +namespace Allen::Rich::Detector { + struct PhotonDetector { + __host__ __device__ inline auto pdSmartID() const { return m_pdSmartID; } + + __host__ __device__ inline auto panelDetectionPoint(const Decoding::SmartID id) const noexcept + { + const auto fPixCol = id.pixelCol(); + const auto fPixRow = id.pixelRow(); + const auto xh = (fPixCol - m_numPixColFrac) * m_effectivePixelXSize; + const auto yh = (fPixRow - m_numPixRowFrac) * m_effectivePixelYSize; + return Point {xh, yh, m_localZcoord}; + } + + __host__ __device__ inline auto globalDetectionPoint(const Decoding::SmartID id) const noexcept + { + return transform3DTimesPoint(localToGlobal(), panelDetectionPoint(id)); + } + + __host__ __device__ inline float getEffectivePixelXSize() const noexcept { return m_effectivePixelXSize; } + __host__ __device__ inline float getEffectivePixelYSize() const noexcept { return m_effectivePixelYSize; } + __host__ __device__ inline float getLocalZcoord() const noexcept { return m_localZcoord; } + __host__ __device__ inline float getNumPixRowFrac() const { return m_numPixRowFrac; } + __host__ __device__ inline float getNumPixColFrac() const noexcept { return m_numPixColFrac; } + + /// Effective pixel area + __host__ __device__ inline auto effectivePixelArea() const noexcept { return m_effPixelArea; } + + /// Access the local to global transform + __host__ __device__ inline const std::array& localToGlobal() const noexcept { return m_locToGloM; } + + std::string toString() const + { + std::ostringstream oss; + oss << "PD { " + << "m_locToGloM: " << m_locToGloM[0] << " " << m_locToGloM[1] << " " << m_locToGloM[2] << " " + << m_locToGloM[3] << " " << m_locToGloM[4] << " " << m_locToGloM[5] << " " << m_locToGloM[6] << " " + << m_locToGloM[7] << " " << m_locToGloM[8] << " " << m_locToGloM[9] << " " << m_locToGloM[10] << " " + << m_locToGloM[11] << " " + << ", " + << "m_zeroInPanelFrame: " << m_zeroInPanelFrame[0] << m_zeroInPanelFrame[1] << m_zeroInPanelFrame[2] << ", " + << "effPixelArea: " << m_effPixelArea << ", " + << "numPixels: " << m_numPixels << ", " + << "localZcoord: " << m_localZcoord << ", " + << "NumPixColFrac: " << m_numPixColFrac << ", " + << "NumPixRowFrac: " << m_numPixRowFrac << ", " + << "EffectivePixelXSize: " << m_effectivePixelXSize << ", " + << "EffectivePixelYSize: " << m_effectivePixelYSize << ", " + << "pdSmartID: " << m_pdSmartID << ", " + << "isHType: " << (m_isHType ? "true" : "false") << "isNull: " << (m_isNull ? "true" : "false") << " }"; + return oss.str(); + } + + __host__ __device__ bool getIsNull() const { return m_isNull; } + + __host__ __device__ void setIsNull(bool value) { m_isNull = value; } + + std::array m_locToGloM {}; + std::array m_zeroInPanelFrame {}; + std::uint64_t m_pdSmartID {0}; + float m_effPixelArea {0.f}; + float m_numPixels {0.f}; + float m_localZcoord {0.f}; + float m_numPixColFrac {0.f}; + float m_numPixRowFrac {0.f}; + float m_effectivePixelXSize {0.f}; + float m_effectivePixelYSize {0.f}; + bool m_isHType {false}; + bool m_isNull {false}; + }; +} // namespace Allen::Rich::Detector diff --git a/device/event_model/rich/include/RichPhotonDetectorPanel.cuh b/device/event_model/rich/include/RichPhotonDetectorPanel.cuh new file mode 100644 index 0000000000000000000000000000000000000000..005eb70033d4bdf87a5950fc3a64289a85845bf7 --- /dev/null +++ b/device/event_model/rich/include/RichPhotonDetectorPanel.cuh @@ -0,0 +1,88 @@ +/*****************************************************************************\ + * (c) Copyright 2023 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 "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. * + \*****************************************************************************/ +#pragma once +#include +#include +#include + +namespace Allen::Rich::Detector { + struct PDPanel { + using PDArray = std::array; + using ModuleArray = std::array; + + /// Access the RICH detector type + __host__ __device__ inline auto rich() const noexcept { return m_rich; } + + /// Access the PD panel side + __host__ __device__ inline auto side() const noexcept { return m_side; } + + __host__ __device__ PhotonDetector dePD(const Decoding::SmartID smartID) const + { + PhotonDetector result; + result.setIsNull(true); + + if (!smartID.pdIsSet()) return result; + + const auto localModNum = smartID.pdMod() - m_modNumOffset; + + const bool mod_in_bounds = (localModNum < pdModules().size()); + const auto& mod = + mod_in_bounds ? pdModules()[localModNum] : pdModules()[0]; // Fallback to first (dummy) module if OOB + + const bool pd_in_bounds = (smartID.pdNumInMod() < mod.size()); + result = pd_in_bounds ? mod[smartID.pdNumInMod()] : result; // Keep default if OOB + + return result; + } + + /// Access Panel ID + const std::array& getPanelID() const { return m_panelID; } + + /// Compute the global detection point for a given RichSmartID + __host__ __device__ inline auto globalDetectionPoint(const Decoding::SmartID id) const noexcept + { + const auto pd = dePD(id); + return (!pd.getIsNull() ? pd.globalDetectionPoint(id) : Point {0, 0, 0}); + } + + /// Access the global to local transform + __host__ __device__ inline const auto& globalToPDPanel() const noexcept { return m_gloToPDPanelM; } + + /// Access all owned PD Modules + __host__ __device__ const ModuleArray& pdModules() const noexcept { return m_PDs; } + + std::string toString() const + { + std::stringstream ss; + for (size_t i = 0; i < m_PDs.size(); ++i) { + ss << "RICH " << (int) m_rich << ", Side " << (int) m_side << ", Module No. " << i << ", PDs: "; + for (size_t j = 0; j < m_PDs[i].size(); ++j) { + if (!(m_PDs[i][j].getIsNull())) { + auto pd = m_PDs[i][j]; + ss << pd.toString() << '\n'; + } + else { + ss << " 0\n"; + } + } + ss << "\n"; + } + return ss.str(); + } + + ModuleArray m_PDs {}; + std::array m_gloToPDPanelM {}; // 3D Transform + uint32_t m_modNumOffset {0}; + std::array m_panelID {Type::InvalidDetector, Side::InvalidSide, -1}; + Type m_rich {Type::InvalidDetector}; + Side m_side {Side::InvalidSide}; + }; +} // namespace Allen::Rich::Detector diff --git a/device/event_model/rich/include/RichPixel.cuh b/device/event_model/rich/include/RichPixel.cuh new file mode 100644 index 0000000000000000000000000000000000000000..c2c4c6a5dfc87aa14ad0567ea7922a17df824b06 --- /dev/null +++ b/device/event_model/rich/include/RichPixel.cuh @@ -0,0 +1,93 @@ +/*****************************************************************************\ + * (c) Copyright 2023 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 "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. * + \*****************************************************************************/ +#pragma once +#include +#include + +namespace Allen::Rich::PixelReco { + class Pixel { + public: + // Default constructor + __host__ __device__ Pixel() {} + + // 3D + __host__ __device__ Pixel(const Point& g, const Point& l, const Decoding::SmartID& sid, float area) : + gPos(g), lPos(l), m_smartID(sid), m_effArea(area), m_timeWindow(0.0f) + { + m_flags.setIsNull(false); + m_flags.setIsInnerRegion(!sid.isLargePMT()); + } + + // 4D constructor + __host__ __device__ Pixel(const Point& g, const Point& l, const Decoding::SmartID& sid, float area, float time) : + gPos(g), lPos(l), m_smartID(sid), m_effArea(area), m_timeWindow(time) + { + m_flags.setIsNull(false); + m_flags.setIsInnerRegion(!sid.isLargePMT()); + } + + __host__ __device__ bool isNull() const { return m_flags.isNull(); } + __host__ __device__ bool isInnerRegion() const { return m_flags.isInnerRegion(); } + __host__ __device__ void setIsNull(bool v) { m_flags.setIsNull(v); } + __host__ __device__ void overrideRegions(bool inner) { m_flags.setIsInnerRegion(inner); } + + __host__ __device__ inline auto rich() const { return m_smartID.rich(); } + __host__ __device__ inline auto side() const { return m_smartID.side(); } + + __host__ __device__ auto smartID() const { return m_smartID; } + __host__ __device__ auto gloPos() const { return gPos; } + __host__ __device__ auto locPos() const { return lPos; } + __host__ __device__ auto effArea() const { return m_effArea; } + __host__ __device__ auto timeWindow() const { return m_timeWindow; } + + __host__ std::string toString() + { + std::stringstream ss; + ss << m_smartID.key() << ',' << gPos.x << ',' << gPos.y << ',' << gPos.z << ',' << lPos.x << ',' << lPos.y << ',' + << lPos.z << ',' << m_effArea << ',' << (int) smartID().rich() << ',' << (int) smartID().side() << ',' + << m_timeWindow << ',' << isInnerRegion() << ',' << isNull() << '\n'; + return ss.str(); + } + + private: + Point gPos {0.0f, 0.0f, 0.0f}; + Point lPos {0.0f, 0.0f, 0.0f}; + Decoding::SmartID m_smartID {}; + float m_effArea {0.0f}; + float m_timeWindow {0.0f}; + + struct Flags { + uint8_t value = 3; // all bits set to 1 by default (isInnerRegion=true, isNull=true) + static constexpr uint8_t IsInnerRegionMask = 1 << 0; + static constexpr uint8_t IsNullMask = 1 << 1; + + __host__ __device__ bool isInnerRegion() const { return value & IsInnerRegionMask; } + + __host__ __device__ bool isNull() const { return value & IsNullMask; } + + __host__ __device__ void setIsInnerRegion(bool v) + { + if (v) + value |= IsInnerRegionMask; + else + value &= ~IsInnerRegionMask; + } + + __host__ __device__ void setIsNull(bool v) + { + if (v) + value |= IsNullMask; + else + value &= ~IsNullMask; + } + } m_flags; + }; +} // namespace Allen::Rich::PixelReco diff --git a/device/event_model/rich/include/RichSmartID.cuh b/device/event_model/rich/include/RichSmartID.cuh new file mode 100644 index 0000000000000000000000000000000000000000..e54b71d7bb35ac94db6b8bac2bc887a07440b472 --- /dev/null +++ b/device/event_model/rich/include/RichSmartID.cuh @@ -0,0 +1,336 @@ +/*****************************************************************************\ + * (c) Copyright 2023 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 "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. * + \*****************************************************************************/ +#pragma once +#include +#include +#include +#include "BackendCommon.h" +#include + +namespace Allen::Rich::Decoding { + class SmartID { + KeyType m_key; + + public: + // Set the RICH PD pixel row identifier + __host__ __device__ constexpr void setPixelRow(const DataType row) + { + setData(row, ShiftPixelRow, MaskPixelRow, MaskPixelRowIsSet); + } + + // Set the RICH PD pixel column identifier + __host__ __device__ constexpr void setPixelCol(const DataType col) + { + setData(col, ShiftPixelCol, MaskPixelCol, MaskPixelColIsSet); + } + + __host__ __device__ SmartID() = default; + + __host__ __device__ SmartID(KeyType key) : m_key(key) {} + + // Number of bits for each data field in the word + static constexpr const BitPackType BitsPixelCol = 3; //< Number of bits for MaPMT pixel column + static constexpr const BitPackType BitsPixelRow = 3; //< Number of bits for MaPMT pixel row + static constexpr const BitPackType BitsPDNumInMod = 4; //< Number of bits for MaPMT 'number in module' + static constexpr const BitPackType BitsPDMod = 9; //< Number of bits for MaPMT module + static constexpr const BitPackType BitsPanel = 1; //< Number of bits for MaPMT panel + static constexpr const BitPackType BitsRich = 1; //< Number of bits for RICH detector + static constexpr const BitPackType BitsPixelSubRowIsSet = 1; + static constexpr const BitPackType BitsPixelColIsSet = 1; + static constexpr const BitPackType BitsPixelRowIsSet = 1; + static constexpr const BitPackType BitsPDIsSet = 1; + static constexpr const BitPackType BitsPanelIsSet = 1; + static constexpr const BitPackType BitsRichIsSet = 1; + static constexpr const BitPackType BitsLargePixel = 1; + + // The shifts + static constexpr const BitPackType ShiftPixelCol = 0; + static constexpr const BitPackType ShiftPixelRow = ShiftPixelCol + BitsPixelCol; + static constexpr const BitPackType ShiftPDNumInMod = ShiftPixelRow + BitsPixelRow; + static constexpr const BitPackType ShiftPDMod = ShiftPDNumInMod + BitsPDNumInMod; + static constexpr const BitPackType ShiftPanel = ShiftPDMod + BitsPDMod; + static constexpr const BitPackType ShiftRich = ShiftPanel + BitsPanel; + static constexpr const BitPackType ShiftPixelSubRowIsSet = ShiftRich + BitsRich; + static constexpr const BitPackType ShiftPixelColIsSet = ShiftPixelSubRowIsSet + BitsPixelSubRowIsSet; + static constexpr const BitPackType ShiftPixelRowIsSet = ShiftPixelColIsSet + BitsPixelColIsSet; + static constexpr const BitPackType ShiftPDIsSet = ShiftPixelRowIsSet + BitsPixelRowIsSet; + static constexpr const BitPackType ShiftPanelIsSet = ShiftPDIsSet + BitsPDIsSet; + static constexpr const BitPackType ShiftRichIsSet = ShiftPanelIsSet + BitsPanelIsSet; + static constexpr const BitPackType ShiftLargePixel = ShiftRichIsSet + BitsRichIsSet; + + // The masks + static constexpr const BitPackType MaskPixelCol = (BitPackType)((BitPackType {1} << BitsPixelCol) - BitPackType {1}) + << ShiftPixelCol; + static constexpr const BitPackType MaskPixelRow = (BitPackType)((BitPackType {1} << BitsPixelRow) - BitPackType {1}) + << ShiftPixelRow; + static constexpr const BitPackType MaskPDNumInMod = + (BitPackType)((BitPackType {1} << BitsPDNumInMod) - BitPackType {1}) << ShiftPDNumInMod; + static constexpr const BitPackType MaskPDMod = (BitPackType)((BitPackType {1} << BitsPDMod) - BitPackType {1}) + << ShiftPDMod; + static constexpr const BitPackType MaskPanel = (BitPackType)((BitPackType {1} << BitsPanel) - BitPackType {1}) + << ShiftPanel; + static constexpr const BitPackType MaskRich = (BitPackType)((BitPackType {1} << BitsRich) - BitPackType {1}) + << ShiftRich; + static constexpr const BitPackType MaskPixelSubRowIsSet = + (BitPackType)((BitPackType {1} << BitsPixelSubRowIsSet) - BitPackType {1}) << ShiftPixelSubRowIsSet; + static constexpr const BitPackType MaskPixelColIsSet = + (BitPackType)((BitPackType {1} << BitsPixelColIsSet) - BitPackType {1}) << ShiftPixelColIsSet; + static constexpr const BitPackType MaskPixelRowIsSet = + (BitPackType)((BitPackType {1} << BitsPixelRowIsSet) - BitPackType {1}) << ShiftPixelRowIsSet; + static constexpr const BitPackType MaskPDIsSet = (BitPackType)((BitPackType {1} << BitsPDIsSet) - BitPackType {1}) + << ShiftPDIsSet; + static constexpr const BitPackType MaskPanelIsSet = + (BitPackType)((BitPackType {1} << BitsPanelIsSet) - BitPackType {1}) << ShiftPanelIsSet; + static constexpr const BitPackType MaskRichIsSet = + (BitPackType)((BitPackType {1} << BitsRichIsSet) - BitPackType {1}) << ShiftRichIsSet; + static constexpr const BitPackType MaskLargePixel = + (BitPackType)((BitPackType {1} << BitsLargePixel) - BitPackType {1}) << ShiftLargePixel; + + // Max values + static constexpr const DataType MaxPixelCol = (DataType)(BitPackType {1} << BitsPixelCol) - DataType {1}; + static constexpr const DataType MaxPixelRow = (DataType)(BitPackType {1} << BitsPixelRow) - DataType {1}; + static constexpr const DataType MaxPDNumInMod = (DataType)(BitPackType {1} << BitsPDNumInMod) - DataType {1}; + static constexpr const DataType MaxPDMod = (DataType)(BitPackType {1} << BitsPDMod) - DataType {1}; + static constexpr const DataType MaxPanel = (DataType)(BitPackType {1} << BitsPanel) - DataType {1}; + static constexpr const DataType MaxRich = (DataType)(BitPackType {1} << BitsRich) - DataType {1}; + + // Number of bits for the channel identification (i.e. excluding any time info) + // Currently use the lowest 32 bits for this. + static constexpr const BitPackType NChannelBits = 32; + + __host__ __device__ constexpr inline void + setData(const DataType value, const BitPackType shift, const BitPackType mask) noexcept + { + m_key = ((BitPackType {value} << shift) & mask) | (m_key & ~mask); + } + + __host__ __device__ constexpr inline void + setData(const DataType value, const BitPackType shift, const BitPackType mask, const BitPackType okMask) noexcept + { + m_key = ((BitPackType {value} << shift) & mask) | (m_key & ~mask) | okMask; + } + + __host__ __device__ constexpr inline BitPackType getData(const BitPackType shift, const BitPackType mask) const + noexcept + { + return (m_key & mask) >> shift; + } + + __host__ __device__ constexpr inline BitPackType key() const noexcept { return m_key; } + + __host__ __device__ constexpr inline bool operator==(const SmartID& other) const noexcept + { + return m_key == other.key(); + } + + __host__ __device__ constexpr inline auto isLargePMT() const noexcept + { + return 0 != getData(ShiftLargePixel, MaskLargePixel); + } + + __host__ __device__ constexpr inline auto rich() const noexcept { return getData(ShiftRich, MaskRich); } + + __host__ __device__ constexpr inline auto panel() const noexcept { return getData(ShiftPanel, MaskPanel); } + + __host__ __device__ constexpr inline auto side() const noexcept { return panel(); } + + __host__ __device__ constexpr inline DataType pdMod() const noexcept { return getData(ShiftPDMod, MaskPDMod); } + + __host__ __device__ constexpr inline DataType pdNumInMod() const noexcept + { + return getData(ShiftPDNumInMod, MaskPDNumInMod); + } + + __host__ __device__ constexpr inline DataType panelLocalModuleNum() const noexcept + { + return pdMod() - PanelModuleOffsets()[rich()][panel()]; + } + + __host__ __device__ constexpr inline DataType columnLocalModuleNum() const noexcept + { + return panelLocalModuleNum() % ModulesPerColumn; + } + + __host__ __device__ constexpr inline DataType numPMTsPerEC() const noexcept + { + return isLargePMT() ? HTypePMTsPerEC : RTypePMTsPerEC; + } + + __host__ __device__ constexpr inline DataType elementaryCell() const noexcept + { + return pdNumInMod() / numPMTsPerEC(); + } + + __host__ __device__ constexpr inline DataType pdNumInEC() const noexcept { return pdNumInMod() % numPMTsPerEC(); } + + __host__ __device__ constexpr inline auto pixelColIsSet() const noexcept + { + return 0 != getData(ShiftPixelColIsSet, MaskPixelColIsSet); + } + + __host__ __device__ constexpr inline auto pixelRowIsSet() const noexcept + { + return 0 != getData(ShiftPixelRowIsSet, MaskPixelRowIsSet); + } + + [[nodiscard]] __host__ __device__ constexpr bool pdIsSet() const noexcept + { + return 0 != ((key() & MaskPDIsSet) >> ShiftPDIsSet); + } + + __host__ __device__ constexpr inline DataType pixelCol() const noexcept + { + return getData(ShiftPixelCol, MaskPixelCol); + } + + __host__ __device__ constexpr inline DataType pixelRow() const noexcept + { + return getData(ShiftPixelRow, MaskPixelRow); + } + + __host__ __device__ constexpr inline DataType anodeIndex() const noexcept + { + return pixelRow() * PixelsPerCol + PixelsPerRow - 1 - pixelCol(); + } + + __host__ __device__ constexpr inline bool adcTimeIsSet() const noexcept + { + return (0 != (key() & MaPMT::MaskADCTimeIsSet) >> MaPMT::ShiftADCTimeIsSet); + } + + __host__ __device__ constexpr ADCTimeType adcTime() const + { + return ((key() & MaPMT::MaskADCTime) >> MaPMT::ShiftADCTime); + } + + __host__ __device__ constexpr auto time() const noexcept + { + return (MaPMT::MinTime + (adcTime() * MaPMT::ScaleADCToTime)); + } + + // ostream operator + __host__ friend std::ostream& operator<<(std::ostream& str, const SmartID& id) + { + const auto rich = id.rich(); + const auto panel = id.panel(); + str << "{ "; + str << "PMT"; + str << (id.isLargePMT() ? ":h " : ":r "); + str << (rich == 0 ? "Rich1 " : "Rich2 "); + if (rich == 0) { + str << (panel == 0 ? "Top " : "Bot "); + } + else { + str << (panel == 0 ? "L-A " : "R-C "); + } + str << "PD[Mod,NInMod]: "; + str << std::setfill('0') << std::setw(3) << id.pdMod() << ','; + str << std::setfill('0') << std::setw(2) << id.pdNumInMod() << ' '; + str << "Mod[Col,NInCol] "; + str << std::setfill('0') << std::setw(2) << id.panelLocalModuleNum() << ','; + str << std::setfill('0') << std::setw(2) << id.columnLocalModuleNum() << ' '; + str << " PD[EC,NInEC]: "; + str << id.elementaryCell() << ','; + str << id.pdNumInEC() << ' '; + + const auto pixColSet = id.pixelColIsSet(); + const auto pixRowSet = id.pixelRowIsSet(); + if (pixColSet || pixRowSet) { + if (pixColSet && pixRowSet) { + str << "Pix[Col,Row]: " << std::setfill('0') << std::setw(1) << id.pixelCol() << "," << id.pixelRow(); + } + else { + const auto fSPix = 2; + if (pixColSet) { + str << std::setfill('0') << std::setw(fSPix) << " pixCol" << id.pixelCol(); + } + if (pixRowSet) { + str << std::setfill('0') << std::setw(fSPix) << " pixRow" << id.pixelRow(); + } + } + // Include PMT derived info + if (pixColSet && pixRowSet) { + str << " Anode:" << std::setfill('0') << std::setw(2) << id.anodeIndex(); + } + } + str << " }\n"; + return str; + } + + public: + // Number of PMT pixels per row + static constexpr const DataType PixelsPerRow = 8; + // Number of PMT pixels per column + static constexpr const DataType PixelsPerCol = 8; + // Total number of PMT pixels + static constexpr const DataType TotalPixels = PixelsPerRow * PixelsPerCol; + // Number PMTs per EC for R type PMTs + static constexpr const DataType RTypePMTsPerEC = 4; + // Number PMTs per EC for H type PMTs + static constexpr const DataType HTypePMTsPerEC = 1; + // Number of ECs per module + static constexpr const DataType ECsPerModule = 4; + // Number of modules per column + static constexpr const DataType ModulesPerColumn = 6; + // Number of module columns per panel, in each RICH + static constexpr const std::array ModuleColumnsPerPanel = { +#ifdef USE_DD4HEP + {11, 14} // With dd4hep we have an extra column reserved at the end of each RICH2 panel +#else + {11, 12} +#endif + }; + // Maximum number of module columns in any panel, RICH1 or RICH2 + static constexpr const DataType MaxModuleColumnsAnyPanel = + std::max(ModuleColumnsPerPanel[0], ModuleColumnsPerPanel[1]); + // Number of modules per panel, in each RICH + static constexpr const std::array ModulesPerPanel { + {(ModulesPerColumn * ModuleColumnsPerPanel[0]), (ModulesPerColumn * ModuleColumnsPerPanel[1])}}; + + /// Number of modules in RICH1 + static constexpr const DataType RICH1Modules = 2 * ModulesPerPanel[Rich::Detector::Type::Rich1]; + /// Number of modules in RICH2 + static constexpr const DataType RICH2Modules = 2 * ModulesPerPanel[Rich::Detector::Type::Rich2]; + /// Total number of modules + static constexpr const DataType TotalModules = RICH1Modules + RICH2Modules; + + __host__ __device__ static constexpr std::array, 2> PanelModuleOffsets() + { + return {std::array {0, ModulesPerPanel[0]}, + std::array {2 * ModulesPerPanel[0], (2 * ModulesPerPanel[0]) + ModulesPerPanel[1]}}; + } + + class MaPMT { + public: + static constexpr const DataType MaxPDsPerModule = 16; + static constexpr const DataType MaxModulesPerPanel = 92; + + // Bits for time field. + static constexpr const BitPackType BitsADCTime = 16; + static constexpr const BitPackType BitsADCTimeIsSet = 1; + static constexpr const BitPackType ShiftADCTime = NChannelBits; + static constexpr const BitPackType ShiftADCTimeIsSet = ShiftADCTime + BitsADCTime; + static constexpr const BitPackType MaskADCTime = (BitPackType)((BitPackType {1} << BitsADCTime) - BitPackType {1}) + << ShiftADCTime; + static constexpr const BitPackType MaskADCTimeIsSet = + (BitPackType)((BitPackType {1} << BitsADCTimeIsSet) - BitPackType {1}) << ShiftADCTimeIsSet; + + // Max ADC time that can be stored + static constexpr const ADCTimeType MaxADCTime = static_cast((BitPackType {1} << BitsADCTime) - 1); + + // Parameters for conversion between float and ADC time values + static constexpr const float MinTime = -50.0f; // In nanoseconds + static constexpr const float MaxTime = 150.0f; // In nanoseconds + static constexpr const float ScaleTimeToADC = MaxADCTime / (MaxTime - MinTime); + static constexpr const float ScaleADCToTime = 1.0f / ScaleTimeToADC; + }; + }; // namespace SmartID +} // namespace Allen::Rich::Decoding diff --git a/device/rich/decoding/include/RichTel40CableMapping.cuh b/device/event_model/rich/include/RichTel40CableMapping.cuh similarity index 84% rename from device/rich/decoding/include/RichTel40CableMapping.cuh rename to device/event_model/rich/include/RichTel40CableMapping.cuh index 303d6051224f8fe38cfd6c21a269ae7e3b363a81..f17bed4a943561e887f75d52a81cf472b5622dcb 100644 --- a/device/rich/decoding/include/RichTel40CableMapping.cuh +++ b/device/event_model/rich/include/RichTel40CableMapping.cuh @@ -12,41 +12,34 @@ #include #include +#include -namespace Rich::Future::DAQ::Allen { +namespace Allen::Rich::Decoding { /// Helper class for RICH PMT data format encoding - class Tel40CableMapping final { - - public: + struct Tel40CableMapping final { /// Struct for storing data for each Tel40 Link - class Tel40LinkData final { - public: + struct Tel40LinkData final { /// RICH SmartID - ::Allen::RichSmartID smartID; + SmartID smartID {}; /// Module Number - int32_t moduleNum; + int32_t moduleNum {0}; /// Source ID - int16_t sourceID; + int16_t sourceID {0}; /// Tel40 connector - int8_t connector; + int8_t connector {0}; /// PDMDB number (0,1) - int8_t pdmdbNum; + int8_t pdmdbNum {0}; /// Link number - int8_t linkNum; + int8_t linkNum {0}; /// PMT type bool isHType {false}; /// Is Link Active - bool isActive {false}; - - public: - /// Default constructor - Tel40LinkData() = default; + bool isActive {false}; // TODO: remove (already stored in metadata) }; - class Tel40MetaData final { - public: - uint32_t nActiveLinks; - bool hasInactiveLinks; + struct Tel40MetaData final { + uint32_t validLinkMask {0}; + bool hasInactiveLinks {true}; // TODO: remove }; /// Max number of links(frames) per PDMDB @@ -71,12 +64,11 @@ namespace Rich::Future::DAQ::Allen { using PDMDBData = std::array; /// Tel40 data for each Module - using ModuleTel40Data = std::array; // 300 is totalmodules + using ModuleTel40Data = std::array; using Tel40SourceIDs = std::array, 164>, 2>, 2>; using Tel40SourceMetas = std::array, 2>, 2>; - public: // accessors /// Access the initialisation state @@ -119,12 +111,9 @@ namespace Rich::Future::DAQ::Allen { /// mapping version inline auto version() const { return m_mappingVer; } - private: - // data - /// Tel40 connection mapping data - Tel40SourceIDs m_tel40ConnData; - Tel40SourceMetas m_tel40ConnMeta; + Tel40SourceIDs m_tel40ConnData {}; + Tel40SourceMetas m_tel40ConnMeta {}; /// Tel40 Module Mapping data ModuleTel40Data m_tel40ModuleData; @@ -135,5 +124,4 @@ namespace Rich::Future::DAQ::Allen { /// Mapping version int m_mappingVer {-1}; }; - -} // namespace Rich::Future::DAQ::Allen \ No newline at end of file +} // namespace Allen::Rich::Decoding diff --git a/device/rich/decoding/include/RichDecoding.cuh b/device/rich/decoding/include/RichDecoding.cuh index 715330c90ce6608eff3c53b90b0434bb284d6861..d8367e3c3798a0cc09695bc6d3227779822b9502 100644 --- a/device/rich/decoding/include/RichDecoding.cuh +++ b/device/rich/decoding/include/RichDecoding.cuh @@ -11,7 +11,18 @@ #pragma once #include "AlgorithmTypes.cuh" -#include +#include + +#if defined(TARGET_DEVICE_CUDA) +#define RICH_DECODING_BLOCK_DIM_X 32 +#define RICH_DECODING_BLOCK_DIM_Y 4 +#else +// Default to CPU-like configuration +// Although Allen automatically sets this up for CPU, +// we need a variable for array sizes in compile-time +#define RICH_DECODING_BLOCK_DIM_X 1 +#define RICH_DECODING_BLOCK_DIM_Y 1 +#endif namespace rich_decoding { struct Parameters { @@ -24,7 +35,7 @@ namespace rich_decoding { DEVICE_INPUT(dev_rich_raw_input_types_t, uint) dev_rich_raw_input_types; DEVICE_OUTPUT(dev_rich_hit_offsets_t, unsigned) dev_rich_hit_offsets; HOST_OUTPUT(host_rich_total_number_of_hits_t, unsigned) host_rich_total_number_of_hits; - DEVICE_OUTPUT(dev_smart_ids_t, Allen::RichSmartID) dev_smart_ids; + DEVICE_OUTPUT(dev_smart_ids_t, Allen::Rich::Decoding::SmartID) dev_smart_ids; }; struct rich_decoding_t : public DeviceAlgorithm, Parameters { @@ -37,6 +48,14 @@ namespace rich_decoding { const Allen::Context&) const; private: - Allen::Property m_block_dim_x {this, "block_dim_x", 64, "block dimension x"}; + // Changed to use 2D block dimensions for better parallelism + // threads in x-dimension for warp-level operations + // threads in y-dimension for bank-level parallelism + Allen::Property m_block_dim_x {this, "block_dim_x", {RICH_DECODING_BLOCK_DIM_X}, "block x dimension"}; + Allen::Property m_block_dim_y {this, "block_dim_y", {RICH_DECODING_BLOCK_DIM_Y}, "block y dimension"}; + Allen::Property m_block_dim {this, + "block_dim", + {RICH_DECODING_BLOCK_DIM_X, RICH_DECODING_BLOCK_DIM_Y, 1}, + "block dimensions"}; }; } // namespace rich_decoding diff --git a/device/rich/decoding/include/RichDefinitions.cuh b/device/rich/decoding/include/RichDefinitions.cuh deleted file mode 100644 index 5cb8359fe65ac8f0dd7a43bdaadb3ecd5ec30143..0000000000000000000000000000000000000000 --- a/device/rich/decoding/include/RichDefinitions.cuh +++ /dev/null @@ -1,159 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2023 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 "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. * -\*****************************************************************************/ -#pragma once - -#include -#include -#include "BackendCommon.h" - -namespace Allen { - class RichSmartID { - uint64_t m_key; - - public: - RichSmartID() = default; - - __host__ __device__ RichSmartID(uint64_t key) : m_key(key) {} - - // Number of bits for each data field in the word - static constexpr const unsigned BitsPixelCol = 3; ///< Number of bits for MaPMT pixel column - static constexpr const unsigned BitsPixelRow = 3; ///< Number of bits for MaPMT pixel row - static constexpr const unsigned BitsPDNumInMod = 4; ///< Number of bits for MaPMT 'number in module' - static constexpr const unsigned BitsPDMod = 9; ///< Number of bits for MaPMT module - static constexpr const unsigned BitsPanel = 1; ///< Number of bits for MaPMT panel - static constexpr const unsigned BitsRich = 1; ///< Number of bits for RICH detector - static constexpr const unsigned BitsPixelSubRowIsSet = 1; - static constexpr const unsigned BitsPixelColIsSet = 1; - static constexpr const unsigned BitsPixelRowIsSet = 1; - static constexpr const unsigned BitsPDIsSet = 1; - static constexpr const unsigned BitsPanelIsSet = 1; - static constexpr const unsigned BitsRichIsSet = 1; - static constexpr const unsigned BitsLargePixel = 1; - - // The shifts - static constexpr const unsigned ShiftPixelCol = 0; - static constexpr const unsigned ShiftPixelRow = ShiftPixelCol + BitsPixelCol; - static constexpr const unsigned ShiftPDNumInMod = ShiftPixelRow + BitsPixelRow; - static constexpr const unsigned ShiftPDMod = ShiftPDNumInMod + BitsPDNumInMod; - static constexpr const unsigned ShiftPanel = ShiftPDMod + BitsPDMod; - static constexpr const unsigned ShiftRich = ShiftPanel + BitsPanel; - static constexpr const unsigned ShiftPixelSubRowIsSet = ShiftRich + BitsRich; - static constexpr const unsigned ShiftPixelColIsSet = ShiftPixelSubRowIsSet + BitsPixelSubRowIsSet; - static constexpr const unsigned ShiftPixelRowIsSet = ShiftPixelColIsSet + BitsPixelColIsSet; - static constexpr const unsigned ShiftPDIsSet = ShiftPixelRowIsSet + BitsPixelRowIsSet; - static constexpr const unsigned ShiftPanelIsSet = ShiftPDIsSet + BitsPDIsSet; - static constexpr const unsigned ShiftRichIsSet = ShiftPanelIsSet + BitsPanelIsSet; - static constexpr const unsigned ShiftLargePixel = ShiftRichIsSet + BitsRichIsSet; - - // The masks - static constexpr const unsigned MaskPixelCol = (unsigned) ((1 << BitsPixelCol) - 1) << ShiftPixelCol; - static constexpr const unsigned MaskPixelRow = (unsigned) ((1 << BitsPixelRow) - 1) << ShiftPixelRow; - static constexpr const unsigned MaskPDNumInMod = (unsigned) ((1 << BitsPDNumInMod) - 1) << ShiftPDNumInMod; - static constexpr const unsigned MaskPDMod = (unsigned) ((1 << BitsPDMod) - 1) << ShiftPDMod; - static constexpr const unsigned MaskPanel = (unsigned) ((1 << BitsPanel) - 1) << ShiftPanel; - static constexpr const unsigned MaskRich = (unsigned) ((1 << BitsRich) - 1) << ShiftRich; - static constexpr const unsigned MaskPixelSubRowIsSet = (unsigned) ((1 << BitsPixelSubRowIsSet) - 1) - << ShiftPixelSubRowIsSet; - static constexpr const unsigned MaskPixelColIsSet = (unsigned) ((1 << BitsPixelColIsSet) - 1) << ShiftPixelColIsSet; - static constexpr const unsigned MaskPixelRowIsSet = (unsigned) ((1 << BitsPixelRowIsSet) - 1) << ShiftPixelRowIsSet; - static constexpr const unsigned MaskPDIsSet = (unsigned) ((1 << BitsPDIsSet) - 1) << ShiftPDIsSet; - static constexpr const unsigned MaskPanelIsSet = (unsigned) ((1 << BitsPanelIsSet) - 1) << ShiftPanelIsSet; - static constexpr const unsigned MaskRichIsSet = (unsigned) ((1 << BitsRichIsSet) - 1) << ShiftRichIsSet; - static constexpr const unsigned MaskLargePixel = (unsigned) ((1 << BitsLargePixel) - 1) << ShiftLargePixel; - - // Max Values - static constexpr const unsigned MaxPixelCol = (unsigned) (1 << BitsPixelCol) - 1; - static constexpr const unsigned MaxPixelRow = (unsigned) (1 << BitsPixelRow) - 1; - static constexpr const unsigned MaxPDNumInMod = (unsigned) (1 << BitsPDNumInMod) - 1; - static constexpr const unsigned MaxPDMod = (unsigned) (1 << BitsPDMod) - 1; - static constexpr const unsigned MaxPanel = (unsigned) (1 << BitsPanel) - 1; - static constexpr const unsigned MaxRich = (unsigned) (1 << BitsRich) - 1; - - __host__ __device__ constexpr inline void - setData(const unsigned value, const unsigned shift, const unsigned mask) noexcept - { - m_key = ((static_cast(value) << shift) & mask) | (m_key & ~mask); - } - - __host__ __device__ constexpr inline void setData( - const unsigned value, // - const unsigned shift, // - const unsigned mask, // - const unsigned okMask) noexcept - { - m_key = ((static_cast(value) << shift) & mask) | (m_key & ~mask) | okMask; - } - - __host__ __device__ constexpr inline uint64_t getData(const unsigned shift, const unsigned mask) const noexcept - { - return (m_key & mask) >> shift; - } - - __host__ __device__ constexpr inline auto key() const noexcept { return m_key; } - - __host__ __device__ constexpr inline bool operator==(const RichSmartID& other) const noexcept - { - return m_key == other.key(); - } - - /// ostream operator - __host__ friend std::ostream& operator<<(std::ostream& str, const RichSmartID& id) { return str << id.key(); } - }; // namespace RichSmartID -} // namespace Allen - -namespace Rich::Future::DAQ { - enum DetectorType : std::int8_t { - InvalidDetector = -1, ///< Unspecified Detector - Rich1 = 0, ///< RICH1 detector - Rich2 = 1, ///< RICH2 detector - Rich = 1 ///< Single RICH detector - }; - - class PackedFrameSizes final { - public: - /// Packed type - using IntType = std::uint8_t; - - private: - // Bits for each Size - static const IntType Bits0 = 4; - static const IntType Bits1 = 4; - // shifts - static const IntType Shift0 = 0; - static const IntType Shift1 = Shift0 + Bits0; - // masks - static const IntType Mask0 = (IntType)((1 << Bits0) - 1) << Shift0; - static const IntType Mask1 = (IntType)((1 << Bits1) - 1) << Shift1; - // max values - static const IntType Max0 = (1 << Bits0) - 1; - static const IntType Max1 = (1 << Bits1) - 1; - - public: - /// Contructor from a single word - __host__ __device__ explicit PackedFrameSizes(const IntType d) : m_data(d) {} - - /// Get the overall data - __host__ __device__ inline IntType data() const noexcept { return m_data; } - - /// Get first size word - __host__ __device__ inline IntType size0() const noexcept { return ((data() & Mask0) >> Shift0); } - - /// Get second size word - __host__ __device__ inline IntType size1() const noexcept { return ((data() & Mask1) >> Shift1); } - - /// Get the total size - __host__ __device__ inline auto totalSize() const noexcept { return size0() + size1(); } - - private: - /// The data word - IntType m_data {0}; - }; -} // namespace Rich::Future::DAQ diff --git a/device/rich/decoding/include/RichMakePixels.cuh b/device/rich/decoding/include/RichMakePixels.cuh new file mode 100644 index 0000000000000000000000000000000000000000..efda2205e83a5c7642e23ce129172177840e634a --- /dev/null +++ b/device/rich/decoding/include/RichMakePixels.cuh @@ -0,0 +1,39 @@ +/*****************************************************************************\ + * (c) Copyright 2018-2020 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 "AlgorithmTypes.cuh" +#include +#include + +namespace rich_make_pixels { + struct Parameters { + HOST_INPUT(host_rich_total_number_of_hits_t, unsigned) host_rich_total_number_of_hits; + MASK_INPUT(dev_event_list_t) dev_event_list; + DEVICE_INPUT(dev_smart_ids_t, Allen::Rich::Decoding::SmartID) dev_smart_ids; + DEVICE_INPUT(dev_rich_hit_offsets_t, unsigned) dev_rich_hit_offsets; + DEVICE_OUTPUT(dev_rich_pixels_t, Allen::Rich::PixelReco::Pixel) dev_rich_pixels; + }; + + struct rich_make_pixels_t : public DeviceAlgorithm, Parameters { + void set_arguments_size(ArgumentReferences, const RuntimeOptions&, const Constants&) const; + + void operator()( + const ArgumentReferences&, + const RuntimeOptions&, + const Constants&, + const Allen::Context&) const; + + private: + Allen::Property m_block_dim {this, "block_dim", {128, 1, 1}, "block dimensions"}; + Allen::Property m_current_rich {this, "current_rich", 1, "current rich"}; + }; +} // namespace rich_make_pixels diff --git a/device/rich/decoding/src/RichDecoding.cu b/device/rich/decoding/src/RichDecoding.cu index ec1761452a50cd7b59069e1a9d433d5e1ec94178..1f0f7bb3449a084ec1e935fd9e5d244877ba7887 100644 --- a/device/rich/decoding/src/RichDecoding.cu +++ b/device/rich/decoding/src/RichDecoding.cu @@ -13,391 +13,323 @@ #include #include #include +#include INSTANTIATE_ALGORITHM(rich_decoding::rich_decoding_t) +constexpr unsigned MaxConnectionsPerTel40 = Allen::Rich::Decoding::Tel40CableMapping::MaxConnectionsPerTel40; + __device__ unsigned rich_calculate_number_of_hits_in_raw_bank( const Allen::RawBank& bank, - const Rich::Future::DAQ::Allen::Tel40CableMapping* cable_mapping, - const Rich::Future::DAQ::Allen::PDMDBDecodeMapping* pdmdb_mapping) + const Allen::Rich::Decoding::Tel40CableMapping* cable_mapping, + const Allen::Rich::Decoding::PDMDBDecodeMapping* pdmdb_mapping, + Allen::Rich::Decoding::PackedFrameSizes::IntType* connSizes) { - unsigned number_of_hits = 0; - std::array< - Rich::Future::DAQ::PackedFrameSizes::IntType, - Rich::Future::DAQ::Allen::Tel40CableMapping::MaxConnectionsPerTel40> - connSizes {}; - auto tel40ID = bank.source_id; - const auto& connMeta = cable_mapping->tel40Meta(tel40ID); - const auto nSizeWords = connMeta.nActiveLinks; - const auto nPackedSizeW = (nSizeWords / 2) + (nSizeWords % 2); - - auto dataW = bank.data; - auto bankEnd = bank.data + bank.size; - auto iPayloadWord = 0u; - auto iWord = 0u; - - for (; iWord < nPackedSizeW && dataW != bankEnd; ++dataW, ++iWord, iPayloadWord += 2) { - // Extract the sizes from the packed word - const Rich::Future::DAQ::PackedFrameSizes sizes(*dataW); - // extract sizes for each packed value - connSizes[iPayloadWord] = sizes.size1(); - connSizes[iPayloadWord + 1] = sizes.size0(); + const auto connMeta = cable_mapping->tel40Meta(tel40ID); + const unsigned nActiveLinks = __popc(connMeta.validLinkMask); + const unsigned nPackedSizeW = (nActiveLinks + 1) / 2; + + auto threadRowDataW = bank.data; + auto threadRowBankEnd = threadRowDataW + bank.size; + + // Extract sizes using improved memory coalescing pattern + for (unsigned i = threadIdx.x; i < nPackedSizeW && threadRowDataW + i < threadRowBankEnd; i += blockDim.x) { + const Allen::Rich::Decoding::PackedFrameSizes sizes(threadRowDataW[i]); + connSizes[2 * i] = sizes.size1(); + connSizes[2 * i + 1] = sizes.size0(); } + __syncwarp(); + + auto dataW = threadRowDataW + nPackedSizeW; const auto& connData = cable_mapping->tel40Data(tel40ID); - if (connMeta.hasInactiveLinks) { - for (unsigned iL = 0; iL < connData.size(); ++iL) { - if (!connData[iL].isActive) { - for (auto i = connData.size() - 1; i > iL; --i) { - connSizes[i] = connSizes[i - 1]; - } - connSizes[iL] = 0; - } + + unsigned localHits = 0; // Hits for this specific section of the bank (thread row) + + const unsigned nWords = bank.size - nPackedSizeW; + uint32_t prev_NZS_first_mask = 0; + uint32_t prev_NZS_second_mask = 0; + unsigned iLink = 0; + unsigned sum = 0; + unsigned linkSize = connSizes[0]; + for (unsigned iW = threadIdx.x; iW < nWords; iW += blockDim.x) { + const unsigned loop_mask = ((iW + RICH_DECODING_BLOCK_DIM_X) > nWords) ? + ~(0xFFFFFFFF << (nWords - (iW & ~(RICH_DECODING_BLOCK_DIM_X - 1)))) : + 0xFFFFFFFF; + + uint8_t word = dataW[iW]; + bool isNZS = (word & 0x80) != 0; + + // find which link we are in and the start of the link word range + while (sum + linkSize <= iW) { + sum += linkSize; + iLink++; + linkSize = connSizes[iLink]; } - } - // finally loop over payload words and decode hits - // note iterator starts from where the above header loop ended... - unsigned iLink {0}; - while (dataW != bankEnd && iLink < connSizes.size()) { - // Do we have any words to decode for this link - if (connSizes[iLink] > 0) { - // Get the Tel40 Data for this connection - const auto& cData = connData[iLink]; - - // get the PDMDB data - const auto& frameData = pdmdb_mapping->getFrameData(cData); - - // Loop over the words for this link - uint16_t iW = 0; - while (iW < connSizes[iLink] && dataW != bankEnd) { - - // check MSB for this word - const auto isNZS = (0x80 & *dataW) != 0; - - if (!isNZS) { - // ZS decoding... word is bit index - - // load the anode data for this bit - if ((unsigned) (*dataW) < frameData.size()) { - const auto& aData = frameData[*dataW]; - - // Data 'could' be invalid, e.g. radiation-induced-upsets - // so cannot make this a hard error - if (aData.ec != -1 && aData.pmtInEC != -1 && aData.anode != -1) { - number_of_hits++; - } - } - - // move to next word - ++iW; - ++dataW; - } - else { - // NZS decoding... - - // which half of the payload are we in ? - const bool firstHalf = (0 == iW && connSizes[iLink] > 5); - // Number of words to decode depends on which half of the payload we are in - const auto nNZSwords = (firstHalf ? 6 : 5); - // bit offset per half - const auto halfBitOffset = (firstHalf ? 39 : 0); - - // look forward last NZS word and read backwards to match frame bit order - for (auto iNZS = nNZSwords - 1; iNZS >= 0; --iNZS) { - - // read the NZS word - auto nzsW = *(dataW + iNZS); - // if word zero clear MSB as this is the NZS flag - if (0 == iNZS) { - nzsW &= 0x7F; - } - - // does this word hold any active bits ? - if (nzsW > 0) { - - // Bit offset for this word - const auto bitOffset = halfBitOffset + (8 * (nNZSwords - 1 - iNZS)); - - // word has data so loop over bits to extract - for (auto iLB = 0; iLB < 8; ++iLB) { - // is bit on ? - // if ( isBitOn( nzsW, iLB ) ) { - if ((nzsW & (1 << iLB)) != 0) { - - // form frame bit value - const auto bit = iLB + bitOffset; - - // load the anode data for this bit - if ((size_t)(bit) < frameData.size()) { - const auto& aData = frameData[bit]; - - // Data 'could' be invalid, e.g. radiation-induced-upsets - // so cannot make this a hard error - if (aData.ec != -1 && aData.pmtInEC != -1 && aData.anode != -1) { - number_of_hits++; - } - } - - } // bit is on - } // loop over word bits - - } // word has any data - - } // loop over all NZS words - - // Finally skip the read NZS words - iW += nNZSwords; - dataW += nNZSwords; + // handle first half (6bytes) + bool is_first_half_head = sum == iW && isNZS && linkSize > 5; + uint64_t first_half_mask = __ballot_sync(loop_mask, is_first_half_head); + first_half_mask *= (1 << 6) - 1; // propagate nzs info to the next 6 bits + first_half_mask |= prev_NZS_first_mask; + prev_NZS_first_mask = first_half_mask >> RICH_DECODING_BLOCK_DIM_X; + + bool isFirstHalfNZS = first_half_mask & (1 << threadIdx.x); + + // handle second half (5bytes) + bool is_second_half_head = (sum + linkSize - 5) == iW && !isFirstHalfNZS && isNZS; + uint64_t second_half_mask = __ballot_sync(loop_mask, is_second_half_head); + second_half_mask *= (1 << 5) - 1; // propagate nzs info to the next 5 bits + second_half_mask |= prev_NZS_second_mask; + prev_NZS_second_mask = second_half_mask >> RICH_DECODING_BLOCK_DIM_X; + + // merge propagated masks and extract true isNZS + uint32_t NZS_mask = first_half_mask | second_half_mask; + isNZS = NZS_mask & (1 << threadIdx.x); + + // get frame mask + const auto& cData = connData[__fns(connMeta.validLinkMask, 0, iLink + 1)]; + const auto& frameMask = pdmdb_mapping->getFrameValidMask(cData); + + // decode the word + if (!isNZS) { // ZS decoding + if (static_cast(word) < 86) { + auto b = word + (word >= 39 ? 1 : 0); + if (frameMask[b / 8] & (1 << (b % 8))) { + localHits++; } } + } + else { // NZS decoding + const unsigned iNZS = isFirstHalfNZS ? (iW - sum) : (iW - (sum + linkSize - 5)); + const unsigned nNZSwords = isFirstHalfNZS ? 6 : 5; + const unsigned halfBitOffset = isFirstHalfNZS ? 40 : 0; + const unsigned bit = halfBitOffset + (8 * (nNZSwords - 1 - iNZS)); + localHits += __popc(word & frameMask[bit / 8]); + } + } - } // no data for this link, so just move on - - // move to next Tel40 link - ++iLink; - } // data word loop - - return number_of_hits; + return Allen::warp::reduce_add_warp_sync(localHits); // Return the total hit count for this bank } template __global__ void rich_calculate_number_of_hits( rich_decoding::Parameters parameters, const unsigned event_start, - const Rich::Future::DAQ::Allen::Tel40CableMapping* cable_mapping, - const Rich::Future::DAQ::Allen::PDMDBDecodeMapping* pdmdb_mapping) + const Allen::Rich::Decoding::Tel40CableMapping* cable_mapping, + const Allen::Rich::Decoding::PDMDBDecodeMapping* pdmdb_mapping) { const auto event_number = parameters.dev_event_list[blockIdx.x]; + const auto number_of_events = parameters.dev_event_list.size(); - // Read raw event const auto raw_event = Allen::RawEvent {parameters.dev_rich_raw_input, parameters.dev_rich_raw_input_offsets, parameters.dev_rich_raw_input_sizes, parameters.dev_rich_raw_input_types, event_number + event_start}; - for (unsigned bank_number = threadIdx.x; bank_number < raw_event.number_of_raw_banks; bank_number += blockDim.x) { + __shared__ Allen::Rich::Decoding::PackedFrameSizes::IntType + connSizes_shared[MaxConnectionsPerTel40 * RICH_DECODING_BLOCK_DIM_Y]; + Allen::Rich::Decoding::PackedFrameSizes::IntType* connSizes = connSizes_shared + threadIdx.y * MaxConnectionsPerTel40; + + unsigned hitCountPannel0 = 0; + unsigned hitCountPannel1 = 0; + for (unsigned bank_number = threadIdx.y; bank_number < raw_event.number_of_raw_banks; bank_number += blockDim.y) { const auto bank = raw_event.raw_bank(bank_number); - const auto number_of_hits_in_raw_bank = - rich_calculate_number_of_hits_in_raw_bank(bank, cable_mapping, pdmdb_mapping); - if (number_of_hits_in_raw_bank > 0) { - atomicAdd(parameters.dev_rich_hit_offsets + event_number, number_of_hits_in_raw_bank); + const auto hitCount = rich_calculate_number_of_hits_in_raw_bank(bank, cable_mapping, pdmdb_mapping, connSizes); + const auto side = (bank.source_id >> 10) & 0x1; + if (side == 0) { + hitCountPannel0 += hitCount; + } + else { + hitCountPannel1 += hitCount; } } + if (threadIdx.x == 0) { + atomicAdd(parameters.dev_rich_hit_offsets + event_number, hitCountPannel0); + atomicAdd(parameters.dev_rich_hit_offsets + event_number + number_of_events, hitCountPannel1); + } } +template +__global__ void rich_calculate_number_of_hits_StreamIDs( + rich_decoding::Parameters parameters, + const unsigned event_start) +{ + const auto event_number = parameters.dev_event_list[blockIdx.x]; + const auto number_of_events = parameters.dev_event_list.size(); + + const auto raw_event = Allen::RawEvent {parameters.dev_rich_raw_input, + parameters.dev_rich_raw_input_offsets, + parameters.dev_rich_raw_input_sizes, + parameters.dev_rich_raw_input_types, + event_number + event_start}; + + for (unsigned bank_number = threadIdx.y * blockDim.x + threadIdx.x; bank_number < raw_event.number_of_raw_banks; + bank_number += blockDim.x * blockDim.y) { + const auto bank = raw_event.raw_bank(bank_number); + const auto hitCount = bank.size / sizeof(word_t); + const auto side = (bank.source_id >> 10) & 0x1; + if (side == 0) { + atomicAdd(parameters.dev_rich_hit_offsets + event_number, hitCount); + } + else { + atomicAdd(parameters.dev_rich_hit_offsets + event_number + number_of_events, hitCount); + } + } +} + +// ---------------------------------------------------------------------- + __device__ void rich_decode_bank( const Allen::RawBank& bank, - const Rich::Future::DAQ::Allen::Tel40CableMapping* cable_mapping, - const Rich::Future::DAQ::Allen::PDMDBDecodeMapping* pdmdb_mapping, + const Allen::Rich::Decoding::Tel40CableMapping* cable_mapping, + const Allen::Rich::Decoding::PDMDBDecodeMapping* pdmdb_mapping, + Allen::Rich::Decoding::PackedFrameSizes::IntType* connSizes, unsigned* event_inserted_hits, - Allen::RichSmartID* event_smart_ids) + Allen::Rich::Decoding::SmartID* event_smart_ids) { - std::array< - Rich::Future::DAQ::PackedFrameSizes::IntType, - Rich::Future::DAQ::Allen::Tel40CableMapping::MaxConnectionsPerTel40> - connSizes {}; - auto tel40ID = bank.source_id; - const auto& connMeta = cable_mapping->tel40Meta(tel40ID); - const auto nSizeWords = connMeta.nActiveLinks; - const auto nPackedSizeW = (nSizeWords / 2) + (nSizeWords % 2); - - auto dataW = bank.data; - auto bankEnd = bank.data + bank.size; - auto iPayloadWord = 0u; - auto iWord = 0u; - - for (; iWord < nPackedSizeW && dataW != bankEnd; ++dataW, ++iWord, iPayloadWord += 2) { - // Extract the sizes from the packed word - const Rich::Future::DAQ::PackedFrameSizes sizes(*dataW); - // extract sizes for each packed value - connSizes[iPayloadWord] = sizes.size1(); - connSizes[iPayloadWord + 1] = sizes.size0(); + const auto connMeta = cable_mapping->tel40Meta(tel40ID); + const unsigned nActiveLinks = __popc(connMeta.validLinkMask); + const unsigned nPackedSizeW = (nActiveLinks + 1) / 2; + + auto threadRowDataW = bank.data; + + for (unsigned i = threadIdx.x; i < nPackedSizeW; i += blockDim.x) { + const Allen::Rich::Decoding::PackedFrameSizes sizes(threadRowDataW[i]); + connSizes[2 * i] = sizes.size1(); + connSizes[2 * i + 1] = sizes.size0(); } + __syncwarp(); + + auto dataW = threadRowDataW; + dataW += nPackedSizeW; const auto& connData = cable_mapping->tel40Data(tel40ID); - if (connMeta.hasInactiveLinks) { - for (unsigned iL = 0; iL < connData.size(); ++iL) { - if (!connData[iL].isActive) { - for (auto i = connData.size() - 1; i > iL; --i) { - connSizes[i] = connSizes[i - 1]; - } - connSizes[iL] = 0; - } + + const unsigned nWords = bank.size - nPackedSizeW; + uint32_t prev_NZS_first_mask = 0; + uint32_t prev_NZS_second_mask = 0; + unsigned iLink = 0; + unsigned sum = 0; + unsigned linkSize = connSizes[0]; + for (unsigned iW = threadIdx.x; iW < nWords; iW += blockDim.x) { + const unsigned loop_mask = ((iW + RICH_DECODING_BLOCK_DIM_X) > nWords) ? + ~(0xFFFFFFFF << (nWords - (iW & ~(RICH_DECODING_BLOCK_DIM_X - 1)))) : + 0xFFFFFFFF; + + uint8_t word = dataW[iW]; + bool isNZS = (word & 0x80) != 0; + + // find which link we are in and the start of the link word range + while (sum + linkSize <= iW) { + sum += linkSize; + iLink++; + linkSize = connSizes[iLink]; } - } - // finally loop over payload words and decode hits - // note iterator starts from where the above header loop ended... - unsigned iLink {0}; - while (dataW != bankEnd && iLink < connSizes.size()) { - // Do we have any words to decode for this link - if (connSizes[iLink] > 0) { - // Get the Tel40 Data for this connection - const auto& cData = connData[iLink]; - - // get the PDMDB data - const auto& frameData = pdmdb_mapping->getFrameData(cData); - - // Loop over the words for this link - uint16_t iW = 0; - while (iW < connSizes[iLink] && dataW != bankEnd) { - - // check MSB for this word - const auto isNZS = (0x80 & *dataW) != 0; - - if (!isNZS) { - // ZS decoding... word is bit index - - // load the anode data for this bit - if ((unsigned) (*dataW) < frameData.size()) { - - const auto& aData = frameData[*dataW]; - - // Data 'could' be invalid, e.g. radiation-induced-upsets - // so cannot make this a hard error - if (aData.ec != -1 && aData.pmtInEC != -1 && aData.anode != -1) { - // make a smart ID - auto hitID = Allen::RichSmartID {cData.smartID}; // sets RICH, side, module and PMT type - - // Add the PMT and pixel info - const auto nInMod = (aData.ec * (0 != ((hitID.key() >> 27) & 0x1) ? 1 : 4)) + aData.pmtInEC; - hitID.setData( - cData.moduleNum, - Allen::RichSmartID::ShiftPDMod, - Allen::RichSmartID::MaskPDMod, - Allen::RichSmartID::MaskPDIsSet); - hitID.setData(nInMod, Allen::RichSmartID::ShiftPDNumInMod, Allen::RichSmartID::MaskPDNumInMod); - - const auto row = aData.anode / 8; - const auto col = 8 - 1 - (aData.anode % 8); - hitID.setData( - row, - Allen::RichSmartID::ShiftPixelRow, - Allen::RichSmartID::MaskPixelRow, - Allen::RichSmartID::MaskPixelRowIsSet); - hitID.setData( - col, - Allen::RichSmartID::ShiftPixelCol, - Allen::RichSmartID::MaskPixelCol, - Allen::RichSmartID::MaskPixelColIsSet); - - const auto insert_index = atomicAdd(event_inserted_hits, 1); - event_smart_ids[insert_index] = hitID; - } - } - - // move to next word - ++iW; - ++dataW; + // handle first half (6bytes) + bool is_first_half_head = sum == iW && isNZS && linkSize > 5; + uint64_t first_half_mask = __ballot_sync(loop_mask, is_first_half_head); + first_half_mask *= (1 << 6) - 1; // propagate nzs info to the next 6 bits + first_half_mask |= prev_NZS_first_mask; + prev_NZS_first_mask = first_half_mask >> RICH_DECODING_BLOCK_DIM_X; + + bool isFirstHalfNZS = first_half_mask & (1 << threadIdx.x); + + // handle second half (5bytes) + bool is_second_half_head = (sum + linkSize - 5) == iW && !isFirstHalfNZS && isNZS; + uint64_t second_half_mask = __ballot_sync(loop_mask, is_second_half_head); + second_half_mask *= (1 << 5) - 1; // propagate nzs info to the next 5 bits + second_half_mask |= prev_NZS_second_mask; + prev_NZS_second_mask = second_half_mask >> RICH_DECODING_BLOCK_DIM_X; + + // merge propagated masks and extract true isNZS + uint32_t NZS_mask = first_half_mask | second_half_mask; + isNZS = NZS_mask & (1 << threadIdx.x); + + // get frame mask + const auto& cData = connData[__fns(connMeta.validLinkMask, 0, iLink + 1)]; + const auto& frameMask = pdmdb_mapping->getFrameValidMask(cData); + const auto& frameData = pdmdb_mapping->getFrameData(cData); + + // decode the word + if (!isNZS) { // ZS decoding + if (static_cast(word) < frameData.size()) { + auto b = word + (word >= 39 ? 1 : 0); + if (frameMask[b / 8] & (1 << (b % 8))) { + const auto& aData = frameData[word]; + + // make a smart ID + auto hitID = Allen::Rich::Decoding::SmartID {cData.smartID}; // sets RICH, side, module and PMT type + + // Add the PMT and pixel info + const auto nInMod = (aData.ec * (0 != ((hitID.key() >> 27) & 0x1) ? 1 : 4)) + aData.pmtInEC; + hitID.setData( + nInMod, Allen::Rich::Decoding::SmartID::ShiftPDNumInMod, Allen::Rich::Decoding::SmartID::MaskPDNumInMod); + hitID.setData( + aData.anode, + Allen::Rich::Decoding::SmartID::ShiftPixelCol, + Allen::Rich::Decoding::SmartID::MaskPixelCol | Allen::Rich::Decoding::SmartID::MaskPixelRow, + Allen::Rich::Decoding::SmartID::MaskPixelColIsSet | Allen::Rich::Decoding::SmartID::MaskPixelRowIsSet); + + const auto insert_index = atomicAdd(event_inserted_hits, 1); + event_smart_ids[insert_index] = hitID; } - else { - // NZS decoding... - - // which half of the payload are we in ? - const bool firstHalf = (0 == iW && connSizes[iLink] > 5); - // Number of words to decode depends on which half of the payload we are in - const auto nNZSwords = (firstHalf ? 6 : 5); - // bit offset per half - const auto halfBitOffset = (firstHalf ? 39 : 0); - - // look forward last NZS word and read backwards to match frame bit order - for (auto iNZS = nNZSwords - 1; iNZS >= 0; --iNZS) { - - // read the NZS word - auto nzsW = *(dataW + iNZS); - // if word zero clear MSB as this is the NZS flag - if (0 == iNZS) { - nzsW &= 0x7F; - } - - // does this word hold any active bits ? - if (nzsW > 0) { - - // Bit offset for this word - const auto bitOffset = halfBitOffset + (8 * (nNZSwords - 1 - iNZS)); - - // word has data so loop over bits to extract - for (auto iLB = 0; iLB < 8; ++iLB) { - // is bit on ? - // if ( isBitOn( nzsW, iLB ) ) { - if ((nzsW & (1 << iLB)) != 0) { - - // form frame bit value - const auto bit = iLB + bitOffset; - - // load the anode data for this bit - if ((size_t)(bit) < frameData.size()) { - const auto& aData = frameData[bit]; - // Data 'could' be invalid, e.g. radiation-induced-upsets - // so cannot make this a hard error - if (aData.ec != -1 && aData.pmtInEC != -1 && aData.anode != -1) { - - // make a smart ID - auto hitID = Allen::RichSmartID {cData.smartID}; // sets RICH, side, module and PMT type - - // Add the PMT and pixel info - const auto nInMod = (aData.ec * (0 != ((hitID.key() >> 27) & 0x1) ? 1 : 4)) + aData.pmtInEC; - hitID.setData( - cData.moduleNum, - Allen::RichSmartID::ShiftPDMod, - Allen::RichSmartID::MaskPDMod, - Allen::RichSmartID::MaskPDIsSet); - hitID.setData(nInMod, Allen::RichSmartID::ShiftPDNumInMod, Allen::RichSmartID::MaskPDNumInMod); - - const auto row = aData.anode / 8; - const auto col = 8 - 1 - (aData.anode % 8); - hitID.setData( - row, - Allen::RichSmartID::ShiftPixelRow, - Allen::RichSmartID::MaskPixelRow, - Allen::RichSmartID::MaskPixelRowIsSet); - hitID.setData( - col, - Allen::RichSmartID::ShiftPixelCol, - Allen::RichSmartID::MaskPixelCol, - Allen::RichSmartID::MaskPixelColIsSet); - - const auto insert_index = atomicAdd(event_inserted_hits, 1); - event_smart_ids[insert_index] = hitID; - } - } - - } // bit is on - } // loop over word bits - - } // word has any data - - } // loop over all NZS words - - // Finally skip the read NZS words - iW += nNZSwords; - dataW += nNZSwords; + } + } + else { // NZS decoding + const unsigned iNZS = isFirstHalfNZS ? (iW - sum) : (iW - (sum + linkSize - 5)); + const unsigned nNZSwords = isFirstHalfNZS ? 6 : 5; + const unsigned halfBitOffset = isFirstHalfNZS ? 39 : 0; + const unsigned bitOffset = halfBitOffset + (8 * (nNZSwords - 1 - iNZS)); + + auto b = bitOffset + (bitOffset >= 39 ? 1 : 0); + word &= frameMask[b / 8]; // mask invalid hits + + // Reserve enough space to push hits + auto insert_index = atomicAdd(event_inserted_hits, __popc(word)); + + // Use find first set intrinsic to quickly iterate over valid hits + while (word != 0) { + unsigned iLB = __ffs(word) - 1; // __ffs returns 1-based index + const auto bit = iLB + bitOffset; + + if (bit < frameData.size()) { + const auto& aData = frameData[bit]; + // make a smart ID + auto hitID = Allen::Rich::Decoding::SmartID {cData.smartID}; // sets RICH, side, module and PMT type + + // Add the PMT and pixel info + const auto nInMod = (aData.ec * (0 != ((hitID.key() >> 27) & 0x1) ? 1 : 4)) + aData.pmtInEC; + hitID.setData( + nInMod, Allen::Rich::Decoding::SmartID::ShiftPDNumInMod, Allen::Rich::Decoding::SmartID::MaskPDNumInMod); + hitID.setData( + aData.anode, + Allen::Rich::Decoding::SmartID::ShiftPixelCol, + Allen::Rich::Decoding::SmartID::MaskPixelCol | Allen::Rich::Decoding::SmartID::MaskPixelRow, + Allen::Rich::Decoding::SmartID::MaskPixelColIsSet | Allen::Rich::Decoding::SmartID::MaskPixelRowIsSet); + + event_smart_ids[insert_index++] = hitID; } + word &= ~(1u << iLB); // Clear first set bit } - - } // no data for this link, so just move on - - // move to next Tel40 link - ++iLink; - } // data word loop + } + } + __syncwarp(); } template __global__ void rich_decoding_kernel( rich_decoding::Parameters parameters, const unsigned event_start, - const Rich::Future::DAQ::Allen::Tel40CableMapping* cable_mapping, - const Rich::Future::DAQ::Allen::PDMDBDecodeMapping* pdmdb_mapping, - unsigned* dev_rich_number_of_inserted_hits) + const Allen::Rich::Decoding::Tel40CableMapping* cable_mapping, + const Allen::Rich::Decoding::PDMDBDecodeMapping* pdmdb_mapping) { const auto event_number = parameters.dev_event_list[blockIdx.x]; - - auto event_inserted_hits = dev_rich_number_of_inserted_hits + event_number; - auto event_smart_ids = parameters.dev_smart_ids + parameters.dev_rich_hit_offsets[event_number]; + const auto number_of_events = parameters.dev_event_list.size(); // Read raw event const auto raw_event = Allen::RawEvent {parameters.dev_rich_raw_input, @@ -406,9 +338,76 @@ __global__ void rich_decoding_kernel( parameters.dev_rich_raw_input_types, event_number + event_start}; - for (unsigned bank_number = threadIdx.x; bank_number < raw_event.number_of_raw_banks; bank_number += blockDim.x) { + __shared__ Allen::Rich::Decoding::PackedFrameSizes::IntType + connSizes_shared[MaxConnectionsPerTel40 * RICH_DECODING_BLOCK_DIM_Y]; + Allen::Rich::Decoding::PackedFrameSizes::IntType* connSizes = connSizes_shared + threadIdx.y * MaxConnectionsPerTel40; + + auto event_smart_ids0 = parameters.dev_smart_ids + parameters.dev_rich_hit_offsets[event_number]; + auto event_smart_ids1 = parameters.dev_smart_ids + parameters.dev_rich_hit_offsets[event_number + number_of_events]; + + __shared__ unsigned event_inserted_hits[2]; + if (threadIdx.y == 0) { + for (unsigned i = threadIdx.x; i < 2; i += blockDim.x) { + event_inserted_hits[i] = 0; + } + } + __syncthreads(); + + for (unsigned bank_number = threadIdx.y; bank_number < raw_event.number_of_raw_banks; bank_number += blockDim.y) { + const auto bank = raw_event.raw_bank(bank_number); + const auto side = (bank.source_id >> 10) & 0x1; + rich_decode_bank( + bank, + cable_mapping, + pdmdb_mapping, + connSizes, + &event_inserted_hits[side], + side == 0 ? event_smart_ids0 : event_smart_ids1); + } +} + +template +__global__ void rich_decoding_kernel_StreamIDs(rich_decoding::Parameters parameters, const unsigned event_start) +{ + const auto event_number = parameters.dev_event_list[blockIdx.x]; + const auto number_of_events = parameters.dev_event_list.size(); + + const auto raw_event = Allen::RawEvent {parameters.dev_rich_raw_input, + parameters.dev_rich_raw_input_offsets, + parameters.dev_rich_raw_input_sizes, + parameters.dev_rich_raw_input_types, + event_number + event_start}; + + auto event_smart_ids0 = parameters.dev_smart_ids + parameters.dev_rich_hit_offsets[event_number]; + auto event_smart_ids1 = parameters.dev_smart_ids + parameters.dev_rich_hit_offsets[event_number + number_of_events]; + + __shared__ unsigned event_inserted_hits[2]; + if (threadIdx.y == 0) { + for (unsigned i = threadIdx.x; i < 2; i += blockDim.x) { + event_inserted_hits[i] = 0; + } + } + __syncthreads(); + + for (unsigned bank_number = threadIdx.y; bank_number < raw_event.number_of_raw_banks; bank_number += blockDim.y) { const auto bank = raw_event.raw_bank(bank_number); - rich_decode_bank(bank, cable_mapping, pdmdb_mapping, event_inserted_hits, event_smart_ids); + const auto side = (bank.source_id >> 10) & 0x1; + const word_t* bank_data = reinterpret_cast(bank.data); + + unsigned* bank_inserted_hits = &event_inserted_hits[side]; + Allen::Rich::Decoding::SmartID* event_smart_ids = side == 0 ? event_smart_ids0 : event_smart_ids1; + + for (unsigned iW = threadIdx.x; iW < bank.size / sizeof(word_t); iW += blockDim.x) { + auto hitID = Allen::Rich::Decoding::SmartID {bank_data[iW]}; +#ifdef USE_DD4HEP + if (hitID.rich() == 1) { + const auto pdMod = hitID.pdMod() + (hitID.panel() == 0 ? 6 : 18); + hitID.setData(pdMod, Allen::Rich::Decoding::SmartID::ShiftPDMod, Allen::Rich::Decoding::SmartID::MaskPDMod); + } +#endif + const auto insert_index = atomicAdd(bank_inserted_hits, 1); + event_smart_ids[insert_index] = hitID; + } } } @@ -417,7 +416,7 @@ void rich_decoding::rich_decoding_t::set_arguments_size( const RuntimeOptions&, const Constants&) const { - set_size(arguments, size(arguments) + 1); + set_size(arguments, size(arguments) * 2 + 1); set_size(arguments, 1); } @@ -427,38 +426,62 @@ void rich_decoding::rich_decoding_t::operator()( const Constants& constants, const Allen::Context& context) const { + Allen::memset_async(arguments, 0, context); + const auto bank_version = first(arguments); - if (bank_version != 10) { + if (bank_version == 10) { /// 'Production' PMT version (Used for real data) + Allen::Rich::Decoding::Tel40CableMapping* cable_mapping = constants.dev_rich_cable_mapping; + Allen::Rich::Decoding::PDMDBDecodeMapping* pdmdb_mapping = constants.dev_rich_pdmdb_mapping; + + // Calculate number of hits into dev_rich_hit_offsets_t + global_function( + runtime_options.mep_layout ? rich_calculate_number_of_hits : rich_calculate_number_of_hits)( + dim3(size(arguments)), m_block_dim, context)( + arguments, std::get<0>(runtime_options.event_interval), cable_mapping, pdmdb_mapping); + PrefixSum::prefix_sum(*this, arguments, context); + + // Decode RICH hits + resize(arguments, first(arguments)); + + global_function(runtime_options.mep_layout ? rich_decoding_kernel : rich_decoding_kernel)( + dim3(size(arguments)), m_block_dim, context)( + arguments, std::get<0>(runtime_options.event_interval), cable_mapping, pdmdb_mapping); + } + else if (bank_version == 3) { /// Simple version. Just streams SmartIDs, excluding the higher time bits. (Used for MC + /// data) + global_function( + runtime_options.mep_layout ? rich_calculate_number_of_hits_StreamIDs : + rich_calculate_number_of_hits_StreamIDs)( + dim3(size(arguments)), m_block_dim, context)( + arguments, std::get<0>(runtime_options.event_interval)); + PrefixSum::prefix_sum(*this, arguments, context); + resize(arguments, first(arguments)); + + global_function( + runtime_options.mep_layout ? + rich_decoding_kernel_StreamIDs : + rich_decoding_kernel_StreamIDs)(dim3(size(arguments)), m_block_dim, context)( + arguments, std::get<0>(runtime_options.event_interval)); + } + else if (bank_version == 4) { /// Updated streaming for 64 bit SmartIDs, with adc time information (Used for MC data) + global_function( + runtime_options.mep_layout ? rich_calculate_number_of_hits_StreamIDs : + rich_calculate_number_of_hits_StreamIDs)( + dim3(size(arguments)), m_block_dim, context)( + arguments, std::get<0>(runtime_options.event_interval)); + PrefixSum::prefix_sum(*this, arguments, context); + resize(arguments, first(arguments)); + + global_function( + runtime_options.mep_layout ? + rich_decoding_kernel_StreamIDs : + rich_decoding_kernel_StreamIDs)(dim3(size(arguments)), m_block_dim, context)( + arguments, std::get<0>(runtime_options.event_interval)); + } + else { throw StrException("Rich bank version not supported (" + std::to_string(bank_version) + ")"); } - auto cable_mapping = reinterpret_cast(constants.dev_rich_cable_mapping); - auto pdmdb_mapping = - reinterpret_cast(constants.dev_rich_pdmdb_mapping); - - // Calculate number of hits into dev_rich_hit_offsets_t - Allen::memset_async(arguments, 0, context); - global_function( - runtime_options.mep_layout ? rich_calculate_number_of_hits : rich_calculate_number_of_hits)( - dim3(size(arguments)), dim3(m_block_dim_x), context)( - arguments, std::get<0>(runtime_options.event_interval), cable_mapping, pdmdb_mapping); - - PrefixSum::prefix_sum(*this, arguments, context); - - // Decode RICH hits - auto dev_rich_number_of_inserted_hits = make_device_buffer(arguments, size(arguments)); - Allen::memset_async( - dev_rich_number_of_inserted_hits.data(), 0, dev_rich_number_of_inserted_hits.size_bytes(), context); - resize(arguments, first(arguments)); - - global_function(runtime_options.mep_layout ? rich_decoding_kernel : rich_decoding_kernel)( - dim3(size(arguments)), dim3(m_block_dim_x), context)( - arguments, - std::get<0>(runtime_options.event_interval), - cable_mapping, - pdmdb_mapping, - dev_rich_number_of_inserted_hits.data()); - if (m_verbosity >= logger::debug) { // Print output print(arguments); diff --git a/device/rich/decoding/src/RichMakePixels.cu b/device/rich/decoding/src/RichMakePixels.cu new file mode 100644 index 0000000000000000000000000000000000000000..53c92eded351c4ecef64b0e6d27fcb6dc5a20eb1 --- /dev/null +++ b/device/rich/decoding/src/RichMakePixels.cu @@ -0,0 +1,115 @@ +/*****************************************************************************\ + * (c) Copyright 2018-2020 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 "RichDefinitions.cuh" +#include "Rich.cuh" +#include +#include +#include + +INSTANTIATE_ALGORITHM(rich_make_pixels::rich_make_pixels_t); + +/// Create a Pixel object from a given SmartID +template +__device__ Allen::Rich::PixelReco::Pixel make_pixel( + const Allen::Rich::Detector::PDPanel& panel, + const std::array& g2panel, + const Allen::Rich::Decoding::SmartID& smartID) +{ + bool smartIDSelected = true; + if constexpr (Allen::Rich::m_enable4D[rich]) { + if (smartID.adcTimeIsSet()) { + const auto hitT = smartID.time(); + if (fabsf(hitT - Allen::Rich::m_avHitTime[rich]) > Allen::Rich::m_timeWindow[rich]) { + smartIDSelected = false; + } + } + } + + Allen::Rich::PixelReco::Pixel pixel; + + // if time checks pass, and a non-Null PD matches the SmartID, save the pix + if (smartIDSelected) { + const auto& pixelPD = panel.dePD(smartID); + + // From panel local coordinates to global coordinates. + const auto gPos = panel.globalDetectionPoint(smartID); + // From global coordinates to RICH detector local coordinates. + const auto lPos = Allen::Rich::transform3DTimesPoint(g2panel, gPos); + + const auto isInner = + (Allen::Rich::m_overrideRegions[rich] ? (fabsf(lPos.x) < float(Allen::Rich::m_innerPixX[rich]) && + fabsf(lPos.y) < float(Allen::Rich::m_innerPixY[rich])) : + !smartID.isLargePMT()); + + // 4D specific properties + float timeWindow = 999999; + if constexpr (Allen::Rich::m_enable4D[rich]) { + timeWindow = isInner ? float(Allen::Rich::m_innerTimeWindow[rich]) : float(Allen::Rich::m_outerTimeWindow[rich]); + } + pixel = Allen::Rich::PixelReco::Pixel(gPos, lPos, smartID, pixelPD.effectivePixelArea(), timeWindow); + if constexpr (Allen::Rich::m_overrideRegions[rich]) { + pixel.overrideRegions(isInner); + } + pixel.setIsNull(pixelPD.getIsNull()); + } + return pixel; +} + +/// Kernel function to iterate over Events, SmartIDs, and create Pixels +template +__global__ void rich_make_pixels_kernel( + rich_make_pixels::Parameters parameters, + const Allen::Rich::RichDetector* deRich) +{ + const auto side = blockIdx.y; // 0 or 1 for RICH1 top/bottom or RICH2 left/right + const auto number_of_events = parameters.dev_event_list.size(); + auto offset = parameters.dev_rich_hit_offsets[side * number_of_events]; + auto pixelsSize = parameters.dev_rich_hit_offsets[(side + 1) * number_of_events] - offset; + + const Allen::Rich::Decoding::SmartID* ids = parameters.dev_smart_ids + offset; + Allen::Rich::PixelReco::Pixel* pixels = parameters.dev_rich_pixels + offset; + + const Allen::Rich::Detector::PDPanel& panel = deRich->pdPanels()[side]; + const auto g2panel = panel.globalToPDPanel(); // cache global to pd panel matrix in registers + + // pixel rec loop, grid-stride loop + uint global_tid = blockIdx.x * blockDim.x + threadIdx.x; + uint stride = blockDim.x * gridDim.x; + for (uint i = global_tid; i < pixelsSize; i += stride) { + pixels[i] = make_pixel(panel, g2panel, ids[i]); + } +} + +/// Function to resize Allen sequence output data structure +void rich_make_pixels::rich_make_pixels_t::set_arguments_size( + ArgumentReferences arguments, + const RuntimeOptions&, + const Constants&) const +{ + set_size(arguments, size(arguments)); +} + +void rich_make_pixels::rich_make_pixels_t::operator()( + const ArgumentReferences& arguments, + const RuntimeOptions&, + const Constants& constants, + const Allen::Context& context) const +{ + // Create RICH object from dump + const auto richValue = m_current_rich; + const Allen::Rich::RichDetector* rich = + richValue == 1 ? constants.dev_rich_1_geometry : constants.dev_rich_2_geometry; + const auto& kernel = richValue == 1 ? rich_make_pixels_kernel<0> : rich_make_pixels_kernel<1>; + + // dont launch too many blocks and instead rely on the grid stride loop to amortize launch overhead + global_function(kernel)(dim3(16, 2), m_block_dim, context)(arguments, rich); +} diff --git a/device/utils/warp/include/WarpIntrinsicsTools.cuh b/device/utils/warp/include/WarpIntrinsicsTools.cuh index 2dbe2460f801cdba92ec76a84723af1a4eb44d5d..7cdae0bcfd6d5c2780ac530a715e68ceb0508d25 100644 --- a/device/utils/warp/include/WarpIntrinsicsTools.cuh +++ b/device/utils/warp/include/WarpIntrinsicsTools.cuh @@ -16,6 +16,24 @@ namespace Allen::warp { #if defined(TARGET_DEVICE_CUDA) + __device__ inline unsigned warp_mask() + { + const unsigned warp_size = min(blockDim.x, 32); + const unsigned warp_offset = ((threadIdx.y * blockDim.x + threadIdx.x) % 32) / warp_size; + return ((1 << warp_size) - 1) << (warp_size * warp_offset); + } + + __device__ inline unsigned reduce_add_warp_sync(unsigned value) + { + const unsigned warp_size = min(blockDim.x, 32); + const unsigned warp_offset = ((threadIdx.y * blockDim.x + threadIdx.x) % 32) / warp_size; + const unsigned mask = ((1 << warp_size) - 1) << (warp_size * warp_offset); + for (int offset = warp_size / 2; offset > 0; offset /= 2) { + value += __shfl_xor_sync(mask, value, offset); + } + return value; + } + __device__ inline unsigned prefix_sum_and_increase_size(const bool process_element, unsigned& size) { const auto mask = __ballot_sync(0xFFFFFFFF, process_element); @@ -115,6 +133,10 @@ namespace Allen::warp { #elif defined(TARGET_DEVICE_CPU) + __device__ inline unsigned warp_mask() { return 1; } + + __device__ inline unsigned reduce_add_warp_sync(unsigned value) { return value; } + __device__ inline unsigned prefix_sum_and_increase_size(const bool process_element, unsigned& size) { const auto address = size; diff --git a/input/allen_geometries/geometry_run3_2025-v00.01/rich_1_geometry.bin b/input/allen_geometries/geometry_run3_2025-v00.01/rich_1_geometry.bin new file mode 100644 index 0000000000000000000000000000000000000000..207b37c979b357378d6cc349dcd006d079ad741e Binary files /dev/null and b/input/allen_geometries/geometry_run3_2025-v00.01/rich_1_geometry.bin differ diff --git a/input/allen_geometries/geometry_run3_2025-v00.01/rich_2_geometry.bin b/input/allen_geometries/geometry_run3_2025-v00.01/rich_2_geometry.bin new file mode 100644 index 0000000000000000000000000000000000000000..df81a3f049d7c54cc7fa29e305ae3d4f66133d1f Binary files /dev/null and b/input/allen_geometries/geometry_run3_2025-v00.01/rich_2_geometry.bin differ diff --git a/integration/non_event_data/include/Consumers.h b/integration/non_event_data/include/Consumers.h index 619d9329ab59b9456f4355d1690e52fb31f62793..81300672f9bbd650cc86f1dbf90128974b54a02e 100644 --- a/integration/non_event_data/include/Consumers.h +++ b/integration/non_event_data/include/Consumers.h @@ -49,6 +49,50 @@ namespace Consumers { std::reference_wrapper m_constants; }; + struct RichPDMDBDecodeMapping final : public Allen::NonEventData::Consumer { + public: + RichPDMDBDecodeMapping(Constants& constants); + + void consume(std::vector const& dev_rich_pdmdb_mapping) override; + + private: + void initialize(const std::vector& data); + std::reference_wrapper m_constants; + }; + + struct RichTel40CableMapping final : public Allen::NonEventData::Consumer { + public: + RichTel40CableMapping(Constants& constants); + + void consume(std::vector const& dev_rich_cable_mapping) override; + + private: + void initialize(const std::vector& data); + std::reference_wrapper m_constants; + }; + + struct Rich1Geometry final : public Allen::NonEventData::Consumer { + public: + Rich1Geometry(Constants& constants); + + void consume(std::vector const& dev_rich_1_geometry) override; + + private: + void initialize(const std::vector& data); + std::reference_wrapper m_constants; + }; + + struct Rich2Geometry final : public Allen::NonEventData::Consumer { + public: + Rich2Geometry(Constants& constants); + + void consume(std::vector const& dev_rich_2_geometry) override; + + private: + void initialize(const std::vector& data); + std::reference_wrapper m_constants; + }; + struct UTBoards final : public Allen::NonEventData::Consumer { public: UTBoards(Constants& constants); diff --git a/integration/non_event_data/src/RichGeometry.cpp b/integration/non_event_data/src/RichGeometry.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6bac0c00339d489b54ecbae02b52e134c314ce55 --- /dev/null +++ b/integration/non_event_data/src/RichGeometry.cpp @@ -0,0 +1,176 @@ +/*****************************************************************************\ + * * (c) Copyright 2018-2020 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 + +#include +#include +#include +#include "RichPDMDBDecodeMapping.cuh" +#include "RichTel40CableMapping.cuh" +#include "Rich.cuh" + +namespace { + using std::string; + using std::to_string; + using std::vector; +} // namespace + +// PDMDB mapping +Consumers::RichPDMDBDecodeMapping::RichPDMDBDecodeMapping(Constants& constants) : m_constants {constants} {} + +void Consumers::RichPDMDBDecodeMapping::initialize(vector const&) +{ + Allen::malloc((void**) &m_constants.get().dev_rich_pdmdb_mapping, sizeof(Allen::Rich::Decoding::PDMDBDecodeMapping)); +} + +void Consumers::RichPDMDBDecodeMapping::consume(vector const& data) +{ + auto& dev_rich_pdmdb_mapping = m_constants.get().dev_rich_pdmdb_mapping; + if (dev_rich_pdmdb_mapping == nullptr) { + initialize(data); + } + auto& host_rich_pdmdb_mapping = m_constants.get().host_rich_pdmdb_mapping; + host_rich_pdmdb_mapping = data; + host_rich_pdmdb_mapping.resize(sizeof(Allen::Rich::Decoding::PDMDBDecodeMapping)); + + // Patch anode (reverse column) and compute valid masks + Allen::Rich::Decoding::PDMDBDecodeMapping* richPDMDBMapping = + reinterpret_cast(host_rich_pdmdb_mapping.data()); + auto& pdmDataR = richPDMDBMapping->m_pdmDataR; + auto& pdmDataH = richPDMDBMapping->m_pdmDataH; + auto& pdmMaskR = richPDMDBMapping->m_pdmMaskR; + auto& pdmMaskH = richPDMDBMapping->m_pdmMaskH; + + for (unsigned rich = 0; rich < 2; rich++) { + for (unsigned pdmdb = 0; pdmdb < 2; pdmdb++) { + for (unsigned link = 0; link < 6; link++) { + for (unsigned b = 0; b < 11; b++) { + pdmMaskR[rich][pdmdb][link][b] = 0; + if (rich == 0) pdmMaskH[pdmdb][link][b] = 0; + } + for (unsigned bit = 0; bit < 86; bit++) { + auto b = bit + (bit >= 39 ? 1 : 0); + auto& dataR = pdmDataR[rich][pdmdb][link][bit]; + auto& dataH = pdmDataH[pdmdb][link][bit]; + + unsigned validR = dataR.ec != -1 && dataR.pmtInEC != -1 && dataR.anode != -1; + unsigned validH = dataH.ec != -1 && dataH.pmtInEC != -1 && dataH.anode != -1; + + pdmMaskR[rich][pdmdb][link][b / 8] |= validR << (b % 8); + if (rich == 0) pdmMaskH[pdmdb][link][b / 8] |= validH << (b % 8); + + // Patch anode (reverse column) + auto row = dataR.anode / 8; + auto col = 7 - (dataR.anode % 8); + dataR.anode = col + 8 * row; + if (rich == 0) { + row = dataH.anode / 8; + col = 7 - (dataH.anode % 8); + dataH.anode = col + 8 * row; + } + } + } + } + } + + Allen::memcpy( + dev_rich_pdmdb_mapping, + host_rich_pdmdb_mapping.data(), + sizeof(Allen::Rich::Decoding::PDMDBDecodeMapping), + Allen::memcpyHostToDevice); +} + +// Tel40 cable mapping +Consumers::RichTel40CableMapping::RichTel40CableMapping(Constants& constants) : m_constants {constants} {} + +void Consumers::RichTel40CableMapping::initialize(vector const&) +{ + Allen::malloc((void**) &m_constants.get().dev_rich_cable_mapping, sizeof(Allen::Rich::Decoding::Tel40CableMapping)); +} + +void Consumers::RichTel40CableMapping::consume(vector const& data) +{ + auto& dev_rich_cable_mapping = m_constants.get().dev_rich_cable_mapping; + if (dev_rich_cable_mapping == nullptr) { + initialize(data); + } + auto& host_rich_cable_mapping = m_constants.get().host_rich_cable_mapping; + host_rich_cable_mapping = data; + + // patch metadata to hold valid mask, TODO: do this properly in the dumpers + Allen::Rich::Decoding::Tel40CableMapping* richCableMapping = + reinterpret_cast(host_rich_cable_mapping.data()); + auto& metadata = richCableMapping->m_tel40ConnMeta; + auto& sources = richCableMapping->m_tel40ConnData; + for (unsigned i = 0; i < metadata.size(); ++i) { + for (unsigned j = 0; j < metadata[i].size(); ++j) { + for (unsigned k = 0; k < metadata[i][j].size(); ++k) { + uint32_t mask = 0; + for (unsigned l = 0; l < sources[i][j][k].size(); ++l) { + mask |= (sources[i][j][k][l].isActive ? (1 << l) : 0); + + sources[i][j][k][l].smartID.setData( + sources[i][j][k][l].moduleNum, + Allen::Rich::Decoding::SmartID::ShiftPDMod, + Allen::Rich::Decoding::SmartID::MaskPDMod, + Allen::Rich::Decoding::SmartID::MaskPDIsSet); + } + metadata[i][j][k].validLinkMask = mask; + } + } + } + + Allen::memcpy( + dev_rich_cable_mapping, + host_rich_cable_mapping.data(), + sizeof(Allen::Rich::Decoding::Tel40CableMapping), + Allen::memcpyHostToDevice); +} + +// RICH 1 geometry +Consumers::Rich1Geometry::Rich1Geometry(Constants& constants) : m_constants {constants} {} + +void Consumers::Rich1Geometry::initialize(vector const&) +{ + Allen::malloc((void**) &m_constants.get().dev_rich_1_geometry, sizeof(Allen::Rich::RichDetector)); +} + +void Consumers::Rich1Geometry::consume(vector const& data) +{ + auto& dev_rich_1_geometry = m_constants.get().dev_rich_1_geometry; + if (dev_rich_1_geometry == nullptr) { + initialize(data); + } + auto& host_rich_1_geometry = m_constants.get().host_rich_1_geometry; + host_rich_1_geometry = data; + Allen::memcpy( + dev_rich_1_geometry, host_rich_1_geometry.data(), sizeof(Allen::Rich::RichDetector), Allen::memcpyHostToDevice); +} + +// RICH 2 geometry +Consumers::Rich2Geometry::Rich2Geometry(Constants& constants) : m_constants {constants} {} + +void Consumers::Rich2Geometry::initialize(vector const&) +{ + Allen::malloc((void**) &m_constants.get().dev_rich_2_geometry, sizeof(Allen::Rich::RichDetector)); +} + +void Consumers::Rich2Geometry::consume(vector const& data) +{ + auto& dev_rich_2_geometry = m_constants.get().dev_rich_2_geometry; + if (dev_rich_2_geometry == nullptr) { + initialize(data); + } + auto& host_rich_2_geometry = m_constants.get().host_rich_2_geometry; + host_rich_2_geometry = data; + Allen::memcpy( + dev_rich_2_geometry, host_rich_2_geometry.data(), sizeof(Allen::Rich::RichDetector), Allen::memcpyHostToDevice); +} diff --git a/integration/non_event_data/src/Updater.cpp b/integration/non_event_data/src/Updater.cpp index c63b72d4c6b28206bb2b71fcdc3da459cfb0737f..7c37280f129f26839a4cd37c5d2c3aa007e32fd8 100644 --- a/integration/non_event_data/src/Updater.cpp +++ b/integration/non_event_data/src/Updater.cpp @@ -65,7 +65,9 @@ namespace Allen { tuple {NonEventData::MuonGeometry {}, std::string("muon_geometry.bin")}, tuple {NonEventData::MuonLookupTables {}, std::string("muon_tables.bin")}, tuple {NonEventData::RichPDMDBMapping {}, std::string("rich_pdmdbmaps.bin")}, - tuple {NonEventData::RichCableMapping {}, std::string("rich_tel40maps.bin")}}; + tuple {NonEventData::RichCableMapping {}, std::string("rich_tel40maps.bin")}, + tuple {NonEventData::Rich1Geometry {}, std::string("rich_1_geometry.bin")}, + tuple {NonEventData::Rich2Geometry {}, std::string("rich_2_geometry.bin")}}; for_each(producers, [this, &geometry_producer](const auto& p) { using id_t = typename std::remove_reference_t(p))>; diff --git a/main/src/RegisterConsumers.cpp b/main/src/RegisterConsumers.cpp index 4eeee89dfd6b4f8bed2cc80819523194d9f09db6..1e7a6a2e94074fa83226a7144b2e7ce35607b367 100644 --- a/main/src/RegisterConsumers.cpp +++ b/main/src/RegisterConsumers.cpp @@ -77,31 +77,31 @@ void register_consumers( BankTypes::MUON), std::make_tuple( Allen::NonEventData::RichPDMDBMapping {}, - [&constants]() { - return std::make_unique( - constants.host_rich_pdmdb_mapping, constants.dev_rich_pdmdb_mapping); - }, + [&constants]() { return std::make_unique(constants); }, BankTypes::Rich1), std::make_tuple( Allen::NonEventData::RichCableMapping {}, - [&constants]() { - return std::make_unique( - constants.host_rich_cable_mapping, constants.dev_rich_cable_mapping); - }, + [&constants]() { return std::make_unique(constants); }, BankTypes::Rich1), std::make_tuple( Allen::NonEventData::RichPDMDBMapping {}, - [&constants]() { - return std::make_unique( - constants.host_rich_pdmdb_mapping, constants.dev_rich_pdmdb_mapping); - }, + [&constants]() { return std::make_unique(constants); }, + BankTypes::Rich1), + std::make_tuple( + Allen::NonEventData::Rich1Geometry {}, + [&constants]() { return std::make_unique(constants); }, + BankTypes::Rich1), + std::make_tuple( + Allen::NonEventData::RichPDMDBMapping {}, + [&constants]() { return std::make_unique(constants); }, BankTypes::Rich2), std::make_tuple( Allen::NonEventData::RichCableMapping {}, - [&constants]() { - return std::make_unique( - constants.host_rich_cable_mapping, constants.dev_rich_cable_mapping); - }, + [&constants]() { return std::make_unique(constants); }, + BankTypes::Rich2), + std::make_tuple( + Allen::NonEventData::Rich2Geometry {}, + [&constants]() { return std::make_unique(constants); }, BankTypes::Rich2)); const auto unconditional_consumers = std::make_tuple( diff --git a/scripts/ci/test_config.yaml b/scripts/ci/test_config.yaml index a6a85829111ef2ebabfe90b873e5dc96686345a8..cee1a2446173347b7f7a1c03095661f0e45f9124 100644 --- a/scripts/ci/test_config.yaml +++ b/scripts/ci/test_config.yaml @@ -110,6 +110,7 @@ full: - hlt1_pp_forward_then_matching_and_downstream - hlt1_pp_forward_then_matching_with_parkf - hlt1_pp_forward_then_matching_rich #throughput test with RICH decoding -> medium/long term development + - hlt1_pp_forward_then_matching_and_downstream_with_parkf_rich dataset: "MEP_2025_pp_pD2_bu_321834_LHCb_ECEB01_BU_0" # we should keep this one properly updated! geometry: "geometry_run3_2025-v00.01" diff --git a/stream/sequence/include/Constants.cuh b/stream/sequence/include/Constants.cuh index df9da96e33773daf677993955cdf26ed803f4e07..21d74615feab8c56e1f20d1b5b58cd89b7dc20cb 100644 --- a/stream/sequence/include/Constants.cuh +++ b/stream/sequence/include/Constants.cuh @@ -43,11 +43,13 @@ namespace MatchUpstreamMuon { namespace TrackMatchingConsts { struct MagnetParametrization; } -namespace Rich::Future::DAQ::Allen { - class PDMDBDecodeMapping; - class Tel40CableMapping; -} // namespace Rich::Future::DAQ::Allen - +namespace Allen::Rich::Decoding { + struct PDMDBDecodeMapping; + struct Tel40CableMapping; +} // namespace Allen::Rich::Decoding +namespace Allen::Rich { + struct RichDetector; +} // namespace Allen::Rich namespace UT::Constants { struct UTLayerGeometry; } @@ -154,8 +156,13 @@ struct Constants { // Rich std::vector host_rich_pdmdb_mapping; std::vector host_rich_cable_mapping; - char* dev_rich_pdmdb_mapping; - char* dev_rich_cable_mapping; + std::vector host_rich_1_geometry; + std::vector host_rich_2_geometry; + + Allen::Rich::Decoding::PDMDBDecodeMapping* dev_rich_pdmdb_mapping = nullptr; + Allen::Rich::Decoding::Tel40CableMapping* dev_rich_cable_mapping = nullptr; + Allen::Rich::RichDetector* dev_rich_1_geometry = nullptr; + Allen::Rich::RichDetector* dev_rich_2_geometry = nullptr; /** * @brief Reserves and initializes constants.