Logo ROOT  
Reference Guide
TH3.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2// Author: Rene Brun 27/10/95
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 "TROOT.h"
13#include "TBuffer.h"
14#include "TClass.h"
15#include "THashList.h"
16#include "TH3.h"
17#include "TProfile2D.h"
18#include "TH2.h"
19#include "TF3.h"
20#include "TVirtualPad.h"
21#include "TVirtualHistPainter.h"
22#include "THLimitsFinder.h"
23#include "TRandom.h"
24#include "TError.h"
25#include "TMath.h"
26#include "TObjString.h"
27
29
30/** \addtogroup Histograms
31@{
32\class TH3C
33\brief 3-D histogram with a byte per channel (see TH1 documentation)
34\class TH3S
35\brief 3-D histogram with a short per channel (see TH1 documentation)
36\class TH3I
37\brief 3-D histogram with an int per channel (see TH1 documentation)}
38\class TH3F
39\brief 3-D histogram with a float per channel (see TH1 documentation)}
40\class TH3D
41\brief 3-D histogram with a double per channel (see TH1 documentation)}
42@}
43*/
44
45/** \class TH3
46 \ingroup Histograms
47The 3-D histogram classes derived from the 1-D histogram classes.
48All operations are supported (fill, fit).
49Drawing is currently restricted to one single option.
50A cloud of points is drawn. The number of points is proportional to
51cell content.
52
53- TH3C a 3-D histogram with one byte per cell (char)
54- TH3S a 3-D histogram with two bytes per cell (short integer)
55- TH3I a 3-D histogram with four bytes per cell (32 bits integer)
56- TH3F a 3-D histogram with four bytes per cell (float)
57- TH3D a 3-D histogram with eight bytes per cell (double)
58*/
59
60
61////////////////////////////////////////////////////////////////////////////////
62/// Default constructor.
63
65{
66 fDimension = 3;
69}
70
71
72////////////////////////////////////////////////////////////////////////////////
73/// Constructor for fix bin size 3-D histograms.
74/// Creates the main histogram structure.
75///
76/// \param[in] name name of histogram (avoid blanks)
77/// \param[in] title histogram title.
78/// If title is of the form `stringt;stringx;stringy;stringz`,
79/// the histogram title is set to `stringt`,
80/// the x axis title to `stringx`, the y axis title to `stringy`, etc.
81/// \param[in] nbinsx number of bins along the X axis
82/// \param[in] xlow low edge of the X axis first bin
83/// \param[in] xup upper edge of the X axis last bin (not included in last bin)
84/// \param[in] nbinsy number of bins along the Y axis
85/// \param[in] ylow low edge of the Y axis first bin
86/// \param[in] yup upper edge of the Y axis last bin (not included in last bin)
87/// \param[in] nbinsz number of bins along the Z axis
88/// \param[in] zlow low edge of the Z axis first bin
89/// \param[in] zup upper edge of the Z axis last bin (not included in last bin)
90
91TH3::TH3(const char *name,const char *title,Int_t nbinsx,Double_t xlow,Double_t xup
92 ,Int_t nbinsy,Double_t ylow,Double_t yup
93 ,Int_t nbinsz,Double_t zlow,Double_t zup)
94 :TH1(name,title,nbinsx,xlow,xup),
95 TAtt3D()
96{
97 fDimension = 3;
98 if (nbinsy <= 0) {
99 Warning("TH3","nbinsy is <=0 - set to nbinsy = 1");
100 nbinsy = 1;
101 }
102 if (nbinsz <= 0) {
103 Warning("TH3","nbinsz is <=0 - set to nbinsz = 1");
104 nbinsz = 1;
105 }
106 fYaxis.Set(nbinsy,ylow,yup);
107 fZaxis.Set(nbinsz,zlow,zup);
108 fNcells = (nbinsx+2)*(nbinsy+2)*(nbinsz+2);
109 fTsumwy = fTsumwy2 = fTsumwxy = 0;
111}
112
113
114////////////////////////////////////////////////////////////////////////////////
115/// Constructor for variable bin size (along X, Y and Z axis) 3-D histograms using input
116/// arrays of type float.
117///
118/// \param[in] name name of histogram (avoid blanks)
119/// \param[in] title histogram title.
120/// If title is of the form `stringt;stringx;stringy;stringz`
121/// the histogram title is set to `stringt`,
122/// the x axis title to `stringx`, the y axis title to `stringy`, etc.
123/// \param[in] nbinsx number of bins
124/// \param[in] xbins array of low-edges for each bin.
125/// This is an array of type float and size nbinsx+1
126/// \param[in] nbinsy number of bins
127/// \param[in] ybins array of low-edges for each bin.
128/// This is an array of type float and size nbinsy+1
129/// \param[in] nbinsz number of bins
130/// \param[in] zbins array of low-edges for each bin.
131/// This is an array of type float and size nbinsz+1
132
133TH3::TH3(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins
134 ,Int_t nbinsy,const Float_t *ybins
135 ,Int_t nbinsz,const Float_t *zbins)
136 :TH1(name,title,nbinsx,xbins),
137 TAtt3D()
138{
139 fDimension = 3;
140 if (nbinsy <= 0) {Warning("TH3","nbinsy is <=0 - set to nbinsy = 1"); nbinsy = 1; }
141 if (nbinsz <= 0) nbinsz = 1;
142 if (ybins) fYaxis.Set(nbinsy,ybins);
143 else fYaxis.Set(nbinsy,0,1);
144 if (zbins) fZaxis.Set(nbinsz,zbins);
145 else fZaxis.Set(nbinsz,0,1);
146 fNcells = (nbinsx+2)*(nbinsy+2)*(nbinsz+2);
147 fTsumwy = fTsumwy2 = fTsumwxy = 0;
149}
150
151
152////////////////////////////////////////////////////////////////////////////////
153/// Constructor for variable bin size (along X, Y and Z axis) 3-D histograms using input
154/// arrays of type double.
155///
156/// \param[in] name name of histogram (avoid blanks)
157/// \param[in] title histogram title.
158/// If title is of the form `stringt;stringx;stringy;stringz`
159/// the histogram title is set to `stringt`,
160/// the x axis title to `stringx`, the y axis title to `stringy`, etc.
161/// \param[in] nbinsx number of bins
162/// \param[in] xbins array of low-edges for each bin.
163/// This is an array of type double and size nbinsx+1
164/// \param[in] nbinsy number of bins
165/// \param[in] ybins array of low-edges for each bin.
166/// This is an array of type double and size nbinsy+1
167/// \param[in] nbinsz number of bins
168/// \param[in] zbins array of low-edges for each bin.
169/// This is an array of type double and size nbinsz+1
170
171TH3::TH3(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins
172 ,Int_t nbinsy,const Double_t *ybins
173 ,Int_t nbinsz,const Double_t *zbins)
174 :TH1(name,title,nbinsx,xbins),
175 TAtt3D()
176{
177 fDimension = 3;
178 if (nbinsy <= 0) {Warning("TH3","nbinsy is <=0 - set to nbinsy = 1"); nbinsy = 1; }
179 if (nbinsz <= 0) nbinsz = 1;
180 if (ybins) fYaxis.Set(nbinsy,ybins);
181 else fYaxis.Set(nbinsy,0,1);
182 if (zbins) fZaxis.Set(nbinsz,zbins);
183 else fZaxis.Set(nbinsz,0,1);
184 fNcells = (nbinsx+2)*(nbinsy+2)*(nbinsz+2);
185 fTsumwy = fTsumwy2 = fTsumwxy = 0;
187}
188
189
190////////////////////////////////////////////////////////////////////////////////
191/// Private copy constructor.
192/// One should use the copy constructor of the derived classes (e.g. TH3D, TH3F ...).
193/// The list of functions is not copied. (Use Clone() if needed)
194
195TH3::TH3(const TH3 &h) : TH1(), TAtt3D()
196{
197 ((TH3&)h).Copy(*this);
198}
199
200
201////////////////////////////////////////////////////////////////////////////////
202/// Destructor.
203
205{
206}
207
208
209////////////////////////////////////////////////////////////////////////////////
210/// Copy.
211
212void TH3::Copy(TObject &obj) const
213{
214 TH1::Copy(obj);
215 ((TH3&)obj).fTsumwy = fTsumwy;
216 ((TH3&)obj).fTsumwy2 = fTsumwy2;
217 ((TH3&)obj).fTsumwxy = fTsumwxy;
218 ((TH3&)obj).fTsumwz = fTsumwz;
219 ((TH3&)obj).fTsumwz2 = fTsumwz2;
220 ((TH3&)obj).fTsumwxz = fTsumwxz;
221 ((TH3&)obj).fTsumwyz = fTsumwyz;
222}
223
224
225////////////////////////////////////////////////////////////////////////////////
226/// Fill histogram with all entries in the buffer.
227/// action = -1 histogram is reset and refilled from the buffer (called by THistPainter::Paint)
228/// action = 0 histogram is filled from the buffer
229/// action = 1 histogram is filled and buffer is deleted
230/// The buffer is automatically deleted when the number of entries
231/// in the buffer is greater than the number of entries in the histogram
232
234{
235 // do we need to compute the bin size?
236 if (!fBuffer) return 0;
237 Int_t nbentries = (Int_t)fBuffer[0];
238 if (!nbentries) return 0;
239 Double_t *buffer = fBuffer;
240 if (nbentries < 0) {
241 if (action == 0) return 0;
242 nbentries = -nbentries;
243 fBuffer=0;
244 Reset("ICES");
245 fBuffer = buffer;
246 }
247 if (CanExtendAllAxes() || fXaxis.GetXmax() <= fXaxis.GetXmin() ||
248 fYaxis.GetXmax() <= fYaxis.GetXmin() ||
249 fZaxis.GetXmax() <= fZaxis.GetXmin()) {
250 //find min, max of entries in buffer
251 Double_t xmin = fBuffer[2];
253 Double_t ymin = fBuffer[3];
255 Double_t zmin = fBuffer[4];
256 Double_t zmax = zmin;
257 for (Int_t i=1;i<nbentries;i++) {
258 Double_t x = fBuffer[4*i+2];
259 if (x < xmin) xmin = x;
260 if (x > xmax) xmax = x;
261 Double_t y = fBuffer[4*i+3];
262 if (y < ymin) ymin = y;
263 if (y > ymax) ymax = y;
264 Double_t z = fBuffer[4*i+4];
265 if (z < zmin) zmin = z;
266 if (z > zmax) zmax = z;
267 }
270 } else {
271 fBuffer = 0;
272 Int_t keep = fBufferSize; fBufferSize = 0;
277 if (zmin < fZaxis.GetXmin()) ExtendAxis(zmin,&fZaxis);
278 if (zmax >= fZaxis.GetXmax()) ExtendAxis(zmax,&fZaxis);
279 fBuffer = buffer;
280 fBufferSize = keep;
281 }
282 }
283 fBuffer = 0;
284
285 for (Int_t i=0;i<nbentries;i++) {
286 Fill(buffer[4*i+2],buffer[4*i+3],buffer[4*i+4],buffer[4*i+1]);
287 }
288 fBuffer = buffer;
289
290 if (action > 0) { delete [] fBuffer; fBuffer = 0; fBufferSize = 0;}
291 else {
292 if (nbentries == (Int_t)fEntries) fBuffer[0] = -nbentries;
293 else fBuffer[0] = 0;
294 }
295 return nbentries;
296}
297
298
299////////////////////////////////////////////////////////////////////////////////
300/// Accumulate arguments in buffer. When buffer is full, empty the buffer
301///
302/// - `fBuffer[0]` = number of entries in buffer
303/// - `fBuffer[1]` = w of first entry
304/// - `fBuffer[2]` = x of first entry
305/// - `fBuffer[3]` = y of first entry
306/// - `fBuffer[4]` = z of first entry
307
309{
310 if (!fBuffer) return -3;
311 Int_t nbentries = (Int_t)fBuffer[0];
312 if (nbentries < 0) {
313 nbentries = -nbentries;
314 fBuffer[0] = nbentries;
315 if (fEntries > 0) {
316 Double_t *buffer = fBuffer; fBuffer=0;
317 Reset("ICES");
318 fBuffer = buffer;
319 }
320 }
321 if (4*nbentries+4 >= fBufferSize) {
322 BufferEmpty(1);
323 return Fill(x,y,z,w);
324 }
325 fBuffer[4*nbentries+1] = w;
326 fBuffer[4*nbentries+2] = x;
327 fBuffer[4*nbentries+3] = y;
328 fBuffer[4*nbentries+4] = z;
329 fBuffer[0] += 1;
330 return -3;
331}
332
333
334////////////////////////////////////////////////////////////////////////////////
335/// Invalid Fill method
336
338{
339 Error("Fill", "Invalid signature - do nothing");
340 return -1;
341}
342
343
344////////////////////////////////////////////////////////////////////////////////
345/// Increment cell defined by x,y,z by 1 .
346///
347/// The function returns the corresponding global bin number which has its content
348/// incremented by 1
349
351{
352 if (fBuffer) return BufferFill(x,y,z,1);
353
354 Int_t binx, biny, binz, bin;
355 fEntries++;
356 binx = fXaxis.FindBin(x);
357 biny = fYaxis.FindBin(y);
358 binz = fZaxis.FindBin(z);
359 if (binx <0 || biny <0 || binz<0) return -1;
360 bin = binx + (fXaxis.GetNbins()+2)*(biny + (fYaxis.GetNbins()+2)*binz);
361 if (fSumw2.fN) ++fSumw2.fArray[bin];
362 AddBinContent(bin);
363 if (binx == 0 || binx > fXaxis.GetNbins()) {
364 if (!GetStatOverflowsBehaviour()) return -1;
365 }
366
367 if (biny == 0 || biny > fYaxis.GetNbins()) {
368 if (!GetStatOverflowsBehaviour()) return -1;
369 }
370 if (binz == 0 || binz > fZaxis.GetNbins()) {
371 if (!GetStatOverflowsBehaviour()) return -1;
372 }
373 ++fTsumw;
374 ++fTsumw2;
375 fTsumwx += x;
376 fTsumwx2 += x*x;
377 fTsumwy += y;
378 fTsumwy2 += y*y;
379 fTsumwxy += x*y;
380 fTsumwz += z;
381 fTsumwz2 += z*z;
382 fTsumwxz += x*z;
383 fTsumwyz += y*z;
384 return bin;
385}
386
387
388////////////////////////////////////////////////////////////////////////////////
389/// Increment cell defined by x,y,z by a weight w.
390///
391/// If the weight is not equal to 1, the storage of the sum of squares of
392/// weights is automatically triggered and the sum of the squares of weights is incremented
393/// by w^2 in the cell corresponding to x,y,z.
394///
395/// The function returns the corresponding global bin number which has its content
396/// incremented by w
397
399{
400 if (fBuffer) return BufferFill(x,y,z,w);
401
402 Int_t binx, biny, binz, bin;
403 fEntries++;
404 binx = fXaxis.FindBin(x);
405 biny = fYaxis.FindBin(y);
406 binz = fZaxis.FindBin(z);
407 if (binx <0 || biny <0 || binz<0) return -1;
408 bin = binx + (fXaxis.GetNbins()+2)*(biny + (fYaxis.GetNbins()+2)*binz);
409 if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW)) Sumw2(); // must be called before AddBinContent
410 if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
411 AddBinContent(bin,w);
412 if (binx == 0 || binx > fXaxis.GetNbins()) {
413 if (!GetStatOverflowsBehaviour()) return -1;
414 }
415 if (biny == 0 || biny > fYaxis.GetNbins()) {
416 if (!GetStatOverflowsBehaviour()) return -1;
417 }
418 if (binz == 0 || binz > fZaxis.GetNbins()) {
419 if (!GetStatOverflowsBehaviour()) return -1;
420 }
421 fTsumw += w;
422 fTsumw2 += w*w;
423 fTsumwx += w*x;
424 fTsumwx2 += w*x*x;
425 fTsumwy += w*y;
426 fTsumwy2 += w*y*y;
427 fTsumwxy += w*x*y;
428 fTsumwz += w*z;
429 fTsumwz2 += w*z*z;
430 fTsumwxz += w*x*z;
431 fTsumwyz += w*y*z;
432 return bin;
433}
434
435
436////////////////////////////////////////////////////////////////////////////////
437/// Increment cell defined by namex,namey,namez by a weight w
438///
439/// If the weight is not equal to 1, the storage of the sum of squares of
440/// weights is automatically triggered and the sum of the squares of weights is incremented
441/// by w^2 in the corresponding cell.
442/// The function returns the corresponding global bin number which has its content
443/// incremented by w
444
445Int_t TH3::Fill(const char *namex, const char *namey, const char *namez, Double_t w)
446{
447 Int_t binx, biny, binz, bin;
448 fEntries++;
449 binx = fXaxis.FindBin(namex);
450 biny = fYaxis.FindBin(namey);
451 binz = fZaxis.FindBin(namez);
452 if (binx <0 || biny <0 || binz<0) return -1;
453 bin = binx + (fXaxis.GetNbins()+2)*(biny + (fYaxis.GetNbins()+2)*binz);
454 if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW)) Sumw2(); // must be called before AddBinContent
455 if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
456 AddBinContent(bin,w);
457 if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
458 if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
459 if (binz == 0 || binz > fZaxis.GetNbins()) return -1;
460
461 Double_t v = w;
462 fTsumw += v;
463 fTsumw2 += v*v;
464 // skip computation of the statistics along axis that have labels (can be extended and are aphanumeric)
465 UInt_t labelBitMask = GetAxisLabelStatus();
466 if (labelBitMask != TH1::kAllAxes) {
467 Double_t x = (labelBitMask & TH1::kXaxis) ? 0 : fXaxis.GetBinCenter(binx);
468 Double_t y = (labelBitMask & TH1::kYaxis) ? 0 : fYaxis.GetBinCenter(biny);
469 Double_t z = (labelBitMask & TH1::kZaxis) ? 0 : fZaxis.GetBinCenter(binz);
470 fTsumwx += v * x;
471 fTsumwx2 += v * x * x;
472 fTsumwy += v * y;
473 fTsumwy2 += v * y * y;
474 fTsumwxy += v * x * y;
475 fTsumwz += v * z;
476 fTsumwz2 += v * z * z;
477 fTsumwxz += v * x * z;
478 fTsumwyz += v * y * z;
479 }
480 return bin;
481}
482
483
484////////////////////////////////////////////////////////////////////////////////
485/// Increment cell defined by namex,y,namez by a weight w
486///
487/// If the weight is not equal to 1, the storage of the sum of squares of
488/// weights is automatically triggered and the sum of the squares of weights is incremented
489/// by w^2 in the corresponding cell.
490/// The function returns the corresponding global bin number which has its content
491/// incremented by w
492
493Int_t TH3::Fill(const char *namex, Double_t y, const char *namez, Double_t w)
494{
495 Int_t binx, biny, binz, bin;
496 fEntries++;
497 binx = fXaxis.FindBin(namex);
498 biny = fYaxis.FindBin(y);
499 binz = fZaxis.FindBin(namez);
500 if (binx <0 || biny <0 || binz<0) return -1;
501 bin = binx + (fXaxis.GetNbins()+2)*(biny + (fYaxis.GetNbins()+2)*binz);
502 if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW)) Sumw2(); // must be called before AddBinContent
503 if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
504 AddBinContent(bin,w);
505 if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
506 if (biny == 0 || biny > fYaxis.GetNbins()) {
507 if (!GetStatOverflowsBehaviour()) return -1;
508 }
509 if (binz == 0 || binz > fZaxis.GetNbins()) return -1;
510 Double_t v = w;
511 fTsumw += v;
512 fTsumw2 += v*v;
513 fTsumwy += v*y;
514 fTsumwy2 += v*y*y;
515 // skip computation of the statistics along axis that have labels (can be extended and are aphanumeric)
516 UInt_t labelBitMask = GetAxisLabelStatus();
517 if (labelBitMask != (TH1::kXaxis | TH1::kZaxis) ) {
518 Double_t x = (labelBitMask & TH1::kXaxis) ? 0 : fXaxis.GetBinCenter(binx);
519 Double_t z = (labelBitMask & TH1::kZaxis) ? 0 : fZaxis.GetBinCenter(binz);
520 fTsumwx += v * x;
521 fTsumwx2 += v * x * x;
522 fTsumwxy += v * x * y;
523 fTsumwz += v * z;
524 fTsumwz2 += v * z * z;
525 fTsumwxz += v * x * z;
526 fTsumwyz += v * y * z;
527 }
528 return bin;
529}
530
531
532////////////////////////////////////////////////////////////////////////////////
533/// Increment cell defined by namex,namey,z by a weight w
534///
535/// If the weight is not equal to 1, the storage of the sum of squares of
536/// weights is automatically triggered and the sum of the squares of weights is incremented
537/// by w^2 in the corresponding cell.
538/// The function returns the corresponding global bin number which has its content
539/// incremented by w
540
541Int_t TH3::Fill(const char *namex, const char *namey, Double_t z, Double_t w)
542{
543 Int_t binx, biny, binz, bin;
544 fEntries++;
545 binx = fXaxis.FindBin(namex);
546 biny = fYaxis.FindBin(namey);
547 binz = fZaxis.FindBin(z);
548 if (binx <0 || biny <0 || binz<0) return -1;
549 bin = binx + (fXaxis.GetNbins()+2)*(biny + (fYaxis.GetNbins()+2)*binz);
550 if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW)) Sumw2(); // must be called before AddBinContent
551 if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
552 AddBinContent(bin,w);
553 if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
554 if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
555 if (binz == 0 || binz > fZaxis.GetNbins()) {
556 if (!GetStatOverflowsBehaviour()) return -1;
557 }
558 Double_t v = w;
559 fTsumw += v;
560 fTsumw2 += v*v;
561 fTsumwz += v*z;
562 fTsumwz2 += v*z*z;
563 // skip computation of the statistics along axis that have labels (can be extended and are aphanumeric)
564 UInt_t labelBitMask = GetAxisLabelStatus();
565 if (labelBitMask != (TH1::kXaxis | TH1::kYaxis)) {
566 Double_t x = (labelBitMask & TH1::kXaxis) ? 0 : fXaxis.GetBinCenter(binx);
567 Double_t y = (labelBitMask & TH1::kYaxis) ? 0 : fYaxis.GetBinCenter(biny);
568 fTsumwx += v * x;
569 fTsumwx2 += v * x * x;
570 fTsumwy += v * y;
571 fTsumwy2 += v * y * y;
572 fTsumwxy += v * x * y;
573 fTsumwxz += v * x * z;
574 fTsumwyz += v * y * z;
575 }
576 return bin;
577}
578
579
580////////////////////////////////////////////////////////////////////////////////
581/// Increment cell defined by x,namey,namez by a weight w
582///
583/// If the weight is not equal to 1, the storage of the sum of squares of
584/// weights is automatically triggered and the sum of the squares of weights is incremented
585/// by w^2 in the corresponding cell.
586/// The function returns the corresponding global bin number which has its content
587/// incremented by w
588
589Int_t TH3::Fill(Double_t x, const char *namey, const char *namez, Double_t w)
590{
591 Int_t binx, biny, binz, bin;
592 fEntries++;
593 binx = fXaxis.FindBin(x);
594 biny = fYaxis.FindBin(namey);
595 binz = fZaxis.FindBin(namez);
596 if (binx <0 || biny <0 || binz<0) return -1;
597 bin = binx + (fXaxis.GetNbins()+2)*(biny + (fYaxis.GetNbins()+2)*binz);
598 if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW)) Sumw2(); // must be called before AddBinContent
599 if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
600 AddBinContent(bin,w);
601 if (binx == 0 || binx > fXaxis.GetNbins()) {
602 if (!GetStatOverflowsBehaviour()) return -1;
603 }
604 if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
605 if (binz == 0 || binz > fZaxis.GetNbins()) return -1;
606
607 // skip computation of the statistics along axis that have labels (can be extended and are aphanumeric)
608 UInt_t labelBitMask = GetAxisLabelStatus();
609 Double_t y = (labelBitMask & TH1::kYaxis) ? 0 : fYaxis.GetBinCenter(biny);
610 Double_t z = (labelBitMask & TH1::kZaxis) ? 0 : fZaxis.GetBinCenter(binz);
611 Double_t v = w;
612 fTsumw += v;
613 fTsumw2 += v * v;
614 fTsumwx += v * x;
615 fTsumwx2 += v * x * x;
616 fTsumwy += v * y;
617 fTsumwy2 += v * y * y;
618 fTsumwxy += v * x * y;
619 fTsumwz += v * z;
620 fTsumwz2 += v * z * z;
621 fTsumwxz += v * x * z;
622 fTsumwyz += v * y * z;
623 return bin;
624}
625
626////////////////////////////////////////////////////////////////////////////////
627/// Increment cell defined by namex , y ,z by a weight w
628///
629/// If the weight is not equal to 1, the storage of the sum of squares of
630/// weights is automatically triggered and the sum of the squares of weights is incremented
631/// by w^2 in the corresponding cell.
632/// The function returns the corresponding global bin number which has its content
633/// incremented by w
634
635Int_t TH3::Fill(const char * namex, Double_t y, Double_t z, Double_t w)
636{
637 Int_t binx, biny, binz, bin;
638 fEntries++;
639 binx = fXaxis.FindBin(namex);
640 biny = fYaxis.FindBin(y);
641 binz = fZaxis.FindBin(z);
642 if (binx < 0 || biny < 0 || binz < 0)
643 return -1;
644 bin = binx + (fXaxis.GetNbins() + 2) * (biny + (fYaxis.GetNbins() + 2) * binz);
645 if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW))
646 Sumw2(); // must be called before AddBinContent
647 if (fSumw2.fN)
648 fSumw2.fArray[bin] += w * w;
649 AddBinContent(bin, w);
650 if (binx == 0 || binx > fXaxis.GetNbins()) {
651 return -1;
652 }
653 if (biny == 0 || biny > fYaxis.GetNbins()) {
655 return -1;
656 }
657 if (binz == 0 || binz > fZaxis.GetNbins()) {
659 return -1;
660 }
661 UInt_t labelBitMask = GetAxisLabelStatus();
662 Double_t x = (labelBitMask & TH1::kXaxis) ? 0 : fXaxis.GetBinCenter(binx);
663 Double_t v = w;
664 fTsumw += v;
665 fTsumw2 += v * v;
666 fTsumwx += v * x;
667 fTsumwx2 += v * x * x;
668 fTsumwy += v * y;
669 fTsumwy2 += v * y * y;
670 fTsumwxy += v * x * y;
671 fTsumwz += v * z;
672 fTsumwz2 += v * z * z;
673 fTsumwxz += v * x * z;
674 fTsumwyz += v * y * z;
675 return bin;
676}
677
678////////////////////////////////////////////////////////////////////////////////
679/// Increment cell defined by x,namey,z by a weight w
680///
681/// If the weight is not equal to 1, the storage of the sum of squares of
682/// weights is automatically triggered and the sum of the squares of weights is incremented
683/// by w^2 in the corresponding cell.
684/// The function returns the corresponding global bin number which has its content
685/// incremented by w
686
687Int_t TH3::Fill(Double_t x, const char *namey, Double_t z, Double_t w)
688{
689 Int_t binx, biny, binz, bin;
690 fEntries++;
691 binx = fXaxis.FindBin(x);
692 biny = fYaxis.FindBin(namey);
693 binz = fZaxis.FindBin(z);
694 if (binx <0 || biny <0 || binz<0) return -1;
695 bin = binx + (fXaxis.GetNbins()+2)*(biny + (fYaxis.GetNbins()+2)*binz);
696 if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW)) Sumw2(); // must be called before AddBinContent
697 if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
698 AddBinContent(bin,w);
699 if (binx == 0 || binx > fXaxis.GetNbins()) {
700 if (!GetStatOverflowsBehaviour()) return -1;
701 }
702 if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
703 if (binz == 0 || binz > fZaxis.GetNbins()) {
704 if (!GetStatOverflowsBehaviour()) return -1;
705 }
706 UInt_t labelBitMask = GetAxisLabelStatus();
707 Double_t y = (labelBitMask & TH1::kYaxis) ? 0 : fYaxis.GetBinCenter(biny);
708 Double_t v = w;
709 fTsumw += v;
710 fTsumw2 += v*v;
711 fTsumwx += v*x;
712 fTsumwx2 += v*x*x;
713 fTsumwy += v*y;
714 fTsumwy2 += v*y*y;
715 fTsumwxy += v*x*y;
716 fTsumwz += v*z;
717 fTsumwz2 += v*z*z;
718 fTsumwxz += v*x*z;
719 fTsumwyz += v*y*z;
720 return bin;
721}
722
723
724////////////////////////////////////////////////////////////////////////////////
725/// Increment cell defined by x,y,namez by a weight w
726///
727/// If the weight is not equal to 1, the storage of the sum of squares of
728/// weights is automatically triggered and the sum of the squares of weights is incremented
729/// by w^2 in the corresponding cell.
730/// The function returns the corresponding global bin number which has its content
731/// incremented by w
732
733Int_t TH3::Fill(Double_t x, Double_t y, const char *namez, Double_t w)
734{
735 Int_t binx, biny, binz, bin;
736 fEntries++;
737 binx = fXaxis.FindBin(x);
738 biny = fYaxis.FindBin(y);
739 binz = fZaxis.FindBin(namez);
740 if (binx <0 || biny <0 || binz<0) return -1;
741 bin = binx + (fXaxis.GetNbins()+2)*(biny + (fYaxis.GetNbins()+2)*binz);
742 if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW)) Sumw2(); // must be called before AddBinContent
743 if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
744 AddBinContent(bin,w);
745 if (binx == 0 || binx > fXaxis.GetNbins()) {
746 if (!GetStatOverflowsBehaviour()) return -1;
747 }
748 if (biny == 0 || biny > fYaxis.GetNbins()) {
749 if (!GetStatOverflowsBehaviour()) return -1;
750 }
751 if (binz == 0 || binz > fZaxis.GetNbins()) return -1;
752 UInt_t labelBitMask = GetAxisLabelStatus();
753 Double_t z = (labelBitMask & TH1::kZaxis) ? 0 : fZaxis.GetBinCenter(binz);
754 Double_t v = w;
755 fTsumw += v;
756 fTsumw2 += v*v;
757 fTsumwx += v*x;
758 fTsumwx2 += v*x*x;
759 fTsumwy += v*y;
760 fTsumwy2 += v*y*y;
761 fTsumwxy += v*x*y;
762 fTsumwz += v*z;
763 fTsumwz2 += v*z*z;
764 fTsumwxz += v*x*z;
765 fTsumwyz += v*y*z;
766 return bin;
767}
768
769
770////////////////////////////////////////////////////////////////////////////////
771/// Fill histogram following distribution in function fname.
772///
773/// @param fname : Function name used for filling the historam
774/// @param ntimes : number of times the histogram is filled
775/// @param rng : (optional) Random number generator used to sample
776///
777/// The distribution contained in the function fname (TF1) is integrated
778/// over the channel contents.
779/// It is normalized to 1.
780/// Getting one random number implies:
781/// - Generating a random number between 0 and 1 (say r1)
782/// - Look in which bin in the normalized integral r1 corresponds to
783/// - Fill histogram channel
784/// ntimes random numbers are generated
785///
786/// N.B. By dfault this methods approximates the integral of the function in each bin with the
787/// function value at the center of the bin, mutiplied by the bin width
788///
789/// One can also call TF1::GetRandom to get a random variate from a function.
790
791void TH3::FillRandom(const char *fname, Int_t ntimes, TRandom * rng)
792{
793 Int_t bin, binx, biny, binz, ibin, loop;
794 Double_t r1, x, y,z, xv[3];
795 // Search for fname in the list of ROOT defined functions
796 TObject *fobj = gROOT->GetFunction(fname);
797 if (!fobj) { Error("FillRandom", "Unknown function: %s",fname); return; }
798 TF3 *f1 = dynamic_cast<TF3*>( fobj );
799 if (!f1) { Error("FillRandom", "Function: %s is not a TF3, is a %s",fname,fobj->IsA()->GetName()); return; }
800
801 TAxis & xAxis = fXaxis;
802 TAxis & yAxis = fYaxis;
803 TAxis & zAxis = fZaxis;
804
805 // in case axes of histogram are not defined use the function axis
807 Double_t xmin,xmax,ymin,ymax,zmin,zmax;
808 f1->GetRange(xmin,ymin,zmin,xmax,ymax,zmax);
809 Info("FillRandom","Using function axis and range ([%g,%g],[%g,%g],[%g,%g])",xmin, xmax,ymin,ymax,zmin,zmax);
810 xAxis = *(f1->GetHistogram()->GetXaxis());
811 yAxis = *(f1->GetHistogram()->GetYaxis());
812 zAxis = *(f1->GetHistogram()->GetZaxis());
813 }
814
815 // Allocate temporary space to store the integral and compute integral
816 Int_t nbinsx = xAxis.GetNbins();
817 Int_t nbinsy = yAxis.GetNbins();
818 Int_t nbinsz = zAxis.GetNbins();
819 Int_t nxy = nbinsx*nbinsy;
820 Int_t nbins = nbinsx*nbinsy*nbinsz;
821
822 Double_t *integral = new Double_t[nbins+1];
823 ibin = 0;
824 integral[ibin] = 0;
825 // approximate integral with function value at bin center
826 for (binz=1;binz<=nbinsz;binz++) {
827 xv[2] = zAxis.GetBinCenter(binz);
828 for (biny=1;biny<=nbinsy;biny++) {
829 xv[1] = yAxis.GetBinCenter(biny);
830 for (binx=1;binx<=nbinsx;binx++) {
831 xv[0] = xAxis.GetBinCenter(binx);
832 ibin++;
833 Double_t fint = f1->EvalPar(xv, nullptr);
834 // uncomment this line to have the integral computation in a bin
835 // Double_t fint = f1->Integral(xAxis.GetBinLowEdge(binx), xAxis.GetBinUpEdge(binx),
836 // yAxis.GetBinLowEdge(biny), yAxis.GetBinUpEdge(biny),
837 // zAxis.GetBinLowEdge(binz), zAxis.GetBinUpEdge(binz));
838 integral[ibin] = integral[ibin-1] + fint;
839 }
840 }
841 }
842
843 // Normalize integral to 1
844 if (integral[nbins] == 0 ) {
845 delete [] integral;
846 Error("FillRandom", "Integral = zero"); return;
847 }
848 for (bin=1;bin<=nbins;bin++) integral[bin] /= integral[nbins];
849
850 // Start main loop ntimes
851 if (fDimension < 2) nbinsy = -1;
852 if (fDimension < 3) nbinsz = -1;
853 for (loop=0;loop<ntimes;loop++) {
854 r1 = (rng) ? rng->Rndm() : gRandom->Rndm();
855 ibin = TMath::BinarySearch(nbins,&integral[0],r1);
856 binz = ibin/nxy;
857 biny = (ibin - nxy*binz)/nbinsx;
858 binx = 1 + ibin - nbinsx*(biny + nbinsy*binz);
859 if (nbinsz) binz++;
860 if (nbinsy) biny++;
861 x = xAxis.GetBinCenter(binx);
862 y = yAxis.GetBinCenter(biny);
863 z = zAxis.GetBinCenter(binz);
864 Fill(x,y,z, 1.);
865 }
866 delete [] integral;
867}
868
869
870////////////////////////////////////////////////////////////////////////////////
871/// Fill histogram following distribution in histogram h.
872///
873/// @param h : Histogram pointer used for smpling random number
874/// @param ntimes : number of times the histogram is filled
875/// @param rng : (optional) Random number generator used for sampling
876///
877/// The distribution contained in the histogram h (TH3) is integrated
878/// over the channel contents.
879/// It is normalized to 1.
880/// Getting one random number implies:
881/// - Generating a random number between 0 and 1 (say r1)
882/// - Look in which bin in the normalized integral r1 corresponds to
883/// - Fill histogram channel
884/// ntimes random numbers are generated
885
886void TH3::FillRandom(TH1 *h, Int_t ntimes, TRandom * rng)
887{
888 if (!h) { Error("FillRandom", "Null histogram"); return; }
889 if (fDimension != h->GetDimension()) {
890 Error("FillRandom", "Histograms with different dimensions"); return;
891 }
892
893 if (h->ComputeIntegral() == 0) return;
894
895 TH3 *h3 = (TH3*)h;
896 Int_t loop;
897 Double_t x,y,z;
898 for (loop=0;loop<ntimes;loop++) {
899 h3->GetRandom3(x,y,z,rng);
900 Fill(x,y,z);
901 }
902}
903
904////////////////////////////////////////////////////////////////////////////////
905/// Project slices along Z in case of a 3-D histogram, then fit each slice
906/// with function f1 and make a 2-d histogram for each fit parameter
907/// Only cells in the bin range [binminx,binmaxx] and [binminy,binmaxy] are considered.
908/// if f1=0, a gaussian is assumed
909/// Before invoking this function, one can set a subrange to be fitted along Z
910/// via f1->SetRange(zmin,zmax)
911/// The argument option (default="QNR") can be used to change the fit options.
912/// "Q" means Quiet mode
913/// "N" means do not show the result of the fit
914/// "R" means fit the function in the specified function range
915///
916/// Note that the generated histograms are added to the list of objects
917/// in the current directory. It is the user's responsibility to delete
918/// these histograms.
919///
920/// Example: Assume a 3-d histogram h3
921/// Root > h3->FitSlicesZ(); produces 4 TH2D histograms
922/// with h3_0 containing parameter 0(Constant) for a Gaus fit
923/// of each cell in X,Y projected along Z
924/// with h3_1 containing parameter 1(Mean) for a gaus fit
925/// with h3_2 containing parameter 2(StdDev) for a gaus fit
926/// with h3_chi2 containing the chisquare/number of degrees of freedom for a gaus fit
927///
928/// Root > h3->Fit(0,15,22,0,0,10);
929/// same as above, but only for bins 15 to 22 along X
930/// and only for cells in X,Y for which the corresponding projection
931/// along Z has more than cut bins filled.
932///
933/// NOTE: To access the generated histograms in the current directory, do eg:
934/// TH2D *h3_1 = (TH2D*)gDirectory->Get("h3_1");
935
936void TH3::FitSlicesZ(TF1 *f1, Int_t binminx, Int_t binmaxx, Int_t binminy, Int_t binmaxy, Int_t cut, Option_t *option)
937{
938 //Int_t nbinsz = fZaxis.GetNbins();
939
940 // get correct first and last bins for outer axes used in the loop doing the slices
941 // when using default values (0,-1) check if an axis range is set in outer axis
942 // do same as in DoProjection for inner axis
943 auto computeFirstAndLastBin = [](const TAxis & outerAxis, Int_t &firstbin, Int_t &lastbin) {
944 Int_t nbins = outerAxis.GetNbins();
945 if ( lastbin < firstbin && outerAxis.TestBit(TAxis::kAxisRange) ) {
946 firstbin = outerAxis.GetFirst();
947 lastbin = outerAxis.GetLast();
948 // For special case of TAxis::SetRange, when first == 1 and last
949 // = N and the range bit has been set, the TAxis will return 0
950 // for both.
951 if (firstbin == 0 && lastbin == 0) {
952 firstbin = 1;
953 lastbin = nbins;
954 }
955 }
956 if (firstbin < 0) firstbin = 0;
957 if (lastbin < 0 || lastbin > nbins + 1) lastbin = nbins + 1;
958 if (lastbin < firstbin) {firstbin = 0; lastbin = nbins + 1;}
959 };
960
961 computeFirstAndLastBin(fXaxis, binminx, binmaxx);
962 computeFirstAndLastBin(fYaxis, binminy, binmaxy);
963
964 // limits for the axis of the fit results histograms are different
965 auto computeAxisLimits = [](const TAxis & outerAxis, Int_t firstbin, Int_t lastbin,
966 Int_t &nBins, Double_t &xMin, Double_t & xMax) {
967 Int_t firstOutBin = std::max(firstbin,1);
968 Int_t lastOutBin = std::min(lastbin,outerAxis.GetNbins() ) ;
969 nBins = lastOutBin-firstOutBin+1;
970 xMin = outerAxis.GetBinLowEdge(firstOutBin);
971 xMax = outerAxis.GetBinUpEdge(lastOutBin);
972 // return first bin that is used in case of variable bin size axis
973 return firstOutBin;
974 };
975 Int_t nbinsX = 0;
976 Double_t xMin, xMax = 0;
977 Int_t firstBinXaxis = computeAxisLimits(fXaxis, binminx, binmaxx, nbinsX, xMin, xMax);
978 Int_t nbinsY = 0;
979 Double_t yMin, yMax = 0;
980 Int_t firstBinYaxis = computeAxisLimits(fYaxis, binminy, binmaxy, nbinsY, yMin, yMax);
981
982 //default is to fit with a gaussian
983 if (f1 == 0) {
984 f1 = (TF1*)gROOT->GetFunction("gaus");
985 if (f1 == 0) f1 = new TF1("gaus","gaus",fZaxis.GetXmin(),fZaxis.GetXmax());
987 }
988 const char *fname = f1->GetName();
989 Int_t npar = f1->GetNpar();
990 Double_t *parsave = new Double_t[npar];
991 f1->GetParameters(parsave);
992
993 //Create one 2-d histogram for each function parameter
994 Int_t ipar;
996 TString title;
997 std::vector<TH1*> hlist(npar+1); // include also chi2 histogram
998 const TArrayD *xbins = fXaxis.GetXbins();
999 const TArrayD *ybins = fYaxis.GetXbins();
1000 for (ipar=0;ipar<= npar;ipar++) {
1001 if (ipar < npar) {
1002 // fitted parameter histograms
1003 name = TString::Format("%s_%d",GetName(),ipar);
1004 title = TString::Format("Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar));
1005 } else {
1006 // chi2 histograms
1007 name = TString::Format("%s_chi2",GetName());
1008 title = "chisquare";
1009 }
1010 if (xbins->fN == 0 && ybins->fN == 0) {
1011 hlist[ipar] = new TH2D(name, title,
1012 nbinsX, xMin, xMax,
1013 nbinsY, yMin, yMax);
1014 } else if (xbins->fN > 0 && ybins->fN > 0 ) {
1015 hlist[ipar] = new TH2D(name, title,
1016 nbinsX, &xbins->fArray[firstBinXaxis],
1017 nbinsY, &ybins->fArray[firstBinYaxis]);
1018 }
1019 // mixed case do not exist for TH3
1020 R__ASSERT(hlist[ipar]);
1021
1022 hlist[ipar]->GetXaxis()->SetTitle(fXaxis.GetTitle());
1023 hlist[ipar]->GetYaxis()->SetTitle(fYaxis.GetTitle());
1024 }
1025 TH1 * hchi2 = hlist.back();
1026
1027 //Loop on all cells in X,Y generate a projection along Z
1028 TH1D *hpz = nullptr;
1029 TString opt(option);
1030 // add option "N" when fitting the 2D histograms
1031 opt += " N ";
1032
1033 for (Int_t biny=binminy; biny<=binmaxy; biny++) {
1034 for (Int_t binx=binminx; binx<=binmaxx; binx++) {
1035 // use TH3::ProjectionZ
1036 hpz = ProjectionZ("R_temp",binx,binx,biny,biny);
1037
1038 Double_t nentries = hpz->GetEntries();
1039 if ( nentries <= 0 || nentries < cut) {
1040 if (!opt.Contains("Q"))
1041 Info("FitSlicesZ","Slice (%d,%d) skipped, the number of entries is zero or smaller than the given cut value, n=%f",binx,biny,nentries);
1042 continue;
1043 }
1044 f1->SetParameters(parsave);
1045 Int_t bin = hlist[0]->FindBin( fXaxis.GetBinCenter(binx), fYaxis.GetBinCenter(biny) );
1046 if (!opt.Contains("Q")) {
1047 int ibx,iby,ibz = 0;
1048 hlist[0]->GetBinXYZ(bin,ibx,iby,ibz);
1049 Info("DoFitSlices","Slice fit [(%f,%f),(%f,%f)]",hlist[0]->GetXaxis()->GetBinLowEdge(ibx), hlist[0]->GetXaxis()->GetBinUpEdge(ibx),
1050 hlist[0]->GetYaxis()->GetBinLowEdge(iby), hlist[0]->GetYaxis()->GetBinUpEdge(iby));
1051 }
1052 hpz->Fit(fname,opt.Data());
1053 Int_t npfits = f1->GetNumberFitPoints();
1054 if (npfits > npar && npfits >= cut) {
1055 for (ipar=0;ipar<npar;ipar++) {
1056 hlist[ipar]->SetBinContent(bin,f1->GetParameter(ipar));
1057 hlist[ipar]->SetBinError(bin,f1->GetParError(ipar));
1058 }
1059 hchi2->SetBinContent(bin,f1->GetChisquare()/(npfits-npar));
1060 }
1061 else {
1062 if (!opt.Contains("Q"))
1063 Info("FitSlicesZ","Fitted slice (%d,%d) skipped, the number of fitted points is too small, n=%d",binx,biny,npfits);
1064 }
1065 }
1066 }
1067 delete [] parsave;
1068 delete hpz;
1069}
1070
1071
1072////////////////////////////////////////////////////////////////////////////////
1073/// See comments in TH1::GetBin
1074
1075Int_t TH3::GetBin(Int_t binx, Int_t biny, Int_t binz) const
1076{
1077 Int_t ofy = fYaxis.GetNbins() + 1; // code duplication unavoidable because TH3 does not inherit from TH2
1078 if (biny < 0) biny = 0;
1079 if (biny > ofy) biny = ofy;
1080
1081 Int_t ofz = fZaxis.GetNbins() + 1; // overflow bin
1082 if (binz < 0) binz = 0;
1083 if (binz > ofz) binz = ofz;
1084
1085 return TH1::GetBin(binx) + (fXaxis.GetNbins() + 2) * (biny + (fYaxis.GetNbins() + 2) * binz);
1086}
1087
1088
1089////////////////////////////////////////////////////////////////////////////////
1090/// Compute first cell (binx,biny,binz) in the range [firstx,lastx](firsty,lasty][firstz,lastz] for which
1091/// diff = abs(cell_content-c) <= maxdiff
1092/// In case several cells in the specified range with diff=0 are found
1093/// the first cell found is returned in binx,biny,binz.
1094/// In case several cells in the specified range satisfy diff <=maxdiff
1095/// the cell with the smallest difference is returned in binx,biny,binz.
1096/// In all cases the function returns the smallest difference.
1097///
1098/// NOTE1: if firstx <= 0, firstx is set to bin 1
1099/// if (lastx < firstx then firstx is set to the number of bins in X
1100/// ie if firstx=0 and lastx=0 (default) the search is on all bins in X.
1101/// if firsty <= 0, firsty is set to bin 1
1102/// if (lasty < firsty then firsty is set to the number of bins in Y
1103/// ie if firsty=0 and lasty=0 (default) the search is on all bins in Y.
1104/// if firstz <= 0, firstz is set to bin 1
1105/// if (lastz < firstz then firstz is set to the number of bins in Z
1106/// ie if firstz=0 and lastz=0 (default) the search is on all bins in Z.
1107/// NOTE2: if maxdiff=0 (default), the first cell with content=c is returned.
1108
1110 Int_t firstx, Int_t lastx,
1111 Int_t firsty, Int_t lasty,
1112 Int_t firstz, Int_t lastz,
1113 Double_t maxdiff) const
1114{
1115 if (fDimension != 3) {
1116 binx = 0;
1117 biny = 0;
1118 binz = 0;
1119 Error("GetBinWithContent3","function is only valid for 3-D histograms");
1120 return 0;
1121 }
1122 if (firstx <= 0) firstx = 1;
1123 if (lastx < firstx) lastx = fXaxis.GetNbins();
1124 if (firsty <= 0) firsty = 1;
1125 if (lasty < firsty) lasty = fYaxis.GetNbins();
1126 if (firstz <= 0) firstz = 1;
1127 if (lastz < firstz) lastz = fZaxis.GetNbins();
1128 Int_t binminx = 0, binminy=0, binminz=0;
1129 Double_t diff, curmax = 1.e240;
1130 for (Int_t k=firstz;k<=lastz;k++) {
1131 for (Int_t j=firsty;j<=lasty;j++) {
1132 for (Int_t i=firstx;i<=lastx;i++) {
1133 diff = TMath::Abs(GetBinContent(i,j,k)-c);
1134 if (diff <= 0) {binx = i; biny=j; binz=k; return diff;}
1135 if (diff < curmax && diff <= maxdiff) {curmax = diff, binminx=i; binminy=j;binminz=k;}
1136 }
1137 }
1138 }
1139 binx = binminx;
1140 biny = binminy;
1141 binz = binminz;
1142 return curmax;
1143}
1144
1145
1146////////////////////////////////////////////////////////////////////////////////
1147/// Return correlation factor between axis1 and axis2.
1148
1150{
1151 if (axis1 < 1 || axis2 < 1 || axis1 > 3 || axis2 > 3) {
1152 Error("GetCorrelationFactor","Wrong parameters");
1153 return 0;
1154 }
1155 if (axis1 == axis2) return 1;
1156 Double_t stddev1 = GetStdDev(axis1);
1157 if (stddev1 == 0) return 0;
1158 Double_t stddev2 = GetStdDev(axis2);
1159 if (stddev2 == 0) return 0;
1160 return GetCovariance(axis1,axis2)/stddev1/stddev2;
1161}
1162
1163
1164////////////////////////////////////////////////////////////////////////////////
1165/// Return covariance between axis1 and axis2.
1166
1168{
1169 if (axis1 < 1 || axis2 < 1 || axis1 > 3 || axis2 > 3) {
1170 Error("GetCovariance","Wrong parameters");
1171 return 0;
1172 }
1173 Double_t stats[kNstat];
1174 GetStats(stats);
1175 Double_t sumw = stats[0];
1176 Double_t sumwx = stats[2];
1177 Double_t sumwx2 = stats[3];
1178 Double_t sumwy = stats[4];
1179 Double_t sumwy2 = stats[5];
1180 Double_t sumwxy = stats[6];
1181 Double_t sumwz = stats[7];
1182 Double_t sumwz2 = stats[8];
1183 Double_t sumwxz = stats[9];
1184 Double_t sumwyz = stats[10];
1185
1186 if (sumw == 0) return 0;
1187 if (axis1 == 1 && axis2 == 1) {
1188 return TMath::Abs(sumwx2/sumw - sumwx*sumwx/(sumw*sumw));
1189 }
1190 if (axis1 == 2 && axis2 == 2) {
1191 return TMath::Abs(sumwy2/sumw - sumwy*sumwy/(sumw*sumw));
1192 }
1193 if (axis1 == 3 && axis2 == 3) {
1194 return TMath::Abs(sumwz2/sumw - sumwz*sumwz/(sumw*sumw));
1195 }
1196 if ((axis1 == 1 && axis2 == 2) || (axis1 == 2 && axis2 == 1)) {
1197 return sumwxy/sumw - sumwx*sumwy/(sumw*sumw);
1198 }
1199 if ((axis1 == 1 && axis2 == 3) || (axis1 == 3 && axis2 == 1)) {
1200 return sumwxz/sumw - sumwx*sumwz/(sumw*sumw);
1201 }
1202 if ((axis1 == 2 && axis2 == 3) || (axis1 == 3 && axis2 == 2)) {
1203 return sumwyz/sumw - sumwy*sumwz/(sumw*sumw);
1204 }
1205 return 0;
1206}
1207
1208
1209////////////////////////////////////////////////////////////////////////////////
1210/// Return 3 random numbers along axis x , y and z distributed according
1211/// to the cell-contents of this 3-dim histogram
1212/// @param[out] x reference to random generated x value
1213/// @param[out] y reference to random generated y value
1214/// @param[out] z reference to random generated z value
1215/// @param[in] rng (optional) Random number generator pointer used (default is gRandom)
1216
1218{
1219 Int_t nbinsx = GetNbinsX();
1220 Int_t nbinsy = GetNbinsY();
1221 Int_t nbinsz = GetNbinsZ();
1222 Int_t nxy = nbinsx*nbinsy;
1223 Int_t nbins = nxy*nbinsz;
1224 Double_t integral;
1225 // compute integral checking that all bins have positive content (see ROOT-5894)
1226 if (fIntegral) {
1227 if (fIntegral[nbins+1] != fEntries) integral = ComputeIntegral(true);
1228 else integral = fIntegral[nbins];
1229 } else {
1230 integral = ComputeIntegral(true);
1231 }
1232 if (integral == 0 ) { x = 0; y = 0; z = 0; return;}
1233 // case histogram has negative bins
1234 if (integral == TMath::QuietNaN() ) { x = TMath::QuietNaN(); y = TMath::QuietNaN(); z = TMath::QuietNaN(); return;}
1235
1236 if (!rng) rng = gRandom;
1237 Double_t r1 = rng->Rndm();
1238 Int_t ibin = TMath::BinarySearch(nbins,fIntegral,(Double_t) r1);
1239 Int_t binz = ibin/nxy;
1240 Int_t biny = (ibin - nxy*binz)/nbinsx;
1241 Int_t binx = ibin - nbinsx*(biny + nbinsy*binz);
1242 x = fXaxis.GetBinLowEdge(binx+1);
1243 if (r1 > fIntegral[ibin]) x +=
1244 fXaxis.GetBinWidth(binx+1)*(r1-fIntegral[ibin])/(fIntegral[ibin+1] - fIntegral[ibin]);
1245 y = fYaxis.GetBinLowEdge(biny+1) + fYaxis.GetBinWidth(biny+1)*rng->Rndm();
1246 z = fZaxis.GetBinLowEdge(binz+1) + fZaxis.GetBinWidth(binz+1)*rng->Rndm();
1247}
1248
1249
1250////////////////////////////////////////////////////////////////////////////////
1251/// Fill the array stats from the contents of this histogram
1252/// The array stats must be correctly dimensioned in the calling program.
1253/// stats[0] = sumw
1254/// stats[1] = sumw2
1255/// stats[2] = sumwx
1256/// stats[3] = sumwx2
1257/// stats[4] = sumwy
1258/// stats[5] = sumwy2
1259/// stats[6] = sumwxy
1260/// stats[7] = sumwz
1261/// stats[8] = sumwz2
1262/// stats[9] = sumwxz
1263/// stats[10]= sumwyz
1264
1265void TH3::GetStats(Double_t *stats) const
1266{
1267 if (fBuffer) ((TH3*)this)->BufferEmpty();
1268
1269 Int_t bin, binx, biny, binz;
1270 Double_t w,err;
1271 Double_t x,y,z;
1273 for (bin=0;bin<11;bin++) stats[bin] = 0;
1274
1275 Int_t firstBinX = fXaxis.GetFirst();
1276 Int_t lastBinX = fXaxis.GetLast();
1277 Int_t firstBinY = fYaxis.GetFirst();
1278 Int_t lastBinY = fYaxis.GetLast();
1279 Int_t firstBinZ = fZaxis.GetFirst();
1280 Int_t lastBinZ = fZaxis.GetLast();
1281 // include underflow/overflow if TH1::StatOverflows(kTRUE) in case no range is set on the axis
1284 if (firstBinX == 1) firstBinX = 0;
1285 if (lastBinX == fXaxis.GetNbins() ) lastBinX += 1;
1286 }
1288 if (firstBinY == 1) firstBinY = 0;
1289 if (lastBinY == fYaxis.GetNbins() ) lastBinY += 1;
1290 }
1292 if (firstBinZ == 1) firstBinZ = 0;
1293 if (lastBinZ == fZaxis.GetNbins() ) lastBinZ += 1;
1294 }
1295 }
1296
1297 // check for labels axis . In that case corresponsing statistics do not make sense and it is set to zero
1298 Bool_t labelXaxis = ((const_cast<TAxis&>(fXaxis)).GetLabels() && fXaxis.CanExtend() );
1299 Bool_t labelYaxis = ((const_cast<TAxis&>(fYaxis)).GetLabels() && fYaxis.CanExtend() );
1300 Bool_t labelZaxis = ((const_cast<TAxis&>(fZaxis)).GetLabels() && fZaxis.CanExtend() );
1301
1302 for (binz = firstBinZ; binz <= lastBinZ; binz++) {
1303 z = (!labelZaxis) ? fZaxis.GetBinCenter(binz) : 0;
1304 for (biny = firstBinY; biny <= lastBinY; biny++) {
1305 y = (!labelYaxis) ? fYaxis.GetBinCenter(biny) : 0;
1306 for (binx = firstBinX; binx <= lastBinX; binx++) {
1307 bin = GetBin(binx,biny,binz);
1308 x = (!labelXaxis) ? fXaxis.GetBinCenter(binx) : 0;
1309 //w = TMath::Abs(GetBinContent(bin));
1310 w = RetrieveBinContent(bin);
1311 err = TMath::Abs(GetBinError(bin));
1312 stats[0] += w;
1313 stats[1] += err*err;
1314 stats[2] += w*x;
1315 stats[3] += w*x*x;
1316 stats[4] += w*y;
1317 stats[5] += w*y*y;
1318 stats[6] += w*x*y;
1319 stats[7] += w*z;
1320 stats[8] += w*z*z;
1321 stats[9] += w*x*z;
1322 stats[10]+= w*y*z;
1323 }
1324 }
1325 }
1326 } else {
1327 stats[0] = fTsumw;
1328 stats[1] = fTsumw2;
1329 stats[2] = fTsumwx;
1330 stats[3] = fTsumwx2;
1331 stats[4] = fTsumwy;
1332 stats[5] = fTsumwy2;
1333 stats[6] = fTsumwxy;
1334 stats[7] = fTsumwz;
1335 stats[8] = fTsumwz2;
1336 stats[9] = fTsumwxz;
1337 stats[10]= fTsumwyz;
1338 }
1339}
1340
1341
1342////////////////////////////////////////////////////////////////////////////////
1343/// Return integral of bin contents. Only bins in the bins range are considered.
1344/// By default the integral is computed as the sum of bin contents in the range.
1345/// if option "width" is specified, the integral is the sum of
1346/// the bin contents multiplied by the bin width in x, y and in z.
1347
1349{
1352 fZaxis.GetFirst(),fZaxis.GetLast(),option);
1353}
1354
1355
1356////////////////////////////////////////////////////////////////////////////////
1357/// Return integral of bin contents in range [binx1,binx2],[biny1,biny2],[binz1,binz2]
1358/// for a 3-D histogram
1359/// By default the integral is computed as the sum of bin contents in the range.
1360/// if option "width" is specified, the integral is the sum of
1361/// the bin contents multiplied by the bin width in x, y and in z.
1362
1363Double_t TH3::Integral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2,
1364 Int_t binz1, Int_t binz2, Option_t *option) const
1365{
1366 Double_t err = 0;
1367 return DoIntegral(binx1,binx2,biny1,biny2,binz1,binz2,err,option);
1368}
1369
1370
1371////////////////////////////////////////////////////////////////////////////////
1372/// Return integral of bin contents in range [binx1,binx2],[biny1,biny2],[binz1,binz2]
1373/// for a 3-D histogram. Calculates also the integral error using error propagation
1374/// from the bin errors assuming that all the bins are uncorrelated.
1375/// By default the integral is computed as the sum of bin contents in the range.
1376/// if option "width" is specified, the integral is the sum of
1377/// the bin contents multiplied by the bin width in x, y and in z.
1378
1380 Int_t binz1, Int_t binz2,
1381 Double_t & error, Option_t *option) const
1382{
1383 return DoIntegral(binx1,binx2,biny1,biny2,binz1,binz2,error,option,kTRUE);
1384}
1385
1386////////////////////////////////////////////////////////////////////////////////
1387///Not yet implemented
1388
1390{
1391 Error("Interpolate","This function must be called with 3 arguments for a TH3");
1392 return 0;
1393}
1394
1395
1396////////////////////////////////////////////////////////////////////////////////
1397///Not yet implemented
1398
1400{
1401 Error("Interpolate","This function must be called with 3 arguments for a TH3");
1402 return 0;
1403}
1404
1405
1406////////////////////////////////////////////////////////////////////////////////
1407/// Given a point P(x,y,z), Interpolate approximates the value via trilinear interpolation
1408/// based on the 8 nearest bin center points (corner of the cube surrounding the points)
1409/// The Algorithm is described in http://en.wikipedia.org/wiki/Trilinear_interpolation
1410/// The given values (x,y,z) must be between first bin center and last bin center for each coordinate:
1411///
1412/// fXAxis.GetBinCenter(1) < x < fXaxis.GetBinCenter(nbinX) AND
1413/// fYAxis.GetBinCenter(1) < y < fYaxis.GetBinCenter(nbinY) AND
1414/// fZAxis.GetBinCenter(1) < z < fZaxis.GetBinCenter(nbinZ)
1415
1417{
1418 Int_t ubx = fXaxis.FindFixBin(x);
1419 if ( x < fXaxis.GetBinCenter(ubx) ) ubx -= 1;
1420 Int_t obx = ubx + 1;
1421
1422 Int_t uby = fYaxis.FindFixBin(y);
1423 if ( y < fYaxis.GetBinCenter(uby) ) uby -= 1;
1424 Int_t oby = uby + 1;
1425
1426 Int_t ubz = fZaxis.FindFixBin(z);
1427 if ( z < fZaxis.GetBinCenter(ubz) ) ubz -= 1;
1428 Int_t obz = ubz + 1;
1429
1430
1431// if ( IsBinUnderflow(GetBin(ubx, uby, ubz)) ||
1432// IsBinOverflow (GetBin(obx, oby, obz)) ) {
1433 if (ubx <=0 || uby <=0 || ubz <= 0 ||
1434 obx > fXaxis.GetNbins() || oby > fYaxis.GetNbins() || obz > fZaxis.GetNbins() ) {
1435 Error("Interpolate","Cannot interpolate outside histogram domain.");
1436 return 0;
1437 }
1438
1442
1443 Double_t xd = (x - fXaxis.GetBinCenter(ubx)) / xw;
1444 Double_t yd = (y - fYaxis.GetBinCenter(uby)) / yw;
1445 Double_t zd = (z - fZaxis.GetBinCenter(ubz)) / zw;
1446
1447
1448 Double_t v[] = { GetBinContent( ubx, uby, ubz ), GetBinContent( ubx, uby, obz ),
1449 GetBinContent( ubx, oby, ubz ), GetBinContent( ubx, oby, obz ),
1450 GetBinContent( obx, uby, ubz ), GetBinContent( obx, uby, obz ),
1451 GetBinContent( obx, oby, ubz ), GetBinContent( obx, oby, obz ) };
1452
1453
1454 Double_t i1 = v[0] * (1 - zd) + v[1] * zd;
1455 Double_t i2 = v[2] * (1 - zd) + v[3] * zd;
1456 Double_t j1 = v[4] * (1 - zd) + v[5] * zd;
1457 Double_t j2 = v[6] * (1 - zd) + v[7] * zd;
1458
1459
1460 Double_t w1 = i1 * (1 - yd) + i2 * yd;
1461 Double_t w2 = j1 * (1 - yd) + j2 * yd;
1462
1463
1464 Double_t result = w1 * (1 - xd) + w2 * xd;
1465
1466 return result;
1467}
1468
1469
1470////////////////////////////////////////////////////////////////////////////////
1471/// Statistical test of compatibility in shape between
1472/// THIS histogram and h2, using Kolmogorov test.
1473/// Default: Ignore under- and overflow bins in comparison
1474///
1475/// option is a character string to specify options
1476/// "U" include Underflows in test
1477/// "O" include Overflows
1478/// "N" include comparison of normalizations
1479/// "D" Put out a line of "Debug" printout
1480/// "M" Return the Maximum Kolmogorov distance instead of prob
1481///
1482/// The returned function value is the probability of test
1483/// (much less than one means NOT compatible)
1484///
1485/// The KS test uses the distance between the pseudo-CDF's obtained
1486/// from the histogram. Since in more than 1D the order for generating the pseudo-CDF is
1487/// arbitrary, we use the pseudo-CDF's obtained from all the possible 6 combinations of the 3 axis.
1488/// The average of all the maximum distances obtained is used in the tests.
1489
1490Double_t TH3::KolmogorovTest(const TH1 *h2, Option_t *option) const
1491{
1492 TString opt = option;
1493 opt.ToUpper();
1494
1495 Double_t prb = 0;
1496 TH1 *h1 = (TH1*)this;
1497 if (h2 == 0) return 0;
1498 const TAxis *xaxis1 = h1->GetXaxis();
1499 const TAxis *xaxis2 = h2->GetXaxis();
1500 const TAxis *yaxis1 = h1->GetYaxis();
1501 const TAxis *yaxis2 = h2->GetYaxis();
1502 const TAxis *zaxis1 = h1->GetZaxis();
1503 const TAxis *zaxis2 = h2->GetZaxis();
1504 Int_t ncx1 = xaxis1->GetNbins();
1505 Int_t ncx2 = xaxis2->GetNbins();
1506 Int_t ncy1 = yaxis1->GetNbins();
1507 Int_t ncy2 = yaxis2->GetNbins();
1508 Int_t ncz1 = zaxis1->GetNbins();
1509 Int_t ncz2 = zaxis2->GetNbins();
1510
1511 // Check consistency of dimensions
1512 if (h1->GetDimension() != 3 || h2->GetDimension() != 3) {
1513 Error("KolmogorovTest","Histograms must be 3-D\n");
1514 return 0;
1515 }
1516
1517 // Check consistency in number of channels
1518 if (ncx1 != ncx2) {
1519 Error("KolmogorovTest","Number of channels in X is different, %d and %d\n",ncx1,ncx2);
1520 return 0;
1521 }
1522 if (ncy1 != ncy2) {
1523 Error("KolmogorovTest","Number of channels in Y is different, %d and %d\n",ncy1,ncy2);
1524 return 0;
1525 }
1526 if (ncz1 != ncz2) {
1527 Error("KolmogorovTest","Number of channels in Z is different, %d and %d\n",ncz1,ncz2);
1528 return 0;
1529 }
1530
1531 // Check consistency in channel edges
1532 Bool_t afunc1 = kFALSE;
1533 Bool_t afunc2 = kFALSE;
1534 Double_t difprec = 1e-5;
1535 Double_t diff1 = TMath::Abs(xaxis1->GetXmin() - xaxis2->GetXmin());
1536 Double_t diff2 = TMath::Abs(xaxis1->GetXmax() - xaxis2->GetXmax());
1537 if (diff1 > difprec || diff2 > difprec) {
1538 Error("KolmogorovTest","histograms with different binning along X");
1539 return 0;
1540 }
1541 diff1 = TMath::Abs(yaxis1->GetXmin() - yaxis2->GetXmin());
1542 diff2 = TMath::Abs(yaxis1->GetXmax() - yaxis2->GetXmax());
1543 if (diff1 > difprec || diff2 > difprec) {
1544 Error("KolmogorovTest","histograms with different binning along Y");
1545 return 0;
1546 }
1547 diff1 = TMath::Abs(zaxis1->GetXmin() - zaxis2->GetXmin());
1548 diff2 = TMath::Abs(zaxis1->GetXmax() - zaxis2->GetXmax());
1549 if (diff1 > difprec || diff2 > difprec) {
1550 Error("KolmogorovTest","histograms with different binning along Z");
1551 return 0;
1552 }
1553
1554 // Should we include Uflows, Oflows?
1555 Int_t ibeg = 1, jbeg = 1, kbeg = 1;
1556 Int_t iend = ncx1, jend = ncy1, kend = ncz1;
1557 if (opt.Contains("U")) {ibeg = 0; jbeg = 0; kbeg = 0;}
1558 if (opt.Contains("O")) {iend = ncx1+1; jend = ncy1+1; kend = ncz1+1;}
1559
1560 Int_t i,j,k,bin;
1561 Double_t sum1 = 0;
1562 Double_t sum2 = 0;
1563 Double_t w1 = 0;
1564 Double_t w2 = 0;
1565 for (i = ibeg; i <= iend; i++) {
1566 for (j = jbeg; j <= jend; j++) {
1567 for (k = kbeg; k <= kend; k++) {
1568 bin = h1->GetBin(i,j,k);
1569 sum1 += h1->GetBinContent(bin);
1570 sum2 += h2->GetBinContent(bin);
1571 Double_t ew1 = h1->GetBinError(bin);
1572 Double_t ew2 = h2->GetBinError(bin);
1573 w1 += ew1*ew1;
1574 w2 += ew2*ew2;
1575 }
1576 }
1577 }
1578
1579
1580 // Check that both scatterplots contain events
1581 if (sum1 == 0) {
1582 Error("KolmogorovTest","Integral is zero for h1=%s\n",h1->GetName());
1583 return 0;
1584 }
1585 if (sum2 == 0) {
1586 Error("KolmogorovTest","Integral is zero for h2=%s\n",h2->GetName());
1587 return 0;
1588 }
1589 // calculate the effective entries.
1590 // the case when errors are zero (w1 == 0 or w2 ==0) are equivalent to
1591 // compare to a function. In that case the rescaling is done only on sqrt(esum2) or sqrt(esum1)
1592 Double_t esum1 = 0, esum2 = 0;
1593 if (w1 > 0)
1594 esum1 = sum1 * sum1 / w1;
1595 else
1596 afunc1 = kTRUE; // use later for calculating z
1597
1598 if (w2 > 0)
1599 esum2 = sum2 * sum2 / w2;
1600 else
1601 afunc2 = kTRUE; // use later for calculating z
1602
1603 if (afunc2 && afunc1) {
1604 Error("KolmogorovTest","Errors are zero for both histograms\n");
1605 return 0;
1606 }
1607
1608 // Find Kolmogorov distance
1609 // order is arbitrary take average of all possible 6 starting orders x,y,z
1610 int order[3] = {0,1,2};
1611 int binbeg[3];
1612 int binend[3];
1613 int ibin[3];
1614 binbeg[0] = ibeg; binbeg[1] = jbeg; binbeg[2] = kbeg;
1615 binend[0] = iend; binend[1] = jend; binend[2] = kend;
1616 Double_t vdfmax[6]; // there are in total 6 combinations
1617 int icomb = 0;
1618 Double_t s1 = 1./(6.*sum1);
1619 Double_t s2 = 1./(6.*sum2);
1620 Double_t rsum1=0, rsum2=0;
1621 do {
1622 // loop on bins
1623 Double_t dmax = 0;
1624 for (i = binbeg[order[0] ]; i <= binend[order[0] ]; i++) {
1625 for ( j = binbeg[order[1] ]; j <= binend[order[1] ]; j++) {
1626 for ( k = binbeg[order[2] ]; k <= binend[order[2] ]; k++) {
1627 ibin[ order[0] ] = i;
1628 ibin[ order[1] ] = j;
1629 ibin[ order[2] ] = k;
1630 bin = h1->GetBin(ibin[0],ibin[1],ibin[2]);
1631 rsum1 += s1*h1->GetBinContent(bin);
1632 rsum2 += s2*h2->GetBinContent(bin);
1633 dmax = TMath::Max(dmax, TMath::Abs(rsum1-rsum2));
1634 }
1635 }
1636 }
1637 vdfmax[icomb] = dmax;
1638 icomb++;
1639 } while (TMath::Permute(3,order) );
1640
1641
1642 // get average of distances
1643 Double_t dfmax = TMath::Mean(6,vdfmax);
1644
1645 // Get Kolmogorov probability
1646 Double_t factnm;
1647 if (afunc1) factnm = TMath::Sqrt(sum2);
1648 else if (afunc2) factnm = TMath::Sqrt(sum1);
1649 else factnm = TMath::Sqrt(sum1*sum2/(sum1+sum2));
1650 Double_t z = dfmax*factnm;
1651
1652 prb = TMath::KolmogorovProb(z);
1653
1654 Double_t prb1 = 0, prb2 = 0;
1655 // option N to combine normalization makes sense if both afunc1 and afunc2 are false
1656 if (opt.Contains("N") && !(afunc1 || afunc2 ) ) {
1657 // Combine probabilities for shape and normalization
1658 prb1 = prb;
1659 Double_t d12 = esum1-esum2;
1660 Double_t chi2 = d12*d12/(esum1+esum2);
1661 prb2 = TMath::Prob(chi2,1);
1662 // see Eadie et al., section 11.6.2
1663 if (prb > 0 && prb2 > 0) prb = prb*prb2*(1-TMath::Log(prb*prb2));
1664 else prb = 0;
1665 }
1666
1667 // debug printout
1668 if (opt.Contains("D")) {
1669 printf(" Kolmo Prob h1 = %s, sum1=%g\n",h1->GetName(),sum1);
1670 printf(" Kolmo Prob h2 = %s, sum2=%g\n",h2->GetName(),sum2);
1671 printf(" Kolmo Probabil = %f, Max Dist = %g\n",prb,dfmax);
1672 if (opt.Contains("N"))
1673 printf(" Kolmo Probabil = %f for shape alone, =%f for normalisation alone\n",prb1,prb2);
1674 }
1675 // This numerical error condition should never occur:
1676 if (TMath::Abs(rsum1-1) > 0.002) Warning("KolmogorovTest","Numerical problems with h1=%s\n",h1->GetName());
1677 if (TMath::Abs(rsum2-1) > 0.002) Warning("KolmogorovTest","Numerical problems with h2=%s\n",h2->GetName());
1678
1679 if (opt.Contains("M")) return dfmax; // return average of max distance
1680
1681 return prb;
1682}
1683
1684
1685////////////////////////////////////////////////////////////////////////////////
1686/// Project a 3-D histogram into a 1-D histogram along X.
1687///
1688/// The projection is always of the type TH1D.
1689/// The projection is made from the cells along the X axis
1690/// ranging from iymin to iymax and izmin to izmax included.
1691/// By default, underflow and overflows are included in both the Y and Z axis.
1692/// By Setting iymin=1 and iymax=NbinsY the underflow and/or overflow in Y will be excluded
1693/// By setting izmin=1 and izmax=NbinsZ the underflow and/or overflow in Z will be excluded
1694///
1695/// if option "e" is specified, the errors are computed.
1696/// if option "d" is specified, the projection is drawn in the current pad.
1697/// if option "o" original axis range of the target axes will be
1698/// kept, but only bins inside the selected range will be filled.
1699///
1700/// NOTE that if a TH1D named "name" exists in the current directory or pad
1701/// the histogram is reset and filled again with the projected contents of the TH3.
1702///
1703/// implemented using Project3D
1704
1705TH1D *TH3::ProjectionX(const char *name, Int_t iymin, Int_t iymax,
1706 Int_t izmin, Int_t izmax, Option_t *option) const
1707{
1708 // in case of default name append the parent name
1709 TString hname = name;
1710 if (hname == "_px") hname = TString::Format("%s%s", GetName(), name);
1711 TString title = TString::Format("%s ( Projection X )",GetTitle());
1712
1713 // when projecting in Z outer axis are Y and Z (order is important. It is defined in the DoProject1D function)
1714 return DoProject1D(hname, title, iymin, iymax, izmin, izmax, &fXaxis, &fYaxis, &fZaxis, option);
1715}
1716
1717
1718////////////////////////////////////////////////////////////////////////////////
1719/// Project a 3-D histogram into a 1-D histogram along Y.
1720///
1721/// The projection is always of the type TH1D.
1722/// The projection is made from the cells along the Y axis
1723/// ranging from ixmin to ixmax and izmin to izmax included.
1724/// By default, underflow and overflow are included in both the X and Z axis.
1725/// By setting ixmin=1 and ixmax=NbinsX the underflow and/or overflow in X will be excluded
1726/// By setting izmin=1 and izmax=NbinsZ the underflow and/or overflow in Z will be excluded
1727///
1728/// if option "e" is specified, the errors are computed.
1729/// if option "d" is specified, the projection is drawn in the current pad.
1730/// if option "o" original axis range of the target axes will be
1731/// kept, but only bins inside the selected range will be filled.
1732///
1733/// NOTE that if a TH1D named "name" exists in the current directory or pad,
1734/// the histogram is reset and filled again with the projected contents of the TH3.
1735///
1736/// implemented using Project3D
1737
1738TH1D *TH3::ProjectionY(const char *name, Int_t ixmin, Int_t ixmax,
1739 Int_t izmin, Int_t izmax, Option_t *option) const
1740{
1741 TString hname = name;
1742 if (hname == "_py") hname = TString::Format("%s%s", GetName(), name);
1743 TString title = TString::Format("%s ( Projection Y )",GetTitle());
1744
1745 // when projecting in Z outer axis are X and Y (order is important. It is defined in the DoProject1D function)
1746 return DoProject1D(hname, title, ixmin, ixmax, izmin, izmax, &fYaxis, &fXaxis, &fZaxis, option);
1747}
1748
1749////////////////////////////////////////////////////////////////////////////////
1750/// Project a 3-D histogram into a 1-D histogram along Z.
1751///
1752/// The projection is always of the type TH1D.
1753/// The projection is made from the cells along the Z axis
1754/// ranging from ixmin to ixmax and iymin to iymax included.
1755/// By default, bins 1 to nx and 1 to ny are included
1756/// By default, underflow and overflow are included in both the X and Y axis.
1757/// By Setting ixmin=1 and ixmax=NbinsX the underflow and/or overflow in X will be excluded
1758/// By setting iymin=1 and/or iymax=NbinsY the underflow and/or overflow in Y will be excluded
1759///
1760/// if option "e" is specified, the errors are computed.
1761/// if option "d" is specified, the projection is drawn in the current pad.
1762/// if option "o" original axis range of the target axes will be
1763/// kept, but only bins inside the selected range will be filled.
1764///
1765/// NOTE that if a TH1D named "name" exists in the current directory or pad,
1766/// the histogram is reset and filled again with the projected contents of the TH3.
1767///
1768/// implemented using Project3D
1769
1770TH1D *TH3::ProjectionZ(const char *name, Int_t ixmin, Int_t ixmax,
1771 Int_t iymin, Int_t iymax, Option_t *option) const
1772{
1773
1774 TString hname = name;
1775 if (hname == "_pz") hname = TString::Format("%s%s", GetName(), name);
1776 TString title = TString::Format("%s ( Projection Z )",GetTitle());
1777
1778 // when projecting in Z outer axis are X and Y (order is important. It is defined in the DoProject1D function)
1779 return DoProject1D(hname, title, ixmin, ixmax, iymin, iymax, &fZaxis, &fXaxis, &fYaxis, option);
1780}
1781
1782
1783////////////////////////////////////////////////////////////////////////////////
1784/// internal method performing the projection to 1D histogram
1785/// called from TH3::Project3D
1786
1787TH1D *TH3::DoProject1D(const char* name, const char * title, int imin1, int imax1, int imin2, int imax2,
1788 const TAxis* projAxis, const TAxis * axis1, const TAxis * axis2, Option_t * option) const
1789{
1790
1791 TString opt = option;
1792 opt.ToLower();
1793
1794 // save previous axis range and bits
1795 // Int_t iminOld1 = axis1->GetFirst();
1796 // Int_t imaxOld1 = axis1->GetLast();
1797 // Int_t iminOld2 = axis2->GetFirst();
1798 // Int_t imaxOld2 = axis2->GetLast();
1799 // Bool_t hadRange1 = axis1->TestBit(TAxis::kAxisRange);
1800 // Bool_t hadRange2 = axis2->TestBit(TAxis::kAxisRange);
1801
1802 // need to cast-away constness to set range
1803 TAxis out1(*axis1);
1804 TAxis out2(*axis2);
1805 // const_cast<TAxis *>(axis1)->SetRange(imin1, imax1);
1806 // const_cast<TAxis*>(axis2)->SetRange(imin2,imax2);
1807 out1.SetRange(imin1, imax1);
1808 out2.SetRange(imin2, imax2);
1809
1810 Bool_t computeErrors = GetSumw2N();
1811 if (opt.Contains("e") ) {
1812 computeErrors = kTRUE;
1813 opt.Remove(opt.First("e"),1);
1814 }
1815 Bool_t originalRange = kFALSE;
1816 if (opt.Contains('o') ) {
1817 originalRange = kTRUE;
1818 opt.Remove(opt.First("o"),1);
1819 }
1820
1821 TH1D * h1 = DoProject1D(name, title, projAxis, &out1, &out2, computeErrors, originalRange,true,true);
1822
1823 // // restore original range
1824 // if (axis1->TestBit(TAxis::kAxisRange)) {
1825 // if (hadRange1) const_cast<TAxis*>(axis1)->SetRange(iminOld1,imaxOld1);
1826 // if (axis2->TestBit(TAxis::kAxisRange)) const_cast<TAxis*>(axis2)->SetRange(iminOld2,imaxOld2);
1827 // // we need also to restore the original bits
1828
1829 // draw in current pad
1830 if (h1 && opt.Contains("d")) {
1831 opt.Remove(opt.First("d"),1);
1832 TVirtualPad *padsav = gPad;
1833 TVirtualPad *pad = gROOT->GetSelectedPad();
1834 if (pad) pad->cd();
1835 if (!gPad || !gPad->FindObject(h1)) {
1836 h1->Draw(opt);
1837 } else {
1838 h1->Paint(opt);
1839 }
1840 if (padsav) padsav->cd();
1841 }
1842
1843 return h1;
1844}
1845
1846////////////////////////////////////////////////////////////////////////////////
1847/// internal methdod performing the projection to 1D histogram
1848/// called from other TH3::DoProject1D
1849
1850TH1D *TH3::DoProject1D(const char* name, const char * title, const TAxis* projX,
1851 const TAxis * out1, const TAxis * out2,
1852 bool computeErrors, bool originalRange,
1853 bool useUF, bool useOF) const
1854{
1855 // Create the projection histogram
1856 TH1D *h1 = 0;
1857
1858 // Get range to use as well as bin limits
1859 // Projected range must be inside and not outside original one (ROOT-8781)
1860 Int_t ixmin = std::max(projX->GetFirst(),1);
1861 Int_t ixmax = std::min(projX->GetLast(),projX->GetNbins());
1862 Int_t nx = ixmax-ixmin+1;
1863
1864 // Create the histogram, either reseting a preexisting one
1865 TObject *h1obj = gROOT->FindObject(name);
1866 if (h1obj && h1obj->InheritsFrom(TH1::Class())) {
1867 if (h1obj->IsA() != TH1D::Class() ) {
1868 Error("DoProject1D","Histogram with name %s must be a TH1D and is a %s",name,h1obj->ClassName());
1869 return 0;
1870 }
1871 h1 = (TH1D*)h1obj;
1872 // reset histogram and re-set the axis in any case
1873 h1->Reset();
1874 const TArrayD *bins = projX->GetXbins();
1875 if ( originalRange )
1876 {
1877 if (bins->fN == 0) {
1878 h1->SetBins(projX->GetNbins(),projX->GetXmin(),projX->GetXmax());
1879 } else {
1880 h1->SetBins(projX->GetNbins(),bins->fArray);
1881 }
1882 } else {
1883 if (bins->fN == 0) {
1884 h1->SetBins(nx,projX->GetBinLowEdge(ixmin),projX->GetBinUpEdge(ixmax));
1885 } else {
1886 h1->SetBins(nx,&bins->fArray[ixmin-1]);
1887 }
1888 }
1889 }
1890
1891 if (!h1) {
1892 const TArrayD *bins = projX->GetXbins();
1893 if ( originalRange )
1894 {
1895 if (bins->fN == 0) {
1896 h1 = new TH1D(name,title,projX->GetNbins(),projX->GetXmin(),projX->GetXmax());
1897 } else {
1898 h1 = new TH1D(name,title,projX->GetNbins(),bins->fArray);
1899 }
1900 } else {
1901 if (bins->fN == 0) {
1902 h1 = new TH1D(name,title,nx,projX->GetBinLowEdge(ixmin),projX->GetBinUpEdge(ixmax));
1903 } else {
1904 h1 = new TH1D(name,title,nx,&bins->fArray[ixmin-1]);
1905 }
1906 }
1907 }
1908
1909 // Copy the axis attributes and the axis labels if needed.
1910 h1->GetXaxis()->ImportAttributes(projX);
1911 THashList* labels = projX->GetLabels();
1912 if (labels) {
1913 TIter iL(labels);
1914 TObjString* lb;
1915 Int_t i = 1;
1916 while ((lb=(TObjString*)iL())) {
1917 h1->GetXaxis()->SetBinLabel(i,lb->String().Data());
1918 i++;
1919 }
1920 }
1921 h1->SetLineColor(this->GetLineColor());
1922 h1->SetFillColor(this->GetFillColor());
1923 h1->SetMarkerColor(this->GetMarkerColor());
1924 h1->SetMarkerStyle(this->GetMarkerStyle());
1925
1926 // Activate errors
1927 if ( computeErrors && (h1->GetSumw2N() != h1->GetNcells() ) ) h1->Sumw2();
1928
1929 // Set references to the axies in case out1 or out2 ar enot provided
1930 // and one can use the histogram axis given projX
1931 if (out1 == nullptr && out2 == nullptr) {
1932 if (projX == GetXaxis()) {
1933 out1 = GetYaxis();
1934 out2 = GetZaxis();
1935 } else if (projX == GetYaxis()) {
1936 out1 = GetXaxis();
1937 out2 = GetZaxis();
1938 } else {
1939 out1 = GetXaxis();
1940 out2 = GetYaxis();
1941 }
1942 }
1943 R__ASSERT(out1 != nullptr && out2 != nullptr);
1944
1945 Int_t *refX = 0, *refY = 0, *refZ = 0;
1946 Int_t ixbin, out1bin, out2bin;
1947 if (projX == GetXaxis()) {
1948 refX = &ixbin;
1949 refY = &out1bin;
1950 refZ = &out2bin;
1951 }
1952 if (projX == GetYaxis()) {
1953 refX = &out1bin;
1954 refY = &ixbin;
1955 refZ = &out2bin;
1956 }
1957 if (projX == GetZaxis()) {
1958 refX = &out1bin;
1959 refY = &out2bin;
1960 refZ = &ixbin;
1961 }
1962 R__ASSERT (refX != 0 && refY != 0 && refZ != 0);
1963
1964 // Fill the projected histogram excluding underflow/overflows if considered in the option
1965 // if specified in the option (by default they considered)
1966 Double_t totcont = 0;
1967
1968 Int_t out1min = out1->GetFirst();
1969 Int_t out1max = out1->GetLast();
1970 // GetFirst(), GetLast() can return (0,0) when the range bit is set artificially (see TAxis::SetRange)
1971 //if (out1min == 0 && out1max == 0) { out1min = 1; out1max = out1->GetNbins(); }
1972 // correct for underflow/overflows
1973 if (useUF && !out1->TestBit(TAxis::kAxisRange) ) out1min -= 1;
1974 if (useOF && !out1->TestBit(TAxis::kAxisRange) ) out1max += 1;
1975 Int_t out2min = out2->GetFirst();
1976 Int_t out2max = out2->GetLast();
1977// if (out2min == 0 && out2max == 0) { out2min = 1; out2max = out2->GetNbins(); }
1978 if (useUF && !out2->TestBit(TAxis::kAxisRange) ) out2min -= 1;
1979 if (useOF && !out2->TestBit(TAxis::kAxisRange) ) out2max += 1;
1980
1981 for (ixbin=0;ixbin<=1+projX->GetNbins();ixbin++) {
1982 if ( projX->TestBit(TAxis::kAxisRange) && ( ixbin < ixmin || ixbin > ixmax )) continue;
1983
1984 Double_t cont = 0;
1985 Double_t err2 = 0;
1986
1987 // loop on the bins to be integrated (outbin should be called inbin)
1988 for (out1bin = out1min; out1bin <= out1max; out1bin++) {
1989 for (out2bin = out2min; out2bin <= out2max; out2bin++) {
1990
1991 Int_t bin = GetBin(*refX, *refY, *refZ);
1992
1993 // sum the bin contents and errors if needed
1994 cont += RetrieveBinContent(bin);
1995 if (computeErrors) {
1996 Double_t exyz = GetBinError(bin);
1997 err2 += exyz*exyz;
1998 }
1999 }
2000 }
2001 Int_t ix = h1->FindBin( projX->GetBinCenter(ixbin) );
2002 h1->SetBinContent(ix ,cont);
2003 if (computeErrors) h1->SetBinError(ix, TMath::Sqrt(err2) );
2004 // sum all content
2005 totcont += cont;
2006
2007 }
2008
2009 // since we use a combination of fill and SetBinError we need to reset and recalculate the statistics
2010 // for weighted histograms otherwise sumw2 will be wrong.
2011 // We can keep the original statistics from the TH3 if the projected sumw is consistent with original one
2012 // i.e. when no events are thrown away
2013 bool resetStats = true;
2014 double eps = 1.E-12;
2015 if (IsA() == TH3F::Class() ) eps = 1.E-6;
2016 if (fTsumw != 0 && TMath::Abs( fTsumw - totcont) < TMath::Abs(fTsumw) * eps) resetStats = false;
2017
2018 bool resetEntries = resetStats;
2019 // entries are calculated using underflow/overflow. If excluded entries must be reset
2020 resetEntries |= !useUF || !useOF;
2021
2022
2023 if (!resetStats) {
2024 Double_t stats[kNstat];
2025 GetStats(stats);
2026 if ( projX == GetYaxis() ) {
2027 stats[2] = stats[4];
2028 stats[3] = stats[5];
2029 }
2030 else if ( projX == GetZaxis() ) {
2031 stats[2] = stats[7];
2032 stats[3] = stats[8];
2033 }
2034 h1->PutStats(stats);
2035 }
2036 else {
2037 // reset statistics
2038 h1->ResetStats();
2039 }
2040 if (resetEntries) {
2041 // in case of error calculation (i.e. when Sumw2() is set)
2042 // use the effective entries for the entries
2043 // since this is the only way to estimate them
2044 Double_t entries = TMath::Floor( totcont + 0.5); // to avoid numerical rounding
2045 if (computeErrors) entries = h1->GetEffectiveEntries();
2046 h1->SetEntries( entries );
2047 }
2048 else {
2050 }
2051
2052 return h1;
2053}
2054
2055
2056////////////////////////////////////////////////////////////////////////////////
2057/// internal method performing the projection to a 2D histogram
2058/// called from TH3::Project3D
2059
2060TH2D *TH3::DoProject2D(const char* name, const char * title, const TAxis* projX, const TAxis* projY,
2061 bool computeErrors, bool originalRange,
2062 bool useUF, bool useOF) const
2063{
2064 TH2D *h2 = 0;
2065
2066 // Get range to use as well as bin limits
2067 Int_t ixmin = std::max(projX->GetFirst(),1);
2068 Int_t ixmax = std::min(projX->GetLast(),projX->GetNbins());
2069 Int_t iymin = std::max(projY->GetFirst(),1);
2070 Int_t iymax = std::min(projY->GetLast(),projY->GetNbins());
2071
2072 Int_t nx = ixmax-ixmin+1;
2073 Int_t ny = iymax-iymin+1;
2074
2075 // Create the histogram, either reseting a preexisting one
2076 // or creating one from scratch.
2077 // Does an object with the same name exists?
2078 TObject *h2obj = gROOT->FindObject(name);
2079 if (h2obj && h2obj->InheritsFrom(TH1::Class())) {
2080 if ( h2obj->IsA() != TH2D::Class() ) {
2081 Error("DoProject2D","Histogram with name %s must be a TH2D and is a %s",name,h2obj->ClassName());
2082 return 0;
2083 }
2084 h2 = (TH2D*)h2obj;
2085 // reset histogram and its axes
2086 h2->Reset();
2087 const TArrayD *xbins = projX->GetXbins();
2088 const TArrayD *ybins = projY->GetXbins();
2089 if ( originalRange ) {
2090 h2->SetBins(projY->GetNbins(),projY->GetXmin(),projY->GetXmax()
2091 ,projX->GetNbins(),projX->GetXmin(),projX->GetXmax());
2092 // set bins for mixed axis do not exists - need to set afterwards the variable bins
2093 if (ybins->fN != 0)
2094 h2->GetXaxis()->Set(projY->GetNbins(),&ybins->fArray[iymin-1]);
2095 if (xbins->fN != 0)
2096 h2->GetYaxis()->Set(projX->GetNbins(),&xbins->fArray[ixmin-1]);
2097 } else {
2098 h2->SetBins(ny,projY->GetBinLowEdge(iymin),projY->GetBinUpEdge(iymax)
2099 ,nx,projX->GetBinLowEdge(ixmin),projX->GetBinUpEdge(ixmax));
2100 if (ybins->fN != 0)
2101 h2->GetXaxis()->Set(ny,&ybins->fArray[iymin-1]);
2102 if (xbins->fN != 0)
2103 h2->GetYaxis()->Set(nx,&xbins->fArray[ixmin-1]);
2104 }
2105 }
2106
2107
2108 if (!h2) {
2109 const TArrayD *xbins = projX->GetXbins();
2110 const TArrayD *ybins = projY->GetXbins();
2111 if ( originalRange )
2112 {
2113 if (xbins->fN == 0 && ybins->fN == 0) {
2114 h2 = new TH2D(name,title,projY->GetNbins(),projY->GetXmin(),projY->GetXmax()
2115 ,projX->GetNbins(),projX->GetXmin(),projX->GetXmax());
2116 } else if (ybins->fN == 0) {
2117 h2 = new TH2D(name,title,projY->GetNbins(),projY->GetXmin(),projY->GetXmax()
2118 ,projX->GetNbins(),&xbins->fArray[ixmin-1]);
2119 } else if (xbins->fN == 0) {
2120 h2 = new TH2D(name,title,projY->GetNbins(),&ybins->fArray[iymin-1]
2121 ,projX->GetNbins(),projX->GetXmin(),projX->GetXmax());
2122 } else {
2123 h2 = new TH2D(name,title,projY->GetNbins(),&ybins->fArray[iymin-1],projX->GetNbins(),&xbins->fArray[ixmin-1]);
2124 }
2125 } else {
2126 if (xbins->fN == 0 && ybins->fN == 0) {
2127 h2 = new TH2D(name,title,ny,projY->GetBinLowEdge(iymin),projY->GetBinUpEdge(iymax)
2128 ,nx,projX->GetBinLowEdge(ixmin),projX->GetBinUpEdge(ixmax));
2129 } else if (ybins->fN == 0) {
2130 h2 = new TH2D(name,title,ny,projY->GetBinLowEdge(iymin),projY->GetBinUpEdge(iymax)
2131 ,nx,&xbins->fArray[ixmin-1]);
2132 } else if (xbins->fN == 0) {
2133 h2 = new TH2D(name,title,ny,&ybins->fArray[iymin-1]
2134 ,nx,projX->GetBinLowEdge(ixmin),projX->GetBinUpEdge(ixmax));
2135 } else {
2136 h2 = new TH2D(name,title,ny,&ybins->fArray[iymin-1],nx,&xbins->fArray[ixmin-1]);
2137 }
2138 }
2139 }
2140
2141 // Copy the axis attributes and the axis labels if needed.
2142 THashList* labels1 = 0;
2143 THashList* labels2 = 0;
2144 // "xy"
2145 h2->GetXaxis()->ImportAttributes(projY);
2146 h2->GetYaxis()->ImportAttributes(projX);
2147 labels1 = projY->GetLabels();
2148 labels2 = projX->GetLabels();
2149 if (labels1) {
2150 TIter iL(labels1);
2151 TObjString* lb;
2152 Int_t i = 1;
2153 while ((lb=(TObjString*)iL())) {
2154 h2->GetXaxis()->SetBinLabel(i,lb->String().Data());
2155 i++;
2156 }
2157 }
2158 if (labels2) {
2159 TIter iL(labels2);
2160 TObjString* lb;
2161 Int_t i = 1;
2162 while ((lb=(TObjString*)iL())) {
2163 h2->GetYaxis()->SetBinLabel(i,lb->String().Data());
2164 i++;
2165 }
2166 }
2167 h2->SetLineColor(this->GetLineColor());
2168 h2->SetFillColor(this->GetFillColor());
2169 h2->SetMarkerColor(this->GetMarkerColor());
2170 h2->SetMarkerStyle(this->GetMarkerStyle());
2171
2172 // Activate errors
2173 if ( computeErrors && (h2->GetSumw2N() != h2->GetNcells()) ) h2->Sumw2();
2174
2175 // Set references to the axis, so that the bucle has no branches.
2176 const TAxis* out = 0;
2177 if ( projX != GetXaxis() && projY != GetXaxis() ) {
2178 out = GetXaxis();
2179 } else if ( projX != GetYaxis() && projY != GetYaxis() ) {
2180 out = GetYaxis();
2181 } else {
2182 out = GetZaxis();
2183 }
2184
2185 Int_t *refX = 0, *refY = 0, *refZ = 0;
2186 Int_t ixbin, iybin, outbin;
2187 if ( projX == GetXaxis() && projY == GetYaxis() ) { refX = &ixbin; refY = &iybin; refZ = &outbin; }
2188 if ( projX == GetYaxis() && projY == GetXaxis() ) { refX = &iybin; refY = &ixbin; refZ = &outbin; }
2189 if ( projX == GetXaxis() && projY == GetZaxis() ) { refX = &ixbin; refY = &outbin; refZ = &iybin; }
2190 if ( projX == GetZaxis() && projY == GetXaxis() ) { refX = &iybin; refY = &outbin; refZ = &ixbin; }
2191 if ( projX == GetYaxis() && projY == GetZaxis() ) { refX = &outbin; refY = &ixbin; refZ = &iybin; }
2192 if ( projX == GetZaxis() && projY == GetYaxis() ) { refX = &outbin; refY = &iybin; refZ = &ixbin; }
2193 R__ASSERT (refX != 0 && refY != 0 && refZ != 0);
2194
2195 // Fill the projected histogram excluding underflow/overflows if considered in the option
2196 // if specified in the option (by default they considered)
2197 Double_t totcont = 0;
2198
2199 Int_t outmin = out->GetFirst();
2200 Int_t outmax = out->GetLast();
2201 // GetFirst(), GetLast() can return (0,0) when the range bit is set artificially (see TAxis::SetRange)
2202 if (outmin == 0 && outmax == 0) { outmin = 1; outmax = out->GetNbins(); }
2203 // correct for underflow/overflows
2204 if (useUF && !out->TestBit(TAxis::kAxisRange) ) outmin -= 1;
2205 if (useOF && !out->TestBit(TAxis::kAxisRange) ) outmax += 1;
2206
2207 for (ixbin=0;ixbin<=1+projX->GetNbins();ixbin++) {
2208 if ( projX->TestBit(TAxis::kAxisRange) && ( ixbin < ixmin || ixbin > ixmax )) continue;
2209 Int_t ix = h2->GetYaxis()->FindBin( projX->GetBinCenter(ixbin) );
2210
2211 for (iybin=0;iybin<=1+projY->GetNbins();iybin++) {
2212 if ( projY->TestBit(TAxis::kAxisRange) && ( iybin < iymin || iybin > iymax )) continue;
2213 Int_t iy = h2->GetXaxis()->FindBin( projY->GetBinCenter(iybin) );
2214
2215 Double_t cont = 0;
2216 Double_t err2 = 0;
2217
2218 // loop on the bins to be integrated (outbin should be called inbin)
2219 for (outbin = outmin; outbin <= outmax; outbin++) {
2220
2221 Int_t bin = GetBin(*refX,*refY,*refZ);
2222
2223 // sum the bin contents and errors if needed
2224 cont += RetrieveBinContent(bin);
2225 if (computeErrors) {
2226 Double_t exyz = GetBinError(bin);
2227 err2 += exyz*exyz;
2228 }
2229
2230 }
2231
2232 // remember axis are inverted
2233 h2->SetBinContent(iy , ix, cont);
2234 if (computeErrors) h2->SetBinError(iy, ix, TMath::Sqrt(err2) );
2235 // sum all content
2236 totcont += cont;
2237
2238 }
2239 }
2240
2241 // since we use fill we need to reset and recalculate the statistics (see comment in DoProject1D )
2242 // or keep original statistics if consistent sumw2
2243 bool resetStats = true;
2244 double eps = 1.E-12;
2245 if (IsA() == TH3F::Class() ) eps = 1.E-6;
2246 if (fTsumw != 0 && TMath::Abs( fTsumw - totcont) < TMath::Abs(fTsumw) * eps) resetStats = false;
2247
2248 bool resetEntries = resetStats;
2249 // entries are calculated using underflow/overflow. If excluded entries must be reset
2250 resetEntries |= !useUF || !useOF;
2251
2252 if (!resetStats) {
2253 Double_t stats[kNstat];
2254 Double_t oldst[kNstat]; // old statistics
2255 for (Int_t i = 0; i < kNstat; ++i) { oldst[i] = 0; }
2256 GetStats(oldst);
2257 std::copy(oldst,oldst+kNstat,stats);
2258 // not that projX refer to Y axis and projX refer to the X axis of projected histogram
2259 // nothing to do for projection in Y vs X
2260 if ( projY == GetXaxis() && projX == GetZaxis() ) { // case XZ
2261 stats[4] = oldst[7];
2262 stats[5] = oldst[8];
2263 stats[6] = oldst[9];
2264 }
2265 if ( projY == GetYaxis() ) {
2266 stats[2] = oldst[4];
2267 stats[3] = oldst[5];
2268 if ( projX == GetXaxis() ) { // case YX
2269 stats[4] = oldst[2];
2270 stats[5] = oldst[3];
2271 }
2272 if ( projX == GetZaxis() ) { // case YZ
2273 stats[4] = oldst[7];
2274 stats[5] = oldst[8];
2275 stats[6] = oldst[10];
2276 }
2277 }
2278 else if ( projY == GetZaxis() ) {
2279 stats[2] = oldst[7];
2280 stats[3] = oldst[8];
2281 if ( projX == GetXaxis() ) { // case ZX
2282 stats[4] = oldst[2];
2283 stats[5] = oldst[3];
2284 stats[6] = oldst[9];
2285 }
2286 if ( projX == GetYaxis() ) { // case ZY
2287 stats[4] = oldst[4];
2288 stats[5] = oldst[5];
2289 stats[6] = oldst[10];
2290 }
2291 }
2292 // set the new statistics
2293 h2->PutStats(stats);
2294 }
2295 else {
2296 // recalculate the statistics
2297 h2->ResetStats();
2298 }
2299
2300 if (resetEntries) {
2301 // use the effective entries for the entries
2302 // since this is the only way to estimate them
2303 Double_t entries = h2->GetEffectiveEntries();
2304 if (!computeErrors) entries = TMath::Floor( entries + 0.5); // to avoid numerical rounding
2305 h2->SetEntries( entries );
2306 }
2307 else {
2308 h2->SetEntries( fEntries );
2309 }
2310
2311
2312 return h2;
2313}
2314
2315
2316////////////////////////////////////////////////////////////////////////////////
2317/// Project a 3-d histogram into 1 or 2-d histograms depending on the
2318/// option parameter, which may contain a combination of the characters x,y,z,e
2319/// - option = "x" return the x projection into a TH1D histogram
2320/// - option = "y" return the y projection into a TH1D histogram
2321/// - option = "z" return the z projection into a TH1D histogram
2322/// - option = "xy" return the x versus y projection into a TH2D histogram
2323/// - option = "yx" return the y versus x projection into a TH2D histogram
2324/// - option = "xz" return the x versus z projection into a TH2D histogram
2325/// - option = "zx" return the z versus x projection into a TH2D histogram
2326/// - option = "yz" return the y versus z projection into a TH2D histogram
2327/// - option = "zy" return the z versus y projection into a TH2D histogram
2328///
2329/// NB: the notation "a vs b" means "a" vertical and "b" horizontal
2330///
2331/// option = "o" original axis range of the target axes will be
2332/// kept, but only bins inside the selected range will be filled.
2333///
2334/// If option contains the string "e", errors are computed
2335///
2336/// The projection is made for the selected bins only.
2337/// To select a bin range along an axis, use TAxis::SetRange, eg
2338/// h3.GetYaxis()->SetRange(23,56);
2339///
2340/// NOTE 1: The generated histogram is named th3name + option
2341/// eg if the TH3* h histogram is named "myhist", then
2342/// h->Project3D("xy"); produces a TH2D histogram named "myhist_xy"
2343/// if a histogram of the same type already exists, it is overwritten.
2344/// The following sequence
2345/// h->Project3D("xy");
2346/// h->Project3D("xy2");
2347/// will generate two TH2D histograms named "myhist_xy" and "myhist_xy2"
2348/// A different name can be generated by attaching a string to the option
2349/// For example h->Project3D("name_xy") will generate an histogram with the name: h3dname_name_xy.
2350///
2351/// NOTE 2: If an histogram of the same type already exists,
2352/// the histogram is reset and filled again with the projected contents of the TH3.
2353///
2354/// NOTE 3: The number of entries in the projected histogram is estimated from the number of
2355/// effective entries for all the cells included in the projection.
2356///
2357/// NOTE 4: underflow/overflow are included by default in the projection
2358/// To exclude underflow and/or overflow (for both axis in case of a projection to a 1D histogram) use option "NUF" and/or "NOF"
2359/// With SetRange() you can have all bins except underflow/overflow only if you set the axis bit range as
2360/// following after having called SetRange: axis->SetRange(1, axis->GetNbins());
2361
2363{
2364 TString opt = option; opt.ToLower();
2365 Int_t pcase = 0;
2366 TString ptype;
2367 if (opt.Contains("x")) { pcase = 1; ptype = "x"; }
2368 if (opt.Contains("y")) { pcase = 2; ptype = "y"; }
2369 if (opt.Contains("z")) { pcase = 3; ptype = "z"; }
2370 if (opt.Contains("xy")) { pcase = 4; ptype = "xy"; }
2371 if (opt.Contains("yx")) { pcase = 5; ptype = "yx"; }
2372 if (opt.Contains("xz")) { pcase = 6; ptype = "xz"; }
2373 if (opt.Contains("zx")) { pcase = 7; ptype = "zx"; }
2374 if (opt.Contains("yz")) { pcase = 8; ptype = "yz"; }
2375 if (opt.Contains("zy")) { pcase = 9; ptype = "zy"; }
2376
2377 if (pcase == 0) {
2378 Error("Project3D","No projection axis specified - return a NULL pointer");
2379 return 0;
2380 }
2381 // do not remove ptype from opt to use later in the projected histo name
2382
2383 Bool_t computeErrors = GetSumw2N();
2384 if (opt.Contains("e") ) {
2385 computeErrors = kTRUE;
2386 opt.Remove(opt.First("e"),1);
2387 }
2388
2389 Bool_t useUF = kTRUE;
2390 Bool_t useOF = kTRUE;
2391 if (opt.Contains("nuf") ) {
2392 useUF = kFALSE;
2393 opt.Remove(opt.Index("nuf"),3);
2394 }
2395 if (opt.Contains("nof") ) {
2396 useOF = kFALSE;
2397 opt.Remove(opt.Index("nof"),3);
2398 }
2399
2400 Bool_t originalRange = kFALSE;
2401 if (opt.Contains('o') ) {
2402 originalRange = kTRUE;
2403 opt.Remove(opt.First("o"),1);
2404 }
2405
2406
2407 // Create the projection histogram
2408 TH1 *h = 0;
2409
2410 TString name = GetName();
2411 TString title = GetTitle();
2412 name += "_"; name += opt; // opt may include a user defined name
2413 title += " "; title += ptype; title += " projection";
2414
2415 switch (pcase) {
2416 case 1:
2417 // "x"
2418 h = DoProject1D(name, title, this->GetXaxis(), nullptr, nullptr,
2419 computeErrors, originalRange, useUF, useOF);
2420 break;
2421
2422 case 2:
2423 // "y"
2424 h = DoProject1D(name, title, this->GetYaxis(), nullptr, nullptr,
2425 computeErrors, originalRange, useUF, useOF);
2426 break;
2427
2428 case 3:
2429 // "z"
2430 h = DoProject1D(name, title, this->GetZaxis(), nullptr, nullptr,
2431 computeErrors, originalRange, useUF, useOF);
2432 break;
2433
2434 case 4:
2435 // "xy"
2436 h = DoProject2D(name, title, this->GetXaxis(),this->GetYaxis(),
2437 computeErrors, originalRange, useUF, useOF);
2438 break;
2439
2440 case 5:
2441 // "yx"
2442 h = DoProject2D(name, title, this->GetYaxis(),this->GetXaxis(),
2443 computeErrors, originalRange, useUF, useOF);
2444 break;
2445
2446 case 6:
2447 // "xz"
2448 h = DoProject2D(name, title, this->GetXaxis(),this->GetZaxis(),
2449 computeErrors, originalRange, useUF, useOF);
2450 break;
2451
2452 case 7:
2453 // "zx"
2454 h = DoProject2D(name, title, this->GetZaxis(),this->GetXaxis(),
2455 computeErrors, originalRange, useUF, useOF);
2456 break;
2457
2458 case 8:
2459 // "yz"
2460 h = DoProject2D(name, title, this->GetYaxis(),this->GetZaxis(),
2461 computeErrors, originalRange, useUF, useOF);
2462 break;
2463
2464 case 9:
2465 // "zy"
2466 h = DoProject2D(name, title, this->GetZaxis(),this->GetYaxis(),
2467 computeErrors, originalRange, useUF, useOF);
2468 break;
2469
2470 }
2471
2472 // draw in current pad
2473 if (h && opt.Contains("d")) {
2474 opt.Remove(opt.First("d"),1);
2475 TVirtualPad *padsav = gPad;
2476 TVirtualPad *pad = gROOT->GetSelectedPad();
2477 if (pad) pad->cd();
2478 if (!gPad || !gPad->FindObject(h)) {
2479 h->Draw(opt);
2480 } else {
2481 h->Paint(opt);
2482 }
2483 if (padsav) padsav->cd();
2484 }
2485
2486 return h;
2487}
2488
2489
2490////////////////////////////////////////////////////////////////////////////////
2491/// internal function to fill the bins of the projected profile 2D histogram
2492/// called from DoProjectProfile2D
2493
2495 const TAxis & a1, const TAxis & a2, const TAxis & a3,
2496 Int_t bin1, Int_t bin2, Int_t bin3,
2497 Int_t inBin, Bool_t useWeights ) const {
2498 Double_t cont = GetBinContent(inBin);
2499 if (!cont) return;
2500 TArrayD & binSumw2 = *(p2->GetBinSumw2());
2501 if (useWeights && binSumw2.fN <= 0) useWeights = false;
2502 if (!useWeights) p2->SetBit(TH1::kIsNotW); // to use Fill for setting the bin contents of the Profile
2503 // the following fill update wrongly the fBinSumw2- need to save it before
2504 Double_t u = a1.GetBinCenter(bin1);
2505 Double_t v = a2.GetBinCenter(bin2);
2506 Double_t w = a3.GetBinCenter(bin3);
2507 Int_t outBin = p2->FindBin(u, v);
2508 if (outBin <0) return;
2509 Double_t tmp = 0;
2510 if ( useWeights ) tmp = binSumw2.fArray[outBin];
2511 p2->Fill( u , v, w, cont);
2512 if (useWeights ) binSumw2.fArray[outBin] = tmp + fSumw2.fArray[inBin];
2513}
2514
2515
2516////////////////////////////////////////////////////////////////////////////////
2517/// internal method to project to a 2D Profile
2518/// called from TH3::Project3DProfile
2519
2520TProfile2D *TH3::DoProjectProfile2D(const char* name, const char * title, const TAxis* projX, const TAxis* projY,
2521 bool originalRange, bool useUF, bool useOF) const
2522{
2523 // Get the ranges where we will work.
2524 Int_t ixmin = std::max(projX->GetFirst(),1);
2525 Int_t ixmax = std::min(projX->GetLast(),projX->GetNbins());
2526 Int_t iymin = std::max(projY->GetFirst(),1);
2527 Int_t iymax = std::min(projY->GetLast(),projY->GetNbins());
2528
2529 Int_t nx = ixmax-ixmin+1;
2530 Int_t ny = iymax-iymin+1;
2531
2532 // Create the projected profiles
2533 TProfile2D *p2 = 0;
2534
2535 // Create the histogram, either reseting a preexisting one
2536 // Does an object with the same name exists?
2537 TObject *p2obj = gROOT->FindObject(name);
2538 if (p2obj && p2obj->InheritsFrom(TH1::Class())) {
2539 if (p2obj->IsA() != TProfile2D::Class() ) {
2540 Error("DoProjectProfile2D","Histogram with name %s must be a TProfile2D and is a %s",name,p2obj->ClassName());
2541 return 0;
2542 }
2543 p2 = (TProfile2D*)p2obj;
2544 // reset existing profile and re-set bins
2545 p2->Reset();
2546 const TArrayD *xbins = projX->GetXbins();
2547 const TArrayD *ybins = projY->GetXbins();
2548 if ( originalRange ) {
2549 p2->SetBins(projY->GetNbins(),projY->GetXmin(),projY->GetXmax()
2550 ,projX->GetNbins(),projX->GetXmin(),projX->GetXmax());
2551 // set bins for mixed axis do not exists - need to set afterwards the variable bins
2552 if (ybins->fN != 0)
2553 p2->GetXaxis()->Set(projY->GetNbins(),&ybins->fArray[iymin-1]);
2554 if (xbins->fN != 0)
2555 p2->GetYaxis()->Set(projX->GetNbins(),&xbins->fArray[ixmin-1]);
2556 } else {
2557 p2->SetBins(ny,projY->GetBinLowEdge(iymin),projY->GetBinUpEdge(iymax)
2558 ,nx,projX->GetBinLowEdge(ixmin),projX->GetBinUpEdge(ixmax));
2559 if (ybins->fN != 0)
2560 p2->GetXaxis()->Set(ny,&ybins->fArray[iymin-1]);
2561 if (xbins->fN != 0)
2562 p2->GetYaxis()->Set(nx,&xbins->fArray[ixmin-1]);
2563 }
2564 }
2565
2566 if (!p2) {
2567 const TArrayD *xbins = projX->GetXbins();
2568 const TArrayD *ybins = projY->GetXbins();
2569 if ( originalRange ) {
2570 if (xbins->fN == 0 && ybins->fN == 0) {
2571 p2 = new TProfile2D(name,title,projY->GetNbins(),projY->GetXmin(),projY->GetXmax()
2572 ,projX->GetNbins(),projX->GetXmin(),projX->GetXmax());
2573 } else if (ybins->fN == 0) {
2574 p2 = new TProfile2D(name,title,projY->GetNbins(),projY->GetXmin(),projY->GetXmax()
2575 ,projX->GetNbins(),&xbins->fArray[ixmin-1]);
2576 } else if (xbins->fN == 0) {
2577 p2 = new TProfile2D(name,title,projY->GetNbins(),&ybins->fArray[iymin-1]
2578 ,projX->GetNbins(),projX->GetXmin(),projX->GetXmax());
2579 } else {
2580 p2 = new TProfile2D(name,title,projY->GetNbins(),&ybins->fArray[iymin-1],projX->GetNbins(),&xbins->fArray[ixmin-1]);
2581 }
2582 } else {
2583 if (xbins->fN == 0 && ybins->fN == 0) {
2584 p2 = new TProfile2D(name,title,ny,projY->GetBinLowEdge(iymin),projY->GetBinUpEdge(iymax)
2585 ,nx,projX->GetBinLowEdge(ixmin),projX->GetBinUpEdge(ixmax));
2586 } else if (ybins->fN == 0) {
2587 p2 = new TProfile2D(name,title,ny,projY->GetBinLowEdge(iymin),projY->GetBinUpEdge(iymax)
2588 ,nx,&xbins->fArray[ixmin-1]);
2589 } else if (xbins->fN == 0) {
2590 p2 = new TProfile2D(name,title,ny,&ybins->fArray[iymin-1]
2591 ,nx,projX->GetBinLowEdge(ixmin),projX->GetBinUpEdge(ixmax));
2592 } else {
2593 p2 = new TProfile2D(name,title,ny,&ybins->fArray[iymin-1],nx,&xbins->fArray[ixmin-1]);
2594 }
2595 }
2596 }
2597
2598 // Set references to the axis, so that the loop has no branches.
2599 const TAxis* outAxis = 0;
2600 if ( projX != GetXaxis() && projY != GetXaxis() ) {
2601 outAxis = GetXaxis();
2602 } else if ( projX != GetYaxis() && projY != GetYaxis() ) {
2603 outAxis = GetYaxis();
2604 } else {
2605 outAxis = GetZaxis();
2606 }
2607
2608 // Weights management
2609 bool useWeights = (GetSumw2N() > 0);
2610 // store sum of w2 in profile if histo is weighted
2611 if (useWeights && (p2->GetBinSumw2()->fN != p2->GetNcells() ) ) p2->Sumw2();
2612
2613 // Set references to the bins, so that the loop has no branches.
2614 Int_t *refX = 0, *refY = 0, *refZ = 0;
2615 Int_t ixbin, iybin, outbin;
2616 if ( projX == GetXaxis() && projY == GetYaxis() ) { refX = &ixbin; refY = &iybin; refZ = &outbin; }
2617 if ( projX == GetYaxis() && projY == GetXaxis() ) { refX = &iybin; refY = &ixbin; refZ = &outbin; }
2618 if ( projX == GetXaxis() && projY == GetZaxis() ) { refX = &ixbin; refY = &outbin; refZ = &iybin; }
2619 if ( projX == GetZaxis() && projY == GetXaxis() ) { refX = &iybin; refY = &outbin; refZ = &ixbin; }
2620 if ( projX == GetYaxis() && projY == GetZaxis() ) { refX = &outbin; refY = &ixbin; refZ = &iybin; }
2621 if ( projX == GetZaxis() && projY == GetYaxis() ) { refX = &outbin; refY = &iybin; refZ = &ixbin; }
2622 R__ASSERT (refX != 0 && refY != 0 && refZ != 0);
2623
2624 Int_t outmin = outAxis->GetFirst();
2625 Int_t outmax = outAxis->GetLast();
2626 // GetFirst, GetLast can return underflow or overflow bins
2627 // correct for underflow/overflows
2628 if (useUF && !outAxis->TestBit(TAxis::kAxisRange) ) outmin -= 1;
2629 if (useOF && !outAxis->TestBit(TAxis::kAxisRange) ) outmax += 1;
2630
2631 TArrayD & binSumw2 = *(p2->GetBinSumw2());
2632 if (useWeights && binSumw2.fN <= 0) useWeights = false;
2633 if (!useWeights) p2->SetBit(TH1::kIsNotW);
2634
2635 // Call specific method for the projection
2636 for (ixbin=0;ixbin<=1+projX->GetNbins();ixbin++) {
2637 if ( (ixbin < ixmin || ixbin > ixmax) && projX->TestBit(TAxis::kAxisRange)) continue;
2638 for ( iybin=0;iybin<=1+projY->GetNbins();iybin++) {
2639 if ( (iybin < iymin || iybin > iymax) && projX->TestBit(TAxis::kAxisRange)) continue;
2640
2641 // profile output bin
2642 Int_t poutBin = p2->FindBin(projY->GetBinCenter(iybin), projX->GetBinCenter(ixbin));
2643 if (poutBin <0) continue;
2644 // loop on the bins to be integrated (outbin should be called inbin)
2645 for (outbin = outmin; outbin <= outmax; outbin++) {
2646
2647 Int_t bin = GetBin(*refX,*refY,*refZ);
2648
2649 //DoFillProfileProjection(p2, *projY, *projX, *outAxis, iybin, ixbin, outbin, bin, useWeights);
2650
2651 Double_t cont = RetrieveBinContent(bin);
2652 if (!cont) continue;
2653
2654 Double_t tmp = 0;
2655 // the following fill update wrongly the fBinSumw2- need to save it before
2656 if ( useWeights ) tmp = binSumw2.fArray[poutBin];
2657 p2->Fill( projY->GetBinCenter(iybin) , projX->GetBinCenter(ixbin), outAxis->GetBinCenter(outbin), cont);
2658 if (useWeights ) binSumw2.fArray[poutBin] = tmp + fSumw2.fArray[bin];
2659
2660 }
2661 }
2662 }
2663
2664 // recompute statistics for the projected profiles
2665 // forget about preserving old statistics
2666 bool resetStats = true;
2667 Double_t stats[kNstat];
2668 // reset statistics
2669 if (resetStats)
2670 for (Int_t i=0;i<kNstat;i++) stats[i] = 0;
2671
2672 p2->PutStats(stats);
2673 Double_t entries = fEntries;
2674 // recalculate the statistics
2675 if (resetStats) {
2676 entries = p2->GetEffectiveEntries();
2677 if (!useWeights) entries = TMath::Floor( entries + 0.5); // to avoid numerical rounding
2678 p2->SetEntries( entries );
2679 }
2680
2681 p2->SetEntries(entries);
2682
2683 return p2;
2684}
2685
2686
2687////////////////////////////////////////////////////////////////////////////////
2688/// Project a 3-d histogram into a 2-d profile histograms depending
2689/// on the option parameter
2690/// option may contain a combination of the characters x,y,z
2691/// option = "xy" return the x versus y projection into a TProfile2D histogram
2692/// option = "yx" return the y versus x projection into a TProfile2D histogram
2693/// option = "xz" return the x versus z projection into a TProfile2D histogram
2694/// option = "zx" return the z versus x projection into a TProfile2D histogram
2695/// option = "yz" return the y versus z projection into a TProfile2D histogram
2696/// option = "zy" return the z versus y projection into a TProfile2D histogram
2697/// NB: the notation "a vs b" means "a" vertical and "b" horizontal
2698///
2699/// option = "o" original axis range of the target axes will be
2700/// kept, but only bins inside the selected range will be filled.
2701///
2702/// The projection is made for the selected bins only.
2703/// To select a bin range along an axis, use TAxis::SetRange, eg
2704/// h3.GetYaxis()->SetRange(23,56);
2705///
2706/// NOTE 1: The generated histogram is named th3name + "_p" + option
2707/// eg if the TH3* h histogram is named "myhist", then
2708/// h->Project3D("xy"); produces a TProfile2D histogram named "myhist_pxy".
2709/// The following sequence
2710/// h->Project3DProfile("xy");
2711/// h->Project3DProfile("xy2");
2712/// will generate two TProfile2D histograms named "myhist_pxy" and "myhist_pxy2"
2713/// So, passing additional characters in the option string one can customize the name.
2714///
2715/// NOTE 2: If a profile of the same type already exists with compatible axes,
2716/// the profile is reset and filled again with the projected contents of the TH3.
2717/// In the case of axes incompatibility, an error is reported and a NULL pointer is returned.
2718///
2719/// NOTE 3: The number of entries in the projected profile is estimated from the number of
2720/// effective entries for all the cells included in the projection.
2721///
2722/// NOTE 4: underflow/overflow are by default excluded from the projection
2723/// (Note that this is a different default behavior compared to the projection to an histogram)
2724/// To include the underflow and/or overflow use option "UF" and/or "OF"
2725
2727{
2728 TString opt = option; opt.ToLower();
2729 Int_t pcase = 0;
2730 TString ptype;
2731 if (opt.Contains("xy")) { pcase = 4; ptype = "xy"; }
2732 if (opt.Contains("yx")) { pcase = 5; ptype = "yx"; }
2733 if (opt.Contains("xz")) { pcase = 6; ptype = "xz"; }
2734 if (opt.Contains("zx")) { pcase = 7; ptype = "zx"; }
2735 if (opt.Contains("yz")) { pcase = 8; ptype = "yz"; }
2736 if (opt.Contains("zy")) { pcase = 9; ptype = "zy"; }
2737
2738 if (pcase == 0) {
2739 Error("Project3D","No projection axis specified - return a NULL pointer");
2740 return 0;
2741 }
2742 // do not remove ptype from opt to use later in the projected histo name
2743
2744 Bool_t useUF = kFALSE;
2745 if (opt.Contains("uf") ) {
2746 useUF = kTRUE;
2747 opt.Remove(opt.Index("uf"),2);
2748 }
2749 Bool_t useOF = kFALSE;
2750 if (opt.Contains("of") ) {
2751 useOF = kTRUE;
2752 opt.Remove(opt.Index("of"),2);
2753 }
2754
2755 Bool_t originalRange = kFALSE;
2756 if (opt.Contains('o') ) {
2757 originalRange = kTRUE;
2758 opt.Remove(opt.First("o"),1);
2759 }
2760
2761 // Create the projected profile
2762 TProfile2D *p2 = 0;
2763 TString name = GetName();
2764 TString title = GetTitle();
2765 name += "_p"; name += opt; // opt may include a user defined name
2766 title += " profile "; title += ptype; title += " projection";
2767
2768 // Call the method with the specific projected axes.
2769 switch (pcase) {
2770 case 4:
2771 // "xy"
2772 p2 = DoProjectProfile2D(name, title, GetXaxis(), GetYaxis(), originalRange, useUF, useOF);
2773 break;
2774
2775 case 5:
2776 // "yx"
2777 p2 = DoProjectProfile2D(name, title, GetYaxis(), GetXaxis(), originalRange, useUF, useOF);
2778 break;
2779
2780 case 6:
2781 // "xz"
2782 p2 = DoProjectProfile2D(name, title, GetXaxis(), GetZaxis(), originalRange, useUF, useOF);
2783 break;
2784
2785 case 7:
2786 // "zx"
2787 p2 = DoProjectProfile2D(name, title, GetZaxis(), GetXaxis(), originalRange, useUF, useOF);
2788 break;
2789
2790 case 8:
2791 // "yz"
2792 p2 = DoProjectProfile2D(name, title, GetYaxis(), GetZaxis(), originalRange, useUF, useOF);
2793 break;
2794
2795 case 9:
2796 // "zy"
2797 p2 = DoProjectProfile2D(name, title, GetZaxis(), GetYaxis(), originalRange, useUF, useOF);
2798 break;
2799
2800 }
2801
2802 return p2;
2803}
2804
2805
2806////////////////////////////////////////////////////////////////////////////////
2807/// Replace current statistics with the values in array stats
2808
2810{
2811 TH1::PutStats(stats);
2812 fTsumwy = stats[4];
2813 fTsumwy2 = stats[5];
2814 fTsumwxy = stats[6];
2815 fTsumwz = stats[7];
2816 fTsumwz2 = stats[8];
2817 fTsumwxz = stats[9];
2818 fTsumwyz = stats[10];
2819}
2820
2821
2822////////////////////////////////////////////////////////////////////////////////
2823/// Rebin only the X axis
2824/// see Rebin3D
2825
2826TH3 *TH3::RebinX(Int_t ngroup, const char *newname)
2827{
2828 return Rebin3D(ngroup, 1, 1, newname);
2829}
2830
2831
2832////////////////////////////////////////////////////////////////////////////////
2833/// Rebin only the Y axis
2834/// see Rebin3D
2835
2836TH3 *TH3::RebinY(Int_t ngroup, const char *newname)
2837{
2838 return Rebin3D(1, ngroup, 1, newname);
2839}
2840
2841
2842////////////////////////////////////////////////////////////////////////////////
2843/// Rebin only the Z axis
2844/// see Rebin3D
2845
2846TH3 *TH3::RebinZ(Int_t ngroup, const char *newname)
2847{
2848 return Rebin3D(1, 1, ngroup, newname);
2849
2850}
2851
2852
2853////////////////////////////////////////////////////////////////////////////////
2854/// Rebin this histogram grouping nxgroup/nygroup/nzgroup bins along the xaxis/yaxis/zaxis together.
2855///
2856/// if newname is not blank a new temporary histogram hnew is created.
2857/// else the current histogram is modified (default)
2858/// The parameter nxgroup/nygroup indicate how many bins along the xaxis/yaxis of this
2859/// have to me merged into one bin of hnew
2860/// If the original histogram has errors stored (via Sumw2), the resulting
2861/// histograms has new errors correctly calculated.
2862///
2863/// examples: if hpxpy is an existing TH3 histogram with 40 x 40 x 40 bins
2864/// hpxpypz->Rebin3D(); // merges two bins along the xaxis and yaxis in one in hpxpypz
2865/// // Carefull: previous contents of hpxpy are lost
2866/// hpxpypz->RebinX(5); //merges five bins along the xaxis in one in hpxpypz
2867/// TH3 *hnew = hpxpypz->RebinY(5,"hnew"); // creates a new histogram hnew
2868/// // merging 5 bins of h1 along the yaxis in one bin
2869///
2870/// NOTE : If nxgroup/nygroup is not an exact divider of the number of bins,
2871/// along the xaxis/yaxis the top limit(s) of the rebinned histogram
2872/// is changed to the upper edge of the xbin=newxbins*nxgroup resp.
2873/// ybin=newybins*nygroup and the corresponding bins are added to
2874/// the overflow bin.
2875/// Statistics will be recomputed from the new bin contents.
2876
2877TH3 *TH3::Rebin3D(Int_t nxgroup, Int_t nygroup, Int_t nzgroup, const char *newname)
2878{
2879 Int_t i,j,k,xbin,ybin,zbin;
2880 Int_t nxbins = fXaxis.GetNbins();
2881 Int_t nybins = fYaxis.GetNbins();
2882 Int_t nzbins = fZaxis.GetNbins();
2887 Double_t zmin = fZaxis.GetXmin();
2888 Double_t zmax = fZaxis.GetXmax();
2889 if ((nxgroup <= 0) || (nxgroup > nxbins)) {
2890 Error("Rebin", "Illegal value of nxgroup=%d",nxgroup);
2891 return 0;
2892 }
2893 if ((nygroup <= 0) || (nygroup > nybins)) {
2894 Error("Rebin", "Illegal value of nygroup=%d",nygroup);
2895 return 0;
2896 }
2897 if ((nzgroup <= 0) || (nzgroup > nzbins)) {
2898 Error("Rebin", "Illegal value of nzgroup=%d",nzgroup);
2899 return 0;
2900 }
2901
2902 Int_t newxbins = nxbins/nxgroup;
2903 Int_t newybins = nybins/nygroup;
2904 Int_t newzbins = nzbins/nzgroup;
2905
2906 // Save old bin contents into a new array
2907 Double_t entries = fEntries;
2908 Double_t *oldBins = new Double_t[fNcells];
2909 for (Int_t ibin = 0; ibin < fNcells; ibin++) {
2910 oldBins[ibin] = RetrieveBinContent(ibin);
2911 }
2912 Double_t *oldSumw2 = 0;
2913 if (fSumw2.fN != 0) {
2914 oldSumw2 = new Double_t[fNcells];
2915 for (Int_t ibin = 0; ibin < fNcells; ibin++) {
2916 oldSumw2[ibin] = fSumw2.fArray[ibin];
2917 }
2918 }
2919
2920 // create a clone of the old histogram if newname is specified
2921 TH3 *hnew = this;
2922 if (newname && strlen(newname)) {
2923 hnew = (TH3*)Clone();
2924 hnew->SetName(newname);
2925 }
2926
2927 // save original statistics
2928 Double_t stat[kNstat];
2929 GetStats(stat);
2930 bool resetStat = false;
2931
2932
2933 // change axis specs and rebuild bin contents array
2934 if (newxbins*nxgroup != nxbins) {
2935 xmax = fXaxis.GetBinUpEdge(newxbins*nxgroup);
2936 resetStat = true; //stats must be reset because top bins will be moved to overflow bin
2937 }
2938 if (newybins*nygroup != nybins) {
2939 ymax = fYaxis.GetBinUpEdge(newybins*nygroup);
2940 resetStat = true; //stats must be reset because top bins will be moved to overflow bin
2941 }
2942 if (newzbins*nzgroup != nzbins) {
2943 zmax = fZaxis.GetBinUpEdge(newzbins*nzgroup);
2944 resetStat = true; //stats must be reset because top bins will be moved to overflow bin
2945 }
2946 // save the TAttAxis members (reset by SetBins) for x axis
2947 Int_t nXdivisions = fXaxis.GetNdivisions();
2948 Color_t xAxisColor = fXaxis.GetAxisColor();
2949 Color_t xLabelColor = fXaxis.GetLabelColor();
2950 Style_t xLabelFont = fXaxis.GetLabelFont();
2951 Float_t xLabelOffset = fXaxis.GetLabelOffset();
2952 Float_t xLabelSize = fXaxis.GetLabelSize();
2953 Float_t xTickLength = fXaxis.GetTickLength();
2954 Float_t xTitleOffset = fXaxis.GetTitleOffset();
2955 Float_t xTitleSize = fXaxis.GetTitleSize();
2956 Color_t xTitleColor = fXaxis.GetTitleColor();
2957 Style_t xTitleFont = fXaxis.GetTitleFont();
2958 // save the TAttAxis members (reset by SetBins) for y axis
2959 Int_t nYdivisions = fYaxis.GetNdivisions();
2960 Color_t yAxisColor = fYaxis.GetAxisColor();
2961 Color_t yLabelColor = fYaxis.GetLabelColor();
2962 Style_t yLabelFont = fYaxis.GetLabelFont();
2963 Float_t yLabelOffset = fYaxis.GetLabelOffset();
2964 Float_t yLabelSize = fYaxis.GetLabelSize();
2965 Float_t yTickLength = fYaxis.GetTickLength();
2966 Float_t yTitleOffset = fYaxis.GetTitleOffset();
2967 Float_t yTitleSize = fYaxis.GetTitleSize();
2968 Color_t yTitleColor = fYaxis.GetTitleColor();
2969 Style_t yTitleFont = fYaxis.GetTitleFont();
2970 // save the TAttAxis members (reset by SetBins) for z axis
2971 Int_t nZdivisions = fZaxis.GetNdivisions();
2972 Color_t zAxisColor = fZaxis.GetAxisColor();
2973 Color_t zLabelColor = fZaxis.GetLabelColor();
2974 Style_t zLabelFont = fZaxis.GetLabelFont();
2975 Float_t zLabelOffset = fZaxis.GetLabelOffset();
2976 Float_t zLabelSize = fZaxis.GetLabelSize();
2977 Float_t zTickLength = fZaxis.GetTickLength();
2978 Float_t zTitleOffset = fZaxis.GetTitleOffset();
2979 Float_t zTitleSize = fZaxis.GetTitleSize();
2980 Color_t zTitleColor = fZaxis.GetTitleColor();
2981 Style_t zTitleFont = fZaxis.GetTitleFont();
2982
2983 // copy merged bin contents (ignore under/overflows)
2984 if (nxgroup != 1 || nygroup != 1 || nzgroup != 1) {
2985 if (fXaxis.GetXbins()->GetSize() > 0 || fYaxis.GetXbins()->GetSize() > 0 || fZaxis.GetXbins()->GetSize() > 0) {
2986 // variable bin sizes in x or y, don't treat both cases separately
2987 Double_t *xbins = new Double_t[newxbins+1];
2988 for (i = 0; i <= newxbins; ++i) xbins[i] = fXaxis.GetBinLowEdge(1+i*nxgroup);
2989 Double_t *ybins = new Double_t[newybins+1];
2990 for (i = 0; i <= newybins; ++i) ybins[i] = fYaxis.GetBinLowEdge(1+i*nygroup);
2991 Double_t *zbins = new Double_t[newzbins+1];
2992 for (i = 0; i <= newzbins; ++i) zbins[i] = fZaxis.GetBinLowEdge(1+i*nzgroup);
2993 hnew->SetBins(newxbins,xbins, newybins, ybins, newzbins, zbins);//changes also errors array (if any)
2994 delete [] xbins;
2995 delete [] ybins;
2996 delete [] zbins;
2997 } else {
2998 hnew->SetBins(newxbins, xmin, xmax, newybins, ymin, ymax, newzbins, zmin, zmax);//changes also errors array
2999 }
3000
3001 Double_t binContent, binSumw2;
3002 Int_t oldxbin = 1;
3003 Int_t oldybin = 1;
3004 Int_t oldzbin = 1;
3005 Int_t bin;
3006 for (xbin = 1; xbin <= newxbins; xbin++) {
3007 oldybin=1;
3008 oldzbin=1;
3009 for (ybin = 1; ybin <= newybins; ybin++) {
3010 oldzbin=1;
3011 for (zbin = 1; zbin <= newzbins; zbin++) {
3012 binContent = 0;
3013 binSumw2 = 0;
3014 for (i = 0; i < nxgroup; i++) {
3015 if (oldxbin+i > nxbins) break;
3016 for (j =0; j < nygroup; j++) {
3017 if (oldybin+j > nybins) break;
3018 for (k =0; k < nzgroup; k++) {
3019 if (oldzbin+k > nzbins) break;
3020 //get global bin (same conventions as in TH1::GetBin(xbin,ybin)
3021 bin = oldxbin + i + (oldybin + j)*(nxbins + 2) + (oldzbin + k)*(nxbins + 2)*(nybins + 2);
3022 binContent += oldBins[bin];
3023 if (oldSumw2) binSumw2 += oldSumw2[bin];
3024 }
3025 }
3026 }
3027 Int_t ibin = hnew->GetBin(xbin,ybin,zbin); // new bin number
3028 hnew->SetBinContent(ibin, binContent);
3029 if (oldSumw2) hnew->fSumw2.fArray[ibin] = binSumw2;
3030 oldzbin += nzgroup;
3031 }
3032 oldybin += nygroup;
3033 }
3034 oldxbin += nxgroup;
3035 }
3036
3037 // compute new underflow/overflows for the 8 vertices
3038 for (Int_t xover = 0; xover <= 1; xover++) {
3039 for (Int_t yover = 0; yover <= 1; yover++) {
3040 for (Int_t zover = 0; zover <= 1; zover++) {
3041 binContent = 0;
3042 binSumw2 = 0;
3043 // make loop in case of only underflow/overflow
3044 for (xbin = xover*oldxbin; xbin <= xover*(nxbins+1); xbin++) {
3045 for (ybin = yover*oldybin; ybin <= yover*(nybins+1); ybin++) {
3046 for (zbin = zover*oldzbin; zbin <= zover*(nzbins+1); zbin++) {
3047 bin = GetBin(xbin,ybin,zbin);
3048 binContent += oldBins[bin];
3049 if (oldSumw2) binSumw2 += oldSumw2[bin];
3050 }
3051 }
3052 }
3053 Int_t binNew = hnew->GetBin( xover *(newxbins+1),
3054 yover*(newybins+1), zover*(newzbins+1) );
3055 hnew->SetBinContent(binNew,binContent);
3056 if (oldSumw2) hnew->fSumw2.fArray[binNew] = binSumw2;
3057 }
3058 }
3059 }
3060
3061 Double_t binContent0, binContent2, binContent3, binContent4;
3062 Double_t binError0, binError2, binError3, binError4;
3063 Int_t oldxbin2, oldybin2, oldzbin2;
3064 Int_t ufbin, ofbin, ofbin2, ofbin3, ofbin4;
3065
3066 // recompute under/overflow contents in y for the new x and z bins
3067 oldxbin2 = 1;
3068 oldybin2 = 1;
3069 oldzbin2 = 1;
3070 for (xbin = 1; xbin<=newxbins; xbin++) {
3071 oldzbin2 = 1;
3072 for (zbin = 1; zbin<=newzbins; zbin++) {
3073 binContent0 = binContent2 = 0;
3074 binError0 = binError2 = 0;
3075 for (i=0; i<nxgroup; i++) {
3076 if (oldxbin2+i > nxbins) break;
3077 for (k=0; k<nzgroup; k++) {
3078 if (oldzbin2+k > nzbins) break;
3079 //old underflow bin (in y)
3080 ufbin = oldxbin2 + i + (nxbins+2)*(nybins+2)*(oldzbin2+k);
3081 binContent0 += oldBins[ufbin];
3082 if (oldSumw2) binError0 += oldSumw2[ufbin];
3083 for (ybin = oldybin; ybin <= nybins + 1; ybin++) {
3084 //old overflow bin (in y)
3085 ofbin = ufbin + ybin*(nxbins+2);
3086 binContent2 += oldBins[ofbin];
3087 if (oldSumw2) binError2 += oldSumw2[ofbin];
3088 }
3089 }
3090 }
3091 hnew->SetBinContent(xbin,0,zbin,binContent0);
3092 hnew->SetBinContent(xbin,newybins+1,zbin,binContent2);
3093 if (oldSumw2) {
3094 hnew->SetBinError(xbin,0,zbin,TMath::Sqrt(binError0));
3095 hnew->SetBinError(xbin,newybins+1,zbin,TMath::Sqrt(binError2) );
3096 }
3097 oldzbin2 += nzgroup;
3098 }
3099 oldxbin2 += nxgroup;
3100 }
3101
3102 // recompute under/overflow contents in x for the new y and z bins
3103 oldxbin2 = 1;
3104 oldybin2 = 1;
3105 oldzbin2 = 1;
3106 for (ybin = 1; ybin<=newybins; ybin++) {
3107 oldzbin2 = 1;
3108 for (zbin = 1; zbin<=newzbins; zbin++) {
3109 binContent0 = binContent2 = 0;
3110 binError0 = binError2 = 0;
3111 for (j=0; j<nygroup; j++) {
3112 if (oldybin2+j > nybins) break;
3113 for (k=0; k<nzgroup; k++) {
3114 if (oldzbin2+k > nzbins) break;
3115 //old underflow bin (in y)
3116 ufbin = (oldybin2 + j)*(nxbins+2) + (nxbins+2)*(nybins+2)*(oldzbin2+k);
3117 binContent0 += oldBins[ufbin];
3118 if (oldSumw2) binError0 += oldSumw2[ufbin];
3119 for (xbin = oldxbin; xbin <= nxbins + 1; xbin++) {
3120 //old overflow bin (in x)
3121 ofbin = ufbin + xbin;
3122 binContent2 += oldBins[ofbin];
3123 if (oldSumw2) binError2 += oldSumw2[ofbin];
3124 }
3125 }
3126 }
3127 hnew->SetBinContent(0,ybin,zbin,binContent0);
3128 hnew->SetBinContent(newxbins+1,ybin,zbin,binContent2);
3129 if (oldSumw2) {
3130 hnew->SetBinError(0,ybin,zbin,TMath::Sqrt(binError0));
3131 hnew->SetBinError(newxbins+1,ybin,zbin,TMath::Sqrt(binError2) );
3132 }
3133 oldzbin2 += nzgroup;
3134 }
3135 oldybin2 += nygroup;
3136 }
3137
3138 // recompute under/overflow contents in z for the new x and y bins
3139 oldxbin2 = 1;
3140 oldybin2 = 1;
3141 oldzbin2 = 1;
3142 for (xbin = 1; xbin<=newxbins; xbin++) {
3143 oldybin2 = 1;
3144 for (ybin = 1; ybin<=newybins; ybin++) {
3145 binContent0 = binContent2 = 0;
3146 binError0 = binError2 = 0;
3147 for (i=0; i<nxgroup; i++) {
3148 if (oldxbin2+i > nxbins) break;
3149 for (j=0; j<nygroup; j++) {
3150 if (oldybin2+j > nybins) break;
3151 //old underflow bin (in z)
3152 ufbin = oldxbin2 + i + (nxbins+2)*(oldybin2+j);
3153 binContent0 += oldBins[ufbin];
3154 if (oldSumw2) binError0 += oldSumw2[ufbin];
3155 for (zbin = oldzbin; zbin <= nzbins + 1; zbin++) {
3156 //old overflow bin (in z)
3157 ofbin = ufbin + (nxbins+2)*(nybins+2)*zbin;
3158 binContent2 += oldBins[ofbin];
3159 if (oldSumw2) binError2 += oldSumw2[ofbin];
3160 }
3161 }
3162 }
3163 hnew->SetBinContent(xbin,ybin,0,binContent0);
3164 hnew->SetBinContent(xbin,ybin,newzbins+1,binContent2);
3165 if (oldSumw2) {
3166 hnew->SetBinError(xbin,ybin,0,TMath::Sqrt(binError0));
3167 hnew->SetBinError(xbin,ybin,newzbins+1,TMath::Sqrt(binError2) );
3168 }
3169 oldybin2 += nygroup;
3170 }
3171 oldxbin2 += nxgroup;
3172 }
3173
3174 // recompute under/overflow contents in y, z for the new x
3175 oldxbin2 = 1;
3176 oldybin2 = 1;
3177 oldzbin2 = 1;
3178 for (xbin = 1; xbin<=newxbins; xbin++) {
3179 binContent0 = 0;
3180 binContent2 = 0;
3181 binContent3 = 0;
3182 binContent4 = 0;
3183 binError0 = 0;
3184 binError2 = 0;
3185 binError3 = 0;
3186 binError4 = 0;
3187 for (i=0; i<nxgroup; i++) {
3188 if (oldxbin2+i > nxbins) break;
3189 ufbin = oldxbin2 + i; //
3190 binContent0 += oldBins[ufbin];
3191 if (oldSumw2) binError0 += oldSumw2[ufbin];
3192 for (ybin = oldybin; ybin <= nybins + 1; ybin++) {
3193 ofbin3 = ufbin+ybin*(nxbins+2);
3194 binContent3 += oldBins[ ofbin3 ];
3195 if (oldSumw2) binError3 += oldSumw2[ofbin3];
3196 for (zbin = oldzbin; zbin <= nzbins + 1; zbin++) {
3197 //old overflow bin (in z)
3198 ofbin4 = oldxbin2 + i + ybin*(nxbins+2) + (nxbins+2)*(nybins+2)*zbin;
3199 binContent4 += oldBins[ofbin4];
3200 if (oldSumw2) binError4 += oldSumw2[ofbin4];
3201 }
3202 }
3203 for (zbin = oldzbin; zbin <= nzbins + 1; zbin++) {
3204 ofbin2 = ufbin+zbin*(nxbins+2)*(nybins+2);
3205 binContent2 += oldBins[ ofbin2 ];
3206 if (oldSumw2) binError2 += oldSumw2[ofbin2];
3207 }
3208 }
3209 hnew->SetBinContent(xbin,0,0,binContent0);
3210 hnew->SetBinContent(xbin,0,newzbins+1,binContent2);
3211 hnew->SetBinContent(xbin,newybins+1,0,binContent3);
3212 hnew->SetBinContent(xbin,newybins+1,newzbins+1,binContent4);
3213 if (oldSumw2) {
3214 hnew->SetBinError(xbin,0,0,TMath::Sqrt(binError0));
3215 hnew->SetBinError(xbin,0,newzbins+1,TMath::Sqrt(binError2) );
3216 hnew->SetBinError(xbin,newybins+1,0,TMath::Sqrt(binError3) );
3217 hnew->SetBinError(xbin,newybins+1,newzbins+1,TMath::Sqrt(binError4) );
3218 }
3219 oldxbin2 += nxgroup;
3220 }
3221
3222 // recompute under/overflow contents in x, y for the new z
3223 oldxbin2 = 1;
3224 oldybin2 = 1;
3225 oldzbin2 = 1;
3226 for (zbin = 1; zbin<=newzbins; zbin++) {
3227 binContent0 = 0;
3228 binContent2 = 0;
3229 binContent3 = 0;
3230 binContent4 = 0;
3231 binError0 = 0;
3232 binError2 = 0;
3233 binError3 = 0;
3234 binError4 = 0;
3235 for (i=0; i<nzgroup; i++) {
3236 if (oldzbin2+i > nzbins) break;
3237 ufbin = (oldzbin2 + i)*(nxbins+2)*(nybins+2); //
3238 binContent0 += oldBins[ufbin];
3239 if (oldSumw2) binError0 += oldSumw2[ufbin];
3240 for (ybin = oldybin; ybin <= nybins + 1; ybin++) {
3241 ofbin3 = ufbin+ybin*(nxbins+2);
3242 binContent3 += oldBins[ ofbin3 ];
3243 if (oldSumw2) binError3 += oldSumw2[ofbin3];
3244 for (xbin = oldxbin; xbin <= nxbins + 1; xbin++) {
3245 //old overflow bin (in z)
3246 ofbin4 = ufbin + xbin + ybin*(nxbins+2);
3247 binContent4 += oldBins[ofbin4];
3248 if (oldSumw2) binError4 += oldSumw2[ofbin4];
3249 }
3250 }
3251 for (xbin = oldxbin; xbin <= nxbins + 1; xbin++) {
3252 ofbin2 = xbin +(oldzbin2+i)*(nxbins+2)*(nybins+2);
3253 binContent2 += oldBins[ ofbin2 ];
3254 if (oldSumw2) binError2 += oldSumw2[ofbin2];
3255 }
3256 }
3257 hnew->SetBinContent(0,0,zbin,binContent0);
3258 hnew->SetBinContent(0,newybins+1,zbin,binContent3);
3259 hnew->SetBinContent(newxbins+1,0,zbin,binContent2);
3260 hnew->SetBinContent(newxbins+1,newybins+1,zbin,binContent4);
3261 if (oldSumw2) {
3262 hnew->SetBinError(0,0,zbin,TMath::Sqrt(binError0));
3263 hnew->SetBinError(0,newybins+1,zbin,TMath::Sqrt(binError3) );
3264 hnew->SetBinError(newxbins+1,0,zbin,TMath::Sqrt(binError2) );
3265 hnew->SetBinError(newxbins+1,newybins+1,zbin,TMath::Sqrt(binError4) );
3266 }
3267 oldzbin2 += nzgroup;
3268 }
3269
3270 // recompute under/overflow contents in x, z for the new y
3271 oldxbin2 = 1;
3272 oldybin2 = 1;
3273 oldzbin2 = 1;
3274 for (ybin = 1; ybin<=newybins; ybin++) {
3275 binContent0 = 0;
3276 binContent2 = 0;
3277 binContent3 = 0;
3278 binContent4 = 0;
3279 binError0 = 0;
3280 binError2 = 0;
3281 binError3 = 0;
3282 binError4 = 0;
3283 for (i=0; i<nygroup; i++) {
3284 if (oldybin2+i > nybins) break;
3285 ufbin = (oldybin2 + i)*(nxbins+2); //
3286 binContent0 += oldBins[ufbin];
3287 if (oldSumw2) binError0 += oldSumw2[ufbin];
3288 for (xbin = oldxbin; xbin <= nxbins + 1; xbin++) {
3289 ofbin3 = ufbin+xbin;
3290 binContent3 += oldBins[ ofbin3 ];
3291 if (oldSumw2) binError3 += oldSumw2[ofbin3];
3292 for (zbin = oldzbin; zbin <= nzbins + 1; zbin++) {
3293 //old overflow bin (in z)
3294 ofbin4 = xbin + (nxbins+2)*(nybins+2)*zbin+(oldybin2+i)*(nxbins+2);
3295 binContent4 += oldBins[ofbin4];
3296 if (oldSumw2) binError4 += oldSumw2[ofbin4];
3297 }
3298 }
3299 for (zbin = oldzbin; zbin <= nzbins + 1; zbin++) {
3300 ofbin2 = (oldybin2+i)*(nxbins+2)+zbin*(nxbins+2)*(nybins+2);
3301 binContent2 += oldBins[ ofbin2 ];
3302 if (oldSumw2) binError2 += oldSumw2[ofbin2];
3303 }
3304 }
3305 hnew->SetBinContent(0,ybin,0,binContent0);
3306 hnew->SetBinContent(0,ybin,newzbins+1,binContent2);
3307 hnew->SetBinContent(newxbins+1,ybin,0,binContent3);
3308 hnew->SetBinContent(newxbins+1,ybin,newzbins+1,binContent4);
3309 if (oldSumw2) {
3310 hnew->SetBinError(0,ybin,0,TMath::Sqrt(binError0));
3311 hnew->SetBinError(0,ybin,newzbins+1,TMath::Sqrt(binError2) );
3312 hnew->SetBinError(newxbins+1,ybin,0,TMath::Sqrt(binError3) );
3313 hnew->SetBinError(newxbins+1,ybin,newzbins+1,TMath::Sqrt(binError4) );
3314 }
3315 oldybin2 += nygroup;
3316 }
3317 }
3318
3319 // Restore x axis attributes
3320 fXaxis.SetNdivisions(nXdivisions);
3321 fXaxis.SetAxisColor(xAxisColor);
3322 fXaxis.SetLabelColor(xLabelColor);
3323 fXaxis.SetLabelFont(xLabelFont);
3324 fXaxis.SetLabelOffset(xLabelOffset);
3325 fXaxis.SetLabelSize(xLabelSize);
3326 fXaxis.SetTickLength(xTickLength);
3327 fXaxis.SetTitleOffset(xTitleOffset);
3328 fXaxis.SetTitleSize(xTitleSize);
3329 fXaxis.SetTitleColor(xTitleColor);
3330 fXaxis.SetTitleFont(xTitleFont);
3331 // Restore y axis attributes
3332 fYaxis.SetNdivisions(nYdivisions);
3333 fYaxis.SetAxisColor(yAxisColor);
3334 fYaxis.SetLabelColor(yLabelColor);
3335 fYaxis.SetLabelFont(yLabelFont);
3336 fYaxis.SetLabelOffset(yLabelOffset);
3337 fYaxis.SetLabelSize(yLabelSize);
3338 fYaxis.SetTickLength(yTickLength);
3339 fYaxis.SetTitleOffset(yTitleOffset);
3340 fYaxis.SetTitleSize(yTitleSize);
3341 fYaxis.SetTitleColor(yTitleColor);
3342 fYaxis.SetTitleFont(yTitleFont);
3343 // Restore z axis attributes
3344 fZaxis.SetNdivisions(nZdivisions);
3345 fZaxis.SetAxisColor(zAxisColor);
3346 fZaxis.SetLabelColor(zLabelColor);
3347 fZaxis.SetLabelFont(zLabelFont);
3348 fZaxis.SetLabelOffset(zLabelOffset);
3349 fZaxis.SetLabelSize(zLabelSize);
3350 fZaxis.SetTickLength(zTickLength);
3351 fZaxis.SetTitleOffset(zTitleOffset);
3352 fZaxis.SetTitleSize(zTitleSize);
3353 fZaxis.SetTitleColor(zTitleColor);
3354 fZaxis.SetTitleFont(zTitleFont);
3355
3356 //restore statistics and entries modified by SetBinContent
3357 hnew->SetEntries(entries);
3358 if (!resetStat) hnew->PutStats(stat);
3359
3360 delete [] oldBins;
3361 if (oldSumw2) delete [] oldSumw2;
3362 return hnew;
3363}
3364
3365
3366////////////////////////////////////////////////////////////////////////////////
3367/// Reset this histogram: contents, errors, etc.
3368
3370{
3371 TH1::Reset(option);
3372 TString opt = option;
3373 opt.ToUpper();
3374 if (opt.Contains("ICE") && !opt.Contains("S")) return;
3375 fTsumwy = 0;
3376 fTsumwy2 = 0;
3377 fTsumwxy = 0;
3378 fTsumwz = 0;
3379 fTsumwz2 = 0;
3380 fTsumwxz = 0;
3381 fTsumwyz = 0;
3382}
3383
3384
3385////////////////////////////////////////////////////////////////////////////////
3386/// Set bin content.
3387
3389{
3390 fEntries++;
3391 fTsumw = 0;
3392 if (bin < 0) return;
3393 if (bin >= fNcells) return;
3394 UpdateBinContent(bin, content);
3395}
3396
3397
3398////////////////////////////////////////////////////////////////////////////////
3399/// Stream an object of class TH3.
3400
3401void TH3::Streamer(TBuffer &R__b)
3402{
3403 if (R__b.IsReading()) {
3404 UInt_t R__s, R__c;
3405 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
3406 if (R__v > 2) {
3407 R__b.ReadClassBuffer(TH3::Class(), this, R__v, R__s, R__c);
3408 return;
3409 }
3410 //====process old versions before automatic schema evolution
3411 TH1::Streamer(R__b);
3412 TAtt3D::Streamer(R__b);
3413 R__b.CheckByteCount(R__s, R__c, TH3::IsA());
3414 //====end of old versions
3415
3416 } else {
3417 R__b.WriteClassBuffer(TH3::Class(),this);
3418 }
3419}
3420
3421
3422//______________________________________________________________________________
3423// TH3C methods
3424// TH3C a 3-D histogram with one byte per cell (char)
3425//______________________________________________________________________________
3426
3427ClassImp(TH3C);
3428
3429
3430////////////////////////////////////////////////////////////////////////////////
3431/// Constructor.
3432
3434{
3435 SetBinsLength(27);
3436 if (fgDefaultSumw2) Sumw2();
3437}
3438
3439
3440////////////////////////////////////////////////////////////////////////////////
3441/// Destructor.
3442
3444{
3445}
3446
3447
3448////////////////////////////////////////////////////////////////////////////////
3449/// Constructor for fix bin size 3-D histograms
3450/// (see TH3::TH3 for explanation of parameters)
3451
3452TH3C::TH3C(const char *name,const char *title,Int_t nbinsx,Double_t xlow,Double_t xup
3453 ,Int_t nbinsy,Double_t ylow,Double_t yup
3454 ,Int_t nbinsz,Double_t zlow,Double_t zup)
3455 :TH3(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup,nbinsz,zlow,zup)
3456{
3458 if (fgDefaultSumw2) Sumw2();
3459
3460 if (xlow >= xup || ylow >= yup || zlow >= zup) SetBuffer(fgBufferSize);
3461}
3462
3463
3464////////////////////////////////////////////////////////////////////////////////
3465/// Constructor for variable bin size 3-D histograms.
3466/// (see TH3::TH3 for explanation of parameters)
3467
3468TH3C::TH3C(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins
3469 ,Int_t nbinsy,const Float_t *ybins
3470 ,Int_t nbinsz,const Float_t *zbins)
3471 :TH3(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins)
3472{
3474 if (fgDefaultSumw2) Sumw2();
3475}
3476
3477
3478////////////////////////////////////////////////////////////////////////////////
3479/// Constructor for variable bin size 3-D histograms.
3480/// (see TH3::TH3 for explanation of parameters)
3481
3482TH3C::TH3C(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins
3483 ,Int_t nbinsy,const Double_t *ybins
3484 ,Int_t nbinsz,const Double_t *zbins)
3485 :TH3(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins)
3486{
3488 if (fgDefaultSumw2) Sumw2();
3489}
3490
3491
3492////////////////////////////////////////////////////////////////////////////////
3493/// Copy constructor.
3494/// The list of functions is not copied. (Use Clone() if needed)
3495
3496TH3C::TH3C(const TH3C &h3c) : TH3(), TArrayC()
3497{
3498 ((TH3C&)h3c).Copy(*this);
3499}
3500
3501
3502////////////////////////////////////////////////////////////////////////////////
3503/// Increment bin content by 1.
3504
3506{
3507 if (fArray[bin] < 127) fArray[bin]++;
3508}
3509
3510
3511////////////////////////////////////////////////////////////////////////////////
3512/// Increment bin content by w.
3513
3515{
3516 Int_t newval = fArray[bin] + Int_t(w);
3517 if (newval > -128 && newval < 128) {fArray[bin] = Char_t(newval); return;}
3518 if (newval < -127) fArray[bin] = -127;
3519 if (newval > 127) fArray[bin] = 127;
3520}
3521
3522
3523////////////////////////////////////////////////////////////////////////////////
3524/// Copy this 3-D histogram structure to newth3.
3525
3526void TH3C::Copy(TObject &newth3) const
3527{
3528 TH3::Copy((TH3C&)newth3);
3529}
3530
3531
3532////////////////////////////////////////////////////////////////////////////////
3533/// Reset this histogram: contents, errors, etc.
3534
3536{
3537 TH3::Reset(option);
3539 // should also reset statistics once statistics are implemented for TH3
3540}
3541
3542
3543////////////////////////////////////////////////////////////////////////////////
3544/// Set total number of bins including under/overflow
3545/// Reallocate bin contents array
3546
3548{
3549 if (n < 0) n = (fXaxis.GetNbins()+2)*(fYaxis.GetNbins()+2)*(fZaxis.GetNbins()+2);
3550 fNcells = n;
3551 TArrayC::Set(n);
3552}
3553
3554
3555////////////////////////////////////////////////////////////////////////////////
3556/// When the mouse is moved in a pad containing a 3-d view of this histogram
3557/// a second canvas shows a projection type given as option.
3558/// To stop the generation of the projections, delete the canvas
3559/// containing the projection.
3560/// option may contain a combination of the characters x,y,z,e
3561/// option = "x" return the x projection into a TH1D histogram
3562/// option = "y" return the y projection into a TH1D histogram
3563/// option = "z" return the z projection into a TH1D histogram
3564/// option = "xy" return the x versus y projection into a TH2D histogram
3565/// option = "yx" return the y versus x projection into a TH2D histogram
3566/// option = "xz" return the x versus z projection into a TH2D histogram
3567/// option = "zx" return the z versus x projection into a TH2D histogram
3568/// option = "yz" return the y versus z projection into a TH2D histogram
3569/// option = "zy" return the z versus y projection into a TH2D histogram
3570/// option can also include the drawing option for the projection, eg to draw
3571/// the xy projection using the draw option "box" do
3572/// myhist.SetShowProjection("xy box");
3573/// This function is typically called from the context menu.
3574/// NB: the notation "a vs b" means "a" vertical and "b" horizontal
3575
3576void TH3::SetShowProjection(const char *option,Int_t nbins)
3577{
3578 GetPainter();
3579
3580 if (fPainter) fPainter->SetShowProjection(option,nbins);
3581}
3582
3583
3584////////////////////////////////////////////////////////////////////////////////
3585/// Stream an object of class TH3C.
3586
3587void TH3C::Streamer(TBuffer &R__b)
3588{
3589 if (R__b.IsReading()) {
3590 UInt_t R__s, R__c;
3591 if (R__b.GetParent() && R__b.GetVersionOwner() < 22300) return;
3592 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
3593 if (R__v > 2) {
3594 R__b.ReadClassBuffer(TH3C::Class(), this, R__v, R__s, R__c);
3595 return;
3596 }
3597 //====process old versions before automatic schema evolution
3598 if (R__v < 2) {
3599 R__b.ReadVersion();
3600 TH1::Streamer(R__b);
3601 TArrayC::Streamer(R__b);
3602 R__b.ReadVersion(&R__s, &R__c);
3603 TAtt3D::Streamer(R__b);
3604 } else {
3605 TH3::Streamer(R__b);
3606 TArrayC::Streamer(R__b);
3607 R__b.CheckByteCount(R__s, R__c, TH3C::IsA());
3608 }
3609 //====end of old versions
3610
3611 } else {
3612 R__b.WriteClassBuffer(TH3C::Class(),this);
3613 }
3614}
3615
3616
3617////////////////////////////////////////////////////////////////////////////////
3618/// Operator =
3619
3621{
3622 if (this != &h1) ((TH3C&)h1).Copy(*this);
3623 return *this;
3624}
3625
3626
3627////////////////////////////////////////////////////////////////////////////////
3628/// Operator *
3629
3631{
3632 TH3C hnew = h1;
3633 hnew.Scale(c1);
3634 hnew.SetDirectory(0);
3635 return hnew;
3636}
3637
3638
3639////////////////////////////////////////////////////////////////////////////////
3640/// Operator +
3641
3643{
3644 TH3C hnew = h1;
3645 hnew.Add(&h2,1);
3646 hnew.SetDirectory(0);
3647 return hnew;
3648}
3649
3650
3651////////////////////////////////////////////////////////////////////////////////
3652/// Operator -
3653
3655{
3656 TH3C hnew = h1;
3657 hnew.Add(&h2,-1);
3658 hnew.SetDirectory(0);
3659 return hnew;
3660}
3661
3662
3663////////////////////////////////////////////////////////////////////////////////
3664/// Operator *
3665
3667{
3668 TH3C hnew = h1;
3669 hnew.Multiply(&h2);
3670 hnew.SetDirectory(0);
3671 return hnew;
3672}
3673
3674
3675////////////////////////////////////////////////////////////////////////////////
3676/// Operator /
3677
3679{
3680 TH3C hnew = h1;
3681 hnew.Divide(&h2);
3682 hnew.SetDirectory(0);
3683 return hnew;
3684}
3685
3686
3687//______________________________________________________________________________
3688// TH3S methods
3689// TH3S a 3-D histogram with two bytes per cell (short integer)
3690//______________________________________________________________________________
3691
3692ClassImp(TH3S);
3693
3694
3695////////////////////////////////////////////////////////////////////////////////
3696/// Constructor.
3697
3699{
3700 SetBinsLength(27);
3701 if (fgDefaultSumw2) Sumw2();
3702}
3703
3704
3705////////////////////////////////////////////////////////////////////////////////
3706/// Destructor.
3707
3709{
3710}
3711
3712
3713////////////////////////////////////////////////////////////////////////////////
3714/// Constructor for fix bin size 3-D histograms.
3715/// (see TH3::TH3 for explanation of parameters)
3716
3717TH3S::TH3S(const char *name,const char *title,Int_t nbinsx,Double_t xlow,Double_t xup
3718 ,Int_t nbinsy,Double_t ylow,Double_t yup
3719 ,Int_t nbinsz,Double_t zlow,Double_t zup)
3720 :TH3(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup,nbinsz,zlow,zup)
3721{
3723 if (fgDefaultSumw2) Sumw2();
3724
3725 if (xlow >= xup || ylow >= yup || zlow >= zup) SetBuffer(fgBufferSize);
3726}
3727
3728
3729////////////////////////////////////////////////////////////////////////////////
3730/// Constructor for variable bin size 3-D histograms.
3731/// (see TH3::TH3 for explanation of parameters)
3732
3733TH3S::TH3S(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins
3734 ,Int_t nbinsy,const Float_t *ybins
3735 ,Int_t nbinsz,const Float_t *zbins)
3736 :TH3(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins)
3737{
3739 if (fgDefaultSumw2) Sumw2();
3740}
3741
3742
3743////////////////////////////////////////////////////////////////////////////////
3744/// Constructor for variable bin size 3-D histograms.
3745/// (see TH3::TH3 for explanation of parameters)
3746
3747TH3S::TH3S(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins
3748 ,Int_t nbinsy,const Double_t *ybins
3749 ,Int_t nbinsz,const Double_t *zbins)
3750 :TH3(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins)
3751{
3753 if (fgDefaultSumw2) Sumw2();
3754}
3755
3756
3757////////////////////////////////////////////////////////////////////////////////
3758/// Copy Constructor.
3759/// The list of functions is not copied. (Use Clone() if needed)
3760
3761TH3S::TH3S(const TH3S &h3s) : TH3(), TArrayS()
3762{
3763 ((TH3S&)h3s).Copy(*this);
3764}
3765
3766
3767////////////////////////////////////////////////////////////////////////////////
3768/// Increment bin content by 1.
3769
3771{
3772 if (fArray[bin] < 32767) fArray[bin]++;
3773}
3774
3775
3776////////////////////////////////////////////////////////////////////////////////
3777/// Increment bin content by w.
3778
3780{
3781 Int_t newval = fArray[bin] + Int_t(w);
3782 if (newval > -32768 && newval < 32768) {fArray[bin] = Short_t(newval); return;}
3783 if (newval < -32767) fArray[bin] = -32767;
3784 if (newval > 32767) fArray[bin] = 32767;
3785}
3786
3787
3788////////////////////////////////////////////////////////////////////////////////
3789/// Copy this 3-D histogram structure to newth3.
3790
3791void TH3S::Copy(TObject &newth3) const
3792{
3793 TH3::Copy((TH3S&)newth3);
3794}
3795
3796
3797////////////////////////////////////////////////////////////////////////////////
3798/// Reset this histogram: contents, errors, etc.
3799
3801{
3802 TH3::Reset(option);
3804 // should also reset statistics once statistics are implemented for TH3
3805}
3806
3807
3808////////////////////////////////////////////////////////////////////////////////
3809/// Set total number of bins including under/overflow
3810/// Reallocate bin contents array
3811
3813{
3814 if (n < 0) n = (fXaxis.GetNbins()+2)*(fYaxis.GetNbins()+2)*(fZaxis.GetNbins()+2);
3815 fNcells = n;
3816 TArrayS::Set(n);
3817}
3818
3819
3820////////////////////////////////////////////////////////////////////////////////
3821/// Stream an object of class TH3S.
3822
3823void TH3S::Streamer(TBuffer &R__b)
3824{
3825 if (R__b.IsReading()) {
3826 UInt_t R__s, R__c;
3827 if (R__b.GetParent() && R__b.GetVersionOwner() < 22300) return;
3828 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
3829 if (R__v > 2) {
3830 R__b.ReadClassBuffer(TH3S::Class(), this, R__v, R__s, R__c);
3831 return;
3832 }
3833 //====process old versions before automatic schema evolution
3834 if (R__v < 2) {
3835 R__b.ReadVersion();
3836 TH1::Streamer(R__b);
3837 TArrayS::Streamer(R__b);
3838 R__b.ReadVersion(&R__s, &R__c);
3839 TAtt3D::Streamer(R__b);
3840 } else {
3841 TH3::Streamer(R__b);
3842 TArrayS::Streamer(R__b);
3843 R__b.CheckByteCount(R__s, R__c, TH3S::IsA());
3844 }
3845 //====end of old versions
3846
3847 } else {
3848 R__b.WriteClassBuffer(TH3S::Class(),this);
3849 }
3850}
3851
3852
3853////////////////////////////////////////////////////////////////////////////////
3854/// Operator =
3855
3857{
3858 if (this != &h1) ((TH3S&)h1).Copy(*this);
3859 return *this;
3860}
3861
3862
3863////////////////////////////////////////////////////////////////////////////////
3864/// Operator *
3865
3867{
3868 TH3S hnew = h1;
3869 hnew.Scale(c1);
3870 hnew.SetDirectory(0);
3871 return hnew;
3872}
3873
3874
3875////////////////////////////////////////////////////////////////////////////////
3876/// Operator +
3877
3879{
3880 TH3S hnew = h1;
3881 hnew.Add(&h2,1);
3882 hnew.SetDirectory(0);
3883 return hnew;
3884}
3885
3886
3887////////////////////////////////////////////////////////////////////////////////
3888/// Operator -
3889
3891{
3892 TH3S hnew = h1;
3893 hnew.Add(&h2,-1);
3894 hnew.SetDirectory(0);
3895 return hnew;
3896}
3897
3898
3899////////////////////////////////////////////////////////////////////////////////
3900/// Operator *
3901
3903{
3904 TH3S hnew = h1;
3905 hnew.Multiply(&h2);
3906 hnew.SetDirectory(0);
3907 return hnew;
3908}
3909
3910
3911////////////////////////////////////////////////////////////////////////////////
3912/// Operator /
3913
3915{
3916 TH3S hnew = h1;
3917 hnew.Divide(&h2);
3918 hnew.SetDirectory(0);
3919 return hnew;
3920}
3921
3922
3923//______________________________________________________________________________
3924// TH3I methods
3925// TH3I a 3-D histogram with four bytes per cell (32 bits integer)
3926//______________________________________________________________________________
3927
3928ClassImp(TH3I);
3929
3930
3931////////////////////////////////////////////////////////////////////////////////
3932/// Constructor.
3933
3935{
3936 SetBinsLength(27);
3937 if (fgDefaultSumw2) Sumw2();
3938}
3939
3940
3941////////////////////////////////////////////////////////////////////////////////
3942/// Destructor.
3943
3945{
3946}
3947
3948
3949////////////////////////////////////////////////////////////////////////////////
3950/// Constructor for fix bin size 3-D histograms
3951/// (see TH3::TH3 for explanation of parameters)
3952
3953TH3I::TH3I(const char *name,const char *title,Int_t nbinsx,Double_t xlow,Double_t xup
3954 ,Int_t nbinsy,Double_t ylow,Double_t yup
3955 ,Int_t nbinsz,Double_t zlow,Double_t zup)
3956 :TH3(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup,nbinsz,zlow,zup)
3957{
3959 if (fgDefaultSumw2) Sumw2();
3960
3961 if (xlow >= xup || ylow >= yup || zlow >= zup) SetBuffer(fgBufferSize);
3962}
3963
3964
3965////////////////////////////////////////////////////////////////////////////////
3966/// Constructor for variable bin size 3-D histograms
3967/// (see TH3::TH3 for explanation of parameters)
3968
3969TH3I::TH3I(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins
3970 ,Int_t nbinsy,const Float_t *ybins
3971 ,Int_t nbinsz,const Float_t *zbins)
3972 :TH3(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins)
3973{
3975 if (fgDefaultSumw2) Sumw2();
3976}
3977
3978
3979////////////////////////////////////////////////////////////////////////////////
3980/// Constructor for variable bin size 3-D histograms
3981/// (see TH3::TH3 for explanation of parameters)
3982
3983TH3I::TH3I(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins
3984 ,Int_t nbinsy,const Double_t *ybins
3985 ,Int_t nbinsz,const Double_t *zbins)
3986 :TH3(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins)
3987{
3989 if (fgDefaultSumw2) Sumw2();
3990}
3991
3992
3993////////////////////////////////////////////////////////////////////////////////
3994/// Copy constructor.
3995/// The list of functions is not copied. (Use Clone() if needed)
3996
3997TH3I::TH3I(const TH3I &h3i) : TH3(), TArrayI()
3998{
3999 ((TH3I&)h3i).Copy(*this);
4000}
4001
4002
4003////////////////////////////////////////////////////////////////////////////////
4004/// Increment bin content by 1.
4005
4007{
4008 if (fArray[bin] < INT_MAX) fArray[bin]++;
4009}
4010
4011
4012////////////////////////////////////////////////////////////////////////////////
4013/// Increment bin content by w.
4014
4016{
4017 Long64_t newval = fArray[bin] + Long64_t(w);
4018 if (newval > -INT_MAX && newval < INT_MAX) {fArray[bin] = Int_t(newval); return;}
4019 if (newval < -INT_MAX) fArray[bin] = -INT_MAX;
4020 if (newval > INT_MAX) fArray[bin] = INT_MAX;
4021}
4022
4023
4024////////////////////////////////////////////////////////////////////////////////
4025/// Copy this 3-D histogram structure to newth3.
4026
4027void TH3I::Copy(TObject &newth3) const
4028{
4029 TH3::Copy((TH3I&)newth3);
4030}
4031
4032
4033////////////////////////////////////////////////////////////////////////////////
4034/// Reset this histogram: contents, errors, etc.
4035
4037{
4038 TH3::Reset(option);
4040 // should also reset statistics once statistics are implemented for TH3
4041}
4042
4043
4044////////////////////////////////////////////////////////////////////////////////
4045/// Set total number of bins including under/overflow
4046/// Reallocate bin contents array
4047
4049{
4050 if (n < 0) n = (fXaxis.GetNbins()+2)*(fYaxis.GetNbins()+2)*(fZaxis.GetNbins()+2);
4051 fNcells = n;
4052 TArrayI::Set(n);
4053}
4054
4055
4056////////////////////////////////////////////////////////////////////////////////
4057/// Operator =
4058
4060{
4061 if (this != &h1) ((TH3I&)h1).Copy(*this);
4062 return *this;
4063}
4064
4065
4066////////////////////////////////////////////////////////////////////////////////
4067/// Operator *
4068
4070{
4071 TH3I hnew = h1;
4072 hnew.Scale(c1);
4073 hnew.SetDirectory(0);
4074 return hnew;
4075}
4076
4077
4078////////////////////////////////////////////////////////////////////////////////
4079/// Operator +
4080
4082{
4083 TH3I hnew = h1;
4084 hnew.Add(&h2,1);
4085 hnew.SetDirectory(0);
4086 return hnew;
4087}
4088
4089
4090////////////////////////////////////////////////////////////////////////////////
4091/// Operator _
4092
4094{
4095 TH3I hnew = h1;
4096 hnew.Add(&h2,-1);
4097 hnew.SetDirectory(0);
4098 return hnew;
4099}
4100
4101
4102////////////////////////////////////////////////////////////////////////////////
4103/// Operator *
4104
4106{
4107 TH3I hnew = h1;
4108 hnew.Multiply(&h2);
4109 hnew.SetDirectory(0);
4110 return hnew;
4111}
4112
4113
4114////////////////////////////////////////////////////////////////////////////////
4115/// Operator /
4116
4118{
4119 TH3I hnew = h1;
4120 hnew.Divide(&h2);
4121 hnew.SetDirectory(0);
4122 return hnew;
4123}
4124
4125
4126//______________________________________________________________________________
4127// TH3F methods
4128// TH3F a 3-D histogram with four bytes per cell (float)
4129//______________________________________________________________________________
4130
4131ClassImp(TH3F);
4132
4133
4134////////////////////////////////////////////////////////////////////////////////
4135/// Constructor.
4136
4138{
4139 SetBinsLength(27);
4140 if (fgDefaultSumw2) Sumw2();
4141}
4142
4143
4144////////////////////////////////////////////////////////////////////////////////
4145/// Destructor.
4146
4148{
4149}
4150
4151
4152////////////////////////////////////////////////////////////////////////////////
4153/// Constructor for fix bin size 3-D histograms
4154/// (see TH3::TH3 for explanation of parameters)
4155
4156TH3F::TH3F(const char *name,const char *title,Int_t nbinsx,Double_t xlow,Double_t xup
4157 ,Int_t nbinsy,Double_t ylow,Double_t yup
4158 ,Int_t nbinsz,Double_t zlow,Double_t zup)
4159 :TH3(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup,nbinsz,zlow,zup)
4160{
4162 if (fgDefaultSumw2) Sumw2();
4163
4164 if (xlow >= xup || ylow >= yup || zlow >= zup) SetBuffer(fgBufferSize);
4165}
4166
4167
4168////////////////////////////////////////////////////////////////////////////////
4169/// Constructor for variable bin size 3-D histograms
4170/// (see TH3::TH3 for explanation of parameters)
4171
4172TH3F::TH3F(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins
4173 ,Int_t nbinsy,const Float_t *ybins
4174 ,Int_t nbinsz,const Float_t *zbins)
4175 :TH3(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins)
4176{
4178 if (fgDefaultSumw2) Sumw2();
4179}
4180
4181
4182////////////////////////////////////////////////////////////////////////////////
4183/// Constructor for variable bin size 3-D histograms
4184/// (see TH3::TH3 for explanation of parameters)
4185
4186TH3F::TH3F(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins
4187 ,Int_t nbinsy,const Double_t *ybins
4188 ,Int_t nbinsz,const Double_t *zbins)
4189 :TH3(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins)
4190{
4192 if (fgDefaultSumw2) Sumw2();
4193}
4194
4195
4196////////////////////////////////////////////////////////////////////////////////
4197/// Copy constructor.
4198/// The list of functions is not copied. (Use Clone() if needed)
4199
4200TH3F::TH3F(const TH3F &h3f) : TH3(), TArrayF()
4201{
4202 ((TH3F&)h3f).Copy(*this);
4203}
4204
4205
4206////////////////////////////////////////////////////////////////////////////////
4207/// Copy this 3-D histogram structure to newth3.
4208
4209void TH3F::Copy(TObject &newth3) const
4210{
4211 TH3::Copy((TH3F&)newth3);
4212}
4213
4214
4215////////////////////////////////////////////////////////////////////////////////
4216/// Reset this histogram: contents, errors, etc.
4217
4219{
4220 TH3::Reset(option);
4222 // should also reset statistics once statistics are implemented for TH3
4223}
4224
4225
4226////////////////////////////////////////////////////////////////////////////////
4227/// Set total number of bins including under/overflow
4228/// Reallocate bin contents array
4229
4231{
4232 if (n < 0) n = (fXaxis.GetNbins()+2)*(fYaxis.GetNbins()+2)*(fZaxis.GetNbins()+2);
4233 fNcells = n;
4234 TArrayF::Set(n);
4235}
4236
4237
4238////////////////////////////////////////////////////////////////////////////////
4239/// Stream an object of class TH3F.
4240
4241void TH3F::Streamer(TBuffer &R__b)
4242{
4243 if (R__b.IsReading()) {
4244 UInt_t R__s, R__c;
4245 if (R__b.GetParent() && R__b.GetVersionOwner() < 22300) return;
4246 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
4247 if (R__v > 2) {
4248 R__b.ReadClassBuffer(TH3F::Class(), this, R__v, R__s, R__c);
4249 return;
4250 }
4251 //====process old versions before automatic schema evolution
4252 if (R__v < 2) {
4253 R__b.ReadVersion();
4254 TH1::Streamer(R__b);
4255 TArrayF::Streamer(R__b);
4256 R__b.ReadVersion(&R__s, &R__c);
4257 TAtt3D::Streamer(R__b);
4258 } else {
4259 TH3::Streamer(R__b);
4260 TArrayF::Streamer(R__b);
4261 R__b.CheckByteCount(R__s, R__c, TH3F::IsA());
4262 }
4263 //====end of old versions
4264
4265 } else {
4266 R__b.WriteClassBuffer(TH3F::Class(),this);
4267 }
4268}
4269
4270
4271////////////////////////////////////////////////////////////////////////////////
4272/// Operator =
4273
4275{
4276 if (this != &h1) ((TH3F&)h1).Copy(*this);
4277 return *this;
4278}
4279
4280
4281////////////////////////////////////////////////////////////////////////////////
4282/// Operator *
4283
4285{
4286 TH3F hnew = h1;
4287 hnew.Scale(c1);
4288 hnew.SetDirectory(0);
4289 return hnew;
4290}
4291
4292
4293////////////////////////////////////////////////////////////////////////////////
4294/// Operator +
4295
4297{
4298 TH3F hnew = h1;
4299 hnew.Add(&h2,1);
4300 hnew.SetDirectory(0);
4301 return hnew;
4302}
4303
4304
4305////////////////////////////////////////////////////////////////////////////////
4306/// Operator -
4307
4309{
4310 TH3F hnew = h1;
4311 hnew.Add(&h2,-1);
4312 hnew.SetDirectory(0);
4313 return hnew;
4314}
4315
4316
4317////////////////////////////////////////////////////////////////////////////////
4318/// Operator *
4319
4321{
4322 TH3F hnew = h1;
4323 hnew.Multiply(&h2);
4324 hnew.SetDirectory(0);
4325 return hnew;
4326}
4327
4328
4329////////////////////////////////////////////////////////////////////////////////
4330/// Operator /
4331
4333{
4334 TH3F hnew = h1;
4335 hnew.Divide(&h2);
4336 hnew.SetDirectory(0);
4337 return hnew;
4338}
4339
4340
4341//______________________________________________________________________________
4342// TH3D methods
4343// TH3D a 3-D histogram with eight bytes per cell (double)
4344//______________________________________________________________________________
4345
4346ClassImp(TH3D);
4347
4348
4349////////////////////////////////////////////////////////////////////////////////
4350/// Constructor.
4351
4353{
4354 SetBinsLength(27);
4355 if (fgDefaultSumw2) Sumw2();
4356}
4357
4358
4359////////////////////////////////////////////////////////////////////////////////
4360/// Destructor.
4361
4363{
4364}
4365
4366
4367////////////////////////////////////////////////////////////////////////////////
4368/// Constructor for fix bin size 3-D histograms
4369/// (see TH3::TH3 for explanation of parameters)
4370
4371TH3D::TH3D(const char *name,const char *title,Int_t nbinsx,Double_t xlow,Double_t xup
4372 ,Int_t nbinsy,Double_t ylow,Double_t yup
4373 ,Int_t nbinsz,Double_t zlow,Double_t zup)
4374 :TH3(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup,nbinsz,zlow,zup)
4375{
4377 if (fgDefaultSumw2) Sumw2();
4378
4379 if (xlow >= xup || ylow >= yup || zlow >= zup) SetBuffer(fgBufferSize);
4380}
4381
4382
4383////////////////////////////////////////////////////////////////////////////////
4384/// Constructor for variable bin size 3-D histograms
4385/// (see TH3::TH3 for explanation of parameters)
4386
4387TH3D::TH3D(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins
4388 ,Int_t nbinsy,const Float_t *ybins
4389 ,Int_t nbinsz,const Float_t *zbins)
4390 :TH3(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins)
4391{
4393 if (fgDefaultSumw2) Sumw2();
4394}
4395
4396
4397////////////////////////////////////////////////////////////////////////////////
4398/// Constructor for variable bin size 3-D histograms
4399/// (see TH3::TH3 for explanation of parameters)
4400
4401TH3D::TH3D(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins
4402 ,Int_t nbinsy,const Double_t *ybins
4403 ,Int_t nbinsz,const Double_t *zbins)
4404 :TH3(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins)
4405{
4407 if (fgDefaultSumw2) Sumw2();
4408}
4409
4410
4411////////////////////////////////////////////////////////////////////////////////
4412/// Copy constructor.
4413/// The list of functions is not copied. (Use Clone() if needed)
4414
4415TH3D::TH3D(const TH3D &h3d) : TH3(), TArrayD()
4416{
4417 ((TH3D&)h3d).Copy(*this);
4418}
4419
4420
4421////////////////////////////////////////////////////////////////////////////////
4422/// Copy this 3-D histogram structure to newth3.
4423
4424void TH3D::Copy(TObject &newth3) const
4425{
4426 TH3::Copy((TH3D&)newth3);
4427}
4428
4429
4430////////////////////////////////////////////////////////////////////////////////
4431/// Reset this histogram: contents, errors, etc.
4432
4434{
4435 TH3::Reset(option);
4437 // should also reset statistics once statistics are implemented for TH3
4438}
4439
4440
4441////////////////////////////////////////////////////////////////////////////////
4442/// Set total number of bins including under/overflow
4443/// Reallocate bin contents array
4444
4446{
4447 if (n < 0) n = (fXaxis.GetNbins()+2)*(fYaxis.GetNbins()+2)*(fZaxis.GetNbins()+2);
4448 fNcells = n;
4449 TArrayD::Set(n);
4450}
4451
4452
4453////////////////////////////////////////////////////////////////////////////////
4454/// Stream an object of class TH3D.
4455
4456void TH3D::Streamer(TBuffer &R__b)
4457{
4458 if (R__b.IsReading()) {
4459 UInt_t R__s, R__c;
4460 if (R__b.GetParent() && R__b.GetVersionOwner() < 22300) return;
4461 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
4462 if (R__v > 2) {
4463 R__b.ReadClassBuffer(TH3D::Class(), this, R__v, R__s, R__c);
4464 return;
4465 }
4466 //====process old versions before automatic schema evolution
4467 if (R__v < 2) {
4468 R__b.ReadVersion();
4469 TH1::Streamer(R__b);
4470 TArrayD::Streamer(R__b);
4471 R__b.ReadVersion(&R__s, &R__c);
4472 TAtt3D::Streamer(R__b);
4473 } else {
4474 TH3::Streamer(R__b);
4475 TArrayD::Streamer(R__b);
4476 R__b.CheckByteCount(R__s, R__c, TH3D::IsA());
4477 }
4478 //====end of old versions
4479
4480 } else {
4481 R__b.WriteClassBuffer(TH3D::Class(),this);
4482 }
4483}
4484
4485
4486////////////////////////////////////////////////////////////////////////////////
4487/// Operator =
4488
4490{
4491 if (this != &h1) ((TH3D&)h1).Copy(*this);
4492 return *this;
4493}
4494
4495
4496////////////////////////////////////////////////////////////////////////////////
4497/// Operator *
4498
4500{
4501 TH3D hnew = h1;
4502 hnew.Scale(c1);
4503 hnew.SetDirectory(0);
4504 return hnew;
4505}
4506
4507
4508////////////////////////////////////////////////////////////////////////////////
4509/// Operator +
4510
4512{
4513 TH3D hnew = h1;
4514 hnew.Add(&h2,1);
4515 hnew.SetDirectory(0);
4516 return hnew;
4517}
4518
4519
4520////////////////////////////////////////////////////////////////////////////////
4521/// Operator -
4522
4524{
4525 TH3D hnew = h1;
4526 hnew.Add(&h2,-1);
4527 hnew.SetDirectory(0);
4528 return hnew;
4529}
4530
4531
4532////////////////////////////////////////////////////////////////////////////////
4533/// Operator *
4534
4536{
4537 TH3D hnew = h1;
4538 hnew.Multiply(&h2);
4539 hnew.SetDirectory(0);
4540 return hnew;
4541}
4542
4543
4544////////////////////////////////////////////////////////////////////////////////
4545/// Operator /
4546
4548{
4549 TH3D hnew = h1;
4550 hnew.Divide(&h2);
4551 hnew.SetDirectory(0);
4552 return hnew;
4553}
void Class()
Definition: Class.C:29
#define c(i)
Definition: RSha256.hxx:101
#define s1(x)
Definition: RSha256.hxx:91
#define h(i)
Definition: RSha256.hxx:106
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:45
short Version_t
Definition: RtypesCore.h:65
char Char_t
Definition: RtypesCore.h:33
unsigned int UInt_t
Definition: RtypesCore.h:46
const Bool_t kFALSE
Definition: RtypesCore.h:101
bool Bool_t
Definition: RtypesCore.h:63
short Short_t
Definition: RtypesCore.h:39
double Double_t
Definition: RtypesCore.h:59
short Color_t
Definition: RtypesCore.h:92
long long Long64_t
Definition: RtypesCore.h:80
short Style_t
Definition: RtypesCore.h:89
float Float_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:100
const char Option_t
Definition: RtypesCore.h:66
#define ClassImp(name)
Definition: Rtypes.h:364
#define R__ASSERT(e)
Definition: TError.h:118
char name[80]
Definition: TGX11.cxx:110
TH3C operator/(TH3C &h1, TH3C &h2)
Operator /.
Definition: TH3.cxx:3678
TH3C operator*(Float_t c1, TH3C &h1)
Operator *.
Definition: TH3.cxx:3630
TH3C operator-(TH3C &h1, TH3C &h2)
Operator -.
Definition: TH3.cxx:3654
TH3C operator+(TH3C &h1, TH3C &h2)
Operator +.
Definition: TH3.cxx:3642
float xmin
Definition: THbookFile.cxx:95
int nentries
Definition: THbookFile.cxx:91
float ymin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
float ymax
Definition: THbookFile.cxx:95
#define gROOT
Definition: TROOT.h:404
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
#define gPad
Definition: TVirtualPad.h:287
Array of chars or bytes (8 bits per element).
Definition: TArrayC.h:27
void Set(Int_t n)
Set size of this array to n chars.
Definition: TArrayC.cxx:105
Char_t * fArray
Definition: TArrayC.h:30
void Copy(TArrayC &array) const
Definition: TArrayC.h:42
void Reset(Char_t val=0)
Definition: TArrayC.h:47
Array of doubles (64 bits per element).
Definition: TArrayD.h:27
Double_t * fArray
Definition: TArrayD.h:30
void Copy(TArrayD &array) const
Definition: TArrayD.h:42
void Set(Int_t n)
Set size of this array to n doubles.
Definition: TArrayD.cxx:106
void Reset()
Definition: TArrayD.h:47
Array of floats (32 bits per element).
Definition: TArrayF.h:27
void Copy(TArrayF &array) const
Definition: TArrayF.h:42
void Reset()
Definition: TArrayF.h:47
void Set(Int_t n)
Set size of this array to n floats.
Definition: TArrayF.cxx:105
Array of integers (32 bits per element).
Definition: TArrayI.h:27
Int_t * fArray
Definition: TArrayI.h:30
void Set(Int_t n)
Set size of this array to n ints.
Definition: TArrayI.cxx:105
void Reset()
Definition: TArrayI.h:47
void Copy(TArrayI &array) const
Definition: TArrayI.h:42
Array of shorts (16 bits per element).
Definition: TArrayS.h:27
void Set(Int_t n)
Set size of this array to n shorts.
Definition: TArrayS.cxx:105
void Reset()
Definition: TArrayS.h:47
void Copy(TArrayS &array) const
Definition: TArrayS.h:42
Short_t * fArray
Definition: TArrayS.h:30
Int_t fN
Definition: TArray.h:38
Int_t GetSize() const
Definition: TArray.h:47
Use this attribute class when an object should have 3D capabilities.
Definition: TAtt3D.h:19
virtual Color_t GetTitleColor() const
Definition: TAttAxis.h:46
virtual Color_t GetLabelColor() const
Definition: