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 
32 public:
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 
70 template <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 
129 template <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 
142 template <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 
166 template <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 
429 template <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 
449  Double_t xmin, xmax;
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 (inclusing 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 
511 template <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 
528 template <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 
559 template <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 
625 template <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 
692 template <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 
703 template <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 
767 template <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
TProfileHelper::Scale
static void Scale(T *p, Double_t c1, Option_t *option)
Definition: TProfileHelper.h:512
TProfileHelper::Merge
static Long64_t Merge(T *p, TCollection *list)
Definition: TProfileHelper.h:167
TProfileHelper::SetErrorOption
static void SetErrorOption(T *p, Option_t *opt)
Definition: TProfileHelper.h:693
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:768
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
TProfileHelper::BuildArray
static void BuildArray(T *p)
Definition: TProfileHelper.h:130
TList::AddFirst
virtual void AddFirst(TObject *obj)
Add object at the beginning of the list.
Definition: TList.cxx:100
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:704
xmax
float xmax
Definition: THbookFile.cxx:95
sum
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345
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:626
TH1::GetDefaultSumw2
static Bool_t GetDefaultSumw2()
Return kTRUE if TH1::Sumw2 must be called when creating new histograms.
Definition: TH1.cxx:4374
Int_t
int Int_t
Definition: RtypesCore.h:45
TProfileHelper::Sumw2
static void Sumw2(T *p, Bool_t flag)
Definition: TProfileHelper.h:529
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
TH1::kAllAxes
@ kAllAxes
Definition: TH1.h:75
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
TH1::kNoAxis
@ kNoAxis
NOTE: Must always be 0 !!!
Definition: TH1.h:71
TProfileHelper::Add
static Bool_t Add(T *p, const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2=1)
Definition: TProfileHelper.h:71
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
TH1Merger
Definition: TH1Merger.h:20
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:430
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:1398
TMath::IsNaN
Bool_t IsNaN(Double_t x)
Definition: TMath.h:892
TProfileHelper::LabelsDeflate
static void LabelsDeflate(T *p, Option_t *)
Definition: TProfileHelper.h:560
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
s1
#define s1(x)
Definition: RSha256.hxx:91
gDebug
Int_t gDebug
Definition: TROOT.cxx:590
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:143
THashList.h
Double_t
double Double_t
Definition: RtypesCore.h:59
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:822
Info
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition: TError.cxx:220
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:30
c2
return c2
Definition: legend2.C:14
TH1Merger.h
kERRORSPREADG
@ kERRORSPREADG
Definition: TProfile.h:28
TIter
Definition: TCollection.h:233
TCollection
Collection abstract base class.
Definition: TCollection.h:63
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
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
TH1::kNstat
@ kNstat
Definition: TH1.h:183
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