Skip to content
Snippets Groups Projects
Commit d7d2f5ac authored by Sebastien Ponce's avatar Sebastien Ponce
Browse files

Moved most of cellid based histo sink code to LHCb

parent 198ce53d
No related branches found
No related tags found
1 merge request!3107Modernized Calo monitoring and dropped MoniAlg
...@@ -9,208 +9,12 @@ ...@@ -9,208 +9,12 @@
* or submit itself to any jurisdiction. * * or submit itself to any jurisdiction. *
\***********************************************************************************/ \***********************************************************************************/
#include <CaloFutureUtils/CellIDHistogram.h> #include <CaloFutureUtils/CellIDHistogramUtils.h>
#include <Detector/Calo/CaloCellID.h>
#include <HistogramPersistencySvc/RootHistogramSinkBase.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 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 { namespace LHCb::Calo::Histograming::Sink {
using namespace std::string_literals; using namespace std::string_literals;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment