Logo ROOT  
Reference Guide
TProfileHelper.h
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2// Author: David Gonzalez Maline 18/01/2008
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#ifndef ROOT_TProfileHelper
13#define ROOT_TProfileHelper
14
15
16//////////////////////////////////////////////////////////////////////////
17// //
18// TProfileHelper //
19// //
20// Profile helper class. //
21// //
22//////////////////////////////////////////////////////////////////////////
23
24#include "TH1.h"
25#include "TError.h"
26#include "THashList.h"
27#include "TMath.h"
28#include "TH1Merger.h"
29
31
32public:
33 template <typename T>
34 static Bool_t Add(T* p, const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2=1);
35
36 template <typename T>
37 static void BuildArray(T* p);
38
39 template <typename T>
40 static Double_t GetBinEffectiveEntries(T* p, Int_t bin);
41
42 template <typename T>
43 static Long64_t Merge(T* p, TCollection *list);
44
45 template <typename T>
46 static T* ExtendAxis(T* p, Double_t x, TAxis *axis);
47
48 template <typename T>
49 static void Scale(T* p, Double_t c1, Option_t * option);
50
51 template <typename T>
52 static void Sumw2(T* p, Bool_t flag );
53
54 template <typename T>
55 static void LabelsDeflate(T* p, Option_t *);
56
57 template <typename T>
58 static void LabelsInflate(T* p, Option_t *);
59
60 template <typename T>
61 static Double_t GetBinError(T* p, Int_t bin);
62
63 template <typename T>
64 static void SetBinEntries(T* p, Int_t bin, Double_t w);
65
66 template <typename T>
67 static void SetErrorOption(T* p, Option_t * opt);
68};
69
70template <typename T>
72{
73 // Performs the operation: this = c1*h1 + c2*h2
74
75 T *p1 = (T*)h1;
76 T *p2 = (T*)h2;
77
78 // delete buffer if it is there since it will become invalid
79 if (p->fBuffer) p->BufferEmpty(1);
80
81// Check profile compatibility
82 Int_t nx = p->GetNbinsX();
83 Int_t ny = p->GetNbinsY();
84 Int_t nz = p->GetNbinsZ();
85
86 if ( nx != p1->GetNbinsX() || nx != p2->GetNbinsX() ||
87 ny != p1->GetNbinsY() || ny != p2->GetNbinsY() ||
88 nz != p1->GetNbinsZ() || nz != p2->GetNbinsZ() ) {
89 Error("TProfileHelper::Add","Attempt to add profiles with different number of bins");
90 return kFALSE;
91 }
92
93// Add statistics
94 Double_t ac1 = TMath::Abs(c1);
95 Double_t ac2 = TMath::Abs(c2);
96 p->fEntries = ac1*p1->GetEntries() + ac2*p2->GetEntries();
98 Int_t i;
99 for (i=0;i<TH1::kNstat;i++) {s0[i] = s1[i] = s2[i] = 0;}
100 p->GetStats(s0);
101 p1->GetStats(s1);
102 p2->GetStats(s2);
103 for (i=0;i<TH1::kNstat;i++) {
104 if (i == 1) s0[i] = c1*c1*s1[i] + c2*c2*s2[i];
105 else s0[i] = ac1*s1[i] + ac2*s2[i];
106 }
107 p->PutStats(s0);
108
109 // create sumw2 per bin if not set
110 if (p->fBinSumw2.fN == 0 && (p1->fBinSumw2.fN != 0 || p2->fBinSumw2.fN != 0)) p->Sumw2();
111
112 // Make the loop over the bins to calculate the Addition
113 Double_t *cu1 = p1->GetW(); Double_t *cu2 = p2->GetW();
114 Double_t *er1 = p1->GetW2(); Double_t *er2 = p2->GetW2();
115 Double_t *en1 = p1->GetB(); Double_t *en2 = p2->GetB();
116 Double_t *ew1 = p1->GetB2(); Double_t *ew2 = p2->GetB2();
117 // if p1 has not the sum of weight squared/bin stored use just the sum of weights
118 if (ew1 == nullptr) ew1 = en1;
119 if (ew2 == nullptr) ew2 = en2;
120 for (Int_t bin = 0; bin < p->fN; bin++) {
121 p->fArray[bin] = c1*cu1[bin] + c2*cu2[bin];
122 p->fSumw2.fArray[bin] = ac1*er1[bin] + ac2*er2[bin];
123 p->fBinEntries.fArray[bin] = ac1*en1[bin] + ac2*en2[bin];
124 if (p->fBinSumw2.fN ) p->fBinSumw2.fArray[bin] = ac1*ac1*ew1[bin] + ac2*ac2*ew2[bin];
125 }
126 return kTRUE;
127}
128
129template <typename T>
131 // Build the extra profile data structure in addition to the histograms
132 // this are: array of bin entries: fBinEntries
133 // array of sum of profiled observable value - squared
134 // stored in TH1::fSumw2
135 // array of some of weight squared (optional) in TProfile::fBinSumw2
136 p->fBinEntries.Set(p->fNcells);
137 p->fSumw2.Set(p->fNcells);
138 if (TH1::GetDefaultSumw2() || p->fBinSumw2.fN > 0 ) p->fBinSumw2.Set(p->fNcells);
139}
140
141
142template <typename T>
144{
145// Return bin effective entries for a weighted filled Profile histogram.
146// In case of an unweighted profile, it is equivalent to the number of entries per bin
147// The effective entries is defined as the square of the sum of the weights divided by the
148// sum of the weights square.
149// TProfile::Sumw2() must be called before filling the profile with weights.
150// Only by calling this method the sum of the square of the weights per bin is stored.
151//
152
153 if (p->fBuffer) p->BufferEmpty();
154
155 if (bin < 0 || bin >= p->fNcells) return 0;
156 double sumOfWeights = p->fBinEntries.fArray[bin];
157 if ( p->fBinSumw2.fN == 0 || p->fBinSumw2.fN != p->fNcells) {
158 // this can happen when reading an old file
159 p->fBinSumw2.Set(0);
160 return sumOfWeights;
161 }
162 double sumOfWeightsSquare = p->fBinSumw2.fArray[bin];
163 return ( sumOfWeightsSquare > 0 ? sumOfWeights * sumOfWeights / sumOfWeightsSquare : 0 );
164}
165
166template <typename T>
168 //Merge all histograms in the collection in this histogram.
169 //This function computes the min/max for the axes,
170 //compute a new number of bins, if necessary,
171 //add bin contents, errors and statistics.
172 //If overflows are present and limits are different the function will fail.
173 //The function returns the total number of entries in the result histogram
174 //if the merge is successful, -1 otherwise.
175 //
176 //IMPORTANT remark. The 2 axis x and y may have different number
177 //of bins and different limits, BUT the largest bin width must be
178 //a multiple of the smallest bin width and the upper limit must also
179 //be a multiple of the bin width.
180
181 if (!li) return 0;
182 if (li->IsEmpty()) return (Int_t) p->GetEntries();
183
184 TList inlist;
185 inlist.AddAll(li);
186
187 // use TH1Merger class
188 TH1Merger merger(*p, *li, "");
189 Bool_t ret = merger();
190
191 return (ret) ? p->GetEntries() : -1;
192
193#ifdef OLD_PROFILE_MERGE
194
195 TAxis newXAxis;
196 TAxis newYAxis;
197 TAxis newZAxis;
198 Bool_t initialLimitsFound = kFALSE;
199 Bool_t allSameLimits = kTRUE;
200 Bool_t allHaveLimits = kTRUE;
201// Bool_t firstNonEmptyHist = kTRUE;
202
203 TIter next(&inlist);
204 T* h = p;
205
206 do {
207
208 // skip empty histograms
209 // if (h->fTsumw == 0 && h->GetEntries() == 0) continue;
210
211 Bool_t hasLimits = h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax();
212 allHaveLimits = allHaveLimits && hasLimits;
213
214 if (hasLimits) {
215 h->BufferEmpty();
216
217#ifdef LATER
218 // this is done in case the first histograms are empty and
219 // the histogram have different limits
220 if (firstNonEmptyHist ) {
221 // set axis limits in the case the first histogram was empty
222 if (h != p ) {
223 if (!p->SameLimitsAndNBins(p->fXaxis, *(h->GetXaxis())) )
224 p->fXaxis.Set(h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax());
225 if (!p->SameLimitsAndNBins(p->fYaxis, *(h->GetYaxis())) )
226 p->fYaxis.Set(h->GetYaxis()->GetNbins(), h->GetYaxis()->GetXmin(),h->GetYaxis()->GetXmax());
227 if (!p->SameLimitsAndNBins(p->fZaxis, *(h->GetZaxis())) )
228 p->fZaxis.Set(h->GetZaxis()->GetNbins(), h->GetZaxis()->GetXmin(),h->GetZaxis()->GetXmax());
229 }
230 firstNonEmptyHist = kFALSE;
231 }
232#endif
233
234 // this is executed the first time an histogram with limits is found
235 // to set some initial values on the new axis
236 if (!initialLimitsFound) {
237 initialLimitsFound = kTRUE;
238 newXAxis.Set(h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),
239 h->GetXaxis()->GetXmax());
240 if ( p->GetDimension() >= 2 )
241 newYAxis.Set(h->GetYaxis()->GetNbins(), h->GetYaxis()->GetXmin(),
242 h->GetYaxis()->GetXmax());
243 if ( p->GetDimension() >= 3 )
244 newZAxis.Set(h->GetZaxis()->GetNbins(), h->GetZaxis()->GetXmin(),
245 h->GetZaxis()->GetXmax());
246 }
247 else {
248 // check first if histograms have same bins
249 if (!p->SameLimitsAndNBins(newXAxis, *(h->GetXaxis())) ||
250 !p->SameLimitsAndNBins(newYAxis, *(h->GetYaxis())) ||
251 !p->SameLimitsAndNBins(newZAxis, *(h->GetZaxis())) ) {
252
253 allSameLimits = kFALSE;
254 // recompute in this case the optimal limits
255 // The condition to works is that the histogram have same bin with
256 // and one common bin edge
257
258 if (!p->RecomputeAxisLimits(newXAxis, *(h->GetXaxis()))) {
259 Error("TProfileHelper::Merge", "Cannot merge profiles %d dim - limits are inconsistent:\n "
260 "first: (%d, %f, %f), second: (%d, %f, %f)",p->GetDimension(),
261 newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax(),
262 h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),
263 h->GetXaxis()->GetXmax());
264 return -1;
265 }
266 if (p->GetDimension() >= 2 && !p->RecomputeAxisLimits(newYAxis, *(h->GetYaxis()))) {
267 Error("TProfileHelper::Merge", "Cannot merge profiles %d dim - limits are inconsistent:\n "
268 "first: (%d, %f, %f), second: (%d, %f, %f)",p->GetDimension(),
269 newYAxis.GetNbins(), newYAxis.GetXmin(), newYAxis.GetXmax(),
270 h->GetYaxis()->GetNbins(), h->GetYaxis()->GetXmin(),
271 h->GetYaxis()->GetXmax());
272 return -1;
273 }
274 if (p->GetDimension() >= 3 && !p->RecomputeAxisLimits(newZAxis, *(h->GetZaxis()))) {
275 Error("TProfileHelper::Merge", "Cannot merge profiles %d dim - limits are inconsistent:\n "
276 "first: (%d, %f, %f), second: (%d, %f, %f)",p->GetDimension(),
277 newZAxis.GetNbins(), newZAxis.GetXmin(), newZAxis.GetXmax(),
278 h->GetZaxis()->GetNbins(), h->GetZaxis()->GetXmin(),
279 h->GetZaxis()->GetXmax());
280 return -1;
281 }
282 }
283 }
284 }
285 } while ( ( h = dynamic_cast<T*> ( next() ) ) != NULL );
286 if (!h && (*next) ) {
287 Error("TProfileHelper::Merge","Attempt to merge object of class: %s to a %s",
288 (*next)->ClassName(),p->ClassName());
289 return -1;
290 }
291
292 next.Reset();
293
294 // In the case of histogram with different limits
295 // newX(Y)Axis will now have the new found limits
296 // but one needs first to clone this histogram to perform the merge
297 // The clone is not needed when all histograms have the same limits
298 T * hclone = 0;
299 if (!allSameLimits) {
300 // We don't want to add the clone to gDirectory,
301 // so remove our kMustCleanup bit temporarily
302 Bool_t mustCleanup = p->TestBit(kMustCleanup);
303 if (mustCleanup) p->ResetBit(kMustCleanup);
304 hclone = (T*)p->IsA()->New();
305 R__ASSERT(hclone);
306 hclone->SetDirectory(0);
307 p->Copy(*hclone);
308 if (mustCleanup) p->SetBit(kMustCleanup);
309 p->BufferEmpty(1); // To remove buffer.
310 p->Reset(); // BufferEmpty sets limits so we can't use it later.
311 p->SetEntries(0);
312 inlist.AddFirst(hclone);
313 }
314
315 if (!allSameLimits && initialLimitsFound) {
316 Int_t b[] = { newXAxis.GetNbins(), newYAxis.GetNbins(), newZAxis.GetNbins() };
317 Double_t v[] = { newXAxis.GetXmin(), newXAxis.GetXmax(),
318 newYAxis.GetXmin(), newYAxis.GetXmax(),
319 newZAxis.GetXmin(), newZAxis.GetXmax() };
320 p->SetBins(b, v);
321 }
322
323 if (!allHaveLimits) {
324 // fill this histogram with all the data from buffers of histograms without limits
325 while ( (h = dynamic_cast<T*> (next()) ) ) {
326 if (h->GetXaxis()->GetXmin() >= h->GetXaxis()->GetXmax() && h->fBuffer) {
327 // no limits
328 Int_t nbentries = (Int_t)h->fBuffer[0];
329 Double_t v[5];
330 for (Int_t i = 0; i < nbentries; i++)
331 if ( p->GetDimension() == 3 ) {
332 v[0] = h->fBuffer[5*i + 2];
333 v[1] = h->fBuffer[5*i + 3];
334 v[2] = h->fBuffer[5*i + 4];
335 v[3] = h->fBuffer[5*i + 5];
336 v[4] = h->fBuffer[5*i + 1];
337 p->Fill(v);
338 } else if ( p->GetDimension() == 2 ) {
339 v[0] = h->fBuffer[4*i + 2];
340 v[1] = h->fBuffer[4*i + 3];
341 v[2] = h->fBuffer[4*i + 4];
342 v[3] = h->fBuffer[4*i + 1];
343 v[4] = 0;
344 p->Fill(v);
345 }
346 else if ( p->GetDimension() == 1 ) {
347 v[0] = h->fBuffer[3*i + 2];
348 v[1] = h->fBuffer[3*i + 3];
349 v[2] = h->fBuffer[3*i + 1];
350 v[3] = v[4] = 0;
351 p->Fill(v);
352 }
353 }
354 }
355 if (!initialLimitsFound) {
356 if (hclone) {
357 inlist.Remove(hclone);
358 delete hclone;
359 }
360 return (Int_t) p->GetEntries(); // all histograms have been processed
361 }
362 next.Reset();
363 }
364
365 //merge bin contents and errors
366 Double_t stats[TH1::kNstat], totstats[TH1::kNstat];
367 for (Int_t i=0;i<TH1::kNstat;i++) {totstats[i] = stats[i] = 0;}
368 p->GetStats(totstats);
369 Double_t nentries = p->GetEntries();
370 Bool_t canExtend = p->CanExtendAllAxes();
371 p->SetCanExtend(TH1::kNoAxis); // reset, otherwise setting the under/overflow will extend the axis
372
373 while ( (h=static_cast<T*>(next())) ) {
374 // process only if the histogram has limits; otherwise it was processed before
375
376 if (h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax()) {
377 // import statistics
378 h->GetStats(stats);
379 for (Int_t i = 0; i < TH1::kNstat; i++)
380 totstats[i] += stats[i];
381 nentries += h->GetEntries();
382
383 for ( Int_t hbin = 0; hbin < h->fN; ++hbin ) {
384 Int_t pbin = hbin;
385 if (!allSameLimits) {
386 // histogram have different limits:
387 // find global bin number in p given the x,y,z axis bin numbers in h
388 // in case of non equal axes
389 // we can use FindBin on p axes because SetCanExtend(TH1::kNoAxis) has been called
390 if ( h->GetW()[hbin] != 0 && (h->IsBinUnderflow(hbin) || h->IsBinOverflow(hbin)) ) {
391 // reject cases where underflow/overflow are there and bin content is not zero
392 Error("TProfileHelper::Merge", "Cannot merge profiles - they have"
393 " different limits and underflows/overflows are present."
394 " The initial profile is now broken!");
395 return -1;
396 }
397 Int_t hbinx, hbiny, hbinz;
398 h->GetBinXYZ(hbin, hbinx, hbiny, hbinz);
399
400 pbin = p->GetBin( p->fXaxis.FindBin( h->GetXaxis()->GetBinCenter(hbinx) ),
401 p->fYaxis.FindBin( h->GetYaxis()->GetBinCenter(hbiny) ),
402 p->fZaxis.FindBin( h->GetZaxis()->GetBinCenter(hbinz) ) );
403 }
404
405
406 p->fArray[pbin] += h->GetW()[hbin];
407 p->fSumw2.fArray[pbin] += h->GetW2()[hbin];
408 p->fBinEntries.fArray[pbin] += h->GetB()[hbin];
409 if (p->fBinSumw2.fN) {
410 if ( h->GetB2() ) p->fBinSumw2.fArray[pbin] += h->GetB2()[hbin];
411 else p->fBinSumw2.fArray[pbin] += h->GetB()[hbin];
412 }
413 }
414 }
415 }
416 if (canExtend) p->SetCanExtend(TH1::kAllAxes);
417
418 //copy merged stats
419 p->PutStats(totstats);
420 p->SetEntries(nentries);
421 if (hclone) {
422 inlist.Remove(hclone);
423 delete hclone;
424 }
425 return (Long64_t)nentries;
426#endif
427}
428
429template <typename T>
431{
432// Profile histogram is resized along axis such that x is in the axis range.
433// The new axis limits are recomputed by doubling iteratively
434// the current axis range until the specified value x is within the limits.
435// The algorithm makes a copy of the histogram, then loops on all bins
436// of the old histogram to fill the extended histogram.
437// Takes into account errors (Sumw2) if any.
438// The axis must be extendable before invoking this function.
439// Ex: h->GetXaxis()->SetCanExtend(kTRUE)
440
441
442 if (!axis->CanExtend()) return 0;
443 if (axis->GetXmin() >= axis->GetXmax()) return 0;
444 if (axis->GetNbins() <= 0) return 0;
445 if (TMath::IsNaN(x)) { // x may be a NaN
446 return 0;
447 }
448
450 if (!p->FindNewAxisLimits(axis, x, xmin, xmax))
451 return 0;
452
453 //save a copy of this histogram
454 T* hold = (T*)p->IsA()->New();
455 R__ASSERT(hold);
456 hold->SetDirectory(0);
457 p->Copy(*hold);
458 //set new axis limits but keep same number of bins
459 axis->SetLimits(xmin,xmax);
460 if (p->fBinSumw2.fN) hold->Sumw2();
461
462 // total bins (including underflow /overflow)
463 Int_t nx = p->fXaxis.GetNbins() + 2;
464 Int_t ny = (p->GetDimension() > 1) ? p->fYaxis.GetNbins() + 2 : 1;
465 Int_t nz = (p->GetDimension() > 2) ? p->fZaxis.GetNbins() + 2 : 1;
466
467 Int_t iaxis = 0;
468 if (axis == p->GetXaxis()) iaxis = 1;
469 if (axis == p->GetYaxis()) iaxis = 2;
470 if (axis == p->GetZaxis()) iaxis = 3;
471 Bool_t firstw = kTRUE;
472
473 //now loop on all bins and refill
474 p->Reset("ICE"); //reset only Integral, contents and Errors
475
476 // need to consider also underflow/overflow in the non-extending axes
477 Double_t xc,yc,zc;
478 Int_t ix, iy, iz, binx, biny, binz;
479 for (binz=0;binz< nz;binz++) {
480 zc = hold->GetZaxis()->GetBinCenter(binz);
481 iz = p->fZaxis.FindFixBin(zc);
482 for (biny=0;biny<ny;biny++) {
483 yc = hold->GetYaxis()->GetBinCenter(biny);
484 iy = p->fYaxis.FindFixBin(yc);
485 for (binx=0;binx<nx;binx++) {
486 xc = hold->GetXaxis()->GetBinCenter(binx);
487 ix = p->fXaxis.FindFixBin(xc);
488 Int_t sourceBin = hold->GetBin(binx,biny,binz);
489 // skip empty bins
490 if (hold->fBinEntries.fArray[sourceBin] == 0) continue;
491 if (hold->IsBinUnderflow(sourceBin, iaxis) || hold->IsBinOverflow(sourceBin, iaxis)) {
492 if (firstw) {
493 Warning("ExtendAxis",
494 "Histogram %s has underflow or overflow in the %s that is extendable"
495 " their content will be lost",p->GetName(),axis->GetName());
496 firstw = kFALSE;
497 }
498 continue;
499 }
500 Int_t destinationBin = p->GetBin(ix,iy,iz);
501 p->AddBinContent(destinationBin, hold->fArray[sourceBin]);
502 p->fBinEntries.fArray[destinationBin] += hold->fBinEntries.fArray[sourceBin];
503 p->fSumw2.fArray[destinationBin] += hold->fSumw2.fArray[sourceBin];
504 if (p->fBinSumw2.fN) p->fBinSumw2.fArray[destinationBin] += hold->fBinSumw2.fArray[sourceBin];
505 }
506 }
507 }
508 return hold;
509}
510
511template <typename T>
513{
514 Double_t ac1 = TMath::Abs(c1);
515
516 // Make the loop over the bins to calculate the Addition
517 Int_t bin;
518 Double_t *cu1 = p->GetW();
519 Double_t *er1 = p->GetW2();
520 Double_t *en1 = p->GetB();
521 for (bin=0;bin<p->fN;bin++) {
522 p->fArray[bin] = c1*cu1[bin];
523 p->fSumw2.fArray[bin] = ac1*ac1*er1[bin];
524 p->fBinEntries.fArray[bin] = en1[bin];
525 }
526}
527
528template <typename T>
530{
531 // Create/Delete structure to store sum of squares of weights per bin *-*-*-*-*-*-*-*
532 // This is needed to compute the correct statistical quantities
533 // of a profile filled with weights
534 //
535 //
536 // This function is automatically called when the histogram is created
537 // if the static function TH1::SetDefaultSumw2 has been called before.
538
539 if (!flag) {
540 // clear array if existing or do nothing
541 if (p->fBinSumw2.fN > 0 ) p->fBinSumw2.Set(0);
542 return;
543 }
544
545 if ( p->fBinSumw2.fN == p->fNcells) {
546 if (!p->fgDefaultSumw2)
547 Warning("Sumw2","Sum of squares of profile bin weights structure already created");
548 return;
549 }
550
551 p->fBinSumw2.Set(p->fNcells);
552
553 // by default fill with the sum of weights which are stored in fBinEntries
554 for (Int_t bin=0; bin<p->fNcells; bin++) {
555 p->fBinSumw2.fArray[bin] = p->fBinEntries.fArray[bin];
556 }
557}
558
559template <typename T>
561{
562 // Reduce the number of bins for this axis to the number of bins having a label.
563 // Works only for the given axis passed in the option
564 // The method will remove only the extra bins existing after the last "labeled" bin.
565 // Note that if there are "un-labeled" bins present between "labeled" bins they will not be removed
566
567
568 TAxis *axis = p->GetXaxis();
569 if (ax[0] == 'y' || ax[0] == 'Y') axis = p->GetYaxis();
570 if (ax[0] == 'z' || ax[0] == 'Z') axis = p->GetZaxis();
571 if (!axis) {
572 Error("TProfileHelper::LabelsDeflate","Invalid axis option %s",ax);
573 return;
574 }
575 if (!axis->GetLabels()) return;
576
577 // find bin with last labels
578 // bin number is object ID in list of labels
579 // therefore max bin number is number of bins of the deflated histograms
580 TIter next(axis->GetLabels());
581 TObject *obj;
582 Int_t nbins = 0;
583 while ((obj = next())) {
584 Int_t ibin = obj->GetUniqueID();
585 if (ibin > nbins) nbins = ibin;
586 }
587 if (nbins < 1) nbins = 1;
588
589 // do nothing in case it was the last bin
590 if (nbins==axis->GetNbins()) return;
591
592 T *hold = (T*)p->IsA()->New();;
593 hold->SetDirectory(0);
594 p->Copy(*hold);
595
596
597 Double_t xmin = axis->GetXmin();
598 Double_t xmax = axis->GetBinUpEdge(nbins);
599 axis->SetRange(0,0);
600 // set the new bins and range
601 axis->Set(nbins,xmin,xmax);
602 p->SetBinsLength(-1); // reset the number of cells
603 p->fBinEntries.Set(p->fN);
604 p->fSumw2.Set(p->fN);
605 if (p->fBinSumw2.fN) p->fBinSumw2.Set(p->fN);
606
607 // reset the content
608 p->Reset("ICE");
609
610 //now loop on all bins and refill
611 Int_t bin,binx,biny,binz;
612 for (bin =0; bin < hold->fN; ++bin)
613 {
614 hold->GetBinXYZ(bin, binx, biny, binz);
615 Int_t ibin = p->GetBin(binx, biny, binz);
616 p->fArray[ibin] += hold->fArray[bin];
617 p->fBinEntries.fArray[ibin] += hold->fBinEntries.fArray[bin];
618 p->fSumw2.fArray[ibin] += hold->fSumw2.fArray[bin];
619 if (p->fBinSumw2.fN) p->fBinSumw2.fArray[ibin] += hold->fBinSumw2.fArray[bin];
620 }
621
622 delete hold;
623}
624
625template <typename T>
627{
628// Double the number of bins for axis.
629// Refill histogram
630// This function is called by TAxis::FindBin(const char *label)
631// Works only for the given axis
632
633 if (gDebug) Info("LabelsInflate","Inflate label for axis %s of profile %s",ax,p->GetName());
634
635 Int_t iaxis = p->AxisChoice(ax);
636 TAxis *axis = 0;
637 if (iaxis == 1) axis = p->GetXaxis();
638 if (iaxis == 2) axis = p->GetYaxis();
639 if (iaxis == 3) axis = p->GetZaxis();
640 if (!axis) return;
641 // TAxis *axis = p->GetXaxis();
642 // if (ax[0] == 'y' || ax[0] == 'Y') axis = p->GetYaxis();
643
644 T *hold = (T*)p->IsA()->New();;
645 hold->SetDirectory(0);
646 p->Copy(*hold);
647
648
649// Int_t nbxold = p->fXaxis.GetNbins();
650// Int_t nbyold = p->fYaxis.GetNbins();
651 Int_t nbins = axis->GetNbins();
652 Double_t xmin = axis->GetXmin();
653 Double_t xmax = axis->GetXmax();
654 xmax = xmin + 2*(xmax-xmin);
655 axis->SetRange(0,0);
656 // double the number of bins
657 nbins *= 2;
658 axis->Set(nbins,xmin,xmax);
659 // reset the array of content according to the axis
660 p->SetBinsLength(-1);
661 Int_t ncells = p->fN;
662 p->fBinEntries.Set(ncells);
663 p->fSumw2.Set(ncells);
664 if (p->fBinSumw2.fN) p->fBinSumw2.Set(ncells);
665
666 p->Reset("ICE"); // reset content and error
667
668 //now loop on all old bins and refill excluding underflow/overflow in
669 // the axis that has the bin doubled
670 Int_t binx, biny, binz = 0;
671 for (Int_t ibin =0; ibin < hold->fNcells; ibin++) {
672 // get the binx,y,z values . The x-y-z (axis) bin values will stay the same between new-old after the expanding
673 hold->GetBinXYZ(ibin,binx,biny,binz);
674 Int_t bin = p->GetBin(binx,biny,binz);
675
676 // underflow and overflow will be cleaned up because their meaning has been altered
677 if (hold->IsBinUnderflow(ibin,iaxis) || hold->IsBinOverflow(ibin,iaxis)) {
678 if (gDebug && hold->fBinEntries.fArray[ibin] > 0) Info("LabelsInflate","Content for underflow/overflow of bin (%d,%d,%d) will be lost",binx,biny,binz);
679 continue;
680 }
681 else {
682 p->fArray[bin] = hold->fArray[ibin];
683 p->fBinEntries.fArray[bin] = hold->fBinEntries.fArray[ibin];
684 p->fSumw2.fArray[bin] = hold->fSumw2.fArray[ibin];
685 if (p->fBinSumw2.fN) p->fBinSumw2.fArray[bin] = hold->fBinSumw2.fArray[ibin];
686 if (gDebug) Info("LabelsInflate","Copy Content from bin (%d,%d,%d) from %d in %d (%f,%f)",binx,biny,binz, ibin, bin, hold->fArray[ibin],hold->fBinEntries.fArray[ibin] );
687 }
688 }
689 delete hold;
690}
691
692template <typename T>
694 // set the profile option
695 TString opt = option;
696 opt.ToLower();
697 p->fErrorMode = kERRORMEAN;
698 if (opt.Contains("s")) p->fErrorMode = kERRORSPREAD;
699 if (opt.Contains("i")) p->fErrorMode = kERRORSPREADI;
700 if (opt.Contains("g")) p->fErrorMode = kERRORSPREADG;
701}
702
703template <typename T>
705{
706 // compute bin error of profile histograms
707
708 if (p->fBuffer) p->BufferEmpty();
709
710 if (bin < 0 || bin >= p->fNcells) return 0;
711 Double_t cont = p->fArray[bin]; // sum of bin w *y
712 Double_t sum = p->fBinEntries.fArray[bin]; // sum of bin weights
713 Double_t err2 = p->fSumw2.fArray[bin]; // sum of bin w * y^2
714 Double_t neff = p->GetBinEffectiveEntries(bin); // (sum of w)^2 / (sum of w^2)
715 if (sum == 0) return 0; // for empty bins
716 // case the values y are gaussian distributed y +/- sigma and w = 1/sigma^2
717 if (p->fErrorMode == kERRORSPREADG) {
718 return 1./TMath::Sqrt(sum);
719 }
720 // compute variance in y (eprim2) and standard deviation in y (eprim)
721 Double_t contsum = cont/sum;
722 Double_t eprim2 = TMath::Abs(err2/sum - contsum*contsum);
723 Double_t eprim = TMath::Sqrt(eprim2);
724
725 if (p->fErrorMode == kERRORSPREADI) {
726 if (eprim != 0) return eprim/TMath::Sqrt(neff);
727 // in case content y is an integer (so each my has an error +/- 1/sqrt(12)
728 // when the std(y) is zero
729 return 1/TMath::Sqrt(12*neff);
730 }
731
732 // if approximate compute the sums (of w, wy and wy2) using all the bins
733 // when the variance in y is zero
734 Double_t test = 1;
735 if (err2 != 0 && neff < 5) test = eprim2*sum/err2;
736 //Int_t cellLimit = (p->GetDimension() == 3)?1000404:10404;
737 if (p->fgApproximate && (test < 1.e-4 || eprim2 <= 0)) {
738 Double_t stats[TH1::kNstat] = {0};
739 p->GetStats(stats);
740 Double_t ssum = stats[0];
741 // for 1D profile
742 int index = 4; // index in the stats array for 1D
743 if (p->GetDimension() == 2) index = 7; // for 2D
744 if (p->GetDimension() == 3) index = 11; // for 3D
745 Double_t scont = stats[index];
746 Double_t serr2 = stats[index+1];
747
748 // compute mean and variance in y
749 Double_t scontsum = scont/ssum; // global mean
750 Double_t seprim2 = TMath::Abs(serr2/ssum - scontsum*scontsum); // global variance
751 eprim = 2*TMath::Sqrt(seprim2); // global std (why factor of 2 ??)
752 sum = ssum;
753 }
754 sum = TMath::Abs(sum);
755
756 // case option "S" return standard deviation in y
757 if (p->fErrorMode == kERRORSPREAD) return eprim;
758
759 // default case : fErrorMode = kERRORMEAN
760 // return standard error on the mean of y
761 //if (neff == 0) std::cerr << "NEFF = 0 for bin " << bin << " " << eprim << " " << neff << " " << std::endl;
762 return eprim/TMath::Sqrt(neff);
763
764}
765
766
767template <typename T>
769// Set the number of entries in bin for the profile
770// In case the profile stores the sum of weight squares - set the sum of weight square to the number entries
771// (i.e. assume the entries have been filled all with a weight == 1 )
772
773 if (bin < 0 || bin >= p->fNcells) return;
774 p->fBinEntries.fArray[bin] = w;
775 if (p->fBinSumw2.fN) p->fBinSumw2.fArray[bin] = w;
776
777}
778
779#endif
#define b(i)
Definition: RSha256.hxx:100
#define s0(x)
Definition: RSha256.hxx:90
#define s1(x)
Definition: RSha256.hxx:91
#define h(i)
Definition: RSha256.hxx:106
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
bool Bool_t
Definition: RtypesCore.h:63
double Double_t
Definition: RtypesCore.h:59
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 R__ASSERT(e)
Definition: TError.h:118
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition: TError.cxx:220
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition: TError.cxx:231
float xmin
Definition: THbookFile.cxx:95
int nentries
Definition: THbookFile.cxx:91
float xmax
Definition: THbookFile.cxx:95
@ kMustCleanup
Definition: TObject.h:355
@ kERRORSPREAD
Definition: TProfile.h:28
@ kERRORSPREADG
Definition: TProfile.h:28
@ kERRORSPREADI
Definition: TProfile.h:28
@ kERRORMEAN
Definition: TProfile.h:28
Int_t gDebug
Definition: TROOT.cxx:592
Class to manage histogram axis.
Definition: TAxis.h:30
Bool_t CanExtend() const
Definition: TAxis.h:82
Double_t GetXmax() const
Definition: TAxis.h:134
virtual void Set(Int_t nbins, Double_t xmin, Double_t xmax)
Initialize axis with fix bins.
Definition: TAxis.cxx:731
virtual void SetLimits(Double_t xmin, Double_t xmax)
Definition: TAxis.h:154
Double_t GetXmin() const
Definition: TAxis.h:133
Int_t GetNbins() const
Definition: TAxis.h:121
virtual void SetRange(Int_t first=0, Int_t last=0)
Set the viewing range for the axis using bin numbers.
Definition: TAxis.cxx:920
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition: TAxis.cxx:528
THashList * GetLabels() const
Definition: TAxis.h:117
Collection abstract base class.
Definition: TCollection.h:65
virtual void AddAll(const TCollection *col)
Add all objects from collection col to this collection.
virtual Bool_t IsEmpty() const
Definition: TCollection.h:188
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
@ kAllAxes
Definition: TH1.h:75
@ kNoAxis
NOTE: Must always be 0 !!!
Definition: TH1.h:71
static Bool_t GetDefaultSumw2()
Return kTRUE if TH1::Sumw2 must be called when creating new histograms.
Definition: TH1.cxx:4379
@ kNstat
Size of statistics data (up to TProfile3D)
Definition: TH1.h:183
virtual Int_t BufferEmpty(Int_t action=0)
Fill histogram with all entries in the buffer.
Definition: TH1.cxx:1403
void Reset()
Definition: TCollection.h:254
A doubly linked list.
Definition: TList.h:44
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:822
virtual void AddFirst(TObject *obj)
Add object at the beginning of the list.
Definition: TList.cxx:100
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Mother of all ROOT objects.
Definition: TObject.h:37
virtual UInt_t GetUniqueID() const
Return the unique object id.
Definition: TObject.cxx:377
static void LabelsInflate(T *p, Option_t *)
static Double_t GetBinError(T *p, Int_t bin)
static T * ExtendAxis(T *p, Double_t x, TAxis *axis)
static void Sumw2(T *p, Bool_t flag)
static void SetBinEntries(T *p, Int_t bin, Double_t w)
static void Scale(T *p, Double_t c1, Option_t *option)
static void SetErrorOption(T *p, Option_t *opt)
static Long64_t Merge(T *p, TCollection *list)
static void BuildArray(T *p)
static Bool_t Add(T *p, const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2=1)
static Double_t GetBinEffectiveEntries(T *p, Int_t bin)
static void LabelsDeflate(T *p, Option_t *)
Basic string class.
Definition: TString.h:136
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1150
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
return c1
Definition: legend1.C:41
Double_t x[n]
Definition: legend1.C:17
for(Int_t i=0;i< n;i++)
Definition: legend1.C:18
TH1F * h1
Definition: legend1.C:5
return c2
Definition: legend2.C:14
double T(double x)
Definition: ChebyshevPol.h:34
Bool_t IsNaN(Double_t x)
Definition: TMath.h:842
Double_t Sqrt(Double_t x)
Definition: TMath.h:641
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
Definition: test.py:1
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345