From a8f075c5cea7d91b450a4a172085d77e9f0657e8 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Tue, 28 Nov 2023 10:33:20 +0000 Subject: [PATCH 01/23] Add new MC checker RichMCDetectorHits --- Rich/RichFutureRecCheckers/CMakeLists.txt | 1 + .../src/RichMCDetectorHits.cpp | 162 ++++++++++++++++++ .../ConfiguredRecoMonitors.py | 25 ++- 3 files changed, 181 insertions(+), 7 deletions(-) create mode 100644 Rich/RichFutureRecCheckers/src/RichMCDetectorHits.cpp diff --git a/Rich/RichFutureRecCheckers/CMakeLists.txt b/Rich/RichFutureRecCheckers/CMakeLists.txt index 9d664de5879..e558c30fd52 100644 --- a/Rich/RichFutureRecCheckers/CMakeLists.txt +++ b/Rich/RichFutureRecCheckers/CMakeLists.txt @@ -22,6 +22,7 @@ gaudi_add_module(RichFutureRecCheckers src/RichMCTrackResolution.cpp src/RichPIDQC.cpp src/RichMCOpticalPhotons.cpp + src/RichMCDetectorHits.cpp LINK Boost::headers Gaudi::GaudiAlgLib diff --git a/Rich/RichFutureRecCheckers/src/RichMCDetectorHits.cpp b/Rich/RichFutureRecCheckers/src/RichMCDetectorHits.cpp new file mode 100644 index 00000000000..97c30f0beb7 --- /dev/null +++ b/Rich/RichFutureRecCheckers/src/RichMCDetectorHits.cpp @@ -0,0 +1,162 @@ +/*****************************************************************************\ +* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration * +* * +* This software is distributed under the terms of the GNU General Public * +* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". * +* * +* In applying this licence, CERN does not waive the privileges and immunities * +* granted to it by virtue of its status as an Intergovernmental Organization * +* or submit itself to any jurisdiction. * +\*****************************************************************************/ + +// base class +#include "RichFutureRecBase/RichRecHistoAlgBase.h" + +// Gaudi Functional +#include "LHCbAlgs/Consumer.h" + +// Rich Utils +#include "RichUtils/RichDAQDefinitions.h" + +// Relations +#include "RichFutureMCUtils/RichSmartIDMCUtils.h" + +// event data +#include "Event/MCRichDigitSummary.h" + +// Rec Event +#include "RichFutureRecEvent/RichRecSIMDPixels.h" + +// STD +#include +#include +#include +#include + +namespace Rich::Future::Rec::MC::Moni { + + /** @class DetectorHits RichDetectorHits.h + * + * Monitors the data sizes in the RICH detectors using MC + * + * @author Chris Jones + * @date 2023-11-15 + */ + + class DetectorHits final : public LHCb::Algorithm::Consumer> { + + public: + /// Standard constructor + DetectorHits( const std::string& name, ISvcLocator* pSvcLocator ) + : Consumer( name, pSvcLocator, + // input data + {KeyValue{"RichSIMDPixelSummariesLocation", SIMDPixelSummariesLocation::Default}, + KeyValue{"RichDigitSummariesLocation", LHCb::MCRichDigitSummaryLocation::Default}} ) {} + + public: + /// Functional operator + void operator()( const SIMDPixelSummaries& pixels, // + const LHCb::MCRichDigitSummarys& digitSums ) const override { + + using namespace Rich::Future::MC::Relations; + + // Count hits per channel ID + std::unordered_map hitsPerID; + + // MC helper + SmartIDUtils smartIDMCHelper( digitSums ); + + // Loop over RICHes + for ( const auto rich : activeDetectors() ) { + // Loop over pixels for this RICH + for ( const auto& pix : pixels.range( rich ) ) { + // scalar loop + for ( std::size_t i = 0; i < SIMDPixel::SIMDFP::Size; ++i ) { + if ( !pix.validMask()[i] ) { continue; } + // Hit ID + const auto id = pix.smartID()[i]; + // count hits per channel (strip time info) + ++hitsPerID[id.channelDataOnly()]; + } // scalar loop + } // SIMD pixels + } // RICHes + + // the lock + std::lock_guard lock( m_updateLock ); + + // loop again over collected hits + for ( const auto& [id, count] : hitsPerID ) { + + // Get MC truth for this ID + const auto histCodes = smartIDMCHelper.mcDigitHistoryCodes( id ); + + // which RICH + const auto rich = id.rich(); + + // Count signal hits in this channel + const auto nSignal = + std::count_if( histCodes.begin(), histCodes.end(), // + [&rich]( const auto& code ) { + return ( code.signalEvent() && ( ( Rich::Rich1 == rich && code.c4f10Hit() ) || + ( Rich::Rich2 == rich && code.cf4Hit() ) ) ); + } ); + const auto sigFrac = ( !histCodes.empty() ? (double)nSignal / (double)histCodes.size() : 0.0 ); + + // fill histos + fillHisto( h_hitsPerID[rich], histCodes.size() ); + fillHisto( h_digitsPerID[rich], count ); + fillHisto( h_signalPerID[rich], nSignal ); + fillHisto( h_sigFracVOcc[rich], histCodes.size(), sigFrac ); + } + } + + protected: + /// Pre-Book all histograms + StatusCode prebookHistograms() override { + bool ok = true; + + const std::size_t nHitMax = 30; + + for ( const auto rich : activeDetectors() ) { + ok &= saveAndCheck( h_hitsPerID[rich], // + richHisto1D( Rich::HistogramID( "nHitsPerChanID", rich ), // + "# hits per channel ID (total nHits>0)", // + -0.5, nHitMax + 0.5, nHitMax + 1 ) ); + ok &= saveAndCheck( h_signalPerID[rich], // + richHisto1D( Rich::HistogramID( "nSignalHitsPerChanID", rich ), // + "# signal hits per channel ID (total nHits>0)", // + -0.5, nHitMax + 0.5, nHitMax + 1 ) ); + ok &= saveAndCheck( h_digitsPerID[rich], // + richHisto1D( Rich::HistogramID( "nDigitsPerChanID", rich ), // + "# digitised hits per channel ID (total nHits>0)", // + -0.5, nHitMax + 0.5, nHitMax + 1 ) ); + ok &= saveAndCheck( h_sigFracVOcc[rich], // + richProfile1D( Rich::HistogramID( "sigFracVOcc", rich ), // + "Signal fraction versus occupancy (total nHits>0)", // + -0.5, nHitMax + 0.5, nHitMax + 1 ) ); + } + + return StatusCode{ok}; + } + + private: + // data + + /// mutex lock + mutable std::mutex m_updateLock; + + private: + // cached histograms + + DetectorArray h_hitsPerID = {{}}; + DetectorArray h_signalPerID = {{}}; + DetectorArray h_digitsPerID = {{}}; + DetectorArray h_sigFracVOcc = {{}}; + }; + + // Declaration of the Algorithm Factory + DECLARE_COMPONENT( DetectorHits ) + +} // namespace Rich::Future::Rec::MC::Moni diff --git a/Rich/RichFutureRecMonitors/python/RichFutureRecMonitors/ConfiguredRecoMonitors.py b/Rich/RichFutureRecMonitors/python/RichFutureRecMonitors/ConfiguredRecoMonitors.py index 5ec67d801af..4f17a382bb8 100644 --- a/Rich/RichFutureRecMonitors/python/RichFutureRecMonitors/ConfiguredRecoMonitors.py +++ b/Rich/RichFutureRecMonitors/python/RichFutureRecMonitors/ConfiguredRecoMonitors.py @@ -141,6 +141,8 @@ def RichRecMonitors( from Configurables import Rich__Future__Rec__Moni__DetectorHits as DetectorHits richHits = makeRichAlg(DetectorHits, "RichRecPixelQC" + GroupName, algprops) + richHits.RichSIMDPixelSummariesLocation = cLocs[ + "RichSIMDPixelSummariesLocation"] all.Members += [richHits] # clusters @@ -386,21 +388,21 @@ def RichRecCheckers( # List of monitors to run checkers={ "Expert": [ - "PIDPerformance", "PhotonCherenkovAngles", + "RichHits", "PIDPerformance", "PhotonCherenkovAngles", "CherenkovResolution", "CKResParameterisation", "TrackResolution", "MCOpticalPhotons", "Time" ], "OfflineFull": [ - "PIDPerformance", "PhotonCherenkovAngles", "TrackResolution", - "Time" + "RichHits", "PIDPerformance", "PhotonCherenkovAngles", + "TrackResolution", "Time" ], "OfflineExpress": [ - "PIDPerformance", "PhotonCherenkovAngles", "TrackResolution", - "Time" + "RichHits", "PIDPerformance", "PhotonCherenkovAngles", + "TrackResolution", "Time" ], "Online": [ - "PIDPerformance", "PhotonCherenkovAngles", "TrackResolution", - "Time" + "RichHits", "PIDPerformance", "PhotonCherenkovAngles", + "TrackResolution", "Time" ], "None": [] }, @@ -494,6 +496,15 @@ def RichRecCheckers( checks = checkers[histograms] + extra_checkers + # Rich Hits + if "RichHits" in checks: + from Configurables import Rich__Future__Rec__MC__Moni__DetectorHits as DetectorHits + richHits = makeRichAlg(DetectorHits, "RichMCHits" + GroupName, + algprops) + richHits.RichSIMDPixelSummariesLocation = cLocs[ + "RichSIMDPixelSummariesLocation"] + all.Members += [richHits] + # Loop over tracks for tktype, trackLocation in sorted(inputTrackLocations.items()): -- GitLab From 065bed21c240de60b4156c3c1692cae780b29fae Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Tue, 28 Nov 2023 10:34:35 +0000 Subject: [PATCH 02/23] RichRecSIMDPixels: handle maintaining count internally via emplace_back --- .../RichFutureRecEvent/RichRecSIMDPixels.h | 33 +++++++++++-------- .../src/RichSIMDSummaryPixels.cpp | 7 ++-- 2 files changed, 25 insertions(+), 15 deletions(-) diff --git a/Rich/RichFutureRecEvent/include/RichFutureRecEvent/RichRecSIMDPixels.h b/Rich/RichFutureRecEvent/include/RichFutureRecEvent/RichRecSIMDPixels.h index 7aa1417b3e3..462c07bc894 100644 --- a/Rich/RichFutureRecEvent/include/RichFutureRecEvent/RichRecSIMDPixels.h +++ b/Rich/RichFutureRecEvent/include/RichFutureRecEvent/RichRecSIMDPixels.h @@ -11,13 +11,6 @@ #pragma once -// STL -#include -#include -#include -#include -#include - // Utils #include "RichUtils/RichPixelCluster.h" #include "RichUtils/RichSIMDTypes.h" @@ -34,6 +27,13 @@ // Event #include "RichFutureRecEvent/RichRecSpacePoints.h" +// STL +#include +#include +#include +#include +#include + namespace Rich::Future::Rec { /** @class SIMDPixel RichFutureRecEvent/RichRecSIMDPixels.h @@ -68,8 +68,7 @@ namespace Rich::Future::Rec { /// Constructor from RICH and panel ID SIMDPixel( const Rich::DetectorType rich, // const Rich::Side side ) - : m_rich( rich ) // - , m_side( side ) {} + : m_rich( rich ), m_side( side ) {} /// Constructor from full data SIMDPixel( const Rich::DetectorType rich, // @@ -315,10 +314,12 @@ namespace Rich::Future::Rec { public: // setters - /// Count hits in each RICH and panel - inline void addHit( const Rich::DetectorType rich, // - const Rich::Side side ) noexcept { - ++m_nHitsSIMD[rich][side]; + /// emplace back a SIMD pixel into the container + template + inline auto& emplace_back( ARGS&&... args ) { + auto& pix = SIMD::STDVector::emplace_back( std::forward( args )... ); + countHit( pix.rich(), pix.side() ); + return pix; } /// Initialise the ranges once filled @@ -341,6 +342,12 @@ namespace Rich::Future::Rec { } private: + /// Count hits in each RICH and panel + inline void countHit( const Rich::DetectorType rich, // + const Rich::Side side ) noexcept { + ++m_nHitsSIMD[rich][side]; + } + /// Set RICH range inline void setRange( const Rich::DetectorType rich, // const std::size_t beginit, // diff --git a/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDSummaryPixels.cpp b/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDSummaryPixels.cpp index 367027b0206..d2657bbe21c 100644 --- a/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDSummaryPixels.cpp +++ b/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDSummaryPixels.cpp @@ -149,8 +149,6 @@ SIMDSummaryPixels::operator()( const Rich::PDPixelCluster::Vector& clusters, // summaries.emplace_back( lastRich, lastSide, smartIDs, // gPos, lPos, // effArea, scClusIn, mask ); - // count - summaries.addHit( lastRich, lastSide ); // reset index = 0; }; @@ -181,6 +179,8 @@ SIMDSummaryPixels::operator()( const Rich::PDPixelCluster::Vector& clusters, // // Scalar cluster index std::size_t clusIn = 0; + _ri_debug << "Found " << clusters.size() << " clusters" << endmsg; + // Note relying on fact they are sorted by panel and rich here. for ( const auto&& [cluster, gloPos] : Ranges::ConstZip( clusters, gPoints ) ) { @@ -258,6 +258,9 @@ SIMDSummaryPixels::operator()( const Rich::PDPixelCluster::Vector& clusters, // // Initialise the summaries, sets ranges etc. summaries.init(); + _ri_debug << "Created for RICH1/RICH2 " << summaries.nHitsScalar( Rich::Rich1 ) << " / " + << summaries.nHitsScalar( Rich::Rich2 ) << " hits" << endmsg; + // ------------------------------------------------------------------------------------------- // Some verbose bugging printout -- GitLab From 5aefb2747ada81dec481ed89d65fb56673d462c9 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Tue, 28 Nov 2023 10:35:40 +0000 Subject: [PATCH 03/23] RichDetectorHits: Use SIMD pixels rather than raw decoded data --- .../src/RichDetectorHits.cpp | 203 ++++++++---------- 1 file changed, 95 insertions(+), 108 deletions(-) diff --git a/Rich/RichFutureRecMonitors/src/RichDetectorHits.cpp b/Rich/RichFutureRecMonitors/src/RichDetectorHits.cpp index 6b3f230bbee..f5338aea5fc 100644 --- a/Rich/RichFutureRecMonitors/src/RichDetectorHits.cpp +++ b/Rich/RichFutureRecMonitors/src/RichDetectorHits.cpp @@ -16,13 +16,17 @@ #include "LHCbAlgs/Consumer.h" // Rich Utils -#include "RichFutureUtils/RichDecodedData.h" -#include "RichFutureUtils/RichSmartIDs.h" #include "RichUtils/RichDAQDefinitions.h" +// Rec Event +#include "RichFutureRecEvent/RichRecSIMDPixels.h" + // AIDA #include "AIDA/IAxis.h" +// boost +#include "boost/container/small_vector.hpp" + // STD #include #include @@ -39,39 +43,23 @@ namespace Rich::Future::Rec::Moni { * @date 2016-12-06 */ - class DetectorHits final : public LHCb::Algorithm::Consumer< - void( const DAQ::DecodedData&, // - const Rich::Utils::RichSmartIDs& ), // - LHCb::DetDesc::usesBaseAndConditions> { + class DetectorHits final : public LHCb::Algorithm::Consumer> { public: /// Standard constructor DetectorHits( const std::string& name, ISvcLocator* pSvcLocator ) : Consumer( name, pSvcLocator, // input data - {KeyValue{"DecodedDataLocation", DAQ::DecodedDataLocation::Default}, - // input conditions data - KeyValue{"RichSmartIDs", Rich::Utils::RichSmartIDs::DefaultConditionKey}} ) { + {KeyValue{"RichSIMDPixelSummariesLocation", SIMDPixelSummariesLocation::Default}} ) { // print some stats on the final plots setProperty( "HistoPrint", true ).ignore(); setProperty( "NBins2DHistos", 50 ).ignore(); } - /// Initialize - StatusCode initialize() override { - return Consumer::initialize().andThen( [&] { Rich::Utils::RichSmartIDs::addConditionDerivation( this ); } ); - } - public: /// Functional operator - void operator()( const DAQ::DecodedData& data, const Rich::Utils::RichSmartIDs& smartIDsHelper ) const override { - - // Count active PDs in each RICH - DetectorArray activePDs = {{}}; - // Count hits in each RICH - DetectorArray richHits = {{}}; - // Count active modules - DetectorArray richModules = {{}}; + void operator()( const SIMDPixelSummaries& pixels ) const override { // the lock std::lock_guard lock( m_updateLock ); @@ -79,85 +67,97 @@ namespace Rich::Future::Rec::Moni { // Loop over RICHes for ( const auto rich : activeDetectors() ) { - // count for occupancy in each PD - using VD = std::vector; - std::vector pdOccXY{nBins2D(), VD( nBins2D(), 0.0 )}; - - // data for this RICH - const auto& rD = data[rich]; - - // sides per RICH - for ( const auto& pD : rD ) { - - // PD modules per side - richModules[rich] += pD.size(); - for ( const auto& mD : pD ) { - - // PDs per module - for ( const auto& PD : mD ) { - - // PD ID - const auto pdID = PD.pdID(); - if ( pdID.isValid() ) { - - // Vector of SmartIDs - const auto& rawIDs = PD.smartIDs(); - - // Do we have any hits - if ( !rawIDs.empty() ) { - - // Side - const auto side = pdID.panel(); - - // Fill average HPD occ plot - fillHisto( h_nTotalPixsPerPD[rich], rawIDs.size() ); - if ( Rich::Rich2 == rich ) { - fillHisto( ( pdID.isHTypePMT() ? h_nTotalPixsPerPDR2H : h_nTotalPixsPerPDR2R ), rawIDs.size() ); - } + // ID for last seen PD + LHCb::RichSmartID lastPD; - // count active PDs - ++activePDs[rich]; + // count hits in current PD + std::size_t curPDhits{0}; - // count hits - richHits[rich] += rawIDs.size(); + // cache iX and iY indices for hits in current PD + boost::container::small_vector, LHCb::RichSmartID::MaPMT::TotalPixels> hitiXiY; - // loop over hits - for ( const auto id : rawIDs ) { - // get global/local coords - const auto gPos = smartIDsHelper.globalPosition( id ); - const auto lPos = smartIDsHelper.globalToPDPanel( gPos ); - // fill plots - fillHisto( h_pixXYGlo[rich][side], gPos.X(), gPos.Y() ); - fillHisto( h_pixXGlo[rich][side], gPos.X() ); - fillHisto( h_pixYGlo[rich][side], gPos.Y() ); - fillHisto( h_pixXYLoc[rich], lPos.X(), lPos.Y() ); - fillHisto( h_pixXLoc[rich], lPos.X() ); - fillHisto( h_pixYLoc[rich], lPos.Y() ); - // cache occupancy data - const auto iX = - std::clamp( h_pdXYOcc[rich]->xAxis().coordToIndex( lPos.X() ), 0, (int)nBins2D() - 1 ); - const auto iY = - std::clamp( h_pdXYOcc[rich]->yAxis().coordToIndex( lPos.Y() ), 0, (int)nBins2D() - 1 ); - pdOccXY[iX][iY] = ( 100.0 * rawIDs.size() ) / (double)LHCb::RichSmartID::MaPMT::TotalPixels; - } - - } // PD has hits - - } // PDID is valid - - } // PDs + // count for occupancy in each PD + using VD = std::vector; + std::vector pdOccXY{nBins2D(), VD( nBins2D(), 0.0 )}; - } // modules + // Count active PDs in each RICH + std::uint32_t activePDs = 0; + + // functor to fill data for a given PD + auto newPD = [&]() { + // If last PD ID is valid fill plots etc. for the last PD + if ( lastPD.isValid() ) { + // Fill average HPD occ plot + fillHisto( h_nTotalPixsPerPD[rich], curPDhits ); + if ( Rich::Rich2 == rich ) { + fillHisto( ( lastPD.isHTypePMT() ? h_nTotalPixsPerPDR2H : h_nTotalPixsPerPDR2R ), curPDhits ); + } + // PD occ + for ( const auto& [iX, iY] : hitiXiY ) { + pdOccXY[iX][iY] = ( 100.0 * curPDhits ) / (double)LHCb::RichSmartID::MaPMT::TotalPixels; + } + } + // reset caches + hitiXiY.clear(); + curPDhits = 0; + }; - } // panels + // Fill RICH hits plots + const auto nHits = pixels.nHitsScalar( rich ); + if ( nHits > 0 ) { fillHisto( h_nTotalPixs[rich], nHits ); } + + // Loop over pixels for this RICH + for ( const auto& pix : pixels.range( rich ) ) { + // scalar loop + for ( std::size_t i = 0; i < SIMDPixel::SIMDFP::Size; ++i ) { + if ( !pix.validMask()[i] ) { continue; } + + // Hit data + const auto id = pix.smartID()[i]; + const auto pdID = id.pdID(); + const auto side = pdID.panel(); + + // If new PD fill plots and caches + if ( lastPD != pdID ) { + newPD(); + lastPD = pdID; + // count active PDs + ++activePDs; + } + + // Count hits in current PD + ++curPDhits; + + // get global/local coords + const auto gPos = pix.gloPos( i ); + const auto lPos = pix.locPos( i ); + + // fill plots + fillHisto( h_pixXYGlo[rich][side], gPos.X(), gPos.Y() ); + fillHisto( h_pixXGlo[rich][side], gPos.X() ); + fillHisto( h_pixYGlo[rich][side], gPos.Y() ); + fillHisto( h_pixXYLoc[rich], lPos.X(), lPos.Y() ); + fillHisto( h_pixXLoc[rich], lPos.X() ); + fillHisto( h_pixYLoc[rich], lPos.Y() ); + + // cache occupancy data + const auto iX = std::clamp( h_pdXYOcc[rich]->xAxis().coordToIndex( lPos.X() ), 0, (int)nBins2D() - 1 ); + const auto iY = std::clamp( h_pdXYOcc[rich]->yAxis().coordToIndex( lPos.Y() ), 0, (int)nBins2D() - 1 ); + hitiXiY.emplace_back( iX, iY ); + + } // scalar loop + } // SIMD pixel loop + + // call once for last hit + newPD(); // Fill occupancies for ( auto iX = 0u; iX < nBins2D(); ++iX ) { for ( auto iY = 0u; iY < nBins2D(); ++iY ) { - const auto X = - 0.5 * ( h_pdXYOcc[rich]->xAxis().binLowerEdge( iX ) + h_pdXYOcc[rich]->xAxis().binUpperEdge( iX ) ); - const auto Y = - 0.5 * ( h_pdXYOcc[rich]->yAxis().binLowerEdge( iY ) + h_pdXYOcc[rich]->yAxis().binUpperEdge( iY ) ); + const auto X = 0.5 * ( h_pdXYOcc[rich]->xAxis().binLowerEdge( iX ) + // + h_pdXYOcc[rich]->xAxis().binUpperEdge( iX ) ); + const auto Y = 0.5 * ( h_pdXYOcc[rich]->yAxis().binLowerEdge( iY ) + // + h_pdXYOcc[rich]->yAxis().binUpperEdge( iY ) ); const auto occ = pdOccXY[iX][iY]; fillHisto( h_pdXYOcc[rich], X, Y, occ ); fillHisto( h_pdXOcc[rich], X, occ ); @@ -165,18 +165,10 @@ namespace Rich::Future::Rec::Moni { } } - } // RICHes - - // Loop over RICHes - for ( const auto rich : activeDetectors() ) { // Fill active PD plots - fillHisto( h_nActivePDs[rich], activePDs[rich] ); - // Fill RICH hits plots - if ( richHits[rich] > 0 ) { - fillHisto( h_nTotalPixs[rich], richHits[rich] ); - fillHisto( h_nActiveModules[rich], richModules[rich] ); - } - } + fillHisto( h_nActivePDs[rich], activePDs ); + + } // RICHes } protected: @@ -212,9 +204,6 @@ namespace Rich::Future::Rec::Moni { ok &= saveAndCheck( h_nActivePDs[rich], // richHisto1D( Rich::HistogramID( "nActivePDs", rich ), // "# Active PDs (nHits>0)", -0.5, 2000.5, 2001 ) ); - ok &= saveAndCheck( h_nActiveModules[rich], // - richHisto1D( Rich::HistogramID( "nActiveModules", rich ), // - "# Active PD Modules (nHits>0)", -0.5, 150.5, 151 ) ); // hit points (global) for ( const auto side : Rich::sides() ) { const auto minX = @@ -309,8 +298,6 @@ namespace Rich::Future::Rec::Moni { DetectorArray h_nTotalPixs = {{}}; /// Number active PDs histograms DetectorArray h_nActivePDs = {{}}; - /// Number active module histograms - DetectorArray h_nActiveModules = {{}}; DetectorArray> h_pixXYGlo = {{}}; DetectorArray> h_pixXGlo = {{}}; -- GitLab From 50d3b470a14eab0ea5dd0c6f1e1150a507646368 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Tue, 28 Nov 2023 10:36:23 +0000 Subject: [PATCH 04/23] ConfiguredRichReco.py: Set default max PD occupancy to 'large number' --- .../python/RichFutureRecSys/ConfiguredRichReco.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py b/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py index e86a778efe8..9efec26fd88 100644 --- a/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py +++ b/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py @@ -304,7 +304,7 @@ def RichRecoSequence( maxClusters=200000, # Maximum PD occupancy - maxPDOccupancy=(64, 64), + maxPDOccupancy=(999999, 999999), # Should pixel clustering be run, for either ( RICH1, RICH2 ) applyPixelClustering=(False, False), -- GitLab From e12ae322b8eff8734708687efdba222ffe12fe0e Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Tue, 28 Nov 2023 10:37:09 +0000 Subject: [PATCH 05/23] RichPixelUseMCInfo: Add support for emulating arbitrary pixel sizes --- .../src/RichPixelUseMCInfo.cpp | 125 ++++++++++++++++-- 1 file changed, 111 insertions(+), 14 deletions(-) diff --git a/Rich/RichFutureRecMCAlgorithms/src/RichPixelUseMCInfo.cpp b/Rich/RichFutureRecMCAlgorithms/src/RichPixelUseMCInfo.cpp index d85c14b4630..239d3f81894 100644 --- a/Rich/RichFutureRecMCAlgorithms/src/RichPixelUseMCInfo.cpp +++ b/Rich/RichFutureRecMCAlgorithms/src/RichPixelUseMCInfo.cpp @@ -27,8 +27,15 @@ // Utils #include "RichFutureUtils/RichSmartIDs.h" +// Boost +#include "boost/container/small_vector.hpp" + // STD +#include +#include #include +#include +#include namespace Rich::Future::Rec::MC { @@ -64,10 +71,28 @@ namespace Rich::Future::Rec::MC { return Transformer::initialize().andThen( [&] { // create the RICH smartID helper instance Rich::Utils::RichSmartIDs::addConditionDerivation( this ); + if ( m_mcPixQuantizeInner > 0 || m_mcPixQuantizeOuter > 0 ) { + info() << "Will quantise inner/outer pixels to " << m_mcPixQuantizeInner.value() << "/" + << m_mcPixQuantizeOuter.value() << " mm square" << endmsg; + } // setProperty( "OutputLevel", MSG::VERBOSE ).ignore(); } ); } + private: + // types + + struct XYData { + double x{0}, y{0}; + SIMDPixel* pix{nullptr}; + std::size_t index{0}; + auto time() const { return ( pix ? pix->hitTime()[index] : -999.0f ); } + void disable() { + if ( pix ) { pix->validMask()[index] = false; } + } + using Vector = boost::container::small_vector; + }; + public: /// Functional operator SIMDPixelSummaries operator()( SIMDPixelSummaries const& pixels, // @@ -141,9 +166,19 @@ namespace Rich::Future::Rec::MC { // Clone the pixels to update them SIMDPixelSummaries out_pixs = pixels; + // local quantized (x,y) hits made for each PD + std::map madeXY; + + // disabled hits + DetectorArray disHits{0, 0}; + // Loop over summary data and update positions for ( auto& pix : out_pixs ) { + // Which rich and side + const auto rich = pix.rich(); + const auto side = pix.side(); + // flag to indicate if positions have been updated bool hasPosUpdate = false; @@ -163,6 +198,7 @@ namespace Rich::Future::Rec::MC { // Find best MCHit to use const auto* mch = &pix_mchits.front(); if ( pix_mchits.size() > 1 ) { + // If ID has time set use that to select MC hit if ( id.adcTimeIsSet() ) { double minTimeDiff = 999999.9; @@ -187,7 +223,7 @@ namespace Rich::Future::Rec::MC { } } } - } + } // more than one MC hit if ( m_useMCPos ) { // Update hit position using MC @@ -196,18 +232,57 @@ namespace Rich::Future::Rec::MC { auto mc_lpos = smartIDsHelper.globalToPDPanel( mc_gpos ); // Apply Z correction const DetectorArray zCorr{-10.8, -10.8}; - mc_lpos.SetZ( mc_lpos.Z() + zCorr[pix.rich()] ); - // Shift back to global - mc_gpos = smartIDsHelper.globalPosition( mc_lpos, pix.rich(), pix.side() ); - _ri_verbo << " -> original LPos " << i << " " << pix.locPos( i ) << endmsg; - _ri_verbo << " -> MC update LPos " << i << " " << mc_lpos << endmsg; - _ri_verbo << " -> original GPos " << i << " " << pix.gloPos( i ) << endmsg; - _ri_verbo << " -> MC update GPos " << i << " " << mc_gpos << endmsg; - // Update position for this scalar - pix.gloPos().X()[i] = mc_gpos.X(); - pix.gloPos().Y()[i] = mc_gpos.Y(); - pix.gloPos().Z()[i] = mc_gpos.Z(); - hasPosUpdate = true; + mc_lpos.SetZ( mc_lpos.Z() + zCorr[rich] ); + // If requested quantize the pixel position inside the pixel, + // to emulate smaller pixels + const bool isInner = ( fabs( mc_lpos.X() ) < m_innerPixX[rich] && // + fabs( mc_lpos.Y() ) < m_innerPixY[rich] ); + const auto pixSize = ( isInner ? m_mcPixQuantizeInner : m_mcPixQuantizeOuter ); + if ( pixSize > 0.0 ) { + auto quantize = [&pixSize]( const auto x ) { return std::round( x / pixSize ) * pixSize; }; + mc_lpos.SetX( quantize( mc_lpos.X() ) ); + mc_lpos.SetY( quantize( mc_lpos.Y() ) ); + // have we already made a hit at this quantized (x,y) for this PD + auto& xys = madeXY[id.pdID()]; + const XYData new_xyD{mc_lpos.X(), mc_lpos.Y(), &pix, i}; + auto old_xyD = std::find_if( xys.begin(), xys.end(), [&new_xyD]( const auto& d ) { + const double tol = 1.0e-4; + return ( std::fabs( d.x - new_xyD.x ) < tol && std::fabs( d.y - new_xyD.y ) < tol ); + } ); + if ( old_xyD != xys.end() ) { + // which hit is the earilest ? + _ri_verbo << "Found existing hit for (x,y) = " << mc_lpos.X() << "," << mc_lpos.Y() << " time " + << pix.hitTime()[i] << " prev time " << old_xyD->time() << endmsg; + if ( pix.hitTime()[i] < old_xyD->time() ) { + _ri_verbo << " -> early, keep and disable old" << endmsg; + // keep this hit, disable the previous one + old_xyD->disable(); + *old_xyD = new_xyD; + ++disHits[rich]; + } else { + _ri_verbo << " -> late, disable" << endmsg; + // disable this hit, keep the previous one + pix.validMask()[i] = false; + ++disHits[rich]; + } + } else { + // Not seen this (x,y) before so just add to the list + xys.push_back( new_xyD ); + } + } + if ( pix.validMask()[i] ) { + // Shift back to global + mc_gpos = smartIDsHelper.globalPosition( mc_lpos, rich, side ); + _ri_verbo << " -> Original LPos " << i << " " << pix.locPos( i ) << endmsg; + _ri_verbo << " -> MC update LPos " << i << " " << mc_lpos << endmsg; + _ri_verbo << " -> Original GPos " << i << " " << pix.gloPos( i ) << endmsg; + _ri_verbo << " -> MC update GPos " << i << " " << mc_gpos << endmsg; + // Update position for this scalar + pix.gloPos().X()[i] = mc_gpos.X(); + pix.gloPos().Y()[i] = mc_gpos.Y(); + pix.gloPos().Z()[i] = mc_gpos.Z(); + hasPosUpdate = true; + } } if ( m_useMCTime ) { @@ -225,13 +300,23 @@ namespace Rich::Future::Rec::MC { } // scalar loop // finally update SIMD local positions if needed - if ( hasPosUpdate ) { pix.locPos() = smartIDsHelper.globalToPDPanel( pix.rich(), pix.side(), pix.gloPos() ); } + if ( hasPosUpdate ) { pix.locPos() = smartIDsHelper.globalToPDPanel( rich, side, pix.gloPos() ); } } // SIMD pixel loop + // Count disabled hits per event + for ( const auto rich : {Rich::Rich1, Rich::Rich2} ) { m_disabledHits[rich] += disHits[rich]; } + return out_pixs; } + private: + // counters + + /// # Disabled Hits + mutable DetectorArray> m_disabledHits{ + {{this, "# Rich1 Disabled Hits"}, {this, "# Rich2 Disabled Hits"}}}; + private: // properties @@ -253,6 +338,18 @@ namespace Rich::Future::Rec::MC { Gaudi::Property> m_timeShift{ this, "TimeCalib", {0., 40.}, "Global time shift for each RICH, to get both to same calibrated point"}; + /// Inner region spatial quantisation to apply to MC local coordinate hits (<0 means do not apply) + Gaudi::Property m_mcPixQuantizeInner{this, "InnerPixelQuantization", -1.0 * Gaudi::Units::mm}; + + /// Outer region spatial quantisation to apply to MC local coordinate hits (<0 means do not apply) + Gaudi::Property m_mcPixQuantizeOuter{this, "OuterPixelQuantization", -1.0 * Gaudi::Units::mm}; + + /// Size in X defining the inner pixels for each RICH + Gaudi::Property> m_innerPixX{this, "InnerPixelsX", {99999.9, 99999.9}}; + + /// Size in Y defining the inner pixels for each RICH + Gaudi::Property> m_innerPixY{this, "InnerPixelsY", {99999.9, 250.0}}; + private: // Handles to access MCRichHits other than Signal, as they are not necessary there DataObjectReadHandle m_prevHits{this, "PrevLocation", "Prev/" + LHCb::MCRichHitLocation::Default}; -- GitLab From a10eac0ea930df7a2574d97df1c6bbc19973ce9e Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Tue, 28 Nov 2023 10:37:47 +0000 Subject: [PATCH 06/23] RichSIMDQuarticPhotonReco: Prefer std::distance to manual subtraction --- .../src/RichSIMDQuarticPhotonReco.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp index e7ff4a8f346..9aba7e5a2aa 100644 --- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp +++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp @@ -9,9 +9,6 @@ * or submit itself to any jurisdiction. * \*****************************************************************************/ -// STL -#include - // Base class #include "RichCommonQuarticPhotonReco.h" @@ -29,6 +26,13 @@ #include "Kernel/RichDetectorType.h" #include "Kernel/RichSide.h" +// STL +#include +#include +#include +#include +#include + namespace Rich::Future::Rec { namespace { @@ -324,7 +328,7 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& seg : pixels.range( Rich::Rich2, Rich::right ) ) ); // SIMD Pixel index in container (start index) for this range - Relations::PhotonToParents::PixelIndex pixIndex = pixR.begin() - pixels.begin(); + Relations::PhotonToParents::PixelIndex pixIndex = std::distance( pixels.begin(), pixR.begin() ); // loop over selected pixels for ( auto ipix = pixR.begin(); ipix != pixR.end(); ++ipix, ++pixIndex ) { -- GitLab From 282718d178c56cd4a364eac733cfe5a73e148858 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Tue, 28 Nov 2023 10:39:55 +0000 Subject: [PATCH 07/23] RichFuture.py: Add example options to control emulated pixel sizes --- Rich/RichFutureRecSys/examples/RichFuture.py | 110 ++++++++++--------- 1 file changed, 59 insertions(+), 51 deletions(-) diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py index 8711b0bb10b..6e9e85825bd 100644 --- a/Rich/RichFutureRecSys/examples/RichFuture.py +++ b/Rich/RichFutureRecSys/examples/RichFuture.py @@ -58,9 +58,9 @@ rootFileBaseName = "RichFuture-" + myBuild + "-" + myConf + "-" + histos HistogramPersistencySvc().OutputFile = rootFileBaseName + "-Histos.root" # Event numbers -nEvents = (1000 if not batchMode else 99999999) -EventSelector().PrintFreq = (100 if not batchMode else 250) -#LHCbApp().SkipEvents = 2 +nEvents = (100 if not batchMode else 99999999) +EventSelector().PrintFreq = (10 if not batchMode else 250) +#LHCbApp().SkipEvents = 2 # Just to initialise DDDBConf() @@ -78,9 +78,9 @@ SequencerTimerTool("ToolSvc.SequencerTimerTool").NameSize = 40 # Auditors AuditorSvc().Auditors += ["FPEAuditor"] FPEAuditor().ActivateAt = ["Execute"] -#AuditorSvc().Auditors += [ "NameAuditor" ] +#AuditorSvc().Auditors += ["NameAuditor"] -# The overall sequence. Filled below. +# The overall sequence.Filled below. all = GaudiSequencer("All", MeasureTime=True) # -------------------------------------------------------------------------------------- @@ -109,18 +109,18 @@ useMCTracks = True useMCHits = True # Get time windows from environment vars if set -# otherwise use default values as passed. -PixWin = float(os.getenv("PIXWIN", "1.0")) -PhotWin = float(os.getenv("PHOTWIN", "0.1")) +# otherwise use default values as passed. +PixWin = float(os.getenv("PIXWIN", "1.000")) +PhotWin = float(os.getenv("PHOTWIN", "0.075")) # Enable 4D reco in each RICH is4D = "ENABLE4D" in os.environ enable4D = (is4D, is4D) # Average expected time for signal avSignalTime = (13.03, 52.94) -# Pixel time window size in each RICH (in ns) +# Pixel time window size in each RICH( in ns ) pixTimeWindow = (PixWin, PixWin) -# Photon reco time window size in each RICH (ns) +# Photon reco time window size in each RICH( ns ) photTimeWindow = (PhotWin, PhotWin) # Min Momentum cuts @@ -128,19 +128,19 @@ minP = (1.0 * GeV, 1.0 * GeV, 1.0 * GeV) minPt = 0.2 * GeV # Add time to reco tracks using MC -# Only used if useMCTracks=False and 4D enabled +# Only used if useMCTracks = False and 4D enabled tkSegAddTimeFromMC = ((enable4D[0] or enable4D[1]) and not useMCTracks) # Add time info to Rich hits -# Only if useMCHits=False and 4D enabled +# Only if useMCHits = False and 4D enabled usePixelMCTime = ((enable4D[0] or enable4D[1]) and not useMCHits) # Cheat pixel positions using MC info -usePixelMCPos = False +usePixelMCPos = True # Filter 'real' tracks based on MC info. # Designed to help match those made from MCRichTracks, where certain -# track types cannot be made (e.g. ghosts, or BT tracks). +# track types cannot be made( e.g.ghosts, or BT tracks ). filterTksUsingMC = False # Unpacking sequence before reco is run @@ -176,6 +176,13 @@ if UseDD4Hep: if tkSegAddTimeFromMC or useMCHits or useMCTracks or filterTksUsingMC: preUnpackSeq.Members += [UnpackMCVertex(), UnpackMCParticle()] +if usePixelMCPos: + from Configurables import Rich__Future__Rec__MC__RichPixelUseMCInfo as PixUseMCInfo + innerPixQ = float(os.getenv("INNERPIXQUANT", "2.8")) + outerPixQ = float(os.getenv("OUTERPIXQUANT", "5.6")) + PixUseMCInfo("RichSIMDPixelsUseMCInfo").InnerPixelQuantization = innerPixQ + PixUseMCInfo("RichSIMDPixelsUseMCInfo").OuterPixelQuantization = outerPixQ + if usePixelMCPos or usePixelMCTime or useMCHits: fetcher.DataKeys += ['pSim/Rich/Hits'] preUnpackSeq.Members += [MCRichHitUnpacker()] @@ -217,10 +224,11 @@ else: richMCDecode.RejectAllBackgrounds = False richMCDecode.RejectScintillation = False richMCDecode.IncludeTimeInfo = True - richMCDecode.AllowMultipleHits = False + richMCDecode.AllowMultipleHits = True + richMCDecode.MaxHitsPerChannel = (99999, 99999) all.Members += [richMCDecode] -# Now get the RICH sequence +# Now get the RICH sequence if filterTksUsingMC or useMCTracks: # Build track info from RICH extended MC info @@ -272,27 +280,27 @@ if not useMCTracks: # Input tracks #tkLocs = { - # "Long": tkFilt.LongTracks, - # "Down": tkFilt.DownstreamTracks, - # "Up": tkFilt.UpstreamTracks, - # "Seed": tkFilt.Ttracks + #"Long" : tkFilt.LongTracks, + #"Down" : tkFilt.DownstreamTracks, + #"Up" : tkFilt.UpstreamTracks, + #"Seed" : tkFilt.Ttracks #} - #tkLocs = { "All" : tkFilt.InputTracks } + #tkLocs = {"All" : tkFilt.InputTracks } tkLocs = {"Long": tkFilt.LongTracks} - #tkLocs = { "Up" : tkFilt.UpstreamTracks } - #tkLocs = { "Down" : tkFilt.DownstreamTracks } + #tkLocs = {"Up" : tkFilt.UpstreamTracks } + #tkLocs = {"Down" : tkFilt.DownstreamTracks } # Output PID #pidLocs = { - # "Long": "Rec/Rich/LongPIDs", - # "Down": "Rec/Rich/DownPIDs", - # "Up": "Rec/Rich/UpPIDs", - # "Seed": "Rec/Rich/SeedPIDs" + #"Long" : "Rec/Rich/LongPIDs", + #"Down" : "Rec/Rich/DownPIDs", + #"Up" : "Rec/Rich/UpPIDs", + #"Seed" : "Rec/Rich/SeedPIDs" #} - #pidLocs = { "All" : "Rec/Rich/PIDs" } + #pidLocs = {"All" : "Rec/Rich/PIDs" } pidLocs = {"Long": "Rec/Rich/LongPIDs"} - #pidLocs = { "Up" : "Rec/Rich/UpPIDs" } - #pidLocs = { "Down" : "Rec/Rich/DownPIDs" } + #pidLocs = {"Up" : "Rec/Rich/UpPIDs" } + #pidLocs = {"Down" : "Rec/Rich/DownPIDs" } else: @@ -335,7 +343,7 @@ rads = ["Rich1Gas", "Rich2Gas"] parts = [ "electron", "muon", "pion", "kaon", "proton", "deuteron", "belowThreshold" ] -#parts = ["pion","kaon","proton","deuteron","belowThreshold"] +#parts = ["pion", "kaon", "proton", "deuteron", "belowThreshold"] photonRecoType = "Quartic" #photonRecoType = "FastQuartic" @@ -345,13 +353,13 @@ photonSel = "Nominal" #photonSel = "Online" #photonSel = "None" -# Number ring points for ray tracing (scaled by CK theta) +# Number ring points for ray tracing( scaled by CK theta ) ringPointsMin = (16, 16, 16) ringPointsMax = (96, 96, 96) # Compute the ring share CK theta values -#rSTol = 0.075 # as fraction of sat. CK theta -#newCKRingTol = (rSTol, rSTol, rSTol) +#rSTol = 0.075 #as fraction of sat.CK theta +#newCKRingTol = ( rSTol, rSTol, rSTol ) newCKRingTol = (0.0, 0.05, 0.1) # Detectable Yield calculation @@ -368,7 +376,7 @@ WR2 = float(os.getenv("PIXBCKWEIGHTRICH2", ("1.0" if enable4D[0] else "1.0"))) pdBckWeights = [(WR1, WR2), (WR1, WR2), (WR1, WR2), (WR1, WR2)] # Hit treatment in background alg -# Need to try and understand the difference here between best 3D and 4D tunings +#Need to try and understand the difference here between best 3D and 4D tunings PDBackIgnoreHitData = [is4D, False, False, False] # background thresholds @@ -380,16 +388,16 @@ from RichFutureRecSys.ConfiguredRichReco import defaultNSigmaCuts ckResSigma = defaultNSigmaCuts() #S = 4.5 -#for tk in tkLocs.keys(): -# # Override presel cuts -# ckResSigma[photonSel][tkCKResTreatment][tk][0] = (S, S, S) -# # Override sel cuts -# ckResSigma[photonSel][tkCKResTreatment][tk][1] = (99,99,99) +# for tk in tkLocs.keys(): +# #Override presel cuts +#ckResSigma[photonSel][tkCKResTreatment][tk][0] = ( S, S, S ) +# #Override sel cuts +#ckResSigma[photonSel][tkCKResTreatment][tk][1] = ( 99, 99, 99 ) # Track CK res scale factors tkCKResScaleF = None # Use nominal values -#scF = 1.2 -#tkCKResScaleF = (scF, scF, scF) +#scF = 1.2 +#tkCKResScaleF = ( scF, scF, scF ) # Merged output finalPIDLoc = "Rec/Rich/PIDs" @@ -409,24 +417,24 @@ pixelClustering = (False, False) #pixelClustering = ( True, True ) # Truncate CK angles ? -#truncateAngles = (False, False, False) +#truncateAngles = ( False, False, False ) truncateAngles = (True, True, True) # Max PD occupancy -maxPDOcc = (int(os.getenv("MAXPDOCCR1", "64")), - int(os.getenv("MAXPDOCCR2", "64"))) +maxPDOcc = (int(os.getenv("MAXPDOCCR1", "99999")), + int(os.getenv("MAXPDOCCR2", "99999"))) # Minimum photon probability minPhotonProb = None # Nominal # Likelihood Settings -# Number Iterations +#Number Iterations nIts = 2 # likelihood thresholds likeThres = [-1e-2, -1e-3, -1e-4, -1e-5] # nominal #likeThres = [-1e-3, -1e-3, -1e-4, -1e-5] -# OutputLevel (-1 means no setting) +# OutputLevel( -1 means no setting ) outLevel = -1 # -------------------------------------------------------------------------------------- @@ -542,7 +550,7 @@ if enableMCChecks: if enablePIDTuples: from Configurables import NTupleSvc app.ExtSvc.append(NTupleSvc()) - # ProtoParticle monitoring. For detailed PID monitoring. + # ProtoParticle monitoring.For detailed PID monitoring. protoMoniSeq = GaudiSequencer("ProtoMonitoring", MeasureTime=True) all.Members += [protoMoniSeq] # Remake the protos with the new Rich PID objects @@ -570,17 +578,17 @@ if enablePIDTuples: GlobalRecoChecks().ANNTupleName = rootFileBaseName + "-ANNPIDTuple.root" # ANNPID monitoring GlobalRecoChecks().AddANNPIDInfo = False - #GlobalRecoChecks().AddANNPIDInfo = True + # GlobalRecoChecks().AddANNPIDInfo = True # -------------------------------------------------------------------------------------- # -------------------------------------------------------------------------------------- # Example command lines # -------------------------------------------------------------------------------------- -# Normal running (old data format) +# Normal running( old data format ) # gaudirun.py ~/LHCbCMake/Feature/Rec/Rich/RichFutureRecSys/examples/{RichFuture.py,data/PMTs/MCFormat/MCXDstFiles.py} 2>&1 | tee RichFuture-${User_release_area##/*/}-${CMTCONFIG}.log -# Normal running (realistic data format) +# Normal running (realistic data format) # gaudirun.py ~/LHCbCMake/Feature/Rec/Rich/RichFutureRecSys/examples/{RichFuture.py,data/PMTs/RealTel40Format/MCLDstFiles.py} 2>&1 | tee RichFuture-${User_release_area##/*/}-${CMTCONFIG}.log # Run on XSIM files -- GitLab From 7b88ef5670494ac49e8f110a1937b9cc7240394d Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Tue, 28 Nov 2023 10:49:35 +0000 Subject: [PATCH 08/23] GlobalReco: Small updates to scripts to make RICH performance plots --- Rec/GlobalReco/root/GlobalPID.C | 15 ++- Rec/GlobalReco/root/MakeRichPlots.C | 6 +- Rec/GlobalReco/root/RichKaonIDCompareFiles.C | 125 +++++++++++++------ 3 files changed, 102 insertions(+), 44 deletions(-) diff --git a/Rec/GlobalReco/root/GlobalPID.C b/Rec/GlobalReco/root/GlobalPID.C index 0cdcf3a66d9..e9702ae8991 100755 --- a/Rec/GlobalReco/root/GlobalPID.C +++ b/Rec/GlobalReco/root/GlobalPID.C @@ -606,11 +606,15 @@ void GlobalPID::saveFigures( const std::string& type ) { void GlobalPID::loadTTree( const std::string& filename ) { resetFile(); - std::cout << "Opening file : " << filename << std::endl; - m_file = std::make_unique( filename.c_str() ); - m_file->cd( ( filename + ":/ChargedProtoTuple" ).c_str() ); - fChain = (TTree*)gDirectory->Get( "protoPtuple" ); - Init( fChain ); + if ( boost::filesystem::exists( filename ) ) { + std::cout << "Opening file '" << filename << "'" << std::endl; + m_file = std::make_unique( filename.c_str() ); + m_file->cd( ( filename + ":/ChargedProtoTuple" ).c_str() ); + fChain = (TTree*)gDirectory->Get( "protoPtuple" ); + Init( fChain ); + } else { + std::cout << "File '" << filename << "' does not exist" << std::endl; + } } void GlobalPID::resetFile() { @@ -618,6 +622,7 @@ void GlobalPID::resetFile() { m_file->Close(); m_file.reset( nullptr ); } + fChain = nullptr; } GlobalPID::GlobalPID() {} diff --git a/Rec/GlobalReco/root/MakeRichPlots.C b/Rec/GlobalReco/root/MakeRichPlots.C index 2c824906eca..51c0e57ad13 100755 --- a/Rec/GlobalReco/root/MakeRichPlots.C +++ b/Rec/GlobalReco/root/MakeRichPlots.C @@ -31,9 +31,9 @@ void MakeRichPlots() { const std::string fName = "RichFuture-Feature-x86_64_v2-centos7-gcc12+detdesc-opt-Expert-ProtoTuple.root"; - const std::vector nametags = {"FromMCHits/3D/lumi-2.0e33", "FromMCHits/3D/lumi-1.2e34", - "FromMCHits/4D/lumi-2.0e33/PixWin-3.000/PhotWin-0.100", - "FromMCHits/4D/lumi-1.2e34/PixWin-3.000/PhotWin-0.100"}; + const std::vector nametags = {"FromMCHits-V3/3D/lumi-2.0e33", "FromMCHits-V3/3D/lumi-1.2e34", + "FromMCHits-V3/4D/lumi-2.0e33/PixWin-1.000/PhotWin-0.100", + "FromMCHits-V3/4D/lumi-1.2e34/PixWin-1.000/PhotWin-0.100"}; for ( auto nametag : nametags ) { diff --git a/Rec/GlobalReco/root/RichKaonIDCompareFiles.C b/Rec/GlobalReco/root/RichKaonIDCompareFiles.C index c13ba7cbb0d..72adca4e672 100755 --- a/Rec/GlobalReco/root/RichKaonIDCompareFiles.C +++ b/Rec/GlobalReco/root/RichKaonIDCompareFiles.C @@ -9,10 +9,14 @@ * or submit itself to any jurisdiction. * \*****************************************************************************/ +#include +#include + #include #include #include #include +#include #include #include "GlobalPID.C" @@ -30,45 +34,91 @@ void RichKaonIDCompareFiles() { // dataSets["Run3Tracking-3D"].push_back( "DSTs-MCTracks" ); // dataSets["Run3Tracking-3D"].push_back( "DSTs-RecoTracks" ); - // const std::string bName = "FromMCHits-V2"; - // for ( const std::string lumi : {"1.5e34", "1.2e34", "1.0e34", "3.0e33", "2.0e33", "2.0e32"} ) { - // dataSets[bName + "-3D"].push_back( bName + "/3D/lumi-" + lumi ); - // for ( const std::string pixWin : {"3.000"} ) { - // for ( const std::string photWin : {"1.000", "0.500", "0.250", "0.100", "0.050", "0.010"} ) { - // const auto n = bName + "/4D/lumi-" + lumi + "/PixWin-" + pixWin + "/PhotWin-" + photWin; - // dataSets[bName + "-4D-PixW_" + pixWin + "-PhoW_" + photWin].push_back( n ); - // dataSets[bName + "-4D-PixW_" + pixWin + "-Lumi_" + lumi].push_back( n ); + auto nA = []( const std::string& label, const std::string& value ) { + return ( value.empty() ? "" : label + "-" + value + "/" ); + }; + + auto add = [&]( const auto& tag, const auto& name ) { + const auto full_name = dir + "/" + name + "/" + fName; + if ( boost::filesystem::exists( full_name ) ) { + dataSets[tag].push_back( name ); + } else { + std::cout << "WARNING: '" << full_name << "' does not exist" << std::endl; + } + }; + + // for ( const std::string bName : {"PixSize-R1-V1", "PixSize-R2-V1", "PixSize-R1R2-V1"} ) { + // for ( const std::string lumi : {"1.2e34"} ) { + // for ( const std::string innerQ : {"1.00", "2.00", "3.00"} ) { + // for ( const std::string outerQ : {"", "2.00", "4.00", "6.00"} ) { + // const auto d1 = nA( "lumi", lumi ) + nA( "InnerPixQ", innerQ ) + nA( "OuterPixQ", outerQ ); + // add( bName + "/3D", bName + "/3D/" + d1 ); + // for ( const std::string pixWin : {"1.000"} ) { + // for ( const std::string photWin : {"0.300", "0.150", "0.075", "0.050"} ) { + // const auto data = bName + "/4D/" + d1 + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin ); + // add( bName + "/4D-PixW_" + pixWin + "-PhoW_" + photWin, data ); + // // add( bName + "/4D-PixW_" + pixWin + "-Lumi_" + lumi, data ); + // } + // } + // } // } // } // } - for ( const std::string lumi : {"1.2e34", "2.0e33", "2.0e32"} ) { - const std::string bName = "BkgTuneRich1"; - for ( const std::string w : {"0.600", "0.800", "1.000", "1.200", "1.400"} ) { - dataSets["3D/lumi-" + lumi + "/" + bName].push_back( bName + "/3D/lumi-" + lumi + "/R1PixBckW-" + w ); + const std::string batch = "FromMCHits-V5"; + const std::string innerQ = "2.80"; + const std::string outerQ = "5.60"; + const bool addEmulated = true; + const auto QS = nA( "InnerPixQ", innerQ ) + nA( "OuterPixQ", outerQ ); + for ( const std::string vari : {"Nominal", "MultiHits", "MultiHits-PerfectMCPos"} ) { + const auto bName = batch + "-" + vari; + for ( const std::string lumi : {"1.5e34", "1.2e34", "1.0e34", "3.0e33", "2.0e33"} ) { + add( vari + "/3D", bName + "/3D/" + nA( "lumi", lumi ) ); + if ( addEmulated ) { add( vari + "/3D", "PixSize-R1R2-V1/3D/" + nA( "lumi", lumi ) + QS ); } for ( const std::string pixWin : {"1.000"} ) { - for ( const std::string photWin : {"0.250", "0.100", "0.050"} ) { - const auto name = "4D/lumi-" + lumi + "/PixWin-" + pixWin + "/" + bName + "_PhotWin-" + photWin; - dataSets[name].push_back( bName + "/4D/lumi-" + lumi + "/R1PixBckW-" + w + "/PixWin-" + pixWin + "/PhotWin-" + - photWin ); + for ( const std::string photWin : {"0.300", "0.150", "0.075", "0.050"} ) { + auto data = bName + "/4D/" + nA( "lumi", lumi ) + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin ); + add( vari + "/4D-PixW_" + pixWin + "-PhoW_" + photWin, data ); + add( vari + "/4D-PixW_" + pixWin + "-Lumi_" + lumi, data ); + if ( addEmulated ) { + data = "PixSize-R1R2-V1/4D/" + nA( "lumi", lumi ) + QS + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin ); + add( vari + "/4D-PixW_" + pixWin + "-PhoW_" + photWin, data ); + add( vari + "/4D-PixW_" + pixWin + "-Lumi_" + lumi, data ); + } } } } } - for ( const std::string lumi : {"1.2e34", "2.0e33", "2.0e32"} ) { - const std::string bName = "BkgTuneRich2"; - for ( const std::string w : {"0.600", "0.800", "1.000", "1.200", "1.400"} ) { - dataSets["3D/lumi-" + lumi + "/" + bName].push_back( bName + "/3D/lumi-" + lumi + "/R2PixBckW-" + w ); - for ( const std::string pixWin : {"1.000"} ) { - for ( const std::string photWin : {"0.250", "0.100", "0.050"} ) { - const auto name = "4D/lumi-" + lumi + "/PixWin-" + pixWin + "/" + bName + "_PhotWin-" + photWin; - dataSets[name].push_back( bName + "/4D/lumi-" + lumi + "/R2PixBckW-" + w + "/PixWin-" + pixWin + "/PhotWin-" + - photWin ); - } - } - } - } + // for ( const std::string lumi : {"1.2e34", "2.0e33", "2.0e32"} ) { + // const std::string bName = "BkgTuneRich1"; + // for ( const std::string w : {"0.600", "0.800", "1.000", "1.200", "1.400"} ) { + // dataSets["3D/lumi-" + lumi + "/" + bName].push_back( bName + "/3D/lumi-" + lumi + "/R1PixBckW-" + w ); + // for ( const std::string pixWin : {"1.000"} ) { + // for ( const std::string photWin : {"0.250", "0.100", "0.050"} ) { + // const auto name = "4D/lumi-" + lumi + "/PixWin-" + pixWin + "/" + bName + "_PhotWin-" + photWin; + // dataSets[name].push_back( bName + "/4D/lumi-" + lumi + "/R1PixBckW-" + w + "/PixWin-" + pixWin + + // "/PhotWin-" + + // photWin ); + // } + // } + // } + // } + + // for ( const std::string lumi : {"1.2e34", "2.0e33", "2.0e32"} ) { + // const std::string bName = "BkgTuneRich2"; + // for ( const std::string w : {"0.600", "0.800", "1.000", "1.200", "1.400"} ) { + // dataSets["3D/lumi-" + lumi + "/" + bName].push_back( bName + "/3D/lumi-" + lumi + "/R2PixBckW-" + w ); + // for ( const std::string pixWin : {"1.000"} ) { + // for ( const std::string photWin : {"0.250", "0.100", "0.050"} ) { + // const auto name = "4D/lumi-" + lumi + "/PixWin-" + pixWin + "/" + bName + "_PhotWin-" + photWin; + // dataSets[name].push_back( bName + "/4D/lumi-" + lumi + "/R2PixBckW-" + w + "/PixWin-" + pixWin + + // "/PhotWin-" + + // photWin ); + // } + // } + // } + // } const Long64_t nTracks = 1e6; @@ -93,8 +143,8 @@ void RichKaonIDCompareFiles() { gConf.useFixedGraphRange = true; gConf.minGraphX = 80; gConf.maxGraphX = 100; - gConf.minGraphY = 0.05; - gConf.maxGraphY = 100; + gConf.minGraphY = 0.08; + gConf.maxGraphY = 60; } else if ( GlobalPID::Upstream == tktype ) { gConf.useFixedGraphRange = true; gConf.minGraphX = 40; @@ -155,17 +205,20 @@ void RichKaonIDCompareFiles() { gConf.writeCutValues = false; - using PlotData = std::vector>; + using PlotData = std::vector>; // colours ... https://root.cern.ch/doc/master/classTColor.html - const std::array colors{kYellow - 1, kRed + 1, kGreen + 2, kBlue + 1, kMagenta + 2, kBlack, - kRed - 6, kBlue - 1, kCyan + 2, kGreen - 5, kGray + 2}; + const auto colors = + std::array{kRed + 1, kBlue + 1, kGreen + 2, kMagenta + 2, kBlack + 0, kOrange + 7, kCyan - 3, + kMagenta - 7, kGreen - 8, kRed - 6, kCyan + 2, kGray + 2, kAzure + 8}; auto lastColor = colors.begin(); auto genTuple = [&]( const auto& tag ) { if ( colors.end() == lastColor ) { lastColor = colors.begin(); } - return std::make_tuple( dir + "/" + tag + "/" + fName, tag, *( lastColor++ ) ); + auto fullname = dir + "/" + tag + "/" + fName; + boost::replace_all( fullname, "//", "/" ); + return std::make_tuple( std::move( fullname ), tag, *( lastColor++ ) ); }; PlotData plotdata; @@ -180,7 +233,7 @@ void RichKaonIDCompareFiles() { pid->config.title += " | " + dname; pid->config.subtitle = title; pid->config.superImpose = ( iPlot++ != 0 ); - pid->config.color = color; + pid->config.color = (Color_t)color; // create the plot pid->makeCurve( nTracks ); } -- GitLab From 32c87b9d4cadb609471de93bcb0d91be630667de Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Tue, 28 Nov 2023 10:49:56 +0000 Subject: [PATCH 09/23] Add support for emulated pixel sizes to RICH example job submission scripts --- .../RichFutureRecSys/examples/jobs/RunJobs.py | 105 +++++++++++------- Rich/RichFutureRecSys/examples/jobs/submit.sh | 4 +- 2 files changed, 69 insertions(+), 40 deletions(-) diff --git a/Rich/RichFutureRecSys/examples/jobs/RunJobs.py b/Rich/RichFutureRecSys/examples/jobs/RunJobs.py index 0d654bca92b..1548abd88c3 100755 --- a/Rich/RichFutureRecSys/examples/jobs/RunJobs.py +++ b/Rich/RichFutureRecSys/examples/jobs/RunJobs.py @@ -30,7 +30,7 @@ parser.add_argument( parser.add_argument( "--photwins", help="Photon time window(s) (only for 4D)", - default="1.0,0.5,0.25,0.1,0.05,0.01", + default="0.300,0.150,0.075,0.050", type=str) parser.add_argument( "--pixbckwR1", @@ -42,6 +42,10 @@ parser.add_argument( help="Pixel background weights for RICH2", default="", type=str) +parser.add_argument( + "--innerPixQ", help="Inner Pixel Quantisation", default="", type=str) +parser.add_argument( + "--outerPixQ", help="Outer Pixel Quantisation", default="", type=str) parser.add_argument("-N", "--name", help="Job Name", default="Test", type=str) args = parser.parse_args() config = vars(args) @@ -76,13 +80,66 @@ for Lumi in config["lumis"].split(","): os.environ["LUMI"] = Lumi - for r1PixW in config["pixbckwR1"].split(","): - for r2PixW in config["pixbckwR2"].split(","): - if enable4D: - for PixWin in config["pixwins"].split(","): - for PhotWin in config["photwins"].split(","): - - jobName = config["name"] + "/4D/lumi-" + Lumi + for innerPixQ in config["innerPixQ"].split(","): + for outerPixQ in config["outerPixQ"].split(","): + + for r1PixW in config["pixbckwR1"].split(","): + for r2PixW in config["pixbckwR2"].split(","): + if enable4D: + for PixWin in config["pixwins"].split(","): + for PhotWin in config["photwins"].split(","): + + jobName = config["name"] + "/4D/lumi-" + Lumi + if len(innerPixQ) != 0: + jobName += "/InnerPixQ-" + "{:.2f}".format( + float(innerPixQ)) + os.environ["INNERPIXQUANT"] = innerPixQ + if len(outerPixQ) != 0: + jobName += "/OuterPixQ-" + "{:.2f}".format( + float(outerPixQ)) + os.environ["OUTERPIXQUANT"] = outerPixQ + if len(r1PixW) != 0: + jobName += "/R1PixBckW-" + "{:.3f}".format( + float(r1PixW)) + os.environ["PIXBCKWEIGHTRICH1"] = r1PixW + if len(r2PixW) != 0: + jobName += "/R2PixBckW-" + "{:.3f}".format( + float(r2PixW)) + os.environ["PIXBCKWEIGHTRICH2"] = r2PixW + if len(PixWin) != 0: + jobName += "/PixWin-" + "{:.3f}".format( + float(PixWin)) + os.environ["PIXWIN"] = PixWin + if len(PhotWin) != 0: + jobName += "/PhotWin-" + "{:.3f}".format( + float(PhotWin)) + os.environ["PHOTWIN"] = PhotWin + + print("Submitting jobs for", jobName) + subprocess.run([ + optsroot + "jobs/submit.sh", jobName, + nFiles(float(Lumi)), + nFilesPerJob(float(Lumi)) + ]) + + unsetenv("INNERPIXQUANT") + unsetenv("OUTERPIXQUANT") + unsetenv("PIXBCKWEIGHTRICH1") + unsetenv("PIXBCKWEIGHTRICH2") + unsetenv("PIXWIN") + unsetenv("PHOTWIN") + + else: + + jobName = config["name"] + "/3D/lumi-" + Lumi + if len(innerPixQ) != 0: + jobName += "/InnerPixQ-" + "{:.2f}".format( + float(innerPixQ)) + os.environ["INNERPIXQUANT"] = innerPixQ + if len(outerPixQ) != 0: + jobName += "/OuterPixQ-" + "{:.2f}".format( + float(outerPixQ)) + os.environ["OUTERPIXQUANT"] = outerPixQ if len(r1PixW) != 0: jobName += "/R1PixBckW-" + "{:.3f}".format( float(r1PixW)) @@ -91,14 +148,6 @@ for Lumi in config["lumis"].split(","): jobName += "/R2PixBckW-" + "{:.3f}".format( float(r2PixW)) os.environ["PIXBCKWEIGHTRICH2"] = r2PixW - if len(PixWin) != 0: - jobName += "/PixWin-" + "{:.3f}".format( - float(PixWin)) - os.environ["PIXWIN"] = PixWin - if len(PhotWin) != 0: - jobName += "/PhotWin-" + "{:.3f}".format( - float(PhotWin)) - os.environ["PHOTWIN"] = PhotWin print("Submitting jobs for", jobName) subprocess.run([ @@ -107,29 +156,9 @@ for Lumi in config["lumis"].split(","): nFilesPerJob(float(Lumi)) ]) + unsetenv("INNERPIXQUANT") + unsetenv("OUTERPIXQUANT") unsetenv("PIXBCKWEIGHTRICH1") unsetenv("PIXBCKWEIGHTRICH2") - unsetenv("PIXWIN") - unsetenv("PHOTWIN") - - else: - - jobName = config["name"] + "/3D/lumi-" + Lumi - if len(r1PixW) != 0: - jobName += "/R1PixBckW-" + "{:.3f}".format(float(r1PixW)) - os.environ["PIXBCKWEIGHTRICH1"] = r1PixW - if len(r2PixW) != 0: - jobName += "/R2PixBckW-" + "{:.3f}".format(float(r2PixW)) - os.environ["PIXBCKWEIGHTRICH2"] = r2PixW - - print("Submitting jobs for", jobName) - subprocess.run([ - optsroot + "jobs/submit.sh", jobName, - nFiles(float(Lumi)), - nFilesPerJob(float(Lumi)) - ]) - - unsetenv("PIXBCKWEIGHTRICH1") - unsetenv("PIXBCKWEIGHTRICH2") unsetenv("LUMI") diff --git a/Rich/RichFutureRecSys/examples/jobs/submit.sh b/Rich/RichFutureRecSys/examples/jobs/submit.sh index 35f3225d554..91b7a33ec0d 100755 --- a/Rich/RichFutureRecSys/examples/jobs/submit.sh +++ b/Rich/RichFutureRecSys/examples/jobs/submit.sh @@ -111,8 +111,8 @@ environment = CONDOR_ID=\$(Cluster).\$(Process) JobBatchName = ${JobName//} # Requirements -#Requirements = (POOL != "GENERAL" && HAS_r02 && OpSysAndVer == "CentOS7" && TARGET.has_avx2) -Requirements = (POOL != "GENERAL" && HAS_r02 && OpSysAndVer == "CentOS7") +#Requirements = (POOL != "GENERAL" && HAS_r01 && HAS_r02 && OpSysAndVer == "CentOS7" && TARGET.has_avx2) +Requirements = (POOL != "GENERAL" && HAS_r01 && HAS_r02 && OpSysAndVer == "CentOS7") # Rank hosts according to floating point speed Rank = kflops -- GitLab From 4aea3096f7aee4087fcf2ce3daa4f486f67b7079 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Tue, 5 Dec 2023 14:45:52 +0000 Subject: [PATCH 10/23] RichSIMDPhotonPredictedPixelSignal: Remove unused includes --- .../src/RichSIMDPhotonPredictedPixelSignal.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp index 2f32d8c00ee..a740a34ebdf 100644 --- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp +++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp @@ -37,9 +37,7 @@ // STL #include -#include #include -#include namespace Rich::Future::Rec { -- GitLab From 8c81e4b613ae72fa7d778aaf12dac77cd0112cd1 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Mon, 11 Dec 2023 11:35:24 +0000 Subject: [PATCH 11/23] RichFutureRecMCAlgorithms: Add two new MC cheating algorithms for photon information --- Rich/RichFutureRecMCAlgorithms/CMakeLists.txt | 2 + .../src/RichPhotonUseMCInfo.cpp | 216 ++++++++++++++++++ .../src/RichTrackCKAnglesUseMCInfo.cpp | 151 ++++++++++++ .../RichFutureRecSys/ConfiguredRichReco.py | 68 +++++- 4 files changed, 432 insertions(+), 5 deletions(-) create mode 100644 Rich/RichFutureRecMCAlgorithms/src/RichPhotonUseMCInfo.cpp create mode 100644 Rich/RichFutureRecMCAlgorithms/src/RichTrackCKAnglesUseMCInfo.cpp diff --git a/Rich/RichFutureRecMCAlgorithms/CMakeLists.txt b/Rich/RichFutureRecMCAlgorithms/CMakeLists.txt index 10067b4cb6c..42548752b43 100644 --- a/Rich/RichFutureRecMCAlgorithms/CMakeLists.txt +++ b/Rich/RichFutureRecMCAlgorithms/CMakeLists.txt @@ -18,6 +18,8 @@ gaudi_add_module(RichFutureRecMCAlgorithms src/RichTrSegMakerFromMCRichTracks.cpp src/RichSegmentAddTimeFromMC.cpp src/RichPixelUseMCInfo.cpp + src/RichPhotonUseMCInfo.cpp + src/RichTrackCKAnglesUseMCInfo.cpp src/RichMCTrackFilter.cpp LINK Gaudi::GaudiAlgLib diff --git a/Rich/RichFutureRecMCAlgorithms/src/RichPhotonUseMCInfo.cpp b/Rich/RichFutureRecMCAlgorithms/src/RichPhotonUseMCInfo.cpp new file mode 100644 index 00000000000..99d58111233 --- /dev/null +++ b/Rich/RichFutureRecMCAlgorithms/src/RichPhotonUseMCInfo.cpp @@ -0,0 +1,216 @@ +/*****************************************************************************\ +* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration * +* * +* This software is distributed under the terms of the GNU General Public * +* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". * +* * +* In applying this licence, CERN does not waive the privileges and immunities * +* granted to it by virtue of its status as an Intergovernmental Organization * +* or submit itself to any jurisdiction. * +\*****************************************************************************/ + +// Gaudi +#include "GaudiKernel/ParsersFactory.h" +#include "GaudiKernel/PhysicalConstants.h" +#include "GaudiKernel/RndmGenerators.h" +#include "GaudiKernel/StdArrayAsProperty.h" + +// base class +#include "RichFutureRecBase/RichRecAlgBase.h" + +// Gaudi Functional +#include "LHCbAlgs/Transformer.h" + +// Event Model +#include "Event/MCRichHit.h" +#include "Event/MCRichOpticalPhoton.h" +#include "Event/Track.h" +#include "RichFutureRecEvent/RichRecCherenkovAngles.h" +#include "RichFutureRecEvent/RichRecCherenkovPhotons.h" +#include "RichFutureRecEvent/RichRecRelations.h" +#include "RichFutureRecEvent/RichRecSIMDPixels.h" + +// MC Relations +#include "RichFutureMCUtils/RichMCHitUtils.h" +#include "RichFutureMCUtils/RichMCOpticalPhotonUtils.h" +#include "RichFutureMCUtils/RichRecMCHelper.h" + +// Rich Utils +#include "RichUtils/ZipRange.h" + +// ROOT +#include + +// STD +#include +#include +#include + +namespace Rich::Future::Rec::MC { + + namespace { + using Data = SIMDCherenkovPhoton::Vector; + } + + /** @class RichPhotonUseMCInfo RichPhotonUseMCInfo.h + * + * Use MC truth to cheat RICH Cherenkov photons + * + * @author Chris Jones + * @date 2023-09-20 + */ + + class RichPhotonUseMCInfo final + : public LHCb::Algorithm::Transformer>> { + + public: + /// Standard constructor + RichPhotonUseMCInfo( const std::string& name, ISvcLocator* pSvcLocator ) + : Transformer( name, pSvcLocator, + // data inputs + {KeyValue{"InCherenkovPhotonLocation", SIMDCherenkovPhotonLocation::Default}, + KeyValue{"RichSIMDPixelSummariesLocation", SIMDPixelSummariesLocation::Default}, + KeyValue{"TracksLocation", LHCb::TrackLocation::Default}, + KeyValue{"SignalCherenkovAnglesLocation", CherenkovAnglesLocation::Signal}, + KeyValue{"PhotonToParentsLocation", Relations::PhotonToParentsLocation::Default}, + KeyValue{"SegmentToTrackLocation", Relations::SegmentToTrackLocation::Default}, + KeyValue{"TrackToMCParticlesRelations", Rich::Future::MC::Relations::TrackToMCParticles}, + KeyValue{"RichDigitSummariesLocation", LHCb::MCRichDigitSummaryLocation::Default}, + KeyValue{"MCRichHitsLocation", LHCb::MCRichHitLocation::Default}, + KeyValue{"MCRichOpticalPhotonsLocation", LHCb::MCRichOpticalPhotonLocation::Default}}, + // output data + {KeyValue{"OutCherenkovPhotonLocation", SIMDCherenkovPhotonLocation::Default + "Out"}} ) { + // setProperty( "OutputLevel", MSG::VERBOSE ).ignore(); + } + + /// Initialize + StatusCode initialize() override { + return Transformer::initialize().andThen( [&] { + if ( !m_useMcPhotCKT ) { + for ( const auto rich : activeDetectors() ) { + info() << rich << " CK theta resolution func: " << m_ckThetaSmearVP[rich] << endmsg; + m_ckVP[rich].emplace( Rich::text( rich ).c_str(), m_ckThetaSmearVP[rich].c_str(), 0.0, 120000.0 ); + } + m_unitGauss = Rndm::Numbers( randSvc(), Rndm::Gauss( 0.0, 1.0 ) ); + } + } ); + } + + /// Functional operator + Data operator()( const Data& photons, // + const SIMDPixelSummaries& pixels, // + const LHCb::Track::Range& tracks, // + const CherenkovAngles::Vector& ckAngles, // + const Relations::PhotonToParents::Vector& photRels, // + const Relations::SegmentToTrackVector& segToTkRel, // + const Rich::Future::MC::Relations::TkToMCPRels& tkmcrels, // + const LHCb::MCRichDigitSummarys& digitSums, // + const LHCb::MCRichHits& mchits, // + const LHCb::MCRichOpticalPhotons& mcphotons ) const override { + + // Make local MC helper objects + const Helper mcHelper( tkmcrels, digitSums ); + const MCHitUtils hitHelper( mchits ); + const MCOpticalPhotonUtils photHelper( mcphotons ); + + // Clone input photons before modifying them + auto new_photons = photons; + + // Loop over new photons to update them using MC + for ( auto&& [phot, rels] : Ranges::Zip( new_photons, photRels ) ) { + + // Get the track for this photon + const auto tkIn = segToTkRel[rels.segmentIndex()]; + const auto* track = tracks[tkIn]; + // Get pixel for this photon + const auto& pixel = pixels[rels.pixelIndex()]; + + // Get the MCParticle(s) for this track + const auto mcPs = mcHelper.mcParticles( *track, false, 0.75 ); + if ( mcPs.empty() ) { continue; } + const auto mcP = mcPs.front(); + + // Get true PID type from first MCParticle + const auto pid = mcHelper.mcParticleType( mcP ); + if ( Rich::Unknown == pid ) { continue; } + + // Loop over scalar entries in SIMD photon + for ( std::size_t i = 0; i < SIMDCherenkovPhoton::SIMDFP::Size; ++i ) { + // Select valid entries + if ( !phot.validityMask()[i] ) { continue; } + + // ID for this hit + const auto id = pixel.smartID()[i]; + + // get the MCHits for this ID + const auto mcHits = hitHelper.mcRichHits( id ); + // ... and then the MC Optical Photons also linked to this track's MCParticle(s) + const auto mcPhots = photHelper.mcOpticalPhotons( mcHits, mcPs ); + _ri_verbo << "Found " << mcPhots.size() << " MC photons" << endmsg; + if ( !mcPhots.empty() ) { + // This is a true CK photon so update the photon angles + if ( m_useMcPhotCKT ) { + // Use Cherenkov angles from first MCPhoton. Should really only be 1 + const auto mcPhot = mcPhots.front(); + if ( mcPhot ) { // should never fail but just in case ... + phot.CherenkovTheta()[i] = mcPhot->cherenkovTheta(); + phot.CherenkovPhi()[i] = mcPhot->cherenkovPhi(); + _ri_verbo << " -> Updated from MCPhoton " << phot.CherenkovTheta()[i] << " " << phot.CherenkovPhi()[i] + << endmsg; + } + } else { + // Use expected values for true PID type with smearing + const auto ckTheta = ckAngles[rels.segmentIndex()][pid]; + // resolution for this momentum + const auto ptot = mcP->p(); + const auto ckRes = m_ckVP[phot.rich()]->Eval( ptot ); + // New CK theta value with smearing + const auto smeared_ckT = ckTheta + ( m_unitGauss.shoot() * ckRes ); + // update the value in the photon + _ri_verbo << " -> Updated from TK expected: " << phot.rich() << " theta=" << ckTheta << " Ptot=" << ptot + << " CKRes=" << ckRes << " | New CKtheta=" << smeared_ckT << endmsg; + phot.CherenkovTheta()[i] = smeared_ckT; + } + } + } + } + + return new_photons; + } + + private: + // data + + /// Unit Gaussian for CK theta smearing + Rndm::Numbers m_unitGauss{}; + + // TF1 functions for CK theta smear value as a function of P + DetectorArray> m_ckVP{{}}; + + private: + // properties + + /// Use true MC optical photon values + Gaudi::Property m_useMcPhotCKT{this, "UseMCPhotonCKThetaValues", false}; + + /// CK theta smearing function versus P for each RICH + Gaudi::Property> m_ckThetaSmearVP{ + this, + "CKThetaSmearFunc", + {"0.00078+0.0012*(std::tanh(-x/5000.0)+1.0)", "0.00065+0.0011*(std::tanh(-x/5000.0)+1.0)"}}; + }; + + // Declaration of the Algorithm Factory + DECLARE_COMPONENT( RichPhotonUseMCInfo ) + +} // namespace Rich::Future::Rec::MC diff --git a/Rich/RichFutureRecMCAlgorithms/src/RichTrackCKAnglesUseMCInfo.cpp b/Rich/RichFutureRecMCAlgorithms/src/RichTrackCKAnglesUseMCInfo.cpp new file mode 100644 index 00000000000..ff4680e5c16 --- /dev/null +++ b/Rich/RichFutureRecMCAlgorithms/src/RichTrackCKAnglesUseMCInfo.cpp @@ -0,0 +1,151 @@ +/*****************************************************************************\ +* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration * +* * +* This software is distributed under the terms of the GNU General Public * +* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". * +* * +* In applying this licence, CERN does not waive the privileges and immunities * +* granted to it by virtue of its status as an Intergovernmental Organization * +* or submit itself to any jurisdiction. * +\*****************************************************************************/ + +// Gaudi +#include "GaudiKernel/ParsersFactory.h" +#include "GaudiKernel/PhysicalConstants.h" +#include "GaudiKernel/StdArrayAsProperty.h" + +// base class +#include "RichFutureRecBase/RichRecAlgBase.h" + +// Gaudi Functional +#include "LHCbAlgs/Transformer.h" + +// Event Model +#include "Event/MCParticle.h" +#include "Event/MCRichHit.h" +#include "Event/MCRichOpticalPhoton.h" +#include "Event/MCRichSegment.h" +#include "Event/MCRichTrack.h" +#include "Event/Track.h" +#include "RichFutureRecEvent/RichRecCherenkovAngles.h" +#include "RichFutureRecEvent/RichRecRelations.h" + +// MC Relations +#include "RichFutureMCUtils/TrackToMCParticle.h" + +// Rich Utils +#include "RichUtils/ZipRange.h" + +// STD +#include +#include + +namespace Rich::Future::Rec::MC { + + namespace { + using Data = CherenkovAngles::Vector; + } + + /** @class RichTrackCKAnglesUseMCInfo RichTrackCKAnglesUseMCInfo.h + * + * Use MC truth to cheat RICH track Cherenkov Angles + * + * @author Chris Jones + * @date 2023-09-20 + */ + + class RichTrackCKAnglesUseMCInfo final + : public LHCb::Algorithm::Transformer>> { + + public: + /// Standard constructor + RichTrackCKAnglesUseMCInfo( const std::string& name, ISvcLocator* pSvcLocator ) + : Transformer( name, pSvcLocator, + // data inputs + {KeyValue{"InCherenkovAnglesLocation", CherenkovAnglesLocation::Signal}, + KeyValue{"TracksLocation", LHCb::TrackLocation::Default}, + KeyValue{"TrackSegmentsLocation", LHCb::RichTrackSegmentLocation::Default}, + KeyValue{"SegmentToTrackLocation", Relations::SegmentToTrackLocation::Default}, + KeyValue{"TrackToMCParticlesRelations", Rich::Future::MC::Relations::TrackToMCParticles}, + KeyValue{"MCRichTracksLocation", LHCb::MCRichTrackLocation::Default}}, + // output data + {KeyValue{"OutCherenkovAnglesLocation", CherenkovAnglesLocation::Signal + "Out"}} ) { + // setProperty( "OutputLevel", MSG::VERBOSE ).ignore(); + } + + public: + /// Functional operator + Data operator()( const Data& in_angles, // + const LHCb::Track::Range& tracks, // + const LHCb::RichTrackSegment::Vector& segments, // + const Relations::SegmentToTrackVector& segToTkRel, // + const Rich::Future::MC::Relations::TkToMCPRels& tkmcrels, // + const LHCb::MCRichTracks& mcRichTks ) const override { + + // MC utils + const Rich::Future::MC::Relations::TrackToMCParticle mcHelper( tkmcrels ); + + // create map from MCP -> MCRichTrack + std::unordered_map mcpToMCRT; + for ( const auto mcRichTk : mcRichTks ) { + const auto mcp = ( mcRichTk ? mcRichTk->mcParticle() : nullptr ); + if ( mcp ) { mcpToMCRT[mcp] = mcRichTk; } + } + + // clone input before modifying + auto new_angles = in_angles; + + // loop over angles to modify + for ( auto&& [segment, angles, tkIn] : Ranges::Zip( segments, new_angles, segToTkRel ) ) { + + // associated track object + const auto* track = tracks[tkIn]; + + // Get the MCParticle for this track + const auto mcPs = mcHelper.mcParticles( *track, false, 0.75 ); + if ( mcPs.empty() ) { continue; } + // just use the first + const auto mcP = mcPs.front(); + + // The True MCParticle type + const auto pid = mcHelper.mcParticleType( mcP ); + if ( Rich::Unknown == pid ) { continue; } + + // Get the MC Rich Track if available + const auto iMCRT = mcpToMCRT.find( mcP ); + if ( iMCRT == mcpToMCRT.end() ) { continue; } + const auto mcRT = iMCRT->second; + // then the MC segment + const auto mcSeg = ( mcRT ? mcRT->segmentInRad( segment.radiator() ) : nullptr ); + if ( !mcSeg ) { continue; } + // .. finally the MC photons + const auto& mcPhots = mcSeg->mcRichOpticalPhotons(); + if ( mcPhots.empty() ) { continue; } + + // Loop over photons and form average CK theta for signal + float theta{0}; + for ( const auto& mcPhot : mcPhots ) { theta += mcPhot->cherenkovTheta(); } + theta /= (float)mcPhots.size(); + + // Update the expected value for true type + _ri_verbo << "Updating " << angles[pid] << " PID=" << pid << " to " << theta << endmsg; + angles[pid] = theta; + } + + return new_angles; + } + + private: + // properties + }; + + // Declaration of the Algorithm Factory + DECLARE_COMPONENT( RichTrackCKAnglesUseMCInfo ) + +} // namespace Rich::Future::Rec::MC diff --git a/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py b/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py index 9efec26fd88..23b9441978b 100644 --- a/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py +++ b/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py @@ -371,6 +371,9 @@ def RichRecoSequence( # Add 4D origin vertex info to track segments using MC tkSegAddTimeFromMC=False, + # Use MC to update expected CK theta values for each segment + tkSegUseMCExpCKTheta=False, + #=========================================================== # Photon treatment options #=========================================================== @@ -404,6 +407,15 @@ def RichRecoSequence( # Only used if 4D reconstruction is active photonTimeWindow=(0.25, 0.25), + # MC cheating options + + # Use MC to cheat photon parameters + photUseMC=False, + + # Function to use when emulating a given CK resolution for each RICH + photEmulatedResMC=("0.00078+0.0012*(std::tanh(-x/5000.0)+1.0)", + "0.00065+0.0011*(std::tanh(-x/5000.0)+1.0)"), + #=========================================================== # Settings for the global PID minimisation #=========================================================== @@ -691,6 +703,15 @@ def RichRecoSequence( segCr.MCParticlesLinkLocation = "Link/" + trackLocation segs.Members += [segCr] + if photUseMC or tkSegUseMCExpCKTheta: + from Configurables import Rich__Future__MC__TrackToMCParticleRelations as TkToMCPRels + mcRels = makeRichAlg(TkToMCPRels, "RichUseMCRels" + name, algprops, + tktype) + segs.Members += [mcRels] + mcRels.TrackToMCParticlesRelations = "Rec/Track/RichMCRels" + tktype + mcRels.TracksLocation = trackLocation + mcRels.MCParticlesLinkLocation = "Link/" + trackLocation + # Track global impact points on PD plane tkGloPnts = makeRichAlg(TrackPanelGlobalPoints, "RichTrackGloPoints" + name, algprops, tktype) @@ -819,9 +840,22 @@ def RichRecoSequence( sigChAngles.SignalPhotonYieldLocation = locs[ "SignalPhotonYieldLocation"] # Outputs - sigChAngles.SignalCherenkovAnglesLocation = locs[ - "SignalCherenkovAnglesLocation"] + sigLoc = locs["SignalCherenkovAnglesLocation"] + sigChAngles.SignalCherenkovAnglesLocation = sigLoc segs.Members += [sigChAngles] + if tkSegUseMCExpCKTheta: + from Configurables import Rich__Future__Rec__MC__RichTrackCKAnglesUseMCInfo as UseMCTKAngles + useMCAngles = makeRichAlg(UseMCTKAngles, + "RichSignalCKAnglesUseMC" + name, + algprops, tktype) + segs.Members += [useMCAngles] + sigChAngles.SignalCherenkovAnglesLocation = sigLoc + "BeforeMCCheat" + useMCAngles.InCherenkovAnglesLocation = sigLoc + "BeforeMCCheat" + useMCAngles.TracksLocation = trackLocation + useMCAngles.TrackSegmentsLocation = locs["TrackSegmentsLocation"] + useMCAngles.SegmentToTrackLocation = locs["SegmentToTrackLocation"] + useMCAngles.TrackToMCParticlesRelations = "Rec/Track/RichMCRels" + tktype + useMCAngles.OutCherenkovAnglesLocation = sigLoc # Track Resolutions tkRes = makeRichAlg(TrackCKResolutions, "RichCKResolutions" + name, @@ -848,6 +882,8 @@ def RichRecoSequence( # RICH Photon Treatment. # ================================================================== + photAlgs = [] + if photonReco == "Quartic" or photonReco == "FastQuartic": from Configurables import Rich__Future__Rec__SIMDQuarticPhotonReco as PhotonReco elif photonReco == "CKEstiFromRadius": @@ -859,6 +895,7 @@ def RichRecoSequence( # The photon reconstruction photReco = makeRichAlg(PhotonReco, "RichPhotonReco" + name, algprops, tktype) + photAlgs += [photReco] # Inputs photReco.TrackSegmentsLocation = locs["TrackSegmentsLocation"] photReco.CherenkovAnglesLocation = locs[ @@ -925,9 +962,32 @@ def RichRecoSequence( raise ValueError("Unknown photon reconstruction algorithm '" + photonReco + "'") + if photUseMC: + from Configurables import Rich__Future__Rec__MC__RichPhotonUseMCInfo as PhotonUseMCInfo + photMCInfo = makeRichAlg(PhotonUseMCInfo, + "RichPhotonUseMCInfo" + name, algprops, + tktype) + photAlgs += [photMCInfo] + pLoc = locs["CherenkovPhotonLocation"] + photMCInfo.OutCherenkovPhotonLocation = pLoc + photMCInfo.InCherenkovPhotonLocation = pLoc + "BeforeMC" + photReco.CherenkovPhotonLocation = pLoc + "BeforeMC" + photMCInfo.TrackToMCParticlesRelations = mcRels.TrackToMCParticlesRelations + photMCInfo.RichSIMDPixelSummariesLocation = cLocs[ + "RichSIMDPixelSummariesLocation"] + photMCInfo.TracksLocation = trackLocation + photMCInfo.PhotonToParentsLocation = locs[ + "PhotonToParentsLocation"] + photMCInfo.SegmentToTrackLocation = locs["SegmentToTrackLocation"] + photMCInfo.SignalCherenkovAnglesLocation = locs[ + "SignalCherenkovAnglesLocation"] + photMCInfo.CKThetaSmearFunc = photEmulatedResMC + photMCInfo.UseMCPhotonCKThetaValues = False + # Predicted pixel signals photPredSig = makeRichAlg( PhotonPredSignal, "RichPredPixelSignal" + name, algprops, tktype) + photAlgs += [photPredSig] # Input photPredSig.RichSIMDPixelSummariesLocation = cLocs[ "RichSIMDPixelSummariesLocation"] @@ -950,9 +1010,7 @@ def RichRecoSequence( # The photon reco sequence phots = GaudiSequencer( - "RichPhotons" + name, - MeasureTime=MeasureTime, - Members=[photReco, photPredSig]) + "RichPhotons" + name, MeasureTime=MeasureTime, Members=photAlgs) # Add to final sequence tkSeq.Members += [phots] -- GitLab From 86c8f9a26c41a525bbc47b64e0e9130bc31cc00c Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Mon, 11 Dec 2023 11:36:25 +0000 Subject: [PATCH 12/23] RichTrackCherenkovAngles: Tidy up emplace_back call --- .../src/RichTrackCherenkovAngles.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.cpp index 7c73d29d661..4cc93f6fea4 100644 --- a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.cpp +++ b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.cpp @@ -62,9 +62,9 @@ namespace Rich::Future::Rec { // iterate over segments and spectra together. for ( auto&& [segment, spectra, yields] : Ranges::ConstZip( segments, tkSpectra, tkYields ) ) { + // make a new object for the cherenkov angles - anglesV.emplace_back(); - auto& angles = anglesV.back(); + auto& angles = anglesV.emplace_back(); //_ri_verbo << "Tk" << endmsg; -- GitLab From 8eb473aa1d9866b660227b3080b72b9d6d46fd5b Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Mon, 11 Dec 2023 11:37:28 +0000 Subject: [PATCH 13/23] RichMCSIMDPhotonCherenkovAngles: Use fabs in res versus P plots --- .../src/RichMCSIMDPhotonCherenkovAngles.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Rich/RichFutureRecCheckers/src/RichMCSIMDPhotonCherenkovAngles.cpp b/Rich/RichFutureRecCheckers/src/RichMCSIMDPhotonCherenkovAngles.cpp index bbad7038385..c207ffd6c14 100644 --- a/Rich/RichFutureRecCheckers/src/RichMCSIMDPhotonCherenkovAngles.cpp +++ b/Rich/RichFutureRecCheckers/src/RichMCSIMDPhotonCherenkovAngles.cpp @@ -239,11 +239,11 @@ StatusCode SIMDPhotonCherenkovAngles::prebookHistograms() { "Reconstructed CKPhi - MC fake photons", // 0.0, 2.0 * Gaudi::Units::pi, nBins1D(), // "Cherenkov Phi / rad" ) ); - ok &= saveAndCheck( h_ckResTrueVP[rad], // - richProfile1D( HID( "ckResTrueVP", rad ), // - " V P - MC true photons - True Type", // - m_minP[rad], m_maxP[rad], nBins1D(), // - "Track Momentum (MeV/c)", " / rad" ) ); + ok &= saveAndCheck( h_ckResTrueVP[rad], // + richProfile1D( HID( "ckResTrueVP", rad ), // + "<|Rec-Exp CKTheta|> V P - MC true photons - True Type", // + m_minP[rad], m_maxP[rad], nBins1D(), // + "Track Momentum (MeV/c)", "<|Rec-Exp CKTheta|> / rad" ) ); // Plots assuming pion hypo @@ -430,7 +430,7 @@ void SIMDPhotonCherenkovAngles::operator()( const Summary::Track::Vector& if ( betaC ) { fillHisto( h_ckResTrueBetaCut[rad], deltaTheta, mcPW ); } fillHisto( h_thetaRecTrue[rad], thetaRec, mcPW ); fillHisto( h_phiRecTrue[rad], phiRec, mcPW ); - fillHisto( h_ckResTrueVP[rad], pTot, deltaTheta, mcPW ); + fillHisto( h_ckResTrueVP[rad], pTot, fabs( deltaTheta ), mcPW ); if ( Rich::Rich2Gas == rad ) { fillHisto( ( id.isHTypePMT() ? h_ckResTrueRich2HType : h_ckResTrueRich2RType ), deltaTheta, mcPW ); } -- GitLab From 3f7f8ff692b01b0f0989aa3c0329067e97bfe6ed Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Mon, 11 Dec 2023 11:38:03 +0000 Subject: [PATCH 14/23] Updat example RICH options --- Rich/RichFutureRecSys/examples/RichFuture.py | 39 ++++++++++++++------ 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py index 6e9e85825bd..d4b7b010d72 100644 --- a/Rich/RichFutureRecSys/examples/RichFuture.py +++ b/Rich/RichFutureRecSys/examples/RichFuture.py @@ -58,8 +58,8 @@ rootFileBaseName = "RichFuture-" + myBuild + "-" + myConf + "-" + histos HistogramPersistencySvc().OutputFile = rootFileBaseName + "-Histos.root" # Event numbers -nEvents = (100 if not batchMode else 99999999) -EventSelector().PrintFreq = (10 if not batchMode else 250) +nEvents = (1000 if not batchMode else 99999999) +EventSelector().PrintFreq = (100 if not batchMode else 250) #LHCbApp().SkipEvents = 2 # Just to initialise @@ -109,7 +109,7 @@ useMCTracks = True useMCHits = True # Get time windows from environment vars if set -# otherwise use default values as passed. +# otherwise use default values as passed. PixWin = float(os.getenv("PIXWIN", "1.000")) PhotWin = float(os.getenv("PHOTWIN", "0.075")) @@ -118,9 +118,9 @@ is4D = "ENABLE4D" in os.environ enable4D = (is4D, is4D) # Average expected time for signal avSignalTime = (13.03, 52.94) -# Pixel time window size in each RICH( in ns ) +# Pixel time window size in each RICH (ns) pixTimeWindow = (PixWin, PixWin) -# Photon reco time window size in each RICH( ns ) +# Photon reco time window size in each RICH (ns) photTimeWindow = (PhotWin, PhotWin) # Min Momentum cuts @@ -136,7 +136,13 @@ tkSegAddTimeFromMC = ((enable4D[0] or enable4D[1]) and not useMCTracks) usePixelMCTime = ((enable4D[0] or enable4D[1]) and not useMCHits) # Cheat pixel positions using MC info -usePixelMCPos = True +usePixelMCPos = False + +# Cheat photon parameters from MC +photUseMC = True + +# Cheat expected track CK theta valus from MC +tkSegUseMCExpCKTheta = False # Filter 'real' tracks based on MC info. # Designed to help match those made from MCRichTracks, where certain @@ -173,7 +179,7 @@ if UseDD4Hep: from Configurables import LHCb__Det__LbDD4hep__IOVProducer as IOVProducer all.Members += [IOVProducer("ReserveIOVDD4hep", ODIN=decodeODIN.ODIN)] -if tkSegAddTimeFromMC or useMCHits or useMCTracks or filterTksUsingMC: +if tkSegAddTimeFromMC or useMCHits or useMCTracks or filterTksUsingMC or tkSegUseMCExpCKTheta: preUnpackSeq.Members += [UnpackMCVertex(), UnpackMCParticle()] if usePixelMCPos: @@ -224,13 +230,13 @@ else: richMCDecode.RejectAllBackgrounds = False richMCDecode.RejectScintillation = False richMCDecode.IncludeTimeInfo = True - richMCDecode.AllowMultipleHits = True + richMCDecode.AllowMultipleHits = False richMCDecode.MaxHitsPerChannel = (99999, 99999) all.Members += [richMCDecode] # Now get the RICH sequence -if filterTksUsingMC or useMCTracks: +if tkSegUseMCExpCKTheta or photUseMC or filterTksUsingMC or useMCTracks: # Build track info from RICH extended MC info fetcher.DataKeys += [ 'pSim/Rich/Segments', 'pSim/Rich/Tracks', 'pSim/MCVertices', @@ -395,9 +401,15 @@ ckResSigma = defaultNSigmaCuts() #ckResSigma[photonSel][tkCKResTreatment][tk][1] = ( 99, 99, 99 ) # Track CK res scale factors -tkCKResScaleF = None # Use nominal values -#scF = 1.2 -#tkCKResScaleF = ( scF, scF, scF ) +#tkCKResScaleF = None # Use nominal values +scF = 0.25 +tkCKResScaleF = ( scF, scF, scF ) + +# CK theta resolution funcs to use for the emulation when enabled +photEmulatedResMC = (os.getenv("R1CKRESFUNC", + str(scF)+"*(0.00078+0.0012*(std::tanh(-x/5000.0)+1.0))"), + os.getenv("R2CKRESFUNC", + str(scF)+"*(0.00065+0.0011*(std::tanh(-x/5000.0)+1.0))")) # Merged output finalPIDLoc = "Rec/Rich/PIDs" @@ -447,6 +459,9 @@ RichRec = RichRecoSequence( enable4DReco=enable4D, tkSegAddTimeFromMC=tkSegAddTimeFromMC, pixUsePosFromMC=usePixelMCPos, + photUseMC=photUseMC, + photEmulatedResMC=photEmulatedResMC, + tkSegUseMCExpCKTheta=tkSegUseMCExpCKTheta, averageHitTime=avSignalTime, pixelTimeWindow=pixTimeWindow, photonTimeWindow=photTimeWindow, -- GitLab From 8bb0ec4f4ccde72fc03bcacca2c8ffa692f347df Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Mon, 11 Dec 2023 11:38:26 +0000 Subject: [PATCH 15/23] Update RICH condor job submission script --- .../RichFutureRecSys/examples/jobs/RunJobs.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/Rich/RichFutureRecSys/examples/jobs/RunJobs.py b/Rich/RichFutureRecSys/examples/jobs/RunJobs.py index 1548abd88c3..893bfbaff8b 100755 --- a/Rich/RichFutureRecSys/examples/jobs/RunJobs.py +++ b/Rich/RichFutureRecSys/examples/jobs/RunJobs.py @@ -46,6 +46,17 @@ parser.add_argument( "--innerPixQ", help="Inner Pixel Quantisation", default="", type=str) parser.add_argument( "--outerPixQ", help="Outer Pixel Quantisation", default="", type=str) +parser.add_argument( + "--ckResFuncR1", + help="RICH1 CK theta resolution function for emulation", + default="", + type=str) +parser.add_argument( + "--ckResFuncR2", + help="RICH1 CK theta resolution function for emulation", + default="", + type=str) + parser.add_argument("-N", "--name", help="Job Name", default="Test", type=str) args = parser.parse_args() config = vars(args) @@ -75,6 +86,10 @@ def unsetenv(var): if var in os.environ: del os.environ[var] +if len(config["ckResFuncR1"]) != 0: + os.environ["R1CKRESFUNC"] = config["ckResFuncR1"] +if len(config["ckResFuncR2"]) != 0: + os.environ["R2CKRESFUNC"] = config["ckResFuncR2"] for Lumi in config["lumis"].split(","): @@ -162,3 +177,7 @@ for Lumi in config["lumis"].split(","): unsetenv("PIXBCKWEIGHTRICH2") unsetenv("LUMI") + + +unsetenv("R1CKRESFUNC") +unsetenv("R2CKRESFUNC") -- GitLab From d48dd59f34a5fade2b1d7f5d9d7f9f9d16803efd Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Mon, 11 Dec 2023 11:38:54 +0000 Subject: [PATCH 16/23] Update GlobalReco/root/RichKaonIDCompareFiles.C --- Rec/GlobalReco/root/RichKaonIDCompareFiles.C | 71 ++++++++++---------- 1 file changed, 36 insertions(+), 35 deletions(-) diff --git a/Rec/GlobalReco/root/RichKaonIDCompareFiles.C b/Rec/GlobalReco/root/RichKaonIDCompareFiles.C index 72adca4e672..512e78ff9be 100755 --- a/Rec/GlobalReco/root/RichKaonIDCompareFiles.C +++ b/Rec/GlobalReco/root/RichKaonIDCompareFiles.C @@ -47,49 +47,50 @@ void RichKaonIDCompareFiles() { } }; - // for ( const std::string bName : {"PixSize-R1-V1", "PixSize-R2-V1", "PixSize-R1R2-V1"} ) { - // for ( const std::string lumi : {"1.2e34"} ) { - // for ( const std::string innerQ : {"1.00", "2.00", "3.00"} ) { - // for ( const std::string outerQ : {"", "2.00", "4.00", "6.00"} ) { - // const auto d1 = nA( "lumi", lumi ) + nA( "InnerPixQ", innerQ ) + nA( "OuterPixQ", outerQ ); - // add( bName + "/3D", bName + "/3D/" + d1 ); - // for ( const std::string pixWin : {"1.000"} ) { - // for ( const std::string photWin : {"0.300", "0.150", "0.075", "0.050"} ) { - // const auto data = bName + "/4D/" + d1 + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin ); - // add( bName + "/4D-PixW_" + pixWin + "-PhoW_" + photWin, data ); - // // add( bName + "/4D-PixW_" + pixWin + "-Lumi_" + lumi, data ); - // } - // } - // } - // } - // } - // } - - const std::string batch = "FromMCHits-V5"; - const std::string innerQ = "2.80"; - const std::string outerQ = "5.60"; - const bool addEmulated = true; - const auto QS = nA( "InnerPixQ", innerQ ) + nA( "OuterPixQ", outerQ ); - for ( const std::string vari : {"Nominal", "MultiHits", "MultiHits-PerfectMCPos"} ) { - const auto bName = batch + "-" + vari; + for ( const std::string bName : {"CKRes/V2/Nominal", "CKRes/V2/Emulated-100pc", "CKRes/V2/Emulated-75pc", + "CKRes/V2/Emulated-50pc", "CKRes/V2/Emulated-25pc"} ) { for ( const std::string lumi : {"1.5e34", "1.2e34", "1.0e34", "3.0e33", "2.0e33"} ) { - add( vari + "/3D", bName + "/3D/" + nA( "lumi", lumi ) ); - if ( addEmulated ) { add( vari + "/3D", "PixSize-R1R2-V1/3D/" + nA( "lumi", lumi ) + QS ); } + // for ( const std::string innerQ : {"1.00", "2.00", "3.00"} ) { + // for ( const std::string outerQ : {"", "2.00", "4.00", "6.00"} ) { + const auto d1 = nA( "lumi", lumi ); // + nA( "InnerPixQ", innerQ ) + nA( "OuterPixQ", outerQ ); + add( bName + "/3D", bName + "/3D/" + d1 ); for ( const std::string pixWin : {"1.000"} ) { for ( const std::string photWin : {"0.300", "0.150", "0.075", "0.050"} ) { - auto data = bName + "/4D/" + nA( "lumi", lumi ) + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin ); - add( vari + "/4D-PixW_" + pixWin + "-PhoW_" + photWin, data ); - add( vari + "/4D-PixW_" + pixWin + "-Lumi_" + lumi, data ); - if ( addEmulated ) { - data = "PixSize-R1R2-V1/4D/" + nA( "lumi", lumi ) + QS + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin ); - add( vari + "/4D-PixW_" + pixWin + "-PhoW_" + photWin, data ); - add( vari + "/4D-PixW_" + pixWin + "-Lumi_" + lumi, data ); - } + const auto data = bName + "/4D/" + d1 + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin ); + add( bName + "/4D-PixW_" + pixWin + "-PhoW_" + photWin, data ); + add( bName + "/4D-PixW_" + pixWin + "-Lumi_" + lumi, data ); } } + // } + // } } } + // const std::string batch = "FromMCHits-V5"; + // const std::string innerQ = "2.80"; + // const std::string outerQ = "5.60"; + // const bool addEmulated = true; + // const auto QS = nA( "InnerPixQ", innerQ ) + nA( "OuterPixQ", outerQ ); + // for ( const std::string vari : {"Nominal", "MultiHits", "MultiHits-PerfectMCPos"} ) { + // const auto bName = batch + "-" + vari; + // for ( const std::string lumi : {"1.5e34", "1.2e34", "1.0e34", "3.0e33", "2.0e33"} ) { + // add( vari + "/3D", bName + "/3D/" + nA( "lumi", lumi ) ); + // if ( addEmulated ) { add( vari + "/3D", "PixSize-R1R2-V1/3D/" + nA( "lumi", lumi ) + QS ); } + // for ( const std::string pixWin : {"1.000"} ) { + // for ( const std::string photWin : {"0.300", "0.150", "0.075", "0.050"} ) { + // auto data = bName + "/4D/" + nA( "lumi", lumi ) + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin ); + // add( vari + "/4D-PixW_" + pixWin + "-PhoW_" + photWin, data ); + // add( vari + "/4D-PixW_" + pixWin + "-Lumi_" + lumi, data ); + // if ( addEmulated ) { + // data = "PixSize-R1R2-V1/4D/" + nA( "lumi", lumi ) + QS + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin + // ); add( vari + "/4D-PixW_" + pixWin + "-PhoW_" + photWin, data ); add( vari + "/4D-PixW_" + pixWin + + // "-Lumi_" + lumi, data ); + // } + // } + // } + // } + // } + // for ( const std::string lumi : {"1.2e34", "2.0e33", "2.0e32"} ) { // const std::string bName = "BkgTuneRich1"; // for ( const std::string w : {"0.600", "0.800", "1.000", "1.200", "1.400"} ) { -- GitLab From d09fdbacea5fb89861e822a53ce5f0648fdba208 Mon Sep 17 00:00:00 2001 From: Gitlab CI Date: Mon, 11 Dec 2023 11:40:01 +0000 Subject: [PATCH 17/23] Fixed formatting patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/34633301 --- Rich/RichFutureRecSys/examples/RichFuture.py | 13 +++++++------ Rich/RichFutureRecSys/examples/jobs/RunJobs.py | 2 +- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py index d4b7b010d72..34ba1de9624 100644 --- a/Rich/RichFutureRecSys/examples/RichFuture.py +++ b/Rich/RichFutureRecSys/examples/RichFuture.py @@ -402,14 +402,15 @@ ckResSigma = defaultNSigmaCuts() # Track CK res scale factors #tkCKResScaleF = None # Use nominal values -scF = 0.25 -tkCKResScaleF = ( scF, scF, scF ) +scF = 0.25 +tkCKResScaleF = (scF, scF, scF) # CK theta resolution funcs to use for the emulation when enabled -photEmulatedResMC = (os.getenv("R1CKRESFUNC", - str(scF)+"*(0.00078+0.0012*(std::tanh(-x/5000.0)+1.0))"), - os.getenv("R2CKRESFUNC", - str(scF)+"*(0.00065+0.0011*(std::tanh(-x/5000.0)+1.0))")) +photEmulatedResMC = ( + os.getenv("R1CKRESFUNC", + str(scF) + "*(0.00078+0.0012*(std::tanh(-x/5000.0)+1.0))"), + os.getenv("R2CKRESFUNC", + str(scF) + "*(0.00065+0.0011*(std::tanh(-x/5000.0)+1.0))")) # Merged output finalPIDLoc = "Rec/Rich/PIDs" diff --git a/Rich/RichFutureRecSys/examples/jobs/RunJobs.py b/Rich/RichFutureRecSys/examples/jobs/RunJobs.py index 893bfbaff8b..a867f5d7c6e 100755 --- a/Rich/RichFutureRecSys/examples/jobs/RunJobs.py +++ b/Rich/RichFutureRecSys/examples/jobs/RunJobs.py @@ -86,6 +86,7 @@ def unsetenv(var): if var in os.environ: del os.environ[var] + if len(config["ckResFuncR1"]) != 0: os.environ["R1CKRESFUNC"] = config["ckResFuncR1"] if len(config["ckResFuncR2"]) != 0: @@ -178,6 +179,5 @@ for Lumi in config["lumis"].split(","): unsetenv("LUMI") - unsetenv("R1CKRESFUNC") unsetenv("R2CKRESFUNC") -- GitLab From 915cc5479d2bd183c21655431c7b107007d2fa29 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Thu, 11 Jan 2024 13:07:29 +0000 Subject: [PATCH 18/23] RichMCOpticalPhotons: Add some new histograms --- .../src/RichMCOpticalPhotons.cpp | 196 ++++++++++++++---- 1 file changed, 161 insertions(+), 35 deletions(-) diff --git a/Rich/RichFutureRecCheckers/src/RichMCOpticalPhotons.cpp b/Rich/RichFutureRecCheckers/src/RichMCOpticalPhotons.cpp index 5fe08f92e6a..2df52cf9c4a 100644 --- a/Rich/RichFutureRecCheckers/src/RichMCOpticalPhotons.cpp +++ b/Rich/RichFutureRecCheckers/src/RichMCOpticalPhotons.cpp @@ -48,6 +48,7 @@ // STD #include #include +#include #include namespace Rich::Future::Rec::MC::Moni { @@ -122,6 +123,8 @@ namespace Rich::Future::Rec::MC::Moni { const DetectorArray panelXsizes{{650, 800}}; const DetectorArray panelYsizes{{680, 750}}; + const DetectorArray richXsizes{{300, 1200}}; + const DetectorArray richYsizes{{300, 1200}}; const DetectorArray shiftXrange{{10, 100}}; const DetectorArray shiftYrange{{10, 25}}; @@ -209,26 +212,80 @@ namespace Rich::Future::Rec::MC::Moni { "Track Momentum (MeV/c)", " / mrad" ) ); // loop over active mass hypos for ( const auto pid : activeParticlesNoBT() ) { - // book yield histos - ok &= saveAndCheck( h_yield[rad][pid], // - richHisto1D( HID( "yield", rad, pid ), // - "Photon Yield (>0)", // - 0, m_maxYield[rad], nBins1D(), // + const auto det = ( Rich::Rich2Gas == rad ? Rich::Rich2 : Rich::Rich1 ); + // Book MC yield histos + ok &= saveAndCheck( h_mcYield[rad][pid], // + richHisto1D( HID( "mcYield", rad, pid ), // + "MC Photon Yield (>0)", // + 0.5, m_maxYield[rad] + 0.5, m_maxYield[rad], // "Photon Yield (>0)" ) ); - ok &= saveAndCheck( h_yieldVp[rad][pid], // - richProfile1D( HID( "yieldVp", rad, pid ), // - "Photon Yield (>0) V P (MeV/c)", // - 1.0 * GeV, 100.0 * GeV, nBins1D(), // + ok &= saveAndCheck( h_mcYieldVp[rad][pid], // + richProfile1D( HID( "mcYieldVp", rad, pid ), // + "MC Photon Yield (>0) V P (MeV/c)", // + m_minP[rad], m_maxP[rad], nBins1D(), // "Track Momentum (MeV/c)", "Photon Yield (>0)" ) ); - ok &= saveAndCheck( h_yieldDiff[rad][pid], // - richHisto1D( HID( "yieldDiff", rad, pid ), // - "Expected-MC Photon Yield (>0)", // - -50.0, 50.0, nBins1D(), // + ok &= saveAndCheck( h_mcYieldVx[rad][pid], // + richProfile1D( HID( "mcYieldVx", rad, pid ), // + "MC Photon Yield (>0) V Panel Local X", // + -richXsizes[det], richXsizes[det], nBins1D(), // + "Panel Local X / mm", "Photon Yield (>0)" ) ); + ok &= saveAndCheck( h_mcYieldVy[rad][pid], // + richProfile1D( HID( "mcYieldVy", rad, pid ), // + "MC Photon Yield (>0) V Panel Local Y", // + -richYsizes[det], richYsizes[det], nBins1D(), // + "Panel Local Y / mm", "Photon Yield (>0)" ) ); + ok &= saveAndCheck( h_mcYieldVxy[rad][pid], // + richProfile2D( HID( "mcYieldVxy", rad, pid ), // + "MC Photon Yield (>0) V Panel Local (X,Y)", // + -richXsizes[det], richXsizes[det], nBins2D(), // + -richYsizes[det], richYsizes[det], nBins2D(), // + "Panel Local X / mm", "Panel Local Y / mm", "Photon Yield (>0)" ) ); + ok &= saveAndCheck( h_mcYieldDiff[rad][pid], // + richHisto1D( HID( "mcYieldDiff", rad, pid ), // + "MC Expected-MC Photon Yield (>0)", // + -50.0, 50.0, nBins1D(), // "Exp-MC Photon Yield" ) ); - ok &= saveAndCheck( h_yieldDiffVp[rad][pid], // - richProfile1D( HID( "yieldDiffVp", rad, pid ), // - "Expected-MC Photon Yield (>0) V P (MeV/c)", // - 1.0 * GeV, 100.0 * GeV, nBins1D(), // + ok &= saveAndCheck( h_mcYieldDiffVp[rad][pid], // + richProfile1D( HID( "mcYieldDiffVp", rad, pid ), // + "MC Expected-MC Photon Yield (>0) V P (MeV/c)", // + 1.0 * GeV, 100.0 * GeV, nBins1D(), // + "Track Momentum (MeV/c)", "Exp-MC Photon Yield" ) ); + // book reco yield histos + ok &= saveAndCheck( h_recoYield[rad][pid], // + richHisto1D( HID( "recoYield", rad, pid ), // + "Reco Photon Yield (>0)", // + 0.5, m_maxYield[rad] + 0.5, m_maxYield[rad], // + "Photon Yield (>0)" ) ); + ok &= saveAndCheck( h_recoYieldVp[rad][pid], // + richProfile1D( HID( "recoYieldVp", rad, pid ), // + "Reco Photon Yield (>0) V P (MeV/c)", // + m_minP[rad], m_maxP[rad], nBins1D(), // + "Track Momentum (MeV/c)", "Photon Yield (>0)" ) ); + ok &= saveAndCheck( h_recoYieldVx[rad][pid], // + richProfile1D( HID( "recoYieldVx", rad, pid ), // + "Reco Photon Yield (>0) V Panel Local X", // + -richXsizes[det], richXsizes[det], nBins1D(), // + "Panel Local X / mm", "Photon Yield (>0)" ) ); + ok &= saveAndCheck( h_recoYieldVy[rad][pid], // + richProfile1D( HID( "recoYieldVy", rad, pid ), // + "Reco Photon Yield (>0) V Panel Local Y", // + -richYsizes[det], richYsizes[det], nBins1D(), // + "Panel Local Y / mm", "Photon Yield (>0)" ) ); + ok &= saveAndCheck( h_recoYieldVxy[rad][pid], // + richProfile2D( HID( "recoYieldVxy", rad, pid ), // + "Reco Photon Yield (>0) V Panel Local (X,Y)", // + -richXsizes[det], richXsizes[det], nBins2D(), // + -richYsizes[det], richYsizes[det], nBins2D(), // + "Panel Local X / mm", "Panel Local Y / mm", "Photon Yield (>0)" ) ); + ok &= saveAndCheck( h_recoYieldDiff[rad][pid], // + richHisto1D( HID( "recoYieldDiff", rad, pid ), // + "Reco Expected-MC Photon Yield (>0)", // + -50.0, 50.0, nBins1D(), // + "Exp-MC Photon Yield" ) ); + ok &= saveAndCheck( h_recoYieldDiffVp[rad][pid], // + richProfile1D( HID( "recoYieldDiffVp", rad, pid ), // + "Reco Expected-MC Photon Yield (>0) V P (MeV/c)", // + 1.0 * GeV, 100.0 * GeV, nBins1D(), // "Track Momentum (MeV/c)", "Exp-MC Photon Yield" ) ); } } @@ -270,6 +327,10 @@ namespace Rich::Future::Rec::MC::Moni { const auto rad = seg.radiator(); if ( !radiatorIsActive( rad ) ) { continue; } + // Segment momentum + const auto pTot = seg.bestMomentumMag(); + if ( pTot < m_minP[rad] || pTot > m_maxP[rad] ) { continue; } + // Track pointer const auto tk = tracks.at( tkIndex ); @@ -281,20 +342,30 @@ namespace Rich::Future::Rec::MC::Moni { // The True MCParticle type const auto pid = mcHelper.mcParticleType( mcP ); if ( mcP && pid != Rich::Unknown ) { + // beta cut for true MC type + const auto mcbeta = richPartProps()->beta( pTot, pid ); + const auto betaC = ( mcbeta >= m_minBeta[rad] && mcbeta <= m_maxBeta[rad] ); + if ( !betaC ) { continue; } if ( yield[pid] > 0 ) { // MC Optical Photons const auto mcPhots = photHelper.mcOpticalPhotons( mcP ); - // Count signal photons in this radiator - const auto nPhots = std::count_if( mcPhots.begin(), mcPhots.end(), // - [&rad]( const auto mcPhot ) { - const auto mcH = ( mcPhot ? mcPhot->mcRichHit() : nullptr ); - return ( mcH && mcH->radiator() == rad && mcH->isSignal() ); - } ); - fillHisto( h_yield[rad][pid], nPhots ); - fillHisto( h_yieldVp[rad][pid], seg.bestMomentumMag(), nPhots ); - const auto yield_diff = yield[pid] - (double)nPhots; - fillHisto( h_yieldDiff[rad][pid], yield_diff ); - fillHisto( h_yieldDiffVp[rad][pid], seg.bestMomentumMag(), yield_diff ); + // Count MC signal photons in this radiator + const auto nPhots = // + std::count_if( mcPhots.begin(), mcPhots.end(), // + [&rad]( const auto mcPhot ) { + const auto mcH = ( mcPhot ? mcPhot->mcRichHit() : nullptr ); + return ( mcH && mcH->radiator() == rad && mcH->isSignal() ); + } ); + if ( nPhots > 0 ) { + fillHisto( h_mcYield[rad][pid], nPhots ); + fillHisto( h_mcYieldVp[rad][pid], pTot, nPhots ); + fillHisto( h_mcYieldVx[rad][pid], seg.entryPoint().X(), nPhots ); + fillHisto( h_mcYieldVy[rad][pid], seg.entryPoint().Y(), nPhots ); + fillHisto( h_mcYieldVxy[rad][pid], seg.entryPoint().X(), seg.entryPoint().Y(), nPhots ); + const auto yield_diff = yield[pid] - (double)nPhots; + fillHisto( h_mcYieldDiff[rad][pid], yield_diff ); + fillHisto( h_mcYieldDiffVp[rad][pid], pTot, yield_diff ); + } } } } // MCPs @@ -310,6 +381,10 @@ namespace Rich::Future::Rec::MC::Moni { // Need MC info so skip tracks without any if ( mcPs.empty() ) { continue; } + // True MC Signal counts for this track + RadiatorArray nSigPhots{{}}; + RadiatorArray segIndex{-1, -1, -1}; + // loop over photons for this track for ( const auto photIn : sumTk.photonIndices() ) { @@ -326,6 +401,7 @@ namespace Rich::Future::Rec::MC::Moni { // Rich / Radiator info const auto rich = seg.rich(); const auto rad = seg.radiator(); + segIndex[rad] = rels.segmentIndex(); if ( !radiatorIsActive( rad ) ) { continue; } // Segment momentum @@ -352,6 +428,9 @@ namespace Rich::Future::Rec::MC::Moni { // Weight per MCPhoton const auto mcPW = ( !mcPhots.empty() ? 1.0 / (double)mcPhots.size() : 1.0 ); + // Count true signal reco photons for each segment + if ( !mcPhots.empty() ) { ++nSigPhots[rad]; } + // Loop over MC photons for ( const auto mcPhot : mcPhots ) { @@ -395,7 +474,37 @@ namespace Rich::Future::Rec::MC::Moni { } // SIMD scalar loop } // photons - } // tracks + + // Reco yield plots + for ( const auto rad : activeRadiators() ) { + if ( segIndex[rad] >= 0 ) { + const auto pid = mcHelper.mcParticleType( mcPs.front() ); + if ( pid != Rich::Unknown ) { + const auto& yield = yields[segIndex[rad]]; + if ( yield[pid] > 0 ) { + const auto& seg = segments[segIndex[rad]]; + const auto pTot = seg.bestMomentumMag(); + // beta cut for true MC type + const auto mcbeta = richPartProps()->beta( pTot, pid ); + const auto betaC = ( mcbeta >= m_minBeta[rad] && mcbeta <= m_maxBeta[rad] ); + if ( !betaC ) { continue; } + const auto nPhots = nSigPhots[rad]; + if ( nPhots > 0 ) { + fillHisto( h_recoYield[rad][pid], nPhots ); + fillHisto( h_recoYieldVp[rad][pid], pTot, nPhots ); + fillHisto( h_recoYieldVx[rad][pid], seg.entryPoint().X(), nPhots ); + fillHisto( h_recoYieldVy[rad][pid], seg.entryPoint().Y(), nPhots ); + fillHisto( h_recoYieldVxy[rad][pid], seg.entryPoint().X(), seg.entryPoint().Y(), nPhots ); + const auto yield_diff = yield[pid] - (double)nPhots; + fillHisto( h_recoYieldDiff[rad][pid], yield_diff ); + fillHisto( h_recoYieldDiffVp[rad][pid], pTot, yield_diff ); + } + } + } + } + } + + } // tracks } private: @@ -421,21 +530,38 @@ namespace Rich::Future::Rec::MC::Moni { DetectorArray h_pixXLocShift2D = {{}}; DetectorArray h_pixYLocShift2D = {{}}; - RadiatorArray> h_yield = {{}}; - RadiatorArray> h_yieldVp = {{}}; - RadiatorArray> h_yieldDiff = {{}}; - RadiatorArray> h_yieldDiffVp = {{}}; + RadiatorArray> h_mcYield = {{}}; + RadiatorArray> h_mcYieldVp = {{}}; + RadiatorArray> h_mcYieldVx = {{}}; + RadiatorArray> h_mcYieldVy = {{}}; + RadiatorArray> h_mcYieldVxy = {{}}; + RadiatorArray> h_mcYieldDiff = {{}}; + RadiatorArray> h_mcYieldDiffVp = {{}}; + + RadiatorArray> h_recoYield = {{}}; + RadiatorArray> h_recoYieldVp = {{}}; + RadiatorArray> h_recoYieldVx = {{}}; + RadiatorArray> h_recoYieldVy = {{}}; + RadiatorArray> h_recoYieldVxy = {{}}; + RadiatorArray> h_recoYieldDiff = {{}}; + RadiatorArray> h_recoYieldDiffVp = {{}}; private: // JOs /// Min momentum by radiator (MeV/c) Gaudi::Property> m_minP{ - this, "MinP", {2.0 * Gaudi::Units::GeV, 3.0 * Gaudi::Units::GeV, 30.0 * Gaudi::Units::GeV}}; + this, "MinP", {2.0 * Gaudi::Units::GeV, 3.0 * Gaudi::Units::GeV, 16.0 * Gaudi::Units::GeV}}; /// Max momentum by radiator (MeV/c) Gaudi::Property> m_maxP{ - this, "MaxP", {20.0 * Gaudi::Units::GeV, 70.0 * Gaudi::Units::GeV, 120.0 * Gaudi::Units::GeV}}; + this, "MaxP", {20.0 * Gaudi::Units::GeV, 80.0 * Gaudi::Units::GeV, 120.0 * Gaudi::Units::GeV}}; + + /// minimum beta value for tracks + Gaudi::Property> m_minBeta{this, "MinBeta", {0.9999f, 0.9999f, 0.9999f}}; + + /// maximum beta value for tracks + Gaudi::Property> m_maxBeta{this, "MaxBeta", {999.99f, 999.99f, 999.99f}}; /// Maximum photon yield Gaudi::Property> m_maxYield{this, "MaximumYields", {80, 80, 80}}; -- GitLab From 01cff13727871e2dbdb14ac3942d8958c7aedf3d Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Thu, 11 Jan 2024 13:08:03 +0000 Subject: [PATCH 19/23] RichTrSegMakerFromMCRichTracks: Add option to enable some rate of ghost track creation --- .../src/RichTrSegMakerFromMCRichTracks.cpp | 37 ++++++++++++++++++- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/Rich/RichFutureRecMCAlgorithms/src/RichTrSegMakerFromMCRichTracks.cpp b/Rich/RichFutureRecMCAlgorithms/src/RichTrSegMakerFromMCRichTracks.cpp index 7fc85fd97bb..ac8bb6dfb49 100644 --- a/Rich/RichFutureRecMCAlgorithms/src/RichTrSegMakerFromMCRichTracks.cpp +++ b/Rich/RichFutureRecMCAlgorithms/src/RichTrSegMakerFromMCRichTracks.cpp @@ -157,6 +157,7 @@ namespace Rich::Future::Rec::MC { } if ( !radiatorIsActive( Rich::Rich1Gas ) ) { _ri_debug << "Track segments for Rich1Gas are disabled" << endmsg; } if ( !radiatorIsActive( Rich::Rich2Gas ) ) { _ri_debug << "Track segments for Rich2Gas are disabled" << endmsg; } + // derived condition data Detector::Rich1::addConditionDerivation( this ); Detector::Rich2::addConditionDerivation( this ); @@ -172,6 +173,8 @@ namespace Rich::Future::Rec::MC { m_randZeroOne = Rndm::Numbers( randSvc(), Rndm::Flat( 0.0, 1.0 ) ); + if ( m_ghostFrac > 0.0 ) { info() << "Will create " << 100.0 * m_ghostFrac << "% ghost tracks" << endmsg; } + // return return sc; } @@ -299,6 +302,12 @@ namespace Rich::Future::Rec::MC { } } + // Should we use this track to create a fake ghost track ? + const bool asGhost = ( m_randZeroOne.shoot() <= m_ghostFrac ); + // Flip track parameters in X and/or Y (must be at least one) + const bool ghostXflip = ( asGhost && m_randZeroOne.shoot() < 0.5 ); + const bool ghostYflip = ( asGhost && ( !ghostXflip || m_randZeroOne.shoot() < 0.5 ) ); + // Loop over all radiators for ( const auto radiator : detInfo.radiators ) { @@ -392,6 +401,27 @@ namespace Rich::Future::Rec::MC { continue; } + // If we are using this track as a ghost, flip the data in X and/or Y + // to break any correlations with RICH hits. + if ( asGhost ) { + if ( ghostXflip ) { + entryStateMomentum.SetX( -entryStateMomentum.X() ); + exitStateMomentum.SetX( -exitStateMomentum.X() ); + midStateMomentum.SetX( -midStateMomentum.X() ); + entryPoint.SetX( -entryPoint.X() ); + exitPoint.SetX( -exitPoint.X() ); + midPoint.SetX( -midPoint.X() ); + } + if ( ghostYflip ) { + entryStateMomentum.SetY( -entryStateMomentum.Y() ); + exitStateMomentum.SetY( -exitStateMomentum.Y() ); + midStateMomentum.SetY( -midStateMomentum.Y() ); + entryPoint.SetY( -entryPoint.Y() ); + exitPoint.SetY( -exitPoint.Y() ); + midPoint.SetY( -midPoint.Y() ); + } + } + _ri_verbo << " -> Passed selection cuts " << endmsg; _ri_verbo << " -> Rad Points | entry=" << entryPoint << " mid=" << midPoint << " exit=" << exitPoint << endmsg; @@ -432,8 +462,8 @@ namespace Rich::Future::Rec::MC { tkToSegsRel.emplace_back( tkKey, tkIndex ); auto& tkRels = tkToSegsRel.back(); tkRels.segmentIndices = std::move( segList ); - // fill Track <-> MCParticle link table - tkToMCPs[mcp].push_back( tk.get() ); + // fill Track <-> MCParticle link table (if not being made as a ghost) + if ( !asGhost ) { tkToMCPs[mcp].push_back( tk.get() ); } // Finally pass ownership to container _ri_verbo << "Created Track : " << *tk << endmsg; tks.insert( tk.release(), tkKey ); @@ -522,6 +552,9 @@ namespace Rich::Future::Rec::MC { /// Allow use of ideal state creator when RICH data is missing Gaudi::Property> m_useIdealStates{this, "UseIdealStates", {true, true, true}}; + /// Fraction of tracks to use to make 'ghosts' + Gaudi::Property m_ghostFrac{this, "GhostFraction", 0.0}; + /// Nominal z positions of states at RICHes Gaudi::Property> m_nomZstates{ this, -- GitLab From 102be43ce24b18ae4019537cf98285f2c48e8b16 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Thu, 11 Jan 2024 13:08:39 +0000 Subject: [PATCH 20/23] RichTrackParameterisedCherenkovResolutions: clean up --- ...TrackParameterisedCherenkovResolutions.cpp | 217 +++++++++--------- 1 file changed, 103 insertions(+), 114 deletions(-) diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackParameterisedCherenkovResolutions.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackParameterisedCherenkovResolutions.cpp index bd7d2542d9a..32090051e11 100644 --- a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackParameterisedCherenkovResolutions.cpp +++ b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackParameterisedCherenkovResolutions.cpp @@ -9,11 +9,6 @@ * or submit itself to any jurisdiction. * \*****************************************************************************/ -// STL -#include -#include -#include - // Gaudi Kernel #include "GaudiKernel/ParsersFactory.h" #include "GaudiKernel/PhysicalConstants.h" @@ -41,8 +36,17 @@ // Resolution parameters #include "ResolutionParameters.h" +// STL +#include +#include +#include + namespace Rich::Future::Rec { + namespace { + using OutD = CherenkovResolutions::Vector; + } + /** @class TrackParameterisedCherenkovResolutions * * Computes the expected Cherenkov resolutions for the given track segments @@ -51,8 +55,8 @@ namespace Rich::Future::Rec { * @date 2016-09-30 */ class TrackParameterisedCherenkovResolutions final - : public LHCb::Algorithm::Transformer>> { public: @@ -69,13 +73,95 @@ namespace Rich::Future::Rec { } /// Initialization after creation - StatusCode initialize() override; + StatusCode initialize() override { + + // Sets up various tools and services + auto sc = Transformer::initialize(); + if ( !sc ) return sc; + + // Track type + auto tkType = m_tkType.value(); + // handle minor name variations... + if ( "Upstream" == tkType ) { tkType = "Up"; } + if ( "Downstream" == tkType ) { tkType = "Down"; } + // If 'MC', i.e. we are faking tracks using the RICH extended MC info, + // then just use parameters for Long tracks + if ( "MC" == tkType ) { tkType = "Long"; } + + // PID type loop + for ( const auto pid : activeParticles() ) { + + // Load the resolution parameters for the configured track type + const std::string mode = ( m_modes[pid] == Momentum ? "Momentum" : "CKTheta" ); + if ( tkType.empty() ) { + error() << "Track type must be configured" << endmsg; + return StatusCode::FAILURE; + } + const auto params = resolutionParameters( tkType, mode ); + + // Initialise the interpolators for each radiator + for ( const auto rad : activeRadiators() ) { + if ( params[rad][pid].empty() || !m_resInterp[rad][pid].initInterpolator( params[rad][pid] ) ) { + error() << "Failed to initialise interpolator for " << rad << " " << pid << " " << m_tkType << endmsg; + return StatusCode::FAILURE; + } + } + + } // pids + + // return + return sc; + } public: /// Algorithm execution via transform - CherenkovResolutions::Vector operator()( const LHCb::RichTrackSegment::Vector& segments, // - const CherenkovAngles::Vector& ckAngles // - ) const override; + OutD operator()( const LHCb::RichTrackSegment::Vector& segments, // + const CherenkovAngles::Vector& ckAngles ) const override { + + using namespace Gaudi::Units; + + // data to return + OutD resV; + resV.reserve( segments.size() ); + + // Loop over input segment data + for ( auto&& [segment, angles] : Ranges::ConstZip( segments, ckAngles ) ) { + + // Add a resolution entry for this segment + auto& res = resV.emplace_back(); + + // Is the lowest mass hypo for this track above threshold + const bool aboveThreshold = ( angles[lightestActiveHypo()] > 0 ); + + // check if segment is valid + if ( aboveThreshold ) { + + // Which radiator + const auto rad = segment.radiator(); + + // Loop over active hypos ( note including BT here ) + for ( const auto pid : activeParticles() ) { + + // get the variable for the interpolation mode, clamped to valid range + const double var = std::clamp( ( m_modes[pid] == Momentum ? segment.bestMomentumMag() : angles[pid] ), + m_resInterp[rad][pid].minX(), m_resInterp[rad][pid].maxX() ); + + // Compute the final resolution, including global scale factor + const auto ckRes = + std::clamp( m_scale[rad] * m_resInterp[rad][pid].value( var ), m_minRes[rad], m_maxRes[rad] ); + //_ri_verbo << rad << " " << pid << " " << segment.bestMomentumMag() << " " << ckRes << endmsg; + + // save value + res.setData( pid, ckRes ); + + } // hypo loop + + } // above threshold + + } // track loop + + return resV; + } private: // types @@ -108,108 +194,11 @@ namespace Rich::Future::Rec { /// Absolute max CK theta resolution per radiator Gaudi::Property> m_maxRes{this, "MaxCKThetaRes", {0.003f, 0.0025f, 0.001f}}; - }; // namespace Rich::Future::Rec - -} // namespace Rich::Future::Rec + /// Absolute min CK theta resolution per radiator + Gaudi::Property> m_minRes{this, "MinCKThetaRes", {0.0f, 0.0f, 0.0f}}; + }; -// All code is in general Rich reconstruction namespace -using namespace Rich::Future::Rec; + // Declaration of the Algorithm Factory + DECLARE_COMPONENT( TrackParameterisedCherenkovResolutions ) -//============================================================================= - -StatusCode TrackParameterisedCherenkovResolutions::initialize() { - - // Sets up various tools and services - auto sc = Transformer::initialize(); - if ( !sc ) return sc; - - // Track type - auto tkType = m_tkType.value(); - // handle minor name variations... - if ( "Upstream" == tkType ) { tkType = "Up"; } - if ( "Downstream" == tkType ) { tkType = "Down"; } - // If 'MC', i.e. we are faking tracks using the RICH extended MC info, - // then just use parameters for Long tracks - if ( "MC" == tkType ) { tkType = "Long"; } - - // PID type loop - for ( const auto pid : activeParticles() ) { - - // Load the resolution parameters for the configured track type - const std::string mode = ( m_modes[pid] == Momentum ? "Momentum" : "CKTheta" ); - if ( tkType.empty() ) { - error() << "Track type must be configured" << endmsg; - return StatusCode::FAILURE; - } - const auto params = resolutionParameters( tkType, mode ); - - // Initialise the interpolators for each radiator - for ( const auto rad : activeRadiators() ) { - if ( params[rad][pid].empty() || !m_resInterp[rad][pid].initInterpolator( params[rad][pid] ) ) { - error() << "Failed to initialise interpolator for " << rad << " " << pid << " " << m_tkType << endmsg; - return StatusCode::FAILURE; - } - } - - } // pids - - // return - return sc; -} - -//============================================================================= - -CherenkovResolutions::Vector // -TrackParameterisedCherenkovResolutions::operator()( const LHCb::RichTrackSegment::Vector& segments, // - const CherenkovAngles::Vector& ckAngles ) const { - - using namespace Gaudi::Units; - - // data to return - CherenkovResolutions::Vector resV; - resV.reserve( segments.size() ); - - // Loop over input segment data - for ( auto&& [segment, angles] : Ranges::ConstZip( segments, ckAngles ) ) { - - // Add a resolution entry for this segment - auto& res = resV.emplace_back(); - - // Is the lowest mass hypo for this track above threshold - const bool aboveThreshold = ( angles[lightestActiveHypo()] > 0 ); - - // check if segment is valid - if ( aboveThreshold ) { - - // Which radiator - const auto rad = segment.radiator(); - - // Loop over active hypos ( note including BT here ) - for ( const auto pid : activeParticles() ) { - - // get the variable for the interpolation mode, clamped to valid range - const double var = std::clamp( ( m_modes[pid] == Momentum ? segment.bestMomentumMag() : angles[pid] ), - m_resInterp[rad][pid].minX(), m_resInterp[rad][pid].maxX() ); - - // Compute the final resolution, including global scale factor - const auto ckRes = std::min( m_scale[rad] * m_resInterp[rad][pid].value( var ), m_maxRes[rad] ); - //_ri_verbo << rad << " " << pid << " " << segment.bestMomentumMag() << " " << ckRes << endmsg; - - // save value - res.setData( pid, ckRes ); - - } // hypo loop - - } // above threshold - - } // track loop - - return resV; -} - -//============================================================================= - -// Declaration of the Algorithm Factory -DECLARE_COMPONENT( TrackParameterisedCherenkovResolutions ) - -//============================================================================= +} // namespace Rich::Future::Rec -- GitLab From 36df59edc68c3f3c09779cc32228af2571ee4585 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Thu, 11 Jan 2024 13:09:17 +0000 Subject: [PATCH 21/23] Fix typo in submit script option comment --- Rich/RichFutureRecSys/examples/jobs/RunJobs.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Rich/RichFutureRecSys/examples/jobs/RunJobs.py b/Rich/RichFutureRecSys/examples/jobs/RunJobs.py index a867f5d7c6e..884b9286f7b 100755 --- a/Rich/RichFutureRecSys/examples/jobs/RunJobs.py +++ b/Rich/RichFutureRecSys/examples/jobs/RunJobs.py @@ -53,10 +53,9 @@ parser.add_argument( type=str) parser.add_argument( "--ckResFuncR2", - help="RICH1 CK theta resolution function for emulation", + help="RICH2 CK theta resolution function for emulation", default="", type=str) - parser.add_argument("-N", "--name", help="Job Name", default="Test", type=str) args = parser.parse_args() config = vars(args) -- GitLab From a00dce47de9d654cd51707bf957280857f725388 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Thu, 11 Jan 2024 13:09:52 +0000 Subject: [PATCH 22/23] RichFutureRecSys: Update example options --- Rich/RichFutureRecSys/examples/RichFuture.py | 23 +++++++++++++------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py index 34ba1de9624..6db0452864a 100644 --- a/Rich/RichFutureRecSys/examples/RichFuture.py +++ b/Rich/RichFutureRecSys/examples/RichFuture.py @@ -58,7 +58,7 @@ rootFileBaseName = "RichFuture-" + myBuild + "-" + myConf + "-" + histos HistogramPersistencySvc().OutputFile = rootFileBaseName + "-Histos.root" # Event numbers -nEvents = (1000 if not batchMode else 99999999) +nEvents = (10000 if not batchMode else 99999999) EventSelector().PrintFreq = (100 if not batchMode else 250) #LHCbApp().SkipEvents = 2 @@ -105,8 +105,8 @@ if UseDD4Hep: # -------------------------------------------------------------------------------------- # Reco options -useMCTracks = True -useMCHits = True +useMCTracks = False +useMCHits = False # Get time windows from environment vars if set # otherwise use default values as passed. @@ -139,7 +139,7 @@ usePixelMCTime = ((enable4D[0] or enable4D[1]) and not useMCHits) usePixelMCPos = False # Cheat photon parameters from MC -photUseMC = True +photUseMC = False # Cheat expected track CK theta valus from MC tkSegUseMCExpCKTheta = False @@ -230,8 +230,12 @@ else: richMCDecode.RejectAllBackgrounds = False richMCDecode.RejectScintillation = False richMCDecode.IncludeTimeInfo = True - richMCDecode.AllowMultipleHits = False + richMCDecode.AllowMultipleHits = True richMCDecode.MaxHitsPerChannel = (99999, 99999) + # Nominal (0.95, 0.95) + # 50% resolution (0.75, 0.85) + # 25% resolution (0.55, 0.75) + richMCDecode.ReadoutEfficiency = (0.95, 0.95) all.Members += [richMCDecode] # Now get the RICH sequence @@ -401,9 +405,7 @@ ckResSigma = defaultNSigmaCuts() #ckResSigma[photonSel][tkCKResTreatment][tk][1] = ( 99, 99, 99 ) # Track CK res scale factors -#tkCKResScaleF = None # Use nominal values -scF = 0.25 -tkCKResScaleF = (scF, scF, scF) +scF = 1.0 # CK theta resolution funcs to use for the emulation when enabled photEmulatedResMC = ( @@ -412,6 +414,11 @@ photEmulatedResMC = ( os.getenv("R2CKRESFUNC", str(scF) + "*(0.00065+0.0011*(std::tanh(-x/5000.0)+1.0))")) +# Track CK res scale factors +#tkCKResScaleF = None # Use nominal values (1.0) +if scF < 0.5: scF = 0.5 +tkCKResScaleF = (scF, scF, scF) + # Merged output finalPIDLoc = "Rec/Rich/PIDs" -- GitLab From 555fce7ae07571cbf3a764b03a4bf5d530a9f8dd Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Thu, 11 Jan 2024 13:10:17 +0000 Subject: [PATCH 23/23] GlobalReco: Update ROOT script for RICH plots --- Rec/GlobalReco/root/RichKaonIDCompareFiles.C | 59 +++++++++++--------- 1 file changed, 32 insertions(+), 27 deletions(-) diff --git a/Rec/GlobalReco/root/RichKaonIDCompareFiles.C b/Rec/GlobalReco/root/RichKaonIDCompareFiles.C index 512e78ff9be..786b3088d34 100755 --- a/Rec/GlobalReco/root/RichKaonIDCompareFiles.C +++ b/Rec/GlobalReco/root/RichKaonIDCompareFiles.C @@ -23,23 +23,21 @@ void RichKaonIDCompareFiles() { - const std::string dir = "/usera/jonesc/LHCb/output/U2/InclB/WithSpill"; - // const std::string dir = "/usera/jonesc/LHCbCMake/Feature/Rec/output/test"; + // const std::string dir = "/usera/jonesc/LHCb/output/U2/InclB/WithSpill"; + const std::string dir = "/usera/jonesc/LHCbCMake/Feature/Rec/output/test/RecoHits"; - const std::string fName = "RichFuture-Feature-x86_64_v2-centos7-gcc12+detdesc-opt-Expert-ProtoTuple.root"; - // const std::string fName = "RichFuture-Feature-x86_64_v3-centos7-gcc12+detdesc-opt-Expert-ProtoTuple.root"; + // const std::string fName = "RichFuture-Feature-x86_64_v2-centos7-gcc12+detdesc-opt-Expert-ProtoTuple.root"; + const std::string fName = "RichFuture-Feature-x86_64_v3-centos7-gcc12+detdesc-opt-Expert-ProtoTuple.root"; std::map> dataSets; - // dataSets["Run3Tracking-3D"].push_back( "DSTs-MCTracks" ); - // dataSets["Run3Tracking-3D"].push_back( "DSTs-RecoTracks" ); - auto nA = []( const std::string& label, const std::string& value ) { return ( value.empty() ? "" : label + "-" + value + "/" ); }; auto add = [&]( const auto& tag, const auto& name ) { - const auto full_name = dir + "/" + name + "/" + fName; + auto full_name = dir + "/" + name + "/" + fName; + boost::replace_all( full_name, "//", "/" ); if ( boost::filesystem::exists( full_name ) ) { dataSets[tag].push_back( name ); } else { @@ -47,24 +45,30 @@ void RichKaonIDCompareFiles() { } }; - for ( const std::string bName : {"CKRes/V2/Nominal", "CKRes/V2/Emulated-100pc", "CKRes/V2/Emulated-75pc", - "CKRes/V2/Emulated-50pc", "CKRes/V2/Emulated-25pc"} ) { - for ( const std::string lumi : {"1.5e34", "1.2e34", "1.0e34", "3.0e33", "2.0e33"} ) { - // for ( const std::string innerQ : {"1.00", "2.00", "3.00"} ) { - // for ( const std::string outerQ : {"", "2.00", "4.00", "6.00"} ) { - const auto d1 = nA( "lumi", lumi ); // + nA( "InnerPixQ", innerQ ) + nA( "OuterPixQ", outerQ ); - add( bName + "/3D", bName + "/3D/" + d1 ); - for ( const std::string pixWin : {"1.000"} ) { - for ( const std::string photWin : {"0.300", "0.150", "0.075", "0.050"} ) { - const auto data = bName + "/4D/" + d1 + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin ); - add( bName + "/4D-PixW_" + pixWin + "-PhoW_" + photWin, data ); - add( bName + "/4D-PixW_" + pixWin + "-Lumi_" + lumi, data ); - } - } - // } - // } - } - } + add( "Run3Tracking-3D", "RecoTracks" ); + add( "Run3Tracking-3D", "MCTracks-NoGhosts" ); + add( "Run3Tracking-3D", "MCTracks-5pcGhosts" ); + + // for ( const std::string batch : {"CKRes/V2/"} ) { + // for ( const std::string N : {"Nominal", "Emulated-100pc", "Emulated-75pc", "Emulated-50pc", "Emulated-25pc"} ) { + // const auto bName = batch + N; + // for ( const std::string lumi : {"1.5e34", "1.2e34", "1.0e34", "3.0e33", "2.0e33"} ) { + // for ( const std::string innerQ : {"", "1.00"} ) { + // for ( const std::string outerQ : {"", "1.00"} ) { + // const auto d1 = nA( "lumi", lumi ) + nA( "InnerPixQ", innerQ ) + nA( "OuterPixQ", outerQ ); + // add( bName + "/3D", bName + "/3D/" + d1 ); + // const std::string pixWin = "1.000"; + // for ( const std::string photWin : {"0.300", "0.150", "0.075", "0.050"} ) { + // const auto data = bName + "/4D/" + d1 + nA( "PixWin", pixWin ) + nA( "PhotWin", photWin ); + // add( bName + "/4D-PhoW_" + photWin, data ); + // add( bName + "/4D-Lumi_" + lumi, data ); + // add( batch + "4D-Lumi_" + lumi + "-PhoW_" + photWin, data ); + // } + // } + // } + // } + // } + // } // const std::string batch = "FromMCHits-V5"; // const std::string innerQ = "2.80"; @@ -123,7 +127,8 @@ void RichKaonIDCompareFiles() { const Long64_t nTracks = 1e6; - const double GeV = 1000; + const double MeV = 1.0; + const double GeV = 1000 * MeV; for ( const auto& [dname, tags] : dataSets ) { -- GitLab