Logo ROOT   master
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 
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  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 
420 template <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 
437  Double_t xmin, xmax;
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 
481 template <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 
498 template <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 
529 template <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 
595 template <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 
650 template <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 
661 template <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 
725 template <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
for(Int_t i=0;i< n;i++)
Definition: legend1.C:18
static void SetBinEntries(T *p, Int_t bin, Double_t w)
static void Scale(T *p, Double_t c1, Option_t *option)
static long int sum(long int i)
Definition: Factory.cxx:2274
float xmin
Definition: THbookFile.cxx:93
long long Long64_t
Definition: RtypesCore.h:69
virtual void SetLimits(Double_t xmin, Double_t xmax)
Definition: TAxis.h:154
static Double_t GetBinError(T *p, Int_t bin)
const char Option_t
Definition: RtypesCore.h:62
#define s0(x)
Definition: RSha256.hxx:90
return c1
Definition: legend1.C:41
double T(double x)
Definition: ChebyshevPol.h:34
static void SetErrorOption(T *p, Option_t *opt)
static void LabelsInflate(T *p, Option_t *)
virtual void AddFirst(TObject *obj)
Add object at the beginning of the list.
Definition: TList.cxx:97
virtual void AddAll(const TCollection *col)
Add all objects from collection col to this collection.
#define R__ASSERT(e)
Definition: TError.h:96
Definition: test.py:1
Basic string class.
Definition: TString.h:131
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
int Int_t
Definition: RtypesCore.h:41
Bool_t CanExtend() const
Definition: TAxis.h:82
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
void Reset()
Definition: TCollection.h:252
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition: TAxis.cxx:525
static void BuildArray(T *p)
Double_t GetXmin() const
Definition: TAxis.h:133
Double_t x[n]
Definition: legend1.C:17
static Bool_t GetDefaultSumw2()
Return kTRUE if TH1::Sumw2 must be called when creating new histograms.
Definition: TH1.cxx:4276
THashList * GetLabels() const
Definition: TAxis.h:117
static void LabelsDeflate(T *p, Option_t *)
static T * ExtendAxis(T *p, Double_t x, TAxis *axis)
TH1F * h1
Definition: legend1.C:5
void Error(const char *location, const char *msgfmt,...)
A doubly linked list.
Definition: TList.h:44
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
Class to manage histogram axis.
Definition: TAxis.h:30
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:819
Collection abstract base class.
Definition: TCollection.h:63
#define s1(x)
Definition: RSha256.hxx:91
float xmax
Definition: THbookFile.cxx:93
void Warning(const char *location, const char *msgfmt,...)
#define h(i)
Definition: RSha256.hxx:106
const Bool_t kFALSE
Definition: RtypesCore.h:88
static Bool_t Add(T *p, const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2=1)
virtual Bool_t IsEmpty() const
Definition: TCollection.h:186
return c2
Definition: legend2.C:14
static Double_t GetBinEffectiveEntries(T *p, Int_t bin)
int nentries
Definition: THbookFile.cxx:89
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
The TH1 histogram class.
Definition: TH1.h:56
Mother of all ROOT objects.
Definition: TObject.h:37
NOTE: Must always be 0 !!!
Definition: TH1.h:69
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual Int_t BufferEmpty(Int_t action=0)
Fill histogram with all entries in the buffer.
Definition: TH1.cxx:1345
static void Sumw2(T *p, Bool_t flag)
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Int_t GetNbins() const
Definition: TAxis.h:121
virtual void Set(Int_t nbins, Double_t xmin, Double_t xmax)
Initialize axis with fix bins.
Definition: TAxis.cxx:728
const Bool_t kTRUE
Definition: RtypesCore.h:87
Double_t GetXmax() const
Definition: TAxis.h:134
static Long64_t Merge(T *p, TCollection *list)