source: trunk/MagicSoft/Mars/mhbase/MH.cc@ 9464

Last change on this file since 9464 was 9414, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 56.5 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: MH.cc,v 1.49 2009-03-30 08:06:41 tbretz Exp $
3! --------------------------------------------------------------------------
4!
5! *
6! * This file is part of MARS, the MAGIC Analysis and Reconstruction
7! * Software. It is distributed to you in the hope that it can be a useful
8! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
9! * It is distributed WITHOUT ANY WARRANTY.
10! *
11! * Permission to use, copy, modify and distribute this software and its
12! * documentation for any purpose is hereby granted without fee,
13! * provided that the above copyright notice appear in all copies and
14! * that both that copyright notice and this permission notice appear
15! * in supporting documentation. It is provided "as is" without express
16! * or implied warranty.
17! *
18!
19!
20! Author(s): Thomas Bretz 07/2001 <mailto:tbretz@astro.uni-wuerzburg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2008
23!
24!
25\* ======================================================================== */
26
27//////////////////////////////////////////////////////////////////////////////
28// //
29// MH //
30// //
31// This is a base tasks for mars histograms. It defines a common interface //
32// for filling the histograms with events (MH::Fill) which is used by a //
33// common 'filler' And a SetupFill member function which may be used //
34// by MFillH. The idea is: //
35// 1) If your Histogram can become filled by one single container //
36// (like MHHillas) you overload MH::Fill and it gets called with //
37// a pointer to the container with which it should be filled. //
38// //
39// 2) You histogram needs several containers to get filled. Than you //
40// have to overload MH::SetupFill and get the necessary objects from //
41// the parameter list. Use this objects in Fill to fill your //
42// histogram. //
43// //
44// If you want to create your own histogram class the new class must be //
45// derived from MH (instead of the base MParContainer) and you must //
46// the fill function of MH. This is the function which is called to fill //
47// the histogram(s) by the data of a corresponding parameter container. //
48// //
49// Remark: the static member function (eg MakeDefCanvas) can be called //
50// from everywhere using: MH::MakeDefCanvas(...) //
51// //
52//////////////////////////////////////////////////////////////////////////////
53
54#include "MH.h"
55
56#include <TH1.h>
57#include <TH2.h>
58#include <TH3.h>
59#include <TColor.h>
60#include <TMath.h>
61#include <TClass.h>
62#include <TStyle.h> // TStyle::GetScreenFactor
63#include <TGaxis.h>
64#include <TCanvas.h>
65#include <TLegend.h>
66#include <TPaveStats.h>
67#include <TBaseClass.h>
68#include <THashList.h>
69#if ROOT_VERSION_CODE > ROOT_VERSION(3,04,01)
70#include <THLimitsFinder.h>
71#endif
72#include <TProfile2D.h> // root > 5.18
73
74#include "MLog.h"
75#include "MLogManip.h"
76
77#include "MString.h"
78
79#include "MParList.h"
80#include "MParContainer.h"
81
82#include "MBinning.h"
83
84#include "MArrayD.h"
85#include "MArrayF.h"
86
87ClassImp(MH);
88
89using namespace std;
90
91// --------------------------------------------------------------------------
92//
93// Default Constructor. It sets name and title only. Typically you won't
94// need to change this.
95//
96MH::MH(const char *name, const char *title)
97 : fSerialNumber(0), fNumExecutions(0)
98
99{
100 //
101 // set the name and title of this object
102 //
103 fName = name ? name : "MH";
104 fTitle = title ? title : "Base class for Mars histograms";
105}
106
107// --------------------------------------------------------------------------
108//
109// If you want to use the automatic filling of your derived class you
110// must overload this function. If it is not overloaded you cannot use
111// FillH with this class. The argument is a pointer to a container
112// in your paremeter list which is specified in the MFillH constructor.
113// If you are not going to use it you should at least add
114// Bool_t MH::Fill(const MParContainer *) { return kTRUE; }
115// to your class definition.
116//
117Int_t MH::Fill(const MParContainer *par, const Stat_t w)
118{
119 *fLog << warn << GetDescriptor() << ": Fill not overloaded! Can't be used!" << endl;
120 return kERROR;
121}
122
123// --------------------------------------------------------------------------
124//
125// This virtual function is ment as a generalized interface to retrieve
126// a pointer to a root histogram from the MH-derived class.
127//
128TH1 *MH::GetHistByName(const TString name) const
129{
130 *fLog << warn << GetDescriptor() << ": GetHistByName not overloaded! Can't be used!" << endl;
131 return NULL;
132}
133
134// --------------------------------------------------------------------------
135//
136// This is a function which should replace the creation of default
137// canvases like root does. Because this is inconvinient in some aspects.
138// need to change this.
139// You can specify a name for the default canvas and a title. Also
140// width and height can be given.
141// MakeDefCanvas looks for a canvas with the given name. If now name is
142// given the DefCanvasName of root is used. If no such canvas is existing
143// it is created and returned. If such a canvas already exists a new canvas
144// with a name plus anumber is created (the number is calculated by the
145// number of all existing canvases plus one)
146//
147// Normally the canvas size is scaled with gStyle->GetScreenFactor() so
148// that on all screens it looks like the same part of the screen.
149// To suppress this scaling use usescreenfactor=kFALSE. In this case
150// you specify directly the size of the embedded pad.
151//
152TCanvas *MH::MakeDefCanvas(TString name, const char *title,
153 UInt_t w, UInt_t h, Bool_t usescreenfactor)
154{
155 const TList *list = (TList*)gROOT->GetListOfCanvases();
156
157 if (name.IsNull())
158 name = gROOT->GetDefCanvasName();
159
160 if (list->FindObject(name))
161 name += MString::Format(" <%d>", list->GetSize()+1);
162
163 if (!usescreenfactor)
164 {
165 const Float_t cx = gStyle->GetScreenFactor();
166 w += 4;
167 h += 28;
168 w = (int)(w/cx+1);
169 h = (int)(h/cx+1);
170 }
171
172 return new TCanvas(name, title, w, h);
173}
174
175// --------------------------------------------------------------------------
176//
177// This function works like MakeDefCanvas(name, title, w, h) but name
178// and title are retrieved from the given TObject.
179//
180// Normally the canvas size is scaled with gStyle->GetScreenFactor() so
181// that on all screens it looks like the same part of the screen.
182// To suppress this scaling use usescreenfactor=kFALSE. In this case
183// you specify directly the size of the embedded pad.
184//
185TCanvas *MH::MakeDefCanvas(const TObject *obj,
186 UInt_t w, UInt_t h, Bool_t usescreenfactor)
187{
188 if (!usescreenfactor)
189 {
190 const Float_t cx = gStyle->GetScreenFactor();
191 w += 4;
192 h += 28;
193 h = (int)(h/cx+1);
194 w = (int)(w/cx+1);
195 }
196
197 return MakeDefCanvas(obj->GetName(), obj->GetTitle(), w, h);
198}
199
200// --------------------------------------------------------------------------
201//
202// Search in gPad for all objects with the name name and remove all of them
203// (TList::Remove)
204//
205void MH::RemoveFromPad(const char *name)
206{
207 if (!gPad)
208 return;
209
210 TList *list = gPad->GetListOfPrimitives();
211 if (!list)
212 return;
213
214 TObject *obj = 0;
215 while ((obj = gPad->FindObject(name)))
216 list->Remove(obj);
217}
218
219// --------------------------------------------------------------------------
220//
221// If labels are set for this axis the correct MBinning corresponding
222// to the existing label range is returned (this is necessary to
223// maintain the correct number of bins in the histogram)
224// otherwise the given binning is returned.
225//
226MBinning MH::GetBinningForLabels(TAxis &x, const MBinning *bins)
227{
228 if (!x.GetLabels())
229 return *bins;
230
231 const Int_t n = TMath::Max(x.GetLabels()->GetEntries(), 1);
232 return MBinning(n, 0, n);
233}
234
235// --------------------------------------------------------------------------
236//
237// If Labels are set this function deletes the fXbins Array from
238// the axis (which makes the axis a variable bin-size axis)
239// and sets the Nbins, Xmin and Xmax according to the number of labels.
240//
241void MH::RestoreBinningForLabels(TAxis &x)
242{
243 if (!x.GetLabels())
244 return;
245
246 const Int_t n = TMath::Max(x.GetLabels()->GetEntries(), 1);
247 x.Set(n, 0, n);
248
249 const_cast<TArrayD*>(x.GetXbins())->Set(0);
250}
251
252// --------------------------------------------------------------------------
253//
254// Applies a given binning to a 1D-histogram. In case the axis has labels
255// (e.g. GetXaxis()->GetLabels()) the binning is set according to the
256// labels.
257//
258void MH::SetBinning(TH1 *h, const MBinning *binsx)
259{
260 //
261 // Another strange behaviour: TAxis::Set deletes the axis title!
262 //
263 TAxis &x = *h->GetXaxis();
264
265#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
266 TString xtitle = x.GetTitle();
267#endif
268
269#if ROOT_VERSION_CODE < ROOT_VERSION(5,12,00)
270 // All this is reset by TAxis::Set
271 const TAttAxis att(x);
272 const Bool_t tm(x.GetTimeDisplay());
273 const TString tf(x.GetTimeFormat());
274
275 //
276 // This is a necessary workaround if one wants to set
277 // non-equidistant bins after the initialization
278 // TH1D::fNcells must be set correctly.
279 //
280 h->SetBins(binsx->GetNumBins(), 0, 1);
281
282 //
283 // Set the binning of the current histogram to the binning
284 // in one of the two given histograms
285 //
286 x.Set(binsx->GetNumBins(), binsx->GetEdges());
287
288 // All this is reset by TAxis::Set
289 att.Copy(x);
290 x.SetTimeDisplay(tm);
291 x.SetTimeFormat(tf);
292#else
293 if (!x.GetLabels())
294 h->SetBins(binsx->GetNumBins(), binsx->GetEdges());
295#endif
296
297
298#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
299 x.SetTitle(xtitle);
300#endif
301}
302
303// --------------------------------------------------------------------------
304//
305// Applies given binnings to the two axis of a 2D-histogram.
306// In case the axis has labels (e.g. GetXaxis()->GetLabels())
307// the binning is set according to the labels.
308//
309void MH::SetBinning(TH2 *h, const MBinning *binsx, const MBinning *binsy)
310{
311 TAxis &x = *h->GetXaxis();
312 TAxis &y = *h->GetYaxis();
313
314 const MBinning bx(GetBinningForLabels(x, binsx));
315 const MBinning by(GetBinningForLabels(y, binsy));
316
317 //
318 // Another strange behaviour: TAxis::Set deletes the axis title!
319 //
320#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
321 TString xtitle = x.GetTitle();
322 TString ytitle = y.GetTitle();
323#endif
324
325#if ROOT_VERSION_CODE < ROOT_VERSION(5,12,00)
326 // All this is reset by TAxis::Set
327 const TAttAxis attx(x);
328 const TAttAxis atty(y);
329 const Bool_t tmx(x.GetTimeDisplay());
330 const Bool_t tmy(y.GetTimeDisplay());
331 const TString tfx(x.GetTimeFormat());
332 const TString tfy(y.GetTimeFormat());
333
334 //
335 // This is a necessary workaround if one wants to set
336 // non-equidistant bins after the initialization
337 // TH1D::fNcells must be set correctly.
338 //
339 h->SetBins(bx.GetNumBins(), 0, 1,
340 by.GetNumBins(), 0, 1);
341
342 //
343 // Set the binning of the current histogram to the binning
344 // in one of the two given histograms
345 //
346 x.Set(bx.GetNumBins(), bx.GetEdges());
347 y.Set(by.GetNumBins(), by.GetEdges());
348
349 // All this is reset by TAxis::Set
350 attx.Copy(x);
351 atty.Copy(y);
352 x.SetTimeDisplay(tmx);
353 y.SetTimeDisplay(tmy);
354 x.SetTimeFormat(tfx);
355 y.SetTimeFormat(tfy);
356#else
357 if (h->InheritsFrom(TProfile2D::Class()))
358 {
359 h->SetBins(bx.GetNumBins(), 0, 1,
360 by.GetNumBins(), 0, 1);
361
362 h->SetBinsLength();
363
364 x.Set(bx.GetNumBins(), bx.GetEdges());
365 y.Set(by.GetNumBins(), by.GetEdges());
366 }
367 else
368 h->SetBins(bx.GetNumBins(), bx.GetEdges(),
369 by.GetNumBins(), by.GetEdges());
370#endif
371
372 RestoreBinningForLabels(x);
373 RestoreBinningForLabels(y);
374
375#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
376 x.SetTitle(xtitle);
377 y.SetTitle(ytitle);
378#endif
379}
380
381// --------------------------------------------------------------------------
382//
383// Applies given binnings to the three axis of a 3D-histogram
384// In case the axis has labels (e.g. GetXaxis()->GetLabels())
385// the binning is set according to the labels.
386//
387void MH::SetBinning(TH3 *h, const MBinning *binsx, const MBinning *binsy, const MBinning *binsz)
388{
389 //
390 // Another strange behaviour: TAxis::Set deletes the axis title!
391 //
392 TAxis &x = *h->GetXaxis();
393 TAxis &y = *h->GetYaxis();
394 TAxis &z = *h->GetZaxis();
395
396 const MBinning bx(GetBinningForLabels(x, binsx));
397 const MBinning by(GetBinningForLabels(y, binsy));
398 const MBinning bz(GetBinningForLabels(z, binsz));
399
400#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
401 TString xtitle = x.GetTitle();
402 TString ytitle = y.GetTitle();
403 TString ztitle = z.GetTitle();
404#endif
405
406#if ROOT_VERSION_CODE < ROOT_VERSION(5,12,00)
407 // All this is reset by TAxis::Set
408 const TAttAxis attx(x);
409 const TAttAxis atty(y);
410 const TAttAxis attz(z);
411 const Bool_t tmx(x.GetTimeDisplay());
412 const Bool_t tmy(y.GetTimeDisplay());
413 const Bool_t tmz(z.GetTimeDisplay());
414 const TString tfx(x.GetTimeFormat());
415 const TString tfy(y.GetTimeFormat());
416 const TString tfz(z.GetTimeFormat());
417#endif
418
419 //
420 // This is a necessary workaround if one wants to set
421 // non-equidistant bins after the initialization
422 // TH1D::fNcells must be set correctly.
423 //
424 h->SetBins(bx.GetNumBins(), 0, 1,
425 by.GetNumBins(), 0, 1,
426 bz.GetNumBins(), 0, 1);
427
428 //
429 // Set the binning of the current histogram to the binning
430 // in one of the two given histograms
431 //
432 x.Set(bx.GetNumBins(), bx.GetEdges());
433 y.Set(by.GetNumBins(), by.GetEdges());
434 z.Set(bz.GetNumBins(), bz.GetEdges());
435
436 RestoreBinningForLabels(x);
437 RestoreBinningForLabels(y);
438 RestoreBinningForLabels(z);
439
440#if ROOT_VERSION_CODE < ROOT_VERSION(5,12,00)
441 // All this is reset by TAxis::Set
442 attx.Copy(x);
443 atty.Copy(y);
444 attz.Copy(z);
445 x.SetTimeDisplay(tmx);
446 y.SetTimeDisplay(tmy);
447 z.SetTimeDisplay(tmz);
448 x.SetTimeFormat(tfx);
449 y.SetTimeFormat(tfy);
450 z.SetTimeFormat(tfz);
451#endif
452
453#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
454 x.SetTitle(xtitle);
455 y.SetTitle(ytitle);
456 z.SetTitle(ztitle);
457#endif
458}
459
460// --------------------------------------------------------------------------
461//
462// Applies given binning (the n+1 edges) to the axis of a 1D-histogram
463//
464void MH::SetBinning(TH1 *h, const TArrayD &binsx)
465{
466 MBinning bx;
467 bx.SetEdges(binsx);
468 SetBinning(h, &bx);
469}
470
471// --------------------------------------------------------------------------
472//
473// Applies given binning (the n+1 edges) to the two axis of a
474// 2D-histogram
475//
476void MH::SetBinning(TH2 *h, const TArrayD &binsx, const TArrayD &binsy)
477{
478 MBinning bx;
479 MBinning by;
480 bx.SetEdges(binsx);
481 by.SetEdges(binsy);
482 SetBinning(h, &bx, &by);
483}
484
485// --------------------------------------------------------------------------
486//
487// Applies given binning (the n+1 edges) to the three axis of a
488// 3D-histogram
489//
490void MH::SetBinning(TH3 *h, const TArrayD &binsx, const TArrayD &binsy, const TArrayD &binsz)
491{
492 MBinning bx;
493 MBinning by;
494 MBinning bz;
495 bx.SetEdges(binsx);
496 by.SetEdges(binsy);
497 bz.SetEdges(binsz);
498 SetBinning(h, &bx, &by, &bz);
499}
500
501// --------------------------------------------------------------------------
502//
503// Applies the binning of a TAxis (eg from a root histogram) to the axis
504// of a 1D-histogram
505//
506void MH::SetBinning(TH1 *h, const TAxis *binsx)
507{
508 const Int_t nx = binsx->GetNbins();
509
510 TArrayD bx(nx+1);
511 for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
512 bx[nx] = binsx->GetXmax();
513
514 SetBinning(h, bx);
515}
516
517// --------------------------------------------------------------------------
518//
519// Applies the binnings of the TAxis' (eg from a root histogram) to the
520// two axis' of a 2D-histogram
521//
522void MH::SetBinning(TH2 *h, const TAxis *binsx, const TAxis *binsy)
523{
524 const Int_t nx = binsx->GetNbins();
525 const Int_t ny = binsy->GetNbins();
526
527 TArrayD bx(nx+1);
528 TArrayD by(ny+1);
529 for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
530 for (int i=0; i<ny; i++) by[i] = binsy->GetBinLowEdge(i+1);
531 bx[nx] = binsx->GetXmax();
532 by[ny] = binsy->GetXmax();
533
534 SetBinning(h, bx, by);
535}
536
537// --------------------------------------------------------------------------
538//
539// Applies the binnings of the TAxis' (eg from a root histogram) to the
540// three axis' of a 3D-histogram
541//
542void MH::SetBinning(TH3 *h, const TAxis *binsx, const TAxis *binsy, const TAxis *binsz)
543{
544 const Int_t nx = binsx->GetNbins();
545 const Int_t ny = binsy->GetNbins();
546 const Int_t nz = binsz->GetNbins();
547
548 TArrayD bx(nx+1);
549 TArrayD by(ny+1);
550 TArrayD bz(nz+1);
551 for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
552 for (int i=0; i<ny; i++) by[i] = binsy->GetBinLowEdge(i+1);
553 for (int i=0; i<nz; i++) bz[i] = binsz->GetBinLowEdge(i+1);
554 bx[nx] = binsx->GetXmax();
555 by[ny] = binsy->GetXmax();
556 bz[nz] = binsz->GetXmax();
557
558 SetBinning(h, bx, by, bz);
559}
560
561// --------------------------------------------------------------------------
562//
563// Applies the binnings of one root-histogram x to another one h
564// Both histograms must be of the same type: TH1, TH2 or TH3
565//
566void MH::SetBinning(TH1 *h, const TH1 *x)
567{
568 if (h->InheritsFrom(TH3::Class()) && x->InheritsFrom(TH3::Class()))
569 {
570 SetBinning((TH3*)h, x->GetXaxis(), x->GetYaxis(), x->GetZaxis());
571 return;
572 }
573 if (h->InheritsFrom(TH3::Class()) || x->InheritsFrom(TH3::Class()))
574 return;
575 if (h->InheritsFrom(TH2::Class()) && x->InheritsFrom(TH2::Class()))
576 {
577 SetBinning((TH2*)h, x->GetXaxis(), x->GetYaxis());
578 return;
579 }
580 if (h->InheritsFrom(TH2::Class()) || x->InheritsFrom(TH2::Class()))
581 return;
582 if (h->InheritsFrom(TH1::Class()) && x->InheritsFrom(TH1::Class()))
583 {
584 SetBinning(h, x->GetXaxis());
585 return;
586 }
587}
588
589void MH::RemoveFirstBin(TH1 &h)
590{
591 if (h.InheritsFrom(TH2::Class()) || h.InheritsFrom(TH3::Class()))
592 return;
593
594 const Bool_t haserr = h.GetSumw2N()>0;
595
596 const Int_t n0 = h.GetNbinsX();
597 if (n0<2)
598 return;
599
600 TArrayD val(n0-1);
601 TArrayD er(haserr ? n0-1 : 0);
602 for (int i=1; i<n0; i++)
603 {
604 val[i-1] = h.GetBinContent(i+1);
605 if (haserr)
606 er[i-1] = h.GetBinError(i+1);
607 }
608
609 MBinning bins;
610 bins.SetEdges(h, 'x');
611 bins.RemoveFirstEdge();
612 bins.Apply(h);
613
614 h.Reset();
615
616 for (int i=1; i<n0; i++)
617 {
618 h.SetBinContent(i, val[i-1]);
619 if (haserr)
620 h.SetBinError(i, er[i-1]);
621 }
622}
623
624// --------------------------------------------------------------------------
625//
626// Multiplies all entries in a TArrayD by a float f
627//
628void MH::ScaleArray(TArrayD &bins, Double_t f)
629{
630 for (int i=0; i<bins.GetSize(); i++)
631 bins[i] *= f;
632}
633
634// --------------------------------------------------------------------------
635//
636// Scales the binning of a TAxis by a float f
637//
638TArrayD MH::ScaleAxis(TAxis &axe, Double_t f)
639{
640 TArrayD arr(axe.GetNbins()+1);
641
642 for (int i=1; i<=axe.GetNbins()+1; i++)
643 arr[i-1] = axe.GetBinLowEdge(i);
644
645 ScaleArray(arr, f);
646
647 return arr;
648}
649
650// --------------------------------------------------------------------------
651//
652// Scales the binning of one, two or three axis of a histogram by a float f
653//
654void MH::ScaleAxis(TH1 *h, Double_t fx, Double_t fy, Double_t fz)
655{
656 if (h->InheritsFrom(TH3::Class()))
657 {
658 SetBinning((TH3*)h,
659 ScaleAxis(*h->GetXaxis(), fx),
660 ScaleAxis(*h->GetYaxis(), fy),
661 ScaleAxis(*h->GetZaxis(), fz));
662 return;
663 }
664
665 if (h->InheritsFrom(TH2::Class()))
666 {
667 SetBinning((TH2*)h,
668 ScaleAxis(*h->GetXaxis(), fx),
669 ScaleAxis(*h->GetYaxis(), fy));
670 return;
671 }
672
673 if (h->InheritsFrom(TH1::Class()))
674 SetBinning(h, ScaleAxis(*h->GetXaxis(), fx));
675}
676
677// --------------------------------------------------------------------------
678//
679// Tries to find a MBinning container with the name "Binning"+name
680// in the given parameter list. If it was found it is applied to the
681// given histogram. This is only valid for 1D-histograms.
682// If the binning is found, but it IsDefault() kTRUE is returned, but
683// no binning is applied.
684//
685Bool_t MH::ApplyBinning(const MParList &plist, TString name, TH1 *h)
686{
687 if (h->InheritsFrom(TH2::Class()) || h->InheritsFrom(TH3::Class()))
688 {
689 gLog << warn << "MH::ApplyBinning: '" << h->GetName() << "' is not a basic TH1 object... no binning applied." << endl;
690 return kFALSE;
691 }
692
693 const MBinning *bins = (MBinning*)plist.FindObject("Binning"+name);
694 if (!bins)
695 {
696 gLog << inf << "Object 'Binning" << name << "' [MBinning] not found... no binning applied." << endl;
697 return kFALSE;
698 }
699
700 if (bins->IsDefault())
701 {
702 gLog << inf << "Object 'Binning" << name << "' [MBinning] is default... binning unchanged." << endl;
703 return kTRUE;
704 }
705
706 SetBinning(h, bins);
707 return kTRUE;
708}
709
710Bool_t MH::ApplyBinning(const MParList &plist, TString x, TString y, TH2 *h)
711{
712 const MBinning *binsx = (MBinning*)plist.FindObject("Binning"+x);
713 const MBinning *binsy = (MBinning*)plist.FindObject("Binning"+y);
714
715 if (!binsx && !binsy)
716 {
717 gLog << inf << "Neither 'Binning" << x << "' nor 'Binning" << y;
718 gLog << "' [MBinning] found... no binning applied." << endl;
719 return kFALSE;
720 }
721
722 if (!binsx)
723 gLog << inf << "Object 'Binning" << x << "' [MBinning] not found... binning unchanged." << endl;
724 if (!binsy)
725 gLog << inf << "Object 'Binning" << y << "' [MBinning] not found... binning unchanged." << endl;
726
727 if (binsx && binsx->IsDefault())
728 {
729 gLog << inf << "Object 'Binning" << x << "' [MBinning] is default... binning unchanged." << endl;
730 binsx = 0;
731 }
732
733 if (binsy && binsy->IsDefault())
734 {
735 gLog << inf << "Object 'Binning" << y << "' [MBinning] is default... binning unchanged." << endl;
736 binsy = 0;
737 }
738
739 MBinning binsxx, binsyy;
740 binsxx.SetEdges(*h, 'x');
741 binsyy.SetEdges(*h, 'y');
742
743 SetBinning(h, binsx?binsx:&binsxx, binsy?binsy:&binsyy);
744
745 return kTRUE;
746}
747
748Bool_t MH::ApplyBinning(const MParList &plist, TString x, TString y, TString z, TH3 *h)
749{
750 const MBinning *binsx = (MBinning*)plist.FindObject("Binning"+x);
751 const MBinning *binsy = (MBinning*)plist.FindObject("Binning"+y);
752 const MBinning *binsz = (MBinning*)plist.FindObject("Binning"+z);
753
754 if (!binsx && !binsy && !binsz)
755 {
756 gLog << inf << "Neither 'Binning" << x << "', 'Binning" << y;
757 gLog << "' nor 'Binning" << z << "' [MBinning] found... ";
758 gLog << "no binning applied." << endl;
759 return kFALSE;
760 }
761
762 if (!binsx)
763 gLog << inf << "Object 'Binning" << x << "' [MBinning] not found... binning unchanged." << endl;
764 if (!binsy)
765 gLog << inf << "Object 'Binning" << y << "' [MBinning] not found... binning unchanged." << endl;
766 if (!binsz)
767 gLog << inf << "Object 'Binning" << z << "' [MBinning] not found... binning unchanged." << endl;
768
769 if (binsx && binsx->IsDefault())
770 {
771 gLog << inf << "Object 'Binning" << x << "' [MBinning] is default... binning unchanged." << endl;
772 binsx = 0;
773 }
774
775 if (binsy && binsy->IsDefault())
776 {
777 gLog << inf << "Object 'Binning" << y << "' [MBinning] is default... binning unchanged." << endl;
778 binsy = 0;
779 }
780
781 if (binsz && binsz->IsDefault())
782 {
783 gLog << inf << "Object 'Binning" << z << "' [MBinning] is default... binning unchanged." << endl;
784 binsy = 0;
785 }
786
787 MBinning binsxx, binsyy, binszz;
788 binsxx.SetEdges(*h, 'x');
789 binsyy.SetEdges(*h, 'y');
790 binszz.SetEdges(*h, 'z');
791
792 SetBinning(h, binsx?binsx:&binsxx, binsy?binsy:&binsyy, binsz?binsz:&binszz);
793
794 return kTRUE;
795}
796
797void MH::FindGoodLimits(Int_t nbins, Int_t &newbins, Double_t &xmin, Double_t &xmax, Bool_t isInteger)
798{
799#if ROOT_VERSION_CODE > ROOT_VERSION(3,04,01)
800 THLimitsFinder::OptimizeLimits(nbins, newbins, xmin, xmax, isInteger);
801#else
802//*-*-*-*-*-*-*-*-*Find reasonable bin values*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
803//*-* ==========================
804
805 Double_t dx = 0.1*(xmax-xmin);
806 Double_t umin = xmin - dx;
807 Double_t umax = xmax + dx;
808
809 if (umin < 0 && xmin >= 0)
810 umin = 0;
811
812 if (umax > 0 && xmax <= 0)
813 umax = 0;
814
815 Double_t binlow =0;
816 Double_t binhigh =0;
817 Double_t binwidth=0;
818
819 TGaxis::Optimize(umin, umax, nbins, binlow, binhigh, nbins, binwidth, "");
820
821 if (binwidth <= 0 || binwidth > 1.e+39)
822 {
823 xmin = -1;
824 xmax = 1;
825 }
826 else
827 {
828 xmin = binlow;
829 xmax = binhigh;
830 }
831
832 if (isInteger)
833 {
834 Int_t ixmin = (Int_t)xmin;
835 Int_t ixmax = (Int_t)xmax;
836 Double_t dxmin = (Double_t)ixmin;
837 Double_t dxmax = (Double_t)ixmax;
838
839 xmin = xmin<0 && xmin!=dxmin ? dxmin - 1 : dxmin;
840 xmax = xmax>0 && xmax!=dxmax ? dxmax + 1 : dxmax;
841
842 if (xmin>=xmax)
843 xmax = xmin+1;
844
845 Int_t bw = 1 + (Int_t)((xmax-xmin)/nbins);
846
847 nbins = (Int_t)((xmax-xmin)/bw);
848
849 if (xmin+nbins*bw < xmax)
850 {
851 nbins++;
852 xmax = xmin +nbins*bw;
853 }
854 }
855
856 newbins = nbins;
857#endif
858}
859
860// --------------------------------------------------------------------------
861//
862// Returns the lowest entry in a histogram which is greater than gt (eg >0)
863//
864Double_t MH::GetMinimumGT(const TH1 &h, Double_t gt)
865{
866 Double_t min = FLT_MAX;
867
868 const Int_t nx = h.GetXaxis()->GetNbins();
869 const Int_t ny = h.GetYaxis()->GetNbins();
870 const Int_t nz = h.GetZaxis()->GetNbins();
871
872 for (int iz=1; iz<=nz; iz++)
873 for (int iy=1; iy<=ny; iy++)
874 for (int ix=1; ix<=nx; ix++)
875 {
876 const Double_t v = h.GetBinContent(h.GetBin(ix, iy, iz));
877 if (gt<v && v<min)
878 min = v;
879 }
880 return min;
881}
882
883// --------------------------------------------------------------------------
884//
885// Returns the bin center in a logarithmic scale. If the given bin
886// number is <1 it is set to 1. If it is =GetNbins() it is set to
887// GetNbins()
888//
889Double_t MH::GetBinCenterLog(const TAxis &axe, Int_t nbin)
890{
891 if (nbin>axe.GetNbins())
892 nbin = axe.GetNbins();
893
894 if (nbin<1)
895 nbin = 1;
896
897 const Double_t lo = axe.GetBinLowEdge(nbin);
898 const Double_t hi = axe.GetBinUpEdge(nbin);
899
900 const Double_t val = log10(lo) + log10(hi);
901
902 return pow(10, val/2);
903}
904
905// --------------------------------------------------------------------------
906//
907// Draws a copy of the two given histograms. Uses title as the pad title.
908// Also layout the two statistic boxes and a legend.
909//
910void MH::DrawSameCopy(const TH1 &hist1, const TH1 &hist2, const TString title)
911{
912 //
913 // Draw first histogram
914 //
915 TH1 *h1 = hist1.DrawCopy();
916 gPad->SetBorderMode(0);
917 gPad->Update();
918
919 // FIXME: Also align max/min with set Maximum/Minimum
920 const Double_t maxbin1 = hist1.GetBinContent(hist1.GetMaximumBin());
921 const Double_t maxbin2 = hist2.GetBinContent(hist2.GetMaximumBin());
922 const Double_t minbin1 = hist1.GetBinContent(hist1.GetMinimumBin());
923 const Double_t minbin2 = hist2.GetBinContent(hist2.GetMinimumBin());
924
925 const Double_t max = TMath::Max(maxbin1, maxbin2);
926 const Double_t min = TMath::Min(minbin1, minbin2);
927
928 h1->SetMaximum(max>0?max*1.05:max*0.95);
929 h1->SetMinimum(max>0?min*0.95:min*1.05);
930
931 TPaveText *t = (TPaveText*)gPad->FindObject("title");
932 if (t)
933 {
934 t->SetName((TString)"MHTitle"); // rename object
935 t->Clear(); // clear old lines
936 t->AddText((TString)" "+title+" "); // add the new title
937 t->SetBit(kCanDelete); // make sure object is deleted
938
939 //
940 // FIXME: This is a stupid workaround to hide the redrawn
941 // (see THistPainter::PaintTitle) title
942 //
943 gPad->Modified(); // indicates a change
944 gPad->Update(); // recreates the original title
945 t->Pop(); // bring our title on top
946 }
947
948 //
949 // Rename first statistics box
950 //
951 TPaveStats *s1 = dynamic_cast<TPaveStats*>(gPad->FindObject("stats"));
952 if (!s1)
953 s1 = dynamic_cast<TPaveStats*>(hist1.GetListOfFunctions()->FindObject("stats"));
954 else
955 s1->SetName((TString)"Stat"+hist1.GetTitle());
956
957 if (s1 && s1->GetX2NDC()>0.95)
958 {
959 const Double_t x1 = s1->GetX1NDC()-0.01;
960 s1->SetX1NDC(x1-(s1->GetX2NDC()-s1->GetX1NDC()));
961 s1->SetX2NDC(x1);
962 }
963
964 //
965 // Draw second histogram
966 //
967 TH1 *h2 = hist2.DrawCopy("sames");
968 gPad->Update();
969
970 //
971 // Draw Legend
972 //
973 TPaveStats *s2 = dynamic_cast<TPaveStats*>(gPad->FindObject("stats"));
974 if (!s2)
975 s2 = dynamic_cast<TPaveStats*>(hist2.GetListOfFunctions()->FindObject("stats"));
976
977 if (s2)
978 {
979 TLegend &l = *new TLegend(s2->GetX1NDC(),
980 s2->GetY1NDC()-0.015-(s2->GetY2NDC()-s2->GetY1NDC())/2,
981 s2->GetX2NDC(),
982 s2->GetY1NDC()-0.01
983 );
984 l.AddEntry(h1, h1->GetTitle());
985 l.AddEntry(h2, h2->GetTitle());
986 l.SetTextSize(s2->GetTextSize());
987 l.SetTextFont(s2->GetTextFont());
988 l.SetBorderSize(s2->GetBorderSize());
989 l.SetBit(kCanDelete);
990 l.Draw();
991 }
992}
993
994// --------------------------------------------------------------------------
995//
996// Draws the two given histograms. Uses title as the pad title.
997// Also layout the two statistic boxes and a legend.
998//
999void MH::DrawSame(TH1 &hist1, TH1 &hist2, const TString title, Bool_t same)
1000{
1001 //
1002 // Draw first histogram
1003 //
1004 hist1.Draw(same?"same":"");
1005 gPad->SetBorderMode(0);
1006 gPad->Update();
1007
1008 if (hist1.GetEntries()>0 && hist2.GetEntries()>0)
1009 {
1010 const Double_t maxbin1 = hist1.GetBinContent(hist1.GetMaximumBin());
1011 const Double_t maxbin2 = hist2.GetBinContent(hist2.GetMaximumBin());
1012 const Double_t minbin1 = hist1.GetBinContent(hist1.GetMinimumBin());
1013 const Double_t minbin2 = hist2.GetBinContent(hist2.GetMinimumBin());
1014
1015 const Double_t max = TMath::Max(maxbin1, maxbin2);
1016 const Double_t min = TMath::Min(minbin1, minbin2);
1017
1018 if (max!=min)
1019 {
1020 hist1.SetMaximum(max>0?max*1.05:max*0.95);
1021 hist1.SetMinimum(max>0?min*0.95:min*1.05);
1022 }
1023 }
1024
1025 TPaveText *t = (TPaveText*)gPad->FindObject("title");
1026 if (t)
1027 {
1028 t->SetName((TString)"MHTitle"); // rename object
1029 t->Clear(); // clear old lines
1030 t->AddText((TString)" "+title+" "); // add the new title
1031 t->SetBit(kCanDelete); // make sure object is deleted
1032
1033 //
1034 // FIXME: This is a stupid workaround to hide the redrawn
1035 // (see THistPainter::PaintTitle) title
1036 //
1037 gPad->Modified(); // indicates a change
1038 gPad->Update(); // recreates the original title
1039 t->Pop(); // bring our title on top
1040 }
1041
1042 //
1043 // Rename first statistics box
1044 //
1045 // Where to get the TPaveStats depends on the root version
1046 TPaveStats *s1 = dynamic_cast<TPaveStats*>(gPad->FindObject("stats"));
1047 if (!s1)
1048 s1 = dynamic_cast<TPaveStats*>(hist1.GetListOfFunctions()->FindObject("stats"));
1049 else
1050 s1->SetName((TString)"Stat"+hist1.GetTitle());
1051
1052 if (s1 && s1->GetX2NDC()>0.95)
1053 {
1054 const Double_t x1 = s1->GetX1NDC()-0.01;
1055 s1->SetX1NDC(x1-(s1->GetX2NDC()-s1->GetX1NDC()));
1056 s1->SetX2NDC(x1);
1057 }
1058
1059 //
1060 // Draw second histogram
1061 //
1062 hist2.Draw("sames");
1063 gPad->Update();
1064
1065 //
1066 // Draw Legend
1067 //
1068 // Where to get the TPaveStats depends on the root version
1069 TPaveStats *s2 = dynamic_cast<TPaveStats*>(gPad->FindObject("stats"));
1070 if (!s2)
1071 s2 = dynamic_cast<TPaveStats*>(hist2.GetListOfFunctions()->FindObject("stats"));
1072
1073 if (s2)
1074 {
1075 TLegend &l = *new TLegend(s2->GetX1NDC(),
1076 s2->GetY1NDC()-0.015-(s2->GetY2NDC()-s2->GetY1NDC())/2,
1077 s2->GetX2NDC(),
1078 s2->GetY1NDC()-0.01
1079 );
1080 l.AddEntry(&hist1, hist1.GetTitle());
1081 l.AddEntry(&hist2, hist2.GetTitle());
1082 l.SetTextSize(s2->GetTextSize());
1083 l.SetTextFont(s2->GetTextFont());
1084 l.SetBorderSize(s2->GetBorderSize());
1085 l.SetBit(kCanDelete);
1086 l.Draw();
1087 }
1088}
1089
1090// --------------------------------------------------------------------------
1091//
1092// If the opt string contains 'nonew' or gPad is not given NULL is returned.
1093// Otherwise the present gPad is returned.
1094//
1095TVirtualPad *MH::GetNewPad(TString &opt)
1096{
1097 opt.ToLower();
1098
1099 if (!opt.Contains("nonew"))
1100 return NULL;
1101
1102 opt.ReplaceAll("nonew", "");
1103
1104 return gPad;
1105}
1106
1107// --------------------------------------------------------------------------
1108//
1109// Encapsulate the TObject::Clone such, that a cloned TH1 (or derived)
1110// object is not added to the current directory, when cloned.
1111//
1112TObject *MH::Clone(const char *name) const
1113{
1114 const Bool_t store = TH1::AddDirectoryStatus();
1115
1116 TH1::AddDirectory(kFALSE);
1117 TObject *o = MParContainer::Clone(name);
1118 TH1::AddDirectory(store);
1119
1120 return o;
1121}
1122
1123// --------------------------------------------------------------------------
1124//
1125// If the opt string contains 'nonew' or gPad is not given a new canvas
1126// with size w/h is created. Otherwise the object is cloned and drawn
1127// to the present pad. The kCanDelete bit is set for the clone.
1128//
1129TObject *MH::DrawClone(Option_t *opt, Int_t w, Int_t h) const
1130{
1131 TString option(opt);
1132
1133 TVirtualPad *p = GetNewPad(option);
1134 if (!p)
1135 p = MakeDefCanvas(this, w, h);
1136 else
1137 if (!option.Contains("same", TString::kIgnoreCase))
1138 p->Clear();
1139
1140 gROOT->SetSelectedPad(NULL);
1141
1142 TObject *o = MParContainer::DrawClone(option);
1143 o->SetBit(kCanDelete);
1144 return o;
1145}
1146
1147// --------------------------------------------------------------------------
1148//
1149// Check whether a class inheriting from MH overwrites the Draw function
1150//
1151Bool_t MH::OverwritesDraw(TClass *cls) const
1152{
1153 if (!cls)
1154 cls = IsA();
1155
1156 //
1157 // Check whether we reached the base class MTask
1158 //
1159 if (TString(cls->GetName())=="MH")
1160 return kFALSE;
1161
1162 //
1163 // Check whether the class cls overwrites Draw
1164 //
1165 if (cls->GetMethodAny("Draw"))
1166 return kTRUE;
1167
1168 //
1169 // If the class itself doesn't overload it check all it's base classes
1170 //
1171 TBaseClass *base=NULL;
1172 TIter NextBase(cls->GetListOfBases());
1173 while ((base=(TBaseClass*)NextBase()))
1174 {
1175 if (OverwritesDraw(base->GetClassPointer()))
1176 return kTRUE;
1177 }
1178
1179 return kFALSE;
1180}
1181
1182// --------------------------------------------------------------------------
1183//
1184// Cuts the bins containing only zeros at the edges.
1185//
1186// A new number of bins can be defined with nbins != 0
1187// In the case of nbins == 0, no rebinning will take place
1188//
1189// Returns the new (real) number of bins
1190//
1191Int_t MH::StripZeros(TH1 *h, Int_t nbins)
1192{
1193 TAxis &axe = *h->GetXaxis();
1194
1195 const Int_t min1 = axe.GetFirst();
1196 const Int_t max1 = axe.GetLast();
1197 const Int_t range1 = max1-min1;
1198
1199 //
1200 // Check for useless zeros
1201 //
1202 if (range1 == 0)
1203 return 0;
1204
1205 Int_t min2 = 0;
1206 for (int i=min1; i<=max1; i++)
1207 if (h->GetBinContent(i) != 0)
1208 {
1209 min2 = i;
1210 break;
1211 }
1212
1213 //
1214 // If the histogram consists of zeros only
1215 //
1216 if (min2 == max1)
1217 return 0;
1218
1219 Int_t max2 = 0;
1220 for (int i=max1; i>=min2; i--)
1221 if (h->GetBinContent(i) != 0)
1222 {
1223 max2 = i;
1224 break;
1225 }
1226
1227 //
1228 // Appying TAxis->SetRange before ReBin does not work ...
1229 // But this workaround helps quite fine
1230 //
1231 Axis_t min = h->GetBinLowEdge(min2);
1232 Axis_t max = h->GetBinLowEdge(max2)+h->GetBinWidth(max2);
1233
1234 Int_t nbins2 = max2-min2;
1235 //
1236 // Check for rebinning
1237 //
1238 if (nbins > 0)
1239 {
1240 const Int_t ngroup = (Int_t)(nbins2*h->GetNbinsX()/nbins/(max1-min1));
1241 if (ngroup > 1)
1242 {
1243 h->Rebin(ngroup);
1244 nbins2 /= ngroup;
1245 }
1246 }
1247
1248 Int_t newbins = 0;
1249 FindGoodLimits(nbins2, newbins, min, max, kFALSE);
1250 axe.SetRangeUser(min,max);
1251 return axe.GetLast()-axe.GetFirst();
1252}
1253
1254// --------------------------------------------------------------------------
1255//
1256// In contradiction to TPad::FindObject this function searches recursively
1257// in a pad for an object. gPad is the default.
1258//
1259TObject *MH::FindObjectInPad(const char *name, TVirtualPad *pad)
1260{
1261 if (!pad)
1262 pad = gPad;
1263
1264 if (!pad)
1265 return NULL;
1266
1267 TObject *o;
1268
1269 TIter Next(pad->GetListOfPrimitives());
1270 while ((o=Next()))
1271 {
1272 if (!strcmp(o->GetName(), name))
1273 return o;
1274
1275 if (o->InheritsFrom("TPad"))
1276 if ((o = FindObjectInPad(name, (TVirtualPad*)o)))
1277 return o;
1278 }
1279 return NULL;
1280}
1281
1282// --------------------------------------------------------------------------
1283//
1284//
1285//
1286TH1I* MH::ProjectArray(const TArrayF &array, Int_t nbins,
1287 const char* name, const char* title)
1288{
1289 const Int_t size = array.GetSize();
1290
1291 TH1I *h1=0;
1292
1293 //check if histogram with identical name exist
1294 TObject *h1obj = gROOT->FindObject(name);
1295 if (h1obj && h1obj->InheritsFrom("TH1I"))
1296 {
1297 h1 = (TH1I*)h1obj;
1298 h1->Reset();
1299 }
1300
1301 Double_t min = size>0 ? array[0] : 0;
1302 Double_t max = size>0 ? array[0] : 1;
1303
1304 // first loop over array to find the min and max
1305 for (Int_t i=1; i<size;i++)
1306 {
1307 max = TMath::Max((Double_t)array[i], max);
1308 min = TMath::Min((Double_t)array[i], min);
1309 }
1310
1311 Int_t newbins = 0;
1312 FindGoodLimits(nbins, newbins, min, max, kFALSE);
1313
1314 if (!h1)
1315 {
1316 h1 = new TH1I(name, title, nbins, min, max);
1317 h1->SetXTitle("");
1318 h1->SetYTitle("Counts");
1319 h1->SetDirectory(NULL);
1320 }
1321
1322 // Second loop to fill the histogram
1323 for (Int_t i=0;i<size;i++)
1324 h1->Fill(array[i]);
1325
1326 return h1;
1327}
1328
1329// --------------------------------------------------------------------------
1330//
1331//
1332//
1333TH1I* MH::ProjectArray(const TArrayD &array, Int_t nbins, const char* name, const char* title)
1334{
1335 const Int_t size = array.GetSize();
1336
1337 Double_t min = size>0 ? array[0] : 0;
1338 Double_t max = size>0 ? array[0] : 1;
1339
1340 TH1I *h1=0;
1341
1342 //check if histogram with identical name exist
1343 TObject *h1obj = gROOT->FindObject(name);
1344 if (h1obj && h1obj->InheritsFrom("TH1I"))
1345 {
1346 h1 = (TH1I*)h1obj;
1347 h1->Reset();
1348 }
1349
1350 // first loop over array to find the min and max
1351 for (Int_t i=1; i<size;i++)
1352 {
1353 max = TMath::Max(array[i], max);
1354 min = TMath::Min(array[i], min);
1355 }
1356
1357 Int_t newbins = 0;
1358 FindGoodLimits(nbins, newbins, min, max, kFALSE);
1359
1360 if (!h1)
1361 {
1362 h1 = new TH1I(name, title, newbins, min, max);
1363 h1->SetXTitle("");
1364 h1->SetYTitle("Counts");
1365 h1->SetDirectory(NULL);
1366 }
1367
1368 // Second loop to fill the histogram
1369 for (Int_t i=0;i<size;i++)
1370 h1->Fill(array[i]);
1371
1372 return h1;
1373}
1374
1375// --------------------------------------------------------------------------
1376//
1377//
1378//
1379TH1I* MH::ProjectArray(const MArrayF &array, Int_t nbins,
1380 const char* name, const char* title)
1381{
1382 return ProjectArray(TArrayF(array.GetSize(),array.GetArray()), nbins, name, title);
1383}
1384
1385// --------------------------------------------------------------------------
1386//
1387//
1388//
1389TH1I* MH::ProjectArray(const MArrayD &array, Int_t nbins, const char* name, const char* title)
1390{
1391 return ProjectArray(TArrayD(array.GetSize(),array.GetArray()), nbins, name, title);
1392}
1393
1394// --------------------------------------------------------------------------
1395//
1396// This is a workaround for rrot-version <5.13/04 to get correct
1397// binomial errors even if weights have been used, do:
1398// h->Divide(h1, h2, 5, 2, "b");
1399// MH::SetBinomialErrors(*h, *h1, *h2, 5, 2);
1400//
1401// see http://root.cern.ch/phpBB2/viewtopic.php?p=14818
1402// see http://savannah.cern.ch/bugs/?20722
1403//
1404void MH::SetBinomialErrors(TH1 &hres, const TH1 &h1, const TH1 &h2, Double_t c1, Double_t c2)
1405{
1406 for (Int_t binx=0; binx<=hres.GetNbinsX()+1; binx++)
1407 {
1408 const Double_t b1 = h1.GetBinContent(binx);
1409 const Double_t b2 = h2.GetBinContent(binx);
1410 const Double_t e1 = h1.GetBinError(binx);
1411 const Double_t e2 = h2.GetBinError(binx);
1412
1413 //const Double_t w = c2*b2 ? (c1*b1)/(c2*b2) : 0;
1414 //const Double_t rc = ((1-2*w)*e1*e1+w*w*e2*e2)/(b2*b2);
1415
1416 if (b2==0)
1417 {
1418 hres.SetBinError(binx, 0);
1419 continue;
1420 }
1421
1422 const Double_t c = c2==0 ? 1 : c1/c2;
1423 const Double_t u = b2==0 ? 0 : b1/b2;
1424
1425 const Double_t rc = c*((1-2*u)*e1*e1+u*u*e2*e2)/(b2*b2);
1426
1427 hres.SetBinError(binx, TMath::Sqrt(TMath::Abs(rc)));
1428 }
1429}
1430
1431// --------------------------------------------------------------------------
1432//
1433// Return the first and last bin of the histogram which is not 0
1434//
1435void MH::GetRange(const TH1 &h, Int_t &lo, Int_t &hi)
1436{
1437 lo = 0;
1438 hi = 1;
1439
1440 for (int i=1; i<=h.GetNbinsX(); i++)
1441 {
1442 if (lo==0 && h.GetBinContent(i)>0)
1443 lo = i;
1444
1445 if (h.GetBinContent(i)>0)
1446 hi = i;
1447 }
1448}
1449
1450// --------------------------------------------------------------------------
1451//
1452// Return the lower edge of the first and the upper edge of the last bin
1453// of the histogram which is not 0
1454//
1455void MH::GetRangeUser(const TH1 &h, Axis_t &lo, Axis_t &hi)
1456{
1457 Int_t f, l;
1458 GetRange(h, f, l);
1459
1460 lo = h.GetBinLowEdge(f);
1461 hi = h.GetBinLowEdge(l+1);
1462}
1463
1464// --------------------------------------------------------------------------
1465//
1466// Return the first and last column (x-bin) of the histogram which is not 0.
1467// Therefor a proper projection is produced if the argument is a TH2.
1468//
1469// TH3 are not yet supported
1470//
1471void MH::GetRangeX(const TH1 &hist, Int_t &lo, Int_t &hi)
1472{
1473 if (hist.InheritsFrom(TH3::Class()))
1474 return;
1475
1476 if (hist.InheritsFrom(TH2::Class()))
1477 {
1478 TH1 *pro = static_cast<const TH2&>(hist).ProjectionX();
1479 GetRange(*pro, lo, hi);
1480 delete pro;
1481 return;
1482 }
1483
1484 GetRange(hist, lo, hi);
1485}
1486
1487// --------------------------------------------------------------------------
1488//
1489// Return the first and last row (y-bin) of the histogram which is not 0.
1490// Therefor a proper projection is produced if the argument is a TH2.
1491//
1492// TH3 are not yet supported
1493//
1494void MH::GetRangeY(const TH1 &hist, Int_t &lo, Int_t &hi)
1495{
1496 if (hist.InheritsFrom(TH3::Class()))
1497 return;
1498
1499 if (hist.InheritsFrom(TH2::Class()))
1500 {
1501 TH1 *pro = static_cast<const TH2&>(hist).ProjectionY();
1502 GetRange(*pro, lo, hi);
1503 delete pro;
1504 }
1505}
1506
1507// --------------------------------------------------------------------------
1508//
1509// Return the lower edge of the first and the upper edge of the last bin
1510// of the histogram h returned by GetRangeX
1511//
1512void MH::GetRangeUserX(const TH1 &h, Axis_t &lo, Axis_t &hi)
1513{
1514 Int_t f, l;
1515 GetRangeX(h, f, l);
1516
1517 lo = h.GetXaxis()->GetBinLowEdge(f);
1518 hi = h.GetXaxis()->GetBinLowEdge(l+1);
1519}
1520
1521// --------------------------------------------------------------------------
1522//
1523// Return the lower edge of the first and the upper edge of the last bin
1524// of the histogram h returned by GetRangeY
1525//
1526void MH::GetRangeUserY(const TH1 &h, Axis_t &lo, Axis_t &hi)
1527{
1528 Int_t f, l;
1529 GetRangeY(h, f, l);
1530
1531 lo = h.GetYaxis()->GetBinLowEdge(f);
1532 hi = h.GetYaxis()->GetBinLowEdge(l+1);
1533}
1534
1535// --------------------------------------------------------------------------
1536//
1537// See MTask::PrintSkipped
1538//
1539void MH::PrintSkipped(UInt_t n, const char *str)
1540{
1541 *fLog << " " << setw(7) << n;
1542 if (GetNumExecutions()>0)
1543 {
1544 *fLog << " (" << MString::Format("%5.1f", 100.*n/GetNumExecutions());
1545 *fLog << "%)" << str << endl;
1546 }
1547 *fLog << " Evts skipped: " << str << endl;
1548}
1549
1550#ifdef CreateGradientColorTable
1551#error CreateGradientColorTable already defined
1552#endif
1553
1554#if ROOT_VERSION_CODE < ROOT_VERSION(5,18,00)
1555#define CreateGradientColorTable gStyle->CreateGradientColorTable
1556#else
1557#define CreateGradientColorTable TColor::CreateGradientColorTable
1558#endif
1559
1560// --------------------------------------------------------------------------
1561//
1562// Calls gStyle->SetPalette. Allowed palettes are:
1563// pretty
1564// deepblue: darkblue -> lightblue
1565// lightblue: black -> blue -> white
1566// greyscale: black -> white
1567// glow1: black -> darkred -> orange -> yellow -> white
1568// glow2:
1569// glowsym: lightblue -> blue -> black -> darkred -> orange -> yellow -> white
1570// redish: darkred -> lightred
1571// bluish: darkblue -> lightblue
1572// small1:
1573//
1574// If the palette name contains 'inv' the order of the colors is inverted.
1575//
1576// The second argument determines the number of colors for the palette.
1577// The default is 50. 'pretty' always has 50 colors.
1578//
1579// (Remark: Drawing 3D object like TH2D with surf3 allows a maximum
1580// of 99 colors)
1581//
1582void MH::SetPalette(TString paletteName, Int_t ncol)
1583{
1584 Bool_t found=kFALSE;
1585
1586 paletteName.ToLower();
1587
1588 const Bool_t inverse = paletteName.Contains("inv");
1589
1590 if (paletteName.Contains("pretty"))
1591 {
1592 gStyle->SetPalette(1, 0);
1593 ncol=50;
1594 found=kTRUE;
1595 }
1596
1597 if (paletteName.Contains("deepblue"))
1598 {
1599 Double_t s[5] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
1600 Double_t r[5] = { 0.00, 0.09, 0.18, 0.09, 0.00 };
1601 Double_t g[5] = { 0.01, 0.02, 0.39, 0.68, 0.97 };
1602 Double_t b[5] = { 0.17, 0.39, 0.62, 0.79, 0.97 };
1603 CreateGradientColorTable(5, s, r, g, b, ncol);
1604 found=kTRUE;
1605 }
1606
1607 if (paletteName.Contains("lightblue"))
1608 {
1609 Double_t s[5] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
1610 Double_t r[5] = { 0.00, 0.09, 0.18, 0.09, 0.00 };
1611 Double_t g[5] = { 0.00, 0.02, 0.40, 0.70, 1.00 };
1612 Double_t b[5] = { 0.00, 0.27, 0.51, 0.81, 1.00 };
1613 CreateGradientColorTable(5, s, r, g, b, ncol);
1614 found=kTRUE;
1615 }
1616
1617 if (paletteName.Contains("greyscale"))
1618 {
1619 double s[2] = {0.00, 1.00};
1620 double r[2] = {0.00, 1.00};
1621 double g[2] = {0.00, 1.00};
1622 double b[2] = {0.00, 1.00};
1623 CreateGradientColorTable(2, s, r, g, b, ncol);
1624 found=kTRUE;
1625 }
1626
1627 if (paletteName.Contains("glow1"))
1628 {
1629 double s[5] = {0., 0.10, 0.45, 0.75, 1.00};
1630 double r[5] = {0., 0.35, 0.85, 1.00, 1.00};
1631 double g[5] = {0., 0.10, 0.20, 0.73, 1.00};
1632 double b[5] = {0., 0.03, 0.06, 0.00, 1.00};
1633 CreateGradientColorTable(5, s, r, g, b, ncol);
1634 found=kTRUE;
1635 }
1636
1637 if (paletteName.Contains("glow2"))
1638 {
1639 double s[4] = {0.00, 0.50, 0.75, 1.00};
1640 double r[4] = {0.24, 0.67, 1.00, 1.00};
1641 double g[4] = {0.03, 0.04, 0.80, 1.00};
1642 double b[4] = {0.03, 0.04, 0.00, 1.00};
1643 CreateGradientColorTable(4, s, r, g, b, ncol);
1644 found=kTRUE;
1645 }
1646
1647 if (paletteName.Contains("glowsym"))
1648 {
1649 double s[8] = {0.00, 0.17, 0.39, 0.50, 0.55, 0.72, 0.88, 1.00};
1650 double r[8] = {0.09, 0.18, 0.09, 0.00, 0.35, 0.85, 1.00, 1.00};
1651 double g[8] = {0.70, 0.40, 0.02, 0.00, 0.10, 0.20, 0.73, 1.00};
1652 double b[8] = {0.81, 0.51, 0.27, 0.00, 0.03, 0.06, 0.00, 1.00};
1653 CreateGradientColorTable(8, s, r, g, b, ncol);
1654 found=kTRUE;
1655 }/*
1656 if (paletteName.Contains("glows2"))
1657 {
1658 double s[10] = {0.00, 0.17, 0.35, 0.50, 0.65, 0.73, 0.77, 0.85, 0.92, 1.00};
1659 double r[10] = {0.09, 0.18, 0.09, 0.00, 0.00, 0.20, 0.55, 0.85, 1.00, 1.00};
1660 double g[10] = {0.81, 0.51, 0.27, 0.00, 0.00, 0.05, 0.10, 0.20, 0.73, 1.00};
1661 double b[10] = {0.70, 0.40, 0.02, 0.00, 0.27, 0.40, 0.35, 0.16, 0.03, 1.00};
1662 gStyle->CreateGradientColorTable(10, s, r, g, b, ncol);
1663 found=kTRUE;
1664 }*/
1665
1666 if (paletteName.Contains("redish"))
1667 {
1668 double s[3] = {0., 0.5, 1.};
1669 double r[3] = {0., 1.0, 1.};
1670 double g[3] = {0., 0.0, 1.};
1671 double b[3] = {0., 0.0, 1.};
1672 CreateGradientColorTable(3, s, r, g, b, ncol);
1673 found=kTRUE;
1674 }
1675
1676 if (paletteName.Contains("bluish"))
1677 {
1678 double s[3] = {0., 0.5, 1.};
1679 double r[3] = {0., 0.0, 1.};
1680 double g[3] = {0., 0.0, 1.};
1681 double b[3] = {0., 1.0, 1.};
1682 CreateGradientColorTable(3, s, r, g, b, ncol);
1683 found=kTRUE;
1684 }
1685
1686 if (paletteName.Contains("small1"))
1687 {
1688 double s[4] = {0.00, 0.50, 0.95, 1.};
1689 double r[4] = {0.04, 0.28, 0.98, 1.};
1690 double g[4] = {0.28, 0.93, 0.03, 1.};
1691 double b[4] = {0.79, 0.11, 0.03, 1.};
1692 CreateGradientColorTable(4, s, r, g, b, ncol);
1693 found=kTRUE;
1694 }
1695 if (paletteName.Contains("pepe"))
1696 {
1697 double s[5] = {0.0, 0.6, 0.7, 0.9, 1.0 };
1698 double r[5] = {0.1, 0.1, 1.0, 1.0, 1.0 };
1699 double g[5] = {0.1, 0.1, 0.0, 1.0, 1.0 };
1700 double b[5] = {0.2, 0.7, 0.0, 0.3, 0.9 };
1701 CreateGradientColorTable(5, s, r, g, b, ncol);
1702 found=kTRUE;
1703 }
1704 if (paletteName.Contains("temp"))
1705 {
1706 Double_t s[7] = { 0.00, 0.30, 0.45, 0.75, 0.88, 0.95, 1.00 };
1707 Double_t r[7] = { 0.00, 0.00, 0.30, 0.60, 1.00, 1.00, 1.00 };
1708 Double_t g[7] = { 0.00, 0.00, 0.00, 0.00, 0.20, 1.00, 1.00 };
1709 Double_t b[7] = { 0.10, 0.60, 0.30, 0.00, 0.20, 0.00, 1.00 };
1710 CreateGradientColorTable(7, s, r, g, b, ncol);
1711 found=kTRUE;
1712 }
1713
1714 if (inverse)
1715 {
1716 TArrayI c(ncol);
1717 for (int i=0; i<ncol; i++)
1718 c[ncol-i-1] = gStyle->GetColorPalette(i);
1719 gStyle->SetPalette(ncol, c.GetArray());
1720 }
1721
1722 if (!found)
1723 gLog << warn << "MH::SetPalette: Palette " << paletteName << " unknown... ignored." << endl;
1724}
1725
1726// --------------------------------------------------------------------------
1727//
1728// Unfortunately in TH1::GetObjectInfo the buffer is just 64 characters
1729// which is sometimes to small. This is just a copy of the code but the
1730// buffer has been increased to 128 which should fairly be enough.
1731//
1732// Necessary for root <= 5.22/00
1733//
1734char *MH::GetObjectInfoH(Int_t px, Int_t py, const TH1 &h)
1735{
1736 const TH1 *fH = &h;
1737 const TAxis *fXaxis = h.GetXaxis();
1738 const TAxis *fYaxis = h.GetYaxis();
1739
1740 // Redefines TObject::GetObjectInfo.
1741 // Displays the histogram info (bin number, contents, integral up to bin
1742 // corresponding to cursor position px,py
1743
1744 if (!gPad) return (char*)"";
1745
1746 static char info[128];
1747 Double_t x = gPad->PadtoX(gPad->AbsPixeltoX(px));
1748 Double_t y = gPad->PadtoY(gPad->AbsPixeltoY(py));
1749 Double_t x1 = gPad->PadtoX(gPad->AbsPixeltoX(px+1));
1750 const char *drawOption = fH->GetDrawOption();
1751 Double_t xmin, xmax, uxmin,uxmax;
1752 Double_t ymin, ymax, uymin,uymax;
1753 if (fH->GetDimension() == 2) {
1754 if (gPad->GetView() || strncmp(drawOption,"cont",4) == 0
1755 || strncmp(drawOption,"CONT",4) == 0) {
1756 uxmin=gPad->GetUxmin();
1757 uxmax=gPad->GetUxmax();
1758 xmin = fXaxis->GetBinLowEdge(fXaxis->GetFirst());
1759 xmax = fXaxis->GetBinUpEdge(fXaxis->GetLast());
1760 x = xmin +(xmax-xmin)*(x-uxmin)/(uxmax-uxmin);
1761 uymin=gPad->GetUymin();
1762 uymax=gPad->GetUymax();
1763 ymin = fYaxis->GetBinLowEdge(fYaxis->GetFirst());
1764 ymax = fYaxis->GetBinUpEdge(fYaxis->GetLast());
1765 y = ymin +(ymax-ymin)*(y-uymin)/(uymax-uymin);
1766 }
1767 }
1768 Int_t binx,biny,binmin,binx1;
1769 if (gPad->IsVertical()) {
1770 binx = fXaxis->FindFixBin(x);
1771 binmin = fXaxis->GetFirst();
1772 binx1 = fXaxis->FindFixBin(x1);
1773 // special case if more than 1 bin in x per pixel
1774 if (binx1-binx>1 && fH->GetDimension() == 1) {
1775 Double_t binval=fH->GetBinContent(binx);
1776 Int_t binnear=binx;
1777 for (Int_t ibin=binx+1; ibin<binx1; ibin++) {
1778 Double_t binvaltmp = fH->GetBinContent(ibin);
1779 if (TMath::Abs(y-binvaltmp) < TMath::Abs(y-binval)) {
1780 binval=binvaltmp;
1781 binnear=ibin;
1782 }
1783 }
1784 binx = binnear;
1785 }
1786 } else {
1787 x1 = gPad->PadtoY(gPad->AbsPixeltoY(py+1));
1788 binx = fXaxis->FindFixBin(y);
1789 binmin = fXaxis->GetFirst();
1790 binx1 = fXaxis->FindFixBin(x1);
1791 // special case if more than 1 bin in x per pixel
1792 if (binx1-binx>1 && fH->GetDimension() == 1) {
1793 Double_t binval=fH->GetBinContent(binx);
1794 Int_t binnear=binx;
1795 for (Int_t ibin=binx+1; ibin<binx1; ibin++) {
1796 Double_t binvaltmp = fH->GetBinContent(ibin);
1797 if (TMath::Abs(x-binvaltmp) < TMath::Abs(x-binval)) {
1798 binval=binvaltmp;
1799 binnear=ibin;
1800 }
1801 }
1802 binx = binnear;
1803 }
1804 }
1805 if (fH->GetDimension() == 1) {
1806 Double_t integ = 0;
1807 for (Int_t bin=binmin;bin<=binx;bin++) {integ += fH->GetBinContent(bin);}
1808 sprintf(info,"(x=%g, y=%g, binx=%d, binc=%g, Sum=%g)",x,y,binx,fH->GetBinContent(binx),integ);
1809 } else {
1810 biny = fYaxis->FindFixBin(y);
1811 sprintf(info,"(x=%g, y=%g, binx=%d, biny=%d, binc=%g)",x,y,binx,biny,fH->GetCellContent(binx,biny));
1812 }
1813 return info;
1814}
1815
1816// --------------------------------------------------------------------------
1817//
1818// Unfortunately in TProfile::GetObjectInfo the buffer is just 64 characters
1819// which is sometimes to small. This is just a copy of the code but the
1820// buffer has been increased to 128 which should fairly be enough.
1821//
1822// Necessary for root <= 5.22/00
1823//
1824char *MH::GetObjectInfoP(Int_t px, Int_t py, const TProfile &p)
1825{
1826 if (!gPad) return (char*)"";
1827 static char info[128];
1828 Double_t x = gPad->PadtoX(gPad->AbsPixeltoX(px));
1829 Double_t y = gPad->PadtoY(gPad->AbsPixeltoY(py));
1830 Int_t binx = p.GetXaxis()->FindFixBin(x);
1831 sprintf(info,"(x=%g, y=%g, binx=%d, binc=%g, bine=%g, binn=%d)", x, y, binx, p.GetBinContent(binx), p.GetBinError(binx), (Int_t)p.GetBinEntries(binx));
1832 return info;
1833}
1834
1835// --------------------------------------------------------------------------
1836//
1837// Unfortunately TH1::GetObjectInfo and TProfile::GetObjectInfo can
1838// result in buffer ovwerflows therefor we have to re-implement these
1839// function by our own.
1840//
1841// Necessary for root <= 5.22/00
1842//
1843char *MH::GetObjectInfo(Int_t px, Int_t py, const TObject &o)
1844{
1845 if (!o.InheritsFrom(TH1::Class()))
1846 return o.GetObjectInfo(px, py);
1847
1848 if (o.InheritsFrom(TProfile::Class()))
1849 return GetObjectInfoP(px, py, static_cast<const TProfile&>(o));
1850
1851 if (o.InheritsFrom("MHCamera"))
1852 return o.GetObjectInfo(px, py);
1853
1854 if (o.InheritsFrom(TH1::Class()))
1855 return GetObjectInfoH(px, py, static_cast<const TH1&>(o));
1856
1857 return const_cast<char*>("MH::GetObjectInfo: unknown class.");
1858}
1859
1860// --------------------------------------------------------------------------
1861//
1862// Set the pad-range such that at the smallest width it is still
1863// two times max and that the displayed "pixels" (i.e. its coordinate
1864// system) have the width to height ratio of aspect.
1865//
1866void MH::SetPadRange(Float_t max, Float_t aspect)
1867{
1868 if (!gPad)
1869 return;
1870
1871 const Float_t w = gPad->GetWw();
1872 const Float_t h = gPad->GetWh();
1873
1874 if (w>aspect*h)
1875 {
1876 const Double_t dx = ((w/h-aspect)/2+1)*max;
1877 gPad->Range(-dx, -max, dx, max);
1878 }
1879 else
1880 {
1881 const Double_t dy = ((h/w-1./aspect)/2+1)*max;
1882 gPad->Range(-max, -dy, max, dy);
1883 }
1884}
1885
1886// --------------------------------------------------------------------------
1887//
1888// Set the range of a pad in a way that the coordinates fit into the pad
1889// without abberation.
1890//
1891void MH::SetPadRange(Float_t x0, Float_t y0, Float_t x1, Float_t y1)
1892{
1893 if (!gPad)
1894 return;
1895
1896 const Float_t w = x1-x0; // Width in user coordinates
1897 const Float_t h = y1-y0; // Hieght in user coordinates
1898
1899 const Float_t ww = gPad->GetWw()*gPad->GetAbsWNDC(); // Width of pad in pixels
1900 const Float_t hh = gPad->GetWh()*gPad->GetAbsHNDC(); // Height of pad in pixels
1901
1902 if (ww/hh > w/h)
1903 {
1904 const Double_t dx = (ww/hh-w/h)/2*h;
1905
1906 gPad->Range(x0-dx, y0, x1+dx, y1);
1907 }
1908 else
1909 {
1910 const Double_t dy = (hh/ww-h/w)/2*w;
1911
1912 gPad->Range(x0, y0-dy, x1, y1+dy);
1913 }
1914}
Note: See TracBrowser for help on using the repository browser.