Extended CaloHistogramSink to Profile histograms

Also allows to either store them as is, or convert to a pair of regular histogrms for mean and rms

CaloFuture/CaloFutureMoniDst/src/CaloHistogramSink.cpp

@@ -17,10 +17,12 @@
 
 #include <TH2D.h>
 #include <TH3D.h>
+#include <TProfile2D.h>
+#include <TProfile3D.h>
 
-namespace {
+namespace GHSink = Gaudi::Histograming::Sink;
 
-  using namespace Gaudi::Histograming::Sink;
+namespace {
 
   /// description of Calorimeter geometry -- FIXME: READ FROM GEOMETRY
   struct CaloParams {
@@ -34,10 +36,55 @@ namespace {
   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      = TH2D;
-    using WeightType = double;
+    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 ) {
@@ -46,7 +93,8 @@ namespace {
       int         nbRow = ( cp.centre - cp.firstRow ) * 2;
       auto        xmax  = cp.maxCellSizes * nbCol / 2.;
       auto        ymax  = cp.maxCellSizes * nbRow / 2.;
-      return TH2D{name.c_str(), title.c_str(), nbCol * cp.reg, -xmax, xmax, nbRow * cp.reg, -ymax, ymax};
+      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 ) {
@@ -67,21 +115,26 @@ namespace {
         for ( unsigned int ic = 0; ic < ibox; ic++ ) {
           int    iic = ( id.col() - cp.centre ) * ibox + ic;
           double x   = cellSize * ( iic + 0.5 );
-          histo.SetBinContent( histo.FindBin( x, y ), weight );
+          auto   bin = histo.FindBin( x, y );
+          updateHistoBin<flavour>( histo, bin, weight );
         }
       }
     }
 
-    static void fillMetaData( TH2D& histo, nlohmann::json const&, unsigned int nentries ) {
+    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      = TH3D;
-    using WeightType = double;
+    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 ) {
@@ -90,11 +143,20 @@ namespace {
       int         nbRow = ( cp.centre - cp.firstRow ) * 2;
       auto        xmax  = cp.maxCellSizes * nbCol / 2.;
       auto        ymax  = cp.maxCellSizes * nbRow / 2.;
-      return TH3D{name.c_str(), title.c_str(), nbCol * cp.reg,   -xmax,          xmax,          nbRow * cp.reg,
-                  -ymax,        ymax,          (int)axis1.nBins, axis1.minValue, axis1.maxValue};
+      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( TH3D& histo, unsigned int i, const WeightType& weight ) {
+    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();
@@ -118,39 +180,69 @@ namespace {
           int    iic    = ( id.col() - cp.centre ) * ibox + ic;
           double x      = cellSize * ( iic + 0.5 );
           int    xIndex = histo.GetXaxis()->FindBin( x );
-          histo.SetBinContent( histo.GetBin( xIndex, yIndex, index3rdD ), weight );
+          auto   bin    = histo.GetBin( xIndex, yIndex, index3rdD );
+          updateHistoBin<flavour>( histo, bin, weight );
         }
       }
     }
 
-    static void fillMetaData( TH3D& histo, nlohmann::json const&, unsigned int nentries ) {
+    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;
 
-  struct Root : public Gaudi::Histograming::Sink::Base {
-
+  struct Root : public GHSink::Base {
     using Base::Base;
-
-    HistoRegistry const registry = {
-        {{"histogram:CellIDHistogram"s, 1}, &Gaudi::Histograming::Sink::saveRootHisto<Traits1D>},
-        {{"histogram:WeightedCellIDHistogram"s, 1}, &Gaudi::Histograming::Sink::saveRootHisto<Traits1D>},
-        {{"histogram:CellIDHistogram"s, 2}, &Gaudi::Histograming::Sink::saveRootHisto<Traits2D>},
-        {{"histogram:WeightedCellIDHistogram"s, 2}, &Gaudi::Histograming::Sink::saveRootHisto<Traits2D>},
+    template <int ND>
+    struct saveProfileRootHistoWrapper {
+      saveProfileRootHistoWrapper( Root& root ) : m_root( root ) {}
+      void operator()( TFile& file, std::string dir, std::string name, nlohmann::json const& j ) {
+        return saveProfileRootHisto<ND>( file, dir, name, j, m_root.m_convertProfiles.value() );
+      }
+      Root& m_root;
     };
-
+    HistoRegistry const registry = {
+        {{"histogram:CellIDHistogram"s, 1}, &GHSink::saveRootHisto<Traits1D<HistoFlavour::Regular>>},
+        {{"histogram:WeightedCellIDHistogram"s, 1}, &GHSink::saveRootHisto<Traits1D<HistoFlavour::Regular>>},
+        {{"histogram:ProfileCellIDHistogram"s, 1}, saveProfileRootHistoWrapper<1>( *this )},
+        {{"histogram:WeightedProfileCellIDHistogram"s, 1}, saveProfileRootHistoWrapper<1>( *this )},
+        {{"histogram:CellIDHistogram"s, 2}, &GHSink::saveRootHisto<Traits2D<HistoFlavour::Regular>>},
+        {{"histogram:WeightedCellIDHistogram"s, 2}, &GHSink::saveRootHisto<Traits2D<HistoFlavour::Regular>>},
+        {{"histogram:ProfileCellIDHistogram"s, 2}, saveProfileRootHistoWrapper<2>( *this )},
+        {{"histogram:WeightedProfileCellIDHistogram"s, 2}, saveProfileRootHistoWrapper<2>( *this )}};
     StatusCode initialize() override {
       return Base::initialize().andThen( [&] {
         for ( auto& [id, func] : registry ) { registerHandler( id, func ); }
       } );
     }
+    Gaudi::Property<bool> m_convertProfiles{
+        this, "ConvertProfileHistos", true,
+        "If true, no profile histos are saved. These are converted to a pair of regular histograms, one plotting the "
+        "mean and the other the rms of the original data"};
   };
 
   DECLARE_COMPONENT( Root )