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

Last change on this file since 4905 was 4889, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 13.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, 8/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-2004
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28// MHEffectiveOnTime
29//
30// Fills a 2-D histogram Time-difference vs. Theta
31//
32// From this histogram the effective on-time is determined by a fit.
33// The result of the fit (see Fit()) and the fit-parameters (like chi^2)
34// are stored in corresponding histograms
35//
36//////////////////////////////////////////////////////////////////////////////
37
38#include "MHEffectiveOnTime.h"
39
40#include <TF1.h>
41#include <TLatex.h>
42#include <TCanvas.h>
43#include <TMinuit.h>
44#include <TRandom.h>
45
46#include "MTime.h"
47#include "MParameters.h"
48#include "MPointingPos.h"
49
50#include "MBinning.h"
51#include "MParList.h"
52
53#include "MLog.h"
54#include "MLogManip.h"
55
56ClassImp(MHEffectiveOnTime);
57
58using namespace std;
59
60// --------------------------------------------------------------------------
61//
62// Default Constructor. It sets name and title of the histogram.
63//
64MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title)
65 : fPointPos(0), fTime(0), fParam(0), fIsFinalized(kFALSE), fInterval(60)
66{
67 //
68 // set the name and title of this object
69 //
70 fName = name ? name : "MHEffectiveOnTime";
71 fTitle = title ? title : "2-D histogram in Theta and time difference";
72
73 // Main histogram
74 fHTimeDiff.SetXTitle("\\Delta t [s]");
75 fHTimeDiff.SetYTitle("\\Theta [\\circ]");
76 fHTimeDiff.UseCurrentStyle();
77 fHTimeDiff.SetDirectory(NULL);
78
79 // effective on time versus theta
80 fHEffOn.SetName("EffOnTime");
81 fHEffOn.SetTitle("Effective On Time T_{eff}");
82 fHEffOn.SetXTitle("\\Theta [\\circ]");
83 fHEffOn.SetYTitle("T_{eff} [s]");
84 fHEffOn.UseCurrentStyle();
85 fHEffOn.SetDirectory(NULL);
86 fHEffOn.GetYaxis()->SetTitleOffset(1.2);
87 //fHEffOn.Sumw2();
88
89 // chi2/NDF versus theta
90 fHChi2.SetName("Chi2/NDF");
91 fHChi2.SetTitle("\\chi^{2}/NDF of Effective On Time Fit");
92 fHChi2.SetXTitle("\\Theta [\\circ]");
93 fHChi2.SetYTitle("\\chi^{2}/NDF");
94 fHChi2.UseCurrentStyle();
95 fHChi2.SetDirectory(NULL);
96 fHChi2.GetYaxis()->SetTitleOffset(1.2);
97 //fHChi2.Sumw2();
98
99 // chi2 probability
100 fHProb.SetName("Prob");
101 fHProb.SetTitle("\\chi^{2} Probability of Fit");
102 fHProb.SetXTitle("\\Theta [\\circ]");
103 fHProb.SetYTitle("p [%]");
104 fHProb.UseCurrentStyle();
105 fHProb.SetDirectory(NULL);
106 fHProb.GetYaxis()->SetTitleOffset(1.2);
107 fHProb.SetMaximum(101);
108 //fHChi2.Sumw2();
109
110 // lambda versus theta
111 fHLambda.SetName("lambda");
112 fHLambda.SetTitle("lambda of Effective On Time Fit");
113 fHLambda.SetXTitle("\\Theta [\\circ]");
114 fHLambda.SetYTitle("\\lambda [Hz]");
115 fHLambda.UseCurrentStyle();
116 fHLambda.SetDirectory(NULL);
117 // fHLambda.Sumw2();
118
119 // N0 versus theta
120 fHN0.SetName("N0");
121 fHN0.SetTitle("Ideal number of events N_{0}");
122 fHN0.SetXTitle("\\Theta [\\circ]");
123 fHN0.SetYTitle("N_{0}");
124 fHN0.UseCurrentStyle();
125 fHN0.SetDirectory(NULL);
126 //fHN0del.Sumw2();
127
128 // setup binning
129 MBinning btheta("BinningTheta");
130 btheta.SetEdgesCos(51, 0, 60);
131
132 MBinning btime("BinningTimeDiff");
133 btime.SetEdges(50, 0, 0.1);
134
135 MH::SetBinning(&fHTimeDiff, &btime, &btheta);
136
137 btheta.Apply(fHEffOn);
138 btheta.Apply(fHChi2);
139 btheta.Apply(fHLambda);
140 btheta.Apply(fHN0);
141 btheta.Apply(fHProb);
142}
143
144// FIXME: Just for a preliminary check
145static Double_t testval = 0;
146static Double_t testerr = 0;
147
148// --------------------------------------------------------------------------
149//
150// Set the binnings and prepare the filling of the histogram
151//
152Bool_t MHEffectiveOnTime::SetupFill(const MParList *plist)
153{
154 fPointPos = (MPointingPos*)plist->FindObject("MPointingPos");
155 if (!fPointPos)
156 {
157 *fLog << err << dbginf << "MPointingPos not found... aborting." << endl;
158 return kFALSE;
159 }
160
161 // FIXME: Remove const-qualifier from base-class!
162 fTime = (MTime*)const_cast<MParList*>(plist)->FindCreateObj("MTime", "MTimeEffectiveOnTime");
163 if (!fTime)
164 return kFALSE;
165 fParam = (MParameterDerr*)const_cast<MParList*>(plist)->FindCreateObj("MParameterDerr", "MEffectiveOnTime");
166 if (!fParam)
167 return kFALSE;
168
169 const MBinning* binsdtime = (MBinning*)plist->FindObject("BinningTimeDiff");
170 const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta");
171 if (!binstheta || !binsdtime)
172 *fLog << warn << dbginf << "At least one MBinning not found... ignored." << endl;
173 else
174 {
175 SetBinning(&fHTimeDiff, binsdtime, binstheta);
176
177 binstheta->Apply(fHEffOn);
178 binstheta->Apply(fHChi2);
179 binstheta->Apply(fHLambda);
180 binstheta->Apply(fHN0);
181 binstheta->Apply(fHProb);
182 }
183
184 return kTRUE;
185}
186
187Bool_t MHEffectiveOnTime::Finalize()
188{
189 FitThetaBins();
190 Calc();
191
192 fIsFinalized = kTRUE;
193
194 return kTRUE;
195}
196
197void MHEffectiveOnTime::FitThetaBins()
198{
199 const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999));
200
201 // nbins = number of Theta bins
202 const Int_t nbins = fHTimeDiff.GetNbinsY();
203
204 TH1D *h=0;
205 for (int i=1; i<=nbins; i++)
206 {
207 // TH1D &h = *hist->ProjectionX("Calc-theta", i, i);
208 h = fHTimeDiff.ProjectionX(name, i, i, "E");
209
210 Double_t res[7];
211 if (!FitH(h, res))
212 continue;
213
214 // the effective on time is Nm/lambda
215 fHEffOn.SetBinContent(i, res[0]);
216 fHEffOn.SetBinError (i, res[1]);
217
218 // plot chi2-probability of fit
219 fHProb.SetBinContent(i, res[2]);
220
221 // plot chi2/NDF of fit
222 fHChi2.SetBinContent(i, res[3]);
223
224 // lambda of fit
225 fHLambda.SetBinContent(i, res[4]);
226 fHLambda.SetBinError (i, res[5]);
227
228 // N0 of fit
229 fHN0.SetBinContent(i, res[6]);
230
231 // Rdead (from fit) is the fraction from real time lost by the dead time
232 //fHRdead.SetBinContent(i, Rdead);
233 //fHRdead.SetBinError (i,dRdead);
234 }
235
236 // Histogram is reused via gROOT->FindObject()
237 // Need to be deleted only once
238 if (h)
239 delete h;
240}
241
242Bool_t MHEffectiveOnTime::FitH(TH1D *h, Double_t *res) const
243{
244 const Double_t Nm = h->Integral();
245
246 // FIXME: Do fit only if contents of bin has changed
247 if (Nm<=0)
248 return kFALSE;
249
250 // determine range (yq[0], yq[1]) of time differences
251 // where fit should be performed;
252 // require a fraction >=xq[0] of all entries to lie below yq[0]
253 // and a fraction <=xq[1] of all entries to lie below yq[1];
254 // within the range (yq[0], yq[1]) there must be no empty bin;
255 // choose pedestrian approach as long as GetQuantiles is not available
256 Double_t xq[2] = { 0.05, 0.95 };
257 Double_t yq[2];
258 h->GetQuantiles(2, yq, xq);
259
260 // Nmdel = Nm * binwidth, with Nm = number of observed events
261 const Double_t Nmdel = h->Integral("width");
262
263 //
264 // Setup Poisson function for the fit:
265 // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
266 //
267 // parameter 0 = lambda
268 // parameter 1 = N0*del
269 //
270 TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
271 func.SetParNames("lambda", "N0", "del");
272
273 func.SetParameter(0, 100); // Hz
274 func.SetParameter(1, Nm);
275 func.FixParameter(2, Nmdel/Nm);
276
277 // options : 0 do not plot the function
278 // I use integral of function in bin rather than value at bin center
279 // R use the range specified in the function range
280 // Q quiet mode
281 h->Fit(&func, "0IRQ");
282
283 const Double_t chi2 = func.GetChisquare();
284 const Int_t NDF = func.GetNDF();
285
286 // was fit successful ?
287 if (NDF<=0 || chi2>=2.5*NDF)
288 return kFALSE;
289
290 const Double_t lambda = func.GetParameter(0);
291 const Double_t N0 = func.GetParameter(1);
292 const Double_t prob = func.GetProb();
293
294 /*
295 *fLog << all << "Nm/lambda=" << Nm/lambda << " chi2/NDF=";
296 *fLog << (NDF ? chi2/NDF : 0.0) << " lambda=";
297 *fLog << lambda << " N0=" << N0 << endl;
298 */
299
300 Double_t emat[2][2];
301 gMinuit->mnemat((Double_t*)emat, 2);
302
303 const Double_t dldl = emat[0][0];
304 //const Double_t dN0dN0 = emat[1][1];
305
306 const Double_t teff = Nm/lambda;
307 const Double_t dteff = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm);
308
309 const Double_t dl = TMath::Sqrt(dldl);
310
311 //const Double_t kappa = Nm/N0;
312 //const Double_t Rdead = 1.0 - kappa;
313 //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
314
315 // the effective on time is Nm/lambda
316 res[0] = teff;
317 res[1] = dteff;
318
319 // plot chi2-probability of fit
320 res[2] = prob*100;
321
322 // plot chi2/NDF of fit
323 res[3] = NDF ? chi2/NDF : 0.0;
324
325 // lambda of fit
326 res[4] = lambda;
327 res[5] = dl;
328
329 // N0 of fit
330 res[6] = N0;
331
332 // Rdead (from fit) is the fraction from real time lost by the dead time
333 //fHRdead.SetBinContent(i, Rdead);
334 //fHRdead.SetBinError (i,dRdead);
335
336 return kTRUE;
337}
338
339void MHEffectiveOnTime::Paint(Option_t *opt)
340{
341 TVirtualPad *padsave = gPad;
342
343 TString o(opt);
344 if (o==(TString)"fit")
345 {
346 TH1D *h0=0;
347
348 padsave->cd(1);
349 if ((h0 = (TH1D*)gPad->FindObject("TimeDiff")))
350 {
351 h0 = fHTimeDiff.ProjectionX("TimeDiff", -1, 9999, "E");
352 if (h0->GetEntries()>0)
353 gPad->SetLogy();
354 }
355
356 padsave->cd(2);
357 if ((h0 = (TH1D*)gPad->FindObject("Theta")))
358 fHTimeDiff.ProjectionY("Theta", -1, 9999, "E");
359
360 if (!fIsFinalized)
361 FitThetaBins();
362 }
363 if (o==(TString)"paint")
364 {
365 Double_t error = 0;
366 for (int i=0; i<fHEffOn.GetXaxis()->GetNbins(); i++)
367 error += fHEffOn.GetBinError(i);
368
369 TLatex text(0.45, 0.94, Form("T_{eff} = %.1fs \\pm %.1fs", fHEffOn.Integral(), error));
370 text.SetBit(TLatex::kTextNDC);
371 text.SetTextSize(0.04);
372 text.Paint();
373 }
374 gPad = padsave;
375}
376
377// --------------------------------------------------------------------------
378//
379// Draw the histogram
380//
381void MHEffectiveOnTime::Draw(Option_t *opt)
382{
383 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
384 pad->SetBorderMode(0);
385
386 AppendPad("fit");
387
388 pad->Divide(2,2);
389
390 TH1 *h;
391
392 pad->cd(1);
393 gPad->SetBorderMode(0);
394 h = fHTimeDiff.ProjectionX("TimeDiff", -1, 9999, "E");
395 h->SetTitle("Distribution of \\Delta t [s]");
396 h->SetXTitle("\\Delta t [s]");
397 h->SetYTitle("Counts");
398 h->SetDirectory(NULL);
399 h->SetMarkerStyle(kFullDotMedium);
400 h->SetBit(kCanDelete);
401 h->Draw();
402
403 pad->cd(2);
404 gPad->SetBorderMode(0);
405 h = fHTimeDiff.ProjectionY("Theta", -1, 9999, "E");
406 h->SetTitle("Distribution of \\Theta [\\circ]");
407 h->SetXTitle("\\Theta [\\circ]");
408 h->SetYTitle("Counts");
409 h->SetDirectory(NULL);
410 h->SetMarkerStyle(kFullDotMedium);
411 h->SetBit(kCanDelete);
412 h->GetYaxis()->SetTitleOffset(1.1);
413 h->Draw();
414
415 pad->cd(3);
416 gPad->SetBorderMode(0);
417 fHEffOn.Draw();
418 AppendPad("paint");
419
420 pad->cd(4);
421 gPad->SetBorderMode(0);
422 fHProb.Draw();
423}
424
425void MHEffectiveOnTime::Calc()
426{
427 TH1D *h = fHTimeDiff.ProjectionX("", -1, 99999, "E");
428 h->SetDirectory(0);
429
430 Double_t res[7];
431 Bool_t rc = FitH(h, res);
432 delete h;
433
434 if (!rc)
435 return;
436
437 const Double_t val = res[0];
438 const Double_t error = res[1];
439
440 fParam->SetVal(val-fEffOnTime0, error-fEffOnErr0);
441 fParam->SetReadyToSave();
442
443 testval += fParam->GetVal();
444 testerr += fParam->GetErr();
445
446 fEffOnTime0 = val;
447 fEffOnErr0 = error;
448
449 MTime now(*fTime);
450 now.AddMilliSeconds(fInterval*1000);
451
452 *fLog <<all << now << " - ";// << fLastTime-fTime;
453 *fLog << Form("T_{eff} = %.1fs \\pm %.1fs",
454 fParam->GetVal(), fParam->GetErr());
455 *fLog << Form(" %.1f %.1f %.1f %.1f",
456 val, testval, error, testerr) << endl;
457
458 fTime->AddMilliSeconds(fInterval*1000);
459}
460
461// --------------------------------------------------------------------------
462//
463// Fill the histogram
464//
465Bool_t MHEffectiveOnTime::Fill(const MParContainer *par, const Stat_t w)
466{
467 const MTime *time = dynamic_cast<const MTime*>(par);
468 if (!time)
469 {
470 *fLog << err << "ERROR - MHEffectiveOnTimeTime::Fill without argument or container doesn't inherit from MTime... abort." << endl;
471 return kFALSE;
472 }
473
474 if (*fTime==MTime())
475 {
476 *fLog << all << *time << " - " << *fTime << " " << TMath::Floor(*time) << endl;
477 fEffOnTime0 = 0;
478 fEffOnErr0 = 0;
479
480 fParam->SetVal(0, 0);
481 fParam->SetReadyToSave();
482
483 *fTime = *time;
484
485 // Make this just a ns before the first event
486 fTime->Minus1ns();
487 }
488
489 if (*fTime+fInterval<*time)
490 {
491 FitThetaBins();
492 Calc();
493 }
494
495 fHTimeDiff.Fill(*time-fLastTime, fPointPos->GetZd(), w);
496 fLastTime = *time;
497
498 return kTRUE;
499}
500
Note: See TracBrowser for help on using the repository browser.