From 1a78b9cebf20130f317821a4be4be8f54e95e011 Mon Sep 17 00:00:00 2001
From: Lorenzo Moneta <Lorenzo.Moneta@cern.ch>
Date: Wed, 7 May 2014 16:36:27 +0200
Subject: [PATCH] Implement correctly the projections from a TProfile2D to
TProfile and from a TProfile3D to a TProfile2D. Add corresponding tests in
stressHistogram. This fixes bug ROOT-6220
---
hist/hist/inc/TH3.h | 23 ++--
hist/hist/inc/TProfile2D.h | 4 +
hist/hist/inc/TProfile3D.h | 7 +-
hist/hist/src/TH2.cxx | 2 +-
hist/hist/src/TProfile.cxx | 16 ++-
hist/hist/src/TProfile2D.cxx | 152 ++++++++++++++++++++++----
hist/hist/src/TProfile3D.cxx | 190 +++++++++++++++++++++++++++++++--
test/stressHistogram.cxx | 200 ++++++++++++++++++++++++++++++++---
8 files changed, 530 insertions(+), 64 deletions(-)
diff --git a/hist/hist/inc/TH3.h b/hist/hist/inc/TH3.h
index 9dfaa629e8c..266ad5a46c8 100644
--- a/hist/hist/inc/TH3.h
+++ b/hist/hist/inc/TH3.h
@@ -121,7 +121,7 @@ public:
TH1D *ProjectionZ(const char *name="_pz", Int_t ixmin=0, Int_t ixmax=-1, Int_t iymin=0,
Int_t iymax=-1, Option_t *option="") const; // *MENU*
TH1 *Project3D(Option_t *option="x") const; // *MENU*
- TProfile2D *Project3DProfile(Option_t *option="xy") const; // *MENU*
+ virtual TProfile2D *Project3DProfile(Option_t *option="xy") const; // *MENU*
virtual void PutStats(Double_t *stats);
virtual TH3 *RebinX(Int_t ngroup = 2, const char *newname = "");
virtual TH3 *RebinY(Int_t ngroup = 2, const char *newname = "");
@@ -134,15 +134,26 @@ public:
virtual void SetShowProjection(const char *option="xy",Int_t nbins=1); // *MENU*
protected:
- TH1D *DoProject1D(const char* name, const char * title, TAxis* projX,
+
+ virtual TH1D *DoProject1D(const char* name, const char * title, TAxis* projX,
bool computeErrors, bool originalRange,
bool useUF, bool useOF) const;
- TH2D *DoProject2D(const char* name, const char * title, TAxis* projX, TAxis* projY,
+ virtual TH2D *DoProject2D(const char* name, const char * title, TAxis* projX, TAxis* projY,
bool computeErrors, bool originalRange,
bool useUF, bool useOF) const;
- TProfile2D *DoProjectProfile2D(const char* name, const char * title, TAxis* projX, TAxis* projY,
- bool originalRange, bool useUF, bool useOF) const;
-
+ virtual TProfile2D *DoProjectProfile2D(const char* name, const char * title, TAxis* projX, TAxis* projY,
+ bool originalRange, bool useUF, bool useOF) const;
+
+ // these functions are need to be used inside TProfile3D::DoProjectProfile2D
+ static TH1D *DoProject1D(const TH3 & h, const char* name, const char * title, TAxis* projX,
+ bool computeErrors, bool originalRange, bool useUF, bool useOF) {
+ return h.DoProject1D(name, title, projX, computeErrors, originalRange, useUF, useOF);
+ }
+ static TH2D *DoProject2D(const TH3 & h, const char* name, const char * title, TAxis* projX, TAxis* projY,
+ bool computeErrors, bool originalRange, bool useUF, bool useOF) {
+ return h.DoProject2D(name, title, projX,projY, computeErrors, originalRange, useUF, useOF);
+ }
+
ClassDef(TH3,5) //3-Dim histogram base class
};
diff --git a/hist/hist/inc/TProfile2D.h b/hist/hist/inc/TProfile2D.h
index d7aa86ade7a..8d680c11ceb 100644
--- a/hist/hist/inc/TProfile2D.h
+++ b/hist/hist/inc/TProfile2D.h
@@ -53,6 +53,8 @@ protected:
nbins[1], range[2], range[3]); };
Int_t Fill(const Double_t* v) { return Fill(v[0], v[1], v[2], v[3]); };
+ virtual TProfile *DoProfile(bool onX, const char *name, Int_t firstbin, Int_t lastbin, Option_t *option) const;
+
using TH2::Fill;
Int_t Fill(Double_t, Double_t) {return TH2::Fill(0); } //MayNotUse
@@ -127,6 +129,8 @@ public:
virtual Bool_t Multiply(const TH1 *h1);
virtual Bool_t Multiply(const TH1 *h1, const TH1 *h2, Double_t c1=1, Double_t c2=1, Option_t *option=""); // *MENU*
TH2D *ProjectionXY(const char *name="_pxy", Option_t *option="e") const;
+ TProfile *ProfileX(const char *name="_pfx", Int_t firstybin=0, Int_t lastybin=-1, Option_t *option="") const; // *MENU*
+ TProfile *ProfileY(const char *name="_pfy", Int_t firstxbin=0, Int_t lastxbin=-1, Option_t *option="") const; // *MENU*
virtual void PutStats(Double_t *stats);
virtual void Reset(Option_t *option="");
virtual TProfile2D *Rebin2D(Int_t nxgroup=2, Int_t nygroup=2, const char *newname="");
diff --git a/hist/hist/inc/TProfile3D.h b/hist/hist/inc/TProfile3D.h
index 926a5a5d231..86031f3e1d9 100644
--- a/hist/hist/inc/TProfile3D.h
+++ b/hist/hist/inc/TProfile3D.h
@@ -69,7 +69,9 @@ protected:
//virtual void UpdateBinContent(Int_t bin, Double_t content);
virtual Double_t GetBinErrorSqUnchecked(Int_t bin) const { Double_t err = GetBinError(bin); return err*err; }
-
+ virtual TProfile2D *DoProjectProfile2D(const char* name, const char * title, TAxis* projX, TAxis* projY,
+ bool originalRange, bool useUF, bool useOF) const;
+
private:
Double_t *GetB() {return &fBinEntries.fArray[0];}
Double_t *GetB2() {return (fBinSumw2.fN ? &fBinSumw2.fArray[0] : 0 ); }
@@ -127,7 +129,8 @@ public:
virtual Bool_t Multiply(TF1 *h1, Double_t c1=1);
virtual Bool_t Multiply(const TH1 *h1);
virtual Bool_t Multiply(const TH1 *h1, const TH1 *h2, Double_t c1=1, Double_t c2=1, Option_t *option=""); // *MENU*
- TH3D *ProjectionXYZ(const char *name="_pxyz", Option_t *option="e") const;
+ virtual TH3D *ProjectionXYZ(const char *name="_pxyz", Option_t *option="e") const;
+ virtual TProfile2D *Project3DProfile(Option_t *option="xy") const; // *MENU*
virtual void PutStats(Double_t *stats);
virtual void Reset(Option_t *option="");
virtual void SavePrimitive(std::ostream &out, Option_t *option = "");
diff --git a/hist/hist/src/TH2.cxx b/hist/hist/src/TH2.cxx
index 721123c438b..8aa096b9673 100644
--- a/hist/hist/src/TH2.cxx
+++ b/hist/hist/src/TH2.cxx
@@ -2123,7 +2123,7 @@ TProfile *TH2::ProfileX(const char *name, Int_t firstybin, Int_t lastybin, Optio
//
// if option "d" is specified, the profile is drawn in the current pad.
//
- // if option "o" original axis range of the taget axes will be
+ // if option "o" original axis range of the target axes will be
// kept, but only bins inside the selected range will be filled.
//
// The option can also be used to specify the projected profile error type.
diff --git a/hist/hist/src/TProfile.cxx b/hist/hist/src/TProfile.cxx
index c13761dc237..950dd8a84b0 100644
--- a/hist/hist/src/TProfile.cxx
+++ b/hist/hist/src/TProfile.cxx
@@ -1238,7 +1238,7 @@ TH1D *TProfile::ProjectionX(const char *name, Option_t *option) const
if (opt.Contains("e")) computeErrors = kTRUE;
if (opt.Contains("w")) binWeight = kTRUE;
if (opt.Contains("c=e")) {cequalErrors = kTRUE; computeErrors=kFALSE;}
- if (computeErrors || binWeight ) h1->Sumw2();
+ if (computeErrors || binWeight || (binEntries && fBinSumw2.fN) ) h1->Sumw2();
// Fill the projected histogram
Double_t cont;
@@ -1255,15 +1255,11 @@ TH1D *TProfile::ProjectionX(const char *name, Option_t *option) const
if (computeErrors ) h1->SetBinError(bin , GetBinError(bin) );
// in case of option W bin error is deduced from bin sum of z**2 values of profile
// this is correct only if the profile is filled with weights =1
- if (binWeight) h1->SetBinError(bin , TMath::Sqrt(fSumw2.fArray[bin] ) );
- // in case of bin entries and h1 has sumw2 set, we need to set also the bin error
- if (binEntries && h1->GetSumw2() ) {
- Double_t err2;
- if (fBinSumw2.fN)
- err2 = fBinSumw2.fArray[bin];
- else
- err2 = cont; // this is fBinEntries.fArray[bin]
- h1->SetBinError(bin, TMath::Sqrt(err2 ) );
+ if (binWeight) h1->GetSumw2()->fArray[bin] = fSumw2.fArray[bin];
+ // in case of bin entries and profile is weighted, we need to set also the bin error
+ if (binEntries && fBinSumw2.fN ) {
+ R__ASSERT( h1->GetSumw2() );
+ h1->GetSumw2()->fArray[bin] = fBinSumw2.fArray[bin];
}
}
diff --git a/hist/hist/src/TProfile2D.cxx b/hist/hist/src/TProfile2D.cxx
index f26fe2c11e5..0d9b9569494 100644
--- a/hist/hist/src/TProfile2D.cxx
+++ b/hist/hist/src/TProfile2D.cxx
@@ -1192,21 +1192,21 @@ TH2D *TProfile2D::ProjectionXY(const char *name, Option_t *option) const
// This option makes sense only for profile filled with all weights =1.
// When the profile is weighted (filled with weights different than 1) the
// bin error of the projected histogram (obtained using this option "W") cannot be
-// correctly computed from the information stored in the profile.
+// correctly computed from the information stored in the profile. In that case the
+// obtained histogram contains as bin error square the weighted sum of the square of the
+// profiled observable (TProfile2D::fSumw2[bin] )
TString opt = option;
opt.ToLower();
+ // Create the projection histogram
+ // name of projected histogram is by default name of orginal histogram + _pxy
+ TString pname(name);
+ if (pname.IsNull() || pname == "_pxy")
+ pname = TString(GetName() ) + TString("_pxy");
- // Create the projection histogram
- char *pname = (char*)name;
- if (strcmp(name,"_px") == 0) {
- Int_t nch = strlen(GetName()) + 4;
- pname = new char[nch];
- snprintf(pname,nch,"%s%s",GetName(),name);
- }
Int_t nx = fXaxis.GetNbins();
Int_t ny = fYaxis.GetNbins();
const TArrayD *xbins = fXaxis.GetXbins();
@@ -1229,8 +1229,7 @@ TH2D *TProfile2D::ProjectionXY(const char *name, Option_t *option) const
if (opt.Contains("e")) computeErrors = kTRUE;
if (opt.Contains("w")) binWeight = kTRUE;
if (opt.Contains("c=e")) {cequalErrors = kTRUE; computeErrors=kFALSE;}
- if (computeErrors || binWeight ) h1->Sumw2();
- if (pname != name) delete [] pname;
+ if (computeErrors || binWeight || (binEntries && fBinSumw2.fN) ) h1->Sumw2();
// Fill the projected histogram
Int_t bin,binx, biny;
@@ -1250,15 +1249,11 @@ TH2D *TProfile2D::ProjectionXY(const char *name, Option_t *option) const
if (computeErrors ) h1->SetBinError(bin , GetBinError(bin) );
// in case of option W bin error is deduced from bin sum of z**2 values of profile
// this is correct only if the profile is unweighted
- if (binWeight) h1->SetBinError(bin , TMath::Sqrt(fSumw2.fArray[bin] ) );
- // in case of bin entries and h1 has sumw2 set set, we need to set the also bin error
- if (binEntries && h1->GetSumw2() ) {
- Double_t err2;
- if (fBinSumw2.fN)
- err2 = fBinSumw2.fArray[bin];
- else
- err2 = cont; // this is fBinEntries.fArray[bin]
- h1->SetBinError(bin, TMath::Sqrt(err2 ) );
+ if (binWeight) h1->GetSumw2()->fArray[bin] = fSumw2.fArray[bin];
+ // in case of bin entries and profile is weighted, we need to set also the bin error
+ if (binEntries && fBinSumw2.fN ) {
+ R__ASSERT( h1->GetSumw2() );
+ h1->GetSumw2()->fArray[bin] = fBinSumw2.fArray[bin];
}
}
}
@@ -1266,6 +1261,125 @@ TH2D *TProfile2D::ProjectionXY(const char *name, Option_t *option) const
return h1;
}
+//______________________________________________________________________________
+TProfile *TProfile2D::ProfileX(const char *name, Int_t firstybin, Int_t lastybin, Option_t *option) const
+{
+ // *-*-*-*-*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.
+ // The result is a 1D profile which contains the combination of all the considered bins along Y
+ // By default, bins 1 to ny are included
+ // When all bins are included, the number of entries in the projection
+ // is set to the number of entries of the 2-D histogram, otherwise
+ // the number of entries is incremented by 1 for all non empty cells.
+ //
+ // The option can also be used to specify the projected profile error type.
+ // Values which can be used are 's', 'i', or 'g'. See TProfile::BuildOptions for details
+ //
+ //
+ return DoProfile(true, name, firstybin, lastybin, option);
+}
+
+//______________________________________________________________________________
+TProfile *TProfile2D::ProfileY(const char *name, Int_t firstxbin, Int_t lastxbin, Option_t *option) const
+{
+ // *-*-*-*-*Project a 2-D histogram into a profile histogram along X*-*-*-*-*-*
+ // *-* ========================================================
+ //
+ // The projection is made from the channels along the X axis
+ // ranging from firstybin to lastybin included.
+ // The result is a 1D profile which contains the combination of all the considered bins along X
+ // By default, bins 1 to ny are included
+ // When all bins are included, the number of entries in the projection
+ // is set to the number of entries of the 2-D histogram, otherwise
+ // the number of entries is incremented by 1 for all non empty cells.
+ //
+ // The option can also be used to specify the projected profile error type.
+ // Values which can be used are 's', 'i', or 'g'. See TProfile::BuildOptions for details
+ //
+ //
+ //
+ return DoProfile(false, name, firstxbin, lastxbin, option);
+}
+
+//______________________________________________________________________________
+TProfile * TProfile2D::DoProfile(bool onX, const char *name, Int_t firstbin, Int_t lastbin, Option_t *option) const {
+ // implementation of ProfileX or ProfileY for a TProfile2D
+ // Do correctly the combination of the bin averages when doing the projection
+
+ TString opt = option;
+ opt.ToLower();
+ bool originalRange = opt.Contains("o");
+
+ TString expectedName = ( onX ? "_pfx" : "_pfy" );
+
+ TString pname(name);
+ if (pname.IsNull() || name == expectedName)
+ pname = TString(GetName() ) + expectedName;
+
+ const TAxis& outAxis = ( onX ? fXaxis : fYaxis );
+ const TArrayD *bins = outAxis.GetXbins();
+ Int_t firstOutBin = outAxis.GetFirst();
+ Int_t lastOutBin = outAxis.GetLast();
+
+ TProfile * p1 = 0;
+ // case of fixed bins
+ if (bins->fN == 0) {
+ if (originalRange)
+ p1 = new TProfile(pname,GetTitle(), outAxis.GetNbins(), outAxis.GetXmin(), outAxis.GetXmax(), opt );
+ else
+ p1 = new TProfile(pname,GetTitle(), lastOutBin-firstOutBin+1,
+ outAxis.GetBinLowEdge(firstOutBin),outAxis.GetBinUpEdge(lastOutBin), opt);
+ } else {
+ // case of variable bins
+ if (originalRange )
+ p1 = new TProfile(pname,GetTitle(),outAxis.GetNbins(),bins->fArray,opt);
+ else
+ p1 = new TProfile(pname,GetTitle(),lastOutBin-firstOutBin+1,&bins->fArray[firstOutBin-1],opt);
+
+ }
+
+ if (fBinSumw2.fN) p1->Sumw2();
+
+ // make projection in a 2D first
+ TH2D * h2dW = ProjectionXY("h2temp-W","W");
+ TH2D * h2dN = ProjectionXY("h2temp-N","B");
+
+ h2dW->SetDirectory(0); h2dN->SetDirectory(0);
+
+
+ TString opt1 = (originalRange) ? "o" : "";
+ TH1D * h1W = (onX) ? h2dW->ProjectionX("h1temp-W",firstbin,lastbin,opt1) : h2dW->ProjectionY("h1temp-W",firstbin,lastbin,opt1);
+ TH1D * h1N = (onX) ? h2dN->ProjectionX("h1temp-N",firstbin,lastbin,opt1) : h2dN->ProjectionY("h1temp-N",firstbin,lastbin,opt1);
+ h1W->SetDirectory(0); h1N->SetDirectory(0);
+
+
+ // fill the bin content
+ R__ASSERT( h1W->fN == p1->fN );
+ R__ASSERT( h1N->fN == p1->fN );
+ R__ASSERT( h1W->GetSumw2()->fN != 0); // h1W should always be a weighted histogram since h2dW is
+ for (int i = 0; i < p1->fN ; ++i) {
+ p1->fArray[i] = h1W->GetBinContent(i); // array of profile is sum of all values
+ p1->GetSumw2()->fArray[i] = h1W->GetSumw2()->fArray[i]; // array of content square of profile is weight square of the W projected histogram
+ p1->SetBinEntries(i, h1N->GetBinContent(i) );
+ if (fBinSumw2.fN) p1->GetBinSumw2()->fArray[i] = h1N->GetSumw2()->fArray[i]; // sum of weight squares are stored to compute errors in h1N histogram
+ }
+ // delete the created histograms
+ delete h2dW;
+ delete h2dN;
+ delete h1W;
+ delete h1N;
+
+ // Also we need to set the entries since they have not been correctly calculated during the projection
+ // we can only set them to the effective entries
+ p1->SetEntries( p1->GetEffectiveEntries() );
+
+ return p1;
+}
+
+
//______________________________________________________________________________
void TProfile2D::PutStats(Double_t *stats)
{
diff --git a/hist/hist/src/TProfile3D.cxx b/hist/hist/src/TProfile3D.cxx
index aee1a4fb901..4085feb4379 100644
--- a/hist/hist/src/TProfile3D.cxx
+++ b/hist/hist/src/TProfile3D.cxx
@@ -10,6 +10,7 @@
*************************************************************************/
#include "TProfile3D.h"
+#include "TProfile2D.h"
#include "THashList.h"
#include "TMath.h"
#include "THLimitsFinder.h"
@@ -1064,47 +1065,214 @@ TH3D *TProfile3D::ProjectionXYZ(const char *name, Option_t *option) const
// will be equal to the GetBinEntries(bin) of the profile,
// if option "C=E" the bin contents of the projection are set to the
// bin errors of the profile
-//
+// if option "E" is specified the errors of the projected histogram are computed and set
+// to be equal to the errors of the profile.
+// Option "E" is defined as the default one in the header file.
+// if option "" is specified the histogram errors are simply the sqrt of its content
+// if option "B" is specified, the content of bin of the returned histogram
+// will be equal to the GetBinEntries(bin) of the profile,
+// if option "C=E" the bin contents of the projection are set to the
+// bin errors of the profile
+// if option "W" is specified the bin content of the projected histogram is set to the
+// product of the bin content of the profile and the entries.
+// With this option the returned histogram will be equivalent to the one obtained by
+// filling directly a TH2D using the 3-rd value as a weight.
+// This option makes sense only for profile filled with all weights =1.
+// When the profile is weighted (filled with weights different than 1) the
+// bin error of the projected histogram (obtained using this option "W") cannot be
+// correctly computed from the information stored in the profile. In that case the
+// obtained histogram contains as bin error square the weighted sum of the square of the
+// profiled observable (TProfile2D::fSumw2[bin] )
+
TString opt = option;
opt.ToLower();
Int_t nx = fXaxis.GetNbins();
Int_t ny = fYaxis.GetNbins();
Int_t nz = fZaxis.GetNbins();
+ const TArrayD *xbins = fXaxis.GetXbins();
+ const TArrayD *ybins = fYaxis.GetXbins();
+ const TArrayD *zbins = fZaxis.GetXbins();
// Create the projection histogram
TString pname = name;
if (pname == "_px") {
- pname = GetName(); pname.Append("_px");
+ pname = GetName(); pname.Append("_pxyz");
+ }
+ TH3D *h1 = 0 ;
+ if (xbins->fN == 0 && ybins->fN == 0 && zbins->fN == 0)
+ h1 = new TH3D(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax(),ny,fYaxis.GetXmin(),fYaxis.GetXmax(),nz,fZaxis.GetXmin(),fZaxis.GetXmax());
+ else if ( xbins->fN != 0 && ybins->fN != 0 && zbins->fN != 0)
+ h1 = new TH3D(pname,GetTitle(),nx,xbins->GetArray(),ny,ybins->GetArray(), nz,zbins->GetArray() );
+ else {
+ Error("ProjectionXYZ","Histogram has an axis with variable bins and an axis with fixed bins. This case is not cupported - return a null pointer");
+ return 0;
}
- TH3D *h1 = new TH3D(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax(),ny,fYaxis.GetXmin(),fYaxis.GetXmax(),nz,fZaxis.GetXmin(),fZaxis.GetXmax());
+
+
Bool_t computeErrors = kFALSE;
Bool_t cequalErrors = kFALSE;
Bool_t binEntries = kFALSE;
+ Bool_t binWeight = kFALSE;
+
if (opt.Contains("b")) binEntries = kTRUE;
if (opt.Contains("e")) computeErrors = kTRUE;
+ if (opt.Contains("w")) binWeight = kTRUE;
if (opt.Contains("c=e")) {cequalErrors = kTRUE; computeErrors=kFALSE;}
- if (computeErrors) h1->Sumw2();
+ if (computeErrors || binWeight || (binEntries && fBinSumw2.fN) ) h1->Sumw2();
// Fill the projected histogram
Int_t bin,binx,biny,binz;
- Double_t cont,err;
+ Double_t cont;
for (binx =0;binx<=nx+1;binx++) {
for (biny =0;biny<=ny+1;biny++) {
for (binz =0;binz<=nz+1;binz++) {
bin = GetBin(binx,biny,binz);
- if (binEntries) cont = GetBinEntries(bin);
- else cont = GetBinContent(bin);
- err = GetBinError(bin);
- if (cequalErrors) h1->SetBinContent(binx,biny,binz, err);
- else h1->SetBinContent(binx,biny,binz,cont);
- if (computeErrors) h1->SetBinError(binx,biny,binz,err);
+
+ if (binEntries) cont = GetBinEntries(bin);
+ else if (cequalErrors) cont = GetBinError(bin);
+ else if (binWeight) cont = GetBinContent(bin) * GetBinEntries(bin);
+ else cont = GetBinContent(bin); // default case
+
+ h1->SetBinContent(bin ,cont);
+
+ // if option E projected histogram errors are same as profile
+ if (computeErrors ) h1->SetBinError(bin , GetBinError(bin) );
+ // in case of option W bin error is deduced from bin sum of z**2 values of profile
+ // this is correct only if the profile is unweighted
+ if (binWeight) h1->GetSumw2()->fArray[bin] = fSumw2.fArray[bin];
+ // in case of bin entries and profile is weighted, we need to set also the bin error
+ if (binEntries && fBinSumw2.fN ) {
+ R__ASSERT( h1->GetSumw2() );
+ h1->GetSumw2()->fArray[bin] = fBinSumw2.fArray[bin];
+ }
}
}
}
h1->SetEntries(fEntries);
return h1;
}
+//______________________________________________________________________________
+TProfile2D *TProfile3D::Project3DProfile(Option_t *option) const
+{
+ // *-*-*-*-*Project a 3-D profile into a 2D-profile histogram depending
+ // on the option parameter
+ // option may contain a combination of the characters x,y,z
+ // option = "xy" return the x versus y projection into a TProfile2D histogram
+ // option = "yx" return the y versus x projection into a TProfile2D histogram
+ // option = "xz" return the x versus z projection into a TProfile2D histogram
+ // option = "zx" return the z versus x projection into a TProfile2D histogram
+ // option = "yz" return the y versus z projection into a TProfile2D histogram
+ // option = "zy" return the z versus y projection into a TProfile2D histogram
+ // NB: the notation "a vs b" means "a" vertical and "b" horizontalalong X*-*-*-*-*-*
+ //
+ // The resulting profile contains the combination of all the considered bins along X
+ // By default, all bins are included considering also underflow/overflows
+ //
+ // The option can also be used to specify the projected profile error type.
+ // Values which can be used are 's', 'i', or 'g'. See TProfile::BuildOptions for details
+ //
+ // To select a bin range along an axis, use TAxis::SetRange, eg
+ // h3.GetYaxis()->SetRange(23,56);
+ //
+ //
+
+ // can call TH3 method which will call the virtual method :DoProjectProfile2D reimplented below
+ // but need to add underflow/overflow
+ TString opt(option);
+ opt.Append(" UF OF");
+ return TH3::Project3DProfile(opt);
+}
+
+//______________________________________________________________________________
+TProfile2D *TProfile3D::DoProjectProfile2D(const char* name, const char * title, TAxis* projX, TAxis* projY,
+ bool originalRange, bool useUF, bool useOF) const
+{
+ // internal method to project to a 2D Profile
+ // called from TH3::Project3DProfile but re-implemented in case of the TPRofile3D since what is done is different
+
+ // Get the ranges where we will work.
+ Int_t ixmin = projX->GetFirst();
+ Int_t ixmax = projX->GetLast();
+ Int_t iymin = projY->GetFirst();
+ Int_t iymax = projY->GetLast();
+ if (ixmin == 0 && ixmax == 0) { ixmin = 1; ixmax = projX->GetNbins(); }
+ if (iymin == 0 && iymax == 0) { iymin = 1; iymax = projY->GetNbins(); }
+ Int_t nx = ixmax-ixmin+1;
+ Int_t ny = iymax-iymin+1;
+
+ // Create the projected profiles
+ TProfile2D *p2 = 0;
+ // Create always a new TProfile2D (not as in the case of TH3 projection)
+
+ const TArrayD *xbins = projX->GetXbins();
+ const TArrayD *ybins = projY->GetXbins();
+ // assume all axis have variable bins or have fixed bins
+ if ( originalRange ) {
+ if (xbins->fN == 0 && ybins->fN == 0) {
+ p2 = new TProfile2D(name,title,projY->GetNbins(),projY->GetXmin(),projY->GetXmax()
+ ,projX->GetNbins(),projX->GetXmin(),projX->GetXmax());
+ } else {
+ p2 = new TProfile2D(name,title,projY->GetNbins(),&ybins->fArray[iymin-1],projX->GetNbins(),&xbins->fArray[ixmin-1]);
+ }
+ } else {
+ if (xbins->fN == 0 && ybins->fN == 0) {
+ p2 = new TProfile2D(name,title,ny,projY->GetBinLowEdge(iymin),projY->GetBinUpEdge(iymax)
+ ,nx,projX->GetBinLowEdge(ixmin),projX->GetBinUpEdge(ixmax));
+ } else {
+ p2 = new TProfile2D(name,title,ny,&ybins->fArray[iymin-1],nx,&xbins->fArray[ixmin-1]);
+ }
+ }
+
+ // weights
+ bool useWeights = (fBinSumw2.fN != 0);
+ if (useWeights) p2->Sumw2();
+
+ // make projection in a 3D first
+ TH3D * h3dW = ProjectionXYZ("h3temp-W","W");
+ TH3D * h3dN = ProjectionXYZ("h3temp-N","B");
+
+ h3dW->SetDirectory(0); h3dN->SetDirectory(0);
+
+ // note that h3dW is always a weighted histogram - so we need to compute error in the projection
+ TAxis * projX_hW = h3dW->GetXaxis();
+ TAxis * projX_hN = h3dN->GetXaxis();
+ if (projX == GetYaxis() ) { projX_hW = h3dW->GetYaxis(); projX_hN = h3dN->GetYaxis(); }
+ if (projX == GetZaxis() ) { projX_hW = h3dW->GetZaxis(); projX_hN = h3dN->GetZaxis(); }
+ TAxis * projY_hW = h3dW->GetYaxis();
+ TAxis * projY_hN = h3dN->GetYaxis();
+ if (projY == GetXaxis() ) { projY_hW = h3dW->GetXaxis(); projY_hN = h3dN->GetXaxis(); }
+ if (projY == GetZaxis() ) { projY_hW = h3dW->GetZaxis(); projY_hN = h3dN->GetZaxis(); }
+
+ TH2D * h2W = TH3::DoProject2D(*h3dW,"htemp-W","",projX_hW, projY_hW, true, originalRange, useUF, useOF);
+ TH2D * h2N = TH3::DoProject2D(*h3dN,"htemp-N","",projX_hN, projY_hN, useWeights, originalRange, useUF, useOF);
+ h2W->SetDirectory(0); h2N->SetDirectory(0);
+
+
+ // fill the bin content
+ R__ASSERT( h2W->fN == p2->fN );
+ R__ASSERT( h2N->fN == p2->fN );
+ R__ASSERT( h2W->GetSumw2()->fN != 0); // h2W should always be a weighted histogram since h3dW is weighted
+ for (int i = 0; i < p2->fN ; ++i) {
+ //std::cout << " proj bin " << i << " " << h2W->fArray[i] << " " << h2N->fArray[i] << std::endl;
+ p2->fArray[i] = h2W->fArray[i]; // array of profile is sum of all values
+ p2->GetSumw2()->fArray[i] = h2W->GetSumw2()->fArray[i]; // array of content square of profile is weight square of the W projected histogram
+ p2->SetBinEntries(i, h2N->fArray[i] );
+ if (useWeights) p2->GetBinSumw2()->fArray[i] = h2N->GetSumw2()->fArray[i]; // sum of weight squares are stored to compute errors in h1N histogram
+ }
+ // delete the created histograms
+ delete h3dW;
+ delete h3dN;
+ delete h2W;
+ delete h2N;
+
+ // Also we need to set the entries since they have not been correctly calculated during the projection
+ // we can only set them to the effective entries
+ p2->SetEntries( p2->GetEffectiveEntries() );
+
+ return p2;
+
+}
//______________________________________________________________________________
void TProfile3D::PutStats(Double_t *stats)
diff --git a/test/stressHistogram.cxx b/test/stressHistogram.cxx
index b0d54baf14d..cd72fb34a82 100644
--- a/test/stressHistogram.cxx
+++ b/test/stressHistogram.cxx
@@ -8499,18 +8499,23 @@ bool testTH3toTH2()
static const double centre_deviation = 0.3;
-class ProjectionTester {
+struct ProjectionTester {
// This class implements the tests for all types of projections of
// all the classes tested in this file.
-private:
+//public:
static const unsigned int binsizeX = 8;
static const unsigned int binsizeY = 10;
static const unsigned int binsizeZ = 12;
static const int lower_limit = 0;
static const int upper_limit = 10;
-
+ static const int lower_limitX = 0;
+ static const int upper_limitX = 10;
+ static const int lower_limitY = -5;
+ static const int upper_limitY = 10;
+ static const int lower_limitZ = -10;
+ static const int upper_limitZ = 10;
TH3D* h3;
TH2D* h2XY;
@@ -8554,7 +8559,21 @@ private:
TH1D* hw1YZ;
TH1D* hw1ZX;
TH1D* hw1ZY;
-
+
+ TProfile3D* p3;
+
+ TProfile2D* p2XY;
+ TProfile2D* p2XZ;
+ TProfile2D* p2YX;
+ TProfile2D* p2YZ;
+ TProfile2D* p2ZX;
+ TProfile2D* p2ZY;
+
+ TProfile* p1X;
+ TProfile* p1Y;
+ TProfile* p1Z;
+
+
THnSparseD* s3;
THnD* n3;
@@ -8562,12 +8581,14 @@ private:
public:
- ProjectionTester()
+
+ ProjectionTester(bool useWeights = false)
{
+ buildWithWeights = useWeights;
+ CreateProfiles();
CreateHistograms();
- buildWithWeights = false;
}
-
+
void CreateHistograms()
{
h3 = new TH3D("h3","h3", binsizeX, lower_limit, upper_limit,
@@ -8609,7 +8630,7 @@ public:
pe2ZX = new TProfile2D("pe2ZX", "pe2ZX", binsizeZ, lower_limit, upper_limit,
binsizeX, lower_limit, upper_limit);
pe2ZY = new TProfile2D("pe2ZY", "pe2ZY", binsizeZ, lower_limit, upper_limit,
- binsizeY, lower_limit, upper_limit);
+ binsizeY, lower_limit, upper_limit);
h2wXY = new TH2D("h2wXY", "h2wXY", binsizeX, lower_limit, upper_limit,
binsizeY, lower_limit, upper_limit);
@@ -8659,6 +8680,32 @@ public:
n3 = new THnD("n3","n3", 3, bsize, xmin, xmax);
}
+
+ void CreateProfiles() {
+
+ // create Profile histograms
+ p3 = new TProfile3D("p3","p3", binsizeX, lower_limitX, upper_limitX,
+ binsizeY, lower_limitY, upper_limitY,
+ binsizeZ, lower_limitZ, upper_limitZ);
+
+ p2XY = new TProfile2D("p2XY", "p2XY", binsizeX, lower_limitX, upper_limitX,
+ binsizeY, lower_limitY, upper_limitY);
+ p2XZ = new TProfile2D("p2XZ", "p2XZ", binsizeX, lower_limitX, upper_limitX,
+ binsizeZ, lower_limitZ, upper_limitZ);
+ p2YX = new TProfile2D("p2YX", "p2YX", binsizeY, lower_limitY, upper_limitY,
+ binsizeX, lower_limitX, upper_limitX);
+ p2YZ = new TProfile2D("p2YZ", "p2YZ", binsizeY, lower_limitY, upper_limitY,
+ binsizeZ, lower_limitZ, upper_limitZ);
+ p2ZX = new TProfile2D("p2ZX", "p2ZX", binsizeZ, lower_limitZ, upper_limitZ,
+ binsizeX, lower_limitX, upper_limitX);
+ p2ZY = new TProfile2D("p2ZY", "p2ZY", binsizeZ, lower_limitZ, upper_limitZ,
+ binsizeY, lower_limitY, upper_limitY);
+
+ p1X = new TProfile("p1X", "pe1X", binsizeX, lower_limitX, upper_limitX);
+ p1Y = new TProfile("p1Y", "pe1Y", binsizeY, lower_limitY, upper_limitY);
+ p1Z = new TProfile("p1Z", "pe1Z", binsizeZ, lower_limitZ, upper_limitZ);
+
+ }
void DeleteHistograms()
{
@@ -8710,6 +8757,20 @@ public:
delete s3;
delete n3;
+ // profiles
+ delete p3;
+
+ delete p2XY;
+ delete p2XZ;
+ delete p2YX;
+ delete p2YZ;
+ delete p2ZX;
+ delete p2ZY;
+
+ delete p1X;
+ delete p1Y;
+ delete p1Z;
+
// delete all histogram in gROOT
TList * l = gROOT->GetList();
TIter next(l);
@@ -8727,6 +8788,7 @@ public:
void buildHistograms()
{
+
if (h3->GetSumw2N() ) {s3->Sumw2(); n3->Sumw2();}
for ( int ix = 0; ix <= h3->GetXaxis()->GetNbins() + 1; ++ix ) {
@@ -8807,6 +8869,7 @@ public:
void buildHistogramsWithWeights()
{
+
s3->Sumw2();
n3->Sumw2();
@@ -8994,6 +9057,7 @@ public:
buildWithWeights = false;
}
+
int compareHistograms()
{
int status = 0;
@@ -9262,7 +9326,97 @@ public:
return status;
}
-
+
+
+ void buildProfiles() {
+
+ if (buildWithWeights) {
+ p3->Sumw2();
+ p2XY->Sumw2(); p2YX->Sumw2(); p2YZ->Sumw2();
+ p2XZ->Sumw2(); p2ZX->Sumw2(); p2ZY->Sumw2();
+ p1X->Sumw2(); p1Y->Sumw2(); p1Z->Sumw2();
+ }
+
+
+ // use a different way to fill the histogram
+ for (int i = 0; i < 100000; ++i) {
+
+ // use in range in X but only overflow in Y and underflow/overflow in Z
+ double x = gRandom->Uniform(lower_limitX, upper_limitX );
+ double y = gRandom->Uniform(lower_limitY, upper_limitY+2.);
+ double z = gRandom->Uniform(lower_limitZ-1, upper_limitZ+1);
+ double u = TMath::Gaus(x,0,3)*TMath::Gaus(y,3,5)*TMath::Gaus(z,-3,10);
+
+ double w = 1;
+ if (buildWithWeights) w += x*x + (y-2)*(y-2) + (z+2)*(z+2);
+
+ p3->Fill(x,y,z,u,w);
+
+ p2XY->Fill(x,y,u,w);
+ p2YX->Fill(y,x,u,w);
+ p2XZ->Fill(x,z,u,w);
+ p2ZX->Fill(z,x,u,w);
+ p2YZ->Fill(y,z,u,w);
+ p2ZY->Fill(z,y,u,w);
+
+ p1X->Fill(x,u,w);
+ p1Y->Fill(y,u,w);
+ p1Z->Fill(z,u,w);
+
+ }
+
+ // reset the statistics to get same statistics computed from bin centers
+ p1X->ResetStats();
+ p1Y->ResetStats();
+ p1Z->ResetStats();
+
+ p2XY->ResetStats();
+ p2YX->ResetStats();
+ p2XZ->ResetStats();
+ p2ZX->ResetStats();
+ p2YZ->ResetStats();
+ p2ZY->ResetStats();
+ }
+
+
+ // actual test of profile projections
+ int compareProfiles()
+ {
+ int status = 0;
+ int options = 0;
+
+ // TProfile2d derived from TProfile3d
+ options = cmpOptStats;
+ //options = cmpOptPrint;
+ status += equals("TProfile3D -> XY", p2XY, p3->Project3DProfile("yx"), options);
+ status += equals("TProfile3D -> YX", p2YX, p3->Project3DProfile("xy"), options);
+ status += equals("TProfile3D -> XZ", p2XZ, p3->Project3DProfile("zx"), options);
+ status += equals("TProfile3D -> ZX", p2ZX, p3->Project3DProfile("xz"), options);
+ status += equals("TProfile3D -> YZ", p2YZ, p3->Project3DProfile("zy"), options);
+ status += equals("TProfile3D -> ZY", p2ZY, p3->Project3DProfile("yz"), options);
+ options = 0;
+ if ( defaultEqualOptions & cmpOptPrint )
+ cout << "----------------------------------------------" << endl;
+
+ // TProfile1 derived from TProfile2D from TProfile3D
+ options = cmpOptStats;
+ //options = cmpOptDebug;
+ TProfile2D* tmp1 = 0;
+ status += equals("TProfile2D -> X", p1X, p2XY->ProfileX(), options);
+ tmp1 = p3->Project3DProfile("xz");
+ status += equals("TProfile3D -> X", p1X, tmp1->ProfileY(), options);
+ delete tmp1; tmp1 = 0;
+ status += equals("TProfile2D -> Y", p1Y, p2ZY->ProfileY(), options);
+ tmp1 = p3->Project3DProfile("xy");
+ status += equals("TProfile3D -> X", p1Y, tmp1->ProfileX(), options);
+ delete tmp1; tmp1 = 0;
+ status += equals("TProfile2D -> Z", p1Z, p2ZX->ProfileX(), options);
+ tmp1 = p3->Project3DProfile("zy");
+ status += equals("TProfile3D -> Z", p1Z, tmp1->ProfileY(), options);
+ delete tmp1; tmp1 = 0;
+
+ return status;
+ }
};
int stressHistogram()
@@ -9310,10 +9464,18 @@ int stressHistogram()
//Ht->buildHistograms(2,4,5,6,8,10);
status = ht->compareHistograms();
GlobalStatus += status;
- printResult("Testing Projections without weights..............................", status);
delete ht;
+ printResult("Testing Histogram Projections without weights....................", status);
+
+ ProjectionTester* htp = new ProjectionTester();
+ htp->buildProfiles();
+ status = htp->compareProfiles();
+ GlobalStatus += status;
+ delete htp;
+
+ printResult("Testing Profile Projections without weights......................", status);
- // Test 2
+ // Test 3-4
if ( defaultEqualOptions & cmpOptPrint )
std::cout << "**********************************\n"
<< " Test with weights \n"
@@ -9324,8 +9486,16 @@ int stressHistogram()
ht2->buildHistogramsWithWeights();
status = ht2->compareHistograms();
GlobalStatus += status;
- printResult("Testing Projections with weights.................................", status);
+ printResult("Testing Histogram Projections with weights.......................", status);
delete ht2;
+
+
+ ProjectionTester* htp2 = new ProjectionTester(true);
+ htp2->buildProfiles();
+ status = htp2->compareProfiles();
+ GlobalStatus += status;
+ printResult("Testing Profile Projections with weights.......................", status);
+ delete htp2;
// Test 3
// Range Tests
@@ -9925,17 +10095,17 @@ int equals(const char* msg, TH1D* h1, TH1D* h2, int options, double ERRORLIMIT)
differents += (bool) equals(h1->GetXaxis()->GetNbins() , h2->GetXaxis()->GetNbins() );
if (debug) {
- cout << "Nbins = " << h1->GetXaxis()->GetNbins() << " , " << h2->GetXaxis()->GetNbins() << " | " << differents << std::endl;
+ cout << "Nbins = " << h1->GetXaxis()->GetNbins() << " | " << h2->GetXaxis()->GetNbins() << " | " << differents << std::endl;
}
differents += (bool) equals(h1->GetXaxis()->GetXmin() , h2->GetXaxis()->GetXmin() );
if (debug) {
- cout << "Xmin = " << h1->GetXaxis()->GetXmin() << " , " << h2->GetXaxis()->GetXmin() << " | " << differents << std::endl;
+ cout << "Xmin = " << h1->GetXaxis()->GetXmin() << " | " << h2->GetXaxis()->GetXmin() << " | " << differents << std::endl;
}
differents += (bool) equals(h1->GetXaxis()->GetXmax() , h2->GetXaxis()->GetXmax() );
if (debug) {
- cout << "Xmax = " << h1->GetXaxis()->GetXmax() << " , " << h2->GetXaxis()->GetXmax() << endl;
+ cout << "Xmax = " << h1->GetXaxis()->GetXmax() << " | " << h2->GetXaxis()->GetXmax() << endl;
}
for ( int i = 0; i <= h1->GetNbinsX() + 1; ++i )
--
GitLab