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