#include "TProfile.h"
#include "TMath.h"
#include "THashList.h"
#include "TF1.h"
#include "THLimitsFinder.h"
#include "Riostream.h"
#include "TVirtualPad.h"
#include "TError.h"
#include "TClass.h"
Bool_t TProfile::fgApproximate = kFALSE;
ClassImp(TProfile)
/*
<img src="gif/profile.gif">
*/
//End_Html
TProfile::TProfile() : TH1D()
{
BuildOptions(0,0,"");
}
TProfile::~TProfile()
{
}
TProfile::TProfile(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup,Option_t *option)
: TH1D(name,title,nbins,xlow,xup)
{
BuildOptions(0,0,option);
}
TProfile::TProfile(const char *name,const char *title,Int_t nbins,const Float_t *xbins,Option_t *option)
: TH1D(name,title,nbins,xbins)
{
BuildOptions(0,0,option);
}
TProfile::TProfile(const char *name,const char *title,Int_t nbins,const Double_t *xbins,Option_t *option)
: TH1D(name,title,nbins,xbins)
{
BuildOptions(0,0,option);
}
TProfile::TProfile(const char *name,const char *title,Int_t nbins,const Double_t *xbins,Double_t ylow,Double_t yup,Option_t *option)
: TH1D(name,title,nbins,xbins)
{
BuildOptions(ylow,yup,option);
}
TProfile::TProfile(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup,Double_t ylow,Double_t yup,Option_t *option)
: TH1D(name,title,nbins,xlow,xup)
{
BuildOptions(ylow,yup,option);
}
void TProfile::BuildOptions(Double_t ymin, Double_t ymax, Option_t *option)
{
SetErrorOption(option);
fBinEntries.Set(fNcells);
Sumw2();
fYmin = ymin;
fYmax = ymax;
fScaling = kFALSE;
fTsumwy = fTsumwy2 = 0;
}
TProfile::TProfile(const TProfile &profile) : TH1D()
{
((TProfile&)profile).Copy(*this);
}
void TProfile::Add(TF1 *, Double_t, Option_t * )
{
Error("Add","Function not implemented for TProfile");
return;
}
void TProfile::Add(const TH1 *h1, Double_t c1)
{
if (!h1) {
Error("Add","Attempt to add a non-existing profile");
return;
}
if (!h1->InheritsFrom("TProfile")) {
Error("Add","Attempt to add a non-profile object");
return;
}
TProfile *p1 = (TProfile*)h1;
Int_t nbinsx = GetNbinsX();
if (nbinsx != p1->GetNbinsX()) {
Error("Add","Attempt to add profiles with different number of bins");
return;
}
Double_t ac1 = TMath::Abs(c1);
fEntries += ac1*p1->GetEntries();
fTsumw += ac1*p1->fTsumw;
fTsumw2 += c1*c1*p1->fTsumw2;
fTsumwx += ac1*p1->fTsumwx;
fTsumwx2 += ac1*p1->fTsumwx2;
fTsumwy += ac1*p1->fTsumwy;
fTsumwy2 += ac1*p1->fTsumwy2;
Int_t bin;
Double_t *cu1 = p1->GetW();
Double_t *er1 = p1->GetW2();
Double_t *en1 = p1->GetB();
for (bin=0;bin<=nbinsx+1;bin++) {
fArray[bin] += c1*cu1[bin];
fSumw2.fArray[bin] += ac1*er1[bin];
fBinEntries.fArray[bin] += ac1*en1[bin];
}
}
void TProfile::Add(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2)
{
if (!h1 || !h2) {
Error("Add","Attempt to add a non-existing profile");
return;
}
if (!h1->InheritsFrom("TProfile")) {
Error("Add","Attempt to add a non-profile object");
return;
}
TProfile *p1 = (TProfile*)h1;
if (!h2->InheritsFrom("TProfile")) {
Error("Add","Attempt to add a non-profile object");
return;
}
TProfile *p2 = (TProfile*)h2;
Int_t nbinsx = GetNbinsX();
if (nbinsx != p1->GetNbinsX() || nbinsx != p2->GetNbinsX()) {
Error("Add","Attempt to add profiles with different number of bins");
return;
}
Double_t ac1 = TMath::Abs(c1);
Double_t ac2 = TMath::Abs(c2);
fEntries = ac1*p1->GetEntries() + ac2*p2->GetEntries();
fTsumw = ac1*p1->fTsumw + ac2*p2->fTsumw;
fTsumw2 = c1*c1*p1->fTsumw2 + c2*c2*p2->fTsumw2;
fTsumwx = ac1*p1->fTsumwx + ac2*p2->fTsumwx;
fTsumwx2 = ac1*p1->fTsumwx2 + ac2*p2->fTsumwx2;
fTsumwy = ac1*p1->fTsumwy + ac2*p2->fTsumwy;
fTsumwy2 = ac1*p1->fTsumwy2 + ac2*p2->fTsumwy2;
Int_t bin;
Double_t *cu1 = p1->GetW();
Double_t *cu2 = p2->GetW();
Double_t *er1 = p1->GetW2();
Double_t *er2 = p2->GetW2();
Double_t *en1 = p1->GetB();
Double_t *en2 = p2->GetB();
for (bin=0;bin<=nbinsx+1;bin++) {
fArray[bin] = c1*cu1[bin] + c2*cu2[bin];
if (fScaling) {
fSumw2.fArray[bin] = ac1*ac1*er1[bin] + ac2*ac2*er2[bin];
fBinEntries.fArray[bin] = en1[bin];
} else {
fSumw2.fArray[bin] = ac1*er1[bin] + ac2*er2[bin];
fBinEntries.fArray[bin] = ac1*en1[bin] + ac2*en2[bin];
}
}
}
void TProfile::Approximate(Bool_t approx)
{
fgApproximate = approx;
}
Int_t TProfile::BufferEmpty(Int_t action)
{
if (!fBuffer) return 0;
Int_t nbentries = (Int_t)fBuffer[0];
if (!nbentries) return 0;
Double_t *buffer = fBuffer;
if (nbentries < 0) {
if (action == 0) return 0;
nbentries = -nbentries;
fBuffer=0;
Reset();
fBuffer = buffer;
}
if (TestBit(kCanRebin) || fXaxis.GetXmax() <= fXaxis.GetXmin()) {
Double_t xmin = fBuffer[2];
Double_t xmax = xmin;
for (Int_t i=1;i<nbentries;i++) {
Double_t x = fBuffer[3*i+2];
if (x < xmin) xmin = x;
if (x > xmax) xmax = x;
}
if (fXaxis.GetXmax() <= fXaxis.GetXmin()) {
THLimitsFinder::GetLimitsFinder()->FindGoodLimits(this,xmin,xmax);
} else {
fBuffer = 0;
Int_t keep = fBufferSize; fBufferSize = 0;
if (xmin < fXaxis.GetXmin()) RebinAxis(xmin,&fXaxis);
if (xmax >= fXaxis.GetXmax()) RebinAxis(xmax,&fXaxis);
fBuffer = buffer;
fBufferSize = keep;
}
}
fBuffer = 0;
for (Int_t i=0;i<nbentries;i++) {
Fill(buffer[3*i+2],buffer[3*i+3],buffer[3*i+1]);
}
fBuffer = buffer;
if (action > 0) { delete [] fBuffer; fBuffer = 0; fBufferSize = 0;}
else {
if (nbentries == (Int_t)fEntries) fBuffer[0] = -nbentries;
else fBuffer[0] = 0;
}
return nbentries;
}
Int_t TProfile::BufferFill(Double_t x, Double_t y, Double_t w)
{
if (!fBuffer) return -2;
Int_t nbentries = (Int_t)fBuffer[0];
if (nbentries < 0) {
nbentries = -nbentries;
fBuffer[0] = nbentries;
if (fEntries > 0) {
Double_t *buffer = fBuffer; fBuffer=0;
Reset();
fBuffer = buffer;
}
}
if (3*nbentries+3 >= fBufferSize) {
BufferEmpty(1);
return Fill(x,y,w);
}
fBuffer[3*nbentries+1] = w;
fBuffer[3*nbentries+2] = x;
fBuffer[3*nbentries+3] = y;
fBuffer[0] += 1;
return -2;
}
void TProfile::Copy(TObject &obj) const
{
TH1D::Copy(((TProfile&)obj));
fBinEntries.Copy(((TProfile&)obj).fBinEntries);
for (int bin=0;bin<fNcells;bin++) {
((TProfile&)obj).fArray[bin] = fArray[bin];
((TProfile&)obj).fSumw2.fArray[bin] = fSumw2.fArray[bin];
}
((TProfile&)obj).fYmin = fYmin;
((TProfile&)obj).fYmax = fYmax;
((TProfile&)obj).fScaling = fScaling;
((TProfile&)obj).fErrorMode = fErrorMode;
((TProfile&)obj).fTsumwy = fTsumwy;
((TProfile&)obj).fTsumwy2 = fTsumwy2;
}
void TProfile::Divide(TF1 *, Double_t )
{
Error("Divide","Function not implemented for TProfile");
return;
}
void TProfile::Divide(const TH1 *h1)
{
if (!h1) {
Error("Divide","Attempt to divide a non-existing profile");
return;
}
if (!h1->InheritsFrom("TH1")) {
Error("Divide","Attempt to divide by a non-profile or non-histogram object");
return;
}
TProfile *p1 = (TProfile*)h1;
Int_t nbinsx = GetNbinsX();
if (nbinsx != p1->GetNbinsX()) {
Error("Divide","Attempt to divide profiles with different number of bins");
return;
}
fEntries = fTsumw = fTsumw2 = fTsumwx = fTsumwx2 = fTsumwy = fTsumwy2 = 0;
Int_t bin;
Double_t *cu1=0, *er1=0, *en1=0;
Double_t e0,e1,c12;
if (h1->InheritsFrom("TProfile")) {
cu1 = p1->GetW();
er1 = p1->GetW2();
en1 = p1->GetB();
}
Double_t c0,c1,w,z,x;
for (bin=0;bin<=nbinsx+1;bin++) {
c0 = fArray[bin];
if (cu1) c1 = cu1[bin];
else c1 = h1->GetBinContent(bin);
if (c1) w = c0/c1;
else w = 0;
fArray[bin] = w;
z = TMath::Abs(w);
x = fXaxis.GetBinCenter(bin);
fEntries++;
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
fTsumwy += z*c1;
fTsumwx2 += z*c1*c1;
e0 = fSumw2.fArray[bin];
if (er1) e1 = er1[bin];
else {e1 = h1->GetBinError(bin); e1*=e1;}
c12= c1*c1;
if (!c1) fSumw2.fArray[bin] = 0;
else fSumw2.fArray[bin] = (e0*c1*c1 + e1*c0*c0)/(c12*c12);
if (!en1) continue;
if (!en1[bin]) fBinEntries.fArray[bin] = 0;
else fBinEntries.fArray[bin] /= en1[bin];
}
}
void TProfile::Divide(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
{
TString opt = option;
opt.ToLower();
Bool_t binomial = kFALSE;
if (opt.Contains("b")) binomial = kTRUE;
if (!h1 || !h2) {
Error("Divide","Attempt to divide a non-existing profile");
return;
}
if (!h1->InheritsFrom("TProfile")) {
Error("Divide","Attempt to divide a non-profile object");
return;
}
TProfile *p1 = (TProfile*)h1;
if (!h2->InheritsFrom("TProfile")) {
Error("Divide","Attempt to divide by a non-profile object");
return;
}
TProfile *p2 = (TProfile*)h2;
Int_t nbinsx = GetNbinsX();
if (nbinsx != p1->GetNbinsX() || nbinsx != p2->GetNbinsX()) {
Error("Divide","Attempt to divide profiles with different number of bins");
return;
}
if (!c2) {
Error("Divide","Coefficient of dividing profile cannot be zero");
return;
}
printf("WARNING!!: The algorithm in TProfile::Divide computing the errors is not accurate\n");
printf(" Instead of Divide(TProfile *h1, TProfile *h2, do:\n");
printf(" TH1D *p1 = h1->ProjectionX();\n");
printf(" TH1D *p2 = h2->ProjectionX();\n");
printf(" p1->Divide(p2);\n");
fEntries = fTsumw = fTsumw2 = fTsumwx = fTsumwx2 = 0;
Int_t bin;
Double_t *cu1 = p1->GetW();
Double_t *cu2 = p2->GetW();
Double_t *er1 = p1->GetW2();
Double_t *er2 = p2->GetW2();
Double_t *en1 = p1->GetB();
Double_t *en2 = p2->GetB();
Double_t b1,b2,w,z,x,ac1,ac2;
ac1 = TMath::Abs(c1);
ac2 = TMath::Abs(c2);
for (bin=0;bin<=nbinsx+1;bin++) {
b1 = cu1[bin];
b2 = cu2[bin];
if (b2) w = c1*b1/(c2*b2);
else w = 0;
fArray[bin] = w;
z = TMath::Abs(w);
x = fXaxis.GetBinCenter(bin);
fEntries++;
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
Double_t e1 = er1[bin];
Double_t e2 = er2[bin];
Double_t b22= b2*b2*TMath::Abs(c2);
if (!b2) fSumw2.fArray[bin] = 0;
else {
if (binomial) {
fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/b2);
} else {
fSumw2.fArray[bin] = ac1*ac2*(e1*b2*b2 + e2*b1*b1)/(b22*b22);
}
}
if (en2[bin]) fBinEntries.fArray[bin] = en1[bin]/en2[bin];
else fBinEntries.fArray[bin] = 0;
}
}
TH1 *TProfile::DrawCopy(Option_t *option) const
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
TProfile *newpf = (TProfile*)Clone();
newpf->SetDirectory(0);
newpf->SetBit(kCanDelete);
newpf->AppendPad(option);
return newpf;
}
Int_t TProfile::Fill(Double_t x, Double_t y)
{
if (fBuffer) return BufferFill(x,y,1);
Int_t bin;
if (fYmin != fYmax) {
if (y <fYmin || y> fYmax) return -1;
}
fEntries++;
bin =fXaxis.FindBin(x);
AddBinContent(bin, y);
fSumw2.fArray[bin] += (Double_t)y*y;
fBinEntries.fArray[bin] += 1;
if (bin == 0 || bin > fXaxis.GetNbins()) {
if (!fgStatOverflows) return -1;
}
fTsumw++;
fTsumw2++;
fTsumwx += x;
fTsumwx2 += x*x;
fTsumwy += y;
fTsumwy2 += y*y;
return bin;
}
Int_t TProfile::Fill(const char *namex, Double_t y)
{
Int_t bin;
if (fYmin != fYmax) {
if (y <fYmin || y> fYmax) return -1;
}
fEntries++;
bin =fXaxis.FindBin(namex);
AddBinContent(bin, y);
fSumw2.fArray[bin] += (Double_t)y*y;
fBinEntries.fArray[bin] += 1;
if (bin == 0 || bin > fXaxis.GetNbins()) {
if (!fgStatOverflows) return -1;
}
Double_t x = fXaxis.GetBinCenter(bin);
fTsumw++;
fTsumw2++;
fTsumwx += x;
fTsumwx2 += x*x;
fTsumwy += y;
fTsumwy2 += y*y;
return bin;
}
Int_t TProfile::Fill(Double_t x, Double_t y, Double_t w)
{
if (fBuffer) return BufferFill(x,y,w);
Int_t bin;
if (fYmin != fYmax) {
if (y <fYmin || y> fYmax) return -1;
}
Double_t z= (w > 0 ? w : -w);
fEntries++;
bin =fXaxis.FindBin(x);
AddBinContent(bin, z*y);
fSumw2.fArray[bin] += z*y*y;
fBinEntries.fArray[bin] += z;
if (bin == 0 || bin > fXaxis.GetNbins()) {
if (!fgStatOverflows) return -1;
}
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
fTsumwy += z*y;
fTsumwy2 += z*y*y;
return bin;
}
Int_t TProfile::Fill(const char *namex, Double_t y, Double_t w)
{
Int_t bin;
if (fYmin != fYmax) {
if (y <fYmin || y> fYmax) return -1;
}
Double_t z= (w > 0 ? w : -w);
fEntries++;
bin =fXaxis.FindBin(namex);
AddBinContent(bin, z*y);
fSumw2.fArray[bin] += z*y*y;
fBinEntries.fArray[bin] += z;
if (bin == 0 || bin > fXaxis.GetNbins()) {
if (!fgStatOverflows) return -1;
}
Double_t x = fXaxis.GetBinCenter(bin);
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x;
fTsumwx2 += z*x*x;
fTsumwy += z*y;
fTsumwy2 += z*y*y;
return bin;
}
void TProfile::FillN(Int_t ntimes, const Double_t *x, const Double_t *y, const Double_t *w, Int_t stride)
{
Int_t bin,i;
ntimes *= stride;
for (i=0;i<ntimes;i+=stride) {
if (fYmin != fYmax) {
if (y[i] <fYmin || y[i]> fYmax) continue;
}
Double_t z= (w[i] > 0 ? w[i] : -w[i]);
fEntries++;
bin =fXaxis.FindBin(x[i]);
AddBinContent(bin, z*y[i]);
fSumw2.fArray[bin] += z*y[i]*y[i];
fBinEntries.fArray[bin] += z;
if (bin == 0 || bin > fXaxis.GetNbins()) {
if (!fgStatOverflows) continue;
}
fTsumw += z;
fTsumw2 += z*z;
fTsumwx += z*x[i];
fTsumwx2 += z*x[i]*x[i];
fTsumwy += z*y[i];
fTsumwy2 += z*y[i]*y[i];
}
}
Double_t TProfile::GetBinContent(Int_t bin) const
{
if (fBuffer) ((TProfile*)this)->BufferEmpty();
if (bin < 0 || bin >= fNcells) return 0;
if (fBinEntries.fArray[bin] == 0) return 0;
if (!fArray) return 0;
return fArray[bin]/fBinEntries.fArray[bin];
}
Double_t TProfile::GetBinEntries(Int_t bin) const
{
if (fBuffer) ((TProfile*)this)->BufferEmpty();
if (bin < 0 || bin >= fNcells) return 0;
return fBinEntries.fArray[bin];
}
Double_t TProfile::GetBinError(Int_t bin) const
{
if (fBuffer) ((TProfile*)this)->BufferEmpty();
if (bin < 0 || bin >= fNcells) return 0;
Double_t cont = fArray[bin];
Double_t sum = fBinEntries.fArray[bin];
Double_t err2 = fSumw2.fArray[bin];
if (sum == 0) return 0;
Double_t eprim;
Double_t contsum = cont/sum;
Double_t eprim2 = TMath::Abs(err2/sum - contsum*contsum);
eprim = TMath::Sqrt(eprim2);
Double_t test = 1;
if (err2 != 0 && sum < 5) test = eprim2*sum/err2;
if (fgApproximate && fNcells <=1002 && (test < 1.e-4 || eprim2 < 1e-6)) {
Double_t scont, ssum, serr2;
scont = ssum = serr2 = 0;
for (Int_t i=1;i<fNcells;i++) {
if (fSumw2.fArray[i] <= 0) continue;
scont += fArray[i];
ssum += fBinEntries.fArray[i];
serr2 += fSumw2.fArray[i];
}
Double_t scontsum = scont/ssum;
Double_t seprim2 = TMath::Abs(serr2/ssum - scontsum*scontsum);
eprim = 2*TMath::Sqrt(seprim2);
sum = ssum;
}
sum = TMath::Abs(sum);
if (fErrorMode == kERRORMEAN) return eprim/TMath::Sqrt(sum);
else if (fErrorMode == kERRORSPREAD) return eprim;
else if (fErrorMode == kERRORSPREADI) {
if (eprim != 0) return eprim/TMath::Sqrt(sum);
return 1/TMath::Sqrt(12*sum);
}
else if (fErrorMode == kERRORSPREADG) {
return 1./TMath::Sqrt(sum);
}
else return eprim;
}
Option_t *TProfile::GetErrorOption() const
{
if (fErrorMode == kERRORSPREAD) return "s";
if (fErrorMode == kERRORSPREADI) return "i";
if (fErrorMode == kERRORSPREADG) return "g";
return "";
}
char* TProfile::GetObjectInfo(Int_t px, Int_t py) const
{
if (!gPad) return (char*)"";
static char info[64];
Double_t x = gPad->PadtoX(gPad->AbsPixeltoX(px));
Double_t y = gPad->PadtoY(gPad->AbsPixeltoY(py));
Int_t binx = GetXaxis()->FindFixBin(x);
sprintf(info,"(x=%g, y=%g, binx=%d, binc=%g, bine=%g, binn=%d)", x, y, binx, GetBinContent(binx), GetBinError(binx), (Int_t)GetBinEntries(binx));
return info;
}
void TProfile::GetStats(Double_t *stats) const
{
if (fBuffer) ((TProfile*)this)->BufferEmpty();
Int_t bin, binx;
if (fTsumw == 0 || fXaxis.TestBit(TAxis::kAxisRange)) {
Double_t w;
Double_t x;
for (bin=0;bin<6;bin++) stats[bin] = 0;
if (!fBinEntries.fArray) return;
for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
w = fBinEntries.fArray[binx];
x = fXaxis.GetBinCenter(binx);
stats[0] += w;
stats[1] += w*w;
stats[2] += w*x;
stats[3] += w*x*x;
stats[4] += fArray[binx];
stats[5] += fSumw2.fArray[binx];
}
} else {
if (fTsumwy == 0 && fTsumwy2 == 0) {
TProfile *p = (TProfile*)this;
for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
p->fTsumwy += fArray[binx];
p->fTsumwy2 += fSumw2.fArray[binx];
}
}
stats[0] = fTsumw;
stats[1] = fTsumw2;
stats[2] = fTsumwx;
stats[3] = fTsumwx2;
stats[4] = fTsumwy;
stats[5] = fTsumwy2;
}
}
void TProfile::LabelsDeflate(Option_t *)
{
if (!fXaxis.GetLabels()) return;
TIter next(fXaxis.GetLabels());
TObject *obj;
Int_t nbins = 0;
while ((obj = next())) {
if (obj->GetUniqueID()) nbins++;
}
if (nbins < 2) nbins = 2;
TProfile *hold = (TProfile*)Clone();
hold->SetDirectory(0);
Double_t xmin = fXaxis.GetXmin();
Double_t xmax = fXaxis.GetBinUpEdge(nbins);
fXaxis.SetRange(0,0);
fXaxis.Set(nbins,xmin,xmax);
Int_t ncells = nbins+2;
SetBinsLength(ncells);
fBinEntries.Set(ncells);
fSumw2.Set(ncells);
Int_t bin;
for (bin=1;bin<=nbins;bin++) {
fArray[bin] = hold->fArray[bin];
fBinEntries.fArray[bin] = hold->fBinEntries.fArray[bin];
fSumw2.fArray[bin] = hold->fSumw2.fArray[bin];
}
delete hold;
}
void TProfile::LabelsInflate(Option_t *)
{
TProfile *hold = (TProfile*)Clone();
hold->SetDirectory(0);
Int_t nbold = fXaxis.GetNbins();
Int_t nbins = nbold;
Double_t xmin = fXaxis.GetXmin();
Double_t xmax = fXaxis.GetXmax();
xmax = xmin + 2*(xmax-xmin);
fXaxis.SetRange(0,0);
fXaxis.Set(2*nbins,xmin,xmax);
nbins *= 2;
Int_t ncells = nbins+2;
SetBinsLength(ncells);
fBinEntries.Set(ncells);
fSumw2.Set(ncells);
Int_t bin;
for (bin=1;bin<=nbins;bin++) {
if (bin <= nbold) {
fArray[bin] = hold->fArray[bin];
fBinEntries.fArray[bin] = hold->fBinEntries.fArray[bin];
fSumw2.fArray[bin] = hold->fSumw2.fArray[bin];
} else {
fArray[bin] = 0;
fBinEntries.fArray[bin] = 0;
fSumw2.fArray[bin] = 0;
}
}
delete hold;
}
void TProfile::LabelsOption(Option_t *option, Option_t * )
{
THashList *labels = fXaxis.GetLabels();
if (!labels) {
Warning("LabelsOption","Cannot sort. No labels");
return;
}
TString opt = option;
opt.ToLower();
if (opt.Contains("h")) {
fXaxis.SetBit(TAxis::kLabelsHori);
fXaxis.ResetBit(TAxis::kLabelsVert);
fXaxis.ResetBit(TAxis::kLabelsDown);
fXaxis.ResetBit(TAxis::kLabelsUp);
}
if (opt.Contains("v")) {
fXaxis.SetBit(TAxis::kLabelsVert);
fXaxis.ResetBit(TAxis::kLabelsHori);
fXaxis.ResetBit(TAxis::kLabelsDown);
fXaxis.ResetBit(TAxis::kLabelsUp);
}
if (opt.Contains("u")) {
fXaxis.SetBit(TAxis::kLabelsUp);
fXaxis.ResetBit(TAxis::kLabelsVert);
fXaxis.ResetBit(TAxis::kLabelsDown);
fXaxis.ResetBit(TAxis::kLabelsHori);
}
if (opt.Contains("d")) {
fXaxis.SetBit(TAxis::kLabelsDown);
fXaxis.ResetBit(TAxis::kLabelsVert);
fXaxis.ResetBit(TAxis::kLabelsHori);
fXaxis.ResetBit(TAxis::kLabelsUp);
}
Int_t sort = -1;
if (opt.Contains("a")) sort = 0;
if (opt.Contains(">")) sort = 1;
if (opt.Contains("<")) sort = 2;
if (sort < 0) return;
Int_t n = TMath::Min(fXaxis.GetNbins(), labels->GetSize());
Int_t *a = new Int_t[n+2];
Int_t i,j;
Double_t *cont = new Double_t[n+2];
Double_t *sumw = new Double_t[n+2];
Double_t *errors = new Double_t[n+2];
Double_t *ent = new Double_t[n+2];
THashList *labold = new THashList(labels->GetSize(),1);
TIter nextold(labels);
TObject *obj;
while ((obj=nextold())) {
labold->Add(obj);
}
labels->Clear();
if (sort > 0) {
for (i=1;i<=n;i++) {
sumw[i-1] = fArray[i];
errors[i-1] = fSumw2.fArray[i];
ent[i-1] = fBinEntries.fArray[i];
if (fBinEntries.fArray[i] == 0) cont[i-1] = 0;
else cont[i-1] = fArray[i]/fBinEntries.fArray[i];
}
if (sort ==1) TMath::Sort(n,cont,a,kTRUE);
else TMath::Sort(n,cont,a,kFALSE);
for (i=1;i<=n;i++) {
fArray[i] = sumw[a[i-1]];
fSumw2.fArray[i] = errors[a[i-1]];
fBinEntries.fArray[i] = ent[a[i-1]];
}
for (i=1;i<=n;i++) {
obj = labold->At(a[i-1]);
labels->Add(obj);
obj->SetUniqueID(i);
}
} else {
const UInt_t kUsed = 1<<18;
TObject *objk=0;
a[0] = 0;
a[n+1] = n+1;
for (i=1;i<=n;i++) {
const char *label = "zzzzzzzzzzzz";
for (j=1;j<=n;j++) {
obj = labold->At(j-1);
if (!obj) continue;
if (obj->TestBit(kUsed)) continue;
if (strcmp(label,obj->GetName()) < 0) continue;
objk = obj;
a[i] = j;
label = obj->GetName();
}
if (objk) {
objk->SetUniqueID(i);
labels->Add(objk);
objk->SetBit(kUsed);
}
}
for (i=1;i<=n;i++) {
obj = labels->At(i-1);
if (!obj) continue;
obj->ResetBit(kUsed);
}
for (i=1;i<=n;i++) {
sumw[i] = fArray[a[i]];
errors[i] = fSumw2.fArray[a[i]];
ent[i] = fBinEntries.fArray[a[i]];
}
for (i=1;i<=n;i++) {
fArray[i] = sumw[i];
fSumw2.fArray[i] = errors[i];
fBinEntries.fArray[i] = ent[i];
}
}
delete labold;
if (a) delete [] a;
if (sumw) delete [] sumw;
if (cont) delete [] cont;
if (errors) delete [] errors;
if (ent) delete [] ent;
}
Long64_t TProfile::Merge(TCollection *li)
{
if (!li) return 0;
if (li->IsEmpty()) return (Int_t) GetEntries();
TList inlist;
TH1* hclone = (TH1*)Clone("FirstClone");
R__ASSERT(hclone);
BufferEmpty(1);
Reset();
SetEntries(0);
inlist.Add(hclone);
inlist.AddAll(li);
TAxis newXAxis;
Bool_t initialLimitsFound = kFALSE;
Bool_t same = kTRUE;
Bool_t allHaveLimits = kTRUE;
TIter next(&inlist);
while (TObject *o = next()) {
TProfile* h = dynamic_cast<TProfile*> (o);
if (!h) {
Error("Add","Attempt to add object of class: %s to a %s",
o->ClassName(),this->ClassName());
return -1;
}
Bool_t hasLimits = h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax();
allHaveLimits = allHaveLimits && hasLimits;
if (hasLimits) {
h->BufferEmpty();
if (!initialLimitsFound) {
initialLimitsFound = kTRUE;
newXAxis.Set(h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),
h->GetXaxis()->GetXmax());
}
else {
if (!RecomputeAxisLimits(newXAxis, *(h->GetXaxis()))) {
Error("Merge", "Cannot merge histograms - limits are inconsistent:\n "
"first: (%d, %f, %f), second: (%d, %f, %f)",
newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax(),
h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),
h->GetXaxis()->GetXmax());
}
}
}
}
next.Reset();
same = same && SameLimitsAndNBins(newXAxis, *GetXaxis());
if (!same && initialLimitsFound)
SetBins(newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax());
if (!allHaveLimits) {
while (TProfile* h = (TProfile*)next()) {
if (h->GetXaxis()->GetXmin() >= h->GetXaxis()->GetXmax() && h->fBuffer) {
Int_t nbentries = (Int_t)h->fBuffer[0];
for (Int_t i = 0; i < nbentries; i++)
Fill(h->fBuffer[3*i + 2], h->fBuffer[3*i + 3], h->fBuffer[3*i + 1]);
}
}
if (!initialLimitsFound)
return (Int_t) GetEntries();
next.Reset();
}
Double_t stats[kNstat], totstats[kNstat];
for (Int_t i=0;i<kNstat;i++) {totstats[i] = stats[i] = 0;}
GetStats(totstats);
Double_t nentries = GetEntries();
Int_t binx, ix, nx;
Bool_t canRebin=TestBit(kCanRebin);
ResetBit(kCanRebin);
while (TProfile* h=(TProfile*)next()) {
if (h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax()) {
h->GetStats(stats);
for (Int_t i=0;i<kNstat;i++)
totstats[i] += stats[i];
nentries += h->GetEntries();
nx = h->GetXaxis()->GetNbins();
for (binx = 0; binx <= nx + 1; binx++) {
if ((!same) && (binx == 0 || binx == nx + 1)) {
if (h->GetW()[binx] != 0) {
Error("Merge", "Cannot merge histograms - the histograms have"
" different limits and undeflows/overflows are present."
" The initial histogram is now broken!");
return -1;
}
}
ix = fXaxis.FindBin(h->GetBinCenter(binx));
fArray[ix] += h->GetW()[binx];
fSumw2.fArray[ix] += h->GetW2()[binx];
fBinEntries.fArray[ix] += h->GetB()[binx];
}
fEntries += h->GetEntries();
fTsumw += h->fTsumw;
fTsumw2 += h->fTsumw2;
fTsumwx += h->fTsumwx;
fTsumwx2 += h->fTsumwx2;
fTsumwy += h->fTsumwy;
}
}
if (canRebin) SetBit(kCanRebin);
PutStats(totstats);
SetEntries(nentries);
inlist.Remove(hclone);
delete hclone;
return (Long64_t) nentries;
}
void TProfile::Multiply(TF1 *f1, Double_t c1)
{
if (!f1) {
Error("Multiply","Attempt to multiply by a null function");
return;
}
Int_t nbinsx = GetNbinsX();
Double_t xx[1], cf1, ac1 = TMath::Abs(c1);
Double_t s1[10];
Int_t i;
for (i=0;i<10;i++) {s1[i] = 0;}
PutStats(s1);
SetMinimum();
SetMaximum();
Int_t bin;
for (bin=0;bin<=nbinsx+1;bin++) {
xx[0] = fXaxis.GetBinCenter(bin);
if (!f1->IsInside(xx)) continue;
TF1::RejectPoint(kFALSE);
cf1 = f1->EvalPar(xx);
if (TF1::RejectedPoint()) continue;
fArray[bin] *= c1*cf1;
fSumw2.fArray[bin] *= ac1*cf1*cf1;
}
}
void TProfile::Multiply(const TH1 *)
{
Error("Multiply","Multiplication of profile histograms not implemented");
}
void TProfile::Multiply(const TH1 *, const TH1 *, Double_t, Double_t, Option_t *)
{
Error("Multiply","Multiplication of profile histograms not implemented");
}
TH1D *TProfile::ProjectionX(const char *name, Option_t *option) const
{
TString opt = option;
opt.ToLower();
Int_t nx = fXaxis.GetNbins();
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;
const TArrayD *bins = fXaxis.GetXbins();
if (bins->fN == 0) {
h1 = new TH1D(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax());
} else {
h1 = new TH1D(pname,GetTitle(),nx,bins->fArray);
}
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 || binWeight ) h1->Sumw2();
if (pname != name) delete [] pname;
Double_t cont;
for (Int_t bin =0;bin<=nx+1;bin++) {
if (binEntries) cont = GetBinEntries(bin);
else if (cequalErrors) cont = GetBinError(bin);
else if (binWeight) cont = GetBinContent(bin) * GetBinEntries(bin);
else cont = GetBinContent(bin);
h1->SetBinContent(bin ,cont);
if (computeErrors ) h1->SetBinError(bin , GetBinError(bin) );
if (binWeight) h1->SetBinError(bin , TMath::Sqrt(fSumw2.fArray[bin] ) );
}
h1->GetXaxis()->ImportAttributes(this->GetXaxis());
h1->GetYaxis()->ImportAttributes(this->GetYaxis());
THashList* labels=this->GetXaxis()->GetLabels();
if (labels) {
TIter iL(labels);
TObjString* lb;
Int_t i = 1;
while ((lb=(TObjString*)iL())) {
h1->GetXaxis()->SetBinLabel(i,lb->String().Data());
i++;
}
}
h1->SetEntries(fEntries);
return h1;
}
void TProfile::PutStats(Double_t *stats)
{
fTsumw = stats[0];
fTsumw2 = stats[1];
fTsumwx = stats[2];
fTsumwx2 = stats[3];
fTsumwy = stats[4];
fTsumwy2 = stats[5];
}
TH1 *TProfile::Rebin(Int_t ngroup, const char*newname, const Double_t *xbins)
{
Int_t nbins = fXaxis.GetNbins();
Double_t xmin = fXaxis.GetXmin();
Double_t xmax = fXaxis.GetXmax();
if ((ngroup <= 0) || (ngroup > nbins)) {
Error("Rebin", "Illegal value of ngroup=%d",ngroup);
return 0;
}
if (!newname && xbins) {
Error("Rebin","if xbins is specified, newname must be given");
return 0;
}
Int_t newbins = nbins/ngroup;
if (!xbins) {
Int_t nbg = nbins/ngroup;
if (nbg*ngroup != nbins) {
Warning("Rebin", "ngroup=%d must be an exact divider of nbins=%d",ngroup,nbins);
}
}
else {
newbins = ngroup;
ngroup = nbins;
}
Double_t *oldBins = new Double_t[nbins+2];
Double_t *oldCount = new Double_t[nbins+2];
Double_t *oldErrors = new Double_t[nbins+2];
Int_t bin, i;
Double_t *cu1 = GetW();
Double_t *er1 = GetW2();
Double_t *en1 = GetB();
for (bin=0;bin<=nbins+1;bin++) {
oldBins[bin] = cu1[bin];
oldCount[bin] = en1[bin];
oldErrors[bin] = er1[bin];
}
TProfile *hnew = this;
if ((newname && strlen(newname) > 0) || xbins) {
hnew = (TProfile*)Clone(newname);
}
if(!xbins && (newbins*ngroup != nbins)) {
xmax = fXaxis.GetBinUpEdge(newbins*ngroup);
hnew->fTsumw = 0;
}
if(!xbins && (fXaxis.GetXbins()->GetSize() > 0)){
Double_t *bins = new Double_t[newbins+1];
for(i = 0; i <= newbins; ++i) bins[i] = fXaxis.GetBinLowEdge(1+i*ngroup);
hnew->SetBins(newbins,bins);
delete [] bins;
} else if (xbins) {
hnew->SetBins(newbins,xbins);
} else {
hnew->SetBins(newbins,xmin,xmax);
}
Double_t *cu2 = hnew->GetW();
Double_t *er2 = hnew->GetW2();
Double_t *en2 = hnew->GetB();
Int_t oldbin = 1;
Double_t binContent, binCount, binError;
for (bin = 1;bin<=newbins;bin++) {
binContent = 0;
binCount = 0;
binError = 0;
Int_t imax = ngroup;
Double_t xbinmax = hnew->GetXaxis()->GetBinUpEdge(bin);
for (i=0;i<ngroup;i++) {
if( (hnew == this && (oldbin+i > nbins)) ||
(hnew != this && (fXaxis.GetBinCenter(oldbin+i) > xbinmax)) ) {
imax = i;
break;
}
binContent += oldBins[oldbin+i];
binCount += oldCount[oldbin+i];
binError += oldErrors[oldbin+i];
}
cu2[bin] = binContent;
er2[bin] = binError;
en2[bin] = binCount;
oldbin += imax;
}
hnew->fArray[0] = oldBins[0];
hnew->fArray[newbins+1] = oldBins[nbins+1];
hnew->fBinEntries[0] = oldCount[0];
hnew->fBinEntries[newbins+1] = oldCount[nbins+1];
hnew->fSumw2[0] = oldErrors[0];
hnew->fSumw2[newbins+1] = oldErrors[nbins+1];
delete [] oldBins;
delete [] oldCount;
delete [] oldErrors;
return hnew;
}
void TProfile::RebinAxis(Double_t x, TAxis *axis)
{
if (!TestBit(kCanRebin)) return;
if (axis->GetXmin() >= axis->GetXmax()) return;
if (axis->GetNbins() <= 0) return;
Double_t xmin, xmax;
if (!FindNewAxisLimits(axis, x, xmin, xmax))
return;
TProfile *hold = (TProfile*)Clone();
hold->SetDirectory(0);
axis->SetLimits(xmin,xmax);
Int_t nbinsx = fXaxis.GetNbins();
Reset("ICE");
for (Int_t binx = 1; binx <= nbinsx; binx++) {
Int_t destinationBin = fXaxis.FindFixBin(hold->GetXaxis()->GetBinCenter(binx));
Int_t sourceBin = binx;
AddBinContent(destinationBin, hold->fArray[sourceBin]);
fBinEntries.fArray[destinationBin] += hold->fBinEntries.fArray[sourceBin];
fSumw2.fArray[destinationBin] += hold->fSumw2.fArray[sourceBin];
}
fTsumwy = hold->fTsumwy;
fTsumwy2 = hold->fTsumwy2;
delete hold;
}
void TProfile::Reset(Option_t *option)
{
TH1D::Reset(option);
fBinEntries.Reset();
TString opt = option;
opt.ToUpper();
if (opt.Contains("ICE")) return;
fTsumwy = 0;
fTsumwy2 = 0;
}
void TProfile::SavePrimitive(ostream &out, Option_t *option )
{
Bool_t nonEqiX = kFALSE;
Int_t i;
if (GetXaxis()->GetXbins()->fN && GetXaxis()->GetXbins()->fArray) {
nonEqiX = kTRUE;
out << " Double_t xAxis[" << GetXaxis()->GetXbins()->fN
<< "] = {";
for (i = 0; i < GetXaxis()->GetXbins()->fN; i++) {
if (i != 0) out << ", ";
out << GetXaxis()->GetXbins()->fArray[i];
}
out << "}; " << endl;
}
char quote = '"';
out<<" "<<endl;
out<<" "<<ClassName()<<" *";
out<<GetName()<<" = new "<<ClassName()<<"("<<quote<<GetName()<<quote<<","<<quote<<GetTitle()<<quote
<<","<<GetXaxis()->GetNbins();
if (nonEqiX)
out << ", xAxis";
else
out << "," << GetXaxis()->GetXmin()
<< "," << GetXaxis()->GetXmax()
<<","<<quote<<GetErrorOption()<<quote<<");"<<endl;
Int_t bin;
for (bin=0;bin<fNcells;bin++) {
Double_t bi = GetBinEntries(bin);
if (bi) {
out<<" "<<GetName()<<"->SetBinEntries("<<bin<<","<<bi<<");"<<endl;
}
}
for (bin=0;bin<fNcells;bin++) {
Double_t bc = fArray[bin];
if (bc) {
out<<" "<<GetName()<<"->SetBinContent("<<bin<<","<<bc<<");"<<endl;
}
}
if (fSumw2.fN) {
for (bin=0;bin<fNcells;bin++) {
Double_t be = TMath::Sqrt(fSumw2.fArray[bin]);
if (be) {
out<<" "<<GetName()<<"->SetBinError("<<bin<<","<<be<<");"<<endl;
}
}
}
TH1::SavePrimitiveHelp(out, option);
}
void TProfile::Scale(Double_t c1, Option_t *)
{
Double_t ent = fEntries;
fScaling = kTRUE;
Add(this,this,c1,0);
fScaling = kFALSE;
fEntries = ent;
}
void TProfile::SetBinEntries(Int_t bin, Double_t w)
{
if (bin < 0 || bin >= fNcells) return;
fBinEntries.fArray[bin] = w;
}
void TProfile::SetBins(Int_t nx, Double_t xmin, Double_t xmax)
{
fXaxis.Set(nx,xmin,xmax);
fNcells = nx+2;
SetBinsLength(fNcells);
fBinEntries.Set(fNcells);
fSumw2.Set(fNcells);
}
void TProfile::SetBins(Int_t nx, const Double_t *xbins)
{
fXaxis.Set(nx,xbins);
fNcells = nx+2;
SetBinsLength(fNcells);
fBinEntries.Set(fNcells);
fSumw2.Set(fNcells);
}
void TProfile::SetBuffer(Int_t buffersize, Option_t *)
{
if (fBuffer) {
BufferEmpty();
delete [] fBuffer;
fBuffer = 0;
}
if (buffersize <= 0) {
fBufferSize = 0;
return;
}
if (buffersize < 100) buffersize = 100;
fBufferSize = 1 + 3*buffersize;
fBuffer = new Double_t[fBufferSize];
memset(fBuffer,0,8*fBufferSize);
}
void TProfile::SetErrorOption(Option_t *option)
{
TString opt = option;
opt.ToLower();
fErrorMode = kERRORMEAN;
if (opt.Contains("s")) fErrorMode = kERRORSPREAD;
if (opt.Contains("i")) fErrorMode = kERRORSPREADI;
if (opt.Contains("g")) fErrorMode = kERRORSPREADG;
}
void TProfile::Streamer(TBuffer &R__b)
{
if (R__b.IsReading()) {
UInt_t R__s, R__c;
Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v > 2) {
R__b.ReadClassBuffer(TProfile::Class(), this, R__v, R__s, R__c);
return;
}
TH1D::Streamer(R__b);
fBinEntries.Streamer(R__b);
Int_t errorMode;
R__b >> errorMode;
fErrorMode = (EErrorType)errorMode;
if (R__v < 2) {
Float_t ymin,ymax;
R__b >> ymin; fYmin = ymin;
R__b >> ymax; fYmax = ymax;
} else {
R__b >> fYmin;
R__b >> fYmax;
}
R__b.CheckByteCount(R__s, R__c, TProfile::IsA());
} else {
R__b.WriteClassBuffer(TProfile::Class(),this);
}
}
Last change: Fri Dec 5 09:53:13 2008
Last generated: 2008-12-05 09:53
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.