source: trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.cc@ 4905

Last change on this file since 4905 was 2015, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 8.6 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 1/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Wolfgang Wittek 1/2002 <mailto:wittek@mppmu.mpg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2002
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27// //
28// MHEffOnTimeTheta //
29// //
30// calculates the effective on time for each bin in Theta //
31// //
32//////////////////////////////////////////////////////////////////////////////
33
34#include "MHEffOnTimeTheta.h"
35
36#include <TStyle.h>
37
38#include <TF1.h>
39#include <TH2.h>
40#include <TCanvas.h>
41
42#include "MTime.h"
43
44#include "MBinning.h"
45#include "MParList.h"
46
47#include "MLog.h"
48#include "MLogManip.h"
49
50ClassImp(MHEffOnTimeTheta);
51
52// --------------------------------------------------------------------------
53//
54// Default Constructor. It sets name and title of the histograms.
55//
56MHEffOnTimeTheta::MHEffOnTimeTheta(const char *name, const char *title)
57 : fHEffOn()
58{
59 //
60 // set the name and title of this object
61 //
62 fName = name ? name : "MHEffOnTimeTheta";
63 fTitle = title ? title : "1-D histogram of Eff On Time";
64
65 // effective on time versus theta
66 fHEffOn.SetName("EffOn");
67 fHEffOn.SetTitle("Effective On Time vs. Theta");
68
69 fHEffOn.SetDirectory(NULL);
70
71 fHEffOn.SetXTitle("\\Theta [\\circ]");
72 fHEffOn.SetYTitle("t_{eff} [s]");
73
74 // chi2/NDF versus theta
75 fHChi2.SetName("Chi2/NDF");
76 fHChi2.SetTitle("Chi2/NDF of OnTimeFit vs. Theta");
77
78 fHChi2.SetDirectory(NULL);
79
80 fHChi2.SetXTitle("\\Theta [\\circ]");
81 fHChi2.SetYTitle("chi^{2}/NDF");
82
83 // lambda versus theta
84 fHLambda.SetName("lambda");
85 fHLambda.SetTitle("lambda of OnTimeFit vs. Theta");
86
87 fHLambda.SetDirectory(NULL);
88
89 fHLambda.SetXTitle("\\Theta [\\circ]");
90 fHLambda.SetYTitle("\\lambda [Hz]");
91
92 // N0del versus theta
93 fHN0del.SetName("N0del");
94 fHN0del.SetTitle("N0del of OnTimeFit vs. Theta");
95
96 fHN0del.SetDirectory(NULL);
97
98 fHN0del.SetXTitle("\\Theta [\\circ]");
99 fHN0del.SetYTitle("N0del");
100}
101
102// -----------------------------------------------------------------------
103//
104// Calculate the effective on time by fitting the distribution of
105// time differences
106//
107void MHEffOnTimeTheta::Calc(TH2D *hist)
108{
109 // nbins = number of Theta bins
110 const Int_t nbins = hist->GetNbinsY();
111
112 for (int i=1; i<=nbins; i++)
113 {
114 // TH1D &h = *hist->ProjectionX("Calc-theta", i, i);
115 TH1D &h = *hist->ProjectionX("Calc-theta", i, i, "E");
116
117 // Nmdel = Nm * binwidth, with Nm = number of observed events
118 const Double_t Nmdel = h.Integral("width");
119 const Double_t Nm = h.Integral();
120 // Double_t mean = h->GetMean();
121
122 //...................................................
123 // determine range (yq[0], yq[1]) of time differences
124 // where fit should be performed;
125 // require a fraction >=xq[0] of all entries to lie below yq[0]
126 // and a fraction <=xq[1] of all entries to lie below yq[1];
127 // within the range (yq[0], yq[1]) there must be no empty bin;
128 // choose pedestrian approach as long as GetQuantiles is not available
129
130 Double_t xq[2] = { 0.15, 0.95 };
131 Double_t yq[2];
132
133 // GetQuantiles doesn't seem to be available in root 3.01/06
134 // h->GetQuantiles(2,yq,xq);
135
136 const Double_t sumtot = h.Integral();
137 const Int_t jbins = h.GetNbinsX();
138
139 if (sumtot > 0.0)
140 {
141 Double_t sum1 = 0.0;
142 yq[0] = h.GetBinLowEdge(jbins+1);
143 for (int j=1; j<=jbins; j++)
144 {
145 if (sum1 >= xq[0]*sumtot)
146 {
147 yq[0] = h.GetBinLowEdge(j);
148 break;
149 }
150 sum1 += h.GetBinContent(j);
151 }
152
153 Double_t sum2 = 0.0;
154 yq[1] = h.GetBinLowEdge(jbins+1);
155 for (int j=1; j<=jbins; j++)
156 {
157 Double_t content = h.GetBinContent(j);
158 sum2 += content;
159 if (sum2 >= xq[1]*sumtot || content == 0.0)
160 {
161 yq[1] = h.GetBinLowEdge(j);
162 break;
163 }
164 }
165
166 //...................................................
167
168 // parameter 0 = lambda
169 // parameter 1 = N0*del with N0 = ideal number of events
170 // and del = bin width of time difference
171 TF1 func("Poisson", "[1] * [0] * exp(-[0] *x)", yq[0], yq[1]);
172
173 func.SetParameter(0, 100); // [Hz]
174 func.SetParameter(1, Nmdel);
175
176 func.SetParLimits(0, 0, 1000); // [Hz]
177 func.SetParLimits(1, 0, 10*Nmdel);
178
179 func.SetParName(0, "lambda");
180 func.SetParName(1, "Nmdel");
181
182 // options : 0 do not plot the function
183 // I use integral of function in bin rather than value at bin center
184 // R use the range specified in the function range
185 // Q quiet mode
186 h.Fit("Poisson", "0IRQ");
187
188 // gPad->SetLogy();
189 // gStyle->SetOptFit(1011);
190 // h->GetXaxis()->SetTitle("time difference [s]");
191 // h->GetYaxis()->SetTitle("Counts");
192 // h->DrawCopy();
193
194 // func.SetRange(yq[0], yq[1]); // Range of Drawing
195 // func.DrawCopy("same");
196
197 const Double_t lambda = func.GetParameter(0);
198 const Double_t N0del = func.GetParameter(1);
199 const Double_t chi2 = func.GetChisquare();
200 const Int_t NDF = func.GetNDF();
201
202 // was fit successful ?
203 if (NDF>0 && chi2<2.5*NDF)
204 {
205 // the effective on time is Nm/lambda
206 fHEffOn.SetBinContent(i, Nm/lambda);
207
208 // plot chi2/NDF of fit
209 fHChi2.SetBinContent(i, NDF ? chi2/NDF : 0.0);
210
211 // lambda of fit
212 fHLambda.SetBinContent(i, lambda);
213
214 // N0del of fit
215 fHN0del.SetBinContent(i, N0del);
216
217 delete &h;
218 continue;
219 }
220 }
221
222 fHEffOn.SetBinContent (i, 1.e-20);
223 fHChi2.SetBinContent (i, 1.e-20);
224 fHLambda.SetBinContent(i, 1.e-20);
225 fHN0del.SetBinContent (i, 1.e-20);
226
227 delete &h;
228 }
229}
230
231// -------------------------------------------------------------------------
232//
233// Set the binnings and prepare the filling of the histograms
234//
235Bool_t MHEffOnTimeTheta::SetupFill(const MParList *plist)
236{
237 const MBinning* bins = (MBinning*)plist->FindObject("BinningTheta");
238 if (!bins)
239 {
240 *fLog << err << dbginf << "BinningTheta [MBinning] not found... aborting." << endl;
241 return kFALSE;
242 }
243
244 SetBinning(&fHEffOn, bins);
245 SetBinning(&fHChi2, bins);
246 SetBinning(&fHLambda, bins);
247 SetBinning(&fHN0del, bins);
248
249 fHEffOn.Sumw2();
250 fHChi2.Sumw2();
251 fHLambda.Sumw2();
252 fHN0del.Sumw2();
253
254 return kTRUE;
255}
256
257// -------------------------------------------------------------------------
258//
259// Draw the histogram
260//
261void MHEffOnTimeTheta::Draw(Option_t *opt)
262{
263 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
264 pad->SetBorderMode(0);
265
266 pad->Divide(2,2);
267
268 pad->cd(1);
269 gPad->SetBorderMode(0);
270 fHEffOn.Draw(opt);
271
272 pad->cd(2);
273 gPad->SetBorderMode(0);
274 fHChi2.Draw(opt);
275
276 pad->cd(3);
277 gPad->SetBorderMode(0);
278 fHLambda.Draw(opt);
279
280 pad->cd(4);
281 gPad->SetBorderMode(0);
282 fHN0del.Draw(opt);
283
284 pad->Modified();
285 pad->Update();
286}
Note: See TracBrowser for help on using the repository browser.