source: tags/Mars-V2.1/mhistmc/MHMcEnergy.cc

Last change on this file was 2173, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 6.5 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): Javier Lopez 05/2001 <mailto:jlopez@ifae.es>
19! Author(s): Thomas Bretz 05/2001 <mailto:tbretz@uni-sw.gwdg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2001
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MHMcEnergy
29//
30// This class holds the information (histogram and fit function)
31// about the energy threshold for a particular trigger condition.
32//
33////////////////////////////////////////////////////////////////////////////
34#include "MHMcEnergy.h"
35
36#include <stdlib.h>
37#include <iostream>
38
39#include <TH1.h>
40#include <TF1.h>
41#include <TCanvas.h>
42#include <TPaveLabel.h>
43
44#include "MH.h"
45
46#include "MLog.h"
47#include "MLogManip.h"
48
49ClassImp(MHMcEnergy);
50
51using namespace std;
52
53// -------------------------------------------------------------------------
54//
55// Default Constructor.
56//
57MHMcEnergy::MHMcEnergy(const char *name, const char *title)
58{
59 fTitle = title ? title : "Container for an energy distribution histogram";
60
61 // - we initialize the histogram
62 // - we have to give diferent names for the diferent histograms because
63 // root don't allow us to have diferent histograms with the same name
64
65 fHist = new TH1F("", "", 20, 0.5, 4.5);
66
67 fHist->SetDirectory(NULL);
68 fHist->SetXTitle("log(E/GeV)");
69 fHist->SetYTitle("dN/dE");
70
71 SetName(name ? name : "MHMcEnergy");
72}
73
74// -------------------------------------------------------------------------
75//
76// This doesn't only set the name. It tries to get the number from the
77// name and creates also name and title of the histogram.
78//
79// This is necessary for example if a list of such MHMcEnergy histograms
80// is created (s. MParList::CreateObjList)
81//
82void MHMcEnergy::SetName(const char *name)
83{
84 TString cname(name);
85 const char *semicolon = strrchr(cname, ';');
86
87 UInt_t idx = semicolon ? atoi(semicolon+1) : 0;
88
89 fName = cname;
90
91 char text[256];
92 if (idx>0)
93 sprintf(text, "Energy Distribution for trigger condition #%i", idx);
94 else
95 sprintf(text, "Energy Distribution");
96
97 char aux[256];
98 strcpy(aux, "Threshold");
99
100 if (idx>0)
101 sprintf(aux+strlen(aux), " #%i", idx);
102
103 fHist->SetName(aux);
104 fHist->SetTitle(text);
105}
106
107//-------------------------------------------------------------------------
108//
109// Defualt Destructor
110//
111MHMcEnergy::~MHMcEnergy()
112{
113 delete fHist;
114}
115
116//--------------------------------------------------------------------------
117//
118// Fill the histogram with the log10 of the energy for triggered events.
119//
120void MHMcEnergy::Fill(Float_t log10E, Float_t w)
121{
122 fHist->Fill(log10E, w);
123}
124
125// -------------------------------------------------------------------------
126//
127// Fitting function
128//
129void MHMcEnergy::Fit(Axis_t xxmin, Axis_t xxmax)
130{
131 //
132 // 0: don't draw the function (it is drawn together with the histogram)
133 // Q: quiet mode
134 //
135 fHist->Fit("gaus", "Q0", "", xxmin, xxmax);
136
137 TF1 *result = fHist->GetFunction("gaus");
138
139 fThreshold = CalcThreshold(result);
140 fThresholdErr = CalcThresholdErr(result);
141 fGaussPeak = CalcGaussPeak(result);
142 fGaussSigma = CalcGaussSigma(result);
143}
144
145// ------------------------------------------------------------------------
146//
147// Helper function for Draw() and DrawClone() which adds some useful
148// information to the plot.
149//
150void MHMcEnergy::DrawLegend() const
151{
152 char text[256];
153
154 const Float_t min = fHist->GetMinimum();
155 const Float_t max = fHist->GetMaximum();
156 const Float_t sum = min+max;
157
158 sprintf(text, "Energy Threshold = %4.1f +- %4.1f GeV",
159 fThreshold, fThresholdErr);
160
161 TPaveLabel* label = new TPaveLabel(2.2, 0.75*sum, 4.4, 0.90*sum, text);
162 label->SetFillColor(10);
163 label->SetTextSize(0.3);
164 label->SetBit(kCanDelete);
165 label->Draw();
166}
167
168// ------------------------------------------------------------------------
169//
170// Drawing function. It creates its own canvas.
171//
172void MHMcEnergy::Draw(Option_t *option)
173{
174 TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);
175 pad->SetBorderMode(0);
176
177 AppendPad("");
178
179 fHist->Draw(option);
180
181 DrawLegend();
182
183 pad->Modified();
184 pad->Update();
185}
186
187// --------------------------------------------------------------------------
188//
189// Set the number of bins in the histogran.
190//
191void MHMcEnergy::SetNumBins(Int_t nbins)
192{
193 fHist->SetBins(nbins, 0.5, 4.5);
194}
195// --------------------------------------------------------------------------
196//
197// Write the threshold and its error in the standard output
198//
199void MHMcEnergy::Print(Option_t*) const
200{
201 *fLog << all << "Threshold: " << fThreshold << " +- " << fThresholdErr << endl;
202}
203
204// -------------------------------------------------------------------------
205//
206// Return the threshold
207//
208Float_t MHMcEnergy::CalcThreshold(TF1 *gauss)
209{
210 const Float_t p1 = gauss->GetParameter(1);
211
212 return pow(10, p1);
213}
214
215// -------------------------------------------------------------------------
216//
217// Return the error of the threshold.
218//
219Float_t MHMcEnergy::CalcThresholdErr(TF1 *gauss)
220{
221 const Float_t lg10 = log(10.);
222 const Float_t p1 = gauss->GetParameter(1);
223 const Float_t p1err = gauss->GetParError(1);
224
225 // The error has into accuont the error in the fit
226 return pow(10, p1) * p1err * lg10;
227}
228
229// -------------------------------------------------------------------------
230//
231// Return the peak of the fitted gaussan function.
232//
233Float_t MHMcEnergy::CalcGaussPeak(TF1 *gauss)
234{
235 return gauss->GetParameter(1);
236}
237
238// -------------------------------------------------------------------------
239//
240// Return the sigma of the fitted gaussan function.
241//
242Float_t MHMcEnergy::CalcGaussSigma(TF1 *gauss)
243{
244 return gauss->GetParameter(2);
245}
246
247TH1 *MHMcEnergy::GetHistByName(const TString name)
248{
249 return fHist;
250}
Note: See TracBrowser for help on using the repository browser.