ROOT  6.06/09
Reference Guide
TProfile2D.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Author: Rene Brun 16/04/2000
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #include "TProfile2D.h"
13 #include "TMath.h"
14 #include "THLimitsFinder.h"
15 #include "Riostream.h"
16 #include "TVirtualPad.h"
17 #include "TError.h"
18 #include "TClass.h"
19 
20 #include "TProfileHelper.h"
21 
23 
25 
26 /** \class TProfile2D
27  \ingroup Hist
28  Profile2D histograms are used to display the mean
29  value of Z and its RMS for each cell in X,Y.
30  Profile2D histograms are in many cases an
31  elegant replacement of three-dimensional histograms : the inter-relation of three
32  measured quantities X, Y and Z can always be visualized by a three-dimensional
33  histogram or scatter-plot; its representation on the line-printer is not particularly
34  satisfactory, except for sparse data. If Z is an unknown (but single-valued)
35  approximate function of X,Y this function is displayed by a profile2D histogram with
36  much better precision than by a scatter-plot.
37 
38  The following formulae show the cumulated contents (capital letters) and the values
39  displayed by the printing or plotting routines (small letters) of the elements for cell I, J.
40 
41  2
42  H(I,J) = sum Z E(I,J) = sum Z
43  l(I,J) = sum l L(I,J) = sum l
44  h(I,J) = H(I,J)/L(I,J) s(I,J) = sqrt(E(I,J)/L(I,J)- h(I,J)**2)
45  e(I,J) = s(I,J)/sqrt(L(I,J))
46 
47  In the special case where s(I,J) is zero (eg, case of 1 entry only in one cell)
48  the bin error e(I,J) is computed from the average of the s(I,J) for all cells
49  if the static function TProfile2D::Approximate has been called.
50  This simple/crude approximation was suggested in order to keep the cell
51  during a fit operation. But note that this approximation is not the default behaviour.
52 
53  Example of a profile2D histogram
54  ~~~~{.cpp}
55  {
56  TCanvas *c1 = new TCanvas("c1","Profile histogram example",200,10,700,500);
57  hprof2d = new TProfile2D("hprof2d","Profile of pz versus px and py",40,-4,4,40,-4,4,0,20);
58  Float_t px, py, pz;
59  for ( Int_t i=0; i<25000; i++) {
60  gRandom->Rannor(px,py);
61  pz = px*px + py*py;
62  hprof2d->Fill(px,py,pz,1);
63  }
64  hprof2d->Draw();
65  }
66  ~~~~
67 */
68 
69 ////////////////////////////////////////////////////////////////////////////////
70 /// Default constructor for Profile2D histograms
71 
73 {
74  fTsumwz = fTsumwz2 = 0;
75  fScaling = kFALSE;
76  BuildOptions(0,0,"");
77 }
78 
79 ////////////////////////////////////////////////////////////////////////////////
80 /// Default destructor for Profile2D histograms
81 
83 {
84 }
85 
86 ////////////////////////////////////////////////////////////////////////////////
87 /// Normal Constructor for Profile histograms -*
88 ///
89 /// The first eight parameters are similar to TH2D::TH2D.
90 /// All values of z are accepted at filling time.
91 /// To fill a profile2D histogram, one must use TProfile2D::Fill function.
92 ///
93 /// Note that when filling the profile histogram the function Fill
94 /// checks if the variable z is betyween fZmin and fZmax.
95 /// If a minimum or maximum value is set for the Z scale before filling,
96 /// then all values below zmin or above zmax will be discarded.
97 /// Setting the minimum or maximum value for the Z scale before filling
98 /// has the same effect as calling the special TProfile2D constructor below
99 /// where zmin and zmax are specified.
100 ///
101 /// H(I,J) is printed as the cell contents. The errors computed are s(I,J) if CHOPT='S'
102 /// (spread option), or e(I,J) if CHOPT=' ' (error on mean).
103 ///
104 /// See TProfile2D::BuildOptions for explanation of parameters
105 ///
106 /// see other constructors below with all possible combinations of
107 /// fix and variable bin size like in TH2D.
108 
109 TProfile2D::TProfile2D(const char *name,const char *title,Int_t nx,Double_t xlow,Double_t xup,Int_t ny,Double_t ylow,Double_t yup,Option_t *option)
110 : TH2D(name,title,nx,xlow,xup,ny,ylow,yup)
111 {
112  BuildOptions(0,0,option);
113  if (xlow >= xup || ylow >= yup) SetBuffer(fgBufferSize);
114 }
115 
116 ////////////////////////////////////////////////////////////////////////////////
117 /// Create a 2-D Profile with variable bins in X and fix bins in Y
118 
119 TProfile2D::TProfile2D(const char *name,const char *title,Int_t nx,const Double_t *xbins,Int_t ny,Double_t ylow,Double_t yup,Option_t *option)
120 : TH2D(name,title,nx,xbins,ny,ylow,yup)
121 {
122  BuildOptions(0,0,option);
123 }
124 
125 ////////////////////////////////////////////////////////////////////////////////
126 /// Create a 2-D Profile with fix bins in X and variable bins in Y
127 
128 TProfile2D::TProfile2D(const char *name,const char *title,Int_t nx,Double_t xlow,Double_t xup,Int_t ny,const Double_t *ybins,Option_t *option)
129 : TH2D(name,title,nx,xlow,xup,ny,ybins)
130 {
131  BuildOptions(0,0,option);
132 }
133 
134 ////////////////////////////////////////////////////////////////////////////////
135 /// Create a 2-D Profile with variable bins in X and variable bins in Y
136 
137 TProfile2D::TProfile2D(const char *name,const char *title,Int_t nx,const Double_t *xbins,Int_t ny,const Double_t *ybins,Option_t *option)
138 : TH2D(name,title,nx,xbins,ny,ybins)
139 {
140  BuildOptions(0,0,option);
141 }
142 
143 
144 ////////////////////////////////////////////////////////////////////////////////
145 /// Constructor for Profile2D histograms with range in z
146 /// The first eight parameters are similar to TH2D::TH2D.
147 /// Only the values of Z between ZMIN and ZMAX will be considered at filling time.
148 /// zmin and zmax will also be the maximum and minimum values
149 /// on the z scale when drawing the profile2D.
150 ///
151 /// See TProfile2D::BuildOptions for more explanations on errors
152 ///
153 
154 TProfile2D::TProfile2D(const char *name,const char *title,Int_t nx,Double_t xlow,Double_t xup,Int_t ny, Double_t ylow,Double_t yup,Double_t zlow,Double_t zup,Option_t *option)
155 : TH2D(name,title,nx,xlow,xup,ny,ylow,yup)
156 {
157  BuildOptions(zlow,zup,option);
158  if (xlow >= xup || ylow >= yup) SetBuffer(fgBufferSize);
159 }
160 
161 
162 ////////////////////////////////////////////////////////////////////////////////
163 /// Set Profile2D histogram structure and options
164 ///
165 /// zmin: minimum value allowed for z
166 /// zmax: maximum value allowed for z
167 /// if (zmin = zmax = 0) there are no limits on the allowed z values (zmin = -inf, zmax = +inf)
168 ///
169 /// option: this is the option for the computation of the t error of the profile ( TProfile2D::GetBinError )
170 /// possible values for the options are documented in TProfile2D::SetErrorOption
171 
173 {
174  // See TProfile::BuildOptions for a detailed deescription
175  //
176  //
177 
178  SetErrorOption(option);
179 
180  // create extra profile data structire (bin entries/ y^2 and sum of weight square)
182 
183  fZmin = zmin;
184  fZmax = zmax;
185  fScaling = kFALSE;
186  fTsumwz = fTsumwz2 = 0;
187 }
188 
189 ////////////////////////////////////////////////////////////////////////////////
190 /// Copy constructor.
191 
193 {
194  ((TProfile2D&)profile).Copy(*this);
195 }
196 
197 
198 ////////////////////////////////////////////////////////////////////////////////
199 /// Performs the operation: this = this + c1*f1
200 
202 {
203  Error("Add","Function not implemented for TProfile2D");
204  return kFALSE;
205 }
206 
207 
208 ////////////////////////////////////////////////////////////////////////////////
209 /// Performs the operation: this = this + c1*h1
210 
212 {
213  if (!h1) {
214  Error("Add","Attempt to add a non-existing profile");
215  return kFALSE;
216  }
217  if (!h1->InheritsFrom(TProfile2D::Class())) {
218  Error("Add","Attempt to add a non-profile2D object");
219  return kFALSE;
220  }
221 
222  return TProfileHelper::Add(this, this, h1, 1, c1);
223 }
224 
225 ////////////////////////////////////////////////////////////////////////////////
226 /// Replace contents of this profile2D by the addition of h1 and h2*
227 ///
228 /// this = c1*h1 + c2*h2
229 
231 {
232  if (!h1 || !h2) {
233  Error("Add","Attempt to add a non-existing profile");
234  return kFALSE;
235  }
236  if (!h1->InheritsFrom(TProfile2D::Class())) {
237  Error("Add","Attempt to add a non-profile2D object");
238  return kFALSE;
239  }
240  if (!h2->InheritsFrom(TProfile2D::Class())) {
241  Error("Add","Attempt to add a non-profile2D object");
242  return kFALSE;
243  }
244  return TProfileHelper::Add(this, h1, h2, c1, c2);
245 }
246 
247 
248 ////////////////////////////////////////////////////////////////////////////////
249 /// static function
250 /// set the fgApproximate flag. When the flag is true, the function GetBinError
251 /// will approximate the bin error with the average profile error on all bins
252 /// in the following situation only
253 /// - the number of bins in the profile2D is less than 10404 (eg 100x100)
254 /// - the bin number of entries is small ( <5)
255 /// - the estimated bin error is extremely small compared to the bin content
256 /// (see TProfile2D::GetBinError)
257 
259 {
260  fgApproximate = approx;
261 }
262 
263 
264 ////////////////////////////////////////////////////////////////////////////////
265 /// Fill histogram with all entries in the buffer.
266 /// - action = -1 histogram is reset and refilled from the buffer (called by THistPainter::Paint)
267 /// - action = 0 histogram is filled from the buffer
268 /// - action = 1 histogram is filled and buffer is deleted
269 /// The buffer is automatically deleted when the number of entries
270 /// in the buffer is greater than the number of entries in the histogram
271 
273 {
274  // do we need to compute the bin size?
275  if (!fBuffer) return 0;
276  Int_t nbentries = (Int_t)fBuffer[0];
277  if (!nbentries) return 0;
278  Double_t *buffer = fBuffer;
279  if (nbentries < 0) {
280  if (action == 0) return 0;
281  nbentries = -nbentries;
282  fBuffer=0;
283  Reset("ICES"); // reset without deleting the functions
284  fBuffer = buffer;
285  }
287  //find min, max of entries in buffer
288  Double_t xmin = fBuffer[2];
289  Double_t xmax = xmin;
290  Double_t ymin = fBuffer[3];
291  Double_t ymax = ymin;
292  for (Int_t i=1;i<nbentries;i++) {
293  Double_t x = fBuffer[4*i+2];
294  if (x < xmin) xmin = x;
295  if (x > xmax) xmax = x;
296  Double_t y = fBuffer[4*i+3];
297  if (y < ymin) ymin = y;
298  if (y > ymax) ymax = y;
299  }
300  if (fXaxis.GetXmax() <= fXaxis.GetXmin() || fYaxis.GetXmax() <= fYaxis.GetXmin()) {
301  THLimitsFinder::GetLimitsFinder()->FindGoodLimits(this,xmin,xmax,ymin,ymax);
302  } else {
303  fBuffer = 0;
304  Int_t keep = fBufferSize; fBufferSize = 0;
305  if (xmin < fXaxis.GetXmin()) ExtendAxis(xmin,&fXaxis);
306  if (xmax >= fXaxis.GetXmax()) ExtendAxis(xmax,&fXaxis);
307  if (ymin < fYaxis.GetXmin()) ExtendAxis(ymin,&fYaxis);
308  if (ymax >= fYaxis.GetXmax()) ExtendAxis(ymax,&fYaxis);
309  fBuffer = buffer;
310  fBufferSize = keep;
311  }
312  }
313 
314  fBuffer = 0;
315  for (Int_t i=0;i<nbentries;i++) {
316  Fill(buffer[4*i+2],buffer[4*i+3],buffer[4*i+4],buffer[4*i+1]);
317  }
318  fBuffer = buffer;
319 
320  if (action > 0) { delete [] fBuffer; fBuffer = 0; fBufferSize = 0;}
321  else {
322  if (nbentries == (Int_t)fEntries) fBuffer[0] = -nbentries;
323  else fBuffer[0] = 0;
324  }
325  return nbentries;
326 }
327 
328 ////////////////////////////////////////////////////////////////////////////////
329 /// accumulate arguments in buffer. When buffer is full, empty the buffer
330 /// fBuffer[0] = number of entries in buffer
331 /// fBuffer[1] = w of first entry
332 /// fBuffer[2] = x of first entry
333 /// fBuffer[3] = y of first entry
334 /// fBuffer[4] = z of first entry
335 
337 {
338  if (!fBuffer) return -3;
339  Int_t nbentries = (Int_t)fBuffer[0];
340  if (nbentries < 0) {
341  nbentries = -nbentries;
342  fBuffer[0] = nbentries;
343  if (fEntries > 0) {
344  Double_t *buffer = fBuffer; fBuffer=0;
345  Reset("ICES"); // reset without deleting the functions
346  fBuffer = buffer;
347  }
348  }
349  if (4*nbentries+4 >= fBufferSize) {
350  BufferEmpty(1);
351  return Fill(x,y,z,w);
352  }
353  fBuffer[4*nbentries+1] = w;
354  fBuffer[4*nbentries+2] = x;
355  fBuffer[4*nbentries+3] = y;
356  fBuffer[4*nbentries+4] = z;
357  fBuffer[0] += 1;
358  return -2;
359 }
360 
361 ////////////////////////////////////////////////////////////////////////////////
362 /// Copy a Profile2D histogram to a new profile2D histogram*
363 
365 {
366  try {
367  TProfile2D & pobj = dynamic_cast<TProfile2D&>(obj);
368 
369  TH2D::Copy(pobj);
371  fBinSumw2.Copy(pobj.fBinSumw2);
372  for (int bin=0;bin<fNcells;bin++) {
373  pobj.fArray[bin] = fArray[bin];
374  pobj.fSumw2.fArray[bin] = fSumw2.fArray[bin];
375  }
376  pobj.fZmin = fZmin;
377  pobj.fZmax = fZmax;
378  pobj.fScaling = fScaling;
379  pobj.fErrorMode = fErrorMode;
380  pobj.fTsumwz = fTsumwz;
381  pobj.fTsumwz2 = fTsumwz2;
382 
383  } catch(...) {
384  Fatal("Copy","Cannot copy a TProfile2D in a %s",obj.IsA()->GetName());
385  }
386 
387 }
388 
389 
390 ////////////////////////////////////////////////////////////////////////////////
391 /// Performs the operation: this = this/(c1*f1)
392 /// This function is not implemented
393 
395 {
396  Error("Divide","Function not implemented for TProfile2D");
397  return kFALSE;
398 }
399 
400 ////////////////////////////////////////////////////////////////////////////////
401 /// Divide this profile2D by h1
402 ///
403 /// this = this/h1
404 ///
405 /// This function return kFALSE if the divide operation failed
406 
408 {
409 
410  if (!h1) {
411  Error("Divide","Attempt to divide a non-existing profile2D");
412  return kFALSE;
413  }
414  if (!h1->InheritsFrom(TProfile2D::Class())) {
415  Error("Divide","Attempt to divide a non-profile2D object");
416  return kFALSE;
417  }
418  TProfile2D *p1 = (TProfile2D*)h1;
419 
420  // delete buffer if it is there since it will become invalid
421  if (fBuffer) BufferEmpty(1);
422 
423  //*-*- Check profile compatibility
424  Int_t nx = GetNbinsX();
425  if (nx != p1->GetNbinsX()) {
426  Error("Divide","Attempt to divide profiles with different number of bins");
427  return kFALSE;
428  }
429  Int_t ny = GetNbinsY();
430  if (ny != p1->GetNbinsY()) {
431  Error("Divide","Attempt to divide profiles with different number of bins");
432  return kFALSE;
433  }
434 
435  //*-*- Reset statistics
436  fEntries = fTsumw = fTsumw2 = fTsumwx = fTsumwx2 = 0;
437 
438  //*-*- Loop on bins (including underflows/overflows)
439  Int_t bin,binx,biny;
440  Double_t *cu1 = p1->GetW();
441  Double_t *er1 = p1->GetW2();
442  Double_t *en1 = p1->GetB();
443  Double_t c0,c1,w,z,x,y;
444  for (binx =0;binx<=nx+1;binx++) {
445  for (biny =0;biny<=ny+1;biny++) {
446  bin = biny*(fXaxis.GetNbins()+2) + binx;
447  c0 = fArray[bin];
448  c1 = cu1[bin];
449  if (c1) w = c0/c1;
450  else w = 0;
451  fArray[bin] = w;
452  z = TMath::Abs(w);
453  x = fXaxis.GetBinCenter(binx);
454  y = fYaxis.GetBinCenter(biny);
455  fEntries++;
456  fTsumw += z;
457  fTsumw2 += z*z;
458  fTsumwx += z*x;
459  fTsumwx2 += z*x*x;
460  fTsumwy += z*y;
461  fTsumwy2 += z*y*y;
462  fTsumwxy += z*x*y;
463  fTsumwz += z;
464  fTsumwz2 += z*z;
465  Double_t e0 = fSumw2.fArray[bin];
466  Double_t e1 = er1[bin];
467  Double_t c12= c1*c1;
468  if (!c1) fSumw2.fArray[bin] = 0;
469  else fSumw2.fArray[bin] = (e0*c1*c1 + e1*c0*c0)/(c12*c12);
470  if (!en1[bin]) fBinEntries.fArray[bin] = 0;
471  else fBinEntries.fArray[bin] /= en1[bin];
472  }
473  }
474  // mantaining the correct sum of weights square is not supported when dividing
475  // bin error resulting from division of profile needs to be checked
476  if (fBinSumw2.fN) {
477  Warning("Divide","Cannot preserve during the division of profiles the sum of bin weight square");
478  fBinSumw2 = TArrayD();
479  }
480  return kTRUE;
481 }
482 
483 
484 ////////////////////////////////////////////////////////////////////////////////
485 /// Replace contents of this profile2D by the division of h1 by h2
486 ///
487 /// this = c1*h1/(c2*h2)
488 ///
489 /// This function return kFALSE if the divide operation failed
490 
492 {
493  TString opt = option;
494  opt.ToLower();
495  Bool_t binomial = kFALSE;
496  if (opt.Contains("b")) binomial = kTRUE;
497  if (!h1 || !h2) {
498  Error("Divide","Attempt to divide a non-existing profile2D");
499  return kFALSE;
500  }
501  if (!h1->InheritsFrom(TProfile2D::Class())) {
502  Error("Divide","Attempt to divide a non-profile2D object");
503  return kFALSE;
504  }
505  TProfile2D *p1 = (TProfile2D*)h1;
506  if (!h2->InheritsFrom(TProfile2D::Class())) {
507  Error("Divide","Attempt to divide a non-profile2D object");
508  return kFALSE;
509  }
510  TProfile2D *p2 = (TProfile2D*)h2;
511 
512  // delete buffer if it is there since it will become invalid
513  if (fBuffer) BufferEmpty(1);
514 
515  //*-*- Check histogram compatibility
516  Int_t nx = GetNbinsX();
517  if (nx != p1->GetNbinsX() || nx != p2->GetNbinsX()) {
518  Error("Divide","Attempt to divide profiles with different number of bins");
519  return kFALSE;
520  }
521  Int_t ny = GetNbinsY();
522  if (ny != p1->GetNbinsY() || ny != p2->GetNbinsY()) {
523  Error("Divide","Attempt to divide profiles with different number of bins");
524  return kFALSE;
525  }
526  if (!c2) {
527  Error("Divide","Coefficient of dividing profile cannot be zero");
528  return kFALSE;
529  }
530 
531  //*-*- Reset statistics
532  fEntries = fTsumw = fTsumw2 = fTsumwx = fTsumwx2 = 0;
533 
534  //*-*- Loop on bins (including underflows/overflows)
535  Int_t bin,binx,biny;
536  Double_t *cu1 = p1->GetW();
537  Double_t *cu2 = p2->GetW();
538  Double_t *er1 = p1->GetW2();
539  Double_t *er2 = p2->GetW2();
540  Double_t *en1 = p1->GetB();
541  Double_t *en2 = p2->GetB();
542  Double_t b1,b2,w,z,x,y,ac1,ac2;
543  ac1 = TMath::Abs(c1);
544  ac2 = TMath::Abs(c2);
545  for (binx =0;binx<=nx+1;binx++) {
546  for (biny =0;biny<=ny+1;biny++) {
547  bin = biny*(fXaxis.GetNbins()+2) + binx;
548  b1 = cu1[bin];
549  b2 = cu2[bin];
550  if (b2) w = c1*b1/(c2*b2);
551  else w = 0;
552  fArray[bin] = w;
553  z = TMath::Abs(w);
554  x = fXaxis.GetBinCenter(binx);
555  y = fYaxis.GetBinCenter(biny);
556  fEntries++;
557  fTsumw += z;
558  fTsumw2 += z*z;
559  fTsumwx += z*x;
560  fTsumwx2 += z*x*x;
561  fTsumwy += z*y;
562  fTsumwy2 += z*y*y;
563  fTsumwxy += z*x*y;
564  fTsumwz += z;
565  fTsumwz2 += z*z;
566  Double_t e1 = er1[bin];
567  Double_t e2 = er2[bin];
568  //Double_t b22= b2*b2*d2;
569  Double_t b22= b2*b2*TMath::Abs(c2);
570  if (!b2) fSumw2.fArray[bin] = 0;
571  else {
572  if (binomial) {
573  fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/(c2*b2));
574  } else {
575  fSumw2.fArray[bin] = ac1*ac2*(e1*b2*b2 + e2*b1*b1)/(b22*b22);
576  }
577  }
578  if (!en2[bin]) fBinEntries.fArray[bin] = 0;
579  else fBinEntries.fArray[bin] = en1[bin]/en2[bin];
580  }
581  }
582  return kTRUE;
583 }
584 
585 ////////////////////////////////////////////////////////////////////////////////
586 /// Fill a Profile2D histogram (no weights)
587 ///*-* =======================================
588 
590 {
591  if (fBuffer) return BufferFill(x,y,z,1);
592 
593  Int_t bin,binx,biny;
594 
595  if (fZmin != fZmax) {
596  if (z <fZmin || z> fZmax || TMath::IsNaN(z) ) return -1;
597  }
598 
599  fEntries++;
600  binx =fXaxis.FindBin(x);
601  biny =fYaxis.FindBin(y);
602  if (binx <0 || biny <0) return -1;
603  bin = GetBin(binx, biny);
604  fArray[bin] += z;
605  fSumw2.fArray[bin] += z*z;
606  fBinEntries.fArray[bin] += 1;
607  if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1;
608  if (binx == 0 || binx > fXaxis.GetNbins()) {
609  if (!fgStatOverflows) return -1;
610  }
611  if (biny == 0 || biny > fYaxis.GetNbins()) {
612  if (!fgStatOverflows) return -1;
613  }
614  ++fTsumw;
615  ++fTsumw2;
616  fTsumwx += x;
617  fTsumwx2 += x*x;
618  fTsumwy += y;
619  fTsumwy2 += y*y;
620  fTsumwxy += x*y;
621  fTsumwz += z;
622  fTsumwz2 += z*z;
623  return bin;
624 }
625 
626 ////////////////////////////////////////////////////////////////////////////////
627 /// Fill a Profile2D histogram (no weights)
628 ///
629 
630 Int_t TProfile2D::Fill(Double_t x, const char *namey, Double_t z)
631 {
632  Int_t bin,binx,biny;
633 
634  if (fZmin != fZmax) {
635  if (z <fZmin || z> fZmax || TMath::IsNaN(z)) return -1;
636  }
637 
638  fEntries++;
639  binx =fXaxis.FindBin(x);
640  biny =fYaxis.FindBin(namey);
641  if (binx <0 || biny <0) return -1;
642  bin = biny*(fXaxis.GetNbins()+2) + binx;
643  AddBinContent(bin, z);
644  fSumw2.fArray[bin] += (Double_t)z*z;
645  fBinEntries.fArray[bin] += 1;
646  if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1;
647  if (binx == 0 || binx > fXaxis.GetNbins()) {
648  if (!fgStatOverflows) return -1;
649  }
650  if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
651  Double_t y = fYaxis.GetBinCenter(biny);
652  ++fTsumw;
653  ++fTsumw2;
654  fTsumwx += x;
655  fTsumwx2 += x*x;
656  fTsumwy += y;
657  fTsumwy2 += y*y;
658  fTsumwxy += x*y;
659  fTsumwz += z;
660  fTsumwz2 += z*z;
661  return bin;
662 }
663 
664 ////////////////////////////////////////////////////////////////////////////////
665 /// Fill a Profile2D histogram (no weights)
666 ///
667 
668 Int_t TProfile2D::Fill(const char *namex, const char *namey, Double_t z)
669 {
670  Int_t bin,binx,biny;
671 
672  if (fZmin != fZmax) {
673  if (z <fZmin || z> fZmax || TMath::IsNaN(z) ) return -1;
674  }
675 
676  fEntries++;
677  binx =fXaxis.FindBin(namex);
678  biny =fYaxis.FindBin(namey);
679  if (binx <0 || biny <0) return -1;
680  bin = biny*(fXaxis.GetNbins()+2) + binx;
681  AddBinContent(bin, z);
682  fSumw2.fArray[bin] += (Double_t)z*z;
683  fBinEntries.fArray[bin] += 1;
684  if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1;
685  if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
686  if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
687  Double_t x = fYaxis.GetBinCenter(binx);
688  Double_t y = fYaxis.GetBinCenter(biny);
689  ++fTsumw;
690  ++fTsumw2;
691  fTsumwx += x;
692  fTsumwx2 += x*x;
693  fTsumwy += y;
694  fTsumwy2 += y*y;
695  fTsumwxy += x*y;
696  fTsumwz += z;
697  fTsumwz2 += z*z;
698  return bin;
699 }
700 
701 ////////////////////////////////////////////////////////////////////////////////
702 /// Fill a Profile2D histogram (no weights)
703 ///
704 
705 Int_t TProfile2D::Fill(const char *namex, Double_t y, Double_t z)
706 {
707  Int_t bin,binx,biny;
708 
709  if (fZmin != fZmax) {
710  if (z <fZmin || z> fZmax || TMath::IsNaN(z)) return -1;
711  }
712 
713  fEntries++;
714  binx =fXaxis.FindBin(namex);
715  biny =fYaxis.FindBin(y);
716  if (binx <0 || biny <0) return -1;
717  bin = biny*(fXaxis.GetNbins()+2) + binx;
718  AddBinContent(bin, z);
719  fSumw2.fArray[bin] += (Double_t)z*z;
720  fBinEntries.fArray[bin] += 1;
721  if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1;
722  if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
723  if (biny == 0 || biny > fYaxis.GetNbins()) {
724  if (!fgStatOverflows) return -1;
725  }
726  Double_t x = fYaxis.GetBinCenter(binx);
727  ++fTsumw;
728  ++fTsumw2;
729  fTsumwx += x;
730  fTsumwx2 += x*x;
731  fTsumwy += y;
732  fTsumwy2 += y*y;
733  fTsumwxy += x*y;
734  fTsumwz += z;
735  fTsumwz2 += z*z;
736  return bin;
737 }
738 
739 ////////////////////////////////////////////////////////////////////////////////
740 /// Fill a Profile2D histogram with weights
741 
743 {
744  if (fBuffer) return BufferFill(x,y,z,w);
745 
746  Int_t bin,binx,biny;
747 
748  if (fZmin != fZmax) {
749  if (z <fZmin || z> fZmax || TMath::IsNaN(z)) return -1;
750  }
751 
752  Double_t u= w;
753  fEntries++;
754  binx =fXaxis.FindBin(x);
755  biny =fYaxis.FindBin(y);
756  if (binx <0 || biny <0) return -1;
757  bin = biny*(fXaxis.GetNbins()+2) + binx;
758  AddBinContent(bin, u*z);
759  fSumw2.fArray[bin] += u*z*z;
760  if (!fBinSumw2.fN && u != 1.0 && !TestBit(TH1::kIsNotW)) Sumw2(); // must be called before accumulating the entries
761  if (fBinSumw2.fN) fBinSumw2.fArray[bin] += u*u;
762  fBinEntries.fArray[bin] += u;
763  if (binx == 0 || binx > fXaxis.GetNbins()) {
764  if (!fgStatOverflows) return -1;
765  }
766  if (biny == 0 || biny > fYaxis.GetNbins()) {
767  if (!fgStatOverflows) return -1;
768  }
769  fTsumw += u;
770  fTsumw2 += u*u;
771  fTsumwx += u*x;
772  fTsumwx2 += u*x*x;
773  fTsumwy += u*y;
774  fTsumwy2 += u*y*y;
775  fTsumwxy += u*x*y;
776  fTsumwz += u*z;
777  fTsumwz2 += u*z*z;
778  return bin;
779 }
780 
781 ////////////////////////////////////////////////////////////////////////////////
782 /// Return bin content of a Profile2D histogram
783 
785 {
786  if (fBuffer) ((TProfile2D*)this)->BufferEmpty();
787 
788  if (bin < 0 || bin >= fNcells) return 0;
789  if (fBinEntries.fArray[bin] == 0) return 0;
790  if (!fArray) return 0;
791  return fArray[bin]/fBinEntries.fArray[bin];
792 }
793 
794 ////////////////////////////////////////////////////////////////////////////////
795 /// Return bin entries of a Profile2D histogram
796 
798 {
799  if (fBuffer) ((TProfile2D*)this)->BufferEmpty();
800 
801  if (bin < 0 || bin >= fNcells) return 0;
802  return fBinEntries.fArray[bin];
803 }
804 
805 ////////////////////////////////////////////////////////////////////////////////
806 /// Return bin effective entries for a weighted filled Profile histogram.
807 /// In case of an unweighted profile, it is equivalent to the number of entries per bin
808 /// The effective entries is defined as the square of the sum of the weights divided by the
809 /// sum of the weights square.
810 /// TProfile::Sumw2() must be called before filling the profile with weights.
811 /// Only by calling this method the sum of the square of the weights per bin is stored.
812 ///
813 ///*-* =========================================
814 
816 {
817  return TProfileHelper::GetBinEffectiveEntries(this, bin);
818 }
819 
820 ////////////////////////////////////////////////////////////////////////////////
821 /// -*Return bin error of a Profile2D histogram
822 ///
823 /// Computing errors: A moving field
824 /// =================================
825 /// The computation of errors for a TProfile2D has evolved with the versions
826 /// of ROOT. The difficulty is in computing errors for bins with low statistics.
827 /// - prior to version 3.10, we had no special treatment of low statistic bins.
828 /// As a result, these bins had huge errors. The reason is that the
829 /// expression eprim2 is very close to 0 (rounding problems) or 0.
830 /// - The algorithm is modified/protected for the case
831 /// when a TProfile2D is projected (ProjectionX). The previous algorithm
832 /// generated a N^2 problem when projecting a TProfile2D with a large number of
833 /// bins (eg 100000).
834 /// - in version 3.10/02, a new static function TProfile::Approximate
835 /// is introduced to enable or disable (default) the approximation.
836 /// (see also comments in TProfile::GetBinError)
837 
839 {
840  return TProfileHelper::GetBinError((TProfile2D*)this, bin);
841 }
842 
843 ////////////////////////////////////////////////////////////////////////////////
844 /// -*Return option to compute profile2D errors
845 ///*-* =========================================
846 
848 {
849  if (fErrorMode == kERRORSPREAD) return "s";
850  if (fErrorMode == kERRORSPREADI) return "i";
851  if (fErrorMode == kERRORSPREADG) return "g";
852  return "";
853 }
854 
855 ////////////////////////////////////////////////////////////////////////////////
856 /// fill the array stats from the contents of this profile
857 /// The array stats must be correctly dimensionned in the calling program.
858 /// stats[0] = sumw
859 /// stats[1] = sumw2
860 /// stats[2] = sumwx
861 /// stats[3] = sumwx2
862 /// stats[4] = sumwy
863 /// stats[5] = sumwy2
864 /// stats[6] = sumwxy
865 /// stats[7] = sumwz
866 /// stats[8] = sumwz2
867 ///
868 /// If no axis-subrange is specified (via TAxis::SetRange), the array stats
869 /// is simply a copy of the statistics quantities computed at filling time.
870 /// If a sub-range is specified, the function recomputes these quantities
871 /// from the bin contents in the current axis range.
872 
873 void TProfile2D::GetStats(Double_t *stats) const
874 {
875  if (fBuffer) ((TProfile2D*)this)->BufferEmpty();
876 
877  // Loop on bins
879  Int_t bin, binx, biny;
880  Double_t w, w2;
881  Double_t x,y;
882  for (bin=0;bin<9;bin++) stats[bin] = 0;
883  if (!fBinEntries.fArray) return;
884  Int_t firstBinX = fXaxis.GetFirst();
885  Int_t lastBinX = fXaxis.GetLast();
886  Int_t firstBinY = fYaxis.GetFirst();
887  Int_t lastBinY = fYaxis.GetLast();
888  // include underflow/overflow if TH1::StatOverflows(kTRUE) in case no range is set on the axis
889  if (fgStatOverflows) {
890  if ( !fXaxis.TestBit(TAxis::kAxisRange) ) {
891  if (firstBinX == 1) firstBinX = 0;
892  if (lastBinX == fXaxis.GetNbins() ) lastBinX += 1;
893  }
894  if ( !fYaxis.TestBit(TAxis::kAxisRange) ) {
895  if (firstBinY == 1) firstBinY = 0;
896  if (lastBinY == fYaxis.GetNbins() ) lastBinY += 1;
897  }
898  }
899  for (biny = firstBinY; biny <= lastBinY; biny++) {
900  y = fYaxis.GetBinCenter(biny);
901  for (binx = firstBinX; binx <= lastBinX; binx++) {
902  bin = GetBin(binx,biny);
903  w = fBinEntries.fArray[bin];
904  w2 = (fBinSumw2.fN ? fBinSumw2.fArray[bin] : w );
905  x = fXaxis.GetBinCenter(binx);
906  stats[0] += w;
907  stats[1] += w2;
908  stats[2] += w*x;
909  stats[3] += w*x*x;
910  stats[4] += w*y;
911  stats[5] += w*y*y;
912  stats[6] += w*x*y;
913  stats[7] += fArray[bin];
914  stats[8] += fSumw2.fArray[bin];
915  }
916  }
917  } else {
918  stats[0] = fTsumw;
919  stats[1] = fTsumw2;
920  stats[2] = fTsumwx;
921  stats[3] = fTsumwx2;
922  stats[4] = fTsumwy;
923  stats[5] = fTsumwy2;
924  stats[6] = fTsumwxy;
925  stats[7] = fTsumwz;
926  stats[8] = fTsumwz2;
927  }
928 }
929 
930 ////////////////////////////////////////////////////////////////////////////////
931 /// Reduce the number of bins for this axis to the number of bins having a label.
932 
934 {
936 }
937 
938 ////////////////////////////////////////////////////////////////////////////////
939 /// Double the number of bins for axis.
940 /// Refill histogram
941 /// This function is called by TAxis::FindBin(const char *label)
942 
944 {
946 }
947 
948 ////////////////////////////////////////////////////////////////////////////////
949 /// Set option(s) to draw axis with labels
950 /// option = "a" sort by alphabetic order
951 /// = ">" sort by decreasing values
952 /// = "<" sort by increasing values
953 /// = "h" draw labels horizonthal
954 /// = "v" draw labels vertical
955 /// = "u" draw labels up (end of label right adjusted)
956 /// = "d" draw labels down (start of label left adjusted)
957 
959 {
960 
961  TAxis *axis = GetXaxis();
962  if (ax[0] == 'y' || ax[0] == 'Y') axis = GetYaxis();
963  THashList *labels = axis->GetLabels();
964  if (!labels) {
965  Warning("LabelsOption","Cannot sort. No labels");
966  return;
967  }
968  TString opt = option;
969  opt.ToLower();
970  if (opt.Contains("h")) {
971  axis->SetBit(TAxis::kLabelsHori);
974  axis->ResetBit(TAxis::kLabelsUp);
975  }
976  if (opt.Contains("v")) {
977  axis->SetBit(TAxis::kLabelsVert);
980  axis->ResetBit(TAxis::kLabelsUp);
981  }
982  if (opt.Contains("u")) {
983  axis->SetBit(TAxis::kLabelsUp);
987  }
988  if (opt.Contains("d")) {
989  axis->SetBit(TAxis::kLabelsDown);
992  axis->ResetBit(TAxis::kLabelsUp);
993  }
994  Int_t sort = -1;
995  if (opt.Contains("a")) sort = 0;
996  if (opt.Contains(">")) sort = 1;
997  if (opt.Contains("<")) sort = 2;
998  if (sort < 0) return;
999 
1000  Int_t nx = fXaxis.GetNbins()+2;
1001  Int_t ny = fYaxis.GetNbins()+2;
1002  Int_t n = TMath::Min(axis->GetNbins(), labels->GetSize());
1003  Int_t *a = new Int_t[n+2];
1004  Int_t i,j,k,bin;
1005  Double_t *sumw = new Double_t[nx*ny];
1006  Double_t *errors = new Double_t[nx*ny];
1007  Double_t *ent = new Double_t[nx*ny];
1008  THashList *labold = new THashList(labels->GetSize(),1);
1009  TIter nextold(labels);
1010  TObject *obj;
1011  while ((obj=nextold())) {
1012  labold->Add(obj);
1013  }
1014  labels->Clear();
1015  if (sort > 0) {
1016  //---sort by values of bins
1017  Double_t *pcont = new Double_t[n+2];
1018  for (i=0;i<=n;i++) pcont[i] = 0;
1019  for (i=1;i<nx;i++) {
1020  for (j=1;j<ny;j++) {
1021  bin = i+nx*j;
1022  sumw[bin] = fArray[bin];
1023  errors[bin] = fSumw2.fArray[bin];
1024  ent[bin] = fBinEntries.fArray[bin];
1025  if (axis == GetXaxis()) k = i;
1026  else k = j;
1027  if (fBinEntries.fArray[bin] != 0) pcont[k-1] += fArray[bin]/fBinEntries.fArray[bin];
1028  }
1029  }
1030  if (sort ==1) TMath::Sort(n,pcont,a,kTRUE); //sort by decreasing values
1031  else TMath::Sort(n,pcont,a,kFALSE); //sort by increasing values
1032  delete [] pcont;
1033  for (i=0;i<n;i++) {
1034  obj = labold->At(a[i]);
1035  labels->Add(obj);
1036  obj->SetUniqueID(i+1);
1037  }
1038  for (i=1;i<nx;i++) {
1039  for (j=1;j<ny;j++) {
1040  bin = i+nx*j;
1041  if (axis == GetXaxis()) {
1042  fArray[bin] = sumw[a[i-1]+1+nx*j];
1043  fSumw2.fArray[bin] = errors[a[i-1]+1+nx*j];
1044  fBinEntries.fArray[bin] = ent[a[i-1]+1+nx*j];
1045  } else {
1046  fArray[bin] = sumw[i+nx*(a[j-1]+1)];
1047  fSumw2.fArray[bin] = errors[i+nx*(a[j-1]+1)];
1048  fBinEntries.fArray[bin] = ent[i+nx*(a[j-1]+1)];
1049  }
1050  }
1051  }
1052  } else {
1053  //---alphabetic sort
1054  const UInt_t kUsed = 1<<18;
1055  TObject *objk=0;
1056  a[0] = 0;
1057  a[n+1] = n+1;
1058  for (i=1;i<=n;i++) {
1059  const char *label = "zzzzzzzzzzzz";
1060  for (j=1;j<=n;j++) {
1061  obj = labold->At(j-1);
1062  if (!obj) continue;
1063  if (obj->TestBit(kUsed)) continue;
1064  //use strcasecmp for case non-sensitive sort (may be an option)
1065  if (strcmp(label,obj->GetName()) < 0) continue;
1066  objk = obj;
1067  a[i] = j;
1068  label = obj->GetName();
1069  }
1070  if (objk) {
1071  objk->SetUniqueID(i);
1072  labels->Add(objk);
1073  objk->SetBit(kUsed);
1074  }
1075  }
1076  for (i=1;i<=n;i++) {
1077  obj = labels->At(i-1);
1078  if (!obj) continue;
1079  obj->ResetBit(kUsed);
1080  }
1081  for (i=0;i<nx;i++) {
1082  for (j=0;j<ny;j++) {
1083  bin = i+nx*j;
1084  sumw[bin] = fArray[bin];
1085  errors[bin] = fSumw2.fArray[bin];
1086  ent[bin] = fBinEntries.fArray[bin];
1087  }
1088  }
1089  for (i=0;i<nx;i++) {
1090  for (j=0;j<ny;j++) {
1091  bin = i+nx*j;
1092  if (axis == GetXaxis()) {
1093  fArray[bin] = sumw[a[i]+nx*j];
1094  fSumw2.fArray[bin] = errors[a[i]+nx*j];
1095  fBinEntries.fArray[bin] = ent[a[i]+nx*j];
1096  } else {
1097  fArray[bin] = sumw[i+nx*a[j]];
1098  fSumw2.fArray[bin] = errors[i+nx*a[j]];
1099  fBinEntries.fArray[bin] = ent[i+nx*a[j]];
1100  }
1101  }
1102  }
1103  }
1104  delete labold;
1105  if (a) delete [] a;
1106  if (sumw) delete [] sumw;
1107  if (errors) delete [] errors;
1108  if (ent) delete [] ent;
1109 }
1110 
1111 ////////////////////////////////////////////////////////////////////////////////
1112 ///Merge all histograms in the collection in this histogram.
1113 ///This function computes the min/max for the axes,
1114 ///compute a new number of bins, if necessary,
1115 ///add bin contents, errors and statistics.
1116 ///If overflows are present and limits are different the function will fail.
1117 ///The function returns the total number of entries in the result histogram
1118 ///if the merge is successfull, -1 otherwise.
1119 ///
1120 ///IMPORTANT remark. The 2 axis x and y may have different number
1121 ///of bins and different limits, BUT the largest bin width must be
1122 ///a multiple of the smallest bin width and the upper limit must also
1123 ///be a multiple of the bin width.
1124 
1126 {
1127  return TProfileHelper::Merge(this, li);
1128 }
1129 
1130 ////////////////////////////////////////////////////////////////////////////////
1131 /// Performs the operation: this = this*c1*f1
1132 
1134 {
1135  Error("Multiply","Function not implemented for TProfile2D");
1136  return kFALSE;
1137 }
1138 
1139 ////////////////////////////////////////////////////////////////////////////////
1140 /// Multiply this profile2D by h1 -
1141 ///*-* =============================
1142 ///
1143 /// this = this*h1
1144 ///
1145 
1147 {
1148  Error("Multiply","Multiplication of profile2D histograms not implemented");
1149  return kFALSE;
1150 }
1151 
1152 
1153 ////////////////////////////////////////////////////////////////////////////////
1154 ///*-*Replace contents of this profile2D by multiplication of h1 by h2*-*
1155 ///*-* ================================================================
1156 ///
1157 /// this = (c1*h1)*(c2*h2)
1158 ///
1159 
1161 {
1162  Error("Multiply","Multiplication of profile2D histograms not implemented");
1163  return kFALSE;
1164 }
1165 
1166 ////////////////////////////////////////////////////////////////////////////////
1167 ///*-*Project this profile2D into a 2-D histogram along X,Y -*
1168 ///*-* =====================================================
1169 ///
1170 /// The projection is always of the type TH2D.
1171 ///
1172 /// if option "E" is specified the errors of the projected histogram are computed and set
1173 /// to be equal to the errors of the profile.
1174 /// Option "E" is defined as the default one in the header file.
1175 /// if option "" is specified the histogram errors are simply the sqrt of its content
1176 /// if option "B" is specified, the content of bin of the returned histogram
1177 /// will be equal to the GetBinEntries(bin) of the profile,
1178 /// if option "C=E" the bin contents of the projection are set to the
1179 /// bin errors of the profile
1180 /// if option "W" is specified the bin content of the projected histogram is set to the
1181 /// product of the bin content of the profile and the entries.
1182 /// With this option the returned histogram will be equivalent to the one obtained by
1183 /// filling directly a TH2D using the 3-rd value as a weight.
1184 /// This option makes sense only for profile filled with all weights =1.
1185 /// When the profile is weighted (filled with weights different than 1) the
1186 /// bin error of the projected histogram (obtained using this option "W") cannot be
1187 /// correctly computed from the information stored in the profile. In that case the
1188 /// obtained histogram contains as bin error square the weighted sum of the square of the
1189 /// profiled observable (TProfile2D::fSumw2[bin] )
1190 
1191 TH2D *TProfile2D::ProjectionXY(const char *name, Option_t *option) const
1192 {
1193 
1194  TString opt = option;
1195  opt.ToLower();
1196 
1197  // Create the projection histogram
1198  // name of projected histogram is by default name of orginal histogram + _pxy
1199  TString pname(name);
1200  if (pname.IsNull() || pname == "_pxy")
1201  pname = TString(GetName() ) + TString("_pxy");
1202 
1203 
1204  Int_t nx = fXaxis.GetNbins();
1205  Int_t ny = fYaxis.GetNbins();
1206  const TArrayD *xbins = fXaxis.GetXbins();
1207  const TArrayD *ybins = fYaxis.GetXbins();
1208  TH2D * h1 = 0;
1209  if (xbins->fN == 0 && ybins->fN == 0) {
1210  h1 = new TH2D(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax(),ny,fYaxis.GetXmin(),fYaxis.GetXmax());
1211  } else if (xbins->fN == 0) {
1212  h1 = new TH2D(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax(),ny, ybins->GetArray() );
1213  } else if (ybins->fN == 0) {
1214  h1 = new TH2D(pname,GetTitle(),nx,xbins->GetArray(),ny,fYaxis.GetXmin(),fYaxis.GetXmax());
1215  } else {
1216  h1 = new TH2D(pname,GetTitle(),nx,xbins->GetArray(),ny,ybins->GetArray() );
1217  }
1218  Bool_t computeErrors = kFALSE;
1219  Bool_t cequalErrors = kFALSE;
1220  Bool_t binEntries = kFALSE;
1221  Bool_t binWeight = kFALSE;
1222  if (opt.Contains("b")) binEntries = kTRUE;
1223  if (opt.Contains("e")) computeErrors = kTRUE;
1224  if (opt.Contains("w")) binWeight = kTRUE;
1225  if (opt.Contains("c=e")) {cequalErrors = kTRUE; computeErrors=kFALSE;}
1226  if (computeErrors || binWeight || (binEntries && fBinSumw2.fN) ) h1->Sumw2();
1227 
1228  // Fill the projected histogram
1229  Int_t bin,binx, biny;
1230  Double_t cont;
1231  for (binx =0;binx<=nx+1;binx++) {
1232  for (biny =0;biny<=ny+1;biny++) {
1233  bin = GetBin(binx,biny);
1234 
1235  if (binEntries) cont = GetBinEntries(bin);
1236  else if (cequalErrors) cont = GetBinError(bin);
1237  else if (binWeight) cont = GetBinContent(bin) * GetBinEntries(bin);
1238  else cont = GetBinContent(bin); // default case
1239 
1240  h1->SetBinContent(bin ,cont);
1241 
1242  // if option E projected histogram errors are same as profile
1243  if (computeErrors ) h1->SetBinError(bin , GetBinError(bin) );
1244  // in case of option W bin error is deduced from bin sum of z**2 values of profile
1245  // this is correct only if the profile is unweighted
1246  if (binWeight) h1->GetSumw2()->fArray[bin] = fSumw2.fArray[bin];
1247  // in case of bin entries and profile is weighted, we need to set also the bin error
1248  if (binEntries && fBinSumw2.fN ) {
1249  R__ASSERT( h1->GetSumw2() );
1250  h1->GetSumw2()->fArray[bin] = fBinSumw2.fArray[bin];
1251  }
1252  }
1253  }
1254  h1->SetEntries(fEntries);
1255  return h1;
1256 }
1257 
1258 ////////////////////////////////////////////////////////////////////////////////
1259 /// *-*Project a 2-D histogram into a profile histogram along X
1260 /// *-* ========================================================
1261 ///
1262 /// The projection is made from the channels along the Y axis
1263 /// ranging from firstybin to lastybin included.
1264 /// The result is a 1D profile which contains the combination of all the considered bins along Y
1265 /// By default, bins 1 to ny are included
1266 /// When all bins are included, the number of entries in the projection
1267 /// is set to the number of entries of the 2-D histogram, otherwise
1268 /// the number of entries is incremented by 1 for all non empty cells.
1269 ///
1270 /// The option can also be used to specify the projected profile error type.
1271 /// Values which can be used are 's', 'i', or 'g'. See TProfile::BuildOptions for details
1272 ///
1273 ///
1274 
1275 TProfile *TProfile2D::ProfileX(const char *name, Int_t firstybin, Int_t lastybin, Option_t *option) const
1276 {
1277  return DoProfile(true, name, firstybin, lastybin, option);
1278 }
1279 
1280 ////////////////////////////////////////////////////////////////////////////////
1281 /// *-*Project a 2-D histogram into a profile histogram along X
1282 /// *-* ========================================================
1283 ///
1284 /// The projection is made from the channels along the X axis
1285 /// ranging from firstybin to lastybin included.
1286 /// The result is a 1D profile which contains the combination of all the considered bins along X
1287 /// By default, bins 1 to ny are included
1288 /// When all bins are included, the number of entries in the projection
1289 /// is set to the number of entries of the 2-D histogram, otherwise
1290 /// the number of entries is incremented by 1 for all non empty cells.
1291 ///
1292 /// The option can also be used to specify the projected profile error type.
1293 /// Values which can be used are 's', 'i', or 'g'. See TProfile::BuildOptions for details
1294 ///
1295 ///
1296 ///
1297 
1298 TProfile *TProfile2D::ProfileY(const char *name, Int_t firstxbin, Int_t lastxbin, Option_t *option) const
1299 {
1300  return DoProfile(false, name, firstxbin, lastxbin, option);
1301 }
1302 
1303 ////////////////////////////////////////////////////////////////////////////////
1304 /// implementation of ProfileX or ProfileY for a TProfile2D
1305 /// Do correctly the combination of the bin averages when doing the projection
1306 
1307 TProfile * TProfile2D::DoProfile(bool onX, const char *name, Int_t firstbin, Int_t lastbin, Option_t *option) const {
1308  TString opt = option;
1309  opt.ToLower();
1310  bool originalRange = opt.Contains("o");
1311 
1312  TString expectedName = ( onX ? "_pfx" : "_pfy" );
1313 
1314  TString pname(name);
1315  if (pname.IsNull() || name == expectedName)
1316  pname = TString(GetName() ) + expectedName;
1317 
1318  const TAxis& outAxis = ( onX ? fXaxis : fYaxis );
1319  const TArrayD *bins = outAxis.GetXbins();
1320  Int_t firstOutBin = outAxis.GetFirst();
1321  Int_t lastOutBin = outAxis.GetLast();
1322 
1323  TProfile * p1 = 0;
1324  // case of fixed bins
1325  if (bins->fN == 0) {
1326  if (originalRange)
1327  p1 = new TProfile(pname,GetTitle(), outAxis.GetNbins(), outAxis.GetXmin(), outAxis.GetXmax(), opt );
1328  else
1329  p1 = new TProfile(pname,GetTitle(), lastOutBin-firstOutBin+1,
1330  outAxis.GetBinLowEdge(firstOutBin),outAxis.GetBinUpEdge(lastOutBin), opt);
1331  } else {
1332  // case of variable bins
1333  if (originalRange )
1334  p1 = new TProfile(pname,GetTitle(),outAxis.GetNbins(),bins->fArray,opt);
1335  else
1336  p1 = new TProfile(pname,GetTitle(),lastOutBin-firstOutBin+1,&bins->fArray[firstOutBin-1],opt);
1337 
1338  }
1339 
1340  if (fBinSumw2.fN) p1->Sumw2();
1341 
1342  // make projection in a 2D first
1343  TH2D * h2dW = ProjectionXY("h2temp-W","W");
1344  TH2D * h2dN = ProjectionXY("h2temp-N","B");
1345 
1346  h2dW->SetDirectory(0); h2dN->SetDirectory(0);
1347 
1348 
1349  TString opt1 = (originalRange) ? "o" : "";
1350  TH1D * h1W = (onX) ? h2dW->ProjectionX("h1temp-W",firstbin,lastbin,opt1) : h2dW->ProjectionY("h1temp-W",firstbin,lastbin,opt1);
1351  TH1D * h1N = (onX) ? h2dN->ProjectionX("h1temp-N",firstbin,lastbin,opt1) : h2dN->ProjectionY("h1temp-N",firstbin,lastbin,opt1);
1352  h1W->SetDirectory(0); h1N->SetDirectory(0);
1353 
1354 
1355  // fill the bin content
1356  R__ASSERT( h1W->fN == p1->fN );
1357  R__ASSERT( h1N->fN == p1->fN );
1358  R__ASSERT( h1W->GetSumw2()->fN != 0); // h1W should always be a weighted histogram since h2dW is
1359  for (int i = 0; i < p1->fN ; ++i) {
1360  p1->fArray[i] = h1W->GetBinContent(i); // array of profile is sum of all values
1361  p1->GetSumw2()->fArray[i] = h1W->GetSumw2()->fArray[i]; // array of content square of profile is weight square of the W projected histogram
1362  p1->SetBinEntries(i, h1N->GetBinContent(i) );
1363  if (fBinSumw2.fN) p1->GetBinSumw2()->fArray[i] = h1N->GetSumw2()->fArray[i]; // sum of weight squares are stored to compute errors in h1N histogram
1364  }
1365  // delete the created histograms
1366  delete h2dW;
1367  delete h2dN;
1368  delete h1W;
1369  delete h1N;
1370 
1371  // Also we need to set the entries since they have not been correctly calculated during the projection
1372  // we can only set them to the effective entries
1373  p1->SetEntries( p1->GetEffectiveEntries() );
1374 
1375  return p1;
1376 }
1377 
1378 
1379 ////////////////////////////////////////////////////////////////////////////////
1380 /// Replace current statistics with the values in array stats
1381 
1383 {
1384  fTsumw = stats[0];
1385  fTsumw2 = stats[1];
1386  fTsumwx = stats[2];
1387  fTsumwx2 = stats[3];
1388  fTsumwy = stats[4];
1389  fTsumwy2 = stats[5];
1390  fTsumwxy = stats[6];
1391  fTsumwz = stats[7];
1392  fTsumwz2 = stats[8];
1393 }
1394 
1395 ////////////////////////////////////////////////////////////////////////////////
1396 /// -*Reset contents of a Profile2D histogram
1397 ///*-* =======================================
1398 
1400 {
1401  TH2D::Reset(option);
1402  fBinEntries.Reset();
1403  fBinSumw2.Reset();
1404  TString opt = option;
1405  opt.ToUpper();
1406  if (opt.Contains("ICE") && !opt.Contains("S")) return;
1407  fTsumwz = fTsumwz2 = 0;
1408 }
1409 
1410 
1411 ////////////////////////////////////////////////////////////////////////////////
1412 /// Profile histogram is resized along axis such that x is in the axis range.
1413 /// The new axis limits are recomputed by doubling iteratively
1414 /// the current axis range until the specified value x is within the limits.
1415 /// The algorithm makes a copy of the histogram, then loops on all bins
1416 /// of the old histogram to fill the extended histogram.
1417 /// Takes into account errors (Sumw2) if any.
1418 /// The axis must be extendable before invoking this function.
1419 /// Ex: h->GetXaxis()->SetCanExtend(kTRUE)
1420 
1422 {
1423  TProfile2D* hold = TProfileHelper::ExtendAxis(this, x, axis);
1424  if ( hold ) {
1425  fTsumwz = hold->fTsumwz;
1426  fTsumwz2 = hold->fTsumwz2;
1427  delete hold;
1428  }
1429 }
1430 
1431 ////////////////////////////////////////////////////////////////////////////////
1432 /// Rebin this histogram grouping nxgroup/nygroup bins along the xaxis/yaxis together.
1433 /// if newname is not blank a new profile hnew is created.
1434 /// else the current histogram is modified (default)
1435 /// The parameter nxgroup/nygroup indicate how many bins along the xaxis/yaxis of this
1436 /// have to be merged into one bin of hnew
1437 /// If the original profile has errors stored (via Sumw2), the resulting
1438 /// profile has new errors correctly calculated.
1439 ///
1440 /// examples: if hpxpy is an existing TProfile2D profile with 40 x 40 bins
1441 /// hpxpy->Rebin2D(); // merges two bins along the xaxis and yaxis in one
1442 /// // Carefull: previous contents of hpxpy are lost
1443 /// hpxpy->Rebin2D(3,5); // merges 3 bins along the xaxis and 5 bins along the yaxis in one
1444 /// // Carefull: previous contents of hpxpy are lost
1445 /// hpxpy->RebinX(5); //merges five bins along the xaxis in one in hpxpy
1446 /// TProfile2D *hnew = hpxpy->RebinY(5,"hnew"); // creates a new profile hnew
1447 /// // merging 5 bins of hpxpy along the yaxis in one bin
1448 ///
1449 /// NOTE : If nxgroup/nygroup is not an exact divider of the number of bins,
1450 /// along the xaxis/yaxis the top limit(s) of the rebinned profile
1451 /// is changed to the upper edge of the xbin=newxbins*nxgroup resp.
1452 /// ybin=newybins*nygroup and the remaining bins are added to
1453 /// the overflow bin.
1454 /// Statistics will be recomputed from the new bin contents.
1455 
1456 TProfile2D * TProfile2D::Rebin2D(Int_t nxgroup ,Int_t nygroup,const char * newname ) {
1457  //something to do?
1458  if((nxgroup != 1) || (nygroup != 1)){
1459  Int_t nxbins = fXaxis.GetNbins();
1460  Int_t nybins = fYaxis.GetNbins();
1465  if ((nxgroup <= 0) || (nxgroup > nxbins)) {
1466  Error("Rebin", "Illegal value of nxgroup=%d",nxgroup);
1467  return 0;
1468  }
1469  if ((nygroup <= 0) || (nygroup > nybins)) {
1470  Error("Rebin", "Illegal value of nygroup=%d",nygroup);
1471  return 0;
1472  }
1473 
1474  Int_t newxbins = nxbins/nxgroup;
1475  Int_t newybins = nybins/nygroup;
1476 
1477  //warning if bins are added to the overflow bin
1478  if(newxbins*nxgroup != nxbins) {
1479  Warning("Rebin", "nxgroup=%d should be an exact divider of nxbins=%d",nxgroup,nxbins);
1480  }
1481  if(newybins*nygroup != nybins) {
1482  Warning("Rebin", "nygroup=%d should be an exact divider of nybins=%d",nygroup,nybins);
1483  }
1484 
1485  //save old bin contents in new arrays
1486  Double_t *oldBins = new Double_t[(nxbins+2)*(nybins+2)];
1487  Double_t *oldCount = new Double_t[(nxbins+2)*(nybins+2)];
1488  Double_t *oldErrors = new Double_t[(nxbins+2)*(nybins+2)];
1489  Double_t *oldBinw2 = (fBinSumw2.fN ? new Double_t[(nxbins+2)*(nybins+2)] : 0 );
1490  Double_t *cu1 = GetW();
1491  Double_t *er1 = GetW2();
1492  Double_t *en1 = GetB();
1493  Double_t *ew1 = GetB2();
1494  for(Int_t ibin=0; ibin < (nxbins+2)*(nybins+2); ibin++){
1495  oldBins[ibin] = cu1[ibin];
1496  oldCount[ibin] = en1[ibin];
1497  oldErrors[ibin] = er1[ibin];
1498  if (ew1 && fBinSumw2.fN) oldBinw2[ibin] = ew1[ibin];
1499  }
1500 
1501  // create a clone of the old profile if newname is specified
1502  TProfile2D *hnew = this;
1503  if(newname && strlen(newname) > 0) {
1504  hnew = (TProfile2D*)Clone(newname);
1505  }
1506 
1507  // in case of nxgroup/nygroup not an exact divider of nxbins/nybins,
1508  // top limit is changed (see NOTE in method comment)
1509  if(newxbins*nxgroup != nxbins) {
1510  xmax = fXaxis.GetBinUpEdge(newxbins*nxgroup);
1511  hnew->fTsumw = 0; //stats must be reset because top bins will be moved to overflow bin
1512  }
1513  if(newybins*nygroup != nybins) {
1514  ymax = fYaxis.GetBinUpEdge(newybins*nygroup);
1515  hnew->fTsumw = 0; //stats must be reset because top bins will be moved to overflow bin
1516  }
1517 
1518  //rebin the axis
1519  if((fXaxis.GetXbins()->GetSize() > 0) || (fYaxis.GetXbins()->GetSize() > 0)){
1520  Double_t* xbins = new Double_t[newxbins+1];
1521  Double_t* ybins = new Double_t[newybins+1];
1522  for(Int_t i=0; i < newxbins+1; i++)
1523  xbins[i] = fXaxis.GetBinLowEdge(1+i*nxgroup);
1524  for(Int_t j=0; j < newybins+1; j++)
1525  ybins[j] = fYaxis.GetBinLowEdge(1+j*nygroup);
1526  hnew->SetBins(newxbins,xbins,newybins,ybins);
1527  delete [] xbins;
1528  delete [] ybins;
1529  }
1530  //fixed bin size
1531  else{
1532  hnew->SetBins(newxbins,xmin,xmax,newybins,ymin,ymax);
1533  }
1534 
1535  //merge bins
1536  Double_t *cu2 = hnew->GetW();
1537  Double_t *er2 = hnew->GetW2();
1538  Double_t *en2 = hnew->GetB();
1539  Double_t *ew2 = hnew->GetB2();
1540  Double_t binContent, binCount, binError, binSumw2;
1541  //connection between x and y bin number and linear global bin number:
1542  //global bin = xbin + (nxbins+2) * ybin
1543  Int_t oldxbin = 1;
1544  Int_t oldybin = 1;
1545  //global bin number
1546  Int_t bin;
1547  for(Int_t xbin = 1; xbin <= newxbins; xbin++){
1548  oldybin = 1;
1549  for(Int_t ybin = 1; ybin <= newybins; ybin++){
1550  binContent = 0;
1551  binCount = 0;
1552  binError = 0;
1553  binSumw2 = 0;
1554  for(Int_t i=0; i < nxgroup; i++){
1555  if(oldxbin + i > nxbins) break;
1556  for(Int_t j=0; j < nygroup; j++){
1557  if(oldybin + j > nybins) break;
1558  bin = oldxbin + i + (nxbins+2)*(oldybin+j);
1559  binContent += oldBins[bin];
1560  binCount += oldCount[bin];
1561  binError += oldErrors[bin];
1562  if(fBinSumw2.fN) binSumw2 += oldBinw2[bin];
1563  }
1564  }
1565  bin = xbin + (newxbins + 2)*ybin;
1566  cu2[bin] = binContent;
1567  er2[bin] = binError;
1568  en2[bin] = binCount;
1569  if(fBinSumw2.fN) ew2[bin] = binSumw2;
1570  oldybin += nygroup;
1571  }
1572  oldxbin += nxgroup;
1573  }
1574 
1575  //copy the underflow bin in x and y (0,0)
1576  cu2[0] = oldBins[0];
1577  er2[0] = oldErrors[0];
1578  en2[0] = oldCount[0];
1579  if(fBinSumw2.fN) ew2[0] = oldBinw2[0];
1580  //calculate overflow bin in x and y (newxbins+1,newybins+1)
1581  //therefore the oldxbin and oldybin from above are needed!
1582  binContent = 0;
1583  binCount = 0;
1584  binError = 0;
1585  binSumw2 = 0;
1586  for(Int_t i=oldxbin; i <= nxbins+1; i++){
1587  for(Int_t j=oldybin; j <= nybins+1; j++){
1588  //global bin number
1589  bin = i + (nxbins+2)*j;
1590  binContent += oldBins[bin];
1591  binCount += oldCount[bin];
1592  binError += oldErrors[bin];
1593  if(fBinSumw2.fN) binSumw2 += oldBinw2[bin];
1594  }
1595  }
1596  bin = (newxbins+2)*(newybins+2)-1;
1597  cu2[bin] = binContent;
1598  er2[bin] = binError;
1599  en2[bin] = binCount;
1600  if(fBinSumw2.fN) ew2[bin] = binSumw2;
1601  //calculate overflow bin in x and underflow bin in y (newxbins+1,0)
1602  binContent = 0;
1603  binCount = 0;
1604  binError = 0;
1605  binSumw2 = 0;
1606  for(Int_t i=oldxbin; i <= nxbins+1; i++){
1607  bin = i;
1608  binContent += oldBins[bin];
1609  binCount += oldCount[bin];
1610  binError += oldErrors[bin];
1611  if(fBinSumw2.fN) binSumw2 += oldBinw2[bin];
1612  }
1613  bin = newxbins + 1;
1614  cu2[bin] = binContent;
1615  er2[bin] = binError;
1616  en2[bin] = binCount;
1617  if(fBinSumw2.fN) ew2[bin] = binSumw2;
1618  //calculate underflow bin in x and overflow bin in y (0,newybins+1)
1619  binContent = 0;
1620  binCount = 0;
1621  binError = 0;
1622  binSumw2 = 0;
1623  for(Int_t i=oldybin; i <= nybins+1; i++){
1624  bin = i*(nxbins + 2);
1625  binContent += oldBins[bin];
1626  binCount += oldCount[bin];
1627  binError += oldErrors[bin];
1628  if(fBinSumw2.fN) binSumw2 += oldBinw2[bin];
1629  }
1630  bin = (newxbins + 2)*(newybins + 1);
1631  cu2[bin] = binContent;
1632  er2[bin] = binError;
1633  en2[bin] = binCount;
1634  if(fBinSumw2.fN) ew2[bin] = binSumw2;
1635  //calculate under/overflow contents in y for the new x bins
1636  Double_t binContentuf, binCountuf, binErroruf, binSumw2uf;
1637  Double_t binContentof, binCountof, binErrorof, binSumw2of;
1638  Int_t ufbin, ofbin;
1639  Int_t oldxbin2 = 1;
1640  for(Int_t xbin = 1; xbin <= newxbins; xbin++){
1641  binContentuf = 0;
1642  binCountuf = 0;
1643  binErroruf = 0;
1644  binSumw2uf = 0;
1645  binContentof = 0;
1646  binCountof = 0;
1647  binErrorof = 0;
1648  binSumw2of = 0;
1649  for(Int_t i = 0; i < nxgroup; i++){
1650  //index of under/overflow bin for y in old binning
1651  ufbin = (oldxbin2 + i);
1652  binContentuf += oldBins[ufbin];
1653  binCountuf += oldCount[ufbin];
1654  binErroruf += oldErrors[ufbin];
1655  if(fBinSumw2.fN) binSumw2uf += oldBinw2[ufbin];
1656  for(Int_t j = oldybin; j <= nybins+1; j++)
1657  {
1658  ofbin = ufbin + j*(nxbins + 2);
1659  binContentof += oldBins[ofbin];
1660  binCountof += oldCount[ofbin];
1661  binErrorof += oldErrors[ofbin];
1662  if(fBinSumw2.fN) binSumw2of += oldBinw2[ofbin];
1663  }
1664  }
1665  //index of under/overflow bin for y in new binning
1666  ufbin = xbin;
1667  ofbin = ufbin + (newybins + 1)*(newxbins + 2);
1668  cu2[ufbin] = binContentuf;
1669  er2[ufbin] = binErroruf;
1670  en2[ufbin] = binCountuf;
1671  if(fBinSumw2.fN) ew2[ufbin] = binSumw2uf;
1672  cu2[ofbin] = binContentof;
1673  er2[ofbin] = binErrorof;
1674  en2[ofbin] = binCountof;
1675  if(fBinSumw2.fN) ew2[ofbin] = binSumw2of;
1676 
1677  oldxbin2 += nxgroup;
1678  }
1679  //calculate under/overflow contents in x for the new y bins
1680  Int_t oldybin2 = 1;
1681  for(Int_t ybin = 1; ybin <= newybins; ybin++){
1682  binContentuf = 0;
1683  binCountuf = 0;
1684  binErroruf = 0;
1685  binSumw2uf = 0;
1686  binContentof = 0;
1687  binCountof = 0;
1688  binErrorof = 0;
1689  binSumw2of = 0;
1690  for(Int_t i = 0; i < nygroup; i++){
1691  //index of under/overflow bin for x in old binning
1692  ufbin = (oldybin2 + i)*(nxbins+2);
1693  binContentuf += oldBins[ufbin];
1694  binCountuf += oldCount[ufbin];
1695  binErroruf += oldErrors[ufbin];
1696  if(fBinSumw2.fN) binSumw2uf += oldBinw2[ufbin];
1697  for(Int_t j = oldxbin; j <= nxbins+1; j++)
1698  {
1699  ofbin = j + ufbin;
1700  binContentof += oldBins[ofbin];
1701  binCountof += oldCount[ofbin];
1702  binErrorof += oldErrors[ofbin];
1703  if(fBinSumw2.fN) binSumw2of += oldBinw2[ofbin];
1704  }
1705  }
1706  //index of under/overflow bin for x in new binning
1707  ufbin = ybin * (newxbins + 2);
1708  ofbin = newxbins + 1 + ufbin;
1709  cu2[ufbin] = binContentuf;
1710  er2[ufbin] = binErroruf;
1711  en2[ufbin] = binCountuf;
1712  if(fBinSumw2.fN) ew2[ufbin] = binSumw2uf;
1713  cu2[ofbin] = binContentof;
1714  er2[ofbin] = binErrorof;
1715  en2[ofbin] = binCountof;
1716  if(fBinSumw2.fN) ew2[ofbin] = binSumw2of;
1717 
1718  oldybin2 += nygroup;
1719  }
1720 
1721  delete [] oldBins;
1722  delete [] oldCount;
1723  delete [] oldErrors;
1724  if (oldBinw2) delete [] oldBinw2;
1725 
1726  return hnew;
1727  }
1728  //nxgroup == nygroup == 1
1729  else{
1730  if((newname) && (strlen(newname) > 0))
1731  return (TProfile2D*)Clone(newname);
1732  else
1733  return this;
1734  }
1735 }
1736 
1737 ////////////////////////////////////////////////////////////////////////////////
1738 /// Rebin only the X axis
1739 /// see Rebin2D
1740 
1741 TProfile2D * TProfile2D::RebinX(Int_t ngroup,const char * newname ) {
1742  return Rebin2D(ngroup,1,newname);
1743 }
1744 
1745 ////////////////////////////////////////////////////////////////////////////////
1746 /// Rebin only the Y axis
1747 /// see Rebin2D
1748 
1749 TProfile2D * TProfile2D::RebinY(Int_t ngroup,const char * newname ) {
1750  return Rebin2D(1,ngroup,newname);
1751 }
1752 
1753 ////////////////////////////////////////////////////////////////////////////////
1754 /// Save primitive as a C++ statement(s) on output stream out
1755 
1756 void TProfile2D::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
1757 {
1758  //Note the following restrictions in the code generated:
1759  // - variable bin size not implemented
1760  // - SetErrorOption not implemented
1761 
1762 
1763  char quote = '"';
1764  out <<" "<<std::endl;
1765  out <<" "<<ClassName()<<" *";
1766 
1767  out << GetName() << " = new " << ClassName() << "(" << quote
1768  << GetName() << quote << "," << quote<< GetTitle() << quote
1769  << "," << GetXaxis()->GetNbins();
1770  out << "," << GetXaxis()->GetXmin()
1771  << "," << GetXaxis()->GetXmax();
1772  out << "," << GetYaxis()->GetNbins();
1773  out << "," << GetYaxis()->GetXmin()
1774  << "," << GetYaxis()->GetXmax();
1775  out << "," << fZmin
1776  << "," << fZmax;
1777  out << ");" << std::endl;
1778 
1779 
1780  // save bin entries
1781  Int_t bin;
1782  for (bin=0;bin<fNcells;bin++) {
1783  Double_t bi = GetBinEntries(bin);
1784  if (bi) {
1785  out<<" "<<GetName()<<"->SetBinEntries("<<bin<<","<<bi<<");"<<std::endl;
1786  }
1787  }
1788  //save bin contents
1789  for (bin=0;bin<fNcells;bin++) {
1790  Double_t bc = fArray[bin];
1791  if (bc) {
1792  out<<" "<<GetName()<<"->SetBinContent("<<bin<<","<<bc<<");"<<std::endl;
1793  }
1794  }
1795  // save bin errors
1796  if (fSumw2.fN) {
1797  for (bin=0;bin<fNcells;bin++) {
1798  Double_t be = TMath::Sqrt(fSumw2.fArray[bin]);
1799  if (be) {
1800  out<<" "<<GetName()<<"->SetBinError("<<bin<<","<<be<<");"<<std::endl;
1801  }
1802  }
1803  }
1804 
1805  TH1::SavePrimitiveHelp(out, GetName(), option);
1806 }
1807 
1808 ////////////////////////////////////////////////////////////////////////////////
1809 /// Multiply this profile2D by a constant c1
1810 ///
1811 /// this = c1*this
1812 ///
1813 /// This function uses the services of TProfile2D::Add
1814 
1816 {
1817  TProfileHelper::Scale(this, c1, option);
1818 }
1819 
1820 ////////////////////////////////////////////////////////////////////////////////
1821 /// Set the number of entries in bin -
1822 
1824 {
1825  TProfileHelper::SetBinEntries(this, bin, w);
1826 }
1827 
1828 ////////////////////////////////////////////////////////////////////////////////
1829 /// Redefine x and y axis parameters -
1830 
1832 {
1833  TH1::SetBins(nx,xmin, xmax,ny, ymin,ymax);
1836 }
1837 
1838 ////////////////////////////////////////////////////////////////////////////////
1839 /// Redefine x and y axis parameters for variable bin sizes - -*
1840 
1841 void TProfile2D::SetBins(Int_t nx, const Double_t *xbins, Int_t ny, const Double_t *ybins)
1842 {
1843  TH1::SetBins(nx,xbins,ny,ybins);
1846 }
1847 
1848 ////////////////////////////////////////////////////////////////////////////////
1849 /// Set total number of bins including under/overflow
1850 /// Reallocate bin contents array
1851 
1853 {
1856 }
1857 
1858 ////////////////////////////////////////////////////////////////////////////////
1859 /// set the buffer size in units of 8 bytes (double)
1860 
1862 {
1863  if (fBuffer) {
1864  BufferEmpty();
1865  delete [] fBuffer;
1866  fBuffer = 0;
1867  }
1868  if (buffersize <= 0) {
1869  fBufferSize = 0;
1870  return;
1871  }
1872  if (buffersize < 100) buffersize = 100;
1873  fBufferSize = 1 + 4*buffersize;
1874  fBuffer = new Double_t[fBufferSize];
1875  memset(fBuffer,0,sizeof(Double_t)*fBufferSize);
1876 }
1877 
1878 ////////////////////////////////////////////////////////////////////////////////
1879 /// Set option to compute profile2D errors
1880 /// The computation of the bin errors is based on the parameter option:
1881 /// - ' ' (Default) The bin errors are the standard error on the mean of the bin profiled values (Z),
1882 /// i.e. the standard error of the bin contents.
1883 /// Note that if TProfile::Approximate() is called, an approximation is used when
1884 /// the spread in Z is 0 and the number of bin entries is > 0
1885 /// - 's' The bin errors are the standard deviations of the Z bin values
1886 /// Note that if TProfile::Approximate() is called, an approximation is used when
1887 /// the spread in Z is 0 and the number of bin entries is > 0
1888 /// - 'i' Errors are as in default case (standard errors of the bin contents)
1889 /// The only difference is for the case when the spread in Z is zero.
1890 /// In this case for N > 0 the error is 1./SQRT(12.*N)
1891 /// - 'g' Errors are 1./SQRT(W) for W not equal to 0 and 0 for W = 0.
1892 /// W is the sum in the bin of the weights of the profile.
1893 /// This option is for combining measurements z +/- dz,
1894 /// and the profile is filled with values y and weights z = 1/dz**2
1895 ///
1896 /// See TProfile::BuildOptions for a detailed explanation of all options
1897 
1899 {
1900  TProfileHelper::SetErrorOption(this, option);
1901 }
1902 
1903 ////////////////////////////////////////////////////////////////////////////////
1904 /// Stream an object of class TProfile2D.
1905 
1906 void TProfile2D::Streamer(TBuffer &R__b)
1907 {
1908  if (R__b.IsReading()) {
1909  UInt_t R__s, R__c;
1910  Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
1911  if (R__v > 2) {
1912  R__b.ReadClassBuffer(TProfile2D::Class(), this, R__v, R__s, R__c);
1913  return;
1914  }
1915  //====process old versions before automatic schema evolution
1916  TH2D::Streamer(R__b);
1917  fBinEntries.Streamer(R__b);
1918  Int_t errorMode;
1919  R__b >> errorMode;
1920  fErrorMode = (EErrorType)errorMode;
1921  if (R__v < 2) {
1922  Float_t zmin,zmax;
1923  R__b >> zmin; fZmin = zmin;
1924  R__b >> zmax; fZmax = zmax;
1925  } else {
1926  R__b >> fZmin;
1927  R__b >> fZmax;
1928  }
1929  R__b.CheckByteCount(R__s, R__c, TProfile2D::IsA());
1930  //====end of old versions
1931 
1932  } else {
1933  R__b.WriteClassBuffer(TProfile2D::Class(),this);
1934  }
1935 }
1936 
1937 ////////////////////////////////////////////////////////////////////////////////
1938 /// Create/Delete structure to store sum of squares of weights per bin
1939 /// This is needed to compute the correct statistical quantities
1940 /// of a profile filled with weights
1941 ///
1942 /// This function is automatically called when the histogram is created
1943 /// if the static function TH1::SetDefaultSumw2 has been called before.
1944 /// If flag is false the structure is deleted
1945 
1947 {
1948  TProfileHelper::Sumw2(this, flag);
1949 }
const int nx
Definition: kalman.C:16
virtual void SetBinsLength(Int_t n=-1)
Set total number of bins including under/overflow Reallocate bin contents array.
Definition: TH2.cxx:4036
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:429
virtual Long64_t Merge(TCollection *list)
Merge all histograms in the collection in this histogram.
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:52
static void SetBinEntries(T *p, Int_t bin, Double_t w)
static void Scale(T *p, Double_t c1, Option_t *option)
float xmin
Definition: THbookFile.cxx:93
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4629
long long Long64_t
Definition: RtypesCore.h:69
virtual Double_t GetBinEntries(Int_t bin) const
Return bin entries of a Profile2D histogram.
Definition: TProfile2D.cxx:797
Bool_t IsReading() const
Definition: TBuffer.h:81
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:487
virtual void Sumw2(Bool_t flag=kTRUE)
Create/Delete structure to store sum of squares of weights per bin This is needed to compute the corr...
short Version_t
Definition: RtypesCore.h:61
static Double_t GetBinError(T *p, Int_t bin)
float Float_t
Definition: RtypesCore.h:53
virtual void SetDirectory(TDirectory *dir)
By default when an histogram is created, it is added to the list of histogram objects in the current ...
Definition: TH1.cxx:8266
const char Option_t
Definition: RtypesCore.h:62
TCanvas * c1
Definition: legend1.C:2
void Reset()
Definition: TArrayD.h:49
virtual void ExtendAxis(Double_t x, TAxis *axis)
Profile histogram is resized along axis such that x is in the axis range.
float ymin
Definition: THbookFile.cxx:93
TAxis fYaxis
Definition: TH1.h:103
static void SetErrorOption(T *p, Option_t *opt)
virtual void SetBins(Int_t nx, Double_t xmin, Double_t xmax)
Redefine x axis parameters.
Definition: TH1.cxx:8090
static void LabelsInflate(T *p, Option_t *)
static Bool_t fgStatOverflows
flag to add histograms to the directory
Definition: TH1.h:128
virtual void Copy(TObject &hnew) const
Copy a Profile2D histogram to a new profile2D histogram*.
Definition: TProfile2D.cxx:364
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1101
Buffer base class used for serializing objects.
Definition: TBuffer.h:40
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition: TAxis.cxx:489
#define R__ASSERT(e)
Definition: TError.h:98
static THLimitsFinder * GetLimitsFinder()
Return pointer to the current finder.
virtual Int_t CheckByteCount(UInt_t startpos, UInt_t bcnt, const TClass *clss)=0
virtual TProfile2D * RebinY(Int_t ngroup=2, const char *newname="")
Rebin only the Y axis see Rebin2D.
Basic string class.
Definition: TString.h:137
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:170
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1088
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
const Bool_t kFALSE
Definition: Rtypes.h:92
virtual TProfile2D * Rebin2D(Int_t nxgroup=2, Int_t nygroup=2, const char *newname="")
Rebin this histogram grouping nxgroup/nygroup bins along the xaxis/yaxis together.
TArrayD fSumw2
Definition: TH1.h:116
virtual ~TProfile2D()
Default destructor for Profile2D histograms.
Definition: TProfile2D.cxx:82
virtual Int_t GetNbinsX() const
Definition: TH1.h:296
const TKDTreeBinning * bins
Profile Historam.
Definition: TProfile.h:34
void Copy(TArrayD &array) const
Definition: TArrayD.h:44
virtual Int_t FindGoodLimits(TH1 *h, Double_t xmin, Double_t xmax)
compute the best axis limits for the X axis.
Double_t * GetB()
Definition: TProfile2D.h:66
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save primitive as a C++ statement(s) on output stream out.
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:310
virtual Bool_t Multiply(TF1 *h1, Double_t c1=1)
Performs the operation: this = this*c1*f1.
EErrorType
Definition: TProfile.h:30
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:732
Double_t fTsumwx2
Definition: TH1.h:111
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2565
static void BuildArray(T *p)
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:946
static Bool_t fgApproximate
Definition: TProfile2D.h:45
virtual TProfile * 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 avera...
Int_t Fill(const Double_t *v)
Definition: TProfile2D.h:54
Double_t x[n]
Definition: legend1.C:17
virtual Bool_t Divide(TF1 *h1, Double_t c1=1)
Performs the operation: this = this/(c1*f1) This function is not implemented.
Definition: TProfile2D.cxx:394
void Class()
Definition: Class.C:29
THashList implements a hybrid collection class consisting of a hash table and a list to store TObject...
Definition: THashList.h:36
void BuildOptions(Double_t zmin, Double_t zmax, Option_t *option)
Set Profile2D histogram structure and options.
Definition: TProfile2D.cxx:172
const int ny
Definition: kalman.C:17
Double_t * GetB2()
Definition: TProfile2D.h:67
TProfile * ProfileX(const char *name="_pfx", Int_t firstybin=0, Int_t lastybin=-1, Option_t *option="") const
*-*Project a 2-D histogram into a profile histogram along X *-* =====================================...
TArrayD fBinEntries
Definition: TProfile2D.h:37
virtual void SetErrorOption(Option_t *option="")
Set option to compute profile2D errors The computation of the bin errors is based on the parameter op...
virtual TArrayD * GetBinSumw2()
Definition: TProfile.h:111
static double p2(double t, double a, double b, double c)
static void LabelsDeflate(T *p, Option_t *)
static void Approximate(Bool_t approx=kTRUE)
static function set the fgApproximate flag.
Definition: TProfile2D.cxx:258
TH1D * ProjectionX(const char *name="_px", Int_t firstybin=0, Int_t lastybin=-1, Option_t *option="") const
Project a 2-D histogram into a 1-D histogram along X.
Definition: TH2.cxx:2513
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMath.h:1002
Double_t * fArray
Definition: TArrayD.h:32
static T * ExtendAxis(T *p, Double_t x, TAxis *axis)
Double_t fTsumwy2
Definition: TH2.h:41
virtual void LabelsDeflate(Option_t *axis="X")
Reduce the number of bins for this axis to the number of bins having a label.
Definition: TProfile2D.cxx:933
TH1D * ProjectionY(const char *name="_py", Int_t firstxbin=0, Int_t lastxbin=-1, Option_t *option="") const
Project a 2-D histogram into a 1-D histogram along Y.
Definition: TH2.cxx:2553
virtual TArrayD * GetSumw2()
Definition: TH1.h:312
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
TH1F * h1
Definition: legend1.C:5
virtual void SetBinError(Int_t bin, Double_t error)
see convention for numbering bins in TH1::GetBin
Definition: TH1.cxx:8528
Double_t GetXmin() const
Definition: TAxis.h:137
virtual void SetBinsLength(Int_t n=-1)
Set total number of bins including under/overflow Reallocate bin contents array.
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this profile2D by a constant c1.
char * out
Definition: TBase64.cxx:29
Double_t fTsumwx
Definition: TH1.h:110
Double_t fZmin
Definition: TProfile2D.h:39
virtual Double_t GetEffectiveEntries() const
number of effective entries of the histogram, neff = (Sum of weights )^2 / (Sum of weight^2 ) In case...
Definition: TH1.cxx:4073
virtual void SetUniqueID(UInt_t uid)
Set the unique object id.
Definition: TObject.cxx:743
virtual void AddBinContent(Int_t bin)
Increment bin content by 1.
Definition: TH2.h:314
Int_t fN
Definition: TArray.h:40
void Clear(Option_t *option="")
Remove all objects from the list.
Definition: THashList.cxx:167
virtual Bool_t Add(TF1 *h1, Double_t c1=1, Option_t *option="")
Performs the operation: this = this + c1*f1.
Definition: TProfile2D.cxx:201
float ymax
Definition: THbookFile.cxx:93
virtual void Copy(TObject &hnew) const
Copy.
Definition: TH2.cxx:4016
virtual Double_t GetBinContent(Int_t bin) const
Return bin content of a Profile2D histogram.
Definition: TProfile2D.cxx:784
Class to manage histogram axis.
Definition: TAxis.h:36
Int_t GetNbins() const
Definition: TAxis.h:125
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition: TAxis.cxx:499
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:187
Double_t fTsumwxy
Definition: TH2.h:42
Double_t fTsumwy
Definition: TH2.h:40
Collection abstract base class.
Definition: TCollection.h:48
TClass * IsA() const
virtual Int_t BufferEmpty(Int_t action=0)
Fill histogram with all entries in the buffer.
Definition: TProfile2D.cxx:272
static Int_t fgBufferSize
Definition: TH1.h:126
Double_t fZmax
Definition: TProfile2D.h:40
virtual void PutStats(Double_t *stats)
Replace current statistics with the values in array stats.
Double_t fEntries
Definition: TH1.h:107
unsigned int UInt_t
Definition: RtypesCore.h:42
Bool_t TestBit(UInt_t f) const
Definition: TObject.h:173
Double_t fTsumwz2
Definition: TProfile2D.h:43
virtual TProfile2D * RebinX(Int_t ngroup=2, const char *newname="")
Rebin only the X axis see Rebin2D.
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
virtual void LabelsOption(Option_t *option="h", Option_t *axis="X")
Set option(s) to draw axis with labels option = "a" sort by alphabetic order = ">" sort by decreasing...
Definition: TProfile2D.cxx:958
TAxis * GetYaxis()
Definition: TH1.h:320
virtual Double_t GetBinError(Int_t bin) const
-*Return bin error of a Profile2D histogram
Definition: TProfile2D.cxx:838
static double p1(double t, double a, double b)
float xmax
Definition: THbookFile.cxx:93
virtual void SetBuffer(Int_t buffersize, Option_t *option="")
set the buffer size in units of 8 bytes (double)
virtual void LabelsInflate(Option_t *axis="X")
Double the number of bins for axis.
Definition: TProfile2D.cxx:943
Bool_t IsNull() const
Definition: TString.h:387
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:613
virtual void GetStats(Double_t *stats) const
fill the array stats from the contents of this profile The array stats must be correctly dimensionned...
Definition: TProfile2D.cxx:873
virtual Int_t GetBin(Int_t binx, Int_t biny, Int_t binz=0) const
Return Global bin number corresponding to binx,y,z.
Definition: TH2.cxx:956
const Double_t * GetArray() const
Definition: TArrayD.h:45
Int_t GetSize() const
Definition: TArray.h:49
Double_t * GetW()
Definition: TProfile2D.h:68
virtual Int_t FindBin(Double_t x)
Find bin number corresponding to abscissa x.
Definition: TAxis.cxx:264
static Bool_t Add(T *p, const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2=1)
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
Double_t fTsumw2
Definition: TH1.h:109
virtual Int_t BufferFill(Double_t, Double_t)
accumulate arguments in buffer.
Definition: TProfile2D.h:47
virtual Int_t GetSize() const
Definition: TCollection.h:95
return c2
Definition: legend2.C:14
#define ClassImp(name)
Definition: Rtypes.h:279
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:415
double Double_t
Definition: RtypesCore.h:55
virtual void SavePrimitiveHelp(std::ostream &out, const char *hname, Option_t *option="")
helper function for the SavePrimitive functions from TH1 or classes derived from TH1, eg TProfile, TProfile2D.
Definition: TH1.cxx:6864
virtual Double_t GetBinEffectiveEntries(Int_t bin)
Return bin effective entries for a weighted filled Profile histogram.
Definition: TProfile2D.cxx:815
Double_t fTsumw
Definition: TH1.h:108
static Double_t GetBinEffectiveEntries(T *p, Int_t bin)
Double_t GetXmax() const
Definition: TAxis.h:138
Double_t y[n]
Definition: legend1.C:17
TH2D * ProjectionXY(const char *name="_pxy", Option_t *option="e") const
*-*Project this profile2D into a 2-D histogram along X,Y -* *-* =====================================...
The TH1 histogram class.
Definition: TH1.h:80
virtual void SetBinEntries(Int_t bin, Double_t w)
Set the number of entries in bin -.
Profile2D histograms are used to display the mean value of Z and its RMS for each cell in X...
Definition: TProfile2D.h:31
Array of doubles (64 bits per element).
Definition: TArrayD.h:29
#define name(a, b)
Definition: linkTestLib0.cpp:5
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:440
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:449
Mother of all ROOT objects.
Definition: TObject.h:58
virtual void SetBinEntries(Int_t bin, Double_t w)
Set the number of entries in bin.
Definition: TProfile.cxx:1604
virtual Int_t GetNbinsY() const
Definition: TH1.h:297
THashList * GetLabels() const
Definition: TAxis.h:122
virtual void Add(TObject *obj)
Definition: TList.h:81
Int_t fBufferSize
Definition: TH1.h:119
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:567
1-Dim function class
Definition: TF1.h:149
Int_t IsNaN(Double_t x)
Definition: TMath.h:617
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8350
virtual Bool_t CanExtendAllAxes() const
returns true if all axes are extendable
Definition: TH1.cxx:6197
TArrayD fBinSumw2
Definition: TProfile2D.h:44
virtual void SetEntries(Double_t n)
Definition: TH1.h:382
TAxis fXaxis
Definition: TH1.h:102
static void Sumw2(T *p, Bool_t flag)
const TArrayD * GetXbins() const
Definition: TAxis.h:134
Option_t * GetErrorOption() const
-*Return option to compute profile2D errors *-* ========================================= ...
Definition: TProfile2D.cxx:847
void ResetBit(UInt_t f)
Definition: TObject.h:172
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2689
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
virtual void Sumw2(Bool_t flag=kTRUE)
Create/delete structure to store sum of squares of weights per bin — This is needed to compute the c...
Definition: TProfile.cxx:1729
const Bool_t kTRUE
Definition: Rtypes.h:91
TObject * obj
void SetBins(const Int_t *nbins, const Double_t *range)
Definition: TProfile2D.h:52
Double_t * fBuffer
Definition: TH1.h:120
EErrorType fErrorMode
Definition: TProfile2D.h:38
void Set(Int_t n)
Set size of this array to n doubles.
Definition: TArrayD.cxx:104
Double_t fTsumwz
True when TProfile2D::Scale is called.
Definition: TProfile2D.h:42
const Int_t n
Definition: legend1.C:16
static Long64_t Merge(T *p, TCollection *list)
Double_t * GetW2()
Definition: TProfile2D.h:69
TProfile * ProfileY(const char *name="_pfy", Int_t firstxbin=0, Int_t lastxbin=-1, Option_t *option="") const
*-*Project a 2-D histogram into a profile histogram along X *-* =====================================...
Bool_t fScaling
Definition: TProfile2D.h:41
TAxis * GetXaxis()
Definition: TH1.h:319
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
Int_t fNcells
Definition: TH1.h:101
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:297
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:904