//*CMZ : 2.23/08 31/10/99 18.32.25 by Rene Brun
//*CMZ : 2.23/07 25/10/99 08.20.53 by Rene Brun
//*CMZ : 2.23/04 13/10/99 17.58.35 by Rene Brun
//*CMZ : 2.23/03 27/09/99 14.25.41 by Rene Brun
//*CMZ : 2.23/02 07/09/99 14.53.38 by Rene Brun
//*CMZ : 2.23/01 31/08/99 16.59.36 by Rene Brun
//*CMZ : 2.23/00 30/07/99 20.04.37 by Fons Rademakers
//*-- Author : Rene Brun 26/12/94
//*KEEP,CopyRight,T=C.
/*************************************************************************
* Copyright(c) 1995-1999, The ROOT System, All rights reserved. *
* Authors: Rene Brun and Fons Rademakers. *
* *
* For the licensing terms see $ROOTSYS/AA_LICENSE. *
* For the list of contributors see $ROOTSYS/AA_CREDITS. *
*************************************************************************/
//*KEND.
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <ctype.h>
#include <fstream.h>
#include <iostream.h>
//*KEEP,TROOT.
#include "TROOT.h"
//*KEEP,TH1.
#include "TH1.h"
//*KEEP,TH2.
#include "TH2.h"
//*KEEP,TF2.
#include "TF2.h"
//*KEEP,TF3.
#include "TF3.h"
//*KEEP,TVirtualPad.
#include "TVirtualPad.h"
//*KEEP,Foption.
#include "Foption.h"
//*KEEP,TMath.
#include "TMath.h"
//*KEEP,TRandom,T=C++.
#include "TRandom.h"
//*KEEP,TVirtualFitter,T=C++.
#include "TVirtualFitter.h"
//*KEEP,TProfile.
#include "TProfile.h"
//*KEEP,TStyle.
#include "TStyle.h"
//*KEND.
//______________________________________________________________________________
//*-*-*-*-*-*-*-*-*-*-*The H I S T O G R A M Classes*-*-*-*-*-*-*-*-*-*-*-*
//*-* ===============================
//*-*
//*-* ROOT supports the following histogram types:
//*-*
//*-* 1-D histograms:
//*-* TH1C : histograms with one byte per channel. Maximum bin content = 255
//*-* TH1S : histograms with one short per channel. Maximum bin content = 65535
//*-* TH1F : histograms with one float per channel. Maximum precision 7 digits
//*-* TH1D : histograms with one double per channel. Maximum precision 14 digits
//*-*
//*-* 2-D histograms:
//*-* TH2C : histograms with one byte per channel. Maximum bin content = 255
//*-* TH2S : histograms with one short per channel. Maximum bin content = 65535
//*-* TH2F : histograms with one float per channel. Maximum precision 7 digits
//*-* TH2D : histograms with one double per channel. Maximum precision 14 digits
//*-*
//*-* 3-D histograms:
//*-* TH3C : histograms with one byte per channel. Maximum bin content = 255
//*-* TH3S : histograms with one short per channel. Maximum bin content = 65535
//*-* TH3F : histograms with one float per channel. Maximum precision 7 digits
//*-* TH3D : histograms with one double per channel. Maximum precision 14 digits
//*-*
//*-* Profile histograms: See class TProfile
//*-*
//*-*- All histogram classes are derived from the base class TH1
//*-* The TH*C classes also inherit from the array class TArrayC.
//*-* The TH*S classes also inherit from the array class TArrayS.
//*-* The TH*F classes also inherit from the array class TArrayF.
//*-* The TH*D classes also inherit from the array class TArrayD.
//*-*
//*-* Convention for numbering bins
//*-* =============================
//*-* For all histogram types: nbins, xlow, xup
//*-* bin = 0; underflow bin
//*-* bin = 1; first bin with low-edge xlow INCLUDED
//*-* bin = nbins; last bin with upper-edge xup EXCLUDED
//*-* bin = nbins+1; overflow bin
//*-*
//
/*
*/
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Foption_t Foption;
TAxis *xaxis=0;
TAxis *yaxis=0;
TAxis *zaxis=0;
TF1 *gF1=0;
TVirtualFitter *hFitter=0;
static Int_t xfirst,xlast,yfirst,ylast,zfirst,zlast;
static Float_t xmin, xmax, ymin, ymax, zmin, zmax, binwidx, binwidy, binwidz;
const Int_t kCloseDirectory = BIT(7);
const Int_t kNoStats = BIT(9);
const Int_t kUserContour = BIT(10);
const Int_t kAxisRange = BIT(11);
extern void H1InitGaus();
extern void H1InitExpo();
extern void H1InitPolynom();
extern void H1FitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
extern void H1FitLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
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()
{
//*-*-*-*-*-*-*-*-*-*-*Histogram default constructor*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* =============================
fDirectory = 0;
fFunctions = new TList(this);
fIntegral = 0;
fPainter = 0;
}
//______________________________________________________________________________
TH1::~TH1()
{
//*-*-*-*-*-*-*-*-*-*-*Histogram default destructor*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* ============================
if (!TestBit(kNotDeleted)) return;
if (fIntegral) {delete [] fIntegral; fIntegral = 0;}
if (fFunctions) { fFunctions->Delete(); delete fFunctions; }
if (fDirectory) {
if (!fDirectory->TestBit(kCloseDirectory)) fDirectory->GetList()->Remove(this);
}
delete fPainter;
fDirectory = 0;
fFunctions = 0;
}
//______________________________________________________________________________
TH1::TH1(const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
:TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*-*-*Normal constructor for fix bin size histograms*-*-*-*-*-*-*
//*-* ==============================================
//
// Creates the main histogram structure:
// name : name of histogram (avoid blanks)
// title : histogram title
// nbins : number of bins
// xlow : low edge of first bin
// xup : upper edge of last bin (not included in last bin)
//
// When an histogram is created, it is automatically added to the list
// of special objects in the current directory.
// To find the pointer to this histogram in the current directory
// by its name, do:
// TH1F *h1 = (TH1F*)gDirectory->FindObject(name);
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Build();
if (nbins <= 0) nbins = 1;
fXaxis.Set(nbins,xlow,xup);
fNcells = fXaxis.GetNbins()+2;
}
//______________________________________________________________________________
TH1::TH1(const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
:TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*Normal constructor for variable bin size histograms*-*-*-*-*-*-*
//*-* ===================================================
//
// Creates the main histogram structure:
// name : name of histogram (avoid blanks)
// title : histogram title
// nbins : number of bins
// xbins : array of low-edges for each bin
// This is an array of size nbins+1
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Build();
if (nbins <= 0) nbins = 1;
if (xbins) fXaxis.Set(nbins,xbins);
else fXaxis.Set(nbins,0,1);
fNcells = fXaxis.GetNbins()+2;
}
//______________________________________________________________________________
void TH1::Browse(TBrowser *)
{
Draw();
gPad->Update();
}
//______________________________________________________________________________
void TH1::Build()
{
//*-*-*-*-*-*-*-*-*-*Creates histogram basic data structure*-*-*-*-*-*-*-*-*-*
//*-* ======================================
fPainter = 0;
fIntegral = 0;
fEntries = 0;
fNormFactor = 0;
fTsumw = fTsumw2=fTsumwx=fTsumwx2=0;
fMaximum = -1111;
fMinimum = -1111;
fXaxis.SetName("xaxis");
fYaxis.SetName("yaxis");
fZaxis.SetName("zaxis");
fYaxis.Set(1,0.,1.);
fZaxis.Set(1,0.,1.);
UseCurrentStyle();
if (gDirectory) {
TH1 *hold = (TH1*)gDirectory->GetList()->FindObject(GetName());
if (hold) {
Warning("Build","Replacing existing histogram: %s",GetName());
gDirectory->GetList()->Remove(hold);
// delete hold;
}
gDirectory->Append(this);
}
fDirectory = gDirectory;
fFunctions = new TList(this);
}
//______________________________________________________________________________
void TH1::Add(TH1 *h1, Float_t c1)
{
// Performs the operation: this = this + c1*h1
// if errors are defined (see TH1::Sumw2), errors are also recalculated.
// Note that if h1 has Sumw2 set, Sumw2 is automatically called for this
// if not already set.
if (!h1) {
Error("Add","Attempt to add a non-existing histogram");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
//*-*- Check histogram compatibility
if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()) {
Error("Add","Attempt to add histograms with different number of bins");
return;
}
//*-*- Issue a Warning if histogram limits are different
if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
Warning("Add","Attempt to add histograms with different axis limits");
}
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
//*-* Create Sumw2 if h1 has Sumw2 set
if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
//*-*- Add statistics
Float_t ac1 = TMath::Abs(c1);
fEntries += ac1*h1->GetEntries();
fTsumw += ac1*h1->fTsumw;
fTsumw2 += ac1*h1->fTsumw2;
fTsumwx += ac1*h1->fTsumwx;
fTsumwx2 += ac1*h1->fTsumwx2;
//*-*- Loop on bins (including underflows/overflows)
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);
AddBinContent(bin,cu);
if (fSumw2.fN) {
Double_t error1 = h1->GetBinError(bin);
fSumw2.fArray[bin] += c1*c1*error1*error1;
}
}
}
}
}
//______________________________________________________________________________
void TH1::Add(TH1 *h1, TH1 *h2, Float_t c1, Float_t c2)
{
//*-*-*-*-*Replace contents of this histogram by the addition of h1 and h2*-*-*
//*-* ===============================================================
//
// this = c1*h1 + c2*h2
// if errors are defined (see TH1::Sumw2), errors are also recalculated
// Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for this
// if not already set.
//
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();
//*-*- Check histogram compatibility
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;
}
//*-*- Issue a Warning if histogram limits are different
if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
Warning("Add","Attempt to add histograms with different axis limits");
}
if (GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin() ||
GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ||
GetYaxis()->GetXmin() != h2->GetYaxis()->GetXmin() ||
GetYaxis()->GetXmax() != h2->GetYaxis()->GetXmax() ||
GetZaxis()->GetXmin() != h2->GetZaxis()->GetXmin() ||
GetZaxis()->GetXmax() != h2->GetZaxis()->GetXmax()) {
Warning("Add","Attempt to add histograms with different axis limits");
}
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
if (fDimension < 3) nbinsz = -1;
//*-* Create Sumw2 if h1 or h2 have Sumw2 set
if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
//*-*- Add statistics
Float_t ac1 = TMath::Abs(c1);
Float_t ac2 = TMath::Abs(c2);
fEntries = ac1*h1->GetEntries() + ac2*h2->GetEntries();
fTsumw = ac1*h1->fTsumw + ac2*h2->fTsumw;
fTsumw2 = ac1*h1->fTsumw2 + ac2*h2->fTsumw2;
fTsumwx = ac1*h1->fTsumwx + ac2*h2->fTsumwx;
fTsumwx2 = ac1*h1->fTsumwx2 + ac2*h2->fTsumwx2;
//*-*- Loop on bins (including underflows/overflows)
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;
}
}
}
}
}
//______________________________________________________________________________
void TH1::AddBinContent(Int_t)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by 1*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* ==========================
AbstractMethod("AddBinContent");
}
//______________________________________________________________________________
void TH1::AddBinContent(Int_t, Stat_t)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by a weight w*-*-*-*-*-*-*-*-*-*-*
//*-* ===================================
AbstractMethod("AddBinContent");
}
//______________________________________________________________________________
Double_t TH1::ComputeIntegral()
{
// Compute integral (cumulative sum of bins)
// The result stored in fIntegral is used by the GetRandom functions.
// This function is automatically called by GetRandom when the fIntegral
// array does not exist or when the number of entries in the histogram
// has changed since the previous call to GetRandom.
// The resulting integral is normalized to 1
Int_t bin, binx, biny, binz, ibin;
// delete previously computed integral (if any)
if (fIntegral) delete [] fIntegral;
//*-*- Allocate space to store the integral and compute integral
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);
}
}
}
//*-*- Normalize integral to 1
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];
}
//______________________________________________________________________________
void TH1::Copy(TObject &obj)
{
//*-*-*-*-*-*-*Copy this histogram structure to newth1*-*-*-*-*-*-*-*-*-*-*-*
//*-* =======================================
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;
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);
fContour.Copy(((TH1&)obj).fContour);
fSumw2.Copy(((TH1&)obj).fSumw2);
// fFunctions->Copy(((TH1&)obj).fFunctions);
gDirectory->Append(&obj);
// ((TH1&)obj).AppendDirectory();
((TH1&)obj).fDirectory = gDirectory;
}
//______________________________________________________________________________
Int_t TH1::DistancetoPrimitive(Int_t px, Int_t py)
{
//*-*-*-*-*-*-*-*-*-*-*Compute distance from point px,py to a line*-*-*-*-*-*
//*-* ===========================================
//*-* Compute the closest distance of approach from point px,py to elements
//*-* of an histogram.
//*-* The distance is computed in pixels units.
//*-*
//*-* Algorithm:
//*-* Currently, this simple model computes the distance from the mouse
//*-* to the histogram contour only.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
if (!fPainter) return 9999;
return fPainter->DistancetoPrimitive(px,py);
}
//______________________________________________________________________________
void TH1::Divide(TH1 *h1)
{
//*-*-*-*-*-*-*-*-*-*-*Divide this histogram by h1*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* ===========================
//
// this = this/h1
// if errors are defined (see TH1::Sumw2), errors are also recalculated.
// Note that if h1 has Sumw2 set, Sumw2 is automatically called for this
// if not already set.
//
if (!h1) {
Error("Divide","Attempt to divide a non-existing histogram");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
//*-*- Check histogram compatibility
if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()) {
Error("Divide","Attempt to divide histograms with different number of bins");
return;
}
//*-*- Issue a Warning if histogram limits are different
if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
GetZaxis()->GetXmax() != h1->GetZaxis()->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;
//*-* Create Sumw2 if h1 has Sumw2 set
if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
//*-*- Reset statistics
fEntries = fTsumw = fTsumw2 = fTsumwx = fTsumwx2 = 0;
//*-*- Loop on bins (including underflows/overflows)
Int_t bin, binx, biny, binz;
Double_t c0,c1,w,z,x;
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);
z = TMath::Abs(w);
x = fXaxis.GetBinCenter(binx);
fEntries++;
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
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);
}
}
}
}
}
//______________________________________________________________________________
void TH1::Divide(TH1 *h1, TH1 *h2, Float_t c1, Float_t c2, Option_t *option)
{
//*-*-*-*-*Replace contents of this histogram by the division of h1 by h2*-*-*
//*-* ==============================================================
//
// this = c1*h1/(c2*h2)
//
// if errors are defined (see TH1::Sumw2), errors are also recalculated
// Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for this
// if not already set.
//
TString opt = option;
opt.ToLower();
Bool_t binomial = kFALSE;
if (opt.Contains("b")) binomial = kTRUE;
if (!h1 || !h2) {
Error("Divide","Attempt to divide a non-existing histogram");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
//*-*- Check histogram compatibility
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;
}
//*-*- Issue a Warning if histogram limits are different
if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
Warning("Divide","Attempt to divide histograms with different axis limits");
}
if (GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin() ||
GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ||
GetYaxis()->GetXmin() != h2->GetYaxis()->GetXmin() ||
GetYaxis()->GetXmax() != h2->GetYaxis()->GetXmax() ||
GetZaxis()->GetXmin() != h2->GetZaxis()->GetXmin() ||
GetZaxis()->GetXmax() != h2->GetZaxis()->GetXmax()) {
Warning("Divide","Attempt to divide histograms with different axis limits");
}
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
//*-* Create Sumw2 if h1 or h2 have Sumw2 set
if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
//*-*- Reset statistics
fEntries = fTsumw = fTsumw2 = fTsumwx = fTsumwx2 = 0;
//*-*- Loop on bins (including underflows/overflows)
Int_t bin, binx, biny, binz;
Double_t b1,b2,w,z,x,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);
z = TMath::Abs(w);
x = fXaxis.GetBinCenter(binx);
fEntries++;
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
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) {
fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/(c2*b2));
} else {
fSumw2.fArray[bin] = d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1)/(b22*b22);
}
}
}
}
}
}
//______________________________________________________________________________
void TH1::Draw(Option_t *option)
{
//*-*-*-*-*-*-*-*-*-*-*Draw this histogram with options*-*-*-*-*-*-*-*-*-*-*-*
//*-* ================================
//*-*
//*-* This histogram is added in the list of objects to be drawn in the
//*-* current pad.
//*-* It is important to note that only a reference to this histogram is added.
//*-* This means that if the histogram contents change between pad redraws,
//*-* the current histogram contents will be shown.
//*-* Use function DrawCopy to force a copy of the histogram at the time
//*-* of the Draw operation in the Pad list.
//*-*
//*-* The following options are supported on all types:
//*-* "SAME" : Superimpose on previous picture in the same pad
//*-* "CYL" : Use Cylindrical coordinates
//*-* "POL" : Use Polar coordinates
//*-* "SPH" : Use Spherical coordinates
//*-* "PSR" : Use PseudoRapidity/Phi coordinates
//*-* "LEGO" : Draw a lego plot with hidden line removal
//*-* "LEGO1" : Draw a lego plot with hidden surface removal
//*-* "LEGO2" : Draw a lego plot using colors to show the cell contents
//*-* "SURF" : Draw a surface plot with hidden line removal
//*-* "SURF1" : Draw a surface plot with hidden surface removal
//*-* "SURF2" : Draw a surface plot using colors to show the cell contents
//*-* "SURF3" : same as SURF with in addition a contour view drawn on the top
//*-* "SURF4" : Draw a surface using Gouraud shading
//*-*
//*-* The following options are supported for 1-D types:
//*-* "AH" : Draw histogram, but not the axis labels and tick marks
//*-* "B" : Bar chart option
//*-* "C" : Draw a smooth Curve througth the histogram bins
//*-* "E" : Draw error bars
//*-* "E0" : Draw error bars including bins with o contents
//*-* "E1" : Draw error bars with perpendicular lines at the edges
//*-* "E2" : Draw error bars with rectangles
//*-* "E3" : Draw a fill area througth the end points of the vertical error bars
//*-* "E4" : Draw a smoothed filled area through the end points of the error bars
//*-* "L" : Draw a line througth the bin contents
//*-* "P" : Draw current marker at each bin
//*-* "*H" : Draw histogram with a * at each bin
//*-*
//*-*
//*-* The following options are supported for 2-D types:
//*-* "ARR" : arrow mode. Shows gradient between adjacent cells
//*-* "BOX" : a box is drawn for each cell with surface proportional to contents
//*-* "COL" : a box is drawn for each cell with a color scale varying with contents
//*-* "COLZ" : same as "COL". In addition the color mapping is also drawn
//*-* "CONT" : Draw a contour plot (same as CONT0)
//*-* "CONT0" : Draw a contour plot using colors to distinguish contours
//*-* "CONT1" : Draw a contour plot using line styles to distinguish contours
//*-* "CONT2" : Draw a contour plot using the same line style for all contours
//*-* "CONT3" : Draw a contour plot using fill area colors
//*-* "FB" : With LEGO or SURFACE, suppress the Front-Box
//*-* "BB" : With LEGO or SURFACE, suppress the Back-Box
//*-* "SCAT" : Draw a scatter-plot (default)
//*-*
//*-* Note that most options above can be concatenated, example:
//*-* h-Draw("e1same");
//*-* Options are case insensitive
//*-*
//*-* When using the options "BOX", "COL" or "COLZ", the color palette used
//*-* is the one defined in the current style (see TStyle::SetPalette)
//*-*
//*-* When using the "CONT" or "SURF" or "LEGO" options, the number
//*-* of contour levels can be changed via TH1::SetContour.
//*-* (default is 20 equidistant levels)
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) {
//the following statement is necessary in case one attempts to draw
//a temporary histogram already in the current pad
if (TestBit(kCanDelete)) gPad->GetListOfPrimitives()->Remove(this);
gPad->Clear();
}
AppendPad(option);
}
//______________________________________________________________________________
TH1 *TH1::DrawCopy(Option_t *)
{
//*-*-*-*-*-*-*Copy this histogram and Draw in the current pad*-*-*-*-*-*-*-*
//*-* ===============================================
//*-*
//*-* Once the histogram is drawn into the pad, any further modification
//*-* using graphics input will be made on the copy of the histogram,
//*-* and not to the original object.
//*-*
//*-* See Draw for the list of options
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
AbstractMethod("DrawCopy");
return 0;
}
//______________________________________________________________________________
void TH1::DrawPanel()
{
//*-*-*-*-*-*-*Display a panel with all histogram drawing options*-*-*-*-*-*
//*-* ==================================================
//*-*
//*-* See class TDrawPanelHist for example
if (fPainter) fPainter->DrawPanel();
}
//______________________________________________________________________________
void TH1::Eval(TF1 *f1, Option_t *option)
{
//*-*-*-*-*Evaluate function f1 at the center of bins of this histogram-*-*-*-*
//*-* ============================================================
//*-*
//*-* If option "R" is specified, the function is evaluated only
//*-* for the bins included in the function range.
//*-* If option "A" is specified, the value of the function is added to the
//*-* existing bin contents
//*-* If option "S" is specified, the value of the function is used to
//*-* generate an integer value, distributed according to the Poisson
//*-* distribution, with f1 as the mean.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Double_t x[3];
Int_t stat,add,bin,binx,biny,binz,nbinsx, nbinsy, nbinsz;
if (!f1) return;
Double_t fu;
TString opt = option;
opt.ToLower();
if (opt.Contains("a")) add = 1;
else add = 0;
if (opt.Contains("s")) stat = 1;
else stat = 0;
nbinsx = fXaxis.GetNbins();
nbinsy = fYaxis.GetNbins();
nbinsz = fZaxis.GetNbins();
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);
fu = f1->Eval(x[0],x[1],x[2]);
if (stat) fu = gRandom->Poisson(fu);
if (add) AddBinContent(bin,fu);
else SetBinContent(bin,fu);
}
}
}
}
//______________________________________________________________________________
void TH1::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
//*-*-*-*-*-*-*-*-*-*-*Execute action corresponding to one event*-*-*-*
//*-* =========================================
//*-* This member function is called when a histogram is clicked with the locator
//*-*
//*-* If Left button clicked on the bin top value, then the content of this bin
//*-* is modified according to the new position of the mouse when it is released.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
if (fPainter) fPainter->ExecuteEvent(event, px, py);
}
//______________________________________________________________________________
Int_t TH1::Fill(Axis_t x)
{
//*-*-*-*-*-*-*-*-*-*Increment bin with abscissa X by 1*-*-*-*-*-*-*-*-*-*-*
//*-* ==================================
//*-*
//*-* if x is less than the low-edge of the first bin, the Underflow bin is incremented
//*-* if x is greater than the upper edge of last bin, the Overflow bin is incremented
//*-*
//*-* If the storage of the sum of squares of weights has been triggered,
//*-* via the function Sumw2, then the sum of the squares of weights is incremented
//*-* by 1 in the bin corresponding to x.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Int_t bin;
fEntries++;
bin =fXaxis.FindBin(x);
AddBinContent(bin);
if (fSumw2.fN) ++fSumw2.fArray[bin];
if (bin == 0 || bin > fXaxis.GetNbins()) return -1;
++fTsumw;
++fTsumw2;
fTsumwx += x;
fTsumwx2 += x*x;
return bin;
}
//______________________________________________________________________________
Int_t TH1::Fill(Axis_t x, Stat_t w)
{
//*-*-*-*-*-*-*-*Increment bin with abscissa X with a weight w*-*-*-*-*-*-*-*
//*-* =============================================
//*-*
//*-* if x is less than the low-edge of the first bin, the Underflow bin is incremented
//*-* if x is greater than the upper edge of last bin, the Overflow bin is incremented
//*-*
//*-* If the storage of the sum of squares of weights has been triggered,
//*-* via the function Sumw2, then the sum of the squares of weights is incremented
//*-* by w^2 in the bin corresponding to x.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//authorize calling this function for 2-d histograms. assume weight=1
if (fDimension == 2) {
return Fill12(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()) return -1;
Stat_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, Axis_t *x, Double_t *w, Int_t stride)
{
//*-*-*-*-*-*-*-*Fill this histogram with an array x and weights w*-*-*-*-*
//*-* =================================================
//*-*
//*-* ntimes: number of entries in arrays x and w (array size must be ntimes*stride)
//*-* x: array of values to be histogrammed
//*-* w: array of weighs
//*-* stride: step size through arrays x and w
//*-*
//*-* If the storage of the sum of squares of weights has been triggered,
//*-* via the function Sumw2, then the sum of the squares of weights is incremented
//*-* by w[i]^2 in the bin corresponding to x[i].
//*-* if w is NULL each entry is assumed a weight=1
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
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) continue;
Stat_t z= (ww > 0 ? ww : -ww);
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x[i];
fTsumwx2 += z*x[i]*x[i];
}
}
//______________________________________________________________________________
Int_t TH1::Fill12(Axis_t x,Axis_t y)
{
//*-*-*-*-*-*-*-*-*-*-*Increment cell defined by x,y by 1*-*-*-*-*-*-*-*-*-*
//*-* ==================================
//*-*
//*-* if x or/and y is less than the low-edge of the corresponding axis first bin,
//*-* the Underflow cell is incremented.
//*-* if x or/and y is greater than the upper edge of corresponding axis last bin,
//*-* the Overflow cell is incremented.
//*-*
//*-* If the storage of the sum of squares of weights has been triggered,
//*-* via the function Sumw2, then the sum of the squares of weights is incremented
//*-* by 1in the cell corresponding to x,y.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Int_t binx, biny, bin;
fEntries++;
binx = fXaxis.FindBin(x);
biny = fYaxis.FindBin(y);
bin = biny*(fXaxis.GetNbins()+2) + binx;
AddBinContent(bin);
if (fSumw2.fN) ++fSumw2.fArray[bin];
if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
++fTsumw;
++fTsumw2;
fTsumwx += x;
fTsumwx2 += x*x;
return bin;
}
//______________________________________________________________________________
Int_t TH1::Fill(Axis_t x, Axis_t y, Stat_t w)
{
//*-*-*-*-*-*-*-*-*-*-*Increment cell defined by x,y by a weight w*-*-*-*-*-*
//*-* ===========================================
//*-*
//*-* if x or/and y is less than the low-edge of the corresponding axis first bin,
//*-* the Underflow cell is incremented.
//*-* if x or/and y is greater than the upper edge of corresponding axis last bin,
//*-* the Overflow cell is incremented.
//*-*
//*-* If the storage of the sum of squares of weights has been triggered,
//*-* via the function Sumw2, then the sum of the squares of weights is incremented
//*-* by w^2 in the cell corresponding to x,y.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//authorize calling this function for 3-d histograms. assume weight=1
if (fDimension == 3) {
return Fill(x, y, w, 1);
}
Int_t binx, biny, bin;
fEntries++;
binx = fXaxis.FindBin(x);
biny = fYaxis.FindBin(y);
bin = biny*(fXaxis.GetNbins()+2) + binx;
AddBinContent(bin,w);
if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
Stat_t z= (w > 0 ? w : -w);
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
return bin;
}
//______________________________________________________________________________
Int_t TH1::Fill(Axis_t x, Axis_t y, Axis_t z, Stat_t w)
{
//*-*-*-*-*-*-*-*-*-*-*Increment cell defined by x,y,z by a weight w*-*-*-*-*
//*-* =============================================
//*-*
//*-* If the storage of the sum of squares of weights has been triggered,
//*-* via the function Sumw2, then the sum of the squares of weights is incremented
//*-* by w^2 in the cell corresponding to x,y,z.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Int_t binx, biny, binz, bin;
fEntries++;
binx = fXaxis.FindBin(x);
biny = fYaxis.FindBin(y);
binz = fZaxis.FindBin(z);
bin = binx + (fXaxis.GetNbins()+2)*(biny + (fYaxis.GetNbins()+2)*binz);
AddBinContent(bin,w);
if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
if (binz == 0 || binz > fZaxis.GetNbins()) return -1;
return bin;
}
//______________________________________________________________________________
void TH1::FillN(Int_t ntimes, Axis_t *x, Axis_t *y, Double_t *w, Int_t stride)
{
//*-*-*-*-*-*-*Fill a 2-D histogram with an array of values and weights*-*-*-*
//*-* ========================================================
//*-*
//*-* ntimes: number of entries in arrays x and w (array size must be ntimes*stride)
//*-* x: array of x values to be histogrammed
//*-* y: array of y values to be histogrammed
//*-* w: array of weights
//*-* stride: step size through arrays x, y and w
//*-*
//*-* If the storage of the sum of squares of weights has been triggered,
//*-* via the function Sumw2, then the sum of the squares of weights is incremented
//*-* by w[i]^2 in the cell corresponding to x[i],y[i].
//*-* if w is NULL each entry is assumed a weight=1
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Int_t binx, biny, bin, i;
fEntries += ntimes;
Double_t ww = 1;
ntimes *= stride;
for (i=0;i<ntimes;i+=stride) {
binx = fXaxis.FindBin(x[i]);
biny = fYaxis.FindBin(y[i]);
bin = biny*(fXaxis.GetNbins()+2) + binx;
if (w) ww = w[i];
AddBinContent(bin,ww);
if (fSumw2.fN) fSumw2.fArray[bin] += ww*ww;
if (binx == 0 || binx > fXaxis.GetNbins()) continue;
if (biny == 0 || biny > fYaxis.GetNbins()) continue;
Stat_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)
{
//*-*-*-*-*-*-*Fill histogram following distribution in function fname*-*-*-*
//*-* =======================================================
//*-*
//*-* The distribution contained in the function fname (TF1) is integrated
//*-* over the channel contents.
//*-* It is normalized to 1.
//*-* Getting one random number implies:
//*-* - Generating a random number between 0 and 1 (say r1)
//*-* - Look in which bin in the normalized integral r1 corresponds to
//*-* - Fill histogram channel
//*-* ntimes random numbers are generated
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**-*-*-*-*-*-*-*
Int_t bin, binx, biny, binz, ibin, loop;
Double_t r1, x, y,z, xv[3];
//*-*- Search for fname in the list of ROOT defined functions
TF1 *f1 = (TF1*)gROOT->GetFunction(fname);
if (!f1) { Error("FillRandom", "Unknown function: %s",fname); return; }
//*-*- Allocate temporary space to store the integral and compute integral
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
Int_t nxy = nbinsx*nbinsy;
Int_t nbins = nxy*nbinsz;
Double_t *integral = new Double_t[nbins+1];
ibin = 0;
integral[ibin] = 0;
for (binz=1;binz<=nbinsz;binz++) {
xv[2] = fZaxis.GetBinCenter(binz);
for (biny=1;biny<=nbinsy;biny++) {
xv[1] = fYaxis.GetBinCenter(biny);
for (binx=1;binx<=nbinsx;binx++) {
xv[0] = fXaxis.GetBinCenter(binx);
ibin++;
integral[ibin] = integral[ibin-1] + f1->Eval(xv[0],xv[1],xv[2]);
}
}
}
//*-*- Normalize integral to 1
if (integral[nbins] == 0 ) {
Error("FillRandom", "Integral = zero"); return;
}
for (bin=1;bin<=nbins;bin++) integral[bin] /= integral[nbins];
//*-*--------------Start main loop ntimes
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
Bool_t profile = kFALSE;
if (IsA() == TProfile::Class()) profile = kTRUE;
for (loop=0;loop<ntimes;loop++) {
r1 = gRandom->Rndm(loop);
ibin = TMath::BinarySearch(nbins,&integral[0],r1);
binz = ibin/nxy;
biny = (ibin - nxy*binz)/nbinsx;
binx = 1 + ibin - nbinsx*(biny + nbinsy*binz);
if (nbinsz) binz++;
if (nbinsy) biny++;
x = GetBinCenter(binx);
y = GetBinCenter(biny);
z = GetBinCenter(binz);
if (fDimension == 1) {
if (profile) Fill(x,y,1.);
else Fill(x, 1.);
}
else if (fDimension == 2) Fill(x,y, 1.);
else if (fDimension == 3) Fill(x,y,z, 1.);
}
delete [] integral;
}
//______________________________________________________________________________
void TH1::FillRandom(TH1 *h, Int_t ntimes)
{
//*-*-*-*-*-*-*Fill histogram following distribution in histogram h*-*-*-*
//*-* ====================================================
//*-*
//*-* The distribution contained in the histogram h (TH1) is integrated
//*-* over the channel contents.
//*-* It is normalized to 1.
//*-* Getting one random number implies:
//*-* - Generating a random number between 0 and 1 (say r1)
//*-* - Look in which bin in the normalized integral r1 corresponds to
//*-* - Fill histogram channel
//*-* ntimes random numbers are generated
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**-*-*-*-*-*-*-*
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;
Axis_t x,y,z;
if (fDimension == 1) {
for (loop=0;loop<ntimes;loop++) {
x = h->GetRandom();
Fill(x);
}
} else if (fDimension == 2) {
for (loop=0;loop<ntimes;loop++) {
h->GetRandom2(x,y);
Fill(x,y,1.);
}
} else {
for (loop=0;loop<ntimes;loop++) {
h->GetRandom3(x,y,z);
Fill(x,y,z,1.);
}
}
}
//______________________________________________________________________________
Int_t TH1::FindBin(Axis_t x, Axis_t y, Axis_t z)
{
//*-*-*-*-*-*Return Global bin number corresponding to x,y,z*-*-*-*-*-*-*
//*-* ===============================================
//*-*
//*-* 2-D and 3-D histograms are represented with a one dimensional
//*-* structure.
//*-* This has the advantage that all existing functions, such as
//*-* GetBinContent, GetBinError, GetBinFunction work for all dimensions.
//*-* See also TH1::GetBin
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
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;
}
//______________________________________________________________________________
void TH1::Fit(const Text_t *fname ,Option_t *option ,Option_t *goption, Float_t xxmin, Float_t xxmax)
{
//*-*-*-*-*-*-*-*-*-*-*Fit histogram with function fname*-*-*-*-*-*-*-*-*-*-*
//*-* =================================
//*-*
//*-* fname is the name of an already predefined function created by TF1 or TF2
//*-* Predefined functions such as Gaus, Expo and Poln are automatically
//*-* created by ROOT.
//*-*
//*-* The list of fit options is given in parameter option.
//*-* option = "W" Set all errors to 1
//*-* = "I" Use integral of function in bin instead of value at bin center
//*-* = "L" Use Loglikelihood method (default is chisquare method)
//*-* = "U" Use a User specified fitting algorithm (via SetFCN)
//*-* = "Q" Quiet mode (minimum printing)
//*-* = "V" Verbose mode (default is between Q and V)
//*-* = "E" Perform better Errors estimation using Minos technique
//*-* = "R" Use the Range specified in the function range
//*-* = "N" Do not store the graphics function, do not draw
//*-* = "0" Do not plot the result of the fit. By default the fitted function
//*-* is drawn unless the option"N" above is specified.
//*-* = "+" Add this new fitted function to the list of fitted functions
//*-* (by default, any previous function is deleted)
//*-*
//*-* When the fit is drawn (by default), the parameter goption may be used
//*-* to specify a list of graphics options. See TH1::Draw for a complete
//*-* list of these options.
//*-*
//*-* In order to use the Range option, one must first create a function
//*-* with the expression to be fitted. For example, if your histogram
//*-* has a defined range between -4 and 4 and you want to fit a gaussian
//*-* only in the interval 1 to 3, you can do:
//*-* TF1 *f1 = new TF1("f1","gaus",1,3);
//*-* histo->Fit("f1","R");
//*-*
//*-* You can specify boundary limits for some or all parameters via
//*-* f1->SetParLimits(p_number, parmin, parmax);
//*-* if parmin>=parmax, the parameter is fixed
//*-* Note that you are not forced to fix the limits for all parameters.
//*-* For example, if you fit a function with 6 parameters, you can do:
//*-* func->SetParameters(0,3.1,1.e-6,0.1,-8,100);
//*-* func->SetParLimits(4,-10,-4);
//*-* func->SetParLimits(5, 1,1);
//*-* With this setup, parameters 0->3 can vary freely
//*-* Parameter 4 has boundaries [-10,-4] with initial value -8
//*-* Parameter 5 is fixed to 100.
//*-*
//*-* Note that option "I" gives better results but is slower.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Int_t i, npar,nvpar,nparx;
Double_t par, we, al, bl;
Double_t eplus,eminus,eparab,globcc,amin,edm,errdef,werr;
Double_t params[100], arglist[100];
TF1 *fnew1;
TF2 *fnew2;
TF3 *fnew3;
xfirst = fXaxis.GetFirst();
xlast = fXaxis.GetLast();
binwidx = fXaxis.GetBinWidth(xlast);
xmin = fXaxis.GetBinLowEdge(xfirst);
xmax = fXaxis.GetBinLowEdge(xlast) +binwidx;
yfirst = fYaxis.GetFirst();
ylast = fYaxis.GetLast();
binwidy = fYaxis.GetBinWidth(ylast);
ymin = fYaxis.GetBinLowEdge(yfirst);
ymax = fYaxis.GetBinLowEdge(ylast) +binwidy;
zfirst = fZaxis.GetFirst();
zlast = fZaxis.GetLast();
binwidz = fZaxis.GetBinWidth(zlast);
zmin = fZaxis.GetBinLowEdge(zfirst);
zmax = fZaxis.GetBinLowEdge(zlast) +binwidz;
xaxis = &fXaxis;
yaxis = &fYaxis;
zaxis = &fZaxis;
//*-*- Check if Minuit is initialized and create special functions
hFitter = TVirtualFitter::Fitter(this);
hFitter->Clear();
//*-*- Get pointer to the function by searching in the list of functions in ROOT
gF1 = (TF1*)gROOT->GetFunction(fname);
if (!gF1) { Error("Fit", "Unknown function: %s",fname); return; }
npar = gF1->GetNpar();
if (npar <=0) { Error("Fit", "Illegal number of parameters = %d",npar); return; }
//*-*- Check that function has same dimension as histogram
if (gF1->GetNdim() == 1 && GetDimension() > 1) {
Error("Fit", "Function %s is not 2-D",fname); return; }
if (gF1->GetNdim() == 2 && GetDimension() < 2) {
Error("Fit", "Function %s is not 1-D",fname); return; }
if (xxmin != xxmax) gF1->SetRange(xxmin,ymin,zmin,xxmax,ymax,zmax);
//*-*- Decode list of options into Foption
if (!FitOptionsMake(option)) return;
if (xxmin != xxmax) {
gF1->SetRange(xxmin,ymin,zmin,xxmax,ymax,zmax);
Foption.Range = 1;
}
//*-*- Is a Fit range specified?
Float_t fxmin, fymin, fzmin, fxmax, fymax, fzmax;
if (Foption.Range) {
gF1->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;
xfirst = fXaxis.FindBin(xmin); if (xfirst < 1) xfirst = 1;
xlast = fXaxis.FindBin(xmax); if (xlast > fXaxis.GetLast()) xlast = fXaxis.GetLast();
yfirst = fYaxis.FindBin(ymin); if (yfirst < 1) yfirst = 1;
ylast = fYaxis.FindBin(ymax); if (ylast > fYaxis.GetLast()) ylast = fYaxis.GetLast();
zfirst = fZaxis.FindBin(zmin); if (zfirst < 1) zfirst = 1;
zlast = fZaxis.FindBin(zmax); if (zlast > fZaxis.GetLast()) zlast = fZaxis.GetLast();
} else {
gF1->SetRange(xmin,ymin,zmin,xmax,ymax,zmax);
}
//*-*- If case of a predefined function, then compute initial values of parameters
Int_t special = gF1->GetNumber();
if (special == 100) H1InitGaus();
else if (special == 400) H1InitGaus();
else if (special == 200) H1InitExpo();
else if (special == 299+npar) H1InitPolynom();
//*-*- Some initialisations
if (!Foption.Verbose) {
arglist[0] = -1;
hFitter->ExecuteCommand("SET PRINT", arglist,1);
arglist[0] = 0;
hFitter->ExecuteCommand("SET NOW", arglist,0);
}
//*-*- Set error criterion for chisquare or likelihood methods
//*-*- MINUIT ERRDEF should not be set to 0.5 in case of loglikelihood fit.
//*-*- because the FCN is already multiplied by 2 in H1FitLikelihood
//*-*- if Hoption.User is specified, assume that the user has already set
//*-*- his minimization function via SetFCN.
arglist[0] = 1;
if (Foption.Like) {
hFitter->SetFCN(H1FitLikelihood);
} else {
if (!Foption.User) hFitter->SetFCN(H1FitChisquare);
}
hFitter->ExecuteCommand("SET ERR",arglist,1);
//*-*- Transfer names and initial values of parameters to Minuit
Int_t nfixed = 0;
for (i=0;i<npar;i++) {
par = gF1->GetParameter(i);
gF1->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,gF1->GetParName(i),par,we,al,bl);
}
if(nfixed > 0)hFitter->ExecuteCommand("FIX",arglist,nfixed); // Otto
//*-*- Set Gradient
if (Foption.Gradient) {
if (Foption.Gradient == 1) arglist[0] = 1;
else arglist[0] = 0;
hFitter->ExecuteCommand("SET GRAD",arglist,1);
}
//*-*- Reset Print level
if (Foption.Verbose) {
arglist[0] = 0; hFitter->ExecuteCommand("SET PRINT", arglist,1);
}
//*-*- Perform minimization
arglist[0] = hFitter->GetMaxIterations();
arglist[1] = 1; //epsilon
hFitter->ExecuteCommand("MIGRAD",arglist,2);
if (Foption.Errors) {
hFitter->ExecuteCommand("HESSE",arglist,0);
hFitter->ExecuteCommand("MINOS",arglist,0);
}
//*-*- Get return status
char parName[50];
for (i=0;i<npar;i++) {
hFitter->GetParameter(i,parName, par,we,al,bl);
if (Foption.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;
gF1->SetParameter(i,par);
gF1->SetParError(i,werr);
}
hFitter->GetStats(amin,edm,errdef,nvpar,nparx);
//*-*- Print final values of parameters.
if (!Foption.Quiet) {
if (Foption.Errors) hFitter->PrintResults(4,amin);
else hFitter->PrintResults(3,amin);
}
//*-* If Log Likelihood, compute an equivalent chisquare
if (Foption.Like) H1FitChisquare(npar, params, amin, params, 1);
gF1->SetChisquare(amin);
//*-*- Store fitted function in histogram functions list and draw
if (!Foption.Nostore) {
if (!Foption.Plus) fFunctions->Delete();
if (GetDimension() < 2) {
fnew1 = new TF1();
gF1->Copy(*fnew1);
fFunctions->Add(fnew1);
fnew1->SetParent(this);
fnew1->Save(xmin,xmax);
} else if (GetDimension() < 3) {
fnew2 = new TF2();
gF1->Copy(*fnew2);
fFunctions->Add(fnew2);
fnew2->SetParent(this);
} else {
fnew3 = new TF3();
gF1->Copy(*fnew3);
fFunctions->Add(fnew3);
fnew3->SetParent(this);
}
if (TestBit(kCanDelete)) return;
if (!Foption.Nograph && GetDimension() < 3) Draw(goption);
}
}
//______________________________________________________________________________
void TH1::FitPanel()
{
//*-*-*-*-*-*-*Display a panel with all histogram fit options*-*-*-*-*-*
//*-* ==============================================
//*-*
//*-* See class TFitPanel for example
if (fPainter) fPainter->FitPanel();
}
//______________________________________________________________________________
void TH1::FitSlicesX(TF1 *f1, Int_t binmin, Int_t binmax, Int_t cut, Option_t *option)
{
// Project slices along X in case of a 2-D histogram, then fit each slice
// with function f1 and make a histogram for each fit parameter
// Only bins along Y between binmin and binmax are considered.
// if f1=0, a gaussian is assumed
// Before invoking this function, one can set a subrange to be fitted along X
// via f1->SetRange(xmin,xmax)
// The argument option (default="QNR") can be used to change the fit options.
// "Q" means Quiet mode
// "N" means do not show the result of the fit
// "R" means fit the function in the specified function range
//
// Note that the generated histograms are added to the list of objects
// in the current directory. It is the user's responsability to delete
// these histograms.
//
// Example: Assume a 2-d histogram h2
// Root > h2->FitSlicesX(); produces 4 TH1D histograms
// with h2_0 containing parameter 0(Constant) for a Gaus fit
// of each bin in Y projected along X
// with h2_1 containing parameter 1(Mean) for a gaus fit
// with h2_2 containing parameter 2(RMS) for a gaus fit
// with h2_chi2 containing the chisquare/number of degrees of freedom for a gaus fit
//
// Root > h2->FitSlicesX(0,15,22,10);
// same as above, but only for bins 15 to 22 along Y
// and only for bins in Y for which the corresponding projection
// along X has more than cut bins filled.
//
//only valid for 2-D histogram
if (fDimension != 2) {
Error("FitSlicesX", "Function only valid for 2-D histograms");
return;
}
Int_t nbins = fYaxis.GetNbins();
if (binmin < 1) binmin = 1;
if (binmax > nbins) binmax = nbins;
if (binmax < binmin) {binmin = 1; binmax = nbins;}
//default is to fit with a gaussian
if (f1 == 0) {
f1 = (TF1*)gROOT->GetFunction("gaus");
if (f1 == 0) f1 = new TF1("gaus","gaus",fXaxis.GetXmin(),fXaxis.GetXmax());
else f1->SetRange(fXaxis.GetXmin(),fXaxis.GetXmax());
}
const char *fname = f1->GetName();
Int_t npar = f1->GetNpar();
Double_t *parsave = new Double_t[npar];
f1->GetParameters(parsave);
//Create one histogram for each function parameter
Int_t ipar;
char name[80], title[80];
TH1D *hlist[25];
for (ipar=0;ipar<npar;ipar++) {
sprintf(name,"%s_%d",GetName(),ipar);
sprintf(title,"Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar));
hlist[ipar] = new TH1D(name,title, nbins, fYaxis.GetXmin(), fYaxis.GetXmax());
}
sprintf(name,"%s_chi2",GetName());
TH1D *hchi2 = new TH1D(name,"chisquare", nbins, fYaxis.GetXmin(), fYaxis.GetXmax());
//Loop on all bins in Y, generate a projection along X
Int_t bin;
Int_t nentries;
for (bin=binmin;bin<=binmax;bin++) {
TH1D *hpx = ProjectionX("_temp",bin,bin,"e");
if (hpx == 0) continue;
nentries = Int_t(hpx->GetEntries());
if (nentries == 0 || nentries < cut) {delete hpx; continue;}
f1->SetParameters(parsave);
hpx->Fit(fname,option);
Int_t npfits = f1->GetNumberFitPoints();
if (npfits > npar && npfits >= cut) {
for (ipar=0;ipar<npar;ipar++) {
hlist[ipar]->Fill(GetBinCenter(bin),f1->GetParameter(ipar));
hlist[ipar]->SetBinError(bin,f1->GetParError(ipar));
}
hchi2->Fill(GetBinCenter(bin),f1->GetChisquare()/(npfits-npar));
}
delete hpx;
}
delete [] parsave;
}
//______________________________________________________________________________
void TH1::FitSlicesY(TF1 *f1, Int_t binmin, Int_t binmax, Int_t cut, Option_t *option)
{
// Project slices along Y in case of a 2-D histogram, then fit each slice
// with function f1 and make a histogram for each fit parameter
// Only bins along X between binmin and binmax are considered.
// if f1=0, a gaussian is assumed
// Before invoking this function, one can set a subrange to be fitted along Y
// via f1->SetRange(ymin,ymax)
// The argument option (default="QNR") can be used to change the fit options.
// "Q" means Quiet mode
// "N" means do not show the result of the fit
// "R" means fit the function in the specified function range
//
// Note that the generated histograms are added to the list of objects
// in the current directory. It is the user's responsability to delete
// these histograms.
//
// Example: Assume a 2-d histogram h2
// Root > h2->FitSlicesY(); produces 4 TH1D histograms
// with h2_0 containing parameter 0(Constant) for a Gaus fit
// of each bin in X projected along Y
// with h2_1 containing parameter 1(Mean) for a gaus fit
// with h2_2 containing parameter 2(RMS) for a gaus fit
// with h2_chi2 containing the chisquare/number of degrees of freedom for a gaus fit
//
// Root > h2->FitSlicesY(0,15,22,10);
// same as above, but only for bins 15 to 22 along X
// and only for bins in X for which the corresponding projection
// along Y has more than cut bins filled.
//
// A complete example of this function is given in tutorial:fitslicesy.C
// with the following output:
//
/*
*/
//
//only valid for 2-D histogram
if (fDimension != 2) {
Error("FitSlicesY", "Function only valid for 2-D histograms");
return;
}
Int_t nbins = fXaxis.GetNbins();
if (binmin < 1) binmin = 1;
if (binmax > nbins) binmax = nbins;
if (binmax < binmin) {binmin = 1; binmax = nbins;}
//default is to fit with a gaussian
if (f1 == 0) {
f1 = (TF1*)gROOT->GetFunction("gaus");
if (f1 == 0) f1 = new TF1("gaus","gaus",fYaxis.GetXmin(),fYaxis.GetXmax());
else f1->SetRange(fYaxis.GetXmin(),fYaxis.GetXmax());
}
const char *fname = f1->GetName();
Int_t npar = f1->GetNpar();
Double_t *parsave = new Double_t[npar];
f1->GetParameters(parsave);
//Create one histogram for each function parameter
Int_t ipar;
char name[80], title[80];
TH1D *hlist[25];
for (ipar=0;ipar<npar;ipar++) {
sprintf(name,"%s_%d",GetName(),ipar);
sprintf(title,"Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar));
hlist[ipar] = new TH1D(name,title, nbins, fXaxis.GetXmin(), fXaxis.GetXmax());
}
sprintf(name,"%s_chi2",GetName());
TH1D *hchi2 = new TH1D(name,"chisquare", nbins, fXaxis.GetXmin(), fXaxis.GetXmax());
//Loop on all bins in X, generate a projection along Y
Int_t bin;
Int_t nentries;
for (bin=binmin;bin<=binmax;bin++) {
TH1D *hpy = ProjectionY("_temp",bin,bin,"e");
if (hpy == 0) continue;
nentries = Int_t(hpy->GetEntries());
if (nentries == 0 || nentries < cut) {delete hpy; continue;}
f1->SetParameters(parsave);
hpy->Fit(fname,option);
Int_t npfits = f1->GetNumberFitPoints();
if (npfits > npar && npfits >= cut) {
for (ipar=0;ipar<npar;ipar++) {
hlist[ipar]->Fill(GetBinCenter(bin),f1->GetParameter(ipar));
hlist[ipar]->SetBinError(bin,f1->GetParError(ipar));
}
hchi2->Fill(GetBinCenter(bin),f1->GetChisquare()/(npfits-npar));
}
delete hpy;
}
delete [] parsave;
}
//______________________________________________________________________________
void TH1::FitSlicesZ(TF1 *f1, Int_t binminx, Int_t binmaxx, Int_t binminy, Int_t binmaxy, Int_t cut, Option_t *option)
{
// Project slices along Z in case of a 3-D histogram, then fit each slice
// with function f1 and make a 2-d histogram for each fit parameter
// Only cells in the bin range [binminx,binmaxx] and [binminy,binmaxy] are considered.
// if f1=0, a gaussian is assumed
// Before invoking this function, one can set a subrange to be fitted along Z
// via f1->SetRange(zmin,zmax)
// The argument option (default="QNR") can be used to change the fit options.
// "Q" means Quiet mode
// "N" means do not show the result of the fit
// "R" means fit the function in the specified function range
//
//
// Example: Assume a 3-d histogram h3
// Root > h3->FitSlicesZ(); produces 4 TH2D histograms
// with h3_0 containing parameter 0(Constant) for a Gaus fit
// of each cell in X,Y projected along Z
// with h3_1 containing parameter 1(Mean) for a gaus fit
// with h3_2 containing parameter 2(RMS) for a gaus fit
// with h3_chi2 containing the chisquare/number of degrees of freedom for a gaus fit
//
// Root > h3->Fit(0,15,22,0,0,10);
// same as above, but only for bins 15 to 22 along X
// and only for cells in X,Y for which the corresponding projection
// along Z has more than cut bins filled.
//only valid for 3-D histogram
if (fDimension != 3) {
Error("FitSlicesZ", "Function only valid for 3-D histograms");
return;
}
Int_t nbinsx = fXaxis.GetNbins();
Int_t nbinsy = fYaxis.GetNbins();
Int_t nbinsz = fZaxis.GetNbins();
if (binminx < 1) binminx = 1;
if (binmaxx > nbinsx) binmaxx = nbinsx;
if (binmaxx < binminx) {binminx = 1; binmaxx = nbinsx;}
if (binminy < 1) binminy = 1;
if (binmaxy > nbinsy) binmaxy = nbinsy;
if (binmaxy < binminy) {binminy = 1; binmaxy = nbinsy;}
//default is to fit with a gaussian
if (f1 == 0) {
f1 = (TF1*)gROOT->GetFunction("gaus");
if (f1 == 0) f1 = new TF1("gaus","gaus",fZaxis.GetXmin(),fZaxis.GetXmax());
else f1->SetRange(fZaxis.GetXmin(),fZaxis.GetXmax());
}
const char *fname = f1->GetName();
Int_t npar = f1->GetNpar();
Double_t *parsave = new Double_t[npar];
f1->GetParameters(parsave);
//Create one 2-d histogram for each function parameter
Int_t ipar;
char name[80], title[80];
TH2D *hlist[25];
for (ipar=0;ipar<npar;ipar++) {
sprintf(name,"%s_%d",GetName(),ipar);
sprintf(title,"Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar));
hlist[ipar] = new TH2D(name,title, nbinsx, fXaxis.GetXmin(), fXaxis.GetXmax()
, nbinsy, fYaxis.GetXmin(), fYaxis.GetXmax());
}
sprintf(name,"%s_chi2",GetName());
TH2D *hchi2 = new TH2D(name,"chisquare", nbinsx, fXaxis.GetXmin(), fXaxis.GetXmax()
, nbinsy, fYaxis.GetXmin(), fYaxis.GetXmax());
//Loop on all cells in X,Y generate a projection along Z
TH1D *hpz = new TH1D("R_temp","_temp",nbinsz, fZaxis.GetXmin(), fZaxis.GetXmax());
Int_t bin,binx,biny,binz;
for (biny=binminy;biny<=binmaxy;biny++) {
Float_t y = fYaxis.GetBinCenter(biny);
for (binx=binminx;binx<=binmaxx;binx++) {
Float_t x = fXaxis.GetBinCenter(binx);
hpz->Reset();
Int_t nfill = 0;
for (binz=1;binz<=nbinsz;binz++) {
bin = GetBin(binx,biny,binz);
Float_t w = GetBinContent(bin);
if (w == 0) continue;
hpz->Fill(fZaxis.GetBinCenter(binz),w);
hpz->SetBinError(binz,GetBinError(bin));
nfill++;
}
if (nfill < cut) continue;
f1->SetParameters(parsave);
hpz->Fit(fname,option);
Int_t npfits = f1->GetNumberFitPoints();
if (npfits > npar && npfits >= cut) {
for (ipar=0;ipar<npar;ipar++) {
hlist[ipar]->Fill(x,y,f1->GetParameter(ipar));
hlist[ipar]->SetCellError(binx,biny,f1->GetParError(ipar));
}
hchi2->Fill(x,y,f1->GetChisquare()/(npfits-npar));
}
}
}
delete [] parsave;
delete hpz;
}
//______________________________________________________________________________
Text_t *TH1::GetObjectInfo(Int_t px, Int_t py)
{
// Redefines TObject::GetObjectInfo.
// Displays the histogram info (bin number, contents, integral up to bin
// corresponding to cursor position px,py
//
return fPainter->GetObjectInfo(px,py);
}
//______________________________________________________________________________
Int_t TH1::FitOptionsMake(Option_t *choptin)
{
//*-*-*-*-*-*-*-*-*Decode string choptin and fill Foption structure*-*-*-*-*-*
//*-* ================================================
Foption.Quiet = 0;
Foption.Verbose = 0;
Foption.Bound = 0;
Foption.Like = 0;
Foption.User = 0;
Foption.W1 = 0;
Foption.Errors = 0;
Foption.Range = 0;
Foption.Gradient= 0;
Foption.Nograph = 0;
Foption.Nostore = 0;
Foption.Plus = 0;
Foption.Integral= 0;
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")) Foption.Quiet = 1;
if (strstr(chopt,"V")){Foption.Verbose = 1; Foption.Quiet = 0;}
if (strstr(chopt,"L")) Foption.Like = 1;
if (strstr(chopt,"W")) Foption.W1 = 1;
if (strstr(chopt,"E")) Foption.Errors = 1;
if (strstr(chopt,"R")) Foption.Range = 1;
if (strstr(chopt,"G")) Foption.Gradient= 1;
if (strstr(chopt,"N")) Foption.Nostore = 1;
if (strstr(chopt,"0")) Foption.Nograph = 1;
if (strstr(chopt,"+")) Foption.Plus = 1;
if (strstr(chopt,"I")) Foption.Integral= 1;
if (strstr(chopt,"U")){Foption.User = 1; Foption.Like = 0;}
return 1;
}
//______________________________________________________________________________
void H1FitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
{
//*-*-*-*-*-*Minimization function for H1s using a Chisquare method*-*-*-*-*-*
//*-* ======================================================
Double_t cu,eu,fu,fsum;
Double_t dersum[100], grad[100];
Double_t x[3];
Int_t bin,binx,biny,binz,k;
Float_t binlow, binup, binsize;
Int_t npfits = 0;
npar = gF1->GetNpar();
if (flag == 2) for (k=0;k<npar;k++) dersum[k] = gin[k] = 0;
TH1 *hfit = (TH1*)hFitter->GetObjectFit();
gF1->InitArgs(x,u);
f = 0;
for (binz=zfirst;binz<=zlast;binz++) {
x[2] = zaxis->GetBinCenter(binz);
for (biny=yfirst;biny<=ylast;biny++) {
x[1] = yaxis->GetBinCenter(biny);
for (binx=xfirst;binx<=xlast;binx++) {
bin = hfit->GetBin(binx,biny,binz);
cu = hfit->GetBinContent(bin);
x[0] = xaxis->GetBinCenter(binx);
if (Foption.Integral) {
binlow = xaxis->GetBinLowEdge(binx);
binsize = xaxis->GetBinWidth(binx);
binup = binlow + binsize;
fu = gF1->Integral(binlow,binup,u)/binsize;
} else {
fu = gF1->EvalPar(x,u);
}
if (Foption.W1) {
if (cu == 0) continue;
eu = 1;
} else {
eu = hfit->GetBinError(bin);
if (eu <= 0) continue;
}
if (flag == 2) {
for (k=0;k<npar;k++) dersum[k] += 1; //should be the derivative
}
npfits++;
if (flag == 2) {
for (k=0;k<npar;k++) grad[k] += dersum[k]*(fu-cu)/eu; dersum[k] = 0;
}
fsum = (cu-fu)/eu;
f += fsum*fsum;
//printf("binx=%d, x=%g, cu=%g, fu=%g, eu=%g, fsum=%g, f=%g, u1=%g, u2=%g, u3=%gn",binx,x[0],cu,fu,eu,fsum,f,u[0],u[1],u[2]);
}
}
}
gF1->SetNumberFitPoints(npfits);
}
//______________________________________________________________________________
void H1FitLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
{
//*-*-*-*-*-*Minimization function for H1s using a Likelihood method*-*-*-*-*-*
//*-* =======================================================
// Basically, it forms the likelihood by determining the Poisson
// probability that given a number of entries in a particualar bin,
// the fit would predict it's value. This is then done for each bin,
// and the sum of the logs is taken as the likelihood.
Double_t cu,fu,fobs,fsub;
Double_t dersum[100];
Double_t x[3];
Int_t bin,binx,biny,binz,k,icu;
Float_t binlow, binup, binsize;
Int_t npfits = 0;
npar = gF1->GetNpar();
if (flag == 2) for (k=0;k<npar;k++) dersum[k] = gin[k] = 0;
TH1 *hfit = (TH1*)hFitter->GetObjectFit();
gF1->InitArgs(x,u);
f = 0;
for (binz=zfirst;binz<=zlast;binz++) {
x[2] = zaxis->GetBinCenter(binz);
for (biny=yfirst;biny<=ylast;biny++) {
x[1] = yaxis->GetBinCenter(biny);
for (binx=xfirst;binx<=xlast;binx++) {
bin = hfit->GetBin(binx,biny,binz);
cu = hfit->GetBinContent(bin);
x[0] = xaxis->GetBinCenter(binx);
if (Foption.Integral) {
binlow = xaxis->GetBinLowEdge(binx);
binsize = xaxis->GetBinWidth(binx);
binup = binlow + binsize;
fu = gF1->Integral(binlow,binup,u)/binsize;
} else {
fu = gF1->EvalPar(x,u);
}
npfits++;
if (flag == 2) {
for (k=0;k<npar;k++) {
dersum[k] += 1; //should be the derivative
//grad[k] += dersum[k]*(fu-cu)/eu; dersum[k] = 0;
}
}
if (fu < 1.e-9) fu = 1.e-9;
icu = Int_t(cu);
fsub = -fu +icu*TMath::Log(fu);
fobs = hFitter->GetSumLog(icu);
fsub -= fobs;
f -= fsub;
}
}
}
f *= 2;
gF1->SetNumberFitPoints(npfits);
}
//______________________________________________________________________________
void H1InitGaus()
{
//*-*-*-*-*-*Compute Initial values of parameters for a gaussian*-*-*-*-*-*-*
//*-* ===================================================
Double_t allcha, sumx, sumx2, x, val, rms, mean;
Int_t bin;
static Double_t sqrtpi = 2.506628;
//*-*- Compute mean value and RMS of the histogram in the given range
TH1 *curHist = (TH1*)hFitter->GetObjectFit();
Double_t valmax = curHist->GetBinContent(xfirst);
allcha = sumx = sumx2 = 0;
for (bin=xfirst;bin<=xlast;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 = TMath::Sqrt(sumx2/allcha - mean*mean);
if (rms == 0) rms = binwidx*(xlast-xfirst+1)/4;
//if the distribution is really gaussian, the best approximation
//is binwidx*allcha/(sqrtpi*rms)
//However, in case of non-gaussian tails, this underestimates
//the normalisation constant. In this case the maximum value
//is a better approximation.
//We take the average of both quantities
Double_t constant = 0.5*(valmax+binwidx*allcha/(sqrtpi*rms));
//In case the mean value is outside the histo limits and
//the RMS is bigger than the range, we take
// mean = center of bins
// rms = half range
Float_t xmin = curHist->GetXaxis()->GetXmin();
Float_t xmax = curHist->GetXaxis()->GetXmax();
if ((mean < xmin || mean > xmax) && rms > (xmax-xmin)) {
mean = 0.5*(xmax+xmin);
rms = 0.5*(xmax-xmin);
}
gF1->SetParameter(0,constant);
gF1->SetParameter(1,mean);
gF1->SetParameter(2,rms);
gF1->SetParLimits(2,0,10*rms);
}
//______________________________________________________________________________
void H1InitExpo()
{
//*-*-*-*-*-*Compute Initial values of parameters for an exponential*-*-*-*-*
//*-* =======================================================
Double_t constant, slope;
Int_t ifail;
Int_t nchanx = xlast - xfirst + 1;
H1LeastSquareLinearFit(-nchanx, constant, slope, ifail);
gF1->SetParameter(0,constant);
gF1->SetParameter(1,slope);
}
//______________________________________________________________________________
void H1InitPolynom()
{
//*-*-*-*-*-*Compute Initial values of parameters for a polynom*-*-*-*-*-*-*
//*-* ===================================================
Double_t fitpar[25];
Int_t nchanx = xlast - xfirst + 1;
Int_t npar = gF1->GetNpar();
if (nchanx <=1 || npar == 1) {
TH1 *curHist = (TH1*)hFitter->GetObjectFit();
fitpar[0] = curHist->GetSumOfWeights()/Float_t(nchanx);
} else {
H1LeastSquareFit( nchanx, npar, fitpar);
}
for (Int_t i=0;i<npar;i++) gF1->SetParameter(i, fitpar[i]);
}
//______________________________________________________________________________
void H1LeastSquareFit(Int_t n, Int_t m, Double_t *a)
{
//*-*-*-*-*-*-*-*Least squares lpolynomial fitting without weights*-*-*-*-*-*-*
//*-* =================================================
//*-*
//*-* n number of points to fit
//*-* m number of parameters
//*-* a array of parameters
//*-*
//*-* based on CERNLIB routine LSQ: Translated to C++ by Rene Brun
//*-* (E.Keil. revised by B.Schorr, 23.10.1981.)
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
static Double_t zero = 0.;
static Double_t one = 1.;
static Int_t idim = 20;
static Double_t b[400] /* was [20][20] */;
static Int_t i, k, l, ifail;
static Double_t power;
static 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] = Float_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;
}
TH1 *curHist = (TH1*)hFitter->GetObjectFit();
for (k = xfirst; k <= xlast; ++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)
{
//*-*-*-*-*-*-*-*-*-*Least square linear fit without weights*-*-*-*-*-*-*-*-*
//*-* =======================================
//*-*
//*-* extracted from CERNLIB LLSQ: Translated to C++ by Rene Brun
//*-* (added to LSQ by B. Schorr, 15.02.1982.)
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
static Double_t xbar, ybar, x2bar;
static Int_t i, n;
static Double_t xybar;
static Float_t fn, xk, yk;
static Double_t det;
n = TMath::Abs(ndata);
ifail = -2;
xbar = ybar = x2bar = xybar = 0;
TH1 *curHist = (TH1*)hFitter->GetObjectFit();
for (i = xfirst; i <= xlast; ++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 = Float_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)
{
//*-*-*-*-*-*-*-*Extracted from CERN Program library routine DSEQN*-*-*-*-*-*
//*-* =================================================
//*-*
//*-* : Translated to C++ by Rene Brun
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Int_t a_dim1, a_offset, b_dim1, b_offset;
static Int_t nmjp1, i, j, l;
static Int_t im1, jp1, nm1, nmi;
static Double_t s1, s21, s22;
static Double_t one = 1.;
/* Parameter adjustments */
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)
{
//*-*-*-*-*-*Return Global bin number corresponding to binx,y,z*-*-*-*-*-*-*
//*-* ==================================================
//*-*
//*-* 2-D and 3-D histograms are represented with a one dimensional
//*-* structure.
//*-* This has the advantage that all existing functions, such as
//*-* GetBinContent, GetBinError, GetBinFunction work for all dimensions.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
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;
}
//______________________________________________________________________________
Axis_t TH1::GetRandom()
{
// return a random number distributed according the histogram bin contents.
// This function checks if the bins integral exists. If not, the integral
// is evaluated, normalized to one.
// The integral is automatically recomputed if the number of entries
// is not the same then when the integral was computed.
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 = ComputeIntegral();
} else {
integral = ComputeIntegral();
if (integral == 0 || fIntegral == 0) return 0;
}
Float_t r1 = gRandom->Rndm();
Int_t ibin = TMath::BinarySearch(nbinsx,&fIntegral[0],r1);
return GetBinLowEdge(ibin+1)
+GetBinWidth(ibin+1)*(fIntegral[ibin+1]-r1)/(fIntegral[ibin+1] - fIntegral[ibin]);
}
//______________________________________________________________________________
void TH1::GetRandom2(Axis_t &x, Axis_t &y)
{
// return 2 random numbers along axis x and y distributed according
// the cellcontents of a 2-dim histogram
if (fDimension != 2) {
Error("GetRandom","Function only valid for 2-d histograms");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbins = nbinsx*nbinsy;
Double_t integral;
if (fIntegral) {
if (fIntegral[nbins+1] != fEntries) integral = ComputeIntegral();
} else {
integral = ComputeIntegral();
if (integral == 0 || fIntegral == 0) return;
}
Float_t r1 = gRandom->Rndm();
Int_t ibin = TMath::BinarySearch(nbins,&fIntegral[0],r1);
Int_t biny = ibin/nbinsx;
Int_t binx = ibin - nbinsx*biny;
x = fXaxis.GetBinLowEdge(binx+1)
+fXaxis.GetBinWidth(binx+1)*(fIntegral[ibin+1]-r1)/(fIntegral[ibin+1] - fIntegral[ibin]);
y = fYaxis.GetBinLowEdge(biny+1) + fYaxis.GetBinWidth(biny+1)*gRandom->Rndm();
}
//______________________________________________________________________________
void TH1::GetRandom3(Axis_t &x, Axis_t &y, Axis_t &z)
{
// return 3 random numbers along axis x , y and z distributed according
// the cellcontents of a 3-dim histogram
if (fDimension != 3) {
Error("GetRandom","Function only valid for 3-d histograms");
return;
}
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
Int_t nxy = nbinsx*nbinsy;
Int_t nbins = nxy*nbinsz;
Double_t integral;
if (fIntegral) {
if (fIntegral[nbins+1] != fEntries) integral = ComputeIntegral();
} else {
integral = ComputeIntegral();
if (integral == 0 || fIntegral == 0) return;
}
Float_t r1 = gRandom->Rndm();
Int_t ibin = TMath::BinarySearch(nbins,&fIntegral[0],r1);
Int_t binz = ibin/nxy;
Int_t biny = (ibin - nxy*binz)/nbinsx;
Int_t binx = ibin - nbinsx*(biny + nbinsy*binz);
x = fXaxis.GetBinLowEdge(binx+1)
+fXaxis.GetBinWidth(binx+1)*(fIntegral[ibin+1]-r1)/(fIntegral[ibin+1] - fIntegral[ibin]);
y = fYaxis.GetBinLowEdge(biny+1) + fYaxis.GetBinWidth(biny+1)*gRandom->Rndm();
z = fZaxis.GetBinLowEdge(binz+1) + fZaxis.GetBinWidth(binz+1)*gRandom->Rndm();
}
//______________________________________________________________________________
Stat_t TH1::GetBinContent(Int_t)
{
//*-*-*-*-*-*-*Return content of global bin number bin*-*-*-*-*-*-*-*-*-*-*-*
//*-* =======================================
AbstractMethod("GetBinContent");
return 0;
}
//______________________________________________________________________________
void TH1::Multiply(TH1 *h1)
{
//*-*-*-*-*-*-*-*-*-*-*Multiply this histogram by h1*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* =============================
//
// this = this*h1
//
// If errors of this are available (TH1::Sumw2), errors are recalculated.
// Note that if h1 has Sumw2 set, Sumw2 is automatically called for this
// if not already set.
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();
//*-*- Check histogram compatibility
if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()) {
Error("Multiply","Attempt to multiply histograms with different number of bins");
return;
}
//*-*- Issue a Warning if histogram limits are different
if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
GetZaxis()->GetXmax() != h1->GetZaxis()->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;
//*-* Create Sumw2 if h1 has Sumw2 set
if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
//*-*- Reset statistics
fEntries = fTsumw = fTsumw2 = fTsumwx = fTsumwx2 = 0;
//*-*- Loop on bins (including underflows/overflows)
Int_t bin, binx, biny, binz;
Double_t c0,c1,w,z,x;
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);
z = TMath::Abs(w);
x = fXaxis.GetBinCenter(binx);
fEntries++;
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
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);
}
}
}
}
}
//______________________________________________________________________________
void TH1::Multiply(TH1 *h1, TH1 *h2, Float_t c1, Float_t c2, Option_t *option)
{
//*-*-*-*-*Replace contents of this histogram by multiplication of h1 by h2*-*
//*-* ================================================================
//
// this = (c1*h1)*(c2*h2)
//
// If errors of this are available (TH1::Sumw2), errors are recalculated.
// Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for this
// if not already set.
//
TString opt = option;
opt.ToLower();
// Bool_t binomial = kFALSE;
// if (opt.Contains("b")) binomial = kTRUE;
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();
//*-*- Check histogram compatibility
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;
}
//*-*- Issue a Warning if histogram limits are different
if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
Warning("Multiply","Attempt to multiply histograms with different axis limits");
}
if (GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin() ||
GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ||
GetYaxis()->GetXmin() != h2->GetYaxis()->GetXmin() ||
GetYaxis()->GetXmax() != h2->GetYaxis()->GetXmax() ||
GetZaxis()->GetXmin() != h2->GetZaxis()->GetXmin() ||
GetZaxis()->GetXmax() != h2->GetZaxis()->GetXmax()) {
Warning("Multiply","Attempt to multiply histograms with different axis limits");
}
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
//*-* Create Sumw2 if h1 or h2 have Sumw2 set
if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
//*-*- Reset statistics
fEntries = fTsumw = fTsumw2 = fTsumwx = fTsumwx2 = 0;
//*-*- Loop on bins (including underflows/overflows)
Int_t bin, binx, biny, binz;
Double_t b1,b2,w,z,x,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);
z = TMath::Abs(w);
x = fXaxis.GetBinCenter(binx);
fEntries++;
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
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);
}
}
}
}
}
//______________________________________________________________________________
void TH1::Paint(Option_t *option)
{
//*-*-*-*-*-*-*-*-*Control routine to paint any kind of histograms*-*-*-*-*-*-*
//*-* ===============================================
//
// This function is automatically called by TCanvas::Update.
// (see TH1::Draw for the list of options)
if (!fPainter) fPainter = TVirtualHistPainter::HistPainter(this);
if (fPainter) fPainter->Paint(option);
}
//______________________________________________________________________________
TProfile *TH1::ProfileX(const char *name, Int_t firstybin, Int_t lastybin, Option_t *option)
{
//*-*-*-*-*Project a 2-D histogram into a profile histogram along X*-*-*-*-*-*
//*-* ========================================================
//
// The projection is made from the channels along the Y axis
// ranging from firstybin to lastybin included.
//
// Check if this is a 2-d histogram
if (GetDimension() != 2) {
Error("ProfileX","Not a 2-D histogram");
return 0;
}
TString opt = option;
opt.ToLower();
Int_t nx = fXaxis.GetNbins();
Int_t ny = fYaxis.GetNbins();
if (firstybin < 0) firstybin = 0;
if (lastybin > ny) lastybin = ny;
// Create the profile histogram
char *pname = (char*)name;
if (strcmp(name,"_pfx") == 0) {
Int_t nch = strlen(GetName()) + 5;
pname = new char[nch];
sprintf(pname,"%s%s",GetName(),name);
}
TProfile *h1 = new TProfile(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax(),option);
if (pname != name) delete [] pname;
// Fill the profile histogram
Double_t cont;
Double_t entries = 0;
for (Int_t binx =0;binx<=nx;binx++) {
for (Int_t biny=firstybin;biny<=lastybin;biny++) {
cont = GetCellContent(binx,biny);
if (cont) {
h1->Fill(fXaxis.GetBinCenter(binx),fYaxis.GetBinCenter(biny), cont);
entries += cont;
}
}
}
h1->SetEntries(entries);
return h1;
}
//______________________________________________________________________________
TProfile *TH1::ProfileY(const char *name, Int_t firstxbin, Int_t lastxbin, Option_t *option)
{
//*-*-*-*-*Project a 2-D histogram into a profile histogram along Y*-*-*-*-*-*
//*-* ========================================================
//
// The projection is made from the channels along the X axis
// ranging from firstxbin to lastxbin included.
//
// Check if this is a 2-d histogram
if (GetDimension() != 2) {
Error("ProfileY","Not a 2-D histogram");
return 0;
}
TString opt = option;
opt.ToLower();
Int_t nx = fXaxis.GetNbins();
Int_t ny = fYaxis.GetNbins();
if (firstxbin < 0) firstxbin = 0;
if (lastxbin > nx) lastxbin = nx;
// Create the projection histogram
char *pname = (char*)name;
if (strcmp(name,"_pfy") == 0) {
Int_t nch = strlen(GetName()) + 5;
pname = new char[nch];
sprintf(pname,"%s%s",GetName(),name);
}
TProfile *h1 = new TProfile(pname,GetTitle(),nx,fYaxis.GetXmin(),fYaxis.GetXmax(),option);
if (pname != name) delete [] pname;
// Fill the profile histogram
Double_t cont;
Double_t entries = 0;
for (Int_t biny =0;biny<=ny;biny++) {
for (Int_t binx=firstxbin;binx<=lastxbin;binx++) {
cont = GetCellContent(binx,biny);
if (cont) {
h1->Fill(fYaxis.GetBinCenter(biny),fXaxis.GetBinCenter(binx), cont);
entries += cont;
}
}
}
h1->SetEntries(entries);
return h1;
}
//______________________________________________________________________________
TH1D *TH1::ProjectionX(const char *name, Int_t firstybin, Int_t lastybin, Option_t *option)
{
//*-*-*-*-*Project a 2-D histogram into a 1-D histogram along X*-*-*-*-*-*-*
//*-* ====================================================
//
// The projection is always of the type TH1D.
// The projection is made from the channels along the Y axis
// ranging from firstybin to lastybin included.
//
// if option "E" is specified, the errors are computed.
//
// Check if this is a 2-d histogram
if (GetDimension() != 2) {
Error("ProjectionX","Not a 2-D histogram");
return 0;
}
TString opt = option;
opt.ToLower();
Int_t nx = fXaxis.GetNbins();
Int_t ny = fYaxis.GetNbins();
if (firstybin < 0) firstybin = 0;
if (lastybin > ny) lastybin = ny;
// Create the projection histogram
char *pname = (char*)name;
if (strcmp(name,"_px") == 0) {
Int_t nch = strlen(GetName()) + 4;
pname = new char[nch];
sprintf(pname,"%s%s",GetName(),name);
}
TH1D *h1 = new TH1D(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax());
Bool_t computeErrors = kFALSE;
if (opt.Contains("e")) {h1->Sumw2(); computeErrors = kTRUE;}
if (pname != name) delete [] pname;
// Fill the projected histogram
Double_t cont,err,err2;
Double_t entries = 0;
for (Int_t binx =0;binx<=nx+1;binx++) {
err2 = 0;
for (Int_t biny=firstybin;biny<=lastybin;biny++) {
cont = GetCellContent(binx,biny);
err = GetCellError(binx,biny);
err2 += err*err;
if (cont) {
h1->Fill(fXaxis.GetBinCenter(binx), cont);
entries += cont;
}
}
if (computeErrors) h1->SetBinError(binx,TMath::Sqrt(err2));
}
h1->SetEntries(entries);
return h1;
}
//______________________________________________________________________________
TH1D *TH1::ProjectionY(const char *name, Int_t firstxbin, Int_t lastxbin, Option_t *option)
{
//*-*-*-*-*Project a 2-D histogram into a 1-D histogram along Y*-*-*-*-*-*-*
//*-* ====================================================
//
// The projection is always of the type TH1D.
// The projection is made from the channels along the X axis
// ranging from firstxbin to lastxbin included.
//
// if option "E" is specified, the errors are computed.
//
// Check if this is a 2-d histogram
if (GetDimension() != 2) {
Error("ProjectionY","Not a 2-D histogram");
return 0;
}
TString opt = option;
opt.ToLower();
Int_t nx = fXaxis.GetNbins();
Int_t ny = fYaxis.GetNbins();
if (firstxbin < 0) firstxbin = 0;
if (lastxbin > nx) lastxbin = nx;
// Create the projection histogram
char *pname = (char*)name;
if (strcmp(name,"_py") == 0) {
Int_t nch = strlen(GetName()) + 4;
pname = new char[nch];
sprintf(pname,"%s%s",GetName(),name);
}
TH1D *h1 = new TH1D(pname,GetTitle(),ny,fYaxis.GetXmin(),fYaxis.GetXmax());
Bool_t computeErrors = kFALSE;
if (opt.Contains("e")) {h1->Sumw2(); computeErrors = kTRUE;}
if (pname != name) delete [] pname;
// Fill the projected histogram
Double_t cont,err,err2;
Double_t entries = 0;
for (Int_t biny =0;biny<=ny+1;biny++) {
err2 = 0;
for (Int_t binx=firstxbin;binx<=lastxbin;binx++) {
cont = GetCellContent(binx,biny);
err = GetCellError(binx,biny);
err2 += err*err;
if (cont) {
h1->Fill(fYaxis.GetBinCenter(biny), cont);
entries += cont;
}
}
if (computeErrors) h1->SetBinError(biny,TMath::Sqrt(err2));
}
h1->SetEntries(entries);
return h1;
}
//______________________________________________________________________________
TH1D *TH1::ProjectionZ(const char *name, Int_t ixmin, Int_t ixmax, Int_t iymin, Int_t iymax, Option_t *option)
{
//*-*-*-*-*Project a 3-D histogram into a 1-D histogram along Z*-*-*-*-*-*-*
//*-* ====================================================
//
// The projection is always of the type TH1D.
// The projection is made from the cells along the X axis
// ranging from ixmin to ixmax and iymin to iymax included.
//
// if option "E" is specified, the errors are computed.
//
// code from Paola Collins & Hans Dijkstra
if (GetDimension() != 3){
Error("ProjectionZ", "Not a 3-D histogram");
return 0;
}
TString opt = option;
opt.ToLower();
Int_t nx = GetNbinsX();
Int_t ny = GetNbinsY();
Int_t nz = GetNbinsZ();
if (ixmin < 0) ixmin = 0;
if (ixmax > nx) ixmax = nx;
if (iymin < 0) iymin = 0;
if (iymax > nx) iymax = ny;
// Create the projection histogram
char *pname = (char*)name;
if (strcmp(name,"_pz") == 0) {
Int_t nch = strlen(GetName()) + 4;
pname = new char[nch];
sprintf(pname,"%s%s",GetName(),name);
}
TH1D *h1 = new TH1D(pname,GetTitle(),nz,fZaxis.GetXmin(),fZaxis.GetXmax());
Bool_t computeErrors = kFALSE;
if (opt.Contains("e")) {h1->Sumw2(); computeErrors = kTRUE;}
if (pname != name) delete [] pname;
// Fill the projected histogram
Float_t cont,e,e1;
Double_t entries = 0;
Double_t newerror = 0;
for (Int_t ixbin=ixmin;ixbin<=ixmax;ixbin++){
for (Int_t iybin=iymin;iybin<=iymax;iybin++){
for (Int_t binz=0;binz<=(nz+1);binz++){
Int_t bin = GetBin(ixbin,iybin,binz);
cont = GetBinContent(bin);
if (computeErrors) {
e = GetBinError(bin);
e1 = h1->GetBinError(binz);
newerror = TMath::Sqrt(e*e + e1*e1);
}
if (cont) {
h1->Fill(fZaxis.GetBinCenter(binz), cont);
entries += cont;
}
if (computeErrors) h1->SetBinError(binz,newerror);
}
}
}
h1->SetEntries(entries);
return h1;
}
//______________________________________________________________________________
TH1 *TH1::Project3D(Option_t *option)
{
// Project a 3-d histogram into 1 or 2-d histograms depending on the
// option parameter
// option may contain a combination of the characters x,y,z,e
// option = "x" return the x projection into a TH1D histogram
// option = "y" return the y projection into a TH1D histogram
// option = "z" return the z projection into a TH1D histogram
// option = "xy" return the x versus y projection into a TH2D histogram
// option = "yx" return the y versus x projection into a TH2D histogram
// option = "xz" return the x versus z projection into a TH2D histogram
// option = "zx" return the z versus x projection into a TH2D histogram
// option = "yz" return the y versus z projection into a TH2D histogram
// option = "zy" return the z versus y projection into a TH2D histogram
//
// If option contains the string "e", errors are computed
//
// The projection is made for the selected bins only
if (GetDimension() != 3) {
Error("Project3D", "Not a 3-D histogram");
return 0;
}
TString opt = option; opt.ToLower();
Int_t ixmin = fXaxis.GetFirst();
Int_t ixmax = fXaxis.GetLast();
Int_t iymin = fYaxis.GetFirst();
Int_t iymax = fYaxis.GetLast();
Int_t izmin = fZaxis.GetFirst();
Int_t izmax = fZaxis.GetLast();
Int_t nx = ixmax-ixmin+1;
Int_t ny = iymax-iymin+1;
Int_t nz = izmax-izmin+1;
Int_t pcase = 0;
if (opt.Contains("x")) pcase = 1;
if (opt.Contains("y")) pcase = 2;
if (opt.Contains("z")) pcase = 3;
if (opt.Contains("xy")) pcase = 4;
if (opt.Contains("yx")) pcase = 5;
if (opt.Contains("xz")) pcase = 6;
if (opt.Contains("zx")) pcase = 7;
if (opt.Contains("yz")) pcase = 8;
if (opt.Contains("zy")) pcase = 9;
// Create the projection histogram
TH1 *h = 0;
Int_t nch = strlen(GetName()) +opt.Length() +2;
char *name = new char[nch];
sprintf(name,"%s_%s",GetName(),option);
nch = strlen(GetTitle()) +opt.Length() +2;
char *title = new char[nch];
sprintf(title,"%s_%s",GetTitle(),option);
switch (pcase) {
case 1:
h = new TH1D(name,title,nx,fXaxis.GetBinLowEdge(ixmin),fXaxis.GetBinUpEdge(ixmax));
break;
case 2:
h = new TH1D(name,title,ny,fYaxis.GetBinLowEdge(iymin),fYaxis.GetBinUpEdge(iymax));
break;
case 3:
h = new TH1D(name,title,nz,fZaxis.GetBinLowEdge(izmin),fZaxis.GetBinUpEdge(izmax));
break;
case 4:
h = new TH2D(name,title,ny,fYaxis.GetBinLowEdge(iymin),fYaxis.GetBinUpEdge(iymax)
,nx,fXaxis.GetBinLowEdge(ixmin),fXaxis.GetBinUpEdge(ixmax));
break;
case 5:
h = new TH2D(name,title,nx,fXaxis.GetBinLowEdge(ixmin),fXaxis.GetBinUpEdge(ixmax)
,ny,fYaxis.GetBinLowEdge(iymin),fYaxis.GetBinUpEdge(iymax));
break;
case 6:
h = new TH2D(name,title,nz,fZaxis.GetBinLowEdge(izmin),fZaxis.GetBinUpEdge(izmax)
,nx,fXaxis.GetBinLowEdge(ixmin),fXaxis.GetBinUpEdge(ixmax));
break;
case 7:
h = new TH2D(name,title,nx,fXaxis.GetBinLowEdge(ixmin),fXaxis.GetBinUpEdge(ixmax)
,nz,fZaxis.GetBinLowEdge(izmin),fZaxis.GetBinUpEdge(izmax));
break;
case 8:
h = new TH2D(name,title,nz,fZaxis.GetBinLowEdge(izmin),fZaxis.GetBinUpEdge(izmax)
,ny,fYaxis.GetBinLowEdge(iymin),fYaxis.GetBinUpEdge(iymax));
break;
case 9:
h = new TH2D(name,title,ny,fYaxis.GetBinLowEdge(iymin),fYaxis.GetBinUpEdge(iymax)
,nz,fZaxis.GetBinLowEdge(izmin),fZaxis.GetBinUpEdge(izmax));
break;
}
delete [] name;
delete [] title;
if (h == 0) return h;
Bool_t computeErrors = kFALSE;
if (opt.Contains("e")) {h->Sumw2(); computeErrors = kTRUE;}
// Fill the projected histogram
Float_t cont,e,e1;
Double_t entries = 0;
Double_t newerror = 0;
for (Int_t ixbin=ixmin;ixbin<=ixmax;ixbin++){
Int_t ix = ixbin-ixmin+1;
for (Int_t iybin=iymin;iybin<=iymax;iybin++){
Int_t iy = iybin-iymin+1;
for (Int_t izbin=izmin;izbin<=izmax;izbin++){
Int_t iz = izbin-izmin+1;
Int_t bin = GetBin(ixbin,iybin,izbin);
cont = GetBinContent(bin);
switch (pcase) {
case 1:
if (cont) h->Fill(fXaxis.GetBinCenter(ixbin), cont);
if (computeErrors) {
e = GetBinError(bin);
e1 = h->GetBinError(ix);
newerror = TMath::Sqrt(e*e + e1*e1);
h->SetBinError(ix,newerror);
}
break;
case 2:
if (cont) h->Fill(fYaxis.GetBinCenter(iybin), cont);
if (computeErrors) {
e = GetBinError(bin);
e1 = h->GetBinError(iy);
newerror = TMath::Sqrt(e*e + e1*e1);
h->SetBinError(iy,newerror);
}
break;
case 3:
if (cont) h->Fill(fZaxis.GetBinCenter(izbin), cont);
if (computeErrors) {
e = GetBinError(bin);
e1 = h->GetBinError(iz);
newerror = TMath::Sqrt(e*e + e1*e1);
h->SetBinError(iz,newerror);
}
break;
case 4:
if (cont) h->Fill(fYaxis.GetBinCenter(iybin),fXaxis.GetBinCenter(ixbin), cont);
if (computeErrors) {
e = GetBinError(bin);
e1 = h->GetCellError(iy,ix);
newerror = TMath::Sqrt(e*e + e1*e1);
h->SetCellError(iy,ix,newerror);
}
break;
case 5:
if (cont) h->Fill(fXaxis.GetBinCenter(ixbin),fYaxis.GetBinCenter(iybin), cont);
if (computeErrors) {
e = GetBinError(bin);
e1 = h->GetCellError(ix,iy);
newerror = TMath::Sqrt(e*e + e1*e1);
h->SetCellError(ix,iy,newerror);
}
break;
case 6:
if (cont) h->Fill(fZaxis.GetBinCenter(izbin),fXaxis.GetBinCenter(ixbin), cont);
if (computeErrors) {
e = GetBinError(bin);
e1 = h->GetCellError(iz,ix);
newerror = TMath::Sqrt(e*e + e1*e1);
h->SetCellError(iz,ix,newerror);
}
break;
case 7:
if (cont) h->Fill(fXaxis.GetBinCenter(ixbin),fZaxis.GetBinCenter(izbin), cont);
if (computeErrors) {
e = GetBinError(bin);
e1 = h->GetCellError(ix,iz);
newerror = TMath::Sqrt(e*e + e1*e1);
h->SetCellError(ix,iz,newerror);
}
break;
case 8:
if (cont) h->Fill(fZaxis.GetBinCenter(izbin),fYaxis.GetBinCenter(iybin), cont);
if (computeErrors) {
e = GetBinError(bin);
e1 = h->GetCellError(iz,iy);
newerror = TMath::Sqrt(e*e + e1*e1);
h->SetCellError(iz,iy,newerror);
}
break;
case 9:
if (cont) h->Fill(fYaxis.GetBinCenter(iybin),fZaxis.GetBinCenter(izbin), cont);
if (computeErrors) {
e = GetBinError(bin);
e1 = h->GetCellError(iy,iz);
newerror = TMath::Sqrt(e*e + e1*e1);
h->SetCellError(iy,iz,newerror);
}
break;
}
if (cont) {
entries += cont;
}
}
}
}
h->SetEntries(entries);
return h;
}
//______________________________________________________________________________
TH1 *TH1::Rebin(Int_t ngroup, const char*newname)
{
//*-*-*-*-*Rebin this histogram grouping ngroup bins together*-*-*-*-*-*-*-*-*
//*-* ==================================================
// if newname is not blank a new temporary histogram hnew is created.
// else the current histogram is modified (default)
// The parameter ngroup indicates how many bins of this have to me merged
// into one bin of hnew
// If the original histogram has errors stored (via Sumw2), the resulting
// histograms has new errors correctly calculated.
//
// examples: if h1 is an existing TH1F histogram with 100 bins
// h1->Rebin(); //merges two bins in one in h1: previous contents of h1 are lost
// h1->Rebin(5); //merges five bins in one in h1
// TH1F *hnew = h1->Rebin(5,"hnew"); // creates a new histogram hnew
// //merging 5 bins of h1 in one bin
// NOTE: This function is currently implemented only for 1-D histograms
Int_t nbins = fXaxis.GetNbins();
Float_t xmin = fXaxis.GetXmin();
Float_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;
// Save old bin contents into a new array
Double_t *oldBins = new Double_t[nbins];
Int_t bin, i;
for (bin=0;bin<nbins;bin++) oldBins[bin] = GetBinContent(bin+1);
Double_t *oldErrors = 0;
if (fSumw2.fN != 0) {
oldErrors = new Double_t[nbins];
for (bin=0;bin<nbins;bin++) oldErrors[bin] = GetBinError(bin+1);
}
// create a clone of the old histogram if newname is specified
TH1 *hnew = this;
if (strlen(newname) > 0) {
hnew = (TH1*)Clone();
hnew->SetName(newname);
}
// change axis specs and rebuild bin contents array
hnew->SetBins(newbins,xmin,xmax); //this also changes errors array (if any)
// copy merged bin contents (ignore under/overflows)
Int_t oldbin = 0;
Double_t binContent, binError;
for (bin = 0;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+1,binContent);
if (oldErrors) hnew->SetBinError(bin+1,TMath::Sqrt(binError));
oldbin += ngroup;
}
delete [] oldBins;
if (oldErrors) delete [] oldErrors;
return hnew;
}
//______________________________________________________________________________
void TH1::Scale(Float_t c1)
{
//*-*-*-*-*Multiply this histogram by a constant c1*-*-*-*-*-*-*-*-*
//*-* ========================================
//
// this = c1*this
//
// Note that both contents and errors(if any) are scaled.
// This function uses the services of TH1::Add
//
Double_t ent = fEntries;
Add(this,this,c1,0);
fEntries = ent;
//if contours set, must also scale contours
Int_t ncontours = GetContour();
if (ncontours == 0) return;
Float_t *levels = fContour.GetArray();
for (Int_t i=0;i<ncontours;i++) {
levels[i] *= c1;
}
}
// -------------------------------------------------------------------------
void TH1::SmoothArray(Int_t NN, Double_t *XX, Int_t ntimes)
{
// smooth array XX, translation of Hbook routine hsmoof.F
// based on algorithm 353QH twice presented by J. Friedman
// in Proc.of the 1974 CERN School of Computing, Norway, 11-24 August, 1974.
Int_t ii, jj, ik, jk, kk, nn1, nn2;
Double_t hh[5];
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++) {
// first copy original data into temp array
for (ii = 0; ii < NN; ii++) {
YY[ii] = XX[ii];
}
// do 353 i.e. running median 3, 5, and 3 in a single loop
for (kk = 0; kk < 3; kk++) {
ik = 0;
if (kk == 1) ik = 1;
nn1 = ik + 1;
nn2 = NN - ik - 1;
// do all elements beside the first and last point for median 3
// and first two and last 2 for median 5
for (ii = nn1; ii < nn2; ii++) {
for (jj = 0; jj < 3; jj++) {
hh[jj] = YY[ii + jj - 1];
}
ZZ[ii] = TH1::SmoothMedian(3 + 2*ik, hh);
}
if (kk == 0) { // first median 3
// first point
hh[0] = 3*YY[1] - 2*YY[2];
hh[1] = YY[0];
hh[2] = YY[2];
ZZ[0] = TH1::SmoothMedian(3, hh);
// last point
hh[0] = YY[NN - 2];
hh[1] = YY[NN - 1];
hh[2] = 3*YY[NN - 2] - 2*YY[NN - 3];
ZZ[NN - 1] = TH1::SmoothMedian(3, hh);
}
if (kk == 1) { // median 5
// first point remains the same
ZZ[0] = YY[0];
for (ii = 0; ii < 3; ii++) {
hh[ii] = YY[ii];
}
ZZ[1] = TH1::SmoothMedian(3, hh);
// last two points
for (ii = 0; ii < 3; ii++) {
hh[ii] = YY[NN - 3 + ii];
}
ZZ[NN - 2] = TH1::SmoothMedian(3, hh);
ZZ[NN - 1] = YY[NN - 1];
}
}
// quadratic interpolation for flat segments
nn2 = nn2 - 1;
for (ii = nn1 + 1; ii < nn2; 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 = 0;
if ( TMath::Abs(hh[1]) > TMath::Abs(hh[0]) ) jk = -1;
YY[ii] = 0.5*ZZ[ii - 2*jk] + ZZ[ii - jk]/0.75 + ZZ[ii + 2*jk] /6.;
YY[ii + jk] = 0.5*(ZZ[ii + 2*jk] - ZZ[ii - 2*jk]) + ZZ[ii - jk];
}
// running means
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];
// now do the same for residuals
for (ii = 0; ii < NN; ii++) {
YY[ii] = XX[ii] - RR[ii];
}
// do 353 i.e. running median 3, 5, and 3 in a single loop
for (kk = 0; kk < 3; kk++) {
ik = 0;
if (kk == 1) ik = 1;
nn1 = ik + 1;
nn2 = NN - ik - 1;
// do all elements beside the first and last point for median 3
// and first two and last 2 for median 5
for (ii = nn1; ii < nn2; ii++) {
for (jj = 0; jj < 3; jj++) {
hh[jj] = YY[ii + jj - 1];
}
ZZ[ii] = TH1::SmoothMedian(3 + 2*ik, hh);
}
if (kk == 0) { // first median 3
// first point
hh[0] = 3*YY[1] - 2*YY[2];
hh[1] = YY[0];
hh[2] = YY[2];
ZZ[0] = TH1::SmoothMedian(3, hh);
// last point
hh[0] = YY[NN - 2];
hh[1] = YY[NN - 1];
hh[2] = 3*YY[NN - 2] - 2*YY[NN - 3];
ZZ[NN - 1] = TH1::SmoothMedian(3, hh);
}
if (kk == 1) { // median 5
// first point remains the same
ZZ[0] = YY[0];
for (ii = 0; ii < 3; ii++) {
hh[ii] = YY[ii];
}
ZZ[1] = TH1::SmoothMedian(3, hh);
// last two points
for (ii = 0; ii < 3; ii++) {
hh[ii] = YY[NN - 3 + ii];
}
ZZ[NN - 2] = TH1::SmoothMedian(3, hh);
ZZ[NN - 1] = YY[NN - 1];
}
}
// quadratic interpolation for flat segments
nn2 = nn2 - 1;
for (ii = nn1 + 1; ii < nn2; 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 = 0;
if ( TMath::Abs(hh[1]) > TMath::Abs(hh[0]) ) jk = -1;
YY[ii] = 0.5*ZZ[ii - 2*jk] + ZZ[ii - jk]/0.75 + ZZ[ii + 2*jk]/6.;
YY[ii + jk] = 0.5*(ZZ[ii + 2*jk] - ZZ[ii - 2*jk]) + ZZ[ii - jk];
}
// running means
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];
// add smoothed XX and smoothed residuals
for (ii = 0; ii < NN; ii++) {
XX[ii] = RR[ii] + ZZ[ii];
}
}
delete [] YY;
delete [] ZZ;
delete [] RR;
}
// ------------------------------------------------------------------------
Double_t TH1::SmoothMedian(Int_t n, Double_t *a)
{
// return the median of a vector a in monotonic order with length n
// where median is a number which divides sequence of n numbers
// into 2 halves. When n is odd, the median is kth element k = (n + 1) / 2.
// when n is even the median is a mean of the elements k = n/2 and k = n/2 + 1.
Int_t in, imin, imax;
Double_t xm;
if (n%2 == 0) in = n / 2;
else in = n / 2 + 1;
// find array element with maximum content
imax = TMath::LocMax(n,a);
xm = a[imax];
while (in < n) {
imin = TMath::LocMin(n,a); // find array element with minimum content
a[imin] = xm;
in++;
}
imin = TMath::LocMin(n,a);
return a[imin];
}
// ------------------------------------------------------------------------
void TH1::Smooth(Int_t ntimes)
{
// Smooth bin contents of this histogram.
// bin contents are replaced by their smooth values.
// Errors (if any) are not modified.
// algorithm can only be applied to 1-d histograms
if (fDimension != 1) {
Error("Smooth","Smooth only supported for 1-d histograms");
return;
}
Int_t nbins = fXaxis.GetNbins();
Double_t *XX = new Double_t[nbins];
Int_t i;
for (i=0;i<nbins;i++) {
XX[i] = GetBinContent(i+1);
}
TH1::SmoothArray(nbins,XX,ntimes);
for (i=0;i<nbins;i++) {
SetBinContent(i+1,XX[i]);
}
if (gPad) gPad->Modified();
}
//_______________________________________________________________________
void TH1::Streamer(TBuffer &b)
{
//*-*-*-*-*-*-*-*-*Stream a class object*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* =====================
if (b.IsReading()) {
b.ReadVersion(); //Version_t v = b.ReadVersion();
TNamed::Streamer(b);
TAttLine::Streamer(b);
TAttFill::Streamer(b);
TAttMarker::Streamer(b);
b >> fNcells;
fXaxis.Streamer(b);
fYaxis.Streamer(b);
fZaxis.Streamer(b);
b >> fBarOffset;
b >> fBarWidth;
b >> fEntries;
b >> fTsumw;
b >> fTsumw2;
b >> fTsumwx;
b >> fTsumwx2;
b >> fMaximum;
b >> fMinimum;
b >> fNormFactor;
fContour.Streamer(b);
fSumw2.Streamer(b);
fOption.Streamer(b);
fFunctions->Delete();
fFunctions->Streamer(b);
if (!gROOT->ReadingBasket()) {
fDirectory = gDirectory;
if (!gDirectory->GetList()->FindObject(this)) gDirectory->Append(this);
}
} else {
b.WriteVersion(TH1::IsA());
TNamed::Streamer(b);
TAttLine::Streamer(b);
TAttFill::Streamer(b);
TAttMarker::Streamer(b);
b << fNcells;
fXaxis.Streamer(b);
fYaxis.Streamer(b);
fZaxis.Streamer(b);
b << fBarOffset;
b << fBarWidth;
b << fEntries;
b << fTsumw;
b << fTsumw2;
b << fTsumwx;
b << fTsumwx2;
b << fMaximum;
b << fMinimum;
b << fNormFactor;
fContour.Streamer(b);
fSumw2.Streamer(b);
fOption.Streamer(b);
fFunctions->Streamer(b);
}
}
//______________________________________________________________________________
void TH1::Print(Option_t *option)
{
//*-*-*-*-*-*-*Print some global quantities for this histogram*-*-*-*-*-*-*-*
//*-* ===============================================
//
// If option "all" is given, bin contents and errors are also printed.
//
printf( "TH1.Print Name= %s, Entries= %d, Total sum= %gn",GetName(),Int_t(fEntries),GetSumOfWeights());
if (strcasecmp(option,"all")) return;
Int_t bin, binx, biny, binz;
Stat_t w,e;
Float_t x,y,z;
if (fDimension == 1) {
for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
x = fXaxis.GetBinCenter(binx);
w = GetBinContent(binx);
e = GetBinError(binx);
if(fSumw2.fN) printf(" fSumw[%d]=%g, x=%g, error=%gn",binx,w,x,e);
else printf(" fSumw[%d]=%g, x=%gn",binx,w,x);
}
}
if (fDimension == 2) {
for (biny=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) {
y = fYaxis.GetBinCenter(biny);
for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();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=%gn",binx,biny,w,x,y,e);
else printf(" fSumw[%d][%d]=%g, x=%g, y=%gn",binx,biny,w,x,y);
}
}
}
if (fDimension == 3) {
for (binz=fZaxis.GetFirst();binz<=fZaxis.GetLast();binz++) {
z = fZaxis.GetBinCenter(binz);
for (biny=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) {
y = fYaxis.GetBinCenter(biny);
for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();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=%gn",binx,biny,binz,w,x,y,z,e);
else printf(" fSumw[%d][%d][%d]=%g, x=%g, y=%g, z=%gn",binx,biny,binz,w,x,y,z);
}
}
}
}
}
//______________________________________________________________________________
void TH1::Reset(Option_t *)
{
//*-*-*-*-*-*-*-*Reset this histogram: contents, errors, etc*-*-*-*-*-*-*-*
//*-* ===========================================
fTsumw = 0;
fTsumw2 = 0;
fTsumwx = 0;
fTsumwx2 = 0;
fEntries = 0;
fSumw2.Reset();
fFunctions->Delete();
fContour.Set(0);
if (fIntegral) {delete [] fIntegral; fIntegral = 0;}
}
//______________________________________________________________________________
void TH1::SavePrimitive(ofstream &out, Option_t *option)
{
// Save primitive as a C++ statement(s) on output stream out
char quote = '"';
out<<" "<<endl;
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()");
out<<" "<<GetName()<<"->Draw("
<<quote<<option<<quote<<");"<<endl;
}
//______________________________________________________________________________
void TH1::UseCurrentStyle()
{
//*-*-*-*-*-*Replace current attributes by current style*-*-*-*-*
//*-* ===========================================
fXaxis.ResetAttAxis("X");
fYaxis.ResetAttAxis("Y");
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());
}
//______________________________________________________________________________
Stat_t TH1::GetMean(Int_t axis)
{
//*-*-*-*-*-*-*-*Return mean value of this histogram along the X axis*-*-*-*-*
//*-* ====================================================
if (axis <1 || axis > 3) return 0;
Stat_t stats[10];
GetStats(stats);
if (stats[0] == 0) return 0;
Int_t ax = 2*axis;
return stats[ax]/stats[0];
}
//______________________________________________________________________________
Stat_t TH1::GetRMS(Int_t axis)
{
//*-*-*-*-*-*-*-*Return the Root Mean Square value of this histogram*-*-*-*-*
//*-* ===================================================
if (axis <1 || axis > 3) return 0;
Stat_t x, rms2, stats[10];
GetStats(stats);
if (stats[0] == 0) return 0;
Int_t ax = 2*axis;
x = stats[ax]/stats[0];
rms2 = TMath::Abs(stats[ax+1]/stats[0] -x*x);
return TMath::Sqrt(rms2);
}
//______________________________________________________________________________
void TH1::GetStats(Stat_t *stats)
{
// fill the array stats from the contents of this histogram
// The array stats must be correctly dimensionned in the calling program.
// stats[0] = sumw
// stats[1] = sumw2
// stats[2] = sumwx
// stats[3] = sumwx2
// stats[4] = sumwy if 2-d or 3-d
// stats[5] = sumwy2
// stats[6] = sumwz if 3-d
// stats[7] = sumwz2
//
// If no axis-subrange is specified (via TAxis::SetRange), the array stats
// is simply a copy of the statistics quantities computed at filling time.
// If a sub-range is specified, the function recomputes these quantities
// from the bin contents in the current axis range.
// Loop on bins (possibly including underflows/overflows)
Int_t bin, binx, biny, binz;
Stat_t w;
Float_t x,y,z;
if (fDimension == 1) {
if (fTsumw == 0 || fXaxis.TestBit(kAxisRange)) {
for (bin=0;bin<4;bin++) stats[bin] = 0;
for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
x = fXaxis.GetBinCenter(binx);
w = TMath::Abs(GetBinContent(binx));
stats[0] += w;
stats[1] += w*w;
stats[2] += w*x;
stats[3] += w*x*x;
}
} else {
stats[0] = fTsumw;
stats[1] = fTsumw2;
stats[2] = fTsumwx;
stats[3] = fTsumwx2;
}
}
if (fDimension == 2) {
for (bin=0;bin<6;bin++) stats[bin] = 0;
for (biny=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) {
y = fYaxis.GetBinCenter(biny);
for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
bin = GetBin(binx,biny);
x = fXaxis.GetBinCenter(binx);
w = TMath::Abs(GetBinContent(bin));
stats[0] += w;
stats[1] += w*w;
stats[2] += w*x;
stats[3] += w*x*x;
stats[4] += w*y;
stats[5] += w*y*y;
}
}
}
if (fDimension == 3) {
for (bin=0;bin<8;bin++) stats[bin] = 0;
for (binz=fZaxis.GetFirst();binz<=fZaxis.GetLast();binz++) {
z = fZaxis.GetBinCenter(binz);
for (biny=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) {
y = fYaxis.GetBinCenter(biny);
for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
bin = GetBin(binx,biny,binz);
x = fXaxis.GetBinCenter(binx);
w = TMath::Abs(GetBinContent(bin));
stats[0] += w;
stats[1] += w*w;
stats[2] += w*x;
stats[3] += w*x*x;
stats[4] += w*y;
stats[5] += w*y*y;
stats[6] += w*z;
stats[7] += w*z*z;
}
}
}
}
}
//______________________________________________________________________________
Stat_t TH1::GetSumOfWeights()
{
//*-*-*-*-*-*-*-*Return the sum of weights excluding under/overflows*-*-*-*-*
//*-* ===================================================
Int_t bin,binx,biny,binz;
Stat_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;
}
//______________________________________________________________________________
Stat_t TH1::Integral()
{
//Return integral of bin contents. Only bins in the bins range are considered.
Stat_t stats[10];
GetStats(stats);
return stats[0];
}
//______________________________________________________________________________
Stat_t TH1::Integral(Int_t binx1, Int_t binx2)
{
//Return integral of bin contents between binx1 and binx2 for a 1-D histogram
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
if (binx1 < 0) binx1 = 0;
if (binx2 > nbinsx+1) binx2 = nbinsx+1;
if (binx2 < binx1) binx2 = nbinsx;
Stat_t integral = 0;
//*-*- Loop on bins in specified range
Int_t bin, binx, biny, binz;
for (binz=0;binz<=nbinsz+1;binz++) {
for (biny=0;biny<=nbinsy+1;biny++) {
for (binx=binx1;binx<=binx2;binx++) {
bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
integral += GetBinContent(bin);
}
}
}
return integral;
}
//______________________________________________________________________________
Stat_t TH1::Integral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2)
{
//Return integral of bin contents in range [binx1,binx2],[biny1,biny2]
// for a 2-D histogram
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
if (binx1 < 0) binx1 = 0;
if (binx2 > nbinsx+1) binx2 = nbinsx+1;
if (binx2 < binx1) binx2 = nbinsx;
if (biny1 < 0) biny1 = 0;
if (biny2 > nbinsy+1) biny2 = nbinsy+1;
if (biny2 < biny1) biny2 = nbinsy;
Stat_t integral = 0;
//*-*- Loop on bins in specified range
Int_t bin, binx, biny, binz;
for (binz=0;binz<=nbinsz+1;binz++) {
for (biny=biny1;biny<=biny2;biny++) {
for (binx=binx1;binx<=binx2;binx++) {
bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
integral += GetBinContent(bin);
}
}
}
return integral;
}
//______________________________________________________________________________
Stat_t TH1::Integral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Int_t binz1, Int_t binz2)
{
//Return integral of bin contents in range [binx1,binx2],[biny1,biny2],[binz1,binz2]
// for a 3-D histogram
Int_t nbinsx = GetNbinsX();
Int_t nbinsy = GetNbinsY();
Int_t nbinsz = GetNbinsZ();
if (fDimension < 2) nbinsy = -1;
if (fDimension < 3) nbinsz = -1;
if (binx1 < 0) binx1 = 0;
if (binx2 > nbinsx+1) binx2 = nbinsx+1;
if (biny1 < 0) biny1 = 0;
if (biny2 > nbinsy+1) biny2 = nbinsy+1;
if (binz1 < 0) binz1 = 0;
if (binz2 > nbinsz+1) binz2 = nbinsz+1;
Stat_t integral = 0;
//*-*- Loop on bins in specified range
Int_t bin, binx, biny, binz;
for (binz=binz1;binz<=binz2;binz++) {
for (biny=biny1;biny<=biny2;biny++) {
for (binx=binx1;binx<=binx2;binx++) {
bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
integral += GetBinContent(bin);
}
}
}
return integral;
}
//______________________________________________________________________________
void TH1::SetContent(Stat_t *content)
{
//*-*-*-*-*-*-*-*Replace bin contents by the contents of array content*-*-*-*
//*-* =====================================================
Int_t bin;
Stat_t bincontent;
for (bin=0; bin<fNcells; bin++) {
bincontent = *(content + bin);
SetBinContent(bin, bincontent);
}
}
//______________________________________________________________________________
Int_t TH1::GetContour(Float_t *levels)
{
//*-*-*-*-*-*-*-*Return contour values into array levels*-*-*-*-*-*-*-*-*-*
//*-* =======================================
//*-*
//*-* The number of contour levels can be returned by getContourLevel
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
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;
}
//______________________________________________________________________________
Float_t TH1::GetContourLevel(Int_t level)
{
//*-*-*-*-*-*-*-*Return the number of contour levels*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* ===================================
if (level <0 || level >= fContour.fN) return 0;
Float_t zlevel = fContour.fArray[level];
return zlevel;
}
//______________________________________________________________________________
void TH1::SetContour(Int_t nlevels, Float_t *levels)
{
//*-*-*-*-*-*-*-*Set the number and values of contour levels*-*-*-*-*-*-*-*-*
//*-* ===========================================
//
// By default the number of contour levels is set to 20.
//
// if argument levels = 0 or issing, equidistant contours are computed
//
Int_t level;
ResetBit(kUserContour);
if (nlevels <=0 ) {
fContour.Set(0);
return;
}
fContour.Set(nlevels);
//*-*- Contour levels are specified
if (levels) {
SetBit(kUserContour);
for (level=0; level<nlevels; level++) fContour.fArray[level] = levels[level];
} else {
//*-*- contour levels are computed automatically as equidistant contours
Float_t zmin = GetMinimum();
Float_t zmax = GetMaximum();
Float_t dz = (zmax-zmin)/Float_t(nlevels);
if (gPad->GetLogz()) {
if (zmax <= 0) return;
if (zmin <= 0) zmin = 0.001*zmax;
zmin = TMath::Log10(zmin);
zmax = TMath::Log10(zmax);
dz = (zmax-zmin)/Float_t(nlevels);
}
for (level=0; level<nlevels; level++) {
fContour.fArray[level] = zmin + dz*Float_t(level);
}
}
}
//______________________________________________________________________________
void TH1::SetContourLevel(Int_t level, Float_t value)
{
//*-*-*-*-*-*-*-*-*-*-*Set value for one contour level*-*-*-*-*-*-*-*-*-*-*-*
//*-* ===============================
if (level <0 || level >= fContour.fN) return;
SetBit(kUserContour);
fContour.fArray[level] = value;
}
//______________________________________________________________________________
Float_t TH1::GetMaximum()
{
//*-*-*-*-*-*-*-*-*-*-*Return maximum value of all bins*-*-*-*-*-*-*-*-*-*-*-*
//*-* ================================
Int_t binx, biny;
Float_t maximum, value;
if (fMaximum != -1111) return fMaximum;
if (GetDimension() < 2) {
maximum = GetBinContent(1);
for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
value = GetBinContent(binx);
if (value > maximum) maximum = value;
}
} else {
maximum = GetCellContent(1,1);
for(biny=1; biny<=fYaxis.GetNbins(); biny++) {
for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
value = GetCellContent(binx, biny);
if (value > maximum) maximum = value;
}
}
}
return maximum;
}
//______________________________________________________________________________
Float_t TH1::GetMinimum()
{
//*-*-*-*-*-*-*-*-*-*-*Return minimum value of all bins*-*-*-*-*-*-*-*-*-*-*-*
//*-* ================================
Int_t binx, biny;
Float_t minimum, value;
if (fMinimum != -1111) return fMinimum;
if (GetDimension() < 2) {
minimum = GetBinContent(1);
for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
value = GetBinContent(binx);
if (value < minimum) minimum = value;
}
} else {
minimum = GetCellContent(1,1);
for(biny=1; biny<=fYaxis.GetNbins(); biny++) {
for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
value = GetCellContent(binx, biny);
if (value < minimum) minimum = value;
}
}
}
return minimum;
}
//______________________________________________________________________________
void TH1::SetBins(Int_t nx, Float_t xmin, Float_t xmax)
{
//*-*-*-*-*-*-*-*-*Redefine x axis parameters*-*-*-*-*-*-*-*-*-*-*-*
//*-* ===========================
// The X axis parameters are modified.
// The bins content array is resized
// if errors (Sumw2) the errors array is resized
// The previous bin contents are lost
// To change only the axis limits, see TAxis::SetRange
if (GetDimension() != 1) {
Error("SetBins","Operation only valid for 1-d histograms");
return;
}
fXaxis.Set(nx,xmin,xmax);
fNcells = nx+2;
SetBinsLength(fNcells);
if (fSumw2.fN) {
fSumw2.Set(fNcells);
}
}
//______________________________________________________________________________
void TH1::SetBins(Int_t nx, Float_t xmin, Float_t xmax, Int_t ny, Float_t ymin, Float_t ymax)
{
//*-*-*-*-*-*-*-*-*Redefine x and y axis parameters*-*-*-*-*-*-*-*-*-*-*-*
//*-* =================================
// The X and Y axis parameters are modified.
// The bins content array is resized
// if errors (Sumw2) the errors array is resized
// The previous bin contents are lost
// To change only the axis limits, see TAxis::SetRange
if (GetDimension() != 2) {
Error("SetBins","Operation only valid for 2-d histograms");
return;
}
fXaxis.Set(nx,xmin,xmax);
fYaxis.Set(ny,ymin,ymax);
fNcells = (nx+2)*(ny+2);
SetBinsLength(fNcells);
if (fSumw2.fN) {
fSumw2.Set(fNcells);
}
}
//______________________________________________________________________________
void TH1::SetBins(Int_t nx, Float_t xmin, Float_t xmax, Int_t ny, Float_t ymin, Float_t ymax, Int_t nz, Float_t zmin, Float_t zmax)
{
//*-*-*-*-*-*-*-*-*Redefine x, y and z axis parameters*-*-*-*-*-*-*-*-*-*-*-*
//*-* ====================================
// The X, Y and Z axis parameters are modified.
// The bins content array is resized
// if errors (Sumw2) the errors array is resized
// The previous bin contents are lost
// To change only the axis limits, see TAxis::SetRange
if (GetDimension() != 3) {
Error("SetBins","Operation only valid for 3-d histograms");
return;
}
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(Float_t maximum)
{
//*-*-*-*-*-*-*-*-*Set the maximum value for the Y axis*-*-*-*-*-*-*-*-*-*-*-*
//*-* ====================================
// By default the maximum value is automatically set to the maximum
// bin content plyus a margin of 10 per cent.
//
fMaximum = maximum;
}
//______________________________________________________________________________
void TH1::SetMinimum(Float_t minimum)
{
//*-*-*-*-*-*-*-*-*Set the minimum value for the Y axis*-*-*-*-*-*-*-*-*-*-*-*
//*-* ====================================
// By default the minimum value is automatically set to zero if all bin contents
// are positive or the minimum - 10 per cent otherwise.
//
fMinimum = minimum;
}
//______________________________________________________________________________
void TH1::SetDirectory(TDirectory *dir)
{
// Remove reference to this histogram from current directory and add
// reference to new directory dir. dir can be 0 in which case the
// histogram does not belong to any directory.
if (fDirectory == dir) return;
if (fDirectory) fDirectory->GetList()->Remove(this);
fDirectory = dir;
if (fDirectory) fDirectory->GetList()->Add(this);
}
//______________________________________________________________________________
void TH1::SetError(Stat_t *error)
{
//*-*-*-*-*-*-*-*-*Replace bin errors by values in array error*-*-*-*-*-*-*-*-*
//*-* ===========================================
Int_t bin;
Stat_t binerror;
for (bin=0; bin<fNcells; bin++) {
binerror = error[bin];
SetBinError(bin, binerror);
}
}
//______________________________________________________________________________
void TH1::SetName(const Text_t *name)
{
//*-*-*-*-*-*-*-*-*-*-*Change the name of this histogram*-*-*-*-*-*-*-*-*-*-*
//*-* =================================
// Histograms are named objects in a THashList.
// We must update the hashlist if we change the name
if (fDirectory) fDirectory->GetList()->Remove(this);
fName = name;
if (fDirectory) fDirectory->GetList()->Add(this);
}
//______________________________________________________________________________
void TH1::SetStats(Bool_t stats)
{
//*-*-*-*-*-*-*-*-*Set statistics option on/off
//*-* ============================
// By default, the statistics box is drawn.
// The paint options can be selected via gStyle->SetOptStats.
// This function sets/resets the kNoStats bin in the histogram object.
// It has priority over the Style option.
ResetBit(kNoStats);
if (!stats) SetBit(kNoStats);
}
//______________________________________________________________________________
void TH1::Sumw2()
{
//*-*-*-*-*Create structure to store sum of squares of weights*-*-*-*-*-*-*-*
//*-* ===================================================
//*-*
//*-* if histogram is already filled, the sum of squares of weights
//*-* is filled with the existing bin contents
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
if (fSumw2.fN) {
Warning("Sumw2","Sum of squares of weights structure already created");
return;
}
fSumw2.Set(fNcells);
if (fEntries) {
for (Int_t bin=0; bin<fNcells; bin++) {
fSumw2.fArray[bin] = GetBinContent(bin);
}
}
}
//______________________________________________________________________________
TF1 *TH1::GetFunction(const Text_t *name)
{
//*-*-*-*-*Return pointer to function with name*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* ===================================
return (TF1*)fFunctions->FindObject(name);
}
//______________________________________________________________________________
Stat_t TH1::GetBinError(Int_t bin)
{
//*-*-*-*-*-*-*Return value of error associated to bin number bin*-*-*-*-*
//*-* ==================================================
//*-*
//*-* if the sum of squares of weights has been defined (via Sumw2),
//*-* this function returns the sqrt(sum of w2).
//*-* otherwise it returns the sqrt(contents) for this bin.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
if (fSumw2.fN) return TMath::Sqrt(fSumw2.fArray[bin]);
Stat_t error2 = TMath::Abs(GetBinContent(bin));
return TMath::Sqrt(error2);
}
//______________________________________________________________________________
Stat_t TH1::GetCellContent(Int_t binx, Int_t biny)
{
Stat_t binval;
if (GetDimension() > 1) {
binval = GetBinContent(biny*(fXaxis.GetNbins()+2) + binx);
return binval;
}
return GetBinContent(binx);
}
//______________________________________________________________________________
Stat_t TH1::GetCellError(Int_t binx, Int_t biny)
{
if (GetDimension() > 1) return GetBinError(biny*(fXaxis.GetNbins()+2) + binx);
return GetBinError(binx);
}
//______________________________________________________________________________
void TH1::SetBinError(Int_t bin, Stat_t error)
{
if (!fSumw2.fN) Sumw2();
if (bin <0 || bin>= fSumw2.fN) return;
fSumw2.fArray[bin] = error*error;
}
//______________________________________________________________________________
void TH1::SetCellContent(Int_t binx, Int_t biny, Stat_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::SetCellError(Int_t binx, Int_t biny, Stat_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, Stat_t)
{
AbstractMethod("SetBinContent");
}
ClassImp(TH1C)
//______________________________________________________________________________
// TH1C methods
//______________________________________________________________________________
TH1C::TH1C(): TH1(), TArrayC()
{
fDimension = 1;
}
//______________________________________________________________________________
TH1C::TH1C(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
: TH1(name,title,nbins,xlow,xup), TArrayC()
{
fDimension = dim;
}
//______________________________________________________________________________
TH1C::TH1C(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
: TH1(name,title,nbins,xbins), TArrayC()
{
fDimension = dim;
}
//______________________________________________________________________________
TH1C::TH1C(const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
: TH1(name,title,nbins,xlow,xup), TArrayC(nbins+2)
{
//
// Create a 1-Dim histogram with fix bins of type char (one byte per channel)
// ==========================================================================
// (see TH1::TH1 for explanation of parameters)
//
fDimension = 1;
}
//______________________________________________________________________________
TH1C::TH1C(const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
: TH1(name,title,nbins,xbins), TArrayC(nbins+2)
{
//
// Create a 1-Dim histogram with variable bins of type char (one byte per channel)
// ==========================================================================
// (see TH1::TH1 for explanation of parameters)
//
fDimension = 1;
}
//______________________________________________________________________________
TH1C::~TH1C()
{
}
//______________________________________________________________________________
TH1C::TH1C(const TH1C &h1c)
{
((TH1C&)h1c).Copy(*this);
}
//______________________________________________________________________________
void TH1C::AddBinContent(Int_t bin)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by 1*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* ==========================
if (fArray[bin] < 127) fArray[bin]++;
}
//______________________________________________________________________________
void TH1C::AddBinContent(Int_t bin, Stat_t w)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by 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)
{
TH1::Copy(newth1);
TArrayC::Copy((TH1C&)newth1);
}
//______________________________________________________________________________
TH1 *TH1C::DrawCopy(Option_t *option)
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
TH1C *newth1 = (TH1C*)Clone();
newth1->SetDirectory(0);
newth1->SetBit(kCanDelete);
newth1->AppendPad(option);
return newth1;
}
//______________________________________________________________________________
Stat_t TH1C::GetBinContent(Int_t bin)
{
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
return Stat_t (fArray[bin]);
}
//______________________________________________________________________________
void TH1C::Reset(Option_t *option)
{
TH1::Reset(option);
TArrayC::Reset();
}
//______________________________________________________________________________
TH1C& TH1C::operator=(const TH1C &h1)
{
if (this != &h1) ((TH1C&)h1).Copy(*this);
return *this;
}
//______________________________________________________________________________
TH1C operator*(Float_t c1, TH1C &h1)
{
TH1C hnew = h1;
hnew.Scale(c1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1C operator+(TH1C &h1, TH1C &h2)
{
TH1C hnew = h1;
hnew.Add(&h2,1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1C operator-(TH1C &h1, TH1C &h2)
{
TH1C hnew = h1;
hnew.Add(&h2,-1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1C operator*(TH1C &h1, TH1C &h2)
{
TH1C hnew = h1;
hnew.Multiply(&h2);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1C operator/(TH1C &h1, TH1C &h2)
{
TH1C hnew = h1;
hnew.Divide(&h2);
hnew.SetDirectory(0);
return hnew;
}
ClassImp(TH1S)
//______________________________________________________________________________
// TH1S methods
//______________________________________________________________________________
TH1S::TH1S(): TH1(), TArrayS()
{
fDimension = 1;
}
//______________________________________________________________________________
TH1S::TH1S(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
: TH1(name,title,nbins,xlow,xup), TArrayS()
{
fDimension = dim;
}
//______________________________________________________________________________
TH1S::TH1S(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
: TH1(name,title,nbins,xbins), TArrayS()
{
fDimension = dim;
}
//______________________________________________________________________________
TH1S::TH1S(const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
: TH1(name,title,nbins,xlow,xup), TArrayS(nbins+2)
{
//
// Create a 1-Dim histogram with fix bins of type short
// ====================================================
// (see TH1::TH1 for explanation of parameters)
//
fDimension = 1;
}
//______________________________________________________________________________
TH1S::TH1S(const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
: TH1(name,title,nbins,xbins), TArrayS(nbins+2)
{
//
// Create a 1-Dim histogram with variable bins of type short
// =========================================================
// (see TH1::TH1 for explanation of parameters)
//
fDimension = 1;
}
//______________________________________________________________________________
TH1S::~TH1S()
{
}
//______________________________________________________________________________
TH1S::TH1S(const TH1S &h1s)
{
((TH1S&)h1s).Copy(*this);
}
//______________________________________________________________________________
void TH1S::AddBinContent(Int_t bin)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by 1*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* ==========================
if (fArray[bin] < 32767) fArray[bin]++;
}
//______________________________________________________________________________
void TH1S::AddBinContent(Int_t bin, Stat_t w)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by 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)
{
TH1::Copy(newth1);
TArrayS::Copy((TH1S&)newth1);
}
//______________________________________________________________________________
TH1 *TH1S::DrawCopy(Option_t *option)
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
TH1S *newth1 = (TH1S*)Clone();
newth1->SetDirectory(0);
newth1->SetBit(kCanDelete);
newth1->AppendPad(option);
return newth1;
}
//______________________________________________________________________________
Stat_t TH1S::GetBinContent(Int_t bin)
{
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
return Stat_t (fArray[bin]);
}
//______________________________________________________________________________
void TH1S::Reset(Option_t *option)
{
TH1::Reset(option);
TArrayS::Reset();
}
//______________________________________________________________________________
TH1S& TH1S::operator=(const TH1S &h1)
{
if (this != &h1) ((TH1S&)h1).Copy(*this);
return *this;
}
//______________________________________________________________________________
TH1S operator*(Float_t c1, TH1S &h1)
{
TH1S hnew = h1;
hnew.Scale(c1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1S operator+(TH1S &h1, TH1S &h2)
{
TH1S hnew = h1;
hnew.Add(&h2,1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1S operator-(TH1S &h1, TH1S &h2)
{
TH1S hnew = h1;
hnew.Add(&h2,-1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1S operator*(TH1S &h1, TH1S &h2)
{
TH1S hnew = h1;
hnew.Multiply(&h2);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1S operator/(TH1S &h1, TH1S &h2)
{
TH1S hnew = h1;
hnew.Divide(&h2);
hnew.SetDirectory(0);
return hnew;
}
ClassImp(TH1F)
//______________________________________________________________________________
// TH1F methods
//______________________________________________________________________________
TH1F::TH1F(): TH1(), TArrayF()
{
fDimension = 1;
}
//______________________________________________________________________________
TH1F::TH1F(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
: TH1(name,title,nbins,xlow,xup), TArrayF()
{
fDimension = dim;
}
//______________________________________________________________________________
TH1F::TH1F(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
: TH1(name,title,nbins,xbins), TArrayF()
{
fDimension = dim;
}
//______________________________________________________________________________
TH1F::TH1F(const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
: TH1(name,title,nbins,xlow,xup), TArrayF(nbins+2)
{
//
// Create a 1-Dim histogram with fix bins of type float
// ====================================================
// (see TH1::TH1 for explanation of parameters)
//
fDimension = 1;
}
//______________________________________________________________________________
TH1F::TH1F(const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
: TH1(name,title,nbins,xbins), TArrayF(nbins+2)
{
//
// Create a 1-Dim histogram with variable bins of type float
// =========================================================
// (see TH1::TH1 for explanation of parameters)
//
fDimension = 1;
}
//______________________________________________________________________________
TH1F::TH1F(const TH1F &h)
{
((TH1F&)h).Copy(*this);
}
//______________________________________________________________________________
TH1F::~TH1F()
{
}
//______________________________________________________________________________
void TH1F::Copy(TObject &newth1)
{
TH1::Copy(newth1);
TArrayF::Copy((TH1F&)newth1);
}
//______________________________________________________________________________
TH1 *TH1F::DrawCopy(Option_t *option)
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
TH1F *newth1 = (TH1F*)Clone();
newth1->SetDirectory(0);
newth1->SetBit(kCanDelete);
newth1->AppendPad(option);
return newth1;
}
//______________________________________________________________________________
Stat_t TH1F::GetBinContent(Int_t bin)
{
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
return Stat_t (fArray[bin]);
}
//______________________________________________________________________________
void TH1F::Reset(Option_t *option)
{
TH1::Reset(option);
TArrayF::Reset();
}
//______________________________________________________________________________
TH1F& TH1F::operator=(const TH1F &h1)
{
if (this != &h1) ((TH1F&)h1).Copy(*this);
return *this;
}
//______________________________________________________________________________
TH1F operator*(Float_t c1, TH1F &h1)
{
TH1F hnew = h1;
hnew.Scale(c1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1F operator+(TH1F &h1, TH1F &h2)
{
TH1F hnew = h1;
hnew.Add(&h2,1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1F operator-(TH1F &h1, TH1F &h2)
{
TH1F hnew = h1;
hnew.Add(&h2,-1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1F operator*(TH1F &h1, TH1F &h2)
{
TH1F hnew = h1;
hnew.Multiply(&h2);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1F operator/(TH1F &h1, TH1F &h2)
{
TH1F hnew = h1;
hnew.Divide(&h2);
hnew.SetDirectory(0);
return hnew;
}
ClassImp(TH1D)
//______________________________________________________________________________
// TH1D methods
//______________________________________________________________________________
TH1D::TH1D(): TH1(), TArrayD()
{
fDimension = 1;
}
//______________________________________________________________________________
TH1D::TH1D(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
: TH1(name,title,nbins,xlow,xup), TArrayD()
{
fDimension = dim;
}
//______________________________________________________________________________
TH1D::TH1D(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
: TH1(name,title,nbins,xbins), TArrayD()
{
fDimension = dim;
}
//______________________________________________________________________________
TH1D::TH1D(const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
: TH1(name,title,nbins,xlow,xup), TArrayD(nbins+2)
{
//
// Create a 1-Dim histogram with fix bins of type double
// =====================================================
// (see TH1::TH1 for explanation of parameters)
//
fDimension = 1;
}
//______________________________________________________________________________
TH1D::TH1D(const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
: TH1(name,title,nbins,xbins), TArrayD(nbins+2)
{
//
// Create a 1-Dim histogram with variable bins of type double
// =====================================================
// (see TH1::TH1 for explanation of parameters)
//
fDimension = 1;
}
//______________________________________________________________________________
TH1D::~TH1D()
{
}
//______________________________________________________________________________
TH1D::TH1D(const TH1D &h1d)
{
((TH1D&)h1d).Copy(*this);
}
//______________________________________________________________________________
void TH1D::Copy(TObject &newth1)
{
TH1::Copy(newth1);
TArrayD::Copy((TH1D&)newth1);
}
//______________________________________________________________________________
TH1 *TH1D::DrawCopy(Option_t *option)
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
TH1D *newth1 = (TH1D*)Clone();
newth1->SetDirectory(0);
newth1->SetBit(kCanDelete);
newth1->AppendPad(option);
return newth1;
}
//______________________________________________________________________________
Stat_t TH1D::GetBinContent(Int_t bin)
{
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
return Stat_t (fArray[bin]);
}
//______________________________________________________________________________
void TH1D::Reset(Option_t *option)
{
TH1::Reset(option);
TArrayD::Reset();
}
//______________________________________________________________________________
TH1D& TH1D::operator=(const TH1D &h1)
{
if (this != &h1) ((TH1D&)h1).Copy(*this);
return *this;
}
//______________________________________________________________________________
TH1D operator*(Float_t c1, TH1D &h1)
{
TH1D hnew = h1;
hnew.Scale(c1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1D operator+(TH1D &h1, TH1D &h2)
{
TH1D hnew = h1;
hnew.Add(&h2,1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1D operator-(TH1D &h1, TH1D &h2)
{
TH1D hnew = h1;
hnew.Add(&h2,-1);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1D operator*(TH1D &h1, TH1D &h2)
{
TH1D hnew = h1;
hnew.Multiply(&h2);
hnew.SetDirectory(0);
return hnew;
}
//______________________________________________________________________________
TH1D operator/(TH1D &h1, TH1D &h2)
{
TH1D hnew = h1;
hnew.Divide(&h2);
hnew.SetDirectory(0);
return hnew;
}
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.