source: trunk/MagicSoft/Mars/mhist/MH.cc@ 1505

Last change on this file since 1505 was 1505, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 19.2 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//////////////////////////////////////////////////////////////////////////////
48
49#include "MH.h"
50
51#include <TH1.h>
52#include <TH2.h>
53#include <TH3.h>
54#include <TGaxis.h>
55#include <TCanvas.h>
56#include <TLegend.h>
57#include <TPaveStats.h>
58
59#include "MLog.h"
60#include "MLogManip.h"
61
62#include "MParList.h"
63#include "MParContainer.h"
64
65#include "MBinning.h"
66
67ClassImp(MH);
68
69// --------------------------------------------------------------------------
70//
71// Default Constructor. It sets name and title only. Typically you won't
72// need to change this.
73//
74MH::MH(const char *name, const char *title)
75{
76 //
77 // set the name and title of this object
78 //
79 fName = name ? name : "MH";
80 fTitle = title ? title : "Base class for Mars histograms";
81}
82
83// --------------------------------------------------------------------------
84//
85// If you want to use the automatic filling of your derived class you
86// must overload this function. If it is not overloaded you cannot use
87// FillH with this class. The argument is a pointer to a container
88// in your paremeter list which is specified in the MFillH constructor.
89// If you are not going to use it you should at least add
90// Bool_t MH::Fill(const MParContainer *) { return kTRUE; }
91// to your class definition.
92//
93Bool_t MH::Fill(const MParContainer *par)
94{
95 *fLog << GetDescriptor() << ": Fill not overloaded! Can't be used!" << endl;
96 return kFALSE;
97}
98
99// --------------------------------------------------------------------------
100//
101// This is a function which should replace the creation of default
102// canvases like root does. Because this is inconvinient in some aspects.
103// need to change this.
104// You can specify a name for the default canvas and a title. Also
105// width and height can be given.
106// MakeDefCanvas looks for a canvas with the given name. If now name is
107// given the DefCanvasName of root is used. If no such canvas is existing
108// it is created and returned. If such a canvas already exists a new canvas
109// with a name plus anumber is created (the number is calculated by the
110// number of all existing canvases plus one)
111//
112TCanvas *MH::MakeDefCanvas(TString name, const char *title,
113 const UInt_t w, const UInt_t h)
114{
115 const TList *list = (TList*)gROOT->GetListOfCanvases();
116
117 if (name.IsNull())
118 name = gROOT->GetDefCanvasName();
119
120 if (list->FindObject(name))
121 name += Form(" <%d>", list->GetSize()+1);
122
123 return new TCanvas(name, title, w, h);
124}
125
126// --------------------------------------------------------------------------
127//
128// This function works like MakeDefCanvas(name, title, w, h) but name
129// and title are retrieved from the given TObject.
130//
131TCanvas *MH::MakeDefCanvas(const TObject *obj,
132 const UInt_t w, const UInt_t h)
133{
134 return MakeDefCanvas(obj->GetName(), obj->GetTitle(), w, h);
135}
136
137// --------------------------------------------------------------------------
138//
139// Applies a given binning to a 1D-histogram
140//
141void MH::SetBinning(TH1 *h, const MBinning *binsx)
142{
143 //
144 // Another strange behaviour: TAxis::Set deletes the axis title!
145 //
146 TAxis &x = *h->GetXaxis();
147
148#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
149 TString xtitle = x.GetTitle();
150#endif
151
152 //
153 // This is a necessary workaround if one wants to set
154 // non-equidistant bins after the initialization
155 // TH1D::fNcells must be set correctly.
156 //
157 h->SetBins(binsx->GetNumBins(), 0, 1);
158
159 //
160 // Set the binning of the current histogram to the binning
161 // in one of the two given histograms
162 //
163 x.Set(binsx->GetNumBins(), binsx->GetEdges());
164#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
165 x.SetTitle(xtitle);
166#endif
167}
168
169// --------------------------------------------------------------------------
170//
171// Applies given binnings to the two axis of a 2D-histogram
172//
173void MH::SetBinning(TH2 *h, const MBinning *binsx, const MBinning *binsy)
174{
175 TAxis &x = *h->GetXaxis();
176 TAxis &y = *h->GetYaxis();
177
178 //
179 // Another strange behaviour: TAxis::Set deletes the axis title!
180 //
181#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
182 TString xtitle = x.GetTitle();
183 TString ytitle = y.GetTitle();
184#endif
185
186 //
187 // This is a necessary workaround if one wants to set
188 // non-equidistant bins after the initialization
189 // TH1D::fNcells must be set correctly.
190 //
191 h->SetBins(binsx->GetNumBins(), 0, 1,
192 binsy->GetNumBins(), 0, 1);
193
194 //
195 // Set the binning of the current histogram to the binning
196 // in one of the two given histograms
197 //
198 x.Set(binsx->GetNumBins(), binsx->GetEdges());
199 y.Set(binsy->GetNumBins(), binsy->GetEdges());
200#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
201 x.SetTitle(xtitle);
202 y.SetTitle(ytitle);
203#endif
204}
205
206// --------------------------------------------------------------------------
207//
208// Applies given binnings to the three axis of a 3D-histogram
209//
210void MH::SetBinning(TH3 *h, const MBinning *binsx, const MBinning *binsy, const MBinning *binsz)
211{
212 //
213 // Another strange behaviour: TAxis::Set deletes the axis title!
214 //
215 TAxis &x = *h->GetXaxis();
216 TAxis &y = *h->GetYaxis();
217 TAxis &z = *h->GetZaxis();
218
219#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
220 TString xtitle = x.GetTitle();
221 TString ytitle = y.GetTitle();
222 TString ztitle = z.GetTitle();
223#endif
224
225 //
226 // This is a necessary workaround if one wants to set
227 // non-equidistant bins after the initialization
228 // TH1D::fNcells must be set correctly.
229 //
230 h->SetBins(binsx->GetNumBins(), 0, 1,
231 binsy->GetNumBins(), 0, 1,
232 binsz->GetNumBins(), 0, 1);
233
234 //
235 // Set the binning of the current histogram to the binning
236 // in one of the two given histograms
237 //
238 x.Set(binsx->GetNumBins(), binsx->GetEdges());
239 y.Set(binsy->GetNumBins(), binsy->GetEdges());
240 z.Set(binsz->GetNumBins(), binsz->GetEdges());
241#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
242 x.SetTitle(xtitle);
243 y.SetTitle(ytitle);
244 z.SetTitle(ztitle);
245#endif
246}
247
248// --------------------------------------------------------------------------
249//
250// Applies given binning (the n+1 edges) to the axis of a 1D-histogram
251//
252void MH::SetBinning(TH1 *h, const TArrayD &binsx)
253{
254 MBinning bx;
255 bx.SetEdges(binsx);
256 SetBinning(h, &bx);
257}
258
259// --------------------------------------------------------------------------
260//
261// Applies given binning (the n+1 edges) to the two axis of a
262// 2D-histogram
263//
264void MH::SetBinning(TH2 *h, const TArrayD &binsx, const TArrayD &binsy)
265{
266 MBinning bx;
267 MBinning by;
268 bx.SetEdges(binsx);
269 by.SetEdges(binsy);
270 SetBinning(h, &bx, &by);
271}
272
273// --------------------------------------------------------------------------
274//
275// Applies given binning (the n+1 edges) to the three axis of a
276// 3D-histogram
277//
278void MH::SetBinning(TH3 *h, const TArrayD &binsx, const TArrayD &binsy, const TArrayD &binsz)
279{
280 MBinning bx;
281 MBinning by;
282 MBinning bz;
283 bx.SetEdges(binsx);
284 by.SetEdges(binsy);
285 bz.SetEdges(binsz);
286 SetBinning(h, &bx, &by, &bz);
287}
288
289// --------------------------------------------------------------------------
290//
291// Applies the binning of a TAxis (eg from a root histogram) to the axis
292// of a 1D-histogram
293//
294void MH::SetBinning(TH1 *h, const TAxis *binsx)
295{
296 const Int_t nx = binsx->GetNbins();
297
298 TArrayD bx(nx+1);
299 for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
300 bx[nx] = binsx->GetXmax();
301
302 SetBinning(h, bx);
303}
304
305// --------------------------------------------------------------------------
306//
307// Applies the binnings of the TAxis' (eg from a root histogram) to the
308// two axis' of a 2D-histogram
309//
310void MH::SetBinning(TH2 *h, const TAxis *binsx, const TAxis *binsy)
311{
312 const Int_t nx = binsx->GetNbins();
313 const Int_t ny = binsy->GetNbins();
314
315 TArrayD bx(nx+1);
316 TArrayD by(ny+1);
317 for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
318 for (int i=0; i<ny; i++) by[i] = binsy->GetBinLowEdge(i+1);
319 bx[nx] = binsx->GetXmax();
320 by[ny] = binsy->GetXmax();
321
322 SetBinning(h, bx, by);
323}
324
325// --------------------------------------------------------------------------
326//
327// Applies the binnings of the TAxis' (eg from a root histogram) to the
328// three axis' of a 3D-histogram
329//
330void MH::SetBinning(TH3 *h, const TAxis *binsx, const TAxis *binsy, const TAxis *binsz)
331{
332 const Int_t nx = binsx->GetNbins();
333 const Int_t ny = binsy->GetNbins();
334 const Int_t nz = binsz->GetNbins();
335
336 TArrayD bx(nx+1);
337 TArrayD by(ny+1);
338 TArrayD bz(nz+1);
339 for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
340 for (int i=0; i<ny; i++) by[i] = binsy->GetBinLowEdge(i+1);
341 for (int i=0; i<nz; i++) bz[i] = binsz->GetBinLowEdge(i+1);
342 bx[nx] = binsx->GetXmax();
343 by[ny] = binsy->GetXmax();
344 bz[nz] = binsz->GetXmax();
345
346 SetBinning(h, bx, by, bz);
347}
348
349// --------------------------------------------------------------------------
350//
351// Applies the binnings of one root-histogram x to another one h
352// Both histograms must be of the same type: TH1, TH2 or TH3
353//
354void MH::SetBinning(TH1 *h, TH1 *x)
355{
356 if (h->InheritsFrom(TH3::Class()) && x->InheritsFrom(TH3::Class()))
357 {
358 SetBinning((TH3*)h, x->GetXaxis(), x->GetYaxis(), x->GetZaxis());
359 return;
360 }
361 if (h->InheritsFrom(TH3::Class()) || x->InheritsFrom(TH3::Class()))
362 return;
363 if (h->InheritsFrom(TH2::Class()) && x->InheritsFrom(TH2::Class()))
364 {
365 SetBinning((TH2*)h, x->GetXaxis(), x->GetYaxis());
366 return;
367 }
368 if (h->InheritsFrom(TH2::Class()) || x->InheritsFrom(TH2::Class()))
369 return;
370 if (h->InheritsFrom(TH1::Class()) && x->InheritsFrom(TH1::Class()))
371 {
372 SetBinning(h, x->GetXaxis());
373 return;
374 }
375}
376
377// --------------------------------------------------------------------------
378//
379// Multiplies all entries in a TArrayD by a float f
380//
381void MH::ScaleArray(TArrayD &bins, Double_t f)
382{
383 for (int i=0; i<bins.GetSize(); i++)
384 bins[i] *= f;
385}
386
387// --------------------------------------------------------------------------
388//
389// Scales the binning of a TAxis by a float f
390//
391TArrayD MH::ScaleAxis(TAxis &axe, Double_t f)
392{
393 TArrayD arr(axe.GetNbins()+1);
394
395 for (int i=1; i<=axe.GetNbins()+1; i++)
396 arr[i-1] = axe.GetBinLowEdge(i);
397
398 ScaleArray(arr, f);
399
400 return arr;
401}
402
403// --------------------------------------------------------------------------
404//
405// Scales the binning of one, two or three axis of a histogram by a float f
406//
407void MH::ScaleAxis(TH1 *h, Double_t fx, Double_t fy, Double_t fz)
408{
409 if (h->InheritsFrom(TH3::Class()))
410 {
411 SetBinning((TH3*)h,
412 ScaleAxis(*h->GetXaxis(), fx),
413 ScaleAxis(*h->GetYaxis(), fy),
414 ScaleAxis(*h->GetZaxis(), fz));
415 return;
416 }
417
418 if (h->InheritsFrom(TH2::Class()))
419 {
420 SetBinning((TH2*)h,
421 ScaleAxis(*h->GetXaxis(), fx),
422 ScaleAxis(*h->GetYaxis(), fy));
423 return;
424 }
425
426 if (h->InheritsFrom(TH1::Class()))
427 SetBinning(h, ScaleAxis(*h->GetXaxis(), fx));
428}
429
430// --------------------------------------------------------------------------
431//
432// Tries to find a MBinning container with the name "Binning"+name
433// in the given parameter list. If it was found it is applied to the
434// given histogram. This is only valid for 1D-histograms
435//
436Bool_t MH::ApplyBinning(const MParList &plist, TString name, TH1 *h)
437{
438 if (h->InheritsFrom(TH2::Class()) || h->InheritsFrom(TH3::Class()))
439 {
440 gLog << warn << "MH::ApplyBinning: '" << h->GetName() << "' is not a basic TH1 object... no binning applied." << endl;
441 return kFALSE;
442 }
443
444 const MBinning *bins = (MBinning*)plist.FindObject("Binning"+name);
445 if (!bins)
446 {
447 gLog << warn << "Object 'Binning" << name << "' [MBinning] not found... no binning applied." << endl;
448 return kFALSE;
449 }
450
451 SetBinning(h, bins);
452 return kTRUE;
453}
454
455void MH::FindGoodLimits(Int_t nbins, Int_t &newbins, Double_t &xmin, Double_t &xmax, Bool_t isInteger)
456{
457//*-*-*-*-*-*-*-*-*Find reasonable bin values*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
458//*-* ==========================
459
460 Double_t dx = 0.1*(xmax-xmin);
461 Double_t umin = xmin - dx;
462 Double_t umax = xmax + dx;
463
464 if (umin < 0 && xmin >= 0)
465 umin = 0;
466
467 if (umax > 0 && xmax <= 0)
468 umax = 0;
469
470 Int_t n=0;
471 Double_t binlow =0;
472 Double_t binhigh =0;
473 Double_t binwidth=0;
474 TGaxis::Optimize(umin, umax, nbins, binlow, binhigh, n, binwidth, "");
475
476 if (binwidth <= 0 || binwidth > 1.e+39)
477 {
478 xmin = -1;
479 xmax = 1;
480 }
481 else
482 {
483 xmin = binlow;
484 xmax = binhigh;
485 }
486
487 if (isInteger)
488 {
489 Int_t ixmin = (Int_t)xmin;
490 Int_t ixmax = (Int_t)xmax;
491 Double_t dxmin = (Double_t)ixmin;
492 Double_t dxmax = (Double_t)ixmax;
493
494 xmin = xmin<0 && xmin!=dxmin ? dxmin - 1 : dxmin;
495 xmax = xmax>0 && xmax!=dxmax ? dxmax + 1 : dxmax;
496
497 if (xmin>=xmax)
498 xmax = xmin+1;
499
500 Int_t bw = 1 + (Int_t)((xmax-xmin)/nbins);
501
502 nbins = (Int_t)((xmax-xmin)/bw);
503
504 if (xmin+nbins*bw < xmax)
505 {
506 nbins++;
507 xmax = xmin +nbins*bw;
508 }
509 }
510
511 newbins = nbins;
512}
513
514// --------------------------------------------------------------------------
515//
516// Returns the lowest entry in a histogram which is greater than gt (eg >0)
517//
518Double_t MH::GetMinimumGT(const TH1 &h, Double_t gt)
519{
520 Double_t min = FLT_MAX;
521
522 const TAxis &axe = *((TH1&)h).GetXaxis();
523
524 for (int i=1; i<=axe.GetNbins(); i++)
525 {
526 Double_t x = h.GetBinContent(i);
527 if (gt<x && x<min)
528 min = x;
529 }
530 return min;
531}
532
533// --------------------------------------------------------------------------
534//
535// Draws a copy of the two given histograms. Uses title as the pad title.
536// Also layout the two statistic boxes and a legend.
537//
538void MH::DrawCopy(const TH1 &hist1, const TH1 &hist2, const TString title)
539{
540 //
541 // Draw first histogram
542 //
543 TH1 *h1 = (TH1*)((TH1&)hist1).DrawCopy();
544 gPad->Update();
545
546 //
547 // Rename first statistics box
548 //
549 TPaveStats &s1 = *(TPaveStats*)gPad->FindObject("stats");
550 s1.SetX1NDC(s1.GetX1NDC()-0.01);
551 s1.SetName((TString)"Stat"+hist1.GetTitle());
552
553 //
554 // Draw second histogram
555 //
556 TH1 *h2 = (TH1*)((TH1&)hist2).DrawCopy("sames");
557
558 gPad->Update();
559
560 //
561 // Set new position of second statistics box
562 //
563 TPaveStats &s2 = *(TPaveStats*)gPad->FindObject("stats");
564 s2.SetX1NDC(s1.GetX1NDC()-(s2.GetX2NDC()-s2.GetX1NDC())-0.01);
565 s2.SetX2NDC(s1.GetX1NDC()-0.01);
566
567 //
568 // Draw Legend
569 //
570 const Int_t n = s1.GetListOfLines()->GetSize();
571 const Double_t h = s1.GetY2NDC()-s1.GetY1NDC();
572 TLegend &l = *new TLegend(s1.GetX1NDC(),
573 s1.GetY1NDC()-0.015-h*2/n,
574 s1.GetX2NDC(),
575 s1.GetY1NDC()-0.01
576 );
577 l.AddEntry(h1, h1->GetTitle());
578 l.AddEntry(h2, h2->GetTitle());
579 l.SetTextSize(s2.GetTextSize());
580 l.SetTextFont(s2.GetTextFont());
581 l.SetBorderSize(s2.GetBorderSize());
582 l.SetBit(kCanDelete);
583
584 // FIXME: Don't change the hist title. Change TPaveText::title
585 h1->SetTitle(" "+title+" ");
586
587 l.Draw();
588}
589
590// --------------------------------------------------------------------------
591//
592// Draws the two given histograms. Uses title as the pad title.
593// Also layout the two statistic boxes and a legend.
594//
595void MH::Draw(TH1 &hist1, TH1 &hist2, const TString title)
596{
597 //
598 // Draw first histogram
599 //
600 hist1.Draw();
601 // FIXME: Don't change the hist title. Change TPaveText::title
602 //h1->SetTitle(" "+title+" ");
603
604 gPad->Update();
605
606 //
607 // Rename first statistics box
608 //
609 TPaveStats &s1 = *(TPaveStats*)gPad->FindObject("stats");
610 s1.SetX1NDC(s1.GetX1NDC()-0.01);
611 s1.SetName((TString)"Stat"+hist1.GetTitle());
612
613 //
614 // Draw second histogram
615 //
616 hist2.Draw("sames");
617
618 gPad->Update();
619
620 //
621 // Set new position of second statistics box
622 //
623 TPaveStats &s2 = *(TPaveStats*)gPad->FindObject("stats");
624 s2.SetX1NDC(s1.GetX1NDC()-(s2.GetX2NDC()-s2.GetX1NDC())-0.01);
625 s2.SetX2NDC(s1.GetX1NDC()-0.01);
626
627 //
628 // Draw Legend
629 //
630 TLegend &l = *new TLegend(s1.GetX1NDC(),
631 s1.GetY1NDC()-0.015-(s1.GetY2NDC()-s1.GetY1NDC())/2,
632 s1.GetX2NDC(),
633 s1.GetY1NDC()-0.01
634 );
635 l.AddEntry(&hist1, hist1.GetTitle());
636 l.AddEntry(&hist2, hist2.GetTitle());
637 l.SetTextSize(s2.GetTextSize());
638 l.SetTextFont(s2.GetTextFont());
639 l.SetBorderSize(s2.GetBorderSize());
640 l.SetBit(kCanDelete);
641 l.Draw();
642}
Note: See TracBrowser for help on using the repository browser.