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 "TCollection.h"
27#include "THashList.h"
28#include "TMath.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// Make the loop over the bins to calculate the Addition
110 Int_t bin;
111 Double_t *cu1 = p1->GetW(); Double_t *cu2 = p2->GetW();
112 Double_t *er1 = p1->GetW2(); Double_t *er2 = p2->GetW2();
113 Double_t *en1 = p1->GetB(); Double_t *en2 = p2->GetB();
114 Double_t *ew1 = p1->GetB2(); Double_t *ew2 = p2->GetB2();
115 // create sumw2 per bin if not set
116 if (p->fBinSumw2.fN == 0 && (p1->fBinSumw2.fN != 0 || p2->fBinSumw2.fN != 0) ) p->Sumw2();
117 // if p1 has not the sum of weight squared/bin stored use just the sum of weights
118 if (ew1 == 0) ew1 = en1;
119 if (ew2 == 0) ew2 = en2;
120 for (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 TAxis newXAxis;
188 TAxis newYAxis;
189 TAxis newZAxis;
190 Bool_t initialLimitsFound = kFALSE;
191 Bool_t allSameLimits = kTRUE;
192 Bool_t allHaveLimits = kTRUE;
193// Bool_t firstNonEmptyHist = kTRUE;
194
195 TIter next(&inlist);
196 T* h = p;
197
198 do {
199
200 // skip empty histograms
201 // if (h->fTsumw == 0 && h->GetEntries() == 0) continue;
202
203 Bool_t hasLimits = h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax();
204 allHaveLimits = allHaveLimits && hasLimits;
205
206 if (hasLimits) {
207 h->BufferEmpty();
208
209#ifdef LATER
210 // this is done in case the first histograms are empty and
211 // the histogram have different limits
212 if (firstNonEmptyHist ) {
213 // set axis limits in the case the first histogram was empty
214 if (h != p ) {
215 if (!p->SameLimitsAndNBins(p->fXaxis, *(h->GetXaxis())) )
216 p->fXaxis.Set(h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax());
217 if (!p->SameLimitsAndNBins(p->fYaxis, *(h->GetYaxis())) )
218 p->fYaxis.Set(h->GetYaxis()->GetNbins(), h->GetYaxis()->GetXmin(),h->GetYaxis()->GetXmax());
219 if (!p->SameLimitsAndNBins(p->fZaxis, *(h->GetZaxis())) )
220 p->fZaxis.Set(h->GetZaxis()->GetNbins(), h->GetZaxis()->GetXmin(),h->GetZaxis()->GetXmax());
221 }
222 firstNonEmptyHist = kFALSE;
223 }
224#endif
225
226 // this is executed the first time an histogram with limits is found
227 // to set some initial values on the new axis
228 if (!initialLimitsFound) {
229 initialLimitsFound = kTRUE;
230 newXAxis.Set(h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),
231 h->GetXaxis()->GetXmax());
232 if ( p->GetDimension() >= 2 )
233 newYAxis.Set(h->GetYaxis()->GetNbins(), h->GetYaxis()->GetXmin(),
234 h->GetYaxis()->GetXmax());
235 if ( p->GetDimension() >= 3 )
236 newZAxis.Set(h->GetZaxis()->GetNbins(), h->GetZaxis()->GetXmin(),
237 h->GetZaxis()->GetXmax());
238 }
239 else {
240 // check first if histograms have same bins
241 if (!p->SameLimitsAndNBins(newXAxis, *(h->GetXaxis())) ||
242 !p->SameLimitsAndNBins(newYAxis, *(h->GetYaxis())) ||
243 !p->SameLimitsAndNBins(newZAxis, *(h->GetZaxis())) ) {
244
245 allSameLimits = kFALSE;
246 // recompute in this case the optimal limits
247 // The condition to works is that the histogram have same bin with
248 // and one common bin edge
249
250 if (!p->RecomputeAxisLimits(newXAxis, *(h->GetXaxis()))) {
251 Error("TProfileHelper::Merge", "Cannot merge profiles %d dim - limits are inconsistent:\n "
252 "first: (%d, %f, %f), second: (%d, %f, %f)",p->GetDimension(),
253 newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax(),
254 h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),
255 h->GetXaxis()->GetXmax());
256 return -1;
257 }
258 if (p->GetDimension() >= 2 && !p->RecomputeAxisLimits(newYAxis, *(h->GetYaxis()))) {
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 newYAxis.GetNbins(), newYAxis.GetXmin(), newYAxis.GetXmax(),
262 h->GetYaxis()->GetNbins(), h->GetYaxis()->GetXmin(),
263 h->GetYaxis()->GetXmax());
264 return -1;
265 }
266 if (p->GetDimension() >= 3 && !p->RecomputeAxisLimits(newZAxis, *(h->GetZaxis()))) {
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 newZAxis.GetNbins(), newZAxis.GetXmin(), newZAxis.GetXmax(),
270 h->GetZaxis()->GetNbins(), h->GetZaxis()->GetXmin(),
271 h->GetZaxis()->GetXmax());
272 return -1;
273 }
274 }
275 }
276 }
277 } while ( ( h = dynamic_cast<T*> ( next() ) ) != NULL );
278 if (!h && (*next) ) {
279 Error("TProfileHelper::Merge","Attempt to merge object of class: %s to a %s",
280 (*next)->ClassName(),p->ClassName());
281 return -1;
282 }
283
284 next.Reset();
285
286 // In the case of histogram with different limits
287 // newX(Y)Axis will now have the new found limits
288 // but one needs first to clone this histogram to perform the merge
289 // The clone is not needed when all histograms have the same limits
290 T * hclone = 0;
291 if (!allSameLimits) {
292 // We don't want to add the clone to gDirectory,
293 // so remove our kMustCleanup bit temporarily
294 Bool_t mustCleanup = p->TestBit(kMustCleanup);
295 if (mustCleanup) p->ResetBit(kMustCleanup);
296 hclone = (T*)p->IsA()->New();
297 R__ASSERT(hclone);
298 hclone->SetDirectory(0);
299 p->Copy(*hclone);
300 if (mustCleanup) p->SetBit(kMustCleanup);
301 p->BufferEmpty(1); // To remove buffer.
302 p->Reset(); // BufferEmpty sets limits so we can't use it later.
303 p->SetEntries(0);
304 inlist.AddFirst(hclone);
305 }
306
307 if (!allSameLimits && initialLimitsFound) {
308 Int_t b[] = { newXAxis.GetNbins(), newYAxis.GetNbins(), newZAxis.GetNbins() };
309 Double_t v[] = { newXAxis.GetXmin(), newXAxis.GetXmax(),
310 newYAxis.GetXmin(), newYAxis.GetXmax(),
311 newZAxis.GetXmin(), newZAxis.GetXmax() };
312 p->SetBins(b, v);
313 }
314
315 if (!allHaveLimits) {
316 // fill this histogram with all the data from buffers of histograms without limits
317 while ( (h = dynamic_cast<T*> (next()) ) ) {
318 if (h->GetXaxis()->GetXmin() >= h->GetXaxis()->GetXmax() && h->fBuffer) {
319 // no limits
320 Int_t nbentries = (Int_t)h->fBuffer[0];
321 Double_t v[5];
322 for (Int_t i = 0; i < nbentries; i++)
323 if ( p->GetDimension() == 3 ) {
324 v[0] = h->fBuffer[5*i + 2];
325 v[1] = h->fBuffer[5*i + 3];
326 v[2] = h->fBuffer[5*i + 4];
327 v[3] = h->fBuffer[5*i + 5];
328 v[4] = h->fBuffer[5*i + 1];
329 p->Fill(v);
330 } else if ( p->GetDimension() == 2 ) {
331 v[0] = h->fBuffer[4*i + 2];
332 v[1] = h->fBuffer[4*i + 3];
333 v[2] = h->fBuffer[4*i + 4];
334 v[3] = h->fBuffer[4*i + 1];
335 v[4] = 0;
336 p->Fill(v);
337 }
338 else if ( p->GetDimension() == 1 ) {
339 v[0] = h->fBuffer[3*i + 2];
340 v[1] = h->fBuffer[3*i + 3];
341 v[2] = h->fBuffer[3*i + 1];
342 v[3] = v[4] = 0;
343 p->Fill(v);
344 }
345 }
346 }
347 if (!initialLimitsFound) {
348 if (hclone) {
349 inlist.Remove(hclone);
350 delete hclone;
351 }
352 return (Int_t) p->GetEntries(); // all histograms have been processed
353 }
354 next.Reset();
355 }
356
357 //merge bin contents and errors
358 Double_t stats[TH1::kNstat], totstats[TH1::kNstat];
359 for (Int_t i=0;i<TH1::kNstat;i++) {totstats[i] = stats[i] = 0;}
360 p->GetStats(totstats);
361 Double_t nentries = p->GetEntries();
362 Bool_t canExtend = p->CanExtendAllAxes();
363 p->SetCanExtend(TH1::kNoAxis); // reset, otherwise setting the under/overflow will extend the axis
364
365 while ( (h=static_cast<T*>(next())) ) {
366 // process only if the histogram has limits; otherwise it was processed before
367
368 if (h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax()) {
369 // import statistics
370 h->GetStats(stats);
371 for (Int_t i = 0; i < TH1::kNstat; i++)
372 totstats[i] += stats[i];
373 nentries += h->GetEntries();
374
375 for ( Int_t hbin = 0; hbin < h->fN; ++hbin ) {
376 Int_t pbin = hbin;
377 if (!allSameLimits) {
378 // histogram have different limits:
379 // find global bin number in p given the x,y,z axis bin numbers in h
380 // in case of non equal axes
381 // we can use FindBin on p axes because SetCanExtend(TH1::kNoAxis) has been called
382 if ( h->GetW()[hbin] != 0 && (h->IsBinUnderflow(hbin) || h->IsBinOverflow(hbin)) ) {
383 // reject cases where underflow/overflow are there and bin content is not zero
384 Error("TProfileHelper::Merge", "Cannot merge profiles - they have"
385 " different limits and underflows/overflows are present."
386 " The initial profile is now broken!");
387 return -1;
388 }
389 Int_t hbinx, hbiny, hbinz;
390 h->GetBinXYZ(hbin, hbinx, hbiny, hbinz);
391
392 pbin = p->GetBin( p->fXaxis.FindBin( h->GetXaxis()->GetBinCenter(hbinx) ),
393 p->fYaxis.FindBin( h->GetYaxis()->GetBinCenter(hbiny) ),
394 p->fZaxis.FindBin( h->GetZaxis()->GetBinCenter(hbinz) ) );
395 }
396
397
398 p->fArray[pbin] += h->GetW()[hbin];
399 p->fSumw2.fArray[pbin] += h->GetW2()[hbin];
400 p->fBinEntries.fArray[pbin] += h->GetB()[hbin];
401 if (p->fBinSumw2.fN) {
402 if ( h->GetB2() ) p->fBinSumw2.fArray[pbin] += h->GetB2()[hbin];
403 else p->fBinSumw2.fArray[pbin] += h->GetB()[hbin];
404 }
405 }
406 }
407 }
408 if (canExtend) p->SetCanExtend(TH1::kAllAxes);
409
410 //copy merged stats
411 p->PutStats(totstats);
412 p->SetEntries(nentries);
413 if (hclone) {
414 inlist.Remove(hclone);
415 delete hclone;
416 }
417 return (Long64_t)nentries;
418}
419
420template <typename T>
422{
423// Profile histogram is resized along axis such that x is in the axis range.
424// The new axis limits are recomputed by doubling iteratively
425// the current axis range until the specified value x is within the limits.
426// The algorithm makes a copy of the histogram, then loops on all bins
427// of the old histogram to fill the extended histogram.
428// Takes into account errors (Sumw2) if any.
429// The axis must be extendable before invoking this function.
430// Ex: h->GetXaxis()->SetCanExtend(kTRUE)
431
432
433 if (!axis->CanExtend()) return 0;
434 if (axis->GetXmin() >= axis->GetXmax()) return 0;
435 if (axis->GetNbins() <= 0) return 0;
436
438 if (!p->FindNewAxisLimits(axis, x, xmin, xmax))
439 return 0;
440
441 //save a copy of this histogram
442 T* hold = (T*)p->IsA()->New();
443 R__ASSERT(hold);
444 hold->SetDirectory(0);
445 p->Copy(*hold);
446 //set new axis limits
447 axis->SetLimits(xmin,xmax);
448 if (p->fBinSumw2.fN) hold->Sumw2();
449
450 Int_t nbinsx = p->fXaxis.GetNbins();
451 Int_t nbinsy = p->fYaxis.GetNbins();
452 Int_t nbinsz = p->fZaxis.GetNbins();
453
454 //now loop on all bins and refill
455 p->Reset("ICE"); //reset only Integral, contents and Errors
456
457 Double_t bx,by,bz;
458 Int_t ix, iy, iz, binx, biny, binz;
459 for (binz=1;binz<=nbinsz;binz++) {
460 bz = hold->GetZaxis()->GetBinCenter(binz);
461 iz = p->fZaxis.FindFixBin(bz);
462 for (biny=1;biny<=nbinsy;biny++) {
463 by = hold->GetYaxis()->GetBinCenter(biny);
464 iy = p->fYaxis.FindFixBin(by);
465 for (binx=1;binx<=nbinsx;binx++) {
466 bx = hold->GetXaxis()->GetBinCenter(binx);
467 ix = p->fXaxis.FindFixBin(bx);
468
469 Int_t sourceBin = hold->GetBin(binx,biny,binz);
470 Int_t destinationBin = p->GetBin(ix,iy,iz);
471 p->AddBinContent(destinationBin, hold->fArray[sourceBin]);
472 p->fBinEntries.fArray[destinationBin] += hold->fBinEntries.fArray[sourceBin];
473 p->fSumw2.fArray[destinationBin] += hold->fSumw2.fArray[sourceBin];
474 if (p->fBinSumw2.fN) p->fBinSumw2.fArray[destinationBin] += hold->fBinSumw2.fArray[sourceBin];
475 }
476 }
477 }
478 return hold;
479}
480
481template <typename T>
483{
484 Double_t ac1 = TMath::Abs(c1);
485
486 // Make the loop over the bins to calculate the Addition
487 Int_t bin;
488 Double_t *cu1 = p->GetW();
489 Double_t *er1 = p->GetW2();
490 Double_t *en1 = p->GetB();
491 for (bin=0;bin<p->fN;bin++) {
492 p->fArray[bin] = c1*cu1[bin];
493 p->fSumw2.fArray[bin] = ac1*ac1*er1[bin];
494 p->fBinEntries.fArray[bin] = en1[bin];
495 }
496}
497
498template <typename T>
500{
501 // Create/Delete structure to store sum of squares of weights per bin *-*-*-*-*-*-*-*
502 // This is needed to compute the correct statistical quantities
503 // of a profile filled with weights
504 //
505 //
506 // This function is automatically called when the histogram is created
507 // if the static function TH1::SetDefaultSumw2 has been called before.
508
509 if (!flag) {
510 // clear array if existing or do nothing
511 if (p->fBinSumw2.fN > 0 ) p->fBinSumw2.Set(0);
512 return;
513 }
514
515 if ( p->fBinSumw2.fN == p->fNcells) {
516 if (!p->fgDefaultSumw2)
517 Warning("Sumw2","Sum of squares of profile bin weights structure already created");
518 return;
519 }
520
521 p->fBinSumw2.Set(p->fNcells);
522
523 // by default fill with the sum of weights which are stored in fBinEntries
524 for (Int_t bin=0; bin<p->fNcells; bin++) {
525 p->fBinSumw2.fArray[bin] = p->fBinEntries.fArray[bin];
526 }
527}
528
529template <typename T>
531{
532 // Reduce the number of bins for this axis to the number of bins having a label.
533 // Works only for the given axis passed in the option
534 // The method will remove only the extra bins existing after the last "labeled" bin.
535 // Note that if there are "un-labeled" bins present between "labeled" bins they will not be removed
536
537
538 TAxis *axis = p->GetXaxis();
539 if (ax[0] == 'y' || ax[0] == 'Y') axis = p->GetYaxis();
540 if (ax[0] == 'z' || ax[0] == 'Z') axis = p->GetZaxis();
541 if (!axis) {
542 Error("TProfileHelper::LabelsDeflate","Invalid axis option %s",ax);
543 return;
544 }
545 if (!axis->GetLabels()) return;
546
547 // find bin with last labels
548 // bin number is object ID in list of labels
549 // therefore max bin number is number of bins of the deflated histograms
550 TIter next(axis->GetLabels());
551 TObject *obj;
552 Int_t nbins = 0;
553 while ((obj = next())) {
554 Int_t ibin = obj->GetUniqueID();
555 if (ibin > nbins) nbins = ibin;
556 }
557 if (nbins < 1) nbins = 1;
558
559 // do nothing in case it was the last bin
560 if (nbins==axis->GetNbins()) return;
561
562 T *hold = (T*)p->IsA()->New();;
563 hold->SetDirectory(0);
564 p->Copy(*hold);
565
566
567 Double_t xmin = axis->GetXmin();
568 Double_t xmax = axis->GetBinUpEdge(nbins);
569 axis->SetRange(0,0);
570 // set the new bins and range
571 axis->Set(nbins,xmin,xmax);
572 p->SetBinsLength(-1); // reset the number of cells
573 p->fBinEntries.Set(p->fN);
574 p->fSumw2.Set(p->fN);
575 if (p->fBinSumw2.fN) p->fBinSumw2.Set(p->fN);
576
577 // reset the content
578 p->Reset("ICE");
579
580 //now loop on all bins and refill
581 Int_t bin,binx,biny,binz;
582 for (bin =0; bin < hold->fN; ++bin)
583 {
584 hold->GetBinXYZ(bin, binx, biny, binz);
585 Int_t ibin = p->GetBin(binx, biny, binz);
586 p->fArray[ibin] += hold->fArray[bin];
587 p->fBinEntries.fArray[ibin] += hold->fBinEntries.fArray[bin];
588 p->fSumw2.fArray[ibin] += hold->fSumw2.fArray[bin];
589 if (p->fBinSumw2.fN) p->fBinSumw2.fArray[ibin] += hold->fBinSumw2.fArray[bin];
590 }
591
592 delete hold;
593}
594
595template <typename T>
597{
598// Double the number of bins for axis.
599// Refill histogram
600// This function is called by TAxis::FindBin(const char *label)
601// Works only for the given axis
602
603
604 TAxis *axis = p->GetXaxis();
605 if (ax[0] == 'y' || ax[0] == 'Y') axis = p->GetYaxis();
606 T *hold = (T*)p->IsA()->New();;
607 hold->SetDirectory(0);
608 p->Copy(*hold);
609
610
611// Int_t nbxold = p->fXaxis.GetNbins();
612// Int_t nbyold = p->fYaxis.GetNbins();
613 Int_t nbins = axis->GetNbins();
614 Double_t xmin = axis->GetXmin();
615 Double_t xmax = axis->GetXmax();
616 xmax = xmin + 2*(xmax-xmin);
617 axis->SetRange(0,0);
618 // double the number of bins
619 nbins *= 2;
620 axis->Set(nbins,xmin,xmax);
621 // reset the array of content according to the axis
622 p->SetBinsLength(-1);
623 Int_t ncells = p->fN;
624 p->fBinEntries.Set(ncells);
625 p->fSumw2.Set(ncells);
626 if (p->fBinSumw2.fN) p->fBinSumw2.Set(ncells);
627
628 //now loop on all bins and refill
629 for (Int_t ibin =0; ibin < p->fN; ibin++)
630 {
631 Int_t binx, biny, binz;
632 p->GetBinXYZ(ibin, binx, biny, binz);
633 Int_t bin = hold->GetBin(binx, biny, binz);
634
635 if (p->IsBinUnderflow(ibin) || p->IsBinOverflow(ibin)) {
636 p->UpdateBinContent(ibin, 0.0);
637 p->fBinEntries.fArray[ibin] = 0.0;
638 p->fSumw2.fArray[ibin] = 0.0;
639 if (p->fBinSumw2.fN) p->fBinSumw2.fArray[ibin] = 0.0;
640 } else {
641 p->fArray[ibin] = hold->fArray[bin];
642 p->fBinEntries.fArray[ibin] = hold->fBinEntries.fArray[bin];
643 p->fSumw2.fArray[ibin] = hold->fSumw2.fArray[bin];
644 if (p->fBinSumw2.fN) p->fBinSumw2.fArray[ibin] = hold->fBinSumw2.fArray[bin];
645 }
646 }
647 delete hold;
648}
649
650template <typename T>
652 // set the profile option
653 TString opt = option;
654 opt.ToLower();
655 p->fErrorMode = kERRORMEAN;
656 if (opt.Contains("s")) p->fErrorMode = kERRORSPREAD;
657 if (opt.Contains("i")) p->fErrorMode = kERRORSPREADI;
658 if (opt.Contains("g")) p->fErrorMode = kERRORSPREADG;
659}
660
661template <typename T>
663{
664 // compute bin error of profile histograms
665
666 if (p->fBuffer) p->BufferEmpty();
667
668 if (bin < 0 || bin >= p->fNcells) return 0;
669 Double_t cont = p->fArray[bin]; // sum of bin w *y
670 Double_t sum = p->fBinEntries.fArray[bin]; // sum of bin weights
671 Double_t err2 = p->fSumw2.fArray[bin]; // sum of bin w * y^2
672 Double_t neff = p->GetBinEffectiveEntries(bin); // (sum of w)^2 / (sum of w^2)
673 if (sum == 0) return 0; // for empty bins
674 // case the values y are gaussian distributed y +/- sigma and w = 1/sigma^2
675 if (p->fErrorMode == kERRORSPREADG) {
676 return 1./TMath::Sqrt(sum);
677 }
678 // compute variance in y (eprim2) and standard deviation in y (eprim)
679 Double_t contsum = cont/sum;
680 Double_t eprim2 = TMath::Abs(err2/sum - contsum*contsum);
681 Double_t eprim = TMath::Sqrt(eprim2);
682
683 if (p->fErrorMode == kERRORSPREADI) {
684 if (eprim != 0) return eprim/TMath::Sqrt(neff);
685 // in case content y is an integer (so each my has an error +/- 1/sqrt(12)
686 // when the std(y) is zero
687 return 1/TMath::Sqrt(12*neff);
688 }
689
690 // if approximate compute the sums (of w, wy and wy2) using all the bins
691 // when the variance in y is zero
692 Double_t test = 1;
693 if (err2 != 0 && neff < 5) test = eprim2*sum/err2;
694 //Int_t cellLimit = (p->GetDimension() == 3)?1000404:10404;
695 if (p->fgApproximate && (test < 1.e-4 || eprim2 <= 0)) {
696 Double_t stats[TH1::kNstat];
697 p->GetStats(stats);
698 Double_t ssum = stats[0];
699 // for 1D profile
700 int index = 4; // index in the stats array for 1D
701 if (p->GetDimension() == 2) index = 7; // for 2D
702 if (p->GetDimension() == 3) index = 11; // for 3D
703 Double_t scont = stats[index];
704 Double_t serr2 = stats[index+1];
705
706 // compute mean and variance in y
707 Double_t scontsum = scont/ssum; // global mean
708 Double_t seprim2 = TMath::Abs(serr2/ssum - scontsum*scontsum); // global variance
709 eprim = 2*TMath::Sqrt(seprim2); // global std (why factor of 2 ??)
710 sum = ssum;
711 }
712 sum = TMath::Abs(sum);
713
714 // case option "S" return standard deviation in y
715 if (p->fErrorMode == kERRORSPREAD) return eprim;
716
717 // default case : fErrorMode = kERRORMEAN
718 // return standard error on the mean of y
719 //if (neff == 0) std::cerr << "NEFF = 0 for bin " << bin << " " << eprim << " " << neff << " " << std::endl;
720 return eprim/TMath::Sqrt(neff);
721
722}
723
724
725template <typename T>
727// Set the number of entries in bin for the profile
728// In case the profile stores the sum of weight squares - set the sum of weight square to the number entries
729// (i.e. assume the entries have been filled all with a weight == 1 )
730
731 if (bin < 0 || bin >= p->fNcells) return;
732 p->fBinEntries.fArray[bin] = w;
733 if (p->fBinSumw2.fN) p->fBinSumw2.fArray[bin] = w;
734
735}
736
737#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:43
const Bool_t kFALSE
Definition: RtypesCore.h:90
double Double_t
Definition: RtypesCore.h:57
long long Long64_t
Definition: RtypesCore.h:71
const Bool_t kTRUE
Definition: RtypesCore.h:89
const char Option_t
Definition: RtypesCore.h:64
#define R__ASSERT(e)
Definition: TError.h:96
void Error(const char *location, const char *msgfmt,...)
void Warning(const char *location, const char *msgfmt,...)
float xmin
Definition: THbookFile.cxx:93
int nentries
Definition: THbookFile.cxx:89
float xmax
Definition: THbookFile.cxx:93
@ 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
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:728
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 from bin first to last.
Definition: TAxis.cxx:914
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition: TAxis.cxx:525
THashList * GetLabels() const
Definition: TAxis.h:117
Collection abstract base class.
Definition: TCollection.h:63
virtual void AddAll(const TCollection *col)
Add all objects from collection col to this collection.
virtual Bool_t IsEmpty() const
Definition: TCollection.h:186
The TH1 histogram class.
Definition: TH1.h:56
static Bool_t GetDefaultSumw2()
Return kTRUE if TH1::Sumw2 must be called when creating new histograms.
Definition: TH1.cxx:4294
@ kAllAxes
Definition: TH1.h:73
@ kNoAxis
NOTE: Must always be 0 !!!
Definition: TH1.h:69
@ kNstat
Definition: TH1.h:179
virtual Int_t BufferEmpty(Int_t action=0)
Fill histogram with all entries in the buffer.
Definition: TH1.cxx:1347
void Reset()
Definition: TCollection.h:252
A doubly linked list.
Definition: TList.h:44
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:821
virtual void AddFirst(TObject *obj)
Add object at the beginning of the list.
Definition: TList.cxx:99
Mother of all ROOT objects.
Definition: TObject.h:37
virtual UInt_t GetUniqueID() const
Return the unique object id.
Definition: TObject.cxx:375
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:131
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
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
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
Definition: test.py:1
static long int sum(long int i)
Definition: Factory.cxx:2275