Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-15 08:57:54 

0001 /***********************************************************************************\
0002 * (c) Copyright 1998-2025 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 <Gaudi/Accumulators/StaticHistogram.h>
0014 #include <Gaudi/MonitoringHub.h>
0015 #include <GaudiKernel/GaudiException.h>
0016 #include <TDirectory.h>
0017 #include <TFile.h>
0018 #include <TH1.h>
0019 #include <TH1D.h>
0020 #include <TH2D.h>
0021 #include <TH3D.h>
0022 #include <TProfile.h>
0023 #include <TProfile2D.h>
0024 #include <TProfile3D.h>
0025 #include <algorithm>
0026 #include <cmath>
0027 #include <functional>
0028 #include <gsl/span>
0029 #include <nlohmann/json.hpp>
0030 #include <range/v3/numeric/accumulate.hpp>
0031 #include <range/v3/range/conversion.hpp>
0032 #include <range/v3/view/split_when.hpp>
0033 #include <range/v3/view/transform.hpp>
0034 #include <string>
0035 #include <type_traits>
0036 #include <vector>
0037 
0038 #include <fmt/format.h>
0039 
0040 // upstream has renamed namespace ranges::view -> ranges::views
0041 #if RANGE_V3_VERSION < 900
0042 namespace ranges::views {
0043   using namespace ranges::view;
0044 }
0045 #endif
0046 
0047 namespace Gaudi::Histograming::Sink::detail {
0048   inline std::string formatTitle( std::string_view title, unsigned int width ) {
0049     if ( title.size() > width ) {
0050       return fmt::format( "\"{:.{}s}...\"", title, width - 3 );
0051     } else {
0052       return fmt::format( "\"{:.{}s}\"", title, width );
0053     }
0054   }
0055   inline std::string formatName( std::string_view name, unsigned int width ) {
0056     if ( name.size() > width ) {
0057       const auto preSize  = ( width / 2 ) - 2;
0058       const auto postSize = width - preSize - 1;
0059       return fmt::format( "{:.{}s}...{:.{}s}", name.substr( 0, preSize ), preSize,
0060                           name.substr( name.size() - postSize, postSize ), postSize );
0061     } else {
0062       return std::string{ name };
0063     }
0064   }
0065   constexpr std::string_view histo1DFormatting =
0066       " | {:{}.{}s} | {:{}s} | {:10} |{:11.5g} | {:<#11.5g}|{:11.5g} |{:11.5g} |";
0067   constexpr std::string_view histo2DFormatting =
0068       " ID={:{}.{}s}  {:{}s}  Ents/All={:>5.0f}/{:<5.0f}<X>/sX={:.5g}/{:<.5g},<Y>/sY={:.5g}/{:<.5g}";
0069   constexpr std::string_view histo3DFormatting =
0070       " ID={:{}.{}s}  {:{}s}  "
0071       "Ents/All={:>5.0f}/{:<5.0f}<X>/sX={:.5g}/{:<.5g},<Y>/sY={:.5g}/{:<.5g},<Z>/sZ={:.5g}/{:<.5g}";
0072 
0073   /// helper struct to print integers with fixed width
0074   struct IntWithFixedWidth {
0075     unsigned int value{ 0 };
0076   };
0077 
0078   /// sqrt or zero
0079   template <typename T>
0080   inline T sqrt_or_zero( const T x ) {
0081     return ( x > 0 ? std::sqrt( x ) : T{ 0 } );
0082   }
0083 
0084 } // namespace Gaudi::Histograming::Sink::detail
0085 
0086 /** fmt dedicated formatter for IntWithFixedWidth
0087  *
0088  *  supports 2 kinds of formatting string :
0089  *   - the empty '{}' which uses a width of 10
0090  *   - '{:n}' where n is an unsigned int, then n is the used width
0091  *  Using this formatter will result is always respecting the given width
0092  *  If the number of digits is less than n, the integer is printed using %d,
0093  *  otherwise it's using scientific notation with n-5 digits. So n has to be
0094  *  at least 6
0095  */
0096 template <>
0097 struct fmt::formatter<Gaudi::Histograming::Sink::detail::IntWithFixedWidth> {
0098   unsigned int width{ 10 };
0099 
0100   constexpr auto parse( format_parse_context& ctx ) {
0101     auto it = ctx.begin();
0102     if ( it == ctx.end() || *it == '}' ) { return it; }
0103     int w = fmt::detail::parse_nonnegative_int( it, ctx.end(), -1 );
0104     if ( w == -1 ) throw fmt::format_error( "bad" );
0105     width = (unsigned int)w;
0106     if ( it == ctx.end() || *it != '}' ) throw fmt::format_error( "bad" );
0107     return it;
0108   }
0109 
0110   auto format( Gaudi::Histograming::Sink::detail::IntWithFixedWidth n, format_context& ctx ) const {
0111     auto s = fmt::format( "{:{}d}", n.value, width );
0112     if ( s.size() > width ) { s = fmt::format( "{:{}.{}g}", (double)n.value, width, width - 5 ); }
0113     return fmt::format_to( ctx.out(), "{:s}", s );
0114   }
0115 };
0116 
0117 namespace Gaudi::Histograming::Sink {
0118 
0119   /**
0120    * templated Traits dealing with Root Histogram filling for standard histograms
0121    *
0122    * Specializatoins of the Traits type must provide the following static functions :
0123    *   - Histo create( std::string& name, std::string& title, Axis& axis... )
0124    *       it should intanciate a new ROOT histogram instance and return it
0125    *   - void fillMetaData( Histo& histo, nlohmann::json const& jsonAxis, unsigned int nentries)
0126    *       it should fill metadata in the ROOT histo histogram provides from the list of axis number of entries
0127    *   - void fill( Histo& histo, unsigned int nbins, const WeightType& weight )
0128    *       it should fill the given number of bin with the given values
0129    * Last point, it should have a static constexpr unsigned int Dimension
0130    */
0131   template <bool isProfile, typename RootHisto, unsigned int N>
0132   struct Traits;
0133 
0134   /**
0135    * generic function to convert json to a ROOT Histogram
0136    *
0137    * returns the Root histogram and the dir where to save it in the Root file
0138    * This may be different from input dir in case name has slashes
0139    */
0140   template <typename Traits>
0141   std::tuple<typename Traits::Histo, std::string> jsonToRootHistogram( std::string& dir, std::string& name,
0142                                                                        nlohmann::json const& j );
0143 
0144   /**
0145    * generic function to convert a ROOT Histogram to json
0146    *
0147    * essentially used for backward compatibility of old HistogramService with MonitoringHub
0148    */
0149   template <typename Histo>
0150   nlohmann::json rootHistogramTojson( Histo const& );
0151 
0152   /**
0153    * generic method to save histograms to files, based on Traits
0154    */
0155   template <typename Traits>
0156   void saveRootHisto( TFile& file, std::string dir, std::string name, nlohmann::json const& j );
0157 
0158   /**
0159    * generic method to save regular histograms to files
0160    *
0161    * Can be used in most cases as the handler function to register into Sink::Base
0162    * contains all the boiler plate code and redirects specific code to the adapted Traits template
0163    */
0164   template <unsigned int N, bool isProfile, typename ROOTHisto>
0165   void saveRootHisto( TFile& file, std::string dir, std::string name, nlohmann::json const& j );
0166 
0167   namespace details {
0168 
0169     /// Small helper struct representing the Axis of an Histogram
0170     struct Axis {
0171       unsigned int nBins{};
0172       double       minValue{};
0173       double       maxValue{};
0174       std::string  title;
0175     };
0176 
0177     /// extract an Axis from json data
0178     inline Axis jsonToAxis( nlohmann::json& jAxis ) {
0179       return { jAxis.at( "nBins" ).get<unsigned int>(), jAxis.at( "minValue" ).get<double>(),
0180                jAxis.at( "maxValue" ).get<double>(),
0181                ";" + jAxis.at( "title" ).get<std::string>() }; // ";" to prepare concatenations of titles
0182     }
0183 
0184     /**
0185      * handles cases where name includes '/' character(s) and move needed part of it
0186      * to dir. Also handle case of absolute names
0187      */
0188     inline void fixNameAndDir( std::string& name, std::string& dir ) {
0189       if ( name[0] == '/' ) {
0190         dir  = "";
0191         name = name.substr( 1 );
0192       }
0193       // take into account the case where name contains '/'s (e.g. "Group/Name") by
0194       // moving the prefix into dir
0195       if ( auto pos = name.rfind( '/' ); pos != std::string::npos ) {
0196         dir += '/' + name.substr( 0, pos );
0197         name = name.substr( pos + 1 );
0198       }
0199     }
0200 
0201     /**
0202      * generic function to convert json to a ROOT Histogram - internal implementation
0203      */
0204     template <typename Traits, std::size_t... index>
0205     std::tuple<typename Traits::Histo, std::string> jsonToRootHistogramInternal( std::string& dir, std::string& name,
0206                                                                                  nlohmann::json const& j,
0207                                                                                  std::index_sequence<index...> ) {
0208       // extract data from json
0209       auto jsonAxis = j.at( "axis" );
0210       auto axis     = std::array{ jsonToAxis( jsonAxis[index] )... };
0211       auto weights  = j.at( "bins" ).get<typename Traits::WeightType>();
0212       auto title    = j.at( "title" ).get<std::string>();
0213       auto nentries = j.at( "nEntries" ).get<unsigned int>();
0214       // weird way ROOT has to give titles to axis
0215       title += ( axis[index].title + ... );
0216       // compute total number of bins, multiplying bins per axis
0217       auto totNBins = ( ( axis[index].nBins + 2 ) * ... );
0218       // fix name and dir if needed
0219       fixNameAndDir( name, dir );
0220       // Create Root histogram calling constructors with the args tuple
0221       auto histo = Traits::create( name, title, axis[index]... );
0222       // fill Histo
0223       Traits::fill( histo, totNBins, weights );
0224       // in case we have sums, overwrite them in the histogram with our more precise values
0225       // FIXME This is only supporting regular histograms. It won't work in case of weighted histograms
0226       if constexpr ( sizeof...( index ) == 1 ) {
0227         if ( j.find( "sum" ) != j.end() ) {
0228           double s[13];
0229           s[0] = j.at( "nTotEntries" ).get<double>();
0230           s[1] = j.at( "nTotEntries" ).get<double>();
0231           s[2] = j.at( "sum" ).get<double>();
0232           s[3] = j.at( "sum2" ).get<double>();
0233           histo.PutStats( s );
0234         }
0235       } else {
0236         if ( j.find( "sumx" ) != j.end() ) {
0237           double s[13];
0238           s[0] = j.at( "nTotEntries" ).get<double>();
0239           s[1] = j.at( "nTotEntries" ).get<double>();
0240           s[2] = j.at( "sumx" ).get<double>();
0241           s[3] = j.at( "sumx2" ).get<double>();
0242           s[4] = j.at( "sumy" ).get<double>();
0243           s[5] = j.at( "sumy2" ).get<double>();
0244           s[6] = j.at( "sumxy" ).get<double>();
0245           if constexpr ( sizeof...( index ) > 2 ) {
0246             s[7]  = j.at( "sumz" ).get<double>();
0247             s[8]  = j.at( "sumz2" ).get<double>();
0248             s[9]  = j.at( "sumxz" ).get<double>();
0249             s[10] = j.at( "sumyz" ).get<double>();
0250           }
0251           histo.PutStats( s );
0252         }
0253       }
0254       // fill histo metadata, e.g. bins and number of entries
0255       Traits::fillMetaData( histo, jsonAxis, nentries );
0256       return { histo, dir };
0257     }
0258 
0259     /**
0260      * Common base for Traits dealing with Histogram conversions to Root
0261      * Provides generic implementation for creating the histogram and filling meta data
0262      * The filling (method fill) is not implemented
0263      */
0264     template <typename RootHisto, unsigned int N>
0265     struct TraitsBase {
0266       static constexpr unsigned int Dimension{ N };
0267       // hleper method for the creation of the Root histogram, using an index_sequence for handling the axis
0268       template <typename AxisArray, std::size_t... index>
0269       static RootHisto create( std::string& name, std::string& title, AxisArray axis, std::index_sequence<index...> ) {
0270         return std::make_from_tuple<RootHisto>(
0271             std::tuple_cat( std::tuple{ name.c_str(), title.c_str() },
0272                             std::tuple{ std::get<index>( axis ).nBins, std::get<index>( axis ).minValue,
0273                                         std::get<index>( axis ).maxValue }... ) );
0274       }
0275       // Generic creation of the Root histogram from name, title and a set of axis
0276       template <typename... Axis>
0277       static RootHisto create( std::string& name, std::string& title, Axis&... axis ) {
0278         return create( name, title, std::make_tuple( axis... ), std::make_index_sequence<sizeof...( Axis )>() );
0279       }
0280       // Generic fill of the metadata of the Root histogram from json data
0281       static void fillMetaData( RootHisto& histo, nlohmann::json const& jsonAxis, unsigned int nentries ) {
0282         auto try_set_bin_labels = [&histo, &jsonAxis]( auto idx ) {
0283           TAxis* axis = nullptr;
0284           switch ( idx ) {
0285           case 0:
0286             axis = histo.GetXaxis();
0287             break;
0288           case 1:
0289             axis = histo.GetYaxis();
0290             break;
0291           case 2:
0292             axis = histo.GetZaxis();
0293             break;
0294           default:
0295             break;
0296           }
0297           if ( axis && jsonAxis[idx].contains( "labels" ) ) {
0298             const auto labels = jsonAxis[idx].at( "labels" );
0299             for ( unsigned int i = 0; i < labels.size(); i++ ) {
0300               axis->SetBinLabel( i + 1, labels[i].template get<std::string>().c_str() );
0301             }
0302           }
0303           if ( axis && jsonAxis[idx].contains( "xbins" ) ) {
0304             const auto xbins = jsonAxis[idx].at( "xbins" );
0305             axis->Set( xbins.size() - 1, xbins.template get<std::vector<double>>().data() );
0306           }
0307         };
0308         for ( unsigned int i = 0; i < jsonAxis.size(); i++ ) try_set_bin_labels( i );
0309         // Fix the number of entries, wrongly computed by ROOT from the numer of calls to fill
0310         histo.SetEntries( nentries );
0311       }
0312     };
0313 
0314     /// helper Wrapper around TProfileX to be able to fill it
0315     template <typename TP>
0316     struct ProfileWrapper : TP {
0317       template <typename... Args>
0318       ProfileWrapper( Args&&... args ) : TP( std::forward<Args>( args )... ) {
0319         // When the static function TH1::SetDefaultSumw2 had been called (passing true), `Sumw2(true)` is automatically
0320         // called by the constructor of TProfile. This is a problem because in the profiles in Gaudi::Accumulators we do
0321         // not keep track of the sum of squares of weights and consequently don't set fBinSumw2 of the TProfile, which
0322         // in turn leads to a self-inconsistent TProfile object. The "fix" is to disable sum of squares explicitly.
0323         this->Sumw2( false );
0324       }
0325       ProfileWrapper( ProfileWrapper const& )             = delete;
0326       ProfileWrapper& operator=( ProfileWrapper const& )  = delete;
0327       ProfileWrapper( ProfileWrapper const&& )            = delete;
0328       ProfileWrapper& operator=( ProfileWrapper const&& ) = delete;
0329       void            setBinNEntries( Int_t i, Int_t n ) { this->fBinEntries.fArray[i] = n; }
0330       void            setBinW2( Int_t i, Double_t v ) { this->fSumw2.fArray[i] = v; }
0331       Double_t        getBinW2( Int_t i ) { return this->fSumw2.fArray[i]; }
0332     };
0333 
0334     /**
0335      * changes to the ROOT directory given in the current ROOT file and returns
0336      * the current directory before the change
0337      */
0338     inline TDirectory* changeDir( TFile& file, const std::string& dir ) {
0339       // remember the current directory
0340       auto previousDir = gDirectory;
0341       // find or create the directory for the histogram
0342       using namespace ranges;
0343       auto is_delimiter        = []( auto c ) { return c == '/' || c == '.'; };
0344       auto transform_to_string = views::transform( []( auto&& rng ) { return rng | to<std::string>; } );
0345       auto currentDir          = accumulate( dir | views::split_when( is_delimiter ) | transform_to_string,
0346                                              file.GetDirectory( "" ), []( auto current, auto&& dir_level ) {
0347                                       if ( current ) {
0348                                         // try to get next level
0349                                         auto nextDir = current->GetDirectory( dir_level.c_str() );
0350                                         // if it does not exist, create it
0351                                         if ( !nextDir ) nextDir = current->mkdir( dir_level.c_str() );
0352                                         // move to next level
0353                                         current = nextDir;
0354                                       }
0355                                       return current;
0356                                     } );
0357       if ( !currentDir )
0358         throw GaudiException( "Could not create directory " + dir, "Histogram::Sink::Root", StatusCode::FAILURE );
0359       // switch to the directory
0360       currentDir->cd();
0361       return previousDir;
0362     }
0363 
0364     template <typename Histo>
0365     nlohmann::json allAxisTojson( Histo const& h ) {
0366       if constexpr ( std::is_base_of_v<TH3D, Histo> ) {
0367         return { *h.GetXaxis(), *h.GetYaxis(), *h.GetZaxis() };
0368       } else if constexpr ( std::is_base_of_v<TProfile2D, Histo> || std::is_base_of_v<TH2D, Histo> ) {
0369         return { *h.GetXaxis(), *h.GetYaxis() };
0370       } else {
0371         return { *h.GetXaxis() };
0372       }
0373     }
0374 
0375     template <typename Histo>
0376     nlohmann::json binsTojson( Histo const& h ) {
0377       if constexpr ( std::is_base_of_v<TProfile, Histo> || std::is_base_of_v<TProfile2D, Histo> ) {
0378         nlohmann::json j;
0379         // ROOT TProfile interface being completely inconsistent, we have to play
0380         // with different ways to access values of the histogram... one per value !
0381         auto* sums  = h.GetArray();
0382         auto* sums2 = h.GetSumw2();
0383         for ( unsigned long n = 0; n < (unsigned long)h.GetSize(); n++ ) {
0384           j.push_back( nlohmann::json{ { h.GetBinEntries( n ), sums[n] }, sums2->At( n ) } );
0385         }
0386         return j;
0387       } else {
0388         return gsl::span{ h.GetArray(), (unsigned long)h.GetSize() };
0389       }
0390     }
0391 
0392     /// automatic translation of Root Histograms to json
0393     template <typename Histo>
0394     nlohmann::json rootHistogramToJson( Histo const& h ) {
0395       bool        isProfile = std::is_base_of_v<TProfile, Histo> || std::is_base_of_v<TProfile2D, Histo>;
0396       std::string type      = isProfile ? "histogram:ProfileHistogram:double" : "histogram:Histogram:double";
0397       auto        j         = nlohmann::json{ { "type", type },
0398                                               { "title", h.GetTitle() },
0399                                               { "dimension", h.GetDimension() },
0400                                               { "empty", (int)h.GetEntries() == 0 },
0401                                               { "nEntries", (int)h.GetEntries() },
0402                                               { "axis", allAxisTojson( h ) },
0403                                               { "bins", binsTojson( h ) } };
0404       if ( !isProfile ) {
0405         double s[13];
0406         h.GetStats( s );
0407         if ( h.GetDimension() == 1 ) {
0408           j["nTotEntries"] = s[0];
0409           j["sum"]         = s[2];
0410           j["sum2"]        = s[3];
0411           j["mean"]        = s[2] / s[0];
0412         } else {
0413           j["nTotEntries"] = s[0];
0414           j["sumx"]        = s[2];
0415           j["sumx2"]       = s[3];
0416           j["meanx"]       = s[2] / s[0];
0417           j["sumy"]        = s[4];
0418           j["sumy2"]       = s[5];
0419           j["sumxy"]       = s[6];
0420           j["meany"]       = s[4] / s[0];
0421           if ( h.GetDimension() >= 3 ) {
0422             j["sumz"]  = s[7];
0423             j["sumz2"] = s[8];
0424             j["sumxz"] = s[9];
0425             j["sumyz"] = s[10];
0426             j["meanz"] = s[7] / s[0];
0427           }
0428         }
0429       }
0430       return j;
0431     }
0432 
0433   } // namespace details
0434 
0435   /**
0436    * Specialization of Traits dealing with non profile Root Histograms
0437    */
0438   template <typename RootHisto, unsigned int N>
0439   struct Traits<false, RootHisto, N> : details::TraitsBase<RootHisto, N> {
0440     using Histo      = RootHisto;
0441     using WeightType = std::vector<double>;
0442     static void fill( Histo& histo, unsigned int nbins, const WeightType& weight ) {
0443       for ( unsigned int i = 0; i < nbins; i++ ) { histo.SetBinContent( i, weight[i] ); }
0444     }
0445   };
0446 
0447   /**
0448    * Specialization of Traits dealing with profile Root Histograms
0449    */
0450   template <typename RootHisto, unsigned int N>
0451   struct Traits<true, RootHisto, N> : details::TraitsBase<details::ProfileWrapper<RootHisto>, N> {
0452     using Histo      = details::ProfileWrapper<RootHisto>;
0453     using WeightType = std::vector<std::tuple<std::tuple<unsigned int, double>, double>>;
0454     static constexpr void fill( Histo& histo, unsigned int nbins, const WeightType& weight ) {
0455       for ( unsigned int i = 0; i < nbins; i++ ) {
0456         auto [c, sumWeight2]       = weight[i];
0457         auto [nEntries, sumWeight] = c;
0458         histo.setBinNEntries( i, nEntries );
0459         histo.SetBinContent( i, sumWeight );
0460         histo.setBinW2( i, sumWeight2 );
0461       }
0462     }
0463   };
0464 
0465   template <typename Traits>
0466   std::tuple<typename Traits::Histo, std::string> jsonToRootHistogram( std::string& dir, std::string& name,
0467                                                                        nlohmann::json const& j ) {
0468     return details::jsonToRootHistogramInternal<Traits>( dir, name, j, std::make_index_sequence<Traits::Dimension>() );
0469   }
0470 
0471   template <typename Traits>
0472   void saveRootHisto( TFile& file, std::string dir, std::string name, nlohmann::json const& j ) {
0473     // get the Root histogram
0474     auto [histo, newDir] = jsonToRootHistogram<Traits>( dir, name, j );
0475     // Change to the proper directory in the ROOT file
0476     auto previousDir = details::changeDir( file, newDir );
0477     // write to file
0478     histo.Write();
0479     // switch back to the previous directory
0480     previousDir->cd();
0481   }
0482 
0483   template <unsigned int N, bool isProfile, typename ROOTHisto>
0484   void saveRootHisto( TFile& file, std::string dir, std::string name, nlohmann::json const& j ) {
0485     saveRootHisto<Traits<isProfile, ROOTHisto, N>>( file, std::move( dir ), std::move( name ), j );
0486   }
0487 
0488   /// Direct conversion of 1D histograms from Gaudi to ROOT format
0489   template <Accumulators::atomicity Atomicity = Accumulators::atomicity::full, typename Arithmetic = double>
0490   details::ProfileWrapper<TProfile> profileHisto1DToRoot( std::string name, Monitoring::Hub::Entity const& ent ) {
0491     // get original histogram from the Entity
0492     auto const& gaudiHisto =
0493         *reinterpret_cast<Gaudi::Accumulators::StaticProfileHistogram<1, Atomicity, Arithmetic>*>( ent.id() );
0494     // convert to Root
0495     auto const&                       gaudiAxis = gaudiHisto.template axis<0>();
0496     details::ProfileWrapper<TProfile> histo{ name.c_str(), gaudiHisto.title().c_str(), gaudiAxis.numBins(),
0497                                              gaudiAxis.minValue(), gaudiAxis.maxValue() };
0498     if ( !gaudiAxis.labels().empty() ) {
0499       auto& axis = *histo.GetXaxis();
0500       for ( unsigned int i = 0; i < gaudiAxis.labels().size(); i++ ) {
0501         axis.SetBinLabel( i + 1, gaudiAxis.labels()[i].c_str() );
0502       }
0503     }
0504     for ( unsigned int i = 0; i < gaudiHisto.totNBins(); i++ ) {
0505       auto const& [tmp, sumw2] = gaudiHisto.binValue( i );
0506       auto const& [nent, sumw] = tmp;
0507       histo.SetBinContent( i, sumw );
0508       histo.setBinNEntries( i, nent );
0509       histo.setBinW2( i, sumw2 );
0510     }
0511     unsigned long totNEntries{ 0 };
0512     for ( unsigned int i = 0; i < gaudiHisto.totNBins(); i++ ) { totNEntries += gaudiHisto.nEntries( i ); }
0513     histo.SetEntries( totNEntries );
0514     return histo;
0515   }
0516 
0517   /// Direct conversion of 2D histograms from Gaudi to ROOT format
0518   template <Accumulators::atomicity Atomicity = Accumulators::atomicity::full, typename Arithmetic = double>
0519   details::ProfileWrapper<TProfile2D> profileHisto2DToRoot( std::string name, Monitoring::Hub::Entity const& ent ) {
0520     // get original histogram from the Entity
0521     auto const& gaudiHisto =
0522         *reinterpret_cast<Gaudi::Accumulators::StaticProfileHistogram<2, Atomicity, Arithmetic>*>( ent.id() );
0523     // convert to Root
0524     auto const&                         gaudiXAxis = gaudiHisto.template axis<0>();
0525     auto const&                         gaudiYAxis = gaudiHisto.template axis<1>();
0526     details::ProfileWrapper<TProfile2D> histo{ name.c_str(),          gaudiHisto.title().c_str(), gaudiXAxis.numBins(),
0527                                                gaudiXAxis.minValue(), gaudiXAxis.maxValue(),      gaudiYAxis.numBins(),
0528                                                gaudiYAxis.minValue(), gaudiYAxis.maxValue() };
0529     if ( !gaudiXAxis.labels().empty() ) {
0530       auto& axis = *histo.GetXaxis();
0531       for ( unsigned int i = 0; i < gaudiXAxis.labels().size(); i++ ) {
0532         axis.SetBinLabel( i + 1, gaudiXAxis.labels()[i].c_str() );
0533       }
0534     }
0535     if ( !gaudiYAxis.labels().empty() ) {
0536       auto& axis = *histo.GetYaxis();
0537       for ( unsigned int i = 0; i < gaudiYAxis.labels().size(); i++ ) {
0538         axis.SetBinLabel( i + 1, gaudiYAxis.labels()[i].c_str() );
0539       }
0540     }
0541     for ( unsigned int i = 0; i < gaudiHisto.totNBins(); i++ ) {
0542       auto const& [tmp, sumw2] = gaudiHisto.binValue( i );
0543       auto const& [nent, sumw] = tmp;
0544       histo.SetBinContent( i, sumw );
0545       histo.setBinNEntries( i, nent );
0546       histo.setBinW2( i, sumw2 );
0547     }
0548     unsigned long totNEntries{ 0 };
0549     for ( unsigned int i = 0; i < gaudiHisto.totNBins(); i++ ) { totNEntries += gaudiHisto.nEntries( i ); }
0550     histo.SetEntries( totNEntries );
0551     return histo;
0552   }
0553 
0554   /// Direct conversion of 3D histograms from Gaudi to ROOT format
0555   template <Accumulators::atomicity Atomicity = Accumulators::atomicity::full, typename Arithmetic = double>
0556   details::ProfileWrapper<TProfile3D> profileHisto3DToRoot( std::string name, Monitoring::Hub::Entity const& ent ) {
0557     // get original histogram from the Entity
0558     auto const& gaudiHisto =
0559         *reinterpret_cast<Gaudi::Accumulators::StaticProfileHistogram<3, Atomicity, Arithmetic>*>( ent.id() );
0560     // convert to Root
0561     auto const&                         gaudiXAxis = gaudiHisto.template axis<0>();
0562     auto const&                         gaudiYAxis = gaudiHisto.template axis<1>();
0563     auto const&                         gaudiZAxis = gaudiHisto.template axis<2>();
0564     details::ProfileWrapper<TProfile3D> histo{ name.c_str(),          gaudiHisto.title().c_str(), gaudiXAxis.numBins(),
0565                                                gaudiXAxis.minValue(), gaudiXAxis.maxValue(),      gaudiYAxis.numBins(),
0566                                                gaudiYAxis.minValue(), gaudiYAxis.maxValue(),      gaudiZAxis.numBins(),
0567                                                gaudiZAxis.minValue(), gaudiZAxis.maxValue() };
0568     if ( !gaudiXAxis.labels().empty() ) {
0569       auto& axis = *histo.GetXaxis();
0570       for ( unsigned int i = 0; i < gaudiXAxis.labels().size(); i++ ) {
0571         axis.SetBinLabel( i + 1, gaudiXAxis.labels()[i].c_str() );
0572       }
0573     }
0574     if ( !gaudiYAxis.labels().empty() ) {
0575       auto& axis = *histo.GetYaxis();
0576       for ( unsigned int i = 0; i < gaudiYAxis.labels().size(); i++ ) {
0577         axis.SetBinLabel( i + 1, gaudiYAxis.labels()[i].c_str() );
0578       }
0579     }
0580     if ( !gaudiZAxis.labels().empty() ) {
0581       auto& axis = *histo.GetZaxis();
0582       for ( unsigned int i = 0; i < gaudiZAxis.labels().size(); i++ ) {
0583         axis.SetBinLabel( i + 1, gaudiZAxis.labels()[i].c_str() );
0584       }
0585     }
0586     for ( unsigned int i = 0; i < gaudiHisto.totNBins(); i++ ) {
0587       auto const& [tmp, sumw2] = gaudiHisto.binValue( i );
0588       auto const& [nent, sumw] = tmp;
0589       histo.SetBinContent( i, sumw );
0590       histo.setBinNEntries( i, nent );
0591       histo.setBinW2( i, sumw2 );
0592     }
0593     unsigned long totNEntries{ 0 };
0594     for ( unsigned int i = 0; i < gaudiHisto.totNBins(); i++ ) { totNEntries += gaudiHisto.nEntries( i ); }
0595     histo.SetEntries( totNEntries );
0596     return histo;
0597   }
0598 
0599   template <unsigned int N, Accumulators::atomicity Atomicity = Accumulators::atomicity::full,
0600             typename Arithmetic = double>
0601   void saveProfileHisto( TFile& file, std::string dir, std::string name, Monitoring::Hub::Entity const& ent ) {
0602     // fix name and dir if needed
0603     details::fixNameAndDir( name, dir );
0604     // Change to the proper directory in the ROOT file
0605     auto previousDir = details::changeDir( file, dir );
0606     // write to file
0607     if constexpr ( N == 1 ) {
0608       profileHisto1DToRoot<Atomicity, Arithmetic>( name, ent ).Write();
0609     } else if constexpr ( N == 2 ) {
0610       profileHisto2DToRoot<Atomicity, Arithmetic>( name, ent ).Write();
0611     } else if constexpr ( N == 3 ) {
0612       profileHisto3DToRoot<Atomicity, Arithmetic>( name, ent ).Write();
0613     }
0614     // switch back to the previous directory
0615     previousDir->cd();
0616   }
0617 
0618   namespace details {
0619     template <typename TYPE>
0620     struct BinAvValue {
0621     private:
0622       TYPE m_min{ 0 }, m_max{ 0 }, m_binSize{ 0 };
0623       bool m_isInt{ false };
0624 
0625     public:
0626       BinAvValue( const TYPE minV, const TYPE maxV, unsigned int nbins, const bool isInt )
0627           : m_min( minV ), m_max( maxV ), m_isInt( isInt ) {
0628         if ( nbins > 0 ) { m_binSize = ( maxV - minV ) / nbins; }
0629       }
0630       auto operator()( const unsigned int n ) const {
0631         TYPE binAv{ 0 };
0632         if ( m_isInt ) {
0633           // Need to find all int values that are in range for this bin and form average
0634           const auto   binMin = m_min + ( n - 1.0 ) * m_binSize;
0635           const auto   binMax = binMin + m_binSize;
0636           unsigned int nAv{ 0 };
0637           for ( auto iv = std::floor( binMin ); iv < std::ceil( binMax ); iv += 1.0 ) {
0638             if ( binMin <= iv && iv <= binMax ) {
0639               binAv += iv;
0640               ++nAv;
0641             }
0642           }
0643           binAv = ( nAv > 0 ? binAv / nAv : 0 );
0644         } else {
0645           // For floats just use bin centre
0646           binAv = m_min + ( n - 0.5 ) * m_binSize;
0647         }
0648         return binAv;
0649       }
0650     };
0651   } // namespace details
0652 
0653   template <Gaudi::Accumulators::atomicity Atomicity = Gaudi::Accumulators::atomicity::full,
0654             typename Arithmetic                      = double>
0655   std::string printProfileHisto1D( std::string_view name, Gaudi::Monitoring::Hub::Entity const& ent,
0656                                    unsigned int stringsWidth = 45 ) {
0657     // Type to use for internal calculation. Use double to avoid precision issues across different builds.
0658     using FPType = double;
0659     // get original histogram from the Entity
0660     auto const& gaudiHisto =
0661         *reinterpret_cast<Gaudi::Accumulators::StaticProfileHistogram<1, Atomicity, Arithmetic>*>( ent.id() );
0662     // compute min and max values, we actually display the axis limits
0663     auto const&        gaudiAxisX = gaudiHisto.template axis<0>();
0664     const FPType       minValueX  = gaudiAxisX.minValue();
0665     const FPType       maxValueX  = gaudiAxisX.maxValue();
0666     const unsigned int nBinsX     = gaudiAxisX.numBins();
0667     // Compute fist and second momenta for normal dimenstion
0668     // Compute 1st and 2nd momenta for the profile dimension
0669     unsigned int              nEntries{ 0 }, nAllEntries{ 0 };
0670     FPType                    totalSumW{ 0 };
0671     FPType                    sumX{ 0 }, sum2X{ 0 }, sum3X{ 0 }, sum4X{ 0 };
0672     const details::BinAvValue xBinAv( minValueX, maxValueX, nBinsX, std::is_integral_v<Arithmetic> );
0673     for ( unsigned int nx = 0; nx <= nBinsX + 1; ++nx ) {
0674       auto const& [tmp, sumw2] = gaudiHisto.binValue( nx );
0675       auto const& [nent, sumw] = tmp;
0676       nAllEntries += nent;
0677       if ( nx > 0 && nx <= nBinsX ) {
0678         const FPType binXValue = xBinAv( nx );
0679         nEntries += nent;
0680         totalSumW += sumw;
0681         auto val = binXValue * sumw;
0682         sumX += val;
0683         val *= binXValue;
0684         sum2X += val;
0685         val *= binXValue;
0686         sum3X += val;
0687         val *= binXValue;
0688         sum4X += val;
0689       }
0690     }
0691     const std::string ftitle = detail::formatTitle( gaudiHisto.title(), stringsWidth - 2 );
0692     if ( nEntries == 0 || !( fabs( totalSumW ) > 0 ) ) { return ""; }
0693     const FPType meanX     = sumX / totalSumW;
0694     const FPType sigmaX2   = ( sum2X / totalSumW ) - std::pow( meanX, 2 );
0695     const FPType stddevX   = detail::sqrt_or_zero( sigmaX2 );
0696     const FPType EX3       = sum3X / totalSumW;
0697     const FPType A         = sigmaX2 * stddevX;
0698     const FPType skewnessX = ( fabs( A ) > 0.0 ? ( EX3 - ( 3 * sigmaX2 + meanX * meanX ) * meanX ) / A : 0.0 );
0699     const FPType B         = sigmaX2 * sigmaX2;
0700     const FPType kurtosisX =
0701         ( fabs( B ) > 0.0 && fabs( A ) > 0.0
0702               ? ( sum4X / totalSumW - meanX * ( 4 * EX3 - meanX * ( 6 * sigmaX2 + 3 * meanX * meanX ) ) ) / B
0703               : 3.0 );
0704     // print
0705     return fmt::format( fmt::runtime( detail::histo1DFormatting ), detail::formatName( name, stringsWidth ),
0706                         stringsWidth, stringsWidth, ftitle, stringsWidth, detail::IntWithFixedWidth{ nEntries }, meanX,
0707                         stddevX, skewnessX, kurtosisX - FPType{ 3.0 } );
0708   }
0709 
0710   template <Gaudi::Accumulators::atomicity Atomicity = Gaudi::Accumulators::atomicity::full,
0711             typename Arithmetic                      = double>
0712   std::string printProfileHisto2D( std::string_view name, Gaudi::Monitoring::Hub::Entity const& ent,
0713                                    unsigned int stringsWidth = 45 ) {
0714     // Type to use for internal calculation. Use double to avoid precision issues across different builds.
0715     using FPType = double;
0716     // get original histogram from the Entity
0717     auto const& gaudiHisto =
0718         *reinterpret_cast<Gaudi::Accumulators::StaticProfileHistogram<2, Atomicity, Arithmetic>*>( ent.id() );
0719     // compute min and max values, we actually display the axis limits
0720     auto const&        gaudiAxisX = gaudiHisto.template axis<0>();
0721     const FPType       minValueX  = gaudiAxisX.minValue();
0722     const FPType       maxValueX  = gaudiAxisX.maxValue();
0723     const unsigned int nBinsX     = gaudiAxisX.numBins();
0724     auto const&        gaudiAxisY = gaudiHisto.template axis<1>();
0725     const FPType       minValueY  = gaudiAxisY.minValue();
0726     const FPType       maxValueY  = gaudiAxisY.maxValue();
0727     const unsigned int nBinsY     = gaudiAxisY.numBins();
0728     // Compute fist and second momenta for normal dimenstion
0729     // Compute 1st and 2nd momenta for the profile dimension
0730     FPType                    nEntries{ 0 }, nAllEntries{ 0 };
0731     FPType                    totalSumW{ 0 };
0732     FPType                    sumX{ 0 }, sumY{ 0 };
0733     FPType                    sum2X{ 0 }, sum2Y{ 0 };
0734     const details::BinAvValue xBinAv( minValueX, maxValueX, nBinsX, std::is_integral_v<Arithmetic> );
0735     const details::BinAvValue yBinAv( minValueY, maxValueY, nBinsY, std::is_integral_v<Arithmetic> );
0736     for ( unsigned int ny = 0; ny <= nBinsY + 1; ++ny ) {
0737       const auto offset = ny * ( nBinsX + 2 );
0738       for ( unsigned int nx = 0; nx <= nBinsX + 1; ++nx ) {
0739         auto const& [tmp, sumw2] = gaudiHisto.binValue( nx + offset );
0740         auto const& [nent, sumw] = tmp;
0741         nAllEntries += nent;
0742         if ( nx > 0 && ny > 0 && nx <= nBinsX && ny <= nBinsY ) {
0743           const FPType binXValue = xBinAv( nx );
0744           const FPType binYValue = yBinAv( ny );
0745           nEntries += nent;
0746           totalSumW += sumw;
0747           sumX += binXValue * sumw;
0748           sum2X += binXValue * binXValue * sumw;
0749           sumY += binYValue * sumw;
0750           sum2Y += binYValue * binYValue * sumw;
0751         }
0752       }
0753     }
0754     const FPType meanX = fabs( totalSumW ) > 0 ? sumX / totalSumW : 0.0;
0755     const FPType stddevX =
0756         fabs( totalSumW ) > 0 ? detail::sqrt_or_zero( ( sum2X - sumX * ( sumX / totalSumW ) ) / totalSumW ) : 0.0;
0757     const FPType meanY = fabs( totalSumW ) > 0 ? sumY / totalSumW : 0.0;
0758     const FPType stddevY =
0759         fabs( totalSumW ) > 0 ? detail::sqrt_or_zero( ( sum2Y - sumY * ( sumY / totalSumW ) ) / totalSumW ) : 0.0;
0760     // print
0761     const std::string ftitle = detail::formatTitle( gaudiHisto.title(), stringsWidth - 2 );
0762     return fmt::format( fmt::runtime( detail::histo2DFormatting ), detail::formatName( name, stringsWidth ),
0763                         stringsWidth, stringsWidth, ftitle, stringsWidth, nEntries, nAllEntries, meanX, stddevX, meanY,
0764                         stddevY );
0765   }
0766 
0767   template <Gaudi::Accumulators::atomicity Atomicity = Gaudi::Accumulators::atomicity::full,
0768             typename Arithmetic                      = double>
0769   std::string printProfileHisto3D( std::string_view name, Gaudi::Monitoring::Hub::Entity const& ent,
0770                                    unsigned int stringsWidth = 45 ) {
0771     // Type to use for internal calculation. Use double to avoid precision issues across different builds.
0772     using FPType = double;
0773     // get original histogram from the Entity
0774     auto const& gaudiHisto =
0775         *reinterpret_cast<Gaudi::Accumulators::StaticProfileHistogram<3, Atomicity, Arithmetic>*>( ent.id() );
0776     // compute min and max values, we actually display the axis limits
0777     auto const&        gaudiAxisX = gaudiHisto.template axis<0>();
0778     const FPType       minValueX  = gaudiAxisX.minValue();
0779     const FPType       maxValueX  = gaudiAxisX.maxValue();
0780     const unsigned int nBinsX     = gaudiAxisX.numBins();
0781     auto const&        gaudiAxisY = gaudiHisto.template axis<1>();
0782     const FPType       minValueY  = gaudiAxisY.minValue();
0783     const FPType       maxValueY  = gaudiAxisY.maxValue();
0784     const unsigned int nBinsY     = gaudiAxisY.numBins();
0785     auto const&        gaudiAxisZ = gaudiHisto.template axis<2>();
0786     const FPType       minValueZ  = gaudiAxisZ.minValue();
0787     const FPType       maxValueZ  = gaudiAxisZ.maxValue();
0788     const unsigned int nBinsZ     = gaudiAxisZ.numBins();
0789     // Compute fist and second momenta for normal dimenstion
0790     // Compute 1st and 2nd momenta for the profile dimension
0791     FPType                    nEntries{ 0 }, nAllEntries{ 0 };
0792     FPType                    totalSumW{};
0793     FPType                    sumX{ 0 }, sumY{ 0 }, sumZ{ 0 };
0794     FPType                    sum2X{ 0 }, sum2Y{ 0 }, sum2Z{ 0 };
0795     const details::BinAvValue xBinAv( minValueX, maxValueX, nBinsX, std::is_integral_v<Arithmetic> );
0796     const details::BinAvValue yBinAv( minValueY, maxValueY, nBinsY, std::is_integral_v<Arithmetic> );
0797     const details::BinAvValue zBinAv( minValueZ, maxValueZ, nBinsZ, std::is_integral_v<Arithmetic> );
0798     for ( unsigned int nz = 0; nz <= nBinsZ + 1; ++nz ) {
0799       const auto offsetz = nz * ( nBinsY + 2 );
0800       for ( unsigned int ny = 0; ny <= nBinsY; ++ny ) {
0801         const auto offset = ( offsetz + ny ) * ( nBinsX + 2 );
0802         for ( unsigned int nx = 0; nx <= nBinsX; ++nx ) {
0803           auto const& [tmp, sumw2] = gaudiHisto.binValue( nx + offset );
0804           auto const& [nent, sumw] = tmp;
0805           nAllEntries += nent;
0806           if ( nx > 0 && ny > 0 && nz > 0 && nx <= nBinsX && ny <= nBinsY && nz <= nBinsZ ) {
0807             const FPType binXValue = xBinAv( nx );
0808             const FPType binYValue = yBinAv( ny );
0809             const FPType binZValue = zBinAv( nz );
0810             nEntries += nent;
0811             totalSumW += sumw;
0812             sumX += binXValue * sumw;
0813             sum2X += binXValue * binXValue * sumw;
0814             sumY += binYValue * sumw;
0815             sum2Y += binYValue * binYValue * sumw;
0816             sumZ += binZValue * sumw;
0817             sum2Z += binZValue * binZValue * sumw;
0818           }
0819         }
0820       }
0821     }
0822     const FPType meanX = fabs( totalSumW ) > 0 ? sumX / totalSumW : 0.0;
0823     const FPType stddevX =
0824         fabs( totalSumW ) > 0 ? detail::sqrt_or_zero( ( sum2X - sumX * ( sumX / totalSumW ) ) / totalSumW ) : 0.0;
0825     const FPType meanY = fabs( totalSumW ) > 0 ? sumY / totalSumW : 0.0;
0826     const FPType stddevY =
0827         fabs( totalSumW ) > 0 ? detail::sqrt_or_zero( ( sum2Y - sumY * ( sumY / totalSumW ) ) / totalSumW ) : 0.0;
0828     const FPType meanZ = fabs( totalSumW ) > 0 ? sumZ / totalSumW : 0.0;
0829     const FPType stddevZ =
0830         fabs( totalSumW ) > 0 ? detail::sqrt_or_zero( ( sum2Z - sumZ * ( sumZ / totalSumW ) ) / totalSumW ) : 0.0;
0831     // print
0832     const std::string ftitle = detail::formatTitle( gaudiHisto.title(), stringsWidth - 2 );
0833     return fmt::format( fmt::runtime( detail::histo3DFormatting ), detail::formatName( name, stringsWidth ),
0834                         stringsWidth, stringsWidth, ftitle, stringsWidth, nEntries, nAllEntries, meanX, stddevX, meanY,
0835                         stddevY, meanZ, stddevZ );
0836   }
0837 
0838   struct BinAccessor {
0839     std::function<double( unsigned int )> m_get;
0840     BinAccessor( std::string_view type, const nlohmann::json& j ) {
0841       auto subtype = std::string_view( type ).substr( 10 ); // remove "histogram:" in front
0842       if ( subtype.substr( 0, 15 ) == "WeightedProfile" || subtype.substr( 0, 7 ) == "Profile" ) {
0843         m_get = [bins = j.at( "bins" ).get<std::vector<std::tuple<std::tuple<unsigned int, double>, double>>>()](
0844                     unsigned int nx ) {
0845           auto [c, sumWeight2] = bins[nx];
0846           auto [n, sumWeight]  = c;
0847           return sumWeight;
0848         };
0849       } else {
0850         m_get = [bins = j.at( "bins" ).get<std::vector<double>>()]( unsigned int nx ) { return bins[nx]; };
0851       }
0852     }
0853     auto operator[]( const unsigned int n ) const { return m_get( n ); }
0854   };
0855 
0856   struct ArthTypeAccessor {
0857     static bool is_integral( const std::string_view type ) {
0858       for ( const std::string s : { ":i", ":l" } ) {
0859         if ( 0 == type.compare( type.length() - s.length(), s.length(), s ) ) { return true; }
0860       }
0861       return false;
0862     }
0863   };
0864 
0865   inline std::string printHistogram1D( std::string_view type, std::string_view name, std::string_view title,
0866                                        const nlohmann::json& j, unsigned int stringsWidth = 45 ) {
0867     const auto& jaxis       = j.at( "axis" );
0868     const auto  minValueX   = jaxis[0].at( "minValue" ).get<double>();
0869     const auto  maxValueX   = jaxis[0].at( "maxValue" ).get<double>();
0870     const auto  nBinsX      = jaxis[0].at( "nBins" ).get<unsigned int>();
0871     const auto  nAllEntries = j.at( "nEntries" ).get<unsigned int>();
0872     // compute the various momenta
0873     BinAccessor               ba{ type, j };
0874     double                    sumX{}, sum2X{}, sum3X{}, sum4X{}, nEntries{ 0 };
0875     const details::BinAvValue xBinAv( minValueX, maxValueX, nBinsX, ArthTypeAccessor::is_integral( type ) );
0876     for ( unsigned int nx = 1; nx <= nBinsX; ++nx ) {
0877       const double binXValue = xBinAv( nx );
0878       const auto   n         = ba[nx];
0879       nEntries += n;
0880       auto val = binXValue * n;
0881       sumX += val;
0882       val *= binXValue;
0883       sum2X += val;
0884       val *= binXValue;
0885       sum3X += val;
0886       val *= binXValue;
0887       sum4X += val;
0888     }
0889     const std::string ftitle = detail::formatTitle( title, stringsWidth - 2 );
0890     if ( !( fabs( nEntries ) > 0.0 ) ) { return ""; }
0891     const double meanX     = sumX / nEntries;
0892     const double sigmaX2   = ( sum2X / nEntries ) - std::pow( meanX, 2 );
0893     const double stddevX   = detail::sqrt_or_zero( sigmaX2 );
0894     const double EX3       = sum3X / nEntries;
0895     const double A         = sigmaX2 * stddevX;
0896     const double skewnessX = ( fabs( A ) > 0.0 ? ( EX3 - ( 3 * sigmaX2 + meanX * meanX ) * meanX ) / A : 0.0 );
0897     const double B         = sigmaX2 * sigmaX2;
0898     const double kurtosisX =
0899         ( fabs( B ) > 0.0 && fabs( A ) > 0.0
0900               ? ( sum4X / nEntries - meanX * ( 4 * EX3 - meanX * ( 6 * sigmaX2 + 3 * meanX * meanX ) ) ) / B
0901               : 3.0 );
0902     // print
0903     return fmt::format( fmt::runtime( detail::histo1DFormatting ), detail::formatName( name, stringsWidth ),
0904                         stringsWidth, stringsWidth, ftitle, stringsWidth, detail::IntWithFixedWidth{ nAllEntries },
0905                         meanX, stddevX, skewnessX, kurtosisX - 3.0 );
0906   }
0907 
0908   inline std::string printHistogram2D( std::string_view type, std::string_view name, std::string_view title,
0909                                        const nlohmann::json& j, unsigned int stringsWidth = 45 ) {
0910     const auto& jaxis       = j.at( "axis" );
0911     const auto  minValueX   = jaxis[0].at( "minValue" ).get<double>();
0912     const auto  maxValueX   = jaxis[0].at( "maxValue" ).get<double>();
0913     const auto  nBinsX      = jaxis[0].at( "nBins" ).get<unsigned int>();
0914     const auto  minValueY   = jaxis[1].at( "minValue" ).get<double>();
0915     const auto  maxValueY   = jaxis[1].at( "maxValue" ).get<double>();
0916     const auto  nBinsY      = jaxis[1].at( "nBins" ).get<unsigned int>();
0917     const auto  nAllEntries = j.at( "nEntries" ).get<double>();
0918     // compute the various memneta
0919     BinAccessor               ba{ type, j };
0920     double                    nEntries{};
0921     double                    sumX{}, sumY{};
0922     double                    sum2X{}, sum2Y{};
0923     const auto                isInt = ArthTypeAccessor::is_integral( type );
0924     const details::BinAvValue xBinAv( minValueX, maxValueX, nBinsX, isInt );
0925     const details::BinAvValue yBinAv( minValueY, maxValueY, nBinsY, isInt );
0926     for ( unsigned int ny = 1; ny <= nBinsY; ++ny ) {
0927       const auto offset = ny * ( nBinsX + 2 );
0928       for ( unsigned int nx = 1; nx <= nBinsX; ++nx ) {
0929         const double binXValue = xBinAv( nx );
0930         const double binYValue = yBinAv( ny );
0931         const auto   n         = ba[offset + nx];
0932         nEntries += n;
0933         sumX += n * binXValue;
0934         sum2X += n * binXValue * binXValue;
0935         sumY += n * binYValue;
0936         sum2Y += n * binYValue * binYValue;
0937       }
0938     }
0939     const double meanX = fabs( nEntries ) > 0 ? sumX / nEntries : 0.0;
0940     const double stddevX =
0941         fabs( nEntries ) > 0 ? detail::sqrt_or_zero( ( sum2X - sumX * ( sumX / nEntries ) ) / nEntries ) : 0.0;
0942     const double meanY = fabs( nEntries ) > 0 ? sumY / nEntries : 0.0;
0943     const double stddevY =
0944         fabs( nEntries ) > 0 ? detail::sqrt_or_zero( ( sum2Y - sumY * ( sumY / nEntries ) ) / nEntries ) : 0.0;
0945     // print
0946     const std::string ftitle = detail::formatTitle( title, stringsWidth - 2 );
0947     return fmt::format( fmt::runtime( detail::histo2DFormatting ), detail::formatName( name, stringsWidth ),
0948                         stringsWidth, stringsWidth, ftitle, stringsWidth, nEntries, nAllEntries, meanX, stddevX, meanY,
0949                         stddevY );
0950   }
0951 
0952   inline std::string printHistogram3D( std::string_view type, std::string_view name, std::string_view title,
0953                                        const nlohmann::json& j, unsigned int stringsWidth = 45 ) {
0954     const auto& jaxis       = j.at( "axis" );
0955     const auto  minValueX   = jaxis[0].at( "minValue" ).get<double>();
0956     const auto  maxValueX   = jaxis[0].at( "maxValue" ).get<double>();
0957     const auto  nBinsX      = jaxis[0].at( "nBins" ).get<unsigned int>();
0958     const auto  minValueY   = jaxis[1].at( "minValue" ).get<double>();
0959     const auto  maxValueY   = jaxis[1].at( "maxValue" ).get<double>();
0960     const auto  nBinsY      = jaxis[1].at( "nBins" ).get<unsigned int>();
0961     const auto  minValueZ   = jaxis[2].at( "minValue" ).get<double>();
0962     const auto  maxValueZ   = jaxis[2].at( "maxValue" ).get<double>();
0963     const auto  nBinsZ      = jaxis[2].at( "nBins" ).get<unsigned int>();
0964     const auto  nAllEntries = j.at( "nEntries" ).get<double>();
0965     // compute the various memneta
0966     BinAccessor               ba{ type, j };
0967     double                    nEntries{};
0968     double                    sumX{}, sumY{}, sumZ{};
0969     double                    sum2X{}, sum2Y{}, sum2Z{};
0970     const auto                isInt = ArthTypeAccessor::is_integral( type );
0971     const details::BinAvValue xBinAv( minValueX, maxValueX, nBinsX, isInt );
0972     const details::BinAvValue yBinAv( minValueY, maxValueY, nBinsY, isInt );
0973     const details::BinAvValue zBinAv( minValueZ, maxValueZ, nBinsZ, isInt );
0974     for ( unsigned int nz = 1; nz <= nBinsZ; ++nz ) {
0975       const auto offsetz = nz * ( nBinsY + 2 );
0976       for ( unsigned int ny = 1; ny <= nBinsY; ++ny ) {
0977         const auto offset = ( offsetz + ny ) * ( nBinsX + 2 );
0978         for ( unsigned int nx = 1; nx <= nBinsX; ++nx ) {
0979           const double binXValue = xBinAv( nx );
0980           const double binYValue = yBinAv( ny );
0981           const double binZValue = zBinAv( nz );
0982           const auto   n         = ba[offset + nx];
0983           nEntries += n;
0984           sumX += n * binXValue;
0985           sum2X += n * binXValue * binXValue;
0986           sumY += n * binYValue;
0987           sum2Y += n * binYValue * binYValue;
0988           sumZ += n * binZValue;
0989           sum2Z += n * binZValue * binZValue;
0990         }
0991       }
0992     }
0993     const double meanX = fabs( nEntries ) > 0 ? sumX / nEntries : 0.0;
0994     const double stddevX =
0995         fabs( nEntries ) > 0 ? detail::sqrt_or_zero( ( sum2X - sumX * ( sumX / nEntries ) ) / nEntries ) : 0.0;
0996     const double meanY = fabs( nEntries ) > 0 ? sumY / nEntries : 0.0;
0997     const double stddevY =
0998         fabs( nEntries ) > 0 ? detail::sqrt_or_zero( ( sum2Y - sumY * ( sumY / nEntries ) ) / nEntries ) : 0.0;
0999     const double meanZ = fabs( nEntries ) > 0 ? sumZ / nEntries : 0.0;
1000     const double stddevZ =
1001         fabs( nEntries ) > 0 ? detail::sqrt_or_zero( ( sum2Z - sumZ * ( sumZ / nEntries ) ) / nEntries ) : 0.0;
1002     // print
1003     const std::string ftitle = detail::formatTitle( title, stringsWidth - 2 );
1004     return fmt::format( fmt::runtime( detail::histo3DFormatting ), detail::formatName( name, stringsWidth ),
1005                         stringsWidth, stringsWidth, ftitle, stringsWidth, nEntries, nAllEntries, meanX, stddevX, meanY,
1006                         stddevY, meanZ, stddevZ );
1007   }
1008 
1009 } // namespace Gaudi::Histograming::Sink
1010 
1011 namespace nlohmann {
1012   /// automatic translation of TAxis to json
1013   inline void to_json( nlohmann::json& j, TAxis const& axis ) {
1014     j = nlohmann::json{ { "nBins", axis.GetNbins() },
1015                         { "minValue", axis.GetXmin() },
1016                         { "maxValue", axis.GetXmax() },
1017                         { "title", axis.GetTitle() } };
1018   }
1019   /// automatic translation of Root histograms to json
1020   inline void to_json( nlohmann::json& j, TH1D const& h ) {
1021     j = Gaudi::Histograming::Sink::details::rootHistogramToJson( h );
1022   }
1023   inline void to_json( nlohmann::json& j, TH2D const& h ) {
1024     j = Gaudi::Histograming::Sink::details::rootHistogramToJson( h );
1025   }
1026   inline void to_json( nlohmann::json& j, TH3D const& h ) {
1027     j = Gaudi::Histograming::Sink::details::rootHistogramToJson( h );
1028   }
1029   inline void to_json( nlohmann::json& j, TProfile const& h ) {
1030     j = Gaudi::Histograming::Sink::details::rootHistogramToJson( h );
1031   }
1032   inline void to_json( nlohmann::json& j, TProfile2D const& h ) {
1033     j = Gaudi::Histograming::Sink::details::rootHistogramToJson( h );
1034   }
1035 
1036 } // namespace nlohmann
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!