diff --git a/UT/UTAlgorithms/CMakeLists.txt b/UT/UTAlgorithms/CMakeLists.txt index dc2c3fede73c6a57e1b18aa014b830e8ddc373b8..eb3ffe933a78aa1ee8df97a650e2dfcd3428db7d 100644 --- a/UT/UTAlgorithms/CMakeLists.txt +++ b/UT/UTAlgorithms/CMakeLists.txt @@ -16,7 +16,6 @@ UT/UTAlgorithms gaudi_add_module(UTAlgorithms SOURCES src/UTClusterContainerCopy.cpp - src/UTClusterCreator_Monitor.cpp src/UTClusterKiller.cpp src/UTClustersToDigits.cpp src/UTClustersToLite.cpp diff --git a/UT/UTAlgorithms/src/UTClusterCreator_Monitor.cpp b/UT/UTAlgorithms/src/UTClusterCreator_Monitor.cpp deleted file mode 100644 index 38ab55ebd2c19f31fe33ec6396fbce092e387bd3..0000000000000000000000000000000000000000 --- a/UT/UTAlgorithms/src/UTClusterCreator_Monitor.cpp +++ /dev/null @@ -1,146 +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. * -\*****************************************************************************/ -#include "DetDesc/GenericConditionAccessorHolder.h" -#include "Detector/UT/ChannelID.h" -#include "Event/UTCluster.h" -#include "Event/UTDigit.h" -#include "Kernel/IUTReadoutTool.h" -#include "Kernel/LHCbConstants.h" -#include "LHCbMath/LHCbMath.h" -#include "TrackInterfaces/IUTClusterPosition.h" -#include "UTDet/DeUTDetector.h" -#include "UTDet/DeUTSector.h" - -#include "LHCbAlgs/Transformer.h" - -/** - * Class for clustering in the UT tracker - * - * @author A.Beiter (based on code by M.Needham) - * @date 2018-09-04 - */ - -class UTClusterCreator_Monitor - : public LHCb::Algorithm::Transformer> { - -public: - // Constructors and destructor - UTClusterCreator_Monitor( const std::string& name, ISvcLocator* pSvcLocator ); - - // IAlgorithm members - StatusCode initialize() override; - LHCb::UTClusters operator()( const LHCb::UTDigits& digitCont, const DeUTDetector& det, - const DeMagnet& magnet ) const override; - -private: - StatusCode createClusters( const DeUTDetector&, const DeMagnet&, const LHCb::UTDigits& digitsCont, - LHCb::UTClusters& clustersCont ) const; - - LHCb::UTCluster::ADCVector strips( const SmartRefVector& clusteredDigits, - const LHCb::Detector::UT::ChannelID closestChan ) const; - - Gaudi::Property m_outputVersion{this, "OutputVersion", 1}; - Gaudi::Property m_maxSize{this, "Size", 4}; - Gaudi::Property m_byBeetle{this, "ByBeetle", true}; - Gaudi::Property m_spillName{this, "Spill", "Central"}; - Gaudi::Property m_conditionLocation{this, "conditionLocation", - "/dd/Conditions/ReadoutConf/UT/ClusteringThresholds"}; - - LHCb::UTCluster::Spill m_spill; - - PublicToolHandle m_positionTool{this, "PositionTool", "UTOnlinePosition"}; - ToolHandle m_readoutTool{this, "ReadoutTool", "UTReadoutTool"}; -}; - -using namespace LHCb; - -DECLARE_COMPONENT( UTClusterCreator_Monitor ) - -UTClusterCreator_Monitor::UTClusterCreator_Monitor( const std::string& name, ISvcLocator* pSvcLocator ) - : Transformer{name, - pSvcLocator, - {{"InputLocation", UTDigitLocation::UTTightDigits}, - {"UTLocation", DeUTDetLocation::location()}, - {"DeMagnetLocation", LHCb::Det::Magnet::det_path}}, - {"OutputLocation", UTClusterLocation::UTClusters}} {} - -StatusCode UTClusterCreator_Monitor::initialize() { - return Transformer::initialize().andThen( [&] { - m_spill = UTCluster::SpillToType( m_spillName ); - return StatusCode::SUCCESS; - } ); -} - -LHCb::UTClusters UTClusterCreator_Monitor::operator()( const UTDigits& digitCont, const DeUTDetector& det, - const DeMagnet& magnet ) const { - UTClusters clusterCont; - clusterCont.setVersion( m_outputVersion ); - createClusters( det, magnet, digitCont, clusterCont ) - .orThrow( "Failed to create clusters", "UTClusterCreator_Monitor" ) - .ignore(); - return clusterCont; -} - -StatusCode UTClusterCreator_Monitor::createClusters( const DeUTDetector& det, const DeMagnet& magnet, - const UTDigits& digitCont, UTClusters& clusterCont ) const { - - auto iterDigit = digitCont.begin(); - while ( iterDigit != digitCont.end() ) { - - auto jterDigit = iterDigit; - ++jterDigit; - - // make a cluster ! - SmartRefVector clusteredDigits; - clusteredDigits.push_back( *iterDigit ); - double totCharge = ( *iterDigit )->depositedCharge(); - - // calculate the interstrip fraction a la Tell1 - IUTClusterPosition::Info measValue = m_positionTool->estimate( det, magnet, clusteredDigits ); - - const double nSum = totCharge; - - // make cluster +set things - const UTLiteCluster clusterLite( measValue.strip, measValue.fractionalPosition, clusteredDigits.size(), true, - true ); - - // link to the DAQ - UTDAQID utdaqID = m_readoutTool->channelIDToDAQID( static_cast( measValue.strip ) ); - if ( utdaqID.id() == UTTell1ID::nullBoard ) { - // screwed - this->err() << "failed to link " << LHCb::UTNames::UniqueSectorToString( measValue.strip ) << endmsg - << "Failed to find Tell1 board" << endmsg; - return StatusCode::FAILURE; - } - - UTCluster* newCluster = new UTCluster( clusterLite, strips( clusteredDigits, measValue.strip ), nSum, utdaqID.id(), - utdaqID.channel(), m_spill ); - // add to container - clusterCont.insert( newCluster, measValue.strip ); - - iterDigit = jterDigit; - } // iterDigit - - return StatusCode::SUCCESS; -} - -UTCluster::ADCVector UTClusterCreator_Monitor::strips( const SmartRefVector& clusteredDigits, - const Detector::UT::ChannelID closestChan ) const { - // make the vector of ADC values stored in the cluster - UTCluster::ADCVector tVec; - for ( const auto& d : clusteredDigits ) { - const int strip = d->channelID().strip() - closestChan.strip(); - tVec.emplace_back( strip, d->depositedCharge() ); - } - return tVec; -} diff --git a/UT/UTAssociators/CMakeLists.txt b/UT/UTAssociators/CMakeLists.txt index 6d81b2e4e32b8a1dab6b4e3f3ab5e536ec17277e..a812808e8cde626547911f67f964f4a5451009ee 100644 --- a/UT/UTAssociators/CMakeLists.txt +++ b/UT/UTAssociators/CMakeLists.txt @@ -20,6 +20,8 @@ gaudi_add_module(UTAssociators src/UTDigit2MCHitLinker.cpp src/UTDigit2MCParticleLinker.cpp src/UTTruthTool.cpp + src/UTTightDigit2MCHitLinker.cpp + src/UTTightDigit2MCParticleLinker.cpp LINK Gaudi::GaudiKernel LHCb::DigiEvent diff --git a/UT/UTAssociators/src/UTDigit2MCHitLinker.cpp b/UT/UTAssociators/src/UTDigit2MCHitLinker.cpp index 678d113f88924e7dc88d9a6c35c5886c9288e6e9..e37ba3b2f24dd91835b484405401eb27ea67a138 100644 --- a/UT/UTAssociators/src/UTDigit2MCHitLinker.cpp +++ b/UT/UTAssociators/src/UTDigit2MCHitLinker.cpp @@ -12,6 +12,7 @@ #include "Event/LinksByKey.h" #include "Event/MCHit.h" #include "Event/MCTruth.h" +#include "Event/MCUTDigit.h" #include "Event/UTDigit.h" #include "LHCbAlgs/Transformer.h" #include "UTTruthTool.h" diff --git a/UT/UTAssociators/src/UTTightDigit2MCHitLinker.cpp b/UT/UTAssociators/src/UTTightDigit2MCHitLinker.cpp new file mode 100644 index 0000000000000000000000000000000000000000..03a6e31ff160657874d9702dce9ee9d0e42aecba --- /dev/null +++ b/UT/UTAssociators/src/UTTightDigit2MCHitLinker.cpp @@ -0,0 +1,112 @@ +/*****************************************************************************\ +* (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. * +\*****************************************************************************/ + +#include "Event/LinksByKey.h" +#include "Event/MCHit.h" +#include "Event/UTCluster.h" +#include "Event/UTDigit.h" +#include "LHCbAlgs/Transformer.h" +#include "Linker/LinkedTo.h" + +using namespace LHCb; + +/** + * @author Andy Beiter (based on code by Matt Needham) + * @date 2018-09-04 + */ + +class UTTightDigit2MCHitLinker : public LHCb::Algorithm::Transformer { +public: + UTTightDigit2MCHitLinker( std::string const& name, ISvcLocator* pSvc ) + : Transformer{name, + pSvc, + {{"InputData", UTDigitLocation::UTTightDigits}, + {"DigitLocation", UTDigitLocation::UTDigits}, + {"hitLocation", "/Event/" + MCHitLocation::UT}, + {"DigitASCTlocation", UTDigitLocation::UTDigits + "2MCHits"}}, + {"OutputData", LHCb::LinksByKey::linkerName( UTClusterLocation::UTClusters + "2MCHits" )}} {} + + LHCb::LinksByKey operator()( const UTDigits&, const UTDigits&, const MCHits&, + const LHCb::LinksByKey& ) const override; + +private: + typedef std::pair HitPair; + typedef std::map HitMap; + std::vector refsToRelate( const HitMap& hitMap, const LHCb::MCHits& hits ) const; + HitMap associateToTruth( const LHCb::UTDigits&, const LHCb::UTDigit&, const LinkedTo& ) const; + + Gaudi::Property m_addSpillOverHits{this, "AddSpillOverHits", false, "Flag to add spill-over to linker table"}; + Gaudi::Property m_minFrac{this, "Minfrac", 0.2, "Minimal charge fraction to link to MCParticle"}; + Gaudi::Property m_oneRef{this, "OneRef", false, "Flag to allow only 1 link for each cluster"}; +}; + +DECLARE_COMPONENT( UTTightDigit2MCHitLinker ) + +LHCb::LinksByKey UTTightDigit2MCHitLinker::operator()( const UTDigits& clusterCont, const UTDigits& digitCont, + const MCHits& mcHits, const LHCb::LinksByKey& links ) const { + // create a linker + auto myLink = LHCb::LinksByKey{std::in_place_type, std::in_place_type, + LHCb::LinksByKey::Order::decreasingWeight}; + + // make link to truth to MCHit from cluster + // TODO: can we really delete this? This algo does require MCHits as input, so they must exist... : + // links->resolveLinks( evtSvc() ); + const auto aTable = LinkedTo{&links}; + + // loop and link UTClusters to MC truth + for ( auto iterClus : clusterCont ) { + // find all hits + // select references to add to table + auto selectedRefs = refsToRelate( associateToTruth( digitCont, *iterClus, aTable ), mcHits ); + if ( !selectedRefs.empty() ) { + if ( !m_oneRef.value() ) { + myLink.link( iterClus, selectedRefs ); + } else { + auto [p, w] = selectedRefs.back(); + myLink.link( iterClus, p, w ); + } + } // refsToRelate ! empty + } // loop iterClus + return myLink; +} + +std::vector UTTightDigit2MCHitLinker::refsToRelate( const HitMap& hitMap, + const MCHits& hits ) const { + std::vector selRefs; + for ( auto [aHit, weight] : hitMap ) { + if ( aHit && weight > m_minFrac ) { + if ( m_addSpillOverHits.value() || &hits == aHit->parent() ) { selRefs.emplace_back( aHit, weight ); } + } + } + return selRefs; +} + +UTTightDigit2MCHitLinker::HitMap +UTTightDigit2MCHitLinker::associateToTruth( const LHCb::UTDigits& digitCont, const UTDigit& aCluster, + const LinkedTo& aTable ) const { + HitMap hitMap; + auto id = aCluster.channelID(); + const UTDigit* aDigit = digitCont.object( id ); + if ( aDigit ) { + double foundCharge = 0.; + for ( const auto& [aHit, weight] : aTable.weightedRange( aDigit ) ) { + hitMap[&aHit] += weight * aDigit->depositedCharge(); + foundCharge += weight * aDigit->depositedCharge(); + } + hitMap[nullptr] += aDigit->depositedCharge() - foundCharge; + } // Digit + + // renormalize to 1 + for ( auto& [_, w] : hitMap ) w /= aCluster.depositedCharge(); + + return hitMap; +} diff --git a/UT/UTAssociators/src/UTTightDigit2MCParticleLinker.cpp b/UT/UTAssociators/src/UTTightDigit2MCParticleLinker.cpp new file mode 100644 index 0000000000000000000000000000000000000000..dd9bd9f819f805ab4ea79717f876e3ad8f0b8643 --- /dev/null +++ b/UT/UTAssociators/src/UTTightDigit2MCParticleLinker.cpp @@ -0,0 +1,108 @@ +/*****************************************************************************\ +* (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. * +\*****************************************************************************/ + +#include "Event/MCHit.h" +#include "Event/MCParticle.h" +#include "Event/UTCluster.h" +#include "Event/UTDigit.h" +#include "Linker/LinkedTo.h" + +#include "LHCbAlgs/Transformer.h" + +/** + * @author Andy Beiter (based on code by Matt Needham) + * @date 2018-09-04 + */ +class UTTightDigit2MCParticleLinker : public LHCb::Algorithm::Transformer { +public: + UTTightDigit2MCParticleLinker( const std::string& name, ISvcLocator* pSvc ) + : Transformer{ + name, + pSvc, + {{"InputData", LHCb::UTDigitLocation::UTTightDigits}, + {"MCParticleLocation", LHCb::MCParticleLocation::Default}, + {"ClusterASCTlocation", LHCb::LinksByKey::linkerName( LHCb::UTClusterLocation::UTClusters + "2MCHits" )}}, + {"OutputData", LHCb::LinksByKey::linkerName( LHCb::UTClusterLocation::UTClusters )}} {} + + LHCb::LinksByKey operator()( const LHCb::UTDigits&, const LHCb::MCParticles&, + const LHCb::LinksByKey& ) const override; + +private: + using PartPair = std::pair; + using ParticleMap = std::map; + std::vector refsToRelate( const ParticleMap&, LHCb::MCParticles const& ) const; + + ParticleMap associateToTruth( const LHCb::UTDigit*, const LinkedTo& ) const; + + Gaudi::Property m_addSpillOverHits{this, "AddSpillOverHits", false, "Flag to add spill-over to linker table"}; + Gaudi::Property m_minFrac{this, "Minfrac", 0.2, "Minimal charge fraction to link to MCParticle"}; + Gaudi::Property m_oneRef{this, "OneRef", false, "Flag to allow only 1 link for each cluster"}; +}; + +DECLARE_COMPONENT( UTTightDigit2MCParticleLinker ) + +LHCb::LinksByKey UTTightDigit2MCParticleLinker:: + operator()( const LHCb::UTDigits& clusterCont, const LHCb::MCParticles& mcParts, const LHCb::LinksByKey& links ) const { + LinkedTo aTable{&links}; + // create a linker + auto aLinker = LHCb::LinksByKey{std::in_place_type, std::in_place_type, + LHCb::LinksByKey::Order::decreasingWeight}; + aLinker.reset(); + + // loop and link UTClusters to MC truth + for ( auto iterClus : clusterCont ) { + // find all particles + ParticleMap partMap = associateToTruth( iterClus, aTable ); + + // select references to add to table + auto selectedRefs = refsToRelate( partMap, mcParts ); + if ( !selectedRefs.empty() ) { + if ( !m_oneRef.value() ) { + aLinker.link( iterClus, selectedRefs ); + } else { + auto [p, w] = selectedRefs.back(); + aLinker.link( iterClus, p, w ); + } + } // refsToRelate != empty + } // loop iterClus + + return aLinker; +} + +std::vector +UTTightDigit2MCParticleLinker::refsToRelate( const ParticleMap& partMap, LHCb::MCParticles const& particles ) const { + std::vector selRefs; + for ( auto const& [aParticle, w] : partMap ) { + if ( aParticle && w > m_minFrac && ( m_addSpillOverHits.value() || &particles == aParticle->parent() ) ) { + selRefs.emplace_back( aParticle, w ); + } + } + return selRefs; +} + +UTTightDigit2MCParticleLinker::ParticleMap +UTTightDigit2MCParticleLinker::associateToTruth( const LHCb::UTDigit* aCluster, + const LinkedTo& aTable ) const { + + ParticleMap partMap; + double foundCharge = 0; + for ( const auto& [hit, weight] : aTable.weightedRange( aCluster ) ) { + const LHCb::MCParticle* aParticle = hit.mcParticle(); + partMap[aParticle] += weight; + foundCharge += weight; + } + + // difference between depEnergy and total cluster charge = noise (due to norm) + partMap[nullptr] += 1.0 - foundCharge; + + return partMap; +} diff --git a/UT/UTCheckers/src/UTEffChecker.cpp b/UT/UTCheckers/src/UTEffChecker.cpp index 6dad8eca4c35b558533b144db4f35a8eeedd85a7..bbefea7dd40bd9e3077ea248dea44e5ffd17db80 100644 --- a/UT/UTCheckers/src/UTEffChecker.cpp +++ b/UT/UTCheckers/src/UTEffChecker.cpp @@ -13,6 +13,8 @@ #include "Event/MCHit.h" #include "Event/MCParticle.h" #include "Event/UTCluster.h" +#include "Event/UTDigit.h" +#include "Kernel/UTIDMapping.h" #include "Kernel/UTLexicalCaster.h" #include "Linker/LinkedFrom.h" #include "Linker/LinkedTo.h" @@ -49,10 +51,13 @@ namespace { * * In the finalize method, a summary of the efficiency per layer is printed. * - * @author A.Beiter based on code by: + * @author Xuhao Yuan based on code by: + + * @author A.Beiter + * @author M.Needham * @author M.Needham * @author J. van Tilburg - * @date 2018-09-04 + * @date 2021-04-02 */ class UTEffChecker final @@ -68,7 +73,7 @@ public: private: void initHistogramsForChannel( unsigned int uniqueID ) const; - void layerEff( LHCb::MCParticle const*, DeUTDetector const&, LinkedFrom const&, + void layerEff( LHCb::MCParticle const*, DeUTDetector const&, LinkedFrom const&, LinkedTo const& ) const; bool isInside( const DeUTLayer*, const LHCb::MCHit*, DeUTDetector const& ) const; @@ -94,6 +99,7 @@ private: mutable HistoMap m_effYLayerHistos; mutable HistoMap m_effXYLayerHistos; mutable std::mutex m_mutex; + const unsigned int m_layerNum = 4u; }; DECLARE_COMPONENT( UTEffChecker ) @@ -114,6 +120,7 @@ StatusCode UTEffChecker::initialize() { } void UTEffChecker::initHistogramsForChannel( unsigned int uniqueID ) const { + // determined side // Create layer name std::string layerName = ", layer " + UT::toString( (int)uniqueID ); // Book x distribution histogram MCHits @@ -122,16 +129,16 @@ void UTEffChecker::initHistogramsForChannel( unsigned int uniqueID ) const { // Book y distribution histogram MCHits m_yLayerHistos.emplace( uniqueID, book1D( (int)uniqueID + 6000, "y distribution hits" + layerName, -1000., 1000., 400 ) ); - // Book x distribution histogram UTClusters + // Book x distribution histogram UTDigits m_effXLayerHistos.emplace( uniqueID, book1D( (int)uniqueID + 7000, "x distribution ass. clusters" + layerName, -1000., 1000., 400 ) ); - // Book y distribution histogram UTClusters + // Book y distribution histogram UTDigits m_effYLayerHistos.emplace( uniqueID, book1D( (int)uniqueID + 8000, "y distribution ass. clusters" + layerName, -1000., 1000., 400 ) ); // Book x vs y distribution histogram MCHits m_xyLayerHistos.emplace( uniqueID, book2D( (int)uniqueID + 9000, "x vs y distribution hits" + layerName, -1000., 1000., 200, -1000., 1000., 200 ) ); - // Book x vs y distribution histogram UTClusters + // Book x vs y distribution histogram UTDigits m_effXYLayerHistos.emplace( uniqueID, book2D( (int)uniqueID + 10000, "x vs y distribution ass. clusters" + layerName, -1000., 1000., 200, -1000., 1000., 200 ) ); } @@ -140,7 +147,7 @@ void UTEffChecker::operator()( const MCParticles& particles, const LHCb::LinksBy DeUTDetector const& det ) const { // Get the UTCluster to MCHit associator auto associator = get( LHCb::LinksByKey::linkerName( m_asctLocation ) ); - auto table = LinkedFrom( associator ); + auto table = LinkedFrom( associator ); if ( !table ) throw GaudiException( "Failed to find table", "UTEffChecker::operator()", StatusCode::FAILURE ); // Get the MCParticle to MCHit associator @@ -170,8 +177,8 @@ StatusCode UTEffChecker::finalize() { err = sqrt( ( eff * ( 100.0 - eff ) ) / (double)nAcc ); if ( !m_pEff ) eff = 1 - eff; } - info() << LHCb::UTNames::UniqueLayerToString( histoId % 100, histoId / 100 ) << " eff " << eff << " +/- " << err - << endmsg; + info() << LHCb::UTNames::UniqueLayerToString( ( histoId / 100 - 1 ) * 2 + histoId % 100 - 1 ) << " eff " << eff + << " +/- " << err << endmsg; } // layer info() << " -----------------------" << endmsg; @@ -190,31 +197,38 @@ StatusCode UTEffChecker::finalize() { } void UTEffChecker::layerEff( MCParticle const* aParticle, DeUTDetector const& det, - LinkedFrom const& table, LinkedTo const& hitTable ) const { + LinkedFrom const& table, LinkedTo const& hitTable ) const { std::scoped_lock lock( m_mutex ); // find all MC hits for this particle - auto hits = hitTable.range( aParticle ); + auto hits = hitTable.range( aParticle ); + const int utIDupper = 0x600000L; if ( hits.empty() ) return; - det.applyToAllLayers( [this, &hits, &table, &det]( DeUTLayer const& layer ) { + std::array, 4> layerHits; + + det.applyToAllLayers( [this, &hits, &table, &det, &layerHits]( DeUTLayer const& layer ) { // look for MCHit in this layer..... - std::vector layerHits; for ( const auto& aHit : hits ) { - auto hitChan = Detector::UT::ChannelID( aHit.sensDetID() ); + auto hitChan = ( aHit.sensDetID() > utIDupper ) + ? Detector::UT::ChannelID( LHCb::UTIDMapping::ReconvertID( aHit.sensDetID() ) ) + : Detector::UT::ChannelID( aHit.sensDetID() ); if ( uniqueHistoID( hitChan ) == uniqueHistoID( layer.elementID() ) && isInside( &layer, &aHit, det ) ) { - layerHits.push_back( &aHit ); + layerHits[layer.elementID().layer()].push_back( &aHit ); } }; - if ( !layerHits.empty() ) { - bool found = std::any_of( layerHits.begin(), layerHits.end(), - [&]( const MCHit* h ) { return !table.range( h ).empty(); } ); + if ( ( ( det.version() == ( DeUTDetector::GeoVersion::v0 ) ) || + ( det.version() == ( DeUTDetector::GeoVersion::v1 ) && layer.elementID().side() == 1 ) ) && + !layerHits[layer.elementID().layer()].empty() ) { + bool found = + std::any_of( layerHits[layer.elementID().layer()].begin(), layerHits[layer.elementID().layer()].end(), + [&]( const MCHit* h ) { return !table.range( h ).empty(); } ); auto histoId = uniqueHistoID( layer.elementID() ); // create histos if needed if ( m_xLayerHistos.find( histoId ) == m_xLayerHistos.end() ) { initHistogramsForChannel( histoId ); } - const MCHit* aHit = layerHits.front(); + const MCHit* aHit = layerHits[layer.elementID().layer()].front(); const Gaudi::XYZPoint midPoint = aHit->midPoint(); // histo vs x @@ -242,7 +256,7 @@ void UTEffChecker::layerEff( MCParticle const* aParticle, DeUTDetector const& de } int UTEffChecker::uniqueHistoID( const Detector::UT::ChannelID aChan ) const { - return aChan.station() * 100 + aChan.layer(); + return aChan.station() * 100 + aChan.layer() % 2 + 1u; } bool UTEffChecker::isInside( const DeUTLayer* aLayer, const MCHit* aHit, DeUTDetector const& det ) const { diff --git a/UT/UTCheckers/src/UTOccupancy.h b/UT/UTCheckers/src/UTOccupancy.h index 2c46b75672764f839ca9a3af8415de9272a0b5cb..ec3514c1d86b2a4c902fc85a7a5d434559c725f6 100644 --- a/UT/UTCheckers/src/UTOccupancy.h +++ b/UT/UTCheckers/src/UTOccupancy.h @@ -129,12 +129,13 @@ void UTOccupancy::fillHistograms( const PBASE* obj, DeUTDetector const& d int offset; - if ( aChan.detRegion() == UTNames::detRegion::RegionC ) { + if ( aChan.detRegion() == static_cast( UTNames::detRegion::RegionC ) ) { offset = utSector->column() - 1; - } else if ( aChan.detRegion() == UTNames::detRegion::RegionB ) { + } else if ( aChan.detRegion() == static_cast( UTNames::detRegion::RegionB ) ) { offset = utSector->column() - 7; - } else if ( aChan.detRegion() == UTNames::detRegion::RegionA ) { - aChan.station() == UTNames::Station::UTa ? offset = utSector->column() - 10 : offset = utSector->column() - 12; + } else if ( aChan.detRegion() == static_cast( UTNames::detRegion::RegionA ) ) { + aChan.station() == static_cast( UTNames::Station::UTa ) ? offset = utSector->column() - 10 + : offset = utSector->column() - 12; } else { this->warning() << "Unknown region " << aChan.detRegion() << endmsg; return; diff --git a/UT/UTMonitors/CMakeLists.txt b/UT/UTMonitors/CMakeLists.txt index 33c75d01c8e5ba002e0ca90a2cfd4b89904d79d6..17da75589ac88e45eefbe12186c8bf55b0a1ece0 100644 --- a/UT/UTMonitors/CMakeLists.txt +++ b/UT/UTMonitors/CMakeLists.txt @@ -15,7 +15,7 @@ UT/UTMonitors gaudi_add_module(UTMonitors SOURCES - src/UTClusterMonitor.cpp + src/UTTightDigitMonitor.cpp src/UTDataSizeMonitor.cpp src/UTErrorMonitor.cpp src/UTFullEventDump.cpp diff --git a/UT/UTMonitors/src/UTClusterMonitor.cpp b/UT/UTMonitors/src/UTTightDigitMonitor.cpp similarity index 54% rename from UT/UTMonitors/src/UTClusterMonitor.cpp rename to UT/UTMonitors/src/UTTightDigitMonitor.cpp index 1299acc1ee8989fd0225812c5b2da075643bd913..4fbbb81c7533496e635ea1dc58def2c56b66e9a5 100644 --- a/UT/UTMonitors/src/UTClusterMonitor.cpp +++ b/UT/UTMonitors/src/UTTightDigitMonitor.cpp @@ -12,9 +12,8 @@ #include "AIDA/IHistogram2D.h" #include "AIDA/IProfile1D.h" #include "AIDA/IProfile2D.h" -#include "Detector/UT/ChannelID.h" #include "Event/ODIN.h" -#include "Event/UTCluster.h" +#include "Event/UTDigit.h" #include "Kernel/IUTReadoutTool.h" #include "Kernel/UTDAQDefinitions.h" #include "Kernel/UTDetectorPlot.h" @@ -38,19 +37,28 @@ #include #include -namespace { - // empty derived condition triggering creation of histograms - struct HistoCreationCond {}; -} // namespace - +//#include // Use root histos to get MPV for sectors (in absence of something better) // Boost accumulator classes to store MPVs of sectors -/** - * @author Andy Beiter (based on code by Mark Tobin) - * @date 2018-09-04 +//----------------------------------------------------------------------------- +// Implementation file for class : UTTightDigitMonitor +// +// 2018-09-04 : Andy Beiter (based on code by Mark Tobin) +//----------------------------------------------------------------------------- + +/** @class UTTightDigitMonitor UTTightDigitMonitor.h + * + * + * @author Xuhao Yuan (based on code by Mark Tobin, Andy Beiter) + * @date 2021-04-02 */ +namespace { + // empty derived condition triggering creation of histograms + struct HistoCreationCond {}; +} // namespace + /// Define some typedefs for boost accumulators namespace UT { @@ -64,14 +72,16 @@ namespace UT { double, boost::accumulators::stats> MeanAccumulator; - class UTClusterMonitor + class UTTightDigitMonitor : public LHCb::Algorithm::Consumer< - void( LHCb::ODIN const&, LHCb::UTClusters const&, DeUTDetector const&, HistoCreationCond const& ), + void( LHCb::ODIN const&, LHCb::UTDigits const&, DeUTDetector const&, HistoCreationCond const& ), LHCb::DetDesc::usesBaseAndConditions> { public: - UTClusterMonitor( const std::string& name, ISvcLocator* pSvcLocator ); - StatusCode initialize() override; - void operator()( LHCb::ODIN const&, LHCb::UTClusters const&, DeUTDetector const&, + /// Standard constructor + UTTightDigitMonitor( const std::string& name, ISvcLocator* pSvcLocator ); + + StatusCode initialize() override; ///< Algorithm initialization + void operator()( LHCb::ODIN const&, LHCb::UTDigits const&, DeUTDetector const&, HistoCreationCond const& ) const override; private: @@ -79,14 +89,14 @@ namespace UT { bool m_histoBooked{false}; /// Fill basic histograms - void fillHistograms( const LHCb::UTCluster* cluster, std::vector& nClustersPerTELL1, + void fillHistograms( const LHCb::UTDigit* digit, std::vector& nDigitsPerTELL, DeUTDetector const& det ) const; /// Fill detailed histograms - void fillDetailedHistograms( const LHCb::UTCluster* cluster, DeUTDetector const& det ) const; + void fillDetailedHistograms( const LHCb::UTDigit* digit, DeUTDetector const& det ) const; /// Fill the hitmaps when required. - void fillClusterMaps( const LHCb::UTCluster* cluster, DeUTDetector const& det ) const; + void fillDigitMaps( const LHCb::UTDigit* digit, DeUTDetector const& det ) const; /// Fill the cluster mpv map per sector void fillMPVMap( DeUTSector const& sector, double charge ) const; @@ -94,48 +104,49 @@ namespace UT { /// Toggles to turn plots off/on Gaudi::Property m_plotBySvcBox{this, "ByServiceBox", false, "Plot by Service Box"}; Gaudi::Property m_plotByDetRegion{this, "ByDetectorRegion", false, "Plot by unique detector region"}; - Gaudi::Property m_plotByPort{this, "ByPort", false, "Plot number of clusters/tell1 by vs port"}; - Gaudi::Property m_hitMaps{this, "HitMaps", false, "True if cluster maps are to be shown"}; + Gaudi::Property m_plotByPort{this, "ByPort", false, "Plot number of digits/tell by vs port"}; + Gaudi::Property m_hitMaps{this, "HitMaps", false, "True if digit maps are to be shown"}; /// Cuts on the data quality Gaudi::Property m_chargeCut{this, "ChargeCut", 0, "Cut on charge"}; - Gaudi::Property m_minNClusters{this, "MinTotalClusters", 0, - "Cut on minimum number of clusters in the event"}; - Gaudi::Property m_maxNClusters{this, "MaxTotalClusters", std::numeric_limits::max(), - "Cut on maximum number of clusters in the event"}; + Gaudi::Property m_minNDigits{this, "MinTotalDigits", 0, + "Cut on minimum number of digits in the event"}; + Gaudi::Property m_maxNDigits{this, "MaxTotalDigits", std::numeric_limits::max(), + "Cut on maximum number of digits in the event"}; Gaudi::Property m_minMPVCharge{this, "MinMPVCharge", 0, "Minimum charge for calculation on MPV"}; Gaudi::Property m_overFlowLimit{this, "OverflowLimit", 5000, - "Overflow limit for number of clusters in event"}; + "Overflow limit for number of digits in event"}; Gaudi::Property m_resetRate{this, "ResetRate", 1000, "Reset rate for MPV calculation"}; mutable Gaudi::Accumulators::Counter<> m_eventsCount{this, "Number of events"}; - unsigned int m_nTELL1s; ///< Number of TELL1 boards expected. + unsigned int m_nTELL40s; ///< Number of TELL40 boards expected. /// Book histograms void bookHistograms( DeUTDetector const& det, IUTReadoutTool::ReadoutInfo const& roInfo ); // filled in monitor clusters - AIDA::IHistogram1D* m_1d_nClusters = nullptr; ///< Number of clusters - AIDA::IHistogram1D* m_1d_nClusters_gt_100 = nullptr; ///< Number of clusters (N > 100) - AIDA::IHistogram1D* m_1d_nClustersVsTELL1 = nullptr; ///< Number of clusters per TELL1 - AIDA::IProfile1D* m_prof_nClustersVsTELL1 = nullptr; ///< Number of clusters per TELL1 - AIDA::IHistogram2D* m_2d_nClustersVsTELL1 = nullptr; ///< Number of clusters per TELL1 + AIDA::IHistogram1D* m_1d_nDigits = nullptr; ///< Number of digits + AIDA::IHistogram1D* m_1d_nDigits_gt_100 = nullptr; ///< Number of digits (N > 100) + AIDA::IHistogram1D* m_1d_nDigitsVsTELL40 = nullptr; ///< Number of digits per TELL40 + AIDA::IProfile1D* m_prof_nDigitsVsTELL40 = nullptr; ///< Number of digits per TELL40 + AIDA::IHistogram2D* m_2d_nDigitsVsTELL40 = nullptr; ///< Number of digits per TELL40 - AIDA::IHistogram1D* m_1d_nClustersVsTELL1Links = nullptr; ///< Number of clusters per TELL1 (split into links) - AIDA::IHistogram1D* m_1d_nClustersVsTELL1Sectors = nullptr; ///< Number of clusters per TELL1 (split into sectors) + AIDA::IHistogram1D* m_1d_nDigitsVsTELL40Links = nullptr; ///< Number of digits per TELL40 (split into links) + AIDA::IHistogram1D* m_1d_nDigitsVsTELL40Sectors = nullptr; ///< Number of digits per TELL40 (split into sectors) - AIDA::IHistogram1D* m_1d_nClusters_overflow = nullptr; ///< Number of clusters where last bin contains overflow info + AIDA::IHistogram1D* m_1d_nDigits_overflow = nullptr; ///< Number of digits where last bin contains overflow info // filled in fillHistograms - AIDA::IHistogram1D* m_1d_ClusterSize = nullptr; ///< Cluster Size - AIDA::IHistogram2D* m_2d_ClusterSizeVsTELL1 = nullptr; ///< Cluster Size vs TELL1 - AIDA::IHistogram2D* m_2d_UTNVsTELL1 = nullptr; ///< Signal to Noise vs TELL1 - AIDA::IHistogram2D* m_2d_ChargeVsTELL1 = nullptr; ///< Cluster Charge vs TELL1 - AIDA::IHistogram2D* m_2d_ClustersPerPortVsTELL1 = nullptr; ///< Clusters per port vs TELL1 + // do we need these two? + AIDA::IHistogram1D* m_1d_DigitSize = nullptr; ///< Digit Size + AIDA::IHistogram2D* m_2d_DigitSizeVsTELL40 = nullptr; ///< Digit Size vs TELL40 + AIDA::IHistogram2D* m_2d_UTNVsTELL40 = nullptr; ///< Signal to Noise vs TELL40 + AIDA::IHistogram2D* m_2d_ChargeVsTELL40 = nullptr; ///< Digit Charge vs TEL401 + AIDA::IHistogram2D* m_2d_DigitsPerPortVsTELL40 = nullptr; ///< Digits per port vs TELL40 - AIDA::IHistogram1D* m_1d_totalCharge = nullptr; ///< Total charge in cluster + AIDA::IHistogram1D* m_1d_totalCharge = nullptr; ///< Total charge in digit /// Plots by service box std::map m_1ds_chargeByServiceBox; @@ -143,7 +154,7 @@ namespace UT { std::map m_1ds_chargeByDetRegion; /// Hitmaps - AIDA::IHistogram2D* m_2d_hitmap = nullptr; ///< Cluster hitmap + AIDA::IHistogram2D* m_2d_hitmap = nullptr; ///< Digit hitmap AIDA::IProfile1D* m_prof_sectorMPVs = nullptr; ///< Sector MPV vs arbitrary sector number AIDA::IProfile1D* m_prof_sectorTruncMean1 = nullptr; ///< Sector MPV calculated using truncated mean losing first and last 15% @@ -158,33 +169,39 @@ namespace UT { unsigned int m_nBeetlePortsPerSector; ///< Number of beetle ports per sector static const unsigned int m_nBeetlePortsPerUTSector = 16u; ///< Number of bins per UT sector in the hitmap (beetle ports) - unsigned int m_nSectorsPerTELL1; + unsigned int m_nSectorsPerTELL40; /// Accumulation of statistics for the MPV per sector std::map m_1ds_chargeBySector; mutable std::map> m_chargeBySector; /// Contains the bin corresponding to a given sectors in the 1d representation of the hitmap std::map m_sectorBins1D; - AIDA::IHistogram1D* m_1d_nClustersVsSector = nullptr; ///< Number of clusters in each sector - AIDA::IHistogram1D* m_1d_nClustersVsBeetlePort = nullptr; ///< Number of clusters in each beetle port + AIDA::IHistogram1D* m_1d_nDigitsVsSector = nullptr; ///< Number of digits in each sector + AIDA::IHistogram1D* m_1d_nDigitsVsBeetlePort = nullptr; ///< Number of digits in each beetle port mutable std::mutex m_mutex; ToolHandle m_readoutTool{this, "ReadoutTool", "UTReadoutTool"}; }; // Declaration of the Algorithm Factory - DECLARE_COMPONENT( UTClusterMonitor ) + DECLARE_COMPONENT( UTTightDigitMonitor ) } // namespace UT -UT::UTClusterMonitor::UTClusterMonitor( const std::string& name, ISvcLocator* pSvcLocator ) +//============================================================================= +// Standard constructor, initializes variables +//============================================================================= +UT::UTTightDigitMonitor::UTTightDigitMonitor( const std::string& name, ISvcLocator* pSvcLocator ) : Consumer{name, pSvcLocator, - {{"ODINLocation", LHCb::ODINLocation::Default}, - {"ClusterLocation", LHCb::UTClusterLocation::UTClusters}, + {KeyValue{"ODINLocation", LHCb::ODINLocation::Default}, + KeyValue{"TightDigitLocation", LHCb::UTDigitLocation::UTTightDigits}, {"UTLocation", DeUTDetLocation::location()}, {"HistoCreationCondLocation", "UTClusterMonitorHistoCreationCond"}}} {} -StatusCode UT::UTClusterMonitor::initialize() { +//============================================================================= +// Initialization +//============================================================================= +StatusCode UT::UTTightDigitMonitor::initialize() { return Consumer::initialize().andThen( [&] { if ( "" == histoTopDir() ) setHistoTopDir( "UT/" ); // Turn all histograms on in expert mode @@ -196,7 +213,7 @@ StatusCode UT::UTClusterMonitor::initialize() { } m_nBeetlePortsPerSector = m_nBeetlePortsPerUTSector; - m_nSectorsPerTELL1 = UTDAQ::noptlinks * UTDAQ::nports / m_nBeetlePortsPerSector; + m_nSectorsPerTELL40 = UTDAQ::noptlinks * UTDAQ::nports / m_nBeetlePortsPerSector; // Triggers booking of histograms when geometry is loaded addConditionDerivation( @@ -211,43 +228,46 @@ StatusCode UT::UTClusterMonitor::initialize() { } ); } -void UT::UTClusterMonitor::operator()( LHCb::ODIN const& odin, LHCb::UTClusters const& clusters, - DeUTDetector const& det, HistoCreationCond const& ) const { +//============================================================================= +// Main execution +//============================================================================= +void UT::UTTightDigitMonitor::operator()( LHCb::ODIN const& odin, LHCb::UTDigits const& digits, DeUTDetector const& det, + HistoCreationCond const& ) const { plot1D( odin.bunchId(), "BCID", "BCID", -0.5, 2807.5, 2808 ); ++m_eventsCount; - /// This plots cluster ADCs vs Sampling time - /// Store number of clusters/TELL1 (48 (42) Tell1s, 1->48 - std::vector nClustersPerTELL1( m_nTELL1s, 0 ); + /// This plots digit ADCs vs Sampling time + /// Store number of digits/TELL40 (48 (42) Tell1s, 1->48 + std::vector nDigitsPerTELL40( m_nTELL40s, 0 ); - const unsigned int nClusters = clusters.size(); - debug() << "Number of clusters in " << inputLocation() << " is " << nClusters << endmsg; + const unsigned int nDigits = digits.size(); + debug() << "Number of digits in " << inputLocation() << " is " << nDigits << endmsg; - if ( nClusters < m_minNClusters || nClusters > m_maxNClusters ) return; + if ( nDigits < m_minNDigits || nDigits > m_maxNDigits ) return; - m_1d_nClusters->fill( nClusters ); - if ( 100 < nClusters ) { m_1d_nClusters_gt_100->fill( nClusters ); } - if ( m_overFlowLimit < nClusters ) - m_1d_nClusters_overflow->fill( m_overFlowLimit ); + m_1d_nDigits->fill( nDigits ); + if ( 100 < nDigits ) { m_1d_nDigits_gt_100->fill( nDigits ); } + if ( m_overFlowLimit < nDigits ) + m_1d_nDigits_overflow->fill( m_overFlowLimit ); else - m_1d_nClusters_overflow->fill( nClusters ); - // Loop over clusters - - for ( const auto& cluster : clusters ) { - fillHistograms( cluster, nClustersPerTELL1, det ); - if ( fullDetail() ) fillDetailedHistograms( cluster, det ); - } // End of cluster iterator - - // Fill histogram for number of clusters/TELL1 - unsigned int TELL1 = 1; - for ( auto itClPerTELL1 = nClustersPerTELL1.begin(); itClPerTELL1 != nClustersPerTELL1.end(); - ++itClPerTELL1, ++TELL1 ) { - verbose() << "TELL1: " << TELL1 << ",clusters: " << ( *itClPerTELL1 ) << endmsg; - unsigned int nClusters = ( *itClPerTELL1 ); - if ( 0 < nClusters ) { - m_2d_nClustersVsTELL1->fill( TELL1, nClusters ); - m_prof_nClustersVsTELL1->fill( TELL1, nClusters ); + m_1d_nDigits_overflow->fill( nDigits ); + // Loop over digits + + for ( const auto& digit : digits ) { + fillHistograms( digit, nDigitsPerTELL40, det ); + if ( fullDetail() ) fillDetailedHistograms( digit, det ); + } // End of digit iterator + + // Fill histogram for number of digits/TELL40 + unsigned int TELL40 = 1; + for ( auto itClPerTELL40 = nDigitsPerTELL40.begin(); itClPerTELL40 != nDigitsPerTELL40.end(); + ++itClPerTELL40, ++TELL40 ) { + verbose() << "TELL40: " << TELL40 << ",digits: " << ( *itClPerTELL40 ) << endmsg; + unsigned int nDigits = ( *itClPerTELL40 ); + if ( 0 < nDigits ) { + m_2d_nDigitsVsTELL40->fill( TELL40, nDigits ); + m_prof_nDigitsVsTELL40->fill( TELL40, nDigits ); } } } @@ -259,7 +279,7 @@ void UT::UTClusterMonitor::operator()( LHCb::ODIN const& odin, LHCb::UTClusters // 2) Truncated mean using only the 1st 70% of the distribution // 3) Binned mean with 20 bins from 0 to 100 ADC counts //================================================================================================================================ -void UT::UTClusterMonitor::fillMPVMap( DeUTSector const& sector, double charge ) const { +void UT::UTTightDigitMonitor::fillMPVMap( DeUTSector const& sector, double charge ) const { unsigned int sectorID = sector.elementID().uniqueSector(); if ( charge > m_minMPVCharge ) { auto lock = std::scoped_lock{m_mutex}; @@ -308,39 +328,43 @@ void UT::UTClusterMonitor::fillMPVMap( DeUTSector const& sector, double charge ) } } -void UT::UTClusterMonitor::bookHistograms( DeUTDetector const& det, IUTReadoutTool::ReadoutInfo const& roInfo ) { +//============================================================================== +// Book histograms +//============================================================================== +void UT::UTTightDigitMonitor::bookHistograms( DeUTDetector const& det, IUTReadoutTool::ReadoutInfo const& roInfo ) { // make sure we fill Histograms only once ! // In case we load 2 different geometries, we'll be called twice. But in that case, our histograms make no sense if ( m_histoBooked ) { throw GaudiException( "Double loading of UT geometry detected. This is not supported", - "UTClusterMonitor::bookHistograms", StatusCode::FAILURE ); + "UTTightDigitMonitor::bookHistograms", StatusCode::FAILURE ); } m_histoBooked = true; - m_nTELL1s = m_readoutTool->nBoard( &roInfo ); + m_nTELL40s = m_readoutTool->nBoard( &roInfo ); + ; - // filled in monitor clusters - m_1d_nClusters = book1D( "Number of clusters", 0., 20000., 2000 ); - m_1d_nClusters_gt_100 = book1D( "Number of clusters (N > 100)", 0., 20000., 2000 ); - m_1d_nClusters_overflow = book1D( "Number of clusters (no overflow)", 0., m_overFlowLimit + 1, 1000 ); - m_2d_nClustersVsTELL1 = book2D( "Number of clusters per TELL1", 0.5, m_nTELL1s + 0.5, m_nTELL1s, 0., 100., 50 ); + // filled in monitor digits + m_1d_nDigits = book1D( "Number of digits", 0., 20000., 2000 ); + m_1d_nDigits_gt_100 = book1D( "Number of digits (N > 100)", 0., 20000., 2000 ); + m_1d_nDigits_overflow = book1D( "Number of digits (no overflow)", 0., m_overFlowLimit + 1, 1000 ); + m_2d_nDigitsVsTELL40 = book2D( "Number of digits per TELL40", 0.5, m_nTELL40s + 0.5, m_nTELL40s, 0., 100., 50 ); - // Number of clusters produced by each TELL1 - m_1d_nClustersVsTELL1 = book1D( "Number of clusters vs TELL1", 0.5, m_nTELL1s + 0.5, m_nTELL1s ); - m_prof_nClustersVsTELL1 = bookProfile1D( "Mean number of clusters per TELL1", 0.5, m_nTELL1s + 0.5, m_nTELL1s ); - m_1d_nClustersVsTELL1Links = - book1D( "Number of clusters vs TELL1 links", 1., m_nTELL1s + 1., m_nTELL1s * UTDAQ::noptlinks ); - m_1d_nClustersVsTELL1Sectors = - book1D( "Number of clusters vs TELL1 sectors", 1., m_nTELL1s + 1., m_nTELL1s * m_nSectorsPerTELL1 ); + // Number of digits produced by each TELL40 + m_1d_nDigitsVsTELL40 = book1D( "Number of digits vs TELL40", 0.5, m_nTELL40s + 0.5, m_nTELL40s ); + m_prof_nDigitsVsTELL40 = bookProfile1D( "Mean number of digits per TELL40", 0.5, m_nTELL40s + 0.5, m_nTELL40s ); + m_1d_nDigitsVsTELL40Links = + book1D( "Number of digits vs TELL40 links", 1., m_nTELL40s + 1., m_nTELL40s * UTDAQ::noptlinks ); + m_1d_nDigitsVsTELL40Sectors = + book1D( "Number of digits vs TELL40 sectors", 1., m_nTELL40s + 1., m_nTELL40s * m_nSectorsPerTELL40 ); // filled in fillHistograms - m_1d_ClusterSize = book1D( "Cluster Size", 0.5, 4.5, 4 ); - m_2d_ClusterSizeVsTELL1 = book2D( "Cluster Size vs TELL1", 0.5, m_nTELL1s + 0.5, m_nTELL1s, 0.5, 4.5, 4 ); - m_2d_UTNVsTELL1 = book2D( "Signal to Noise vs TELL1", 0.5, m_nTELL1s + 0.5, m_nTELL1s, 0., 100., 25 ); - m_2d_ChargeVsTELL1 = book2D( "Cluster Charge vs TELL1", 0.5, m_nTELL1s + 0.5, m_nTELL1s, 0., 60., 60 ); + m_1d_DigitSize = book1D( "Digit Size", 0.5, 4.5, 4 ); + m_2d_DigitSizeVsTELL40 = book2D( "Digit Size vs TELL40", 0.5, m_nTELL40s + 0.5, m_nTELL40s, 0.5, 4.5, 4 ); + m_2d_UTNVsTELL40 = book2D( "Signal to Noise vs TELL40", 0.5, m_nTELL40s + 0.5, m_nTELL40s, 0., 100., 25 ); + m_2d_ChargeVsTELL40 = book2D( "Cluster Charge vs TELL40", 0.5, m_nTELL40s + 0.5, m_nTELL40s, 0., 60., 60 ); if ( m_plotByPort ) { - m_2d_ClustersPerPortVsTELL1 = book2D( "Clusters per port vs TELL1", 0.5, m_nTELL1s + 0.5, m_nTELL1s, 0., 96., 96 ); + m_2d_DigitsPerPortVsTELL40 = book2D( "Digits per port vs TELL40", 0.5, m_nTELL40s + 0.5, m_nTELL40s, 0., 96., 96 ); } m_1d_totalCharge = book1D( "Cluster ADC Values", 0., 200., 200 ); /// list of service boxes @@ -352,13 +376,15 @@ void UT::UTClusterMonitor::bookHistograms( DeUTDetector const& det, IUTReadoutTo } // End of service box loop if ( m_plotByDetRegion ) { - for ( const auto& region : LHCb::UTNames::allDetRegions() ) { +#if 0 + for ( const auto& region : LHCb::UTNames().allDetRegions() ) { m_1ds_chargeByDetRegion[region] = book1D( "Cluster ADC Values " + region, 0., 200., 200 ); }; +#endif } // Book histograms for hitmaps if ( m_hitMaps ) { - std::string idMap = "UT cluster map"; + std::string idMap = "UT digit map"; std::string idMPVMap = "UT Sector MPVs"; UT::UTDetectorPlot hitMap( idMap, idMap, m_nBeetlePortsPerUTSector ); m_2d_hitmap = book2D( hitMap.name(), hitMap.minBinX(), hitMap.maxBinX(), hitMap.nBinX(), hitMap.minBinY(), @@ -376,119 +402,133 @@ void UT::UTClusterMonitor::bookHistograms( DeUTDetector const& det, IUTReadoutTo std::string title = "Total charge in " + sector.nickname(); m_1ds_chargeBySector[sector.elementID().uniqueSector()] = Gaudi::Utils::Aida2ROOT::aida2root( book1D( "BySector/" + idh, title, 0., 100., 50 ) ); + // std::cout << title+":" << (*Sectors)->elementID().uniqueSector() << "," << bin << std::endl; m_sectorBins1D[sector.elementID().uniqueSector()] = bin; ++bin; } ); // 1d representation for online monitoring - m_nSectors = m_sectorBins1D.size(); - m_1d_nClustersVsSector = book1D( "Number of clusters per sector", 0.5, m_nSectors + 0.5, m_nSectors ); - m_1d_nClustersVsBeetlePort = - book1D( "Number of clusters per beetle port", 0.5, m_nSectors + 0.5, m_nSectors * m_nBeetlePortsPerSector ); + m_nSectors = m_sectorBins1D.size(); + m_1d_nDigitsVsSector = book1D( "Number of digits per sector", 0.5, m_nSectors + 0.5, m_nSectors ); + m_1d_nDigitsVsBeetlePort = + book1D( "Number of digits per beetle port", 0.5, m_nSectors + 0.5, m_nSectors * m_nBeetlePortsPerSector ); m_prof_sectorMPVs = bookProfile1D( "Cluster MPV vs Sector", 0.5, m_nSectors + 0.5, m_nSectors ); m_prof_sectorTruncMean1 = bookProfile1D( "Cluster trunc mean vs Sector", 0.5, m_nSectors + 0.5, m_nSectors ); m_prof_sectorTruncMean2 = - bookProfile1D( "Cluster trunc mean (1st 70%) vs Sector", 0.5, m_nSectors + 0.5, m_nSectors ); - m_prof_sectorBinnedMPV = bookProfile1D( "Cluster binned MPV vs Sector", 0.5, m_nSectors + 0.5, m_nSectors ); + bookProfile1D( "Digit trunc mean (1st 70%) vs Sector", 0.5, m_nSectors + 0.5, m_nSectors ); + m_prof_sectorBinnedMPV = bookProfile1D( "Digit binned MPV vs Sector", 0.5, m_nSectors + 0.5, m_nSectors ); } // end of hitmaps } -void UT::UTClusterMonitor::fillHistograms( LHCb::UTCluster const* cluster, std::vector& nClustersPerTELL1, - DeUTDetector const& det ) const { - const double totalCharge = cluster->totalCharge(); +//============================================================================== +// Fill histograms +//============================================================================== +void UT::UTTightDigitMonitor::fillHistograms( const LHCb::UTDigit* digit, std::vector& nDigitsPerTELL40, + DeUTDetector const& det ) const { + + const double totalCharge = digit->depositedCharge(); if ( totalCharge < m_chargeCut ) return; // calculate MPVs // unsigned int sectorID=cluster->channelID().uniqueSector(); - auto sector = det.findSector( cluster->channelID() ); + auto sector = det.findSector( digit->channelID() ); if ( m_hitMaps ) { - fillClusterMaps( cluster, det ); + fillDigitMaps( digit, det ); fillMPVMap( *sector, totalCharge ); unsigned int sectorID = sector->elementID().uniqueSector(); unsigned int bin = m_sectorBins1D.at( sectorID ); - m_1d_nClustersVsSector->fill( bin ); - unsigned int port = ( cluster->channelID().strip() - 1 ) / UTDAQ::nstrips; - m_1d_nClustersVsBeetlePort->fill( ( bin - 0.5 ) + static_cast( port ) / m_nBeetlePortsPerSector ); + m_1d_nDigitsVsSector->fill( bin ); + unsigned int port = ( digit->channelID().strip() - sector->firstStrip() ) / UTDAQ::nstrips; + m_1d_nDigitsVsBeetlePort->fill( ( bin - 0.5 ) + static_cast( port ) / m_nBeetlePortsPerSector ); } - const unsigned int clusterSize = cluster->size(); + const unsigned int digitSize = 1u; - const unsigned int sourceID = cluster->sourceID(); - unsigned int TELL1ID = ( m_readoutTool )->SourceIDToTELLNumber( sourceID ); - nClustersPerTELL1[TELL1ID - 1] += 1; - m_1d_nClustersVsTELL1->fill( TELL1ID ); - unsigned int tell1Link = cluster->tell1Channel() / 128; - m_1d_nClustersVsTELL1Links->fill( ( TELL1ID ) + ( ( 0.5 + tell1Link ) / UTDAQ::noptlinks ) ); - unsigned int tell1Sector = cluster->tell1Channel() / ( m_nBeetlePortsPerSector * UTDAQ::nstrips ); - m_1d_nClustersVsTELL1Sectors->fill( ( TELL1ID ) + ( ( 0.5 + tell1Sector ) / m_nSectorsPerTELL1 ) ); + // const unsigned int sourceID = cluster->sourceID(); + const UTDAQID utdaqID = m_readoutTool->channelIDToDAQID( digit->channelID() ); + const unsigned int tell40channel = UTDAQ::nStripsInUTSector * utdaqID.lane() + utdaqID.channel(); + unsigned int TELL40ID = utdaqID.board(); + nDigitsPerTELL40[TELL40ID] += 1; + // nDigitsPerTELL40[TELL40ID - 1] += 1; + m_1d_nDigitsVsTELL40->fill( TELL40ID ); + unsigned int tell40Link = tell40channel / 128; + m_1d_nDigitsVsTELL40Links->fill( ( TELL40ID ) + ( ( 0.5 + tell40Link ) / UTDAQ::noptlinks ) ); + unsigned int tell40Sector = tell40channel / ( m_nBeetlePortsPerSector * UTDAQ::nstrips ); + m_1d_nDigitsVsTELL40Sectors->fill( ( TELL40ID ) + ( ( 0.5 + tell40Sector ) / m_nSectorsPerTELL40 ) ); - m_1d_ClusterSize->fill( clusterSize ); - m_2d_ClusterSizeVsTELL1->fill( TELL1ID, clusterSize ); + m_1d_DigitSize->fill( digitSize ); + m_2d_DigitSizeVsTELL40->fill( TELL40ID, digitSize ); - double noise = sector->noise( cluster->channelID() ); + double noise = sector->noise( digit->channelID() ); if ( noise > 1e-3 ) { - const double signalToNoise = cluster->totalCharge() / noise; - m_2d_UTNVsTELL1->fill( TELL1ID, signalToNoise ); + const double signalToNoise = digit->depositedCharge() / noise; + m_2d_UTNVsTELL40->fill( TELL40ID, signalToNoise ); } else - Warning( "Zero S/N for some clusters", StatusCode::SUCCESS, 1 ).ignore(); - m_2d_ChargeVsTELL1->fill( TELL1ID, totalCharge ); + Warning( "Zero S/N for some digits", StatusCode::SUCCESS, 1 ).ignore(); + m_2d_ChargeVsTELL40->fill( TELL40ID, totalCharge ); if ( m_plotByPort ) { - const unsigned int tell1Channel = cluster->tell1Channel(); - unsigned int port = tell1Channel / UTDAQ::nstrips; - m_2d_ClustersPerPortVsTELL1->fill( TELL1ID, port ); + // const unsigned int tell1Channel = cluster->tell1Channel(); + unsigned int port = tell40channel / UTDAQ::nstrips; + m_2d_DigitsPerPortVsTELL40->fill( TELL40ID, port ); } // Always fill histograms per readout quadrant // Get service box and set up histogram IDs m_1d_totalCharge->fill( totalCharge ); - std::string svcBox = m_readoutTool->serviceBox( cluster->channelID() ); + std::string svcBox = m_readoutTool->serviceBox( digit->channelID() ); std::string quadrant = svcBox.substr( 0, 2 ); m_1ds_chargeByServiceBox.at( quadrant )->fill( totalCharge ); if ( m_plotBySvcBox ) m_1ds_chargeByServiceBox.at( svcBox )->fill( totalCharge ); - if ( m_plotByDetRegion ) m_1ds_chargeByDetRegion.at( cluster->detRegionName() )->fill( totalCharge ); } - -void UT::UTClusterMonitor::fillDetailedHistograms( const LHCb::UTCluster* cluster, DeUTDetector const& det ) const { +//============================================================================== +/// Fill more detailed histograms +//============================================================================== +void UT::UTTightDigitMonitor::fillDetailedHistograms( const LHCb::UTDigit* digit, DeUTDetector const& det ) const { // high threshold - plot( cluster->highThreshold(), "High threshold", "High threshold", -0.5, 1.5, 2 ); + // plot( digit->highThreshold(), "High threshold", "High threshold", -0.5, 1.5, 2 ); // by strip, modulo a few things.... - const unsigned int strip = cluster->strip(); + const unsigned int strip = ( digit->channelID() ).strip(); plot1D( strip % 8, "strip modulo 8", "strip modulo 8", -0.5, 8.5, 9 ); plot1D( strip % 32, "strip modulo 32", "strip modulo 32", -0.5, 32.5, 33 ); plot1D( strip % 128, "strip modulo 128", "strip modulo 128", -0.5, 128.5, 129 ); // histogram by station - const int station = cluster->station(); - plot1D( station, "Number of clusters per station", "Number of clusters per station", -0.5, 4.5, 5 ); + const int station = ( digit->channelID() ).station(); + plot1D( station, "Number of digits per station", "Number of digits per station", -0.5, 4.5, 5 ); // by layer - const int layer = cluster->layer(); - plot( (double)( 10 * station + layer ), "Number of clusters per layer", "Number of clusters per layer", -0.5, 40.5, - 41 ); + const int layer = digit->layer(); + plot( (double)( 10 * station + layer ), "Number of digits per layer", "Number of digits per layer", -0.5, 40.5, 41 ); +#if 0 const double neighbourSum = cluster->neighbourSum(); plot1D( cluster->pseudoSize(), "Pseudo size of cluster", "Pseudo size of cluster", -0.5, 10.5, 11 ); +#endif - const double totalCharge = cluster->totalCharge(); - plot1D( neighbourSum / totalCharge, "Relative neighbour sum", "Relative neighbour sum", -1.02, 1.02, 51 ); + const double totalCharge = digit->depositedCharge(); + // plot1D( neighbourSum / totalCharge, "Relative neighbour sum", "Relative neighbour sum", -1.02, 1.02, 51 ); - const unsigned int clusterSize = cluster->size(); - if ( 3 > clusterSize ) { +#if 0 + const unsigned int digitSize = 1u; + if ( 3 > digitSize ) { plot1D( neighbourSum, "Neighbour sum (1- and 2-strip clusters)", "Neighbour sum (1- and 2-strip clusters)", -16.5, 26.5, 43 ); plot1D( neighbourSum / totalCharge, "Relative neighbour sum (1- and 2-strip clusters)", "Relative neighbour sum (1- and 2-strip clusters)", -1.02, 1.02, 51 ); } +#endif // charge and S/N - plot1D( totalCharge, cluster->layerName() + "/Charge", "Charge", -0.5, 500.5, 501 ); - auto sector = det.findSector( cluster->channelID() ); - double noise = sector->noise( cluster->channelID() ); + plot1D( totalCharge, digit->layerName() + "/Charge", "Charge", -0.5, 500.5, 501 ); + auto sector = det.findSector( digit->channelID() ); + double noise = sector->noise( digit->channelID() ); if ( noise > 1e-3 ) { - const double signalToNoise = cluster->totalCharge() / noise; - plot1D( signalToNoise, cluster->layerName() + "/Signal to noise", "S/N", 0., 100., 100 ); + const double signalToNoise = digit->depositedCharge() / noise; + plot1D( signalToNoise, LHCb::UTNames().UniqueLayerToString( digit->channelID() ) + "/Signal to noise", "S/N", 0., + 100., 100 ); } else - Warning( "Some cluster have zero noise", StatusCode::SUCCESS, 1 ).ignore(); - // Plot cluster ADCs for 1, 2, 3, 4 strip clusters + Warning( "Some digit have zero noise", StatusCode::SUCCESS, 1 ).ignore(); + // Plot digit ADCs for 1, 2, 3, 4 strip digits +#if 0 std::string svcBox = m_readoutTool->serviceBox( cluster->channelID() ); std::string idhStrip = " (" + std::to_string( clusterSize ) + " strip)"; std::string idh; @@ -501,16 +541,19 @@ void UT::UTClusterMonitor::fillDetailedHistograms( const LHCb::UTCluster* cluste idh = "BySector/s" + std::to_string( cluster->size() ) + "/charge_$sector" + std::to_string( cluster->channelID().uniqueSector() ); plot1D( totalCharge, idh, "Total charge in " + cluster->sectorName(), 0., 200., 200 ); +#endif } - -void UT::UTClusterMonitor::fillClusterMaps( const LHCb::UTCluster* cluster, DeUTDetector const& det ) const { - const std::string idMap = "UT cluster map"; - auto utSector = det.findSector( cluster->channelID() ); +//============================================================================== +/// Fill histograms of digit maps +//============================================================================== +void UT::UTTightDigitMonitor::fillDigitMaps( const LHCb::UTDigit* digit, DeUTDetector const& det ) const { + const std::string idMap = "UT digit map"; + auto utSector = det.findSector( digit->channelID() ); // make the real hit map UT::UTDetectorPlot hitMap( idMap, idMap, m_nBeetlePortsPerUTSector ); UT::UTDetectorPlot::Bins bins = hitMap.toBins( *utSector ); - int port = ( cluster->channelID().strip() - 1 ) / UTDAQ::nstrips; + int port = ( digit->channelID().strip() - utSector->firstStrip() ) / UTDAQ::nstrips; double xBin = bins.xBin - double( port + 1 ) / m_nBeetlePortsPerUTSector + 0.5; int nBins = bins.endBinY - bins.beginBinY; diff --git a/UT/UTTools/src/UTOfflinePosition.cpp b/UT/UTTools/src/UTOfflinePosition.cpp index 1dbeccaf0e19a6a9549b4f94144b906a0c33fabe..39ca912ae92a9399734e3202a3d1e42540bfdaad 100644 --- a/UT/UTTools/src/UTOfflinePosition.cpp +++ b/UT/UTTools/src/UTOfflinePosition.cpp @@ -130,9 +130,7 @@ IUTClusterPosition::Info UTOfflinePosition::estimate( const DeUTDetector& det, c auto info = UTFun::position( adcVector, m_trim ); double stripNum = info.first; - Detector::UT::ChannelID theChan = {firstChan.type(), firstChan.station(), - firstChan.layer(), firstChan.detRegion(), - firstChan.sector(), (unsigned int)stripNum + firstChan.strip()}; + Detector::UT::ChannelID theChan( firstChan.channelID() + ( (unsigned int)stripNum ) ); double stripFrac = stripFraction( stripNum - floor( stripNum ), info.second ); @@ -160,8 +158,7 @@ IUTClusterPosition::Info UTOfflinePosition::estimate( const DeUTDetector& det, c double stripNum = info.first; Detector::UT::ChannelID firstChan = digits.front()->channelID(); - Detector::UT::ChannelID theChan = {firstChan.type(), firstChan.station(), firstChan.layer(), - firstChan.detRegion(), firstChan.sector(), (unsigned int)stripNum}; + Detector::UT::ChannelID theChan( firstChan.channelID() - firstChan.strip() + (unsigned int)stripNum ); double stripFrac = stripFraction( stripNum - floor( stripNum ), info.second ); diff --git a/UT/UTTools/src/UTOnlinePosition.cpp b/UT/UTTools/src/UTOnlinePosition.cpp index 62215634bd5e76de0d3c7fd26cf2b1be844fb94e..4aa741810307c151e1eb2360c23fee2404f87406 100644 --- a/UT/UTTools/src/UTOnlinePosition.cpp +++ b/UT/UTTools/src/UTOnlinePosition.cpp @@ -60,9 +60,7 @@ IUTClusterPosition::Info UTOnlinePosition::estimate( const DeUTDetector&, const } Detector::UT::ChannelID firstChan = aCluster->firstChannel(); - Detector::UT::ChannelID theChan = {firstChan.type(), firstChan.station(), - firstChan.layer(), firstChan.detRegion(), - firstChan.sector(), (unsigned int)stripNum + firstChan.strip()}; + Detector::UT::ChannelID theChan( firstChan.channelID() + (unsigned int)stripNum ); IUTClusterPosition::Info theInfo; theInfo.strip = theChan; @@ -84,8 +82,7 @@ IUTClusterPosition::Info UTOnlinePosition::estimate( const DeUTDetector&, const interStripPos = 0; } - auto theChan = Detector::UT::ChannelID( firstChan.type(), firstChan.station(), firstChan.layer(), - firstChan.detRegion(), firstChan.sector(), (unsigned int)stripNum ); + Detector::UT::ChannelID theChan( firstChan.channelID() - firstChan.strip() + (unsigned int)stripNum ); IUTClusterPosition::Info theInfo; theInfo.strip = theChan; diff --git a/UT/UTTools/src/UTSelectChannelIDByElement.cpp b/UT/UTTools/src/UTSelectChannelIDByElement.cpp index ef0823935dac8932a149c9191a3c41bd80179b83..7a7173f315972fd17ae7e022559bf0dc40d84a17 100644 --- a/UT/UTTools/src/UTSelectChannelIDByElement.cpp +++ b/UT/UTTools/src/UTSelectChannelIDByElement.cpp @@ -16,7 +16,7 @@ #include "UTDet/DeUTDetector.h" #include "UTDet/DeUTLayer.h" #include "UTDet/DeUTSector.h" -#include "UTDet/DeUTStation.h" +#include "UTDet/DeUTSide.h" #include "GaudiAlg/FunctionalTool.h" #include "GaudiKernel/IBinder.h" @@ -26,9 +26,9 @@ namespace { // generic handler for Station or Layer ot Sector #ifdef USE_DD4HEP - using UTElement = std::variant; + using UTElement = std::variant; #else - using UTElement = std::variant; + using UTElement = std::variant; #endif /** @@ -41,11 +41,11 @@ namespace { struct Contains { LHCb::Detector::UT::ChannelID id; #ifdef USE_DD4HEP - bool operator()( DeUTStation const& s ) { return s.contains( id ); } + bool operator()( DeUTSide const& s ) { return s.contains( id ); } bool operator()( DeUTLayer const& l ) { return l.contains( id ); } bool operator()( DeUTSector const& s ) { return s.contains( id ); } #else - bool operator()( DeUTStation const* s ) { return s->contains( id ); } + bool operator()( DeUTSide const* s ) { return s->contains( id ); } bool operator()( DeUTLayer const* l ) { return l->contains( id ); } bool operator()( DeUTSector const* s ) { return s->contains( id ); } #endif @@ -57,7 +57,11 @@ namespace { for ( const auto& name : names ) { if ( name.find( "UT" ) != std::string_view::npos ) { const LHCb::Detector::UT::ChannelID chan = LHCb::UTNames::stringToChannel( name ); +#ifdef USE_DD4HEP if ( chan.sector() != 0 ) { +#else + { +#endif auto s = det.findSector( chan ); if ( s ) { #ifdef USE_DD4HEP @@ -68,7 +72,11 @@ namespace { continue; } } +#ifdef USE_DD4HEP if ( chan.layer() != 0 ) { +#else + { +#endif auto l = det.findLayer( chan ); if ( l ) { #ifdef USE_DD4HEP @@ -79,8 +87,13 @@ namespace { continue; } } +#ifdef USE_DD4HEP if ( chan.station() != 0u ) { auto s = det.findStation( chan ); +#else + { + auto s = det.findSide( chan ); +#endif if ( s ) { #ifdef USE_DD4HEP elements.emplace_back( *s );