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