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