diff --git a/Phys/ParticleMaker/src/ParticleMassDTFMonitor.cpp b/Phys/ParticleMaker/src/ParticleMassDTFMonitor.cpp index ca89cbd70ec0dd8087988f7a8bd3f1cc41ceee6a..6871c2b92d1b187cf3981c9438725cb1f315e6a2 100644 --- a/Phys/ParticleMaker/src/ParticleMassDTFMonitor.cpp +++ b/Phys/ParticleMaker/src/ParticleMassDTFMonitor.cpp @@ -12,14 +12,16 @@ #include "DecayTreeFitter/Fitter.h" #include "DetDesc/DetectorElement.h" #include "Event/Particle.h" -#include "GaudiAlg/GaudiHistoAlg.h" #include "Kernel/IParticlePropertySvc.h" #include "Kernel/ParticleProperty.h" #include "LHCbAlgs/Consumer.h" -using base_type = - LHCb::Algorithm::Consumer>; +#include + +#include + +using base_type = LHCb::Algorithm::Consumer>; class ParticleMassDTFMonitor : public base_type { @@ -29,6 +31,13 @@ public: pSvcLocator, {KeyValue{"Particles", ""}, KeyValue{"StandardGeometryTop", LHCb::standard_geometry_top}}} {} + StatusCode initialize() override { + return base_type::initialize().andThen( [&] { + m_histo.emplace( this, m_histoName, m_histoName, + Gaudi::Accumulators::Axis{m_bins, m_mean - 20 * m_sigma, m_mean + 20 * m_sigma} ); + } ); + } + void operator()( LHCb::Particle::Range const& particles, DetectorElement const& lhcb ) const override { // loop over particles auto b_mass = m_mass.buffer(); @@ -52,7 +61,7 @@ public: b_mass += m; // fill histo with 5sigma range on each side of mean - plot1D( m, m_histoName, m_mean - 20 * m_sigma, m_mean + 20 * m_sigma, m_bins ); + ++( *m_histo )[m]; // count candidates in 1 sigma if ( std::abs( m_mean - m ) / m_sigma < 1 ) b_1sig++; @@ -61,7 +70,7 @@ public: private: // properties - Gaudi::Property m_bins{this, "Bins", 20}; + Gaudi::Property m_bins{this, "Bins", 20}; Gaudi::Property m_mean{this, "MassMean", 5279.65 * Gaudi::Units::MeV}; Gaudi::Property m_sigma{this, "MassSigma", 50 * Gaudi::Units::MeV}; Gaudi::Property m_histoName{this, "histoName", "Invariant Mass"}; @@ -71,6 +80,8 @@ private: mutable Gaudi::Accumulators::StatCounter m_mass{this, "Candidate mass"}; mutable Gaudi::Accumulators::Counter<> m_n1sigma{this, "Candidates in 1 sigma"}; mutable Gaudi::Accumulators::Counter<> m_ndtf{this, "Candidates with succesful DTF fit"}; + // histograms + mutable std::optional> m_histo{}; // property service ServiceHandle m_particlePropertySvc{this, "ParticleProperty", "LHCb::ParticlePropertySvc"}; diff --git a/Pr/PrMCTools/src/VPClusterEfficiency.cpp b/Pr/PrMCTools/src/VPClusterEfficiency.cpp index d463e2027ce30371b8638d0b29d4433c572c9d57..473b57ec11562e35576f2da262f612e5af3b9758 100644 --- a/Pr/PrMCTools/src/VPClusterEfficiency.cpp +++ b/Pr/PrMCTools/src/VPClusterEfficiency.cpp @@ -15,23 +15,23 @@ #include "Event/RawBank.h" #include "Event/VPDigit.h" #include "Event/VPFullCluster.h" -#include "GaudiAlg/GaudiHistoAlg.h" #include "LHCbAlgs/Consumer.h" #include "VPDet/DeVP.h" + +#include "Gaudi/Accumulators/HistogramArray.h" +#include "Gaudi/Accumulators/StaticHistogram.h" + #include -/** @class VPClusterEfficiency VPClusterEfficiency.h - * +/** * Checks the VPCluster efficiency for simulation */ - class VPClusterEfficiency : public LHCb::Algorithm::Consumer&, const LHCb::MCHits&, const LHCb::MCParticles&, const LHCb::LinksByKey&, const LHCb::MCProperty&, const DeVP& ), - LHCb::Algorithm::Traits::usesBaseAndConditions> { + LHCb::Algorithm::Traits::usesConditions> { public: - /// Standard constructor VPClusterEfficiency( const std::string& name, ISvcLocator* pSvcLocator ); /// Consumer operator: takes VP clusters & MCHits to make efficiency plots @@ -47,6 +47,429 @@ private: mutable Gaudi::Accumulators::SigmaCounter<> m_residual_x{this, "Residuals x [mm]"}; mutable Gaudi::Accumulators::SigmaCounter<> m_residual_y{this, "Residuals y [mm]"}; mutable Gaudi::Accumulators::SigmaCounter<> m_purity{this, "Purity"}; + + // MCHit related histograms + mutable Gaudi::Accumulators::HistogramArray, 52> m_mcHitXYPerModule{ + this, + []( int n ) { return fmt::format( "All MCHit x,y VP module {}", n ); }, + []( int n ) { return fmt::format( "All MCHit x,y VP module {}", n ); }, + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_mcHitXY{ + this, "All MCHit x,y VP all modules", "All MCHit x,y VP all modules", {240, -60, 60}, {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_mcHitDistTravel{this, + "Distance travelled by MCHit in VP sensor (mm)", + "Distance travelled by MCHit in VP sensor (mm)", + {100, 0, 1}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_mcHitEnergyDeposit{this, + "Energy depotisted by MCHit in VP sensor (MeV)", + "Energy depotisted by MCHit in VP sensor (MeV)", + {100, 0, 1}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_mcHitEffXYPerModule{ + this, + []( int n ) { return fmt::format( "Efficiency MCHit -> Pixel x,y VP module {}", n ); }, + []( int n ) { return fmt::format( "Efficiency MCHit -> Pixel x,y VP module {}", n ); }, + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<2> m_mcHitEffXY{this, + "Efficiency MCHit -> Pixel x,y VP all modules", + "Efficiency MCHit -> Pixel x,y VP all modules", + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_mcHitEffvDist{ + this, + "Efficiency (MCHit->Pixel) v distance travelled by MCHit in VP sensor (mm)", + "Efficiency (MCHit->Pixel) v distance travelled by MCHit in VP sensor (mm)", + {100, 0, 1}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_mcHitEffvEnergy{ + this, + "Efficiency (MCHit->Pixel) v energy depotisted by MCHit in VP sensor (MeV)", + "Efficiency (MCHit->Pixel) v energy depotisted by MCHit in VP sensor (MeV)", + {100, 0, 1}}; + + // Cluster related histograms + mutable Gaudi::Accumulators::StaticHistogram<1> m_nHitsPerPixel{ + this, "Number of MCHits per VP pixel", "Number of MCHits per VP pixel", {10, -0.5, 9.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_nMCParticlesPerPixel{ + this, "Number of MCParticles per VP pixel", "Number of MCParticles per VP pixel", {10, -0.5, 9.5}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_nMCParticlesPerPixelModule{ + this, + "Number of MCParticles per VP Pixel v module", + "Number of MCParticles per VP Pixel v module", + {10, -0.5, 9.5}, + {52, -0.5, 51.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_nParticlesPerCluster{ + this, "Number of MCParticles per VP Cluster", "Number of MCParticles per VP Cluster", {10, -0.5, 9.5}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_nParticlesPerClusterModule{ + this, + "Number of MCParticles per VP Cluster v module", + "Number of MCParticles per VP Cluster v module", + {10, -0.5, 9.5}, + {52, -0.5, 51.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_spillPerCluster{ + this, "Fraction pixels spill or noise per cluster", "Fraction pixels spill or noise per cluster", {26, 0., 1.04}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_spillPerClusterSize{ + this, + "Fraction pixels spill or noise per cluster v cluster size", + "Fraction pixels spill or noise per cluster v cluster size", + {26, 0., 1.04}, + {50, 0.5, 50.5}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_goodCluterPos{this, + "Good (>70% true) Clusters pos all modules", + "Good (>70% true) Clusters pos all modules", + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_goodCluterPosPerModule{ + this, + []( int n ) { return fmt::format( "Good (>70% true) Clusters pos module {}", n ); }, + []( int n ) { return fmt::format( "Good (>70% true) Clusters pos module {}", n ); }, + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_badCluterPos{this, + "Bad (<70% true) Clusters pos all modules", + "Bad (<70% true) Clusters pos all modules", + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_badCluterPosPerModule{ + this, + []( int n ) { return fmt::format( "Bad (<70% true) Clusters pos module {}", n ); }, + []( int n ) { return fmt::format( "Bad (<70% true) Clusters pos module {}", n ); }, + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_pixelXY{ + this, "Pixel xy all sensors", "Pixel xy all sensors", {240, -60, 60}, {240, -60, 60}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_pixelXYPerModule{ + this, + []( int n ) { return fmt::format( "Pixel xy sensor {}", n ); }, + []( int n ) { return fmt::format( "Pixel xy sensor {}", n ); }, + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_nPixelNotUsed{this, + "Number of pixels not in clusters per event", + "Number of pixels not in clusters per event", + {100, 0., 10000.}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_fracPixelNotInEvent{this, + "Fraction of pixels not in clusters per event", + "Fraction of pixels not in clusters per event", + {100, 0., 1.}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_pixelNotinCluster{this, + "Pixel not in cluster xy all sensors", + "Pixel not in cluster xy all sensors", + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_pixelNotinClusterPerModule{ + this, + []( int n ) { return fmt::format( "Pixel not in cluster xy sensor {}", n ); }, + []( int n ) { return fmt::format( "Pixel not in cluster xy sensor {}", n ); }, + {240, -60, 60}, + {240, -60, 60}}; + + // missed MCHits histograms + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedHitsEta{ + this, "# MCHits not found - eta", "# MCHits not found - eta", {70, -5.5, 5.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedHitsPhi{ + this, "# MCHits not found - phi", "# MCHits not found - phi", {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedHitsPhikk{ + this, "# MCHits not found - phi - phikk", "# MCHits not found - phi - phikk", {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedHitsPhigee{ + this, "# MCHits not found - phi - gammaee", "# MCHits not found - phi - gammaee", {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedHitsR{ + this, "# MCHits not found - r", "# MCHits not found - r", {50, 5.1, 50}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedHitsR7{ + this, "# MCHits not found - r<7mm", "# MCHits not found - r<7mm", {14, 5.1, 7}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedHitsP{ + this, "# MCHits not found - p", "# MCHits not found - p", {120, 0, 100}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedHitsPt{ + this, "# MCHits not found - pt", "# MCHits not found - pt", {120, 0, 10}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedHitsZ{ + this, "# MCHits not found - z", "# MCHits not found - z", {90, -300, 800}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedHitsModule{ + this, "# MCHits not found - module", "# MCHits not found - module", {52, -0.5, 51.5}}; + + // recontructed MCHits histograms + mutable Gaudi::Accumulators::StaticHistogram<1> m_recHitsEta{ + this, "# MCHits reconstructible - eta", "# MCHits reconstructible - eta", {70, -5.5, 5.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recHitsPhi{ + this, "# MCHits reconstructible - phi", "# MCHits reconstructible - phi", {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recHitsPhikk{ + this, "# MCHits reconstructible - phi - phikk", "# MCHits reconstructible - phi - phikk", {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recHitsPhigee{ + this, "# MCHits reconstructible - phi - gammaee", "# MCHits reconstructible - phi - gammaee", {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recHitsR{ + this, "# MCHits reconstructible - r", "# MCHits reconstructible - r", {50, 5.1, 50}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recHitsR7{ + this, "# MCHits reconstructible - r<7mm", "# MCHits reconstructible - r<7mm", {14, 5.1, 7}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recHitsP{ + this, "# MCHits reconstructible - p", "# MCHits reconstructible - p", {120, 0, 100}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recHitsPt{ + this, "# MCHits reconstructible - pt", "# MCHits reconstructible - pt", {120, 0, 10}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recHitsZ{ + this, "# MCHits reconstructible - z", "# MCHits reconstructible - z", {90, -300, 800}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recHitsModule{ + this, "# MCHits reconstructible - module", "# MCHits reconstructible - module", {52, -0.5, 51.5}}; + + // missed MCHits from Velo reconstructible tracks histograms + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedVeloHitsEta{ + this, "# MCHits not found Velo reco - eta", "# MCHits not found Velo reco - eta", {70, -5.5, 5.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedVeloHitsPhi{ + this, "# MCHits not found Velo reco - phi", "# MCHits not found Velo reco - phi", {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedVeloHitsPhikk{this, + "# MCHits not found Velo reco - phi - phikk", + "# MCHits not found Velo reco - phi - phikk", + {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedVeloHitsPhigee{this, + "# MCHits not found Velo reco - phi - gammaee", + "# MCHits not found Velo reco - phi - gammaee", + {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedVeloHitsR{ + this, "# MCHits not found Velo reco - r", "# MCHits not found Velo reco - r", {50, 5.1, 50}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedVeloHitsR7{ + this, "# MCHits not found Velo reco - r<7mm", "# MCHits not found Velo reco - r<7mm", {14, 5.1, 7}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedVeloHitsP{ + this, "# MCHits not found Velo reco - p", "# MCHits not found Velo reco - p", {120, 0, 100}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedVeloHitsPt{ + this, "# MCHits not found Velo reco - pt", "# MCHits not found Velo reco - pt", {120, 0, 10}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedVeloHitsZ{ + this, "# MCHits not found Velo reco - z", "# MCHits not found Velo reco - z", {90, -300, 800}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_missedVeloHitsModule{ + this, "# MCHits not found Velo reco - module", "# MCHits not found Velo reco - module", {52, -0.5, 51.5}}; + + // recontructed MCHits from Velo reconstructible tracks histograms + mutable Gaudi::Accumulators::StaticHistogram<1> m_recVeloHitsEta{ + this, "# MCHits reconstructible Velo reco - eta", "# MCHits reconstructible Velo reco - eta", {70, -5.5, 5.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recVeloHitsPhi{ + this, "# MCHits reconstructible Velo reco - phi", "# MCHits reconstructible Velo reco - phi", {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recVeloHitsPhikk{this, + "# MCHits reconstructible Velo reco - phi - phikk", + "# MCHits reconstructible Velo reco - phi - phikk", + {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recVeloHitsPhigee{ + this, + "# MCHits reconstructible Velo reco - phi - gammaee", + "# MCHits reconstructible Velo reco - phi - gammaee", + {21, -M_PI, M_PI}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recVeloHitsR{ + this, "# MCHits reconstructible Velo reco - r", "# MCHits reconstructible Velo reco - r", {50, 5.1, 50}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recVeloHitsR7{ + this, "# MCHits reconstructible Velo reco - r<7mm", "# MCHits reconstructible Velo reco - r<7mm", {14, 5.1, 7}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recVeloHitsP{ + this, "# MCHits reconstructible Velo reco - p", "# MCHits reconstructible Velo reco - p", {120, 0, 100}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recVeloHitsPt{ + this, "# MCHits reconstructible Velo reco - pt", "# MCHits reconstructible Velo reco - pt", {120, 0, 10}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recVeloHitsZ{ + this, "# MCHits reconstructible Velo reco - z", "# MCHits reconstructible Velo reco - z", {90, -300, 800}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recVeloHitsModule{this, + "# MCHits reconstructible Velo reco - module", + "# MCHits reconstructible Velo reco - module", + {52, -0.5, 51.5}}; + + // recontructed MCHits from Muon histograms + mutable Gaudi::Accumulators::StaticHistogram<1> m_muonR7{this, "# muons r<7mm", "# muons r<7mm", {11, -0.5, 10.5}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_muonEff{ + this, "Efficiency vs # #mu (at least 2 #mu)", "Efficiency vs # #mu (at least 2 #mu)", {11, -0.5, 10.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_muonSemsorsMu{ + this, "Sensors w at least 2 #mu", "Sensors w at least 2 #mu", {208, -0.5, 207.5}}; + + // cluster residual histograms + mutable Gaudi::Accumulators::StaticHistogram<1> m_resX{ + this, "Residuals along x [mm]", "Residuals along x [mm]", {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_resY{ + this, "Residuals along y [mm]", "Residuals along y [mm]", {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_resXPerModule{ + this, + []( int n ) { return fmt::format( "Residuals along x [mm] - module{}", n ); }, + []( int n ) { return fmt::format( "Residuals along x [mm] - module{}", n ); }, + {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_resYPerModule{ + this, + []( int n ) { return fmt::format( "Residuals along y [mm] - module{}", n ); }, + []( int n ) { return fmt::format( "Residuals along y [mm] - module{}", n ); }, + {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_resXVelo{ + this, "Residuals along x - VELO reco [mm]", "Residuals along x - VELO reco [mm]", {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_resYVelo{ + this, "Residuals along y - VELO reco [mm]", "Residuals along y - VELO reco [mm]", {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_resXLong{ + this, "Residuals along x - long [mm]", "Residuals along x - long [mm]", {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_resYLong{ + this, "Residuals along y - long [mm]", "Residuals along y - long [mm]", {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_resXLongPixel{ + this, "Residuals along x - long pixels [mm]", "Residuals along x - long pixels [mm]", {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_resYLongPixel{ + this, "Residuals along y - long pixels [mm]", "Residuals along y - long pixels [mm]", {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_resXEdge{ + this, "Residuals along x - sensor edges [mm]", "Residuals along x - sensor edges [mm]", {200, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_resYEdge{ + this, "Residuals along y - sensor edges [mm]", "Residuals along y - sensor edges [mm]", {200, -0.2, 0.2}}; + + // matrix edge histograms + mutable Gaudi::Accumulators::StaticHistogram<2> m_noRecClosestCluter2D{this, + "Non reconstructed MCHit - closest cluster", + "Non reconstructed MCHit - closest cluster", + {36, -0.99, 0.99}, + {36, -0.99, 0.99}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_noRecClosestCluter1D{ + this, + "Non reconstructed MCHit - closest cluster - 1D", + "Non reconstructed MCHit - closest cluster - 1D", + {48, 0., 0.66}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_noRecClosestRecHit{ + this, + "Non reconstructed MCHit - closest reconstructible MCHit", + "Non reconstructed MCHit - closest reconstructible MCHit", + {48, 0., 0.66}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recClosestRecHit{ + this, + "Reconstructible MCHit - closest reconstructible MCHit", + "Reconstructible MCHit - closest reconstructible MCHit", + {48, 0., 0.66}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recClosestRecHitExt{ + this, + "Reconstructible MCHit - closest reconstructible MCHit - extended", + "Reconstructible MCHit - closest reconstructible MCHit - extended", + {3830, 0., 52.8}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_noRecClosestRecHitphikk{ + this, + "Non reconstructed MCHit - closest reconstructible MCHit - phikk", + "Non reconstructed MCHit - closest reconstructible MCHit - phikk", + {48, 0., 0.66}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recClosestRecHitphikk{ + this, + "Reconstructible MCHit - closest reconstructible MCHit - phikk", + "Reconstructible MCHit - closest reconstructible MCHit - phikk", + {48, 0., 0.66}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_noRecClosestRecHitgee{ + this, + "Non reconstructed MCHit - closest reconstructible MCHit - gammaee", + "Non reconstructed MCHit - closest reconstructible MCHit - gammaee", + {48, 0., 0.66}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_recClosestRecHitgee{ + this, + "Reconstructible MCHit - closest reconstructible MCHit - gammaee", + "Reconstructible MCHit - closest reconstructible MCHit - gammaee", + {48, 0., 0.66}}; + + // distribution number of clusters histograms + mutable Gaudi::Accumulators::StaticHistogram<1> m_nClusters{this, "# clusters", "# clusters", {10000, 0, 10000}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_nUnmatchedClusters{ + this, "# unmatched clusters", "# unmatched clusters", {100, 0, 10000}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_nClustersR{ + this, "# clusters distribution - r", "# clusters distribution - r", {50, 5.1, 50}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_nClustersRPerModule{ + this, + []( int n ) { return fmt::format( "# cluster distribution - r - module{}", n ); }, + []( int n ) { return fmt::format( "# cluster distribution - r - module{}", n ); }, + {50, 5.1, 50}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_nClustersModule{ + this, "# clusters distribution - module", "# clusters distribution - module", {52, -0.5, 51.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_nPixelPerCluster{ + this, "# pixel per cluster distribution", "# pixel per cluster distribution", {51, -0.5, 50.5}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_nPixelPerClusterPerModule{ + this, + []( int n ) { return fmt::format( "# pixel per cluster distribution - module{}", n ); }, + []( int n ) { return fmt::format( "# pixel per cluster distribution - module{}", n ); }, + {51, -0.5, 50.5}}; + + // purity histograms + mutable Gaudi::Accumulators::StaticHistogram<1> m_purityHisto{ + this, "Purity distribution", "Purity distribution", {110, 0, 1.1}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_nPixels{ + this, "# pixel distribution - seen", "# pixel distribution - seen", {20, 0, 100}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_purityPixels{ + this, "Purity vs # pixels", "Purity vs # pixels", {20, 0, 100}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_nPixelsInCluster{ + this, "Number of pixels in cluster", "Number of pixels in cluster", {33, 0, 32}}; + + // pixel stat histograms + mutable Gaudi::Accumulators::StaticHistogram<1> m_cluterTohit{ + this, "# associated cluster to hit", "# associated cluster to hit", {11, -0.5, 10.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_pixels{ + this, "# pixel distribution", "# pixel distribution", {10, 0, 100}}; + + // efficiency vs occupancy histograms + mutable Gaudi::Accumulators::HistogramArray, 208> m_recMCHitPerSensor{ + this, + []( int n ) { return fmt::format( "Reconstructible MCHit x,y VP sensor {}", n ); }, + []( int n ) { return fmt::format( "Reconstructible MCHit x,y VP sensor {}", n ); }, + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_recMCHit{this, + "Reconstructible MCHit x,y VP all sensors", + "Reconstructible MCHit x,y VP all sensors", + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::HistogramArray, 208> m_recMCHitVeloPerSensor{ + this, + []( int n ) { return fmt::format( "Reconstructible MCHit x,y VP sensor - VELO {}", n ); }, + []( int n ) { return fmt::format( "Reconstructible MCHit x,y VP sensor - VELO {}", n ); }, + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_recMCHitVelo{this, + "Reconstructible MCHit x,y VP all sensors - VELO", + "Reconstructible MCHit x,y VP all sensors - VELO", + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::HistogramArray, 208> m_missedMCHitPerSensor{ + this, + []( int n ) { return fmt::format( "Missed MCHit x,y VP sensor {}", n ); }, + []( int n ) { return fmt::format( "Missed MCHit x,y VP sensor {}", n ); }, + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_missedMCHit{ + this, "Missed MCHit x,y VP all sensors", "Missed MCHit x,y VP all sensors", {240, -60, 60}, {240, -60, 60}}; + mutable Gaudi::Accumulators::HistogramArray, 208> m_missedMCHitVeloPerSensor{ + this, + []( int n ) { return fmt::format( "Missed MCHit x,y VP sensor - VELO {}", n ); }, + []( int n ) { return fmt::format( "Missed MCHit x,y VP sensor - VELO {}", n ); }, + {240, -60, 60}, + {240, -60, 60}}; + mutable Gaudi::Accumulators::StaticHistogram<2> m_missedMCHitVelo{this, + "Missed MCHit x,y VP all sensors - VELO", + "Missed MCHit x,y VP all sensors - VELO", + {240, -60, 60}, + {240, -60, 60}}; + + mutable Gaudi::Accumulators::HistogramArray, 52> m_recHitsRPerModule{ + this, + []( int n ) { return fmt::format( "# MCHits reconstructible - r - module{}", n ); }, + []( int n ) { return fmt::format( "# MCHits reconstructible - r - module{}", n ); }, + {50, 5.1, 50}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_recHitsR7PerModule{ + this, + []( int n ) { return fmt::format( "# MCHits reconstructible - r<7mm - module{}", n ); }, + []( int n ) { return fmt::format( "# MCHits reconstructible - r<7mm - module{}", n ); }, + {14, 5.1, 7}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_recHitsVeloRPerModule{ + this, + []( int n ) { return fmt::format( "# MCHits reconstructible Velo reco - r - module{}", n ); }, + []( int n ) { return fmt::format( "# MCHits reconstructible Velo reco - r - module{}", n ); }, + {50, 5.1, 50}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_recHitsVeloR7PerModule{ + this, + []( int n ) { return fmt::format( "# MCHits reconstructible Velo reco - r<7mm - module{}", n ); }, + []( int n ) { return fmt::format( "# MCHits reconstructible Velo reco - r<7mm - module{}", n ); }, + {14, 5.1, 7}}; + + mutable Gaudi::Accumulators::HistogramArray, 52> m_missedHitsRPerModule{ + this, + []( int n ) { return fmt::format( "# MCHits not found - r - module{}", n ); }, + []( int n ) { return fmt::format( "# MCHits not found - r - module{}", n ); }, + {50, 5.1, 50}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_missedHitsR7PerModule{ + this, + []( int n ) { return fmt::format( "# MCHits not found - r<7mm - module{}", n ); }, + []( int n ) { return fmt::format( "# MCHits not found - r<7mm - module{}", n ); }, + {14, 5.1, 7}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_missedHitsVeloRPerModule{ + this, + []( int n ) { return fmt::format( "# MCHits not found Velo reco - r - module{}", n ); }, + []( int n ) { return fmt::format( "# MCHits not found Velo reco - r - module{}", n ); }, + {50, 5.1, 50}}; + mutable Gaudi::Accumulators::HistogramArray, 52> m_missedHitsVeloR7PerModule{ + this, + []( int n ) { return fmt::format( "# MCHits not found Velo reco - r<7mm - module{}", n ); }, + []( int n ) { return fmt::format( "# MCHits not found Velo reco - r<7mm - module{}", n ); }, + {14, 5.1, 7}}; }; DECLARE_COMPONENT( VPClusterEfficiency ) @@ -157,33 +580,26 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra } } + unsigned int modNum = 0; for ( const auto& hitsInModule : hitsInModules ) { - std::string strModNum; - if ( hitsInModule.size() != 0 ) strModNum = std::to_string( hitsInModule[0]->sensDetID() / 4 ); for ( const auto& mcH : hitsInModule ) { // all MCHits - plot2D( mcH->midPoint().x(), mcH->midPoint().y(), "All MCHit x,y VP module " + strModNum, -60, 60, -60, 60, 240, - 240 ); - plot2D( mcH->midPoint().x(), mcH->midPoint().y(), "All MCHit x,y VP all modules ", -60, 60, -60, 60, 240, 240 ); + ++m_mcHitXYPerModule[modNum][{mcH->midPoint().x(), mcH->midPoint().y()}]; + ++m_mcHitXY[{mcH->midPoint().x(), mcH->midPoint().y()}]; // distance particle traveled in sensor & energy auto dist = mcH->pathLength(); - plot1D( dist, "Distance travelled by MCHit in VP sensor (mm)", 0, 1, 100 ); + ++m_mcHitDistTravel[dist]; // energy auto energy = mcH->energy(); - plot1D( energy, "Energy depotisted by MCHit in VP sensor (MeV)", 0, 1, 100 ); + ++m_mcHitEnergyDeposit[energy]; // were MCHits found double efficiency = ( channelIdForMCHit.find( mcH ) != channelIdForMCHit.end() ); - - profile2D( mcH->midPoint().x(), mcH->midPoint().y(), efficiency, - "Efficiency MCHit -> Pixel x,y VP module " + strModNum, -60, 60, -60, 60, 240, 240 ); - profile2D( mcH->midPoint().x(), mcH->midPoint().y(), efficiency, "Efficiency MCHit -> Pixel x,y VP all modules", - -60, 60, -60, 60, 240, 240 ); - - profile1D( dist, efficiency, "Efficiency (MCHit->Pixel) v distance travelled by MCHit in VP sensor (mm)", 0, 1, - 100 ); - profile1D( energy, efficiency, "Efficiency (MCHit->Pixel) v energy depotisted by MCHit in VP sensor (MeV)", 0, 1, - 100 ); + m_mcHitEffXYPerModule[modNum][{mcH->midPoint().x(), mcH->midPoint().y()}] += efficiency; + m_mcHitEffXY[{mcH->midPoint().x(), mcH->midPoint().y()}] += efficiency; + m_mcHitEffvDist[dist] += efficiency; + m_mcHitEffvEnergy[energy] += efficiency; } + modNum++; } // what pixels are in the event -- copy them all now, then delete as they are found in clusters @@ -197,7 +613,7 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra for ( auto& channelID : clus.pixels() ) { unsigned int numMC = ( MCHitForchannelId.find( channelID.channelID() )->second.size() ); - plot1D( (double)numMC, "Number of MCHits per VP pixel", -0.5, 9.5, 10 ); + ++m_nHitsPerPixel[numMC]; // count MCParticles contributing to Cluster std::set mcKeysPix; for ( auto& mcHit : ( *MCHitForchannelId.find( channelID.channelID() ) ).second ) { @@ -210,10 +626,8 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra mcKeysPix.insert( particle->key() ); } unsigned int numMCP = mcKeysPix.size(); - plot1D( (double)numMCP, "Number of MCParticles per VP pixel", -0.5, 9.5, 10 ); - plot2D( (double)numMCP, (double)channelID.module(), "Number of MCParticles per VP Pixel v module", -0.5, 9.5, - -0.5, 51.5, 10, 52 ); - + ++m_nMCParticlesPerPixel[numMCP]; + ++m_nMCParticlesPerPixelModule[{numMCP, channelID.module()}]; // check if pixel is from MCHit if ( numMC != 0 ) { ++mainPix; @@ -224,23 +638,19 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra for ( auto& key : mcKeysPix ) mcKeysClus.insert( key ); // add to list of cluster MCParticles } double numMCParticles = (double)mcKeysClus.size(); - plot1D( numMCParticles, "Number of MCParticles per VP Cluster", -0.5, 9.5, 10 ); - plot2D( numMCParticles, (double)clus.channelID().module(), "Number of MCParticles per VP Cluster v module", -0.5, - 9.5, -0.5, 51.5, 10, 52 ); + ++m_nParticlesPerCluster[numMCParticles]; + ++m_nParticlesPerClusterModule[{numMCParticles, clus.channelID().module()}]; double fracOther = ( (double)otherPix ) / ( (double)( otherPix + mainPix ) ); - plot1D( fracOther, "Fraction pixels spill or noise per cluster", 0., 1.04, 26 ); - plot2D( fracOther, clus.pixels().size(), "Fraction pixels spill or noise per cluster v cluster size", 0., 1.04, 0.5, - 50.5, 26, 50 ); - + ++m_spillPerCluster[fracOther]; + ++m_spillPerClusterSize[{fracOther, clus.pixels().size()}]; // plot positions of "good" clusters i.e. > 70% from MCHits - unsigned int module = clus.channelID().module(); - auto strModNum = std::to_string( module ); + unsigned int module = clus.channelID().module(); if ( fracOther < 0.3 ) { - plot2D( clus.x(), clus.y(), "Good (>70% true) Clusters pos all modules", -60, 60, -60, 60, 240, 240 ); - plot2D( clus.x(), clus.y(), "Good (>70% true) Clusters pos module " + strModNum, -60, 60, -60, 60, 240, 240 ); + ++m_goodCluterPos[{clus.x(), clus.y()}]; + ++m_goodCluterPosPerModule[module][{clus.x(), clus.y()}]; } else { - plot2D( clus.x(), clus.y(), "Bad (<70% true) Clusters pos all modules", -60, 60, -60, 60, 240, 240 ); - plot2D( clus.x(), clus.y(), "Bad (<70% true) Clusters pos module " + strModNum, -60, 60, -60, 60, 240, 240 ); + ++m_badCluterPos[{clus.x(), clus.y()}]; + ++m_badCluterPosPerModule[module][{clus.x(), clus.y()}]; } } @@ -248,26 +658,21 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra for ( auto& channelID : digits ) { const DeVPSensor& sensor = det.sensor( channelID.sensor() ); Gaudi::XYZPoint pointGlobal = sensor.channelToGlobalPoint( channelID ); - plot2D( pointGlobal.x(), pointGlobal.y(), "Pixel xy all sensors", -60, 60, -60, 60, 240, 240 ); - auto strModNum = std::to_string( channelID.module() ); - plot2D( pointGlobal.x(), pointGlobal.y(), "Pixel xy sensor " + strModNum, -60, 60, -60, 60, 240, 240 ); + ++m_pixelXY[{pointGlobal.x(), pointGlobal.y()}]; + ++m_pixelXYPerModule[channelID.module()][{pointGlobal.x(), pointGlobal.y()}]; } // lost pixels double nPixNotUsed = (double)( pixUsed.size() ); if ( msgLevel( MSG::DEBUG ) ) { debug() << "Found " << nPixNotUsed << " VP Digits unused in clusters" << endmsg; } - plot1D( nPixNotUsed, "Number of pixels not in clusters per event", 0., 10000, 100 ); - if ( digits.size() > 0 ) { - plot1D( nPixNotUsed / (double)digits.size(), "Fraction of pixels not in clusters per event", 0., 1., 100 ); - } + ++m_nPixelNotUsed[nPixNotUsed]; + if ( digits.size() > 0 ) { ++m_fracPixelNotInEvent[nPixNotUsed / (double)digits.size()]; } for ( auto& channelID : pixUsed ) { const DeVPSensor& sensor = det.sensor( channelID.sensor() ); Gaudi::XYZPoint pointGlobal = sensor.channelToGlobalPoint( channelID ); - plot2D( pointGlobal.x(), pointGlobal.y(), "Pixel not in cluster xy all sensors", -60, 60, -60, 60, 240, 240 ); - auto strModNum = std::to_string( channelID.module() ); - plot2D( pointGlobal.x(), pointGlobal.y(), "Pixel not in cluster xy sensor " + strModNum, -60, 60, -60, 60, 240, - 240 ); + ++m_pixelNotinCluster[{pointGlobal.x(), pointGlobal.y()}]; + ++m_pixelNotinClusterPerModule[channelID.module()][{pointGlobal.x(), pointGlobal.y()}]; } // cluster efficiency @@ -334,105 +739,99 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra // plots for missed MCHits for ( auto& mcH : hitsMissed ) { - plot1D( mcH->mcParticle()->momentum().Eta(), "# MCHits not found - eta", -5.5, 5.5, 70 ); - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits not found - phi", -M_PI, M_PI, 21 ); + ++m_missedHitsEta[mcH->mcParticle()->momentum().Eta()]; + ++m_missedHitsPhi[mcH->mcParticle()->momentum().phi()]; if ( nullptr != mcH->mcParticle()->originVertex()->mother() ) { if ( abs( mcH->mcParticle()->particleID().pid() ) == 321 && abs( mcH->mcParticle()->originVertex()->mother()->particleID().pid() ) == 333 ) { - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits not found - phi - phikk", -M_PI, M_PI, 21 ); + ++m_missedHitsPhikk[mcH->mcParticle()->momentum().phi()]; } if ( abs( mcH->mcParticle()->particleID().pid() ) == 11 && abs( mcH->mcParticle()->originVertex()->mother()->particleID().pid() ) == 22 ) { - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits not found - phi - gammaee", -M_PI, M_PI, 21 ); + ++m_missedHitsPhigee[mcH->mcParticle()->momentum().phi()]; } } - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits not found - r", 5.1, 50, 50 ); - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits not found - r<7mm", 5.1, 7, 14 ); - plot1D( mcH->mcParticle()->p() / 1000.0, "# MCHits not found - p", 0, 100, 120 ); - plot1D( mcH->mcParticle()->pt() / 1000.0, "# MCHits not found - pt", 0, 10, 120 ); - plot1D( mcH->midPoint().z(), "# MCHits not found - z", -300, 800, 90 ); - plot1D( mcH->sensDetID() / 4, "# MCHits not found - module", -0.5, 51.5, 52 ); + auto r = sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ); + ++m_missedHitsR[r]; + ++m_missedHitsR7[r]; + ++m_missedHitsP[mcH->mcParticle()->p() / 1000.0]; + ++m_missedHitsPt[mcH->mcParticle()->pt() / 1000.0]; + ++m_missedHitsZ[mcH->midPoint().z()]; + ++m_missedHitsModule[mcH->sensDetID() / 4]; } // plots for reconstructed MCHits for ( auto& mcH : hitsReconstructible ) { - plot1D( mcH->mcParticle()->momentum().Eta(), "# MCHits reconstructible - eta", -5.5, 5.5, 70 ); - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits reconstructible - phi", -M_PI, M_PI, 21 ); + ++m_recHitsEta[mcH->mcParticle()->momentum().Eta()]; + ++m_recHitsPhi[mcH->mcParticle()->momentum().phi()]; if ( nullptr != mcH->mcParticle()->originVertex()->mother() ) { if ( abs( mcH->mcParticle()->particleID().pid() ) == 321 && abs( mcH->mcParticle()->originVertex()->mother()->particleID().pid() ) == 333 ) { - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits reconstructible - phi - phikk", -M_PI, M_PI, 21 ); + ++m_recHitsPhikk[mcH->mcParticle()->momentum().phi()]; } if ( abs( mcH->mcParticle()->particleID().pid() ) == 11 && abs( mcH->mcParticle()->originVertex()->mother()->particleID().pid() ) == 22 ) { - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits reconstructible - phi - gammaee", -M_PI, M_PI, 21 ); + ++m_recHitsPhigee[mcH->mcParticle()->momentum().phi()]; } } - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits reconstructible - r", 5.1, 50, 50 ); - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits reconstructible - r<7mm", 5.1, 7, 14 ); - plot1D( mcH->mcParticle()->p() / 1000.0, "# MCHits reconstructible - p", 0, 100, 120 ); - plot1D( mcH->mcParticle()->pt() / 1000.0, "# MCHits reconstructible - pt", 0, 10, 120 ); - plot1D( mcH->midPoint().z(), "# MCHits reconstructible - z", -300, 800, 90 ); - plot1D( mcH->sensDetID() / 4, "# MCHits reconstructible - module", -0.5, 51.5, 52 ); + auto r = sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ); + ++m_recHitsR[r]; + ++m_recHitsR7[r]; + ++m_recHitsP[mcH->mcParticle()->p() / 1000.0]; + ++m_recHitsPt[mcH->mcParticle()->pt() / 1000.0]; + ++m_recHitsZ[mcH->midPoint().z()]; + ++m_recHitsModule[mcH->sensDetID() / 4]; } // plots for missed MCHits from Velo reconstructible track for ( auto& mcH : hitsMissedVeloReco ) { - plot1D( mcH->mcParticle()->momentum().Eta(), "# MCHits not found Velo reco - eta", -5.5, 5.5, 70 ); - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits not found Velo reco - phi", -M_PI, M_PI, 21 ); + ++m_missedVeloHitsEta[mcH->mcParticle()->momentum().Eta()]; + ++m_missedVeloHitsPhi[mcH->mcParticle()->momentum().phi()]; if ( nullptr != mcH->mcParticle()->originVertex()->mother() ) { if ( abs( mcH->mcParticle()->particleID().pid() ) == 321 && abs( mcH->mcParticle()->originVertex()->mother()->particleID().pid() ) == 333 ) { - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits not found Velo reco - phi - phikk", -M_PI, M_PI, 21 ); + ++m_missedVeloHitsPhikk[mcH->mcParticle()->momentum().phi()]; } if ( abs( mcH->mcParticle()->particleID().pid() ) == 11 && abs( mcH->mcParticle()->originVertex()->mother()->particleID().pid() ) == 22 ) { - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits not found Velo reco - phi - gammaee", -M_PI, M_PI, 21 ); + ++m_missedVeloHitsPhigee[mcH->mcParticle()->momentum().phi()]; } } - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits not found Velo reco - r", 5.1, 50, 50 ); - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits not found Velo reco - r<7mm", 5.1, 7, 14 ); - plot1D( mcH->mcParticle()->p() / 1000.0, "# MCHits not found Velo reco - p", 0, 100, 120 ); - plot1D( mcH->mcParticle()->pt() / 1000.0, "# MCHits not found Velo reco - pt", 0, 10, 120 ); - plot1D( mcH->midPoint().z(), "# MCHits not found Velo reco - z", -300, 800, 90 ); - plot1D( mcH->sensDetID() / 4, "# MCHits not found Velo reco - module", -0.5, 51.5, 52 ); + auto r = sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ); + ++m_missedVeloHitsR[r]; + ++m_missedVeloHitsR7[r]; + ++m_missedVeloHitsP[mcH->mcParticle()->p() / 1000.0]; + ++m_missedVeloHitsPt[mcH->mcParticle()->pt() / 1000.0]; + ++m_missedVeloHitsZ[mcH->midPoint().z()]; + ++m_missedVeloHitsModule[mcH->sensDetID() / 4]; } // plots for reconstructed MCHits from Velo reconstructible track for ( auto& mcH : hitsReconstructibleVeloReco ) { - plot1D( mcH->mcParticle()->momentum().Eta(), "# MCHits reconstructible Velo reco - eta", -5.5, 5.5, 70 ); - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits reconstructible Velo reco - phi", -M_PI, M_PI, 21 ); + ++m_recVeloHitsEta[mcH->mcParticle()->momentum().Eta()]; + ++m_recVeloHitsPhi[mcH->mcParticle()->momentum().phi()]; if ( nullptr != mcH->mcParticle()->originVertex()->mother() ) { if ( abs( mcH->mcParticle()->particleID().pid() ) == 321 && abs( mcH->mcParticle()->originVertex()->mother()->particleID().pid() ) == 333 ) { - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits reconstructible Velo reco - phi - phikk", -M_PI, M_PI, - 21 ); + ++m_recVeloHitsPhikk[mcH->mcParticle()->momentum().phi()]; } if ( abs( mcH->mcParticle()->particleID().pid() ) == 11 && abs( mcH->mcParticle()->originVertex()->mother()->particleID().pid() ) == 22 ) { - plot1D( mcH->mcParticle()->momentum().phi(), "# MCHits reconstructible Velo reco - phi - gammaee", -M_PI, M_PI, - 21 ); + ++m_recVeloHitsPhigee[mcH->mcParticle()->momentum().phi()]; } } - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits reconstructible Velo reco - r", 5.1, 50, 50 ); - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits reconstructible Velo reco - r<7mm", 5.1, 7, 14 ); - plot1D( mcH->mcParticle()->p() / 1000.0, "# MCHits reconstructible Velo reco - p", 0, 100, 120 ); - plot1D( mcH->mcParticle()->pt() / 1000.0, "# MCHits reconstructible Velo reco - pt", 0, 10, 120 ); - plot1D( mcH->midPoint().z(), "# MCHits reconstructible Velo reco - z", -300, 800, 90 ); - plot1D( mcH->sensDetID() / 4, "# MCHits reconstructible Velo reco - module", -0.5, 51.5, 52 ); + auto r = sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ); + ++m_recVeloHitsR[r]; + ++m_recVeloHitsR7[r]; + ++m_recVeloHitsP[mcH->mcParticle()->p() / 1000.0]; + ++m_recVeloHitsPt[mcH->mcParticle()->pt() / 1000.0]; + ++m_recVeloHitsZ[mcH->midPoint().z()]; + ++m_recVeloHitsModule[mcH->sensDetID() / 4]; } // plots for reconstructed MCHits from muons for ( auto& hitsInSensor : hitsReconstructibleInSensorsMu ) { - plot1D( hitsInSensor.size(), "# muons r<7mm", -0.5, 10.5, 11 ); + ++m_muonR7[hitsInSensor.size()]; if ( hitsInSensor.size() > 1 ) { auto sensor = hitsInSensor[0]->sensDetID(); double efficiency; @@ -441,8 +840,8 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra } else { efficiency = 0.0; } - profile1D( hitsInSensor.size(), efficiency, "Efficiency vs # #mu (at least 2 #mu)", -0.5, 10.5, 11 ); - plot1D( sensor, "Sensors w at least 2 #mu", -0.5, 207.5, 208 ); + m_muonEff[hitsInSensor.size()] += efficiency; + ++m_muonSemsorsMu[sensor]; } } @@ -450,7 +849,7 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra for ( auto& clus : clusters ) { m_num_pix_clu += clus.pixels().size(); std::set associated_hits; - auto strModNum = std::to_string( clus.channelID().module() ); + auto module = clus.channelID().module(); for ( auto& channelID : clus.pixels() ) { for ( auto& mcHit : ( *MCHitForchannelId.find( channelID.channelID() ) ).second ) { if ( channelIdForMCHit.find( mcHit ) != channelIdForMCHit.end() ) { associated_hits.insert( mcHit ); } @@ -463,36 +862,35 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra if ( abs( x_dist ) < 0.1 ) { m_residual_x += x_dist; } if ( abs( y_dist ) < 0.1 ) { m_residual_y += y_dist; } - plot1D( x_dist, "Residuals along x [mm]", -0.2, 0.2, 200 ); - plot1D( y_dist, "Residuals along y [mm]", -0.2, 0.2, 200 ); - plot1D( x_dist, "Residuals along x [mm] - module" + strModNum, -0.2, 0.2, 200 ); - plot1D( y_dist, "Residuals along y [mm] - module" + strModNum, -0.2, 0.2, 200 ); + ++m_resX[x_dist]; + ++m_resY[y_dist]; + ++m_resXPerModule[module][x_dist]; + ++m_resYPerModule[module][y_dist]; if ( trackInfo.hasVelo( mcHit->mcParticle() ) ) { // make residuals for VELO reconstructible tracks - plot1D( x_dist, "Residuals along x - VELO reco [mm]", -0.2, 0.2, 200 ); - plot1D( y_dist, "Residuals along y - VELO reco [mm]", -0.2, 0.2, 200 ); + ++m_resXVelo[x_dist]; + ++m_resYVelo[y_dist]; } if ( trackInfo.hasT( mcHit->mcParticle() ) && trackInfo.hasVelo( mcHit->mcParticle() ) ) { // make residuals for - // long tracks - plot1D( x_dist, "Residuals along x - long reco [mm]", -0.2, 0.2, 200 ); - plot1D( y_dist, "Residuals along y - long reco [mm]", -0.2, 0.2, 200 ); + ++m_resXLong[x_dist]; + ++m_resYLong[y_dist]; } if ( clus.channelID().col() == LHCb::Detector::VPChannelID::ColumnID{255} || clus.channelID().col() == LHCb::Detector::VPChannelID::ColumnID{256} || clus.channelID().col() == LHCb::Detector::VPChannelID::ColumnID{511} || clus.channelID().col() == LHCb::Detector::VPChannelID::ColumnID{512} ) { // make residuals for long pixels - plot1D( x_dist, "Residuals along x - long pixels [mm]", -0.2, 0.2, 200 ); - plot1D( y_dist, "Residuals along y - long pixels [mm]", -0.2, 0.2, 200 ); + ++m_resXLongPixel[x_dist]; + ++m_resYLongPixel[y_dist]; } if ( clus.channelID().col() < LHCb::Detector::VPChannelID::ColumnID{3} || clus.channelID().col() > LHCb::Detector::VPChannelID::ColumnID{764} || clus.channelID().row() < LHCb::Detector::VPChannelID::RowID{3} || clus.channelID().row() > LHCb::Detector::VPChannelID::RowID{252} ) { // make residuals for sensor edges - plot1D( x_dist, "Residuals along x - sensor edges [mm]", -0.2, 0.2, 200 ); - plot1D( y_dist, "Residuals along y - sensor edges [mm]", -0.2, 0.2, 200 ); + ++m_resXEdge[x_dist]; + ++m_resYEdge[y_dist]; } } } @@ -520,8 +918,8 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra } } if ( gotdist ) { - plot2D( x_dist, y_dist, "Non reconstructed MCHit - closest cluster", -0.99, 0.99, -0.99, 0.99, 36, 36 ); - plot1D( dist, "Non reconstructed MCHit - closest cluster - 1D", 0, 0.66, 48 ); + ++m_noRecClosestCluter2D[{x_dist, y_dist}]; + ++m_noRecClosestCluter1D[dist]; } gotdist = false; @@ -540,9 +938,7 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra } } - if ( gotdist && distmcHmcH > 0 ) { - plot1D( distmcHmcH, "Non reconstructed MCHit - closest reconstructible MCHit", 0, 0.66, 48 ); - } + if ( gotdist && distmcHmcH > 0 ) { ++m_noRecClosestRecHit[distmcHmcH]; } } } @@ -564,8 +960,8 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra } } if ( gotdist && distmcHmcH > 0 ) { - plot1D( distmcHmcH, "Reconstructible MCHit - closest reconstructible MCHit", 0, 0.66, 48 ); - plot1D( distmcHmcH, "Reconstructible MCHit - closest reconstructible MCHit - extended", 0, 52.8, 3840 ); + ++m_recClosestRecHit[distmcHmcH]; + ++m_recClosestRecHitExt[distmcHmcH]; } } } @@ -597,9 +993,7 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra } } } - if ( gotdist && dist > 0 ) { - plot1D( dist, "Non reconstructed MCHit - closest reconstructible MCHit - phikk", 0, 0.66, 48 ); - } + if ( gotdist && dist > 0 ) { ++m_noRecClosestRecHitphikk[dist]; } } } } @@ -629,9 +1023,7 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra } } } - if ( gotdist && distmcHmcH > 0 ) { - plot1D( distmcHmcH, "Reconstructible MCHit - closest reconstructible MCHit - phikk", 0, 0.66, 48 ); - } + if ( gotdist && distmcHmcH > 0 ) { ++m_recClosestRecHitphikk[distmcHmcH]; } } } } @@ -664,9 +1056,7 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra } } } - if ( gotdist && dist > 0 ) { - plot1D( dist, "Non reconstructed MCHit - closest reconstructible MCHit - gammaee", 0, 0.66, 48 ); - } + if ( gotdist && dist > 0 ) { ++m_noRecClosestRecHitgee[dist]; } } } } @@ -696,26 +1086,24 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra } } } - if ( gotdist && distmcHmcH > 0 ) { - plot1D( distmcHmcH, "Reconstructible MCHit - closest reconstructible MCHit - gammaee", 0, 0.66, 48 ); - } + if ( gotdist && distmcHmcH > 0 ) { ++m_recClosestRecHitgee[distmcHmcH]; } } } } } // plot distribution number of clusters - m_num_clusters += clusters.size(); - plot1D( clusters.size(), "# clusters", 0, 10000, 10000 ); - plot1D( unmatched_clusters, "# unmatched clusters", 0, 10000, 100 ); + ++m_num_clusters += clusters.size(); + ++m_nClusters[clusters.size()]; + ++m_nUnmatchedClusters[unmatched_clusters]; for ( auto& clus : clusters ) { - auto strModNum = std::to_string( clus.channelID().module() ); - plot1D( sqrt( clus.x() * clus.x() + clus.y() * clus.y() ), "# cluster distribution - r", 5.1, 50, 50 ); - plot1D( sqrt( clus.x() * clus.x() + clus.y() * clus.y() ), "# cluster distribution - r - module" + strModNum, 5.1, - 50, 50 ); - plot1D( clus.channelID().module(), "# cluster distribution - module", -0.5, 51.5, 52 ); - plot1D( clus.pixels().size(), "# pixel per cluster distribution", -0.5, 50.5, 51 ); - plot1D( clus.pixels().size(), "# pixel per cluster distribution - module" + strModNum, -0.5, 50.5, 51 ); + auto module = clus.channelID().module(); + auto r = sqrt( clus.x() * clus.x() + clus.y() * clus.y() ); + ++m_nClustersR[r]; + ++m_nClustersRPerModule[module][r]; + ++m_nClustersModule[module]; + ++m_nPixelPerCluster[clus.pixels().size()]; + ++m_nPixelPerClusterPerModule[module][clus.pixels().size()]; } // purity plots @@ -733,12 +1121,13 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra std::sort( ids_hit.begin(), ids_hit.end() ); std::vector common; set_intersection( ids_clu.begin(), ids_clu.end(), ids_hit.begin(), ids_hit.end(), back_inserter( common ) ); - plot1D( ( (double)common.size() ) / ids_hit.size(), "Purity distribution", 0, 1.1, 110 ); - m_purity += ( (double)common.size() ) / ids_hit.size(); - plot1D( ids_hit.size(), "# pixel distribution - seen", 0, 100, 20 ); - profile1D( ids_hit.size(), ( (double)common.size() ) / ids_hit.size(), "Purity vs # pixels", 0, 100, 20 ); + auto pur = ( (double)common.size() ) / ids_hit.size(); + ++m_purityHisto[pur]; + m_purity += pur; + ++m_nPixels[ids_hit.size()]; + m_purityPixels[ids_hit.size()] += pur; } - plot1D( ids_clu.size(), "Number of pixels in cluster", 0, 32, 33 ); + ++m_nPixelsInCluster[ids_clu.size()]; } // pixel stat plots @@ -753,8 +1142,8 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra } } } - plot1D( associated_clusters.size(), "# associated cluster to hit", -0.5, 10.5, 11 ); - plot1D( channelIdForMCHit.find( mcHit )->second.size(), "# pixel distribution", 0, 100, 20 ); + ++m_cluterTohit[associated_clusters.size()]; + ++m_pixels[channelIdForMCHit.find( mcHit )->second.size()]; m_efficiency += !associated_clusters.empty(); m_num_pix_hit += ( channelIdForMCHit.find( mcHit )->second.size() ); } @@ -764,76 +1153,55 @@ void VPClusterEfficiency::operator()( const LHCb::RawBank::View& ra for ( const auto& mcH : hitsReconstructibleInSensor ) { // reconstructible MCHits const auto sensor = mcH->sensDetID(); - plot2D( mcH->midPoint().x(), mcH->midPoint().y(), - "Reconstructible MCHit x,y VP sensor " + std::to_string( sensor ), -60, 60, -60, 60, 240, 240 ); - plot2D( mcH->midPoint().x(), mcH->midPoint().y(), "Reconstructible MCHit x,y VP all sensors ", -60, 60, -60, 60, - 240, 240 ); + ++m_recMCHitPerSensor[sensor][{mcH->midPoint().x(), mcH->midPoint().y()}]; + ++m_recMCHit[{mcH->midPoint().x(), mcH->midPoint().y()}]; if ( trackInfo.hasVelo( mcH->mcParticle() ) ) { - plot2D( mcH->midPoint().x(), mcH->midPoint().y(), - "Reconstructible MCHit x,y VP sensor - VELO " + std::to_string( sensor ), -60, 60, -60, 60, 240, 240 ); - plot2D( mcH->midPoint().x(), mcH->midPoint().y(), "Reconstructible MCHit x,y VP all sensors - VELO ", -60, 60, - -60, 60, 240, 240 ); + ++m_recMCHitVeloPerSensor[sensor][{mcH->midPoint().x(), mcH->midPoint().y()}]; + ++m_recMCHitVelo[{mcH->midPoint().x(), mcH->midPoint().y()}]; } } } for ( const auto& hitsMissedInSensor : hitsMissedInSensors ) { - - std::string strSenNum; - if ( hitsMissedInSensor.size() != 0 ) strSenNum = std::to_string( hitsMissedInSensor[0]->sensDetID() ); - + int sensor = -1; + if ( hitsMissedInSensor.size() != 0 ) sensor = hitsMissedInSensor[0]->sensDetID(); for ( auto& mcH : hitsMissedInSensor ) { // missed MCHits - plot2D( mcH->midPoint().x(), mcH->midPoint().y(), "Missed MCHit x,y VP sensor " + strSenNum, -60, 60, -60, 60, - 240, 240 ); - plot2D( mcH->midPoint().x(), mcH->midPoint().y(), "Missed MCHit x,y VP all sensors ", -60, 60, -60, 60, 240, - 240 ); + ++m_missedMCHitPerSensor[sensor][{mcH->midPoint().x(), mcH->midPoint().y()}]; + ++m_missedMCHit[{mcH->midPoint().x(), mcH->midPoint().y()}]; if ( trackInfo.hasVelo( mcH->mcParticle() ) ) { - plot2D( mcH->midPoint().x(), mcH->midPoint().y(), "Missed MCHit x,y VP sensor - VELO " + strSenNum, -60, 60, - -60, 60, 240, 240 ); - plot2D( mcH->midPoint().x(), mcH->midPoint().y(), "Missed MCHit x,y VP all sensors - VELO ", -60, 60, -60, 60, - 240, 240 ); + ++m_missedMCHitVeloPerSensor[sensor][{mcH->midPoint().x(), mcH->midPoint().y()}]; + ++m_missedMCHitVelo[{mcH->midPoint().x(), mcH->midPoint().y()}]; } } } for ( const auto& hitsReconstructibleInSensor : hitsReconstructibleInSensors ) { - - std::string strModNum; - if ( hitsReconstructibleInSensor.size() != 0 ) - strModNum = std::to_string( hitsReconstructibleInSensor[0]->sensDetID() / 4 ); - + int module = -1; + if ( hitsReconstructibleInSensor.size() != 0 ) module = hitsReconstructibleInSensor[0]->sensDetID() / 4; for ( const auto& mcH : hitsReconstructibleInSensor ) { // reconstructible MCHits - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits reconstructible - r - module" + strModNum, 5.1, 50, 50 ); - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits reconstructible - r<7mm - module" + strModNum, 5.1, 7, 14 ); + auto r = sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ); + ++m_recHitsRPerModule[module][r]; + ++m_recHitsR7PerModule[module][r]; if ( trackInfo.hasVelo( mcH->mcParticle() ) ) { - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits reconstructible Velo reco - r - module" + strModNum, 5.1, 50, 50 ); - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits reconstructible Velo reco - r<7mm - module" + strModNum, 5.1, 7, 14 ); + ++m_recHitsVeloRPerModule[module][r]; + ++m_recHitsVeloR7PerModule[module][r]; } } } for ( const auto& hitsMissedInSensor : hitsMissedInSensors ) { - - std::string strModNum; - if ( hitsMissedInSensor.size() != 0 ) strModNum = std::to_string( hitsMissedInSensor[0]->sensDetID() / 4 ); - + int module = -1; + if ( hitsMissedInSensor.size() != 0 ) module = hitsMissedInSensor[0]->sensDetID() / 4; for ( auto& mcH : hitsMissedInSensor ) { // missed MCHits - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits not found - r - module" + strModNum, 5.1, 50, 50 ); - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits not found - r<7mm - module" + strModNum, 5.1, 7, 14 ); + auto r = sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ); + ++m_missedHitsRPerModule[module][r]; + ++m_missedHitsR7PerModule[module][r]; if ( trackInfo.hasVelo( mcH->mcParticle() ) ) { - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits not found Velo reco - r - module" + strModNum, 5.1, 50, 50 ); - plot1D( sqrt( mcH->midPoint().x() * mcH->midPoint().x() + mcH->midPoint().y() * mcH->midPoint().y() ), - "# MCHits not found Velo reco - r<7mm - module" + strModNum, 5.1, 7, 14 ); + ++m_missedHitsVeloRPerModule[module][r]; + ++m_missedHitsVeloR7PerModule[module][r]; } } } diff --git a/Rec/RecAlgs/CMakeLists.txt b/Rec/RecAlgs/CMakeLists.txt index 3365f446ec3fb55758f4ba50f40ce7bd0d305ad4..93fdeb5ec75a7051012b27b5e77aaadca240cd54 100644 --- a/Rec/RecAlgs/CMakeLists.txt +++ b/Rec/RecAlgs/CMakeLists.txt @@ -16,9 +16,6 @@ Rec/RecAlgs gaudi_add_module(RecAlgs SOURCES src/AddToProcStatus.cpp - src/EventTimeMonitor.cpp - src/ProcStatAbortMoni.cpp - src/RecProcessingTimeMoni.cpp src/RecSummaryMaker.cpp src/TimingTuple.cpp src/UniqueIDGeneratorAlg.cpp diff --git a/Rec/RecAlgs/src/EventTimeMonitor.cpp b/Rec/RecAlgs/src/EventTimeMonitor.cpp deleted file mode 100644 index de18ba62bca995dd25b1b619afedd5adff324503..0000000000000000000000000000000000000000 --- a/Rec/RecAlgs/src/EventTimeMonitor.cpp +++ /dev/null @@ -1,96 +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 "AIDA/IHistogram1D.h" -#include "Event/ODIN.h" // event & run number -#include "GaudiAlg/GaudiHistoAlg.h" -#include "LHCbAlgs/Consumer.h" - -namespace { - //@FIXME: what about leapyears? - static const std::array month_offsets = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; -} // namespace - -//----------------------------------------------------------------------------- -// Implementation file for class : EventTimeMonitor -// -// 2012-04-19 : Patrick Koppenburg -//----------------------------------------------------------------------------- - -/** @class EventTimeMonitor EventTimeMonitor.h - * - * Creates histograms of event time - * - * @author Patrick Koppenburg - * @date 2012-04-19 - */ - -class EventTimeMonitor final : public LHCb::Algorithm::Consumer> { -public: - /// Standard constructor - EventTimeMonitor( const std::string& name, ISvcLocator* pSvcLocator ); - - StatusCode initialize() override; ///< Algorithm initialization - void operator()( const LHCb::ODIN& ) const override; ///< Algorithm execution - -private: - /// the histogram definition (as property) - Gaudi::Property m_histoS{ - this, "SecondsPlot", {"GPS Seconds", 0, 3600., 3600}, "The parameters of 'delta memory' histogram"}; - Gaudi::Histo1DDef m_histoH{"GPS Hours", -0.5, 23.5, 24}; // the histogram definition (as property) - Gaudi::Histo1DDef m_histoD{"GPS Days", 0.5, 365.5, 366}; // the histogram definition (as property) - Gaudi::Histo1DDef m_histoY{"GPS Year", 2008.5, 2028.5, 20}; // the histogram definition (as property) - mutable AIDA::IHistogram1D* m_plotS = nullptr; // the histogram of seconds - mutable AIDA::IHistogram1D* m_plotH = nullptr; // the histogram of hours - mutable AIDA::IHistogram1D* m_plotD = nullptr; // the histogram of day of year - mutable AIDA::IHistogram1D* m_plotY = nullptr; // the histogram of year -}; - -// Declaration of the Algorithm Factory -DECLARE_COMPONENT( EventTimeMonitor ) - -//============================================================================= -// Standard constructor, initializes variables -//============================================================================= -EventTimeMonitor::EventTimeMonitor( const std::string& name, ISvcLocator* pSvcLocator ) - : Consumer( name, pSvcLocator, KeyValue{"Input", LHCb::ODINLocation::Default} ) {} -//============================================================================= -// Initialization -//============================================================================= -StatusCode EventTimeMonitor::initialize() { - StatusCode sc = Consumer::initialize(); // must be executed first - if ( sc.isFailure() ) return sc; // error printed already by GaudiHistoAlg - - if ( msgLevel( MSG::DEBUG ) ) debug() << "==> Initialize" << endmsg; - if ( produceHistos() ) { - m_plotS = book( m_histoS ); - m_plotH = book( m_histoH ); - m_plotD = book( m_histoD ); - m_plotY = book( m_histoY ); - } - return sc; -} - -//============================================================================= -// Main execution -//============================================================================= -void EventTimeMonitor::operator()( const LHCb::ODIN& odin ) const { - - if ( produceHistos() ) { - const Gaudi::Time gtime = odin.eventTime(); - m_plotY->fill( gtime.year( false ) ); - m_plotD->fill( month_offsets[gtime.month( false )] + gtime.day( false ) ); - m_plotH->fill( gtime.hour( false ) ); - m_plotS->fill( 60 * gtime.minute( false ) + gtime.second( false ) + gtime.nsecond() / 1000000000. ); - } -} - -//============================================================================= diff --git a/Rec/RecAlgs/src/ProcStatAbortMoni.cpp b/Rec/RecAlgs/src/ProcStatAbortMoni.cpp deleted file mode 100644 index e2be17a08e9e53edddbca88fb53c34ff155b1151..0000000000000000000000000000000000000000 --- a/Rec/RecAlgs/src/ProcStatAbortMoni.cpp +++ /dev/null @@ -1,95 +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 "AIDA/IProfile1D.h" -#include "Event/ProcStatus.h" -#include "GaudiAlg/GaudiHistoAlg.h" -#include "GaudiUtils/HistoLabels.h" -#include "LHCbAlgs/Consumer.h" - -//----------------------------------------------------------------------------- -// Implementation file for class : ProcStatAbortMoni -// -// 2010-07-16 : Chris Jones -//----------------------------------------------------------------------------- - -/** @class ProcStatAbortMoni ProcStatAbortMoni.h - * - * Monitor for abort rates in ProcStat - * - * @author Chris Jones - * @date 2010-07-16 - */ - -class ProcStatAbortMoni final - : public LHCb::Algorithm::Consumer> { -public: - /// Standard constructor - ProcStatAbortMoni( const std::string& name, ISvcLocator* pSvcLocator ); - - StatusCode initialize() override; ///< Algorithm initialization - void operator()( const LHCb::ProcStatus& ) const override; ///< Algorithm execution - -private: - /// List of subsystems - Gaudi::Property> m_subSystems{ - this, - "SubSystems", - {"Overall", "Hlt", "VELO", "TT", "IT", "OT", "Tracking", "Vertex", "RICH", "CALO", "MUON", "PROTO"}}; - - /// cache the histogram pointer - AIDA::IProfile1D* m_h = nullptr; -}; - -// Declaration of the Algorithm Factory -DECLARE_COMPONENT( ProcStatAbortMoni ) - -//============================================================================= -// Standard constructor, initializes variables -//============================================================================= -ProcStatAbortMoni::ProcStatAbortMoni( const std::string& name, ISvcLocator* pSvcLocator ) - : Consumer( name, pSvcLocator, {"ProcStatusLocation", LHCb::ProcStatusLocation::Default} ) {} - -//============================================================================= -// Initialization -//============================================================================= -StatusCode ProcStatAbortMoni::initialize() { - StatusCode sc = Consumer::initialize(); - if ( sc.isFailure() ) return sc; - - // book the histo - m_h = bookProfile1D( "aborts", "Processing Abort Rates (%)", -0.5, m_subSystems.size() - 0.5, m_subSystems.size() ); - // Set the bin labels - const bool ok = Gaudi::Utils::Histos::setBinLabels( m_h, m_subSystems ); - - // return - return ok ? sc : StatusCode::FAILURE; -} - -//============================================================================= -// Main execution -//============================================================================= -void ProcStatAbortMoni::operator()( const LHCb::ProcStatus& proc ) const { - // Loop over sub-systems and fill plot - int index( 0 ); - for ( const auto& sys : m_subSystems ) { - ++index; - if ( sys == "Overall" ) { - m_h->fill( index, proc.aborted() ? 100.0 : 0.0 ); - } else { - m_h->fill( index, proc.subSystemAbort( sys ) ? 100.0 : 0.0 ); - } - } - // Debug printout if aborted - if ( proc.aborted() && msgLevel( MSG::DEBUG ) ) { debug() << proc << endmsg; } -} - -//============================================================================= diff --git a/Rec/RecAlgs/src/RecProcessingTimeMoni.cpp b/Rec/RecAlgs/src/RecProcessingTimeMoni.cpp deleted file mode 100644 index d1618d53baa2cef5ca182a50b87033721c69bf30..0000000000000000000000000000000000000000 --- a/Rec/RecAlgs/src/RecProcessingTimeMoni.cpp +++ /dev/null @@ -1,100 +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 "AIDA/IHistogram1D.h" -#include "GaudiAlg/GaudiHistoAlg.h" -#include -//----------------------------------------------------------------------------- -// Implementation file for class : RecProcessingTimeMoni -// -// 2010-07-15 : Chris Jones -//----------------------------------------------------------------------------- - -/** @class RecProcessingTimeMoni RecProcessingTimeMoni.h - * - * Simple monitor making basic processing time plots for the Reconstruction - * - * @author Chris Jones - * @date 2010-07-15 - */ - -class RecProcessingTimeMoni final : public GaudiHistoAlg { - -public: - /// Standard constructor - RecProcessingTimeMoni( const std::string& name, ISvcLocator* pSvcLocator ); - - StatusCode initialize() override; ///< Algorithm initialization - StatusCode execute() override; ///< Algorithm execution - -private: - /// Definition of algorithm name list - typedef std::vector AlgorithmNames; - - AIDA::IHistogram1D* m_hist = nullptr; ///< Pointer to processing time histogram - - Gaudi::Property m_algNames{this, "Algorithms"}; ///< List of algorithm(s) to include in timing - Gaudi::Property m_logMaxTime{this, "LogMaxEventTime", - 8.0}; ///< Job Option for log10(maximum overall processing time) for plots - Gaudi::Property m_logMinTime{this, "LogMinEventTime", - -2.0}; ///< Job Option for log10(minimum overall processing time) for plots -}; - -// Declaration of the Algorithm Factory -DECLARE_COMPONENT( RecProcessingTimeMoni ) - -//============================================================================= -// Standard constructor, initializes variables -//============================================================================= -RecProcessingTimeMoni::RecProcessingTimeMoni( const std::string& name, ISvcLocator* pSvcLocator ) - : GaudiHistoAlg( name, pSvcLocator ) { - setProperty( "HistoTopDir", "Timing/" ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); -} - -//============================================================================= -// Initialization -//============================================================================= -StatusCode RecProcessingTimeMoni::initialize() { - StatusCode sc = GaudiHistoAlg::initialize(); - if ( sc.isFailure() ) return sc; - - // are we properly configured - if ( m_algNames.empty() ) { sc = Warning( "No algorithms to time !" ); } - - // book the histogram at initialization time - m_hist = book( "overallTime", "log10(Event Processing Time / ms)", m_logMinTime, m_logMaxTime, 100 ); - - // return - return sc; -} - -//============================================================================= -// Main execution -//============================================================================= -StatusCode RecProcessingTimeMoni::execute() { - - // Loop over algorithms to include in the timing and add them up - double time = std::accumulate( - m_algNames.begin(), m_algNames.end(), double{0}, [&, chrono = chronoSvc()]( double t, const auto& name ) { - const auto alg_time = chrono->chronoDelta( name + ":Execute", IChronoStatSvc::ELAPSED ) / 1000; - if ( msgLevel( MSG::VERBOSE ) ) verbose() << name << " " << alg_time << endmsg; - return t + alg_time; - } ); - - // only fill if algorithm(s) ran (time>0) - if ( time > 0 ) { - // Take the base 10 log of the time (helps show the large tails) - fill( m_hist, std::log10( time ), 1.0 ); - } - - return StatusCode::SUCCESS; -} - -//============================================================================= diff --git a/Tr/TrackCheckers/CMakeLists.txt b/Tr/TrackCheckers/CMakeLists.txt index 99f4f12b1e1e87537a9ddee8efd031809aedc0e0..7178287a8c2a87a02766f4e599c355cb988e6824 100644 --- a/Tr/TrackCheckers/CMakeLists.txt +++ b/Tr/TrackCheckers/CMakeLists.txt @@ -15,13 +15,9 @@ Tr/TrackCheckers gaudi_add_module(TrackCheckers SOURCES - src/ExtrapolatorChecker.cpp src/Map.cpp src/TrackCheckerBase.cpp - src/TrackCloneChecker.cpp - src/TrackEffChecker.cpp src/TrackIPResolutionCheckerNT.cpp - src/TriggerObjectsCompatibilityChecker.cpp src/TriggerObjectsCompatibilityProfileChecker.cpp src/TrackResChecker.cpp src/VertexChecker.cpp diff --git a/Tr/TrackCheckers/src/ExtrapolatorChecker.cpp b/Tr/TrackCheckers/src/ExtrapolatorChecker.cpp deleted file mode 100644 index 5e0e45c615ccafd991876b80bb937ebba7dcbc00..0000000000000000000000000000000000000000 --- a/Tr/TrackCheckers/src/ExtrapolatorChecker.cpp +++ /dev/null @@ -1,306 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2019 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 "Core/FloatComparison.h" -#include "DetDesc/DetectorElement.h" -#include "DetDesc/GenericConditionAccessorHolder.h" -#include "Event/MCHit.h" -#include "Event/MCParticle.h" -#include "Event/Track.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "LHCbAlgs/Consumer.h" -#include "Linker/LinkedFrom.h" -#include "TrackCheckerBase.h" - -using namespace Gaudi; -using namespace Gaudi::Units; -using namespace LHCb; - -/** @class ExtrapolatorChecker ExtrapolatorChecker.h - * - * This algorithm checks the performance of the TrackExtrapolators. This - * algorithm can be run in Boole or Brunel, since it relies only on the - * presence of MCHits. Algorithm description: - * @li create a State from a MCHit, - * @li extrapolate this State to the next MCHit on the track, - * @li compare extrapolated state with MCHit parameters, - * @li and make the corresponding resolution and pull plots. - * - * The plots are split up in: - * @li first MCHit (extrapolated from true vertex) - * @li from VELO MCHit to VELO MCHit - * @li from TT MCHit to TT MCHit - * @li from IT MCHit to IT MCHit - * @li from OT MCHit to OT MCHit - * @li through RICH1: from VELO MCHit to TT MCHit - * @li through magnet: from TT MCHit to IT/OT MCHit - * - * Note that the extrapolations are done starting from the @b entry @b point - * of the MCHit and extrapolating to the entry point of the next MCHit. A - * correction for the magnetic field is applied on the slopes (@e tx and @e ty) - * of the MCHit. This is required since the displacement vector of the MCHit - * is calculated from the difference between the exit and entry point, and - * therefore this displacement vector already accounts for half the effect of - * the magnetic field. This effect is subtracted in the correction method. - * - * @author Jeroen van Tilburg - * @date 2006-07-06 - */ - -class ExtrapolatorChecker : public LHCb::Algorithm::Consumer< - void( MCParticles const&, DetectorElement const& ), - LHCb::Algorithm::Traits::usesBaseAndConditions> { - -public: - /// Standard constructor - ExtrapolatorChecker( const std::string& name, ISvcLocator* pSvcLocator ); - - /// Algorithm execution - void operator()( const MCParticles&, DetectorElement const& ) const override; - -private: - /// Find the next MCHit belonging to the same MCParticle starting from a z-pos - int findNextHit( const LHCb::MCParticle* mcPart, const double zRec, LHCb::MCHit const*& closestHit, - std::string& detectorName ) const; - - /// Helper function to find the next MCHit - LHCb::MCHit const* findNextXxxHit( const LHCb::MCParticle* mcPart, const double zRec, - LHCb::LinksByKey const* links ) const; - - /// Get the q/p for a given MCHit - double qOverP( const LHCb::MCParticle* mcPart, const LHCb::MCHit* mcHit ) const; - - /// Correct slopes for magnetic field given an MCHit and a MCParticle - void correctSlopes( const LHCb::MCParticle* mcPart, const LHCb::MCHit* mcHit, double& tx, double& ty ) const; - - /// String of the available detectors. - Gaudi::Property> m_dets{this, "Detectors", {"IT", "OT", "TT", "Velo"}}; -}; - -DECLARE_COMPONENT( ExtrapolatorChecker ) - -//============================================================================= -// Standard constructor, initializes variables -//============================================================================= -ExtrapolatorChecker::ExtrapolatorChecker( const std::string& name, ISvcLocator* pSvcLocator ) - : Consumer{name, - pSvcLocator, - {KeyValue{"MCParticleLocation", MCParticleLocation::Default}, - KeyValue{"MCParticleInContainer", LHCb::MCParticleLocation::Default}}} {} - -//============================================================================= -// Main execution -//============================================================================= -void ExtrapolatorChecker::operator()( const MCParticles& particles, DetectorElement const& geometry ) const { - - // Loop over the MCParticles - for ( MCParticle const* mcPart : particles ) { - // Decide whether the MCParticle will be checked - if ( selected( mcPart ) ) { - - // Get the state at the vertex - State stateVtx; - idealStateCreator() - ->createStateVertex( mcPart, stateVtx ) - .ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - - // Find first MCHit - std::string detName; - MCHit const* mcHit = nullptr; - int detID = findNextHit( mcPart, stateVtx.z(), mcHit, detName ); - - // Get the entry point of the MCHit - XYZPoint entryP = mcHit->entry(); - - // Extrapolate through RF foil - extrapolator() - ->propagate( stateVtx, entryP.z(), geometry ) - .ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - TrackVector vec = stateVtx.stateVector(); - TrackSymMatrix cov = stateVtx.covariance(); - - // Correct tx and ty from the MCHit for the magnetic field - double tx, ty; - correctSlopes( mcPart, mcHit, tx, ty ); - - // Determine the resolutions - double dx = vec( 0 ) - entryP.x(); - double dy = vec( 1 ) - entryP.y(); - double dtx = vec( 2 ) - tx; - double dty = vec( 3 ) - ty; - double dqp = vec( 4 ) - qOverP( mcPart, mcHit ); - - // fill the histograms - plot1D( dx, 1, "X resolution at 1st meas", -0.5, 0.5, 100 ); - plot1D( dy, 2, "Y resolution at 1st meas", -0.5, 0.5, 100 ); - plot1D( dtx, 3, "Tx resolution at 1st meas", -0.01, 0.01, 100 ); - plot1D( dty, 4, "Ty resolution at 1st meas", -0.01, 0.01, 100 ); - plot1D( stateVtx.p() - mcHit->p(), 5, "dp at 1st meas", -0.3, 0.3, 100 ); - if ( !LHCb::essentiallyZero( cov( 0, 0 ) ) && !LHCb::essentiallyZero( cov( 1, 1 ) ) && - !LHCb::essentiallyZero( cov( 2, 2 ) ) && !LHCb::essentiallyZero( cov( 3, 3 ) ) ) { - plot1D( dx / sqrt( cov( 0, 0 ) ), 11, "X pull at 1st meas", -5., 5., 100 ); - plot1D( dy / sqrt( cov( 1, 1 ) ), 12, "Y pull at 1st meas", -5., 5., 100 ); - plot1D( dtx / sqrt( cov( 2, 2 ) ), 13, "Tx pull at 1st meas", -5., 5., 100 ); - plot1D( dty / sqrt( cov( 3, 3 ) ), 14, "Ty pull at 1st meas", -5., 5., 100 ); - } - if ( !LHCb::essentiallyZero( cov( 4, 4 ) ) ) - plot1D( dqp / sqrt( cov( 4, 4 ) ), 15, "q/p pull at 1st meas", -5., 5., 100 ); - - const Gaudi::TrackSymMatrix zeroCov; - State state; - state.setState( entryP.x(), entryP.y(), entryP.z(), tx, ty, qOverP( mcPart, mcHit ) ); - state.setCovariance( zeroCov ); - - detID = findNextHit( mcPart, state.z(), mcHit, detName ); - - while ( mcHit ) { - - entryP = mcHit->entry(); - double dz = entryP.z() - state.z(); - - // Extrapolate to next MCHit - extrapolator() - ->propagate( state, entryP.z(), geometry ) - .ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - TrackVector vec = state.stateVector(); - TrackSymMatrix cov = state.covariance(); - - // Correct tx and ty from the MCHit for the magnetic field - correctSlopes( mcPart, mcHit, tx, ty ); - - // Determine the resolutions - double dx = vec( 0 ) - entryP.x(); - double dy = vec( 1 ) - entryP.y(); - double dtx = vec( 2 ) - tx; - double dty = vec( 3 ) - ty; - double dqp = vec( 4 ) - qOverP( mcPart, mcHit ); - - // Define the ranges for the resolution plots - double xr = 0.04; - double tr = 0.001; - double pr = 0.3; - - // Determine the histogram title - int ID = 100 * detID; - std::string title = " at " + detName + " hits"; - if ( dz > 4. * m && ( detName == "OT" || detName == "IT" ) ) { - ID = 1100; - title = " after magnet extrapolation"; - xr = 5.; - tr = 0.005; - pr = 5.; - } else if ( dz > 1. * m && detName == "TT" ) { - ID = 1000; - title = " after RICH1 extrapolation"; - xr = 5.; - tr = 0.005; - pr = 10.; - } - - // fill the histograms - plot1D( dx, ID + 1, "X resolution" + title, -xr, xr, 100 ); - plot1D( dy, ID + 2, "Y resolution" + title, -xr, xr, 100 ); - plot1D( dtx, ID + 3, "Tx resolution" + title, -tr, tr, 100 ); - plot1D( dty, ID + 4, "Ty resolution" + title, -tr, tr, 100 ); - plot1D( state.p() - mcHit->p(), ID + 5, "dp" + title, -pr, pr, 100 ); - if ( !LHCb::essentiallyZero( cov( 0, 0 ) ) && !LHCb::essentiallyZero( cov( 1, 1 ) ) && - !LHCb::essentiallyZero( cov( 2, 2 ) ) && !LHCb::essentiallyZero( cov( 3, 3 ) ) ) { - plot1D( dx / sqrt( cov( 0, 0 ) ), ID + 11, "X pull" + title, -5., 5., 100 ); - plot1D( dy / sqrt( cov( 1, 1 ) ), ID + 12, "Y pull" + title, -5., 5., 100 ); - plot1D( dtx / sqrt( cov( 2, 2 ) ), ID + 13, "Tx pull" + title, -5., 5., 100 ); - plot1D( dty / sqrt( cov( 3, 3 ) ), ID + 14, "Ty pull" + title, -5., 5., 100 ); - } - if ( !LHCb::essentiallyZero( cov( 4, 4 ) ) ) - plot1D( dqp / sqrt( cov( 4, 4 ) ), ID + 15, "q/p pull" + title, -5., 5., 100 ); - - state.setState( entryP.x(), entryP.y(), entryP.z(), tx, ty, qOverP( mcPart, mcHit ) ); - state.setCovariance( zeroCov ); - - detID = findNextHit( mcPart, entryP.z(), mcHit, detName ); - } - } - } // End loop over MCParticles -} - -//============================================================================= -// Find the next MCHit starting from a given z -// looping over the hits in all the tracking detectors -//============================================================================= -int ExtrapolatorChecker::findNextHit( const MCParticle* mcPart, const double zRec, MCHit const*& closestHit, - std::string& detectorName ) const { - detectorName = "Not found!"; - int detID = 0; - closestHit = nullptr; - ; - double closestZ = 100000.; - - for ( auto itDets = m_dets.begin(); itDets != m_dets.end(); ++itDets ) { - auto links = SmartDataPtr{ - evtSvc(), LHCb::LinksByKey::linkerName( MCParticleLocation::Default + "2MC" + *itDets + "Hits" )}; - - MCHit const* tmpClosestHit = findNextXxxHit( mcPart, zRec, links ); - if ( tmpClosestHit && tmpClosestHit->entry().z() > zRec && tmpClosestHit->entry().z() < closestZ ) { - closestHit = tmpClosestHit; - closestZ = tmpClosestHit->entry().z(); - detectorName = *itDets; - detID = 1 + itDets - m_dets.begin(); - } - } - return detID; -} - -//============================================================================= -// Find the next MCHit of type Xxx starting from a given z -//============================================================================= -const MCHit* ExtrapolatorChecker::findNextXxxHit( const MCParticle* mcPart, const double zRec, - const LHCb::LinksByKey* links ) const { - // Retrieve MCParticle to MCHit linker tables - double closestZ = 100000.; - const MCHit* closestHit = nullptr; - for ( const auto& mcHit : LinkedFrom{links}.range( mcPart ) ) { - double zOfHit = mcHit.entry().z(); - // get the closest hit - if ( zOfHit > zRec + 0.1 && zOfHit < closestZ ) { - closestHit = &mcHit; - closestZ = zOfHit; - } - } - - return closestHit; -} - -//============================================================================= -// Determine q/p given an MCHit and a MCParticle -//============================================================================= -double ExtrapolatorChecker::qOverP( const MCParticle* mcPart, const MCHit* mcHit ) const { - double charge = ( mcPart->particleID().threeCharge() ) / 3.; - double momentum = mcPart->p(); - if ( mcHit != NULL && !LHCb::essentiallyZero( mcHit->p() ) ) momentum = mcHit->p(); - return ( momentum > TrackParameters::lowTolerance ) ? charge / momentum : 0.; -} - -//============================================================================= -// Correct slopes for magnetic field given an MCHit and a MCParticle -//============================================================================= -void ExtrapolatorChecker::correctSlopes( const MCParticle* mcPart, const MCHit* mcHit, double& tx, double& ty ) const { - // TODO: I hope this method can be removed as soon as the displacement vector - // in the MCHit is calculated in Gauss using the momentum direction of the - // *entry point*. (JvT: 27/10/2006). - - // Get magnetic field vector - Gaudi::XYZVector B; - fieldSvc()->fieldVector( mcHit->midPoint(), B ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - - // Calculate new displacement vector and tx,ty slopes - Gaudi::XYZVector d = mcHit->displacement(); - Gaudi::XYZVector dNew = d - ( 0.5 * d.R() * qOverP( mcPart, mcHit ) * d.Cross( B ) * eplus * c_light ); - tx = dNew.x() / dNew.z(); - ty = dNew.y() / dNew.z(); -} diff --git a/Tr/TrackCheckers/src/TrackCheckerBase.cpp b/Tr/TrackCheckers/src/TrackCheckerBase.cpp index be07efc3080ba88787c7b648a6f03f993a498864..c8bb7d113c49810d9cb7f00ebde3e7033e021d08 100644 --- a/Tr/TrackCheckers/src/TrackCheckerBase.cpp +++ b/Tr/TrackCheckers/src/TrackCheckerBase.cpp @@ -16,9 +16,8 @@ // Initialization. Check parameters //============================================================================= StatusCode TrackCheckerBase::initialize() { - if ( "" == histoTopDir() ) { setHistoTopDir( "Track/" ); } // Mandatory initialization of GaudiAlgorithm - return GaudiHistoAlg::initialize().andThen( [&] { + return Algorithm::initialize().andThen( [&] { const auto& theMap = TrackMaps::recDescription(); m_recCat = theMap.find( m_selectionCriteria )->second; } ); @@ -31,8 +30,7 @@ const LHCb::MCParticle* TrackCheckerBase::mcTruth( const LHCb::Track& track, con if ( !mcparticle ) { mcparticle = static_cast( mcParts.containedObject( mcPartKey ) ); if ( mcparticle && mcparticle->particleID().threeCharge() == 0 ) { - this->Warning( "Linker table for track contains pointer to particle with zero charge", StatusCode::SUCCESS, 0 ) - .ignore(); + this->warning() << "Linker table for track contains pointer to particle with zero charge" << endmsg; mcparticle = nullptr; } } diff --git a/Tr/TrackCheckers/src/TrackCheckerBase.h b/Tr/TrackCheckers/src/TrackCheckerBase.h index b85647dd92d38605694aba6e67554b041ac183e3..91736cb9b46c84ce4e8c17d0e914e884ac597c35 100644 --- a/Tr/TrackCheckers/src/TrackCheckerBase.h +++ b/Tr/TrackCheckers/src/TrackCheckerBase.h @@ -12,7 +12,7 @@ #pragma once // from Gaudi -#include "GaudiAlg/GaudiHistoAlg.h" +#include "Gaudi/Algorithm.h" // interfaces #include "GaudiKernel/IMagneticFieldSvc.h" @@ -36,11 +36,11 @@ * @date 7-5-2007 */ -class TrackCheckerBase : public GaudiHistoAlg { +class TrackCheckerBase : public Gaudi::Algorithm { public: /** Standard construtor */ - using GaudiHistoAlg::GaudiHistoAlg; + using Algorithm::Algorithm; /** Algorithm initialization */ StatusCode initialize() override; diff --git a/Tr/TrackCheckers/src/TrackCloneChecker.cpp b/Tr/TrackCheckers/src/TrackCloneChecker.cpp deleted file mode 100644 index 9361caa17473b83b4ea769dc6f013d82ece2405b..0000000000000000000000000000000000000000 --- a/Tr/TrackCheckers/src/TrackCloneChecker.cpp +++ /dev/null @@ -1,190 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2019 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 "LHCbAlgs/Consumer.h" -#include "TrackCheckerBase.h" - -#include "boost/format.hpp" - -#include - -//----------------------------------------------------------------------------- -// Implementation file for class : TrackCloneChecker -// -// 2007-09-13 : Chris Jones -//----------------------------------------------------------------------------- - -/** @class TrackCloneChecker TrackCloneChecker.h - * - * Produce some simple plots for the Clone linker information - * - * @author Chris Jones - * @date 2007-09-13 - */ - -class TrackCloneChecker : public LHCb::Algorithm::Consumer> { - -public: - /// Standard constructor - TrackCloneChecker( const std::string& name, ISvcLocator* pSvcLocator ); - - /** Algorithm execute */ - void operator()( const LHCb::Tracks&, const LHCb::MCParticles&, const LHCb::LinksByKey&, - const LHCb::LinksByKey& ) const override; - - StatusCode finalize() override; ///< Algorithm finalize - -private: - /** @class TrackTally TrackCloneChecker.h - * - * Counts track information for clones - * - * @author Chris Jones - * @date 2007-09-13 - */ - struct TrackTally { - unsigned long totalClones{0}; - unsigned long totalNonClones{0}; - unsigned long totalGhosts{0}; - unsigned long rejectedClones{0}; - unsigned long rejectedNonClones{0}; - unsigned long rejectedGhosts{0}; - /// Map for one tally object per track history type - typedef std::map Map; - }; - -private: - /// Get efficiency - inline std::pair getEff1( const double top, const double bot ) const { - return {( bot > 0 ? 100.0 * top / bot : 0 ), - ( bot > 0 ? 100.0 * sqrt( ( top / bot ) * ( 1. - top / bot ) / bot ) : 0 )}; - } - - /// Get efficiency - inline std::pair getEff2( const double top, const double bot ) const { - return {( bot > 0 ? top / bot : 0 ), ( bot > 0 ? sqrt( top ) / bot : 0 )}; - } - -private: - /// Summary map XXXX This is not thread safe. Should be rewritten using Gaudi counters - mutable TrackTally::Map m_trackMap; - - /// KL distance cut - Gaudi::Property m_klCut{this, "CloneCut", 5000}; - - /// Event count - mutable Gaudi::Accumulators::Counter<> m_nEvtsCounter{this, "Nb events"}; -}; - -// Declaration of the Algorithm Factory -DECLARE_COMPONENT( TrackCloneChecker ) - -TrackCloneChecker::TrackCloneChecker( const std::string& name, ISvcLocator* pSvcLocator ) - : Consumer( name, pSvcLocator, - {KeyValue{"TracksInContainer", LHCb::TrackLocation::Default}, - KeyValue{"MCParticleInContainer", LHCb::MCParticleLocation::Default}, - KeyValue{"LinkerInTable", "Link/" + LHCb::TrackLocation::Default}, - KeyValue{"CloneLinkerLocation", "Link/" + LHCb::TrackLocation::Default + "Clones"}} ) {} - -void TrackCloneChecker::operator()( const LHCb::Tracks& tracks, const LHCb::MCParticles& mcParts, - const LHCb::LinksByKey& links, const LHCb::LinksByKey& cloneLinks ) const { - ++m_nEvtsCounter; - - // Builds a lookup table to know which MC particle has clones - // It holds how many reconstructed tracks we have for each MCParticle, identified by key - std::map nbReconstructedTracks; - links.applyToAllLinks( - [&nbReconstructedTracks]( unsigned int, unsigned int mcPartKey, float ) { ++nbReconstructedTracks[mcPartKey]; } ); - - // loop over tracks - for ( const auto* track : tracks ) { - // MCP for main track - const LHCb::MCParticle* mcP = mcTruth( *track, mcParts, links ); - if ( !selected( mcP ) ) continue; - - // pick up the clone info for this track - cloneLinks.applyToLinks( - track->key(), [&tracks, &mcParts, &links, this, &mcP]( unsigned int, unsigned int cloneKey, float weight ) { - const LHCb::Track* cloneTrack = static_cast( tracks.containedObject( cloneKey ) ); - if ( mcP ) { - // MCP for clone track - const LHCb::MCParticle* mcP_clone = this->mcTruth( *cloneTrack, mcParts, links ); - const bool mcSel = ( mcP_clone ? selected( mcP_clone ) : false ); - - // log10(KLdistance) - const double logFLdist = log10( weight ); - // const bool cloneID = ( weight < m_klCut ); - - if ( mcP_clone && mcSel ) { - if ( mcP == mcP_clone ) { - plot1D( logFLdist, "KLDtrueClones", "Log10(KLDistance) | True Clones", -5, 13, 100 ); - } else { - plot1D( logFLdist, "KLDnotClones", "Log10(KLDistance) | Not Clones", -5, 13, 100 ); - } - } else if ( mcP_clone && !mcSel ) { - plot1D( logFLdist, "KLDrejMCPs", "Log10(KLDistance) | Rejected MCParticles", -5, 13, 100 ); - } else { - plot1D( logFLdist, "KLDghosts", "Log10(KLDistance) | Ghosts", -5, 13, 100 ); - } - } - } ); - - // clone ID - const bool cloneID = ( track->info( LHCb::Track::AdditionalInfo::CloneDist, 9e99 ) < m_klCut ); - - // real clone ? - const bool hasMCclones = ( nbReconstructedTracks.at( mcP->key() ) > 1 ); - - // tally object - TrackTally& tally = m_trackMap[track->history()]; - - // tallies - if ( mcP ) { - if ( hasMCclones ) ++tally.totalClones; - if ( !hasMCclones ) ++tally.totalNonClones; - if ( hasMCclones && cloneID ) ++tally.rejectedClones; - if ( !hasMCclones && cloneID ) ++tally.rejectedNonClones; - } else { - ++tally.totalGhosts; - if ( cloneID ) ++tally.rejectedGhosts; - } - - } // track loop -} - -StatusCode TrackCloneChecker::finalize() { - const std::string& lines = - "============================================================================================"; - always() << lines << endmsg; - always() << " Clone summary for '" << inputLocation<3>() << "' IDed clones with KLdist<" << m_klCut << endmsg; - always() << lines << endmsg; - - std::pair r1, r2, r3, r4, r5; - - always() << " Track type NonClones/Evt Clones/Evt ClonesID/% NonClonesID/% GhostsID/%" << endmsg; - for ( auto iM = m_trackMap.begin(); iM != m_trackMap.end(); ++iM ) { - const TrackTally& tally = iM->second; - r1 = getEff1( tally.rejectedClones, tally.totalClones / 2.0 ); - r2 = getEff1( tally.rejectedNonClones, tally.totalNonClones ); - r3 = getEff1( tally.rejectedGhosts, tally.totalGhosts ); - r4 = getEff2( tally.totalNonClones, m_nEvtsCounter.nEntries() ); - r5 = getEff2( tally.totalClones / 2.0, m_nEvtsCounter.nEntries() ); - always() << boost::format( "%15s %6.2f +-%5.2f %6.2f +-%5.2f %6.2f +-%5.2f %6.2f +-%5.2f %6.2f +-%5.2f" ) % - iM->first % r4.first % r4.second % r5.first % r5.second % r1.first % r1.second % r2.first % - r2.second % r3.first % r3.second - << endmsg; - } - - always() << lines << endmsg; - - return Consumer::finalize(); -} diff --git a/Tr/TrackCheckers/src/TrackEffChecker.cpp b/Tr/TrackCheckers/src/TrackEffChecker.cpp deleted file mode 100644 index b34a66781407339723a3dee3f5f30381104aada5..0000000000000000000000000000000000000000 --- a/Tr/TrackCheckers/src/TrackEffChecker.cpp +++ /dev/null @@ -1,743 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2019 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 "AIDA/IHistogram1D.h" -#include "Event/GhostTrackInfo.h" -#include "Event/MCParticle.h" -#include "Event/Track.h" -#include "GaudiKernel/SystemOfUnits.h" -#include "LHCbAlgs/Consumer.h" -#include "Map.h" -#include "TrackCheckerBase.h" - -/** @class TrackEffChecker TrackEffChecker.h - * - * Class for track monitoring - * @author M. Needham. - * @date 6-5-2007 - */ - -class TrackEffChecker : public LHCb::Algorithm::Consumer> { - -public: - /** Standard construtor */ - TrackEffChecker( const std::string& name, ISvcLocator* pSvcLocator ); - - /** Algorithm execute */ - void operator()( const LHCb::Tracks&, const LHCb::MCParticles&, const LHCb::LinksByKey&, - const LHCb::LinksByKey& ) const override; - - /** Algorithm finalize */ - StatusCode finalize() override; - -private: - Gaudi::Property m_requireLong{this, "RequireLongTrack", false}; - - void ghostInfo( const LHCb::Tracks&, const LHCb::MCParticles&, const LHCb::LinksByKey&, unsigned int ) const; - - void effInfo( const LHCb::Tracks&, const LHCb::MCParticles&, const LHCb::LinksByKey&, - const std::vector>& ) const; - - void plots( const std::string& type, const LHCb::Track* track ) const; - - void plots( const std::string& type, const LHCb::MCParticle* part ) const; - - double weightedMeasurementSum( const LHCb::Track* aTrack ) const; - - mutable Gaudi::Accumulators::AveragingCounter m_nTrackCounter{this, "nTrack"}; - mutable Gaudi::Accumulators::StatCounter<> m_nGhostCounter{this, "nGhost"}; - - mutable Gaudi::Accumulators::StatCounter<> m_nClone{this, "nClone"}; - mutable Gaudi::Accumulators::StatCounter<> m_nCloneG5{this, "nCloneG5"}; - mutable Gaudi::Accumulators::StatCounter<> m_nCloneB{this, "nCloneB"}; - mutable Gaudi::Accumulators::StatCounter<> m_nCloneG5B{this, "nCloneG5B"}; - mutable Gaudi::Accumulators::StatCounter<> m_nCloneBAll{this, "nCloneBAll"}; - mutable Gaudi::Accumulators::StatCounter<> m_nCloneG5BAll{this, "nCloneG5BAll"}; - mutable Gaudi::Accumulators::StatCounter<> m_nCloneKsL{this, "nCloneKsL"}; - mutable Gaudi::Accumulators::StatCounter<> m_nCloneG5KsL{this, "nCloneG5KsL"}; - mutable Gaudi::Accumulators::StatCounter<> m_nCloneIP{this, "nCloneIP"}; - mutable Gaudi::Accumulators::StatCounter<> m_nCloneG5IP{this, "nCloneG5IP"}; - - mutable Gaudi::Accumulators::StatCounter<> m_nToFind{this, "nToFind"}; - mutable Gaudi::Accumulators::StatCounter<> m_nFound{this, "nFound"}; - mutable Gaudi::Accumulators::StatCounter<> m_nToFindB{this, "nToFindB"}; - mutable Gaudi::Accumulators::StatCounter<> m_nFoundB{this, "nFoundB"}; - mutable Gaudi::Accumulators::StatCounter<> m_nToFindBAll{this, "nToFindBAll"}; - mutable Gaudi::Accumulators::StatCounter<> m_nFoundBAll{this, "nFoundBAll"}; - mutable Gaudi::Accumulators::StatCounter<> m_nToFindKsL{this, "nToFindKsL"}; - mutable Gaudi::Accumulators::StatCounter<> m_nFoundKsL{this, "nFoundKsL"}; - mutable Gaudi::Accumulators::StatCounter<> m_nToFindIP{this, "nToFindIP"}; - mutable Gaudi::Accumulators::StatCounter<> m_nFoundIP{this, "nFoundIP"}; - - mutable Gaudi::Accumulators::StatCounter<> m_nToFindG5{this, "nToFindG5"}; - mutable Gaudi::Accumulators::StatCounter<> m_nFoundG5{this, "nFoundG5"}; - mutable Gaudi::Accumulators::StatCounter<> m_nToFindG5B{this, "nToFindG5B"}; - mutable Gaudi::Accumulators::StatCounter<> m_nFoundG5B{this, "nFoundG5B"}; - mutable Gaudi::Accumulators::StatCounter<> m_nToFindG5BAll{this, "nToFindG5BAll"}; - mutable Gaudi::Accumulators::StatCounter<> m_nFoundG5BAll{this, "nFoundG5BAll"}; - mutable Gaudi::Accumulators::StatCounter<> m_nToFindG5KsL{this, "nToFindG5KsL"}; - mutable Gaudi::Accumulators::StatCounter<> m_nFoundG5KsL{this, "nFoundG5KsL"}; - mutable Gaudi::Accumulators::StatCounter<> m_nToFindG5IP{this, "nToFindG5IP"}; - mutable Gaudi::Accumulators::StatCounter<> m_nFoundG5IP{this, "nFoundG5IP"}; - - mutable Gaudi::Accumulators::SigmaCounter<> m_efficiency{this, "efficiency"}; - mutable Gaudi::Accumulators::SigmaCounter<> m_efficiencyG5{this, "efficiencyG5"}; - mutable Gaudi::Accumulators::SigmaCounter<> m_efficiencyB{this, "efficiencyB"}; - mutable Gaudi::Accumulators::SigmaCounter<> m_efficiencyG5B{this, "efficiencyG5B"}; - mutable Gaudi::Accumulators::SigmaCounter<> m_efficiencyBAll{this, "efficiencyBAll"}; - mutable Gaudi::Accumulators::SigmaCounter<> m_efficiencyG5BAll{this, "efficiencyG5BAll"}; - mutable Gaudi::Accumulators::SigmaCounter<> m_efficiencyKsL{this, "efficiencyKsL"}; - mutable Gaudi::Accumulators::SigmaCounter<> m_efficiencyG5KsL{this, "efficiencyG5KsL"}; - mutable Gaudi::Accumulators::SigmaCounter<> m_efficiencyIP{this, "efficiencyIP"}; - mutable Gaudi::Accumulators::SigmaCounter<> m_efficiencyG5IP{this, "efficiencyG5IP"}; -}; - -DECLARE_COMPONENT( TrackEffChecker ) - -TrackEffChecker::TrackEffChecker( const std::string& name, ISvcLocator* pSvcLocator ) - : Consumer( name, pSvcLocator, - {KeyValue{"TracksInContainer", LHCb::TrackLocation::Default}, - KeyValue{"MCParticleInContainer", LHCb::MCParticleLocation::Default}, - KeyValue{"LinkerInTable", "Link/" + LHCb::TrackLocation::Default}, - KeyValue{"AllLinksLocation", "Link/Pat/LHCbID"}} ) {} - -void TrackEffChecker::operator()( const LHCb::Tracks& tracks, const LHCb::MCParticles& mcParts, - const LHCb::LinksByKey& links, const LHCb::LinksByKey& allIds ) const { - - std::vector> linkedIds; - allIds.applyToAllLinks( [&linkedIds]( unsigned int srcKey, unsigned int mcPartKey, float ) { - while ( linkedIds.size() <= mcPartKey ) { - std::vector dum; - linkedIds.push_back( dum ); - } - linkedIds[mcPartKey].push_back( srcKey ); - } ); - - unsigned int nTracksInThisEvent = - std::count_if( tracks.begin(), tracks.end(), [requireLong = m_requireLong.value()]( const auto* track ) { - return track->type() == LHCb::Track::Types::Long || !requireLong; - } ); - - m_nTrackCounter += nTracksInThisEvent; - - // we want to count ghosts etc - ghostInfo( tracks, mcParts, links, nTracksInThisEvent ); - - // then we want to look at efficiencies - effInfo( tracks, mcParts, links, linkedIds ); -} - -void TrackEffChecker::ghostInfo( const LHCb::Tracks& tracks, const LHCb::MCParticles& mcParts, - const LHCb::LinksByKey& links, unsigned int nTracksInThisEvent ) const { - - unsigned int nGhost = 0; - std::string type = ""; - for ( const auto* track : tracks ) { - - if ( track->type() != LHCb::Track::Types::Long && m_requireLong.value() ) continue; - - const LHCb::MCParticle* particle = mcTruth( *track, mcParts, links ); - - splitByAlgorithm() == true ? type = Gaudi::Utils::toString( track->history() ) : type = "ALL"; - if ( !particle ) { - ++nGhost; - if ( fullDetail() == true ) { - plots( type + "/ghost", track ); - LHCb::GhostTrackInfo gInfo; - StatusCode sc = ghostClassification()->info( *track, gInfo ); - if ( sc.isSuccess() ) { - const LHCb::GhostTrackInfo::Classification& gtype = gInfo.classification(); - plot( gtype, "ghost classification", -0.5, 30.5, 31 ); - ++counter( Gaudi::Utils::toString( gtype ) ); - } - } - } else { - if ( fullDetail() == true ) plots( type + "/real", track ); - } - } - - // counter for ghost rate - m_nGhostCounter += nGhost; - - // plot the event ghost rate - if ( nTracksInThisEvent != 0 ) plot( nGhost / double( nTracksInThisEvent ), "ghost rate", -0.01, 1.01, 51 ); - if ( fullDetail() == true ) { - // ghost rate versus # interactions - long nVert = visPrimVertTool()->countVertices(); - if ( nTracksInThisEvent != 0 ) - plot2D( nVert, 100 * nGhost / double( nTracksInThisEvent ), "ghost rate vs visible", -0.5, 20.5, -0.5, 100.5, 21, - 101 ); - } -} - -void TrackEffChecker::effInfo( const LHCb::Tracks& tracks, const LHCb::MCParticles& mcParts, - const LHCb::LinksByKey& links, const std::vector>& linkedIds ) const { - - double efficiency = 0; - double efficiencyG5 = 0; - double efficiencyB = 0; - double efficiencyG5B = 0; - double efficiencyBAll = 0; - double efficiencyG5BAll = 0; - double efficiencyKsL = 0; - double efficiencyG5KsL = 0; - double efficiencyIP = 0; - double efficiencyG5IP = 0; - - unsigned int nToFind = 0u; - unsigned int nFound = 0u; - unsigned int nToFindG5 = 0u; - unsigned int nFoundG5 = 0u; - - unsigned int nToFindB = 0u; - unsigned int nFoundB = 0u; - unsigned int nToFindG5B = 0u; - unsigned int nFoundG5B = 0u; - - unsigned int nToFindBAll = 0u; - unsigned int nFoundBAll = 0u; - unsigned int nToFindG5BAll = 0u; - unsigned int nFoundG5BAll = 0u; - - unsigned int nToFindKsL = 0u; - unsigned int nFoundKsL = 0u; - unsigned int nToFindG5KsL = 0u; - unsigned int nFoundG5KsL = 0u; - - unsigned int nToFindIP = 0u; - unsigned int nFoundIP = 0u; - unsigned int nToFindG5IP = 0u; - unsigned int nFoundG5IP = 0u; - - // Build a map of reconstructed tracks for each MCParticle - std::map reconstructedTracks; - links.applyToAllLinks( - [&reconstructedTracks, &tracks]( unsigned int trackKey, unsigned int mcPartKey, float weight ) { - const LHCb::Track* track = static_cast( tracks.containedObject( trackKey ) ); - auto [it, inserted] = reconstructedTracks.try_emplace( mcPartKey, track, 0, weight ); - if ( !inserted ) ++it->second.clone; - } ); - - for ( const auto* mcPart : mcParts ) { - bool reconstructible = false; - bool reconstructibleB = false; - bool reconstructibleBAll = false; - bool reconstructibleKsL = false; - bool largeIP = false; - - reconstructible = selected( mcPart ); - reconstructibleBAll = reconstructible && bAncestorWithReconstructibleDaughters( mcPart ); - reconstructibleB = reconstructible && bAncestor( mcPart ); - reconstructibleKsL = reconstructible && ksLambdaAncestor( mcPart ); - - if ( reconstructible == true ) { - - double x = mcPart->originVertex()->position4vector().x(); - double y = mcPart->originVertex()->position4vector().y(); - double tx = mcPart->momentum().x() / mcPart->momentum().z(); - double ty = mcPart->momentum().y() / mcPart->momentum().z(); - - double IP2 = ( x * x + y * y ) - ( tx * x + ty * y ) * ( tx * x + ty * y ) / ( tx * tx + ty * ty ); - largeIP = reconstructibleB && IP2 > 1. && IP2 < 4.; // 1 to 2 mm - - double nTrue = 0; - - TrackCheckerBase::LinkInfo info = reconstructedTracks[mcPart->key()]; - - if ( info.track ) { - - std::vector ids; - if ( linkedIds.size() > (unsigned int)mcPart->key() ) { - for ( const auto id : linkedIds[mcPart->key()] ) { - LHCb::LHCbID temp; - temp.setDetectorType( LHCb::LHCbID::channelIDtype( id >> 28 ) ); // create LHCbId from int. Clumsy ! - temp.setID( id ); - ids.push_back( temp ); - } - } - - for ( const LHCb::LHCbID& id : ids ) { - if ( id.isVP() ) { - if ( info.track->type() == LHCb::Track::Types::Ttrack || - info.track->type() == LHCb::Track::Types::Downstream ) - continue; - nTrue += 1.; - } else if ( id.isVP() ) { - if ( info.track->type() == LHCb::Track::Types::Ttrack || - info.track->type() == LHCb::Track::Types::Downstream ) - continue; - nTrue += 1.; - } else if ( id.isUT() ) { - if ( info.track->type() != LHCb::Track::Types::Ttrack && info.track->type() != LHCb::Track::Types::Velo ) - nTrue += 1.; - } - } - } - - plots( "mcSelected", mcPart ); - ++nToFind; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) ++nToFindG5; - - if ( reconstructibleB ) { - ++nToFindB; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) ++nToFindG5B; - } - - if ( reconstructibleBAll ) { - ++nToFindBAll; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) ++nToFindG5BAll; - } - - if ( reconstructibleKsL ) { - ++nToFindKsL; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) ++nToFindG5KsL; - } - - if ( largeIP ) { - ++nToFindIP; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) ++nToFindG5IP; - } - - if ( info.track != 0 && ( info.track->type() == LHCb::Track::Types::Long || !m_requireLong.value() ) ) { - ++nFound; - efficiency += info.purity * info.track->lhcbIDs().size() / nTrue; - - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) { - ++nFoundG5; - efficiencyG5 += info.purity * info.track->lhcbIDs().size() / nTrue; - ; - } - - plots( "selected", info.track ); - plots( "mcSelectedAndRec", mcPart ); - - plot( info.purity, "hitpurity", -0.01, 1.01, 51 ); - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) plot( info.purity, "hitpurityG5", -0.01, 1.01, 51 ); - - m_nClone += info.clone; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) m_nCloneG5 += info.clone; - - if ( reconstructibleB ) { - - if ( info.track != 0 && ( info.track->type() == LHCb::Track::Types::Long || !m_requireLong.value() ) ) { - ++nFoundB; - efficiencyB += info.purity * info.track->lhcbIDs().size() / nTrue; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) { - ++nFoundG5B; - efficiencyG5B += info.purity * info.track->lhcbIDs().size() / nTrue; - } - - plot( info.purity, "hitpurityB", -0.01, 1.01, 51 ); - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) plot( info.purity, "hitpurityG5B", -0.01, 1.01, 51 ); - - m_nCloneB += info.clone; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) m_nCloneG5B += info.clone; - } - } - - if ( reconstructibleBAll ) { - - if ( info.track != 0 && ( info.track->type() == LHCb::Track::Types::Long || !m_requireLong.value() ) ) { - ++nFoundBAll; - efficiencyBAll += info.purity * info.track->lhcbIDs().size() / nTrue; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) { - ++nFoundG5BAll; - efficiencyG5BAll += info.purity * info.track->lhcbIDs().size() / nTrue; - } - - plot( info.purity, "hitpurityBAll", -0.01, 1.01, 51 ); - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) plot( info.purity, "hitpurityG5BAll", -0.01, 1.01, 51 ); - - m_nCloneBAll += info.clone; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) - m_nCloneG5BAll += info.clone * info.track->lhcbIDs().size() / nTrue; - } - } - - if ( reconstructibleKsL ) { - - if ( info.track != 0 && ( info.track->type() == LHCb::Track::Types::Long || !m_requireLong.value() ) ) { - ++nFoundKsL; - efficiencyKsL += info.purity * info.track->lhcbIDs().size() / nTrue; - - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) { - ++nFoundG5KsL; - efficiencyG5KsL += info.purity * info.track->lhcbIDs().size() / nTrue; - } - - plot( info.purity, "hitpurityKsL", -0.01, 1.01, 51 ); - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) plot( info.purity, "hitpurityG5KsL", -0.01, 1.01, 51 ); - - m_nCloneKsL += info.clone; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) m_nCloneG5KsL += info.clone; - } - } - - if ( largeIP ) { - - if ( info.track != 0 && ( info.track->type() == LHCb::Track::Types::Long || !m_requireLong.value() ) ) { - ++nFoundIP; - efficiencyIP += info.purity * info.track->lhcbIDs().size() / nTrue; - - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) { - ++nFoundG5IP; - efficiencyG5IP += info.purity * info.track->lhcbIDs().size() / nTrue; - } - - plot( info.purity, "hitpurityIP", -0.01, 1.01, 51 ); - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) plot( info.purity, "hitpurityG5IP", -0.01, 1.01, 51 ); - - m_nCloneIP += info.clone; - if ( mcPart->p() > 5 * Gaudi::Units::GeV ) m_nCloneG5IP += info.clone; - } - } - } - } - } // loop particles - - // update counters - m_nToFind += nToFind; - m_nFound += nFound; - m_nToFindB += nToFindB; - m_nFoundB += nFoundB; - m_nToFindBAll += nToFindBAll; - m_nFoundBAll += nFoundBAll; - m_nToFindKsL += nToFindKsL; - m_nFoundKsL += nFoundKsL; - m_nToFindIP += nToFindIP; - m_nFoundIP += nFoundIP; - - m_nToFindG5 += nToFindG5; - m_nFoundG5 += nFoundG5; - m_nToFindG5B += nToFindG5B; - m_nFoundG5B += nFoundG5B; - m_nToFindG5BAll += nToFindG5BAll; - m_nFoundG5BAll += nFoundG5BAll; - m_nToFindG5KsL += nToFindG5KsL; - m_nFoundG5KsL += nFoundG5KsL; - m_nToFindG5IP += nToFindG5IP; - m_nFoundG5IP += nFoundG5IP; - - m_efficiency += efficiency; - m_efficiencyG5 += efficiencyG5; - m_efficiencyB += efficiencyB; - m_efficiencyG5B += efficiencyG5B; - m_efficiencyBAll += efficiencyBAll; - m_efficiencyG5BAll += efficiencyG5BAll; - m_efficiencyKsL += efficiencyKsL; - m_efficiencyG5KsL += efficiencyG5KsL; - m_efficiencyIP += efficiencyIP; - m_efficiencyG5IP += efficiencyG5IP; - - // event efficiency - if ( nToFind != 0u ) plot( nFound / double( nToFind ), "eff", -0.01, 1.01, 51 ); - if ( nToFindG5 != 0u ) plot( nFoundG5 / double( nToFindG5 ), "effG5", -0.01, 1.01, 51 ); - if ( nToFindB != 0u ) plot( nFoundB / double( nToFindB ), "effB", -0.01, 1.01, 51 ); - if ( nToFindG5B != 0u ) plot( nFoundG5B / double( nToFindG5B ), "effG5B", -0.01, 1.01, 51 ); - if ( nToFindBAll != 0u ) plot( nFoundBAll / double( nToFindBAll ), "effBAll", -0.01, 1.01, 51 ); - if ( nToFindG5BAll != 0u ) plot( nFoundG5BAll / double( nToFindG5BAll ), "effG5BAll", -0.01, 1.01, 51 ); - if ( nToFindKsL != 0u ) plot( nFoundKsL / double( nToFindKsL ), "effKsL", -0.01, 1.01, 51 ); - if ( nToFindG5KsL != 0u ) plot( nFoundG5KsL / double( nToFindG5KsL ), "effG5KsL", -0.01, 1.01, 51 ); - - if ( fullDetail() == true ) { - // ghost rate versus # interactions - long nVert = visPrimVertTool()->countVertices(); - if ( nToFind != 0u ) - plot2D( nVert, 100 * nFound / double( nToFind ), "eff vs visible", -0.5, 20.5, -0.5, 100.5, 21, 101 ); - } -} - -void TrackEffChecker::plots( const std::string& type, const LHCb::MCParticle* part ) const { - - plot( part->p() / Gaudi::Units::GeV, type + "/p", "p", 0., 50., 25 ); - plot( part->pt() / Gaudi::Units::GeV, type + "/pt", "pt", 0., 10., 20 ); - plot( part->pseudoRapidity(), type + "/eta", "eta", 1., 6., 25 ); -} - -void TrackEffChecker::plots( const std::string& type, const LHCb::Track* track ) const { - - // const double nMeas = weightedMeasurementSum(track); - - if ( track->type() != LHCb::Track::Types::Velo && track->history() != LHCb::Track::History::PrPixel ) { - plot( track->pt() / Gaudi::Units::GeV, type + "/pt", "pt", 0., 10., 100 ); - plot( track->p() / Gaudi::Units::GeV, type + "/p", "p", 0., 50., 25 ); - } - plot( track->chi2PerDoF(), type + "/chi2", "chi2", 0., 500., 100 ); - plot( track->probChi2(), type + "/probchi2", "probChi2", 0., 1., 100 ); - plot( track->pseudoRapidity(), type + "/eta", "eta", 1., 6., 25 ); - - // information from the extra info list - const LHCb::Track::ExtraInfo& info = track->extraInfo(); - LHCb::Track::ExtraInfo::const_iterator iterInfo = info.begin(); - for ( ; iterInfo != info.end(); ++iterInfo ) { - LHCb::Track::AdditionalInfo infoName = LHCb::Track::AdditionalInfo( iterInfo->first ); - std::string title = Gaudi::Utils::toString( infoName ); - const TrackMaps::InfoHistMap& histMap = TrackMaps::infoHistDescription(); - TrackMaps::InfoHistMap::const_iterator iterM = histMap.find( infoName ); - if ( iterM != histMap.end() ) { - const TrackMaps::HistoRange range = histMap.find( infoName )->second; - plot( iterInfo->second, type + "/info/" + range.fid, title, range.fxMin, range.fxMax, 100 ); - } - } // iterInfo -} - -double TrackEffChecker::weightedMeasurementSum( const LHCb::Track* aTrack ) const { - - double wSum = 0; - const std::vector& ids = aTrack->lhcbIDs(); - for ( std::vector::const_iterator iter = ids.begin(); iter != ids.end(); ++iter ) { - if ( iter->isVP() == true ) { - wSum += 1; - } else if ( iter->isFT() == true ) { - wSum += 1; - } else { - // nothing - } - } - return wSum; -} - -StatusCode TrackEffChecker::finalize() { - - // ghost rate - std::string histName = "ghost rate"; - AIDA::IHistogram1D* hist = histo1D( histName ); - double eGhost = 0; - if ( hist != 0 ) eGhost = hist->mean(); - const double tGhost = - m_nTrackCounter.sum() == 0 ? 0.0 : double( m_nGhostCounter.sum() ) / double( m_nTrackCounter.sum() ); - - histName = "hitpurity"; - hist = histo1D( histName ); - double pur = 0; - if ( hist != 0 ) pur = hist->mean(); - - histName = "hitpurityG5"; - hist = histo1D( histName ); - double purG5 = 0; - if ( hist != 0 ) purG5 = hist->mean(); - - histName = "hitpurityB"; - hist = histo1D( histName ); - double purB = 0; - if ( hist != 0 ) purB = hist->mean(); - - histName = "hitpurityG5B"; - hist = histo1D( histName ); - double purG5B = 0; - if ( hist != 0 ) purG5B = hist->mean(); - - histName = "hitpurityBAll"; - hist = histo1D( histName ); - double purBAll = 0; - if ( hist != 0 ) purBAll = hist->mean(); - - histName = "hitpurityG5BAll"; - hist = histo1D( histName ); - double purG5BAll = 0; - if ( hist != 0 ) purG5BAll = hist->mean(); - - histName = "hitpurityKsL"; - hist = histo1D( histName ); - double purKsL = 0; - if ( hist != 0 ) purKsL = hist->mean(); - - histName = "hitpurityG5KsL"; - hist = histo1D( histName ); - double purG5KsL = 0; - if ( hist != 0 ) purG5KsL = hist->mean(); - - histName = "hitpurityIP"; - hist = histo1D( histName ); - double purIP = 0; - if ( hist != 0 ) purIP = hist->mean(); - - histName = "hitpurityG5IP"; - hist = histo1D( histName ); - double purG5IP = 0; - if ( hist != 0 ) purG5IP = hist->mean(); - - const double tEff = m_nToFind.nEntries() == 0 ? 0.0 : double( m_nFound.nEntries() ) / double( m_nToFind.nEntries() ); - - const double tEffG5 = - m_nToFindG5.nEntries() == 0 ? 0.0 : double( m_nFoundG5.nEntries() ) / double( m_nToFindG5.nEntries() ); - - const double hEff = - m_nFound.nEntries() == 0 ? 0.0 : double( m_efficiency.nEntries() ) / double( m_nFound.nEntries() ); - - const double hEffG5 = - m_nFoundG5.nEntries() == 0 ? 0.0 : double( m_efficiencyG5.nEntries() ) / double( m_nFoundG5.nEntries() ); - - const double tCloneG5 = - m_nToFindG5.nEntries() == 0 ? 0.0 : double( m_nCloneG5.nEntries() ) / double( m_nToFindG5.nEntries() ); - - const double tClone = - m_nToFind.nEntries() == 0 ? 0.0 : double( m_nClone.nEntries() ) / double( m_nToFind.nEntries() ); - - const double tEffB = - m_nToFindB.nEntries() == 0 ? 0.0 : double( m_nFoundB.nEntries() ) / double( m_nToFindB.nEntries() ); - - const double tEffG5B = - m_nToFindG5B.nEntries() == 0 ? 0.0 : double( m_nFoundG5B.nEntries() ) / double( m_nToFindG5B.nEntries() ); - - const double hEffB = - m_nFoundB.nEntries() == 0 ? 0.0 : double( m_efficiencyB.nEntries() ) / double( m_nFoundB.nEntries() ); - - const double hEffG5B = - m_nFoundG5B.nEntries() == 0 ? 0.0 : double( m_efficiencyG5B.nEntries() ) / double( m_nFoundG5B.nEntries() ); - - const double tCloneG5B = - m_nToFindG5B.nEntries() == 0 ? 0.0 : double( m_nCloneG5B.nEntries() ) / double( m_nToFindG5B.nEntries() ); - - const double tCloneB = - m_nToFindB.nEntries() == 0 ? 0.0 : double( m_nCloneB.nEntries() ) / double( m_nToFindB.nEntries() ); - - const double tEffBAll = - m_nToFindBAll.nEntries() == 0 ? 0.0 : double( m_nFoundBAll.nEntries() ) / double( m_nToFindBAll.nEntries() ); - - const double tEffG5BAll = m_nToFindG5BAll.nEntries() == 0 - ? 0.0 - : double( m_nFoundG5BAll.nEntries() ) / double( m_nToFindG5BAll.nEntries() ); - - const double hEffBAll = - m_nFoundBAll.nEntries() == 0 ? 0.0 : double( m_efficiencyBAll.nEntries() ) / double( m_nFoundBAll.nEntries() ); - - const double hEffG5BAll = m_nFoundG5BAll.nEntries() == 0 - ? 0.0 - : double( m_efficiencyG5BAll.nEntries() ) / double( m_nFoundG5BAll.nEntries() ); - - const double tCloneG5BAll = m_nToFindG5BAll.nEntries() == 0 - ? 0.0 - : double( m_nCloneG5BAll.nEntries() ) / double( m_nToFindG5BAll.nEntries() ); - - const double tCloneBAll = - m_nToFindBAll.nEntries() == 0 ? 0.0 : double( m_nCloneBAll.nEntries() ) / double( m_nToFindBAll.nEntries() ); - - const double tEffKsL = - m_nToFindKsL.nEntries() == 0 ? 0.0 : double( m_nFoundKsL.nEntries() ) / double( m_nToFindKsL.nEntries() ); - - const double tEffG5KsL = - m_nToFindG5KsL.nEntries() == 0 ? 0.0 : double( m_nFoundG5KsL.nEntries() ) / double( m_nToFindG5KsL.nEntries() ); - - const double hEffKsL = - m_nFoundKsL.nEntries() == 0 ? 0.0 : double( m_efficiencyKsL.nEntries() ) / double( m_nFoundKsL.nEntries() ); - - const double hEffG5KsL = - m_nFoundG5KsL.nEntries() == 0 ? 0.0 : double( m_efficiencyG5KsL.nEntries() ) / double( m_nFoundG5KsL.nEntries() ); - - const double tCloneG5KsL = - m_nToFindG5KsL.nEntries() == 0 ? 0.0 : double( m_nCloneG5KsL.nEntries() ) / double( m_nToFindG5KsL.nEntries() ); - - const double tCloneKsL = - m_nToFindKsL.nEntries() == 0 ? 0.0 : double( m_nCloneKsL.nEntries() ) / double( m_nToFindKsL.nEntries() ); - - const double tEffIP = - m_nToFindIP.nEntries() == 0 ? 0.0 : double( m_nFoundIP.nEntries() ) / double( m_nToFindIP.nEntries() ); - - const double tEffG5IP = - m_nToFindG5IP.nEntries() == 0 ? 0.0 : double( m_nFoundG5IP.nEntries() ) / double( m_nToFindG5IP.nEntries() ); - - const double hEffIP = - m_nFoundIP.nEntries() == 0 ? 0.0 : double( m_efficiencyIP.nEntries() ) / double( m_nFoundIP.nEntries() ); - - const double hEffG5IP = - m_nFoundG5IP.nEntries() == 0 ? 0.0 : double( m_efficiencyG5IP.nEntries() ) / double( m_nFoundG5IP.nEntries() ); - - const double tCloneG5IP = - m_nToFindG5IP.nEntries() == 0 ? 0.0 : double( m_nCloneG5IP.nEntries() ) / double( m_nToFindG5IP.nEntries() ); - - const double tCloneIP = - m_nToFindIP.nEntries() == 0 ? 0.0 : double( m_nCloneIP.nEntries() ) / double( m_nToFindIP.nEntries() ); - - info() << " ************** " - << format( " %8.0f tracks including %8.0f ghosts [%5.1f %%] Event average %5.1f ****", - double( m_nTrackCounter.sum() ), double( m_nGhostCounter.sum() ), tGhost * 100, eGhost * 100 ) - << endmsg; - - info() << " all long " - << format( " : %8.0f from %8.0f [%5.1f %%]; %5.1f %% clones; %5.1f %% hit purity; %5.1f %% hit efficiency", - double( m_nFound.nEntries() ), double( m_nToFind.nEntries() ), tEff * 100, tClone * 100, pur * 100, - hEff * 100 ) - << endmsg; - - info() << " long, p > 5 GeV " - << format( " : %8.0f from %8.0f [%5.1f %%]; %5.1f %% clones; %5.1f %% hit purity; %5.1f %% hit efficiecy", - double( m_nFoundG5.nEntries() ), double( m_nToFindG5.nEntries() ), tEffG5 * 100, tCloneG5 * 100, - purG5 * 100, hEffG5 * 100 ) - << endmsg; - - if ( m_nToFindB.nEntries() > 0 ) { - info() << " all long B daughers " - << format( - " : %8.0f from %8.0f [%5.1f %%]; %5.1f %% clones; %5.1f %% hit purity; %5.1f %% hit efficiency", - double( m_nFoundB.nEntries() ), double( m_nToFindB.nEntries() ), tEffB * 100, tCloneB * 100, - purB * 100, hEffB * 100 ) - << endmsg; - - info() << " long B daughters, p > 5 GeV " - << format( - " : %8.0f from %8.0f [%5.1f %%]; %5.1f %% clones; %5.1f %% hit purity; %5.1f %% hit efficiency", - double( m_nFoundG5B.nEntries() ), double( m_nToFindG5B.nEntries() ), tEffG5B * 100, tCloneG5B * 100, - purG5B * 100, hEffG5B * 100 ) - << endmsg; - } - - if ( m_nToFindBAll.nEntries() > 0 ) { - info() << " all long good B daughers " - << format( - " : %8.0f from %8.0f [%5.1f %%]; %5.1f %% clones; %5.1f %% hit purity; %5.1f %% hit efficiency", - double( m_nFoundBAll.nEntries() ), double( m_nToFindBAll.nEntries() ), tEffBAll * 100, - tCloneBAll * 100, purBAll * 100, hEffBAll * 100 ) - << endmsg; - - info() << " long good B daughters , p > 5 GeV " - << format( - " : %8.0f from %8.0f [%5.1f %%]; %5.1f %% clones; %5.1f %% hit purity; %5.1f %% hit efficiency", - double( m_nFoundG5BAll.nEntries() ), double( m_nToFindG5BAll.nEntries() ), tEffG5BAll * 100, - tCloneG5BAll * 100, purG5BAll * 100, hEffG5BAll * 100 ) - << endmsg; - } - - if ( m_nToFindKsL.nEntries() > 0 ) { - info() << " all long Ks/Lambda " - << format( - " : %8.0f from %8.0f [%5.1f %%]; %5.1f %% clones; %5.1f %% hit purity; %5.1f %% hit efficiency", - double( m_nFoundKsL.nEntries() ), double( m_nToFindKsL.nEntries() ), tEffKsL * 100, tCloneKsL * 100, - purKsL * 100, hEffKsL * 100 ) - << endmsg; - - info() << " long Ks/Lambda, p > 5 GeV " - << format( - " : %8.0f from %8.0f [%5.1f %%]; %5.1f %% clones; %5.1f %% hit purity; %5.1f %% hit efficiency", - double( m_nFoundG5KsL.nEntries() ), double( m_nToFindG5KsL.nEntries() ), tEffG5KsL * 100, - tCloneG5KsL * 100, purG5KsL * 100, hEffG5KsL * 100 ) - << endmsg; - } - - if ( m_nToFindIP.nEntries() > 0 ) { - info() << " all long large IP " - << format( - " : %8.0f from %8.0f [%5.1f %%]; %5.1f %% clones; %5.1f %% hit purity; %5.1f %% hit efficiency", - double( m_nFoundIP.nEntries() ), double( m_nToFindIP.nEntries() ), tEffIP * 100, tCloneIP * 100, - purIP * 100, hEffIP * 100 ) - << endmsg; - - info() << " long largeIP, p > 5 GeV " - << format( - " : %8.0f from %8.0f [%5.1f %%]; %5.1f %% clones; %5.1f %% hit purity; %5.1f %% hit efficiency", - double( m_nFoundG5IP.nEntries() ), double( m_nToFindG5IP.nEntries() ), tEffG5IP * 100, - tCloneG5IP * 100, purG5IP * 100, hEffG5IP * 100 ) - << endmsg; - } - - return TrackCheckerBase::finalize(); -} diff --git a/Tr/TrackCheckers/src/TrackResChecker.cpp b/Tr/TrackCheckers/src/TrackResChecker.cpp index cf5d9767e22921baf608b49f71b8148b1a4b640c..59c788786c292fbd5fea9e179fc4ed47beba9a0e 100644 --- a/Tr/TrackCheckers/src/TrackResChecker.cpp +++ b/Tr/TrackCheckers/src/TrackResChecker.cpp @@ -17,8 +17,7 @@ #include "Event/State.h" #include "Event/StateVector.h" #include "Event/Track.h" -#include "GaudiAlg/GaudiHistoTool.h" -#include "GaudiAlg/IHistoTool.h" +#include "Gaudi/Accumulators/StaticHistogram.h" #include "GaudiKernel/SystemOfUnits.h" #include "GaudiUtils/HistoStats.h" #include "LHCbAlgs/Consumer.h" @@ -27,9 +26,52 @@ #include "TrackInterfaces/ITrackProjector.h" #include "TrackInterfaces/ITrackProjectorSelector.h" #include "TrackKernel/TrackFunctors.h" + +#include #include #include +class TrackResChecker; + +namespace { + /// little wrapper around a map of histograms indexed by type/location + /// implementing thread safe insertion of new histograms on demand + template + struct HistoMapImpl : std::map> { + std::mutex m{}; // thread safety of map insertion + TrackResChecker* algo; + std::string name, title; + std::tuple axis; + + HistoMapImpl( TrackResChecker* _algo, std::string _name, std::string _title, Axis... _axis ) + : algo( _algo ), name( _name ), title( _title ), axis( _axis... ) {} + + auto& get( LHCb::Track::Types t ) { return get( toString( t ) ); } + auto& get( LHCb::Track::Types t, std::string const& dir ) { return get( toString( t ) + "/" + dir ); } + auto& get( std::string const& t ) { + if ( !this->contains( t ) ) { + std::scoped_lock lock{m}; + if ( !this->contains( t ) ) { + // now we are alone and there is still nothing, let's create the histogram + this->emplace( std::piecewise_construct, std::forward_as_tuple( t ), + std::forward_as_tuple( algo, fmt::format( "{}/{}", t, name ), title, axis ) ); + } + } + return this->at( t ); + } + }; + template + struct HistoMap; + template <> + struct HistoMap<1> : HistoMapImpl<1, Gaudi::Accumulators::Axis> { + using HistoMapImpl::HistoMapImpl; + }; + template <> + struct HistoMap<2> : HistoMapImpl<2, Gaudi::Accumulators::Axis, Gaudi::Accumulators::Axis> { + using HistoMapImpl::HistoMapImpl; + }; +} // namespace + /** @class TrackResChecker TrackResChecker.h * * Class for track monitoring @@ -46,38 +88,70 @@ public: /** Standard constructor */ TrackResChecker( const std::string& name, ISvcLocator* pSvcLocator ); - /** Algorithm initialize */ - StatusCode initialize() override; - /** Algorithm execute */ void operator()( LHCb::Track::Range const&, LHCb::MCParticles const&, LHCb::LinksByKey const&, DetectorElement const& ) const override; - /** Algorithm finalize */ - StatusCode finalize() override; - private: - void resolutionHistos( IHistoTool const& histotool, LHCb::Track const& track, LHCb::MCParticle const& mcPart, + void resolutionHistos( LHCb::Track::Types type, LHCb::Track const& track, LHCb::MCParticle const& mcPart, IGeometryInfo const& ) const; - void pullplots( const IHistoTool& histotool, const LHCb::State& trueState, const LHCb::State& recState, + void pullplots( LHCb::Track::Types type, const LHCb::State& trueState, const LHCb::State& recState, const std::string& location ) const; - void plotsByMeasType( IHistoTool const& histotool, LHCb::Track const& track, LHCb::MCParticle const& mcPart, + void plotsByMeasType( LHCb::Track::Types type, LHCb::Track const& track, LHCb::MCParticle const& mcPart, IGeometryInfo const& ) const; - const IHistoTool* createHistoTool( const std::string& name ) const; - private: Gaudi::Property m_plotsByMeasType{this, "PlotsByMeasType", false}; + Gaudi::Property m_fullDetail{this, "FullDetail", false}; ToolHandle m_projectorSelector{this, "ProjectorSelector", "TrackProjectorSelector/Projector"}; - std::unordered_map m_histoTools; - mutable Gaudi::Accumulators::MsgCounter m_differentFitHistory{this, "Tracks have different fit history"}; + + mutable HistoMap<1> m_chi2PerDof{this, "chi2PerDof", "chi2PerDof", {1000, 0., 100.}}; + mutable HistoMap<1> m_probChi2{this, "probChi2", "probChi2", {50, 0., 1.}}; + mutable HistoMap<1> m_fitStatus{this, "fitStatus", "fit status", {5, -0.5, 4.5}}; + mutable HistoMap<1> m_truemom{this, "truemom", "true p [GeV]", {100, 0, 50}}; + mutable HistoMap<1> m_truept{this, "truept", "true pT [GeV]", {100, 0, 10}}; + mutable HistoMap<1> m_correctcharge{this, "correctcharge", "correct charge", {2, -0.5, 1.5}}; + + mutable HistoMap<1> m_xres{this, "xres", "x resolution / mm", {101, -0.4, 0.4}}; + mutable HistoMap<1> m_yres{this, "yres", "y resolution / mm", {101, -0.4, 0.4}}; + mutable HistoMap<1> m_txres{this, "txres", "tx resolution", {101, -0.0025, 0.0025}}; + mutable HistoMap<1> m_tyres{this, "tyres", "ty resolution", {101, -0.0025, 0.0025}}; + + mutable HistoMap<1> m_xpull{this, "xpull", "x pull", {101, -5., 5.}}; + mutable HistoMap<1> m_ypull{this, "ypull", "y pull", {101, -5., 5.}}; + mutable HistoMap<1> m_txpull{this, "txpull", "tx pull", {101, -5., 5.}}; + mutable HistoMap<1> m_typull{this, "typull", "ty pull", {101, -5., 5.}}; + + mutable HistoMap<1> m_qopres{this, "qop_res", "qop", {101, -0.02, 0.02}}; + mutable HistoMap<1> m_qoppull{this, "qoppull", "qop pull", {101, -5., 5.}}; + mutable HistoMap<1> m_ppull{this, "ppull", "p pull", {101, -5., 5.}}; + mutable HistoMap<1> m_dpoverp{this, "dpoverp", "dp/p", {101, -0.05, 0.05}}; + mutable HistoMap<1> m_expecteddpoverp{this, "expecteddpoverp", "expected dp/p", {101, 0., 0.01}}; + + mutable HistoMap<1> m_measres{this, "/meas_res", " Measurement resolution", {100, -0.5, 0.5}}; + mutable HistoMap<1> m_measpull{this, "/meas_pull", " Measurement pull", {100, -5., 5.}}; + mutable HistoMap<1> m_measchi2{this, "/meas_chi2", " Measurement chi2", {200, 0., 10.}}; + mutable HistoMap<2> m_measchi2ndf{this, "/meas_chi2ndf", "Measurement chi2 vs ndf", {200, 0., 10.}, {2, 0., 3.}}; + + mutable HistoMap<2> m_vdopp{this, "vertex/dpoverp_vs_p", "dp/p vs p", {25, 0., 50.}, {50, -0.1, 0.1}}; + mutable HistoMap<2> m_vdopeta{this, "vertex/dpoverp_vs_eta", "dp/p vs eta", {20, 2., 5.}, {50, -0.05, 0.05}}; + mutable HistoMap<2> m_vpullp{this, + "vertex/p_pull_vs_p", + "p pull vs p", + { + 25, + 0., + 50., + }, + {50, -10., 10.}}; + mutable HistoMap<2> m_vpulleta{this, "vertex/p_pull_vs_eta", "p pull vs eta", {20, 2., 5.}, {50, -10., 10.}}; }; DECLARE_COMPONENT( TrackResChecker ) @@ -92,29 +166,6 @@ TrackResChecker::TrackResChecker( const std::string& name, ISvcLocator* pSvcLoca KeyValue{"LinkerInTable", "Link/" + LHCb::TrackLocation::Default}, KeyValue{"StandardGeometryTop", LHCb::standard_geometry_top}} ) {} -StatusCode TrackResChecker::initialize() { - return Consumer::initialize().andThen( [&] { - m_histoTools[LHCb::Track::Types::Unknown] = createHistoTool( "ALL" ); - if ( splitByType() ) { - using Type = LHCb::Track::Types; - constexpr auto types = std::array{Type::Velo, Type::VeloBackward, Type::Long, Type::Upstream, - Type::Downstream, Type::Ttrack, Type::Muon}; - for ( auto type : types ) { m_histoTools[type] = createHistoTool( toString( type ) ); } - } - } ); -} - -const IHistoTool* TrackResChecker::createHistoTool( const std::string& name ) const { - IHistoTool* htool = tool( "HistoTool", name, this ); - GaudiHistoTool* ghtool = dynamic_cast( htool ); - ghtool->setHistoTopDir( histoPath() + "/" ); - std::string histodir = ghtool->histoDir(); - size_t pos = histodir.find( '.' ); - if ( pos != std::string::npos ) histodir.erase( 0, pos + 1 ); - ghtool->setHistoDir( histodir ); - return htool; -} - //============================================================================= // Execute //============================================================================= @@ -135,30 +186,24 @@ void TrackResChecker::operator()( LHCb::Track::Range const& tracks, LHCb::MCPart // we actually just want to know if it passes the 'selector' inside the IMCReconstructible // && int(selector()->reconstructible(mcparticle)) > int(IMCReconstructible::NotReconstructible) ) { - // split by type.. - const IHistoTool* histotool( 0 ); - if ( splitByType() ) { - histotool = m_histoTools.at( track->type() ); - } else { - histotool = m_histoTools.at( LHCb::Track::Types::Unknown ); - } + auto type = splitByType() ? track->type() : LHCb::Track::Types::Unknown; // resolutions at predefined z. - resolutionHistos( *histotool, *track, *mcparticle, *lhcb.geometry() ); + resolutionHistos( type, *track, *mcparticle, *lhcb.geometry() ); // prob chi^2 - histotool->plot1D( track->chi2PerDoF(), "chi2PerDof", "chi2PerDof", 0., 100., 1000 ); - histotool->plot1D( track->probChi2(), "probChi2", "probChi2", 0., 1., 50 ); + ++m_chi2PerDof.get( type )[track->chi2PerDoF()]; + ++m_probChi2.get( type )[track->probChi2()]; // fit status - histotool->plot1D( static_cast( track->fitStatus() ), "fitStatus", "fit status", -0.5, 4.5, 5 ); + ++m_fitStatus.get( type )[static_cast( track->fitStatus() )]; - histotool->plot1D( mcparticle->p() / Gaudi::Units::GeV, "truemom", "true p [GeV]", 0, 50, 100 ); - histotool->plot1D( mcparticle->pt() / Gaudi::Units::GeV, "truept", "true pT [GeV]", 0, 10, 100 ); + ++m_truemom.get( type )[mcparticle->p() / Gaudi::Units::GeV]; + ++m_truept.get( type )[mcparticle->pt() / Gaudi::Units::GeV]; // Resolutions and pulls per Measurement type if ( m_plotsByMeasType.value() && nMeasurements( *track ) > 0 ) - plotsByMeasType( *histotool, *track, *mcparticle, *lhcb.geometry() ); + plotsByMeasType( type, *track, *mcparticle, *lhcb.geometry() ); } } } @@ -166,14 +211,14 @@ void TrackResChecker::operator()( LHCb::Track::Range const& tracks, LHCb::MCPart //============================================================================= // //============================================================================= -void TrackResChecker::resolutionHistos( IHistoTool const& htool, LHCb::Track const& track, +void TrackResChecker::resolutionHistos( LHCb::Track::Types type, LHCb::Track const& track, LHCb::MCParticle const& mcPart, IGeometryInfo const& geometry ) const { // pulls at vertex LHCb::State trueStateVertex; idealStateCreator()->createStateVertex( &mcPart, trueStateVertex ).ignore(); LHCb::State vtxState; StatusCode sc = extrapolator()->propagate( track, trueStateVertex.z(), vtxState, geometry ); - if ( sc.isSuccess() ) pullplots( htool, trueStateVertex, vtxState, "vertex" ); + if ( sc.isSuccess() ) pullplots( type, trueStateVertex, vtxState, "vertex" ); // for vertex also make some 2-d plots if ( track.type() == LHCb::Track::Types::Long || track.type() == LHCb::Track::Types::Upstream || @@ -182,23 +227,22 @@ void TrackResChecker::resolutionHistos( IHistoTool const& htool, LHCb::Track con const double ptrue = mcPart.p(); const double eta = mcPart.pseudoRapidity(); // track.pseudoRapidity(); - htool.plot2D( ptrue / Gaudi::Units::GeV, invp * ptrue - 1, "vertex/dpoverp_vs_p", "dp/p vs p", 0., 50., -0.1, 0.1, - 25, 50 ); - htool.plot2D( eta, invp * ptrue - 1, "vertex/dpoverp_vs_eta", "dp/p vs eta", 2., 5., -0.05, 0.05, 20, 50 ); + ++m_vdopp.get( type )[{ptrue / Gaudi::Units::GeV, invp * ptrue - 1}]; + ++m_vdopeta.get( type )[{eta, invp * ptrue - 1}]; const double invperr2 = track.firstState().covariance()( 4, 4 ); if ( invperr2 > 0 ) { const double ppull = ( invp - 1 / ptrue ) / std::sqrt( invperr2 ); - htool.plot2D( ptrue / Gaudi::Units::GeV, ppull, "vertex/p_pull_vs_p", "p pull vs p", 0., 50., -10., 10., 25, 50 ); - htool.plot2D( eta, ppull, "vertex/p_pull_vs_eta", "p pull vs eta", 2., 5., -10., 10., 20, 50 ); + ++m_vpullp.get( type )[{ptrue / Gaudi::Units::GeV, ppull}]; + ++m_vpulleta.get( type )[{eta, ppull}]; } } // fraction of tracks with correct charge bool correctcharge = track.firstState().qOverP() * mcPart.particleID().threeCharge() > 0; - htool.plot1D( correctcharge, "correctcharge", "correct charge", -0.5, 1.5, 2 ); + ++m_correctcharge.get( type )[correctcharge]; - if ( fullDetail() ) { + if ( m_fullDetail.value() ) { for ( const LHCb::State* state : track.states() ) { // skip the closest to beam, since we already have it if ( state->location() == LHCb::State::Location::ClosestToBeam ) continue; @@ -210,7 +254,7 @@ void TrackResChecker::resolutionHistos( IHistoTool const& htool, LHCb::Track con std::string location = state->location() != LHCb::State::Location::LocationUnknown ? toString( state->location() ) : format( "state_%d_mm", int( state_z ) ); - pullplots( htool, trueState, *state, location ); + pullplots( type, trueState, *state, location ); } } } @@ -220,7 +264,7 @@ void TrackResChecker::resolutionHistos( IHistoTool const& htool, LHCb::Track con //============================================================================= // //============================================================================= -void TrackResChecker::pullplots( const IHistoTool& htool, const LHCb::State& trueState, const LHCb::State& recState, +void TrackResChecker::pullplots( LHCb::Track::Types type, const LHCb::State& trueState, const LHCb::State& recState, const std::string& location ) const { // save some typing @@ -234,15 +278,15 @@ void TrackResChecker::pullplots( const IHistoTool& htool, const LHCb::State& tru const double dty = vec( 3 ) - trueVec( 3 ); // fill the histograms - htool.plot1D( dx, location + "/x_res", "x resolution / mm", -0.4, 0.4, 101 ); - htool.plot1D( dy, location + "/y_res", "y resolution / mm", -0.4, 0.4, 101 ); - htool.plot1D( dtx, location + "/tx_res", "tx resolution", -0.0025, 0.0025, 101 ); - htool.plot1D( dty, location + "/ty_res", "ty resolution", -0.0025, 0.0025, 101 ); + ++m_xres.get( type, location )[dx]; + ++m_yres.get( type, location )[dy]; + ++m_txres.get( type, location )[dtx]; + ++m_tyres.get( type, location )[dty]; - htool.plot1D( dx / sqrt( cov( 0, 0 ) + trueCov( 0, 0 ) ), location + "/xpull", "x pull", -5., 5., 101 ); - htool.plot1D( dy / sqrt( cov( 1, 1 ) + trueCov( 1, 1 ) ), location + "/ypull", "y pull", -5., 5., 101 ); - htool.plot1D( dtx / sqrt( cov( 2, 2 ) + trueCov( 2, 2 ) ), location + "/txpull", "tx pull", -5., 5., 101 ); - htool.plot1D( dty / sqrt( cov( 3, 3 ) + trueCov( 3, 3 ) ), location + "/typull", "ty pull", -5., 5., 101 ); + ++m_xpull.get( type, location )[dx / sqrt( cov( 0, 0 ) + trueCov( 0, 0 ) )]; + ++m_ypull.get( type, location )[dy / sqrt( cov( 1, 1 ) + trueCov( 1, 1 ) )]; + ++m_txpull.get( type, location )[dtx / sqrt( cov( 2, 2 ) + trueCov( 2, 2 ) )]; + ++m_typull.get( type, location )[dty / sqrt( cov( 3, 3 ) + trueCov( 3, 3 ) )]; if ( std::abs( cov( 4, 4 ) ) > 1e-20 ) { // test that there was a momentum measurement const double qop = vec( 4 ); @@ -251,19 +295,18 @@ void TrackResChecker::pullplots( const IHistoTool& htool, const LHCb::State& tru const double invptrue = std::abs( qoptrue ); const double qoperr = std::sqrt( cov( 4, 4 ) + trueCov( 4, 4 ) ); // make two pulls, to be sensitive to both a curvature and a momentum bias - htool.plot1D( ( qop - qoptrue ) / invptrue, location + "/qop_res", "qop ", -.02, .02, 101 ); - htool.plot1D( ( qop - qoptrue ) / qoperr, location + "/qoppull", "qop pull", -5., 5., 101 ); - htool.plot1D( ( invp - invptrue ) / qoperr, location + "/ppull", "p pull", -5., 5., 101 ); - htool.plot1D( invp / invptrue - 1, location + "/dpoverp", "dp/p", -0.05, 0.05, 101 ); - if ( invp > 0 ) - htool.plot1D( std::sqrt( cov( 4, 4 ) ) / invp, location + "/expecteddpoverp", "expected dp/p", 0., 0.01, 100 ); + ++m_qopres.get( type, location )[( qop - qoptrue ) / invptrue]; + ++m_qoppull.get( type, location )[( qop - qoptrue ) / qoperr]; + ++m_ppull.get( type, location )[( invp - invptrue ) / qoperr]; + ++m_dpoverp.get( type, location )[invp / invptrue - 1]; + if ( invp > 0 ) ++m_expecteddpoverp.get( type, location )[std::sqrt( cov( 4, 4 ) ) / invp]; } } //============================================================================= // //============================================================================= -void TrackResChecker::plotsByMeasType( IHistoTool const& htool, LHCb::Track const& track, +void TrackResChecker::plotsByMeasType( LHCb::Track::Types type, LHCb::Track const& track, LHCb::MCParticle const& mcPart, IGeometryInfo const& geometry ) const { for ( const auto& measure : measurements( track ) ) { @@ -276,7 +319,7 @@ void TrackResChecker::plotsByMeasType( IHistoTool const& htool, LHCb::Track cons StatusCode sc = extrapolator()->propagate( track, measure.z(), stateAtMeas, geometry ); if ( sc.isSuccess() ) { // make pull plots as before - pullplots( htool, trueStateAtMeas, stateAtMeas, dir ); + pullplots( type, trueStateAtMeas, stateAtMeas, dir ); } // Monitor unbiased measurement resolutions @@ -311,52 +354,14 @@ void TrackResChecker::plotsByMeasType( IHistoTool const& htool, LHCb::Track cons } ), projectResult ); - htool.plot1D( res, dir + "/meas_res", " Measurement resolution", -0.5, 0.5, 100 ); - htool.plot1D( res / errMeasure, dir + "/meas_pull", " Measurement pull", -5., 5., 100 ); - htool.plot1D( chi2, dir + "/meas_chi2", " Measurement chi2", 0., 10., 200 ); - htool.plot2D( chi2, ndf, dir + "/meas_chi2ndf", "Measurement chi2 vs ndf", 0., 10., 0., 3., 200, 2 ); + ++m_measres.get( type, dir )[res]; + ++m_measpull.get( type, dir )[res / errMeasure]; + ++m_measchi2.get( type, dir )[chi2]; + ++m_measchi2ndf.get( type, dir )[{chi2, ndf}]; } } else { - Warning( "could not get projector for measurement", StatusCode::SUCCESS, 0 ).ignore(); + warning() << "could not get projector for measurement" << endmsg; } } } // iterate measurements } - -//============================================================================= -// -//============================================================================= -StatusCode TrackResChecker::finalize() { - info() << " ************************************ " << endmsg; - for ( const auto& ihtool : m_histoTools ) { - const IHistoTool* htool = ihtool.second; - const GaudiHistoTool* ghtool = dynamic_cast( htool ); - for ( const auto& name : - {"vertex/xpull", "vertex/ypull", "vertex/txpull", "vertex/typull", "vertex/ppull", "probChi2"} ) { - const auto pull = htool->histo( HistoID( name ) ); - if ( pull ) - info() << ghtool->histoDir() << "/" << std::setiosflags( std::ios_base::left ) << std::setw( 10 ) - << pull->title() << " " - << format( ": mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", pull->mean(), - Gaudi::Utils::HistoStats::meanErr( pull ), pull->rms(), - Gaudi::Utils::HistoStats::rmsErr( pull ) ) - << endmsg; - } - for ( const auto& name : {"vertex/x_res", "vertex/y_res"} ) { - const auto res = htool->histo( HistoID( name ) ); - if ( res ) - info() << ghtool->histoDir() << "/" << res->title() - << format( ": RMS = %5.3f +/- %5.3f micron", res->rms() * 1000, - Gaudi::Utils::HistoStats::rmsErr( res ) * 1000 ) - << endmsg; - } - const auto dpop = htool->histo( HistoID( "vertex/dpoverp" ) ); - if ( dpop ) - info() << ghtool->histoDir() << "/" << dpop->title() - << format( ": mean = %6.4f +/- %6.4f, RMS = %6.4f +/- %6.4f", dpop->mean(), - Gaudi::Utils::HistoStats::meanErr( dpop ), dpop->rms(), - Gaudi::Utils::HistoStats::rmsErr( dpop ) ) - << endmsg; - } - return TrackCheckerBase::finalize(); -} diff --git a/Tr/TrackCheckers/src/TriggerObjectsCompatibilityChecker.cpp b/Tr/TrackCheckers/src/TriggerObjectsCompatibilityChecker.cpp deleted file mode 100644 index 77a2f009db2d66ba59412da5f0378a3ef527cd7c..0000000000000000000000000000000000000000 --- a/Tr/TrackCheckers/src/TriggerObjectsCompatibilityChecker.cpp +++ /dev/null @@ -1,501 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2019 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/MCParticle.h" -#include "Event/MCVertex.h" -#include "Event/RecVertex.h" -#include "Event/Track.h" -#include "Event/TrackVertexUtils.h" -#include "Gaudi/Accumulators/Histogram.h" -#include "GaudiKernel/ToolHandle.h" -#include "LHCbAlgs/Consumer.h" -#include "TrackCheckerBase.h" -#include "TrackKernel/CubicStateInterpolationTraj.h" -#include "TrackKernel/TrackFunctors.h" -#include - -class TriggerObjectsCompatibilityChecker final - : public LHCb::Algorithm::Consumer> { -public: - TriggerObjectsCompatibilityChecker( const std::string& name, ISvcLocator* pSvcLocator ) - : Consumer( name, pSvcLocator, - {KeyValue{"TrackContainerHLT1", LHCb::TrackLocation::Default}, - KeyValue{"TrackContainerHLT2", LHCb::TrackLocation::Default}, - KeyValue{"MCParticleInput", LHCb::MCParticleLocation::Default}, - KeyValue{"LinkerLocationHLT1", "Link/Pr/LHCbID"}, KeyValue{"LinkerLocationHLT2", "Link/Pr/LHCbID"}, - KeyValue{"PVContainerHLT1", LHCb::RecVertexLocation::Primary}, - KeyValue{"PVContainerHLT2", LHCb::RecVertexLocation::Primary}} ) {} - void operator()( const LHCb::Tracks& hlt1_tracks, const LHCb::Tracks& hlt2_track, - const LHCb::MCParticles& mcparticles, const LHCb::LinksByKey& linker_hlt1, - const LHCb::LinksByKey& linker_hlt2, const LHCb::RecVertices& hlt1_pvs, - const LHCb::RecVertices& hlt2_pvs ) const override; - - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_fitted_track{ - this, "fitted_track", "fitted track", {2, 0, 2}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_hlt1_fitted_track{ - this, "hlt1_fitted_track", "hlt1 fitted track", {2, 0, 2}}; - - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_numPVs_hlt1{ - this, "number_reconstructed_PVs_hlt1", "number_reconstructed_PVs_hlt1", {50, 0, 50}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_numPVs_hlt2{ - this, "number_reconstructed_PVs_hlt2", "number_reconstructed_PVs_hlt2", {50, 0, 50}}; - mutable Gaudi::Accumulators::Histogram<1> m_probChi2_hlt1{this, "probChi2_hlt1", "probChi2_hlt1", {1000, 0, 1.}}; - mutable Gaudi::Accumulators::Histogram<1> m_probChi2_hlt2{this, "probChi2_hlt2", "probChi2_hlt2", {1000, 0, 1.}}; - - mutable Gaudi::Accumulators::Histogram<1> m_chi2_hlt1{this, "chi2_hlt1", "chi2_hlt1", {2000, 0, 200}}; - mutable Gaudi::Accumulators::Histogram<1> m_chi2_hlt2{this, "chi2_hlt2", "chi2_hlt2", {2000, 0, 200}}; - - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_ndof_hlt1{ - this, "ndof_hlt1", "ndof_hlt1", {2, 0, 2}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_ndof_hlt2{ - this, "ndof_hlt2", "ndof_hlt2", {2, 0, 2}}; - - mutable Gaudi::Accumulators::Histogram<1> m_chi2_ndof_hlt1{this, "chi2_ndof_hlt1", "chi2_ndof_hlt1", {500, 0, 50}}; - mutable Gaudi::Accumulators::Histogram<1> m_chi2_ndof_hlt2{this, "chi2_ndof_hlt2", "chi2_ndof_hlt2", {500, 0, 50}}; - - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_type_hlt1{ - this, "track_type_hlt1", "track_type_hlt1", {10, 0, 10}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_type_hlt2{ - this, "track_type_hlt2", "track_type_hlt2", {10, 0, 10}}; - - mutable Gaudi::Accumulators::Histogram<1> m_deltaX{ - this, "deltaX_first_state", "hlt2 - hlt1 x first state", {1000, -0.2, 0.2}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaY{ - this, "deltaY_first_state", "hlt2 - hlt1 y firstate", {1000, -0.2, 0.2}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaZ{ - this, "deltaZ_first_state", "hlt2 - hlt1 z firstate", {1000, -10., 10.}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaTX{ - this, "deltaTX_first_state", "hlt2 - hlt1 tx firstate", {1000, -0.01, 0.01}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaTY{ - this, "deltaTY_first_state", "hlt2 - hlt1 ty firstate", {1000, -0.01, 0.01}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaQOP{ - this, "deltaQOP_first_state", "hlt2 - hlt1 qop firstate", {1000, -0.0001, 0.0001}}; - - mutable Gaudi::Accumulators::Histogram<2> m_2d_X{ - this, "x_first_state_2d", "firstate_x;hlt1;hlt2", {{1000, -1., 1.}, {1000, -1., 1.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_2d_Y{ - this, "y_first_state_2d", "firstate_y;hlt1;hlt2", {{1000, -1., 1.}, {1000, -1., 1.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_2d_Z{ - this, "z_first_state_2d", "firstate_z;hlt1;hlt2", {{1000, -200., 600.}, {1000, -200., 600.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_2d_TX{ - this, "tx_first_state_2d", "firstate_tx;hlt1;hlt2", {{1000, -0.3, 0.3}, {1000, -0.3, 0.3}}}; - mutable Gaudi::Accumulators::Histogram<2> m_2d_TY{ - this, "ty_first_state_2d", "firstate_ty;hlt1;hlt2", {{1000, -0.2, 0.2}, {1000, -0.2, 0.2}}}; - mutable Gaudi::Accumulators::Histogram<2> m_2d_QOP{ - this, "qop_first_state_2d", "firstate_qop;hlt1;hlt2", {{1000, -0.001, 0.001}, {1000, -0.001, 0.001}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_2d_PT{ - this, "pt_first_state_2d", "firstate_pt;hlt1;hlt2", {{1000, 0., 10000.}, {1000, 0., 10000.}}}; - - mutable Gaudi::Accumulators::Histogram<1> m_deltaIP3D{this, "recdeltaIP3D", "delta_recIP3D", {100, -0.5, 0.5}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaIPx{this, "recdeltaIPx", "delta_recIPx", {100, -0.5, 0.5}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaIPy{this, "recdeltaIPy", "delta_recIPy", {100, -0.5, 0.5}}; - - mutable Gaudi::Accumulators::Histogram<2> m_IP3D{ - this, "recIP3D", "recIP3D;hlt1;hlt2", {{100, 0., 1.}, {100, 0., 1.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_IPx{this, "recIPx", "recIPx;hlt1;hlt2", {{100, -1., 1.}, {100, -1., 1.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_IPy{this, "recIPy", "recIPy;hlt1;hlt2", {{100, -1., 1.}, {100, -1., 1.}}}; - - mutable Gaudi::Accumulators::Histogram<1> m_deltaIPxErr{this, "delta_recIPxErr", "delta_recIPxErr", {100, -0.2, 0.2}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaIPyErr{this, "delta_recIPyErr", "delta_recIPyErr", {100, -0.2, 0.2}}; - - mutable Gaudi::Accumulators::Histogram<2> m_IPxErr{ - this, "recIPxErr", "recIPxErr;hlt1;hlt2", {{100, 0., 0.5}, {100, 0., 0.5}}}; - mutable Gaudi::Accumulators::Histogram<2> m_IPyErr{ - this, "recIPyErr", "recIPyErr;hlt1;hlt2", {{100, 0., 0.5}, {100, 0., 0.5}}}; - - mutable Gaudi::Accumulators::Histogram<1> m_deltaIPChi2{this, "delta_recIPChi2", "delta_recIPChi2", {500, -10., 10.}}; - mutable Gaudi::Accumulators::Histogram<2> m_IPChi2{ - this, "recIPChi2", "recIPChi2;hlt1;hlt2", {{500, 0., 20.}, {500, 0., 20.}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_IPChi2_zoom{ - this, "recIPChi2_zoom", "recIPChi2_zoom;hlt1;hlt2", {{500, 0., 2.}, {500, 0., 2.}}}; - - mutable Gaudi::Accumulators::Histogram<2, Gaudi::Accumulators::atomicity::full, int> m_nTracksPerPV{ - this, "nTracksPerPV", "nTracksPV;hlt1;hlt2", {{50, 0, 50}, {50, 0, 50}}}; - mutable Gaudi::Accumulators::Histogram<2, Gaudi::Accumulators::atomicity::full, int> m_nTracks{ - this, "nTracks", "nTracks;hlt1;hlt2", {{100, 0, 100}, {100, 0, 100}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_PV_X{ - this, "x_PV", "PVx;hlt1;hlt2", {{1000, -0.2, 0.2}, {1000, -0.2, 0.2}}}; - mutable Gaudi::Accumulators::Histogram<2> m_PV_Y{ - this, "y_PV", "PVy;hlt1;hlt2", {{1000, -0.2, 0.2}, {1000, -0.2, 0.2}}}; - mutable Gaudi::Accumulators::Histogram<2> m_PV_Z{ - this, "z_PV", "PVz;hlt1;hlt2", {{1000, -200., 200.}, {1000, -200., 200.}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_PV_cov00{ - this, "cov00_PV", "PVcov00;hlt1;hlt2", {{1000, -1., 1.}, {1000, -1., 1.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_PV_cov01{ - this, "cov01_PV", "PVcov01;hlt1;hlt2", {{1000, -1., 1.}, {1000, -1., 1.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_PV_cov02{ - this, "cov02_PV", "PVcov02;hlt1;hlt2", {{1000, -1., 1.}, {1000, -1., 1.}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_PV_cov11{ - this, "cov11_PV", "PVcov11;hlt1;hlt2", {{1000, -1., 1.}, {1000, -1., 1.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_PV_cov12{ - this, "cov12_PV", "PVcov12;hlt1;hlt2", {{1000, -1., 1.}, {1000, -1., 1.}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_PV_cov22{ - this, "cov22_PV", "PVcov12;hlt1;hlt2", {{1000, -1., 1.}, {1000, -1., 1.}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_X{ - this, "x_stAtVtx", "stAtVtx_x;hlt1;hlt2", {{1000, -1., 1.}, {1000, -1., 1.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_Y{ - this, "y_stAtVtx", "stAtVtx_y;hlt1;hlt2", {{1000, -1., 1.}, {1000, -1., 1.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_Z{ - this, "z_stAtVtx", "stAtVtx_z;hlt1;hlt2", {{1000, -200., 200.}, {1000, -200., 200.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_TX{ - this, "tx_stAtVtx", "stAtVtx_tx;hlt1;hlt2", {{1000, -0.3, 0.3}, {1000, -0.3, 0.3}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_TY{ - this, "ty_stAtVtx", "stAtVtx_ty;hlt1;hlt2", {{1000, -0.2, 0.2}, {1000, -0.2, 0.2}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_QOP{ - this, "qop_stAtVtx", "stAtVtx_qop;hlt1;hlt2", {{1000, -0.001, 0.001}, {1000, -0.001, 0.001}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov00{ - this, "cov00_stAtVtx", "stAtVtxcov00;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov01{ - this, "cov01_stAtVtx", "stAtVtxcov01;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov02{ - this, "cov02_stAtVtx", "stAtVtxcov02;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov03{ - this, "cov03_stAtVtx", "stAtVtxcov03;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov04{ - this, "cov04_stAtVtx", "stAtVtxcov04;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov11{ - this, "cov11_stAtVtx", "stAtVtxcov11;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov12{ - this, "cov12_stAtVtx", "stAtVtxcov12;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov13{ - this, "cov13_stAtVtx", "stAtVtxcov13;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov14{ - this, "cov14_stAtVtx", "stAtVtxcov14;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov22{ - this, "cov22_stAtVtx", "stAtVtxcov22;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov23{ - this, "cov23_stAtVtx", "stAtVtxcov23;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov24{ - this, "cov24_stAtVtx", "stAtVtxcov24;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov33{ - this, "cov33_stAtVtx", "stAtVtxcov33;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov34{ - this, "cov34_stAtVtx", "stAtVtxcov34;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - - mutable Gaudi::Accumulators::Histogram<2> m_stAtVtx_cov44{ - this, "cov44_stAtVtx", "stAtVtxcov44;hlt1;hlt2", {{1000, -100., 100.}, {1000, -100., 100.}}}; - - mutable Gaudi::Accumulators::Histogram<2, Gaudi::Accumulators::atomicity::full, int> m_nVPHits{ - this, "nVPHits", "nVPHits;hlt1;hlt2", {{20, 0, 20}, {20, 0, 20}}}; - mutable Gaudi::Accumulators::Histogram<2, Gaudi::Accumulators::atomicity::full, int> m_nUTHits{ - this, "nUTHits", "nUTHits;hlt1;hlt2", {{4, 0, 4}, {4, 0, 4}}}; - mutable Gaudi::Accumulators::Histogram<2, Gaudi::Accumulators::atomicity::full, int> m_nFTHits{ - this, "nFTHits", "nFTHits;hlt1;hlt2", {{4, 8, 12}, {4, 8, 12}}}; -}; - -DECLARE_COMPONENT( TriggerObjectsCompatibilityChecker ) -void TriggerObjectsCompatibilityChecker:: - operator()( const LHCb::Tracks& hlt1_tracks, const LHCb::Tracks& hlt2_tracks, const LHCb::MCParticles& mcparticles, - const LHCb::LinksByKey& linker_hlt1, const LHCb::LinksByKey& linker_hlt2, const LHCb::RecVertices& hlt1_pvs, - const LHCb::RecVertices& hlt2_pvs ) const { - - std::map truepvs; - - for ( const auto& hlt2_track : hlt2_tracks ) { - - const LHCb::MCParticle* hlt2_mcparticle = mcTruth( *hlt2_track, mcparticles, linker_hlt2 ); // CheckerBase - - auto fit = fitResult( *hlt2_track ); - int trackWasFitted = fit && !fit->nodes().empty(); - - ++m_fitted_track[trackWasFitted]; - - if ( !selected( hlt2_mcparticle ) ) continue; - if ( !( hlt2_track->hasVelo() ) ) continue; - if ( int( hlt2_track->type() ) != 3 ) continue; - - const auto& hlt2_state = trackState( *hlt2_track ); - double hlt2_x = hlt2_state.x(); - double hlt2_y = hlt2_state.y(); - double hlt2_z = hlt2_state.z(); - double hlt2_tx = hlt2_state.tx(); - double hlt2_ty = hlt2_state.ty(); - double hlt2_qop = hlt2_state.qOverP(); - double hlt2_pt = hlt2_state.pt(); - - ++m_numPVs_hlt2[int( hlt2_pvs.size() )]; - ++m_numPVs_hlt1[int( hlt1_pvs.size() )]; - ++m_nTracks[{int( hlt1_tracks.size() ), int( hlt2_tracks.size() )}]; - - LHCb::CubicStateInterpolationTraj hlt2_tracktraj( hlt2_state, Gaudi::XYZVector() ); - Gaudi::XYZPoint hlt2_trkpos( hlt2_state.position() ); - Gaudi::XYZVector hlt2_trkdir( hlt2_state.slopes().Unit() ); - - const LHCb::RecVertex* hlt2_recpv( nullptr ); - double hlt2_bestip2( 10000 ); // is this correct? - - for ( const auto& thispv : hlt2_pvs ) { - Gaudi::XYZVector dx = hlt2_trkpos - thispv->position(); - Gaudi::XYZVector delta = dx - hlt2_trkdir * dx.Dot( hlt2_trkdir ); - double ip2 = delta.Mag2(); - if ( hlt2_recpv == 0 || ip2 < hlt2_bestip2 ) { - hlt2_bestip2 = ip2; - hlt2_recpv = thispv; - } - } - - for ( const auto& hlt1_track : hlt1_tracks ) { - const LHCb::MCParticle* hlt1_mcparticle = mcTruth( *hlt1_track, mcparticles, linker_hlt1 ); - - auto fit = fitResult( *hlt1_track ); - int hlt1_trackWasFitted = fit && !fit->nodes().empty(); - - ++m_hlt1_fitted_track[hlt1_trackWasFitted]; - - if ( !selected( hlt1_mcparticle ) ) continue; - - if ( hlt2_mcparticle == hlt1_mcparticle ) { - const auto& hlt1_state = trackState( *hlt1_track ); - double hlt1_x = hlt1_state.x(); - double hlt1_y = hlt1_state.y(); - double hlt1_z = hlt1_state.z(); - double hlt1_tx = hlt1_state.tx(); - double hlt1_ty = hlt1_state.ty(); - double hlt1_qop = hlt1_state.qOverP(); - double hlt1_pt = hlt1_state.pt(); - - LHCb::CubicStateInterpolationTraj hlt1_tracktraj( hlt1_state, Gaudi::XYZVector() ); - Gaudi::XYZPoint hlt1_trkpos( hlt1_state.position() ); - Gaudi::XYZVector hlt1_trkdir( hlt1_state.slopes().Unit() ); - const LHCb::RecVertex* hlt1_recpv( nullptr ); - double hlt1_bestip2( 10000 ); // is this correct? - - for ( const auto& thispv : hlt1_pvs ) { - Gaudi::XYZVector dx = hlt1_trkpos - thispv->position(); - Gaudi::XYZVector delta = dx - hlt1_trkdir * dx.Dot( hlt1_trkdir ); - double ip2 = delta.Mag2(); - if ( hlt1_recpv == 0 || ip2 < hlt1_bestip2 ) { - hlt1_bestip2 = ip2; - hlt1_recpv = thispv; - } - } - - ++m_probChi2_hlt1[hlt1_track->probChi2()]; - ++m_probChi2_hlt2[hlt2_track->probChi2()]; - - ++m_chi2_hlt1[hlt1_track->chi2()]; - ++m_chi2_hlt2[hlt2_track->chi2()]; - - ++m_ndof_hlt1[hlt1_track->nDoF()]; - ++m_ndof_hlt2[hlt2_track->nDoF()]; - - ++m_chi2_ndof_hlt1[hlt1_track->chi2() / hlt1_track->nDoF()]; - ++m_chi2_ndof_hlt2[hlt2_track->chi2() / hlt2_track->nDoF()]; - - ++m_type_hlt1[int( hlt1_track->type() )]; - ++m_type_hlt2[int( hlt2_track->type() )]; - - ++m_deltaX[hlt2_x - hlt1_x]; - ++m_deltaY[hlt2_y - hlt1_y]; - ++m_deltaZ[hlt2_z - hlt1_z]; - ++m_deltaTX[hlt2_tx - hlt1_tx]; - ++m_deltaTY[hlt2_ty - hlt1_ty]; - ++m_deltaQOP[hlt2_qop - hlt1_qop]; - - ++m_2d_X[{hlt1_x, hlt2_x}]; - ++m_2d_Y[{hlt1_y, hlt2_y}]; - ++m_2d_Z[{hlt1_z, hlt2_z}]; - ++m_2d_TX[{hlt1_tx, hlt2_tx}]; - ++m_2d_TY[{hlt1_ty, hlt2_ty}]; - ++m_2d_QOP[{hlt1_qop, hlt2_qop}]; - - ++m_2d_PT[{hlt1_pt, hlt2_pt}]; - ++m_nVPHits[{hlt1_track->nVPHits(), hlt2_track->nVPHits()}]; - ++m_nUTHits[{hlt1_track->nUTHits(), hlt2_track->nUTHits()}]; - ++m_nFTHits[{hlt1_track->nFTHits(), hlt2_track->nFTHits()}]; - - if ( hlt2_recpv && hlt1_recpv ) { - - LHCb::State hlt2_stateAtVtx = hlt2_tracktraj.state( hlt2_recpv->position().z() ); - LHCb::State hlt1_stateAtVtx = hlt1_tracktraj.state( hlt1_recpv->position().z() ); - - double hlt2_tx = hlt2_stateAtVtx.tx(); - double hlt2_recipxerr = std::sqrt( hlt2_state.covariance()( 0, 0 ) + hlt2_recpv->covMatrix()( 0, 0 ) + - 2 * hlt2_tx * hlt2_recpv->covMatrix()( 0, 2 ) + - hlt2_tx * hlt2_tx * hlt2_recpv->covMatrix()( 2, 2 ) ); - double hlt2_ty = hlt2_stateAtVtx.ty(); - double hlt2_recipyerr = std::sqrt( hlt2_state.covariance()( 1, 1 ) + hlt2_recpv->covMatrix()( 1, 1 ) + - 2 * hlt2_ty * hlt2_recpv->covMatrix()( 1, 2 ) + - hlt2_ty * hlt2_ty * hlt2_recpv->covMatrix()( 2, 2 ) ); - - double hlt1_tx = hlt1_stateAtVtx.tx(); - double hlt1_recipxerr = std::sqrt( hlt1_state.covariance()( 0, 0 ) + hlt1_recpv->covMatrix()( 0, 0 ) + - 2 * hlt1_tx * hlt1_recpv->covMatrix()( 0, 2 ) + - hlt1_tx * hlt1_tx * hlt1_recpv->covMatrix()( 2, 2 ) ); - double hlt1_ty = hlt1_stateAtVtx.ty(); - - double hlt1_recipyerr = std::sqrt( hlt1_state.covariance()( 1, 1 ) + hlt1_recpv->covMatrix()( 1, 1 ) + - 2 * hlt1_ty * hlt1_recpv->covMatrix()( 1, 2 ) + - hlt1_ty * hlt1_ty * hlt1_recpv->covMatrix()( 2, 2 ) ); - - ++m_deltaIP3D[std::sqrt( hlt2_bestip2 ) - std::sqrt( hlt1_bestip2 )]; - ++m_deltaIPx[( hlt2_stateAtVtx.x() - hlt2_recpv->position().x() ) - - ( hlt1_stateAtVtx.x() - hlt1_recpv->position().x() )]; - ++m_deltaIPy[( hlt2_stateAtVtx.y() - hlt2_recpv->position().y() ) - - ( hlt1_stateAtVtx.y() - hlt1_recpv->position().y() )]; - - ++m_deltaIPxErr[hlt2_recipxerr - hlt1_recipxerr]; - ++m_deltaIPyErr[hlt2_recipyerr - hlt1_recipyerr]; - - ++m_IP3D[{std::sqrt( hlt1_bestip2 ), std::sqrt( hlt2_bestip2 )}]; - ++m_IPx[{( hlt1_stateAtVtx.x() - hlt1_recpv->position().x() ), - ( hlt2_stateAtVtx.x() - hlt2_recpv->position().x() )}]; - ++m_IPy[{( hlt1_stateAtVtx.y() - hlt1_recpv->position().y() ), - ( hlt2_stateAtVtx.y() - hlt2_recpv->position().y() )}]; - - ++m_IPxErr[{hlt1_recipxerr, hlt2_recipxerr}]; - ++m_IPyErr[{hlt1_recipyerr, hlt2_recipyerr}]; - - auto hlt2_bestipchi2 = - LHCb::TrackVertexUtils::vertexChi2( hlt2_stateAtVtx, hlt2_recpv->position(), hlt2_recpv->covMatrix() ); - auto hlt1_bestipchi2 = - LHCb::TrackVertexUtils::vertexChi2( hlt1_stateAtVtx, hlt1_recpv->position(), hlt1_recpv->covMatrix() ); - - ++m_deltaIPChi2[hlt2_bestipchi2 - hlt1_bestipchi2]; - ++m_IPChi2[{hlt1_bestipchi2, hlt2_bestipchi2}]; - ++m_IPChi2_zoom[{hlt1_bestipchi2, hlt2_bestipchi2}]; - ++m_deltaIPChi2[( hlt1_bestipchi2 - hlt2_bestipchi2 )]; - - ++m_nTracksPerPV[{hlt1_recpv->tracks().size(), hlt2_recpv->tracks().size()}]; - - ++m_PV_X[{hlt1_recpv->position().x(), hlt2_recpv->position().x()}]; - ++m_PV_Y[{hlt1_recpv->position().y(), hlt2_recpv->position().y()}]; - ++m_PV_Z[{hlt1_recpv->position().z(), hlt2_recpv->position().z()}]; - - ++m_PV_cov00[{hlt1_recpv->covMatrix()( 0, 0 ), hlt2_recpv->covMatrix()( 0, 0 )}]; - ++m_PV_cov01[{hlt1_recpv->covMatrix()( 0, 1 ), hlt2_recpv->covMatrix()( 0, 1 )}]; - ++m_PV_cov02[{hlt1_recpv->covMatrix()( 0, 2 ), hlt2_recpv->covMatrix()( 0, 2 )}]; - - ++m_PV_cov11[{hlt1_recpv->covMatrix()( 1, 1 ), hlt2_recpv->covMatrix()( 1, 1 )}]; - ++m_PV_cov12[{hlt1_recpv->covMatrix()( 1, 2 ), hlt2_recpv->covMatrix()( 1, 2 )}]; - - ++m_PV_cov22[{hlt1_recpv->covMatrix()( 2, 2 ), hlt2_recpv->covMatrix()( 2, 2 )}]; - - ++m_stAtVtx_X[{hlt1_stateAtVtx.position().x(), hlt2_stateAtVtx.position().x()}]; - ++m_stAtVtx_Y[{hlt1_stateAtVtx.position().y(), hlt2_stateAtVtx.position().y()}]; - ++m_stAtVtx_Z[{hlt1_stateAtVtx.position().z(), hlt2_stateAtVtx.position().z()}]; - ++m_stAtVtx_TX[{hlt1_stateAtVtx.tx(), hlt2_stateAtVtx.tx()}]; - ++m_stAtVtx_TY[{hlt1_stateAtVtx.ty(), hlt2_stateAtVtx.ty()}]; - ++m_stAtVtx_QOP[{hlt1_stateAtVtx.qOverP(), hlt2_stateAtVtx.qOverP()}]; - - ++m_stAtVtx_cov00[{hlt1_stateAtVtx.posMomCovariance()( 0, 0 ), hlt2_stateAtVtx.posMomCovariance()( 0, 0 )}]; - ++m_stAtVtx_cov01[{hlt1_stateAtVtx.posMomCovariance()( 0, 1 ), hlt2_stateAtVtx.posMomCovariance()( 0, 1 )}]; - ++m_stAtVtx_cov02[{hlt1_stateAtVtx.posMomCovariance()( 0, 2 ), hlt2_stateAtVtx.posMomCovariance()( 0, 2 )}]; - ++m_stAtVtx_cov03[{hlt1_stateAtVtx.posMomCovariance()( 0, 3 ), hlt2_stateAtVtx.posMomCovariance()( 0, 3 )}]; - ++m_stAtVtx_cov04[{hlt1_stateAtVtx.posMomCovariance()( 0, 4 ), hlt2_stateAtVtx.posMomCovariance()( 0, 4 )}]; - - ++m_stAtVtx_cov11[{hlt1_stateAtVtx.posMomCovariance()( 1, 1 ), hlt2_stateAtVtx.posMomCovariance()( 1, 1 )}]; - ++m_stAtVtx_cov12[{hlt1_stateAtVtx.posMomCovariance()( 1, 2 ), hlt2_stateAtVtx.posMomCovariance()( 1, 2 )}]; - ++m_stAtVtx_cov13[{hlt1_stateAtVtx.posMomCovariance()( 1, 3 ), hlt2_stateAtVtx.posMomCovariance()( 1, 3 )}]; - ++m_stAtVtx_cov14[{hlt1_stateAtVtx.posMomCovariance()( 1, 4 ), hlt2_stateAtVtx.posMomCovariance()( 1, 4 )}]; - - ++m_stAtVtx_cov22[{hlt1_stateAtVtx.posMomCovariance()( 2, 2 ), hlt2_stateAtVtx.posMomCovariance()( 2, 2 )}]; - ++m_stAtVtx_cov23[{hlt1_stateAtVtx.posMomCovariance()( 2, 3 ), hlt2_stateAtVtx.posMomCovariance()( 2, 3 )}]; - ++m_stAtVtx_cov24[{hlt1_stateAtVtx.posMomCovariance()( 2, 4 ), hlt2_stateAtVtx.posMomCovariance()( 2, 4 )}]; - - ++m_stAtVtx_cov33[{hlt1_stateAtVtx.posMomCovariance()( 3, 3 ), hlt2_stateAtVtx.posMomCovariance()( 3, 3 )}]; - ++m_stAtVtx_cov34[{hlt1_stateAtVtx.posMomCovariance()( 3, 4 ), hlt2_stateAtVtx.posMomCovariance()( 3, 4 )}]; - - ++m_stAtVtx_cov44[{hlt1_stateAtVtx.posMomCovariance()( 4, 4 ), hlt2_stateAtVtx.posMomCovariance()( 4, 4 )}]; - - } // if reconstructd pv in both cases - else { - - ++m_deltaIP3D[-999999.]; - ++m_deltaIPx[-999999.]; - ++m_deltaIPy[-999999.]; - - ++m_IP3D[{-999999., -999999.}]; - ++m_IPx[{-999999., -999999.}]; - ++m_IPy[{-999999., -999999.}]; - - ++m_IPxErr[{-999999., -999999.}]; - ++m_IPyErr[{-999999., -999999.}]; - - ++m_deltaIPChi2[-999999.]; - ++m_IPChi2[{-999999., -999999.}]; - ++m_IPChi2_zoom[{-999999., -999999.}]; - ++m_deltaIPChi2[-99999.]; - - ++m_nTracksPerPV[{-9999., -9999.}]; - ++m_PV_X[{-9999., -9999.}]; - ++m_PV_Y[{-9999., -9999.}]; - ++m_PV_Z[{-9999., -9999.}]; - - ++m_PV_cov00[{-9999., -9999.}]; - ++m_PV_cov01[{-9999., -9999.}]; - ++m_PV_cov02[{-9999., -9999.}]; - - ++m_PV_cov11[{-9999., -9999.}]; - ++m_PV_cov12[{-9999., -9999.}]; - - ++m_PV_cov22[{-9999., -9999.}]; - - ++m_stAtVtx_X[{-9999., -9999.}]; - ++m_stAtVtx_Y[{-9999., -9999.}]; - ++m_stAtVtx_Z[{-9999., -9999.}]; - ++m_stAtVtx_TX[{-9999., -9999.}]; - ++m_stAtVtx_TY[{-9999., -9999.}]; - ++m_stAtVtx_QOP[{-9999., -9999.}]; - - ++m_stAtVtx_cov00[{-9999., -9999.}]; - ++m_stAtVtx_cov01[{-9999., -9999.}]; - ++m_stAtVtx_cov02[{-9999., -9999.}]; - ++m_stAtVtx_cov03[{-9999., -9999.}]; - ++m_stAtVtx_cov04[{-9999., -9999.}]; - - ++m_stAtVtx_cov11[{-9999., -9999.}]; - ++m_stAtVtx_cov12[{-9999., -9999.}]; - ++m_stAtVtx_cov13[{-9999., -9999.}]; - ++m_stAtVtx_cov14[{-9999., -9999.}]; - - ++m_stAtVtx_cov22[{-9999., -9999.}]; - ++m_stAtVtx_cov23[{-9999., -9999.}]; - ++m_stAtVtx_cov24[{-9999., -9999.}]; - - ++m_stAtVtx_cov33[{-9999., -9999.}]; - ++m_stAtVtx_cov34[{-9999., -9999.}]; - - ++m_stAtVtx_cov44[{-9999., -9999.}]; - - } // if not resonctructed pvs - - // We also want to monitor the reconstructed IP, so the IP with respect to the reconstructed PVs. - - } // if mcparticles are the same for hlt1 and hlt2 - } // loop over hlt1 tracks - - } // loop over hlt2 tracks - - // const LHCb::MCParticle* mcparticle{nullptr}; - // double maxWeight{0}; - // linker.applyToLinks( - // track->key(), [&maxWeight, &mcparticle, &mcparticles]( unsigned int, unsigned int mcPartKey, float weight ) { - // if ( weight > maxWeight ) { - // / maxWeight = weight; - // mcparticle = static_cast( mcparticles.containedObject( mcPartKey ) ); - // } - //} ); -}; diff --git a/Tr/TrackCheckers/src/TriggerObjectsCompatibilityProfileChecker.cpp b/Tr/TrackCheckers/src/TriggerObjectsCompatibilityProfileChecker.cpp index 033c1dab6392b0f36e3d09dc6cca0284adf449a2..f256ee62da05e0c19752de2e5497db1b9d6d1d6c 100644 --- a/Tr/TrackCheckers/src/TriggerObjectsCompatibilityProfileChecker.cpp +++ b/Tr/TrackCheckers/src/TriggerObjectsCompatibilityProfileChecker.cpp @@ -14,7 +14,7 @@ #include "Event/RecVertex.h" #include "Event/Track.h" #include "Event/TrackVertexUtils.h" -#include "Gaudi/Accumulators/Histogram.h" +#include "Gaudi/Accumulators/StaticHistogram.h" #include "LHCbAlgs/Consumer.h" #include "TrackCheckerBase.h" #include "TrackKernel/CubicStateInterpolationTraj.h" @@ -40,154 +40,168 @@ public: const LHCb::LinksByKey& linker_hlt2, const LHCb::RecVertices& hlt1_pvs, const LHCb::RecVertices& hlt2_pvs ) const override; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_hlt2_fitted_track{ + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_hlt2_fitted_track{ this, "hlt2_fitted_track", "hlt2 fitted track", {2, 0, 2}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_hlt1_fitted_track{ + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_hlt1_fitted_track{ this, "hlt1_fitted_track", "hlt1 fitted track", {2, 0, 2}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_numPVs_hlt1{ + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_numPVs_hlt1{ this, "number_reconstructed_PVs_hlt1", "number_reconstructed_PVs_hlt1", {20, 0, 10}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_numPVs_hlt2{ + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_numPVs_hlt2{ this, "number_reconstructed_PVs_hlt2", "number_reconstructed_PVs_hlt2", {20, 0, 10}}; - mutable Gaudi::Accumulators::Histogram<1> m_probChi2_hlt1{this, "probChi2_hlt1", "probChi2_hlt1", {50, 0, 1.}}; - mutable Gaudi::Accumulators::Histogram<1> m_probChi2_hlt2{this, "probChi2_hlt2", "probChi2_hlt2", {50, 0, 1.}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_probChi2_hlt1{this, "probChi2_hlt1", "probChi2_hlt1", {50, 0, 1.}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_probChi2_hlt2{this, "probChi2_hlt2", "probChi2_hlt2", {50, 0, 1.}}; - mutable Gaudi::Accumulators::Histogram<1> m_chi2_hlt1{this, "chi2_hlt1", "chi2_hlt1", {100, 0, 200}}; - mutable Gaudi::Accumulators::Histogram<1> m_chi2_hlt2{this, "chi2_hlt2", "chi2_hlt2", {100, 0, 200}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_chi2_hlt1{this, "chi2_hlt1", "chi2_hlt1", {100, 0, 200}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_chi2_hlt2{this, "chi2_hlt2", "chi2_hlt2", {100, 0, 200}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_ndof_hlt1{ + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_ndof_hlt1{ this, "ndof_hlt1", "ndof_hlt1", {30, 0, 30}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_ndof_hlt2{ + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_ndof_hlt2{ this, "ndof_hlt2", "ndof_hlt2", {30, 0, 30}}; - mutable Gaudi::Accumulators::Histogram<1> m_chi2_ndof_hlt1{this, "chi2_ndof_hlt1", "chi2_ndof_hlt1", {100, 0, 50}}; - mutable Gaudi::Accumulators::Histogram<1> m_chi2_ndof_hlt2{this, "chi2_ndof_hlt2", "chi2_ndof_hlt2", {100, 0, 50}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_chi2_ndof_hlt1{ + this, "chi2_ndof_hlt1", "chi2_ndof_hlt1", {100, 0, 50}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_chi2_ndof_hlt2{ + this, "chi2_ndof_hlt2", "chi2_ndof_hlt2", {100, 0, 50}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_type_hlt1{ + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_type_hlt1{ this, "track_type_hlt1", "track_type_hlt1", {10, 0, 10}}; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_type_hlt2{ + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_type_hlt2{ this, "track_type_hlt2", "track_type_hlt2", {10, 0, 10}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaX{ + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaX{ this, "deltaX_first_state", "hlt2 - hlt1 x firstate", {100, -0.2, 0.2}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaY{ + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaY{ this, "deltaY_first_state", "hlt2 - hlt1 y firstate", {100, -0.2, 0.2}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaZ{ + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaZ{ this, "deltaZ_first_state", "hlt2 - hlt1 z firstate", {100, -2., 2.}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaTX{ + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaTX{ this, "deltaTX_first_state", "hlt2 - hlt1 tx firstate", {100, -0.005, 0.005}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaTY{ + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaTY{ this, "deltaTY_first_state", "hlt2 - hlt1 ty firstate", {100, -0.005, 0.005}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaQOP{ + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaQOP{ this, "deltaQOP_first_state", "hlt2 - hlt1 qop firstate", {100, -0.00005, 0.00005}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_2d_X{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_2d_X{ this, "x_first_state_2d", "firstate_x;hlt1;hlt2", {100, -1., 1.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_2d_Y{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_2d_Y{ this, "y_first_state_2d", "firstate_y;hlt1;hlt2", {100, -1., 1.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_2d_Z{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_2d_Z{ this, "z_first_state_2d", "firstate_z;hlt1;hlt2", {100, -200., 600.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_2d_TX{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_2d_TX{ this, "tx_first_state_2d", "firstate_tx;hlt1;hlt2", {100, -0.3, 0.3}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_2d_TY{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_2d_TY{ this, "ty_first_state_2d", "firstate_ty;hlt1;hlt2", {100, -0.2, 0.2}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_2d_QOP{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_2d_QOP{ this, "qop_first_state_2d", "firstate_qop;hlt1;hlt2", {100, -0.001, 0.001}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_2d_PT{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_2d_PT{ this, "pt_first_state_2d", "firstate_pt;hlt1;hlt2", {100, 0., 10000.}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaIP3D{this, "recdeltaIP3D", "delta_recIP3D", {100, -0.5, 0.5}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaIPx{this, "recdeltaIPx", "delta_recIPx", {100, -0.5, 0.5}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaIPy{this, "recdeltaIPy", "delta_recIPy", {100, -0.5, 0.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaIP3D{this, "recdeltaIP3D", "delta_recIP3D", {100, -0.5, 0.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaIPx{this, "recdeltaIPx", "delta_recIPx", {100, -0.5, 0.5}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaIPy{this, "recdeltaIPy", "delta_recIPy", {100, -0.5, 0.5}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_IP3D{this, "recIP3D", "recIP3D;hlt1;hlt2", {100, 0., 1.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_IPx{this, "recIPx", "recIPx;hlt1;hlt2", {100, -1., 1.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_IPy{this, "recIPy", "recIPy;hlt1;hlt2", {100, -1., 1.}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_IP3D{this, "recIP3D", "recIP3D;hlt1;hlt2", {100, 0., 1.}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_IPx{this, "recIPx", "recIPx;hlt1;hlt2", {100, -1., 1.}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_IPy{this, "recIPy", "recIPy;hlt1;hlt2", {100, -1., 1.}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaIPxErr{this, "delta_recIPxErr", "delta_recIPxErr", {100, -0.2, 0.2}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaIPyErr{this, "delta_recIPyErr", "delta_recIPyErr", {100, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaIPxErr{ + this, "delta_recIPxErr", "delta_recIPxErr", {100, -0.2, 0.2}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaIPyErr{ + this, "delta_recIPyErr", "delta_recIPyErr", {100, -0.2, 0.2}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_IPxErr{this, "recIPxErr", "recIPxErr;hlt1;hlt2", {100, 0., 0.5}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_IPyErr{this, "recIPyErr", "recIPyErr;hlt1;hlt2", {100, 0., 0.5}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_IPxErr{ + this, "recIPxErr", "recIPxErr;hlt1;hlt2", {100, 0., 0.5}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_IPyErr{ + this, "recIPyErr", "recIPyErr;hlt1;hlt2", {100, 0., 0.5}}; - mutable Gaudi::Accumulators::Histogram<1> m_deltaIPChi2{this, "delta_recIPChi2", "delta_recIPChi2", {100, -10., 10.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_IPChi2{this, "recIPChi2", "recIPChi2;hlt1;hlt2", {40, 0., 20.}}; + mutable Gaudi::Accumulators::StaticHistogram<1> m_deltaIPChi2{ + this, "delta_recIPChi2", "delta_recIPChi2", {100, -10., 10.}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_IPChi2{ + this, "recIPChi2", "recIPChi2;hlt1;hlt2", {40, 0., 20.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_IPChi2_zoom{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_IPChi2_zoom{ this, "recIPChi2_zoom", "recIPChi2_zoom;hlt1;hlt2", {20, 0., 2.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_nTracksPerPV{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_nTracksPerPV{ this, "nTracksPerPV", "nTracksPV;hlt1;hlt2", {50, 0, 50}}; - mutable Gaudi::Accumulators::ProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_nTracks{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_nTracks{ this, "nTracks", "nTracks;hlt1;hlt2", {50, 0, 100}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_PV_X{this, "x_PV", "PVx;hlt1;hlt2", {100, -0.5, 0.5}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_PV_Y{this, "y_PV", "PVy;hlt1;hlt2", {100, -0.5, 0.5}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_PV_Z{this, "z_PV", "PVz;hlt1;hlt2", {100, -200., 200.}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_PV_X{this, "x_PV", "PVx;hlt1;hlt2", {100, -0.5, 0.5}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_PV_Y{this, "y_PV", "PVy;hlt1;hlt2", {100, -0.5, 0.5}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_PV_Z{this, "z_PV", "PVz;hlt1;hlt2", {100, -200., 200.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_PV_cov00{this, "cov00_PV", "PVcov00;hlt1;hlt2", {50, 0., 0.0001}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_PV_cov01{this, "cov01_PV", "PVcov01;hlt1;hlt2", {50, -0.05, 0.05}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_PV_cov02{this, "cov02_PV", "PVcov02;hlt1;hlt2", {50, -0.05, 0.05}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_PV_cov00{ + this, "cov00_PV", "PVcov00;hlt1;hlt2", {50, 0., 0.0001}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_PV_cov01{ + this, "cov01_PV", "PVcov01;hlt1;hlt2", {50, -0.05, 0.05}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_PV_cov02{ + this, "cov02_PV", "PVcov02;hlt1;hlt2", {50, -0.05, 0.05}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_PV_cov11{this, "cov11_PV", "PVcov11;hlt1;hlt2", {50, -0.1, 0.1}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_PV_cov12{this, "cov12_PV", "PVcov12;hlt1;hlt2", {50, -0.1, 0.1}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_PV_cov11{ + this, "cov11_PV", "PVcov11;hlt1;hlt2", {50, -0.1, 0.1}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_PV_cov12{ + this, "cov12_PV", "PVcov12;hlt1;hlt2", {50, -0.1, 0.1}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_PV_cov22{this, "cov22_PV", "PVcov12;hlt1;hlt2", {50, 0., 1.}}; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_PV_cov22{ + this, "cov22_PV", "PVcov12;hlt1;hlt2", {50, 0., 1.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_X{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_X{ this, "x_stAtVtx", "stAtVtx_x;hlt1;hlt2", {100, -1., 1.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_Y{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_Y{ this, "y_stAtVtx", "stAtVtx_y;hlt1;hlt2", {100, -1., 1.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_Z{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_Z{ this, "z_stAtVtx", "stAtVtx_z;hlt1;hlt2", {100, -200., 200.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_TX{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_TX{ this, "tx_stAtVtx", "stAtVtx_tx;hlt1;hlt2", {100, -0.3, 0.3}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_TY{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_TY{ this, "ty_stAtVtx", "stAtVtx_ty;hlt1;hlt2", {100, -0.2, 0.2}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_QOP{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_QOP{ this, "qop_stAtVtx", "stAtVtx_qop;hlt1;hlt2", {100, -0.001, 0.001}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov00{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov00{ this, "cov00_stAtVtx", "stAtVtxcov00;hlt1;hlt2", {50, 0., 0.01}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov01{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov01{ this, "cov01_stAtVtx", "stAtVtxcov01;hlt1;hlt2", {50, -1, 10.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov02{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov02{ this, "cov02_stAtVtx", "stAtVtxcov02;hlt1;hlt2", {50, -500., 500.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov03{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov03{ this, "cov03_stAtVtx", "stAtVtxcov03;hlt1;hlt2", {50, -10., 10.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov04{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov04{ this, "cov04_stAtVtx", "stAtVtxcov04;hlt1;hlt2", {50, -10., 10.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov11{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov11{ this, "cov11_stAtVtx", "stAtVtxcov11;hlt1;hlt2", {50, 0., 10.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov12{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov12{ this, "cov12_stAtVtx", "stAtVtxcov12;hlt1;hlt2", {50, -500., 500.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov13{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov13{ this, "cov13_stAtVtx", "stAtVtxcov13;hlt1;hlt2", {50, -1., 1.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov14{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov14{ this, "cov14_stAtVtx", "stAtVtxcov14;hlt1;hlt2", {50, -1., 0.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov22{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov22{ this, "cov22_stAtVtx", "stAtVtxcov22;hlt1;hlt2", {50, -500., 500.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov23{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov23{ this, "cov23_stAtVtx", "stAtVtxcov23;hlt1;hlt2", {50, -500., 500.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov24{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov24{ this, "cov24_stAtVtx", "stAtVtxcov24;hlt1;hlt2", {50, -500., 500.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov33{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov33{ this, "cov33_stAtVtx", "stAtVtxcov33;hlt1;hlt2", {50, 0., 500.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov34{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov34{ this, "cov34_stAtVtx", "stAtVtxcov34;hlt1;hlt2", {50, -1000., 1000.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_stAtVtx_cov44{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_stAtVtx_cov44{ this, "cov44_stAtVtx", "stAtVtxcov44;hlt1;hlt2", {50, 0., 800.}}; - mutable Gaudi::Accumulators::ProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_nVPHits{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_nVPHits{ this, "nVPHits", "nVPHits;hlt1;hlt2", {21, 0, 21}}; - mutable Gaudi::Accumulators::ProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_nUTHits{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_nUTHits{ this, "nUTHits", "nUTHits;hlt1;hlt2", {4, 0, 4}}; - mutable Gaudi::Accumulators::ProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_nFTHits{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_nFTHits{ this, "nFTHits", "nFTHits;hlt1;hlt2", {16, 0, 16}}; }; // Difference between LHCb::Tracks and LHCb::Track::Range? diff --git a/Tr/TrackMonitors/CMakeLists.txt b/Tr/TrackMonitors/CMakeLists.txt index 8e0f8a149733d534455216190f9461b3c9d2b6e5..18ae5e673a8b8b3584d01a079f8b47b070d8af06 100644 --- a/Tr/TrackMonitors/CMakeLists.txt +++ b/Tr/TrackMonitors/CMakeLists.txt @@ -24,16 +24,13 @@ gaudi_add_module(TrackMonitors src/MonitorDetectorCorrelations.cpp src/MonitorECALEnergyRawEventSizeCorrelations.cpp src/TrackCaloMatchMonitor.cpp - src/TrackDiMuonMonitor.cpp src/TrackFitMatchMonitor.cpp src/TrackMonitor.cpp - src/TrackMonitorBase.cpp src/TrackMonitorTupleBase.cpp src/TrackMuonMatchMonitor.cpp src/TrackPV2HalfMonitor.cpp src/TrackParticleMonitor.cpp src/TrackTune.cpp - src/TrackV0Monitor.cpp src/TrackVertexMonitor.cpp src/TrackVPOverlapMonitor.cpp src/UTTrackMonitor.cpp diff --git a/Tr/TrackMonitors/src/TrackDiMuonMonitor.cpp b/Tr/TrackMonitors/src/TrackDiMuonMonitor.cpp deleted file mode 100644 index eecd9e073d93f4c359817f190c1403e9bae9b8eb..0000000000000000000000000000000000000000 --- a/Tr/TrackMonitors/src/TrackDiMuonMonitor.cpp +++ /dev/null @@ -1,220 +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/DetectorElement.h" -#include "DetDesc/GenericConditionAccessorHolder.h" -#include "Event/TwoProngVertex.h" -#include "Kernel/IParticlePropertySvc.h" -#include "Kernel/ITrajPoca.h" -#include "Kernel/ParticleProperty.h" -#include "Magnet/DeMagnet.h" -#include "TrackInterfaces/ITrackVertexer.h" -#include "TrackKernel/TrackTraj.h" - -#include "GaudiAlg/GaudiHistoAlg.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "LHCbAlgs/Transformer.h" - -#include "gsl/gsl_cdf.h" - -namespace { - //============================================================================= - // Check of two tracks do not have two much overlap to be combined as TwoProng - //============================================================================= - - bool overlap( const LHCb::Track& trackA, const LHCb::Track& trackB ) { - // for now, just look at common ancestors - return std::any_of( trackA.ancestors().begin(), trackA.ancestors().end(), [&]( const auto& ancA ) { - return std::any_of( trackB.ancestors().begin(), trackB.ancestors().end(), - [&]( const auto& ancB ) { return ancA == ancB; } ); - } ); - } - -} // namespace -using TwoProngVertices = KeyedContainer; - -class TrackDiMuonMonitor final - : public LHCb::Algorithm::Transformer< - TwoProngVertices( LHCb::Tracks const&, DetectorElement const&, DeMagnet const& ), - LHCb::Algorithm::Traits::usesBaseAndConditions> { -public: - /** Standard construtor */ - TrackDiMuonMonitor( const std::string& name, ISvcLocator* pSvcLocator ) - : Transformer{name, - pSvcLocator, - {KeyValue{"TrackLocation", LHCb::TrackLocation::Muon}, - KeyValue{"StandardGeometryTop", LHCb::standard_geometry_top}, - KeyValue{"Magnet", LHCb::Det::Magnet::det_path}}, - {"DiMuonLocation", "Rec/Vertex/DiMuon"}} {} - - /** Algorithm initialize */ - StatusCode initialize() override; - - /** Algorithm execute */ - TwoProngVertices operator()( LHCb::Tracks const&, DetectorElement const&, DeMagnet const& ) const override; - -private: - bool accept( const LHCb::TwoProngVertex& vertex ) const; - -private: - PublicToolHandle m_pocatool{this, "PocaTool", "TrajPoca"}; - PublicToolHandle m_vertexer{this, "TrackVertexer", "TrackVertexer"}; - double m_maxDistance = 5 * Gaudi::Units::mm; - Gaudi::Property m_minMass{this, "MinMass", 2.6 * Gaudi::Units::GeV}; - Gaudi::Property m_maxMass{this, "MaxMass", 4.0 * Gaudi::Units::GeV}; - double m_minJPsiMass = 3.065 * Gaudi::Units::GeV; - double m_maxJPsiMass = 3.125 * Gaudi::Units::GeV; - Gaudi::Property m_maxChi2TwoProngVertex{this, "MaxChi2TwoProngVertex", 25}; - double m_muonmass = 0.; - - mutable Gaudi::Accumulators::AveragingCounter<> m_num_combinations{this, "numcombinations"}; - mutable Gaudi::Accumulators::AveragingCounter<> m_num_postracks{this, "numpostracks"}; - mutable Gaudi::Accumulators::AveragingCounter<> m_num_negtracks{this, "numnegtracks"}; - mutable Gaudi::Accumulators::AveragingCounter<> m_num_selected{this, "numselected"}; -}; - -// Declaration of the Algorithm Factory -DECLARE_COMPONENT( TrackDiMuonMonitor ) - -StatusCode TrackDiMuonMonitor::initialize() { - return Transformer::initialize().andThen( [&] { - auto propertysvc = service( "LHCb::ParticlePropertySvc", true ); - if ( !propertysvc ) return StatusCode::FAILURE; - const auto* muon = propertysvc->find( "mu+" ); - if ( !muon ) return StatusCode::FAILURE; - m_muonmass = muon->mass(); - setHistoTopDir( "Track/" ); - return StatusCode::SUCCESS; - } ); -} - -TwoProngVertices TrackDiMuonMonitor::operator()( LHCb::Tracks const& muontracks, DetectorElement const& lhcb, - DeMagnet const& magneticfield ) const { - - // Sort them by charge, make some cuts - std::vector postracks, negtracks; - for ( auto track : muontracks ) { - // require tracks with T and (TT or Velo) - bool accepted = track->hasT() && ( track->hasVelo() || track->hasUT() ); - if ( accepted && !track->ancestors().empty() ) { - // muon track: take the first ancestor - track = track->ancestors().front(); - } else { - track = nullptr; - } - if ( track ) { ( track->firstState().qOverP() > 0 ? postracks : negtracks ).push_back( track ); } - } - - // turn them into trajectories - std::vector postrajs; - std::vector negtrajs; - postrajs.reserve( postracks.size() ); - negtrajs.reserve( negtracks.size() ); - for ( const auto& pos : postracks ) postrajs.emplace_back( *pos, &magneticfield ); - for ( const auto& neg : negtracks ) negtrajs.emplace_back( *neg, &magneticfield ); - - m_num_combinations += postracks.size() * negtracks.size(); - m_num_postracks += postracks.size(); - m_num_negtracks += negtracks.size(); - - // Create the output container - TwoProngVertices dimuoncontainer; - - for ( auto ipos = postracks.begin(); ipos != postracks.end(); ++ipos ) - for ( auto ineg = negtracks.begin(); ineg != negtracks.end(); ++ineg ) { - - if ( overlap( **ipos, **ineg ) ) continue; - const LHCb::TrackTraj& postraj = postrajs[std::distance( postracks.begin(), ipos )]; - const LHCb::TrackTraj& negtraj = negtrajs[std::distance( negtracks.begin(), ineg )]; - - double z_seed = postraj.beginRange() + 1 * Gaudi::Units::mm; - double mupos( z_seed ), muneg( z_seed ); - - // Calls pocatool - Gaudi::XYZVector deltaX; - StatusCode sc = m_pocatool->minimize( postraj, mupos, negtraj, muneg, deltaX, 0.001 * Gaudi::Units::mm ); - - if ( !sc.isSuccess() ) { - Warning( "TrajPoca Failure", StatusCode::SUCCESS, 0 ).ignore(); - continue; - } - double distance = deltaX.R(); - double z = 0.5 * ( mupos + muneg ); - if ( distance > m_maxDistance ) continue; - - // now make an invariant mass cut - Gaudi::XYZVector mompos = postraj.momentum( mupos ); - Gaudi::XYZVector momneg = negtraj.momentum( muneg ); - double mumass2 = m_muonmass * m_muonmass; - Gaudi::LorentzVector p4pos( mompos.X(), mompos.Y(), mompos.Z(), std::sqrt( mompos.Mag2() + mumass2 ) ); - Gaudi::LorentzVector p4neg( momneg.X(), momneg.y(), momneg.Z(), std::sqrt( momneg.Mag2() + mumass2 ) ); - double dimuonmass = ( p4pos + p4neg ).M(); - - if ( dimuonmass < m_minMass || m_maxMass < dimuonmass ) continue; - - // finally, create the vertex and cut on the chisquare - LHCb::State posstate = postraj.state( z ); - LHCb::State negstate = negtraj.state( z ); - auto vertex = m_vertexer->fit( posstate, negstate, *lhcb.geometry() ); - if ( !vertex ) continue; - vertex->addToTracks( *ipos ); - vertex->addToTracks( *ineg ); - if ( !accept( *vertex ) ) continue; - dimuoncontainer.add( vertex.release() ); - } - - m_num_selected += dimuoncontainer.size(); - - plot( dimuoncontainer.size(), "multiplicity", "J/psi candidate multiplicity", -0.5, 20.5, 21 ); - - return dimuoncontainer; -} - -bool TrackDiMuonMonitor::accept( const LHCb::TwoProngVertex& vertex ) const { - bool rc( false ); - static double chi2max = gsl_cdf_chisq_Qinv( 1e-6, 1 ); - double chi2 = vertex.chi2(); - double chi2prob = chi2 < chi2max ? gsl_cdf_chisq_Q( chi2, 1 ) : 0; - plot( chi2prob, "chi2prob", "chi2prob", 0, 1 ); - - if ( vertex.chi2() < m_maxChi2TwoProngVertex ) { - Gaudi::LorentzVector p4A = vertex.p4A( m_muonmass ); - Gaudi::LorentzVector p4B = vertex.p4B( m_muonmass ); - Gaudi::LorentzVector p4 = p4A + p4B; - double mass = p4.M(); - plot( mass, "mass", "dimuon mass", m_minMass, m_maxMass ); - - if ( m_minJPsiMass < mass && mass < m_maxJPsiMass ) { - - plot( vertex.position().x(), "vxJPsi", "J/psi candidate vertex x", -2, 2 ); - plot( vertex.position().y(), "vyJPsi", "J/psi candidate vertex y", -2, 2 ); - plot( vertex.position().z(), "vzJPsi", "J/psi candidate vertex z", -200, 200 ); - plot( chi2prob, "chi2probJPsi", "J/psi candidate chi2prob", 0, 1 ); - plot( mass, "massJPsi", "J/psi candidate mass", 3.05 * Gaudi::Units::GeV, 3.15 * Gaudi::Units::GeV ); - const LHCb::Track* postrack = vertex.trackA(); - const LHCb::Track* negtrack = vertex.trackB(); - if ( postrack->firstState().qOverP() < 0 ) std::swap( postrack, negtrack ); - double ppos = postrack->firstState().momentum().R(); - double pneg = negtrack->firstState().momentum().R(); - double pdif = ppos - pneg; - profile1D( pdif, mass, "massVersusMomDif", "dimuon mass versus p_{pos} - p_{neg}", -50 * Gaudi::Units::GeV, - 50 * Gaudi::Units::GeV, 20, "", 3.065 * Gaudi::Units::GeV, 3.125 * Gaudi::Units::GeV ); - plot( p4.P(), "momJPsi", "JPsi candidate momentum", 0, 100 * Gaudi::Units::GeV ); - plot( pdif, "momdifJPsi", "p_{pos} - p_{neg} for JPsi candidates", -50 * Gaudi::Units::GeV, - 50 * Gaudi::Units::GeV ); - plot( p4.Pt(), "ptJPsi", "JPsi candidate Pt", 0, 10 * Gaudi::Units::GeV ); - plot( p4.Eta(), "etaJPsi", "JPsi candidate eta", 2, 5 ); - - rc = true; - } - } - return rc; -} diff --git a/Tr/TrackMonitors/src/TrackMonitorBase.cpp b/Tr/TrackMonitors/src/TrackMonitorBase.cpp deleted file mode 100644 index 450c791d5c0627f53680a835dcf16bd229df5967..0000000000000000000000000000000000000000 --- a/Tr/TrackMonitors/src/TrackMonitorBase.cpp +++ /dev/null @@ -1,21 +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 files -#include "TrackMonitorBase.h" - -//============================================================================= -// Initialization. Check parameters -//============================================================================= -StatusCode TrackMonitorBase::initialize() { - return GaudiHistoAlg::initialize().andThen( [&] { - if ( histoTopDir().empty() ) setHistoTopDir( "Track/" ); - } ); -} diff --git a/Tr/TrackMonitors/src/TrackMonitorBase.h b/Tr/TrackMonitors/src/TrackMonitorBase.h deleted file mode 100644 index 8f5122eb264ae3f265923be67454a8ac71c0406b..0000000000000000000000000000000000000000 --- a/Tr/TrackMonitors/src/TrackMonitorBase.h +++ /dev/null @@ -1,45 +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 -#include "Event/Track.h" -#include "GaudiAlg/GaudiHistoAlg.h" -#include "Kernel/ITrajPoca.h" -#include "TrackInterfaces/ITrackExtrapolator.h" -#include -#include - -/** @class TrackMonitorBase TrackMonitorBase.h "TrackCheckers/TrackMonitorBase" - * - * Base class for track monitoring: essentially a 'box' of common tools - - * @author M. Needham. - * @date 6-5-2007 - */ - -class TrackMonitorBase : public GaudiHistoAlg { - -public: - /** Standard construtor */ - using GaudiHistoAlg::GaudiHistoAlg; - - /** Algorithm initialization */ - StatusCode initialize() override; - -protected: - /** Get a pointer to the track extrapolator - * @return extrapolator - */ - const ITrackExtrapolator* extrapolator() const { return m_extrapolator.get(); } - -private: - ToolHandle m_extrapolator{this, "Extrapolator", - "TrackMasterExtrapolator"}; ///< Pointer to extrapolator -}; diff --git a/Tr/TrackMonitors/src/TrackMuonMatchMonitor.cpp b/Tr/TrackMonitors/src/TrackMuonMatchMonitor.cpp index d21b60f946812714232b7ae97a626acf64750743..a138c7449105a3fea22e4f012527f7b50f5dba77 100644 --- a/Tr/TrackMonitors/src/TrackMuonMatchMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackMuonMatchMonitor.cpp @@ -20,12 +20,10 @@ #include "Event/State.h" #include "Event/Track.h" -#include "GaudiAlg/GaudiHistoAlg.h" +#include "Gaudi/Accumulators/StaticHistogram.h" #include "GaudiKernel/SystemOfUnits.h" #include "LHCbAlgs/Consumer.h" -#include "AIDA/IHistogram1D.h" - namespace { struct Cache { double zM1, MAXsizeX, MAXsizeY; @@ -43,9 +41,9 @@ namespace LHCb { * @date 2010-01-22 */ class TrackMuonMatchMonitor - : public LHCb::Algorithm::Consumer< - void( Track::Range const&, MuonCoords const&, DetectorElement const&, DeMuonDetector const&, Cache const& ), - Algorithm::Traits::usesBaseAndConditions> { + : public LHCb::Algorithm::Consumer> { public: TrackMuonMatchMonitor( const std::string& name, ISvcLocator* pSvcLoc ) : Consumer{name, @@ -65,9 +63,27 @@ namespace LHCb { Gaudi::Property m_maxErrX{this, "MaxErrX", 5 * Gaudi::Units::mm}; Gaudi::Property m_maxErrY{this, "MaxErrY", 20 * Gaudi::Units::mm}; - static constexpr int nREGIONS = 4; - AIDA::IHistogram1D * m_resx_a[nREGIONS], *m_resy_a[nREGIONS], *m_resx_c[nREGIONS], *m_resy_c[nREGIONS]; - double m_hisxmax[nREGIONS]; + mutable std::array, 4> m_resx_a{ + {{this, "resX_ASide_M1R1", "resX_ASide_M1R1", {100, -80, 80}}, + {this, "resX_ASide_M1R2", "resX_ASide_M1R2", {100, -120, 120}}, + {this, "resX_ASide_M1R3", "resX_ASide_M1R3", {100, -160, 160}}, + {this, "resX_ASide_M1R4", "resX_ASide_M1R4", {100, -200, 200}}}}; + mutable std::array, 4> m_resy_a{ + {{this, "resY_ASide_M1R1", "resY_ASide_M1R1", {100, -80, 80}}, + {this, "resY_ASide_M1R2", "resY_ASide_M1R2", {100, -120, 120}}, + {this, "resY_ASide_M1R3", "resY_ASide_M1R3", {100, -160, 160}}, + {this, "resY_ASide_M1R4", "resY_ASide_M1R4", {100, -200, 200}}}}; + mutable std::array, 4> m_resx_c{ + {{this, "resX_CSide_M1R1", "resX_CSide_M1R1", {100, -80, 80}}, + {this, "resX_CSide_M1R2", "resX_CSide_M1R2", {100, -120, 120}}, + {this, "resX_CSide_M1R3", "resX_CSide_M1R3", {100, -160, 160}}, + {this, "resX_CSide_M1R4", "resX_CSide_M1R4", {100, -200, 200}}}}; + mutable std::array, 4> m_resy_c{ + {{this, "resY_CSide_M1R1", "resY_CSide_M1R1", {100, -80, 80}}, + {this, "resY_CSide_M1R2", "resY_CSide_M1R2", {100, -120, 120}}, + {this, "resY_CSide_M1R3", "resY_CSide_M1R3", {100, -160, 160}}, + {this, "resY_CSide_M1R4", "resY_CSide_M1R4", {100, -200, 200}}}}; + static constexpr double m_hisxmax[4] = {80, 120, 160, 200}; }; // Declaration of the Algorithm Factory @@ -82,24 +98,6 @@ StatusCode LHCb::TrackMuonMatchMonitor::initialize() { return {muonDet.getStationZ( m_iMS ), muonDet.getOuterX( m_iMS ), muonDet.getOuterY( m_iMS )}; } ); std::string name; - setHistoTopDir( "Track/" ); - for ( int iR = 0; iR < nREGIONS; ++iR ) { - // in the "signal" region +/- 6 error units respect the track extrapolation point - unsigned int nbin = 100; - double max = 80. + 40. * float( iR ); - double min = -max; - name = "resX_ASide_M1R" + std::to_string( iR + 1 ); - m_resx_a[iR] = book1D( name, name, min, max, nbin ); - name = "resY_ASide_M1R" + std::to_string( iR + 1 ); - m_resy_a[iR] = book1D( name, name, min, max, nbin ); - name = "resX_CSide_M1R" + std::to_string( iR + 1 ); - m_resx_c[iR] = book1D( name, name, min, max, nbin ); - name = "resY_CSide_M1R" + std::to_string( iR + 1 ); - m_resy_c[iR] = book1D( name, name, min, max, nbin ); - m_hisxmax[iR] = max; - } - - return StatusCode::SUCCESS; } ); } @@ -152,14 +150,13 @@ void LHCb::TrackMuonMatchMonitor::operator()( LHCb::Track::Range const& tTracks, double deltaY = hit.y - ( stateAtM1.y() + stateAtM1.ty() * deltaZ ); if ( std::abs( deltaX ) < m_hisxmax[region] && std::abs( deltaY ) < m_hisxmax[region] ) { - - AIDA::IHistogram1D *tempx, *tempy; - - tempx = hit.x > 0 ? m_resx_a[region] : m_resx_c[region]; - tempy = hit.x > 0 ? m_resy_a[region] : m_resy_c[region]; - - tempx->fill( deltaX ); // X residuals on the same Z as the hit - tempy->fill( deltaY ); // Y residuals on the same Z as the hit + if ( hit.x > 0 ) { + ++m_resx_a[region][deltaX]; + ++m_resy_a[region][deltaY]; + } else { + ++m_resx_c[region][deltaX]; + ++m_resy_c[region][deltaY]; + } } } } diff --git a/Tr/TrackMonitors/src/TrackV0Monitor.cpp b/Tr/TrackMonitors/src/TrackV0Monitor.cpp deleted file mode 100644 index 32a677cb84e1a344bc21ea5be24c13c175f4b074..0000000000000000000000000000000000000000 --- a/Tr/TrackMonitors/src/TrackV0Monitor.cpp +++ /dev/null @@ -1,143 +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 "Event/TwoProngVertex.h" -#include "GaudiAlg/GaudiHistoAlg.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "GaudiKernel/ToolHandle.h" -#include "Kernel/IParticlePropertySvc.h" -#include "Kernel/ParticleProperty.h" -#include "LHCbAlgs/Consumer.h" -#include "TrackInterfaces/ITrackVertexer.h" -#include "gsl/gsl_cdf.h" - -class TrackV0Monitor - : public LHCb::Algorithm::Consumer> { -public: - /** Standard construtor */ - TrackV0Monitor( const std::string& name, ISvcLocator* pSvcLocator ); - - /** Algorithm initialize */ - StatusCode initialize() override; - - /** Algorithm execute */ - void operator()( const LHCb::TwoProngVertices&, const LHCb::RecVertex::Range& ) const override; - -private: - void process( const LHCb::TwoProngVertex& vertex, const LHCb::RecVertex::Range& pvs, const std::string& dir ) const; - void computeDecayLength( const LHCb::TwoProngVertex& vertex, const LHCb::RecVertex::Range& pvs, double& ipchi2, - double& decaylength, double& decaylengtherr ) const; - -private: - ToolHandle m_vertexer{this, "TrackVertexer", "TrackVertexer"}; - Gaudi::Property m_maxIPChi2{this, "MaxIPChi2", 25.0}; - double m_pionmass{0}; - double m_protonmass{0}; -}; - -// Declaration of the Algorithm Factory -DECLARE_COMPONENT( TrackV0Monitor ) - -//============================================================================= -// Standard constructor, initializes variables -//============================================================================= -TrackV0Monitor::TrackV0Monitor( const std::string& name, ISvcLocator* pSvcLocator ) - : Consumer{name, - pSvcLocator, - {KeyValue{"V0Location", LHCb::TwoProngVertexLocation::Default}, - KeyValue{"PrimaryLocation", LHCb::RecVertexLocation::Primary}}} {} - -StatusCode TrackV0Monitor::initialize() { - setHistoTopDir( "Track/" ); - StatusCode sc = Consumer::initialize(); - LHCb::IParticlePropertySvc* propertysvc = svc( "LHCb::ParticlePropertySvc", true ); - if ( !propertysvc ) return StatusCode::FAILURE; - - const LHCb::ParticleProperty* pion = propertysvc->find( "pi+" ); - if ( !pion ) return StatusCode::FAILURE; - m_pionmass = pion->mass(); - const LHCb::ParticleProperty* proton = propertysvc->find( "p+" ); - if ( !proton ) return StatusCode::FAILURE; - m_protonmass = proton->mass(); - // const LHCb::ParticleProperty* kshort = propertysvc->find( "KS0" ) ; - // m_k0pid = kshort ? LHCb::ParticleID(kshort->pdgID()) : 0 ; - // const LHCb::ParticleProperty* lambda = propertysvc->find( "Lambda0" ); - // m_lambdapid = lambda ? LHCb::ParticleID(lambda->pdgID()) : 0 ; - // m_lambdabarpid = lambda ? LHCb::ParticleID(lambda->antiParticle()->pdgID()) : 0 ; - - return sc; -} - -void TrackV0Monitor::process( const LHCb::TwoProngVertex& vertex, const LHCb::RecVertex::Range& pvs, - const std::string& dir ) const { - static double chi2max = gsl_cdf_chisq_Qinv( 1e-6, 1 ); - double chi2 = vertex.chi2(); - double chi2prob = chi2 < chi2max ? gsl_cdf_chisq_Q( chi2, 1 ) : 0; - double pipimass = vertex.mass( m_pionmass, m_pionmass ); - double ppimass = vertex.mass( m_protonmass, m_pionmass ); - double pipmass = vertex.mass( m_pionmass, m_protonmass ); - if ( vertex.trackA()->firstState().qOverP() < 0 ) std::swap( pipmass, ppimass ); - - // compute the ipchi2 and the decaylength - double ipchi2( -1 ), decaylength( 0 ), decaylengtherr( 1 ); - computeDecayLength( vertex, pvs, ipchi2, decaylength, decaylengtherr ); - plot( ipchi2, dir + "/ipchi2", "ip chi2", 0, 100 ); - if ( ipchi2 < m_maxIPChi2 ) { - plot( chi2prob, dir + "/chi2prob", "chi2prob", 0, 1 ); - plot( vertex.position().z(), dir + "/z", "z", -100, 2000 ); - plot( pipimass, dir + "/pipimass", "pi+ pi- mass", 450, 550 ); - plot( ppimass, dir + "/ppimass", "p+ pi- mass", 1090, 1140 ); - plot( pipmass, dir + "/pipmass", "p- pi+ mass", 1090, 1140 ); - plot( decaylength / decaylengtherr, dir + "/dls", "decay length significance", -10, 100 ); - } - - if ( fullDetail() ) { - plot( pipimass, dir + "/pipimassall", "pi+ pi- mass", 450, 550 ); - plot( ppimass, dir + "/ppimassall", "p+ pi- mass", 1090, 1140 ); - plot( pipmass, dir + "/pipmassall", "p- pi+ mass", 1090, 1140 ); - plot( decaylength / decaylengtherr, dir + "/dlsall", "decay length significance", -10, 100 ); - plot( vertex.position().x(), dir + "/x", "x", -100, 100 ); - plot( vertex.position().y(), dir + "/y", "y", -100, 100 ); - } -} - -void TrackV0Monitor::operator()( const LHCb::TwoProngVertices& v0container, const LHCb::RecVertex::Range& pvs ) const { - plot( v0container.size(), "v0multiplicity", "v0multiplicity", -0.5, 20.5, 21 ); - for ( const LHCb::TwoProngVertex* vertex : v0container ) { - // determine the type: Long-Long, DS-DS, or DS-Long - const LHCb::Track::Types trackA = vertex->trackA()->type(); - const LHCb::Track::Types trackB = vertex->trackB()->type(); - auto type = trackA == LHCb::Track::Types::Long && trackB == LHCb::Track::Types::Long - ? "LongLong" - : trackA == LHCb::Track::Types::Downstream && trackB == LHCb::Track::Types::Downstream - ? "DownstreamDownstream" - : "Mixed"; - process( *vertex, pvs, type ); - } -} - -void TrackV0Monitor::computeDecayLength( const LHCb::TwoProngVertex& vertex, const LHCb::RecVertex::Range& pvs, - double& ipchi2, double& decaylength, double& decaylengtherr ) const { - // some default values - ipchi2 = -1; - decaylength = 0; - decaylengtherr = 1; - - for ( const auto& pv : pvs ) { - double thisipchi2( -1 ), thisdecaylength( 0 ), thisdecaylengtherr( 1 ); - bool success = m_vertexer->computeDecayLength( vertex, *pv, thisipchi2, thisdecaylength, thisdecaylengtherr ); - if ( success && ( ipchi2 < 0 || thisipchi2 < ipchi2 ) ) { - ipchi2 = thisipchi2; - decaylength = thisdecaylength; - decaylengtherr = thisdecaylengtherr; - } - } -} diff --git a/Tr/TrackMonitors/src/UTTrackMonitor.cpp b/Tr/TrackMonitors/src/UTTrackMonitor.cpp index afe9d3655a18cc24b661fc700f4f3cc913c308d8..529acfc1ae9cce43bfd95db718efd5d6bfbc7dd0 100644 --- a/Tr/TrackMonitors/src/UTTrackMonitor.cpp +++ b/Tr/TrackMonitors/src/UTTrackMonitor.cpp @@ -28,8 +28,6 @@ #include "Kernel/UTNames.h" #include "UTDAQ/UTInfo.h" -#include "TrackMonitorBase.h" - #include "LHCbAlgs/Consumer.h" using namespace LHCb; @@ -42,7 +40,7 @@ namespace { namespace LHCb { class UTTrackMonitor : public LHCb::Algorithm::Consumer> { + LHCb::DetDesc::usesConditions> { public: UTTrackMonitor( const std::string& name, ISvcLocator* pSvcLocator ); diff --git a/Tr/TrackUtils/src/TrackBestTrackCreator.cpp b/Tr/TrackUtils/src/TrackBestTrackCreator.cpp index be6d1186311df530d97a8bc4eaebbab4ae38cc5c..96bee04a1ea1fa8d5164bec8ae551ac60ca7a996 100644 --- a/Tr/TrackUtils/src/TrackBestTrackCreator.cpp +++ b/Tr/TrackUtils/src/TrackBestTrackCreator.cpp @@ -38,12 +38,7 @@ #include #include -#if DEBUGHISTOGRAMS -# include "GaudiAlg/GaudiHistoAlg.h" -using TrackBestTrackCreatorBase = LHCb::DetDesc::ConditionAccessorHolder>>; -#else using TrackBestTrackCreatorBase = LHCb::DetDesc::ConditionAccessorHolder>; -#endif using namespace LHCb::Event::Enum::Track; @@ -429,25 +424,16 @@ void LHCb::TrackBestTrackCreator::fitAndUpdateCounters( TrackData& td, IGeometry bool LHCb::TrackBestTrackCreator::veloOrClones( const TrackData& lhs, const TrackData& rhs ) const { const double f = lhs.overlapFraction( rhs, TrackData::VP ); -#ifdef DEBUGHISTOGRAMS - if ( f > 0 ) plot1D( f, "veloOverlapFractionH1", 0, 1 ); -#endif return ( f > m_maxOverlapFracVelo ); } bool LHCb::TrackBestTrackCreator::TClones( const TrackData& lhs, const TrackData& rhs ) const { const double f = lhs.overlapFraction( rhs, TrackData::T ); -#ifdef DEBUGHISTOGRAMS - if ( f > 0 ) plot1D( f, "TOverlapFractionH1", 0, 1 ); -#endif return f > m_maxOverlapFracT; } bool LHCb::TrackBestTrackCreator::UTClones( const TrackData& lhs, const TrackData& rhs ) const { const double f = lhs.overlapFraction( rhs, TrackData::UT ); -#ifdef DEBUGHISTOGRAMS - if ( f > 0 ) plot1D( f, "UTOverlapFractionH1", 0, 1 ); -#endif return f > m_maxOverlapFracUT; } @@ -460,30 +446,12 @@ bool LHCb::TrackBestTrackCreator::areClones( const TrackData& it, const TrackDat const double dqop = it.qOverP() - jt.qOverP(); switch ( to_index( itrack.type(), jtrack.type() ) ) { case to_index( Track::Types::Long, Track::Types::Long ): -#ifdef DEBUGHISTOGRAMS - if ( TClones( it, jt ) && veloOrClones( it, jt ) ) { - plot( dqop, "LLDqopClones", -1e-5, 1e-5 ); - } else if ( TClones( it, jt ) ) { - plot( dqop, "LLDqopTClones", -1e-5, 1e-5 ); - } else if ( veloOrClones( it, jt ) ) { - plot( dqop, "LLDqopVeloOrClones", -1e-5, 1e-5 ); - } -#endif return TClones( it, jt ) && ( std::abs( dqop ) < m_minLongLongDeltaQoP || veloOrClones( it, jt ) ); case to_index( Track::Types::Long, Track::Types::Downstream ): case to_index( Track::Types::Downstream, Track::Types::Long ): -#ifdef DEBUGHISTOGRAMS - if ( TClones( it, jt ) ) { - plot( dqop, "DLDqop", -2e-5, 2e-5 ); - if ( UTClones( it, jt ) ) plot( dqop, "DLDqopUTClones", -2e-5, 2e-5 ); - } -#endif return TClones( it, jt ) && ( std::abs( dqop ) < m_minLongDownstreamDeltaQoP || UTClones( it, jt ) ); case to_index( Track::Types::Downstream, Track::Types::Downstream ): // it seems that there are no down stream tracks that share T hits ... -#ifdef DEBUGHISTOGRAMS - if ( TClones( it, jt ) ) { plot( dqop, "DDDqop", -1e-4, 1e-4 ); } -#endif return TClones( it, jt ) && UTClones( it, jt ); case to_index( Track::Types::Long, Track::Types::Upstream ): case to_index( Track::Types::Upstream, Track::Types::Long ): diff --git a/Tr/TrackUtils/src/TrackCloneKiller.cpp b/Tr/TrackUtils/src/TrackCloneKiller.cpp index 30056eb632266fcecac64741a7655b3a8594151d..8a82f0fd722a4bd1084968399a4f8d8ededa3cb8 100644 --- a/Tr/TrackUtils/src/TrackCloneKiller.cpp +++ b/Tr/TrackUtils/src/TrackCloneKiller.cpp @@ -29,12 +29,6 @@ #include #include #include -#if DEBUGHISTOGRAMS -# include "GaudiAlg/GaudiHistoAlg.h" -using TrackCloneKillerBase = GaudiHistoAlg; -#else -using TrackCloneKillerBase = GaudiAlgorithm; -#endif namespace { @@ -71,20 +65,8 @@ namespace { * - initial release, largely copied from TrackBestTrackCreator */ class TrackCloneKiller final - : public LHCb::Algorithm::Transformer> { + : public LHCb::Algorithm::Transformer { public: - /// Standard constructor - - using base_class_t = LHCb::Algorithm::Transformer>; - // using base_class_t::addConditionDerivation; - using base_class_t::debug; - using base_class_t::error; - using base_class_t::info; - using base_class_t::inputLocation; - using base_class_t::msgLevel; - TrackCloneKiller( const std::string& name, ISvcLocator* pSvcLocator ); virtual StatusCode initialize() override; ///< Algorithm initialization @@ -231,26 +213,16 @@ LHCb::Tracks TrackCloneKiller::operator()( const LHCb::Tracks& inTracks, const L bool TrackCloneKiller::veloOrClones( const TrackData& lhs, const TrackData& rhs ) const { const double f = lhs.overlapFraction( rhs, TrackData::VP ); -#ifdef DEBUGHISTOGRAMS - if ( f > 0 ) plot1D( fR, "veloOverlapFractionH1", 0, 1 ); -#else return f > m_maxOverlapFracVelo; -#endif } bool TrackCloneKiller::TClones( const TrackData& lhs, const TrackData& rhs ) const { const double f = lhs.overlapFraction( rhs, TrackData::T ); -#ifdef DEBUGHISTOGRAMS - if ( f > 0 ) plot1D( f, "TOverlapFractionH1", 0, 1 ); -#endif return f > m_maxOverlapFracT; } bool TrackCloneKiller::UTClones( const TrackData& lhs, const TrackData& rhs ) const { const double f = lhs.overlapFraction( rhs, TrackData::UT ); -#ifdef DEBUGHISTOGRAMS - if ( f > 0 ) plot1D( f, "UTOverlapFractionH1", 0, 1 ); -#endif return f > m_maxOverlapFracUT; } @@ -264,30 +236,12 @@ bool TrackCloneKiller::areClones( const TrackData& it, const TrackData& jt ) con const double dqop = it.qOverP() - jt.qOverP(); switch ( to_index( itrack.type(), jtrack.type() ) ) { case to_index( LHCb::Track::Types::Long, LHCb::Track::Types::Long ): -#ifdef DEBUGHISTOGRAMS - if ( TClones( it, jt ) && veloOrClones( it, jt ) ) { - plot( dqop, "LLDqopClones", -1e-5, 1e-5 ); - } else if ( TClones( it, jt ) ) { - plot( dqop, "LLDqopTClones", -1e-5, 1e-5 ); - } else if ( veloOrClones( it, jt ) ) { - plot( dqop, "LLDqopVeloOrClones", -1e-5, 1e-5 ); - } -#endif return TClones( it, jt ) && ( std::abs( dqop ) < m_minLongLongDeltaQoP || veloOrClones( it, jt ) ); case to_index( LHCb::Track::Types::Long, LHCb::Track::Types::Downstream ): case to_index( LHCb::Track::Types::Downstream, LHCb::Track::Types::Long ): -#ifdef DEBUGHISTOGRAMS - if ( TClones( it, jt ) ) { - plot( dqop, "DLDqop", -2e-5, 2e-5 ); - if ( UTClones( it, jt ) ) plot( dqop, "DLDqopUTClones", -2e-5, 2e-5 ); - } -#endif return TClones( it, jt ) && ( std::abs( dqop ) < m_minLongDownstreamDeltaQoP || UTClones( it, jt ) ); case to_index( LHCb::Track::Types::Downstream, LHCb::Track::Types::Downstream ): // it seems that there are no down stream tracks that share T hits ... -#ifdef DEBUGHISTOGRAMS - if ( TClones( it, jt ) ) { plot( dqop, "DDDqop", -1e-4, 1e-4 ); } -#endif return TClones( it, jt ) && UTClones( it, jt ); case to_index( LHCb::Track::Types::Long, LHCb::Track::Types::Upstream ): case to_index( LHCb::Track::Types::Upstream, LHCb::Track::Types::Long ):