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

Last change on this file since 5250 was 5143, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 37.9 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 07/2001 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2002
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26// //
27// MH //
28// //
29// This is a base tasks for mars histograms. It defines a common interface //
30// for filling the histograms with events (MH::Fill) which is used by a //
31// common 'filler' And a SetupFill member function which may be used //
32// by MFillH. The idea is: //
33// 1) If your Histogram can become filled by one single container //
34// (like MHHillas) you overload MH::Fill and it gets called with //
35// a pointer to the container with which it should be filled. //
36// //
37// 2) You histogram needs several containers to get filled. Than you //
38// have to overload MH::SetupFill and get the necessary objects from //
39// the parameter list. Use this objects in Fill to fill your //
40// histogram. //
41// //
42// If you want to create your own histogram class the new class must be //
43// derived from MH (instead of the base MParContainer) and you must //
44// the fill function of MH. This is the function which is called to fill //
45// the histogram(s) by the data of a corresponding parameter container. //
46// //
47// Remark: the static member function (eg MakeDefCanvas) can be called //
48// from everywhere using: MH::MakeDefCanvas(...) //
49// //
50//////////////////////////////////////////////////////////////////////////////
51
52#include "MH.h"
53
54#include <TH1.h>
55#include <TH2.h>
56#include <TH3.h>
57#include <TStyle.h> // TStyle::GetScreenFactor
58#include <TGaxis.h>
59#include <TCanvas.h>
60#include <TLegend.h>
61#include <TPaveStats.h>
62#include <TBaseClass.h>
63#if ROOT_VERSION_CODE > ROOT_VERSION(3,04,01)
64#include <THLimitsFinder.h>
65#endif
66
67#include "MLog.h"
68#include "MLogManip.h"
69
70#include "MParList.h"
71#include "MParContainer.h"
72
73#include "MBinning.h"
74
75#include "MArrayD.h"
76#include "MArrayF.h"
77
78ClassImp(MH);
79
80using namespace std;
81
82// --------------------------------------------------------------------------
83//
84// Default Constructor. It sets name and title only. Typically you won't
85// need to change this.
86//
87MH::MH(const char *name, const char *title)
88{
89 //
90 // set the name and title of this object
91 //
92 fName = name ? name : "MH";
93 fTitle = title ? title : "Base class for Mars histograms";
94}
95
96// --------------------------------------------------------------------------
97//
98// If you want to use the automatic filling of your derived class you
99// must overload this function. If it is not overloaded you cannot use
100// FillH with this class. The argument is a pointer to a container
101// in your paremeter list which is specified in the MFillH constructor.
102// If you are not going to use it you should at least add
103// Bool_t MH::Fill(const MParContainer *) { return kTRUE; }
104// to your class definition.
105//
106Bool_t MH::Fill(const MParContainer *par, const Stat_t w)
107{
108 *fLog << warn << GetDescriptor() << ": Fill not overloaded! Can't be used!" << endl;
109 return kFALSE;
110}
111
112// --------------------------------------------------------------------------
113//
114// This virtual function is ment as a generalized interface to retrieve
115// a pointer to a root histogram from the MH-derived class.
116//
117TH1 *MH::GetHistByName(const TString name)
118{
119 *fLog << warn << GetDescriptor() << ": GetHistByName not overloaded! Can't be used!" << endl;
120 return NULL;
121}
122
123// --------------------------------------------------------------------------
124//
125// This is a function which should replace the creation of default
126// canvases like root does. Because this is inconvinient in some aspects.
127// need to change this.
128// You can specify a name for the default canvas and a title. Also
129// width and height can be given.
130// MakeDefCanvas looks for a canvas with the given name. If now name is
131// given the DefCanvasName of root is used. If no such canvas is existing
132// it is created and returned. If such a canvas already exists a new canvas
133// with a name plus anumber is created (the number is calculated by the
134// number of all existing canvases plus one)
135//
136// Normally the canvas size is scaled with gStyle->GetScreenFactor() so
137// that on all screens it looks like the same part of the screen.
138// To suppress this scaling use usescreenfactor=kFALSE. In this case
139// you specify directly the size of the embedded pad.
140//
141TCanvas *MH::MakeDefCanvas(TString name, const char *title,
142 UInt_t w, UInt_t h, Bool_t usescreenfactor)
143{
144 const TList *list = (TList*)gROOT->GetListOfCanvases();
145
146 if (name.IsNull())
147 name = gROOT->GetDefCanvasName();
148
149 if (list->FindObject(name))
150 name += Form(" <%d>", list->GetSize()+1);
151
152 if (!usescreenfactor)
153 {
154 const Float_t cx = gStyle->GetScreenFactor();
155 w += 4;
156 h += 28;
157 w = (int)(w/cx+1);
158 h = (int)(h/cx+1);
159 }
160
161 return new TCanvas(name, title, w, h);
162}
163
164// --------------------------------------------------------------------------
165//
166// This function works like MakeDefCanvas(name, title, w, h) but name
167// and title are retrieved from the given TObject.
168//
169// Normally the canvas size is scaled with gStyle->GetScreenFactor() so
170// that on all screens it looks like the same part of the screen.
171// To suppress this scaling use usescreenfactor=kFALSE. In this case
172// you specify directly the size of the embedded pad.
173//
174TCanvas *MH::MakeDefCanvas(const TObject *obj,
175 UInt_t w, UInt_t h, Bool_t usescreenfactor)
176{
177 if (!usescreenfactor)
178 {
179 const Float_t cx = gStyle->GetScreenFactor();
180 w += 4;
181 h += 28;
182 h = (int)(h/cx+1);
183 w = (int)(w/cx+1);
184 }
185
186 return MakeDefCanvas(obj->GetName(), obj->GetTitle(), w, h);
187}
188
189// --------------------------------------------------------------------------
190//
191// Applies a given binning to a 1D-histogram
192//
193void MH::SetBinning(TH1 *h, const MBinning *binsx)
194{
195 //
196 // Another strange behaviour: TAxis::Set deletes the axis title!
197 //
198 TAxis &x = *h->GetXaxis();
199
200#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
201 TString xtitle = x.GetTitle();
202#endif
203
204 // All this is reset by TAxis::Set
205 const TAttAxis att(x);
206 const Bool_t tm(x.GetTimeDisplay());
207 const TString tf(x.GetTimeFormat());
208
209 //
210 // This is a necessary workaround if one wants to set
211 // non-equidistant bins after the initialization
212 // TH1D::fNcells must be set correctly.
213 //
214 h->SetBins(binsx->GetNumBins(), 0, 1);
215
216 //
217 // Set the binning of the current histogram to the binning
218 // in one of the two given histograms
219 //
220 x.Set(binsx->GetNumBins(), binsx->GetEdges());
221
222 // All this is reset by TAxis::Set
223 att.Copy(x);
224 x.SetTimeDisplay(tm);
225 x.SetTimeFormat(tf);
226
227#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
228 x.SetTitle(xtitle);
229#endif
230}
231
232// --------------------------------------------------------------------------
233//
234// Applies given binnings to the two axis of a 2D-histogram
235//
236void MH::SetBinning(TH2 *h, const MBinning *binsx, const MBinning *binsy)
237{
238 TAxis &x = *h->GetXaxis();
239 TAxis &y = *h->GetYaxis();
240
241 //
242 // Another strange behaviour: TAxis::Set deletes the axis title!
243 //
244#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
245 TString xtitle = x.GetTitle();
246 TString ytitle = y.GetTitle();
247#endif
248
249 // All this is reset by TAxis::Set
250 const TAttAxis attx(x);
251 const TAttAxis atty(y);
252 const Bool_t tmx(x.GetTimeDisplay());
253 const Bool_t tmy(y.GetTimeDisplay());
254 const TString tfx(x.GetTimeFormat());
255 const TString tfy(y.GetTimeFormat());
256
257 //
258 // This is a necessary workaround if one wants to set
259 // non-equidistant bins after the initialization
260 // TH1D::fNcells must be set correctly.
261 //
262 h->SetBins(binsx->GetNumBins(), 0, 1,
263 binsy->GetNumBins(), 0, 1);
264
265 //
266 // Set the binning of the current histogram to the binning
267 // in one of the two given histograms
268 //
269 x.Set(binsx->GetNumBins(), binsx->GetEdges());
270 y.Set(binsy->GetNumBins(), binsy->GetEdges());
271
272 // All this is reset by TAxis::Set
273 attx.Copy(x);
274 atty.Copy(y);
275 x.SetTimeDisplay(tmx);
276 y.SetTimeDisplay(tmy);
277 x.SetTimeFormat(tfx);
278 y.SetTimeFormat(tfy);
279
280#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
281 x.SetTitle(xtitle);
282 y.SetTitle(ytitle);
283#endif
284}
285
286// --------------------------------------------------------------------------
287//
288// Applies given binnings to the three axis of a 3D-histogram
289//
290void MH::SetBinning(TH3 *h, const MBinning *binsx, const MBinning *binsy, const MBinning *binsz)
291{
292 //
293 // Another strange behaviour: TAxis::Set deletes the axis title!
294 //
295 TAxis &x = *h->GetXaxis();
296 TAxis &y = *h->GetYaxis();
297 TAxis &z = *h->GetZaxis();
298
299#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
300 TString xtitle = x.GetTitle();
301 TString ytitle = y.GetTitle();
302 TString ztitle = z.GetTitle();
303#endif
304
305 // All this is reset by TAxis::Set
306 const TAttAxis attx(x);
307 const TAttAxis atty(y);
308 const TAttAxis attz(z);
309 const Bool_t tmx(x.GetTimeDisplay());
310 const Bool_t tmy(y.GetTimeDisplay());
311 const Bool_t tmz(z.GetTimeDisplay());
312 const TString tfx(x.GetTimeFormat());
313 const TString tfy(y.GetTimeFormat());
314 const TString tfz(z.GetTimeFormat());
315
316 //
317 // This is a necessary workaround if one wants to set
318 // non-equidistant bins after the initialization
319 // TH1D::fNcells must be set correctly.
320 //
321 h->SetBins(binsx->GetNumBins(), 0, 1,
322 binsy->GetNumBins(), 0, 1,
323 binsz->GetNumBins(), 0, 1);
324
325 //
326 // Set the binning of the current histogram to the binning
327 // in one of the two given histograms
328 //
329 x.Set(binsx->GetNumBins(), binsx->GetEdges());
330 y.Set(binsy->GetNumBins(), binsy->GetEdges());
331 z.Set(binsz->GetNumBins(), binsz->GetEdges());
332
333 // All this is reset by TAxis::Set
334 attx.Copy(x);
335 atty.Copy(y);
336 attz.Copy(z);
337 x.SetTimeDisplay(tmx);
338 y.SetTimeDisplay(tmy);
339 z.SetTimeDisplay(tmz);
340 x.SetTimeFormat(tfx);
341 y.SetTimeFormat(tfy);
342 z.SetTimeFormat(tfz);
343
344#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
345 x.SetTitle(xtitle);
346 y.SetTitle(ytitle);
347 z.SetTitle(ztitle);
348#endif
349}
350
351// --------------------------------------------------------------------------
352//
353// Applies given binning (the n+1 edges) to the axis of a 1D-histogram
354//
355void MH::SetBinning(TH1 *h, const TArrayD &binsx)
356{
357 MBinning bx;
358 bx.SetEdges(binsx);
359 SetBinning(h, &bx);
360}
361
362// --------------------------------------------------------------------------
363//
364// Applies given binning (the n+1 edges) to the two axis of a
365// 2D-histogram
366//
367void MH::SetBinning(TH2 *h, const TArrayD &binsx, const TArrayD &binsy)
368{
369 MBinning bx;
370 MBinning by;
371 bx.SetEdges(binsx);
372 by.SetEdges(binsy);
373 SetBinning(h, &bx, &by);
374}
375
376// --------------------------------------------------------------------------
377//
378// Applies given binning (the n+1 edges) to the three axis of a
379// 3D-histogram
380//
381void MH::SetBinning(TH3 *h, const TArrayD &binsx, const TArrayD &binsy, const TArrayD &binsz)
382{
383 MBinning bx;
384 MBinning by;
385 MBinning bz;
386 bx.SetEdges(binsx);
387 by.SetEdges(binsy);
388 bz.SetEdges(binsz);
389 SetBinning(h, &bx, &by, &bz);
390}
391
392// --------------------------------------------------------------------------
393//
394// Applies the binning of a TAxis (eg from a root histogram) to the axis
395// of a 1D-histogram
396//
397void MH::SetBinning(TH1 *h, const TAxis *binsx)
398{
399 const Int_t nx = binsx->GetNbins();
400
401 TArrayD bx(nx+1);
402 for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
403 bx[nx] = binsx->GetXmax();
404
405 SetBinning(h, bx);
406}
407
408// --------------------------------------------------------------------------
409//
410// Applies the binnings of the TAxis' (eg from a root histogram) to the
411// two axis' of a 2D-histogram
412//
413void MH::SetBinning(TH2 *h, const TAxis *binsx, const TAxis *binsy)
414{
415 const Int_t nx = binsx->GetNbins();
416 const Int_t ny = binsy->GetNbins();
417
418 TArrayD bx(nx+1);
419 TArrayD by(ny+1);
420 for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
421 for (int i=0; i<ny; i++) by[i] = binsy->GetBinLowEdge(i+1);
422 bx[nx] = binsx->GetXmax();
423 by[ny] = binsy->GetXmax();
424
425 SetBinning(h, bx, by);
426}
427
428// --------------------------------------------------------------------------
429//
430// Applies the binnings of the TAxis' (eg from a root histogram) to the
431// three axis' of a 3D-histogram
432//
433void MH::SetBinning(TH3 *h, const TAxis *binsx, const TAxis *binsy, const TAxis *binsz)
434{
435 const Int_t nx = binsx->GetNbins();
436 const Int_t ny = binsy->GetNbins();
437 const Int_t nz = binsz->GetNbins();
438
439 TArrayD bx(nx+1);
440 TArrayD by(ny+1);
441 TArrayD bz(nz+1);
442 for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
443 for (int i=0; i<ny; i++) by[i] = binsy->GetBinLowEdge(i+1);
444 for (int i=0; i<nz; i++) bz[i] = binsz->GetBinLowEdge(i+1);
445 bx[nx] = binsx->GetXmax();
446 by[ny] = binsy->GetXmax();
447 bz[nz] = binsz->GetXmax();
448
449 SetBinning(h, bx, by, bz);
450}
451
452// --------------------------------------------------------------------------
453//
454// Applies the binnings of one root-histogram x to another one h
455// Both histograms must be of the same type: TH1, TH2 or TH3
456//
457void MH::SetBinning(TH1 *h, const TH1 *x)
458{
459 if (h->InheritsFrom(TH3::Class()) && x->InheritsFrom(TH3::Class()))
460 {
461 SetBinning((TH3*)h, ((TH1*)x)->GetXaxis(), ((TH1*)x)->GetYaxis(), ((TH1*)x)->GetZaxis());
462 return;
463 }
464 if (h->InheritsFrom(TH3::Class()) || x->InheritsFrom(TH3::Class()))
465 return;
466 if (h->InheritsFrom(TH2::Class()) && x->InheritsFrom(TH2::Class()))
467 {
468 SetBinning((TH2*)h, ((TH1*)x)->GetXaxis(), ((TH1*)x)->GetYaxis());
469 return;
470 }
471 if (h->InheritsFrom(TH2::Class()) || x->InheritsFrom(TH2::Class()))
472 return;
473 if (h->InheritsFrom(TH1::Class()) && x->InheritsFrom(TH1::Class()))
474 {
475 SetBinning(h, ((TH1*)x)->GetXaxis());
476 return;
477 }
478}
479
480void MH::RemoveFirstBin(TH1 &h)
481{
482 if (h.InheritsFrom(TH2::Class()) || h.InheritsFrom(TH3::Class()))
483 return;
484
485 const Bool_t haserr = h.GetSumw2N()>0;
486
487 const Int_t n0 = h.GetNbinsX();
488 if (n0<2)
489 return;
490
491 TArrayD val(n0-1);
492 TArrayD err(haserr ? n0-1 : 0);
493 for (int i=1; i<n0; i++)
494 {
495 val[i-1] = h.GetBinContent(i+1);
496 if (haserr)
497 err[i-1] = h.GetBinError(i+1);
498 }
499
500 MBinning bins;
501 bins.SetEdges(h, 'x');
502 bins.RemoveFirstEdge();
503 bins.Apply(h);
504
505 h.Reset();
506
507 for (int i=1; i<n0; i++)
508 {
509 h.SetBinContent(i, val[i-1]);
510 if (haserr)
511 h.SetBinError(i, err[i-1]);
512 }
513}
514
515// --------------------------------------------------------------------------
516//
517// Multiplies all entries in a TArrayD by a float f
518//
519void MH::ScaleArray(TArrayD &bins, Double_t f)
520{
521 for (int i=0; i<bins.GetSize(); i++)
522 bins[i] *= f;
523}
524
525// --------------------------------------------------------------------------
526//
527// Scales the binning of a TAxis by a float f
528//
529TArrayD MH::ScaleAxis(TAxis &axe, Double_t f)
530{
531 TArrayD arr(axe.GetNbins()+1);
532
533 for (int i=1; i<=axe.GetNbins()+1; i++)
534 arr[i-1] = axe.GetBinLowEdge(i);
535
536 ScaleArray(arr, f);
537
538 return arr;
539}
540
541// --------------------------------------------------------------------------
542//
543// Scales the binning of one, two or three axis of a histogram by a float f
544//
545void MH::ScaleAxis(TH1 *h, Double_t fx, Double_t fy, Double_t fz)
546{
547 if (h->InheritsFrom(TH3::Class()))
548 {
549 SetBinning((TH3*)h,
550 ScaleAxis(*h->GetXaxis(), fx),
551 ScaleAxis(*h->GetYaxis(), fy),
552 ScaleAxis(*h->GetZaxis(), fz));
553 return;
554 }
555
556 if (h->InheritsFrom(TH2::Class()))
557 {
558 SetBinning((TH2*)h,
559 ScaleAxis(*h->GetXaxis(), fx),
560 ScaleAxis(*h->GetYaxis(), fy));
561 return;
562 }
563
564 if (h->InheritsFrom(TH1::Class()))
565 SetBinning(h, ScaleAxis(*h->GetXaxis(), fx));
566}
567
568// --------------------------------------------------------------------------
569//
570// Tries to find a MBinning container with the name "Binning"+name
571// in the given parameter list. If it was found it is applied to the
572// given histogram. This is only valid for 1D-histograms
573//
574Bool_t MH::ApplyBinning(const MParList &plist, TString name, TH1 *h)
575{
576 if (h->InheritsFrom(TH2::Class()) || h->InheritsFrom(TH3::Class()))
577 {
578 gLog << warn << "MH::ApplyBinning: '" << h->GetName() << "' is not a basic TH1 object... no binning applied." << endl;
579 return kFALSE;
580 }
581
582 const MBinning *bins = (MBinning*)plist.FindObject("Binning"+name);
583 if (!bins)
584 {
585 gLog << inf << "Object 'Binning" << name << "' [MBinning] not found... no binning applied." << endl;
586 return kFALSE;
587 }
588
589 SetBinning(h, bins);
590 return kTRUE;
591}
592
593void MH::FindGoodLimits(Int_t nbins, Int_t &newbins, Double_t &xmin, Double_t &xmax, Bool_t isInteger)
594{
595#if ROOT_VERSION_CODE > ROOT_VERSION(3,04,01)
596 THLimitsFinder::OptimizeLimits(nbins, newbins, xmin, xmax, isInteger);
597#else
598//*-*-*-*-*-*-*-*-*Find reasonable bin values*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
599//*-* ==========================
600
601 Double_t dx = 0.1*(xmax-xmin);
602 Double_t umin = xmin - dx;
603 Double_t umax = xmax + dx;
604
605 if (umin < 0 && xmin >= 0)
606 umin = 0;
607
608 if (umax > 0 && xmax <= 0)
609 umax = 0;
610
611 Double_t binlow =0;
612 Double_t binhigh =0;
613 Double_t binwidth=0;
614
615 TGaxis::Optimize(umin, umax, nbins, binlow, binhigh, nbins, binwidth, "");
616
617 if (binwidth <= 0 || binwidth > 1.e+39)
618 {
619 xmin = -1;
620 xmax = 1;
621 }
622 else
623 {
624 xmin = binlow;
625 xmax = binhigh;
626 }
627
628 if (isInteger)
629 {
630 Int_t ixmin = (Int_t)xmin;
631 Int_t ixmax = (Int_t)xmax;
632 Double_t dxmin = (Double_t)ixmin;
633 Double_t dxmax = (Double_t)ixmax;
634
635 xmin = xmin<0 && xmin!=dxmin ? dxmin - 1 : dxmin;
636 xmax = xmax>0 && xmax!=dxmax ? dxmax + 1 : dxmax;
637
638 if (xmin>=xmax)
639 xmax = xmin+1;
640
641 Int_t bw = 1 + (Int_t)((xmax-xmin)/nbins);
642
643 nbins = (Int_t)((xmax-xmin)/bw);
644
645 if (xmin+nbins*bw < xmax)
646 {
647 nbins++;
648 xmax = xmin +nbins*bw;
649 }
650 }
651
652 newbins = nbins;
653#endif
654}
655
656// --------------------------------------------------------------------------
657//
658// Returns the lowest entry in a histogram which is greater than gt (eg >0)
659//
660Double_t MH::GetMinimumGT(const TH1 &h, Double_t gt)
661{
662 Double_t min = FLT_MAX;
663
664 const TAxis &axex = *((TH1&)h).GetXaxis();
665 const TAxis &axey = *((TH1&)h).GetYaxis();
666 const TAxis &axez = *((TH1&)h).GetZaxis();
667
668 for (int iz=1; iz<=axez.GetNbins(); iz++)
669 for (int iy=1; iy<=axey.GetNbins(); iy++)
670 for (int ix=1; ix<=axex.GetNbins(); ix++)
671 {
672 const Double_t v = h.GetBinContent(h.GetBin(ix, iy, iz));
673 if (gt<v && v<min)
674 min = v;
675 }
676 return min;
677}
678
679// --------------------------------------------------------------------------
680//
681// Returns the bin center in a logarithmic scale. If the given bin
682// number is <1 it is set to 1. If it is =GetNbins() it is set to
683// GetNbins()
684//
685Double_t MH::GetBinCenterLog(const TAxis &axe, Int_t nbin)
686{
687 if (nbin>axe.GetNbins())
688 nbin = axe.GetNbins();
689
690 if (nbin<1)
691 nbin = 1;
692
693 const Double_t lo = axe.GetBinLowEdge(nbin);
694 const Double_t hi = axe.GetBinUpEdge(nbin);
695
696 const Double_t val = log10(lo) + log10(hi);
697
698 return pow(10, val/2);
699}
700
701// --------------------------------------------------------------------------
702//
703// Draws a copy of the two given histograms. Uses title as the pad title.
704// Also layout the two statistic boxes and a legend.
705//
706void MH::DrawSameCopy(const TH1 &hist1, const TH1 &hist2, const TString title)
707{
708 //
709 // Draw first histogram
710 //
711 TH1 *h1 = ((TH1&)hist1).DrawCopy();
712 gPad->SetBorderMode(0);
713 gPad->Update();
714
715 // FIXME: Also align max/min with set Maximum/Minimum
716 const Double_t maxbin1 = hist1.GetBinContent(hist1.GetMaximumBin());
717 const Double_t maxbin2 = hist2.GetBinContent(hist2.GetMaximumBin());
718 const Double_t minbin1 = hist1.GetBinContent(hist1.GetMinimumBin());
719 const Double_t minbin2 = hist2.GetBinContent(hist2.GetMinimumBin());
720
721 const Double_t max = TMath::Max(maxbin1, maxbin2);
722 const Double_t min = TMath::Min(minbin1, minbin2);
723
724 h1->SetMaximum(max>0?max*1.05:max*0.95);
725 h1->SetMinimum(max>0?min*0.95:min*1.05);
726
727 TPaveText *t = (TPaveText*)gPad->FindObject("title");
728 if (t)
729 {
730 t->SetName((TString)"MHTitle"); // rename object
731 t->Clear(); // clear old lines
732 t->AddText((TString)" "+title+" "); // add the new title
733 t->SetBit(kCanDelete); // make sure object is deleted
734
735 //
736 // FIXME: This is a stupid workaround to hide the redrawn
737 // (see THistPainter::PaintTitle) title
738 //
739 gPad->Modified(); // indicates a change
740 gPad->Update(); // recreates the original title
741 t->Pop(); // bring our title on top
742 }
743
744 //
745 // Rename first statistics box
746 //
747 TPaveStats *s1 = (TPaveStats*)gPad->FindObject("stats");
748 if (!s1)
749 s1 = (TPaveStats*)hist1.GetListOfFunctions()->FindObject("stats");
750 else
751 s1->SetName((TString)"Stat"+hist1.GetTitle());
752
753 if (s1 && s1->GetX2NDC()>0.95)
754 {
755 const Double_t x1 = s1->GetX1NDC()-0.01;
756 s1->SetX1NDC(x1-(s1->GetX2NDC()-s1->GetX1NDC()));
757 s1->SetX2NDC(x1);
758 }
759
760 //
761 // Draw second histogram
762 //
763 TH1 *h2 = ((TH1&)hist2).DrawCopy("sames");
764 gPad->Update();
765
766 //
767 // Draw Legend
768 //
769 TPaveStats *s2 = (TPaveStats*)gPad->FindObject("stats");
770 if (!s2)
771 s2 = (TPaveStats*)hist2.GetListOfFunctions()->FindObject("stats");
772
773 if (s2)
774 {
775 TLegend &l = *new TLegend(s2->GetX1NDC(),
776 s2->GetY1NDC()-0.015-(s2->GetY2NDC()-s2->GetY1NDC())/2,
777 s2->GetX2NDC(),
778 s2->GetY1NDC()-0.01
779 );
780 l.AddEntry(h1, h1->GetTitle());
781 l.AddEntry(h2, h2->GetTitle());
782 l.SetTextSize(s2->GetTextSize());
783 l.SetTextFont(s2->GetTextFont());
784 l.SetBorderSize(s2->GetBorderSize());
785 l.SetBit(kCanDelete);
786 l.Draw();
787 }
788}
789
790// --------------------------------------------------------------------------
791//
792// Draws the two given histograms. Uses title as the pad title.
793// Also layout the two statistic boxes and a legend.
794//
795void MH::DrawSame(TH1 &hist1, TH1 &hist2, const TString title)
796{
797 //
798 // Draw first histogram
799 //
800 hist1.Draw();
801 gPad->SetBorderMode(0);
802 gPad->Update();
803
804 if (hist1.GetEntries()>0 && hist2.GetEntries()>0)
805 {
806 const Double_t maxbin1 = hist1.GetBinContent(hist1.GetMaximumBin());
807 const Double_t maxbin2 = hist2.GetBinContent(hist2.GetMaximumBin());
808 const Double_t minbin1 = hist1.GetBinContent(hist1.GetMinimumBin());
809 const Double_t minbin2 = hist2.GetBinContent(hist2.GetMinimumBin());
810
811 const Double_t max = TMath::Max(maxbin1, maxbin2);
812 const Double_t min = TMath::Min(minbin1, minbin2);
813
814 if (max!=min)
815 {
816 hist1.SetMaximum(max>0?max*1.05:max*0.95);
817 hist1.SetMinimum(max>0?min*0.95:min*1.05);
818 }
819 }
820
821 TPaveText *t = (TPaveText*)gPad->FindObject("title");
822 if (t)
823 {
824 t->SetName((TString)"MHTitle"); // rename object
825 t->Clear(); // clear old lines
826 t->AddText((TString)" "+title+" "); // add the new title
827 t->SetBit(kCanDelete); // make sure object is deleted
828
829 //
830 // FIXME: This is a stupid workaround to hide the redrawn
831 // (see THistPainter::PaintTitle) title
832 //
833 gPad->Modified(); // indicates a change
834 gPad->Update(); // recreates the original title
835 t->Pop(); // bring our title on top
836 }
837
838 //
839 // Rename first statistics box
840 //
841 // Where to get the TPaveStats depends on the root version
842 TPaveStats *s1 = (TPaveStats*)gPad->FindObject("stats");
843 if (!s1)
844 s1 = (TPaveStats*)hist1.GetListOfFunctions()->FindObject("stats");
845 else
846 s1->SetName((TString)"Stat"+hist1.GetTitle());
847
848 if (s1 && s1->GetX2NDC()>0.95)
849 {
850 const Double_t x1 = s1->GetX1NDC()-0.01;
851 s1->SetX1NDC(x1-(s1->GetX2NDC()-s1->GetX1NDC()));
852 s1->SetX2NDC(x1);
853 }
854
855 //
856 // Draw second histogram
857 //
858 hist2.Draw("sames");
859 gPad->Update();
860
861 //
862 // Draw Legend
863 //
864 // Where to get the TPaveStats depends on the root version
865 TPaveStats *s2 = (TPaveStats*)gPad->FindObject("stats");
866 if (!s2)
867 s2 = (TPaveStats*)hist2.GetListOfFunctions()->FindObject("stats");
868
869 if (s2)
870 {
871 TLegend &l = *new TLegend(s2->GetX1NDC(),
872 s2->GetY1NDC()-0.015-(s2->GetY2NDC()-s2->GetY1NDC())/2,
873 s2->GetX2NDC(),
874 s2->GetY1NDC()-0.01
875 );
876 l.AddEntry(&hist1, hist1.GetTitle());
877 l.AddEntry(&hist2, hist2.GetTitle());
878 l.SetTextSize(s2->GetTextSize());
879 l.SetTextFont(s2->GetTextFont());
880 l.SetBorderSize(s2->GetBorderSize());
881 l.SetBit(kCanDelete);
882 l.Draw();
883 }
884}
885
886// --------------------------------------------------------------------------
887//
888// If the opt string contains 'nonew' or gPad is not given NULL is returned.
889// Otherwise the present gPad is returned.
890//
891TVirtualPad *MH::GetNewPad(TString &opt)
892{
893 opt.ToLower();
894
895 if (!opt.Contains("nonew"))
896 return NULL;
897
898 opt.ReplaceAll("nonew", "");
899
900 return gPad;
901}
902
903// --------------------------------------------------------------------------
904//
905// Encapsulate the TObject::Clone such, that a cloned TH1 (or derived)
906// object is not added to the current directory, when cloned.
907//
908TObject *MH::Clone(const char *name) const
909{
910 const Bool_t store = TH1::AddDirectoryStatus();
911
912 TH1::AddDirectory(kFALSE);
913 TObject *o = MParContainer::Clone(name);
914 TH1::AddDirectory(store);
915
916 return o;
917}
918
919// --------------------------------------------------------------------------
920//
921// If the opt string contains 'nonew' or gPad is not given a new canvas
922// with size w/h is created. Otherwise the object is cloned and drawn
923// to the present pad. The kCanDelete bit is set for the clone.
924//
925TObject *MH::DrawClone(Option_t *opt, Int_t w, Int_t h) const
926{
927 TString option(opt);
928
929 TVirtualPad *p = GetNewPad(option);
930 if (!p)
931 p = MakeDefCanvas(this, w, h);
932 else
933 if (!option.Contains("same", TString::kIgnoreCase))
934 p->Clear();
935
936 gROOT->SetSelectedPad(NULL);
937
938 TObject *o = MParContainer::DrawClone(option);
939 o->SetBit(kCanDelete);
940 return o;
941}
942
943// --------------------------------------------------------------------------
944//
945// Check whether a class inheriting from MH overwrites the Draw function
946//
947Bool_t MH::OverwritesDraw(TClass *cls) const
948{
949 if (!cls)
950 cls = IsA();
951
952 //
953 // Check whether we reached the base class MTask
954 //
955 if (TString(cls->GetName())=="MH")
956 return kFALSE;
957
958 //
959 // Check whether the class cls overwrites Draw
960 //
961 if (cls->GetMethodAny("Draw"))
962 return kTRUE;
963
964 //
965 // If the class itself doesn't overload it check all it's base classes
966 //
967 TBaseClass *base=NULL;
968 TIter NextBase(cls->GetListOfBases());
969 while ((base=(TBaseClass*)NextBase()))
970 {
971 if (OverwritesDraw(base->GetClassPointer()))
972 return kTRUE;
973 }
974
975 return kFALSE;
976}
977
978// --------------------------------------------------------------------------
979//
980// Cuts the bins containing only zeros at the edges.
981//
982// A new number of bins can be defined with nbins != 0
983// In the case of nbins == 0, no rebinning will take place
984//
985// Returns the new (real) number of bins
986//
987Int_t MH::StripZeros(TH1 *h, Int_t nbins)
988{
989 TAxis &axe = *h->GetXaxis();
990
991 const Int_t min1 = axe.GetFirst();
992 const Int_t max1 = axe.GetLast();
993 const Int_t range1 = max1-min1;
994
995 //
996 // Check for useless zeros
997 //
998 if (range1 == 0)
999 return 0;
1000
1001 Int_t min2 = 0;
1002 for (int i=min1; i<=max1; i++)
1003 if (h->GetBinContent(i) != 0)
1004 {
1005 min2 = i;
1006 break;
1007 }
1008
1009 //
1010 // If the histogram consists of zeros only
1011 //
1012 if (min2 == max1)
1013 return 0;
1014
1015 Int_t max2 = 0;
1016 for (int i=max1; i>=min2; i--)
1017 if (h->GetBinContent(i) != 0)
1018 {
1019 max2 = i;
1020 break;
1021 }
1022
1023 //
1024 // Appying TAxis->SetRange before ReBin does not work ...
1025 // But this workaround helps quite fine
1026 //
1027 Axis_t min = h->GetBinLowEdge(min2);
1028 Axis_t max = h->GetBinLowEdge(max2)+h->GetBinWidth(max2);
1029
1030 Int_t nbins2 = max2-min2;
1031 //
1032 // Check for rebinning
1033 //
1034 if (nbins > 0)
1035 {
1036 const Int_t ngroup = (Int_t)(nbins2*h->GetNbinsX()/nbins/(max1-min1));
1037 if (ngroup > 1)
1038 {
1039 h->Rebin(ngroup);
1040 nbins2 /= ngroup;
1041 }
1042 }
1043
1044 Int_t newbins = 0;
1045 FindGoodLimits(nbins2, newbins, min, max, kFALSE);
1046 axe.SetRangeUser(min,max);
1047 return axe.GetLast()-axe.GetFirst();
1048}
1049
1050void MH::ProjectionX(TH1D &dest, const TH2 &src, Int_t firstybin, Int_t lastybin)
1051{
1052 //*-*-*-*-*Project a 2-D histogram into a 1-D histogram along X*-*-*-*-*-*-*
1053 //*-* ====================================================
1054 //
1055 // The projection dest is always of the type TH1D.
1056 // The projection is made from the channels along the Y axis
1057 // ranging from firstybin to lastybin included.
1058 // By default, bins 1 to ny are included
1059 // When all bins are included, the number of entries in the projection
1060 // is set to the number of entries of the 2-D histogram, otherwise
1061 // the number of entries is incremented by 1 for all non empty cells.
1062 //
1063 // if Sumw2() was called for dest, the errors are computed.
1064 //
1065 TAxis &axex = *((TH2&)src).GetXaxis();
1066 TAxis &axey = *((TH2&)src).GetYaxis();
1067
1068 const Int_t nx = axex.GetNbins();
1069 const Int_t ny = axey.GetNbins();
1070 if (firstybin < 0)
1071 firstybin = 1;
1072 if (lastybin > ny)
1073 lastybin = ny;
1074
1075 dest.Reset();
1076 SetBinning(&dest, &axex);
1077
1078 // Create the projection histogram
1079 const Bool_t computeErrors = dest.GetSumw2N() ? 1 : 0;
1080
1081 // Fill the projected histogram
1082 for (Int_t binx=0; binx<=nx+1; binx++)
1083 {
1084 Double_t err2 = 0;
1085 for (Int_t biny=firstybin; biny<=lastybin; biny++)
1086 {
1087 const Double_t cont = src.GetCellContent(binx,biny);
1088 const Double_t err1 = src.GetCellError(binx,biny);
1089 err2 += err1*err1;
1090 if (cont)
1091 dest.Fill(axex.GetBinCenter(binx), cont);
1092 }
1093 if (computeErrors)
1094 dest.SetBinError(binx, TMath::Sqrt(err2));
1095 }
1096 if (firstybin <=1 && lastybin >= ny)
1097 dest.SetEntries(src.GetEntries());
1098}
1099
1100void MH::ProjectionY(TH1D &dest, const TH2 &src, Int_t firstxbin, Int_t lastxbin)
1101{
1102 //*-*-*-*-*Project a 2-D histogram into a 1-D histogram along X*-*-*-*-*-*-*
1103 //*-* ====================================================
1104 //
1105 // The projection dest is always of the type TH1D.
1106 // The projection is made from the channels along the Y axis
1107 // ranging from firstybin to lastybin included.
1108 // By default, bins 1 to ny are included
1109 // When all bins are included, the number of entries in the projection
1110 // is set to the number of entries of the 2-D histogram, otherwise
1111 // the number of entries is incremented by 1 for all non empty cells.
1112 //
1113 // if Sumw2() was called for dest, the errors are computed.
1114 //
1115 TAxis &axex = *((TH2&)src).GetXaxis();
1116 TAxis &axey = *((TH2&)src).GetYaxis();
1117
1118 const Int_t nx = axex.GetNbins();
1119 const Int_t ny = axey.GetNbins();
1120 if (firstxbin < 0)
1121 firstxbin = 1;
1122 if (lastxbin > nx)
1123 lastxbin = nx;
1124
1125 dest.Reset();
1126 SetBinning(&dest, &axey);
1127
1128 // Create the projection histogram
1129 const Bool_t computeErrors = dest.GetSumw2N() ? 1 : 0;
1130
1131 // Fill the projected histogram
1132 for (Int_t biny=0; biny<=ny+1; biny++)
1133 {
1134 Double_t err2 = 0;
1135 for (Int_t binx=firstxbin; binx<=lastxbin; binx++)
1136 {
1137 const Double_t cont = src.GetCellContent(binx,biny);
1138 const Double_t err1 = src.GetCellError(binx,biny);
1139 err2 += err1*err1;
1140 if (cont)
1141 dest.Fill(axey.GetBinCenter(biny), cont);
1142 }
1143 if (computeErrors)
1144 dest.SetBinError(biny, TMath::Sqrt(err2));
1145 }
1146 if (firstxbin <=1 && lastxbin >= nx)
1147 dest.SetEntries(src.GetEntries());
1148}
1149
1150// --------------------------------------------------------------------------
1151//
1152// In contradiction to TPad::FindObject this function searches recursively
1153// in a pad for an object. gPad is the default.
1154//
1155TObject *MH::FindObjectInPad(const char *name, TVirtualPad *pad)
1156{
1157 if (!pad)
1158 pad = gPad;
1159
1160 if (!pad)
1161 return NULL;
1162
1163 TObject *o;
1164
1165 TIter Next(pad->GetListOfPrimitives());
1166 while ((o=Next()))
1167 {
1168 if (!strcmp(o->GetName(), name))
1169 return o;
1170
1171 if (o->InheritsFrom("TPad"))
1172 if ((o = FindObjectInPad(name, (TVirtualPad*)o)))
1173 return o;
1174 }
1175 return NULL;
1176}
1177
1178// --------------------------------------------------------------------------
1179//
1180// M.Gaug added this withouz Documentation
1181//
1182TH1I* MH::ProjectArray(const TArrayF &array, Int_t nbins,
1183 const char* name, const char* title)
1184{
1185 const Int_t size = array.GetSize();
1186
1187 TH1I *h1=0;
1188
1189 //check if histogram with identical name exist
1190 TObject *h1obj = gROOT->FindObject(name);
1191 if (h1obj && h1obj->InheritsFrom("TH1I"))
1192 {
1193 h1 = (TH1I*)h1obj;
1194 h1->Reset();
1195 }
1196
1197 Double_t min = size>0 ? array[0] : 0;
1198 Double_t max = size>0 ? array[0] : 1;
1199
1200 // first loop over array to find the min and max
1201 for (Int_t i=1; i<size;i++)
1202 {
1203 max = TMath::Max((Double_t)array[i], max);
1204 min = TMath::Min((Double_t)array[i], min);
1205 }
1206
1207 Int_t newbins = 0;
1208 FindGoodLimits(nbins, newbins, min, max, kFALSE);
1209
1210 if (!h1)
1211 {
1212 h1 = new TH1I(name, title, nbins, min, max);
1213 h1->SetXTitle("");
1214 h1->SetYTitle("Counts");
1215 h1->SetDirectory(NULL);
1216 }
1217
1218 // Second loop to fill the histogram
1219 for (Int_t i=0;i<size;i++)
1220 h1->Fill(array[i]);
1221
1222 return h1;
1223}
1224
1225// --------------------------------------------------------------------------
1226//
1227// M.Gaug added this withouz Documentation
1228//
1229TH1I* MH::ProjectArray(const TArrayD &array, Int_t nbins, const char* name, const char* title)
1230{
1231 const Int_t size = array.GetSize();
1232
1233 Double_t min = size>0 ? array[0] : 0;
1234 Double_t max = size>0 ? array[0] : 1;
1235
1236 TH1I *h1=0;
1237
1238 //check if histogram with identical name exist
1239 TObject *h1obj = gROOT->FindObject(name);
1240 if (h1obj && h1obj->InheritsFrom("TH1I"))
1241 {
1242 h1 = (TH1I*)h1obj;
1243 h1->Reset();
1244 }
1245
1246 // first loop over array to find the min and max
1247 for (Int_t i=1; i<size;i++)
1248 {
1249 max = TMath::Max(array[i], max);
1250 min = TMath::Min(array[i], min);
1251 }
1252
1253 Int_t newbins = 0;
1254 FindGoodLimits(nbins, newbins, min, max, kFALSE);
1255
1256 if (!h1)
1257 {
1258 h1 = new TH1I(name, title, newbins, min, max);
1259 h1->SetXTitle("");
1260 h1->SetYTitle("Counts");
1261 h1->SetDirectory(NULL);
1262 }
1263
1264 // Second loop to fill the histogram
1265 for (Int_t i=0;i<size;i++)
1266 h1->Fill(array[i]);
1267
1268 return h1;
1269}
1270
1271// --------------------------------------------------------------------------
1272//
1273// M.Gaug added this withouz Documentation
1274//
1275TH1I* MH::ProjectArray(const MArrayF &array, Int_t nbins,
1276 const char* name, const char* title)
1277{
1278 return ProjectArray(TArrayF(array.GetSize(),array.GetArray()), nbins, name, title);
1279}
1280
1281// --------------------------------------------------------------------------
1282//
1283// M.Gaug added this withouz Documentation
1284//
1285TH1I* MH::ProjectArray(const MArrayD &array, Int_t nbins, const char* name, const char* title)
1286{
1287 return ProjectArray(TArrayD(array.GetSize(),array.GetArray()), nbins, name, title);
1288}
1289
1290// --------------------------------------------------------------------------
1291//
1292// See MTask::PrintSkipped
1293//
1294void MH::PrintSkipped(UInt_t n, const char *str)
1295{
1296 *fLog << " " << setw(7) << n << " (";
1297 *fLog << Form("%5.1f", 100.*n/GetNumExecutions());
1298 *fLog << "%) Evts skipped: " << str << endl;
1299}
Note: See TracBrowser for help on using the repository browser.