Histograming/Sink/Utils.h

0001 /***********************************************************************************\
0002 * (c) Copyright 1998-2023 CERN for the benefit of the LHCb and ATLAS collaborations *
0003 *                                                                                   *
0004 * This software is distributed under the terms of the Apache version 2 licence,     *
0005 * copied verbatim in the file "LICENSE".                                            *
0006 *                                                                                   *
0007 * In applying this licence, CERN does not waive the privileges and immunities       *
0008 * granted to it by virtue of its status as an Intergovernmental Organization        *
0009 * or submit itself to any jurisdiction.                                             *
0010 \***********************************************************************************/
0011 #pragma once
0012
0013 #include <GaudiKernel/GaudiException.h>
0014 #include <TDirectory.h>
0015 #include <TFile.h>
0016 #include <TH1.h>
0017 #include <TH1D.h>
0018 #include <TH2D.h>
0019 #include <TH3D.h>
0020 #include <TProfile.h>
0021 #include <TProfile2D.h>
0022 #include <algorithm>
0023 #include <functional>
0024 #include <gsl/span>
0025 #include <nlohmann/json.hpp>
0026 #include <range/v3/numeric/accumulate.hpp>
0027 #include <range/v3/range/conversion.hpp>
0028 #include <range/v3/view/split_when.hpp>
0029 #include <range/v3/view/transform.hpp>
0030 #include <string>
0031 #include <vector>
0032
0033 // upstream has renamed namespace ranges::view -> ranges::views
0034 #if RANGE_V3_VERSION < 900
0035 namespace ranges::views {
0036   using namespace ranges::view;
0037 }
0038 #endif
0039
0040 namespace Gaudi::Histograming::Sink {
0041
0042   /**
0043    * templated Traits dealing with Root Histogram filling for standard histograms
0044    *
0045    * Specializatoins of the Traits type must provide the following static functions :
0046    *   - Histo create( std::string& name, std::string& title, Axis& axis... )
0047    *       it should intanciate a new ROOT histogram instance and return it
0048    *   - void fillMetaData( Histo& histo, nlohmann::json const& jsonAxis, unsigned int nentries)
0049    *       it should fill metadata in the ROOT histo histogram provides from the list of axis number of entries
0050    *   - void fill( Histo& histo, unsigned int i, const WeightType& weight )
0051    *       it should fill the given bin with the given value
0052    * Last point, it should have a static constexpr unsigned int Dimension
0053    */
0054   template <bool isProfile, typename RootHisto, unsigned int N>
0055   struct Traits;
0056
0057   /**
0058    * generic function to convert json to a ROOT Histogram
0059    *
0060    * returns the Root histogram and the dir where to save it in the Root file
0061    * This may be different from input dir in case name has slashes
0062    */
0063   template <typename Traits>
0064   std::tuple<typename Traits::Histo, std::string> jsonToRootHistogram( std::string& dir, std::string& name,
0065                                                                        nlohmann::json const& j );
0066
0067   /**
0068    * generic function to convert a ROOT Histogram to json
0069    *
0070    * essentially used for backward compatibility of old HistogramService with MonitoringHub
0071    */
0072   template <typename Histo>
0073   nlohmann::json rootHistogramTojson( Histo const& );
0074
0075   /**
0076    * generic method to save histograms to files, based on Traits
0077    */
0078   template <typename Traits>
0079   void saveRootHisto( TFile& file, std::string dir, std::string name, nlohmann::json const& j );
0080
0081   /**
0082    * generic method to save regular histograms to files
0083    *
0084    * Can be used in most cases as the handler function to register into Sink::Base
0085    * contains all the boiler plate code and redirects specific code to the adapted Traits template
0086    */
0087   template <unsigned int N, bool isProfile, typename ROOTHisto>
0088   void saveRootHisto( TFile& file, std::string dir, std::string name, nlohmann::json const& j );
0089
0090   namespace details {
0091
0092     /// Small helper struct representing the Axis of an Histogram
0093     struct Axis {
0094       unsigned int nBins;
0095       double       minValue;
0096       double       maxValue;
0097       std::string  title;
0098     };
0099
0100     /// extract an Axis from json data
0101     inline Axis jsonToAxis( nlohmann::json& jAxis ) {
0102       return { jAxis.at( "nBins" ).get<unsigned int>(), jAxis.at( "minValue" ).get<double>(),
0103                jAxis.at( "maxValue" ).get<double>(),
0104                ";" + jAxis.at( "title" ).get<std::string>() }; // ";" to prepare concatenations of titles
0105     }
0106
0107     /**
0108      * generic function to convert json to a ROOT Histogram - internal implementation
0109      */
0110     template <typename Traits, std::size_t... index>
0111     std::tuple<typename Traits::Histo, std::string> jsonToRootHistogramInternal( std::string& dir, std::string& name,
0112                                                                                  nlohmann::json const& j,
0113                                                                                  std::index_sequence<index...> ) {
0114       // extract data from json
0115       auto jsonAxis = j.at( "axis" );
0116       auto axis     = std::array{ jsonToAxis( jsonAxis[index] )... };
0117       auto weights  = j.at( "bins" ).get<std::vector<typename Traits::WeightType>>();
0118       auto title    = j.at( "title" ).get<std::string>();
0119       auto nentries = j.at( "nEntries" ).get<unsigned int>();
0120       // weird way ROOT has to give titles to axis
0121       title += ( axis[index].title + ... );
0122       // compute total number of bins, multiplying bins per axis
0123       auto totNBins = ( ( axis[index].nBins + 2 ) * ... );
0124       assert( weights.size() == totNBins );
0125
0126       if ( name[0] == '/' ) {
0127         dir  = "";
0128         name = name.substr( 1 );
0129       }
0130
0131       // take into account the case where name contains '/'s (e.g. "Group/Name") by
0132       // moving the prefix into dir
0133       if ( auto pos = name.rfind( '/' ); pos != std::string::npos ) {
0134         dir += '/' + name.substr( 0, pos );
0135         name = name.substr( pos + 1 );
0136       }
0137
0138       // Create Root histogram calling constructors with the args tuple
0139       auto histo = Traits::create( name, title, axis[index]... );
0140
0141       // fill Histo
0142       for ( unsigned int i = 0; i < totNBins; i++ ) Traits::fill( histo, i, weights[i] );
0143       // in case we have sums, overwrite them in the histogram with our more precise values
0144       // FIXME This is only supporting regular histograms. It won't work in case of weighted histograms
0145       if constexpr ( sizeof...( index ) == 1 ) {
0146         if ( j.find( "sum" ) != j.end() ) {
0147           double s[13];
0148           s[0] = j.at( "nTotEntries" ).get<double>();
0149           s[1] = j.at( "nTotEntries" ).get<double>();
0150           s[2] = j.at( "sum" ).get<double>();
0151           s[3] = j.at( "sum2" ).get<double>();
0152           histo.PutStats( s );
0153         }
0154       } else {
0155         if ( j.find( "sumx" ) != j.end() ) {
0156           double s[13];
0157           s[0] = j.at( "nTotEntries" ).get<double>();
0158           s[1] = j.at( "nTotEntries" ).get<double>();
0159           s[2] = j.at( "sumx" ).get<double>();
0160           s[3] = j.at( "sumx2" ).get<double>();
0161           s[4] = j.at( "sumy" ).get<double>();
0162           s[5] = j.at( "sumy2" ).get<double>();
0163           s[6] = j.at( "sumxy" ).get<double>();
0164           if constexpr ( sizeof...( index ) > 2 ) {
0165             s[7]  = j.at( "sumz" ).get<double>();
0166             s[8]  = j.at( "sumz2" ).get<double>();
0167             s[9]  = j.at( "sumxz" ).get<double>();
0168             s[10] = j.at( "sumyz" ).get<double>();
0169           }
0170           histo.PutStats( s );
0171         }
0172       }
0173       // fill histo metadata, e.g. bins and number of entries
0174       Traits::fillMetaData( histo, jsonAxis, nentries );
0175       return { histo, dir };
0176     }
0177
0178     /**
0179      * Common base for Traits dealing with Histogram conversions to Root
0180      * Provides generic implementation for creating the histogram and filling meta data
0181      * The filling (method fill) is not implemented
0182      */
0183     template <typename RootHisto, unsigned int N>
0184     struct TraitsBase {
0185       static constexpr unsigned int Dimension{ N };
0186       // hleper method for the creation of the Root histogram, using an index_sequence for handling the axis
0187       template <typename AxisArray, std::size_t... index>
0188       static RootHisto create( std::string& name, std::string& title, AxisArray axis, std::index_sequence<index...> ) {
0189         return std::make_from_tuple<RootHisto>(
0190             std::tuple_cat( std::tuple{ name.c_str(), title.c_str() },
0191                             std::tuple{ std::get<index>( axis ).nBins, std::get<index>( axis ).minValue,
0192                                         std::get<index>( axis ).maxValue }... ) );
0193       }
0194       // Generic creation of the Root histogram from name, title and a set of axis
0195       template <typename... Axis>
0196       static RootHisto create( std::string& name, std::string& title, Axis&... axis ) {
0197         return create( name, title, std::make_tuple( axis... ), std::make_index_sequence<sizeof...( Axis )>() );
0198       }
0199       // Generic fill of the metadata of the Root histogram from json data
0200       static void fillMetaData( RootHisto& histo, nlohmann::json const& jsonAxis, unsigned int nentries ) {
0201         auto try_set_bin_labels = [&histo, &jsonAxis]( auto idx ) {
0202           TAxis* axis = nullptr;
0203           switch ( idx ) {
0204           case 0:
0205             axis = histo.GetXaxis();
0206             break;
0207           case 1:
0208             axis = histo.GetYaxis();
0209             break;
0210           case 2:
0211             axis = histo.GetZaxis();
0212             break;
0213           default:
0214             break;
0215           }
0216           if ( axis && jsonAxis[idx].contains( "labels" ) ) {
0217             const auto labels = jsonAxis[idx].at( "labels" );
0218             for ( unsigned int i = 0; i < labels.size(); i++ ) {
0219               axis->SetBinLabel( i + 1, labels[i].template get<std::string>().c_str() );
0220             }
0221           }
0222           if ( axis && jsonAxis[idx].contains( "xbins" ) ) {
0223             const auto xbins = jsonAxis[idx].at( "xbins" );
0224             axis->Set( xbins.size() - 1, xbins.template get<std::vector<double>>().data() );
0225           }
0226         };
0227         for ( unsigned int i = 0; i < jsonAxis.size(); i++ ) try_set_bin_labels( i );
0228         // Fix the number of entries, wrongly computed by ROOT from the numer of calls to fill
0229         histo.SetEntries( nentries );
0230       }
0231     };
0232
0233     /// helper Wrapper around TProfileX to be able to fill it
0234     template <typename TP>
0235     struct ProfileWrapper : TP {
0236       template <typename... Args>
0237       ProfileWrapper( Args&&... args ) : TP( std::forward<Args>( args )... ) {
0238         // When the static function TH1::SetDefaultSumw2 had been called (passing true), `Sumw2(true)` is automatically
0239         // called by the constructor of TProfile. This is a problem because in the profiles in Gaudi::Accumulators we do
0240         // not keep track of the sum of squares of weights and consequently don't set fBinSumw2 of the TProfile, which
0241         // in turn leads to a self-inconsistent TProfile object. The "fix" is to disable sum of squares explicitly.
0242         this->Sumw2( false );
0243       }
0244       ProfileWrapper( ProfileWrapper const& )             = delete;
0245       ProfileWrapper& operator=( ProfileWrapper const& )  = delete;
0246       ProfileWrapper( ProfileWrapper const&& )            = delete;
0247       ProfileWrapper& operator=( ProfileWrapper const&& ) = delete;
0248       void            setBinNEntries( Int_t i, Int_t n ) { this->fBinEntries.fArray[i] = n; }
0249       void            setBinW2( Int_t i, Double_t v ) { this->fSumw2.fArray[i] = v; }
0250     };
0251
0252     /**
0253      * changes to the ROOT directory given in the current ROOT file and returns
0254      * the current directory before the change
0255      */
0256     inline TDirectory* changeDir( TFile& file, std::string dir ) {
0257       // remember the current directory
0258       auto previousDir = gDirectory;
0259       // find or create the directory for the histogram
0260       using namespace ranges;
0261       auto is_delimiter        = []( auto c ) { return c == '/' || c == '.'; };
0262       auto transform_to_string = views::transform( []( auto&& rng ) { return rng | to<std::string>; } );
0263       auto currentDir          = accumulate( dir | views::split_when( is_delimiter ) | transform_to_string,
0264                                              file.GetDirectory( "" ), []( auto current, auto&& dir_level ) {
0265                                       if ( current ) {
0266                                         // try to get next level
0267                                         auto nextDir = current->GetDirectory( dir_level.c_str() );
0268                                         // if it does not exist, create it
0269                                         if ( !nextDir ) nextDir = current->mkdir( dir_level.c_str() );
0270                                         // move to next level
0271                                         current = nextDir;
0272                                       }
0273                                       return current;
0274                                     } );
0275       if ( !currentDir )
0276         throw GaudiException( "Could not create directory " + dir, "Histogram::Sink::Root", StatusCode::FAILURE );
0277       // switch to the directory
0278       currentDir->cd();
0279       return previousDir;
0280     }
0281
0282     template <typename Histo>
0283     nlohmann::json allAxisTojson( Histo const& h ) {
0284       if constexpr ( std::is_base_of_v<TH3D, Histo> ) {
0285         return { *h.GetXaxis(), *h.GetYaxis(), *h.GetZaxis() };
0286       } else if constexpr ( std::is_base_of_v<TProfile2D, Histo> || std::is_base_of_v<TH2D, Histo> ) {
0287         return { *h.GetXaxis(), *h.GetYaxis() };
0288       } else {
0289         return { *h.GetXaxis() };
0290       }
0291     }
0292
0293     template <typename Histo>
0294     nlohmann::json binsTojson( Histo const& h ) {
0295       if constexpr ( std::is_base_of_v<TProfile, Histo> || std::is_base_of_v<TProfile2D, Histo> ) {
0296         nlohmann::json j;
0297         // ROOT TProfile interface being completely inconsistent, we have to play
0298         // with different ways to access values of the histogram... one per value !
0299         auto* sums  = h.GetArray();
0300         auto* sums2 = h.GetSumw2();
0301         for ( unsigned long n = 0; n < (unsigned long)h.GetSize(); n++ ) {
0302           j.push_back( nlohmann::json{ { h.GetBinEntries( n ), sums[n] }, sums2->At( n ) } );
0303         }
0304         return j;
0305       } else {
0306         return gsl::span{ h.GetArray(), (unsigned long)h.GetSize() };
0307       }
0308     }
0309
0310     /// automatic translation of Root Histograms to json
0311     template <typename Histo>
0312     nlohmann::json rootHistogramToJson( Histo const& h ) {
0313       bool        isProfile = std::is_base_of_v<TProfile, Histo> || std::is_base_of_v<TProfile2D, Histo>;
0314       std::string type      = isProfile ? "histogram:ProfileHistogram:double" : "histogram:Histogram:double";
0315       auto        j         = nlohmann::json{ { "type", type },
0316                                               { "title", h.GetTitle() },
0317                                               { "dimension", h.GetDimension() },
0318                                               { "empty", (int)h.GetEntries() == 0 },
0319                                               { "nEntries", (int)h.GetEntries() },
0320                                               { "axis", allAxisTojson( h ) },
0321                                               { "bins", binsTojson( h ) } };
0322       if ( !isProfile ) {
0323         double s[13];
0324         h.GetStats( s );
0325         if ( h.GetDimension() == 1 ) {
0326           j["nTotEntries"] = s[0];
0327           j["sum"]         = s[2];
0328           j["sum2"]        = s[3];
0329           j["mean"]        = s[2] / s[0];
0330         } else {
0331           j["nTotEntries"] = s[0];
0332           j["sumx"]        = s[2];
0333           j["sumx2"]       = s[3];
0334           j["meanx"]       = s[2] / s[0];
0335           j["sumy"]        = s[4];
0336           j["sumy2"]       = s[5];
0337           j["sumxy"]       = s[6];
0338           j["meany"]       = s[4] / s[0];
0339           if ( h.GetDimension() >= 3 ) {
0340             j["sumz"]  = s[7];
0341             j["sumz2"] = s[8];
0342             j["sumxz"] = s[9];
0343             j["sumyz"] = s[10];
0344             j["meanz"] = s[7] / s[0];
0345           }
0346         }
0347       }
0348       return j;
0349     }
0350
0351   } // namespace details
0352
0353   /**
0354    * Specialization of Traits dealing with non profile Root Histograms
0355    */
0356   template <typename RootHisto, unsigned int N>
0357   struct Traits<false, RootHisto, N> : details::TraitsBase<RootHisto, N> {
0358     using Histo      = RootHisto;
0359     using WeightType = double;
0360     static void fill( Histo& histo, unsigned int i, const WeightType& weight ) { histo.SetBinContent( i, weight ); }
0361   };
0362
0363   /**
0364    * Specialization of Traits dealing with profile Root Histograms
0365    */
0366   template <typename RootHisto, unsigned int N>
0367   struct Traits<true, RootHisto, N> : details::TraitsBase<details::ProfileWrapper<RootHisto>, N> {
0368     using Histo      = details::ProfileWrapper<RootHisto>;
0369     using WeightType = std::tuple<std::tuple<unsigned int, double>, double>;
0370     static constexpr void fill( Histo& histo, unsigned int i, const WeightType& weight ) {
0371       auto [c, sumWeight2]       = weight;
0372       auto [nEntries, sumWeight] = c;
0373       histo.setBinNEntries( i, nEntries );
0374       histo.SetBinContent( i, sumWeight );
0375       histo.setBinW2( i, sumWeight2 );
0376     };
0377   };
0378
0379   template <typename Traits>
0380   std::tuple<typename Traits::Histo, std::string> jsonToRootHistogram( std::string& dir, std::string& name,
0381                                                                        nlohmann::json const& j ) {
0382     return details::jsonToRootHistogramInternal<Traits>( dir, name, j, std::make_index_sequence<Traits::Dimension>() );
0383   }
0384
0385   template <typename Traits>
0386   void saveRootHisto( TFile& file, std::string dir, std::string name, nlohmann::json const& j ) {
0387     // get the Root histogram
0388     auto [histo, newDir] = jsonToRootHistogram<Traits>( dir, name, j );
0389     // Change to the proper directory in the ROOT file
0390     auto previousDir = details::changeDir( file, newDir );
0391     // write to file
0392     histo.Write();
0393     // switch back to the previous directory
0394     previousDir->cd();
0395   }
0396
0397   template <unsigned int N, bool isProfile, typename ROOTHisto>
0398   void saveRootHisto( TFile& file, std::string dir, std::string name, nlohmann::json const& j ) {
0399     saveRootHisto<Traits<isProfile, ROOTHisto, N>>( file, std::move( dir ), std::move( name ), j );
0400   }
0401
0402 } // namespace Gaudi::Histograming::Sink
0403
0404 namespace nlohmann {
0405   /// automatic translation of TAxis to json
0406   inline void to_json( nlohmann::json& j, TAxis const& axis ) {
0407     j = nlohmann::json{ { "nBins", axis.GetNbins() },
0408                         { "minValue", axis.GetXmin() },
0409                         { "maxValue", axis.GetXmax() },
0410                         { "title", axis.GetTitle() } };
0411   }
0412   /// automatic translation of Root histograms to json
0413   inline void to_json( nlohmann::json& j, TH1D const& h ) {
0414     j = Gaudi::Histograming::Sink::details::rootHistogramToJson( h );
0415   }
0416   inline void to_json( nlohmann::json& j, TH2D const& h ) {
0417     j = Gaudi::Histograming::Sink::details::rootHistogramToJson( h );
0418   }
0419   inline void to_json( nlohmann::json& j, TH3D const& h ) {
0420     j = Gaudi::Histograming::Sink::details::rootHistogramToJson( h );
0421   }
0422   inline void to_json( nlohmann::json& j, TProfile const& h ) {
0423     j = Gaudi::Histograming::Sink::details::rootHistogramToJson( h );
0424   }
0425   inline void to_json( nlohmann::json& j, TProfile2D const& h ) {
0426     j = Gaudi::Histograming::Sink::details::rootHistogramToJson( h );
0427   }
0428
0429 } // namespace nlohmann