diff --git a/CaloFuture/CaloFutureMoniDst/src/CaloHistogramSink.cpp b/CaloFuture/CaloFutureMoniDst/src/CaloHistogramSink.cpp
index 5a9634acd8fe68c7b1c47fd46d96954540db74b9..31af8adf5af2eb97ef719d1caea717bd1ec05fd6 100644
--- a/CaloFuture/CaloFutureMoniDst/src/CaloHistogramSink.cpp
+++ b/CaloFuture/CaloFutureMoniDst/src/CaloHistogramSink.cpp
@@ -9,208 +9,12 @@
* or submit itself to any jurisdiction. *
\***********************************************************************************/
-#include <CaloFutureUtils/CellIDHistogram.h>
-#include <Detector/Calo/CaloCellID.h>
+#include <CaloFutureUtils/CellIDHistogramUtils.h>
#include <HistogramPersistencySvc/RootHistogramSinkBase.h>
-#include <HistogramPersistencySvc/RootHistogramUtils.h>
-
-#include <TH2D.h>
-#include <TH3D.h>
-#include <TProfile2D.h>
-#include <TProfile3D.h>
namespace GHSink = Gaudi::Histograming::Sink;
-namespace {
-
- /// description of Calorimeter geometry -- FIXME: READ FROM GEOMETRY
- struct CaloParams {
- int centre; // centre of the Calorimeter in nb rows/cols
- int reg; // nb divisions of biggest cell used in the histogram
- int firstRow; // nb empty rows at the top/bottom, as the calorimeters are not squares
- float maxCellSizes; // size in mm of the biggest cell
- };
- constexpr CaloParams EcalParams = {32, 6, 6, 121.2};
- constexpr CaloParams HcalParams = {16, 2, 3, 262.6};
- constexpr unsigned int NEcalRootHistBins = ( EcalParams.centre * 2 * EcalParams.reg + 2 ) *
- ( ( EcalParams.centre - EcalParams.firstRow ) * 2 * EcalParams.reg + 2 );
-
- /// Small enum for the different type of histos a Trait can deal with
- enum class HistoFlavour {
- Regular, // plain, regular histogram, meaning non profile
- Profile, // profile histogram, to be kept as such
- ProfileMean, // profile histogram, to be converted to a regular one hosting the mean values
- ProfileRMS // profile histogram, to be converted to a regular one hosting the RMS values
- };
-
- /// patches a name depending on the flavour
- template <HistoFlavour flavour>
- std::string patchName( std::string& name ) {
- if constexpr ( flavour == HistoFlavour::ProfileMean ) return name + "_Mean";
- if constexpr ( flavour == HistoFlavour::ProfileRMS ) return name + "_RMS";
- return name;
- }
-
- /// update a histogram bin, according to the given flavour
- template <HistoFlavour flavour, typename Histo, typename WeightType>
- void updateHistoBin( Histo& histo, int bin, WeightType weight ) {
- if constexpr ( flavour == HistoFlavour::Regular ) {
- histo.SetBinContent( bin, weight );
- } else {
- auto [c, sumWeight2] = weight;
- auto [nEntries, sumWeight] = c;
- if constexpr ( flavour == HistoFlavour::ProfileRMS ) {
- auto rms = 0;
- if ( nEntries > 0 ) {
- auto variance = ( sumWeight2 - sumWeight * ( sumWeight / nEntries ) ) / nEntries;
- rms = ( variance > 0 ) ? sqrt( variance ) : 0;
- }
- histo.SetBinContent( bin, rms );
- } else if constexpr ( flavour == HistoFlavour::ProfileMean ) {
- histo.SetBinContent( bin, nEntries > 0 ? sumWeight / nEntries : 0 );
- } else if constexpr ( flavour == HistoFlavour::Profile ) {
- histo.setBinNEntries( bin, nEntries );
- histo.SetBinContent( bin, sumWeight );
- histo.setBinW2( bin, sumWeight2 );
- }
- }
- }
-
- /// Traits handling 1D CellID based histograms, all flavours
- template <HistoFlavour flavour = HistoFlavour::Regular>
- struct Traits1D {
- static constexpr unsigned int Dimension{1};
- using Histo =
- std::conditional_t<flavour == HistoFlavour::Profile, GHSink::details::ProfileWrapper<TProfile2D>, TH2D>;
- using WeightType = std::conditional_t<flavour != HistoFlavour::Regular,
- std::tuple<std::tuple<unsigned int, double>, double>, double>;
-
- template <typename Axis>
- static Histo create( std::string& name, std::string& title, Axis& axis ) {
- auto const& cp = axis.minValue == 0 ? EcalParams : HcalParams;
- int nbCol = cp.centre * 2;
- int nbRow = ( cp.centre - cp.firstRow ) * 2;
- auto xmax = cp.maxCellSizes * nbCol / 2.;
- auto ymax = cp.maxCellSizes * nbRow / 2.;
- return Histo{
- patchName<flavour>( name ).c_str(), title.c_str(), nbCol * cp.reg, -xmax, xmax, nbRow * cp.reg, -ymax, ymax};
- }
-
- static auto fill( Histo& histo, unsigned int i, const WeightType& weight ) {
- bool isEcal = histo.GetNcells() == NEcalRootHistBins;
- auto const& cp = isEcal ? EcalParams : HcalParams;
- unsigned int const nbCells =
- isEcal ? LHCb::Detector::Calo::Index::nbEcalCells() : LHCb::Detector::Calo::Index::nbHcalCells();
- // ignore under / overflow bin, they are not used by construction
- if ( i == 0 || i == nbCells ) return;
- int cellNumber = i - 1 + ( isEcal ? 0 : LHCb::Detector::Calo::Index::nbEcalCells() );
- LHCb::Detector::Calo::CellID id = LHCb::Detector::Calo::DenseIndex::details::toCellID( cellNumber );
- // loop over cell area ( ibox * ibox bins depending on CaloFuture/Area)
- const unsigned int ibox = cp.reg / ( id.area() + 1 );
- auto cellSize = cp.maxCellSizes / cp.reg;
- for ( unsigned int ir = 0; ir < ibox; ir++ ) {
- int iir = ( id.row() - cp.centre ) * ibox + ir;
- double y = cellSize * ( iir + 0.5 );
- for ( unsigned int ic = 0; ic < ibox; ic++ ) {
- int iic = ( id.col() - cp.centre ) * ibox + ic;
- double x = cellSize * ( iic + 0.5 );
- auto bin = histo.FindBin( x, y );
- updateHistoBin<flavour>( histo, bin, weight );
- }
- }
- }
-
- static void fillMetaData( Histo& histo, nlohmann::json const&, unsigned int nentries ) {
- auto const& cp = histo.GetNcells() == NEcalRootHistBins ? EcalParams : HcalParams;
- histo.SetEntries( nentries * cp.reg * cp.reg );
- }
- };
-
- /// Traits handling 2D CellID based histograms, all flavours
- template <HistoFlavour flavour = HistoFlavour::Regular>
- struct Traits2D {
- static constexpr unsigned int Dimension{2};
- using Histo =
- std::conditional_t<flavour == HistoFlavour::Profile, GHSink::details::ProfileWrapper<TProfile3D>, TH3D>;
- using WeightType = std::conditional_t<flavour != HistoFlavour::Regular,
- std::tuple<std::tuple<unsigned int, double>, double>, double>;
-
- template <typename Axis0, typename Axis1>
- static Histo create( std::string& name, std::string& title, Axis0& axis0, Axis1& axis1 ) {
- auto const& cp = axis0.minValue == 0 ? EcalParams : HcalParams;
- int nbCol = cp.centre * 2;
- int nbRow = ( cp.centre - cp.firstRow ) * 2;
- auto xmax = cp.maxCellSizes * nbCol / 2.;
- auto ymax = cp.maxCellSizes * nbRow / 2.;
- return Histo{patchName<flavour>( name ).c_str(),
- title.c_str(),
- nbCol * cp.reg,
- -xmax,
- xmax,
- nbRow * cp.reg,
- -ymax,
- ymax,
- (int)axis1.nBins,
- axis1.minValue,
- axis1.maxValue};
- }
-
- static auto fill( Histo& histo, unsigned int i, const WeightType& weight ) {
- bool isEcal = histo.GetNcells() == NEcalRootHistBins;
- unsigned int const nbCells =
- isEcal ? LHCb::Detector::Calo::Index::nbEcalCells() : LHCb::Detector::Calo::Index::nbHcalCells();
- // first extract value in the non CellID dimension (so 3rd dimension of the resulting histo)
- auto index3rdD = i / ( nbCells + 2 );
- // then extract value for the cellID dimension
- i = i % ( nbCells + 2 );
- // ignore under / overflow bin, they are not used by construction for cellIDs
- if ( i == 0 || i == nbCells ) return;
- int cellNumber = i - 1 + ( isEcal ? 0 : LHCb::Detector::Calo::Index::nbEcalCells() );
- LHCb::Detector::Calo::CellID id = LHCb::Detector::Calo::DenseIndex::details::toCellID( cellNumber );
- // loop over cell area ( ibox * ibox bins depending on CaloFuture/Area)
- auto const& cp = isEcal ? EcalParams : HcalParams;
- const unsigned int ibox = cp.reg / ( id.area() + 1 );
- auto cellSize = cp.maxCellSizes / cp.reg;
- for ( unsigned int ir = 0; ir < ibox; ir++ ) {
- int iir = ( id.row() - cp.centre ) * ibox + ir;
- double y = cellSize * ( iir + 0.5 );
- int yIndex = histo.GetYaxis()->FindBin( y );
- for ( unsigned int ic = 0; ic < ibox; ic++ ) {
- int iic = ( id.col() - cp.centre ) * ibox + ic;
- double x = cellSize * ( iic + 0.5 );
- int xIndex = histo.GetXaxis()->FindBin( x );
- auto bin = histo.GetBin( xIndex, yIndex, index3rdD );
- updateHistoBin<flavour>( histo, bin, weight );
- }
- }
- }
-
- static void fillMetaData( Histo& histo, nlohmann::json const&, unsigned int nentries ) {
- auto const& cp = histo.GetNcells() == NEcalRootHistBins ? EcalParams : HcalParams;
- histo.SetEntries( nentries * cp.reg * cp.reg );
- }
- };
-
- /// Dedicated handler for Profile Histograms, respecting conversion to standard ones when required
- template <int ND>
- void saveProfileRootHisto( TFile& file, std::string dir, std::string name, nlohmann::json const& j,
- bool convertProfiles ) {
- if ( convertProfiles ) {
- // when converting, we save 2 histograms, one for mean, one for RMS
- using TrMean =
- std::conditional_t<ND == 1, Traits1D<HistoFlavour::ProfileMean>, Traits2D<HistoFlavour::ProfileMean>>;
- using TrRMS = std::conditional_t<ND == 1, Traits1D<HistoFlavour::ProfileRMS>, Traits2D<HistoFlavour::ProfileRMS>>;
- GHSink::saveRootHisto<TrMean>( file, dir, name, j );
- GHSink::saveRootHisto<TrRMS>( file, dir, name, j );
- } else {
- using Tr = std::conditional_t<ND == 1, Traits1D<HistoFlavour::Profile>, Traits2D<HistoFlavour::Profile>>;
- GHSink::saveRootHisto<Tr>( file, dir, name, j );
- }
- }
-
-} // namespace
-
namespace LHCb::Calo::Histograming::Sink {
using namespace std::string_literals;