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

Last change on this file since 4921 was 4920, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 16.1 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),
66 fNumEvents(200*60), fNameProjDeltaT(Form("DeltaT_%p", this)),
67 fNameProjTheta(Form("Theta_%p", this))
68{
69 //
70 // set the name and title of this object
71 //
72 fName = name ? name : "MHEffectiveOnTime";
73 fTitle = title ? title : "2-D histogram in Theta and time difference";
74
75 // Main histogram
76 fH2DeltaT.SetName("DeltaT");
77 fH2DeltaT.SetXTitle("\\Delta t [s]");
78 fH2DeltaT.SetYTitle("\\Theta [\\circ]");
79 fH2DeltaT.SetZTitle("Count");
80 fH2DeltaT.UseCurrentStyle();
81 fH2DeltaT.SetDirectory(NULL);
82
83 // Main histogram
84 fH1DeltaT.SetName("DeltaT");
85 fH1DeltaT.SetXTitle("\\Delta t [s]");
86 fH1DeltaT.SetYTitle("Counts");
87 fH1DeltaT.UseCurrentStyle();
88 fH1DeltaT.SetDirectory(NULL);
89
90 // effective on time versus theta
91 fHEffOnTheta.SetName("EffOnTheta");
92 fHEffOnTheta.SetTitle("Effective On Time T_{eff}");
93 fHEffOnTheta.SetXTitle("\\Theta [\\circ]");
94 fHEffOnTheta.SetYTitle("T_{eff} [s]");
95 fHEffOnTheta.UseCurrentStyle();
96 fHEffOnTheta.SetDirectory(NULL);
97 fHEffOnTheta.GetYaxis()->SetTitleOffset(1.2);
98 //fHEffOn.Sumw2();
99
100 // effective on time versus time
101 fHEffOnTime.SetName("EffOnTime");
102 fHEffOnTime.SetTitle("Effective On Time T_{eff}");
103 fHEffOnTime.SetXTitle("Time");
104 fHEffOnTime.SetYTitle("T_{eff} [s]");
105 fHEffOnTime.UseCurrentStyle();
106 fHEffOnTime.SetDirectory(NULL);
107 fHEffOnTime.GetYaxis()->SetTitleOffset(1.2);
108 fHEffOnTime.GetXaxis()->SetLabelSize(0.033);
109 fHEffOnTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
110 fHEffOnTime.GetXaxis()->SetTimeDisplay(1);
111 fHEffOnTime.Sumw2();
112
113 // chi2 probability
114 fHProbTheta.SetName("ProbTheta");
115 fHProbTheta.SetTitle("\\chi^{2} Probability of Fit");
116 fHProbTheta.SetXTitle("\\Theta [\\circ]");
117 fHProbTheta.SetYTitle("p [%]");
118 fHProbTheta.UseCurrentStyle();
119 fHProbTheta.SetDirectory(NULL);
120 fHProbTheta.GetYaxis()->SetTitleOffset(1.2);
121 fHProbTheta.SetMaximum(101);
122
123 // chi2 probability
124 fHProbTime.SetName("ProbTime");
125 fHProbTime.SetTitle("\\chi^{2} Probability of Fit");
126 fHProbTime.SetXTitle("Time");
127 fHProbTime.SetYTitle("p [%]");
128 fHProbTime.UseCurrentStyle();
129 fHProbTime.SetDirectory(NULL);
130 fHProbTime.GetYaxis()->SetTitleOffset(1.2);
131 fHProbTime.GetXaxis()->SetLabelSize(0.033);
132 fHProbTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
133 fHProbTime.GetXaxis()->SetTimeDisplay(1);
134 fHProbTime.SetMaximum(101);
135
136 // lambda versus theta
137 fHLambda.SetName("lambda");
138 fHLambda.SetTitle("lambda of Effective On Time Fit");
139 fHLambda.SetXTitle("\\Theta [\\circ]");
140 fHLambda.SetYTitle("\\lambda [Hz]");
141 fHLambda.UseCurrentStyle();
142 fHLambda.SetDirectory(NULL);
143 // fHLambda.Sumw2();
144
145 // N0 versus theta
146 fHN0.SetName("N0");
147 fHN0.SetTitle("Ideal number of events N_{0}");
148 fHN0.SetXTitle("\\Theta [\\circ]");
149 fHN0.SetYTitle("N_{0}");
150 fHN0.UseCurrentStyle();
151 fHN0.SetDirectory(NULL);
152 //fHN0del.Sumw2();
153
154 // chi2/NDF versus theta
155 /*
156 fHChi2.SetName("Chi2/NDF");
157 fHChi2.SetTitle("\\chi^{2}/NDF of Effective On Time Fit");
158 fHChi2.SetXTitle("\\Theta [\\circ]");
159 fHChi2.SetYTitle("\\chi^{2}/NDF");
160 fHChi2.UseCurrentStyle();
161 fHChi2.SetDirectory(NULL);
162 fHChi2.GetYaxis()->SetTitleOffset(1.2);
163 //fHChi2.Sumw2();
164 */
165
166 // setup binning
167 MBinning btheta("BinningTheta");
168 btheta.SetEdgesCos(101, 0, 60);
169
170 MBinning btime("BinningDeltaT");
171 btime.SetEdges(50, 0, 0.1);
172
173 MH::SetBinning(&fH2DeltaT, &btime, &btheta);
174
175 btime.Apply(fH1DeltaT);
176
177 btheta.Apply(fHEffOnTheta);
178 btheta.Apply(fHLambda);
179 btheta.Apply(fHN0);
180 btheta.Apply(fHProbTheta);
181 //btheta.Apply(fHChi2);
182}
183
184// --------------------------------------------------------------------------
185//
186// Set the binnings and prepare the filling of the histogram
187//
188Bool_t MHEffectiveOnTime::SetupFill(const MParList *plist)
189{
190 fPointPos = (MPointingPos*)plist->FindObject("MPointingPos");
191 if (!fPointPos)
192 {
193 *fLog << err << dbginf << "MPointingPos not found... aborting." << endl;
194 return kFALSE;
195 }
196
197 // FIXME: Remove const-qualifier from base-class!
198 fTime = (MTime*)const_cast<MParList*>(plist)->FindCreateObj("MTime", "MTimeEffectiveOnTime");
199 if (!fTime)
200 return kFALSE;
201 fParam = (MParameterDerr*)const_cast<MParList*>(plist)->FindCreateObj("MParameterDerr", "MEffectiveOnTime");
202 if (!fParam)
203 return kFALSE;
204
205 const MBinning* binsdtime = (MBinning*)plist->FindObject("BinningDeltaT");
206 const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta");
207 if (binsdtime)
208 binsdtime->Apply(fH1DeltaT);
209 if (binstheta)
210 {
211 binstheta->Apply(fHEffOnTheta);
212 binstheta->Apply(fHLambda);
213 binstheta->Apply(fHN0);
214 binstheta->Apply(fHProbTheta);
215 //binstheta->Apply(fHChi2);
216 }
217 if (binstheta && binsdtime)
218 SetBinning(&fH2DeltaT, binsdtime, binstheta);
219
220 return kTRUE;
221}
222
223Bool_t MHEffectiveOnTime::FitH(TH1D *h, Double_t *res, Bool_t paint) const
224{
225 const Double_t Nm = h->Integral();
226
227 // FIXME: Do fit only if contents of bin has changed
228 if (Nm<=0)
229 return kFALSE;
230
231 // determine range (yq[0], yq[1]) of time differences
232 // where fit should be performed;
233 // require a fraction >=xq[0] of all entries to lie below yq[0]
234 // and a fraction <=xq[1] of all entries to lie below yq[1];
235 // within the range (yq[0], yq[1]) there must be no empty bin;
236 // choose pedestrian approach as long as GetQuantiles is not available
237 Double_t xq[2] = { 0.05, 0.95 };
238 Double_t yq[2];
239 h->GetQuantiles(2, yq, xq);
240
241 // Nmdel = Nm * binwidth, with Nm = number of observed events
242 const Double_t Nmdel = h->Integral("width");
243
244 //
245 // Setup Poisson function for the fit:
246 // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
247 //
248 // parameter 0 = lambda
249 // parameter 1 = N0*del
250 //
251 TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
252 func.SetParNames("lambda", "N0", "del");
253
254 func.SetParameter(0, 100); // Hz
255 func.SetParameter(1, Nm);
256 func.FixParameter(2, Nmdel/Nm);
257
258 // options : 0 do not plot the function
259 // I use integral of function in bin rather than value at bin center
260 // R use the range specified in the function range
261 // Q quiet mode
262 h->Fit(&func, "0IRQ");
263
264 const Double_t chi2 = func.GetChisquare();
265 const Int_t NDF = func.GetNDF();
266
267 // was fit successful ?
268 const Bool_t ok = NDF>0 && chi2<2.5*NDF;
269
270 if (paint)
271 {
272 func.SetLineWidth(2);
273 func.SetLineColor(ok ? kGreen : kRed);
274 func.Paint("same");
275 }
276
277 if (!ok)
278 return kFALSE;
279
280 const Double_t lambda = func.GetParameter(0);
281 const Double_t N0 = func.GetParameter(1);
282 const Double_t prob = func.GetProb();
283
284 /*
285 *fLog << all << "Nm/lambda=" << Nm/lambda << " chi2/NDF=";
286 *fLog << (NDF ? chi2/NDF : 0.0) << " lambda=";
287 *fLog << lambda << " N0=" << N0 << endl;
288 */
289
290 Double_t emat[2][2];
291 gMinuit->mnemat((Double_t*)emat, 2);
292
293 const Double_t dldl = emat[0][0];
294 //const Double_t dN0dN0 = emat[1][1];
295
296 const Double_t teff = Nm/lambda;
297 const Double_t dteff = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm);
298
299 const Double_t dl = TMath::Sqrt(dldl);
300
301 //const Double_t kappa = Nm/N0;
302 //const Double_t Rdead = 1.0 - kappa;
303 //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
304
305 // the effective on time is Nm/lambda
306 res[0] = teff;
307 res[1] = dteff;
308
309 // plot chi2-probability of fit
310 res[2] = prob*100;
311
312 // plot chi2/NDF of fit
313 //res[3] = NDF ? chi2/NDF : 0.0;
314
315 // lambda of fit
316 res[3] = lambda;
317 res[4] = dl;
318
319 // N0 of fit
320 res[5] = N0;
321
322 // Rdead (from fit) is the fraction from real time lost by the dead time
323 //fHRdead.SetBinContent(i, Rdead);
324 //fHRdead.SetBinError (i,dRdead);
325
326 return kTRUE;
327}
328
329void MHEffectiveOnTime::FitThetaBins()
330{
331 const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999));
332
333 // nbins = number of Theta bins
334 const Int_t nbins = fH2DeltaT.GetNbinsY();
335
336 TH1D *h=0;
337 for (int i=1; i<=nbins; i++)
338 {
339 // TH1D &h = *hist->ProjectionX("Calc-theta", i, i);
340 h = fH2DeltaT.ProjectionX(name, i, i, "E");
341
342 Double_t res[6];
343 if (!FitH(h, res))
344 continue;
345
346 // the effective on time is Nm/lambda
347 fHEffOnTheta.SetBinContent(i, res[0]);
348 fHEffOnTheta.SetBinError (i, res[1]);
349
350 // plot chi2-probability of fit
351 fHProbTheta.SetBinContent(i, res[2]);
352
353 // plot chi2/NDF of fit
354 //fHChi2.SetBinContent(i, res[3]);
355
356 // lambda of fit
357 fHLambda.SetBinContent(i, res[3]);
358 fHLambda.SetBinError (i, res[4]);
359
360 // N0 of fit
361 fHN0.SetBinContent(i, res[5]);
362
363 // Rdead (from fit) is the fraction from real time lost by the dead time
364 //fHRdead.SetBinContent(i, Rdead);
365 //fHRdead.SetBinError (i,dRdead);
366 }
367
368 // Histogram is reused via gROOT->FindObject()
369 // Need to be deleted only once
370 if (h)
371 delete h;
372}
373
374void MHEffectiveOnTime::FitTimeBin()
375{
376 //
377 // Fit histogram
378 //
379 Double_t res[6];
380 if (!FitH(&fH1DeltaT, res))
381 return;
382
383 // Reset Histogram
384 fH1DeltaT.Reset();
385
386 //
387 // Prepare Histogram
388 //
389
390 // Get x-axis
391 TAxis &x = *fHEffOnTime.GetXaxis();
392
393 // Get number of bins
394 const Int_t n = x.GetNbins();
395
396 // Enhance binning
397 MBinning bins;
398 bins.SetEdges(x);
399 bins.AddEdge(fLastTime.GetAxisTime());
400 bins.Apply(fHEffOnTime);
401 bins.Apply(fHProbTime);
402
403 //
404 // Fill histogram
405 //
406 fHEffOnTime.SetBinContent(n, res[0]);
407 fHEffOnTime.SetBinError(n, res[1]);
408
409 fHProbTime.SetBinContent(n, res[2]);
410
411 //
412 // Now prepare output
413 //
414 fParam->SetVal(res[0], res[1]);
415 fParam->SetReadyToSave();
416
417 *fTime = fLastTime;
418
419 // Include the current event
420 fTime->Plus1ns();
421
422 *fLog << fLastTime << ": Val=" << res[0] << " Err=" << res[1] << endl;
423}
424
425// --------------------------------------------------------------------------
426//
427// Fill the histogram
428//
429Bool_t MHEffectiveOnTime::Fill(const MParContainer *par, const Stat_t w)
430{
431 const MTime *time = dynamic_cast<const MTime*>(par);
432 if (!time)
433 {
434 *fLog << err << "ERROR - MHEffectiveOnTime::Fill without argument or container doesn't inherit from MTime... abort." << endl;
435 return kFALSE;
436 }
437
438 if (fLastTime==MTime())
439 {
440 MBinning bins;
441 bins.SetEdges(2, time->GetAxisTime()-fNumEvents/200, time->GetAxisTime());
442 bins.Apply(fHEffOnTime);
443 bins.Apply(fHProbTime);
444
445 fParam->SetVal(0, 0);
446 fParam->SetReadyToSave();
447
448 *fTime = *time;
449
450 // Make this 1ns before the first event!
451 fTime->Minus1ns();
452 }
453
454 const Double_t dt = *time-fLastTime;
455
456 fH2DeltaT.Fill(dt, fPointPos->GetZd(), w);
457 fH1DeltaT.Fill(dt, w);
458
459 fLastTime = *time;
460
461 if (fH1DeltaT.GetEntries()>=fNumEvents)
462 FitTimeBin();
463
464 return kTRUE;
465}
466
467Bool_t MHEffectiveOnTime::Finalize()
468{
469 FitThetaBins();
470 FitTimeBin();
471
472 fIsFinalized = kTRUE;
473
474 return kTRUE;
475}
476
477void MHEffectiveOnTime::PaintText(Double_t val, Double_t error) const
478{
479 TLatex text(0.45, 0.94, Form("T_{eff} = %.1fs \\pm %.1fs", val, error));
480 text.SetBit(TLatex::kTextNDC);
481 text.SetTextSize(0.04);
482 text.Paint();
483}
484
485void MHEffectiveOnTime::Paint(Option_t *opt)
486{
487 TH1D *h=0;
488
489 TString o(opt);
490 if (o==(TString)"fit")
491 {
492 TVirtualPad *pad = gPad;
493
494 pad->GetPad(1)->cd(1);
495 if ((h = (TH1D*)gPad->FindObject(fNameProjDeltaT)))
496 {
497 h = fH2DeltaT.ProjectionX(fNameProjDeltaT, -1, 9999, "E");
498 if (h->GetEntries()>0)
499 gPad->SetLogy();
500 }
501
502 pad->GetPad(2)->cd(1);
503 if ((h = (TH1D*)gPad->FindObject(fNameProjTheta)))
504 fH2DeltaT.ProjectionY(fNameProjTheta, -1, 9999, "E");
505
506 if (!fIsFinalized)
507 FitThetaBins();
508 return;
509 }
510 if (o==(TString)"paint")
511 {
512 if ((h = (TH1D*)gPad->FindObject(fNameProjDeltaT)))
513 {
514 Double_t res[6];
515 if (FitH(h, res, kTRUE))
516 PaintText(res[0], res[1]);
517 }
518 return;
519 }
520
521 h=0;
522 if (o==(TString)"theta")
523 h = &fHEffOnTheta;
524 if (o==(TString)"time")
525 h = &fHEffOnTime;
526
527 if (!h)
528 return;
529
530 Double_t error = 0;
531 for (int i=0; i<h->GetXaxis()->GetNbins(); i++)
532 error += h->GetBinError(i);
533
534 PaintText(h->Integral(), error);
535}
536
537// --------------------------------------------------------------------------
538//
539// Draw the histogram
540//
541void MHEffectiveOnTime::Draw(Option_t *opt)
542{
543 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
544 pad->SetBorderMode(0);
545
546 AppendPad("fit");
547
548 pad->Divide(2, 1, 0, 0);
549
550 TH1 *h;
551
552 pad->cd(1);
553 gPad->SetBorderMode(0);
554 gPad->Divide(1, 3, 0, 0);
555 pad->GetPad(1)->cd(1);
556 gPad->SetBorderMode(0);
557 h = fH2DeltaT.ProjectionX(fNameProjDeltaT, -1, 9999, "E");
558 h->SetTitle("Distribution of \\Delta t [s]");
559 h->SetXTitle("\\Delta t [s]");
560 h->SetYTitle("Counts");
561 h->SetDirectory(NULL);
562 h->SetMarkerStyle(kFullDotMedium);
563 h->SetBit(kCanDelete);
564 h->Draw();
565 AppendPad("paint");
566
567 pad->GetPad(1)->cd(2);
568 gPad->SetBorderMode(0);
569 fHProbTime.Draw();
570
571 pad->GetPad(1)->cd(3);
572 gPad->SetBorderMode(0);
573 fHEffOnTime.Draw();
574 AppendPad("time");
575
576 pad->cd(2);
577 gPad->SetBorderMode(0);
578 gPad->Divide(1, 3, 0, 0);
579
580 pad->GetPad(2)->cd(1);
581 gPad->SetBorderMode(0);
582 h = fH2DeltaT.ProjectionY(fNameProjTheta, -1, 9999, "E");
583 h->SetTitle("Distribution of \\Theta [\\circ]");
584 h->SetXTitle("\\Theta [\\circ]");
585 h->SetYTitle("Counts");
586 h->SetDirectory(NULL);
587 h->SetMarkerStyle(kFullDotMedium);
588 h->SetBit(kCanDelete);
589 h->GetYaxis()->SetTitleOffset(1.1);
590 h->Draw();
591
592 pad->GetPad(2)->cd(2);
593 gPad->SetBorderMode(0);
594 fHProbTheta.Draw();
595
596 pad->GetPad(2)->cd(3);
597 gPad->SetBorderMode(0);
598 fHEffOnTheta.Draw();
599 AppendPad("theta");
600}
Note: See TracBrowser for help on using the repository browser.