From 33a8cf16e2eef02bd4dd9139a6c874cbdd76a7a0 Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Mon, 17 Oct 2022 11:09:12 +0200 Subject: [PATCH 01/15] Dropped unused RichDigiAlgMoni class --- Rich/RichReadout/CMakeLists.txt | 1 - .../src/component/RichDigiAlgMoni.cpp | 428 ------------------ .../src/component/RichDigiAlgMoni.h | 155 ------- 3 files changed, 584 deletions(-) delete mode 100644 Rich/RichReadout/src/component/RichDigiAlgMoni.cpp delete mode 100755 Rich/RichReadout/src/component/RichDigiAlgMoni.h diff --git a/Rich/RichReadout/CMakeLists.txt b/Rich/RichReadout/CMakeLists.txt index 260efc860..2288b3fac 100644 --- a/Rich/RichReadout/CMakeLists.txt +++ b/Rich/RichReadout/CMakeLists.txt @@ -28,7 +28,6 @@ gaudi_add_module(RichReadout src/component/RichCopySummedDepositsToDigits.cpp src/component/RichDetailedFrontEndResponse.cpp src/component/RichDetailedFrontEndResponsePMT.cpp - src/component/RichDigiAlgMoni.cpp src/component/RichPixelReadout.cpp src/component/RichRegistry.cpp src/component/RichSignal.cpp diff --git a/Rich/RichReadout/src/component/RichDigiAlgMoni.cpp b/Rich/RichReadout/src/component/RichDigiAlgMoni.cpp deleted file mode 100644 index bba538a9a..000000000 --- a/Rich/RichReadout/src/component/RichDigiAlgMoni.cpp +++ /dev/null @@ -1,428 +0,0 @@ -/*****************************************************************************\ -* (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. * -\*****************************************************************************/ -// Units -#include "GaudiKernel/SystemOfUnits.h" - -// From GaudiKernel -#include "GaudiKernel/SmartDataPtr.h" - -// From PartProp -#include "Kernel/IParticlePropertySvc.h" -#include "Kernel/ParticleProperty.h" - -// local -#include "RichDigiAlgMoni.h" - -//----------------------------------------------------------------------------- -// Implementation file for class : RichDigiAlgMoni -// -// 2003-09-08 : Chris Jones (Christopher.Rob.Jones@cern.ch) -//----------------------------------------------------------------------------- - -using namespace Rich::MC::Digi; - -DECLARE_COMPONENT( AlgMoni ) - -// Standard constructor, initializes variables -AlgMoni::AlgMoni( const std::string& name, ISvcLocator* pSvcLocator ) - : HistoAlgBase( name, pSvcLocator ), m_ID( 0 ), m_maxID( 500 ) { - // Declare job options - declareProperty( "InputMCRichDigits", m_mcdigitTES = LHCb::MCRichDigitLocation::Default ); - declareProperty( "InputMCRichDeposits", m_mcdepTES = LHCb::MCRichDepositLocation::Default ); - declareProperty( "InputMCRichHits", m_mchitTES = LHCb::MCRichHitLocation::Default ); - declareProperty( "MaxHPDHistos", m_maxID ); -} - -// Initialisation -StatusCode AlgMoni::initialize() { - // Initialize base class - const StatusCode sc = HistoAlgBase::initialize(); - if ( sc.isFailure() ) { return sc; } - - // get tools - acquireTool( "RichSmartIDTool", m_smartIDTool, nullptr, true ); - acquireTool( "RichMCTruthTool", m_mcTool, nullptr, true ); - - // RichDet - m_richSys = getDet( DeRichLocations::RichSystem ); - - // Retrieve particle property service - LHCb::IParticlePropertySvc* ppSvc = svc( "LHCb::ParticlePropertySvc", true ); - - // Retrieve particle masses - m_particleMass[Rich::Unknown] = 0; - m_particleMass[Rich::Electron] = ppSvc->find( "e+" )->mass() / Gaudi::Units::MeV; - m_particleMass[Rich::Muon] = ppSvc->find( "mu+" )->mass() / Gaudi::Units::MeV; - m_particleMass[Rich::Pion] = ppSvc->find( "pi+" )->mass() / Gaudi::Units::MeV; - m_particleMass[Rich::Kaon] = ppSvc->find( "K+" )->mass() / Gaudi::Units::MeV; - m_particleMass[Rich::Proton] = ppSvc->find( "p+" )->mass() / Gaudi::Units::MeV; - - // release particle property service - release( ppSvc ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - - return sc; -} - -// Main execution -StatusCode AlgMoni::execute() { - _ri_debug << "Execute" << endmsg; - - // Maps for number of pe's - PhotMap ckPhotMapHit; - PhotMap ckPhotMapDig; - PhotMap ckPhotMapDep; - - // Locate MCRichDigits - // ================================================================================ - auto* mcRichDigits = get( m_mcdigitTES ); - _ri_debug << "Successfully located " << mcRichDigits->size() << " MCRichDigits at " << m_mcdigitTES << endmsg; - - // PD occupancy - PDMulti pdMult; - - // Initialise mult counts - std::vector digMult( Rich::NRiches, 0 ); - std::vector hpdCKMult( Rich::NRiches, 0 ); - - Rich::Map hasHPDQuartzCKBkg; - - // Histogramming - PD_GLOBAL_POSITIONS; - const double maxMult[] = {5000, 2000}; - const Rich::HistoID hid; - - // Loop over all MCRichDigits - for ( const auto* McDig : *mcRichDigits ) { - - // Smart ID - const auto id = McDig->key(); - - // position for this ID - Gaudi::XYZPoint point; - const auto sc = m_smartIDTool->globalPosition( id, point ); - if ( !sc ) { continue; } - - // RICH - const auto rich = id.rich(); - - // increment digit count - ++digMult[rich]; - - // increment PD multiplicity count - ++pdMult[id.pdID()]; - if ( McDig->history().hpdQuartzCK() ) { - ++hpdCKMult[rich]; - hasHPDQuartzCKBkg[id.pdID()] = true; - } - - // Position plots - plot1D( point.x(), hid( rich, "digitXglo" ), "digits x global", xMinPDGlo[rich], xMaxPDGlo[rich] ); - plot1D( point.x(), hid( rich, "digitYglo" ), "digits y global", yMinPDGlo[rich], yMaxPDGlo[rich] ); - plot1D( point.x(), hid( rich, "digitZglo" ), "digits z global", zMinPDGlo[rich], zMaxPDGlo[rich] ); - plot2D( point.x(), point.y(), hid( rich, "digitXYglo" ), "digits yVx global", xMinPDGlo[rich], xMaxPDGlo[rich], - yMinPDGlo[rich], yMaxPDGlo[rich] ); - plot2D( point.z(), point.x(), hid( rich, "digitZXglo" ), "digits xVz global", zMinPDGlo[rich], zMaxPDGlo[rich], - xMinPDGlo[rich], xMaxPDGlo[rich] ); - plot2D( point.z(), point.y(), hid( rich, "digitZYglo" ), "digits yVz global", zMinPDGlo[rich], zMaxPDGlo[rich], - yMinPDGlo[rich], yMaxPDGlo[rich] ); - - // loop over all hits associated to the digit - const auto& mcHits = McDig->hits(); - if ( mcHits.empty() ) { warning() << "MCRichDigit " << id << " has no MCRichHits..." << endmsg; } - - // hit mult per digit - plot1D( mcHits.size(), hid( rich, "hitsPerDig" ), "Hit mult. per Digit", -0.5, 10.5, 11 ); - - bool thisDigCounted = false; - for ( const auto& mchit : mcHits ) { - const auto* hit = mchit.mcRichHit(); - - // Compare digit/hit err for signal hits only - if ( !hit->isBackground() ) { - - plot1D( point.x() - hit->entry().x(), hid( rich, "digErrX" ), "digit-hit diff X", -10, 10 ); - plot1D( point.y() - hit->entry().y(), hid( rich, "digErrY" ), "digit-hit diff Y", -10, 10 ); - plot1D( point.z() - hit->entry().z(), hid( rich, "digErrZ" ), "digit-hit diff Z", -10, 10 ); - plot1D( std::sqrt( ( point - hit->entry() ).Mag2() ), hid( rich, "digErrR" ), "digit-hit diff R", -10, 10 ); - - // Get MCRichOptical Photon - const auto* mcPhot = m_mcTool->mcOpticalPhoton( hit ); - if ( mcPhot ) { - // Compare position on HPD entrance window - const Gaudi::XYZPoint& mcPoint = mcPhot->pdIncidencePoint(); - plot1D( point.X() - mcPoint.X(), hid( rich, "pdImpX" ), "dX on HPD entry window : CK signal only", -30, 30 ); - plot1D( point.Y() - mcPoint.Y(), hid( rich, "pdImpY" ), "dY on HPD entry window : CK signal only", -30, 30 ); - plot1D( point.Z() - mcPoint.Z(), hid( rich, "pdImpZ" ), "dZ on HPD entry window : CK signal only", -30, 30 ); - plot1D( std::sqrt( ( point - mcPoint ).Mag2() ), hid( rich, "pdImpR" ), - "dR on HPD entry window : CK signal only", 0, 10 ); - } - } - - // Count beta=1 PEs - countNPE( ckPhotMapDig, hit ); - - // count digits from charged tracks - if ( !thisDigCounted && hit->chargedTrack() ) { thisDigCounted = true; } - - } // end hits loop - - } // MCRichDigits Loop - - // Fill multiplicity plots - - plot1D( digMult[Rich::Rich1], hid( Rich::Rich1, "digitMult" ), "Overall Digit RICH1 Occupancy", 0, - maxMult[Rich::Rich1] ); - plot1D( digMult[Rich::Rich2], hid( Rich::Rich2, "digitMult" ), "Overall Digit RICH2 Occupancy", 0, - maxMult[Rich::Rich2] ); - - plot1D( hpdCKMult[Rich::Rich1], hid( Rich::Rich1, "hpdQCKDigMult" ), "HPD Quartz CK Digit RICH1 Occupancy", 0, - maxMult[Rich::Rich1] ); - plot1D( hpdCKMult[Rich::Rich2], hid( Rich::Rich2, "hpdQCKDigMult" ), "HPD Quartz CK Digit RICH2 Occupancy", 0, - maxMult[Rich::Rich2] ); - - // Fill HPD occupancy plots - for ( const auto& occ : pdMult ) { - plot1D( occ.second, hid( ( occ.first ).rich(), "digPerPD" ), "Digits per HPD (nDigits>0)", -0.5, 200.5, 201 ); - std::ostringstream HPD; - HPD << occ.first; - plot1D( occ.second, hid( occ.first.rich(), "hpdOcc/" + (std::string)m_richSys->hardwareID( occ.first ) ), - HPD.str() + " Digits (nDigits>0)", -0.5, 200.5, 201 ); - if ( hasHPDQuartzCKBkg[occ.first] ) { - plot1D( occ.second, - hid( occ.first.rich(), "hpdOccWithHPDQCK/" + (std::string)m_richSys->hardwareID( occ.first ) ), - HPD.str() + " Digits (nDigits>0 and nHPDQCK>0)", -0.5, 200.5, 201 ); - } - } - - // Locate MCRichDeposits - // ================================================================================ - SmartDataPtr deps( eventSvc(), m_mcdepTES ); - if ( !deps ) { - _ri_debug << "Cannot locate MCRichDeposits at " << m_mcdepTES << endmsg; - } else { - _ri_debug << "Successfully located " << deps->size() << " MCRichDeposits at " << m_mcdepTES << endmsg; - - std::vector nChargedTracks( Rich::NRiches, 0 ); - - // map of backgrounds in each HPD - PartMap pmap1, pmap2, pmap3; - - // Loop over deposits - for ( const auto* dep : *deps ) { - - // Which RICH ? - const auto rich = dep->parentHit()->rich(); - - // Count beta=1 PEs - countNPE( ckPhotMapDep, dep->parentHit() ); - - // count hits from charged tracks - if ( dep->parentHit()->chargedTrack() ) ++nChargedTracks[rich]; - - // study hpd quartz window backgrounds - if ( dep->history().hpdQuartzCK() ) { - // key for this hit HPD and MCParticle pair - const PartKey key( dep->smartID().pdID(), dep->parentHit()->mcParticle() ); - // add pointer to deposit to vector for this key - pmap1[key].push_back( dep ); - if ( msgLevel( MSG::VERBOSE ) ) { - const auto locpoint = m_smartIDTool->globalToPDPanel( dep->parentHit()->entry() ); - verbose() << dep->smartID() << " glo=" << dep->parentHit()->entry() << " loc=" << locpoint - << " mcp=" << dep->parentHit()->mcParticle() << endmsg; - } - } - - // study gas quartz window backgrounds - if ( dep->history().gasQuartzCK() ) { - // key for this hit HPD and MCParticle pair - const PartKey key( dep->smartID().pdID(), dep->parentHit()->mcParticle() ); - pmap2[key].push_back( dep ); - } - - // HPD Reflections - if ( dep->history().hpdReflection() ) { - // key for this hit HPD - const PartKey key( dep->smartID().pdID(), nullptr ); - pmap3[key].push_back( dep ); - } - - } // MCRichDeposits loop - - // fill background plots - fillHPDPlots( pmap1, "HPDQuartzBkg", "HPD Quartz Window CK" ); - fillHPDPlots( pmap2, "GasQuartzBkg", "Gas Quartz Window CK" ); - fillHPDPlots( pmap3, "HPDReflectionBkg", "HPD Internal reflections" ); - - } // MCRichDeposits exist - - // Locate MCRichHits - // ================================================================================ - SmartDataPtr hits( eventSvc(), m_mchitTES ); - if ( !hits ) { - _ri_debug << "Cannot locate MCRichHits at " << m_mchitTES << endmsg; - } else { - - // Hit mult counters - std::vector hitMult( Rich::NRiches, 0 ); - std::vector hpdQCKMult( Rich::NRiches, 0 ); - - // Loop over hits - for ( const auto* hit : *hits ) { - - // Which RICH ? - const auto rich = hit->rich(); - - if ( trueCKHit( hit ) ) { - plot1D( hit->timeOfFlight(), hid( rich, "tofHitSignal" ), "Signal MCRichHit TOF", -100, 100 ); - plot1D( hit->energy(), hid( rich, "enHitSignal" ), "Signal MCRichHit Energy", 0, 100 ); - } else { - plot1D( hit->timeOfFlight(), hid( rich, "tofHitBkgrd" ), "Background MCRichHit TOF", -100, 100 ); - plot1D( hit->energy(), hid( rich, "enHitBkgrd" ), "Background MCRichHit Energy", 0, 100 ); - } - - // increment hit count - ++hitMult[rich]; - if ( hit->hpdQuartzCK() ) { ++hpdQCKMult[rich]; } - - // Increment PES count for high beta tracks - countNPE( ckPhotMapHit, hit ); - - // Plot hit positions - const auto& point = hit->entry(); - plot1D( point.x(), hid( rich, "hitXGlo" ), "mchits x global", xMinPDGlo[rich], xMaxPDGlo[rich] ); - plot1D( point.y(), hid( rich, "hitYGlo" ), "mchits y global", yMinPDGlo[rich], yMaxPDGlo[rich] ); - plot1D( point.z(), hid( rich, "hitZGlo" ), "mchits z global", zMinPDGlo[rich], zMaxPDGlo[rich] ); - plot2D( point.x(), point.y(), hid( rich, "hitXYglo" ), "mchits yVx global", xMinPDGlo[rich], xMaxPDGlo[rich], - yMinPDGlo[rich], yMaxPDGlo[rich] ); - plot2D( point.z(), point.x(), hid( rich, "hitZXglo" ), "mchits xVz global", zMinPDGlo[rich], zMaxPDGlo[rich], - xMinPDGlo[rich], xMaxPDGlo[rich] ); - plot2D( point.z(), point.y(), hid( rich, "hitZYglo" ), "mchits yVz global", zMinPDGlo[rich], zMaxPDGlo[rich], - yMinPDGlo[rich], yMaxPDGlo[rich] ); - } - - // Fill multiplicity plots - - plot1D( hitMult[Rich::Rich1], hid( Rich::Rich1, "hitMult" ), "Overall MCRichHit RICH1 Occupancy", 0, - maxMult[Rich::Rich1] ); - plot1D( hitMult[Rich::Rich2], hid( Rich::Rich2, "hitMult" ), "Overall MCRichHit RICH2 Occupancy", 0, - maxMult[Rich::Rich2] ); - - plot1D( hpdQCKMult[Rich::Rich1], hid( Rich::Rich1, "hpdQCKHitMult" ), "HPD Quartz CK MCRichHit RICH1 Occupancy", 0, - maxMult[Rich::Rich1] ); - plot1D( hpdQCKMult[Rich::Rich2], hid( Rich::Rich2, "hpdQCKHitMult" ), "HPD Quartz CK MCRichHit RICH2 Occupancy", 0, - maxMult[Rich::Rich2] ); - - } // MCRichHits exist - - // Loop over counted PES for MCRichHits - for ( const auto& phot : ckPhotMapHit ) { - plot1D( phot.second, hid( phot.first.second, "npesHit" ), "# MCRichHit p.e.s : beta=1", -0.5, 100.5, 101 ); - } - - // Loop over counted PES for MCRichDeposits - for ( const auto& phot : ckPhotMapDep ) { - plot1D( phot.second, hid( phot.first.second, "npesDep" ), "# MCRichDeposit p.e.s : beta=1", -0.5, 100.5, 101 ); - } - - // Loop over counted PES for MCRichDigits - for ( const auto& phot : ckPhotMapDig ) { - const auto digCount = phot.second; - plot1D( digCount, hid( phot.first.second, "npesDig" ), "# MCRichDigit p.e.s : beta=1", -0.5, 100.5, 101 ); - - // locate the entry for the MCRichDeps - const auto depCount = ckPhotMapDep[phot.first]; - - const auto frac = ( depCount != 0 ? digCount / depCount : 0 ); - plot1D( frac, hid( phot.first.second, "npesRetained" ), "Retained p.e.s : beta=1", 0, 1 ); - } - - return StatusCode::SUCCESS; -} - -void AlgMoni::fillHPDPlots( const PartMap& pmap, const std::string& plotsDir, const std::string& plotsName ) { - // histogramming - PD_LOCAL_POSITIONS_X; - PD_LOCAL_POSITIONS_Y; - const Rich::HistoID hid; - - // count for each RICH - std::vector nHPDs( Rich::NRiches, 0 ); - - _ri_debug << "Filling plots for " << plotsName << endmsg; - - // plots for each HPD - for ( const auto& P : pmap ) { - // SmartID for this HPD - const auto hpdID = P.first.first; - // increment count for RICH - ++nHPDs[hpdID.rich()]; - // Centre point of HPD in local coords - Gaudi::XYZPoint hpdGlo; - const auto sc = m_smartIDTool->pdPosition( hpdID, hpdGlo ); - if ( !sc ) continue; - const auto hpdLoc = m_smartIDTool->globalToPDPanel( hpdGlo ); - // create histo title - std::ostringstream HPD; - HPD << plotsName << " " << hpdID << " " << hpdGlo; - // histo ID - std::string Hid = plotsDir + "/" + boost::lexical_cast( m_ID++ ); - // loop over deposits - Rich::Map uniqueIDs; - _ri_debug << "Found " << P.second.size() << " Deposits for " << Hid << " " << HPD.str() << endmsg; - for ( const auto* D : P.second ) { - uniqueIDs[D->smartID()] = true; - // hit point on silicon in local coords - const auto hitP = m_smartIDTool->globalToPDPanel( D->parentHit()->entry() ) - hpdLoc; - // fill single event HPD plot - if ( m_ID < m_maxID ) { - plot2D( hitP.X(), hitP.Y(), Hid, HPD.str(), -8 * Gaudi::Units::mm, 8 * Gaudi::Units::mm, -8 * Gaudi::Units::mm, - 8 * Gaudi::Units::mm, 32, 32 ); - } - // Intergrated plot - plot2D( hitP.X(), hitP.Y(), plotsDir + "/" + "intBkg", "Integrated " + plotsName, -8 * Gaudi::Units::mm, - 8 * Gaudi::Units::mm, -8 * Gaudi::Units::mm, 8 * Gaudi::Units::mm, 32, 32 ); - } - // plot number of hits - plot1D( P.second.size(), plotsDir + "/" + hid( hpdID.rich(), "/nHitsPerHPD" ).id(), - "# hits per HPD from " + plotsName, 0.5, 200.5, 100 ); - // plot number of digits - plot1D( uniqueIDs.size(), plotsDir + "/" + hid( hpdID.rich(), "/nDigsPerHPD" ).id(), - "# digits per HPD from " + plotsName, 0.5, 200.5, 100 ); - // position of HPD - plot2D( hpdLoc.X(), hpdLoc.Y(), plotsDir + "/" + hid( hpdID.rich(), "hitHPDs" ).id(), "HPD Hit Map", - xMinPDLoc[hpdID.rich()], xMaxPDLoc[hpdID.rich()], yMinPDLoc[hpdID.rich()], yMaxPDLoc[hpdID.rich()] ); - } - - // number of affected HPD plots - plot1D( nHPDs[Rich::Rich1], plotsDir + "/" + hid( Rich::Rich1, "/nHPDsPerEv" ).id(), - "# RICH1 HPDs per event with " + plotsName, -0.5, 20.5, 21 ); - plot1D( nHPDs[Rich::Rich2], plotsDir + "/" + hid( Rich::Rich2, "/nHPDsPerEv" ).id(), - "# RICH2 HPDs per event with " + plotsName, -0.5, 20.5, 21 ); -} - -void AlgMoni::countNPE( PhotMap& photMap, const LHCb::MCRichHit* hit ) { - - // MCParticle momentum - const auto tkPtot = momentum( hit->mcParticle() ); - - // Parent particle hypothesis - const auto mcid = m_mcTool->mcParticleType( hit->mcParticle() ); - - // Which radiator - const auto rad = ( Rich::RadiatorType )( hit->radiator() ); - - // Increment PES count for high beta tracks - if ( tkPtot > 1 * Gaudi::Units::GeV && mcid != Rich::Unknown && mcid != Rich::Electron && hit->isSignal() && - mcBeta( hit->mcParticle() ) > 0.99 ) { - PhotPair pairC( hit->mcParticle(), rad ); - ++photMap[pairC]; - } -} diff --git a/Rich/RichReadout/src/component/RichDigiAlgMoni.h b/Rich/RichReadout/src/component/RichDigiAlgMoni.h deleted file mode 100755 index b6127ebea..000000000 --- a/Rich/RichReadout/src/component/RichDigiAlgMoni.h +++ /dev/null @@ -1,155 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration * -* * -* This software is distributed under the terms of the GNU General Public * -* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". * -* * -* In applying this licence, CERN does not waive the privileges and immunities * -* granted to it by virtue of its status as an Intergovernmental Organization * -* or submit itself to any jurisdiction. * -\*****************************************************************************/ -#ifndef RICHMONITOR_RICHDIGIALGMONI_H -#define RICHMONITOR_RICHDIGIALGMONI_H 1 - -// STD -#include - -// base class -#include "RichKernel/RichHistoAlgBase.h" - -// from Gaudi -#include "GaudiKernel/Vector3DTypes.h" - -// Event model -#include "Event/MCRichDeposit.h" -#include "Event/MCRichDigit.h" -#include "Event/MCRichHit.h" -#include "Event/MCRichOpticalPhoton.h" -#include "Event/MCRichSummedDeposit.h" -#include "Kernel/RichParticleIDType.h" - -// Utils -#include "RichUtils/RichHashMap.h" -#include "RichUtils/RichMap.h" - -// Kernel -#include "Kernel/RichDetectorType.h" -#include "Kernel/RichParticleIDType.h" -#include "Kernel/RichRadiatorType.h" - -// GSL -#include "gsl/gsl_math.h" - -// interfaces -#include "MCInterfaces/IRichMCTruthTool.h" -#include "RichInterfaces/IRichSmartIDTool.h" - -// boost -#include "boost/lexical_cast.hpp" - -// RichDet -#include "RichDet/DeRichSystem.h" - -// temporary histogramming numbers -#include "RichDetParams.h" - -namespace Rich { - namespace MC { - namespace Digi { - - /** @class AlgMoni RichDigiAlgMoni.h - * - * Monitor for Rich Digitisation algorithm performance - * - * @author Chris Jones (Christopher.Rob.Jones@cern.ch) - * @date 2003-09-08 - */ - - class AlgMoni final : public HistoAlgBase { - - public: - /// Standard constructor - AlgMoni( const std::string& name, ISvcLocator* pSvcLocator ); - - StatusCode initialize() override final; ///< Algorithm initialization - StatusCode execute() override final; ///< Algorithm execution - - private: // definitions - /// key for HPD+MCParticle hits - typedef std::pair PartKey; - /// mapping given hits for each HPD+MCParticle Pair - typedef std::map> PartMap; - - private: // methods - // Map to count cherenkov photons for each radiator - typedef std::pair PhotPair; - typedef Rich::Map PhotMap; - - // PD occupancies - typedef Rich::Map PDMulti; - - /// Particle mass - double mass( const LHCb::MCParticle* mcPart ); - - /// Returns beta for a given MCParticle - double mcBeta( const LHCb::MCParticle* mcPart ); - - /// Returns the momentum for a given MCParticle - double momentum( const LHCb::MCParticle* mcPart ); - - /// Count the number of photo electrons - void countNPE( PhotMap& photMap, const LHCb::MCRichHit* hit ); - - /// IS this a true cherenkov signal hit ? - bool trueCKHit( const LHCb::MCRichHit* hit ); - - /// Fill histograms for each HPD silicon wafer background plots - void fillHPDPlots( const PartMap& pmap, const std::string& plotsDir, const std::string& plotsName ); - - private: // data - /// Pointer to RichSmartID tool - const Rich::ISmartIDTool* m_smartIDTool = nullptr; - - /// Pointer to MC truth tool - const Rich::MC::IMCTruthTool* m_mcTool = nullptr; - - /// Pointer to RICH system detector element - const DeRichSystem* m_richSys = nullptr; - - // job options - std::string m_mcdigitTES; ///< Location of MCRichDigits in TES - std::string m_mcdepTES; ///< Location of MCRichDeposits in TES - std::string m_mchitTES; ///< Location of MCRichHits in TES - - /// Particle masses - Rich::Map m_particleMass; - - // histo ID - unsigned int m_ID; - - // max number of Quartz CK HPD event histos to make - unsigned int m_maxID; - }; - - inline double AlgMoni::mass( const LHCb::MCParticle* mcPart ) { - return m_particleMass[m_mcTool->mcParticleType( mcPart )]; - } - - inline double AlgMoni::momentum( const LHCb::MCParticle* mcPart ) { - return ( mcPart ? mcPart->momentum().P() : 0 ); - } - - inline double AlgMoni::mcBeta( const LHCb::MCParticle* mcPart ) { - if ( !mcPart ) return 0; - const double pTot = momentum( mcPart ); - const double Esquare = pTot * pTot + gsl_pow_2( mass( mcPart ) ); - return ( Esquare > 0 ? pTot / sqrt( Esquare ) : 0 ); - } - - inline bool AlgMoni::trueCKHit( const LHCb::MCRichHit* hit ) { return hit->isSignal(); } - - } // namespace Digi - } // namespace MC -} // namespace Rich - -#endif // RICHMONITOR_RICHDIGIALGMONI_H -- GitLab From b42bc1e1c13c7424789194aa1f01c1994b72893d Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Tue, 4 Oct 2022 08:24:32 +0200 Subject: [PATCH 02/15] Converted UTTightDigitCreator to use derived conditions --- .../src/UTTightDigitCreator.cpp | 75 ++++++++----------- 1 file changed, 32 insertions(+), 43 deletions(-) diff --git a/UT/UTDigiAlgorithms/src/UTTightDigitCreator.cpp b/UT/UTDigiAlgorithms/src/UTTightDigitCreator.cpp index a9884a1cc..ac3c1c097 100644 --- a/UT/UTDigiAlgorithms/src/UTTightDigitCreator.cpp +++ b/UT/UTDigiAlgorithms/src/UTTightDigitCreator.cpp @@ -18,10 +18,15 @@ #include "LHCbAlgs/Transformer.h" -#ifdef USE_DD4HEP -// needed by getDet/existDet API. Should be dropped when made properly functional -# include "DetDesc/Condition.h" -#endif +#include + +namespace { + struct Conditions { + double digitSig2NoiseThreshold{3.0}; + double clusterSig2NoiseThreshold{4.0}; + double highThreshold{10.0}; + }; +} // namespace /** * Class for clustering in the UT tracker @@ -29,27 +34,22 @@ * @author A.Beiter (based on code by M.Needham) * @date 2018-09-04 */ - -class UTTightDigitCreator - : public LHCb::Algorithm::Transformer> { +class UTTightDigitCreator : public LHCb::Algorithm::Transformer< + LHCb::UTDigits( const LHCb::UTDigits&, DeUTDetector const&, Conditions const& ), + LHCb::DetDesc::usesBaseAndConditions> { public: UTTightDigitCreator( const std::string& name, ISvcLocator* pSvcLocator ); StatusCode initialize() override; - LHCb::UTDigits operator()( const LHCb::UTDigits&, DeUTDetector const& ) const override; + LHCb::UTDigits operator()( const LHCb::UTDigits&, DeUTDetector const&, Conditions const& ) const override; private: - bool aboveDigitSignalToNoise( const LHCb::UTDigit* aDigit, const DeUTSector& aSector ) const; + bool aboveDigitSignalToNoise( const LHCb::UTDigit* aDigit, const DeUTSector& aSector, + double digitSig2NoiseThreshold ) const; StatusCode loadCutsFromConditions(); - Gaudi::Property m_forceOptions{this, "forceOptions", false}; - Gaudi::Property m_digitSig2NoiseThreshold{this, "DigitSignal2Noise", 3.0}; - Gaudi::Property m_clusterSig2NoiseThreshold{this, "ClusterSignal2Noise", 4.0}; - Gaudi::Property m_highThreshold{this, "HighSignal2Noise", 10.0}; - Gaudi::Property m_maxSize{this, "Size", 4}; - + Gaudi::Property m_maxSize{this, "Size", 4}; Gaudi::Property m_conditionLocation{this, "conditionLocation", "/dd/Conditions/ReadoutConf/UT/ClusteringThresholds"}; @@ -67,27 +67,29 @@ DECLARE_COMPONENT( UTTightDigitCreator ) UTTightDigitCreator::UTTightDigitCreator( const std::string& name, ISvcLocator* pSvcLocator ) : Transformer{name, pSvcLocator, - {{"InputLocation", UTDigitLocation::UTDigits}, {"UTLocation", DeUTDetLocation::location()}}, + {{"InputLocation", UTDigitLocation::UTDigits}, + {"UTLocation", DeUTDetLocation::location()}, + {"Conditions", name + "_Conditions"}}, {"OutputLocation", UTDigitLocation::UTTightDigits}} {} StatusCode UTTightDigitCreator::initialize() { return Transformer::initialize().andThen( [&]() -> StatusCode { - // get the cuts from condition if present - if ( !existDet( m_conditionLocation ) || m_forceOptions ) { - info() << "Default to jobOptions for cuts" << endmsg; - return StatusCode::SUCCESS; - } else { - registerCondition( m_conditionLocation, &UTTightDigitCreator::loadCutsFromConditions ); - return runUpdate(); - } + addConditionDerivation( {m_conditionLocation}, inputLocation(), + [&]( YAML::Node const& cInfo ) -> Conditions { + info() << "Loading cuts tagged as " << cInfo["tag"].as() << endmsg; + return {cInfo["HitThreshold"].as(), cInfo["ConfirmationThreshold"].as(), + cInfo["SpilloverThreshold"].as()}; + } ); + return StatusCode::SUCCESS; } ); } -LHCb::UTDigits UTTightDigitCreator::operator()( const UTDigits& digitCont, DeUTDetector const& det ) const { +LHCb::UTDigits UTTightDigitCreator::operator()( const UTDigits& digitCont, DeUTDetector const& det, + Conditions const& conditions ) const { UTDigits digitCont_ADCcut; for ( const auto& digit : digitCont ) { auto aSector = det.findSector( digit->channelID() ); - if ( !aboveDigitSignalToNoise( digit, *aSector ) ) continue; + if ( !aboveDigitSignalToNoise( digit, *aSector, conditions.digitSig2NoiseThreshold ) ) continue; // link to the DAQ UTDAQID utdaqID = m_readoutTool->channelIDToDAQID( (const Detector::UT::ChannelID)digit->channelID() ); @@ -102,24 +104,11 @@ LHCb::UTDigits UTTightDigitCreator::operator()( const UTDigits& digitCont, DeUTD digitCont_ADCcut.insert( new UTDigit( digit->channelID(), digit->depositedCharge(), utdaqID.id() ), digit->channelID() ); } - return digitCont_ADCcut; } -bool UTTightDigitCreator::aboveDigitSignalToNoise( const UTDigit* aDigit, const DeUTSector& aSector ) const { - const unsigned int threshold = LHCb::Math::round( m_digitSig2NoiseThreshold * aSector.noise( aDigit->channelID() ) ); +bool UTTightDigitCreator::aboveDigitSignalToNoise( const UTDigit* aDigit, const DeUTSector& aSector, + double digitSig2NoiseThreshold ) const { + const unsigned int threshold = LHCb::Math::round( digitSig2NoiseThreshold * aSector.noise( aDigit->channelID() ) ); return aDigit->depositedCharge() > threshold; } - -StatusCode UTTightDigitCreator::loadCutsFromConditions() { - - // load conditions - auto* cInfo = getDet( m_conditionLocation ); - info() << "Loading cuts tagged as " << cInfo->param( "tag" ) << endmsg; - - m_digitSig2NoiseThreshold = cInfo->param( "HitThreshold" ); - m_clusterSig2NoiseThreshold = cInfo->param( "ConfirmationThreshold" ); - m_highThreshold = cInfo->param( "SpilloverThreshold" ); - - return StatusCode::SUCCESS; -} -- GitLab From 64e873186b7b387f3a516133120ecbdc927ae3fa Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Tue, 4 Oct 2022 08:25:08 +0200 Subject: [PATCH 03/15] Modernized MCFTG4AttenuationInterpolationTool and switched to derived conditions --- .../MCFTG4AttenuationInterpolationTool.cpp | 227 ++++++++++-------- .../src/MCFTG4AttenuationInterpolationTool.h | 64 ----- 2 files changed, 126 insertions(+), 165 deletions(-) delete mode 100644 FT/FTDigitisation/src/MCFTG4AttenuationInterpolationTool.h diff --git a/FT/FTDigitisation/src/MCFTG4AttenuationInterpolationTool.cpp b/FT/FTDigitisation/src/MCFTG4AttenuationInterpolationTool.cpp index 140285af2..66a24209c 100644 --- a/FT/FTDigitisation/src/MCFTG4AttenuationInterpolationTool.cpp +++ b/FT/FTDigitisation/src/MCFTG4AttenuationInterpolationTool.cpp @@ -8,15 +8,79 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ -#include "MCFTG4AttenuationInterpolationTool.h" -#include + +#include "IMCFTAttenuationTool.h" // Base Class + +#include "Kernel/STLExtensions.h" +#include "Kernel/TaggedBool.h" +#include + +#include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/GaudiException.h" +#include "GaudiKernel/SystemOfUnits.h" + #include -//----------------------------------------------------------------------------- -// Implementation file for class : MCFTG4AttenuationInterpolationTool -// -// 2018-28-08 : Martin Bieker -//----------------------------------------------------------------------------- +#include + +namespace { + enum class ParameterName { mask, a_refl, a_dir, b_refl, b_dir, c_refl, c_dir, d_refl, d_dir, y_pos }; + ParameterName convert( std::string const& key ) { + if ( key == "mask" ) return ParameterName::mask; + if ( key == "a_refl" ) return ParameterName::a_refl; + if ( key == "a_dir" ) return ParameterName::a_dir; + if ( key == "b_refl" ) return ParameterName::b_refl; + if ( key == "b_dir" ) return ParameterName::b_dir; + if ( key == "c_refl" ) return ParameterName::c_refl; + if ( key == "c_dir" ) return ParameterName::c_dir; + if ( key == "d_refl" ) return ParameterName::d_refl; + if ( key == "d_dir" ) return ParameterName::d_dir; + if ( key == "y_pos" ) + return ParameterName::y_pos; + else + throw std::runtime_error( "unknown ParameterName " + key ); + } + + struct Conditions { + unsigned int nBinsX, nBinsY; + std::vector xEdges, yEdges, effDir, effRef; + std::map> parameters; + }; +} // namespace + +/** + * Tool that reads LYAM parameters from the CondDB and creates a interpolated attenuation map + * for a given irradiation dose and (time) age of the fibres. + * + * @author Martin Bieker based on implementation of M. Demmer and J. Wishahi + * @date 2018-28-08 + */ +class MCFTG4AttenuationInterpolationTool + : public extends, IMCFTAttenuationTool> { + +public: + using extends::extends; + StatusCode initialize() override; + + /// Calculate the direct attenuation and the attenuation with reflection + Attenuation attenuation( double x, double y ) const override; + +protected: + using ReflectedPhotons = LHCb::tagged_bool; + double calculateAttenuation( int binID, MCFTG4AttenuationInterpolationTool::ReflectedPhotons isRefl, + Conditions const& conds ); + +private: + Conditions processParameterMap( YAML::Node const& cond, double mirrorReflectivity ); + int findBin( LHCb::span axis, double position ) const; + Gaudi::Property m_parameterLocation{this, "ParameterLocation", + "/dd/Conditions/Calibration/FT/AttenuationInterpolation", + "Location of ParameterFile in CondDB"}; + Gaudi::Property m_fibreAge{this, "FibreAge", 0, "Age of fibres in months"}; + Gaudi::Property m_irrad{this, "IrradiationLevel", 0, "Irradiation level in 1/fb"}; + Gaudi::Property m_mirrorReflectivity{this, "MirrorReflectivity", 0.75, "Reflectivity of the Mirror"}; + ConditionAccessor m_conditions{this, name() + "_Conditions"}; +}; // names of the parameters describing the binning of the LYAMs namespace { @@ -26,96 +90,71 @@ namespace { // Declaration of the Tool Factory DECLARE_COMPONENT( MCFTG4AttenuationInterpolationTool ) -//========================================================================= -// Retrieve the Attenuation map from the CondDB -//========================================================================= StatusCode MCFTG4AttenuationInterpolationTool::initialize() { - auto sc = base_class::initialize(); - if ( sc.isFailure() ) { return sc; } - info() << "Using the FT attenuation interpolation tool with L= " << std::to_string( m_irrad ) << "/fb " << endmsg; - // When conditions are not found, return failure - if ( existDet( m_parameterLocation ) == false ) { - return Error( "No parameter maps found in SIMCOND: update to a newer tag", StatusCode::FAILURE ); - } - - processParameterMap( m_parameterLocation ); - - if ( !validateMap() ) { - return Error( "Consitency check of parameter map failed: Please check SIMCOND tag", StatusCode::FAILURE ); - } - fillMaps(); - - if ( msgLevel( MSG::DEBUG ) ) { debug() << "nBinsX :" << m_nBinsX << " nBinsY: " << m_nBinsY << std::endl; } - - return sc; + return extends::initialize().andThen( [&] { + info() << "Using the FT attenuation interpolation tool with L= " << std::to_string( m_irrad ) << "/fb " << endmsg; + addConditionDerivation( {m_parameterLocation}, m_conditions.key(), [&]( YAML::Node const& cond ) -> Conditions { + return this->processParameterMap( cond, m_mirrorReflectivity ); + } ); + return StatusCode::SUCCESS; + } ); } -//============================================================================= // Return calculated attenution to MCFTDepositCreator -//============================================================================= IMCFTAttenuationTool::Attenuation MCFTG4AttenuationInterpolationTool::attenuation( double hitXPosition, double hitYPosition ) const { hitXPosition = std::abs( hitXPosition ); hitYPosition = std::abs( hitYPosition ); - int binX = findBin( m_xEdges, hitXPosition ); - int binY = findBin( m_yEdges, hitYPosition ); + auto& conds = m_conditions.get(); + int binX = findBin( conds.xEdges, hitXPosition ); + int binY = findBin( conds.yEdges, hitYPosition ); // Hit outside quadrant (should never happen!) if ( binX == -1 || binY == -1 ) { return {0.0, 0.0}; } else { - return {m_effDir[binY * m_nBinsX + binX], m_effRef[binY * m_nBinsX + binX]}; + return {conds.effDir[binY * conds.nBinsX + binX], conds.effRef[binY * conds.nBinsX + binX]}; } } -MCFTG4AttenuationInterpolationTool::ParameterName -MCFTG4AttenuationInterpolationTool::convert( std::string const& key ) const { - if ( key == "mask" ) return ParameterName::mask; - if ( key == "a_refl" ) return ParameterName::a_refl; - if ( key == "a_dir" ) return ParameterName::a_dir; - if ( key == "b_refl" ) return ParameterName::b_refl; - if ( key == "b_dir" ) return ParameterName::b_dir; - if ( key == "c_refl" ) return ParameterName::c_refl; - if ( key == "c_dir" ) return ParameterName::c_dir; - if ( key == "d_refl" ) return ParameterName::d_refl; - if ( key == "d_dir" ) return ParameterName::d_dir; - if ( key == "y_pos" ) - return ParameterName::y_pos; - else - throw std::runtime_error( "unknown ParameterName " + key ); -} - -//============================================================================= // Read in and validate parameter map from CondDB -//============================================================================= -void MCFTG4AttenuationInterpolationTool::processParameterMap( const std::string& conditionsLocation ) { - - // Get the conditions - Condition* cond = getDet( conditionsLocation ); - +Conditions MCFTG4AttenuationInterpolationTool::processParameterMap( YAML::Node const& cond, + double mirrorReflectivity ) { + Conditions outConds; // Process binning information - m_nBinsX = cond->param( "x_n_bins" ); - m_nBinsY = cond->param( "y_n_bins" ); - m_xEdges = cond->param>( "x_edges" ); - m_yEdges = cond->param>( "y_edges" ); - std::sort( m_xEdges.begin(), m_xEdges.end() ); - std::sort( m_yEdges.begin(), m_yEdges.end() ); + outConds.nBinsX = cond["x_n_bins"].as(); + outConds.nBinsY = cond["y_n_bins"].as(); + outConds.xEdges = cond["x_edges"].as>(); + outConds.yEdges = cond["y_edges"].as>(); + std::sort( begin( outConds.xEdges ), end( outConds.xEdges ) ); + std::sort( begin( outConds.yEdges ), end( outConds.yEdges ) ); // Loop over maps of parameters - for ( auto parameterName : cond->paramNames() ) { - if ( s_binParams.count( parameterName ) != 0 ) { continue; } - m_parameters[convert( parameterName )] = cond->param>( parameterName ); - } // for -} + for ( auto it = cond.begin(); it != cond.end(); it++ ) { + auto name = it->first.as(); + if ( s_binParams.count( name ) != 0 ) { continue; } + outConds.parameters[convert( name )] = it->second.as>(); + } -bool MCFTG4AttenuationInterpolationTool::validateMap() { - if ( m_nBinsX != m_xEdges.size() - 1 || m_nBinsY != m_yEdges.size() - 1 ) { - return false; - } else { - return std::none_of( m_parameters.begin(), m_parameters.end(), - [nb = m_nBinsX * m_nBinsY]( const auto& i ) { return nb != std::get<1>( i ).size(); } ); + // check for consistency + if ( outConds.nBinsX != outConds.xEdges.size() - 1 || outConds.nBinsY != outConds.yEdges.size() - 1 || + std::any_of( + outConds.parameters.begin(), outConds.parameters.end(), + [nb = outConds.nBinsX * outConds.nBinsY]( const auto& i ) { return nb != std::get<1>( i ).size(); } ) ) { + throw GaudiException( "MCFTG4AttenuationInterpolationTool::processParameterMap", + "Consitency check of parameter map failed: Please check SIMCOND tag", StatusCode::FAILURE ); + } + + // fill maps + outConds.effDir.reserve( outConds.nBinsX * outConds.nBinsY ); + outConds.effRef.reserve( outConds.nBinsX * outConds.nBinsY ); + for ( unsigned int bin = 0; bin < outConds.nBinsX * outConds.nBinsY; ++bin ) { + outConds.effDir.push_back( calculateAttenuation( bin, ReflectedPhotons{false}, outConds ) ); + outConds.effRef.push_back( calculateAttenuation( bin, ReflectedPhotons{true}, outConds ) * mirrorReflectivity ); } + + return outConds; } int MCFTG4AttenuationInterpolationTool::findBin( LHCb::span axis, double position ) const { @@ -131,36 +170,22 @@ int MCFTG4AttenuationInterpolationTool::findBin( LHCb::span axis, return iterator - axis.begin() - 1; } -void MCFTG4AttenuationInterpolationTool::fillMaps() { - // Map can only be filled once - assert( m_effDir.empty() ); - assert( m_effRef.empty() ); - - m_effDir.reserve( m_nBinsX * m_nBinsY ); - m_effRef.reserve( m_nBinsX * m_nBinsY ); - - for ( unsigned int bin = 0; bin < m_nBinsX * m_nBinsY; ++bin ) { - m_effDir.push_back( calculateAttenuation( bin, ReflectedPhotons{false} ) ); - m_effRef.push_back( calculateAttenuation( bin, ReflectedPhotons{true} ) * m_mirrorReflectivity ); - } -} - double MCFTG4AttenuationInterpolationTool::calculateAttenuation( - int binID, MCFTG4AttenuationInterpolationTool::ReflectedPhotons isRefl ) { + int binID, MCFTG4AttenuationInterpolationTool::ReflectedPhotons isRefl, Conditions const& conds ) { if ( isRefl ) { - return m_parameters[ParameterName::mask][binID] * - ( ( m_parameters[ParameterName::a_refl][binID] * m_irrad + - +m_parameters[ParameterName::b_refl][binID] * - exp( -m_parameters[ParameterName::c_refl][binID] * m_irrad ) + - m_parameters[ParameterName::d_refl][binID] ) * - exp( -4.2e-4 * ( m_parameters[ParameterName::y_pos][binID] + 2.5 ) * m_fibreAge ) ); + return conds.parameters.at( ParameterName::mask )[binID] * + ( ( conds.parameters.at( ParameterName::a_refl )[binID] * m_irrad + + +conds.parameters.at( ParameterName::b_refl )[binID] * + exp( -conds.parameters.at( ParameterName::c_refl )[binID] * m_irrad ) + + conds.parameters.at( ParameterName::d_refl )[binID] ) * + exp( -4.2e-4 * ( conds.parameters.at( ParameterName::y_pos )[binID] + 2.5 ) * m_fibreAge ) ); } else { - return m_parameters[ParameterName::mask][binID] * - ( ( m_parameters[ParameterName::a_dir][binID] * m_irrad + - +m_parameters[ParameterName::b_dir][binID] * - exp( -m_parameters[ParameterName::c_dir][binID] * m_irrad ) + - m_parameters[ParameterName::d_dir][binID] ) * - exp( -4.2e-4 * ( 2.5 - m_parameters[ParameterName::y_pos][binID] ) * m_fibreAge ) ); + return conds.parameters.at( ParameterName::mask )[binID] * + ( ( conds.parameters.at( ParameterName::a_dir )[binID] * m_irrad + + +conds.parameters.at( ParameterName::b_dir )[binID] * + exp( -conds.parameters.at( ParameterName::c_dir )[binID] * m_irrad ) + + conds.parameters.at( ParameterName::d_dir )[binID] ) * + exp( -4.2e-4 * ( 2.5 - conds.parameters.at( ParameterName::y_pos )[binID] ) * m_fibreAge ) ); } } diff --git a/FT/FTDigitisation/src/MCFTG4AttenuationInterpolationTool.h b/FT/FTDigitisation/src/MCFTG4AttenuationInterpolationTool.h deleted file mode 100644 index fba51fe4e..000000000 --- a/FT/FTDigitisation/src/MCFTG4AttenuationInterpolationTool.h +++ /dev/null @@ -1,64 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2018 CERN for the benefit of the LHCb Collaboration * -* * -* This software is distributed under the terms of the GNU General Public * -* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". * -* * -* In applying this licence, CERN does not waive the privileges and immunities * -* granted to it by virtue of its status as an Intergovernmental Organization * -* or submit itself to any jurisdiction. * -\*****************************************************************************/ -#ifndef MCFTG4ATTENUATIONINTERPOLATIONTOOL_H -#define MCFTG4ATTENUATIONINTERPOLATIONTOOL_H 1 - -// from Gaudi -#include "GaudiAlg/GaudiTool.h" -#include "GaudiKernel/SystemOfUnits.h" -#include "Kernel/TaggedBool.h" - -#include "IMCFTAttenuationTool.h" // Base Class - -#include "Kernel/STLExtensions.h" - -/** @class MCFTG4AttenuationInterpolationTool MCFTG4AttenuationInterpolationTool.h - * - * Tool that reads LYAM parameters from the CondDB and creates a interpolated attenuation map - * for a given irradiation dose and (time) age of the fibres. - * - * @author Martin Bieker based on implementation of M. Demmer and J. Wishahi - * @date 2018-28-08 - */ - -class MCFTG4AttenuationInterpolationTool : public extends { - -public: - using extends::extends; - - /// Initialize the LYAM - StatusCode initialize() override; - - /// Calculate the direct attenuation and the attenuation with reflection - Attenuation attenuation( double x, double y ) const override; - -protected: - using ReflectedPhotons = LHCb::tagged_bool; - double calculateAttenuation( int binID, MCFTG4AttenuationInterpolationTool::ReflectedPhotons isRefl ); - unsigned int m_nBinsX, m_nBinsY; - std::vector m_xEdges, m_yEdges, m_effDir, m_effRef; - enum class ParameterName { mask, a_refl, a_dir, b_refl, b_dir, c_refl, c_dir, d_refl, d_dir, y_pos }; - ParameterName convert( std::string const& key ) const; - std::map> m_parameters; - -private: - void processParameterMap( const std::string& conditionsLocation ); - bool validateMap(); - void fillMaps(); - int findBin( LHCb::span axis, double position ) const; - Gaudi::Property m_parameterLocation{this, "ParameterLocation", - "/dd/Conditions/Calibration/FT/AttenuationInterpolation", - "Location of ParameterFile in CondDB"}; - Gaudi::Property m_fibreAge{this, "FibreAge", 0, "Age of fibres in months"}; - Gaudi::Property m_irrad{this, "IrradiationLevel", 0, "Irradiation level in 1/fb"}; - Gaudi::Property m_mirrorReflectivity{this, "MirrorReflectivity", 0.75, "Reflectivity of the Mirror"}; -}; -#endif // MCFG4TATTENUATIONTOOL_H -- GitLab From 19d9005ddd9d58ba4468bf355ca4399242bf7c4f Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Tue, 4 Oct 2022 08:44:56 +0200 Subject: [PATCH 04/15] Modernized slightly MCFTAttenuationTool --- FT/FTDigitisation/src/IMCFTAttenuationTool.h | 10 +- FT/FTDigitisation/src/MCFTAttenuationTool.cpp | 247 +++++++++++------- FT/FTDigitisation/src/MCFTAttenuationTool.h | 107 -------- 3 files changed, 161 insertions(+), 203 deletions(-) delete mode 100644 FT/FTDigitisation/src/MCFTAttenuationTool.h diff --git a/FT/FTDigitisation/src/IMCFTAttenuationTool.h b/FT/FTDigitisation/src/IMCFTAttenuationTool.h index cbc4f9bb0..f45bc0222 100644 --- a/FT/FTDigitisation/src/IMCFTAttenuationTool.h +++ b/FT/FTDigitisation/src/IMCFTAttenuationTool.h @@ -8,22 +8,17 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ -#ifndef IMCFTATTENUATIONTOOL_H -#define IMCFTATTENUATIONTOOL_H 1 +#pragma once -// from Gaudi #include "GaudiKernel/IAlgTool.h" -/** @class IMCFTAttenuationTool IMCFTAttenuationTool.h - * +/** * Interface for the tool that calculates the light attenuation for the MCFTDepositCreator * * @author Michel De Cian * @date 2015-01-15 */ - struct IMCFTAttenuationTool : extend_interfaces { - // Return the interface ID DeclareInterfaceID( IMCFTAttenuationTool, 1, 0 ); @@ -32,4 +27,3 @@ struct IMCFTAttenuationTool : extend_interfaces { }; virtual Attenuation attenuation( double hitXPosition, double hitYPosition ) const = 0; }; -#endif // IMCFTATTENUATIONTOOL_H diff --git a/FT/FTDigitisation/src/MCFTAttenuationTool.cpp b/FT/FTDigitisation/src/MCFTAttenuationTool.cpp index d0c943c1b..2b5fdb8b7 100644 --- a/FT/FTDigitisation/src/MCFTAttenuationTool.cpp +++ b/FT/FTDigitisation/src/MCFTAttenuationTool.cpp @@ -8,120 +8,191 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ -// local -#include "MCFTAttenuationTool.h" -//----------------------------------------------------------------------------- -// Implementation file for class : MCFTAttenuationTool -// -// 2015-01-15 : Michel De Cian -//----------------------------------------------------------------------------- +#include "IMCFTAttenuationTool.h" + +#include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/SystemOfUnits.h" + +/** + * Tool to calculate the attenuation length for direct and reflectec light in a scintillating fibre. + * + * - ShortAttenuationLength: Distance along the fibre to divide the light amplitude by a factor e : short component + * - LongAttenuationLength: Distance along the fibre to divide the light amplitude by a factor e : long component + * - FractionShort: Fraction of short attenuation at SiPM + * - XMaxIrradiatedZone: Size in x of the zone where fibres are irradiated + * - YMaxIrradiatedZone: Size in y of the zone where fibres are irradiated + * - IrradiatedAttenuationLength: Attenuation length by steps in the zone + * - ReflectionCoefficient: Reflection coefficient of the fibre mirrored side, from 0 to 1 + * - BeginReflectionLossY: Beginning of zone where reflection is too late and lost + * - EndReflectionLossY: End of zone where reflection is too late and lost + * - XStepOfMap: Step along x-axis of the FullAttenuationMap (in mm) + * - YStepOfMap: Step along y-axis of the FullAttenuationMap (in mm) + * + * @author Michel De Cian, based on Eric Cogneras' implementation + * @date 2015-01-15 + */ + +class MCFTAttenuationTool : public extends { + +public: + using base_class::base_class; + StatusCode initialize() override; + + /// Calculate the direct attenuation and the attenuation with reflection + IMCFTAttenuationTool::Attenuation attenuation( double x, double y ) const override; + +private: + inline double calcAtt( const double fracX, const double fracY, const int kx, const int ky ) const { + return fracX * ( fracY * m_transmissionMap[m_nYSteps * ( kx + 1 ) + ky + 1] + + ( 1 - fracY ) * m_transmissionMap[m_nYSteps * ( kx + 1 ) + ky] ) + + ( 1 - fracX ) * ( fracY * m_transmissionMap[m_nYSteps * kx + ky + 1] + + ( 1 - fracY ) * m_transmissionMap[m_nYSteps * kx + ky] ); + } + + inline double calcAttRef( const double fracX, const double fracY, const int kx, const int ky ) const { + return fracX * ( fracY * m_transmissionRefMap[m_nYSteps * ( kx + 1 ) + ky + 1] + + ( 1 - fracY ) * m_transmissionRefMap[m_nYSteps * ( kx + 1 ) + ky] ) + + ( 1 - fracX ) * ( fracY * m_transmissionRefMap[m_nYSteps * kx + ky + 1] + + ( 1 - fracY ) * m_transmissionRefMap[m_nYSteps * kx + ky] ); + } + + Gaudi::Property m_shortAttenuationLength{this, "ShortAttenuationLength", 200 * Gaudi::Units::mm, + "Attenuation length of the light along the fibre: short component"}; + Gaudi::Property m_longAttenuationLength{this, "LongAttenuationLength", 4700 * Gaudi::Units::mm, + "Attenuation length of the light along the fibre: long component"}; + Gaudi::Property m_fractionShort{this, "FractionShort", 0.18, "Fraction of short attenuation at SiPM"}; + Gaudi::Property m_xMaxIrradiatedZone{this, "XMaxIrradiatedZone", 2000. * Gaudi::Units::mm, + "Size in x of the zone where fibres are irradiated"}; + Gaudi::Property m_yMaxIrradiatedZone{this, "YMaxIrradiatedZone", 500. * Gaudi::Units::mm, + "Size in y of the zone where fibres are irradiated"}; + Gaudi::Property> m_irradiatedAttenuationLength{ + this, + "IrradiatedAttenuationLength", + {3000. * Gaudi::Units::mm, 3000. * Gaudi::Units::mm, 3000. * Gaudi::Units::mm, 1500. * Gaudi::Units::mm, + 300. * Gaudi::Units::mm}, + "Attenuation length by steps in the zone"}; + Gaudi::Property m_reflectionCoefficient{this, "ReflectionCoefficient", 0.7, + "Reflection coefficient of the fibre mirrored side, from 0 to 1"}; + Gaudi::Property m_beginReflectionLossY{this, "BeginReflectionLossY", 1000. * Gaudi::Units::mm, + "begin of zone where reflection is too late and lost"}; + Gaudi::Property m_endReflectionLossY{this, "EndReflectionLossY", 1500. * Gaudi::Units::mm, + "end of zone where reflection is too late and lost"}; + + Gaudi::Property m_xStepOfMap{this, "XStepOfMap", 200. * Gaudi::Units::mm, + "Step along X-axis of the FullAttenuationMap(in mm)"}; + Gaudi::Property m_yStepOfMap{this, "YStepOfMap", 100. * Gaudi::Units::mm, + "Step along Y-axis of the FullAttenuationMap(in mm)"}; + + Gaudi::Property m_xMax{this, "XMax", 3164.4 * Gaudi::Units::mm, + "Length of the FullAttenuationMap(in mm) in X-axis"}; + + Gaudi::Property m_yMax{this, "YMax", 2426.0 * Gaudi::Units::mm, + "Length of the FullAttenuationMap(in mm) in Y-axis"}; + + int m_nXSteps; + int m_nYSteps; + + std::vector m_transmissionMap; ///< Maps hits to transmitted energy from the direct pulse + std::vector m_transmissionRefMap; ///< Maps hits to transmitted energy from the reflected pulse +}; // Declaration of the Tool Factory DECLARE_COMPONENT( MCFTAttenuationTool ) -//============================================================================= -// Calculate the attenuation -//============================================================================= IMCFTAttenuationTool::Attenuation MCFTAttenuationTool::attenuation( double hitXPosition, double hitYPosition ) const { - //== Compute attenuation by interpolation in the table //== No difference made here between x-layers and (u,v)-layers : all layers assumed to be x-layers const int kx = std::min( m_nXSteps - 2, (int)( fabs( hitXPosition ) / m_xStepOfMap ) ); const int ky = std::min( m_nYSteps - 2, (int)( fabs( hitYPosition ) / m_yStepOfMap ) ); const double fracX = fabs( hitXPosition ) / m_xStepOfMap - kx; const double fracY = fabs( hitYPosition ) / m_yStepOfMap - ky; - return {calcAtt( fracX, fracY, kx, ky ), calcAttRef( fracX, fracY, kx, ky )}; } -//========================================================================= -// Calculate the transmission map -//========================================================================= StatusCode MCFTAttenuationTool::initialize() { - auto sc = base_class::initialize(); - if ( sc.isFailure() ) return sc; - - //== Generate the transmission map. - //== Two attenuation length, then in the radiation area a different attenuation length - //== Compute also the transmission of the reflected signal, attenuated by the two length - // and the irradiated area - - // Setup maps size - m_nXSteps = m_xMax / m_xStepOfMap + 2; // add x=0 and x > max position - m_nYSteps = m_yMax / m_yStepOfMap + 2; // same for y - - // const double xMax = (m_nXSteps - 1) * m_xStepOfMap; - const double yMax = ( m_nYSteps - 1 ) * m_yStepOfMap; - - // m_transmissionMap set at 1E-10 initialisatino value to avoid 'division by zero' bug - m_transmissionMap.resize( m_nXSteps * m_nYSteps, 1E-10 ); - m_transmissionRefMap.resize( m_nXSteps * m_nYSteps, 1E-10 ); - - // Loop on the x axis - for ( int kx = 0; m_nXSteps > kx; ++kx ) { - const double x = kx * m_xStepOfMap; - - // -- is this a model, a fact, an urban legend? - double radZoneSize = 2 * m_yMaxIrradiatedZone * ( 1 - x / m_xMaxIrradiatedZone ); - if ( 0. > radZoneSize ) radZoneSize = 0.; - const double yBoundaryRadZone = .5 * radZoneSize; - - // Loop on the y axis and calculate transmission map - for ( int ky = 0; m_nYSteps > ky; ++ky ) { - const double y = yMax - ky * m_yStepOfMap; - - // Compute attenuation according to the crossed fibre length - double att = ( m_fractionShort * exp( -( yMax - y ) / m_shortAttenuationLength ) + - ( 1 - m_fractionShort ) * exp( -( yMax - y ) / m_longAttenuationLength ) ); - - // plot2D(x,y,"FibreAttenuationMap","Attenuation coefficient as a function of the position (Quarter 3); x [mm];y - // [mm]", - // 0., m_nXSteps*m_xStepOfMap, 0.,m_nYSteps*m_yStepOfMap, m_nXSteps, m_nYSteps, att); - - if ( y < yBoundaryRadZone ) { - att = ( m_fractionShort * exp( -( yMax - yBoundaryRadZone ) / m_shortAttenuationLength ) + - ( 1 - m_fractionShort ) * exp( -( yMax - yBoundaryRadZone ) / m_longAttenuationLength ) ); - - double lInRadiation = yBoundaryRadZone - y; - - for ( unsigned int kz = 0; m_irradiatedAttenuationLength.size() > kz; ++kz ) { - if ( lInRadiation > m_yStepOfMap ) { - att *= exp( -m_yStepOfMap / m_irradiatedAttenuationLength[kz] ); - } else if ( lInRadiation > 0. ) { - att *= exp( -lInRadiation / m_irradiatedAttenuationLength[kz] ); + return extends::initialize().andThen( [&] { + //== Generate the transmission map. + //== Two attenuation length, then in the radiation area a different attenuation length + //== Compute also the transmission of the reflected signal, attenuated by the two length + // and the irradiated area + + // Setup maps size + m_nXSteps = m_xMax / m_xStepOfMap + 2; // add x=0 and x > max position + m_nYSteps = m_yMax / m_yStepOfMap + 2; // same for y + + // const double xMax = (m_nXSteps - 1) * m_xStepOfMap; + const double yMax = ( m_nYSteps - 1 ) * m_yStepOfMap; + + // m_transmissionMap set at 1E-10 initialisatino value to avoid 'division by zero' bug + m_transmissionMap.resize( m_nXSteps * m_nYSteps, 1E-10 ); + m_transmissionRefMap.resize( m_nXSteps * m_nYSteps, 1E-10 ); + + // Loop on the x axis + for ( int kx = 0; m_nXSteps > kx; ++kx ) { + const double x = kx * m_xStepOfMap; + + // -- is this a model, a fact, an urban legend? + double radZoneSize = 2 * m_yMaxIrradiatedZone * ( 1 - x / m_xMaxIrradiatedZone ); + if ( 0. > radZoneSize ) radZoneSize = 0.; + const double yBoundaryRadZone = .5 * radZoneSize; + + // Loop on the y axis and calculate transmission map + for ( int ky = 0; m_nYSteps > ky; ++ky ) { + const double y = yMax - ky * m_yStepOfMap; + + // Compute attenuation according to the crossed fibre length + double att = ( m_fractionShort * exp( -( yMax - y ) / m_shortAttenuationLength ) + + ( 1 - m_fractionShort ) * exp( -( yMax - y ) / m_longAttenuationLength ) ); + + // plot2D(x,y,"FibreAttenuationMap","Attenuation coefficient as a function of the position (Quarter 3); x [mm];y + // [mm]", + // 0., m_nXSteps*m_xStepOfMap, 0.,m_nYSteps*m_yStepOfMap, m_nXSteps, m_nYSteps, att); + + if ( y < yBoundaryRadZone ) { + att = ( m_fractionShort * exp( -( yMax - yBoundaryRadZone ) / m_shortAttenuationLength ) + + ( 1 - m_fractionShort ) * exp( -( yMax - yBoundaryRadZone ) / m_longAttenuationLength ) ); + + double lInRadiation = yBoundaryRadZone - y; + + for ( unsigned int kz = 0; m_irradiatedAttenuationLength.size() > kz; ++kz ) { + if ( lInRadiation > m_yStepOfMap ) { + att *= exp( -m_yStepOfMap / m_irradiatedAttenuationLength[kz] ); + } else if ( lInRadiation > 0. ) { + att *= exp( -lInRadiation / m_irradiatedAttenuationLength[kz] ); + } + lInRadiation -= m_yStepOfMap; } - lInRadiation -= m_yStepOfMap; } + m_transmissionMap[m_nYSteps * ( kx + 1 ) - ky - 1] = att; } - m_transmissionMap[m_nYSteps * ( kx + 1 ) - ky - 1] = att; - } - // Loop on the y axis and calculate transmission of the reflected pulse - // The energy from the reflected pulse is calculated as the energy coming from y = 0 and attenuated - // by the reflection coefficient, further attenuated depending on the y. - if ( m_reflectionCoefficient > 0. ) { + // Loop on the y axis and calculate transmission of the reflected pulse + // The energy from the reflected pulse is calculated as the energy coming from y = 0 and attenuated + // by the reflection coefficient, further attenuated depending on the y. + if ( m_reflectionCoefficient > 0. ) { - // Reflected energy from y = 0 - double reflected = m_transmissionMap[m_nYSteps * kx] * m_reflectionCoefficient; + // Reflected energy from y = 0 + double reflected = m_transmissionMap[m_nYSteps * kx] * m_reflectionCoefficient; - for ( int kk = 0; m_nYSteps > kk; ++kk ) { - // const double y = kk * m_yStepOfMap; + for ( int kk = 0; m_nYSteps > kk; ++kk ) { + // const double y = kk * m_yStepOfMap; - // Evaluate the attenuation factor. - // Rather than calculating it use information from the transmission map in the same point - double att = 0; + // Evaluate the attenuation factor. + // Rather than calculating it use information from the transmission map in the same point + double att = 0; - if ( m_nYSteps > ( kk + 1 ) ) - att = m_transmissionMap[m_nYSteps * kx + kk] / m_transmissionMap[m_nYSteps * kx + kk + 1]; + if ( m_nYSteps > ( kk + 1 ) ) + att = m_transmissionMap[m_nYSteps * kx + kk] / m_transmissionMap[m_nYSteps * kx + kk + 1]; - // Fill the map - m_transmissionRefMap[m_nYSteps * kx + kk] = reflected; + // Fill the map + m_transmissionRefMap[m_nYSteps * kx + kk] = reflected; - // Prepare the reflection ratio for next iteration - reflected *= att; + // Prepare the reflection ratio for next iteration + reflected *= att; + } } } - } - return sc; + return StatusCode::SUCCESS; + } ); } diff --git a/FT/FTDigitisation/src/MCFTAttenuationTool.h b/FT/FTDigitisation/src/MCFTAttenuationTool.h deleted file mode 100644 index 567d5690e..000000000 --- a/FT/FTDigitisation/src/MCFTAttenuationTool.h +++ /dev/null @@ -1,107 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration * -* * -* This software is distributed under the terms of the GNU General Public * -* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". * -* * -* In applying this licence, CERN does not waive the privileges and immunities * -* granted to it by virtue of its status as an Intergovernmental Organization * -* or submit itself to any jurisdiction. * -\*****************************************************************************/ -#ifndef MCFTATTENUATIONTOOL_H -#define MCFTATTENUATIONTOOL_H 1 - -// from Gaudi -#include "GaudiAlg/GaudiTool.h" -#include "GaudiKernel/SystemOfUnits.h" - -#include "IMCFTAttenuationTool.h" // Interface - -/** @class MCFTAttenuationTool MCFTAttenuationTool.h - * - * Tool to calculate the attenuation length for direct and reflectec light in a scintillating fibre. - * - * - ShortAttenuationLength: Distance along the fibre to divide the light amplitude by a factor e : short component - * - LongAttenuationLength: Distance along the fibre to divide the light amplitude by a factor e : long component - * - FractionShort: Fraction of short attenuation at SiPM - * - XMaxIrradiatedZone: Size in x of the zone where fibres are irradiated - * - YMaxIrradiatedZone: Size in y of the zone where fibres are irradiated - * - IrradiatedAttenuationLength: Attenuation length by steps in the zone - * - ReflectionCoefficient: Reflection coefficient of the fibre mirrored side, from 0 to 1 - * - BeginReflectionLossY: Beginning of zone where reflection is too late and lost - * - EndReflectionLossY: End of zone where reflection is too late and lost - * - XStepOfMap: Step along x-axis of the FullAttenuationMap (in mm) - * - YStepOfMap: Step along y-axis of the FullAttenuationMap (in mm) - * - * - * - * @author Michel De Cian, based on Eric Cogneras' implementation - * @date 2015-01-15 - */ - -class MCFTAttenuationTool : public extends { - -public: - using base_class::base_class; - - /// Initialize the transmission map - StatusCode initialize() override; - - /// Calculate the direct attenuation and the attenuation with reflection - IMCFTAttenuationTool::Attenuation attenuation( double x, double y ) const override; - -private: - inline double calcAtt( const double fracX, const double fracY, const int kx, const int ky ) const { - return fracX * ( fracY * m_transmissionMap[m_nYSteps * ( kx + 1 ) + ky + 1] + - ( 1 - fracY ) * m_transmissionMap[m_nYSteps * ( kx + 1 ) + ky] ) + - ( 1 - fracX ) * ( fracY * m_transmissionMap[m_nYSteps * kx + ky + 1] + - ( 1 - fracY ) * m_transmissionMap[m_nYSteps * kx + ky] ); - } - - inline double calcAttRef( const double fracX, const double fracY, const int kx, const int ky ) const { - return fracX * ( fracY * m_transmissionRefMap[m_nYSteps * ( kx + 1 ) + ky + 1] + - ( 1 - fracY ) * m_transmissionRefMap[m_nYSteps * ( kx + 1 ) + ky] ) + - ( 1 - fracX ) * ( fracY * m_transmissionRefMap[m_nYSteps * kx + ky + 1] + - ( 1 - fracY ) * m_transmissionRefMap[m_nYSteps * kx + ky] ); - } - - Gaudi::Property m_shortAttenuationLength{this, "ShortAttenuationLength", 200 * Gaudi::Units::mm, - "Attenuation length of the light along the fibre: short component"}; - Gaudi::Property m_longAttenuationLength{this, "LongAttenuationLength", 4700 * Gaudi::Units::mm, - "Attenuation length of the light along the fibre: long component"}; - Gaudi::Property m_fractionShort{this, "FractionShort", 0.18, "Fraction of short attenuation at SiPM"}; - Gaudi::Property m_xMaxIrradiatedZone{this, "XMaxIrradiatedZone", 2000. * Gaudi::Units::mm, - "Size in x of the zone where fibres are irradiated"}; - Gaudi::Property m_yMaxIrradiatedZone{this, "YMaxIrradiatedZone", 500. * Gaudi::Units::mm, - "Size in y of the zone where fibres are irradiated"}; - Gaudi::Property> m_irradiatedAttenuationLength{ - this, - "IrradiatedAttenuationLength", - {3000. * Gaudi::Units::mm, 3000. * Gaudi::Units::mm, 3000. * Gaudi::Units::mm, 1500. * Gaudi::Units::mm, - 300. * Gaudi::Units::mm}, - "Attenuation length by steps in the zone"}; - Gaudi::Property m_reflectionCoefficient{this, "ReflectionCoefficient", 0.7, - "Reflection coefficient of the fibre mirrored side, from 0 to 1"}; - Gaudi::Property m_beginReflectionLossY{this, "BeginReflectionLossY", 1000. * Gaudi::Units::mm, - "begin of zone where reflection is too late and lost"}; - Gaudi::Property m_endReflectionLossY{this, "EndReflectionLossY", 1500. * Gaudi::Units::mm, - "end of zone where reflection is too late and lost"}; - - Gaudi::Property m_xStepOfMap{this, "XStepOfMap", 200. * Gaudi::Units::mm, - "Step along X-axis of the FullAttenuationMap(in mm)"}; - Gaudi::Property m_yStepOfMap{this, "YStepOfMap", 100. * Gaudi::Units::mm, - "Step along Y-axis of the FullAttenuationMap(in mm)"}; - - Gaudi::Property m_xMax{this, "XMax", 3164.4 * Gaudi::Units::mm, - "Length of the FullAttenuationMap(in mm) in X-axis"}; - - Gaudi::Property m_yMax{this, "YMax", 2426.0 * Gaudi::Units::mm, - "Length of the FullAttenuationMap(in mm) in Y-axis"}; - - int m_nXSteps; - int m_nYSteps; - - std::vector m_transmissionMap; ///< Maps hits to transmitted energy from the direct pulse - std::vector m_transmissionRefMap; ///< Maps hits to transmitted energy from the reflected pulse -}; -#endif // MCFTATTENUATIONTOOL_H -- GitLab From 141536a0f8fe9ecdca9f0f7718276d42b5dcc7a2 Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Tue, 4 Oct 2022 10:20:18 +0200 Subject: [PATCH 05/15] Modernized SiPMResponse and switched to derived conditions --- .../options/Digitisation_BasicOptions.py | 2 - FT/FTDigitisation/src/SiPMResponse.cpp | 164 +++++++----------- FT/FTDigitisation/src/SiPMResponse.h | 41 ++--- 3 files changed, 75 insertions(+), 132 deletions(-) diff --git a/FT/FTDigitisation/options/Digitisation_BasicOptions.py b/FT/FTDigitisation/options/Digitisation_BasicOptions.py index 8b7dceec7..943700d3a 100644 --- a/FT/FTDigitisation/options/Digitisation_BasicOptions.py +++ b/FT/FTDigitisation/options/Digitisation_BasicOptions.py @@ -105,8 +105,6 @@ def doIt(): digitCreator.addTool(SiPMResponse, name="SiPMResponse") digitCreator.SiPMResponse.ElectronicsResponse = "pacific5q_pz5" - #digitCreator.SiPMResponse.Tshift = 5 - #digitCreator.SiPMGainShift = 1.0 depositCreator = MCFTDepositCreator("MCFTDepositCreator") depositCreator.addTool(FTSiPMTool("FTSiPMTool")) diff --git a/FT/FTDigitisation/src/SiPMResponse.cpp b/FT/FTDigitisation/src/SiPMResponse.cpp index 15d68c72c..7f625159d 100644 --- a/FT/FTDigitisation/src/SiPMResponse.cpp +++ b/FT/FTDigitisation/src/SiPMResponse.cpp @@ -8,12 +8,11 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ -// from Gaudi +#include "SiPMResponse.h" + #include "DetDesc/Condition.h" -#include "GaudiKernel/SystemOfUnits.h" -// local -#include "SiPMResponse.h" +#include //----------------------------------------------------------------------------- // Implementation file for class : SiPMResponse @@ -26,113 +25,68 @@ // Declaration of the Tool Factory DECLARE_COMPONENT( SiPMResponse ) -//============================================================================= -// Standard constructor, initializes variables -//============================================================================= +namespace { + GaudiMath::Interpolation::Type typeFromString( std::string const& splineType ) { + GaudiMath::Interpolation::Type aType; + using namespace GaudiMath::Interpolation; + if ( splineType == "Cspline" ) + aType = Cspline; + else if ( splineType == "Linear" ) + aType = Linear; + else if ( splineType == "Polynomial" ) + aType = Polynomial; + else if ( splineType == "Akima" ) + aType = Akima; + else if ( splineType == "Akima_Periodic" ) + aType = Akima_Periodic; + else if ( splineType == "Cspline_Periodic" ) + aType = Cspline_Periodic; + else + aType = Cspline; // default to cubic + return aType; + } +} // namespace + SiPMResponse::SiPMResponse( const std::string& type, const std::string& name, const IInterface* parent ) - : GaudiTool( type, name, parent ), m_responseSpline( 0 ) { + : ConditionAccessorHolder( type, name, parent ) { declareInterface( this ); } -//============================================================================= -// Destructor -//============================================================================= -SiPMResponse::~SiPMResponse() { - if ( 0 != m_responseSpline ) delete m_responseSpline; -} -//============================================================================= -// Initialize the pulse shape -//============================================================================= StatusCode SiPMResponse::initialize() { - StatusCode sc = GaudiTool::initialize(); - if ( sc.isFailure() ) return Error( "Failed to initialize", sc ); - - // initialise SiPMResponse as zero in case conditions DB fails - std::vector times; - std::vector values; - if ( !existDet( detSvc(), m_conditionLocation.value() ) ) { - return Error( "SiPMResponse found not in conditions DB", StatusCode::FAILURE ); - } else { - Condition* rInfo = getDet( m_conditionLocation ); - times = rInfo->param>( "times" ); - if ( m_electronicsResponse == "pacific5q_pz5" ) { - values = rInfo->param>( "values_P5_pz5" ); - m_sipmGainShift = rInfo->param( "gainshift_P5_pz5" ); - m_tshift = rInfo->param( "timeshift_P5_pz5" ); - } - - if ( m_electronicsResponse == "pacific5q_pz6" ) { - values = rInfo->param>( "values_P5_pz6" ); - m_sipmGainShift = rInfo->param( "gainshift_P5_pz6" ); - m_tshift = rInfo->param( "timeshift_P5_pz6" ); - } - } - - // shift times by tshift - std::for_each( times.begin(), times.end(), [this]( double& n ) { n += m_tshift; } ); - - // Check if the times and values are of equal length - if ( times.size() != values.size() ) { return Error( "inconsistent data !", StatusCode::FAILURE ); } - - // Store the first and last entry of the vector of times - m_tMin = times.front(); - m_tMax = times.back(); - - // Normalize function. Set maximum to 1. - auto maxVal = *max_element( std::begin( values ), std::end( values ) ); - std::transform( values.begin(), values.end(), values.begin(), - std::bind1st( std::multiplies(), 1. / maxVal ) ); - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "SiPMResponse max = " << maxVal << endmsg; - debug() << "SiPMResponse after reweighting:" << endmsg; - debug() << "----------------------------" << endmsg; - for ( unsigned int i = 0; i < values.size(); ++i ) - debug() << "\t t = " << times[i] << "\t val = " << values[i] << endmsg; - } - - // Fit the spline to the data - m_responseSpline = new GaudiMath::SimpleSpline( times, values, typeFromString() ); - - return sc; + return ConditionAccessorHolder::initialize().andThen( [&] { + addConditionDerivation( {m_conditionLocation}, m_conditions.key(), [&]( YAML::Node const& rInfo ) -> Conditions { + // Get data from condition DB + auto times = rInfo["times"].as>(); + std::vector values; + double tshift{0}; // Time shift (t0) of the response shape [ns] + double sipmGainShift{1}; + if ( m_electronicsResponse == "pacific5q_pz5" ) { + values = rInfo["values_P5_pz5"].as>(); + sipmGainShift = rInfo["gainshift_P5_pz5"].as(); + tshift = rInfo["timeshift_P5_pz5"].as(); + } else if ( m_electronicsResponse == "pacific5q_pz6" ) { + values = rInfo["values_P5_pz6"].as>(); + sipmGainShift = rInfo["gainshift_P5_pz6"].as(); + tshift = rInfo["timeshift_P5_pz6"].as(); + } + // shift times by tshift + std::for_each( times.begin(), times.end(), [tshift]( double& n ) { n += tshift; } ); + // Check if the times and values are of equal length + if ( times.size() != values.size() ) { + throw GaudiException( "SiPMResponse::initialize", "inconsistent data !", StatusCode::FAILURE ); + } + // Normalize function. Set maximum to 1. + double maxVal = *std::max_element( std::begin( values ), std::end( values ) ); + std::transform( values.begin(), values.end(), values.begin(), [maxVal]( double item ) { return item / maxVal; } ); + // Fit the spline to the data + return {{times, values, typeFromString( m_splineType )}, times.front(), times.back(), sipmGainShift}; + } ); + return StatusCode::SUCCESS; + } ); } double SiPMResponse::response( const double time ) const { - return ( ( time > m_tMin ) && ( time < m_tMax ) ? m_responseSpline->eval( time ) * m_sipmGainShift : 0.0 ); -} - -std::string SiPMResponse::responsefunction() const { return ( m_electronicsResponse ); } - -double SiPMResponse::gainshift() const { return ( m_sipmGainShift ); } - -GaudiMath::Interpolation::Type SiPMResponse::typeFromString() const { - GaudiMath::Interpolation::Type aType; - using namespace GaudiMath::Interpolation; - if ( m_splineType == "Cspline" ) - aType = Cspline; - else if ( m_splineType == "Linear" ) - aType = Linear; - else if ( m_splineType == "Polynomial" ) - aType = Polynomial; - else if ( m_splineType == "Akima" ) - aType = Akima; - else if ( m_splineType == "Akima_Periodic" ) - aType = Akima_Periodic; - else if ( m_splineType == "Cspline_Periodic" ) - aType = Cspline_Periodic; - else - aType = Cspline; // default to cubic - - return aType; + auto& conds = m_conditions.get(); + return ( ( time > conds.tMin ) && ( time < conds.tMax ) ? conds.responseSpline.eval( time ) * conds.sipmGainShift + : 0.0 ); } - -void SiPMResponse::sample( std::vector& times, std::vector& values ) const { - double t = m_tMin; - while ( t < m_tMax ) { - times.push_back( t ); - values.push_back( response( t ) ); - t += m_samplingDt; - } // loop times -} - -//============================================================================= diff --git a/FT/FTDigitisation/src/SiPMResponse.h b/FT/FTDigitisation/src/SiPMResponse.h index fc0bec78f..f90592d1f 100644 --- a/FT/FTDigitisation/src/SiPMResponse.h +++ b/FT/FTDigitisation/src/SiPMResponse.h @@ -8,55 +8,46 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ -#ifndef SIPMRESPONSE_H -#define SIPMRESPONSE_H 1 +#pragma once + +#include "DetDesc/GenericConditionAccessorHolder.h" -// from Gaudi #include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/SystemOfUnits.h" #include "GaudiMath/GaudiMath.h" static const InterfaceID IID_SiPMResponse( "SiPMResponse", 1, 0 ); -/** @class SiPMResponse SiPMResponse.h Boole_v26r9/SiPMResponse.h - * +/** * This class describes the SiPM response to a single pulse * * @author Maurizio Martinelli * @date 2013-11-12 */ -class SiPMResponse : public GaudiTool { +class SiPMResponse : public LHCb::DetDesc::ConditionAccessorHolder { public: - // Return the interface ID static const InterfaceID& interfaceID() { return IID_SiPMResponse; } - - /// Standard constructor SiPMResponse( const std::string& type, const std::string& name, const IInterface* parent ); - - virtual ~SiPMResponse(); ///< Destructor - StatusCode initialize() override; virtual double response( const double time ) const; - virtual std::string responsefunction() const; - virtual double gainshift() const; - -protected: - GaudiMath::Interpolation::Type typeFromString() const; - void sample( std::vector& times, std::vector& val ) const; + virtual std::string responsefunction() const { return ( m_electronicsResponse ); }; + virtual double gainshift() const { return ( m_conditions.get().sipmGainShift ); }; private: - GaudiMath::SimpleSpline* m_responseSpline; - double m_tMin; - double m_tMax; + struct Conditions { + GaudiMath::SimpleSpline responseSpline; + double tMin; + double tMax; + /// calibration of the gain + double sipmGainShift; + }; + ConditionAccessor m_conditions{this, name() + "_Conditions"}; - // properties Gaudi::Property m_samplingDt{this, "samplingDt", 0.1 * Gaudi::Units::ns, "Sampling time step"}; Gaudi::Property m_splineType{this, "splineType", "Cspline", "The spline type"}; Gaudi::Property m_electronicsResponse{ this, "ElectronicsResponse", "pacific5q_pz5", "Select which electronics response (pacific) to use {pacific4,pacific5q_pz5,pacific5q_pz5}"}; - Gaudi::Property m_tshift{this, "Tshift", 0, "Time shift (t0) of the response shape [ns]"}; - Gaudi::Property m_sipmGainShift{this, "SiPMGainShift", 1., "calibration of the gain"}; Gaudi::Property m_conditionLocation{this, "conditionLocation", "/dd/Conditions/Sim/FT/SiPMResponse"}; }; -#endif // SIPMRESPONSE_H -- GitLab From c7b3983bee024caabd7d9db59ee68692864ec6ef Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Tue, 4 Oct 2022 10:20:37 +0200 Subject: [PATCH 06/15] Modernized MCFTG4AttenuationTool and switched to derived conditions --- .../src/MCFTG4AttenuationTool.cpp | 170 ++++++++++-------- FT/FTDigitisation/src/MCFTG4AttenuationTool.h | 58 ------ 2 files changed, 95 insertions(+), 133 deletions(-) delete mode 100644 FT/FTDigitisation/src/MCFTG4AttenuationTool.h diff --git a/FT/FTDigitisation/src/MCFTG4AttenuationTool.cpp b/FT/FTDigitisation/src/MCFTG4AttenuationTool.cpp index b305b8e90..0f3847f4d 100644 --- a/FT/FTDigitisation/src/MCFTG4AttenuationTool.cpp +++ b/FT/FTDigitisation/src/MCFTG4AttenuationTool.cpp @@ -8,98 +8,118 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ -#include "MCFTG4AttenuationTool.h" -#include -//----------------------------------------------------------------------------- -// Implementation file for class : MCFTG4ttenuationTool -// -// 2017-07-11 : Martin Bieker -//----------------------------------------------------------------------------- +#include "IMCFTAttenuationTool.h" + +#include "Kernel/STLExtensions.h" +#include + +#include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/SystemOfUnits.h" + +#include + +namespace { + struct Conditions { + unsigned int nBinsX, nBinsY; + std::vector xEdges, yEdges, effDir, effRef; + }; + + int findBin( LHCb::span axis, double position ) { + auto iterator = std::upper_bound( axis.begin(), axis.end(), position ); + // If position is outside of map: use first or last bin. + if ( iterator == axis.end() ) + --iterator; + else if ( iterator == axis.begin() ) + ++iterator; + return iterator - axis.begin() - 1; + } +} // namespace + +/** + * Tool that reads attenutation maps from ConDB + * + * @author Martin Bieker based on implementation of M. Demmer and J. Wishahi + * @date 2015-01-15 + * @date 2017-07-11 + */ +class MCFTG4AttenuationTool : public extends, IMCFTAttenuationTool> { +public: + using extends::extends; + + /// Initialize the transmission map + StatusCode initialize() override; + + /// Calculate the direct attenuation and the attenuation with reflection + IMCFTAttenuationTool::Attenuation attenuation( double x, double y ) const override; + +private: + Gaudi::Property m_mirrorReflectivity{this, "MirrorReflectivity", 0.75, // from LHCb-PUB-2014-020 + "Reflectivity of mirror at the end of the fibre mat (0-1)"}; + Gaudi::Property m_irradiation{this, "Irradiation", 50, "Simulated radiation damage"}; + Gaudi::Property m_irradiationLinearModel{this, "IrradiationLinearModel", false, + "Radiation damage model (default is power-law model)"}; + Gaudi::Property m_agingFibres{this, "AgingFibres", false, "Simulate time aging of the fibres"}; + Gaudi::Property m_replaceModules{this, "ReplaceModules", false, "Simulate replacement of inner modules"}; + + ConditionAccessor m_conditions{this, name() + "_Conditions"}; +}; // Declaration of the Tool Factory DECLARE_COMPONENT( MCFTG4AttenuationTool ) -//============================================================================= // Calculate the attenuation -//============================================================================= IMCFTAttenuationTool::Attenuation MCFTG4AttenuationTool::attenuation( double hitXPosition, double hitYPosition ) const { hitXPosition = std::abs( hitXPosition ); hitYPosition = std::abs( hitYPosition ); - int binX = findBin( m_xEdges, hitXPosition ); - int binY = findBin( m_yEdges, hitYPosition ); + auto& conds = m_conditions.get(); + int binX = findBin( conds.xEdges, hitXPosition ); + int binY = findBin( conds.yEdges, hitYPosition ); // Hit outside quadrant (should never happen!) if ( binX == -1 || binY == -1 ) { return {0.0, 0.0}; } else { - return {m_effDir[binY * m_nBinsX + binX], m_effRef[binY * m_nBinsX + binX]}; + return {conds.effDir[binY * conds.nBinsX + binX], conds.effRef[binY * conds.nBinsX + binX]}; } } -//========================================================================= -// Retrieve the Attenuation map from the CondDB -//========================================================================= StatusCode MCFTG4AttenuationTool::initialize() { - auto sc = base_class::initialize(); - if ( sc.isFailure() ) { return sc; } - - // Select the type of radiation map - std::string conditionsLocation; - conditionsLocation = - "/dd/Conditions/Calibration/FT/AttenuationMap" + std::to_string( m_irradiation ) + "fb_z_920-940"; - if ( m_irradiationLinearModel ) conditionsLocation += "_linearmodel"; - if ( m_agingFibres ) conditionsLocation += "_aging"; - if ( m_replaceModules ) conditionsLocation += "_replacedmodule"; - - info() << "Selecting " << std::to_string( m_irradiation ) << "/fb attenuation map using " - << std::string( m_irradiationLinearModel ? "linear" : "power-law" ) << " radiation damage model " - << std::string( m_agingFibres ? "with" : "without" ) << " aging fibres and " - << std::string( m_replaceModules ? "with" : "without" ) << " replaced center module. " << endmsg; - - // When conditions are not found, return failure - if ( existDet( conditionsLocation ) == false ) { - return Error( "No attenuation maps found in SIMCOND: update to a newer tag", StatusCode::FAILURE ); - } - - // Get the conditions - Condition* cond = getDet( conditionsLocation ); - - m_nBinsX = cond->param( "x_n_bins" ); - m_nBinsY = cond->param( "y_n_bins" ); - - m_xEdges = cond->param>( "x_edges" ); - m_yEdges = cond->param>( "y_edges" ); - std::sort( m_xEdges.begin(), m_xEdges.end() ); - std::sort( m_yEdges.begin(), m_yEdges.end() ); - - m_effDir = cond->param>( "eff_dir" ); - m_effRef = cond->param>( "eff_refl" ); - - // Validate Attenuatom Map - if ( !validateMap() ) return Error( "Attenuation Map does not match bin counts" ); - - std::transform( m_effRef.begin(), m_effRef.end(), m_effRef.begin(), - [this]( double val ) { return val * m_mirrorReflectivity; } ); - - if ( msgLevel( MSG::DEBUG ) ) { debug() << "nBinsX :" << m_nBinsX << " nBinsY: " << m_nBinsY << std::endl; } - - return sc; -} - -bool MCFTG4AttenuationTool::validateMap() { - if ( m_nBinsX != m_xEdges.size() - 1 || m_nBinsY != m_yEdges.size() - 1 ) { return false; } - if ( m_nBinsX * m_nBinsY != m_effDir.size() || m_nBinsX * m_nBinsY != m_effRef.size() ) { return false; } - return true; -} - -int MCFTG4AttenuationTool::findBin( LHCb::span axis, double position ) const { - auto iterator = std::upper_bound( axis.begin(), axis.end(), position ); - // If position is outside of map: use first or last bin. - if ( iterator == axis.end() ) - --iterator; - else if ( iterator == axis.begin() ) - ++iterator; - return iterator - axis.begin() - 1; + return extends::initialize().andThen( [&] { + std::string conditionsLocation = + "/dd/Conditions/Calibration/FT/AttenuationMap" + std::to_string( m_irradiation ) + "fb_z_920-940"; + // Select the type of radiation map + if ( m_irradiationLinearModel ) conditionsLocation += "_linearmodel"; + if ( m_agingFibres ) conditionsLocation += "_aging"; + if ( m_replaceModules ) conditionsLocation += "_replacedmodule"; + info() << "Selecting " << std::to_string( m_irradiation ) << "/fb attenuation map using " + << std::string( m_irradiationLinearModel ? "linear" : "power-law" ) << " radiation damage model " + << std::string( m_agingFibres ? "with" : "without" ) << " aging fibres and " + << std::string( m_replaceModules ? "with" : "without" ) << " replaced center module. " << endmsg; + + addConditionDerivation( {conditionsLocation}, m_conditions.key(), [&]( YAML::Node const& cond ) -> Conditions { + // get parameters from condition DB + auto nBinsX = cond["x_n_bins"].as(); + auto nBinsY = cond["y_n_bins"].as(); + auto xEdges = cond["x_edges"].as>(); + auto yEdges = cond["y_edges"].as>(); + std::sort( xEdges.begin(), xEdges.end() ); + std::sort( yEdges.begin(), yEdges.end() ); + auto effDir = cond["eff_dir"].as>(); + auto effRef = cond["eff_refl"].as>(); + + // validate parameters + if ( nBinsX != xEdges.size() - 1 || nBinsY != yEdges.size() - 1 || nBinsX * nBinsY != effDir.size() || + nBinsX * nBinsY != effRef.size() ) { + throw GaudiException( "MCFTG4AttenuationTool::initialize", "Attenuation Map does not match bin counts", + StatusCode::FAILURE ); + } + + std::transform( effRef.begin(), effRef.end(), effRef.begin(), + [this]( double val ) { return val * m_mirrorReflectivity; } ); + return {nBinsX, nBinsY, xEdges, yEdges, effDir, effRef}; + } ); + return StatusCode::SUCCESS; + } ); } diff --git a/FT/FTDigitisation/src/MCFTG4AttenuationTool.h b/FT/FTDigitisation/src/MCFTG4AttenuationTool.h deleted file mode 100644 index 37625114a..000000000 --- a/FT/FTDigitisation/src/MCFTG4AttenuationTool.h +++ /dev/null @@ -1,58 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration * -* * -* This software is distributed under the terms of the GNU General Public * -* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". * -* * -* In applying this licence, CERN does not waive the privileges and immunities * -* granted to it by virtue of its status as an Intergovernmental Organization * -* or submit itself to any jurisdiction. * -\*****************************************************************************/ -#ifndef MCFTG4ATTENUATIONTOOL_H -#define MCFTG4ATTENUATIONTOOL_H 1 - -// from Gaudi -#include "GaudiAlg/GaudiTool.h" -#include "GaudiKernel/SystemOfUnits.h" - -#include "Kernel/STLExtensions.h" - -#include "IMCFTAttenuationTool.h" // Interface - -/** @class MCFTAttenuationTool MCFTAttenuationTool.h - * - * Tool that reads attenutation maps from ConDB - * - * - * @author Martin Bieker based on implementation of M. Demmer and J. Wishahi - * @date 2015-01-15 - */ - -class MCFTG4AttenuationTool : public extends { - -public: - using base_class::base_class; - - /// Initialize the transmission map - StatusCode initialize() override; - - /// Calculate the direct attenuation and the attenuation with reflection - IMCFTAttenuationTool::Attenuation attenuation( double x, double y ) const override; - -private: - bool validateMap(); - int findBin( LHCb::span axis, double position ) const; - - unsigned int m_nBinsX, m_nBinsY; - std::vector m_xEdges, m_yEdges, m_effDir, m_effRef; - - // properties - Gaudi::Property m_mirrorReflectivity{this, "MirrorReflectivity", 0.75, // from LHCb-PUB-2014-020 - "Reflectivity of mirror at the end of the fibre mat (0-1)"}; - Gaudi::Property m_irradiation{this, "Irradiation", 50, "Simulated radiation damage"}; - Gaudi::Property m_irradiationLinearModel{this, "IrradiationLinearModel", false, - "Radiation damage model (default is power-law model)"}; - Gaudi::Property m_agingFibres{this, "AgingFibres", false, "Simulate time aging of the fibres"}; - Gaudi::Property m_replaceModules{this, "ReplaceModules", false, "Simulate replacement of inner modules"}; -}; -#endif // MCFG4TATTENUATIONTOOL_H -- GitLab From 907a432324146807c948c6c31506cc0385ee4268 Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Tue, 4 Oct 2022 10:22:57 +0200 Subject: [PATCH 07/15] Modernized VPDepositCreator and switched to derived conditions --- VP/VPDigitisation/src/VPDepositCreator.cpp | 118 ++++++++------------- 1 file changed, 42 insertions(+), 76 deletions(-) diff --git a/VP/VPDigitisation/src/VPDepositCreator.cpp b/VP/VPDigitisation/src/VPDepositCreator.cpp index b3ec75e7f..42183cd90 100644 --- a/VP/VPDigitisation/src/VPDepositCreator.cpp +++ b/VP/VPDigitisation/src/VPDepositCreator.cpp @@ -8,32 +8,27 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ -#include "DetDesc/ConditionAccessor.h" -#include "DetDesc/ConditionAccessorHolder.h" #include "DetDesc/DetectorElement.h" #include "Detector/VP/VPChannelID.h" #include "Event/GenHeader.h" #include "Event/MCHit.h" #include "Event/MCVPDigit.h" -#include "GaudiAlg/GaudiTupleAlg.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "GaudiKernel/RndmGenerators.h" #include "Kernel/VPDataFunctor.h" +#include "LHCbAlgs/Transformer.h" #include "LHCbMath/LHCbMath.h" #include "VPDet/DeVP.h" #include "VPRadiationDamageTool.h" + +#include "GaudiAlg/GaudiTupleAlg.h" +#include "GaudiKernel/PhysicalConstants.h" +#include "GaudiKernel/RndmGenerators.h" + #include #include -/** @class VPDepositCreator - * Using MCHits as input, this Functional::Transformer simulates the spatial ionisation - * pattern as well as the drift and diffusion of charge carriers in the sensor - * and produces an MCVPDigit for each pixel with a charge deposit. - * - * @author Marcin Kucharczyk - * @date 20/09/09 - */ +#ifdef USE_DD4HEP +# include +#endif namespace VPDepositCreatorConditions { struct ConditionsCache { @@ -47,15 +42,22 @@ namespace VPDepositCreatorConditions { }; } // namespace VPDepositCreatorConditions +/** + * Using MCHits as input, this Functional::Transformer simulates the spatial ionisation + * pattern as well as the drift and diffusion of charge carriers in the sensor + * and produces an MCVPDigit for each pixel with a charge deposit. + * + * @author Marcin Kucharczyk + * @date 20/09/09 + */ class VPDepositCreator - : public Gaudi::Functional::Transformer< + : public LHCb::Algorithm::Transformer< LHCb::MCVPDigits( const LHCb::MCHits& mcPrevPrev, const LHCb::MCHits& mcPrev, const LHCb::MCHits& mcMain, const LHCb::MCHits& mcNext, const LHCb::GenHeader& genHeader, const DeVP&, const VPDepositCreatorConditions::ConditionsCache& ), LHCb::DetDesc::usesBaseAndConditions> { public: - /// Standard constructor VPDepositCreator( const std::string& name, ISvcLocator* pSvcLocator ); StatusCode initialize() override; ///< Algorithm initialization LHCb::MCVPDigits operator()( const LHCb::MCHits& mcPrevPrev, const LHCb::MCHits& mcPrev, const LHCb::MCHits& mcMain, @@ -169,9 +171,6 @@ private: DECLARE_COMPONENT( VPDepositCreator ) -//============================================================================= -// Constructor -//============================================================================= VPDepositCreator::VPDepositCreator( const std::string& name, ISvcLocator* pSvcLocator ) : Transformer{name, pSvcLocator, @@ -184,11 +183,7 @@ VPDepositCreator::VPDepositCreator( const std::string& name, ISvcLocator* pSvcLo KeyValue{"ConditionsCache", "AlgorithmSpecific-" + name + "-ConditionsCache"}}, {"OutputLocation", LHCb::MCVPDigitLocation::Default}} {} -//============================================================================= -// Initialisation -//============================================================================= StatusCode VPDepositCreator::initialize() { - return Transformer::initialize() .andThen( [&] { if ( !m_irradiated ) { @@ -208,70 +203,41 @@ StatusCode VPDepositCreator::initialize() { .andThen( [&] { return m_uniform.initialize( randSvc(), Rndm::Flat( 0.0, 1.0 ) ); } ) .andThen( [&] { addConditionDerivation( +#ifdef USE_DD4HEP + {inputLocation(), "Conditions/Calibration/VP/VPDeadTimeParam"}, + inputLocation(), + [=]( DeVP const& det, YAML::Node const& deadTimeParam ) -> VPDepositCreatorConditions::ConditionsCache { +#else {inputLocation()}, inputLocation(), - [=]( const DeVP& det ) { + [=]( DeVP const& det ) -> VPDepositCreatorConditions::ConditionsCache { +#endif // Calculate the diffusion coefficient. - const double kt = m_temperature * Gaudi::Units::k_Boltzmann / Gaudi::Units::eV; + const double kt = 2. * m_temperature * Gaudi::Units::k_Boltzmann / Gaudi::Units::eV; std::array tmpDiffusionCoefficient{}; det.runOnAllSensors( [&]( const DeVPSensor& sensor ) { - auto module = sensor.module(); - tmpDiffusionCoefficient[module] = sqrt( 2. * kt * sensor.siliconThickness() / ( sensor.sensorHV() ) ); - if ( msgLevel( MSG::VERBOSE ) ) { - verbose() << "Diffusion Coefficient module " << module << " = " << tmpDiffusionCoefficient[module] - << endmsg; - } + tmpDiffusionCoefficient[sensor.module()] = + sqrt( kt * sensor.siliconThickness() / ( sensor.sensorHV() ) ); } ); +#ifdef USE_DD4HEP + return {tmpDiffusionCoefficient, + deadTimeParam["FunctionParam"].as, 52>>()}; +#else std::array, 52> tmpDeadTimeParam{}; - if ( !existDet( detSvc(), "Conditions/Calibration/VP/VPDeadTimeParam" ) ) { - Warning( "VP deadtime parameters not in conditions DB", StatusCode::SUCCESS ).ignore(); - // expand data storage to match number of modules - for ( unsigned int module = 0; module < 52; ++module ) { - tmpDeadTimeParam[module][0] = s_defaultVals[module] * s_defaultDeadParam[0]; - tmpDeadTimeParam[module][1] = s_defaultDeadParam[1]; - tmpDeadTimeParam[module][2] = s_defaultDeadParam[2]; - tmpDeadTimeParam[module][3] = s_defaultDeadParam[3]; - } - } else { - // load from conditions DB - det.runOnAllSensors( [&]( const DeVPSensor& sensor ) { - auto module = sensor.module(); - std::string conditionName = - format( "Conditions/Calibration/VP/VPDeadTimeParam/Module%02i-DeadTime", module ); - if ( !existDet( detSvc(), conditionName ) ) { - Warning( "Condition directory exists but missing condition " + conditionName, StatusCode::FAILURE ) - .ignore(); - tmpDeadTimeParam[module][0] = s_defaultVals[module] * s_defaultDeadParam[0]; - tmpDeadTimeParam[module][1] = s_defaultDeadParam[1]; - tmpDeadTimeParam[module][2] = s_defaultDeadParam[2]; - tmpDeadTimeParam[module][3] = s_defaultDeadParam[3]; - } else { - auto cond = get( detSvc(), conditionName ); - if ( cond->param>( "FunctionParam" ).size() != 4 ) - throw GaudiException( - "Must be exactly 4 parameters in " + conditionName + "/FunctionParam found " + - std::to_string( cond->param>( "FunctionParam" ).size() ), - "VPDepositCreator", StatusCode::FAILURE ); - std::copy_n( cond->param>( "FunctionParam" ).begin(), 4, - tmpDeadTimeParam[module].begin() ); - } - if ( msgLevel( MSG::VERBOSE ) ) - verbose() << "Module " << module << " has deadtime parameters " - << format( "%f (%f exp(%f r) + (1-%f) exp(%f r))", tmpDeadTimeParam[module][0], - tmpDeadTimeParam[module][1], tmpDeadTimeParam[module][2], - tmpDeadTimeParam[module][1], tmpDeadTimeParam[module][3] ) - << endmsg; - } ); - } - return VPDepositCreatorConditions::ConditionsCache( std::move( tmpDiffusionCoefficient ), - std::move( tmpDeadTimeParam ) ); + det.runOnAllSensors( [&]( const DeVPSensor& sensor ) { + auto module = sensor.module(); + std::string conditionName = + format( "Conditions/Calibration/VP/VPDeadTimeParam/Module%02i-DeadTime", module ); + auto cond = get( detSvc(), conditionName ); + std::copy_n( cond->param>( "FunctionParam" ).begin(), 4, + tmpDeadTimeParam[module].begin() ); + } ); + return {tmpDiffusionCoefficient, tmpDeadTimeParam}; +#endif } ); setHistoTopDir( "VP/" ); } ); } -//========================================================================= -// Main execution -//========================================================================= LHCb::MCVPDigits VPDepositCreator::operator()( const LHCb::MCHits& mcPrevPrev, const LHCb::MCHits& mcPrev, const LHCb::MCHits& mcMain, const LHCb::MCHits& mcNext, const LHCb::GenHeader& genHeader, const DeVP& det, -- GitLab From 1cb61555b5212675c041dd8ea501eb2479103104 Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Tue, 4 Oct 2022 10:47:30 +0200 Subject: [PATCH 08/15] Modernized VPDigitCreator and switched to derived conditions --- VP/VPDigitisation/src/VPDigitCreator.cpp | 72 ++++++++---------------- 1 file changed, 23 insertions(+), 49 deletions(-) diff --git a/VP/VPDigitisation/src/VPDigitCreator.cpp b/VP/VPDigitisation/src/VPDigitCreator.cpp index ecfe6451d..d96aaf45c 100644 --- a/VP/VPDigitisation/src/VPDigitCreator.cpp +++ b/VP/VPDigitisation/src/VPDigitCreator.cpp @@ -8,36 +8,43 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ -#include "DetDesc/Condition.h" #include "Event/MCVPDigit.h" #include "Event/VPDigit.h" -#include "GaudiAlg/GaudiTupleAlg.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" #include "Kernel/VPConstants.h" #include "Kernel/VPDataFunctor.h" +#include "LHCbAlgs/Transformer.h" + +#include "GaudiAlg/GaudiTupleAlg.h" +#include "GaudiKernel/RndmGenerators.h" + #include #include +#include + using namespace LHCb; -/** @class VPDigitCreator VPDigitCreator.h +namespace { + using Threshold = std::array; +} // namespace + +/** * Using MCVPDigits as input, this algorithm simulates the response of the * VeloPix ASIC and produces a VPDigit for each MCVPDigit above threshold. * * @author Thomas Britton * @date 2010-07-07 */ - -class VPDigitCreator : public Gaudi::Functional::Transformer> { +class VPDigitCreator + : public LHCb::Algorithm::Transformer> { public: /// Standard constructor VPDigitCreator( const std::string& name, ISvcLocator* pSvcLocator ); StatusCode initialize() override; ///< Algorithm initialization - LHCb::VPDigits operator()( const LHCb::MCVPDigits& mcdigits ) const override; + LHCb::VPDigits operator()( LHCb::MCVPDigits const&, Threshold const& ) const override; private: /// Noise in number of electrons @@ -84,12 +91,6 @@ private: Gaudi::Property m_thresholdConditionName{this, "ThresholdCondition", "Conditions/Calibration/VP/VPDigitisationParam/Thresholds"}; - /// Discrimination threshold in number of electrons - std::array m_threshold; - - /// Default Discrimination threshold in number of electrons (if Condition missing) - static constexpr double s_defaultThreshold = 1000.0; - /// Random number generators mutable Rndm::Numbers m_gauss; mutable Rndm::Numbers m_uniform; @@ -98,18 +99,12 @@ private: DECLARE_COMPONENT( VPDigitCreator ) -//============================================================================= -// Standard constructor -//============================================================================= VPDigitCreator::VPDigitCreator( const std::string& name, ISvcLocator* pSvcLocator ) : Transformer{name, pSvcLocator, - {"InputLocation", LHCb::MCVPDigitLocation::Default}, + {{"InputLocation", LHCb::MCVPDigitLocation::Default}, {"Threshold", name + "_Threshold"}}, {"OutputLocation", LHCb::VPDigitLocation::Default}} {} -//============================================================================= -// Initialization -//============================================================================= StatusCode VPDigitCreator::initialize() { return Transformer::initialize() .andThen( [&] { return m_gauss.initialize( randSvc(), Rndm::Gauss( 0., 1. ) ); } ) @@ -124,37 +119,16 @@ StatusCode VPDigitCreator::initialize() { return m_poisson.initialize( randSvc(), Rndm::Poisson( averageNoisy ) ); } ) .andThen( [&] { - // set to default - for ( auto& t : m_threshold ) t = s_defaultThreshold; - if ( !existDet( detSvc(), m_thresholdConditionName.value() ) ) { - Warning( "VP threshold parameters not in conditions DB", StatusCode::SUCCESS ).ignore(); - } else { - std::string paramName{"ThresholdPerSensor"}; - auto cond = get( detSvc(), m_thresholdConditionName.value() ); - if ( cond->param>( paramName ).size() != 208 ) - throw GaudiException( "Must be exactly 208 values in " + m_thresholdConditionName.value() + '/' + paramName, - "VPDigitCreator", StatusCode::FAILURE ); - std::copy_n( cond->param>( paramName ).begin(), 208, m_threshold.begin() ); - } - if ( msgLevel( MSG::VERBOSE ) ) { - verbose() << "Sensors have a threshold in electrons of" << endmsg; - for ( unsigned int module = 0; module < 52; ++module ) { - verbose() << format( "Module %2i :", module ); - for ( unsigned int sens = 4 * module; sens < 4 * ( module + 1 ); ++sens ) { - verbose() << format( " %3i:%7.1f", sens, m_threshold[sens] ); - } - verbose() << endmsg; - } - } + addConditionDerivation( {m_thresholdConditionName}, inputLocation(), + [&]( YAML::Node const& cond ) -> Threshold { + return cond["ThresholdPerSensor"].as>(); + } ); return StatusCode( StatusCode::SUCCESS ); } ) .andThen( [&] { setHistoTopDir( "VP/" ); } ); } -//============================================================================= -// Main execution -//============================================================================= -LHCb::VPDigits VPDigitCreator::operator()( const LHCb::MCVPDigits& mcdigits ) const { +LHCb::VPDigits VPDigitCreator::operator()( LHCb::MCVPDigits const& mcdigits, Threshold const& threshold ) const { // Create a container for the digits and transfer ownership to the TES. LHCb::VPDigits digits; unsigned int numSpillover = 0; @@ -188,7 +162,7 @@ LHCb::VPDigits VPDigitCreator::operator()( const LHCb::MCVPDigits& mcdigits ) co 207.5, 0, 50000, 208, 100 ); } // Check if the collected charge exceeds the threshold. - if ( chargeMain < m_threshold[to_unsigned( channel.sensor() )] + m_electronicNoise * m_gauss() ) continue; + if ( chargeMain < threshold[to_unsigned( channel.sensor() )] + m_electronicNoise * m_gauss() ) continue; if ( m_monitoring ) { plot1D( chargeMain, "MainChargeAfter", "Pixel harge in main BX, e- after threshold", 0, 50000, 100 ); -- GitLab From 5be5d4b7bfd363525595981accbce9a2bde6e1e3 Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Mon, 17 Oct 2022 11:11:30 +0200 Subject: [PATCH 09/15] Fixed Muon code for dd4hep case --- Muon/MuonAlgs/src/MuonDigitization.cpp | 7 ++++++- Muon/MuonBackground/src/MuonBackground.cpp | 4 ++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/Muon/MuonAlgs/src/MuonDigitization.cpp b/Muon/MuonAlgs/src/MuonDigitization.cpp index 03764107e..fc7052e7a 100644 --- a/Muon/MuonAlgs/src/MuonDigitization.cpp +++ b/Muon/MuonAlgs/src/MuonDigitization.cpp @@ -253,14 +253,19 @@ void MuonDigitization::addChamberNoise() { // the readout number is essentially meaningless... // it is chamber noise..... if ( numberOfNoiseHit == 0 ) continue; +#ifdef USE_DD4HEP + DeMuonChamber p_chamber = m_muonDetector->getChamber( k, s, chamber ); +#else DeMuonChamber& p_chamber = m_muonDetector->getChamber( k, s, chamber ); +#endif for ( int hit = 1; hit <= numberOfNoiseHit; ++hit ) { // define position of hit const double x = m_flatDist() * m_muonDetector->getSensAreaX( k, s ); const double y = m_flatDist() * m_muonDetector->getSensAreaY( k, s ); const double time = m_flatDist() * m_BXTime; - const int gap = std::min( (int)( m_flatDist() * gapPerChamber ), m_muonDetector->gapsInRegion( k, s ) - 1 ); + const int gap = + std::min( (int)( m_flatDist() * gapPerChamber ), (int)m_muonDetector->gapsInRegion( k, s ) - 1 ); #ifdef USE_DD4HEP auto [point1, point2] = p_chamber.getGapPoints( gap, x, y ); #else diff --git a/Muon/MuonBackground/src/MuonBackground.cpp b/Muon/MuonBackground/src/MuonBackground.cpp index 2f83418b0..f91e95122 100755 --- a/Muon/MuonBackground/src/MuonBackground.cpp +++ b/Muon/MuonBackground/src/MuonBackground.cpp @@ -808,7 +808,11 @@ MuonBackground::getRadialPosition( int index, int station ) { float ypos = float( r * sin( globalPhi ) * Gaudi::Units::mm * m_unitLength ); try { +#ifdef USE_DD4HEP + auto chamber = m_muonDetector->pos2StChamber( xpos, ypos, station ); +#else auto& chamber = m_muonDetector->pos2StChamber( xpos, ypos, station ); +#endif // check n ot only that hit is inside chamber but also gap... if ( chamber.checkHitAndGapInChamber( xpos, ypos ) ) { // remember this is the chamber number inside a region.... -- GitLab From 8b5c028203a4249866d2b823713584005758f498 Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Tue, 18 Oct 2022 14:30:01 +0200 Subject: [PATCH 10/15] Started modernizing RichDigitQC did not finish as it's rather a complicated beast --- Rich/RichDigiQC/src/RichDigitQC.cpp | 194 +++++++++++++++++++--------- Rich/RichDigiQC/src/RichDigitQC.h | 165 ----------------------- 2 files changed, 136 insertions(+), 223 deletions(-) delete mode 100755 Rich/RichDigiQC/src/RichDigitQC.h diff --git a/Rich/RichDigiQC/src/RichDigitQC.cpp b/Rich/RichDigiQC/src/RichDigitQC.cpp index a6dc321b3..651a7879d 100644 --- a/Rich/RichDigiQC/src/RichDigitQC.cpp +++ b/Rich/RichDigiQC/src/RichDigitQC.cpp @@ -9,58 +9,159 @@ * or submit itself to any jurisdiction. * \*****************************************************************************/ -//------------------------------------------------------------------------------- -/** @file RichDigitQC.cpp - * +/** * Implementation file for RICH Digitisation Quality Control algorithm : RichDigitQC * * @author Chris Jones Christopher.Rob.Jones@cern.ch * @date 2003-09-08 */ -//------------------------------------------------------------------------------- -// local -#include "RichDigitQC.h" +#include "RichKernel/RichHistoAlgBase.h" + +#include "Event/MCRichDigit.h" + +#include "Kernel/RichDetectorType.h" +#include "Kernel/RichSmartID.h" + +#include "RichUtils/RichHashMap.h" +#include "RichUtils/RichMap.h" +#include "RichUtils/RichPoissonEffFunctor.h" +#include "RichUtils/RichStatDivFunctor.h" + +#include "RichDet/DeRichSystem.h" + +#include "MCInterfaces/IRichMCTruthTool.h" +#include "RichInterfaces/IRichSmartIDTool.h" + +#include "boost/format.hpp" +#include "boost/lexical_cast.hpp" + +#include + +namespace Rich::MC::Digi { + + /** @class DigitQC RichDigitQC.h + * + * Monitor for Rich digitisation and DAQ simulation + * + * @author Chris Jones (Christopher.Rob.Jones@cern.ch) + * @date 2003-09-08 + */ + class DigitQC final : public Rich::HistoAlgBase { + + public: + DigitQC( const std::string& name, ISvcLocator* pSvcLocator ); + + StatusCode initialize() override final; // Algorithm initialization + StatusCode execute() override final; // Algorithm execution + StatusCode finalize() override final; // Algorithm finalization + + private: + /// Returns the location of container for the MCRichHit associated to the given digit + std::string mchitLocation( const LHCb::MCRichDigit* digit ) const; + + /// Pointer to RICH system detector element + const DeRichSystem* m_richSys; + + ToolHandle m_smartIDs{this, "RichSmartIDTool", "Rich::SmartIDTool/RichSmartIDTool"}; + ToolHandle m_mcTool{this, "RichMCTruthTool", "Rich::MC::MCTruthTool/RichMCTruthTool"}; + + Gaudi::Property m_digitTDS{this, "InputDigits", LHCb::MCRichDigitLocation::Default, + "Location of MCRichDigits in TES"}; + + /// Number of events processed + unsigned long long int m_evtC{0}; + + /// Counter for hits in each HPD + typedef Rich::HashMap HPDCounter; + Rich::Map m_nHPDSignal; ///< Tally for HPD occupancy, in each RICH + + typedef Rich::HashMap SpillCount; + + typedef Rich::Map SpillDetCount; + + /// Number of digitised hits per RICH detector and event location + SpillDetCount m_spillDigits; + + /// Number of signal digitised hits per RICH detector and event location + SpillDetCount m_spillDigitsSignal; + + /// Number of raw MC hits hits per RICH detector and event location + SpillDetCount m_totalSpills; + + /// Counter for hit tallies + typedef std::map HitTally; + + /// Number of digitised hits in each RICH + HitTally m_allDigits; -//------------------------------------------------------------------------------- + /// Number of rayleigh scattered hits in each RICH + HitTally m_scattHits; -using namespace Rich::MC::Digi; + /// Number of charged track hits in each RICH + HitTally m_chrgTkHits; -DECLARE_COMPONENT( DigitQC ) + /// Number of gas quartz window CK hits in each RICH + HitTally m_gasQCK; + + /// Number of HPD quartz window CK hits in each RICH + HitTally m_hpdQCK; + + /// Number of nitrogen CK hits in each RICH + HitTally m_nitroQCK; + + /// Number of nitrogen CK hits in each RICH + HitTally m_aeroFiltQCK; + + /// Number of background hits in each RICH + HitTally m_bkgHits; + + /// Number of charge shared hits in each RICH + HitTally m_chrgShrHits; + + /// Number of HPD Reflection hits in each RICH + HitTally m_hpdReflHits; + + /// Number of silicon backscatter hits in each RICH + HitTally m_siBackScatt; + + /// Number of signal-induced noise hits in each RICH + HitTally m_signalInducedNoiseHits; + + /// List of event locations to look for MCRichHits in + typedef Rich::HashMap EventLocations; + EventLocations m_evtLocs; + }; + + inline std::string DigitQC::mchitLocation( const LHCb::MCRichDigit* digit ) const { + // Always just use the first hit, since this will always be signal if the digit has a signal contribution + // In case of signal-induced noise digit, return proper information ( no associated hits ) + return ( digit->hits().empty() ? ( digit->history().signalInducedNoise() ? "Digit due to SIN" : "UNKNOWN" ) + : objectLocation( digit->hits().front().mcRichHit()->parent() ) ); + } + + DECLARE_COMPONENT( DigitQC ) + +} // namespace Rich::MC::Digi // Standard constructor, initializes variables -DigitQC::DigitQC( const std::string& name, ISvcLocator* pSvcLocator ) - : Rich::HistoAlgBase( name, pSvcLocator ), m_richSys( 0 ), m_smartIDs( 0 ), m_mcTool( 0 ), m_evtC( 0 ) { +Rich::MC::Digi::DigitQC::DigitQC( const std::string& name, ISvcLocator* pSvcLocator ) + : Rich::HistoAlgBase( name, pSvcLocator ), m_richSys( 0 ), m_evtC( 0 ) { setProperty( "HistoDir", "DIGI/DIGITQC" ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); setProperty( "HistoOffSet", 10000 ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - // Declare job options - declareProperty( "InputDigits", m_digitTDS = LHCb::MCRichDigitLocation::Default ); - declareProperty( "ExtraHistos", m_extraHists = false ); } // Initialisation -StatusCode DigitQC::initialize() { - // Initialize base class - const StatusCode sc = Rich::HistoAlgBase::initialize(); - if ( sc.isFailure() ) { return sc; } - - // acquire tools - acquireTool( "RichSmartIDTool", m_smartIDs, NULL, true ); - acquireTool( "RichMCTruthTool", m_mcTool, NULL, true ); - - // RichDet - m_richSys = getDet( DeRichLocations::RichSystem ); - - // Warn if extra histos are enabled - if ( m_extraHists ) Warning( "Extra histograms are enabled", StatusCode::SUCCESS ).ignore(); - - return sc; +StatusCode Rich::MC::Digi::DigitQC::initialize() { + return Rich::HistoAlgBase::initialize().andThen( [&] { + // RichDet + m_richSys = getDet( DeRichLocations::RichSystem ); + return StatusCode::SUCCESS; + } ); } // Main execution -StatusCode DigitQC::execute() { - if ( msgLevel( MSG::DEBUG ) ) debug() << "Execute" << endmsg; - +StatusCode Rich::MC::Digi::DigitQC::execute() { // Locate MCRichDigits LHCb::MCRichDigits* richDigits = get( m_digitTDS ); @@ -72,7 +173,6 @@ StatusCode DigitQC::execute() { SpillDetCount spills; for ( LHCb::MCRichDigits::const_iterator iDigit = richDigits->begin(); iDigit != richDigits->end(); ++iDigit ) { const LHCb::MCRichDigit* mcDig = *iDigit; - // Get Rich ID const Rich::DetectorType rich = mcDig->key().rich(); @@ -151,12 +251,6 @@ StatusCode DigitQC::execute() { totDet += ( *iHPD ).second; } plot1D( totDet, RICH + " : Detector occupancy", 0, 15000, 150 ); - if ( m_extraHists ) { - plot1D( backs[rich], RICH + " : # Background hits", 0, 5000, 50 ); - for ( auto iC = spills[rich].begin(); iC != spills[rich].end(); ++iC ) { - plot1D( iC->second, RICH + " : # Spillover hits " + iC->first, 0, 5000, 100 ); - } - } // Event printout //------------------------------------------------------------------------------ @@ -173,7 +267,7 @@ StatusCode DigitQC::execute() { } // Finalize -StatusCode DigitQC::finalize() { +StatusCode Rich::MC::Digi::DigitQC::finalize() { // Statistical calculators const Rich::StatDivFunctor occ( "%8.2f +-%5.2f" ); @@ -195,11 +289,7 @@ StatusCode DigitQC::finalize() { unsigned int maxOcc( 0 ), minOcc( 999999 ); Rich::DAQ::PDHardwareID maxHPD( 0 ), minHPD( 0 ); for ( auto iHPD = m_nHPDSignal[rich].begin(); iHPD != m_nHPDSignal[rich].end(); ++iHPD ) { - const auto hID( m_richSys->hardwareID( ( *iHPD ).first ) ); - Gaudi::XYZPoint hpdGlo; - const auto ok = m_smartIDs->pdPosition( ( *iHPD ).first, hpdGlo ); - if ( !ok ) continue; - const auto hpdLoc( m_smartIDs->globalToPDPanel( hpdGlo ) ); + const auto hID( m_richSys->hardwareID( ( *iHPD ).first ) ); totDetSignal += ( *iHPD ).second; // make sure that the choice of maxHPD and minHPD is independent of the iteration order if ( iHPD->second >= maxOcc ) { @@ -210,18 +300,6 @@ StatusCode DigitQC::finalize() { if ( iHPD->second < minOcc || hID < minHPD ) { minHPD = hID; } minOcc = iHPD->second; } - - if ( m_extraHists ) { - plot2D( hpdLoc.x(), hpdLoc.y(), RICH + " : HPD hardware ID layout", -800, 800, -600, 600, 100, 100, - hID.data() ); - plot2D( hpdLoc.x(), hpdLoc.y(), RICH + " : SmartID Row layout", -800, 800, -600, 600, 100, 100, - ( *iHPD ).first.pdNumInCol() ); - plot2D( hpdLoc.x(), hpdLoc.y(), RICH + " : SmartID Col layout", -800, 800, -600, 600, 100, 100, - ( *iHPD ).first.pdCol() ); - } - debug() << " HPD " << ( *iHPD ).first; - debug() << " Global position : " << hpdGlo << endmsg << " Local position : " << hpdLoc << endmsg - << " Hit occupancy : " << occ( ( *iHPD ).second, m_evtC ) << " hits/event" << endmsg; } info() << " " << RICH << " : Av. total hit occupancy " << occ( m_allDigits[rich], m_evtC ) << " hits/event" diff --git a/Rich/RichDigiQC/src/RichDigitQC.h b/Rich/RichDigiQC/src/RichDigitQC.h deleted file mode 100755 index 445dde374..000000000 --- a/Rich/RichDigiQC/src/RichDigitQC.h +++ /dev/null @@ -1,165 +0,0 @@ -/*****************************************************************************\ -* (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. * -\*****************************************************************************/ - -//------------------------------------------------------------------------------------ -/** @file RichDigitQC.h - * - * Header file for RICH Digitisation Quality Control algorithm : Rich::MC::Digi::DigitQC - * - * @author Chris Jones Christopher.Rob.Jones@cern.ch - * @date 2003-09-08 - */ -//------------------------------------------------------------------------------------ - -#pragma once - -// STD -#include - -// Boost -#include "boost/lexical_cast.hpp" - -// base class -#include "RichKernel/RichHistoAlgBase.h" - -// Event model -#include "Event/MCRichDigit.h" - -// Kernel -#include "Kernel/RichDetectorType.h" -#include "Kernel/RichSmartID.h" - -// Rich Utils -#include "RichUtils/RichHashMap.h" -#include "RichUtils/RichMap.h" -#include "RichUtils/RichPoissonEffFunctor.h" -#include "RichUtils/RichStatDivFunctor.h" - -// RichDet -#include "RichDet/DeRichSystem.h" - -// RICH Interfaces -#include "MCInterfaces/IRichMCTruthTool.h" -#include "RichInterfaces/IRichSmartIDTool.h" - -// Boost -#include "boost/format.hpp" - -namespace Rich::MC::Digi { - - /** @class DigitQC RichDigitQC.h - * - * Monitor for Rich digitisation and DAQ simulation - * - * @author Chris Jones (Christopher.Rob.Jones@cern.ch) - * @date 2003-09-08 - */ - - class DigitQC final : public Rich::HistoAlgBase { - - public: - /// Standard constructor - DigitQC( const std::string& name, ISvcLocator* pSvcLocator ); - - StatusCode initialize() override final; // Algorithm initialization - StatusCode execute() override final; // Algorithm execution - StatusCode finalize() override final; // Algorithm finalization - - private: // methods - /// Returns the location of container for the MCRichHit associated to the given digit - std::string mchitLocation( const LHCb::MCRichDigit* digit ) const; - - private: // data - /// Pointer to RICH system detector element - const DeRichSystem* m_richSys; - - /// Pointer to RichSmartID tool - const Rich::ISmartIDTool* m_smartIDs; - - /// Pointer to MC truth tool - const Rich::MC::IMCTruthTool* m_mcTool = nullptr; - - // job options - std::string m_digitTDS; ///< Location of MCRichDigits in TES - bool m_extraHists; ///< Flag to turn on the production of additional histograms - - /// Number of events processed - unsigned long long int m_evtC{0}; - - /// Counter for hits in each HPD - typedef Rich::HashMap HPDCounter; - Rich::Map m_nHPDSignal; ///< Tally for HPD occupancy, in each RICH - - typedef Rich::HashMap SpillCount; - - typedef Rich::Map SpillDetCount; - - /// Number of digitised hits per RICH detector and event location - SpillDetCount m_spillDigits; - - /// Number of signal digitised hits per RICH detector and event location - SpillDetCount m_spillDigitsSignal; - - /// Number of raw MC hits hits per RICH detector and event location - SpillDetCount m_totalSpills; - - /// Counter for hit tallies - typedef std::map HitTally; - - /// Number of digitised hits in each RICH - HitTally m_allDigits; - - /// Number of rayleigh scattered hits in each RICH - HitTally m_scattHits; - - /// Number of charged track hits in each RICH - HitTally m_chrgTkHits; - - /// Number of gas quartz window CK hits in each RICH - HitTally m_gasQCK; - - /// Number of HPD quartz window CK hits in each RICH - HitTally m_hpdQCK; - - /// Number of nitrogen CK hits in each RICH - HitTally m_nitroQCK; - - /// Number of nitrogen CK hits in each RICH - HitTally m_aeroFiltQCK; - - /// Number of background hits in each RICH - HitTally m_bkgHits; - - /// Number of charge shared hits in each RICH - HitTally m_chrgShrHits; - - /// Number of HPD Reflection hits in each RICH - HitTally m_hpdReflHits; - - /// Number of silicon backscatter hits in each RICH - HitTally m_siBackScatt; - - /// Number of signal-induced noise hits in each RICH - HitTally m_signalInducedNoiseHits; - - /// List of event locations to look for MCRichHits in - typedef Rich::HashMap EventLocations; - EventLocations m_evtLocs; - }; - - inline std::string DigitQC::mchitLocation( const LHCb::MCRichDigit* digit ) const { - // Always just use the first hit, since this will always be signal if the digit has a signal contribution - // In case of signal-induced noise digit, return proper information ( no associated hits ) - return ( digit->hits().empty() ? ( digit->history().signalInducedNoise() ? "Digit due to SIN" : "UNKNOWN" ) - : objectLocation( digit->hits().front().mcRichHit()->parent() ) ); - } - -} // namespace Rich::MC::Digi \ No newline at end of file -- GitLab From 693d5d71d72db84e17e193a6330cf9e8a93025bd Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Tue, 18 Oct 2022 14:36:26 +0200 Subject: [PATCH 11/15] Do not use Rich__MC__Digi__DigitQC in DD4hep mode as it's not DD4hep compliant at this stage --- Digi/Boole/python/Boole/Configuration.py | 2 +- Rich/RichDigiQC/CMakeLists.txt | 8 +++++++- Rich/RichDigiQC/src/{ => legacy}/RichDigitQC.cpp | 0 3 files changed, 8 insertions(+), 2 deletions(-) rename Rich/RichDigiQC/src/{ => legacy}/RichDigitQC.cpp (100%) diff --git a/Digi/Boole/python/Boole/Configuration.py b/Digi/Boole/python/Boole/Configuration.py index fd4a51541..bcfe6322a 100755 --- a/Digi/Boole/python/Boole/Configuration.py +++ b/Digi/Boole/python/Boole/Configuration.py @@ -895,7 +895,7 @@ class Boole(LHCbConfigurableUser): FTLiteClusterMonitor() ] - if "Rich" in moniDets: + if "Rich" in moniDets and not UseDD4Hep: from Configurables import Rich__MC__Digi__DigitQC GaudiSequencer("MoniRichSeq").Members += [ Rich__MC__Digi__DigitQC("RiDigitQC") diff --git a/Rich/RichDigiQC/CMakeLists.txt b/Rich/RichDigiQC/CMakeLists.txt index ee919bd9b..d9e3071cc 100644 --- a/Rich/RichDigiQC/CMakeLists.txt +++ b/Rich/RichDigiQC/CMakeLists.txt @@ -13,10 +13,16 @@ Rich/RichDigiQC --------------- #]=======================================================================] +if (NOT USE_DD4HEP) + set(DetDescOnlyFiles + src/legacy/RichDigitQC.cpp + ) +endif() + gaudi_add_module(RichDigiQC SOURCES src/RichDigiDataObjVerifier.cpp - src/RichDigitQC.cpp + ${DetDescOnlyFiles} LINK Boost::headers LHCb::LHCbKernel diff --git a/Rich/RichDigiQC/src/RichDigitQC.cpp b/Rich/RichDigiQC/src/legacy/RichDigitQC.cpp similarity index 100% rename from Rich/RichDigiQC/src/RichDigitQC.cpp rename to Rich/RichDigiQC/src/legacy/RichDigitQC.cpp -- GitLab From 4dd815eb955a47599bfde3d5c5bcce8cf30cb940 Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Tue, 18 Oct 2022 14:42:13 +0200 Subject: [PATCH 12/15] Modernized RichDetailedFrontEndResponsePMT and used derived conditions to be DD4hep compatible --- .../python/RichDigiSys/Configuration.py | 2 +- .../RichDetailedFrontEndResponsePMT.cpp | 845 ++++++++++++------ .../RichDetailedFrontEndResponsePMT.h | 372 -------- 3 files changed, 568 insertions(+), 651 deletions(-) delete mode 100755 Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.h diff --git a/Rich/RichDigiSys/python/RichDigiSys/Configuration.py b/Rich/RichDigiSys/python/RichDigiSys/Configuration.py index 589a3b837..386ee0ab0 100755 --- a/Rich/RichDigiSys/python/RichDigiSys/Configuration.py +++ b/Rich/RichDigiSys/python/RichDigiSys/Configuration.py @@ -54,7 +54,7 @@ class RichDigiSysConf(LHCbConfigurableUser): "SinNrOfTimeSlots": 3, "SinOccScaleToCurrentNu": True, "TimeGate": True, - "TimeGateWindowBegin": [18.00, 48.00], + "TimeGateWindowBegin": (18.00, 48.00), "TimeGateLookupTableR1": [], "TimeGateLookupTableR2": [], "TestRawFormatDecoding": False diff --git a/Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.cpp b/Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.cpp index 60ddb3b42..f68210ed5 100755 --- a/Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.cpp +++ b/Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.cpp @@ -9,182 +9,503 @@ * or submit itself to any jurisdiction. * \*****************************************************************************/ -// Event +#include "RichChannelPropertiesPMT.h" + +#include "Kernel/RichDetectorType.h" +#include "Kernel/RichSmartID.h" +#include "RichKernel/RichHistoAlgBase.h" + +#ifdef USE_DD4HEP +# include "Detector/Rich/DeRich.h" +#else +# include "RichDet/DeRichPMT.h" +# include "RichDet/DeRichSystem.h" +#endif + +#include "RichInterfaces/IRichSmartIDTool.h" + +#include "RichUtils/RichHashMap.h" +#include "RichUtils/RichMap.h" + #include "Event/GenHeader.h" +#include "Event/MCRichDeposit.h" +#include "Event/MCRichDigit.h" +#include "Event/MCRichHit.h" +#include "Event/MCRichSummedDeposit.h" + +#include "Kernel/STLExtensions.h" +#include "LHCbAlgs/Transformer.h" + +#include "GaudiKernel/RndmGenerators.h" +#include "GaudiKernel/StdArrayAsProperty.h" + +#include +#include +#include + +#include + +namespace { + // properties for the time gate (all times in ns) + constexpr float TimeGateWindowRange = 25.0f; + constexpr float TimeGateSlotWidth = 3.125f; + constexpr unsigned int TimeGateNrOfSlots = 8; + typedef std::bitset<8> TimeGateSlots; + typedef std::vector TimeGateLookupTable; + + // SIN ratio given in the DB is measured for a certain time window range (SinRatioTimeRange); SIN hits are + // added with probability based on this measured ratio in a number of required time slots (separately in each + // time slot); SinNrOfTimeSlots corresponds to the number of time slots before the main time-gate window where + // the SIN should be simulated (including the main time slot) - TimeGate simulation needs to be activated to + // enable this (otherwise SIN will be simulated only for the main time slot) + constexpr float SinRatioTimeRange{25.0f}; + + // occupancies used for the SIN simulation are given in the DB for a reference nu (m_sinOccRefNu); they are + // scaled accordingly (sinOccScaleFactor) to the current (run-time) nu + constexpr float SinOccRefNu{7.60022f}; + + // parameters for the monitoring histograms + constexpr float HistInputSignalMax = 5.0; + constexpr float HistReadoutDelayMax = 20.0; + constexpr float HistTimeOverThresholdMax = 35.0; + constexpr float HistTimeOfAcquisitionMax = 100.0; + + constexpr unsigned int m_nrOfBinsInHistogram1D = 100; + constexpr unsigned int m_nrOfBinsInHistogram2D = 50; + + struct TimeInfoProperties { + TimeInfoProperties( IRndmGenSvc* randSvc, float _pdPhotoelectronTimeOfFlight, float _pdTransitTimeMeanTypeR, + float _pdTransitTimeSpreadTypeR, float _pdTransitTimeMeanTypeH, float _pdTransitTimeSpreadTypeH, + float _timeInfoIndexIntervalInElectrons, unsigned int _readoutTimeInfoIndexMin, + unsigned int _readoutTimeInfoIndexMax ) + : pdPhotoelectronTimeOfFlight{_pdPhotoelectronTimeOfFlight} + , timeInfoIndexIntervalInElectrons{_timeInfoIndexIntervalInElectrons} + , readoutTimeInfoIndexMin{_readoutTimeInfoIndexMin} + , readoutTimeInfoIndexMax{_readoutTimeInfoIndexMax} { + pdTransitTimeTypeR.initialize( randSvc, Rndm::Gauss( _pdTransitTimeMeanTypeR, _pdTransitTimeSpreadTypeR ) ) + .orThrow( "Failed to initialize Rndm", "TimeInfoProperties::TimeInfoProperties" ) + .andThen( [&] { + return pdTransitTimeTypeH.initialize( randSvc, + Rndm::Gauss( _pdTransitTimeMeanTypeH, _pdTransitTimeSpreadTypeH ) ); + } ) + .orThrow( "Failed to initialize Rndm", "TimeInfoProperties::TimeInfoProperties" ) + .andThen( [&] { return gaussForGain.initialize( randSvc, Rndm::Gauss( 0., 1. ) ); } ) + .orThrow( "Failed to initialize Rndm", "DetailedFrontEndResponsePMT::initialize" ) + .andThen( [&] { return flatForSinProb.initialize( randSvc, Rndm::Flat( 0., 1. ) ); } ) + .orThrow( "Failed to initialize Rndm", "DetailedFrontEndResponsePMT::initialize" ) + .andThen( [&] { return flatForSinTime.initialize( randSvc, Rndm::Flat( 0., SinRatioTimeRange ) ); } ) + .orThrow( "Failed to initialize Rndm", "DetailedFrontEndResponsePMT::initialize" ); + } + ~TimeInfoProperties() { + pdTransitTimeTypeR.finalize().ignore(); + pdTransitTimeTypeH.finalize().ignore(); + gaussForGain.finalize().ignore(); + flatForSinProb.finalize().ignore(); + flatForSinTime.finalize().ignore(); + } + // time of flight of a photoelectron in a PD (from photocathode to anode) in Gauss [ns] + float pdPhotoelectronTimeOfFlight{0.}; + float timeInfoIndexIntervalInElectrons{0.}; + unsigned int readoutTimeInfoIndexMin{0}; + unsigned int readoutTimeInfoIndexMax{0}; + +#ifndef USE_DD4HEP + std::vector> readoutDelay{}; + std::vector> readoutTimeOverThreshold{}; +#endif + Rndm::Numbers pdTransitTimeTypeR{}; + Rndm::Numbers pdTransitTimeTypeH{}; + Rndm::Numbers gaussForGain{}; + Rndm::Numbers flatForSinProb{}; + Rndm::Numbers flatForSinTime{}; + }; + +#ifdef USE_DD4HEP + using DeRichSystem = LHCb::Detector::DeRich; +#endif + +} // namespace + +namespace Rich::MC::Digi { + + /** + * @author Marcin Kucharczyk, Mariusz Witek + * @date 2015-10-08 + */ + class DetailedFrontEndResponsePMT final + : public LHCb::Algorithm::Transformer< + LHCb::MCRichDigits( LHCb::MCRichSummedDeposits const&, LHCb::GenHeader const&, DeRichSystem const&, + LHCb::RichSmartID::Vector const&, TimeInfoProperties const& ), + LHCb::DetDesc::usesBaseAndConditions> { + + public: + DetailedFrontEndResponsePMT( const std::string& name, ISvcLocator* pSvcLocator ) + : Transformer( name, pSvcLocator, + {{"MCRichSummedDepositsLocation", LHCb::MCRichSummedDepositLocation::Default}, + {"GenHeaderLocation", LHCb::GenHeaderLocation::Default}, + {"RichSystemLocation", DeRichLocations::RichSystem}, + {"ReadoutChannelsLocation", name + "_ReadoutChannelsLocation"}, + {"TimeInfoPropertiesLocation", name + "_TIPLocation"}}, + {"MCRichDigitsLocation", LHCb::MCRichDigitLocation::Default} ) { + setProperty( "HistoDir", "DIGI/PDRESPONSE" ).ignore(); + } -// Local -#include "RichDetailedFrontEndResponsePMT.h" + StatusCode initialize() override; + LHCb::MCRichDigits operator()( LHCb::MCRichSummedDeposits const&, LHCb::GenHeader const&, DeRichSystem const&, + LHCb::RichSmartID::Vector const&, TimeInfoProperties const& ) const override; + + private: + class PMTChannel; + typedef Rich::Map PMTChannelContainer; + PMTChannelContainer Analog( LHCb::MCRichSummedDeposits const&, LHCb::GenHeader const&, DeRichSystem const&, + LHCb::RichSmartID::Vector const&, TimeInfoProperties const& ) const; + LHCb::MCRichDigits Digital( DeRichSystem const&, TimeInfoProperties const&, PMTChannelContainer& ) const; + + // get gain based on channel properties (gain mean/rms); require that gain is higher than certain value + // (optionally) + float getChannelGain( DeRichSystem const&, TimeInfoProperties const&, LHCb::RichSmartID const&, + float minValue = 0. ) const; + // read data needed for time info simulation from the DB + TimeInfoProperties readTimeInfoPropertiesFromDB( YAML::Node const&, IRndmGenSvc* ) const; + // get time from signal arrival at the anode to it's acquisition at the readout output [ns] + float getReadoutDelay( TimeInfoProperties const&, float, float ) const; + // get time-over-threshold of the readout output signal after acquisition [ns] + float getReadoutTimeOverThreshold( TimeInfoProperties const&, float, float ) const; + // read the single photoelectron energy deposit from the DB ( must be the same number as used in Gauss to create + // the deposit ) + StatusCode readPhotoelectronEnergyFromDB(); + + private: + ServiceHandle m_randSvc{this, "RndmGenSvc", "RndmGenSvc"}; + ToolHandle m_smartIDs{this, "RichSmartIDTool", "Rich::SmartIDTool/RichSmartIDTool"}; + ToolHandle m_propertyTool{this, "RichChannelPropertiesPMT", "RichChannelPropertiesPMT"}; + + Gaudi::Property m_readParametersFromDB{this, "ReadParametersFromDB", true}; + Gaudi::Property m_sin{this, "Sin", false}; + Gaudi::Property m_spillover{this, "SpilloverModel", "None"}; + Gaudi::Property m_gainMean{this, "GainMean", 0.9}; + Gaudi::Property m_gainRms{this, "GainRms", 0.2}; + Gaudi::Property m_threshold{this, "Threshold", 0.125}; + + // time gate is applied to 'raw' deposit times with t0 at the collision time (no offsets) + // the beginning of the time-gate window and the lookup table are given separately for R1 / R2 + // if the lookup table is not provided, a 25 ns time window is used by default + Gaudi::Property m_timeGate{this, "TimeGate", true}; + Gaudi::Property> m_timeGateWindowBegin{this, "TimeGateWindowBegin", {18.00, 48.00}}; + Gaudi::Property> m_timeGateLookupTableR1{this, "TimeGateLookupTableR1", {}}; + Gaudi::Property> m_timeGateLookupTableR2{this, "TimeGateLookupTableR2", {}}; + + // end of the time-gate window [ns] (based on the begin and range) + std::array m_timeGateWindowEnd{0., 0.}; + + // time-gate lookup tables (parsed to a proper format from a string from the Gaudi property) + TimeGateLookupTable m_lookupTableR1{}; + TimeGateLookupTable m_lookupTableR2{}; + + // SIN ratio given in the DB is measured for a certain time window range (SinRatioTimeRange); SIN hits are + // added with probability based on this measured ratio in a number of required time slots (separately in each + // time slot); SinNrOfTimeSlots corresponds to the number of time slots before the main time-gate window where + // the SIN should be simulated (including the main time slot) - TimeGate simulation needs to be activated to + // enable this (otherwise SIN will be simulated only for the main time slot) + Gaudi::Property m_sinNrOfTimeSlots{this, "SinNrOfTimeSlots", 3}; + + // occupancies used for the SIN simulation are given in the DB for a reference nu (m_sinOccRefNu); they are + // scaled accordingly (sinOccScaleFactor) to the current (run-time) nu + Gaudi::Property m_sinOccScaleToCurrentNu{this, "SinOccScaleToCurrentNu", true}; + mutable float m_current_sinOccScaleFactor{0.}; // only used to display values nly once in the log + mutable std::mutex m_sinOccScaleFactorMutex; + + // energy of a single photoelectron deposit from Gauss [MeV] + float m_photoelectronEnergy{0.}; + + // structure for deposits (not the same as RichDeposit, it is what is effectively seen by the Claro - can be a + // SIN hit, or two deposits close in time which are seen as a single one) + struct PMTDeposit { + float inputSignal{0.}; + float timeOfAcquisition{0.}; + float timeOverThreshold{0.}; + }; + + // helper class to process all deposits within a single PMT channel (including SIN) + class PMTChannel { + + private: + typedef std::vector PMTDepositContainer; + PMTDepositContainer m_deposits{}; + TimeGateSlots m_outputTimeSampling{}; + bool m_hasSin{false}; + LHCb::MCRichSummedDeposit* m_summedDeposit{nullptr}; + typedef std::vector TimeGatePatternMatches; + TimeGatePatternMatches m_timeGatePatternMatches{}; + + public: + inline void addDeposit( const float inputSignal, // + const float timeOfAcquisition, // + const float timeOverThreshold ) { + PMTDeposit deposit; + deposit.inputSignal = inputSignal; + deposit.timeOfAcquisition = timeOfAcquisition; + deposit.timeOverThreshold = timeOverThreshold; + m_deposits.push_back( deposit ); + } -using namespace Rich::MC::Digi; + // add the given deposit (with a certain time-of-acquisition and time-over-threshold) to the Claro output + // pattern + inline void processDeposit( const float timeOfAcquisition, // + const float timeOverThreshold, // + const float timeGateMin, // + const float timeGateMax, // + const float timeGateBin ) { -//----------------------------------------------------------------------------- -// Implementation file for class : RichDetailedFrontEndResponsePMT -// 2015-10-08 : Marcin Kucharczyk, Mariusz Witek -// 2017-07-22: Bartosz Malecki -//----------------------------------------------------------------------------- -DECLARE_COMPONENT( DetailedFrontEndResponsePMT ) + // check if the deposit has any overlap with the time-sampling window + if ( isOverlapBetweenSignalAndTimeSampling( timeOfAcquisition, timeOverThreshold, timeGateMin, timeGateMax ) ) { -DetailedFrontEndResponsePMT::DetailedFrontEndResponsePMT( const std::string& name, ISvcLocator* pSvcLocator ) - : Rich::HistoAlgBase( name, pSvcLocator ) { - setProperty( "HistoDir", "DIGI/PDRESPONSE" ).ignore(); -} + const auto firstTimeSlotIndex = getTimeSlotIndex( timeOfAcquisition, timeGateMin, timeGateBin ); -StatusCode DetailedFrontEndResponsePMT::initialize() { - StatusCode sc = Rich::HistoAlgBase::initialize(); - if ( sc.isFailure() ) { return sc; } - // Create a collection of all valid readout channels (pixels) - const Rich::ISmartIDTool* smartIDs( nullptr ); - acquireTool( "RichSmartIDTool", smartIDs, 0, true ); - m_readoutChannelsList = smartIDs->readoutChannelList(); - _ri_debug << "Retrieved " << m_readoutChannelsList.size() << " pixels in active list" << endmsg; - releaseTool( smartIDs ); - // Load parameters from database - m_propertyTool = tool( "RichChannelPropertiesPMT", this ); - sc = m_propertyTool->ReadChannelParameters( m_readoutChannelsList ); - if ( !sc ) return Error( "Failed to read channel parameters", sc ); - - // Get DeRichSystem with global RICH info - m_deRichSystem = getDetIfExists( DeRichLocations::RichSystem ); - if ( !m_deRichSystem ) { - error() << "Can't load the RichSystem detector element from location: " << DeRichLocations::RichSystem << endmsg; - return StatusCode::FAILURE; - } + const auto lastTimeSlotIndex = + getTimeSlotIndex( timeOfAcquisition + timeOverThreshold, timeGateMin, timeGateBin ); - // load time info from DB - sc = readTimeInfoPropertiesFromDB(); - if ( !sc ) { return Error( "Failed to read time info properties from the DB.", sc ); } + // set time slots status within the time-sampling window; time slots ordering is reversed wrt to the + // 'bitset' format + for ( auto i = firstTimeSlotIndex; i < lastTimeSlotIndex + 1; ++i ) { + if ( i >= 0 && (std::size_t)i < m_outputTimeSampling.size() ) { + m_outputTimeSampling.set( m_outputTimeSampling.size() - 1 - i ); + } + } + } + } - // load single photoelectron deposit energy from DB - sc = readPhotoelectronEnergyFromDB(); - if ( !sc ) { return Error( "Failed to read photoelectron deposit energy from the DB.", sc ); } + // apply a time-gate on the Claro output + inline bool applyTimeGate( const TimeGateLookupTable& timeGateLookupTable ); - // check if the channel properties can be read from the DB - if ( m_readParametersFromDB.value() ) { + inline const TimeGateSlots& getOutputTimeSampling() const { return m_outputTimeSampling; } + inline const PMTDepositContainer& getDeposits() const { return m_deposits; } - const auto* firstDePMT = getDePmtFromSmartId( m_readoutChannelsList[0] ); - if ( !firstDePMT->hasCondition( DeRichLocations::PMTPropertiesCondName ) ) { + inline LHCb::MCRichSummedDeposit* getSummedDeposit() const { return m_summedDeposit; } + inline void setSummedDeposit( LHCb::MCRichSummedDeposit* summedDeposit ) { m_summedDeposit = summedDeposit; } - // force to use default values (instead of reading them from DB) if the relevant condition in DB doesn't exist - m_readParametersFromDB.setValue( false ); - _ri_debug << "Channel parameters not available in this DB version. Falling back to use the default values." - << endmsg; - } - // in particular, check if the SIN info is available (it is done here to avoid reading info from a particular PMT - // detector element again in the section relevant for m_sin) - else if ( !( firstDePMT->condition( DeRichLocations::PMTPropertiesCondName ) - ->exists( DeRichLocations::PMTSinInfoParamName ) ) ) { - m_sin.setValue( false ); - _ri_debug << "SIN info not available in this DB version. SIN will not be simulated." << endmsg; - } - } + inline bool hasSin() const noexcept { return m_hasSin; } + inline void setHasSin() { m_hasSin = true; } - // validate input for time gate - if ( m_timeGate.value() ) { + inline const TimeGatePatternMatches& getTimeGatePatternMatches() const { return m_timeGatePatternMatches; } - if ( !m_deRichSystem->hasCondition( DeRichLocations::ReadoutTimeInfoCondName ) ) { - m_timeGate.setValue( false ); - _ri_debug << "Can't simulate time gate without access to the time info in the DB. Unsetting this option." - << endmsg; - } + private: + // check if the given signal has any overlap (in time) with the time-sampling window + inline bool isOverlapBetweenSignalAndTimeSampling( const float timeOfAcquisition, // + const float timeOverThreshold, // + const float timeGateMin, // + const float timeGateMax ) const { + return ( ( timeOfAcquisition <= timeGateMax ) && ( timeOfAcquisition + timeOverThreshold >= timeGateMin ) ); + } - // validate input params - if ( m_timeGateWindowBegin.value().size() != 2 ) { - m_timeGate.setValue( false ); - _ri_debug << "Unexpected size of the vectors with time gate properties. The vectors should have 2 elements (one " - "param for each of R1/R2). Time gate will not be simulated." - << endmsg; - } else { + // get index of the time slot to which the given time corresponds + inline int getTimeSlotIndex( const float time, const float timeGateMin, const float timeGateBin ) const { + const auto var = ( 0 != timeGateBin ? ( time - timeGateMin ) / timeGateBin : 0 ); + // mimics std::floor, but avoid uses it as throws strange FPE with clang.. + auto j = static_cast( var ); + if ( var < 0 ) { --j; } + return j; + } - // set the end of the time-gate window (based on the beginning and range) - m_timeGateWindowEnd[0] = m_timeGateWindowBegin[0] + m_timeGateWindowRange; - m_timeGateWindowEnd[1] = m_timeGateWindowBegin[1] + m_timeGateWindowRange; + // process a single line of the lookup table + inline bool processLookupTableLine( const TimeGateSlots& line, const TimeGateSlots& outputTimeSampling ) const { + return ( outputTimeSampling == line ); + } + }; - // parse and validate the lookup tables - if ( !parseLookupTable( m_timeGateLookupTableR1.value(), m_lookupTableR1 ) || - !parseLookupTable( m_timeGateLookupTableR2.value(), m_lookupTableR2 ) ) { - m_timeGate.setValue( false ); - _ri_debug - << "Time-gate configuration is not right (the number of characters in each lookup table pattern should " - "match the number of slots: " - << m_timeGateNrOfSlots - << "). Time gate " - "will not be simulated." - << endmsg; + // helper functions + inline const DeRichPMT* getDePmtFromSmartId( DeRichSystem const& deRich, LHCb::RichSmartID const smartID ) const { + return (DeRichPMT*)deRich.dePD( smartID ); + } + inline float getChannelGainMean( DeRichSystem const& deRich, LHCb::RichSmartID const smartID ) const { + return m_readParametersFromDB.value() ? getDePmtFromSmartId( deRich, smartID )->PmtChannelGainMean( smartID ) + : m_gainMean.value(); + } + inline float getChannelGainRms( DeRichSystem const& deRich, LHCb::RichSmartID const smartID ) const { + return m_readParametersFromDB.value() ? getDePmtFromSmartId( deRich, smartID )->PmtChannelGainRms( smartID ) + : m_gainRms.value(); + } + inline float getChannelThreshold( DeRichSystem const& deRich, LHCb::RichSmartID const smartID ) const { + return m_readParametersFromDB.value() ? getDePmtFromSmartId( deRich, smartID )->PmtChannelThreshold( smartID ) + : m_threshold.value(); + } + inline float getChannelSinRatio( DeRichSystem const& deRich, LHCb::RichSmartID const smartID ) const { + return m_readParametersFromDB.value() ? getDePmtFromSmartId( deRich, smartID )->PmtChannelSinRatio( smartID ) + : 0.; + } + inline float getPmtAverageOccupancy( DeRichSystem const& deRich, LHCb::RichSmartID const smartID, + float sinOccScaleFactor ) const { + return m_readParametersFromDB.value() + ? ( m_sinOccScaleToCurrentNu + ? getDePmtFromSmartId( deRich, smartID )->PmtAverageOccupancy() * sinOccScaleFactor + : getDePmtFromSmartId( deRich, smartID )->PmtAverageOccupancy() ) + : 0.; + } + + inline void updateHistoryWithSin( LHCb::MCRichDigit*& digit ) const { + auto digitHistory = digit->history(); + digitHistory.setSignalInducedNoise( true ); + digit->setHistory( digitHistory ); + } + + // provide a valid index for time info input set (take the closest one from the DB) + inline unsigned int validTimeInfoInputIndex( TimeInfoProperties const& tip, float thresholdInElectrons ) const { + auto inputIndex = round( thresholdInElectrons / tip.timeInfoIndexIntervalInElectrons ); + if ( inputIndex < tip.readoutTimeInfoIndexMin ) { + warning() << "No reliable time info for this threshold setting. The closest time info input will be used." + << endmsg; + inputIndex = tip.readoutTimeInfoIndexMin; + } else if ( inputIndex > tip.readoutTimeInfoIndexMax ) { + warning() << "No reliable time info for this threshold setting. The closest time info input will be used." + << endmsg; + inputIndex = tip.readoutTimeInfoIndexMax; } + // offset the index ( should start from 0 ) + inputIndex -= tip.readoutTimeInfoIndexMin; + return inputIndex; } - } - if ( m_sin.value() ) { + // provide time of acquisition at the Claro output (taking all offsets into account) + inline float timeOfAcquisition( TimeInfoProperties const& tip, float timeOfDeposit, float inputSignal, + float threshold, bool isLargePMT ) const { + const auto pdTransitTime = isLargePMT ? tip.pdTransitTimeTypeH() : tip.pdTransitTimeTypeR(); + return ( timeOfDeposit - tip.pdPhotoelectronTimeOfFlight + pdTransitTime + + getReadoutDelay( tip, inputSignal, threshold ) ); + } - // don't simulate SIN if the channel properties from DB are not available - if ( !m_readParametersFromDB.value() ) { + // provide time of acquisition at the readout output for SIN hits ( random time within the time window of + // interest ) + inline float timeOfAcquisitionForSin( TimeInfoProperties const& tip, const float timeWindowMin ) const { + return timeWindowMin + tip.flatForSinTime(); + } - m_sin.setValue( false ); - _ri_debug << "Can't simulate SIN without access to the channel properties in the DB. Unsetting this option." - << endmsg; + // return a proper number if given timeGate property is different for R1/R2 + inline float timeGateParam( LHCb::RichSmartID const smartID, LHCb::span timeGateParamProperty ) const { + const auto rich = smartID.rich(); + if ( rich == Rich::Rich1 ) { + return timeGateParamProperty[0]; + } else if ( rich == Rich::Rich2 ) { + return timeGateParamProperty[1]; + } else { + warning() << "Given smartID not in R1/R2." << endmsg; + return 0.; + } } - // sanitise SIN time simulation settings - else if ( !m_timeGate.value() ) { - m_sinNrOfTimeSlots.setValue( 1 ); - _ri_debug << "Time-gate simulation is not enabled. SIN will be simulated only in the main time slot." << endmsg; + + // parse & validate the patterns in a lookup table + inline bool parseLookupTable( const std::vector& lookupTableFrom, + TimeGateLookupTable& lookupTableTo ) const { + bool isValid = true; + for ( auto const& line : lookupTableFrom ) { + if ( line.length() != TimeGateNrOfSlots ) { + isValid = false; + break; + } else { + lookupTableTo.emplace_back( line ); + } + } + return isValid; } - } + }; - sc = m_gaussForGain.initialize( randSvc(), Rndm::Gauss( 0., 1. ) ); - if ( sc.isFailure() ) return Error( "Failed to initialize Rndm", sc ); - sc = m_flatForSinProb.initialize( randSvc(), Rndm::Flat( 0., 1. ) ); - if ( sc.isFailure() ) return Error( "Failed to initialize Rndm", sc ); - sc = m_flatForSinTime.initialize( randSvc(), Rndm::Flat( 0., m_sinRatioTimeRange ) ); - if ( sc.isFailure() ) return Error( "Failed to initialize Rndm", sc ); - sc = - m_pdTransitTimeTypeR.initialize( randSvc(), Rndm::Gauss( m_pdTransitTimeMeanTypeR, m_pdTransitTimeSpreadTypeR ) ); - if ( sc.isFailure() ) return Error( "Failed to initialize Rndm", sc ); - sc = - m_pdTransitTimeTypeH.initialize( randSvc(), Rndm::Gauss( m_pdTransitTimeMeanTypeH, m_pdTransitTimeSpreadTypeH ) ); - if ( sc.isFailure() ) return Error( "Failed to initialize Rndm", sc ); - - return sc; -} + DECLARE_COMPONENT( DetailedFrontEndResponsePMT ) -StatusCode DetailedFrontEndResponsePMT::finalize() { - const auto sc = ( m_gaussForGain.finalize() && m_flatForSinProb.finalize() && m_flatForSinTime.finalize() && - m_pdTransitTimeTypeR.finalize() && m_pdTransitTimeTypeH.finalize() ); - if ( !sc ) { Warning( "Failed to finalize random number generators" ).ignore(); } - return Rich::HistoAlgBase::finalize(); -} +} // namespace Rich::MC::Digi -StatusCode DetailedFrontEndResponsePMT::execute() { - _ri_debug << "Execute" << endmsg; - m_sDeposits = get( m_mcRichSumDepsLocation ); - _ri_debug << "Successfully located " << m_sDeposits->size() << " MCRichSummedDeposits at " << m_mcRichSumDepsLocation - << endmsg; +//----------------------------------------------------------------------------- +// Implementation file for class : RichDetailedFrontEndResponsePMT +// 2015-10-08 : Marcin Kucharczyk, Mariusz Witek +// 2017-07-22: Bartosz Malecki +//----------------------------------------------------------------------------- - // estimate factor of occupancy increase wrt to the one provided for a reference nu (to be tried once) - if ( m_sinOccScaleToCurrentNu && !m_sinOccScaleFactorIsInit ) { - if ( initSinOccScaleFactor().isFailure() ) { - return Error( "Failed to initialize scale factor for the SIN occupancy." ); - } else { - m_sinOccScaleFactorIsInit = true; +StatusCode Rich::MC::Digi::DetailedFrontEndResponsePMT::initialize() { + return Transformer::initialize().andThen( [&]() -> StatusCode { + addConditionDerivation( {std::string{m_smartIDs->getPDPanelsPath()}}, inputLocation(), + [&]( RichPDPanels const& panels ) { // unused but smartIDs depends on them + // Create a collection of all valid readout channels (pixels) + LHCb::RichSmartID::Vector readoutChannelsList = m_smartIDs->readoutChannelList( panels ); + // Load parameters from database + m_propertyTool->ReadChannelParameters( readoutChannelsList ) + .orThrow( "Failed to read channel parameters", + "DetailedFrontEndResponsePMT::initialize" ); + return readoutChannelsList; + } ); + + addConditionDerivation( + {DeRichLocations::ReadoutTimeInfoCondPath}, inputLocation(), + [&]( YAML::Node const& cond ) { return this->readTimeInfoPropertiesFromDB( cond, &*m_randSvc ); } ); + + // load single photoelectron deposit energy from DB + readPhotoelectronEnergyFromDB().orThrow( "Failed to read photoelectron deposit energy from the DB.", + "DetailedFrontEndResponsePMT::initialize" ); + + // validate input for time gate + if ( m_timeGate.value() ) { + // set the end of the time-gate window (based on the beginning and range) + m_timeGateWindowEnd[0] = m_timeGateWindowBegin[0] + TimeGateWindowRange; + m_timeGateWindowEnd[1] = m_timeGateWindowBegin[1] + TimeGateWindowRange; + // parse and validate the lookup tables + if ( !parseLookupTable( m_timeGateLookupTableR1.value(), m_lookupTableR1 ) || + !parseLookupTable( m_timeGateLookupTableR2.value(), m_lookupTableR2 ) ) { + throw GaudiException( + "DetailedFrontEndResponsePMT::initialize", + fmt::format( + "Time-gate configuration is not right (the number of characters in each lookup table pattern should " + "match the number of slots: {})", + TimeGateNrOfSlots ), + StatusCode::FAILURE ); + } } - } - // container for processing channels with any deposits (signal, SIN, etc.) - PMTChannelContainer channels; + if ( m_sin.value() ) { + // check channel properties from DB are available if we simulate SIN + if ( !m_readParametersFromDB.value() ) { + throw GaudiException( + "DetailedFrontEndResponsePMT::initialize", + "Can't simulate SIN without access to the channel properties in the DB. Unsetting this option.", + StatusCode::FAILURE ); + } + // sanitise SIN time simulation settings + if ( !m_timeGate.value() ) { + m_sinNrOfTimeSlots.setValue( 1 ); + _ri_debug << "Time-gate simulation is not enabled. SIN will be simulated only in the main time slot." << endmsg; + } + } + return StatusCode::SUCCESS; + } ); +} +LHCb::MCRichDigits Rich::MC::Digi::DetailedFrontEndResponsePMT:: + operator()( LHCb::MCRichSummedDeposits const& deposits, LHCb::GenHeader const& genHeader, DeRichSystem const& deRich, + LHCb::RichSmartID::Vector const& readoutChannelsList, TimeInfoProperties const& tip ) const { // Analog simulation - StatusCode sc = Analog( channels ); - if ( sc.isFailure() ) { return sc; } + PMTChannelContainer channels = Analog( deposits, genHeader, deRich, readoutChannelsList, tip ); // Digital simulation - sc = Digital( channels ); - - return sc; + return Digital( deRich, tip, channels ); } -StatusCode DetailedFrontEndResponsePMT::Analog( PMTChannelContainer& channels ) { - _ri_debug << "Analog Simulation" << endmsg; - +Rich::MC::Digi::DetailedFrontEndResponsePMT::PMTChannelContainer Rich::MC::Digi::DetailedFrontEndResponsePMT::Analog( + LHCb::MCRichSummedDeposits const& deposits, LHCb::GenHeader const& genHeader, DeRichSystem const& deRich, + LHCb::RichSmartID::Vector const& readoutChannelsList, TimeInfoProperties const& tip ) const { + PMTChannelContainer channels; + // estimate factor of occupancy increase wrt to the one provided for a reference nu (to be tried once) + float sinOccScaleFactor{0.}; + if ( m_sinOccScaleToCurrentNu ) { + sinOccScaleFactor = static_cast( genHeader.beamParameters()->nu() ) / SinOccRefNu; + if ( m_current_sinOccScaleFactor != sinOccScaleFactor ) { + std::scoped_lock l( m_sinOccScaleFactorMutex ); + if ( m_current_sinOccScaleFactor != sinOccScaleFactor ) { + m_current_sinOccScaleFactor = sinOccScaleFactor; + info() << "Occupancy for the SIN simulation will be scaled by a factor of " << std::fixed + << std::setprecision( 4 ) << sinOccScaleFactor + << " wrt the reference one provided for nu = " << std::fixed << std::setprecision( 5 ) << SinOccRefNu + << "." << endmsg; + } + } + } // ANALOG RESPONSE - for ( auto const& summedDeposit : *m_sDeposits ) { + for ( auto const& summedDeposit : deposits ) { // for backward-compatibility with the spillover hits, don't process purely spillover summed deposits if the time // processing is not applied (it would add many irrelevant deposits) it correspond to the early spillover studies, @@ -200,7 +521,7 @@ StatusCode DetailedFrontEndResponsePMT::Analog( PMTChannelContainer& channels ) // values common for the channel const auto channelSmartID = channel.getSummedDeposit()->key().pixelID(); - const auto channelThreshold = getChannelThreshold( channelSmartID ); + const auto channelThreshold = getChannelThreshold( deRich, channelSmartID ); // process each deposit from the summed deposit for ( auto const& deposit : channel.getSummedDeposit()->deposits() ) { @@ -210,7 +531,7 @@ StatusCode DetailedFrontEndResponsePMT::Analog( PMTChannelContainer& channels ) // Claro input signal (number of electrons from the PMT [Me-]) // in principle the gain value can be different for each photoelectron (but leave this 'averaged' approach, unless // a more detailed description is available - those gain values will not be independent of each other) - const float inputSignal = nrOfPhotoelectrons * getChannelGain( channelSmartID ); + const float inputSignal = nrOfPhotoelectrons * getChannelGain( deRich, tip, channelSmartID ); // save only the deposits that will be visible at the Claro output (are above threshold) if ( inputSignal < channelThreshold ) continue; @@ -218,10 +539,10 @@ StatusCode DetailedFrontEndResponsePMT::Analog( PMTChannelContainer& channels ) // time properties of the deposit const auto depositTimeOfAcquisition = m_timeGate.value() - ? timeOfAcquisition( deposit->time(), inputSignal, channelThreshold, channelSmartID.isLargePMT() ) + ? timeOfAcquisition( tip, deposit->time(), inputSignal, channelThreshold, channelSmartID.isLargePMT() ) : 0.; const auto depositTimeOverThreshold = - m_timeGate.value() ? getReadoutTimeOverThreshold( inputSignal, channelThreshold ) : 0.; + m_timeGate.value() ? getReadoutTimeOverThreshold( tip, inputSignal, channelThreshold ) : 0.; channel.addDeposit( inputSignal, depositTimeOfAcquisition, depositTimeOverThreshold ); } @@ -232,24 +553,24 @@ StatusCode DetailedFrontEndResponsePMT::Analog( PMTChannelContainer& channels ) // add SIN deposits if ( m_sin.value() ) { - for ( auto const& readoutChannel : m_readoutChannelsList ) { + for ( auto const& readoutChannel : readoutChannelsList ) { // get SIN probability for the given channel - const float sinProbInChannel = - ( 1 - getPmtAverageOccupancy( readoutChannel ) ) * - ( 1 - exp( -getChannelSinRatio( readoutChannel ) * getPmtAverageOccupancy( readoutChannel ) ) ); + const float sinProbInChannel = ( 1 - getPmtAverageOccupancy( deRich, readoutChannel, sinOccScaleFactor ) ) * + ( 1 - exp( -getChannelSinRatio( deRich, readoutChannel ) * + getPmtAverageOccupancy( deRich, readoutChannel, sinOccScaleFactor ) ) ); for ( unsigned int i = 0; i < m_sinNrOfTimeSlots.value(); ++i ) { // add SIN with the givenl probability (all values are expected to be within 0-1 range) - if ( m_flatForSinProb() < sinProbInChannel ) { + if ( tip.flatForSinProb() < sinProbInChannel ) { // get properties of the SIN deposit - const auto channelThreshold = getChannelThreshold( readoutChannel ); + const auto channelThreshold = getChannelThreshold( deRich, readoutChannel ); // make sure that the inputSignal will pass the threshold (assume that we add SIN with some probability, which // should not be reduced by the threshold) - const auto inputSignal = getChannelGain( readoutChannel, channelThreshold ); + const auto inputSignal = getChannelGain( deRich, tip, readoutChannel, channelThreshold ); // if ( inputSignal < channelThreshold ) continue; // remember to apply threshold if the SIN inputSignal will // not be anymore a random number above threshold value from the definition @@ -257,11 +578,11 @@ StatusCode DetailedFrontEndResponsePMT::Analog( PMTChannelContainer& channels ) // TOA is a random number within time range of the given time slot const auto timeOfAcquisition = m_timeGate.value() - ? timeOfAcquisitionForSin( timeGateParam( readoutChannel, m_timeGateWindowBegin.value() ) - - i * m_sinRatioTimeRange ) + ? timeOfAcquisitionForSin( tip, timeGateParam( readoutChannel, m_timeGateWindowBegin.value() ) - + i * SinRatioTimeRange ) : 0.; const auto timeOverThreshold = - m_timeGate.value() ? getReadoutTimeOverThreshold( inputSignal, channelThreshold ) : 0.; + m_timeGate.value() ? getReadoutTimeOverThreshold( tip, inputSignal, channelThreshold ) : 0.; // add the SIN deposit to the existing channel or create a new one (and get an iterator to it) const auto channelIteratorAndIfIsNew = channels.emplace( readoutChannel, PMTChannel() ); @@ -271,16 +592,13 @@ StatusCode DetailedFrontEndResponsePMT::Analog( PMTChannelContainer& channels ) } } } - - return StatusCode::SUCCESS; + return channels; } -StatusCode DetailedFrontEndResponsePMT::Digital( PMTChannelContainer& channels ) { - - _ri_debug << "Digital simulation" << endmsg; - - auto mcRichDigits = new LHCb::MCRichDigits(); - put( mcRichDigits, m_mcRichDigitsLocation ); +LHCb::MCRichDigits Rich::MC::Digi::DetailedFrontEndResponsePMT::Digital( DeRichSystem const& deRich, + TimeInfoProperties const& tip, + PMTChannelContainer& channels ) const { + LHCb::MCRichDigits mcRichDigits{}; // Loop over all channels for ( auto& curChannel : channels ) { @@ -306,13 +624,16 @@ StatusCode DetailedFrontEndResponsePMT::Digital( PMTChannelContainer& channels ) if ( m_timeGate.value() ) { // values common for the channel - const auto threshold = getChannelThreshold( curChannel.first ); + const auto threshold = getChannelThreshold( deRich, curChannel.first ); const auto curTimeWindowBegin = timeGateParam( curChannel.first, m_timeGateWindowBegin.value() ); const auto curTimeWindowEnd = timeGateParam( curChannel.first, m_timeGateWindowEnd ); const auto rich = curChannel.first.rich(); // check to make sure for the following choices on R1 / R2 lookup table / histograms (is this necessary?) - if ( !( rich == Rich::Rich1 || rich == Rich::Rich2 ) ) { return Error( "Unrecognized RICH detector" ); } + if ( !( rich == Rich::Rich1 || rich == Rich::Rich2 ) ) { + throw GaudiException( "DetailedFrontEndResponsePMT::Digital", "Unrecognized RICH detector", + StatusCode::FAILURE ); + } const auto& curLookupTable = ( rich == Rich::Rich1 ? m_lookupTableR1 : m_lookupTableR2 ); @@ -322,25 +643,24 @@ StatusCode DetailedFrontEndResponsePMT::Digital( PMTChannelContainer& channels ) if ( produceHistos() ) { const auto richAsText = Rich::text( rich ); plot1D( deposit.inputSignal, richAsText + " : Input signal", richAsText + " : Claro input signal [Me-]", 0., - m_histInputSignalMax, m_nrOfBinsInHistogram1D ); + HistInputSignalMax, m_nrOfBinsInHistogram1D ); plot1D( deposit.timeOfAcquisition, richAsText + " : Time of acquisition", - richAsText + " : Deposit time of acquisition [ns]", 0., m_histTimeOfAcquisitionMax, + richAsText + " : Deposit time of acquisition [ns]", 0., HistTimeOfAcquisitionMax, m_nrOfBinsInHistogram1D ); plot1D( deposit.timeOverThreshold, richAsText + " : Time over threshold", - richAsText + " : Deposit time over threshold [ns]", 0., m_histTimeOverThresholdMax, + richAsText + " : Deposit time over threshold [ns]", 0., HistTimeOverThresholdMax, m_nrOfBinsInHistogram1D ); - plot2D( deposit.inputSignal, getReadoutDelay( deposit.inputSignal, threshold ), + plot2D( deposit.inputSignal, getReadoutDelay( tip, deposit.inputSignal, threshold ), richAsText + " : Delay vs Input signal", richAsText + " : Readout delay [ns] vs Input signal [Me-]", - 0., m_histInputSignalMax, 0., m_histReadoutDelayMax, m_nrOfBinsInHistogram2D, - m_nrOfBinsInHistogram2D ); + 0., HistInputSignalMax, 0., HistReadoutDelayMax, m_nrOfBinsInHistogram2D, m_nrOfBinsInHistogram2D ); plot2D( deposit.inputSignal, deposit.timeOverThreshold, richAsText + " : TOT vs Input signal", - richAsText + " : TOT [ns] vs Input signal [Me-]", 0., m_histInputSignalMax, 0., - m_histTimeOverThresholdMax, m_nrOfBinsInHistogram2D, m_nrOfBinsInHistogram2D ); + richAsText + " : TOT [ns] vs Input signal [Me-]", 0., HistInputSignalMax, 0., + HistTimeOverThresholdMax, m_nrOfBinsInHistogram2D, m_nrOfBinsInHistogram2D ); } // add the contribution from the given deposit to the Claro output pattern curChannel.second.processDeposit( deposit.timeOfAcquisition, deposit.timeOverThreshold, curTimeWindowBegin, - curTimeWindowEnd, m_timeGateSlotWidth ); + curTimeWindowEnd, TimeGateSlotWidth ); _ri_debug << "Readout output pattern: " << curChannel.second.getOutputTimeSampling() << endmsg; } @@ -366,7 +686,7 @@ StatusCode DetailedFrontEndResponsePMT::Digital( PMTChannelContainer& channels ) // create the MCRichDigit auto newDigit = new LHCb::MCRichDigit(); - mcRichDigits->insert( newDigit, curChannel.first ); + mcRichDigits.insert( newDigit, curChannel.first ); // check if summedDeposit exists (can be not the case for purely SIN channels) if ( curChannel.second.getSummedDeposit() != nullptr ) { @@ -384,23 +704,26 @@ StatusCode DetailedFrontEndResponsePMT::Digital( PMTChannelContainer& channels ) _ri_debug << newDigit->history() << " MCRichDigit" << endmsg; } - return StatusCode::SUCCESS; + return mcRichDigits; } //============================================================================= // Return a valid gain for pixel //============================================================================= -float DetailedFrontEndResponsePMT::getChannelGain( const LHCb::RichSmartID& smartID, const float& minValue ) { +float Rich::MC::Digi::DetailedFrontEndResponsePMT::getChannelGain( DeRichSystem const& deRich, + TimeInfoProperties const& tip, + LHCb::RichSmartID const& smartID, + float minValue ) const { - const auto channelGainMean = getChannelGainMean( smartID ); - const auto channelGainRms = getChannelGainRms( smartID ); + const auto channelGainMean = getChannelGainMean( deRich, smartID ); + const auto channelGainRms = getChannelGainRms( deRich, smartID ); // try to get valid value (maxNrTries times) float tmpGain; const int maxNrTries = 100; for ( int i = 0; i < maxNrTries; ++i ) { - tmpGain = m_gaussForGain() * channelGainRms + channelGainMean; + tmpGain = tip.gaussForGain() * channelGainRms + channelGainMean; if ( tmpGain > minValue ) break; @@ -417,102 +740,88 @@ float DetailedFrontEndResponsePMT::getChannelGain( const LHCb::RichSmartID& smar //============================================================================= // Readout time info //============================================================================= -StatusCode DetailedFrontEndResponsePMT::readTimeInfoPropertiesFromDB() { - - // check if the condition exists - if ( m_deRichSystem->hasCondition( DeRichLocations::ReadoutTimeInfoCondName ) ) { - - const auto cond = m_deRichSystem->condition( DeRichLocations::ReadoutTimeInfoCondName ); - m_pdPhotoelectronTimeOfFlight = cond->paramAsDouble( "PDPhotoelectronTimeOfFlight" ); - m_pdTransitTimeMeanTypeR = cond->paramAsDouble( "PDTransitTimeMeanTypeR" ); - m_pdTransitTimeSpreadTypeR = cond->paramAsDouble( "PDTransitTimeSpreadTypeR" ); - m_pdTransitTimeMeanTypeH = cond->paramAsDouble( "PDTransitTimeMeanTypeH" ); - m_pdTransitTimeSpreadTypeH = cond->paramAsDouble( "PDTransitTimeSpreadTypeH" ); - - m_timeInfoIndexIntervalInElectrons = cond->paramAsDouble( "TimeInfoIndexIntervalInElectrons" ); - - m_readoutTimeInfoIndexMin = cond->paramAsInt( "ReadoutTimeInfo_IndexMin" ); - m_readoutTimeInfoIndexMax = cond->paramAsInt( "ReadoutTimeInfo_IndexMax" ); - - const auto readoutDelayPropBaseName = cond->paramAsString( "ReadoutDelayCommonLoc" ); - const auto readoutToTPropBaseName = cond->paramAsString( "ReadoutTimeOverThresholdCommonLoc" ); - - for ( unsigned int i = m_readoutTimeInfoIndexMin; i < m_readoutTimeInfoIndexMax + 1; ++i ) { - - const auto iAsString = std::to_string( i ); - - SmartDataPtr readoutDelay( detSvc(), readoutDelayPropBaseName + iAsString ); - SmartDataPtr readoutToT( detSvc(), readoutToTPropBaseName + iAsString ); - - m_readoutDelay.push_back( - std::make_shared( readoutDelay, false, gsl_interp_cspline ) ); - m_readoutTimeOverThreshold.push_back( - std::make_shared( readoutToT, false, gsl_interp_cspline ) ); - } - - // check validity of the interpolators - for ( auto const& el : m_readoutDelay ) { - - if ( !el->valid() ) { - error() << "Incorrect time info format. Check properties with readout delay/time-over-threshold in the DB." - << endmsg; - return StatusCode::FAILURE; - } +TimeInfoProperties +Rich::MC::Digi::DetailedFrontEndResponsePMT::readTimeInfoPropertiesFromDB( YAML::Node const& cond, + IRndmGenSvc* randSvc ) const { + TimeInfoProperties tip{randSvc, + cond["PDPhotoelectronTimeOfFlight"].as(), + cond["PDTransitTimeMeanTypeR"].as(), + cond["PDTransitTimeSpreadTypeR"].as(), + cond["PDTransitTimeMeanTypeH"].as(), + cond["PDTransitTimeSpreadTypeH"].as(), + cond["TimeInfoIndexIntervalInElectrons"].as(), + cond["ReadoutTimeInfo_IndexMin"].as(), + cond["ReadoutTimeInfo_IndexMax"].as()}; + +#ifndef USE_DD4HEP + // FIXME tabulated properties have to be replaced by plain json arrays in DD4hep and this code adapted + const auto readoutDelayPropBaseName = cond["ReadoutDelayCommonLoc"].as(); + const auto readoutToTPropBaseName = cond["ReadoutTimeOverThresholdCommonLoc"].as(); + for ( unsigned int i = tip.readoutTimeInfoIndexMin; i < tip.readoutTimeInfoIndexMax + 1; ++i ) { + SmartDataPtr readoutDelay( detSvc(), fmt::format( "{}{}", readoutDelayPropBaseName, i ) ); + SmartDataPtr readoutToT( detSvc(), fmt::format( "{}{}", readoutToTPropBaseName, i ) ); + tip.readoutDelay.push_back( + std::make_shared( readoutDelay, false, gsl_interp_cspline ) ); + tip.readoutTimeOverThreshold.push_back( + std::make_shared( readoutToT, false, gsl_interp_cspline ) ); + } + // check validity of the interpolators + for ( auto const& el : tip.readoutDelay ) { + if ( !el->valid() ) { + throw GaudiException( + "DetailedFrontEndResponsePMT::readTimeInfoPropertiesFromDB", + "Incorrect time info format. Check properties with readout delay/time-over-threshold in the DB.", + StatusCode::FAILURE ); } - - for ( auto const& el : m_readoutTimeOverThreshold ) { - - if ( !el->valid() ) { - error() << "Incorrect time info format. Check properties with readout delay/time-over-threshold in the DB." - << endmsg; - return StatusCode::FAILURE; - } + } + for ( auto const& el : tip.readoutTimeOverThreshold ) { + if ( !el->valid() ) { + throw GaudiException( + "DetailedFrontEndResponsePMT::readTimeInfoPropertiesFromDB", + "Incorrect time info format. Check properties with readout delay/time-over-threshold in the DB.", + StatusCode::FAILURE ); } - - return StatusCode::SUCCESS; - - } else { - - _ri_debug - << "Condition with time info not available in this DB version. Any time information will be initialized to " - "default values ('0' / nullptr)." - << endmsg; - return StatusCode::SUCCESS; } +#endif + return tip; } -float DetailedFrontEndResponsePMT::getReadoutDelay( const float inputSignal, const float threshold ) const { +float Rich::MC::Digi::DetailedFrontEndResponsePMT::getReadoutDelay( TimeInfoProperties const& tip, float inputSignal, + float threshold ) const { - if ( m_readoutDelay.empty() ) { + if ( tip.readoutDelay.empty() ) { return 0.; } else { - const auto readoutTimeInfoIndex = validTimeInfoInputIndex( threshold ); + const auto readoutTimeInfoIndex = validTimeInfoInputIndex( tip, threshold ); if ( msgLevel( MSG::DEBUG ) ) { - if ( !( m_readoutDelay[readoutTimeInfoIndex]->withinInputRange( inputSignal ) ) ) + if ( !( tip.readoutDelay[readoutTimeInfoIndex]->withinInputRange( inputSignal ) ) ) debug() << "No reliable time info for this input signal. Returned value will be from an extrapolation." << endmsg; } - return m_readoutDelay[readoutTimeInfoIndex]->value( inputSignal ); + return tip.readoutDelay[readoutTimeInfoIndex]->value( inputSignal ); } } -float DetailedFrontEndResponsePMT::getReadoutTimeOverThreshold( const float inputSignal, const float threshold ) const { +float Rich::MC::Digi::DetailedFrontEndResponsePMT::getReadoutTimeOverThreshold( TimeInfoProperties const& tip, + float inputSignal, + float threshold ) const { - if ( m_readoutTimeOverThreshold.empty() ) { + if ( tip.readoutTimeOverThreshold.empty() ) { return 0.; } else { - const auto readoutTimeInfoIndex = validTimeInfoInputIndex( threshold ); + const auto readoutTimeInfoIndex = validTimeInfoInputIndex( tip, threshold ); if ( msgLevel( MSG::DEBUG ) ) { - if ( !( m_readoutDelay[readoutTimeInfoIndex]->withinInputRange( inputSignal ) ) ) + if ( !( tip.readoutDelay[readoutTimeInfoIndex]->withinInputRange( inputSignal ) ) ) debug() << "No reliable time info for this input signal. Returned value will be from an extrapolation." << endmsg; } - return m_readoutTimeOverThreshold[readoutTimeInfoIndex]->value( inputSignal ); + return tip.readoutTimeOverThreshold[readoutTimeInfoIndex]->value( inputSignal ); } } // apply a time-gate on the time-sampled readout output -inline bool DetailedFrontEndResponsePMT::PMTChannel::applyTimeGate( const TimeGateLookupTable& timeGateLookupTable ) { +inline bool Rich::MC::Digi::DetailedFrontEndResponsePMT::PMTChannel::applyTimeGate( + const TimeGateLookupTable& timeGateLookupTable ) { auto passesTimeGate = false; @@ -541,42 +850,22 @@ inline bool DetailedFrontEndResponsePMT::PMTChannel::applyTimeGate( const TimeGa //============================================================================= // Global RICH properties //============================================================================= -StatusCode DetailedFrontEndResponsePMT::readPhotoelectronEnergyFromDB() { - +StatusCode Rich::MC::Digi::DetailedFrontEndResponsePMT::readPhotoelectronEnergyFromDB() { +#ifndef USE_DD4HEP SmartDataPtr photoelectronEnergyTab( detSvc(), DeRichLocations::PMTHighVoltageTabPropLoc ); - if ( !photoelectronEnergyTab ) { return Error( "No info on single photoelectron deposit energy found in the DB in the following location: " + DeRichLocations::PMTHighVoltageTabPropLoc ); } else { m_photoelectronEnergy = photoelectronEnergyTab->table().begin()->second; - if ( !( m_photoelectronEnergy > 0. ) ) { - error() << "Invalid value of single photoelectron deposit energy ( must be a positive number ): " - << m_photoelectronEnergy << endmsg; - return StatusCode::FAILURE; - } else { - _ri_debug << "Initialized single photoelectron deposit energy: " << m_photoelectronEnergy << endmsg; - return StatusCode::SUCCESS; - } } -} - -StatusCode DetailedFrontEndResponsePMT::initSinOccScaleFactor() { - - SmartRef genHeader{get( LHCb::GenHeaderLocation::Default )}; - if ( !genHeader ) { - m_sinOccScaleFactor = 1.0f; - - info() << "Gen header not found. Occupancy for the SIN simulation will not be scaled to the current nu. The one " - "provided for a reference nu will be used (nu = " - << std::fixed << std::setprecision( 5 ) << m_sinOccRefNu << ")." << endmsg; - return StatusCode::SUCCESS; +#endif + if ( !( m_photoelectronEnergy > 0. ) ) { + error() << "Invalid value of single photoelectron deposit energy ( must be a positive number ): " + << m_photoelectronEnergy << endmsg; + return StatusCode::FAILURE; } else { - m_sinOccScaleFactor = static_cast( genHeader->beamParameters()->nu() ) / m_sinOccRefNu; - - info() << "Occupancy for the SIN simulation will be scaled by a factor of " << std::fixed << std::setprecision( 4 ) - << m_sinOccScaleFactor << " wrt the reference one provided for nu = " << std::fixed << std::setprecision( 5 ) - << m_sinOccRefNu << "." << endmsg; + _ri_debug << "Initialized single photoelectron deposit energy: " << m_photoelectronEnergy << endmsg; return StatusCode::SUCCESS; } } diff --git a/Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.h b/Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.h deleted file mode 100755 index e6f318253..000000000 --- a/Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.h +++ /dev/null @@ -1,372 +0,0 @@ -/*****************************************************************************\ -* (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. * -\*****************************************************************************/ - -#pragma once - -// STL -#include -#include - -// from Kernel -#include "Kernel/RichDetectorType.h" -#include "Kernel/RichSmartID.h" -#include "RichKernel/RichHistoAlgBase.h" - -// Interfaces -#include "RichInterfaces/IRichSmartIDTool.h" - -// Utils -#include "RichUtils/RichHashMap.h" -#include "RichUtils/RichMap.h" - -// from Gaudi -#include "GaudiKernel/RndmGenerators.h" - -// from Event -#include "Event/MCRichDeposit.h" -#include "Event/MCRichDigit.h" -#include "Event/MCRichHit.h" -#include "Event/MCRichSummedDeposit.h" - -// from RICHDet -#include "RichDet/DeRichPMT.h" -#include "RichDet/DeRichSystem.h" - -// Local -#include "RichChannelPropertiesPMT.h" - -/** @class DetailedFrontEndResponsePMT DetailedFrontEndResponsePMT.h - * - * @author Marcin Kucharczyk, Mariusz Witek - * @date 2015-10-08 - */ - -using namespace LHCb; - -namespace Rich { - namespace MC { - namespace Digi { - - class DetailedFrontEndResponsePMT final : public Rich::HistoAlgBase { - - public: - DetailedFrontEndResponsePMT( const std::string& name, ISvcLocator* pSvcLocator ); - - StatusCode initialize() override final; - StatusCode execute() override final; - StatusCode finalize() override final; - - private: - class PMTChannel; - typedef Rich::Map PMTChannelContainer; - StatusCode Analog( PMTChannelContainer& channels ); - StatusCode Digital( PMTChannelContainer& channels ); - - // get gain based on channel properties (gain mean/rms); require that gain is higher than certain value - // (optionally) - float getChannelGain( const LHCb::RichSmartID& RichSmartID, const float& minValue = 0. ); - // read data needed for time info simulation from the DB - StatusCode readTimeInfoPropertiesFromDB(); - // get time from signal arrival at the anode to it's acquisition at the readout output [ns] - float getReadoutDelay( const float inputSignal, const float threshold ) const; - // get time-over-threshold of the readout output signal after acquisition [ns] - float getReadoutTimeOverThreshold( const float inputSignal, const float threshold ) const; - // read the single photoelectron energy deposit from the DB ( must be the same number as used in Gauss to create - // the deposit ) - StatusCode readPhotoelectronEnergyFromDB(); - - private: - Gaudi::Property m_readParametersFromDB{this, "ReadParametersFromDB", true}; - Gaudi::Property m_sin{this, "Sin", false}; - Gaudi::Property m_spillover{this, "SpilloverModel", "None"}; - Gaudi::Property m_gainMean{this, "GainMean", 0.9}; - Gaudi::Property m_gainRms{this, "GainRms", 0.2}; - Gaudi::Property m_threshold{this, "Threshold", 0.125}; - - // properties for the time gate (all times in ns) - const float m_timeGateWindowRange = 25.0f; - const float m_timeGateSlotWidth = 3.125f; - const unsigned int m_timeGateNrOfSlots = 8; - typedef std::bitset<8> TimeGateSlots; - typedef std::vector TimeGateLookupTable; - - // time gate is applied to 'raw' deposit times with t0 at the collision time (no offsets) - // the beginning of the time-gate window and the lookup table are given separately for R1 / R2 - // if the lookup table is not provided, a 25 ns time window is used by default - Gaudi::Property m_timeGate{this, "TimeGate", true}; - Gaudi::Property> m_timeGateWindowBegin{this, "TimeGateWindowBegin", {18.00, 48.00}}; - Gaudi::Property> m_timeGateLookupTableR1{this, "TimeGateLookupTableR1", {}}; - Gaudi::Property> m_timeGateLookupTableR2{this, "TimeGateLookupTableR2", {}}; - - // end of the time-gate window [ns] (based on the begin and range) - std::vector m_timeGateWindowEnd{0., 0.}; - - // time-gate lookup tables (parsed to a proper format from a string from the Gaudi property) - TimeGateLookupTable m_lookupTableR1{}; - TimeGateLookupTable m_lookupTableR2{}; - - // SIN ratio given in the DB is measured for a certain time window range (m_sinRatioTimeRange); SIN hits are - // added with probability based on this measured ratio in a number of required time slots (separately in each - // time slot); SinNrOfTimeSlots corresponds to the number of time slots before the main time-gate window where - // the SIN should be simulated (including the main time slot) - TimeGate simulation needs to be activated to - // enable this (otherwise SIN will be simulated only for the main time slot) - Gaudi::Property m_sinNrOfTimeSlots{this, "SinNrOfTimeSlots", 3}; - const float m_sinRatioTimeRange{25.0f}; - - // occupancies used for the SIN simulation are given in the DB for a reference nu (m_sinOccRefNu); they are - // scaled accordingly (m_sinOccScaleFactor) to the current (run-time) nu - Gaudi::Property m_sinOccScaleToCurrentNu{this, "SinOccScaleToCurrentNu", true}; - float m_sinOccScaleFactor{0.}; - const float m_sinOccRefNu{7.60022f}; - bool m_sinOccScaleFactorIsInit{false}; - StatusCode initSinOccScaleFactor(); - - Gaudi::Property m_mcRichSumDepsLocation{this, "MCRichSummedDepositsLocation", - LHCb::MCRichSummedDepositLocation::Default}; - Gaudi::Property m_mcRichDigitsLocation{this, "MCRichDigitsLocation", - LHCb::MCRichDigitLocation::Default}; - - // parameters for the monitoring histograms - const float m_histInputSignalMax = 5.0; - const float m_histReadoutDelayMax = 20.0; - const float m_histTimeOverThresholdMax = 35.0; - const float m_histTimeOfAcquisitionMax = 100.0; - - const unsigned int m_nrOfBinsInHistogram1D = 100; - const unsigned int m_nrOfBinsInHistogram2D = 50; - - RichChannelPropertiesPMT* m_propertyTool = nullptr; - DeRichSystem* m_deRichSystem = nullptr; - LHCb::RichSmartID::Vector m_readoutChannelsList{}; - LHCb::MCRichSummedDeposits* m_sDeposits = nullptr; - mutable Rndm::Numbers m_gaussForGain{}; - mutable Rndm::Numbers m_flatForSinProb{}; - mutable Rndm::Numbers m_flatForSinTime{}; - mutable Rndm::Numbers m_pdTransitTimeTypeR{}; - mutable Rndm::Numbers m_pdTransitTimeTypeH{}; - - // energy of a single photoelectron deposit from Gauss [MeV] - float m_photoelectronEnergy{0.}; - - float m_timeInfoIndexIntervalInElectrons{0.}; - unsigned int m_readoutTimeInfoIndexMin{0}; - unsigned int m_readoutTimeInfoIndexMax{0}; - - std::vector> m_readoutDelay{}; - std::vector> m_readoutTimeOverThreshold{}; - - // time of flight of a photoelectron in a PD (from photocathode to anode) in Gauss [ns] - float m_pdPhotoelectronTimeOfFlight{0.}; - // info on PD transit time (mean and spread) for R/H PD types [ns] - float m_pdTransitTimeMeanTypeR{0.}; - float m_pdTransitTimeSpreadTypeR{0.}; - float m_pdTransitTimeMeanTypeH{0.}; - float m_pdTransitTimeSpreadTypeH{0.}; - - // structure for deposits (not the same as RichDeposit, it is what is effectively seen by the Claro - can be a - // SIN hit, or two deposits close in time which are seen as a single one) - struct PMTDeposit { - float inputSignal{0.}; - float timeOfAcquisition{0.}; - float timeOverThreshold{0.}; - }; - - // helper class to process all deposits within a single PMT channel (including SIN) - class PMTChannel { - - private: - typedef std::vector PMTDepositContainer; - PMTDepositContainer m_deposits{}; - TimeGateSlots m_outputTimeSampling{}; - bool m_hasSin{false}; - LHCb::MCRichSummedDeposit* m_summedDeposit{nullptr}; - typedef std::vector TimeGatePatternMatches; - TimeGatePatternMatches m_timeGatePatternMatches{}; - - public: - inline void addDeposit( const float inputSignal, // - const float timeOfAcquisition, // - const float timeOverThreshold ) { - PMTDeposit deposit; - deposit.inputSignal = inputSignal; - deposit.timeOfAcquisition = timeOfAcquisition; - deposit.timeOverThreshold = timeOverThreshold; - m_deposits.push_back( deposit ); - } - - // add the given deposit (with a certain time-of-acquisition and time-over-threshold) to the Claro output - // pattern - inline void processDeposit( const float timeOfAcquisition, // - const float timeOverThreshold, // - const float timeGateMin, // - const float timeGateMax, // - const float timeGateBin ) { - - // check if the deposit has any overlap with the time-sampling window - if ( isOverlapBetweenSignalAndTimeSampling( timeOfAcquisition, timeOverThreshold, timeGateMin, - timeGateMax ) ) { - - const auto firstTimeSlotIndex = getTimeSlotIndex( timeOfAcquisition, timeGateMin, timeGateBin ); - - const auto lastTimeSlotIndex = - getTimeSlotIndex( timeOfAcquisition + timeOverThreshold, timeGateMin, timeGateBin ); - - // set time slots status within the time-sampling window; time slots ordering is reversed wrt to the - // 'bitset' format - for ( auto i = firstTimeSlotIndex; i < lastTimeSlotIndex + 1; ++i ) { - if ( i >= 0 && (std::size_t)i < m_outputTimeSampling.size() ) { - m_outputTimeSampling.set( m_outputTimeSampling.size() - 1 - i ); - } - } - } - } - - // apply a time-gate on the Claro output - inline bool applyTimeGate( const TimeGateLookupTable& timeGateLookupTable ); - - inline const TimeGateSlots& getOutputTimeSampling() const { return m_outputTimeSampling; } - inline const PMTDepositContainer& getDeposits() const { return m_deposits; } - - inline LHCb::MCRichSummedDeposit* getSummedDeposit() { return m_summedDeposit; } - inline void setSummedDeposit( LHCb::MCRichSummedDeposit* summedDeposit ) { m_summedDeposit = summedDeposit; } - - inline bool hasSin() const noexcept { return m_hasSin; } - inline void setHasSin() { m_hasSin = true; } - - inline const TimeGatePatternMatches& getTimeGatePatternMatches() const { return m_timeGatePatternMatches; } - - private: - // check if the given signal has any overlap (in time) with the time-sampling window - inline bool isOverlapBetweenSignalAndTimeSampling( const float timeOfAcquisition, // - const float timeOverThreshold, // - const float timeGateMin, // - const float timeGateMax ) const { - return ( ( timeOfAcquisition <= timeGateMax ) && ( timeOfAcquisition + timeOverThreshold >= timeGateMin ) ); - } - - // get index of the time slot to which the given time corresponds - inline int getTimeSlotIndex( const float time, const float timeGateMin, const float timeGateBin ) const { - const auto var = ( 0 != timeGateBin ? ( time - timeGateMin ) / timeGateBin : 0 ); - // mimics std::floor, but avoid uses it as throws strange FPE with clang.. - auto j = static_cast( var ); - if ( var < 0 ) { --j; } - return j; - } - - // process a single line of the lookup table - inline bool processLookupTableLine( const TimeGateSlots& line, - const TimeGateSlots& outputTimeSampling ) const { - return ( outputTimeSampling == line ); - } - }; - - // helper functions - inline const DeRichPMT* getDePmtFromSmartId( const LHCb::RichSmartID smartID ) const { - return (DeRichPMT*)m_deRichSystem->dePD( smartID ); - } - inline float getChannelGainMean( const LHCb::RichSmartID smartID ) const { - return m_readParametersFromDB.value() ? getDePmtFromSmartId( smartID )->PmtChannelGainMean( smartID ) - : m_gainMean.value(); - } - inline float getChannelGainRms( const LHCb::RichSmartID smartID ) const { - return m_readParametersFromDB.value() ? getDePmtFromSmartId( smartID )->PmtChannelGainRms( smartID ) - : m_gainRms.value(); - } - inline float getChannelThreshold( const LHCb::RichSmartID smartID ) const { - return m_readParametersFromDB.value() ? getDePmtFromSmartId( smartID )->PmtChannelThreshold( smartID ) - : m_threshold.value(); - } - inline float getChannelSinRatio( const LHCb::RichSmartID smartID ) const { - return m_readParametersFromDB.value() ? getDePmtFromSmartId( smartID )->PmtChannelSinRatio( smartID ) : 0.; - } - inline float getPmtAverageOccupancy( const LHCb::RichSmartID smartID ) const { - return m_readParametersFromDB.value() - ? ( m_sinOccScaleToCurrentNu - ? getDePmtFromSmartId( smartID )->PmtAverageOccupancy() * m_sinOccScaleFactor - : getDePmtFromSmartId( smartID )->PmtAverageOccupancy() ) - : 0.; - } - - inline void updateHistoryWithSin( LHCb::MCRichDigit*& digit ) const { - auto digitHistory = digit->history(); - digitHistory.setSignalInducedNoise( true ); - digit->setHistory( digitHistory ); - } - - // provide a valid index for time info input set (take the closest one from the DB) - inline unsigned int validTimeInfoInputIndex( const float thresholdInElectrons ) const { - - auto inputIndex = round( thresholdInElectrons / m_timeInfoIndexIntervalInElectrons ); - - if ( inputIndex < m_readoutTimeInfoIndexMin ) { - warning() << "No reliable time info for this threshold setting. The closest time info input will be used." - << endmsg; - inputIndex = m_readoutTimeInfoIndexMin; - } else if ( inputIndex > m_readoutTimeInfoIndexMax ) { - warning() << "No reliable time info for this threshold setting. The closest time info input will be used." - << endmsg; - inputIndex = m_readoutTimeInfoIndexMax; - } - - // offset the index ( should start from 0 ) - inputIndex -= m_readoutTimeInfoIndexMin; - - return inputIndex; - } - - // provide time of acquisition at the Claro output (taking all offsets into account) - inline float timeOfAcquisition( const float timeOfDeposit, // - const float inputSignal, // - const float threshold, // - const bool isLargePMT ) const { - const auto pdTransitTime = isLargePMT ? m_pdTransitTimeTypeH() : m_pdTransitTimeTypeR(); - return ( timeOfDeposit - m_pdPhotoelectronTimeOfFlight + pdTransitTime + - getReadoutDelay( inputSignal, threshold ) ); - } - - // provide time of acquisition at the readout output for SIN hits ( random time within the time window of - // interest ) - inline float timeOfAcquisitionForSin( const float timeWindowMin ) const { - return timeWindowMin + m_flatForSinTime(); - } - - // return a proper number if given timeGate property is different for R1/R2 - inline float timeGateParam( const LHCb::RichSmartID smartID, std::vector& timeGateParamProperty ) const { - const auto rich = smartID.rich(); - if ( rich == Rich::Rich1 ) { - return timeGateParamProperty[0]; - } else if ( rich == Rich::Rich2 ) { - return timeGateParamProperty[1]; - } else { - warning() << "Given smartID not in R1/R2." << endmsg; - return 0.; - } - } - - // parse & validate the patterns in a lookup table - inline bool parseLookupTable( const std::vector& lookupTableFrom, - TimeGateLookupTable& lookupTableTo ) const { - bool isValid = true; - for ( auto const& line : lookupTableFrom ) { - if ( line.length() != m_timeGateNrOfSlots ) { - isValid = false; - break; - } else { - lookupTableTo.emplace_back( line ); - } - } - return isValid; - } - }; - - } // namespace Digi - } // namespace MC -} // namespace Rich -- GitLab From d9f974923dd150f1876628615e668893ea6b8311 Mon Sep 17 00:00:00 2001 From: Sebastien Ponce Date: Wed, 19 Oct 2022 11:13:51 +0200 Subject: [PATCH 13/15] Do not use RichDetailedFrontEndResponsePMT in DD4hep mode as it's not DD4hep compliant at this stage --- CMakeLists.txt | 11 ++++++++++- Digi/Boole/CMakeLists.txt | 9 ++++++--- Rich/RichDigiSys/python/RichDigiSys/Configuration.py | 3 ++- Rich/RichReadout/CMakeLists.txt | 8 +++++++- .../{ => legacy}/RichDetailedFrontEndResponsePMT.cpp | 2 +- 5 files changed, 26 insertions(+), 7 deletions(-) rename Rich/RichReadout/src/component/{ => legacy}/RichDetailedFrontEndResponsePMT.cpp (99%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 046e90fbc..4f3e3bc2f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,8 +25,18 @@ list(PREPEND CMAKE_MODULE_PATH set(WITH_Boole_PRIVATE_DEPENDENCIES TRUE) include(BooleDependencies) +if (NOT USE_DD4HEP) + set(VarName DetDescOnlySubDirs) +else() + set(VarName LHCB_IGNORE_SUBDIRS) +endif() +set(${VarName} + Rich/RichDigiSys +) + # -- Subdirectories lhcb_add_subdirectories( + ${DetDescOnlySubDirs} Bcm/BcmDigi Calo/CaloDigit Calo/CaloMoniDigi @@ -38,7 +48,6 @@ lhcb_add_subdirectories( Muon/MuonBackground Muon/MuonMoniDigi Rich/RichDigiQC - Rich/RichDigiSys Rich/RichMCAssociators Rich/RichReadout UT/UTDigiAlgorithms diff --git a/Digi/Boole/CMakeLists.txt b/Digi/Boole/CMakeLists.txt index 96e5b6e23..bfb7707d6 100644 --- a/Digi/Boole/CMakeLists.txt +++ b/Digi/Boole/CMakeLists.txt @@ -14,10 +14,13 @@ Digi/Boole #]=======================================================================] gaudi_install(PYTHON) -gaudi_generate_confuserdb() -lhcb_add_confuser_dependencies( + +if (NOT USE_DD4HEP) + gaudi_generate_confuserdb() + lhcb_add_confuser_dependencies( Rich/RichDigiSys -) + ) +endif() lhcb_env(SET BOOLEOPTS ${CMAKE_CURRENT_SOURCE_DIR}/options) diff --git a/Rich/RichDigiSys/python/RichDigiSys/Configuration.py b/Rich/RichDigiSys/python/RichDigiSys/Configuration.py index 386ee0ab0..17e8367b9 100755 --- a/Rich/RichDigiSys/python/RichDigiSys/Configuration.py +++ b/Rich/RichDigiSys/python/RichDigiSys/Configuration.py @@ -24,6 +24,7 @@ from Configurables import ( Rich__MC__Digi__DetailedFrontEndResponsePMT, Rich__MC__Digi__SimpleFrontEndResponse, Rich__MC__Digi__MCRichDigitsToRawBufferAlg, Rich__MC__Digi__SummedDeposits) +from DDDB.CheckDD4Hep import UseDD4Hep # ---------------------------------------------------------------------------------- @@ -127,7 +128,7 @@ class RichDigiSysConf(LHCbConfigurableUser): elif pdModel == "Copy": response = self.makeComponent( Rich__MC__Digi__CopySummedDepositsToDigits, "RichPDResponse") - elif pdModel == "DetailedPMT": + elif pdModel == "DetailedPMT" and not UseDD4Hep: response = self.makeComponent( Rich__MC__Digi__DetailedFrontEndResponsePMT, "RichPDResponse") response.Threshold = self.getProp("Threshold") diff --git a/Rich/RichReadout/CMakeLists.txt b/Rich/RichReadout/CMakeLists.txt index 2288b3fac..526b0261a 100644 --- a/Rich/RichReadout/CMakeLists.txt +++ b/Rich/RichReadout/CMakeLists.txt @@ -20,14 +20,20 @@ gaudi_add_header_only_library(RichReadoutLib LHCb::MCEvent ) +if (NOT USE_DD4HEP) + set(DetDescOnlyFiles + src/component/legacy/RichDetailedFrontEndResponsePMT.cpp + ) +endif() + gaudi_add_module(RichReadout SOURCES + ${DetDescOnlyFiles} src/component/MCRichDigitsToRawBufferAlg.cpp src/component/RichBase.cpp src/component/RichChannelPropertiesPMT.cpp src/component/RichCopySummedDepositsToDigits.cpp src/component/RichDetailedFrontEndResponse.cpp - src/component/RichDetailedFrontEndResponsePMT.cpp src/component/RichPixelReadout.cpp src/component/RichRegistry.cpp src/component/RichSignal.cpp diff --git a/Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.cpp b/Rich/RichReadout/src/component/legacy/RichDetailedFrontEndResponsePMT.cpp similarity index 99% rename from Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.cpp rename to Rich/RichReadout/src/component/legacy/RichDetailedFrontEndResponsePMT.cpp index f68210ed5..df93bd297 100755 --- a/Rich/RichReadout/src/component/RichDetailedFrontEndResponsePMT.cpp +++ b/Rich/RichReadout/src/component/legacy/RichDetailedFrontEndResponsePMT.cpp @@ -9,7 +9,7 @@ * or submit itself to any jurisdiction. * \*****************************************************************************/ -#include "RichChannelPropertiesPMT.h" +#include "../RichChannelPropertiesPMT.h" #include "Kernel/RichDetectorType.h" #include "Kernel/RichSmartID.h" -- GitLab From 9f61b61542bc4de2509385919125e83b8ffad875 Mon Sep 17 00:00:00 2001 From: Christopher Rob Jones Date: Mon, 24 Oct 2022 11:16:50 +0200 Subject: [PATCH 14/15] Default initialise various data members --- FT/FTDigitisation/src/MCFTAttenuationTool.cpp | 4 ++-- FT/FTDigitisation/src/MCFTG4AttenuationTool.cpp | 2 +- FT/FTDigitisation/src/SiPMResponse.h | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/FT/FTDigitisation/src/MCFTAttenuationTool.cpp b/FT/FTDigitisation/src/MCFTAttenuationTool.cpp index 2b5fdb8b7..e7bd4b76f 100644 --- a/FT/FTDigitisation/src/MCFTAttenuationTool.cpp +++ b/FT/FTDigitisation/src/MCFTAttenuationTool.cpp @@ -90,8 +90,8 @@ private: Gaudi::Property m_yMax{this, "YMax", 2426.0 * Gaudi::Units::mm, "Length of the FullAttenuationMap(in mm) in Y-axis"}; - int m_nXSteps; - int m_nYSteps; + int m_nXSteps{0}; + int m_nYSteps{0}; std::vector m_transmissionMap; ///< Maps hits to transmitted energy from the direct pulse std::vector m_transmissionRefMap; ///< Maps hits to transmitted energy from the reflected pulse diff --git a/FT/FTDigitisation/src/MCFTG4AttenuationTool.cpp b/FT/FTDigitisation/src/MCFTG4AttenuationTool.cpp index 0f3847f4d..8da98286a 100644 --- a/FT/FTDigitisation/src/MCFTG4AttenuationTool.cpp +++ b/FT/FTDigitisation/src/MCFTG4AttenuationTool.cpp @@ -21,7 +21,7 @@ namespace { struct Conditions { - unsigned int nBinsX, nBinsY; + unsigned int nBinsX{0}, nBinsY{0}; std::vector xEdges, yEdges, effDir, effRef; }; diff --git a/FT/FTDigitisation/src/SiPMResponse.h b/FT/FTDigitisation/src/SiPMResponse.h index f90592d1f..31c69a4de 100644 --- a/FT/FTDigitisation/src/SiPMResponse.h +++ b/FT/FTDigitisation/src/SiPMResponse.h @@ -37,10 +37,10 @@ public: private: struct Conditions { GaudiMath::SimpleSpline responseSpline; - double tMin; - double tMax; + double tMin{0}; + double tMax{0}; /// calibration of the gain - double sipmGainShift; + double sipmGainShift{0}; }; ConditionAccessor m_conditions{this, name() + "_Conditions"}; -- GitLab From 733b952d9b37bb08e7921617586f4f9647730d9e Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Mon, 24 Oct 2022 10:14:51 +0100 Subject: [PATCH 15/15] Work around floating point exception in RichDetailedFrontEndResponsePMT --- Rich/RichReadout/CMakeLists.txt | 7 +++++++ .../component/legacy/RichDetailedFrontEndResponsePMT.cpp | 6 +----- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/Rich/RichReadout/CMakeLists.txt b/Rich/RichReadout/CMakeLists.txt index 526b0261a..9dfd3dd6e 100644 --- a/Rich/RichReadout/CMakeLists.txt +++ b/Rich/RichReadout/CMakeLists.txt @@ -24,6 +24,13 @@ if (NOT USE_DD4HEP) set(DetDescOnlyFiles src/component/legacy/RichDetailedFrontEndResponsePMT.cpp ) + # This is working around floating point exceptions when compiling RichDetailedFrontEndResponsePMT + # with clang in sse mode. In that mode, the 2 calls to getTimeSlotIndex in processDeposit are + # optimized into a single SSE division. Now the SSE vector register holds 4 floats, and the + # remaining slots are filled with 0s for both operands, leading to a 'invalid' FPE. + if (CMAKE_CXX_COMPILER_ID STREQUAL "Clang" ) + set_source_files_properties( src/component/legacy/RichDetailedFrontEndResponsePMT.cpp PROPERTIES COMPILE_OPTIONS -fno-slp-vectorize ) + endif() endif() gaudi_add_module(RichReadout diff --git a/Rich/RichReadout/src/component/legacy/RichDetailedFrontEndResponsePMT.cpp b/Rich/RichReadout/src/component/legacy/RichDetailedFrontEndResponsePMT.cpp index df93bd297..d46c077cb 100755 --- a/Rich/RichReadout/src/component/legacy/RichDetailedFrontEndResponsePMT.cpp +++ b/Rich/RichReadout/src/component/legacy/RichDetailedFrontEndResponsePMT.cpp @@ -301,11 +301,7 @@ namespace Rich::MC::Digi { // get index of the time slot to which the given time corresponds inline int getTimeSlotIndex( const float time, const float timeGateMin, const float timeGateBin ) const { - const auto var = ( 0 != timeGateBin ? ( time - timeGateMin ) / timeGateBin : 0 ); - // mimics std::floor, but avoid uses it as throws strange FPE with clang.. - auto j = static_cast( var ); - if ( var < 0 ) { --j; } - return j; + return std::floor( 0 != timeGateBin ? ( time - timeGateMin ) / timeGateBin : 0 ); } // process a single line of the lookup table -- GitLab