diff --git a/Rich/RichFutureRecMonitors/python/RichFutureRecMonitors/ConfiguredRecoMonitors.py b/Rich/RichFutureRecMonitors/python/RichFutureRecMonitors/ConfiguredRecoMonitors.py
index fffb1f773a0218fec782adf670860c6536852b42..ae55a7be4a502abfd310188018fb9dd1eeee6271 100644
--- a/Rich/RichFutureRecMonitors/python/RichFutureRecMonitors/ConfiguredRecoMonitors.py
+++ b/Rich/RichFutureRecMonitors/python/RichFutureRecMonitors/ConfiguredRecoMonitors.py
@@ -266,6 +266,7 @@ def RichRecMonitors(
ckAngles.CherenkovAnglesLocation = locs["SignalCherenkovAnglesLocation"]
ckAngles.SummaryTracksLocation = locs["SummaryTracksLocation"]
ckAngles.PhotonToParentsLocation = locs["PhotonToParentsLocation"]
+ ckAngles.PhotonMirrorDataLocation = locs["PhotonMirrorDataLocation"]
# Options
if onlineBrunelMode:
# Open up the CK res plot range, for the Wide photon selection
@@ -285,6 +286,7 @@ def RichRecMonitors(
ckAnglesT.CherenkovAnglesLocation = locs["SignalCherenkovAnglesLocation"]
ckAnglesT.SummaryTracksLocation = locs["SummaryTracksLocation"]
ckAnglesT.PhotonToParentsLocation = locs["PhotonToParentsLocation"]
+ ckAnglesT.PhotonMirrorDataLocation = locs["PhotonMirrorDataLocation"]
# Options
if onlineBrunelMode:
# Clone settings from the nominal monitor
diff --git a/Rich/RichFutureRecMonitors/scripts/ThetaPhiBias.py b/Rich/RichFutureRecMonitors/scripts/ThetaPhiBias.py
new file mode 100755
index 0000000000000000000000000000000000000000..a4940cd40ee4274d0be0ac7e61d371cbc3957344
--- /dev/null
+++ b/Rich/RichFutureRecMonitors/scripts/ThetaPhiBias.py
@@ -0,0 +1,174 @@
+#!/usr/bin/env python
+
+###############################################################################
+# (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. #
+###############################################################################
+
+from array import array
+
+import numpy as np
+from GaudiPython import gbl
+from ROOT import *
+from ROOT import gROOT
+
+gROOT.SetBatch(True)
+
+hfile = TFile.Open("./RichRefIndexCalib.root")
+
+fitter = gbl.Rich.Rec.CKResolutionFitter()
+
+lhcbFont = 62
+lhcbStatFontSize = 0.03
+lhcbStatBoxWidth = 0.15
+lhcbMarkerType = 20
+lhcbMarkerSize = 0.8
+gStyle.SetPadTickX(1)
+gStyle.SetPadTickY(1)
+gStyle.SetOptFit(101)
+gStyle.SetOptStat(1110)
+gStyle.SetStatBorderSize(1)
+gStyle.SetStatFont(lhcbFont)
+gStyle.SetStatFontSize(lhcbStatFontSize)
+gStyle.SetStatW(lhcbStatBoxWidth)
+gStyle.SetMarkerStyle(lhcbMarkerType)
+gStyle.SetMarkerSize(lhcbMarkerSize)
+nPhiBins = 16
+nQuadrants = 4
+phiInc = 2.0 * 3.1415 / nPhiBins
+
+
+def quadString(q):
+ if q == 0:
+ return "Quadrant x>0 y>0"
+ if q == 1:
+ return "Quadrant x<0 y>0"
+ if q == 2:
+ return "Quadrant x>0 y<0"
+ if q == 3:
+ return "Quadrant x<0 y<0"
+ return "UNDEFINED"
+
+
+def padToQuad(p):
+ if p == 0:
+ return 1
+ if p == 1:
+ return 0
+ if p == 2:
+ return 3
+ if p == 3:
+ return 2
+
+
+for rich in ["Rich1", "Rich2"]:
+ canvas = TCanvas("CKFit" + rich, "CKFit" + rich, 1200, 900)
+
+ # Open PDF
+ pdfName = rich + "-CKThetaPhiPlots.pdf"
+ canvas.Print(pdfName + "[")
+
+ corrs = []
+
+ phi_bias = {}
+ phi_bias_err = {}
+ phi_reso = {}
+ phi_reso_err = {}
+ phi_bins = {}
+ phi_bins_err = {}
+
+ for q in range(0, nQuadrants):
+ phi_bias[q] = array("d")
+ phi_bias_err[q] = array("d")
+ phi_reso[q] = array("d")
+ phi_reso_err[q] = array("d")
+ phi_bins[q] = array("d")
+ phi_bins_err[q] = array("d")
+ for phi in range(0, nPhiBins):
+ hpath = (
+ "RICH/RiCKResLong/"
+ + rich
+ + "Gas/Quadrants/Q"
+ + str(q)
+ + "/PhiBins/Bin"
+ + str(phi)
+ + "/ckResAll"
+ )
+ print(hpath)
+
+ h = hfile.Get(hpath)
+
+ res = fitter.fit(
+ h, gbl.Rich.Rich1Gas if rich == "Rich1" else gbl.Rich.Rich2Gas
+ )
+
+ h.Draw("e")
+ res.overallFitFunc.SetLineColor(kBlue + 2)
+ res.bkgFitFunc.SetLineColor(kRed + 3)
+ res.overallFitFunc.Draw("SAME")
+ res.bkgFitFunc.Draw("SAME")
+
+ canvas.Print(pdfName)
+
+ bias = 1000 * res.ckShift if res.fitOK else 0.0
+ bias = float("{:.8f}".format(bias))
+ phi_bias[q].append(bias)
+ phi_bias_err[q].append(1000 * res.ckShiftErr if res.fitOK else 0.0)
+
+ reso = 1000 * res.ckResolution if res.fitOK else 0.0
+ reso = float("{:.8f}".format(reso))
+ phi_reso[q].append(reso)
+ phi_reso_err[q].append(1000 * res.ckResolutionErr if res.fitOK else 0.0)
+
+ phi_bins[q].append((0.5 + phi) * phiInc)
+ phi_bins_err[q].append(0)
+
+ # corrs.append(phi_bias)
+
+ # print(corrs)
+
+ canvas.Clear()
+ canvas.cd()
+ canvas.Divide(2, 2)
+ graph1 = []
+ for p in range(0, nQuadrants):
+ canvas.cd(p + 1)
+ q = padToQuad(p)
+ graph1.append(
+ TGraphErrors(
+ nPhiBins, phi_bins[q], phi_bias[q], phi_bins_err[q], phi_bias_err[q]
+ )
+ )
+ graph1[p].SetTitle(rich + " CK Theta Bias V CK Phi | " + quadString(q))
+ graph1[p].GetXaxis().SetTitle("CK Phi / rad")
+ graph1[p].GetYaxis().SetTitle("CK Theta bias / mrad")
+ graph1[p].Draw("ALP")
+ gPad.Update()
+ canvas.Print(pdfName)
+
+ canvas.Clear()
+ canvas.cd()
+ canvas.Divide(2, 2)
+ graph2 = []
+ for p in range(0, nQuadrants):
+ canvas.cd(p + 1)
+ q = padToQuad(p)
+ graph2.append(
+ TGraphErrors(
+ nPhiBins, phi_bins[q], phi_reso[q], phi_bins_err[q], phi_reso_err[q]
+ )
+ )
+ graph2[p].SetTitle(rich + " CK Theta Resolution V CK Phi | " + quadString(q))
+ graph2[p].GetXaxis().SetTitle("CK Phi / rad")
+ graph2[p].GetYaxis().SetTitle("CK Theta Resolution / mrad")
+ graph2[p].Draw("ALP")
+ canvas.Print(pdfName)
+
+ # Close PDF
+ canvas.Print(pdfName + "]")
diff --git a/Rich/RichFutureRecMonitors/src/RichSIMDPhotonCherenkovAngles.cpp b/Rich/RichFutureRecMonitors/src/RichSIMDPhotonCherenkovAngles.cpp
index 3c9ef62c3da2a5fb7e4b092f9e9fa505cdf96a5b..5e1fa05e7f9d79259b5813f427181714d553e93c 100644
--- a/Rich/RichFutureRecMonitors/src/RichSIMDPhotonCherenkovAngles.cpp
+++ b/Rich/RichFutureRecMonitors/src/RichSIMDPhotonCherenkovAngles.cpp
@@ -29,6 +29,7 @@
#include "RichFutureRecEvent/RichSummaryEventData.h"
// Rich Utils
+#include "RichFutureUtils/RichSIMDMirrorData.h"
#include "RichUtils/RichPDIdentifier.h"
#include "RichUtils/RichTrackSegment.h"
#include "RichUtils/ZipRange.h"
@@ -39,12 +40,42 @@
// STD
#include <algorithm>
#include <array>
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <memory>
#include <mutex>
#include <sstream>
#include <vector>
namespace Rich::Future::Rec::Moni {
+ namespace {
+ /// Inner/Outer regions
+ enum InOutRegion : std::uint8_t { Inner = 0, Outer };
+ inline std::string regionString( const InOutRegion reg ) noexcept { return ( Inner == reg ? "Inner" : "Outer" ); }
+
+ /// arrays for histograms
+ template <typename TYPE>
+ using ColumnHists = Hist::Array<TYPE, LHCb::RichSmartID::MaPMT::MaxModuleColumnsAnyPanel>;
+ template <typename TYPE>
+ using InOutArray = Hist::Array<TYPE, 2, InOutRegion>;
+
+ /// For Phi binning
+ static constexpr std::size_t NPhiBins = 16;
+ static constexpr double phiInc = ( 2.0 * M_PI / NPhiBins );
+ template <typename T>
+ inline auto phiBin( const T& phi ) {
+ const std::size_t bin = std::floor( phi / phiInc );
+ return std::min( NPhiBins - 1, bin );
+ }
+ inline auto binMinPhi( const std::size_t bin ) noexcept { return bin * phiInc; }
+ inline auto binMaxPhi( const std::size_t bin ) noexcept { return ( bin + 1 ) * phiInc; }
+ inline auto binAvgPhi( const std::size_t bin ) noexcept { return ( bin + 0.5 ) * phiInc; }
+ template <typename TYPE>
+ using PhiHists = Hist::Array<TYPE, NPhiBins>;
+ } // namespace
+
/** @class SIMDPhotonCherenkovAngles RichSIMDPhotonCherenkovAngles.h
*
* Monitors the reconstructed cherenkov angles.
@@ -60,7 +91,8 @@ namespace Rich::Future::Rec::Moni {
const Relations::PhotonToParents::Vector&, //
const LHCb::RichTrackSegment::Vector&, //
const CherenkovAngles::Vector&, //
- const SIMDCherenkovPhoton::Vector& ),
+ const SIMDCherenkovPhoton::Vector&, //
+ const SIMDMirrorData::Vector& ),
Gaudi::Functional::Traits::BaseClass_t<HistoAlgBase>> {
public:
@@ -73,7 +105,8 @@ namespace Rich::Future::Rec::Moni {
KeyValue{ "PhotonToParentsLocation", Relations::PhotonToParentsLocation::Default },
KeyValue{ "TrackSegmentsLocation", LHCb::RichTrackSegmentLocation::Default },
KeyValue{ "CherenkovAnglesLocation", CherenkovAnglesLocation::Signal },
- KeyValue{ "CherenkovPhotonLocation", SIMDCherenkovPhotonLocation::Default } } ) {
+ KeyValue{ "CherenkovPhotonLocation", SIMDCherenkovPhotonLocation::Default },
+ KeyValue{ "PhotonMirrorDataLocation", SIMDMirrorDataLocation::Default } } ) {
// setProperty( "OutputLevel", MSG::VERBOSE ).ignore();
// print some stats on the final plots
setProperty( "HistoPrint", true ).ignore();
@@ -88,41 +121,57 @@ namespace Rich::Future::Rec::Moni {
const Relations::PhotonToParents::Vector& photToSegPix, //
const LHCb::RichTrackSegment::Vector& segments, //
const CherenkovAngles::Vector& expTkCKThetas, //
- const SIMDCherenkovPhoton::Vector& photons ) const override;
+ const SIMDCherenkovPhoton::Vector& photons, //
+ const SIMDMirrorData::Vector& mirrorData ) const override;
protected:
/// Pre-Book all histograms
StatusCode prebookHistograms() override;
- private:
- /// Get the per PD resolution histogram ID
- inline std::string pdResPlotID( const LHCb::RichSmartID hpd ) const {
- const Rich::DAQ::PDIdentifier hid( hpd );
- std::ostringstream id;
- id << "PDs/pd-" << hid.number();
- return id.str();
- }
-
private:
// cached data
// CK theta histograms
- mutable Hist::DetArray<Hist::H1D<>> h_thetaRec;
- mutable Hist::DetArray<Hist::H1D<>> h_phiRec;
- mutable Hist::DetArray<Hist::H1D<>> h_ckResAll;
- mutable Hist::DetArray<Hist::PanelArray<Hist::H1D<>>> h_ckResAllPerPanel;
- mutable Hist::DetArray<Hist::H1D<>> h_ckResAllBetaCut;
- mutable Hist::DetArray<Hist::PanelArray<Hist::H1D<>>> h_ckResAllPerPanelBetaCut;
- mutable Hist::DetArray<Hist::H1D<>> h_ckResAllInner;
- mutable Hist::DetArray<Hist::H1D<>> h_ckResAllOuter;
- mutable Hist::DetArray<Hist::PanelArray<Hist::H1D<>>> h_ckResAllPerPanelInner;
- mutable Hist::DetArray<Hist::PanelArray<Hist::H1D<>>> h_ckResAllPerPanelOuter;
-
- using ColumnHists =
- Hist::DetArray<Hist::PanelArray<Hist::Array<Hist::H1D<>, LHCb::RichSmartID::MaPMT::MaxModuleColumnsAnyPanel>>>;
- mutable ColumnHists h_ckResAllPerCol;
- mutable ColumnHists h_ckResAllPerColInner;
- mutable ColumnHists h_ckResAllPerColOuter;
+ mutable Hist::DetArray<Hist::H1D<>> h_thetaRec;
+ mutable Hist::DetArray<Hist::H1D<>> h_phiRec;
+ mutable Hist::DetArray<Hist::H1D<>> h_ckResAll;
+ mutable Hist::DetArray<Hist::QuadArray<Hist::H1D<>>> h_ckResAllPerQuad;
+ mutable Hist::DetArray<Hist::PanelArray<Hist::H1D<>>> h_ckResAllPerPanel;
+ mutable Hist::DetArray<Hist::H1D<>> h_ckResAllBetaCut;
+ mutable Hist::DetArray<Hist::H1D<>> h_ckResAllDiffSide;
+ mutable Hist::DetArray<Hist::PanelArray<Hist::H1D<>>> h_ckResAllPerPanelBetaCut;
+ mutable InOutArray<Hist::H1D<>> h_ckResAllPerRegion;
+ mutable Hist::PanelArray<InOutArray<Hist::H1D<>>> h_ckResAllPerPanelPerRegion;
+ mutable Hist::DetArray<Hist::PanelArray<ColumnHists<Hist::H1D<>>>> h_ckResAllPerCol;
+ mutable Hist::PanelArray<ColumnHists<InOutArray<Hist::H1D<>>>> h_ckResAllPerColPerRegion;
+ mutable Hist::DetArray<Hist::H2D<>> h_ckResAllVersusPhi;
+ mutable Hist::DetArray<Hist::QuadArray<Hist::H2D<>>> h_ckResAllVersusPhiPerQuad;
+
+ // Phi region histograms
+ struct PhiBinHists {
+ Hist::DetArray<PhiHists<Hist::H1D<>>> ckResAllPhi;
+ Hist::DetArray<Hist::QuadArray<PhiHists<Hist::H1D<>>>> ckResAllPhiPerQuad;
+ };
+ mutable std::unique_ptr<PhiBinHists> h_phiBinHists;
+
+ // Types for mirror specific histograms
+ // Intentionally wrap as a unique pointer and only allocated when we know they are possible
+ // during the event loop, as filling these requires the additional (optional) mirror data from
+ // the photon reconstruction. Use atomicity::none as we handle thread locking ourselves here.
+ template <std::size_t NMIRRORS, Gaudi::Accumulators::atomicity ATOMICITY = Gaudi::Accumulators::atomicity::none>
+ struct MirrorHists {
+ template <typename TYPE>
+ using HistArray = Hist::DetArray<Hist::PanelArray<Hist::Array<TYPE, NMIRRORS>>>;
+ HistArray<Hist::H1D<Hist::DefaultArithmeticType, ATOMICITY>> thetaRec;
+ HistArray<Hist::H1D<Hist::DefaultArithmeticType, ATOMICITY>> phiRec;
+ HistArray<Hist::H1D<Hist::DefaultArithmeticType, ATOMICITY>> ckResAll;
+ HistArray<Hist::H2D<Hist::DefaultArithmeticType, ATOMICITY>> ckResAllVersusPhi;
+ };
+ using PrimMirrorHists = MirrorHists<56u>;
+ using SecMirrorHists = MirrorHists<40u>;
+ mutable std::unique_ptr<PrimMirrorHists> h_primMirrHists;
+ mutable std::unique_ptr<SecMirrorHists> h_secMirrHists;
+ mutable std::mutex m_mirrorLock; ///< Manuyally handle locking for mirror hists
private:
// properties and tools
@@ -142,10 +191,89 @@ namespace Rich::Future::Rec::Moni {
/// Histogram ranges for CK resolution plots
Gaudi::Property<DetectorArray<float>> m_ckResRange{ this, "CKResHistoRange", { 0.0026f, 0.002f } };
+ /// Enable Detailed Histograms
+ Gaudi::Property<bool> m_detailedHists{ this, "DetailedHistograms", false };
+
/** Track selector.
* Longer term should get rid of this and pre-filter the input data instead.
*/
ToolHandle<const ITrackSelector> m_tkSel{ this, "TrackSelector", "TrackSelector" };
+
+ private:
+ void createMirrorHistArrays() const {
+ // As both primary and secondary hists are allocated at the same time only
+ // need to check for one of them.
+ if ( !h_primMirrHists.get() ) {
+ std::scoped_lock lock( m_mirrorLock );
+ if ( !h_primMirrHists.get() ) {
+ h_primMirrHists = std::make_unique<PrimMirrorHists>();
+ h_secMirrHists = std::make_unique<SecMirrorHists>();
+ }
+ }
+ }
+ void initPrimMirrorHists( const Rich::DetectorType rich, const Rich::Side side, const std::size_t mirrN ) const {
+ assert( h_primMirrHists.get() );
+ auto& mirrHists = ( *h_primMirrHists.get() );
+ if ( !mirrHists.ckResAll[rich][side].at( mirrN ).has_value() ) {
+ const auto MirrS = std::to_string( mirrN );
+ const auto rad = radType( rich );
+ bool ok = true;
+ ok &= initHist( mirrHists.thetaRec[rich][side].at( mirrN ),
+ HID( "Mirrors/Primary/Mirror" + MirrS + "/thetaRec", rad ), //
+ "Reconstructed CKTheta - All photons - Primary Mirror " + MirrS, //
+ m_ckThetaMin[rich], m_ckThetaMax[rich], nBins1D(), //
+ "Cherenkov Theta / rad", "Entries" );
+ ok &= initHist( mirrHists.phiRec[rich][side].at( mirrN ), //
+ HID( "Mirrors/Primary/Mirror" + MirrS + "/phiRec", rad ), //
+ "Reconstructed CKPhi - All photons - Primary Mirror " + MirrS, //
+ 0.0, 2.0 * Gaudi::Units::pi, nBins1D(), //
+ "Cherenkov Phi / rad", "Entries" );
+ ok &= initHist( mirrHists.ckResAll[rich][side].at( mirrN ), //
+ HID( "Mirrors/Primary/Mirror" + MirrS + "/ckResAll", rad ), //
+ "Rec-Exp CKTheta - All photons - Primary Mirror " + MirrS, //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
+ "delta(Cherenkov theta) / rad", "Entries" );
+ ok &= initHist( mirrHists.ckResAllVersusPhi[rich][side].at( mirrN ), //
+ HID( "Mirrors/Primary/Mirror" + MirrS + "/ckResAllVersusPhi", rad ), //
+ "Rec-Exp CKTheta Versus CKPhi - All photons - Primary Mirror " + MirrS, //
+ 0, 2.0 * M_PI, nBins2D(), //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins2D(), //
+ "Cherenkov phi / rad", "delta(Cherenkov theta) / rad", "Entries" );
+ if ( !ok ) { throw std::runtime_error( "Failed to initialise primary mirror histograms" ); }
+ }
+ }
+
+ void initSecMirrorHists( const Rich::DetectorType rich, const Rich::Side side, const std::size_t mirrN ) const {
+ assert( h_secMirrHists.get() );
+ auto& mirrHists = ( *h_secMirrHists.get() );
+ if ( !mirrHists.ckResAll[rich][side].at( mirrN ).has_value() ) {
+ const auto MirrS = std::to_string( mirrN );
+ const auto rad = radType( rich );
+ bool ok = true;
+ ok &= initHist( mirrHists.thetaRec[rich][side].at( mirrN ),
+ HID( "Mirrors/Secondary/Mirror" + MirrS + "/thetaRec", rad ), //
+ "Reconstructed CKTheta - All photons - Primary Mirror " + MirrS, //
+ m_ckThetaMin[rich], m_ckThetaMax[rich], nBins1D(), //
+ "Cherenkov Theta / rad", "Entries" );
+ ok &= initHist( mirrHists.phiRec[rich][side].at( mirrN ), //
+ HID( "Mirrors/Secondary/Mirror" + MirrS + "/phiRec", rad ), //
+ "Reconstructed CKPhi - All photons - Primary Mirror " + MirrS, //
+ 0.0, 2.0 * Gaudi::Units::pi, nBins1D(), //
+ "Cherenkov Phi / rad", "Entries" );
+ ok &= initHist( mirrHists.ckResAll[rich][side].at( mirrN ), //
+ HID( "Mirrors/Secondary/Mirror" + MirrS + "/ckResAll", rad ), //
+ "Rec-Exp CKTheta - All photons - Secondary Mirror " + MirrS, //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
+ "delta(Cherenkov theta) / rad", "Entries" );
+ ok &= initHist( mirrHists.ckResAllVersusPhi[rich][side].at( mirrN ), //
+ HID( "Mirrors/Secondary/Mirror" + MirrS + "/ckResAllVersusPhi", rad ), //
+ "Rec-Exp CKTheta Versus CKPhi - All photons - Secondary Mirror " + MirrS, //
+ 0, 2.0 * M_PI, nBins2D(), //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins2D(), //
+ "Cherenkov phi / rad", "delta(Cherenkov theta) / rad", "Entries" );
+ if ( !ok ) { throw std::runtime_error( "Failed to initialise secondary mirror histograms" ); }
+ }
+ }
};
} // namespace Rich::Future::Rec::Moni
@@ -158,9 +286,11 @@ StatusCode SIMDPhotonCherenkovAngles::prebookHistograms() {
bool ok = true;
- // Loop over radiators
- for ( const auto rad : activeRadiators() ) {
- const auto rich = richType( rad );
+ if ( m_detailedHists ) { h_phiBinHists = std::make_unique<PhiBinHists>(); }
+
+ // Loop over RICHes
+ for ( const auto rich : activeDetectors() ) {
+ const auto rad = radType( rich );
// beta cut string
std::string betaS = "beta";
@@ -171,75 +301,118 @@ StatusCode SIMDPhotonCherenkovAngles::prebookHistograms() {
ok &= initHist( h_thetaRec[rich], HID( "thetaRec", rad ),
"Reconstructed CKTheta - All photons", //
m_ckThetaMin[rich], m_ckThetaMax[rich], nBins1D(), //
- "Cherenkov Theta / rad" );
+ "Cherenkov Theta / rad", "Entries" );
ok &= initHist( h_phiRec[rich], HID( "phiRec", rad ),
"Reconstructed CKPhi - All photons", //
0.0, 2.0 * Gaudi::Units::pi, nBins1D(), //
- "Cherenkov Phi / rad" );
+ "Cherenkov Phi / rad", "Entries" );
ok &= initHist( h_ckResAll[rich], HID( "ckResAll", rad ),
"Rec-Exp CKTheta - All photons", //
-m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
+ "delta(Cherenkov theta) / rad", "Entries" );
+ ok &= initHist( h_ckResAllDiffSide[rich], HID( "ckResAllDiffSide", rad ),
+ "Rec-Exp CKTheta - All photons - Cross Sides", //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
+ "delta(Cherenkov theta) / rad", "Entries" );
// with beta cut
ok &= initHist( h_ckResAllBetaCut[rich], HID( "ckResAllBetaCut", rad ), //
"Rec-Exp CKTheta - All photons - " + betaS, //
-m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
+ "delta(Cherenkov theta) / rad", "Entries" );
+
+ // 2D theta versus phi
+ ok &= initHist( h_ckResAllVersusPhi[rich], HID( "ckResAllVersusPhi", rad ), //
+ "Rec-Exp CKTheta Versus CKPhi - All photons", //
+ 0, 2.0 * M_PI, nBins2D(), //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins2D(), //
+ "Cherenkov phi / rad", "delta(Cherenkov theta) / rad", "Entries" );
+
+ // Quadrants
+ for ( const auto q : Rich::Utils::quadrants() ) {
+ ok &= initHist( h_ckResAllPerQuad[rich][q], //
+ HID( "Quadrants/Q" + std::to_string( q ) + "/ckResAll", rad ), //
+ "Rec-Exp CKTheta - All photons" + Rich::Utils::quadString( q ), //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D() );
+ ok &= initHist( h_ckResAllVersusPhiPerQuad[rich][q],
+ HID( "Quadrants/Q" + std::to_string( q ) + "/ckResAllVersusPhi", rad ), //
+ "Rec-Exp CKTheta Versus CKPhi - All photons" + Rich::Utils::quadString( q ), //
+ 0, 2.0 * M_PI, nBins2D(), //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins2D(), //
+ "Cherenkov phi / rad", "delta(Cherenkov theta) / rad", "Entries" );
+ }
+
+ if ( m_detailedHists ) {
+ // Phi Bins
+ for ( std::size_t bin = 0; bin < h_phiBinHists->ckResAllPhi[rich].size(); ++bin ) {
+ const auto binS = std::to_string( bin );
+ const auto binMin = binMinPhi( bin );
+ const auto binMax = binMaxPhi( bin );
+ ok &= ( phiBin( binAvgPhi( bin ) ) == bin );
+ ok &= ( binMin <= binAvgPhi( bin ) && binAvgPhi( bin ) <= binMax );
+ ok &= initHist( h_phiBinHists->ckResAllPhi[rich][bin], HID( "PhiBins/Bin" + binS + "/ckResAll", rad ),
+ "Rec-Exp CKTheta - PhiBin " + binS + " " + std::to_string( binMin ) + "->" +
+ std::to_string( binMax ) + " - All photons",
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), "delta(Cherenkov theta) / rad", "Entries" );
+ for ( const auto q : Rich::Utils::quadrants() ) {
+ ok &=
+ initHist( h_phiBinHists->ckResAllPhiPerQuad[rich][q][bin], //
+ HID( "Quadrants/Q" + std::to_string( q ) + "/PhiBins/Bin" + binS + "/ckResAll", rad ), //
+ "Rec-Exp CKTheta - PhiBin " + binS + " " + std::to_string( binMin ) + "->" +
+ std::to_string( binMax ) + " - All photons" + Rich::Utils::quadString( q ),
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), "delta(Cherenkov theta) / rad", "Entries" );
+ }
+ }
+ }
// loop over detector sides
for ( const auto side : Rich::sides() ) {
ok &= initHist( h_ckResAllPerPanel[rich][side], HID( "ckResAllPerPanel", side, rad ), //
"Rec-Exp CKTheta - All photons", //
-m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
+ "delta(Cherenkov theta) / rad", "Entries" );
// beta cut
- ok &= initHist( h_ckResAllPerPanelBetaCut[rich][side], HID( "ckResAllPerPanelBetaCut", side, rad ), //
- "Rec-Exp CKTheta - All photons - " + betaS, //
- -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
+ ok &= initHist( h_ckResAllPerPanelBetaCut[rich][side], //
+ HID( "ckResAllPerPanelBetaCut", side, rad ), //
+ "Rec-Exp CKTheta - All photons - " + betaS, //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
+ "delta(Cherenkov theta) / rad", "Entries" );
// Per Column resolutions
- const auto nCols = LHCb::RichSmartID::MaPMT::ModuleColumnsPerPanel[rich];
- for ( std::size_t iCol = 0; iCol < nCols; ++iCol ) {
+ for ( std::size_t iCol = 0; iCol < LHCb::RichSmartID::MaPMT::ModuleColumnsPerPanel[rich]; ++iCol ) {
const auto sCol = std::to_string( iCol );
- ok &= initHist( h_ckResAllPerCol[rich][side].at( iCol ), HID( "cols/ckResAllCol" + sCol, side, rad ), //
- "Rec-Exp CKTheta - All photons - Column " + sCol, //
- -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
- // Plots for Inner and outer regions
- ok &= initHist( h_ckResAllPerColInner[rich][side].at( iCol ), //
- HID( "cols/Inner/ckResAllCol" + sCol, side, rad ), //
- "Rec-Exp CKTheta - All photons - Column " + sCol + " - Inner Region", //
- -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
- ok &= initHist( h_ckResAllPerColOuter[rich][side].at( iCol ), //
- HID( "cols/Outer/ckResAllCol" + sCol, side, rad ), //
- "Rec-Exp CKTheta - All photons - Column " + sCol + " - Outer Region", //
- -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
+ ok &= initHist( h_ckResAllPerCol[rich][side].at( iCol ), //
+ HID( "Columns/Col" + sCol + "/ckResAll", side, rad ), //
+ "Rec-Exp CKTheta - All photons - Column " + sCol, //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
+ "delta(Cherenkov theta) / rad", "Entries" );
}
}
- // Make plots for inner and outer regions
- ok &= initHist( h_ckResAllInner[rich], HID( "ckResAllInner", rad ), //
- "Rec-Exp CKTheta - All photons - Inner Region", //
- -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
- ok &= initHist( h_ckResAllOuter[rich], HID( "ckResAllOuter", rad ), //
- "Rec-Exp CKTheta - All photons - Outer Region", //
- -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
- for ( const auto side : Rich::sides() ) {
- ok &= initHist( h_ckResAllPerPanelInner[rich][side], HID( "ckResAllPerPanelInner", side, rad ), //
- "Rec-Exp CKTheta - All photons - Inner Region", //
- -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
- ok &= initHist( h_ckResAllPerPanelOuter[rich][side], HID( "ckResAllPerPanelOuter", side, rad ), //
- "Rec-Exp CKTheta - All photons - Outer Region", //
- -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
- "delta(Cherenkov theta) / rad" );
+ // Make plots for inner and outer regions for RICH2 only (for now, might change in U2)
+ if ( rich == Rich::Rich2 ) {
+ for ( const auto reg : { Inner, Outer } ) {
+ const auto regS = "Region/" + regionString( reg );
+ ok &= initHist( h_ckResAllPerRegion[reg], HID( regS + "/ckResAll", rad ), //
+ "Rec-Exp CKTheta - All photons - " + regS + " Region", //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
+ "delta(Cherenkov theta) / rad", "Entries" );
+ for ( const auto side : Rich::sides() ) {
+ ok &= initHist( h_ckResAllPerPanelPerRegion[side][reg], HID( regS + "/ckResAllPerPanel", side, rad ), //
+ "Rec-Exp CKTheta - All photons - " + regS + " Region", //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
+ "delta(Cherenkov theta) / rad", "Entries" );
+ for ( std::size_t iCol = 0; iCol < LHCb::RichSmartID::MaPMT::ModuleColumnsPerPanel[rich]; ++iCol ) {
+ const auto sCol = std::to_string( iCol );
+ ok &= initHist( h_ckResAllPerColPerRegion[side].at( iCol )[reg], //
+ HID( regS + "/Columns/Col" + sCol + "/ckResAll", side, rad ), //
+ "Rec-Exp CKTheta - All photons - Column " + sCol + " - " + regS + " Region", //
+ -m_ckResRange[rich], m_ckResRange[rich], nBins1D(), //
+ "delta(Cherenkov theta) / rad", "Entries" );
+ }
+ }
+ }
}
- } // active rad loop
+ } // active detector loop
return StatusCode{ ok };
}
@@ -252,27 +425,32 @@ void SIMDPhotonCherenkovAngles::operator()( const LHCb::Track::Range&
const Relations::PhotonToParents::Vector& photToSegPix, //
const LHCb::RichTrackSegment::Vector& segments, //
const CherenkovAngles::Vector& expTkCKThetas, //
- const SIMDCherenkovPhoton::Vector& photons ) const {
+ const SIMDCherenkovPhoton::Vector& photons, //
+ const SIMDMirrorData::Vector& mirrorData ) const {
+
+ const bool hasMirrorData = ( !mirrorData.empty() && mirrorData.size() == photons.size() );
+ if ( m_detailedHists && hasMirrorData ) { createMirrorHistArrays(); }
// local histo buffers
- auto hb_thetaRec = h_thetaRec.buffer();
- auto hb_phiRec = h_phiRec.buffer();
- auto hb_ckResAll = h_ckResAll.buffer();
- auto hb_ckResAllPerPanel = h_ckResAllPerPanel.buffer();
- auto hb_ckResAllInner = h_ckResAllInner.buffer();
- auto hb_ckResAllOuter = h_ckResAllOuter.buffer();
- auto hb_ckResAllPerPanelInner = h_ckResAllPerPanelInner.buffer();
- auto hb_ckResAllPerPanelOuter = h_ckResAllPerPanelOuter.buffer();
- auto hb_ckResAllBetaCut = h_ckResAllBetaCut.buffer();
- auto hb_ckResAllPerPanelBetaCut = h_ckResAllPerPanelBetaCut.buffer();
- auto hb_ckResAllPerCol = h_ckResAllPerCol.buffer();
- auto hb_ckResAllPerColInner = h_ckResAllPerColInner.buffer();
- auto hb_ckResAllPerColOuter = h_ckResAllPerColOuter.buffer();
+ auto hb_thetaRec = h_thetaRec.buffer();
+ auto hb_phiRec = h_phiRec.buffer();
+ auto hb_ckResAll = h_ckResAll.buffer();
+ auto hb_ckResAllPerPanel = h_ckResAllPerPanel.buffer();
+ auto hb_ckResAllPerRegion = h_ckResAllPerRegion.buffer();
+ auto hb_ckResAllPerPanelPerRegion = h_ckResAllPerPanelPerRegion.buffer();
+ auto hb_ckResAllBetaCut = h_ckResAllBetaCut.buffer();
+ auto hb_ckResAllPerPanelBetaCut = h_ckResAllPerPanelBetaCut.buffer();
+ auto hb_ckResAllPerCol = h_ckResAllPerCol.buffer();
+ auto hb_ckResAllPerColPerRegion = h_ckResAllPerColPerRegion.buffer();
+ auto hb_ckResAllVersusPhi = h_ckResAllVersusPhi.buffer();
+ auto hb_ckResAllPerQuad = h_ckResAllPerQuad.buffer();
+ auto hb_ckResAllVersusPhiPerQuad = h_ckResAllVersusPhiPerQuad.buffer();
+ auto hb_ckResAllDiffSide = h_ckResAllDiffSide.buffer();
// loop over the track containers
for ( const auto&& [tk, sumTk] : Ranges::ConstZip( tracks, sumTracks ) ) {
// Is this track selected ?
- if ( !m_tkSel.get()->accept( *tk ) ) continue;
+ if ( !m_tkSel.get()->accept( *tk ) ) { continue; }
// loop over photons for this track
for ( const auto photIn : sumTk.photonIndices() ) {
@@ -306,6 +484,9 @@ void SIMDPhotonCherenkovAngles::operator()( const LHCb::Track::Range&
// expected CK theta
const auto thetaExp = expCKangles[pid];
+ // Quadrant
+ const auto q = Rich::Utils::quadrant( seg.bestPoint().x(), seg.bestPoint().y() );
+
// Loop over scalar entries in SIMD photon
for ( std::size_t i = 0; i < SIMDCherenkovPhoton::SIMDFP::Size; ++i ) {
// Select valid entries
@@ -325,19 +506,54 @@ void SIMDPhotonCherenkovAngles::operator()( const LHCb::Track::Range&
const auto side = id.panel();
// Inner or outer region ?
- const auto isInner = simdPix.isInnerRegion()[i];
+ const auto region = ( simdPix.isInnerRegion()[i] ? Inner : Outer );
- // fill some plots
+ // Segment and photon in same RICH side ?
+ const bool sameSide = ( Rich::Rich1 == rich ? seg.bestPoint().y() * simdPix.gloPos().y()[i] > 0.0
+ : seg.bestPoint().x() * simdPix.gloPos().x()[i] > 0.0 );
+
+ // fill some general plots
++hb_thetaRec[rich][thetaRec];
++hb_phiRec[rich][phiRec];
++hb_ckResAll[rich][deltaTheta];
+ if ( !sameSide ) { ++hb_ckResAllDiffSide[rich][deltaTheta]; }
+ ++hb_ckResAllPerQuad[rich][q][deltaTheta];
+ ++hb_ckResAllVersusPhi[rich][{ phiRec, deltaTheta }];
+ ++hb_ckResAllVersusPhiPerQuad[rich][q][{ phiRec, deltaTheta }];
+ if ( m_detailedHists ) {
+ const auto phiB = phiBin( phiRec );
+ ++( h_phiBinHists->ckResAllPhi[rich].at( phiB )[deltaTheta] );
+ ++( h_phiBinHists->ckResAllPhiPerQuad[rich][q].at( phiB )[deltaTheta] );
+ }
++hb_ckResAllPerPanel[rich][side][deltaTheta];
- ++( isInner ? hb_ckResAllInner[rich] : hb_ckResAllOuter[rich] )[deltaTheta];
- ++( isInner ? hb_ckResAllPerPanelInner[rich][side] : hb_ckResAllPerPanelOuter[rich][side] )[deltaTheta];
+ if ( rich == Rich::Rich2 ) {
+ ++hb_ckResAllPerRegion[region][deltaTheta];
+ ++hb_ckResAllPerPanelPerRegion[side][region][deltaTheta];
+ }
if ( betaC ) {
++hb_ckResAllBetaCut[rich][deltaTheta];
++hb_ckResAllPerPanelBetaCut[rich][side][deltaTheta];
}
+
+ // Mirror plots
+ if ( m_detailedHists && hasMirrorData ) {
+ // Manually thread locking for these hists...
+ std::scoped_lock lock( m_mirrorLock );
+ const auto& mirrordata = mirrorData[photIn];
+ const std::size_t primMirrN = mirrordata.primaryMirrors()[i]->mirrorNumber();
+ const std::size_t secMirrN = mirrordata.secondaryMirrors()[i]->mirrorNumber();
+ initPrimMirrorHists( rich, side, primMirrN );
+ initSecMirrorHists( rich, side, secMirrN );
+ ++h_primMirrHists->thetaRec[rich][side][primMirrN][thetaRec];
+ ++h_secMirrHists->thetaRec[rich][side][secMirrN][thetaRec];
+ ++h_primMirrHists->phiRec[rich][side][primMirrN][phiRec];
+ ++h_secMirrHists->phiRec[rich][side][secMirrN][phiRec];
+ ++h_primMirrHists->ckResAll[rich][side][primMirrN][deltaTheta];
+ ++h_secMirrHists->ckResAll[rich][side][secMirrN][deltaTheta];
+ ++h_primMirrHists->ckResAllVersusPhi[rich][side][primMirrN][{ phiRec, deltaTheta }];
+ ++h_secMirrHists->ckResAllVersusPhi[rich][side][secMirrN][{ phiRec, deltaTheta }];
+ }
+
// FIXME
// Note very old MC samples have a different number of columns.
// Handle here by just silently rejecting those out-of-range.
@@ -346,8 +562,7 @@ void SIMDPhotonCherenkovAngles::operator()( const LHCb::Track::Range&
const std::size_t iCol = id.panelLocalModuleColumn();
if ( iCol < hb_ckResAllPerCol[rich][side].size() ) {
++hb_ckResAllPerCol[rich][side][iCol][deltaTheta];
- ++( isInner ? hb_ckResAllPerColInner[rich][side][iCol]
- : hb_ckResAllPerColOuter[rich][side][iCol] )[deltaTheta];
+ if ( rich == Rich::Rich2 ) { ++hb_ckResAllPerColPerRegion[side][iCol][region][deltaTheta]; }
}
} // valid scalars
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
index 1e229a105d4d6965761e4b9fdf357eaa24a1b9e7..3124c9c9c744b1bce9440e8541eb963a460f9054 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
@@ -67,7 +67,7 @@ namespace Rich::Future::Rec {
#ifdef USE_DD4HEP
m_ckBiasCorrs = { 2.15e-4, 1.8e-5 };
#else
- m_ckBiasCorrs = { 1.2e-4, 5.0e-5 };
+ m_ckBiasCorrs = { 0.6e-4, 0.5e-5 };
#endif
}
@@ -170,6 +170,12 @@ namespace Rich::Future::Rec {
/// Flag to enable the truncation of the output Cherenkov theta and phi values, for each radiator
Gaudi::Property<DetectorArray<bool>> m_ckValueTruncate{ this, "TruncateCKAngles", { true, true } };
+ /** Radiator middle point z fraction. Value gives the fraction of the trajectory from the entry
+ * to the exit point to use for the 'middle' point. So 0.5 means mid point, 0.1 is closer to
+ * the entry point and 0.9 closer to the exit point. */
+ Gaudi::Property<DetectorArray<double>> m_midPointZFrac{
+ this, "MiddlePointZFraction", { 0.5, 0.5 }, "middle point z fraction." };
+
/// Enabled 4D reconstruction
Gaudi::Property<DetectorArray<bool>> m_enable4D{ this, "Enable4D", { false, false }, "Enable 4D reconstruction" };
};
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
index 3943d0726d42b21d5c769f5f3e47f5387ee4db85..957c0f74cfa0759eb2030d291f9110a551a5bd42 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
@@ -810,15 +810,16 @@ SIMDQuarticPhotonReco::getBestGasEmissionPoint( const Rich::RadiatorType radiato
) const {
using namespace LHCb::SIMD;
- // Default point for emission point is the middle
- SIMDFP alongTkFrac = HALF;
-
// Default to all OK
auto ok = SIMDFP::mask_type( true );
// Default emission point is in the middle
emissionPoint = segment.bestPoint();
+ // Default point for emission point
+ const SIMDFP defaultFrac( m_midPointZFrac[radiator] );
+ SIMDFP alongTkFrac = defaultFrac;
+
// which radiator ?
if ( radiator == Rich::Rich1Gas ) {
// First reflection and hit point on same y side ?
@@ -835,8 +836,8 @@ SIMDQuarticPhotonReco::getBestGasEmissionPoint( const Rich::RadiatorType radiato
const auto f = SIMDFP::One() / diff;
fraction( m1 ) = abs( sphReflPoint1.Y() * f );
fraction( m2 ) = abs( sphReflPoint2.Y() * f );
- alongTkFrac( m1 ) = HALF * fraction;
- alongTkFrac( m2 ) = SIMDFP::One() - ( HALF * fraction );
+ alongTkFrac( m1 ) = defaultFrac * fraction;
+ alongTkFrac( m2 ) = SIMDFP::One() - ( defaultFrac * fraction );
emissionPoint = segment.bestPoint( alongTkFrac );
}
} else {
@@ -855,8 +856,8 @@ SIMDQuarticPhotonReco::getBestGasEmissionPoint( const Rich::RadiatorType radiato
const auto f = SIMDFP::One() / diff;
fraction( m1 ) = abs( sphReflPoint1.X() * f );
fraction( m2 ) = abs( sphReflPoint2.X() * f );
- alongTkFrac( m1 ) = HALF * fraction;
- alongTkFrac( m2 ) = SIMDFP::One() - ( HALF * fraction );
+ alongTkFrac( m1 ) = defaultFrac * fraction;
+ alongTkFrac( m2 ) = SIMDFP::One() - ( defaultFrac * fraction );
emissionPoint = segment.bestPoint( alongTkFrac );
}
}
diff --git a/Rich/RichFutureRecSys/tests/options/rich-dst-reco.py b/Rich/RichFutureRecSys/tests/options/rich-dst-reco.py
index 74a2ec3473d1bb1c62b50ea990d2cba483182df2..1d7690cbfd882000a1062fb83e5d9cbc8557cd30 100644
--- a/Rich/RichFutureRecSys/tests/options/rich-dst-reco.py
+++ b/Rich/RichFutureRecSys/tests/options/rich-dst-reco.py
@@ -99,7 +99,7 @@ app = ApplicationMgr(
ExtSvc=[
ToolSvc(),
AuditorSvc(),
- MessageSvcSink(),
+ MessageSvcSink(HistoStringsWidth=80),
RootHistoSink(FileName=myName + "-hists.root"),
],
AuditAlgorithms=True,