//////////////////////////////////////////////////////////////////////////
// //
// GFHistManip //
//
// Author: Gero Flucke
// Date: April 2nd, 2003
// last update: $Date: 2005/12/15 11:21:32 $
// by: $Author: flucke $
// //
// Encapsulate Gero's routines for histogram manipulation. //
// //
//////////////////////////////////////////////////////////////////////////
#include "GFUtils/GFHistManip.h"
#include <iostream>
#include <string.h>
#include <limits.h>
#include <stdlib.h>
#include "TMath.h"
#include "TString.h" // for Form ?
#include "TError.h"
#include "TObjArray.h"
#include "TArrayD.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TGraphAsymmErrors.h"
#include "TAxis.h"
#include "TClass.h"
#include "TDirectory.h"
#include "TROOT.h"
#include "GFUtils/GFHistArray.h"
#include "GFUtils/GFMath.h"
using namespace std;
class TMatrixD;
ClassImp(GFHistManip)
//__________________________________________________________________
TH1* GFHistManip::CreateMaxRelDiff(TObjArray* hists, const TH1* reference)
{
if(!CheckContent(hists, TH1::Class()) || !reference) return NULL;
// check binning??? also for multi dimensional hists?
TH1 *reldiffs = CreateAnalog(reference, "relDiff", kFALSE); // no sumw2 structure, no fill
// FIXME: IsSameBinning!
for(Int_t bin = 1; bin <= reference->GetNbinsX(); ++bin){
Double_t refNum = reference->GetBinContent(bin);
if(!refNum){
cerr << "<GFHistManip::CreateMaxRelDiff>: reference histogram with bin " <<
bin << " = " << refNum << endl;
continue;
}
for(Int_t nHist = 0; nHist < hists->GetEntriesFast(); ++nHist){
if(!hists->At(nHist)) continue;
TH1* h = static_cast<TH1*>(hists->At(nHist));
Double_t curNum = h->GetBinContent(bin);
Double_t relDiff = TMath::Abs(refNum-curNum)/curNum;
if(relDiff > reldiffs->GetBinContent(bin)){
reldiffs->SetBinContent(bin, relDiff);
}
}
}
return reldiffs;
}
//__________________________________________________________________
TH1* GFHistManip::CreateHistDevSquareSum(const GFHistArray *hists, const TH1 *hRef,
Int_t flagUpLow)
{
// quadaratic adding up of deviations of hists to hRef,
// hists must be above (below) hRef for flagUpLow > 0 (< 0)
// 1-D only
if(flagUpLow == 0 || !hists || !hists->GetEntriesFast() || !hRef || hRef->GetDimension() !=1){
return NULL;
}
for(Int_t i = 0; i < hists->GetEntriesFast(); ++i){ // check same number of bins
if(!IsSameBinning(hists->At(i), hRef)){
::Error("GFHistManip::CreateHistDevSquareSum", "hist %d other bins than hRef", i);
return NULL;
}
}
const Int_t sign = TMath::Sign(1, flagUpLow);
TString name(Form("%s%sSqrErr", hRef->GetName(), (flagUpLow>0 ? "Upp" : "Low")));
GFHistManip::MakeUniqueName(name);
TH1 *hResult = CreateAnalog(hRef, name, kFALSE); // no sumw2 structure, no fill
for(Int_t iBin = 1; iBin <= hRef->GetNbinsX(); ++iBin){
const Double_t centerCont = hRef->GetBinContent(iBin);
Double_t sumDeviation2 = 0.;
for(Int_t iHist = 0; iHist < hists->GetEntriesFast(); ++iHist){
const Double_t histCont = hists->At(iHist)->GetBinContent(iBin);
if((histCont - centerCont) * sign < 0.){
::Error("GFHistManip::CreateHistDevSquareSum", "hist %d is %f in bin %d, center is %f",
iHist, histCont, iBin, centerCont);
delete hResult;
return NULL;
}
sumDeviation2 += ((histCont - centerCont)*(histCont - centerCont));
}
hResult->SetBinContent(iBin, TMath::Sqrt(sumDeviation2)*sign + centerCont);
}
return hResult;
}
//__________________________________________________________________
TH1* GFHistManip::CreateHistMinMax(const GFHistArray *hists, Int_t minMax)
{
// only 1D!
if(!hists || hists->GetEntriesFast() < 1 || minMax == 0) return NULL;
TIter iter(hists);
TH1 *h1 = hists->First();
while(TH1 *h = static_cast<TH1*>(iter.Next())){
if(h->GetDimension() != 1 || !IsSameBinning(h, h1)){
::Error("GFHistManip::CreateHistMinMax", "%s has different binning or > 1D",h->GetName());
return NULL;
}
}
const Int_t sign = TMath::Sign(1, minMax); // no sumw2!
TH1* result = CreateAnalog(h1, Form("%s%s", h1->GetName(), (sign>0 ? "Max" : "Min")), kFALSE);
for(Int_t iBin = 0; iBin <= h1->GetNbinsX(); ++iBin){
Double_t content = -sign * DBL_MAX;
iter.Reset();
while(TH1 *h = static_cast<TH1*>(iter.Next())){
if(sign * (h->GetBinContent(iBin)-content) > 0.){
content = h->GetBinContent(iBin);
}
}
result->SetBinContent(iBin, content);
}
return result;
}
//__________________________________________________________________
GFHistArray* GFHistManip::CreateHistArray1D(TH2* hist2D, const char* varName,
const char* histName, Option_t *option)
{
// create array of 1D hists projected from hist2D
// var name is for title of 1D hists which tells about bin borders
// if histName is beginning of names of 1D hists (hist2D->GetName() if NULL)
// if 'option' contains 'o', include overflow, if 'u' underflow
// if 'option' contains 'm' followed by a number, merge bins
// (do not add 'm' several times, the first is taken, even if there is
// no number after it!)
// by default projects along x, but if option contains 'y' do it along y
TString opt(option);
opt.ToLower();
Int_t nMerge = 1;
if (opt.Contains("m")) {
char *endPtr = NULL;
nMerge = strtol(opt.Data()+opt.Index("m")+1, &endPtr, 0);
if (endPtr == opt.Data()+opt.Index("m")+1) {
::Error("CreateHistArray1D", "miss a number after option 'm' in %s", option);
nMerge = 1;
}
}
return CreateHistArray1D(hist2D, varName, histName, !opt.Contains("u"), !opt.Contains("o"),
!opt.Contains("y"), nMerge);
}
//__________________________________________________________________
GFHistArray* GFHistManip::CreateHistArray1D(TH2* hist2D, const char* varName,
const char* histName, Bool_t ignoreOverUnderFlow,
Bool_t projectX, UInt_t mergeNBins)
{
// outdated, use version with 'option'
return CreateHistArray1D(hist2D, varName, histName, ignoreOverUnderFlow,
ignoreOverUnderFlow, projectX, mergeNBins);
}
//__________________________________________________________________
GFHistArray* GFHistManip::CreateHistArray1D(TH2* hist2D, const char* varName,
const char* histName, Bool_t ignoreUnderFlow,
Bool_t ignoreOverFlow,
Bool_t projectX, UInt_t mergeNBins)
{
// private helper method
TString nameStr(histName ? histName : hist2D->GetName());
if(mergeNBins < 1) mergeNBins = 1;
const Int_t numBins = projectX ? hist2D->GetNbinsY() : hist2D->GetNbinsX();
GFHistArray* array = new GFHistArray(numBins+2);//+2 maybe for over-/underflow
TAxis *axis = projectX ? hist2D->GetYaxis() : hist2D->GetXaxis();
if (!ignoreUnderFlow) {
TString curName(Form("%s_0", nameStr.Data()));
GFHistManip::MakeUniqueName(curName);
TH1D* hist = GFHistManip::CreateProjection(hist2D, curName, 0, 0, projectX);
hist->SetEntries(hist->Integral());// needed???
hist->SetTitle(Form("%s < %f", varName, axis->GetBinUpEdge(0)));
array->Add(hist);
}
for(Long_t i = 1; i <= numBins; i += mergeNBins){
Int_t lastBin = i + mergeNBins - 1;
if(lastBin > numBins) lastBin = numBins;
TString title(Form("%f #leq %s < %f", axis->GetBinLowEdge(i), varName,
axis->GetBinUpEdge(lastBin)));
TString curName = nameStr + Form("_%d", i);
GFHistManip::MakeUniqueName(curName);//ROOT refills existing hist in Projection!
TH1D* hist = GFHistManip::CreateProjection(hist2D, curName, i, lastBin, projectX);
hist->SetEntries(hist->Integral());// needed???
hist->SetTitle(title);
array->Add(hist);
}
if (!ignoreOverFlow) {
TString curName(Form("%s_%d", nameStr.Data(), numBins+1));
GFHistManip::MakeUniqueName(curName);
TH1D* hist = GFHistManip::CreateProjection(hist2D, curName,numBins+1,numBins+1, projectX);
hist->SetEntries(hist->Integral());// needed???
hist->SetTitle(Form("%f #leq %s", axis->GetBinLowEdge(numBins+1), varName));
array->Add(hist);
}
return array;
}
//________________________________________________
TH1D* GFHistManip::CreateProjection(const TH2* h2D, const char *name,
Int_t bin1, Int_t bin2, Bool_t projectX)
{
// Create projection of 'h2D' along x/y from 'bin1' to 'bin2', name of new hit is 'name'.
// Error option is selected if h2D has errors stored,
// sumw2 is fixed as in ROOT >= 4.01.
if(!h2D) return NULL;
TH1D* result = NULL;
if(name == NULL) name = "";// ProjectionX/Y crashes for NULL...
if(h2D->GetSumw2N()){ // if special errors exist: compute new errors!
result = projectX ?
h2D->ProjectionX(name, bin1, bin2,"E") : h2D->ProjectionY(name, bin1, bin2, "E");
#if ROOT_VERSION_CODE < ROOT_VERSION(4,4,0) // (4,2,0) ??
// due to bug in ROOT (still present in 4.00_08) sumw2 has to be fixed here 'by hand'
Stat_t s[11];
result->GetStats(s);
s[1] = 0.;
Int_t firstBin = result->GetXaxis()->GetFirst();// what if we wnat to include underflow?
Int_t lastBin = result->GetXaxis()->GetLast(); // and overflow?
for(Int_t i = firstBin; i <= lastBin; ++i){
Stat_t err = result->GetBinError(i);
s[1] += (err*err);
}
result->PutStats(s);
#endif
} else {
result = projectX ?
h2D->ProjectionX(name, bin1, bin2) : h2D->ProjectionY(name, bin1, bin2);
}
return result;
}
//________________________________________________
TH1* GFHistManip::CreateHistYBins(TH2* hist2D, const char* titleAdd, UInt_t mergeNBins)
{
// merge 'mergeNBins'
TString name(hist2D->GetName());
name.Append("yAxis");
GFHistManip::MakeUniqueName(name);
TString title(hist2D->GetTitle());
if(titleAdd) title += titleAdd;
else title += ", y-Axis";
title += Form(";%s", hist2D->GetYaxis()->GetTitle());
if(mergeNBins < 2){ // 0 makes no sense, treat as 1
const TArrayD* bins = hist2D->GetYaxis()->GetXbins();
return (0 == bins->GetSize())
? new TH1F(name, title, hist2D->GetNbinsY(),
hist2D->GetYaxis()->GetXmin(), hist2D->GetYaxis()->GetXmax())
: new TH1F(name, title, bins->GetSize()-1, bins->GetArray());
} else {
TArrayD newBins(hist2D->GetNbinsY()/mergeNBins + 2);// +2, not +1: possible 'left over' bin
Int_t j = 0;
for(Int_t i = 1; i <= hist2D->GetNbinsY(); i+=mergeNBins){
newBins[j] = hist2D->GetYaxis()->GetBinLowEdge(i);
++j;
}
newBins[j] = hist2D->GetYaxis()->GetXmax();
return new TH1F(name, title, j, newBins.GetArray());
}
}
//________________________________________________
TH1* GFHistManip::CreateHistXBins(TH2* hist2D, const char* titleAdd, UInt_t mergeNBins)
{
TString name(hist2D->GetName());
name.Append("xAxis");
GFHistManip::MakeUniqueName(name);
TString title(hist2D->GetTitle());
if(titleAdd) title += titleAdd;
else title += ", x-Axis";
title += Form(";%s", hist2D->GetXaxis()->GetTitle());
if(mergeNBins < 2){ // 0 makes no sense, treat as 1
const TArrayD* bins = hist2D->GetXaxis()->GetXbins();
return (0 == bins->GetSize())
? new TH1F(name, title, hist2D->GetNbinsX(),
hist2D->GetXaxis()->GetXmin(), hist2D->GetXaxis()->GetXmax())
: new TH1F(name, title, bins->GetSize()-1, bins->GetArray());
} else {
TArrayD newBins(hist2D->GetNbinsX()/mergeNBins + 2);// +2, not +1: possible 'left over' bin
Int_t j = 0;
for(Int_t i = 1; i <= hist2D->GetNbinsX(); i+=mergeNBins){
newBins[j] = hist2D->GetXaxis()->GetBinLowEdge(i);
++j;
}
newBins[j] = hist2D->GetXaxis()->GetXmax();
return new TH1F(name, title, j, newBins.GetArray());
}
}
//_________________________________________________
TH2* GFHistManip::CreateHistFlipXY(TH2* h)
{
// create a hist with flipped X and Y bin contents (ignore overflow bins...)
// add 'Flip' to name, always use arrays for bins
// stats might be slightly different...
if(!h) return NULL;
TArrayD newXbins(h->GetNbinsY()+1);
for(Int_t i = 1; i < newXbins.GetSize(); ++i){
newXbins[i-1] = h->GetYaxis()->GetBinLowEdge(i);
}
newXbins[newXbins.GetSize()-1] = h->GetYaxis()->GetXmax();
TArrayD newYbins(h->GetNbinsX()+1);
for(Int_t i = 1; i < newYbins.GetSize(); ++i){
newYbins[i-1] = h->GetXaxis()->GetBinLowEdge(i);
}
newYbins[newYbins.GetSize()-1] = h->GetXaxis()->GetXmax();
TH2* result = NULL;
if(h->InheritsFrom(TH2D::Class())){
result = new TH2D("tmpname", h->GetTitle(), newXbins.GetSize()-1, newXbins.GetArray(),
newYbins.GetSize()-1, newYbins.GetArray());
} else {
result = new TH2F("tmpname", h->GetTitle(), newXbins.GetSize()-1, newXbins.GetArray(),
newYbins.GetSize()-1, newYbins.GetArray());
}
result->SetName(Form("%sFlip", h->GetName()));
result->SetXTitle(h->GetYaxis()->GetTitle());
result->SetYTitle(h->GetXaxis()->GetTitle());
// including over and underflow:
for(Int_t iy = 0; iy <= result->GetNbinsY()+1; ++iy){
for(Int_t ix = 0; ix <= result->GetNbinsX()+1; ++ix){
result->SetBinContent(ix, iy, h->GetBinContent(iy, ix));
if(h->GetSumw2N()) result->SetBinError(ix, iy, h->GetBinError(iy, ix));
}
}
result->SetEntries(h->GetEntries());
Stat_t s[11];
h->GetStats(s);
Stat_t sumwx = s[4];
Stat_t sumwx2 = s[5];
Stat_t sumwy = s[2];
Stat_t sumwy2 = s[3];
s[2] = sumwx;
s[3] = sumwx2;
s[4] = sumwy;
s[5] = sumwy2;
result->PutStats(s);
return result;
}
//_________________________________________________
void GFHistManip::PrintTable(const TH1* hist, const char* title)
{
// formated printing the bin content and error, starting with title
const char line[]
= "------------------------------------------------------------------------------n";
if(!hist || hist->GetDimension() > 1) {
printf("%s| no hist or > 1D for %-54s |n%s", line, title, line);
return;
}
printf("%s| %-74s |n%s", line, title, line);
Int_t writtenBins = 0;
Int_t step = 5;
while(writtenBins < hist->GetNbinsX()){
printf("| Bin |");
Int_t lastToWrite = TMath::Min(hist->GetNbinsX(), writtenBins+step);
for(Int_t i = writtenBins+1; i <= lastToWrite; ++i){
printf("% 5.3g -% 5.3g |", hist->GetBinLowEdge(i), hist->GetBinLowEdge(i+1));
}
printf("n| Value|");
for(Int_t i = writtenBins+1; i <= lastToWrite; ++i){
printf(" % 10.3g |", hist->GetBinContent(i));
}
printf("n| Error|");
for(Int_t i = writtenBins+1; i <= lastToWrite; ++i){
printf(" % 10.3g |", hist->GetBinError(i));
}
printf("n%s", line);
writtenBins += step;
}
}
//_____________________________________________________
void GFHistManip::AddToHistsName(const char* nameAdd, GFHistArray* hists,
const char* title)
{
// namAdd is added to hist name, if title is given, adds it to title
if(!hists) return;
TString add(nameAdd);
TString titleAdd(title);
for(Int_t i = 0; i < hists->GetEntriesFast(); ++i){
if(!(*hists)[i]) continue;
(*hists)[i]->SetName((*hists)[i]->GetName()+add);
if(title) (*hists)[i]->SetTitle((*hists)[i]->GetTitle()+titleAdd);
}
}
//_____________________________________________________
void GFHistManip::CallSumw2(GFHistArray* hists)
{
TIter next(hists);
while(TH1* h = static_cast<TH1*>(next())){
h->Sumw2();
}
}
//_____________________________________________________
void GFHistManip::CreateAnalogHists(GFHistArray* input,
GFHistArray*& output,
const char* addName,
const char* addTitle)
{
// creates output if not yet existing
if(!output) ClearCreateArray(output);
for(Int_t i = 0; i <= input->GetLast(); ++i){
TH1* hInput = (*input)[i];
if(!hInput) continue;
TH1* hist = static_cast<TH1*>(hInput->Clone(TString((*input)[i]->GetName())+=addName));
TString title((*input)[i]->GetTitle());
(title += " ")+= (addTitle ? addTitle : addName);
hist->SetTitle(title);
output->AddAtAndExpand(hist, i);
}
}
//_____________________________________________________
TH1* GFHistManip::CreateAnalogHist(const TH1* hBins, const TGraphAsymmErrors* errors,
Bool_t upper)
{
// create a hist with borders bins from hBins, fill with errors->GetY() but add
// GetEYhigh() (if upper = kTRUE) or subtract GetEYlow()
if(!hBins || !errors || hBins->GetNbinsX() != errors->GetN()) return NULL;
// TString name(Form("%s%s", hBins->GetName(), (upper?"upper":"lower")));
TString name(upper? "upper" : "lower");
GFHistManip::MakeUniqueName(name);
TH1* result = CreateAnalog(hBins, name, kFALSE); // no sumw2 structure, no fill
for(Int_t i = 0; i < hBins->GetNbinsX(); ++i){
Double_t binContent = errors->GetY()[i];
if(upper) binContent += errors->GetEYhigh()[i];
else binContent -= errors->GetEYlow()[i];
result->SetBinContent(i+1, binContent);
if(result->GetSumw2N()) result->SetBinError(i+1, FLT_MIN); // just to get it drawn...
}
return result;
}
//_____________________________________________________
TH1* GFHistManip::CreateRelErrorHist(const TH1* hist, Bool_t abs)
{
// create a hist containing the relative errors of hist, if abs = kTRUE only absolute values
if(!hist) return NULL;
if(hist->GetDimension() != 1){
::Error("GFHistManip::CreateRelErrorHist", "not implemented for %d-D",hist->GetDimension());
}
TH1* hRelErr = GFHistManip::CreateAnalog(hist, Form("%sRelErr", hist->GetName()), kFALSE);
for(Int_t bin = 0; bin <= hist->GetNbinsX()+1; ++bin){// include over/underflow
const Stat_t cont = hist->GetBinContent(bin);
const Stat_t err = hist->GetBinError(bin);
if(cont) hRelErr->SetBinContent(bin, err/cont);
if(abs && hRelErr->GetBinContent(bin) < 0.){
hRelErr->SetBinContent(bin, -hRelErr->GetBinContent(bin));
}
}
return hRelErr;
}
//_____________________________________________________
TH1* GFHistManip::CreateAverage(GFHistArray* input, const TArrayD& weights,
const char* name)
{
if(!input || input->GetEntriesFast() != weights.GetSize() ||
input->GetEntriesFast() != input->GetEntries()){
::Error("GFHistManip::CreateAverage", "Problem with input");
return NULL;
}
const Double_t sumOfWeights = weights.GetSum();
if(!sumOfWeights){
::Error("GFHistManip::CreateAverage", "Sum of weights = 0.");
return NULL;
}
TH1* hAverage = static_cast<TH1*>(input->First()->Clone(name));
if(weights.GetSize() == 1) return hAverage;
else Scale(hAverage, weights[0]);
for(Int_t i = 1; i < weights.GetSize(); ++i){
Add2ndTo1st(hAverage, input->At(i), weights[i]);
}
Scale(hAverage, 1./sumOfWeights);
return hAverage;
}
//_____________________________________________________
TH1* GFHistManip::CreateAnalog(const TH1* h, const char* name, Bool_t sumw2,
Bool_t fill)
{
// create a hist with the same bin borders as 'h', setting 'name', copy title from 'h',
// if(sumw2) call Sumw2(), if(fill) fill content (and errors if sumw2) from 'h'
// currently only for 1D hists, TH1D and TH1F,
// axis titles will be copied
// do not create profiles...
if(!h) return NULL;
const Int_t type = h->InheritsFrom(TArrayD::Class()) ? 1 : 0;
if(type == 0 && !h->InheritsFrom(TArrayF::Class())){
::Warning("GFHistManip::CreateAnalog", "Creating TH<n>F not %s", h->IsA()->GetName());
}
TString title(Form("%s;%s;%s;%s",h->GetTitle(), h->GetXaxis()->GetTitle(),
h->GetYaxis()->GetTitle(), h->GetZaxis()->GetTitle()));
TH1* result = NULL;
if(h->GetDimension() == 1){
const Int_t nBins = h->GetNbinsX();
const TAxis* axis = h->GetXaxis();
if(axis->GetXbins()->GetSize()){ // non-equidistant bins
switch(type){
case 1:
result = new TH1D(name, title, nBins, axis->GetXbins()->GetArray());
break;
case 0: default:
result = new TH1F(name, title, nBins, axis->GetXbins()->GetArray());
}
} else { // equidistant bins
switch(type){
case 1:
result = new TH1D(name, title, nBins, axis->GetXmin(), axis->GetXmax());
break;
case 0: default:
result = new TH1F(name, title, nBins, axis->GetXmin(), axis->GetXmax());
}
}
if(sumw2) result->Sumw2();
if(fill){
for(Int_t i = 0; i <= nBins; ++i){ // including under-/overflow
result->SetBinContent(i, h->GetBinContent(i));
if(sumw2) result->SetBinError(i, h->GetBinError(i));
}
Stat_t s[11];
h->GetStats(s);
result->PutStats(s);
}
} else if(h->GetDimension() == 2) {
const Int_t nBinsX = h->GetNbinsX();
const Int_t nBinsY = h->GetNbinsY();
const TAxis* axisX = h->GetXaxis();
const TAxis* axisY = h->GetYaxis();
if (axisX->GetXbins()->GetSize()) {
if (axisY->GetXbins()->GetSize()) {// non-equidistant bins in x/y
switch(type){
case 1:
result = new TH2D(name, title, nBinsX, axisX->GetXbins()->GetArray(),
nBinsY, axisY->GetXbins()->GetArray());
break;
case 0: default:
result = new TH2F(name, title, nBinsX, axisX->GetXbins()->GetArray(),
nBinsY, axisY->GetXbins()->GetArray());
}
} else { // equidistant bins in X, but not Y
::Error("GFHistManip::CreateAnalog","no support for X but not Y equidist");
}
} else { // i.e. not equidistant in X
if (axisY->GetXbins()->GetSize()) { // equidistant bins in Y, but not X
::Error("GFHistManip::CreateAnalog","no support for Y but not X equidist");
} else { // equidistant bins in x/y
switch(type){
case 1:
result = new TH2D(name, title, nBinsX, axisX->GetXmin(), axisX->GetXmax(),
nBinsY, axisY->GetXmin(), axisY->GetXmax());
break;
case 0: default:
result = new TH2F(name, title, nBinsX, axisX->GetXmin(), axisX->GetXmax(),
nBinsY, axisY->GetXmin(), axisY->GetXmax());
}
}
}
if (result) {
if (sumw2) result->Sumw2();
if (fill) {
for (Int_t iX = 0; iX <= nBinsX; ++iX) { // including under-/overflow
for (Int_t iY = 0; iY <= nBinsY; ++iY) { // including under-/overflow
result->SetBinContent(iX, iY, h->GetBinContent(iX, iY));
if (sumw2) result->SetBinError(iX, iY, h->GetBinError(iX, iY));
}
}
Stat_t s[11];
h->GetStats(s);
result->PutStats(s);
}
}
} else {
::Error("GFHistManip::CreateAnalog", "No support for 3D (yet)");
}
return result;
}
//_____________________________________________________
GFHistArray* GFHistManip::ClearCreateArray(GFHistArray*& array)
{
// array is created as new GFHistArray which OWNS all hists that will be filled into it
// array is also returned...
// if array existed: content will be cleared!
// array is set to be owner of all hists in it!
if(array && array->GetEntriesFast()) {
::Warning("GFHistManip::ClearCreateArray",
"Array filled, call its Print and clear content!");
array->Print();
array->Clear();
} else if(!array) {
array = new GFHistArray;
}
array->SetOwner();
return array;
}
//_____________________________________________________
Double_t GFHistManip::Normalise(TH1* hist, Option_t* opt)
{
// scales histogram to be normalised,
// if opt == "width": taking into account the bin width
// returns scale factor
if(!hist) return 0.;
const Double_t integral = hist->Integral(opt);
if(integral) {
const Double_t scale = 1./ integral;
// hist->Scale(scale);
Scale(hist, scale);
if(hist->GetMaximumStored() != -1111.) hist->SetMaximum(hist->GetMaximum() * scale);
if(hist->GetMinimumStored() != -1111.) hist->SetMinimum(hist->GetMinimum() * scale);
return scale;
} else {
return 0.;
}
}
//____________________________________
Double_t GFHistManip::Normalise(TH1* hist, TObjArray* singleHs)
{
// normalise hist 'hist' and multiply all hists in 'singleHs' by the same factor used to
// normalise 'hist'
if(!hist) return 0.;
const Double_t factor = Normalise(hist);
TIter next(singleHs);
while(TObject* obj = next()){
if(!obj->InheritsFrom(TH1::Class())) continue;
TH1* hist = static_cast<TH1*>(obj);
Scale(hist, factor);
if(hist->GetMaximumStored() != -1111) hist->SetMaximum(hist->GetMaximum() * factor);
if(hist->GetMinimumStored() != -1111) hist->SetMinimum(hist->GetMinimum() * factor);
}
return factor;
}
//____________________________________
Double_t GFHistManip::Normalise1stTo2nd(TH1* hist1, TH1* hist2, Option_t* opt)
{
// if opt == "width" (default) integrals take into account the binwidth
const Double_t content = hist1->Integral(opt);
if(content){
const Double_t scale = hist2->Integral(opt) / content;
Scale(hist1, scale);
if(hist1->GetMaximumStored() != -1111.) hist1->SetMaximum(hist1->GetMaximum() * scale);
if(hist1->GetMinimumStored() != -1111.) hist1->SetMinimum(hist1->GetMinimum() * scale);
return scale;
}
return 1.;
}
//____________________________________
void GFHistManip::NormaliseBinsToSum(GFHistArray* hists)
{
// normalises the 'hists' such that for each bin the sum of all gives one,
// must be 1D and must have same number of bins
// errors will be correctly calculated assuming that the input histograms' errors
// are uncorrelated, but you'll not have access to the correlations of the
if(!hists || !hists->GetEntriesFast()) return;
// check to be 1D:
for(Int_t jHist = 0; jHist < hists->GetEntriesFast(); ++jHist){
const Int_t dim = hists->At(jHist)->GetDimension();
if(hists->At(jHist)->GetDimension() != 1) {
::Error("GFHistManip::NormaliseBinsToSum", "hist %d is %d-D, return", jHist, dim);
return;
}
}
if(hists->GetEntriesFast() == 1){ // 1 hist:
for(Int_t i = 1; i <= hists->At(0)->GetNbinsX(); ++i){
hists->At(0)->SetBinContent(i, 1.);
if(hists->At(0)->GetSumw2N()) hists->At(0)->SetBinError(i, 0.);
}
} else if(hists->GetEntriesFast() > 1){
for(Int_t iBin = 1; iBin <= hists->At(0)->GetNbinsX(); ++iBin){
// first get bin content and error this bin for all hists
TArrayD contents(hists->GetEntriesFast());
TArrayD errors(contents.GetSize());
for(Int_t jHist = 0; jHist < contents.GetSize(); ++jHist){
contents[jHist] = hists->At(jHist)->GetBinContent(iBin);
errors[jHist] = hists->At(jHist)->GetBinError(iBin);
}
// call normalisation method to get fractions and their errors
TArrayD fracErrors;
TMatrixD* covar = NULL; // forget about correlations...
TArrayD fractions = GFMath::NormaliseTo1(contents, &errors, &fracErrors, covar);
// fill fractions and errors to this bin of all hists
for(Int_t jHist = 0; jHist < contents.GetSize(); ++jHist){
hists->At(jHist)->SetBinContent(iBin, fractions[jHist]);
hists->At(jHist)->SetBinError(iBin, fracErrors[jHist]);
}
}
}
}
//________________________________________________
void GFHistManip::CorrectForBinWidth(TH1* hist, Bool_t correctError)
{
// divides all bins content and error through the binwidth
// for 1D only
// under-/overflow untouched (what else...?)
if (!hist) return;
if (hist->GetDimension() != 1) {
::Error("CorrectForBinWidth", "cannot deal with %dD hist", hist->GetDimension());
return;
}
for(Int_t i = 1; i <= hist->GetNbinsX(); ++i){
const Double_t binWidth = hist->GetBinWidth(i); // is always != 0
hist->SetBinContent(i, hist->GetBinContent(i) / binWidth);
if(correctError) hist->SetBinError (i, hist->GetBinError(i) / binWidth);
}
}
//________________________________________________
void GFHistManip::ReCorrectForBinWidth(TH1* hist, Bool_t correctError)
{
// multiply all bins' content and error with the binwidth
// for 1D only
// under-/overflow untouched (what else...?)
if (!hist) return;
if (hist->GetDimension() != 1) {
::Error("ReCorrectForBinWidth", "cannot deal with %dD hist", hist->GetDimension());
return;
}
for(Int_t i = 1; i <= hist->GetNbinsX(); ++i){
const Double_t binWidth = hist->GetBinWidth(i);
hist->SetBinContent(i, hist->GetBinContent(i) * binWidth);
if(correctError) hist->SetBinError (i, hist->GetBinError(i) * binWidth);
}
}
//__________________________________________________________________
TH1* GFHistManip::CreateRatioGassnerErrors(const TH1* hCount, const TH1* hDenom)
{
// cf. appendix E of Johannes Gassner's thesis (from 2002)
if(!hCount || !hDenom) return NULL;
if(hCount->GetDimension() != 1 || hDenom->GetDimension() != 1) {
::Error("GFHistManip::CreateRatioGassnerErrors", "currently works for 1-D hists only");
return NULL;
}
const Int_t nBins = hDenom->GetNbinsX();
if(nBins != hCount->GetNbinsX()) return NULL;
TString newName(hDenom->GetName());
GFHistManip::MakeUniqueName(newName+="RatioGassner");
TH1* hRatio = GFHistManip::CreateAnalog(hCount, newName, kTRUE);// sumw2 true
hRatio->Reset();
for(Int_t bin = 0; bin <= nBins+1; ++bin){
const Double_t nCount = hCount->GetBinContent(bin);
const Double_t errCount = hCount->GetBinError(bin);
const Double_t nDenom = hDenom->GetBinContent(bin);
const Double_t errDenom = hDenom->GetBinError(bin);
const TArrayD effErr = GFMath::EfficiencyGassnerError(nCount, nDenom, errCount, errDenom);
hRatio->SetBinContent(bin, effErr[0]);
hRatio->SetBinError (bin, effErr[1]);
// const Double_t eff = hRatio->GetBinContent(bin);
// Double_t relErrEff2 = 0.;
// if(nCount && nDenom) relErrEff2 = errDenom * errDenom / nDenom / nDenom
// + (1. - 2.*eff) * errCount * errCount / nCount / nCount;
// hRatio->SetBinError(bin, TMath::Sqrt(relErrEff2) * eff);
}
return hRatio;
}
//__________________________________________________________________
TH1* GFHistManip::CreateRatioBlobelErrors(const TH1* hCount, const TH1* hDenom)
{
if(!hCount || !hDenom) return NULL;
if(hCount->GetDimension() != 1 || hDenom->GetDimension() != 1) {
::Error("GFHistManip::CreateRatioBlobelErrors", "currently works for 1-D hists only");
return NULL;
}
const Int_t nBins = hDenom->GetNbinsX();
if(nBins != hCount->GetNbinsX()) return NULL;
TString newName(hDenom->GetName());
GFHistManip::MakeUniqueName(newName+="RatioBlobel");
TH1* hRatio = GFHistManip::CreateAnalog(hCount, newName, kTRUE);// sumw2 true
hRatio->Reset();
for(Int_t bin = 0; bin <= nBins+1; ++bin){// include overflow/underflow
const Double_t nCount = hCount->GetBinContent(bin);
const Double_t errCount = hCount->GetBinError(bin);
const Double_t nDenom = hDenom->GetBinContent(bin);
const Double_t errDenom = hDenom->GetBinError(bin);
const TArrayD effErr = GFMath::EfficiencyBlobelError(nCount, nDenom, errCount, errDenom);
hRatio->SetBinContent(bin, effErr[0]);
hRatio->SetBinError (bin, effErr[1]);
}
return hRatio;
}
//__________________________________________________________________
TH1* GFHistManip::CreateRatioBinomErrors(const TH1* hCount, const TH1* hDenom,
TArrayD* upError, TArrayD* lowError,
const TH1* hNoWeightDenom,
TGraphAsymmErrors** graphPtrPtr)
{
// Create a hist hCount/hDenom and set errors as mean of lower and upper error from
// binomial distribution via GFMath::BinomialError.
// If upError or lowError are given, they are filled with the correct errors
// (including over/underflow bins, so lowError[1] is for bin 1, not bin 0).
// If graphPtrPtr is given, it is dereferenced and filled with a pointer to a newly
// created TGraphAysmmErrors that contains the correct asymmetric errors.
//
// If the histograms are weighted, the number of entries (= sum of weights) in each bin
// is scaled to the number of effective entries in the denominater:
//
// N_eff = (sum w_i)^2/sum(w_i^2) (cf. LOOK manual)
//
// If N_eff/N < 0.66 (N number of 'filling actions') this method may not be applicable. So
// if(hNoWeightDenom)
// check for each bin N_eff/hNoWeightDenom->GetBinContent
// else
// check the quantity of the total hDenom (from GetStats and GetEntries())
//
// If any of the checked rations is below 0.66: return NULL
// if any of the given hists
// first lots of checks...
if(!hCount || !hDenom) return NULL;
if(hCount->GetDimension() != 1 || hDenom->GetDimension() != 1) {
::Error("GFHistManip::CreateRatioBinomErrors", "currently works for 1-D hists only");
return NULL;
}
const Int_t nBins = hDenom->GetNbinsX(); // FIXME: all bins!
if(nBins != hCount->GetNbinsX()) return NULL;
const Double_t minNeffRatio = 0.66;
if(!hNoWeightDenom){// in case no detailed check possible, do it roughly overall:
if(!hDenom->GetEntries()) return NULL;
Stat_t s[11];
hDenom->GetStats(s);//now s[1] is sum of squares of weights, s[0] is sum of weights
if(s[0]*s[0]/s[1]/hDenom->GetEntries() < minNeffRatio){
::Error("GFHistManip::CreateRatioBinomErrors", "N_eff %.4g, N %.4g", s[0]*s[0]/s[1],
hDenom->GetEntries());
return NULL;
}
} else if(hNoWeightDenom->GetDimension() != 1 || nBins != hNoWeightDenom->GetNbinsX()){
if(hNoWeightDenom->GetDimension() != 1)
::Error("GFHistManip::CreateRatioBinomErrors", "currently works for 1-D hists only");
return NULL;
}
TH1* hRatio = CreateAnalog(hCount, Form("%sBy%s",hCount->GetName(),hDenom->GetName()), kTRUE);
TArrayD upErr(nBins+2);
TArrayD lowErr(nBins+2);
for(Int_t i = 0; i <= nBins+1; ++i){//include over/underflow
const Double_t count = hCount->GetBinContent(i);
const Double_t denom = hDenom->GetBinContent(i); // sum of weights
const Double_t sqrtDenomSumW2 = hDenom->GetBinError(i); // sqrt(sum of squares of weights)
if(!denom) continue; // no ratio...
if(sqrtDenomSumW2) {
const Double_t nEff = denom*denom/(sqrtDenomSumW2*sqrtDenomSumW2);
if(hNoWeightDenom && nEff/hNoWeightDenom->GetBinContent(i) < minNeffRatio) {
// detailed check of validity requested but failed
::Error("GFHistManip::CreateRatioBinomErrors", "bin %d: N_eff %.4g, N %.4g",
i, nEff, hNoWeightDenom->GetBinContent(i));
delete hRatio;
return NULL;
}
upErr[i]= GFMath::BinomialError(TMath::Nint(count*nEff/denom),TMath::Nint(nEff), kTRUE);
lowErr[i]=GFMath::BinomialError(TMath::Nint(count*nEff/denom),TMath::Nint(nEff),kFALSE);
hRatio->SetBinError(i, (upErr[i]+lowErr[i])/2.);
} else {
hRatio->SetBinError(i, 0.);// to force creation of sumw2 structure, to not get sqrt(ratio)
}
hRatio->SetBinContent(i, count/denom);
}
if(upError) *upError = upErr;
if(lowError) *lowError = lowErr;
if(graphPtrPtr){
TArrayD binCenters(nBins);
TArrayD binValues(nBins);
for(Int_t i = 0; i < nBins; ++i){
binCenters[i] = hRatio->GetXaxis()->GetBinCenter(i+1);
binValues[i] = hRatio->GetBinContent(i+1);
}
*graphPtrPtr = new TGraphAsymmErrors(nBins, binCenters.GetArray(),
binValues.GetArray(), NULL, NULL, // nothing in X
lowErr.GetArray()+1, upErr.GetArray()+1);//underflow +1
}
return hRatio;
}
//____________________________________
TH1* GFHistManip::CreateEff(const GFHistArray* histsCount,const GFHistArray* histsDenom,
const TArrayD& lumi, TGraphAsymmErrors** graphPtrPtr)
{
// add histsCount and histsDenom lumi weighted and create efficiency hist via
// CreateRatioBinomErrors (so check for NULL!),
// graphPtrPtr is also passed to that routine.
if(!histsCount || ! histsDenom || histsCount->GetEntriesFast() != histsDenom->GetEntriesFast()
|| histsDenom->GetEntriesFast() != lumi.GetSize()){
return NULL;
}
TH1* hCountAll = NULL;
TH1* hDenomAll = NULL;
Double_t sumw2 = 0.;
Double_t entries = 0.;
Stat_t s[11];
for(Int_t i = 0; i < lumi.GetSize(); ++i){
if(!lumi[i]) continue; // skip if no lumi
TH1* hCount = histsCount->At(i);
if(!hCountAll){
hCountAll = CreateAnalog(hCount, Form("%sMerge", hCount->GetName()), kTRUE, kTRUE);
Scale(hCountAll, 1./lumi[i]);
} else {
GFHistManip::Add2ndTo1st(hCountAll, hCount, 1./lumi[i]);
}
// hCountAll = static_cast<TH1*>(hCount->Clone());
// hCountAll->Reset();
// }
// hCountAll->Add(hCount, 1./lumi[i]);
TH1* hDenom = histsDenom->At(i);
hDenom->GetStats(s);
sumw2 += (s[1]/(lumi[i]*lumi[i]));
entries += hDenom->GetEntries();
if(!hDenomAll){
hDenomAll = CreateAnalog(hDenom, Form("%sMerge", hDenom->GetName()), kTRUE, kTRUE);
Scale(hDenomAll, 1./lumi[i]);
} else {
GFHistManip::Add2ndTo1st(hDenomAll, hDenom, 1./lumi[i]);
}
// hDenomAll = static_cast<TH1*>(hDenom->Clone());
// hDenomAll->Reset();
// }
// hDenomAll->Add(hDenom, 1./lumi[i]); // GFHistManip::Add2ndTo1st(...) ?
}
hDenomAll->GetStats(s);
s[1] = sumw2;
hDenomAll->PutStats(s);
hDenomAll->SetEntries(entries);
TH1* eff = CreateRatioBinomErrors(hCountAll, hDenomAll, NULL, NULL, NULL, graphPtrPtr);
delete hDenomAll;
delete hCountAll;
return eff;
}
//__________________________________________________________________
Bool_t GFHistManip::Scale(TH1* hist, Double_t factor)
{
// due to a bug in ROOT (fixed in 4.00_0?), the histogram's sum of squares of weights
// is not properly adjusted when scaling a histogram
if(!hist) return kFALSE;
#if ROOT_VERSION_CODE < ROOT_VERSION(4,2,0)
Stat_t s[11];
hist->GetStats(s);
const Stat_t sumW2 = s[1] * factor * factor;
hist->Scale(factor);
hist->GetStats(s);
s[1] = sumW2;
hist->PutStats(s);
#else
hist->Scale(factor);
#endif
return kTRUE;
}
//__________________________________________________________________
Bool_t GFHistManip::Add2ndTo1st(TH1* h1, const TH1* h2, Double_t factor)
{
// similar to GFHistManip::Scale
if(!h1 || !h2) return kFALSE;
#if ROOT_VERSION_CODE < ROOT_VERSION(5,4,0)
if (factor != 1. && h1->InheritsFrom(TProfile::Class())) {
::Warning("Add2ndTo1st", "use ROOT >= 5.04.00 to add TProfile with factor != 1");
}
#endif
#if ROOT_VERSION_CODE < ROOT_VERSION(4,2,0)
Stat_t s1[11], s2[11];
h1->GetStats(s1);
h2->GetStats(s2);
const Stat_t sumW2 = s1[1] + s2[1] * factor * factor;
h1->Add(h2, factor);
h1->GetStats(s1);
s1[1] = sumW2;
h1->PutStats(s1);
#else
h1->Add(h2, factor);
#endif
return kTRUE;
}
//__________________________________________________________________
Bool_t GFHistManip::AddProfiles(TProfile *p1, const TProfile *p2,
Double_t f1, Double_t f2)
{
// there is a bug in TProfile::Add(...) which makes the error bars wrong,
// fixed by the ROOT team (Rene Brun) on August 9th, 2005
// this fix is not yet tested or does not seem to work...
::Error("AddProfiles", "method does not seem to work as expected!");
if(!p1 || !p2 || p1->GetNbinsX() != p2->GetNbinsX()) return kFALSE;
#if ROOT_VERSION_CODE < ROOT_VERSION(5,4,0)
Stat_t s1[11], s2[11];
p1->GetStats(s1);
p2->GetStats(s2);
const Stat_t sumW2 = s1[1]*f1*f1 + s2[1]*f2*f2;
const Int_t nBins = p1->GetNbinsX() + 2;// incl. under- and overflow
TArrayD sumwy2(nBins);
for(Int_t i = 0; i < nBins; ++i){
const Double_t e1 = p1->TH1::GetBinError(i);// returns sqrt(fSumw2[i])
const Double_t e2 = p2->TH1::GetBinError(i);
cout << i << "con1 " << p1->GetBinContent(i) << " e1 " << e1 << " " << p1->GetBinError(i)
<< "con2 " << p2->GetBinContent(i) << " e2 " << e2 << " " << p2->GetBinError(i)
<< endl;
sumwy2[i] = f1*e1*e1 + f2*e2*e2;// TProfile::fSumw2 is used to store weight*y^2
}
p1->Add(p1, p2, f1, f2);
p1->GetStats(s1);
s1[1] = sumW2;
p1->PutStats(s1);
for(Int_t i = 0; i < nBins; ++i){
p1->SetBinError(i, TMath::Sqrt(sumwy2[i]));
}
#else
p1->Add(p1, p2, f1, f2);
#endif
return kTRUE;
}
//__________________________________
TH1* GFHistManip::MergeHists(const char *name, const TObjArray *dirs)
{
// dirs is assumed to be an array of TDirectory
// return a hist which is the sum of all hists with 'name'
TH1 *hResult = NULL;
TIter dirIter(dirs);
while (TDirectory *dir = static_cast<TDirectory*>(dirIter())) {
TH1 *h = NULL;
#if ROOT_VERSION_CODE < ROOT_VERSION(4,0,0)
h = static_cast<TH1*>(dir->Get(name));
if (!h && !h->InheritsFrom(TH1::Class())) {
delete h;
h = NULL;
}
#else
dir->GetObject(name, h);
#endif
if (!h) {
::Error("GFHistManip::MergeHists","no TH1 %s in directory %s",
name, dir->GetName());
continue;
}
if (hResult) {
GFHistManip::Add2ndTo1st(hResult, h, 1.);
} else {
hResult = static_cast<TH1*>(h->Clone());
}
delete h;
}
return hResult;
}
//__________________________________________________________________
void GFHistManip::RepairStat1(TH1* hist)
{
// ROOT versions before 4.01_xx calculate a wrong sumw2 if hist is result of a projection...
if(!hist) return;
if(hist->GetDimension() > 1){
::Warning("GFHistManip::RepairStat1", "Doesn't work for %d-D hists", hist->GetDimension());
} else {
Stat_t s[11];
hist->GetStats(s);
s[1] = .0;
for (Int_t binx = hist->GetXaxis()->GetFirst(); binx <= hist->GetXaxis()->GetLast();++binx){
Stat_t err = TMath::Abs(hist->GetBinError(binx));
s[1] += err*err;
}
hist->PutStats(s);
}
}
//__________________________________________________________________
Bool_t GFHistManip::Rebin(const TH1 *input, TH1 *output, Option_t *option)
{
// tries to fill the content of the input histogram into the output
// if possible. If not (inconsistent boundaries) returns false
// if option contains 'e', errors are calculated
// If false is returned, output may be only partially manipulated.
if(!input || !output) {
::Error("GFHistManip::Rebin", "input %p or output %p missing", input, output);
return kFALSE;
}
// check range consistency
if(input->GetXaxis()->GetXmin() > output->GetXaxis()->GetXmin()
&& !GFMath::IsEqualFloat(input->GetXaxis()->GetXmin(), output->GetXaxis()->GetXmin())){
::Error("GFHistManip::Rebin", "input/output minimum inconsistent");
return kFALSE;
}
if(input->GetXaxis()->GetXmax() < output->GetXaxis()->GetXmax()
&& !GFMath::IsEqualFloat(input->GetXaxis()->GetXmax(), output->GetXaxis()->GetXmax())){
::Error("GFHistManip::Rebin", "input/output maximum inconsistent");
return kFALSE;
}
const Bool_t err = TString(option).Contains("E",TString::kIgnoreCase);
// find first bin in input
Int_t inBin = 1;
while(!GFMath::IsEqualFloat(input->GetBinLowEdge(inBin), output->GetBinLowEdge(1))){
++inBin;
if(inBin > input->GetNbinsX()){ // if no more bins...
::Error("GFHistManip::Rebin", "No input bin starts with %f", output->GetBinLowEdge(1));
return kFALSE;
}
}
Int_t firstInBin = inBin;
for(Int_t outBin = 1; outBin <= output->GetNbinsX(); ++outBin){
// sum up input as long as in same output bin
Double_t outEntry = 0.;
Double_t outErr2 = 0.;
while(!GFMath::IsEqualFloat(input->GetBinLowEdge(inBin), output->GetBinLowEdge(outBin + 1))
&& input->GetBinLowEdge(inBin) < output->GetBinLowEdge(outBin + 1)){
outEntry += input->GetBinContent(inBin);
if(err) outErr2 += input->GetBinError(inBin)*input->GetBinError(inBin);
++inBin;
}
// next bin borders consistent?
if(!GFMath::IsEqualFloat(input->GetBinLowEdge(inBin), output->GetBinLowEdge(outBin + 1))){
::Error("GFHistManip::Rebin",
"input (%s) bin %d starts at %f, output (%s) bin %d at %f, diff %f",
input->GetName(), inBin, input->GetBinLowEdge(inBin),
output->GetName(), outBin+1, output->GetBinLowEdge(outBin+1),
input->GetBinLowEdge(inBin) - output->GetBinLowEdge(outBin+1));
return kFALSE;
}
// set sum in output
output->SetBinContent(outBin, outEntry);
if(err) output->SetBinError(outBin, TMath::Sqrt(outErr2));
}
// last check: are the hist integrals the same?
if(!GFMath::IsEqualFloat(output->Integral(), input->Integral(firstInBin, inBin))){
::Error("GFHistManip::Rebin", "sums %s (out) %.3f, %s (in) %.3f, %s all %.3f (ignore)",
output->GetName(), output->Integral(), input->GetName(),
input->Integral(firstInBin, inBin), input->GetName(), input->Integral());
// return kFALSE;
}
return kTRUE;
}
//__________________________________________________________________
Double_t GFHistManip::MinMaxOfHist(const TH1* h, Int_t minMax)
{
// max or min of hist 'h' (max if minMax > 0, min if < 0, else result is 0.)
// If hist's option contains 'e' or it has weights stored and not option HIST,
// error bars are taken into account.
if(minMax == 0) return 0.;
TString option = h->GetOption();
option.ToLower();
option.ReplaceAll("same", 0);
const Int_t exBin = (minMax > 0 ? h->GetMaximumBin() : h->GetMinimumBin());
Double_t result = h->GetBinContent(exBin);
if(option.Contains('e') || (h->GetSumw2N() && !option.Contains("hist"))){
for(Int_t bin = 1; bin <= h->GetNbinsX(); ++bin){ //FIXME: for 2/3D: loop over y/z!
if(minMax > 0) result = TMath::Max(result,(h->GetBinContent(bin) + h->GetBinError(bin)));
else result = TMath::Min(result,(h->GetBinContent(bin) - h->GetBinError(bin)));
}
}
return result;
}
//__________________________________________________________________
Double_t GFHistManip::MinMaxOfHists(const GFHistArray* hists, Int_t minMax)
{
Double_t result = (minMax > 0 ? -DBL_MAX : DBL_MAX);
TIter nextHist(hists);
while(TH1* hist = static_cast<TH1*>(nextHist())){
if(minMax > 0) result = TMath::Max(result, GFHistManip::MinMaxOfHist(hist, minMax));
else result = TMath::Min(result, GFHistManip::MinMaxOfHist(hist, minMax));
}
return result;
}
//__________________________________________________________________
Bool_t GFHistManip::IsSameBinning(const TH1 *h1, const TH1 *h2)
{
if(!h1 || h1->GetDimension() > 1 || !h2 || h2->GetDimension() > 1) return kFALSE;
const Int_t nBins = h1->GetNbinsX();
if(nBins != h2->GetNbinsX()) return kFALSE;
for(Int_t i = 1; i < nBins; ++i){ // till overflow!
if(!GFMath::IsEqualFloat(h1->GetBinLowEdge(i), h2->GetBinLowEdge(i))) return kFALSE;
}
return kTRUE;
}
//__________________________________________________________________
//__________________________________________________________________
//__________________________________________________________________
Bool_t GFHistManip::CheckContent(TCollection *hists, const TClass *cl)
{
// true if all inherit from type cl
TIter iter(hists);
while(TObject* h = iter.Next()){
if(!h->InheritsFrom(cl)) return kFALSE;
}
return kTRUE;
}
//__________________________________________________________________
Bool_t GFHistManip::BinsAreSubRange(const TAxis *axis, Int_t firstBin, Int_t lastBin)
{
if(!axis || firstBin > lastBin || lastBin > axis->GetNbins() || firstBin < 1){
// could be extended to tolerate under/overflow?
return kFALSE;
}
return kTRUE;
}
//__________________________________________________________________
void GFHistManip::MakeUniqueName(TString& name)
{
// remove all special letters
// if 'name' is found by gROOT->FindObject(name), '_<unique_number>' is added
name.ReplaceAll("#", 0);
name.ReplaceAll("{", 0);
name.ReplaceAll("}", 0);
name.ReplaceAll("(", 0);
name.ReplaceAll(")", 0);
name.ReplaceAll("[", 0);
name.ReplaceAll("]", 0);
name.ReplaceAll("*", 0);
name.ReplaceAll("/", 0);
name.ReplaceAll(".", 0);
name.ReplaceAll(",", 0);
name.ReplaceAll(";", 0);
name.ReplaceAll(":", 0);
name.ReplaceAll("<", 0);
name.ReplaceAll(">", 0);
name.ReplaceAll("?", 0);
name.ReplaceAll("+", 0);
name.ReplaceAll("-", 0);
name.ReplaceAll("^", 0);
name.ReplaceAll("!", 0);
name.ReplaceAll(" ", 0);
if(gROOT->FindObject(name)){
name += '_';
Long_t count = 1;
while(gROOT->FindObject(name + count)){
++count;
}
name += count;
}
}
//__________________________________________________________________
Bool_t GFHistManip::RadToDegXaxis(TH1 *&hist, Bool_t perBinWidth)
{
if(!hist) return kFALSE;
if(hist->GetDimension() != 1) {
::Error("GFHistManip::RadToDegXaxis", "%d-D not yet foreseen, bit %s is!",
hist->GetDimension(), hist->GetName());
return kFALSE;
}
Axis_t *bins = new Axis_t[hist->GetNbinsX() + 1];
for(Int_t i = 0; i <= hist->GetNbinsX(); ++i){
bins[i] = hist->GetBinLowEdge(i+1) * TMath::RadToDeg();
}
TDirectory *oldHistsDir = hist->GetDirectory();
hist->SetDirectory(NULL);
TString tit(Form("%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(),
hist->GetYaxis()->GetTitle()));
TH1 *hnew = new TH1F(hist->GetName(), tit, hist->GetNbinsX(), bins);
for(Int_t i = 1; i <= hist->GetNbinsX(); ++i){
Double_t corrFac = (perBinWidth ? TMath::Pi()/180.: 1.);
hnew->SetBinContent(i, hist->GetBinContent(i) * corrFac);
if(hist->GetSumw2N()) hnew->SetBinError(i, hist->GetBinError(i) * corrFac);
}
hnew->SetDirectory(oldHistsDir);
delete hist;
delete [] bins;
hist = hnew;
return kTRUE;
}
//__________________________________________________________________
Bool_t GFHistManip::MergeBins(TH1 *&hist, Int_t bin1, Int_t bin2, Option_t *option)
{
// replacing 'hist' by a histogram where bins 'bin1' to 'bin2' are merged
// If they are both negative, they are counted backwards from the last
// bin (merge last two bins is (-1, -2)
// options: 'width': if hist is 'per bin width'
// 'e' : if the new hist should have calculated errors
if(!hist) return kFALSE;
if(hist->GetDimension() != 1) {
::Error("GFHistManip::MergeBins", "%d-D not yet foreseen, but %s is!",
hist->GetDimension(), hist->GetName());
return kFALSE;
}
if(bin1 * bin2 <= 0) {
::Error("GFHistManip::MergeBins", "both bins must be + or -, but are %d %d", bin1, bin2);
return kFALSE;
}
const Int_t firstBin = bin1 > 0 ? bin1 : bin2;
const Int_t lastBin = bin1 > 0 ? bin2 : bin1;
if(firstBin >= lastBin){
::Error("GFHistManip::MergeBins", "both to be merged %d %d", firstBin, firstBin);
return kFALSE;
}
TString opt(option);
opt.ToLower();
const Bool_t perBinWidth = opt.Contains("width");
const Bool_t errors = opt.Contains("e");
const Int_t nNewBins = hist->GetNbinsX() - (lastBin - firstBin);
Axis_t *bins = new Axis_t[nNewBins + 1];
Int_t diff = 0;
for(Int_t i = 1; i <= hist->GetNbinsX()+1; ++i){
if(i <= firstBin || i > lastBin) bins[i-diff-1] = hist->GetBinLowEdge(i);
else ++diff;
}
TDirectory *oldHistsDir = hist->GetDirectory();
hist->SetDirectory(NULL);
TString tit(Form("%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(),
hist->GetYaxis()->GetTitle()));
TH1 *hNew = new TH1F(hist->GetName(), tit, nNewBins, bins);
Double_t sum = 0., sumW2 = 0.;
diff = 0;
for(Int_t i = 1; i <= hist->GetNbinsX(); ++i){
if(i < firstBin || i > lastBin) {
hNew->SetBinContent(i-diff, hist->GetBinContent(i));
if(errors) hNew->SetBinError(i-diff, hist->GetBinError(i));
} else {
const Double_t corrFac = (perBinWidth ? hist->GetBinWidth(i): 1.);
sum += hist->GetBinContent(i) * corrFac;
if(errors) sumW2 += hist->GetBinError(i) * corrFac;
++diff;
}
if(i == lastBin){
Double_t width = 1.;
if(perBinWidth) width = hist->GetBinLowEdge(lastBin+1) - hist->GetBinLowEdge(firstBin);
hNew->SetBinContent(i-diff+1, sum/width);
if(errors) hNew->SetBinError(i-diff+1, TMath::Sqrt(sumW2)/width);
--diff;
}
}
hNew->SetDirectory(oldHistsDir);
delete hist;
delete [] bins;
hist = hNew;
return kTRUE;
}
//____________________________________
void GFHistManip::MergeOverToUnderFlow(TH1 *hist)
{
if (!hist) return;
if (hist->GetDimension() != 1) {
::Error("MergeOverToUnderFlow", "cannot deal with %dD hist", hist->GetDimension());
return;
}
const Int_t overFlowBin = hist->GetNbinsX() + 1;
hist->SetBinContent(0, hist->GetBinContent(0) + hist->GetBinContent(overFlowBin));
hist->SetBinContent(overFlowBin, 0.);
}
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.