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