Skip to content
Snippets Groups Projects
Commit 3d9152ed authored by Gerardo Ganis's avatar Gerardo Ganis Committed by Danilo Piparo
Browse files

hist/TH1Merger: add support for merging histograms filled in power-of-2 mode 

This requires a special treatment and a new use case
parent 429c109b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
hist/hist/src/TH1Merger.cxx
+ 344
− 26
View file @ 3d9152ed
...@@ -9,7 +9,12 @@ ...@@ -9,7 +9,12 @@
#include "THashList.h" #include "THashList.h"
#include "TClass.h" #include "TClass.h"
#include <iostream> #include <iostream>
#include <limits>
#include <utility>
#define PRINTRANGE(a, b, bn) \
Printf(" base: %f %f %d, %s: %f %f %d", a->GetXmin(), a->GetXmax(), a->GetNbins(), \
bn, b->GetXmin(), b->GetXmax(), b->GetNbins());
Bool_t TH1Merger::AxesHaveLimits(const TH1 * h) { Bool_t TH1Merger::AxesHaveLimits(const TH1 * h) {
Bool_t hasLimits = h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax(); Bool_t hasLimits = h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax();
...@@ -37,6 +42,9 @@ Bool_t TH1Merger::operator() () { ...@@ -37,6 +42,9 @@ Bool_t TH1Merger::operator() () {
if (type == kAllNoLimits) if (type == kAllNoLimits)
return BufferMerge(); return BufferMerge();
if (type == kAutoP2HaveLimits || (type == kAutoP2NeedLimits && AutoP2BufferMerge()))
return AutoP2Merge();
// this is the mixed case - more complicated // this is the mixed case - more complicated
if (type == kHasNewLimits) { if (type == kHasNewLimits) {
// we need to define some new axes // we need to define some new axes
...@@ -52,6 +60,144 @@ Bool_t TH1Merger::operator() () { ...@@ -52,6 +60,144 @@ Bool_t TH1Merger::operator() () {
return kFALSE; return kFALSE;
} }
/////////////////////////////////////////////////////////////////////////////////////////
/// Determine final boundaries and number of bins for histograms created in power-of-2
/// autobin mode.
///
/// Return kTRUE if compatible, updating fNewXaxis accordingly; return kFALSE if something
/// wrong.
///
/// The histograms are not merge-compatible if
///
/// 1. have different variable-size bins
/// 2. larger bin size is not an integer multiple of the smaller one
/// 3. the final estimated range is smalle then the bin size
///
Bool_t TH1Merger::AutoP2BuildAxes(TH1 *h)
{
// They must be both defined
if (!h) {
Error("AutoP2BuildAxes", "undefined histogram: %p", h);
return kFALSE;
}
// They must be created in power-of-2 autobin mode
if (!h->TestBit(TH1::kAutoBinPTwo)) {
Error("AutoP2BuildAxes", "not in autobin-power-of-2 mode!");
return kFALSE;
}
// Point to axes
TAxis *a0 = &fNewXAxis, *a1 = h->GetXaxis();
// This is for future merging of detached ranges (only possible if no over/underflows)
Bool_t canextend = (h->GetBinContent(0) > 0 || h->GetBinContent(a1->GetNbins() + 1) > 0) ? kFALSE : kTRUE;
// The first time we just copy the boundaries and bins
if (a0->GetFirst() == a0->GetLast()) {
a0->Set(a1->GetNbins(), a1->GetXmin(), a1->GetXmax());
// This is for future merging of detached ranges (only possible if no over/underflows)
a0->SetCanExtend(canextend);
return kTRUE;
}
// Bin sizes must be in integer ratio
Double_t bwmax = (a0->GetXmax() - a0->GetXmin()) / a0->GetNbins() ;
Double_t bwmin = (a1->GetXmax() - a1->GetXmin()) / a1->GetNbins() ;
Bool_t b0 = kTRUE;
if (bwmin > bwmax) {
std::swap(bwmax, bwmin);
b0 = kFALSE;
}
if (!(bwmin > 0.)) {
PRINTRANGE(a0, a1, h->GetName());
Error("AutoP2BuildAxes", "minimal bin width negative or null: %f", bwmin);
return kFALSE;
}
Double_t rt;
Double_t re = std::modf(bwmax / bwmin, &rt);
if (re > std::numeric_limits<Double_t>::epsilon()) {
PRINTRANGE(a0, a1, h->GetName());
Error("AutoP2BuildAxes", "bin widths not in integer ratio: %f", re);
return kFALSE;
}
// Range of the merged histogram, taking into account overlaps
Bool_t domax = kFALSE;
Double_t xmax, xmin;
if (a0->GetXmin() < a1->GetXmin()) {
if (a0->GetXmax() < a1->GetXmin()) {
if (!a0->CanExtend() || !canextend) {
PRINTRANGE(a0, a1, h->GetName());
Error("AutoP2BuildAxes", "ranges are disconnected and under/overflows: cannot merge");
return kFALSE;
}
xmax = a1->GetXmax();
xmin = a0->GetXmin();
domax = b0 ? kTRUE : kFALSE;
} else {
if (a0->GetXmax() >= a1->GetXmax()) {
xmax = a1->GetXmax();
xmin = a1->GetXmin();
domax = !b0 ? kTRUE : kFALSE;
} else {
xmax = a0->GetXmax();
xmin = a1->GetXmin();
domax = !b0 ? kTRUE : kFALSE;
}
}
} else {
if (a1->GetXmax() < a0->GetXmin()) {
if (!a0->CanExtend() || !canextend) {
PRINTRANGE(a0, a1, h->GetName());
Error("AutoP2BuildAxes", "ranges are disconnected and under/overflows: cannot merge");
return kFALSE;
}
xmax = a0->GetXmax();
xmin = a1->GetXmin();
domax = !b0 ? kTRUE : kFALSE;
} else {
if (a1->GetXmax() >= a0->GetXmax()) {
xmax = a0->GetXmax();
xmin = a0->GetXmin();
domax = b0 ? kTRUE : kFALSE;
} else {
xmax = a1->GetXmax();
xmin = a0->GetXmin();
domax = b0 ? kTRUE : kFALSE;
}
}
}
Double_t range = xmax - xmin;
re = std::modf(range / bwmax, &rt);
if (rt < 1.) {
PRINTRANGE(a0, a1, h->GetName());
Error("MergeCompatibleHistograms", "range smaller than bin width: %f %f %f", range, bwmax, rt);
return kFALSE;
}
if (re > std::numeric_limits<Double_t>::epsilon()) {
if (domax) {
xmax -= bwmax * re;
} else {
xmin += bwmax * re;
}
}
// Number of bins
Int_t nb = (Int_t)rt;
// Set the result
a0->Set(nb, xmin, xmax);
// This is for future merging of detached ranges (only possible if no over/underflows)
if (!a0->CanExtend()) a0->SetCanExtend(canextend);
// Done
return kTRUE;
}
/** /**
Examine the list of histograms to find out which type of Merge we need to do Examine the list of histograms to find out which type of Merge we need to do
Pass the input list containing the histogram to merge and h0 which is the initial histogram Pass the input list containing the histogram to merge and h0 which is the initial histogram
...@@ -80,6 +226,8 @@ TH1Merger::EMergerType TH1Merger::ExamineHistograms() { ...@@ -80,6 +226,8 @@ TH1Merger::EMergerType TH1Merger::ExamineHistograms() {
Bool_t foundLabelHist = kFALSE; Bool_t foundLabelHist = kFALSE;
Bool_t haveWeights = kFALSE; Bool_t haveWeights = kFALSE;
Bool_t isAutoP2 = kFALSE;
// TAxis newXAxis; // TAxis newXAxis;
// TAxis newYAxis; // TAxis newYAxis;
// TAxis newZAxis; // TAxis newZAxis;
...@@ -89,6 +237,8 @@ TH1Merger::EMergerType TH1Merger::ExamineHistograms() { ...@@ -89,6 +237,8 @@ TH1Merger::EMergerType TH1Merger::ExamineHistograms() {
int dimension = fH0->GetDimension(); int dimension = fH0->GetDimension();
isAutoP2 = fH0->TestBit(TH1::kAutoBinPTwo) ? kTRUE : kFALSE;
// start looping on the histograms // start looping on the histograms
do { do {
...@@ -111,7 +261,16 @@ TH1Merger::EMergerType TH1Merger::ExamineHistograms() { ...@@ -111,7 +261,16 @@ TH1Merger::EMergerType TH1Merger::ExamineHistograms() {
allHaveLimits = allHaveLimits && hasLimits; allHaveLimits = allHaveLimits && hasLimits;
allSameLimits &= allHaveLimits; allSameLimits &= allHaveLimits;
if (isAutoP2 && !h->TestBit(TH1::kAutoBinPTwo)) {
Error("Merge", "Cannot merge histogram - some are in autobin-power-of-2 mode, but not %s!", h->GetName());
return kNotCompatible;
}
if (!isAutoP2 && h->TestBit(TH1::kAutoBinPTwo)) {
Error("Merge", "Cannot merge histogram - %s is in autobin-power-of-2 mode, but not the previous ones",
h->GetName());
return kNotCompatible;
}
if (hasLimits) { if (hasLimits) {
h->BufferEmpty(); h->BufferEmpty();
...@@ -270,8 +429,16 @@ TH1Merger::EMergerType TH1Merger::ExamineHistograms() { ...@@ -270,8 +429,16 @@ TH1Merger::EMergerType TH1Merger::ExamineHistograms() {
} }
// in case of weighted histogram set Sumw2() on fH0 is is not weighted // in case of weighted histogram set Sumw2() on fH0 is is not weighted
if (haveWeights && fH0->GetSumw2N() == 0) fH0->Sumw2(); if (haveWeights && fH0->GetSumw2N() == 0)
fH0->Sumw2();
// AutoP2
if (isAutoP2) {
if (allHaveLimits)
return kAutoP2HaveLimits;
return kAutoP2NeedLimits;
}
// return the type of merge // return the type of merge
if (allHaveLabels) return kAllLabel; if (allHaveLabels) return kAllLabel;
if (allSameLimits) return kAllSameAxes; if (allSameLimits) return kAllSameAxes;
...@@ -369,7 +536,179 @@ void TH1Merger::DefineNewAxes() { ...@@ -369,7 +536,179 @@ void TH1Merger::DefineNewAxes() {
} }
Bool_t TH1Merger::BufferMerge() { void TH1Merger::CopyBuffer(TH1 *hsrc, TH1 *hdes)
{
// Check inputs
if (!hsrc || !hsrc->fBuffer || !hdes || !hdes->fBuffer) {
void *p1 = hsrc ? hsrc->fBuffer : 0;
void *p2 = hdes ? hdes->fBuffer : 0;
Warning("TH1Merger::CopyMerge", "invalid inputs: %p, %p, %p, %p -> do nothing", hsrc, hdes, p1, p2);
}
// Entries from buffers have to be filled one by one
// because FillN doesn't resize histograms.
Int_t nbentries = (Int_t)hsrc->fBuffer[0];
if (hdes->fDimension == 1) {
for (Int_t i = 0; i < nbentries; i++)
hdes->Fill(hsrc->fBuffer[2 * i + 2], hsrc->fBuffer[2 * i + 1]);
}
if (hdes->fDimension == 2) {
auto h2 = dynamic_cast<TH2 *>(hdes);
R__ASSERT(h2);
for (Int_t i = 0; i < nbentries; i++)
h2->Fill(hsrc->fBuffer[3 * i + 2], hsrc->fBuffer[3 * i + 3], hsrc->fBuffer[3 * i + 1]);
}
if (hdes->fDimension == 3) {
auto h3 = dynamic_cast<TH3 *>(hdes);
R__ASSERT(h3);
for (Int_t i = 0; i < nbentries; i++)
h3->Fill(hsrc->fBuffer[4 * i + 2], hsrc->fBuffer[4 * i + 3], hsrc->fBuffer[4 * i + 4],
hsrc->fBuffer[4 * i + 1]);
}
}
Bool_t TH1Merger::AutoP2BufferMerge()
{
TH1 *href = 0, *hist = 0;
TIter nextref(&fInputList);
if (TH1Merger::AxesHaveLimits(fH0)) {
href = fH0;
} else {
while ((hist = (TH1 *)nextref()) && !href) {
if (TH1Merger::AxesHaveLimits(hist))
href = hist;
}
}
Bool_t resetfH0 = kFALSE;
if (!href) {
// Merge all histograms to fH0 and do a final projection
href = fH0;
} else {
if (href != fH0) {
// Temporary add fH0 to the list for buffer merging
fInputList.Add(fH0);
resetfH0 = kTRUE;
}
}
TIter next(&fInputList);
while ((hist = (TH1 *)next())) {
if (!TH1Merger::AxesHaveLimits(hist) && hist->fBuffer) {
if (gDebug)
Info("AutoP2BufferMerge", "merging buffer of %s into %s", hist->GetName(), href->GetName());
CopyBuffer(hist, href);
fInputList.Remove(hist);
}
}
// Final projection
if (href->fBuffer)
href->BufferEmpty(1);
// Reset fH0, if already added, to avoid double counting
if (resetfH0)
fH0->Reset("ICES");
// Done, all histos have been processed
return kTRUE;
}
Bool_t TH1Merger::AutoP2Merge()
{
Double_t stats[TH1::kNstat], totstats[TH1::kNstat];
for (Int_t i = 0; i < TH1::kNstat; i++) {
totstats[i] = stats[i] = 0;
}
TIter next(&fInputList);
TH1 *hist = 0;
// Calculate boundaries and bins
Double_t xmin = 0., xmax = 0.;
if (!(fH0->IsEmpty())) {
hist = fH0;
} else {
while ((hist = (TH1 *)next())) {
if (!hist->IsEmpty())
break;
}
}
if (!hist) {
Warning("TH1Merger::AutoP2Merge", "all histograms look empty!");
return kFALSE;
}
// Start building the axes from the reference histogram
if (!AutoP2BuildAxes(hist)) {
Error("TH1Merger::AutoP2Merge", "cannot create axes from %s", hist->GetName());
return kFALSE;
}
TH1 *h = 0;
while ((h = (TH1 *)next())) {
if (!AutoP2BuildAxes(h)) {
Error("TH1Merger::AutoP2Merge", "cannot merge histogram %s: not merge compatible", h->GetName());
return kFALSE;
}
}
xmin = fNewXAxis.GetXmin();
xmax = fNewXAxis.GetXmax();
Int_t nbins = fNewXAxis.GetNbins();
// Prepare stats
fH0->GetStats(totstats);
// Clone fH0 and add it to the list
if (!fH0->IsEmpty())
fInputList.Add(fH0->Clone());
// reset fH0
fH0->Reset("ICES");
// Set the new boundaries
fH0->SetBins(nbins, xmin, xmax);
next.Reset();
Double_t nentries = 0.;
while ((hist = (TH1 *)next())) {
// process only if the histogram has limits; otherwise it was processed before
// in the case of an existing buffer (see if statement just before)
if (gDebug)
Info("TH1Merger::AutoP2Merge", "merging histogram %s into %s (entries: %f)", hist->GetName(), fH0->GetName(),
hist->GetEntries());
// skip empty histograms
if (hist->IsEmpty())
continue;
// import statistics
hist->GetStats(stats);
for (Int_t i = 0; i < TH1::kNstat; i++)
totstats[i] += stats[i];
nentries += hist->GetEntries();
// Int_t nx = hist->GetXaxis()->GetNbins();
// loop on bins of the histogram and do the merge
for (Int_t ibin = 0; ibin < hist->fNcells; ibin++) {
Double_t cu = hist->RetrieveBinContent(ibin);
Double_t e1sq = TMath::Abs(cu);
if (fH0->fSumw2.fN)
e1sq = hist->GetBinErrorSqUnchecked(ibin);
Double_t xu = hist->GetBinCenter(ibin);
Int_t jbin = fH0->FindBin(xu);
fH0->AddBinContent(jbin, cu);
if (fH0->fSumw2.fN)
fH0->fSumw2.fArray[jbin] += e1sq;
}
}
// copy merged stats
fH0->PutStats(totstats);
fH0->SetEntries(nentries);
return kTRUE;
}
Bool_t TH1Merger::BufferMerge()
{
TIter next(&fInputList); TIter next(&fInputList);
while (TH1* hist = (TH1*)next()) { while (TH1* hist = (TH1*)next()) {
...@@ -378,28 +717,7 @@ Bool_t TH1Merger::BufferMerge() { ...@@ -378,28 +717,7 @@ Bool_t TH1Merger::BufferMerge() {
if (gDebug) if (gDebug)
Info("TH1Merger::BufferMerge","Merging histogram %s into %s",hist->GetName(), fH0->GetName() ); Info("TH1Merger::BufferMerge","Merging histogram %s into %s",hist->GetName(), fH0->GetName() );
CopyBuffer(hist, fH0);
// case of no limits
// Entries from buffers have to be filled one by one
// because FillN doesn't resize histograms.
Int_t nbentries = (Int_t)hist->fBuffer[0];
if (fH0->fDimension == 1) {
for (Int_t i = 0; i < nbentries; i++)
fH0->Fill(hist->fBuffer[2*i + 2], hist->fBuffer[2*i + 1]);
}
if (fH0->fDimension == 2) {
auto h2 = dynamic_cast<TH2*>(fH0);
R__ASSERT(h2);
for (Int_t i = 0; i < nbentries; i++)
h2->Fill(hist->fBuffer[3*i + 2], hist->fBuffer[3*i + 3],hist->fBuffer[3*i + 1] );
}
if (fH0->fDimension == 3) {
auto h3 = dynamic_cast<TH3*>(fH0);
R__ASSERT(h3);
for (Int_t i = 0; i < nbentries; i++)
h3->Fill(hist->fBuffer[4*i + 2], hist->fBuffer[4*i + 3],hist->fBuffer[4*i + 4], hist->fBuffer[4*i + 1] );
}
fInputList.Remove(hist); fInputList.Remove(hist);
} }
} }
......
hist/hist/src/TH1Merger.h
+ 11
− 1
View file @ 3d9152ed
...@@ -23,7 +23,9 @@ public: ...@@ -23,7 +23,9 @@ public:
kAllSameAxes = 0, // histogram have all some axes kAllSameAxes = 0, // histogram have all some axes
kAllNoLimits = 1, // all histogram don't have limits (the buffer is used) kAllNoLimits = 1, // all histogram don't have limits (the buffer is used)
kHasNewLimits = 2, // all histogram don't have limits (the buffer is used) kHasNewLimits = 2, // all histogram don't have limits (the buffer is used)
kAllLabel = 3 // histogram have labels all axis kAllLabel = 3, // histogram have labels all axis
kAutoP2HaveLimits = 4, // P2 (power-of-2) algorithm: all histogram have limits
kAutoP2NeedLimits = 5 // P2 algorithm: some histogram still need projections
}; };
static Bool_t AxesHaveLimits(const TH1 * h); static Bool_t AxesHaveLimits(const TH1 * h);
...@@ -67,12 +69,20 @@ public: ...@@ -67,12 +69,20 @@ public:
private: private:
Bool_t AutoP2BuildAxes(TH1 *);
EMergerType ExamineHistograms(); EMergerType ExamineHistograms();
void DefineNewAxes(); void DefineNewAxes();
void CopyBuffer(TH1 *hsrc, TH1 *hdes);
Bool_t BufferMerge(); Bool_t BufferMerge();
Bool_t AutoP2BufferMerge();
Bool_t AutoP2Merge();
Bool_t SameAxesMerge(); Bool_t SameAxesMerge();
Bool_t DifferentAxesMerge(); Bool_t DifferentAxesMerge();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment