#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <ctype.h>
#include "Riostream.h"
#include "TROOT.h"
#include "TClass.h"
#include "TH1.h"
#include "TH2.h"
#include "TF2.h"
#include "TF3.h"
#include "TPluginManager.h"
#include "TVirtualPad.h"
#include "TRandom.h"
#include "TVirtualFitter.h"
#include "THLimitsFinder.h"
#include "TProfile.h"
#include "TStyle.h"
#include "TVectorF.h"
#include "TVectorD.h"
#include "TBrowser.h"
#include "TObjString.h"
#include "TError.h"
#include "TVirtualFFT.h"
//Begin_Html
/*
<img src="gif/th1_classtree.gif">
*/
//End_Html
TF1 *gF1=0;
const Int_t kNstat = 11;
Int_t TH1::fgBufferSize = 1000;
Bool_t TH1::fgAddDirectory = kTRUE;
Bool_t TH1::fgDefaultSumw2 = kFALSE;
Bool_t TH1::fgStatOverflows= kFALSE;
extern void H1InitGaus();
extern void H1InitExpo();
extern void H1InitPolynom();
extern void H1LeastSquareFit(Int_t n, Int_t m, Double_t *a);
extern void H1LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail);
extern void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b);
ClassImp(TH1)
TH1::TH1(): TNamed(), TAttLine(), TAttFill(), TAttMarker()
{
fDirectory = 0;
fFunctions = new TList;
fNcells = 0;
fIntegral = 0;
fPainter = 0;
fEntries = 0;
fNormFactor = 0;
fTsumw = fTsumw2=fTsumwx=fTsumwx2=0;
fMaximum = -1111;
fMinimum = -1111;
fBufferSize = 0;
fBuffer = 0;
fXaxis.SetName("xaxis");
fYaxis.SetName("yaxis");
fZaxis.SetName("zaxis");
fXaxis.SetParent(this);
fYaxis.SetParent(this);
fZaxis.SetParent(this);
UseCurrentStyle();
}
TH1::~TH1()
{
if (!TestBit(kNotDeleted)) {
return;
}
delete[] fIntegral;
fIntegral = 0;
delete[] fBuffer;
fBuffer = 0;
if (fFunctions) {
fFunctions->SetBit(kInvalidObject);
TObject* obj = 0;
while ((obj = fFunctions->First())) {
while(fFunctions->Remove(obj));
if (!obj->TestBit(kNotDeleted)) {
break;
}
delete obj;
obj = 0;
}
delete fFunctions;
fFunctions = 0;
}
if (fDirectory) {
fDirectory->GetList()->Remove(this);
fDirectory = 0;
}
delete fPainter;
fPainter = 0;
}
TH1::TH1(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
:TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
Build();
if (nbins <= 0) nbins = 1;
fXaxis.Set(nbins,xlow,xup);
fNcells = fXaxis.GetNbins()+2;
}
TH1::TH1(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
:TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
Build();
if (nbins <= 0) nbins = 1;
if (xbins) fXaxis.Set(nbins,xbins);
else fXaxis.Set(nbins,0,1);
fNcells = fXaxis.GetNbins()+2;
if (fgDefaultSumw2) Sumw2();
}
TH1::TH1(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
:TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
Build();
if (nbins <= 0) nbins = 1;
if (xbins) fXaxis.Set(nbins,xbins);
else fXaxis.Set(nbins,0,1);
fNcells = fXaxis.GetNbins()+2;
if (fgDefaultSumw2) Sumw2();
}
TH1::TH1(const TH1 &h) : TNamed(), TAttLine(), TAttFill(), TAttMarker()
{
Copy((TObject&)h);
}
Bool_t TH1::AddDirectoryStatus()
{
return fgAddDirectory;
}
void TH1::Browse(TBrowser *b)
{
Draw(b ? b->GetDrawOption() : "");
gPad->Update();
}
void TH1::Build()
{
fDirectory = 0;
fPainter = 0;
fIntegral = 0;
fEntries = 0;
fNormFactor = 0;
fTsumw = fTsumw2=fTsumwx=fTsumwx2=0;
fMaximum = -1111;
fMinimum = -1111;
fBufferSize = 0;
fBuffer = 0;
fXaxis.SetName("xaxis");
fYaxis.SetName("yaxis");
fZaxis.SetName("zaxis");
fYaxis.Set(1,0.,1.);
fZaxis.Set(1,0.,1.);
fXaxis.SetParent(this);
fYaxis.SetParent(this);
fZaxis.SetParent(this);
SetTitle(fTitle.Data());
fFunctions = new TList;
UseCurrentStyle();
if (fgAddDirectory && gDirectory) {
if (!gDirectory->GetList()) {
Warning("Build","Current directory is not a valid directory");
return;
}
TH1 *hold = (TH1*)gDirectory->GetList()->FindObject(GetName());
if (hold) {
Warning("Build","Replacing existing histogram: %s (Potential memory leak).",GetName());
gDirectory->GetList()->Remove(hold);
hold->SetDirectory(0);
}
gDirectory->Append(this);
fDirectory = gDirectory;
}
}
void TH1::Add(TF1 *f1, Double_t c1, Option_t *option)
{
if (!f1) {
Error("Add","Attempt to add a non-existing function");
return;
}
TString opt = option;
opt.ToLower();
Bool_t integral = kFALSE;
if (opt.Contains("i") && fDimension ==1) integral = kTRUE;
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
Double_t s1[10];
Int_t i;
for (i=0;i<10;i++) {s1[i] = 0;}
PutStats(s1);
SetMinimum();
SetMaximum();
Int_t bin, binx, biny, binz;
Double_t cu=0;
Double_t xx[3];
Double_t *params = 0;
f1->InitArgs(xx,params);
for (binz=0;binz<=nbinsz+1;binz++) {
xx[2] = fZaxis.GetBinCenter(binz);
for (biny=0;biny<=nbinsy+1;biny++) {
xx[1] = fYaxis.GetBinCenter(biny);
for (binx=0;binx<=nbinsx+1;binx++) {
xx[0] = fXaxis.GetBinCenter(binx);
if (!f1->IsInside(xx)) continue;
TF1::RejectPoint(kFALSE);
bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
if (integral) {
xx[0] = fXaxis.GetBinLowEdge(binx);
cu = c1*f1->EvalPar(xx);
cu += c1*f1->Integral(fXaxis.GetBinLowEdge(binx),fXaxis.GetBinUpEdge(binx))*fXaxis.GetBinWidth(binx);
} else {
cu = c1*f1->EvalPar(xx);
}
if (TF1::RejectedPoint()) continue;
Double_t error1 = GetBinError(bin);
AddBinContent(bin,cu);
if (fSumw2.fN) {
fSumw2.fArray[bin] = error1*error1;
}
}
}
}
}
void TH1::Add(const TH1 *h1, Double_t c1)
{
if (!h1) {
Error("Add","Attempt to add a non-existing histogram");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()) {
Error("Add","Attempt to add histograms with different number of bins");
return;
}
if (fXaxis.GetXmin() != h1->fXaxis.GetXmin() ||
fXaxis.GetXmax() != h1->fXaxis.GetXmax() ||
fYaxis.GetXmin() != h1->fYaxis.GetXmin() ||
fYaxis.GetXmax() != h1->fYaxis.GetXmax() ||
fZaxis.GetXmin() != h1->fZaxis.GetXmin() ||
fZaxis.GetXmax() != h1->fZaxis.GetXmax()) {
Warning("Add","Attempt to add histograms with different axis limits");
}
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
fEntries += c1*h1->GetEntries();
Double_t s1[kNstat], s2[kNstat];
Int_t i;
for (i=0;i<kNstat;i++) {s1[i] = s2[i] = 0;}
GetStats(s1);
h1->GetStats(s2);
for (i=0;i<kNstat;i++) {
if (i == 1) s1[i] += c1*c1*s2[i];
else s1[i] += TMath::Abs(c1)*s2[i];
}
PutStats(s1);
SetMinimum();
SetMaximum();
Int_t bin, binx, biny, binz;
Double_t cu;
Double_t factor =1;
if (h1->GetNormFactor() != 0) factor = h1->GetNormFactor()/h1->GetSumOfWeights();;
for (binz=0;binz<=nbinsz+1;binz++) {
for (biny=0;biny<=nbinsy+1;biny++) {
for (binx=0;binx<=nbinsx+1;binx++) {
bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
cu = c1*factor*h1->GetBinContent(bin);
AddBinContent(bin,cu);
if (fSumw2.fN) {
Double_t error1 = factor*h1->GetBinError(bin);
fSumw2.fArray[bin] += c1*c1*error1*error1;
}
}
}
}
}
void TH1::Add(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2)
{
if (!h1 || !h2) {
Error("Add","Attempt to add a non-existing histogram");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()
|| nbinsx != h2->GetNbinsX() || nbinsy != h2->GetNbinsY() || nbinsz != h2->GetNbinsZ()) {
Error("Add","Attempt to add histograms with different number of bins");
return;
}
if (fXaxis.GetXmin() != h1->fXaxis.GetXmin() ||
fXaxis.GetXmax() != h1->fXaxis.GetXmax() ||
fYaxis.GetXmin() != h1->fYaxis.GetXmin() ||
fYaxis.GetXmax() != h1->fYaxis.GetXmax() ||
fZaxis.GetXmin() != h1->fZaxis.GetXmin() ||
fZaxis.GetXmax() != h1->fZaxis.GetXmax()) {
Warning("Add","Attempt to add histograms with different axis limits");
}
if (fXaxis.GetXmin() != h2->fXaxis.GetXmin() ||
fXaxis.GetXmax() != h2->fXaxis.GetXmax() ||
fYaxis.GetXmin() != h2->fYaxis.GetXmin() ||
fYaxis.GetXmax() != h2->fYaxis.GetXmax() ||
fZaxis.GetXmin() != h2->fZaxis.GetXmin() ||
fZaxis.GetXmax() != h2->fZaxis.GetXmax()) {
Warning("Add","Attempt to add histograms::Add with different axis limits");
}
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
if (fDimension < 3) nbinsz = -1;
if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
Double_t nEntries = c1*h1->GetEntries() + c2*h2->GetEntries();
Double_t s1[kNstat], s2[kNstat], s3[kNstat];
Int_t i;
for (i=0;i<kNstat;i++) {s1[i] = s2[i] = s3[i] = 0;}
h1->GetStats(s1);
h2->GetStats(s2);
for (i=0;i<kNstat;i++) {
if (i == 1) s3[i] = c1*c1*s1[i] + c2*c2*s2[i];
else s3[i] = TMath::Abs(c1)*s1[i] + TMath::Abs(c2)*s2[i];
}
PutStats(s3);
SetMinimum();
SetMaximum();
ResetBit(kCanRebin);
Int_t bin, binx, biny, binz;
Double_t cu;
for (binz=0;binz<=nbinsz+1;binz++) {
for (biny=0;biny<=nbinsy+1;biny++) {
for (binx=0;binx<=nbinsx+1;binx++) {
bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
cu = c1*h1->GetBinContent(bin)+ c2*h2->GetBinContent(bin);
SetBinContent(bin,cu);
if (fSumw2.fN) {
Double_t error1 = h1->GetBinError(bin);
Double_t error2 = h2->GetBinError(bin);
fSumw2.fArray[bin] = c1*c1*error1*error1 + c2*c2*error2*error2;
}
}
}
}
SetEntries(nEntries);
}
void TH1::AddBinContent(Int_t)
{
AbstractMethod("AddBinContent");
}
void TH1::AddBinContent(Int_t, Double_t)
{
AbstractMethod("AddBinContent");
}
void TH1::AddDirectory(Bool_t add)
{
fgAddDirectory = add;
}
Int_t TH1::BufferEmpty(Int_t action)
{
if (!fBuffer) return 0;
Int_t nbentries = (Int_t)fBuffer[0];
if (!nbentries) return 0;
Double_t *buffer = fBuffer;
if (nbentries < 0) {
if (action == 0) return 0;
nbentries = -nbentries;
fBuffer=0;
Reset();
fBuffer = buffer;
}
if (TestBit(kCanRebin) || (fXaxis.GetXmax() <= fXaxis.GetXmin())) {
Double_t xmin = fBuffer[2];
Double_t xmax = xmin;
for (Int_t i=1;i<nbentries;i++) {
Double_t x = fBuffer[2*i+2];
if (x < xmin) xmin = x;
if (x > xmax) xmax = x;
}
if (fXaxis.GetXmax() <= fXaxis.GetXmin()) {
THLimitsFinder::GetLimitsFinder()->FindGoodLimits(this,xmin,xmax);
} else {
fBuffer = 0;
Int_t keep = fBufferSize; fBufferSize = 0;
if (xmin < fXaxis.GetXmin()) RebinAxis(xmin,"X");
if (xmax >= fXaxis.GetXmax()) RebinAxis(xmax,"X");
fBuffer = buffer;
fBufferSize = keep;
}
}
FillN(nbentries,&fBuffer[2],&fBuffer[1],2);
if (action > 0) { delete [] fBuffer; fBuffer = 0; fBufferSize = 0;}
else {
if (nbentries == (Int_t)fEntries) fBuffer[0] = -nbentries;
else fBuffer[0] = 0;
}
return nbentries;
}
Int_t TH1::BufferFill(Double_t x, Double_t w)
{
if (!fBuffer) return -2;
Int_t nbentries = (Int_t)fBuffer[0];
if (nbentries < 0) {
nbentries = -nbentries;
fBuffer[0] = nbentries;
if (fEntries > 0) {
Double_t *buffer = fBuffer; fBuffer=0;
Reset();
fBuffer = buffer;
}
}
if (2*nbentries+2 >= fBufferSize) {
BufferEmpty(1);
return Fill(x,w);
}
fBuffer[2*nbentries+1] = w;
fBuffer[2*nbentries+2] = x;
fBuffer[0] += 1;
return -2;
}
Double_t TH1::Chi2Test(const TH1* h2, Option_t *option, Double_t *res) const
{
//Begin_Html <!--
/* -->
<html>
<body>
<h1> <IMG WIDTH="50" HEIGHT="44" ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2.gif" ALT="$\chi^2$"> test for comparing weighted and unweighted histograms</h1>
<p>
Function:
Returns p-value. Other return values are specified by the 3rd parameter <br>
Parameters:
<ul>
<li>h2 - the second histogram</li>
<li>option </li>
<ul>
<li>"UU" = experiment experiment comparison (unweighted-unweighted)</li>
<li>"UW" = experiment MC comparison (unweighted-weighted). Note that the first histogram should be unweighted </li>
<li>"WW" = MC MC comparison (weighted-weighted)</li>
<li>"NORM" = to be used when one or both of the histograms is scaled (unweighted-unweighted)</li>
<li>by default underflows and overlows are not included</li>
<ul>
<li>"OF" = overflows included</li>
<li>"UF" = underflows included</li>
</ul>
<li>"P" = print chi2, ndf, p_value, igood</li>
<li>"CHI2" = returns chi2 instead of p-value</li>
<li>"CHI2/NDF" = returns chi2/ndf</li>
</ul>
<li>res: not empty - computes normalized residuals and returns them in this array</li>
</ul>
</p>
<br>
The current implementation is based on the papers "<IMG WIDTH="25" HEIGHT="22" ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2.gif" ALT="$\chi^2$"> test for comparison
of weighted and unweighted histograms" in Proceedings of PHYSTAT05 and
"Comparison weighted and unweighted histograms", arXiv:physics/0605123 by N.Gagunashvili. This function has been implemented
by Daniel Haertl in August 2006.
<h2>Introduction</h2>
A frequently used technique in data analysis is the comparison of histograms.
First suggested by Pearson [1] the <IMG WIDTH="25" HEIGHT="22" ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2.gif" ALT="$\chi^2$"> test of
homogeneity is used widely for comparing usual (unweighted) histograms.
This paper describes the implementation modified <IMG WIDTH="25" HEIGHT="22" ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2.gif" ALT="$\chi^2$"> tests
for comparison of weighted and unweighted histograms and two weighted
histograms [2] as well as usual Pearson's <IMG WIDTH="25" HEIGHT="22" ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2.gif" ALT="$\chi^2$"> test for
comparison two usual (unweighted) histograms.
<h2>Overview</h2>
Comparison of two histograms expect hypotheses that two histograms
represent the identical distributions. To make a decision <I>p</I>-value should be calculated. The hypotheses of identity is rejected if <I>p</I>-value is lower then
some significance level. Traditionally significance levels 0.1, 0.05 and 0.01 are used.
The comparison procedure should include an analysis of the residuals
which is often helpful in identifying the bins of histograms responsible
for a significant overall <i>X<sup>2</sup></i> value. Residuals are the difference between
bin contents and expected bin contents. Most convenient for analysis are the
normalized residuals. If hypotheses of identity are valid then normalized
residuals are approximately independent and identically distributed
random variables having <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_standard.png" ALT="normal"> distribution. Analysis of
residuals expect test of above mentioned properties of residuals.
Notice that indirectly the analysis of residuals increase the power of <IMG WIDTH="25" HEIGHT="22" ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2.gif" ALT="$\chi^2$"> test.
<h2>Methods of comparison</h2>
<h3><IMG WIDTH="50" HEIGHT="44" ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2.gif" ALT="$\chi^2$"> test for comparison two (unweighted) histograms</h3>
Let us consider two histograms with the same
binning and the number of bins equal to <I>r</I>.
Let us denote the number of events in the <I>i</I>th bin in the first histogram as
<i>n<sub>i</sub></i> and as <i>m<sub>i</sub></i> in the second one. The total number of events in the
first histogram is equal to <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_Nsum.png" ALT="$N=\sum_{i=1}^{r}{n_i}$"> ,
and <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_Msum.png" ALT="$M=\sum_{i=1}^{r}{m_i}$"> in the second histogram.
The hypothesis of identity (homogeneity) [3] is that the
two histograms represent random values with identical distributions.
It is equivalent that there exist <I>r</I> constants
<I>p<sub>1</sub>,...,p<sub>r</sub></I>,
such that <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_p_i_sum.png" ALT=" $\sum_{i=1}^{r} p_i=1$"> ,
and the probability of belonging to the <i>i</i>th bin for some measured value
in both experiments is equal to <i>p<sub>i</sub></i>.
The number of events in the <i>i</i>th bin is a random variable
with a distribution approximated by a Poisson probability distribution
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_Npoisson.png" ALT="$e^{-Np_i}(Np_i)^{n_i}/n_i!$ "> for the first histogram and with
distribution <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_Mpoisson.png" ALT="$e^{-Mp_i}(Mp_i)^{m_i}/m_i!$ "> for the second histogram.
If the hypothesis of homogeneity is valid, then the maximum likelihood
estimator of <i>p<sub>i</sub>, i=1,...,r</i>, is
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_ratio.png" ALT="\hat{p}_i= \frac{n_{i}+m_{i}}{N+M}">
</DIV>
and then
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_m1.png" ALT="X^2=\sum_{i=1}^{r}{\frac{(n_{i}-N\hat{p}_i)^2}{N\hat{p}_i}}
+\sum_{i=1}^{r}{\frac{(m_{i}-M\hat{p}_i)^2}{M\hat{p}_i}} =\frac{1}{MN} \sum_{i=1}^{r}{\frac{(Mn_i-Nm_i)^2}{n_i+m_i}}"><IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_m12.png" ALT="X^2=\sum_{i=1}^{r}{\frac{(n_{i}-N\hat{p}_i)^2}{N\hat{p}_i}}
+\sum_{i=1}^{r}{\frac{(m_{i}-M\hat{p}_i)^2}{M\hat{p}_i}} =\frac{1}{MN} \sum_{i=1}^{r}{\frac{(Mn_i-Nm_i)^2}{n_i+m_i}}">
</DIV>
has approximately a <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2r.png" ALT=" $\chi^2_{(r-1)}$"> distribution [3].
The comparison procedure can include an analysis of the residuals which
is often helpful in identifying the bins of histograms responsible for a
significant overall <i>X<sup>2</sup></i> value. Most convenient for analysis are the
adjusted (normalized) residuals [4]
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_res1.png" ALT="$r_i=\frac{n_{i}-N\hat{p}_i}{\sqrt{N\hat{p}_i}\sqrt{(1-N/(N+M))(1-(n_i+m_i)/(N+M))}}$".
</DIV>
If hypotheses of homogeneity are valid then
residuals <i>r<sub>i</sub></i> are approximately independent and identically distributed
random variables having <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_standard.png" ALT="$\mathcal{N}(0,1)$"> distribution.
The application of the <IMG WIDTH="50" HEIGHT="44" ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2.gif" ALT="$\chi^2$"> test has restrictions related to the
value of the expected frequencies <i>Np<sub>i</sub>, Mp<sub>i</sub>, i=1,...,r</i>.
A conservative rule formulated in [5] is that all
the expectations must be 1 or greater for both histograms. In practical cases when expected frequencies are not known the estimated expected frequencies
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_MpNp.png" ALT=" $M\hat{p}_i$, $N\hat{p}_i, i=1,...,r$"> can be used.
<h3>Unweighted and weighted histograms comparison</h3>
A simple modification of the ideas described above can be used for the
comparison of the usual (unweighted) and
weighted histograms. Let us denote the number of events in the <i>i</i>th bin in the unweighted histogram as
<i>n<sub>i</sub></i> and the common weight of events in the <i>i</i>th bin of the
weighted histogram as <i>w<sub>i</sub></i>. The total number of events in the
unweighted histogram is equal to <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_Nsum.png" ALT="$N=\sum_{i=1}^{r}{n_i}$"> and the total
weight of events in the weighted histogram is equal
to <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_Wsum.png" ALT=" $W=\sum_{i=1}^{r}{w_i}$">.
Let us formulate the hypothesis of identity of an unweighted histogram
to a weighted histogram so that there exist <i>r</i> constants <i>p<sub>1</sub>,...,p<sub>r</sub></i>,
such that <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_p_i_sum.png" ALT="$\sum_{i=1}^{r} p_i=1$>, and the probability of belonging to the <i>i</i>th bin for some measured value
is equal to <i>p<sub>i</sub></i> for the unweighted histogram and expectation value of weight <i>w<sub>i</sub></i> equal to <i>Wp<sub>i</sub></i> for the weighted histogram.
The number of events in the <i>i</i>th bin is a random
variable with distribution approximated by the Poisson probability distribution
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_Npoisson.png" ALT="$e^{-Np_i}(Np_i)^{n_i}/n_i!$ "> for the unweighted histogram.
The weight <i>w<sub>i</sub></i> is a random variable with a distribution approximated by
the normal probability distribution <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_standardw.png" ALT=" $ \mathcal{N}(Wp_i,\sigma_i^2)$ ">, where
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_sigma.png" ALT=" $\sigma_i^2$ "> is the variance of the weight <i>w<sub>i</sub></i>.
If we replace the variance <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_sigma.png" ALT=" $\sigma_i^2$ "> with estimate <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_s.png" ALT=" $s_i^2$ "> (sum of squares of weights of events in the <i>i</i>th bin) and
the hypothesis of identity is valid, then the maximum likelihood
estimator of <i>p<sub>i</sub>,i=1,...,r</i>, is
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_ratio2.png" ALT="\hat{p}_i= \frac{Ww_i-Ns_i^2+\sqrt{(Ww_i-Ns_i^2)^2+4W^2s_i^2n_i}}{2W^2} ">.
</DIV>
We may then use the test statistic
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_m2.png" ALT="X^2=\sum_{i=1}^{r}{\frac{(n_{i}-N\hat{p}_i)^2}{N\hat{p}_i}}
+\sum_{i=1}^{r}{\frac{(w_{i}-W\hat{p}_i)^2}{s_i^2}}">
</DIV>
and it has approximately a <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2r.png" ALT=" $\chi^2_{(r-1)}$"> distribution [2].
This test, as well as the original one [3], has a restriction
on the expected frequencies. The expected frequencies
recommended for the weighted histogram is more than 25.
The value of the minimal expected frequency can be decreased down to 10 for
the case when the weights of the events are close to constant.
In the case of a weighted histogram if the number of events is unknown, then we can apply this recommendation for the equivalent number of events as
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_neq.png" ALT="$n_i^{equiv}={w_i^2}/{s_i^2} \, \text{.}$">.
The minimal expected frequency for an unweighted histogram must be 1.
Notice that any usual (unweighted) histogram can be considered as a weighted histogram with events that have constant weights equal to 1.
The variance <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_z.png" ALT="$z_i^2$"> of the difference between the weight <i>w<sub>i</sub></i> and the estimated expectation value of the weight is approximately equal to:
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_zfor1.png" ALT="$z_i^2=Var(w_{i}-W\hat{p}_i)=N\hat{p}_i(1-N\hat{p}_i)\biggl(\frac{Ws_i^2}
{\sqrt{(Ns_i^2-w_iW)^2+4W^2s_i^2n_i}}\biggr)^2\\
+\frac{s_i^2}{4}\biggl(1+\frac{Ns_i^2-w_iW}
{\sqrt{(Ns_i^2-w_iW)^2+4W^2s_i^2n_i}}\biggr)^2$">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_zfor2.png" ALT="$z_i^2=Var(w_{i}-W\hat{p}_i)=N\hat{p}_i(1-N\hat{p}_i)\biggl(\frac{Ws_i^2}
{\sqrt{(Ns_i^2-w_iW)^2+4W^2s_i^2n_i}}\biggr)^2\\
+\frac{s_i^2}{4}\biggl(1+\frac{Ns_i^2-w_iW}
{\sqrt{(Ns_i^2-w_iW)^2+4W^2s_i^2n_i}}\biggr)^2$">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_zfor3.png" ALT="$z_i^2=Var(w_{i}-W\hat{p}_i)=N\hat{p}_i(1-N\hat{p}_i)\biggl(\frac{Ws_i^2}
{\sqrt{(Ns_i^2-w_iW)^2+4W^2s_i^2n_i}}\biggr)^2\\
+\frac{s_i^2}{4}\biggl(1+\frac{Ns_i^2-w_iW}
{\sqrt{(Ns_i^2-w_iW)^2+4W^2s_i^2n_i}}\biggr)^2$">.
</DIV>
The residuals
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_res2.png" ALT="r_i=\frac{w_{i}-W\hat{p}_i}{z_i}">
</DIV>
have approximately a normal distribution with mean equal to 0 and
standard deviation equal to 1.
<h3>Two weighted histograms comparison</h3>
Let us denote the common weight of events of the <i>i</i>th bin in the first histogram as
<i>w<sub>1i</sub></i> and as <i>w<sub>2i</sub></i> in the second one. The total weight of events in the
first histogram is equal to <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_W1sum.png" ALT="$W_1=\sum_{i=1}^{r}{w_{1i}}$">,
and <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_W2sum.png" ALT="$W_2=\sum_{i=1}^{r}{w_{2i}}$"> in the second histogram.
Let us formulate the hypothesis of
identity of weighted histograms so that there exist <i>r</i> constants <i>p<sub>1</sub>,...,p<sub>r</sub></i>,
such that <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_p_i_sum.png" ALT="$\sum_{i=1}^{r} p_i=1$">, and also expectation value of weight <i>w<sub>1i</sub></i> equal to <i>W<sub>1</sub>p<sub>i</sub></i> and expectation value of weight <i>w<sub>2i</sub></i> equal to <i>W<sub>2</sub>p<sub>i</sub></i>.
Weights in both the histograms are random variables with distributions which
can be
approximated by a normal probability distribution <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_standard_w1.png", ALT="$\mathcal{N}(W_1p_i,\sigma_{1i}^2)$">
for the first histogram and by a distribution
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_standard_w2.png", ALT="$\mathcal{N}(W_2p_i,\sigma_{2i}^2)$"> for the second. Here <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_sigma1.png", ALT="$\sigma_{1i}^2$ "> and
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_sigma2.png", ALT="$\sigma_{2i}^2$ "> are the variances of <i>w<sub>1i</sub></i> and <i>w<sub>2i</sub></i> with estimators <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_s1.png", ALT="$s_{1i}^2$ ">
and <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_s2.png", ALT="$s_{2i}^2$ "> respectively. If the hypothesis of identity is valid,
then the maximum likelihood and Least Square Method estimator
of <i>p<sub>i</sub>,i=1,...,r</i>, is
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_ratio3.png", ALT="\hat{p}_i=\frac{w_{1i}W_1/s_{1i}^2+w_{2i}W_2 /s_{2i}^2}{W_1^2/s_{1i}^2+W_2^2/s_{2i}^2} "> .
</DIV>
We may then use the test statistic
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_m3.png", ALT="X^2=\sum_{i=1}^{r}{\frac{(w_{1i}-W_1\hat{p}_i)^2}{s_{1i}^2}}
+\sum_{i=1}^{r}{\frac{(w_{2i}-W_2\hat{p}_i)^2}{s_{2i}^2}}=\sum _{i=1}^{r}{\frac{(W_1w_{2i}-W_2w_{1i})^2}{W_1^2s_{2i}^2+W_2^2s_{1i}^2}}">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_m32.png", ALT="X^2=\sum_{i=1}^{r}{\frac{(w_{1i}-W_1\hat{p}_i)^2}{s_{1i}^2}}
+\sum_{i=1}^{r}{\frac{(w_{2i}-W_2\hat{p}_i)^2}{s_{2i}^2}}=\sum _{i=1}^{r}{\frac{(W_1w_{2i}-W_2w_{1i})^2}{W_1^2s_{2i}^2+W_2^2s_{1i}^2}}">
</DIV>
and it has approximately a <IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2r.png" ALT=" $\chi^2_{(r-1)}$"> distribution [2]. The normalized or studentised residuals [6]
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_res3.png", ALT=" r_i=\frac{w_{1i}-W_1\hat{p}_i}{s_{1i}\sqrt{1-1/(1+W_2^2s_{1i}^2/W_1^2s_{2i}^2)}} ">
</DIV>
have approximately a normal distribution with mean equal to 0 and
standard deviation 1. A recommended minimal expected frequency is equal to 10 for the proposed test.
<h2>Numerical examples</h2>
The method described herein is now illustrated with an example.
We take a distribution
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_example_formula.png", ALT="\phi(x)=\frac{2}{(x-10)^2+1}+\frac{1}{(x-14)^2+1} ">         (1)
</DIV>
defined on the interval [4,16]. Events distributed
according to the formula (1) are simulated to create the unweighted
histogram.
Uniformly distributed events are simulated for the weighted histogram
with weights calculated by formula (1).
Each histogram has the same number of bins: 20.
Fig. 1 shows the result of comparison of the unweighted histogram with
200 events (minimal expected frequency equal to one) and the weighted histogram with 500 events (minimal expected frequency equal to 25)
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_plot1.jpg", ALT="fig1">
</DIV>
<div>
<caption align=left> Fig 1. An example of comparison of the unweighted histogram with 200 events
and the weighted histogram with 500 events: a) unweighted histogram;
b) weighted
histogram; c) normalized residuals plot; d) normal Q-Q plot of residuals.
</caption>
</div>
<BR><P></P>
The value of the test statistic
<i>X<sup>2</sup></i> is equal to 21.09 with <i>p</i>-value equal to 0.33, therefore the
hypothesis
of identity of the two histograms can be accepted for 0.05 significant level. The behavior of the
normalized residuals plot (see Fig. 1c) and the normal Q-Q plot (see Fig. 1d) of residuals are
regular and we cannot identify the outliers or bins with a big influence on
<i>X<sup>2</sup></i>.<br>
<br>
The second example presented the same two histograms but 17 events was added to
content of bin number 15 in unweighted histogram.
Fig. 2 shows the result of comparison of the unweighted histogram with
217 events (minimal expected frequency equal to one) and the weighted histogram with 500 events (minimal expected frequency equal to 25)
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_plot2.jpg", ALT="fig1_bad">
</DIV>
<div>
<caption align=left> Fig 2. An example of comparison of the unweighted histogram with 217 events
and the weighted histogram with 500 events: a) unweighted histogram;
b) weighted
histogram; c) normalized residuals plot; d) normal Q-Q plot of residuals.
</caption>
</div>
<BR><P></P>
The value of the test statistic
<i>X<sup>2</sup></i> is equal to 32.33 with <i>p</i>-value equal to 0.029, therefore the
hypothesis
of identity of the two histograms is rejected for 0.05 significant level. The behavior of the
normalized residuals plot (see Fig. 2c) and the normal Q-Q plot (see Fig. 2d) of residuals are not
regular and we can identify the outlier or bin with a big influence on
<i>X<sup>2</sup></i>.
<h2>References</h2>
[1] Pearson, K., 1904. On the Theory of Contingency and Its Relation to Association
and Normal Correlation. Drapers' Co. Memoirs, Biometric Series No. 1, London.<br>
<br>
[2] Gagunashvili, N., 2006. <IMG WIDTH="25" HEIGHT="22" ALIGN="MIDDLE" BORDER="0" SRC="gif/chi2_chi2.gif" ALT="$\chi^2$"> test for comparison of weighted and
unweighted histograms.
Statistical Problems in Particle Physics, Astrophysics and Cosmology, Proceedings of PHYSTAT05, Oxford, UK, 12-15 September 2005, Imperial College Press, London, 43-44.<br>
    Gagunashvili,N., Comparison of weighted and unweighted histograms, arXiv:physics/0605123, 2006.<br>
<br>
[3] Cramer, H., 1946. Mathematical methods of statistics. Princeton University Press, Princeton.<br>
<br>
[4] Haberman, S.J., 1973. The analysis of residuals in cross-classified tables. Biometrics 29, 205-220.<br>
<br>
[5] Lewontin, R.C. and Felsenstein, J., 1965. The robustness of homogeneity test
in 2 × N tables. Biometrics 21, 19-33. <br>
<br>
[6] Seber, G.A.F., Lee, A.J., 2003, Linear Regression Analysis. John Wiley & Sons Inc., New York.<br>
<body>
</html>
<!--*/
// -->End_Html
Double_t chi2 = 0;
Int_t ndf = 0, igood = 0;
TString opt = option;
opt.ToUpper();
Double_t prob = Chi2TestX(h2,chi2,ndf,igood,option,res);
if(opt.Contains("P")) {
printf("Chi2 = %f, Prob = %g, NDF = %d, igood = %d\n", chi2,prob,ndf,igood);
}
if(opt.Contains("CHI2/NDF")) {
if (ndf == 0) return 0;
return chi2/ndf;
}
if(opt.Contains("CHI2")) {
return chi2;
}
return prob;
}
Double_t TH1::Chi2TestX(const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood, Option_t *option, Double_t *res) const
{
Int_t i, j, k;
Int_t i_start, i_end;
Int_t j_start, j_end;
Int_t k_start, k_end;
Double_t bin1, bin2;
Double_t err1,err2;
Double_t sum1=0, sum2=0;
chi2 = 0;
ndf = 0;
TString opt = option;
opt.ToUpper();
TAxis *xaxis1 = this->GetXaxis();
TAxis *xaxis2 = h2->GetXaxis();
TAxis *yaxis1 = this->GetYaxis();
TAxis *yaxis2 = h2->GetYaxis();
TAxis *zaxis1 = this->GetZaxis();
TAxis *zaxis2 = h2->GetZaxis();
Int_t nbinx1 = xaxis1->GetNbins();
Int_t nbinx2 = xaxis2->GetNbins();
Int_t nbiny1 = yaxis1->GetNbins();
Int_t nbiny2 = yaxis2->GetNbins();
Int_t nbinz1 = zaxis1->GetNbins();
Int_t nbinz2 = zaxis2->GetNbins();
if (this->GetDimension() != h2->GetDimension() ){
Error("ChistatTestX","Histograms have different dimensions.");
return 0;
}
if (nbinx1 != nbinx2) {
Error("ChistatTestX","different number of x channels");
}
if (nbiny1 != nbiny2) {
Error("ChistatTestX","different number of y channels");
}
if (nbinz1 != nbinz2) {
Error("ChistatTestX","different number of z channels");
}
i_start = j_start = k_start = 1;
i_end = nbinx1;
j_end = nbiny1;
k_end = nbinz1;
if (xaxis1->TestBit(TAxis::kAxisRange)) {
i_start = xaxis1->GetFirst();
i_end = xaxis1->GetLast();
}
if (yaxis1->TestBit(TAxis::kAxisRange)) {
j_start = yaxis1->GetFirst();
j_end = yaxis1->GetLast();
}
if (zaxis1->TestBit(TAxis::kAxisRange)) {
k_start = zaxis1->GetFirst();
k_end = zaxis1->GetLast();
}
ndf = (i_end - i_start + 1)*(j_end - j_start + 1)*(k_end - k_start + 1) - 1;
if (opt.Contains("OF")) {
i_end = ++nbinx1;
j_end = ++nbiny1;
k_end = ++nbinz1;
ndf += nbinx1 + nbiny1 + nbinz1;
}
if (opt.Contains("UF")) {
i_start = j_start = k_start = 0;
ndf += nbinx1 + nbiny1 + nbinz1;
}
for(i=i_start; i<=i_end; i++) {
for (j=j_start; j<=j_end; j++) {
for (k=k_start; k<=k_end; k++) {
bin1 = this->GetBinContent(i,j,k);
bin2 = h2->GetBinContent(i,j,k);
if (!opt.Contains("UU") && bin2 <= 0){
Error("ChistatTestX","Hist2: zero events in bin (%d,%d,%d)\n", i,j,k);
return 0;
}
if (opt.Contains("WW") && bin1 <= 0){
Error("ChistatTestX","Hist1: zero events in bin (%d,%d,%d)\n", i,j,k);
return 0;
}
}
}
}
if (opt.Contains("UU") && opt.Contains("NORM")) {
for (i=i_start; i<=i_end; i++) {
for (j=j_start; j<=j_end; j++) {
for (k=k_start; k<=k_end; k++) {
bin1 = this->GetBinContent(i,j,k);
bin2 = h2->GetBinContent(i,j,k);
err1 = this->GetBinError(i,j,k);
err2 = h2->GetBinError(i,j,k);
if (err1==0) continue;
if (err2==0) continue;
bin1 *= bin1/(err1*err1);
bin2 *= bin2/(err2*err2);
bin1 += 0.5;
bin2 += 0.5;
bin1 = Int_t(bin1);
bin2 = Int_t(bin2);
bin1 = Double_t(bin1);
bin2 = Double_t(bin2);
sum1 += bin1;
sum2 += bin2;
}
}
}
} else {
for (i=i_start; i<=i_end; i++) {
for (j=j_start; j<=j_end; j++) {
for (k=k_start; k<=k_end; k++) {
sum1 += this->GetBinContent(i,j,k);
sum2 += h2->GetBinContent(i,j,k);
}
}
}
}
if (sum1 == 0 || sum2 == 0) {
Error("ChistatTestX","one of the histograms is empty");
return 0;
}
Int_t m=0, n=0;
if (opt.Contains("UU")) {
Double_t sum = sum1 + sum2;
Double_t binsum,temp1,temp2,correc;
for (i=i_start; i<=i_end; i++) {
for (j=j_start; j<=j_end; j++) {
for (k=k_start; k<=k_end; k++) {
bin1 = this->GetBinContent(i,j,k);
bin2 = h2->GetBinContent(i,j,k);
if (bin1 == 0 || bin2 == 0) {
--ndf;
} else {
if (opt.Contains("NORM")) {
err1 = this->GetBinError(i,j,k);
err2 = h2->GetBinError(i,j,k);
bin1 *= bin1/(err1*err1);
bin2 *= bin2/(err2*err2);
bin1 += 0.5;
bin2 += 0.5;
bin1 = Int_t(bin1);
bin2 = Int_t(bin2);
bin1 = Double_t(bin1);
bin2 = Double_t(bin2);
}
binsum = bin1 + bin2;
temp1 = binsum*sum1/sum;
temp2 = binsum*sum2/sum;
if (res)
res[i-i_start] = (bin1-temp1)/TMath::Sqrt(temp1);
if (temp1 < 1) m++;
if (temp2 < 1) n++;
correc = (1-sum1/sum)*(1-binsum/sum);
if (res) {
res[i-i_start] /= TMath::Sqrt(correc);
}
temp1 = sum2*bin1-sum1*bin2;
chi2 += temp1*temp1/binsum;
}
}
}
}
chi2 /= sum1*sum2;
if (m) {
igood = 1;
printf("There is bin in Hist1 with less than 1 exp number of events.\n");
}
if (n) {
igood = 2;
printf("There is bin in Hist2 with less than 1 exp number of events.\n");
}
Double_t prob = TMath::Prob(chi2,ndf);
return prob;
}
if (opt.Contains("UW")) {
Double_t var1,var2;
Double_t probb,temp,temp1,temp2;
for (i=i_start; i<=i_end; i++) {
for (j=j_start; j<=j_end; j++) {
for (k=k_start; k<=k_end; k++) {
Int_t x=0, y=0;
bin1 = this->GetBinContent(i,j,k);
bin2 = h2->GetBinContent(i,j,k);
err2 = h2->GetBinError(i,j,k);
err1 *= err1;
err2 *= err2;
var1 = sum2*bin2 - sum1*err2;
var2 = var1*var1 + 4*sum2*sum2*bin1*err2;
while (var1*var1+bin1 == 0 || var1+var2 == 0) {
sum1++;
bin1++;
x++;
y=1;
var1 = sum2*bin2 - sum1*err2;
var2 = var1*var1 + 4*sum2*sum2*bin1*err2;
}
var2 = TMath::Sqrt(var2);
while (var1+var2 == 0) {
sum1++;
bin1++;
x++;
y=1;
var1 = sum2*bin2 - sum1*err2;
var2 = var1*var1 + 4*sum2*sum2*bin1*err2;
while (var1*var1+bin1 == 0 || var1+var2 == 0) {
sum1++;
bin1++;
x++;
y=1;
var1 = sum2*bin2 - sum1*err2;
var2 = var1*var1 + 4*sum2*sum2*bin1*err2;
}
var2 = TMath::Sqrt(var2);
}
probb = (var1+var2)/(2*sum2*sum2);
temp1 = probb * sum1;
temp2 = probb * sum2;
if (temp1 < 1) m++;
if (bin2*bin2/err2 < 10) n++;
temp = bin1 - temp1;
chi2 += temp*temp/temp1;
temp = bin2 - temp2;
chi2 += temp*temp/err2;
temp1 = sum2*err2/var2;
temp2 = 1 + (sum1*err2 - sum2*bin2)/var2;
temp2 = temp1*temp1*sum1*probb*(1-probb) + temp2*temp2*err2/4;
if (res)
res[i-i_start] = temp/TMath::Sqrt(temp2);
}
}
}
if (m) {
igood = 1;
printf("There is bin in Hist1 with less than 1 exp number of events.\n");
}
if (n) {
igood = 2;
printf("There is bin in Hist2 with less than 10 eff number of events.\n");
}
Double_t prob = TMath::Prob(chi2,ndf);
return prob;
}
if (opt.Contains("WW")) {
Double_t temp,temp1,temp2,temp3;
for (i=i_start; i<=i_end; i++) {
for (j=j_start; j<=j_end; j++) {
for (k=k_start; k<=k_end; k++) {
bin1 = this->GetBinContent(i,j,k);
bin2 = h2->GetBinContent(i,j,k);
err1 = this->GetBinError(i,j,k);
err2 = h2->GetBinError(i,j,k);
err1 *= err1;
err2 *= err2;
temp = sum1*sum1*err2 + sum2*sum2*err1;
temp1 = sum2*bin1 - sum1*bin2;
chi2 += temp1*temp1/temp;
temp2 = bin1*sum1*err2 + bin2*sum2*err1;
temp2 *= sum1/temp;
temp2 = bin1-temp2;
temp3 = sum1*sum1 / (sum1*sum1/err1 + sum2*sum2/err2);
if (res)
res[i-i_start] = temp2/TMath::Sqrt(err1-temp3);
bin1 *= bin1/err1;
bin2 *= bin2/err2;
if (bin1 < 10) m++;
if (bin2 < 10) n++;
}
}
}
if (m) {
igood = 1;
printf("There is bin in Hist1 with less than 10 eff number of events.\n");
}
if (n) {
igood = 2;
printf("There is bin in Hist2 with less than 10 eff number of events.\n");
}
Double_t prob = TMath::Prob(chi2,ndf);
return prob;
}
return 0;
}
Double_t TH1::ComputeIntegral()
{
Int_t bin, binx, biny, binz, ibin;
if (fIntegral) delete [] fIntegral;
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
Int_t nxy = nbinsx*nbinsy;
Int_t nbins = nxy*nbinsz;
fIntegral = new Double_t[nbins+2];
ibin = 0;
fIntegral[ibin] = 0;
for (binz=1;binz<=nbinsz;binz++) {
for (biny=1;biny<=nbinsy;biny++) {
for (binx=1;binx<=nbinsx;binx++) {
ibin++;
bin = GetBin(binx, biny, binz);
fIntegral[ibin] = fIntegral[ibin-1] + GetBinContent(bin);
}
}
}
if (fIntegral[nbins] == 0 ) {
Error("ComputeIntegral", "Integral = zero"); return 0;
}
for (bin=1;bin<=nbins;bin++) fIntegral[bin] /= fIntegral[nbins];
fIntegral[nbins+1] = fEntries;
return fIntegral[nbins];
}
Double_t *TH1::GetIntegral()
{
if (!fIntegral) ComputeIntegral();
return fIntegral;
}
void TH1::Copy(TObject &obj) const
{
if (((TH1&)obj).fDirectory) {
((TH1&)obj).fDirectory->GetList()->Remove(&obj);
((TH1&)obj).fDirectory = 0;
}
TNamed::Copy(obj);
((TH1&)obj).fDimension = fDimension;
((TH1&)obj).fNormFactor= fNormFactor;
((TH1&)obj).fEntries = fEntries;
((TH1&)obj).fNcells = fNcells;
((TH1&)obj).fBarOffset = fBarOffset;
((TH1&)obj).fBarWidth = fBarWidth;
((TH1&)obj).fTsumw = fTsumw;
((TH1&)obj).fTsumw2 = fTsumw2;
((TH1&)obj).fTsumwx = fTsumwx;
((TH1&)obj).fTsumwx2 = fTsumwx2;
((TH1&)obj).fMaximum = fMaximum;
((TH1&)obj).fMinimum = fMinimum;
((TH1&)obj).fOption = fOption;
((TH1&)obj).fBuffer = 0;
((TH1&)obj).fBufferSize= fBufferSize;
if (fBuffer) {
Double_t *buf = new Double_t[fBufferSize];
for (Int_t i=0;i<fBufferSize;i++) buf[i] = fBuffer[i];
((TH1&)obj).fBuffer = buf;
}
TAttLine::Copy(((TH1&)obj));
TAttFill::Copy(((TH1&)obj));
TAttMarker::Copy(((TH1&)obj));
fXaxis.Copy(((TH1&)obj).fXaxis);
fYaxis.Copy(((TH1&)obj).fYaxis);
fZaxis.Copy(((TH1&)obj).fZaxis);
((TH1&)obj).fXaxis.SetParent(&obj);
((TH1&)obj).fYaxis.SetParent(&obj);
((TH1&)obj).fZaxis.SetParent(&obj);
fContour.Copy(((TH1&)obj).fContour);
fSumw2.Copy(((TH1&)obj).fSumw2);
if (fgAddDirectory && gDirectory) {
gDirectory->Append(&obj);
((TH1&)obj).fDirectory = gDirectory;
}
}
Int_t TH1::DistancetoPrimitive(Int_t px, Int_t py)
{
if (!fPainter) return 9999;
return fPainter->DistancetoPrimitive(px,py);
}
void TH1::Divide(TF1 *f1, Double_t c1)
{
if (!f1) {
Error("Add","Attempt to divide by a non-existing function");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
Double_t nEntries = fEntries;
Double_t s1[10];
Int_t i;
for (i=0;i<10;i++) {s1[i] = 0;}
PutStats(s1);
SetMinimum();
SetMaximum();
ResetBit(kCanRebin);
Int_t bin, binx, biny, binz;
Double_t cu,w;
Double_t xx[3];
Double_t *params = 0;
f1->InitArgs(xx,params);
for (binz=0;binz<=nbinsz+1;binz++) {
xx[2] = fZaxis.GetBinCenter(binz);
for (biny=0;biny<=nbinsy+1;biny++) {
xx[1] = fYaxis.GetBinCenter(biny);
for (binx=0;binx<=nbinsx+1;binx++) {
xx[0] = fXaxis.GetBinCenter(binx);
if (!f1->IsInside(xx)) continue;
TF1::RejectPoint(kFALSE);
bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
Double_t error1 = GetBinError(bin);
cu = c1*f1->EvalPar(xx);
if (TF1::RejectedPoint()) continue;
if (cu) w = GetBinContent(bin)/cu;
else w = 0;
SetBinContent(bin,w);
if (fSumw2.fN) {
if (cu != 0) fSumw2.fArray[bin] = error1*error1/(cu*cu);
else fSumw2.fArray[bin] = 0;
}
}
}
}
SetEntries(nEntries);
}
void TH1::Divide(const TH1 *h1)
{
if (!h1) {
Error("Divide","Attempt to divide by a non-existing histogram");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()) {
Error("Divide","Attempt to divide histograms with different number of bins");
return;
}
if (fXaxis.GetXmin() != h1->fXaxis.GetXmin() ||
fXaxis.GetXmax() != h1->fXaxis.GetXmax() ||
fYaxis.GetXmin() != h1->fYaxis.GetXmin() ||
fYaxis.GetXmax() != h1->fYaxis.GetXmax() ||
fZaxis.GetXmin() != h1->fZaxis.GetXmin() ||
fZaxis.GetXmax() != h1->fZaxis.GetXmax()) {
Warning("Divide","Attempt to divide histograms with different axis limits");
}
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
if (fDimension < 3) nbinsz = -1;
if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
Double_t nEntries = fEntries;
fEntries = fTsumw = 0;
ResetBit(kCanRebin);
Int_t bin, binx, biny, binz;
Double_t c0,c1,w;
for (binz=0;binz<=nbinsz+1;binz++) {
for (biny=0;biny<=nbinsy+1;biny++) {
for (binx=0;binx<=nbinsx+1;binx++) {
bin = GetBin(binx,biny,binz);
c0 = GetBinContent(bin);
c1 = h1->GetBinContent(bin);
if (c1) w = c0/c1;
else w = 0;
SetBinContent(bin,w);
fEntries++;
if (fSumw2.fN) {
Double_t e0 = GetBinError(bin);
Double_t e1 = h1->GetBinError(bin);
Double_t c12= c1*c1;
if (!c1) { fSumw2.fArray[bin] = 0; continue;}
fSumw2.fArray[bin] = (e0*e0*c1*c1 + e1*e1*c0*c0)/(c12*c12);
}
}
}
}
Double_t s[kNstat];
GetStats(s);
PutStats(s);
SetEntries(nEntries);
}
void TH1::Divide(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
{
TString opt = option;
opt.ToLower();
Bool_t binomial = kFALSE;
if (opt.Contains("b")) binomial = kTRUE;
if (!h1 || !h2) {
Error("Divide","Attempt to divide by a non-existing histogram");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()
|| nbinsx != h2->GetNbinsX() || nbinsy != h2->GetNbinsY() || nbinsz != h2->GetNbinsZ()) {
Error("Divide","Attempt to divide histograms with different number of bins");
return;
}
if (!c2) {
Error("Divide","Coefficient of dividing histogram cannot be zero");
return;
}
if (fXaxis.GetXmin() != h1->fXaxis.GetXmin() ||
fXaxis.GetXmax() != h1->fXaxis.GetXmax() ||
fYaxis.GetXmin() != h1->fYaxis.GetXmin() ||
fYaxis.GetXmax() != h1->fYaxis.GetXmax() ||
fZaxis.GetXmin() != h1->fZaxis.GetXmin() ||
fZaxis.GetXmax() != h1->fZaxis.GetXmax()) {
Warning("Divide","Attempt to divide histograms with different axis limits");
}
if (fXaxis.GetXmin() != h2->fXaxis.GetXmin() ||
fXaxis.GetXmax() != h2->fXaxis.GetXmax() ||
fYaxis.GetXmin() != h2->fYaxis.GetXmin() ||
fYaxis.GetXmax() != h2->fYaxis.GetXmax() ||
fZaxis.GetXmin() != h2->fZaxis.GetXmin() ||
fZaxis.GetXmax() != h2->fZaxis.GetXmax()) {
Warning("Divide","Attempt to divide histograms with different axis limits");
}
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
Double_t nEntries = fEntries;
fEntries = fTsumw = 0;
SetMinimum();
SetMaximum();
ResetBit(kCanRebin);
Int_t bin, binx, biny, binz;
Double_t b1,b2,w,d1,d2;
d1 = c1*c1;
d2 = c2*c2;
for (binz=0;binz<=nbinsz+1;binz++) {
for (biny=0;biny<=nbinsy+1;biny++) {
for (binx=0;binx<=nbinsx+1;binx++) {
bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
b1 = h1->GetBinContent(bin);
b2 = h2->GetBinContent(bin);
if (b2) w = c1*b1/(c2*b2);
else w = 0;
SetBinContent(bin,w);
fEntries++;
if (fSumw2.fN) {
Double_t e1 = h1->GetBinError(bin);
Double_t e2 = h2->GetBinError(bin);
Double_t b22= b2*b2*d2;
if (!b2) { fSumw2.fArray[bin] = 0; continue;}
if (binomial) {
if (b1 != b2) {
fSumw2.fArray[bin] = TMath::Abs( ( (1.-2.*w)*e1*e1 + w*w*e2*e2 )/(b2*b2) );
} else {
Double_t too_low = 0;
Double_t too_high = 1;
Double_t integral;
Double_t a = b1+1;
Double_t x;
for (Int_t loop=0; loop<20; loop++) {
x = 0.5*(too_high + too_low);
Double_t bt = TMath::Exp(TMath::LnGamma(a+1)-TMath::LnGamma(a)+a*log(x)+log(1-x));
if (x < (a+1.0)/(a+3.0)) integral = 1 - bt*TMath::BetaCf(x,a,1)/a;
else integral = bt*TMath::BetaCf(1-x,1,a);
if (integral > 0.683) too_low = x;
else too_high = x;
}
fSumw2.fArray[bin] = (1-x)*(1-x)/4;
}
} else {
fSumw2.fArray[bin] = d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1)/(b22*b22);
}
}
}
}
}
Double_t s[kNstat];
GetStats(s);
PutStats(s);
if (nEntries == 0) nEntries = h1->GetEntries();
if (nEntries == 0) nEntries = 1;
SetEntries(nEntries);
}
void TH1::Draw(Option_t *option)
{
TString opt = option;
opt.ToLower();
if (gPad) {
if (!gPad->IsEditable()) (gROOT->GetMakeDefCanvas())();
if (opt.Contains("same")) {
if (opt.Contains("same") &&
gPad->GetX1() == 0 && gPad->GetX2() == 1 &&
gPad->GetY1() == 0 && gPad->GetY2() == 1 &&
gPad->GetListOfPrimitives()->GetSize()==0) opt.ReplaceAll("same","");
} else {
if (TestBit(kCanDelete)) gPad->GetListOfPrimitives()->Remove(this);
gPad->Clear();
}
} else {
if (opt.Contains("same")) opt.ReplaceAll("same","");
}
AppendPad(opt.Data());
}
TH1 *TH1::DrawCopy(Option_t *) const
{
AbstractMethod("DrawCopy");
return 0;
}
TH1 *TH1::DrawNormalized(Option_t *option, Double_t norm) const
{
Double_t sum = GetSumOfWeights();
if (sum == 0) {
Error("DrawNormalized","Sum of weights is null. Cannot normalized histogram: %s",GetName());
return 0;
}
Bool_t addStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
TH1 *h = (TH1*)Clone();
h->SetBit(kCanDelete);
h->Scale(norm/sum);
h->Draw(option);
TH1::AddDirectory(addStatus);
return h;
}
void TH1::DrawPanel()
{
if (!fPainter) {Draw(); if (gPad) gPad->Update();}
if (fPainter) fPainter->DrawPanel();
}
void TH1::Eval(TF1 *f1, Option_t *option)
{
Double_t x[3];
Int_t range,stat,add,bin,binx,biny,binz,nbinsx, nbinsy, nbinsz;
if (!f1) return;
Double_t fu;
Double_t e=0;
TString opt = option;
opt.ToLower();
if (opt.Contains("a")) add = 1;
else add = 0;
if (opt.Contains("s")) stat = 1;
else stat = 0;
if (opt.Contains("r")) range = 1;
else range = 0;
nbinsx = fXaxis.GetNbins();
nbinsy = fYaxis.GetNbins();
nbinsz = fZaxis.GetNbins();
if (!add) Reset();
for (binz=1;binz<=nbinsz;binz++) {
x[2] = fZaxis.GetBinCenter(binz);
for (biny=1;biny<=nbinsy;biny++) {
x[1] = fYaxis.GetBinCenter(biny);
for (binx=1;binx<=nbinsx;binx++) {
bin = GetBin(binx,biny,binz);
x[0] = fXaxis.GetBinCenter(binx);
if (range && !f1->IsInside(x)) continue;
fu = f1->Eval(x[0],x[1],x[2]);
if (stat) fu = gRandom->PoissonD(fu);
if (fSumw2.fN) e = fSumw2.fArray[bin];
AddBinContent(bin,fu);
if (fSumw2.fN) fSumw2.fArray[bin] = e+ TMath::Abs(fu);
}
}
}
}
void TH1::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
if (fPainter) fPainter->ExecuteEvent(event, px, py);
}
TH1* TH1::FFT(TH1* h_output, Option_t *option)
{
Int_t ndim[3];
ndim[0] = this->GetNbinsX();
ndim[1] = this->GetNbinsY();
ndim[2] = this->GetNbinsZ();
TVirtualFFT *fft;
TString opt = option;
opt.ToUpper();
if (!opt.Contains("2R")){
if (!opt.Contains("2C") && !opt.Contains("2HC") && !opt.Contains("DHT")) {
opt.Append("R2C");
}
fft = TVirtualFFT::FFT(this->GetDimension(), ndim, opt.Data());
}
else {
Int_t ind = opt.Index("R2R", 3);
Int_t *kind = new Int_t[2];
char t;
t = opt[ind+4];
kind[0] = atoi(&t);
if (h_output->GetDimension()>1) {
t = opt[ind+5];
kind[1] = atoi(&t);
}
fft = TVirtualFFT::SineCosine(this->GetDimension(), ndim, kind, option);
delete [] kind;
}
if (!fft) return 0;
Int_t in=0;
for (Int_t binx = 1; binx<=ndim[0]; binx++) {
for (Int_t biny=1; biny<=ndim[1]; biny++) {
for (Int_t binz=1; binz<=ndim[2]; binz++) {
fft->SetPoint(in, this->GetBinContent(binx, biny, binz));
in++;
}
}
}
fft->Transform();
h_output = TransformHisto(fft, h_output, option);
return h_output;
}
Int_t TH1::Fill(Double_t x)
{
if (fBuffer) return BufferFill(x,1);
Int_t bin;
fEntries++;
bin =fXaxis.FindBin(x);
AddBinContent(bin);
if (fSumw2.fN) ++fSumw2.fArray[bin];
if (bin == 0 || bin > fXaxis.GetNbins()) {
if (!fgStatOverflows) return -1;
}
++fTsumw;
++fTsumw2;
fTsumwx += x;
fTsumwx2 += x*x;
return bin;
}
Int_t TH1::Fill(Double_t x, Double_t w)
{
if (fBuffer) return BufferFill(x,w);
Int_t bin;
fEntries++;
bin =fXaxis.FindBin(x);
AddBinContent(bin, w);
if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
if (bin == 0 || bin > fXaxis.GetNbins()) {
if (!fgStatOverflows) return -1;
}
Double_t z= (w > 0 ? w : -w);
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
return bin;
}
Int_t TH1::Fill(const char *namex, Double_t w)
{
Int_t bin;
fEntries++;
bin =fXaxis.FindBin(namex);
AddBinContent(bin, w);
if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
if (bin == 0 || bin > fXaxis.GetNbins()) return -1;
Double_t x = fXaxis.GetBinCenter(bin);
Double_t z= (w > 0 ? w : -w);
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
return bin;
}
void TH1::FillN(Int_t ntimes, const Double_t *x, const Double_t *w, Int_t stride)
{
Int_t bin,i;
fEntries += ntimes;
Double_t ww = 1;
Int_t nbins = fXaxis.GetNbins();
ntimes *= stride;
for (i=0;i<ntimes;i+=stride) {
bin =fXaxis.FindBin(x[i]);
if (w) ww = w[i];
AddBinContent(bin, ww);
if (fSumw2.fN) fSumw2.fArray[bin] += ww*ww;
if (bin == 0 || bin > nbins) {
if (!fgStatOverflows) continue;
}
Double_t z= (ww > 0 ? ww : -ww);
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x[i];
fTsumwx2 += z*x[i]*x[i];
}
}
void TH1::FillRandom(const char *fname, Int_t ntimes)
{
Int_t bin, binx, ibin, loop;
Double_t r1, x, xv[1];
TF1 *f1 = (TF1*)gROOT->GetFunction(fname);
if (!f1) { Error("FillRandom", "Unknown function: %s",fname); return; }
Int_t nbinsx = GetNbinsX();
Double_t *integral = new Double_t[nbinsx+1];
ibin = 0;
integral[ibin] = 0;
for (binx=1;binx<=nbinsx;binx++) {
xv[0] = fXaxis.GetBinCenter(binx);
ibin++;
Double_t fval = f1->Eval(xv[0]);
integral[ibin] = integral[ibin-1] + TMath::Abs(fval)*fXaxis.GetBinWidth(binx);
}
if (integral[nbinsx] == 0 ) {
Error("FillRandom", "Integral = zero"); return;
}
for (bin=1;bin<=nbinsx;bin++) integral[bin] /= integral[nbinsx];
for (loop=0;loop<ntimes;loop++) {
r1 = gRandom->Rndm(loop);
ibin = TMath::BinarySearch(nbinsx,&integral[0],r1);
binx = 1 + ibin;
x = fXaxis.GetBinLowEdge(ibin+1)
+fXaxis.GetBinWidth(ibin+1)*(r1-integral[ibin])/(integral[ibin+1] - integral[ibin]);
Fill(x, 1.);
}
delete [] integral;
}
void TH1::FillRandom(TH1 *h, Int_t ntimes)
{
if (!h) { Error("FillRandom", "Null histogram"); return; }
if (fDimension != h->GetDimension()) {
Error("FillRandom", "Histograms with different dimensions"); return;
}
if (h->ComputeIntegral() == 0) return;
Int_t loop;
Double_t x;
for (loop=0;loop<ntimes;loop++) {
x = h->GetRandom();
Fill(x);
}
}
Int_t TH1::FindBin(Double_t x, Double_t y, Double_t z)
{
if (GetDimension() < 2) {
return fXaxis.FindBin(x);
}
if (GetDimension() < 3) {
Int_t nx = fXaxis.GetNbins()+2;
Int_t binx = fXaxis.FindBin(x);
Int_t biny = fYaxis.FindBin(y);
return binx + nx*biny;
}
if (GetDimension() < 4) {
Int_t nx = fXaxis.GetNbins()+2;
Int_t ny = fYaxis.GetNbins()+2;
Int_t binx = fXaxis.FindBin(x);
Int_t biny = fYaxis.FindBin(y);
Int_t binz = fZaxis.FindBin(z);
return binx + nx*(biny +ny*binz);
}
return -1;
}
TObject *TH1::FindObject(const char *name) const
{
if (fFunctions) return fFunctions->FindObject(name);
return 0;
}
TObject *TH1::FindObject(const TObject *obj) const
{
if (fFunctions) return fFunctions->FindObject(obj);
return 0;
}
Int_t TH1::Fit(const char *fname ,Option_t *option ,Option_t *goption, Double_t xxmin, Double_t xxmax)
{
char *linear;
linear= (char*)strstr(fname, "++");
TF1 *f1=0;
TF2 *f2=0;
TF3 *f3=0;
Int_t ndim=GetDimension();
if (linear){
if (ndim<2){
f1=new TF1(fname, fname, xxmin, xxmax);
return Fit(f1,option,goption,xxmin,xxmax);
}
else if (ndim<3){
f2=new TF2(fname, fname);
return Fit(f2,option,goption,xxmin,xxmax);
}
else{
f3=new TF3(fname, fname);
return Fit(f3,option,goption,xxmin,xxmax);
}
}
else{
f1 = (TF1*)gROOT->GetFunction(fname);
if (!f1) { Printf("Unknown function: %s",fname); return -1; }
return Fit(f1,option,goption,xxmin,xxmax);
}
}
Int_t TH1::Fit(TF1 *f1 ,Option_t *option ,Option_t *goption, Double_t xxmin, Double_t xxmax)
{
Int_t fitResult = 0;
Int_t i, npar,nvpar,nparx;
Double_t par, we, al, bl;
Double_t eplus,eminus,eparab,globcc,amin,aminref,edm,errdef,werr;
Double_t params[100], arglist[100];
Double_t xmin, xmax, ymin, ymax, zmin, zmax, binwidx, binwidy, binwidz;
Int_t hxfirst, hxlast, hyfirst, hylast, hzfirst, hzlast;
TF1 *fnew1;
TF2 *fnew2;
TF3 *fnew3;
if (!f1) {
Error("Fit", "function may not be null pointer");
return 0;
}
if (f1->IsZombie()) {
Error("Fit", "function is zombie");
return 0;
}
npar = f1->GetNpar();
if (npar <= 0) {
Error("Fit", "function %s has illegal number of parameters = %d", f1->GetName(), npar);
return 0;
}
if (f1->GetNdim() != GetDimension()) {
Error("Fit","function %s dimension, %d, does not match histogram dimension, %d",
f1->GetName(), f1->GetNdim(), GetDimension());
return 0;
}
if (fBuffer) BufferEmpty(1);
hxfirst = fXaxis.GetFirst();
hxlast = fXaxis.GetLast();
binwidx = fXaxis.GetBinWidth(hxlast);
xmin = fXaxis.GetBinLowEdge(hxfirst);
xmax = fXaxis.GetBinLowEdge(hxlast) +binwidx;
hyfirst = fYaxis.GetFirst();
hylast = fYaxis.GetLast();
binwidy = fYaxis.GetBinWidth(hylast);
ymin = fYaxis.GetBinLowEdge(hyfirst);
ymax = fYaxis.GetBinLowEdge(hylast) +binwidy;
hzfirst = fZaxis.GetFirst();
hzlast = fZaxis.GetLast();
binwidz = fZaxis.GetBinWidth(hzlast);
zmin = fZaxis.GetBinLowEdge(hzfirst);
zmax = fZaxis.GetBinLowEdge(hzlast) +binwidz;
Foption_t fitOption;
if (!FitOptionsMake(option,fitOption)) return 0;
if (xxmin != xxmax) {
f1->SetRange(xxmin,ymin,zmin,xxmax,ymax,zmax);
fitOption.Range = 1;
}
Int_t special = f1->GetNumber();
Bool_t linear = f1->IsLinear();
if (special==299+npar)
linear = kTRUE;
if (fitOption.Bound || fitOption.Like || fitOption.Errors || fitOption.Gradient || fitOption.More || fitOption.User|| fitOption.Integral || fitOption.Minuit)
linear = kFALSE;
char l[] ="TLinearFitter";
Int_t strdiff = 0;
Bool_t isSet = kFALSE;
if (TVirtualFitter::GetFitter()){
isSet=kTRUE;
strdiff = strcmp(TVirtualFitter::GetFitter()->IsA()->GetName(), l);
}
if (linear) {
TClass *cl = gROOT->GetClass("TLinearFitter");
if (isSet && strdiff!=0) {
delete TVirtualFitter::GetFitter();
isSet=kFALSE;
}
if (!isSet && cl) {
TVirtualFitter::SetFitter((TVirtualFitter *)cl->New());
}
} else {
if (isSet && strdiff==0){
delete TVirtualFitter::GetFitter();
isSet=kFALSE;
}
if (!isSet)
TVirtualFitter::SetFitter(0);
}
TVirtualFitter *hFitter = TVirtualFitter::Fitter(this, f1->GetNpar());
if (!hFitter) return 0;
hFitter->Clear();
gF1 = f1;
hFitter->SetUserFunc(f1);
if (xxmin != xxmax) f1->SetRange(xxmin,ymin,zmin,xxmax,ymax,zmax);
hFitter->SetFitOption(fitOption);
Double_t fxmin, fymin, fzmin, fxmax, fymax, fzmax;
if (fitOption.Range) {
f1->GetRange(fxmin, fymin, fzmin, fxmax, fymax, fzmax);
if (fxmin > xmin) xmin = fxmin;
if (fymin > ymin) ymin = fymin;
if (fzmin > zmin) zmin = fzmin;
if (fxmax < xmax) xmax = fxmax;
if (fymax < ymax) ymax = fymax;
if (fzmax < zmax) zmax = fzmax;
hxfirst = fXaxis.FindFixBin(xmin); if (hxfirst < 1) hxfirst = 1;
hxlast = fXaxis.FindFixBin(xmax); if (hxlast > fXaxis.GetLast()) hxlast = fXaxis.GetLast();
hyfirst = fYaxis.FindFixBin(ymin); if (hyfirst < 1) hyfirst = 1;
hylast = fYaxis.FindFixBin(ymax); if (hylast > fYaxis.GetLast()) hylast = fYaxis.GetLast();
hzfirst = fZaxis.FindFixBin(zmin); if (hzfirst < 1) hzfirst = 1;
hzlast = fZaxis.FindFixBin(zmax); if (hzlast > fZaxis.GetLast()) hzlast = fZaxis.GetLast();
} else {
f1->SetRange(xmin,ymin,zmin,xmax,ymax,zmax);
}
hFitter->SetXfirst(hxfirst); hFitter->SetXlast(hxlast);
hFitter->SetYfirst(hyfirst); hFitter->SetYlast(hylast);
hFitter->SetZfirst(hzfirst); hFitter->SetZlast(hzlast);
Int_t maxpoints = (hzlast-hzfirst+1)*(hylast-hyfirst+1)*(hxlast-hxfirst+1);
Int_t psize = 2 +fDimension;
if (fitOption.Integral) psize = 2+2*fDimension;
Double_t *cache = hFitter->SetCache(maxpoints,psize);
Int_t np = 0;
for (Int_t binz=hzfirst;binz<=hzlast;binz++) {
for (Int_t biny=hyfirst;biny<=hylast;biny++) {
for (Int_t binx=hxfirst;binx<=hxlast;binx++) {
if (fitOption.Integral) {
if (fDimension > 2) {
cache[6] = fZaxis.GetBinCenter(binz);
cache[7] = fZaxis.GetBinWidth(binz);
}
if (fDimension > 1) {
cache[4] = fYaxis.GetBinCenter(biny);
cache[5] = fYaxis.GetBinWidth(biny);
}
cache[2] = fXaxis.GetBinCenter(binx);
cache[3] = fXaxis.GetBinWidth(binx);
} else {
if (fDimension > 2) {
cache[4] = fZaxis.GetBinCenter(binz);
}
if (fDimension > 1) {
cache[3] = fYaxis.GetBinCenter(biny);
}
cache[2] = fXaxis.GetBinCenter(binx);
}
if (!f1->IsInside(&cache[2])) continue;
Int_t bin = GetBin(binx,biny,binz);
cache[0] = GetBinContent(bin);
cache[1] = GetBinError(bin);
if (fitOption.W1) {
if (fitOption.W1 == 1 && cache[0] == 0) continue;
cache[1] = 1;
}
if (cache[1] == 0) {
if (fitOption.Like) cache[1] = 1;
else continue;
}
np++;
cache += psize;
}
}
}
hFitter->SetCache(np,psize);
if (linear){
fitResult = hFitter->ExecuteCommand("FitHist", 0, 0);
} else {
if (fitOption.Bound) special = 0;
if (special == 100) H1InitGaus();
else if (special == 400) H1InitGaus();
else if (special == 200) H1InitExpo();
else if (special == 299+npar) H1InitPolynom();
if (!fitOption.Verbose) {
arglist[0] = -1;
hFitter->ExecuteCommand("SET PRINT", arglist,1);
arglist[0] = 0;
hFitter->ExecuteCommand("SET NOW", arglist,0);
}
arglist[0] = TVirtualFitter::GetErrorDef();
if (fitOption.Like) {
hFitter->SetFitMethod("H1FitLikelihood");
} else {
if (!fitOption.User) hFitter->SetFitMethod("H1FitChisquare");
}
hFitter->ExecuteCommand("SET Err",arglist,1);
Int_t nfixed = 0;
for (i=0;i<npar;i++) {
par = f1->GetParameter(i);
f1->GetParLimits(i,al,bl);
if (al*bl != 0 && al >= bl) {
al = bl = 0;
arglist[nfixed] = i+1;
nfixed++;
}
we = 0.1*TMath::Abs(bl-al);
if (we == 0) we = 0.3*TMath::Abs(par);
if (we == 0) we = binwidx;
hFitter->SetParameter(i,f1->GetParName(i),par,we,al,bl);
}
if(nfixed > 0)hFitter->ExecuteCommand("FIX",arglist,nfixed);
if (fitOption.Gradient) {
if (fitOption.Gradient == 1) arglist[0] = 1;
else arglist[0] = 0;
hFitter->ExecuteCommand("SET GRAD",arglist,1);
}
if (fitOption.Verbose) {
arglist[0] = 0; hFitter->ExecuteCommand("SET PRINT", arglist,1);
}
Double_t ey, sumw2=0;
for (i=hxfirst;i<=hxlast;i++) {
ey = GetBinError(i);
sumw2 += ey*ey;
}
arglist[0] = TVirtualFitter::GetMaxIterations();
arglist[1] = sumw2*TVirtualFitter::GetPrecision();
fitResult = hFitter->ExecuteCommand("MIGRAD",arglist,2);
if (fitResult != 0) {
if (!fitOption.Quiet) {
Warning("Fit","Abnormal termination of minimization.");
}
}
if (fitOption.More) {
hFitter->ExecuteCommand("IMPROVE",arglist,0);
}
if (fitOption.Errors) {
hFitter->ExecuteCommand("HESSE",arglist,0);
hFitter->ExecuteCommand("MINOS",arglist,0);
}
char parName[50];
for (i=0;i<npar;i++) {
hFitter->GetParameter(i,parName, par,we,al,bl);
if (!fitOption.Errors) werr = we;
else {
hFitter->GetErrors(i,eplus,eminus,eparab,globcc);
if (eplus > 0 && eminus < 0) werr = 0.5*(eplus-eminus);
else werr = we;
}
params[i] = par;
f1->SetParameter(i,par);
f1->SetParError(i,werr);
}
hFitter->GetStats(aminref,edm,errdef,nvpar,nparx);
amin = aminref;
if (fitOption.Like) amin = hFitter->Chisquare(npar, params);
f1->SetChisquare(amin);
f1->SetNDF(f1->GetNumberFitPoints()-npar+nfixed);
}
if (!fitOption.Quiet) {
if (fitOption.Errors) hFitter->PrintResults(4,aminref);
else hFitter->PrintResults(3,aminref);
}
if (!fitOption.Nostore) {
if (!fitOption.Plus) {
TIter next(fFunctions, kIterBackward);
TObject *obj;
while ((obj = next())) {
if (obj->InheritsFrom(TF1::Class())) {
fFunctions->Remove(obj);
delete obj;
}
}
}
if (GetDimension() < 2) {
fnew1 = new TF1();
f1->Copy(*fnew1);
fFunctions->Add(fnew1);
fnew1->SetParent(this);
fnew1->Save(xmin,xmax,0,0,0,0);
if (fitOption.Nograph) fnew1->SetBit(TF1::kNotDraw);
fnew1->SetBit(TFormula::kNotGlobal);
} else if (GetDimension() < 3) {
fnew2 = new TF2();
f1->Copy(*fnew2);
fFunctions->Add(fnew2);
fnew2->SetParent(this);
fnew2->Save(xmin,xmax,ymin,ymax,0,0);
if (fitOption.Nograph) fnew2->SetBit(TF1::kNotDraw);
fnew2->SetBit(TFormula::kNotGlobal);
} else {
fnew3 = new TF3();
f1->Copy(*fnew3);
fFunctions->Add(fnew3);
fnew3->SetParent(this);
fnew3->SetBit(TFormula::kNotGlobal);
}
if (TestBit(kCanDelete)) return fitResult;
if (!fitOption.Nograph && GetDimension() < 3) Draw(goption);
}
return fitResult;
}
void TH1::FitPanel()
{
if (fPainter) fPainter->FitPanel();
}
TH1 *TH1::GetAsymmetry(TH1* h2, Double_t c2, Double_t dc2)
{
Bool_t addStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
TH1 *asym = (TH1*)Clone();
asym->Sumw2();
TH1 *top = (TH1*)asym->Clone();
TH1 *bottom = (TH1*)asym->Clone();
TH1::AddDirectory(addStatus);
TH1 *h1 = this;
top->Add(h1,h2,1,-c2);
bottom->Add(h1,h2,1,c2);
asym->Divide(top,bottom);
Int_t xmax = asym->GetNbinsX();
Int_t ymax = asym->GetNbinsY();
Int_t zmax = asym->GetNbinsZ();
Double_t bot, error, a, b, da, db;
for(Int_t i=1; i<= xmax; i++){
for(Int_t j=1; j<= ymax; j++){
for(Int_t k=1; k<= zmax; k++){
a = h1->GetBinContent(i,j,k);
b = h2->GetBinContent(i,j,k);
bot = bottom->GetBinContent(i,j,k);
if(bot < 1e-6){}
else{
da = h1->GetBinError(i,j,k);
db = h2->GetBinError(i,j,k);
error = 2*TMath::Sqrt(a*a*c2*c2*db*db + c2*c2*b*b*da*da+a*a*b*b*dc2*dc2)/(bot*bot);
asym->SetBinError(i,j,k,error);
}
}
}
}
delete top;
delete bottom;
return asym;
}
Int_t TH1::GetDefaultBufferSize()
{
return fgBufferSize;
}
Bool_t TH1::GetDefaultSumw2()
{
return fgDefaultSumw2;
}
Double_t TH1::GetEntries() const
{
if (fBuffer) ((TH1*)this)->BufferEmpty();
return fEntries;
}
Double_t TH1::GetEffectiveEntries() const
{
Stat_t s[kNstat];
this->GetStats(s);
return (s[1] ? s[0]*s[0]/s[1] : 0.);
}
char *TH1::GetObjectInfo(Int_t px, Int_t py) const
{
return ((TH1*)this)->GetPainter()->GetObjectInfo(px,py);
}
TVirtualHistPainter *TH1::GetPainter(Option_t *option)
{
if (!fPainter) {
TString opt = option;
opt.ToLower();
if (opt.Contains("gl") || gStyle->GetCanvasPreferGL()) {
TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TGLHistPainter");
if (handler && handler->LoadPlugin() != -1)
fPainter = reinterpret_cast<TVirtualHistPainter *>(handler->ExecPlugin(1, this));
}
}
if (!fPainter) fPainter = TVirtualHistPainter::HistPainter(this);
return fPainter;
}
Int_t TH1::GetQuantiles(Int_t nprobSum, Double_t *q, const Double_t *probSum)
{
if (GetDimension() > 1) {
Error("GetQuantiles","Only available for 1-d histograms");
return 0;
}
const Int_t nbins = GetXaxis()->GetNbins();
if (!fIntegral) ComputeIntegral();
if (fIntegral && fIntegral[nbins+1] != fEntries) ComputeIntegral();
Int_t i, ibin;
Double_t *prob = (Double_t*)probSum;
Int_t nq = nprobSum;
if (probSum == 0) {
nq = nbins+1;
prob = new Double_t[nq];
prob[0] = 0;
for (i=1;i<nq;i++) {
prob[i] = fIntegral[i]/fIntegral[nbins];
}
}
for (i = 0; i < nq; i++) {
ibin = TMath::BinarySearch(nbins,fIntegral,prob[i]);
while (ibin < nbins-1 && fIntegral[ibin+1] == prob[i]) {
if (fIntegral[ibin+2] == prob[i]) ibin++;
else break;
}
q[i] = GetBinLowEdge(ibin+1);
const Double_t dint = fIntegral[ibin+1]-fIntegral[ibin];
if (dint > 0) q[i] += GetBinWidth(ibin+1)*(prob[i]-fIntegral[ibin])/dint;
}
if (!probSum) delete [] prob;
return nq;
}
Int_t TH1::FitOptionsMake(Option_t *choptin, Foption_t &fitOption)
{
Int_t nch = strlen(choptin);
if (!nch) return 1;
char chopt[32];
strcpy(chopt,choptin);
for (Int_t i=0;i<nch;i++) chopt[i] = toupper(choptin[i]);
if (strstr(chopt,"Q")) fitOption.Quiet = 1;
if (strstr(chopt,"V")) {fitOption.Verbose = 1; fitOption.Quiet = 0;}
if (strstr(chopt,"L")) fitOption.Like = 1;
if (strstr(chopt,"LL")) fitOption.Like = 2;
if (strstr(chopt,"W")) fitOption.W1 = 1;
if (strstr(chopt,"WW")) fitOption.W1 = 2;
if (strstr(chopt,"E")) fitOption.Errors = 1;
if (strstr(chopt,"M")) fitOption.More = 1;
if (strstr(chopt,"R")) fitOption.Range = 1;
if (strstr(chopt,"G")) fitOption.Gradient= 1;
if (strstr(chopt,"N")) fitOption.Nostore = 1;
if (strstr(chopt,"0")) fitOption.Nograph = 1;
if (strstr(chopt,"+")) fitOption.Plus = 1;
if (strstr(chopt,"I")) fitOption.Integral= 1;
if (strstr(chopt,"B")) fitOption.Bound = 1;
if (strstr(chopt,"U")) {fitOption.User = 1; fitOption.Like = 0;}
if (strstr(chopt,"F")) fitOption.Minuit = 1;
if (strstr(chopt,"C")) fitOption.Nochisq = 1;
return 1;
}
void H1InitGaus()
{
Double_t allcha, sumx, sumx2, x, val, rms, mean;
Int_t bin;
const Double_t sqrtpi = 2.506628;
TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
TH1 *curHist = (TH1*)hFitter->GetObjectFit();
Int_t hxfirst = hFitter->GetXfirst();
Int_t hxlast = hFitter->GetXlast();
Double_t valmax = curHist->GetBinContent(hxfirst);
Double_t binwidx = curHist->GetBinWidth(hxfirst);
allcha = sumx = sumx2 = 0;
for (bin=hxfirst;bin<=hxlast;bin++) {
x = curHist->GetBinCenter(bin);
val = TMath::Abs(curHist->GetBinContent(bin));
if (val > valmax) valmax = val;
sumx += val*x;
sumx2 += val*x*x;
allcha += val;
}
if (allcha == 0) return;
mean = sumx/allcha;
rms = sumx2/allcha - mean*mean;
if (rms > 0) rms = TMath::Sqrt(rms);
else rms = 0;
if (rms == 0) rms = binwidx*(hxlast-hxfirst+1)/4;
Double_t constant = 0.5*(valmax+binwidx*allcha/(sqrtpi*rms));
Double_t xmin = curHist->GetXaxis()->GetXmin();
Double_t xmax = curHist->GetXaxis()->GetXmax();
if ((mean < xmin || mean > xmax) && rms > (xmax-xmin)) {
mean = 0.5*(xmax+xmin);
rms = 0.5*(xmax-xmin);
}
TF1 *f1 = (TF1*)hFitter->GetUserFunc();
f1->SetParameter(0,constant);
f1->SetParameter(1,mean);
f1->SetParameter(2,rms);
f1->SetParLimits(2,0,10*rms);
}
void H1InitExpo()
{
Double_t constant, slope;
Int_t ifail;
TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
Int_t hxfirst = hFitter->GetXfirst();
Int_t hxlast = hFitter->GetXlast();
Int_t nchanx = hxlast - hxfirst + 1;
H1LeastSquareLinearFit(-nchanx, constant, slope, ifail);
TF1 *f1 = (TF1*)hFitter->GetUserFunc();
f1->SetParameter(0,constant);
f1->SetParameter(1,slope);
}
void H1InitPolynom()
{
Double_t fitpar[25];
TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
TF1 *f1 = (TF1*)hFitter->GetUserFunc();
Int_t hxfirst = hFitter->GetXfirst();
Int_t hxlast = hFitter->GetXlast();
Int_t nchanx = hxlast - hxfirst + 1;
Int_t npar = f1->GetNpar();
if (nchanx <=1 || npar == 1) {
TH1 *curHist = (TH1*)hFitter->GetObjectFit();
fitpar[0] = curHist->GetSumOfWeights()/Double_t(nchanx);
} else {
H1LeastSquareFit( nchanx, npar, fitpar);
}
for (Int_t i=0;i<npar;i++) f1->SetParameter(i, fitpar[i]);
}
void H1LeastSquareFit(Int_t n, Int_t m, Double_t *a)
{
const Double_t zero = 0.;
const Double_t one = 1.;
const Int_t idim = 20;
Double_t b[400] ;
Int_t i, k, l, ifail;
Double_t power;
Double_t da[20], xk, yk;
if (m <= 2) {
H1LeastSquareLinearFit(n, a[0], a[1], ifail);
return;
}
if (m > idim || m > n) return;
b[0] = Double_t(n);
da[0] = zero;
for (l = 2; l <= m; ++l) {
b[l-1] = zero;
b[m + l*20 - 21] = zero;
da[l-1] = zero;
}
TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
TH1 *curHist = (TH1*)hFitter->GetObjectFit();
Int_t hxfirst = hFitter->GetXfirst();
Int_t hxlast = hFitter->GetXlast();
for (k = hxfirst; k <= hxlast; ++k) {
xk = curHist->GetBinCenter(k);
yk = curHist->GetBinContent(k);
power = one;
da[0] += yk;
for (l = 2; l <= m; ++l) {
power *= xk;
b[l-1] += power;
da[l-1] += power*yk;
}
for (l = 2; l <= m; ++l) {
power *= xk;
b[m + l*20 - 21] += power;
}
}
for (i = 3; i <= m; ++i) {
for (k = i; k <= m; ++k) {
b[k - 1 + (i-1)*20 - 21] = b[k + (i-2)*20 - 21];
}
}
H1LeastSquareSeqnd(m, b, idim, ifail, 1, da);
for (i=0; i<m; ++i) a[i] = da[i];
}
void H1LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail)
{
Double_t xbar, ybar, x2bar;
Int_t i, n;
Double_t xybar;
Double_t fn, xk, yk;
Double_t det;
n = TMath::Abs(ndata);
ifail = -2;
xbar = ybar = x2bar = xybar = 0;
TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
TH1 *curHist = (TH1*)hFitter->GetObjectFit();
Int_t hxfirst = hFitter->GetXfirst();
Int_t hxlast = hFitter->GetXlast();
for (i = hxfirst; i <= hxlast; ++i) {
xk = curHist->GetBinCenter(i);
yk = curHist->GetBinContent(i);
if (ndata < 0) {
if (yk <= 0) yk = 1e-9;
yk = TMath::Log(yk);
}
xbar += xk;
ybar += yk;
x2bar += xk*xk;
xybar += xk*yk;
}
fn = Double_t(n);
det = fn*x2bar - xbar*xbar;
ifail = -1;
if (det <= 0) {
a0 = ybar/fn;
a1 = 0;
return;
}
ifail = 0;
a0 = (x2bar*ybar - xbar*xybar) / det;
a1 = (fn*xybar - xbar*ybar) / det;
}
void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b)
{
Int_t a_dim1, a_offset, b_dim1, b_offset;
Int_t nmjp1, i, j, l;
Int_t im1, jp1, nm1, nmi;
Double_t s1, s21, s22;
const Double_t one = 1.;
b_dim1 = idim;
b_offset = b_dim1 + 1;
b -= b_offset;
a_dim1 = idim;
a_offset = a_dim1 + 1;
a -= a_offset;
if (idim < n) return;
ifail = 0;
for (j = 1; j <= n; ++j) {
if (a[j + j*a_dim1] <= 0) { ifail = -1; return; }
a[j + j*a_dim1] = one / a[j + j*a_dim1];
if (j == n) continue;
jp1 = j + 1;
for (l = jp1; l <= n; ++l) {
a[j + l*a_dim1] = a[j + j*a_dim1] * a[l + j*a_dim1];
s1 = -a[l + (j+1)*a_dim1];
for (i = 1; i <= j; ++i) { s1 = a[l + i*a_dim1] * a[i + (j+1)*a_dim1] + s1; }
a[l + (j+1)*a_dim1] = -s1;
}
}
if (k <= 0) return;
for (l = 1; l <= k; ++l) {
b[l*b_dim1 + 1] = a[a_dim1 + 1]*b[l*b_dim1 + 1];
}
if (n == 1) return;
for (l = 1; l <= k; ++l) {
for (i = 2; i <= n; ++i) {
im1 = i - 1;
s21 = -b[i + l*b_dim1];
for (j = 1; j <= im1; ++j) {
s21 = a[i + j*a_dim1]*b[j + l*b_dim1] + s21;
}
b[i + l*b_dim1] = -a[i + i*a_dim1]*s21;
}
nm1 = n - 1;
for (i = 1; i <= nm1; ++i) {
nmi = n - i;
s22 = -b[nmi + l*b_dim1];
for (j = 1; j <= i; ++j) {
nmjp1 = n - j + 1;
s22 = a[nmi + nmjp1*a_dim1]*b[nmjp1 + l*b_dim1] + s22;
}
b[nmi + l*b_dim1] = -s22;
}
}
}
Int_t TH1::GetBin(Int_t binx, Int_t biny, Int_t binz) const
{
Int_t nx, ny, nz;
if (GetDimension() < 2) {
nx = fXaxis.GetNbins()+2;
if (binx < 0) binx = 0;
if (binx >= nx) binx = nx-1;
return binx;
}
if (GetDimension() < 3) {
nx = fXaxis.GetNbins()+2;
if (binx < 0) binx = 0;
if (binx >= nx) binx = nx-1;
ny = fYaxis.GetNbins()+2;
if (biny < 0) biny = 0;
if (biny >= ny) biny = ny-1;
return binx + nx*biny;
}
if (GetDimension() < 4) {
nx = fXaxis.GetNbins()+2;
if (binx < 0) binx = 0;
if (binx >= nx) binx = nx-1;
ny = fYaxis.GetNbins()+2;
if (biny < 0) biny = 0;
if (biny >= ny) biny = ny-1;
nz = fZaxis.GetNbins()+2;
if (binz < 0) binz = 0;
if (binz >= nz) binz = nz-1;
return binx + nx*(biny +ny*binz);
}
return -1;
}
Double_t TH1::GetRandom() const
{
if (fDimension > 1) {
Error("GetRandom","Function only valid for 1-d histograms");
return 0;
}
Int_t nbinsx = GetNbinsX();
Double_t integral;
if (fIntegral) {
if (fIntegral[nbinsx+1] != fEntries) integral = ((TH1*)this)->ComputeIntegral();
} else {
integral = ((TH1*)this)->ComputeIntegral();
if (integral == 0 || fIntegral == 0) return 0;
}
Double_t r1 = gRandom->Rndm();
Int_t ibin = TMath::BinarySearch(nbinsx,fIntegral,r1);
Double_t x = GetBinLowEdge(ibin+1);
if (r1 > fIntegral[ibin]) x +=
GetBinWidth(ibin+1)*(r1-fIntegral[ibin])/(fIntegral[ibin+1] - fIntegral[ibin]);
return x;
}
Double_t TH1::GetBinContent(Int_t) const
{
AbstractMethod("GetBinContent");
return 0;
}
Double_t TH1::GetBinContent(Int_t binx, Int_t biny) const
{
Int_t bin = GetBin(binx,biny);
return GetBinContent(bin);
}
Double_t TH1::GetBinContent(Int_t binx, Int_t biny, Int_t binz) const
{
Int_t bin = GetBin(binx,biny,binz);
return GetBinContent(bin);
}
Double_t TH1::GetBinWithContent(Double_t c, Int_t &binx, Int_t firstx, Int_t lastx,Double_t maxdiff) const
{
if (fDimension > 1) {
binx = 0;
Error("GetBinWithContent","function is only valid for 1-D histograms");
return 0;
}
if (firstx <= 0) firstx = 1;
if (lastx < firstx) lastx = fXaxis.GetNbins();
Int_t binminx = 0;
Double_t diff, curmax = 1.e240;
for (Int_t i=firstx;i<=lastx;i++) {
diff = TMath::Abs(GetBinContent(i)-c);
if (diff <= 0) {binx = i; return diff;}
if (diff < curmax && diff <= maxdiff) {curmax = diff, binminx=i;}
}
binx = binminx;
return curmax;
}
TAxis *TH1::GetXaxis() const
{
return &((TH1*)this)->fXaxis;
}
TAxis *TH1::GetYaxis() const
{
return &((TH1*)this)->fYaxis;
}
TAxis *TH1::GetZaxis() const
{
return &((TH1*)this)->fZaxis;
}
void TH1::LabelsDeflate(Option_t *ax)
{
Int_t iaxis = AxisChoice(ax);
TAxis *axis = 0;
if (iaxis == 1) axis = GetXaxis();
if (iaxis == 2) axis = GetYaxis();
if (iaxis == 3) axis = GetZaxis();
if (!axis) return;
if (!axis->GetLabels()) return;
TIter next(axis->GetLabels());
TObject *obj;
Int_t nbins = 0;
while ((obj = next())) {
if (obj->GetUniqueID()) nbins++;
}
if (nbins < 1) nbins = 1;
TH1 *hold = (TH1*)Clone();
hold->SetDirectory(0);
Bool_t timedisp = axis->GetTimeDisplay();
Double_t xmin = axis->GetXmin();
Double_t xmax = axis->GetBinUpEdge(nbins);
if (xmax <= xmin) xmax = xmin +nbins;
axis->SetRange(0,0);
axis->Set(nbins,xmin,xmax);
Int_t nbinsx = hold->GetXaxis()->GetNbins();
Int_t nbinsy = hold->GetYaxis()->GetNbins();
Int_t nbinsz = hold->GetZaxis()->GetNbins();
Int_t ncells = nbinsx+2;
if (GetDimension() > 1) ncells *= nbinsy+2;
if (GetDimension() > 2) ncells *= nbinsz+2;
SetBinsLength(ncells);
Int_t errors = fSumw2.fN;
if (errors) fSumw2.Set(ncells);
axis->SetTimeDisplay(timedisp);
Double_t err,cu;
Int_t bin,ibin,binx,biny,binz;
Double_t oldEntries = fEntries;
for (binz=1;binz<=nbinsz;binz++) {
for (biny=1;biny<=nbinsy;biny++) {
for (binx=1;binx<=nbinsx;binx++) {
bin = hold->GetBin(binx,biny,binz);
ibin= GetBin(binx,biny,binz);
cu = hold->GetBinContent(bin);
SetBinContent(ibin,cu);
if (errors) {
err = hold->GetBinError(bin);
SetBinError(ibin,err);
}
}
}
}
fEntries = oldEntries;
delete hold;
}
void TH1::LabelsInflate(Option_t *ax)
{
Int_t iaxis = AxisChoice(ax);
TAxis *axis = 0;
if (iaxis == 1) axis = GetXaxis();
if (iaxis == 2) axis = GetYaxis();
if (iaxis == 3) axis = GetZaxis();
if (!axis) return;
TH1 *hold = (TH1*)Clone();
hold->SetDirectory(0);
Bool_t timedisp = axis->GetTimeDisplay();
Int_t nbxold = fXaxis.GetNbins();
Int_t nbyold = fYaxis.GetNbins();
Int_t nbzold = fZaxis.GetNbins();
Int_t nbins = axis->GetNbins();
Double_t xmin = axis->GetXmin();
Double_t xmax = axis->GetXmax();
xmax = xmin + 2*(xmax-xmin);
axis->SetRange(0,0);
axis->Set(2*nbins,xmin,xmax);
Int_t nbinsx = fXaxis.GetNbins();
Int_t nbinsy = fYaxis.GetNbins();
Int_t nbinsz = fZaxis.GetNbins();
Int_t ncells = nbinsx+2;
if (GetDimension() > 1) ncells *= nbinsy+2;
if (GetDimension() > 2) ncells *= nbinsz+2;
SetBinsLength(ncells);
Int_t errors = fSumw2.fN;
if (errors) fSumw2.Set(ncells);
axis->SetTimeDisplay(timedisp);
Double_t err,cu;
Double_t oldEntries = fEntries;
Int_t bin,ibin,binx,biny,binz;
for (binz=1;binz<=nbinsz;binz++) {
for (biny=1;biny<=nbinsy;biny++) {
for (binx=1;binx<=nbinsx;binx++) {
bin = hold->GetBin(binx,biny,binz);
ibin= GetBin(binx,biny,binz);
if (binx > nbxold || biny > nbyold || binz > nbzold) bin = -1;
if (bin > 0) cu = hold->GetBinContent(bin);
else cu = 0;
SetBinContent(ibin,cu);
if (errors) {
if (bin > 0) err = hold->GetBinError(bin);
else err = 0;
SetBinError(ibin,err);
}
}
}
}
fEntries = oldEntries;
delete hold;
}
void TH1::LabelsOption(Option_t *option, Option_t *ax)
{
Int_t iaxis = AxisChoice(ax);
TAxis *axis = 0;
if (iaxis == 1) axis = GetXaxis();
if (iaxis == 2) axis = GetYaxis();
if (iaxis == 3) axis = GetZaxis();
if (!axis) return;
THashList *labels = axis->GetLabels();
if (!labels) {
Warning("LabelsOption","Cannot sort. No labels");
return;
}
TString opt = option;
opt.ToLower();
if (opt.Contains("h")) {
axis->SetBit(TAxis::kLabelsHori);
axis->ResetBit(TAxis::kLabelsVert);
axis->ResetBit(TAxis::kLabelsDown);
axis->ResetBit(TAxis::kLabelsUp);
}
if (opt.Contains("v")) {
axis->SetBit(TAxis::kLabelsVert);
axis->ResetBit(TAxis::kLabelsHori);
axis->ResetBit(TAxis::kLabelsDown);
axis->ResetBit(TAxis::kLabelsUp);
}
if (opt.Contains("u")) {
axis->SetBit(TAxis::kLabelsUp);
axis->ResetBit(TAxis::kLabelsVert);
axis->ResetBit(TAxis::kLabelsDown);
axis->ResetBit(TAxis::kLabelsHori);
}
if (opt.Contains("d")) {
axis->SetBit(TAxis::kLabelsDown);
axis->ResetBit(TAxis::kLabelsVert);
axis->ResetBit(TAxis::kLabelsHori);
axis->ResetBit(TAxis::kLabelsUp);
}
Int_t sort = -1;
if (opt.Contains("a")) sort = 0;
if (opt.Contains(">")) sort = 1;
if (opt.Contains("<")) sort = 2;
if (sort < 0) return;
if (sort > 0 && GetDimension() > 2) {
Error("LabelsOption","Sorting by value not implemented for 3-D histograms");
return;
}
Double_t entries = fEntries;
Int_t n = TMath::Min(axis->GetNbins(), labels->GetSize());
Int_t *a = new Int_t[n+2];
Int_t i,j,k;
Double_t *cont = 0;
Double_t *errors = 0;
THashList *labold = new THashList(labels->GetSize(),1);
TIter nextold(labels);
TObject *obj;
while ((obj=nextold())) {
labold->Add(obj);
}
labels->Clear();
if (sort > 0) {
if (GetDimension() == 1) {
cont = new Double_t[n];
if (fSumw2.fN) errors = new Double_t[n];
for (i=1;i<=n;i++) {
cont[i-1] = GetBinContent(i);
if (errors) errors[i-1] = GetBinError(i);
}
if (sort ==1) TMath::Sort(n,cont,a,kTRUE);
else TMath::Sort(n,cont,a,kFALSE);
for (i=1;i<=n;i++) {
SetBinContent(i,cont[a[i-1]]);
if (errors) SetBinError(i,errors[a[i-1]]);
}
for (i=1;i<=n;i++) {
obj = labold->At(a[i-1]);
labels->Add(obj);
obj->SetUniqueID(i);
}
} else if (GetDimension()== 2) {
Double_t *pcont = new Double_t[n+2];
for (i=0;i<=n;i++) pcont[i] = 0;
Int_t nx = fXaxis.GetNbins();
Int_t ny = fYaxis.GetNbins();
cont = new Double_t[(nx+2)*(ny+2)];
if (fSumw2.fN) errors = new Double_t[n];
for (i=1;i<=nx;i++) {
for (j=1;j<=ny;j++) {
cont[i+nx*j] = GetBinContent(i,j);
if (errors) errors[i+nx*j] = GetBinError(i,j);
if (axis == GetXaxis()) k = i;
else k = j;
pcont[k-1] += cont[i+nx*j];
}
}
if (sort ==1) TMath::Sort(n,pcont,a,kTRUE);
else TMath::Sort(n,pcont,a,kFALSE);
for (i=0;i<n;i++) {
obj = labold->At(a[i]);
labels->Add(obj);
obj->SetUniqueID(i+1);
}
delete [] pcont;
for (i=1;i<=nx;i++) {
for (j=1;j<=ny;j++) {
if (axis == GetXaxis()) {
SetBinContent(i,j,cont[a[i-1]+1+nx*j]);
if (errors) SetBinError(i,j,errors[a[i-1]+1+nx*j]);
} else {
SetBinContent(i,j,cont[i+nx*(a[j-1]+1)]);
if (errors) SetBinError(i,j,errors[i+nx*(a[j-1]+1)]);
}
}
}
} else {
}
} else {
const UInt_t kUsed = 1<<18;
TObject *objk=0;
a[0] = 0;
a[n+1] = n+1;
for (i=1;i<=n;i++) {
const char *label = "zzzzzzzzzzzz";
for (j=1;j<=n;j++) {
obj = labold->At(j-1);
if (!obj) continue;
if (obj->TestBit(kUsed)) continue;
if (strcmp(label,obj->GetName()) < 0) continue;
objk = obj;
a[i] = j;
label = obj->GetName();
}
if (objk) {
objk->SetUniqueID(i);
labels->Add(objk);
objk->SetBit(kUsed);
}
}
for (i=1;i<=n;i++) {
obj = labels->At(i-1);
if (!obj) continue;
obj->ResetBit(kUsed);
}
if (GetDimension() == 1) {
cont = new Double_t[n+2];
if (fSumw2.fN) errors = new Double_t[n+2];
for (i=1;i<=n;i++) {
cont[i] = GetBinContent(a[i]);
if (errors) errors[i] = GetBinError(a[i]);
}
for (i=1;i<=n;i++) {
SetBinContent(i,cont[i]);
if (errors) SetBinError(i,errors[i]);
}
} else if (GetDimension()== 2) {
Int_t nx = fXaxis.GetNbins()+2;
Int_t ny = fYaxis.GetNbins()+2;
cont = new Double_t[nx*ny];
if (fSumw2.fN) errors = new Double_t[nx*ny];
for (i=0;i<nx;i++) {
for (j=0;j<ny;j++) {
cont[i+nx*j] = GetBinContent(i,j);
if (errors) errors[i+nx*j] = GetBinError(i,j);
}
}
for (i=0;i<nx;i++) {
for (j=0;j<ny;j++) {
if (axis == GetXaxis()) {
SetBinContent(i,j,cont[a[i]+nx*j]);
if (errors) SetBinError(i,j,errors[a[i]+nx*j]);
} else {
SetBinContent(i,j,cont[i+nx*a[j]]);
if (errors) SetBinError(i,j,errors[i+nx*a[j]]);
}
}
}
} else {
Int_t nx = fXaxis.GetNbins()+2;
Int_t ny = fYaxis.GetNbins()+2;
Int_t nz = fZaxis.GetNbins()+2;
cont = new Double_t[nx*ny*nz];
if (fSumw2.fN) errors = new Double_t[nx*ny*nz];
for (i=0;i<nx;i++) {
for (j=0;j<ny;j++) {
for (k=0;k<nz;k++) {
cont[i+nx*(j+ny*k)] = GetBinContent(i,j,k);
if (errors) errors[i+nx*(j+ny*k)] = GetBinError(i,j,k);
}
}
}
for (i=0;i<nx;i++) {
for (j=0;j<ny;j++) {
for (k=0;k<nz;k++) {
if (axis == GetXaxis()) {
SetBinContent(i,j,k,cont[a[i]+nx*(j+ny*k)]);
if (errors) SetBinError(i,j,k,errors[a[i]+nx*(j+ny*k)]);
} else if (axis == GetYaxis()) {
SetBinContent(i,j,k,cont[i+nx*(a[j]+ny*k)]);
if (errors) SetBinError(i,j,k,errors[i+nx*(a[j]+ny*k)]);
} else {
SetBinContent(i,j,k,cont[i+nx*(j+ny*a[k])]);
if (errors) SetBinError(i,j,k,errors[i+nx*(j+ny*a[k])]);
}
}
}
}
}
}
fEntries = entries;
delete labold;
if (a) delete [] a;
if (cont) delete [] cont;
if (errors) delete [] errors;
}
static Bool_t AlmostEqual(Double_t a, Double_t b, Double_t epsilon = 0.00000001)
{
return TMath::Abs(a - b) < epsilon;
}
static Bool_t AlmostInteger(Double_t a, Double_t epsilon = 0.00000001)
{
return AlmostEqual(a - TMath::Floor(a), 0, epsilon) ||
AlmostEqual(a - TMath::Floor(a), 1, epsilon);
}
Bool_t TH1::SameLimitsAndNBins(const TAxis& axis1, const TAxis& axis2)
{
if ((axis1.GetNbins() == axis2.GetNbins())
&& (axis1.GetXmin() == axis2.GetXmin())
&& (axis1.GetXmax() == axis2.GetXmax()))
return kTRUE;
else
return kFALSE;
}
Bool_t TH1::RecomputeAxisLimits(TAxis& destAxis, const TAxis& anAxis)
{
if (SameLimitsAndNBins(destAxis, anAxis))
return kTRUE;
if (destAxis.GetXbins()->fN || anAxis.GetXbins()->fN)
return kFALSE;
Double_t width1 = destAxis.GetBinWidth(0);
Double_t width2 = anAxis.GetBinWidth(0);
if (width1 == 0 || width2 == 0)
return kFALSE;
Double_t xmin = TMath::Min(destAxis.GetXmin(), anAxis.GetXmin());
Double_t xmax = TMath::Max(destAxis.GetXmax(), anAxis.GetXmax());
Double_t width = TMath::Max(width1, width2);
if (!AlmostInteger(width/width1) || !AlmostInteger(width/width2))
return kFALSE;
Double_t delta;
delta = (destAxis.GetXmin() - xmin)/width1;
if (!AlmostInteger(delta))
xmin -= (TMath::Ceil(delta) - delta)*width1;
delta = (anAxis.GetXmin() - xmin)/width2;
if (!AlmostInteger(delta))
xmin -= (TMath::Ceil(delta) - delta)*width2;
delta = (destAxis.GetXmin() - xmin)/width1;
if (!AlmostInteger(delta))
return kFALSE;
delta = (xmax - destAxis.GetXmax())/width1;
if (!AlmostInteger(delta))
xmax += (TMath::Ceil(delta) - delta)*width1;
delta = (xmax - anAxis.GetXmax())/width2;
if (!AlmostInteger(delta))
xmax += (TMath::Ceil(delta) - delta)*width2;
delta = (xmax - destAxis.GetXmax())/width1;
if (!AlmostInteger(delta))
return kFALSE;
#ifdef DEBUG
if (!AlmostInteger((xmax - xmin) / width)) {
printf("TH1::RecomputeAxisLimits - Impossible\n");
return kFALSE;
}
#endif
destAxis.Set(TMath::Nint((xmax - xmin)/width), xmin, xmax);
return kTRUE;
}
Long64_t TH1::Merge(TCollection *li)
{
if (!li) return 0;
if (li->IsEmpty()) return (Int_t) GetEntries();
Bool_t mustCleanup = TestBit(kMustCleanup);
if (mustCleanup) ResetBit(kMustCleanup);
TList inlist;
TH1* hclone = (TH1*)Clone("FirstClone");
if (mustCleanup) SetBit(kMustCleanup);
R__ASSERT(hclone);
BufferEmpty(1);
Reset();
SetEntries(0);
inlist.Add(hclone);
inlist.AddAll(li);
THashList allLabels;
THashList* labels=GetXaxis()->GetLabels();
Bool_t haveOneLabel=kFALSE;
if (labels) {
TIter iL(labels);
TObjString* lb;
while ((lb=(TObjString*)iL())) {
haveOneLabel |= (lb && lb->String().Length());
if (!allLabels.FindObject(lb))
allLabels.Add(lb);
}
}
TAxis newXAxis;
Bool_t initialLimitsFound = kFALSE;
Bool_t allHaveLabels = haveOneLabel;
Bool_t same = kTRUE;
Bool_t allHaveLimits = kTRUE;
TIter next(&inlist);
while (TObject *o = next()) {
TH1* h = dynamic_cast<TH1*> (o);
if (!h) {
Error("Add","Attempt to add object of class: %s to a %s",
o->ClassName(),this->ClassName());
return -1;
}
Bool_t hasLimits = h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax();
allHaveLimits = allHaveLimits && hasLimits;
if (hasLimits) {
h->BufferEmpty();
if (!initialLimitsFound) {
initialLimitsFound = kTRUE;
newXAxis.Set(h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),
h->GetXaxis()->GetXmax());
}
else {
if (!RecomputeAxisLimits(newXAxis, *(h->GetXaxis()))) {
Error("Merge", "Cannot merge histograms - limits are inconsistent:\n "
"first: (%d, %f, %f), second: (%d, %f, %f)",
newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax(),
h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),
h->GetXaxis()->GetXmax());
}
}
}
if (allHaveLabels) {
THashList* labels=h->GetXaxis()->GetLabels();
Bool_t haveOneLabel=kFALSE;
if (labels) {
TIter iL(labels);
TObjString* lb;
while ((lb=(TObjString*)iL())) {
haveOneLabel |= (lb && lb->String().Length());
if (!allLabels.FindObject(lb)) {
allLabels.Add(lb);
same = kFALSE;
}
}
}
allHaveLabels&=(labels && haveOneLabel);
if (!allHaveLabels)
Warning("Merge","Not all histograms have labels. I will ignore labels,"
" falling back to bin numbering mode.");
}
}
next.Reset();
same = same && SameLimitsAndNBins(newXAxis, *GetXaxis());
if (!same && initialLimitsFound)
SetBins(newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax());
if (!allHaveLimits) {
while (TH1* h = (TH1*)next()) {
if (h->GetXaxis()->GetXmin() >= h->GetXaxis()->GetXmax() && h->fBuffer) {
Int_t nbentries = (Int_t)h->fBuffer[0];
for (Int_t i = 0; i < nbentries; i++)
Fill(h->fBuffer[2*i + 2], h->fBuffer[2*i + 1]);
}
}
if (!initialLimitsFound)
return (Int_t) GetEntries();
next.Reset();
}
Double_t stats[kNstat], totstats[kNstat];
for (Int_t i=0;i<kNstat;i++) {totstats[i] = stats[i] = 0;}
GetStats(totstats);
Double_t nentries = GetEntries();
Int_t binx, ix, nx;
Double_t cu;
Bool_t canRebin=TestBit(kCanRebin);
ResetBit(kCanRebin);
while (TH1* h=(TH1*)next()) {
if (h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax()) {
h->GetStats(stats);
for (Int_t i=0;i<kNstat;i++)
totstats[i] += stats[i];
nentries += h->GetEntries();
nx = h->GetXaxis()->GetNbins();
for (binx = 0; binx <= nx + 1; binx++) {
cu = h->GetBinContent(binx);
if (!allHaveLabels || !binx || binx==nx+1) {
if ((!same) && (binx == 0 || binx == nx + 1)) {
if (cu != 0) {
Error("Merge", "Cannot merge histograms - the histograms have"
" different limits and undeflows/overflows are present."
" The initial histogram is now broken!");
return -1;
}
}
ix = fXaxis.FindBin(h->GetXaxis()->GetBinCenter(binx));
} else {
const char* label=h->GetXaxis()->GetBinLabel(binx);
if (!label) label="";
ix = fXaxis.FindBin(label);
}
if (ix >= 0) AddBinContent(ix,cu);
if (fSumw2.fN) {
Double_t error1 = h->GetBinError(binx);
fSumw2.fArray[ix] += error1*error1;
}
}
}
}
if (canRebin) SetBit(kCanRebin);
PutStats(totstats);
SetEntries(nentries);
inlist.Remove(hclone);
delete hclone;
return (Long64_t)nentries;
}
void TH1::Multiply(TF1 *f1, Double_t c1)
{
if (!f1) {
Error("Add","Attempt to multiply by a non-existing function");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
Double_t nEntries = fEntries;
Double_t s1[10];
Int_t i;
for (i=0;i<10;i++) {s1[i] = 0;}
PutStats(s1);
SetMinimum();
SetMaximum();
ResetBit(kCanRebin);
Int_t bin, binx, biny, binz;
Double_t cu,w;
Double_t xx[3];
Double_t *params = 0;
f1->InitArgs(xx,params);
for (binz=0;binz<=nbinsz+1;binz++) {
xx[2] = fZaxis.GetBinCenter(binz);
for (biny=0;biny<=nbinsy+1;biny++) {
xx[1] = fYaxis.GetBinCenter(biny);
for (binx=0;binx<=nbinsx+1;binx++) {
xx[0] = fXaxis.GetBinCenter(binx);
if (!f1->IsInside(xx)) continue;
TF1::RejectPoint(kFALSE);
bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
Double_t error1 = GetBinError(bin);
cu = c1*f1->EvalPar(xx);
if (TF1::RejectedPoint()) continue;
w = GetBinContent(bin)*cu;
SetBinContent(bin,w);
if (fSumw2.fN) {
fSumw2.fArray[bin] = cu*cu*error1*error1;
}
}
}
}
SetEntries(nEntries);
}
void TH1::Multiply(const TH1 *h1)
{
if (!h1) {
Error("Multiply","Attempt to multiply by a non-existing histogram");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()) {
Error("Multiply","Attempt to multiply histograms with different number of bins");
return;
}
if (fXaxis.GetXmin() != h1->fXaxis.GetXmin() ||
fXaxis.GetXmax() != h1->fXaxis.GetXmax() ||
fYaxis.GetXmin() != h1->fYaxis.GetXmin() ||
fYaxis.GetXmax() != h1->fYaxis.GetXmax() ||
fZaxis.GetXmin() != h1->fZaxis.GetXmin() ||
fZaxis.GetXmax() != h1->fZaxis.GetXmax()) {
Warning("Multiply","Attempt to multiply histograms with different axis limits");
}
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
if (fDimension < 3) nbinsz = -1;
if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
Double_t nEntries = fEntries;
fEntries = fTsumw = 0;
SetMinimum();
SetMaximum();
ResetBit(kCanRebin);
Int_t bin, binx, biny, binz;
Double_t c0,c1,w;
for (binz=0;binz<=nbinsz+1;binz++) {
for (biny=0;biny<=nbinsy+1;biny++) {
for (binx=0;binx<=nbinsx+1;binx++) {
bin = GetBin(binx,biny,binz);
c0 = GetBinContent(bin);
c1 = h1->GetBinContent(bin);
w = c0*c1;
SetBinContent(bin,w);
fEntries++;
if (fSumw2.fN) {
Double_t e0 = GetBinError(bin);
Double_t e1 = h1->GetBinError(bin);
fSumw2.fArray[bin] = (e0*e0*c1*c1 + e1*e1*c0*c0);
}
}
}
}
Double_t s[kNstat];
GetStats(s);
PutStats(s);
SetEntries(nEntries);
}
void TH1::Multiply(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
{
TString opt = option;
opt.ToLower();
if (!h1 || !h2) {
Error("Multiply","Attempt to multiply by a non-existing histogram");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()
|| nbinsx != h2->GetNbinsX() || nbinsy != h2->GetNbinsY() || nbinsz != h2->GetNbinsZ()) {
Error("Multiply","Attempt to multiply histograms with different number of bins");
return;
}
if (fXaxis.GetXmin() != h1->fXaxis.GetXmin() ||
fXaxis.GetXmax() != h1->fXaxis.GetXmax() ||
fYaxis.GetXmin() != h1->fYaxis.GetXmin() ||
fYaxis.GetXmax() != h1->fYaxis.GetXmax() ||
fZaxis.GetXmin() != h1->fZaxis.GetXmin() ||
fZaxis.GetXmax() != h1->fZaxis.GetXmax()) {
Warning("Multiply","Attempt to multiply histograms with different axis limits");
}
if (fXaxis.GetXmin() != h2->fXaxis.GetXmin() ||
fXaxis.GetXmax() != h2->fXaxis.GetXmax() ||
fYaxis.GetXmin() != h2->fYaxis.GetXmin() ||
fYaxis.GetXmax() != h2->fYaxis.GetXmax() ||
fZaxis.GetXmin() != h2->fZaxis.GetXmin() ||
fZaxis.GetXmax() != h2->fZaxis.GetXmax()) {
Warning("Multiply","Attempt to multiply histograms with different axis limits");
}
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
Double_t nEntries = fEntries;
fEntries = fTsumw = fTsumw2 = fTsumwx = fTsumwx2 = 0;
SetMinimum();
SetMaximum();
ResetBit(kCanRebin);
Int_t bin, binx, biny, binz;
Double_t b1,b2,w,d1,d2;
d1 = c1*c1;
d2 = c2*c2;
for (binz=0;binz<=nbinsz+1;binz++) {
for (biny=0;biny<=nbinsy+1;biny++) {
for (binx=0;binx<=nbinsx+1;binx++) {
bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
b1 = h1->GetBinContent(bin);
b2 = h2->GetBinContent(bin);
w = (c1*b1)*(c2*b2);
SetBinContent(bin,w);
fEntries++;
if (fSumw2.fN) {
Double_t e1 = h1->GetBinError(bin);
Double_t e2 = h2->GetBinError(bin);
fSumw2.fArray[bin] = d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1);
}
}
}
}
Double_t s[kNstat];
GetStats(s);
PutStats(s);
SetEntries(nEntries);
}
void TH1::Paint(Option_t *option)
{
GetPainter(option);
if (fPainter) {
if (strlen(option) > 0) fPainter->Paint(option);
else fPainter->Paint(fOption.Data());
}
}
TH1 *TH1::Rebin(Int_t ngroup, const char*newname, const Double_t *xbins)
{
Int_t nbins = fXaxis.GetNbins();
Double_t xmin = fXaxis.GetXmin();
Double_t xmax = fXaxis.GetXmax();
if ((ngroup <= 0) || (ngroup > nbins)) {
Error("Rebin", "Illegal value of ngroup=%d",ngroup);
return 0;
}
if (fDimension > 1 || InheritsFrom("TProfile")) {
Error("Rebin", "Operation valid on 1-D histograms only");
return 0;
}
Int_t newbins = nbins/ngroup;
if (xbins) newbins = ngroup;
Double_t entries = fEntries;
Double_t *oldBins = new Double_t[nbins+2];
Int_t bin, i;
for (bin=0;bin<nbins+2;bin++) oldBins[bin] = GetBinContent(bin);
Double_t *oldErrors = 0;
if (fSumw2.fN != 0) {
oldErrors = new Double_t[nbins+2];
for (bin=0;bin<nbins+2;bin++) oldErrors[bin] = GetBinError(bin);
}
TH1 *hnew = this;
if ((newname && strlen(newname) > 0) || xbins) {
hnew = (TH1*)Clone(newname);
}
if(!xbins && (newbins*ngroup != nbins)) {
xmax = fXaxis.GetBinUpEdge(newbins*ngroup);
hnew->fTsumw = 0;
}
Int_t nDivisions = fXaxis.GetNdivisions();
Color_t axisColor = fXaxis.GetAxisColor();
Color_t labelColor = fXaxis.GetLabelColor();
Style_t labelFont = fXaxis.GetLabelFont();
Float_t labelOffset = fXaxis.GetLabelOffset();
Float_t labelSize = fXaxis.GetLabelSize();
Float_t tickLength = fXaxis.GetTickLength();
Float_t titleOffset = fXaxis.GetTitleOffset();
Float_t titleSize = fXaxis.GetTitleSize();
Color_t titleColor = fXaxis.GetTitleColor();
Style_t titleFont = fXaxis.GetTitleFont();
if(!xbins && (fXaxis.GetXbins()->GetSize() > 0)){
Double_t *bins = new Double_t[newbins+1];
for(Int_t i = 0; i <= newbins; ++i) bins[i] = fXaxis.GetBinLowEdge(1+i*ngroup);
hnew->SetBins(newbins,bins);
delete [] bins;
} else if (xbins) {
ngroup = 1;
hnew->SetBins(newbins,xbins);
} else {
hnew->SetBins(newbins,xmin,xmax);
}
fXaxis.SetNdivisions(nDivisions);
fXaxis.SetAxisColor(axisColor);
fXaxis.SetLabelColor(labelColor);
fXaxis.SetLabelFont(labelFont);
fXaxis.SetLabelOffset(labelOffset);
fXaxis.SetLabelSize(labelSize);
fXaxis.SetTickLength(tickLength);
fXaxis.SetTitleOffset(titleOffset);
fXaxis.SetTitleSize(titleSize);
fXaxis.SetTitleColor(titleColor);
fXaxis.SetTitleFont(titleFont);
Int_t oldbin = 1;
Double_t binContent, binError;
for (bin = 1;bin<=newbins;bin++) {
binContent = 0;
binError = 0;
for (i=0;i<ngroup;i++) {
if (oldbin+i > nbins) break;
binContent += oldBins[oldbin+i];
if (oldErrors) binError += oldErrors[oldbin+i]*oldErrors[oldbin+i];
}
hnew->SetBinContent(bin,binContent);
if (oldErrors) hnew->SetBinError(bin,TMath::Sqrt(binError));
oldbin += ngroup;
}
hnew->SetBinContent(0,oldBins[0]);
hnew->SetBinContent(newbins+1,oldBins[nbins+1]);
hnew->SetEntries(entries);
delete [] oldBins;
if (oldErrors) delete [] oldErrors;
return hnew;
}
Bool_t TH1::FindNewAxisLimits(const TAxis* axis, const Double_t point, Double_t& newMin, Double_t &newMax)
{
Double_t xmin = axis->GetXmin();
Double_t xmax = axis->GetXmax();
if (xmin >= xmax) return kFALSE;
Double_t range = xmax-xmin;
Double_t binsize = range / axis->GetNbins();
Int_t ntimes = 0;
while (point < xmin) {
if (ntimes++ > 64)
return kFALSE;
xmin = xmin - range;
range *= 2;
binsize *= 2;
if ( xmin / binsize - TMath::Floor(xmin / binsize) >= 0.5) {
xmin += 0.5 * binsize;
xmax += 0.5 * binsize;
}
}
while (point >= xmax) {
if (ntimes++ > 64)
return kFALSE;
xmax = xmax + range;
range *= 2;
binsize *= 2;
if ( xmin / binsize - TMath::Floor(xmin / binsize) >= 0.5) {
xmin -= 0.5 * binsize;
xmax -= 0.5 * binsize;
}
}
newMin = xmin;
newMax = xmax;
return kTRUE;
}
void TH1::RebinAxis(Double_t x, const char *ax)
{
if (!TestBit(kCanRebin)) return;
if (TMath::IsNaN(x)) {
ResetBit(kCanRebin);
return;
}
char achoice = toupper(ax[0]);
TAxis *axis = &fXaxis;
if (achoice == 'Y') axis = &fYaxis;
if (achoice == 'Z') axis = &fZaxis;
if (axis->GetXmin() >= axis->GetXmax()) return;
if (axis->GetNbins() <= 0) return;
Double_t xmin, xmax;
if (!FindNewAxisLimits(axis, x, xmin, xmax))
return;
TH1 *hold = (TH1*)Clone();
hold->SetDirectory(0);
axis->SetLimits(xmin,xmax);
Int_t nbinsx = fXaxis.GetNbins();
Int_t nbinsy = fYaxis.GetNbins();
Int_t nbinsz = fZaxis.GetNbins();
Double_t err,cu;
Double_t bx,by,bz;
Int_t errors = GetSumw2N();
Int_t ix,iy,iz,ibin,binx,biny,binz,bin;
Reset("ICE");
for (binz=1;binz<=nbinsz;binz++) {
bz = hold->GetZaxis()->GetBinCenter(binz);
iz = fZaxis.FindFixBin(bz);
for (biny=1;biny<=nbinsy;biny++) {
by = hold->GetYaxis()->GetBinCenter(biny);
iy = fYaxis.FindFixBin(by);
for (binx=1;binx<=nbinsx;binx++) {
bx = hold->GetXaxis()->GetBinCenter(binx);
ix = fXaxis.FindFixBin(bx);
bin = hold->GetBin(binx,biny,binz);
ibin= GetBin(ix,iy,iz);
cu = hold->GetBinContent(bin);
AddBinContent(ibin,cu);
if (errors) {
err = hold->GetBinError(bin);
fSumw2.fArray[ibin] += err*err;
}
}
}
}
delete hold;
}
void TH1::RecursiveRemove(TObject *obj)
{
if (fFunctions) {
if (!fFunctions->TestBit(kInvalidObject)) fFunctions->RecursiveRemove(obj);
}
}
void TH1::Scale(Double_t c1)
{
Double_t ent = fEntries;
Add(this,this,c1,0);
fEntries = ent;
Int_t ncontours = GetContour();
if (ncontours == 0) return;
Double_t *levels = fContour.GetArray();
for (Int_t i=0;i<ncontours;i++) {
levels[i] *= c1;
}
}
void TH1::SetDefaultBufferSize(Int_t buffersize)
{
if (buffersize < 0) buffersize = 0;
fgBufferSize = buffersize;
}
void TH1::SetDefaultSumw2(Bool_t sumw2)
{
fgDefaultSumw2 = sumw2;
}
void TH1::SetTitle(const char *title)
{
fTitle = title;
TString str1 = fTitle, str2;
Int_t isc = str1.Index(";");
Int_t lns = str1.Length();
if (isc >=0 ) {
fTitle = str1(0,isc);
str1 = str1(isc+1, lns);
isc = str1.Index(";");
if (isc >=0 ) {
str2 = str1(0,isc);
fXaxis.SetTitle(str2.Data());
lns = str1.Length();
str1 = str1(isc+1, lns);
isc = str1.Index(";");
if (isc >=0 ) {
str2 = str1(0,isc);
fYaxis.SetTitle(str2.Data());
lns = str1.Length();
str1 = str1(isc+1, lns);
fZaxis.SetTitle(str1.Data());
} else {
fYaxis.SetTitle(str1.Data());
}
} else {
fXaxis.SetTitle(str1.Data());
}
}
if (gPad && TestBit(kMustCleanup)) gPad->Modified();
}
void TH1::SmoothArray(Int_t nn, Double_t *xx, Int_t ntimes)
{
Int_t ii, jj, ik, jk, kk, nn2;
Double_t hh[6] = {0,0,0,0,0,0};
Double_t *yy = new Double_t[nn];
Double_t *zz = new Double_t[nn];
Double_t *rr = new Double_t[nn];
for (Int_t pass=0;pass<ntimes;pass++) {
for (ii = 0; ii < nn; ii++) {
yy[ii] = xx[ii];
}
for (kk = 1; kk <= 3; kk++) {
ik = 0;
if (kk == 2) ik = 1;
nn2 = nn - ik - 1;
for (ii = ik + 1; ii < nn2; ii++) {
for (jj = 0; jj < 3; jj++) {
hh[jj] = yy[ii + jj - 1];
}
zz[ii] = TMath::Median(3 + 2*ik, hh);
}
if (kk == 1) {
hh[0] = 3*yy[1] - 2*yy[2];
hh[1] = yy[0];
hh[2] = yy[1];
zz[0] = TMath::Median(3, hh);
hh[0] = yy[nn - 2];
hh[1] = yy[nn - 1];
hh[2] = 3*yy[nn - 2] - 2*yy[nn - 3];
zz[nn - 1] = TMath::Median(3, hh);
}
if (kk == 2) {
zz[0] = yy[0];
for (ii = 0; ii < 3; ii++) {
hh[ii] = yy[ii];
}
zz[1] = TMath::Median(3, hh);
for (ii = 0; ii < 3; ii++) {
hh[ii] = yy[nn - 3 + ii];
}
zz[nn - 2] = TMath::Median(3, hh);
zz[nn - 1] = yy[nn - 1];
}
}
for (ii = 2; ii < (nn - 2); ii++) {
if (zz[ii - 1] != zz[ii]) continue;
if (zz[ii] != zz[ii + 1]) continue;
hh[0] = zz[ii - 2] - zz[ii];
hh[1] = zz[ii + 2] - zz[ii];
if (hh[0] * hh[1] < 0) continue;
jk = 1;
if ( TMath::Abs(hh[1]) > TMath::Abs(hh[0]) ) jk = -1;
yy[ii] = -0.5*zz[ii - 2*jk] + zz[ii]/0.75 + zz[ii + 2*jk] /6.;
yy[ii + jk] = 0.5*(zz[ii + 2*jk] - zz[ii - 2*jk]) + zz[ii];
}
for (ii = 1; ii < nn - 1; ii++) {
rr[ii] = 0.25*yy[ii - 1] + 0.5*yy[ii] + 0.25*yy[ii + 1];
}
rr[0] = yy[0];
rr[nn - 1] = yy[nn - 1];
for (ii = 0; ii < nn; ii++) {
yy[ii] = xx[ii] - rr[ii];
}
for (kk = 1; kk <= 3; kk++) {
ik = 0;
if (kk == 2) ik = 1;
nn2 = nn - ik - 1;
for (ii = ik + 1; ii < nn2; ii++) {
for (jj = 0; jj < 3; jj++) {
hh[jj] = yy[ii + jj - 1];
}
zz[ii] = TMath::Median(3 + 2*ik, hh);
}
if (kk == 1) {
hh[0] = 3*yy[1] - 2*yy[2];
hh[1] = yy[0];
hh[2] = yy[1];
zz[0] = TMath::Median(3, hh);
hh[0] = yy[nn - 2];
hh[1] = yy[nn - 1];
hh[2] = 3*yy[nn - 2] - 2*yy[nn - 3];
zz[nn - 1] = TMath::Median(3, hh);
}
if (kk == 2) {
zz[0] = yy[0];
for (ii = 0; ii < 3; ii++) {
hh[ii] = yy[ii];
}
zz[1] = TMath::Median(3, hh);
for (ii = 0; ii < 3; ii++) {
hh[ii] = yy[nn - 3 + ii];
}
zz[nn - 2] = TMath::Median(3, hh);
zz[nn - 1] = yy[nn - 1];
}
}
for (ii = 2; ii < (nn - 2); ii++) {
if (zz[ii - 1] != zz[ii]) continue;
if (zz[ii] != zz[ii + 1]) continue;
hh[0] = zz[ii - 2] - zz[ii];
hh[1] = zz[ii + 2] - zz[ii];
if (hh[0] * hh[1] < 0) continue;
jk = 1;
if ( TMath::Abs(hh[1]) > TMath::Abs(hh[0]) ) jk = -1;
yy[ii] = -0.5*zz[ii - 2*jk] + zz[ii]/0.75 + zz[ii + 2*jk]/6.;
yy[ii + jk] = 0.5*(zz[ii + 2*jk] - zz[ii - 2*jk]) + zz[ii];
}
for (ii = 1; ii < (nn - 1); ii++) {
zz[ii] = 0.25*yy[ii - 1] + 0.5*yy[ii] + 0.25*yy[ii + 1];
}
zz[0] = yy[0];
zz[nn - 1] = yy[nn - 1];
for (ii = 0; ii < nn; ii++) {
if (xx[ii] < 0) xx[ii] = rr[ii] + zz[ii];
else xx[ii] = TMath::Abs(rr[ii] + zz[ii]);
}
}
delete [] yy;
delete [] zz;
delete [] rr;
}
void TH1::Smooth(Int_t ntimes, Int_t firstbin, Int_t lastbin)
{
if (fDimension != 1) {
Error("Smooth","Smooth only supported for 1-d histograms");
return;
}
Int_t nbins = fXaxis.GetNbins();
if (firstbin < 0) firstbin = 1;
if (lastbin < 0) lastbin = nbins;
if (lastbin > nbins+1) lastbin = nbins;
nbins = lastbin - firstbin + 1;
Double_t *xx = new Double_t[nbins];
Double_t nent = fEntries;
Int_t i;
for (i=0;i<nbins;i++) {
xx[i] = GetBinContent(i+firstbin);
}
TH1::SmoothArray(nbins,xx,ntimes);
for (i=0;i<nbins;i++) {
SetBinContent(i+firstbin,xx[i]);
}
fEntries = nent;
delete [] xx;
if (gPad) gPad->Modified();
}
void TH1::StatOverflows(Bool_t flag)
{
fgStatOverflows = flag;
}
void TH1::Streamer(TBuffer &b)
{
if (b.IsReading()) {
UInt_t R__s, R__c;
Version_t R__v = b.ReadVersion(&R__s, &R__c);
if (R__v > 2) {
TH1::Class()->ReadBuffer(b, this, R__v, R__s, R__c);
fXaxis.SetParent(this);
fYaxis.SetParent(this);
fZaxis.SetParent(this);
if (fgAddDirectory && !gROOT->ReadingObject()) {
fDirectory = gDirectory;
if (!gDirectory->GetList()->FindObject(this)) gDirectory->Append(this);
}
ResetBit(kCanDelete);
return;
}
TNamed::Streamer(b);
TAttLine::Streamer(b);
TAttFill::Streamer(b);
TAttMarker::Streamer(b);
b >> fNcells;
fXaxis.Streamer(b);
fYaxis.Streamer(b);
fZaxis.Streamer(b);
fXaxis.SetParent(this);
fYaxis.SetParent(this);
fZaxis.SetParent(this);
b >> fBarOffset;
b >> fBarWidth;
b >> fEntries;
b >> fTsumw;
b >> fTsumw2;
b >> fTsumwx;
b >> fTsumwx2;
if (R__v < 2) {
Float_t maximum, minimum, norm;
Float_t *contour=0;
b >> maximum; fMaximum = maximum;
b >> minimum; fMinimum = minimum;
b >> norm; fNormFactor = norm;
Int_t n = b.ReadArray(contour);
fContour.Set(n);
for (Int_t i=0;i<n;i++) fContour.fArray[i] = contour[i];
delete [] contour;
} else {
b >> fMaximum;
b >> fMinimum;
b >> fNormFactor;
fContour.Streamer(b);
}
fSumw2.Streamer(b);
fOption.Streamer(b);
fFunctions->Delete();
fFunctions->Streamer(b);
if (!gROOT->ReadingObject()) {
fDirectory = gDirectory;
if (!gDirectory->GetList()->FindObject(this)) gDirectory->Append(this);
}
b.CheckByteCount(R__s, R__c, TH1::IsA());
} else {
TH1::Class()->WriteBuffer(b,this);
}
}
void TH1::Print(Option_t *option) const
{
printf( "TH1.Print Name = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
TString opt = option;
opt.ToLower();
Int_t all;
if (opt.Contains("all")) all = 0;
else if (opt.Contains("range")) all = 1;
else if (opt.Contains("base")) all = 2;
else return;
Int_t bin, binx, biny, binz;
Int_t firstx=0,lastx=0,firsty=0,lasty=0,firstz=0,lastz=0;
if (all == 0) {
lastx = fXaxis.GetNbins()+1;
if (fDimension > 1) lasty = fYaxis.GetNbins()+1;
if (fDimension > 2) lastz = fZaxis.GetNbins()+1;
} else {
firstx = fXaxis.GetFirst(); lastx = fXaxis.GetLast();
if (fDimension > 1) {firsty = fYaxis.GetFirst(); lasty = fYaxis.GetLast();}
if (fDimension > 2) {firstz = fZaxis.GetFirst(); lastz = fZaxis.GetLast();}
}
if (all== 2) {
printf(" Title = %s\n", GetTitle());
printf(" NbinsX= %d, xmin= %g, xmax=%g", fXaxis.GetNbins(), fXaxis.GetXmin(), fXaxis.GetXmax());
if( fDimension > 1) printf(", NbinsY= %d, ymin= %g, ymax=%g", fYaxis.GetNbins(), fYaxis.GetXmin(), fYaxis.GetXmax());
if( fDimension > 2) printf(", NbinsZ= %d, zmin= %g, zmax=%g", fZaxis.GetNbins(), fZaxis.GetXmin(), fZaxis.GetXmax());
printf("\n");
return;
}
Double_t w,e;
Double_t x,y,z;
if (fDimension == 1) {
for (binx=firstx;binx<=lastx;binx++) {
x = fXaxis.GetBinCenter(binx);
w = GetBinContent(binx);
e = GetBinError(binx);
if(fSumw2.fN) printf(" fSumw[%d]=%g, x=%g, error=%g\n",binx,w,x,e);
else printf(" fSumw[%d]=%g, x=%g\n",binx,w,x);
}
}
if (fDimension == 2) {
for (biny=firsty;biny<=lasty;biny++) {
y = fYaxis.GetBinCenter(biny);
for (binx=firstx;binx<=lastx;binx++) {
bin = GetBin(binx,biny);
x = fXaxis.GetBinCenter(binx);
w = GetBinContent(bin);
e = GetBinError(bin);
if(fSumw2.fN) printf(" fSumw[%d][%d]=%g, x=%g, y=%g, error=%g\n",binx,biny,w,x,y,e);
else printf(" fSumw[%d][%d]=%g, x=%g, y=%g\n",binx,biny,w,x,y);
}
}
}
if (fDimension == 3) {
for (binz=firstz;binz<=lastz;binz++) {
z = fZaxis.GetBinCenter(binz);
for (biny=firsty;biny<=lasty;biny++) {
y = fYaxis.GetBinCenter(biny);
for (binx=firstx;binx<=lastx;binx++) {
bin = GetBin(binx,biny,binz);
x = fXaxis.GetBinCenter(binx);
w = GetBinContent(bin);
e = GetBinError(bin);
if(fSumw2.fN) printf(" fSumw[%d][%d][%d]=%g, x=%g, y=%g, z=%g, error=%g\n",binx,biny,binz,w,x,y,z,e);
else printf(" fSumw[%d][%d][%d]=%g, x=%g, y=%g, z=%g\n",binx,biny,binz,w,x,y,z);
}
}
}
}
}
void TH1::Rebuild(Option_t *)
{
SetBinsLength();
if (fSumw2.fN) {
fSumw2.Set(fNcells);
}
}
void TH1::Reset(Option_t *option)
{
TString opt = option;
opt.ToUpper();
fSumw2.Reset();
if (fIntegral) {delete [] fIntegral; fIntegral = 0;}
if (opt.Contains("ICE")) return;
if (fBuffer) BufferEmpty();
fTsumw = 0;
fTsumw2 = 0;
fTsumwx = 0;
fTsumwx2 = 0;
fEntries = 0;
TObject *stats = fFunctions->FindObject("stats");
fFunctions->Remove(stats);
TObject *obj;
while ((obj = fFunctions->First())) {
while(fFunctions->Remove(obj));
delete obj;
}
if(stats) fFunctions->Add(stats);
fContour.Set(0);
}
void TH1::SavePrimitive(ostream &out, Option_t *option )
{
Bool_t nonEqiX = kFALSE;
Bool_t nonEqiY = kFALSE;
Bool_t nonEqiZ = kFALSE;
Int_t i;
if (GetXaxis()->GetXbins()->fN && GetXaxis()->GetXbins()->fArray) {
nonEqiX = kTRUE;
out << " Double_t xAxis[" << GetXaxis()->GetXbins()->fN
<< "] = {";
for (i = 0; i < GetXaxis()->GetXbins()->fN; i++) {
if (i != 0) out << ", ";
out << GetXaxis()->GetXbins()->fArray[i];
}
out << "}; " << endl;
}
if (fDimension > 1 && GetYaxis()->GetXbins()->fN &&
GetYaxis()->GetXbins()->fArray) {
nonEqiY = kTRUE;
out << " Double_t yAxis[" << GetYaxis()->GetXbins()->fN
<< "] = {";
for (i = 0; i < GetYaxis()->GetXbins()->fN; i++) {
if (i != 0) out << ", ";
out << GetYaxis()->GetXbins()->fArray[i];
}
out << "}; " << endl;
}
if (fDimension > 2 && GetZaxis()->GetXbins()->fN &&
GetZaxis()->GetXbins()->fArray) {
nonEqiZ = kTRUE;
out << " Double_t zAxis[" << GetZaxis()->GetXbins()->fN
<< "] = {";
for (i = 0; i < GetZaxis()->GetXbins()->fN; i++) {
if (i != 0) out << ", ";
out << GetZaxis()->GetXbins()->fArray[i];
}
out << "}; " << endl;
}
char quote = '"';
out <<" "<<endl;
out <<" "<<"TH1"<<" *";
out << GetName() << " = new " << ClassName() << "(" << quote
<< GetName() << quote << "," << quote<< GetTitle() << quote
<< "," << GetXaxis()->GetNbins();
if (nonEqiX)
out << ", xAxis";
else
out << "," << GetXaxis()->GetXmin()
<< "," << GetXaxis()->GetXmax();
if (fDimension > 1) {
out << "," << GetYaxis()->GetNbins();
if (nonEqiY)
out << ", yAxis";
else
out << "," << GetYaxis()->GetXmin()
<< "," << GetYaxis()->GetXmax();
}
if (fDimension > 2) {
out << "," << GetZaxis()->GetNbins();
if (nonEqiZ)
out << ", zAxis";
else
out << "," << GetZaxis()->GetXmin()
<< "," << GetZaxis()->GetXmax();
}
out << ");" << endl;
Int_t bin;
for (bin=0;bin<fNcells;bin++) {
Double_t bc = GetBinContent(bin);
if (bc) {
out<<" "<<GetName()<<"->SetBinContent("<<bin<<","<<bc<<");"<<endl;
}
}
if (fSumw2.fN) {
for (bin=0;bin<fNcells;bin++) {
Double_t be = GetBinError(bin);
if (be) {
out<<" "<<GetName()<<"->SetBinError("<<bin<<","<<be<<");"<<endl;
}
}
}
TH1::SavePrimitiveHelp(out, option);
}
void TH1::SavePrimitiveHelp(ostream &out, Option_t *option )
{
char quote = '"';
if (TMath::Abs(GetBarOffset()) > 1e-5) {
out<<" "<<GetName()<<"->SetBarOffset("<<GetBarOffset()<<");"<<endl;
}
if (TMath::Abs(GetBarWidth()-1) > 1e-5) {
out<<" "<<GetName()<<"->SetBarWidth("<<GetBarWidth()<<");"<<endl;
}
if (fMinimum != -1111) {
out<<" "<<GetName()<<"->SetMinimum("<<fMinimum<<");"<<endl;
}
if (fMaximum != -1111) {
out<<" "<<GetName()<<"->SetMaximum("<<fMaximum<<");"<<endl;
}
if (fNormFactor != 0) {
out<<" "<<GetName()<<"->SetNormFactor("<<fNormFactor<<");"<<endl;
}
if (fEntries != 0) {
out<<" "<<GetName()<<"->SetEntries("<<fEntries<<");"<<endl;
}
if (fDirectory == 0) {
out<<" "<<GetName()<<"->SetDirectory(0);"<<endl;
}
if (TestBit(kNoStats)) {
out<<" "<<GetName()<<"->SetStats(0);"<<endl;
}
if (fOption.Length() != 0) {
out<<" "<<GetName()<<"->SetOption("<<quote<<fOption.Data()<<quote<<");"<<endl;
}
Int_t ncontours = GetContour();
if (ncontours > 0) {
out<<" "<<GetName()<<"->SetContour("<<ncontours<<");"<<endl;
for (Int_t bin=0;bin<ncontours;bin++) {
out<<" "<<GetName()<<"->SetContourLevel("<<bin<<","<<GetContourLevel(bin)<<");"<<endl;
}
}
TObjOptLink *lnk = (TObjOptLink*)fFunctions->FirstLink();
TObject *obj;
while (lnk) {
obj = lnk->GetObject();
obj->SavePrimitive(out,"nodraw");
if (obj->InheritsFrom("TF1")) {
out<<" "<<GetName()<<"->GetListOfFunctions()->Add("<<obj->GetName()<<");"<<endl;
} else if (obj->InheritsFrom("TPaveStats")) {
out<<" "<<GetName()<<"->GetListOfFunctions()->Add(ptstats);"<<endl;
out<<" ptstats->SetParent("<<GetName()<<"->GetListOfFunctions());"<<endl;
} else {
out<<" "<<GetName()<<"->GetListOfFunctions()->Add("<<obj->GetName()<<","<<quote<<lnk->GetOption()<<quote<<");"<<endl;
}
lnk = (TObjOptLink*)lnk->Next();
}
SaveFillAttributes(out,GetName(),0,1001);
SaveLineAttributes(out,GetName(),1,1,1);
SaveMarkerAttributes(out,GetName(),1,1,1);
fXaxis.SaveAttributes(out,GetName(),"->GetXaxis()");
fYaxis.SaveAttributes(out,GetName(),"->GetYaxis()");
fZaxis.SaveAttributes(out,GetName(),"->GetZaxis()");
TString opt = option;
opt.ToLower();
if (!opt.Contains("nodraw")) {
out<<" "<<GetName()<<"->Draw("
<<quote<<option<<quote<<");"<<endl;
}
}
void TH1::UseCurrentStyle()
{
if (gStyle->IsReading()) {
fXaxis.ResetAttAxis("X");
fYaxis.ResetAttAxis("Y");
fZaxis.ResetAttAxis("Z");
SetBarOffset(gStyle->GetBarOffset());
SetBarWidth(gStyle->GetBarWidth());
SetFillColor(gStyle->GetHistFillColor());
SetFillStyle(gStyle->GetHistFillStyle());
SetLineColor(gStyle->GetHistLineColor());
SetLineStyle(gStyle->GetHistLineStyle());
SetLineWidth(gStyle->GetHistLineWidth());
SetMarkerColor(gStyle->GetMarkerColor());
SetMarkerStyle(gStyle->GetMarkerStyle());
SetMarkerSize(gStyle->GetMarkerSize());
Int_t dostat = gStyle->GetOptStat();
if (gStyle->GetOptFit() && !dostat) dostat = 1000000001;
SetStats(dostat);
} else {
gStyle->SetBarOffset(fBarOffset);
gStyle->SetBarWidth(fBarWidth);
gStyle->SetHistFillColor(GetFillColor());
gStyle->SetHistFillStyle(GetFillStyle());
gStyle->SetHistLineColor(GetLineColor());
gStyle->SetHistLineStyle(GetLineStyle());
gStyle->SetHistLineWidth(GetLineWidth());
gStyle->SetMarkerColor(GetMarkerColor());
gStyle->SetMarkerStyle(GetMarkerStyle());
gStyle->SetMarkerSize(GetMarkerSize());
gStyle->SetOptStat(TestBit(kNoStats));
}
TIter next(GetListOfFunctions());
TObject *obj;
while ((obj = next())) {
obj->UseCurrentStyle();
}
}
Double_t TH1::GetMean(Int_t axis) const
{
if (axis<1 || axis>3&&axis<11 || axis>13) return 0;
Double_t stats[kNstat];
for (Int_t i=4;i<kNstat;i++) stats[i] = 0;
GetStats(stats);
if (stats[0] == 0) return 0;
if (axis<10){
Int_t ax[3] = {2,4,7};
return stats[ax[axis-1]]/stats[0];
} else {
Double_t rms = GetRMS(axis-10);
Double_t neff = GetEffectiveEntries();
return ( neff > 0 ? rms/TMath::Sqrt(neff) : 0. );
}
}
Double_t TH1::GetMeanError(Int_t axis) const
{
return GetMean(axis+10);
}
Double_t TH1::GetRMS(Int_t axis) const
{
if (axis<1 || axis>3&&axis<11 || axis>13) return 0;
Double_t x, rms2, stats[kNstat];
for (Int_t i=4;i<kNstat;i++) stats[i] = 0;
GetStats(stats);
if (stats[0] == 0) return 0;
Int_t ax[3] = {2,4,7};
Int_t axm = ax[axis%10 - 1];
x = stats[axm]/stats[0];
rms2 = TMath::Abs(stats[axm+1]/stats[0] -x*x);
if (axis<10)
return TMath::Sqrt(rms2);
else {
Double_t neff = GetEffectiveEntries();
return ( neff > 0 ? TMath::Sqrt(rms2/(2*neff) ) : 0. );
}
}
Double_t TH1::GetRMSError(Int_t axis) const
{
return GetRMS(axis+10);
}
Double_t TH1::GetSkewness(Int_t axis) const
{
const TAxis *ax;
if (axis==1 || axis==11) ax = &fXaxis;
else if (axis==2 || axis==12) ax = &fYaxis;
else if (axis==3 || axis==13) ax = &fZaxis;
else {
Error("GetSkewness", "illegal value of parameter");
return 0;
}
if (axis < 10) {
Double_t x, w, mean, rms, rms3, sum=0;
mean = GetMean(axis);
rms = GetRMS(axis);
rms3 = rms*rms*rms;
Int_t bin;
Double_t np=0;
for (bin=ax->GetFirst(); bin<=ax->GetLast(); bin++){
x = GetBinCenter(bin);
w = GetBinContent(bin);
np+=w;
sum+=w*(x-mean)*(x-mean)*(x-mean);
}
sum/=np*rms3;
return sum;
} else {
Int_t nbins = ax->GetNbins();
return TMath::Sqrt(6./nbins);
}
}
Double_t TH1::GetKurtosis(Int_t axis) const
{
const TAxis *ax;
if (axis==1 || axis==11) ax = &fXaxis;
else if (axis==2 || axis==12) ax = &fYaxis;
else if (axis==3 || axis==13) ax = &fZaxis;
else {
Error("GetKurtosis", "illegal value of parameter");
return 0;
}
if (axis < 10){
Double_t x, w, mean, rms, rms4, sum=0;
mean = GetMean(axis);
rms = GetRMS(axis);
rms4 = rms*rms*rms*rms;
Int_t bin;
Double_t np=0;
for (bin=ax->GetFirst(); bin<=ax->GetLast(); bin++){
x = GetBinCenter(bin);
w = GetBinContent(bin);
np+=w;
sum+=w*(x-mean)*(x-mean)*(x-mean)*(x-mean);
}
sum/=np*rms4;
return sum-3;
} else {
Int_t nbins = ax->GetNbins();
return TMath::Sqrt(24./nbins);
}
}
void TH1::GetStats(Double_t *stats) const
{
if (fBuffer) ((TH1*)this)->BufferEmpty();
Int_t bin, binx;
Double_t w,err;
Double_t x;
if (fTsumw == 0 || fXaxis.TestBit(TAxis::kAxisRange)) {
for (bin=0;bin<4;bin++) stats[bin] = 0;
Int_t firstBinX = fXaxis.GetFirst();
Int_t lastBinX = fXaxis.GetLast();
if (fgStatOverflows && !fXaxis.TestBit(TAxis::kAxisRange)) {
if (firstBinX == 1) firstBinX = 0;
if (lastBinX == fXaxis.GetNbins() ) lastBinX += 1;
}
for (binx = firstBinX; binx <= lastBinX; binx++) {
x = fXaxis.GetBinCenter(binx);
w = TMath::Abs(GetBinContent(binx));
err = TMath::Abs(GetBinError(binx));
stats[0] += w;
stats[1] += err*err;
stats[2] += w*x;
stats[3] += w*x*x;
}
} else {
stats[0] = fTsumw;
stats[1] = fTsumw2;
stats[2] = fTsumwx;
stats[3] = fTsumwx2;
}
}
void TH1::PutStats(Double_t *stats)
{
fTsumw = stats[0];
fTsumw2 = stats[1];
fTsumwx = stats[2];
fTsumwx2 = stats[3];
}
Double_t TH1::GetSumOfWeights() const
{
Int_t bin,binx,biny,binz;
Double_t sum =0;
for(binz=1; binz<=fZaxis.GetNbins(); binz++) {
for(biny=1; biny<=fYaxis.GetNbins(); biny++) {
for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
bin = GetBin(binx,biny,binz);
sum += GetBinContent(bin);
}
}
}
return sum;
}
Double_t TH1::Integral(Option_t *option) const
{
return Integral(fXaxis.GetFirst(),fXaxis.GetLast(),option);
}
Double_t TH1::Integral(Int_t binx1, Int_t binx2, Option_t *option) const
{
Int_t nbinsx = GetNbinsX();
if (binx1 < 0) binx1 = 0;
if (binx2 > nbinsx+1) binx2 = nbinsx+1;
if (binx2 < binx1) binx2 = nbinsx;
Double_t integral = 0;
TString opt = option;
opt.ToLower();
Bool_t width = kFALSE;
if (opt.Contains("width")) width = kTRUE;
Int_t binx;
for (binx=binx1;binx<=binx2;binx++) {
if (width) integral += GetBinContent(binx)*fXaxis.GetBinWidth(binx);
else integral += GetBinContent(binx);
}
return integral;
}
Double_t TH1::KolmogorovTest(const TH1 *h2, Option_t *option) const
{
TString opt = option;
opt.ToUpper();
Double_t prob = 0;
TH1 *h1 = (TH1*)this;
if (h2 == 0) return 0;
TAxis *axis1 = h1->GetXaxis();
TAxis *axis2 = h2->GetXaxis();
Int_t ncx1 = axis1->GetNbins();
Int_t ncx2 = axis2->GetNbins();
if (h1->GetDimension() != 1 || h2->GetDimension() != 1) {
Error("KolmogorovTest","Histograms must be 1-D\n");
return 0;
}
if (ncx1 != ncx2) {
Error("KolmogorovTest","Number of channels is different, %d and %d\n",ncx1,ncx2);
return 0;
}
Double_t difprec = 1e-5;
Double_t diff1 = TMath::Abs(axis1->GetXmin() - axis2->GetXmin());
Double_t diff2 = TMath::Abs(axis1->GetXmax() - axis2->GetXmax());
if (diff1 > difprec || diff2 > difprec) {
Error("KolmogorovTest","histograms with different binning");
return 0;
}
Bool_t afunc1 = kFALSE;
Bool_t afunc2 = kFALSE;
Double_t sum1 = 0, sum2 = 0;
Double_t ew1, ew2, w1 = 0, w2 = 0;
Int_t bin;
for (bin=1;bin<=ncx1;bin++) {
sum1 += h1->GetBinContent(bin);
sum2 += h2->GetBinContent(bin);
ew1 = h1->GetBinError(bin);
ew2 = h2->GetBinError(bin);
w1 += ew1*ew1;
w2 += ew2*ew2;
}
if (sum1 == 0) {
Error("KolmogorovTest","Histogram1 %s integral is zero\n",h1->GetName());
return 0;
}
if (sum2 == 0) {
Error("KolmogorovTest","Histogram2 %s integral is zero\n",h2->GetName());
return 0;
}
Double_t tsum1 = sum1;
Double_t tsum2 = sum2;
if (opt.Contains("U")) {
tsum1 += h1->GetBinContent(0);
tsum2 += h2->GetBinContent(0);
}
if (opt.Contains("O")) {
tsum1 += h1->GetBinContent(ncx1+1);
tsum2 += h2->GetBinContent(ncx1+1);
}
Double_t ne1 = h1->GetEntries();
Double_t ne2 = h2->GetEntries();
Double_t difsum1 = (ne1-tsum1)/tsum1;
Double_t esum1 = sum1;
if (difsum1 > difprec && Int_t(ne1) != ncx1) {
if (opt.Contains("U") || opt.Contains("O")) {
Warning("KolmogorovTest","U/O option with weighted events for hist:%s\n",h1->GetName());
}
if (h1->GetSumw2N() == 0) {
Warning("KolmogorovTest","Weighted events and no Sumw2, hist:%s\n",h1->GetName());
} else {
esum1 = sum1*sum1/w1;
}
}
Double_t difsum2 = (ne2-tsum2)/tsum2;
Double_t esum2 = sum2;
if (difsum2 > difprec && Int_t(ne2) != ncx1) {
if (opt.Contains("U") || opt.Contains("O")) {
Warning("KolmogorovTest","U/O option with weighted events for hist:%s\n",h2->GetName());
}
if (h2->GetSumw2N() == 0) {
Warning("KolmogorovTest","Weighted events and no Sumw2, hist:%s\n",h2->GetName());
} else {
esum2 = sum2*sum2/w2;
}
}
Double_t s1 = 1/tsum1;
Double_t s2 = 1/tsum2;
Double_t dfmax =0, rsum1 = 0, rsum2 = 0;
Int_t first = 1;
Int_t last = ncx1;
if (opt.Contains("U")) first = 0;
if (opt.Contains("O")) last = ncx1+1;
for (bin=first;bin<=last;bin++) {
rsum1 += s1*h1->GetBinContent(bin);
rsum2 += s2*h2->GetBinContent(bin);
dfmax = TMath::Max(dfmax,TMath::Abs(rsum1-rsum2));
}
Double_t z, prb1=0, prb2=0, prb3=0;
if (afunc1) z = dfmax*TMath::Sqrt(esum2);
else if (afunc2) z = dfmax*TMath::Sqrt(esum1);
else z = dfmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
prob = TMath::KolmogorovProb(z);
if (opt.Contains("N")) {
prb1 = prob;
Double_t resum1 = esum1; if (afunc1) resum1 = 0;
Double_t resum2 = esum2; if (afunc2) resum2 = 0;
Double_t d12 = esum1-esum2;
Double_t chi2 = d12*d12/(resum1+resum2);
prb2 = TMath::Prob(chi2,1);
if (prob > 0 && prb2 > 0) prob *= prb2*(1-TMath::Log(prob*prb2));
else prob = 0;
}
const Int_t nEXPT = 1000;
if (opt.Contains("X")) {
Double_t dSEXPT;
Bool_t addStatus = fgAddDirectory;
fgAddDirectory = kFALSE;
TH1F *hDistValues = new TH1F("hDistValues","KS distances",200,0,1);
TH1 *hExpt = (TH1*)Clone();
fgAddDirectory = addStatus;
for (Int_t i=0; i < nEXPT; i++) {
hExpt->Reset();
hExpt->FillRandom(h1,(Int_t)ne2);
dSEXPT = KolmogorovTest(hExpt,"M");
hDistValues->Fill(dSEXPT);
}
prb3 = hDistValues->Integral(hDistValues->FindBin(dfmax),200)/hDistValues->Integral();
delete hDistValues;
delete hExpt;
}
if (opt.Contains("D")) {
printf(" Kolmo Prob h1 = %s, sum1=%g\n",h1->GetName(),sum1);
printf(" Kolmo Prob h2 = %s, sum2=%g\n",h2->GetName(),sum2);
printf(" Kolmo Prob = %g, Max Dist = %g\n",prob,dfmax);
if (opt.Contains("N"))
printf(" Kolmo Prob = %f for shape alone, =%f for normalisation alone\n",prb1,prb2);
if (opt.Contains("X"))
printf(" Kolmo Prob = %f with %d pseudo-experiments\n",prb3,nEXPT);
}
if (TMath::Abs(rsum1-1) > 0.002) Warning("KolmogorovTest","Numerical problems with h1=%s\n",h1->GetName());
if (TMath::Abs(rsum2-1) > 0.002) Warning("KolmogorovTest","Numerical problems with h2=%s\n",h2->GetName());
if(opt.Contains("M")) return dfmax;
else if(opt.Contains("X")) return prb3;
else return prob;
}
void TH1::SetContent(const Double_t *content)
{
Int_t bin;
Double_t bincontent;
for (bin=0; bin<fNcells; bin++) {
bincontent = *(content + bin);
SetBinContent(bin, bincontent);
}
}
Int_t TH1::GetContour(Double_t *levels)
{
Int_t nlevels = fContour.fN;
if (levels) {
if (nlevels == 0) {
nlevels = 20;
SetContour(nlevels);
} else {
if (TestBit(kUserContour) == 0) SetContour(nlevels);
}
for (Int_t level=0; level<nlevels; level++) levels[level] = fContour.fArray[level];
}
return nlevels;
}
Double_t TH1::GetContourLevel(Int_t level) const
{
if (level <0 || level >= fContour.fN) return 0;
Double_t zlevel = fContour.fArray[level];
return zlevel;
}
Double_t TH1::GetContourLevelPad(Int_t level) const
{
if (level <0 || level >= fContour.fN) return 0;
Double_t zlevel = fContour.fArray[level];
if (gPad && gPad->GetLogz() && TestBit(kUserContour)) {
if (zlevel <= 0) return 0;
zlevel = TMath::Log10(zlevel);
}
return zlevel;
}
void TH1::SetBuffer(Int_t buffersize, Option_t * )
{
if (fBuffer) {
BufferEmpty();
delete [] fBuffer;
fBuffer = 0;
}
if (buffersize <= 0) {
fBufferSize = 0;
return;
}
if (buffersize < 100) buffersize = 100;
fBufferSize = 1 + buffersize*(fDimension+1);
fBuffer = new Double_t[fBufferSize];
memset(fBuffer,0,8*fBufferSize);
}
void TH1::SetContour(Int_t nlevels, const Double_t *levels)
{
Int_t level;
ResetBit(kUserContour);
if (nlevels <=0 ) {
fContour.Set(0);
return;
}
fContour.Set(nlevels);
if (levels) {
SetBit(kUserContour);
for (level=0; level<nlevels; level++) fContour.fArray[level] = levels[level];
} else {
Double_t zmin = GetMinimum();
Double_t zmax = GetMaximum();
if ((zmin == zmax) && (zmin != 0)) {
zmax += 0.01*TMath::Abs(zmax);
zmin -= 0.01*TMath::Abs(zmin);
}
Double_t dz = (zmax-zmin)/Double_t(nlevels);
if (gPad && gPad->GetLogz()) {
if (zmax <= 0) return;
if (zmin <= 0) zmin = 0.001*zmax;
zmin = TMath::Log10(zmin);
zmax = TMath::Log10(zmax);
dz = (zmax-zmin)/Double_t(nlevels);
}
for (level=0; level<nlevels; level++) {
fContour.fArray[level] = zmin + dz*Double_t(level);
}
}
}
void TH1::SetContourLevel(Int_t level, Double_t value)
{
if (level <0 || level >= fContour.fN) return;
SetBit(kUserContour);
fContour.fArray[level] = value;
}
Double_t TH1::GetMaximum(Double_t maxval) const
{
if (fMaximum != -1111) return fMaximum;
Int_t bin, binx, biny, binz;
Int_t xfirst = fXaxis.GetFirst();
Int_t xlast = fXaxis.GetLast();
Int_t yfirst = fYaxis.GetFirst();
Int_t ylast = fYaxis.GetLast();
Int_t zfirst = fZaxis.GetFirst();
Int_t zlast = fZaxis.GetLast();
Double_t maximum = -FLT_MAX, value;
for (binz=zfirst;binz<=zlast;binz++) {
for (biny=yfirst;biny<=ylast;biny++) {
for (binx=xfirst;binx<=xlast;binx++) {
bin = GetBin(binx,biny,binz);
value = GetBinContent(bin);
if (value > maximum && value < maxval) maximum = value;
}
}
}
return maximum;
}
Int_t TH1::GetMaximumBin() const
{
Int_t locmax, locmay, locmaz;
return GetMaximumBin(locmax, locmay, locmaz);
}
Int_t TH1::GetMaximumBin(Int_t &locmax, Int_t &locmay, Int_t &locmaz) const
{
Int_t bin, binx, biny, binz;
Int_t locm;
Int_t xfirst = fXaxis.GetFirst();
Int_t xlast = fXaxis.GetLast();
Int_t yfirst = fYaxis.GetFirst();
Int_t ylast = fYaxis.GetLast();
Int_t zfirst = fZaxis.GetFirst();
Int_t zlast = fZaxis.GetLast();
Double_t maximum = -FLT_MAX, value;
locm = locmax = locmay = locmaz = 0;
for (binz=zfirst;binz<=zlast;binz++) {
for (biny=yfirst;biny<=ylast;biny++) {
for (binx=xfirst;binx<=xlast;binx++) {
bin = GetBin(binx,biny,binz);
value = GetBinContent(bin);
if (value > maximum) {
maximum = value;
locm = bin;
locmax = binx;
locmay = biny;
locmaz = binz;
}
}
}
}
return locm;
}
Double_t TH1::GetMinimum(Double_t minval) const
{
if (fMinimum != -1111) return fMinimum;
Int_t bin, binx, biny, binz;
Int_t xfirst = fXaxis.GetFirst();
Int_t xlast = fXaxis.GetLast();
Int_t yfirst = fYaxis.GetFirst();
Int_t ylast = fYaxis.GetLast();
Int_t zfirst = fZaxis.GetFirst();
Int_t zlast = fZaxis.GetLast();
Double_t minimum=FLT_MAX, value;
for (binz=zfirst;binz<=zlast;binz++) {
for (biny=yfirst;biny<=ylast;biny++) {
for (binx=xfirst;binx<=xlast;binx++) {
bin = GetBin(binx,biny,binz);
value = GetBinContent(bin);
if (value < minimum && value > minval) minimum = value;
}
}
}
return minimum;
}
Int_t TH1::GetMinimumBin() const
{
Int_t locmix, locmiy, locmiz;
return GetMinimumBin(locmix, locmiy, locmiz);
}
Int_t TH1::GetMinimumBin(Int_t &locmix, Int_t &locmiy, Int_t &locmiz) const
{
Int_t bin, binx, biny, binz;
Int_t locm;
Int_t xfirst = fXaxis.GetFirst();
Int_t xlast = fXaxis.GetLast();
Int_t yfirst = fYaxis.GetFirst();
Int_t ylast = fYaxis.GetLast();
Int_t zfirst = fZaxis.GetFirst();
Int_t zlast = fZaxis.GetLast();
Double_t minimum = FLT_MAX, value;
locm = locmix = locmiy = locmiz = 0;
for (binz=zfirst;binz<=zlast;binz++) {
for (biny=yfirst;biny<=ylast;biny++) {
for (binx=xfirst;binx<=xlast;binx++) {
bin = GetBin(binx,biny,binz);
value = GetBinContent(bin);
if (value < minimum) {
minimum = value;
locm = bin;
locmix = binx;
locmiy = biny;
locmiz = binz;
}
}
}
}
return locm;
}
void TH1::SetBins(Int_t nx, Double_t xmin, Double_t xmax)
{
if (GetDimension() != 1) {
Error("SetBins","Operation only valid for 1-d histograms");
return;
}
fXaxis.SetRange(0,0);
fXaxis.Set(nx,xmin,xmax);
fYaxis.Set(1,0,1);
fZaxis.Set(1,0,1);
fNcells = nx+2;
SetBinsLength(fNcells);
if (fSumw2.fN) {
fSumw2.Set(fNcells);
}
}
void TH1::SetBins(Int_t nx, const Double_t *xBins)
{
if (GetDimension() != 1) {
Error("SetBins","Operation only valid for 1-d histograms");
return;
}
fXaxis.SetRange(0,0);
fXaxis.Set(nx,xBins);
fYaxis.Set(1,0,1);
fZaxis.Set(1,0,1);
fNcells = nx+2;
SetBinsLength(fNcells);
if (fSumw2.fN) {
fSumw2.Set(fNcells);
}
}
void TH1::SetBins(Int_t nx, Double_t xmin, Double_t xmax, Int_t ny, Double_t ymin, Double_t ymax)
{
if (GetDimension() != 2) {
Error("SetBins","Operation only valid for 2-d histograms");
return;
}
fXaxis.SetRange(0,0);
fYaxis.SetRange(0,0);
fXaxis.Set(nx,xmin,xmax);
fYaxis.Set(ny,ymin,ymax);
fZaxis.Set(1,0,1);
fNcells = (nx+2)*(ny+2);
SetBinsLength(fNcells);
if (fSumw2.fN) {
fSumw2.Set(fNcells);
}
}
void TH1::SetBins(Int_t nx, const Double_t *xBins, Int_t ny, const Double_t *yBins)
{
if (GetDimension() != 2) {
Error("SetBins","Operation only valid for 2-d histograms");
return;
}
fXaxis.SetRange(0,0);
fYaxis.SetRange(0,0);
fXaxis.Set(nx,xBins);
fYaxis.Set(ny,yBins);
fZaxis.Set(1,0,1);
fNcells = (nx+2)*(ny+2);
SetBinsLength(fNcells);
if (fSumw2.fN) {
fSumw2.Set(fNcells);
}
}
void TH1::SetBins(Int_t nx, Double_t xmin, Double_t xmax, Int_t ny, Double_t ymin, Double_t ymax, Int_t nz, Double_t zmin, Double_t zmax)
{
if (GetDimension() != 3) {
Error("SetBins","Operation only valid for 3-d histograms");
return;
}
fXaxis.SetRange(0,0);
fYaxis.SetRange(0,0);
fZaxis.SetRange(0,0);
fXaxis.Set(nx,xmin,xmax);
fYaxis.Set(ny,ymin,ymax);
fZaxis.Set(nz,zmin,zmax);
fNcells = (nx+2)*(ny+2)*(nz+2);
SetBinsLength(fNcells);
if (fSumw2.fN) {
fSumw2.Set(fNcells);
}
}
void TH1::SetMaximum(Double_t maximum)
{
fMaximum = maximum;
}
void TH1::SetMinimum(Double_t minimum)
{
fMinimum = minimum;
}
void TH1::SetDirectory(TDirectory *dir)
{
if (fDirectory == dir) return;
if (fDirectory) fDirectory->GetList()->Remove(this);
fDirectory = dir;
if (fDirectory) fDirectory->GetList()->Add(this);
}
void TH1::SetError(const Double_t *error)
{
Int_t bin;
Double_t binerror;
for (bin=0; bin<fNcells; bin++) {
binerror = error[bin];
SetBinError(bin, binerror);
}
}
void TH1::SetName(const char *name)
{
if (fDirectory) fDirectory->GetList()->Remove(this);
fName = name;
if (fDirectory) fDirectory->GetList()->Add(this);
}
void TH1::SetNameTitle(const char *name, const char *title)
{
if (fDirectory) fDirectory->GetList()->Remove(this);
fName = name;
SetTitle(title);
if (fDirectory) fDirectory->GetList()->Add(this);
}
void TH1::SetStats(Bool_t stats)
{
ResetBit(kNoStats);
if (!stats) {
SetBit(kNoStats);
if (fFunctions) delete fFunctions->FindObject("stats");
}
}
void TH1::Sumw2()
{
if (!fgDefaultSumw2 && fSumw2.fN) {
Warning("Sumw2","Sum of squares of weights structure already created");
return;
}
fSumw2.Set(fNcells);
for (Int_t bin=0; bin<fNcells; bin++) {
fSumw2.fArray[bin] = GetBinContent(bin);
}
}
TF1 *TH1::GetFunction(const char *name) const
{
return (TF1*)fFunctions->FindObject(name);
}
Double_t TH1::GetBinError(Int_t bin) const
{
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
if (fBuffer) ((TH1*)this)->BufferEmpty();
if (fSumw2.fN) return TMath::Sqrt(fSumw2.fArray[bin]);
Double_t error2 = TMath::Abs(GetBinContent(bin));
return TMath::Sqrt(error2);
}
Double_t TH1::GetBinError(Int_t binx, Int_t biny) const
{
Int_t bin = GetBin(binx,biny);
return GetBinError(bin);
}
Double_t TH1::GetBinError(Int_t binx, Int_t biny, Int_t binz) const
{
Int_t bin = GetBin(binx,biny,binz);
return GetBinError(bin);
}
Double_t TH1::GetCellContent(Int_t binx, Int_t biny) const
{
Int_t bin = GetBin(binx,biny);
return GetBinContent(bin);
}
Double_t TH1::GetCellError(Int_t binx, Int_t biny) const
{
Int_t bin = GetBin(binx,biny);
return GetBinError(bin);
}
void TH1::SetBinError(Int_t bin, Double_t error)
{
if (!fSumw2.fN) Sumw2();
if (bin <0 || bin>= fSumw2.fN) return;
fSumw2.fArray[bin] = error*error;
}
void TH1::SetBinContent(Int_t binx, Int_t biny, Double_t content)
{
if (binx <0 || binx>fXaxis.GetNbins()+1) return;
if (biny <0 || biny>fYaxis.GetNbins()+1) return;
SetBinContent(biny*(fXaxis.GetNbins()+2) + binx,content);
}
void TH1::SetBinContent(Int_t binx, Int_t biny, Int_t binz, Double_t content)
{
if (binx <0 || binx>fXaxis.GetNbins()+1) return;
if (biny <0 || biny>fYaxis.GetNbins()+1) return;
if (binz <0 || binz>fZaxis.GetNbins()+1) return;
Int_t bin = GetBin(binx,biny,binz);
SetBinContent(bin,content);
}
void TH1::SetCellContent(Int_t binx, Int_t biny, Double_t content)
{
if (binx <0 || binx>fXaxis.GetNbins()+1) return;
if (biny <0 || biny>fYaxis.GetNbins()+1) return;
SetBinContent(biny*(fXaxis.GetNbins()+2) + binx,content);
}
void TH1::SetBinError(Int_t binx, Int_t biny, Double_t error)
{
if (binx <0 || binx>fXaxis.GetNbins()+1) return;
if (biny <0 || biny>fYaxis.GetNbins()+1) return;
SetBinError(biny*(fXaxis.GetNbins()+2) + binx,error);
}
void TH1::SetBinError(Int_t binx, Int_t biny, Int_t binz, Double_t error)
{
if (binx <0 || binx>fXaxis.GetNbins()+1) return;
if (biny <0 || biny>fYaxis.GetNbins()+1) return;
if (binz <0 || binz>fZaxis.GetNbins()+1) return;
Int_t bin = GetBin(binx,biny,binz);
SetBinError(bin,error);
}
void TH1::SetCellError(Int_t binx, Int_t biny, Double_t error)
{
if (binx <0 || binx>fXaxis.GetNbins()+1) return;
if (biny <0 || biny>fYaxis.GetNbins()+1) return;
if (!fSumw2.fN) Sumw2();
Int_t bin = biny*(fXaxis.GetNbins()+2) + binx;
fSumw2.fArray[bin] = error*error;
}
void TH1::SetBinContent(Int_t, Double_t)
{
AbstractMethod("SetBinContent");
}
TH1 *TH1::ShowBackground(Int_t niter, Option_t *option)
{
return (TH1*)gROOT->ProcessLineFast(Form("TSpectrum::StaticBackground((TH1*)0x%x,%d,\"%s\")",this,niter,option));
}
Int_t TH1::ShowPeaks(Double_t sigma, Option_t *option, Double_t threshold)
{
return (Int_t)gROOT->ProcessLineFast(Form("TSpectrum::StaticSearch((TH1*)0x%x,%g,\"%s\",%g)",this,sigma,option,threshold));
}
TH1* TH1::TransformHisto(TVirtualFFT *fft, TH1* h_output, Option_t *option)
{
if (fft->GetNdim()>2){
printf("Only 1d and 2d\n");
return 0;
}
Int_t binx,biny;
TString opt = option;
opt.ToUpper();
Int_t *n = fft->GetN();
TH1 *hout=0;
if (h_output) hout = h_output;
else {
char name[10];
sprintf(name, "out_%s", opt.Data());
if (fft->GetNdim()==1)
hout = new TH1D(name, name,n[0], 0, n[0]);
else if (fft->GetNdim()==2)
hout = new TH2D(name, name, n[0], 0, n[0], n[1], 0, n[1]);
}
TString type=fft->GetType();
Int_t ind[2];
if (opt.Contains("RE")){
if (type.Contains("2C") || type.Contains("2HC")) {
Double_t re, im;
for (binx = 1; binx<=hout->GetNbinsX(); binx++) {
for (biny=1; biny<=hout->GetNbinsY(); biny++) {
ind[0] = binx-1; ind[1] = biny-1;
fft->GetPointComplex(ind, re, im);
hout->SetBinContent(binx, biny, re);
}
}
} else {
for (binx = 1; binx<=hout->GetNbinsX(); binx++) {
for (biny=1; biny<=hout->GetNbinsY(); biny++) {
ind[0] = binx-1; ind[1] = biny-1;
hout->SetBinContent(binx, biny, fft->GetPointReal(ind));
}
}
}
}
if (opt.Contains("IM")) {
if (type.Contains("2C") || type.Contains("2HC")) {
Double_t re, im;
for (binx = 1; binx<=hout->GetNbinsX(); binx++) {
for (biny=1; biny<=hout->GetNbinsY(); biny++) {
ind[0] = binx-1; ind[1] = biny-1;
fft->GetPointComplex(ind, re, im);
hout->SetBinContent(binx, biny, im);
}
}
} else {
printf("No complex numbers in the output");
return 0;
}
}
if (opt.Contains("MA")) {
if (type.Contains("2C") || type.Contains("2HC")) {
Double_t re, im;
for (binx = 1; binx<=hout->GetNbinsX(); binx++) {
for (biny=1; biny<=hout->GetNbinsY(); biny++) {
ind[0] = binx-1; ind[1] = biny-1;
fft->GetPointComplex(ind, re, im);
hout->SetBinContent(binx, biny, TMath::Sqrt(re*re + im*im));
}
}
} else {
for (binx = 1; binx<=hout->GetNbinsX(); binx++) {
for (biny=1; biny<=hout->GetNbinsY(); biny++) {
ind[0] = binx-1; ind[1] = biny-1;
hout->SetBinContent(binx, biny, TMath::Abs(fft->GetPointReal(ind)));
}
}
}
}
if (opt.Contains("PH")) {
if (type.Contains("2C") || type.Contains("2HC")){
Double_t re, im, ph;
for (binx = 1; binx<=hout->GetNbinsX(); binx++){
for (biny=1; biny<=hout->GetNbinsY(); biny++){
ind[0] = binx-1; ind[1] = biny-1;
fft->GetPointComplex(ind, re, im);
if (TMath::Abs(re) > 1e-13){
ph = TMath::ATan(im/re);
if (re<0 && im<0)
ph -= TMath::Pi();
if (re<0 && im>=0)
ph += TMath::Pi();
} else {
if (TMath::Abs(im) < 1e-13)
ph = 0;
else if (im>0)
ph = TMath::Pi()*0.5;
else
ph = -TMath::Pi()*0.5;
}
hout->SetBinContent(binx, biny, ph);
}
}
} else {
printf("Pure real output, no phase");
return 0;
}
}
return hout;
}
ClassImp(TH1C)
TH1C::TH1C(): TH1(), TArrayC()
{
fDimension = 1;
SetBinsLength(3);
if (fgDefaultSumw2) Sumw2();
}
TH1C::TH1C(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
: TH1(name,title,nbins,xlow,xup)
{
fDimension = 1;
TArrayC::Set(fNcells);
if (xlow >= xup) SetBuffer(fgBufferSize);
if (fgDefaultSumw2) Sumw2();
}
TH1C::TH1C(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
: TH1(name,title,nbins,xbins)
{
fDimension = 1;
TArrayC::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1C::TH1C(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
: TH1(name,title,nbins,xbins)
{
fDimension = 1;
TArrayC::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1C::~TH1C()
{
}
TH1C::TH1C(const TH1C &h1c) : TH1(), TArrayC()
{
((TH1C&)h1c).Copy(*this);
}
void TH1C::AddBinContent(Int_t bin)
{
if (fArray[bin] < 127) fArray[bin]++;
}
void TH1C::AddBinContent(Int_t bin, Double_t w)
{
Int_t newval = fArray[bin] + Int_t(w);
if (newval > -128 && newval < 128) {fArray[bin] = Char_t(newval); return;}
if (newval < -127) fArray[bin] = -127;
if (newval > 127) fArray[bin] = 127;
}
void TH1C::Copy(TObject &newth1) const
{
TH1::Copy(newth1);
TArrayC::Copy((TH1C&)newth1);
}
TH1 *TH1C::DrawCopy(Option_t *option) const
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
TH1C *newth1 = (TH1C*)Clone();
newth1->SetDirectory(0);
newth1->SetBit(kCanDelete);
newth1->AppendPad(opt.Data());
return newth1;
}
Double_t TH1C::GetBinContent(Int_t bin) const
{
if (fBuffer) ((TH1C*)this)->BufferEmpty();
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
if (!fArray) return 0;
return Double_t (fArray[bin]);
}
void TH1C::Reset(Option_t *option)
{
TH1::Reset(option);
TArrayC::Reset();
}
void TH1C::SetBinContent(Int_t bin, Double_t content)
{
if (bin < 0) return;
if (bin >= fNcells-1) {
if (!fXaxis.GetTimeDisplay() && !TestBit(kCanRebin)) {
if (bin == fNcells-1) fArray[bin] = Char_t (content);
return;
}
while (bin >= fNcells-1) LabelsInflate();
}
fArray[bin] = Char_t (content);
fEntries++;
}
void TH1C::SetBinsLength(Int_t n)
{
if (n < 0) n = fXaxis.GetNbins() + 2;
fNcells = n;
TArrayC::Set(n);
}
TH1C& TH1C::operator=(const TH1C &h1)
{
if (this != &h1) ((TH1C&)h1).Copy(*this);
return *this;
}
TH1C operator*(Double_t c1, const TH1C &h1)
{
TH1C hnew = h1;
hnew.Scale(c1);
hnew.SetDirectory(0);
return hnew;
}
TH1C operator+(const TH1C &h1, const TH1C &h2)
{
TH1C hnew = h1;
hnew.Add(&h2,1);
hnew.SetDirectory(0);
return hnew;
}
TH1C operator-(const TH1C &h1, const TH1C &h2)
{
TH1C hnew = h1;
hnew.Add(&h2,-1);
hnew.SetDirectory(0);
return hnew;
}
TH1C operator*(const TH1C &h1, const TH1C &h2)
{
TH1C hnew = h1;
hnew.Multiply(&h2);
hnew.SetDirectory(0);
return hnew;
}
TH1C operator/(const TH1C &h1, const TH1C &h2)
{
TH1C hnew = h1;
hnew.Divide(&h2);
hnew.SetDirectory(0);
return hnew;
}
ClassImp(TH1S)
TH1S::TH1S(): TH1(), TArrayS()
{
fDimension = 1;
SetBinsLength(3);
if (fgDefaultSumw2) Sumw2();
}
TH1S::TH1S(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
: TH1(name,title,nbins,xlow,xup)
{
fDimension = 1;
TArrayS::Set(fNcells);
if (xlow >= xup) SetBuffer(fgBufferSize);
if (fgDefaultSumw2) Sumw2();
}
TH1S::TH1S(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
: TH1(name,title,nbins,xbins)
{
fDimension = 1;
TArrayS::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1S::TH1S(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
: TH1(name,title,nbins,xbins)
{
fDimension = 1;
TArrayS::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1S::~TH1S()
{
}
TH1S::TH1S(const TH1S &h1s) : TH1(), TArrayS()
{
((TH1S&)h1s).Copy(*this);
}
void TH1S::AddBinContent(Int_t bin)
{
if (fArray[bin] < 32767) fArray[bin]++;
}
void TH1S::AddBinContent(Int_t bin, Double_t w)
{
Int_t newval = fArray[bin] + Int_t(w);
if (newval > -32768 && newval < 32768) {fArray[bin] = Short_t(newval); return;}
if (newval < -32767) fArray[bin] = -32767;
if (newval > 32767) fArray[bin] = 32767;
}
void TH1S::Copy(TObject &newth1) const
{
TH1::Copy(newth1);
TArrayS::Copy((TH1S&)newth1);
}
TH1 *TH1S::DrawCopy(Option_t *option) const
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
TH1S *newth1 = (TH1S*)Clone();
newth1->SetDirectory(0);
newth1->SetBit(kCanDelete);
newth1->AppendPad(opt.Data());
return newth1;
}
Double_t TH1S::GetBinContent(Int_t bin) const
{
if (fBuffer) ((TH1S*)this)->BufferEmpty();
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
if (!fArray) return 0;
return Double_t (fArray[bin]);
}
void TH1S::Reset(Option_t *option)
{
TH1::Reset(option);
TArrayS::Reset();
}
void TH1S::SetBinContent(Int_t bin, Double_t content)
{
if (bin < 0) return;
if (bin >= fNcells-1) {
if (!fXaxis.GetTimeDisplay() && !TestBit(kCanRebin)) {
if (bin == fNcells-1) fArray[bin] = Short_t (content);
return;
}
while (bin >= fNcells-1) LabelsInflate();
}
fArray[bin] = Short_t (content);
fEntries++;
}
void TH1S::SetBinsLength(Int_t n)
{
if (n < 0) n = fXaxis.GetNbins() + 2;
fNcells = n;
TArrayS::Set(n);
}
TH1S& TH1S::operator=(const TH1S &h1)
{
if (this != &h1) ((TH1S&)h1).Copy(*this);
return *this;
}
TH1S operator*(Double_t c1, const TH1S &h1)
{
TH1S hnew = h1;
hnew.Scale(c1);
hnew.SetDirectory(0);
return hnew;
}
TH1S operator+(const TH1S &h1, const TH1S &h2)
{
TH1S hnew = h1;
hnew.Add(&h2,1);
hnew.SetDirectory(0);
return hnew;
}
TH1S operator-(const TH1S &h1, const TH1S &h2)
{
TH1S hnew = h1;
hnew.Add(&h2,-1);
hnew.SetDirectory(0);
return hnew;
}
TH1S operator*(const TH1S &h1, const TH1S &h2)
{
TH1S hnew = h1;
hnew.Multiply(&h2);
hnew.SetDirectory(0);
return hnew;
}
TH1S operator/(const TH1S &h1, const TH1S &h2)
{
TH1S hnew = h1;
hnew.Divide(&h2);
hnew.SetDirectory(0);
return hnew;
}
ClassImp(TH1I)
TH1I::TH1I(): TH1(), TArrayI()
{
fDimension = 1;
SetBinsLength(3);
if (fgDefaultSumw2) Sumw2();
}
TH1I::TH1I(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
: TH1(name,title,nbins,xlow,xup)
{
fDimension = 1;
TArrayI::Set(fNcells);
if (xlow >= xup) SetBuffer(fgBufferSize);
if (fgDefaultSumw2) Sumw2();
}
TH1I::TH1I(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
: TH1(name,title,nbins,xbins)
{
fDimension = 1;
TArrayI::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1I::TH1I(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
: TH1(name,title,nbins,xbins)
{
fDimension = 1;
TArrayI::Set(fNcells);
}
TH1I::~TH1I()
{
}
TH1I::TH1I(const TH1I &h1i) : TH1(), TArrayI()
{
((TH1I&)h1i).Copy(*this);
}
void TH1I::AddBinContent(Int_t bin)
{
if (fArray[bin] < 2147483647) fArray[bin]++;
}
void TH1I::AddBinContent(Int_t bin, Double_t w)
{
Int_t newval = fArray[bin] + Int_t(w);
if (newval > -2147483647 && newval < 2147483647) {fArray[bin] = Int_t(newval); return;}
if (newval < -2147483647) fArray[bin] = -2147483647;
if (newval > 2147483647) fArray[bin] = 2147483647;
}
void TH1I::Copy(TObject &newth1) const
{
TH1::Copy(newth1);
TArrayI::Copy((TH1I&)newth1);
}
TH1 *TH1I::DrawCopy(Option_t *option) const
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
TH1I *newth1 = (TH1I*)Clone();
newth1->SetDirectory(0);
newth1->SetBit(kCanDelete);
newth1->AppendPad(opt.Data());
return newth1;
}
Double_t TH1I::GetBinContent(Int_t bin) const
{
if (fBuffer) ((TH1I*)this)->BufferEmpty();
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
if (!fArray) return 0;
return Double_t (fArray[bin]);
}
void TH1I::Reset(Option_t *option)
{
TH1::Reset(option);
TArrayI::Reset();
}
void TH1I::SetBinContent(Int_t bin, Double_t content)
{
if (bin < 0) return;
if (bin >= fNcells-1) {
if (!fXaxis.GetTimeDisplay() && !TestBit(kCanRebin)) {
if (bin == fNcells-1) fArray[bin] = Int_t (content);
return;
}
while (bin >= fNcells-1) LabelsInflate();
}
fArray[bin] = Int_t (content);
fEntries++;
}
void TH1I::SetBinsLength(Int_t n)
{
if (n < 0) n = fXaxis.GetNbins() + 2;
fNcells = n;
TArrayI::Set(n);
}
TH1I& TH1I::operator=(const TH1I &h1)
{
if (this != &h1) ((TH1I&)h1).Copy(*this);
return *this;
}
TH1I operator*(Double_t c1, const TH1I &h1)
{
TH1I hnew = h1;
hnew.Scale(c1);
hnew.SetDirectory(0);
return hnew;
}
TH1I operator+(const TH1I &h1, const TH1I &h2)
{
TH1I hnew = h1;
hnew.Add(&h2,1);
hnew.SetDirectory(0);
return hnew;
}
TH1I operator-(const TH1I &h1, const TH1I &h2)
{
TH1I hnew = h1;
hnew.Add(&h2,-1);
hnew.SetDirectory(0);
return hnew;
}
TH1I operator*(const TH1I &h1, const TH1I &h2)
{
TH1I hnew = h1;
hnew.Multiply(&h2);
hnew.SetDirectory(0);
return hnew;
}
TH1I operator/(const TH1I &h1, const TH1I &h2)
{
TH1I hnew = h1;
hnew.Divide(&h2);
hnew.SetDirectory(0);
return hnew;
}
ClassImp(TH1F)
TH1F::TH1F(): TH1(), TArrayF()
{
fDimension = 1;
SetBinsLength(3);
if (fgDefaultSumw2) Sumw2();
}
TH1F::TH1F(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
: TH1(name,title,nbins,xlow,xup)
{
fDimension = 1;
TArrayF::Set(fNcells);
if (xlow >= xup) SetBuffer(fgBufferSize);
if (fgDefaultSumw2) Sumw2();
}
TH1F::TH1F(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
: TH1(name,title,nbins,xbins)
{
fDimension = 1;
TArrayF::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1F::TH1F(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
: TH1(name,title,nbins,xbins)
{
fDimension = 1;
TArrayF::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1F::TH1F(const TVectorF &v)
: TH1("TVectorF","",v.GetNrows(),0,v.GetNrows())
{
TArrayF::Set(fNcells);
fDimension = 1;
Int_t ivlow = v.GetLwb();
for (Int_t i=0;i<fNcells-2;i++) {
SetBinContent(i+1,v(i+ivlow));
}
TArrayF::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1F::TH1F(const TH1F &h) : TH1(), TArrayF()
{
((TH1F&)h).Copy(*this);
}
TH1F::~TH1F()
{
}
void TH1F::Copy(TObject &newth1) const
{
TH1::Copy(newth1);
TArrayF::Copy((TH1F&)newth1);
}
TH1 *TH1F::DrawCopy(Option_t *option) const
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
TH1F *newth1 = (TH1F*)Clone();
newth1->SetDirectory(0);
newth1->SetBit(kCanDelete);
newth1->AppendPad(opt.Data());
return newth1;
}
Double_t TH1F::GetBinContent(Int_t bin) const
{
if (fBuffer) ((TH1F*)this)->BufferEmpty();
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
if (!fArray) return 0;
return Double_t (fArray[bin]);
}
void TH1F::Reset(Option_t *option)
{
TH1::Reset(option);
TArrayF::Reset();
}
void TH1F::SetBinContent(Int_t bin, Double_t content)
{
if (bin < 0) return;
if (bin >= fNcells-1) {
if (!fXaxis.GetTimeDisplay() && !TestBit(kCanRebin)) {
if (bin == fNcells-1) fArray[bin] = Float_t (content);
return;
}
while (bin >= fNcells-1) LabelsInflate();
}
fArray[bin] = Float_t (content);
fEntries++;
}
void TH1F::SetBinsLength(Int_t n)
{
if (n < 0) n = fXaxis.GetNbins() + 2;
fNcells = n;
TArrayF::Set(n);
}
TH1F& TH1F::operator=(const TH1F &h1)
{
if (this != &h1) ((TH1F&)h1).Copy(*this);
return *this;
}
TH1F operator*(Double_t c1, const TH1F &h1)
{
TH1F hnew = h1;
hnew.Scale(c1);
hnew.SetDirectory(0);
return hnew;
}
TH1F operator+(const TH1F &h1, const TH1F &h2)
{
TH1F hnew = h1;
hnew.Add(&h2,1);
hnew.SetDirectory(0);
return hnew;
}
TH1F operator-(const TH1F &h1, const TH1F &h2)
{
TH1F hnew = h1;
hnew.Add(&h2,-1);
hnew.SetDirectory(0);
return hnew;
}
TH1F operator*(const TH1F &h1, const TH1F &h2)
{
TH1F hnew = h1;
hnew.Multiply(&h2);
hnew.SetDirectory(0);
return hnew;
}
TH1F operator/(const TH1F &h1, const TH1F &h2)
{
TH1F hnew = h1;
hnew.Divide(&h2);
hnew.SetDirectory(0);
return hnew;
}
ClassImp(TH1D)
TH1D::TH1D(): TH1(), TArrayD()
{
fDimension = 1;
SetBinsLength(3);
if (fgDefaultSumw2) Sumw2();
}
TH1D::TH1D(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
: TH1(name,title,nbins,xlow,xup)
{
fDimension = 1;
TArrayD::Set(fNcells);
if (xlow >= xup) SetBuffer(fgBufferSize);
if (fgDefaultSumw2) Sumw2();
}
TH1D::TH1D(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
: TH1(name,title,nbins,xbins)
{
fDimension = 1;
TArrayD::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1D::TH1D(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
: TH1(name,title,nbins,xbins)
{
fDimension = 1;
TArrayD::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1D::TH1D(const TVectorD &v)
: TH1("TVectorD","",v.GetNrows(),0,v.GetNrows())
{
TArrayD::Set(fNcells);
fDimension = 1;
Int_t ivlow = v.GetLwb();
for (Int_t i=0;i<fNcells-2;i++) {
SetBinContent(i+1,v(i+ivlow));
}
TArrayD::Set(fNcells);
if (fgDefaultSumw2) Sumw2();
}
TH1D::~TH1D()
{
}
TH1D::TH1D(const TH1D &h1d) : TH1(), TArrayD()
{
((TH1D&)h1d).Copy(*this);
}
void TH1D::Copy(TObject &newth1) const
{
TH1::Copy(newth1);
TArrayD::Copy((TH1D&)newth1);
}
TH1 *TH1D::DrawCopy(Option_t *option) const
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
TH1D *newth1 = (TH1D*)Clone();
newth1->SetDirectory(0);
newth1->SetBit(kCanDelete);
newth1->AppendPad(opt.Data());
return newth1;
}
Double_t TH1D::GetBinContent(Int_t bin) const
{
if (fBuffer) ((TH1D*)this)->BufferEmpty();
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
if (!fArray) return 0;
return Double_t (fArray[bin]);
}
void TH1D::Reset(Option_t *option)
{
TH1::Reset(option);
TArrayD::Reset();
}
void TH1D::SetBinContent(Int_t bin, Double_t content)
{
if (bin < 0) return;
if (bin >= fNcells-1) {
if (!fXaxis.GetTimeDisplay() && !TestBit(kCanRebin)) {
if (bin == fNcells-1) fArray[bin] = content;
return;
}
while (bin >= fNcells-1) LabelsInflate();
}
fArray[bin] = content;
fEntries++;
}
void TH1D::SetBinsLength(Int_t n)
{
if (n < 0) n = fXaxis.GetNbins() + 2;
fNcells = n;
TArrayD::Set(n);
}
TH1D& TH1D::operator=(const TH1D &h1)
{
if (this != &h1) ((TH1D&)h1).Copy(*this);
return *this;
}
TH1D operator*(Double_t c1, const TH1D &h1)
{
TH1D hnew = h1;
hnew.Scale(c1);
hnew.SetDirectory(0);
return hnew;
}
TH1D operator+(const TH1D &h1, const TH1D &h2)
{
TH1D hnew = h1;
hnew.Add(&h2,1);
hnew.SetDirectory(0);
return hnew;
}
TH1D operator-(const TH1D &h1, const TH1D &h2)
{
TH1D hnew = h1;
hnew.Add(&h2,-1);
hnew.SetDirectory(0);
return hnew;
}
TH1D operator*(const TH1D &h1, const TH1D &h2)
{
TH1D hnew = h1;
hnew.Multiply(&h2);
hnew.SetDirectory(0);
return hnew;
}
TH1D operator/(const TH1D &h1, const TH1D &h2)
{
TH1D hnew = h1;
hnew.Divide(&h2);
hnew.SetDirectory(0);
return hnew;
}
TH1 *R__H(Int_t hid)
{
char hname[20];
if(hid >= 0) sprintf(hname,"h%d",hid);
else sprintf(hname,"h_%d",hid);
return (TH1*)gDirectory->Get(hname);
}
TH1 *R__H(const char * hname)
{
return (TH1*)gDirectory->Get(hname);
}
ROOT page - Class index - Class Hierarchy - 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.