source: trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc@ 5138

Last change on this file since 5138 was 5012, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 23.4 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// Filling this you will get the effective on-time versus theta and
31// observation time.
32//
33// From this histogram the effective on-time is determined by a fit.
34// The result of the fit (see Fit()) and the fit-parameters (like chi^2)
35// are stored in corresponding histograms
36//
37// To determin the efective on time a poisson fit is done. For more details
38// please have a look into the source code of FitH() it should be simple
39// to understand. In this function a Delta-T distribution is fitted, while
40// Delta-T is the time between two consecutive events.
41//
42// The fit is done for projections of a 2D histogram in Theta and Delta T.
43// So you get the effective on time versus theta.
44//
45// To get the effective on-time versus time a histogram is filled with
46// the Delta-T distribution of a number of events set by SetNumEvents().
47// The default is 12000 (roughly 1min at 200Hz)
48//
49// For each "time-bin" the histogram is fitted and the resulting effective
50// on-time is stored in the fHTimeEffOn histogram. Each entry in this
51// histogram is the effective observation time between the upper and
52// lower edges of the bins.
53//
54// In addition the calculated effective on time is stored in a
55// "MEffectiveOnTime [MParameterDerr]" and the corresponding time-stamp
56// (the upper edge of the bin) "MTimeEffectiveOnTime [MTime]"
57//
58// The class takes two binnings from the Parameter list; if these binnings
59// are not available the defaultbinning is used:
60// MBinning("BinningDeltaT"); // Units of seconds
61// MBinning("BinningTheta"); // Units of degrees
62//
63//
64// Usage:
65// ------
66// MFillH fill("MHEffectiveOnTime", "MTime");
67// tlist.AddToList(&fill);
68//
69//
70// Input Container:
71// MPointingPos
72//
73// Output Container:
74// MEffectiveOnTime [MParameterDerr]
75// MTimeEffectiveOnTime [MTime]
76//
77//////////////////////////////////////////////////////////////////////////////
78#include "MHEffectiveOnTime.h"
79
80#include <TF1.h>
81#include <TMinuit.h>
82#include <TRandom.h>
83
84#include <TLatex.h>
85#include <TGaxis.h>
86#include <TCanvas.h>
87#include <TPaveStats.h>
88
89#include "MTime.h"
90#include "MParameters.h"
91#include "MPointingPos.h"
92
93#include "MBinning.h"
94#include "MParList.h"
95
96#include "MLog.h"
97#include "MLogManip.h"
98
99ClassImp(MHEffectiveOnTime);
100
101using namespace std;
102
103// --------------------------------------------------------------------------
104//
105// Default Constructor. It initializes all histograms.
106//
107MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title)
108 : fPointPos(0), fTime(0), fParam(0), fIsFinalized(kFALSE),
109 fNumEvents(200*60), fNameProjDeltaT(Form("DeltaT_%p", this)),
110 fNameProjTheta(Form("Theta_%p", this))
111{
112 //
113 // set the name and title of this object
114 //
115 fName = name ? name : "MHEffectiveOnTime";
116 fTitle = title ? title : "Histogram to determin effective On-Time vs Time and Zenith Angle";
117
118 // Main histogram
119 fH2DeltaT.SetName("DeltaT");
120 fH2DeltaT.SetXTitle("\\Delta t [s]");
121 fH2DeltaT.SetYTitle("\\Theta [\\circ]");
122 fH2DeltaT.SetZTitle("Count");
123 fH2DeltaT.UseCurrentStyle();
124 fH2DeltaT.SetDirectory(NULL);
125
126 // Main histogram
127 fH1DeltaT.SetName("DeltaT");
128 fH1DeltaT.SetXTitle("\\Delta t [s]");
129 fH1DeltaT.SetYTitle("Counts");
130 fH1DeltaT.UseCurrentStyle();
131 fH1DeltaT.SetDirectory(NULL);
132
133 // effective on time versus theta
134 fHThetaEffOn.SetName("EffOnTheta");
135 fHThetaEffOn.SetTitle("Effective On Time T_{eff}");
136 fHThetaEffOn.SetXTitle("\\Theta [\\circ]");
137 fHThetaEffOn.SetYTitle("T_{eff} [s]");
138 fHThetaEffOn.UseCurrentStyle();
139 fHThetaEffOn.SetDirectory(NULL);
140 fHThetaEffOn.GetYaxis()->SetTitleOffset(1.2);
141 fHThetaEffOn.GetYaxis()->SetTitleColor(kBlue);
142 fHThetaEffOn.SetLineColor(kBlue);
143 //fHEffOn.Sumw2();
144
145 // effective on time versus time
146 fHTimeEffOn.SetName("EffOnTime");
147 fHTimeEffOn.SetTitle("Effective On Time T_{eff}");
148 fHTimeEffOn.SetXTitle("Time");
149 fHTimeEffOn.SetYTitle("T_{eff} [s]");
150 fHTimeEffOn.UseCurrentStyle();
151 fHTimeEffOn.SetDirectory(NULL);
152 fHTimeEffOn.GetYaxis()->SetTitleOffset(1.2);
153 fHTimeEffOn.GetXaxis()->SetLabelSize(0.033);
154 fHTimeEffOn.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
155 fHTimeEffOn.GetXaxis()->SetTimeDisplay(1);
156 fHTimeEffOn.GetYaxis()->SetTitleColor(kBlue);
157 fHTimeEffOn.SetLineColor(kBlue);
158
159 // chi2 probability
160 fHThetaProb.SetName("ProbTheta");
161 fHThetaProb.SetTitle("\\chi^{2} Probability of Fit");
162 fHThetaProb.SetXTitle("\\Theta [\\circ]");
163 fHThetaProb.SetYTitle("p [%]");
164 fHThetaProb.UseCurrentStyle();
165 fHThetaProb.SetDirectory(NULL);
166 fHThetaProb.GetYaxis()->SetTitleOffset(1.2);
167 fHThetaProb.SetMaximum(101);
168 fHThetaProb.GetYaxis()->SetTitleColor(kBlue);
169 fHThetaProb.SetLineColor(kBlue);
170
171 // chi2 probability
172 fHTimeProb.SetName("ProbTime");
173 fHTimeProb.SetTitle("\\chi^{2} Probability of Fit");
174 fHTimeProb.SetXTitle("Time");
175 fHTimeProb.SetYTitle("p [%]");
176 fHTimeProb.UseCurrentStyle();
177 fHTimeProb.SetDirectory(NULL);
178 fHTimeProb.GetYaxis()->SetTitleOffset(1.2);
179 fHTimeProb.GetXaxis()->SetLabelSize(0.033);
180 fHTimeProb.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
181 fHTimeProb.GetXaxis()->SetTimeDisplay(1);
182 fHTimeProb.SetMaximum(101);
183 fHTimeProb.GetYaxis()->SetTitleColor(kBlue);
184 fHTimeProb.SetLineColor(kBlue);
185
186 // lambda versus theta
187 fHThetaLambda.SetName("LambdaTheta");
188 fHThetaLambda.SetTitle("Slope (Rate) vs Theta");
189 fHThetaLambda.SetXTitle("\\Theta [\\circ]");
190 fHThetaLambda.SetYTitle("\\lambda [s^{-1}]");
191 fHThetaLambda.UseCurrentStyle();
192 fHThetaLambda.SetDirectory(NULL);
193 fHThetaLambda.SetLineColor(kGreen);
194
195 // lambda versus time
196 fHTimeLambda.SetName("LambdaTime");
197 fHTimeLambda.SetTitle("Slope (Rate) vs Time");
198 fHTimeLambda.SetXTitle("\\Time [\\circ]");
199 fHTimeLambda.SetYTitle("\\lambda [s^{-1}]");
200 fHTimeLambda.UseCurrentStyle();
201 fHTimeLambda.SetDirectory(NULL);
202 fHTimeLambda.GetYaxis()->SetTitleOffset(1.2);
203 fHTimeLambda.GetXaxis()->SetLabelSize(0.033);
204 fHTimeLambda.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
205 fHTimeLambda.GetXaxis()->SetTimeDisplay(1);
206 fHTimeLambda.SetLineColor(kGreen);
207
208 // NDoF versus theta
209 fHThetaNDF.SetName("NDofTheta");
210 fHThetaNDF.SetTitle("Number of Degrees of freedom vs Theta");
211 fHThetaNDF.SetXTitle("\\Theta [\\circ]");
212 fHThetaNDF.SetYTitle("NDoF [#]");
213 fHThetaNDF.UseCurrentStyle();
214 fHThetaNDF.SetDirectory(NULL);
215 fHThetaNDF.SetLineColor(kGreen);
216
217 // NDoF versus time
218 /*
219 fHTimeNDF.SetName("NDofTime");
220 fHTimeNDF.SetTitle("Number of Degrees of freedom vs Time");
221 fHTimeNDF.SetXTitle("Time");
222 fHTimeNDF.SetYTitle("NDoF [#]");
223 fHTimeNDF.UseCurrentStyle();
224 fHTimeNDF.SetDirectory(NULL);
225 fHTimeNDF.GetYaxis()->SetTitleOffset(1.2);
226 fHTimeNDF.GetXaxis()->SetLabelSize(0.033);
227 fHTimeNDF.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
228 fHTimeNDF.GetXaxis()->SetTimeDisplay(1);
229 fHTimeNDF.SetLineColor(kBlue);
230 */
231 // setup binning
232 MBinning btheta("BinningTheta");
233 btheta.SetEdgesCos(100, 0, 60);
234
235 MBinning btime("BinningDeltaT");
236 btime.SetEdges(50, 0, 0.1);
237
238 MH::SetBinning(&fH2DeltaT, &btime, &btheta);
239
240 btime.Apply(fH1DeltaT);
241
242 btheta.Apply(fHThetaEffOn);
243 btheta.Apply(fHThetaLambda);
244 btheta.Apply(fHThetaNDF);
245 btheta.Apply(fHThetaProb);
246 //btheta.Apply(fHChi2);
247}
248
249// --------------------------------------------------------------------------
250//
251// Set the binnings and prepare the filling of the histogram
252//
253Bool_t MHEffectiveOnTime::SetupFill(const MParList *plist)
254{
255 fPointPos = (MPointingPos*)plist->FindObject("MPointingPos");
256 if (!fPointPos)
257 {
258 *fLog << err << dbginf << "MPointingPos not found... aborting." << endl;
259 return kFALSE;
260 }
261
262 // FIXME: Remove const-qualifier from base-class!
263 fTime = (MTime*)const_cast<MParList*>(plist)->FindCreateObj("MTime", "MTimeEffectiveOnTime");
264 if (!fTime)
265 return kFALSE;
266 fParam = (MParameterDerr*)const_cast<MParList*>(plist)->FindCreateObj("MParameterDerr", "MEffectiveOnTime");
267 if (!fParam)
268 return kFALSE;
269
270 const MBinning* binsdtime = (MBinning*)plist->FindObject("BinningDeltaT");
271 const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta");
272 if (binsdtime)
273 binsdtime->Apply(fH1DeltaT);
274 if (binstheta)
275 {
276 binstheta->Apply(fHThetaEffOn);
277 binstheta->Apply(fHThetaLambda);
278 binstheta->Apply(fHThetaNDF);
279 binstheta->Apply(fHThetaProb);
280 //binstheta->Apply(fHChi2);
281 }
282 if (binstheta && binsdtime)
283 SetBinning(&fH2DeltaT, binsdtime, binstheta);
284
285 return kTRUE;
286}
287
288// --------------------------------------------------------------------------
289//
290// Fit a single Delta-T distribution. See source code for more details
291//
292Bool_t MHEffectiveOnTime::FitH(TH1D *h, Double_t *res, Bool_t paint) const
293{
294 const Double_t Nm = h->Integral();
295
296 // FIXME: Do fit only if contents of bin has changed
297 if (Nm<=0)
298 return kFALSE;
299
300 // determine range (yq[0], yq[1]) of time differences
301 // where fit should be performed;
302 // require a fraction >=xq[0] of all entries to lie below yq[0]
303 // and a fraction <=xq[1] of all entries to lie below yq[1];
304 // within the range (yq[0], yq[1]) there must be no empty bin;
305 // choose pedestrian approach as long as GetQuantiles is not available
306 Double_t xq[2] = { 0.05, 0.95 };
307 Double_t yq[2];
308 h->GetQuantiles(2, yq, xq);
309
310 // Nmdel = Nm * binwidth, with Nm = number of observed events
311 const Double_t Nmdel = h->Integral("width");
312
313 //
314 // Setup Poisson function for the fit:
315 // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
316 //
317 // parameter 0 = lambda
318 // parameter 1 = N0*del
319 //
320 TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)");
321 //func.SetParNames("lambda", "N0", "del");
322
323 func.SetParameter(0, 100); // Hz
324 func.SetParameter(1, Nm);
325 func.FixParameter(2, Nmdel/Nm);
326
327 // options : N do not store the function, do not draw
328 // I use integral of function in bin rather than value at bin center
329 // R use the range specified in the function range
330 // Q quiet mode
331 h->Fit(&func, "NIQ", "", yq[0], yq[1]);
332
333 const Double_t chi2 = func.GetChisquare();
334 const Int_t NDF = func.GetNDF();
335
336 // was fit successful ?
337 const Bool_t ok = NDF>0 && chi2<2.5*NDF;
338
339 if (paint)
340 {
341 func.SetLineWidth(2);
342 func.SetLineColor(ok ? kGreen : kRed);
343 func.Paint("same");
344 }
345
346 if (!ok)
347 return kFALSE;
348
349 const Double_t lambda = func.GetParameter(0);
350 //const Double_t N0 = func.GetParameter(1);
351 const Double_t prob = func.GetProb();
352
353 /*
354 *fLog << all << "Nm/lambda=" << Nm/lambda << " chi2/NDF=";
355 *fLog << (NDF ? chi2/NDF : 0.0) << " lambda=";
356 *fLog << lambda << " N0=" << N0 << endl;
357 */
358
359 Double_t emat[2][2];
360 gMinuit->mnemat((Double_t*)emat, 2);
361
362 const Double_t dldl = emat[0][0];
363 //const Double_t dN0dN0 = emat[1][1];
364
365 const Double_t teff = Nm/lambda;
366 const Double_t dteff = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm);
367 const Double_t dl = TMath::Sqrt(dldl);
368
369 //const Double_t kappa = Nm/N0;
370 //const Double_t Rdead = 1.0 - kappa;
371 //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
372
373 // the effective on time is Nm/lambda
374 res[0] = teff;
375 res[1] = dteff;
376
377 // plot chi2-probability of fit
378 res[2] = prob*100;
379
380 // lambda of fit
381 res[3] = lambda;
382 res[4] = dl;
383
384 // NDoF of fit
385 res[5] = NDF;
386
387 // Rdead (from fit) is the fraction from real time lost by the dead time
388 //fHRdead.SetBinContent(i, Rdead);
389 //fHRdead.SetBinError (i,dRdead);
390
391 return kTRUE;
392}
393
394// --------------------------------------------------------------------------
395//
396// Fit a all bins of the distribution in theta. Store the result in the
397// Theta-Histograms
398//
399void MHEffectiveOnTime::FitThetaBins()
400{
401 fHThetaEffOn.Reset();
402 fHThetaProb.Reset();
403 fHThetaLambda.Reset();
404 fHThetaNDF.Reset();
405
406 const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999));
407
408 // nbins = number of Theta bins
409 const Int_t nbins = fH2DeltaT.GetNbinsY();
410
411 TH1D *h=0;
412 for (int i=1; i<=nbins; i++)
413 {
414 // TH1D &h = *hist->ProjectionX("Calc-theta", i, i);
415 h = fH2DeltaT.ProjectionX(name, i, i, "E");
416
417 Double_t res[6];
418 if (!FitH(h, res))
419 continue;
420
421 // the effective on time is Nm/lambda
422 fHThetaEffOn.SetBinContent(i, res[0]);
423 fHThetaEffOn.SetBinError (i, res[1]);
424
425 // plot chi2-probability of fit
426 fHThetaProb.SetBinContent(i, res[2]);
427
428 // plot chi2/NDF of fit
429 //fHChi2.SetBinContent(i, res[3]);
430
431 // lambda of fit
432 fHThetaLambda.SetBinContent(i, res[3]);
433 fHThetaLambda.SetBinError (i, res[4]);
434
435 // NDoF of fit
436 fHThetaNDF.SetBinContent(i, res[5]);
437
438 // Rdead (from fit) is the fraction from real time lost by the dead time
439 //fHRdead.SetBinContent(i, Rdead);
440 //fHRdead.SetBinError (i,dRdead);
441 }
442
443 // Histogram is reused via gROOT->FindObject()
444 // Need to be deleted only once
445 if (h)
446 delete h;
447}
448
449// --------------------------------------------------------------------------
450//
451// Fit the single-time-bin histogram. Store the result in the
452// Time-Histograms
453//
454void MHEffectiveOnTime::FitTimeBin()
455{
456 //
457 // Fit histogram
458 //
459 Double_t res[6];
460 if (!FitH(&fH1DeltaT, res))
461 return;
462
463 // Reset Histogram
464 fH1DeltaT.Reset();
465
466 //
467 // Prepare Histogram
468 //
469
470 // Get number of bins
471 const Int_t n = fHTimeEffOn.GetNbinsX();
472
473 // Enhance binning
474 MBinning bins;
475 bins.SetEdges(fHTimeEffOn, 'x');
476 bins.AddEdge(fLastTime.GetAxisTime());
477 bins.Apply(fHTimeEffOn);
478 bins.Apply(fHTimeProb);
479 bins.Apply(fHTimeLambda);
480 //bins.Apply(fHTimeNDF);
481
482 //
483 // Fill histogram
484 //
485 fHTimeEffOn.SetBinContent(n, res[0]);
486 fHTimeEffOn.SetBinError(n, res[1]);
487
488 fHTimeProb.SetBinContent(n, res[2]);
489
490 fHTimeLambda.SetBinContent(n, res[3]);
491 fHTimeLambda.SetBinError(n, res[4]);
492
493 //fHTimeNDF.SetBinContent(n, res[5]);
494
495 //
496 // Now prepare output
497 //
498 fParam->SetVal(res[0], res[1]);
499 fParam->SetReadyToSave();
500
501 *fTime = fLastTime;
502
503 // Include the current event
504 fTime->Plus1ns();
505
506 *fLog << fLastTime << ": Val=" << res[0] << " Err=" << res[1] << endl;
507}
508
509// --------------------------------------------------------------------------
510//
511// Fill the histogram
512//
513Bool_t MHEffectiveOnTime::Fill(const MParContainer *par, const Stat_t w)
514{
515 const MTime *time = dynamic_cast<const MTime*>(par);
516 if (!time)
517 {
518 *fLog << err << "ERROR - MHEffectiveOnTime::Fill without argument or container doesn't inherit from MTime... abort." << endl;
519 return kFALSE;
520 }
521
522 //
523 // If this is the first call we have to initialize the time-histogram
524 //
525 if (fLastTime==MTime())
526 {
527 MBinning bins;
528 bins.SetEdges(1, time->GetAxisTime()-fNumEvents/200, time->GetAxisTime());
529 bins.Apply(fHTimeEffOn);
530 bins.Apply(fHTimeProb);
531 bins.Apply(fHTimeLambda);
532 //bins.Apply(fHTimeNDF);
533
534 fParam->SetVal(0, 0);
535 fParam->SetReadyToSave();
536
537 *fTime = *time;
538
539 // Make this 1ns before the first event!
540 fTime->Minus1ns();
541 }
542
543 //
544 // Fill time difference into the histograms
545 //
546 const Double_t dt = *time-fLastTime;
547 fLastTime = *time;
548
549 fH2DeltaT.Fill(dt, fPointPos->GetZd(), w);
550 fH1DeltaT.Fill(dt, w);
551
552 //
553 // If we reached the event number limit for the time-bins fit the histogram
554 //
555 if (fH1DeltaT.GetEntries()>=fNumEvents)
556 FitTimeBin();
557
558 return kTRUE;
559}
560
561// --------------------------------------------------------------------------
562//
563// Fit the theta projections of the 2D histogram and the 1D Delta-T
564// distribution
565//
566Bool_t MHEffectiveOnTime::Finalize()
567{
568 FitThetaBins();
569 FitTimeBin();
570 MH::RemoveFirstBin(fHTimeEffOn);
571 MH::RemoveFirstBin(fHTimeProb);
572 MH::RemoveFirstBin(fHTimeLambda);
573 //MH::RemoveFirstBin(fHTimeNDF);
574
575 fIsFinalized = kTRUE;
576
577 return kTRUE;
578}
579
580// --------------------------------------------------------------------------
581//
582// Paint the integral and the error on top of the histogram
583//
584void MHEffectiveOnTime::PaintText(Double_t val, Double_t error) const
585{
586 TLatex text(0.45, 0.94, Form("T_{eff} = %.1fs \\pm %.1fs", val, error));
587 text.SetBit(TLatex::kTextNDC);
588 text.SetTextSize(0.04);
589 text.Paint();
590}
591
592void MHEffectiveOnTime::PaintText(Double_t *res) const
593{
594 TLatex text(0.25, 0.94, Form("T_{eff}=%.1fs\\pm%.1fs \\labda=%.1f\\pm%.1f p=%.1f%% NDF=%d",
595 res[0], res[1], res[3], res[4], res[2], res[5]));
596 text.SetBit(TLatex::kTextNDC);
597 text.SetTextSize(0.04);
598 text.Paint();
599}
600
601void MHEffectiveOnTime::PaintProb(TH1 &h) const
602{
603 Double_t sum = 0;
604 Int_t n = 0;
605 for (int i=0; i<h.GetNbinsX(); i++)
606 if (h.GetBinContent(i+1)>0)
607 {
608 sum += h.GetBinContent(i+1);
609 n++;
610 }
611
612 if (n==0)
613 return;
614
615 TLatex text(0.45, 0.94, Form("\\bar{p} = %.1f%% (n=%d)", sum/n, n));
616 text.SetBit(TLatex::kTextNDC);
617 text.SetTextSize(0.04);
618 text.Paint();
619}
620
621void MHEffectiveOnTime::UpdateRightAxis(TH1 &h)
622{
623 const Double_t max = h.GetMaximum()*1.1;
624 if (max==0)
625 return;
626
627 h.SetNormFactor(h.Integral()*gPad->GetUymax()/max);
628
629 TGaxis *axis = (TGaxis*)gPad->FindObject("RightAxis");
630 if (!axis)
631 return;
632
633 axis->SetX1(gPad->GetUxmax());
634 axis->SetX2(gPad->GetUxmax());
635 axis->SetY1(gPad->GetUymin());
636 axis->SetY2(gPad->GetUymax());
637 axis->SetWmax(max);
638}
639
640// --------------------------------------------------------------------------
641//
642// Prepare painting the histograms
643//
644void MHEffectiveOnTime::Paint(Option_t *opt)
645{
646 TH1D *h=0;
647 TPaveStats *st=0;
648
649 TString o(opt);
650 if (o==(TString)"fit")
651 {
652 TVirtualPad *pad = gPad;
653
654 for (int x=0; x<2; x++)
655 for (int y=0; y<3; y++)
656 {
657 TVirtualPad *p=gPad->GetPad(x+1)->GetPad(y+1);
658 if (!(st = (TPaveStats*)p->GetPrimitive("stats")))
659 continue;
660
661 if (st->GetOptStat()==11)
662 continue;
663
664 const Double_t y1 = st->GetY1NDC();
665 const Double_t y2 = st->GetY2NDC();
666 const Double_t x1 = st->GetX1NDC();
667 const Double_t x2 = st->GetX2NDC();
668
669 st->SetY1NDC((y2-y1)/3+y1);
670 st->SetX1NDC((x2-x1)/3+x1);
671 st->SetOptStat(11);
672 }
673
674 pad->GetPad(1)->cd(1);
675 if ((h = (TH1D*)gPad->FindObject(fNameProjDeltaT)))
676 {
677 h = fH2DeltaT.ProjectionX(fNameProjDeltaT, -1, 9999, "E");
678 if (h->GetEntries()>0)
679 gPad->SetLogy();
680 }
681
682 pad->GetPad(2)->cd(1);
683 if ((h = (TH1D*)gPad->FindObject(fNameProjTheta)))
684 fH2DeltaT.ProjectionY(fNameProjTheta, -1, 9999, "E");
685
686 if (!fIsFinalized)
687 FitThetaBins();
688 return;
689 }
690 if (o==(TString)"paint")
691 {
692 if ((h = (TH1D*)gPad->FindObject(fNameProjDeltaT)))
693 {
694 Double_t res[6];
695 FitH(h, res, kTRUE);
696 PaintText(res);
697 }
698 return;
699 }
700
701 if (o==(TString)"timendf")
702 {
703 // UpdateRightAxis(fHTimeNDF);
704 // FIXME: first bin?
705 PaintProb(fHTimeProb);
706 }
707
708 if (o==(TString)"thetandf")
709 {
710 UpdateRightAxis(fHThetaNDF);
711 // FIXME: first bin?
712 PaintProb(fHThetaProb);
713 }
714
715 h=0;
716 if (o==(TString)"theta")
717 {
718 h = &fHThetaEffOn;
719 UpdateRightAxis(fHThetaLambda);
720 }
721 if (o==(TString)"time")
722 {
723 h = &fHTimeEffOn;
724 UpdateRightAxis(fHTimeLambda);
725 }
726
727 if (!h)
728 return;
729
730 Double_t error = 0;
731 for (int i=0; i<h->GetXaxis()->GetNbins(); i++)
732 error += h->GetBinError(i);
733
734 PaintText(h->Integral(), error);
735}
736
737void MHEffectiveOnTime::DrawRightAxis(const char *title)
738{
739 TGaxis *axis = new TGaxis(gPad->GetUxmax(), gPad->GetUymin(),
740 gPad->GetUxmax(), gPad->GetUymax(),
741 0, 1, 510, "+L");
742 axis->SetName("RightAxis");
743 axis->SetTitle(title);
744 axis->SetTitleOffset(0.9);
745 axis->SetTextColor(kGreen);
746 axis->CenterTitle();
747 axis->SetBit(kCanDelete);
748 axis->Draw();
749}
750
751// --------------------------------------------------------------------------
752//
753// Draw the histogram
754//
755void MHEffectiveOnTime::Draw(Option_t *opt)
756{
757 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
758 pad->SetBorderMode(0);
759
760 AppendPad("fit");
761
762 pad->Divide(2, 1, 0, 0);
763
764 TH1 *h;
765
766 pad->cd(1);
767 gPad->SetBorderMode(0);
768 gPad->Divide(1, 3, 0, 0);
769 pad->GetPad(1)->cd(1);
770 gPad->SetBorderMode(0);
771 h = fH2DeltaT.ProjectionX(fNameProjDeltaT, -1, 9999, "E");
772 h->SetTitle("Distribution of \\Delta t [s]");
773 h->SetXTitle("\\Delta t [s]");
774 h->SetYTitle("Counts");
775 h->SetDirectory(NULL);
776 h->SetMarkerStyle(kFullDotMedium);
777 h->SetBit(kCanDelete);
778 h->Draw();
779 AppendPad("paint");
780
781 pad->GetPad(1)->cd(2);
782 gPad->SetBorderMode(0);
783 fHTimeProb.Draw();
784 AppendPad("timendf");
785 //fHTimeNDF.Draw("same");
786 //DrawRightAxis("NDF");
787
788 pad->GetPad(1)->cd(3);
789 gPad->SetBorderMode(0);
790 fHTimeEffOn.Draw();
791 AppendPad("time");
792 fHTimeLambda.Draw("same");
793 DrawRightAxis("\\lambda [s^{-1}]");
794
795 pad->cd(2);
796 gPad->SetBorderMode(0);
797 gPad->Divide(1, 3, 0, 0);
798
799 pad->GetPad(2)->cd(1);
800 gPad->SetBorderMode(0);
801 h = fH2DeltaT.ProjectionY(fNameProjTheta, -1, 9999, "E");
802 h->SetTitle("Distribution of \\Theta [\\circ]");
803 h->SetXTitle("\\Theta [\\circ]");
804 h->SetYTitle("Counts");
805 h->SetDirectory(NULL);
806 h->SetMarkerStyle(kFullDotMedium);
807 h->SetBit(kCanDelete);
808 h->GetYaxis()->SetTitleOffset(1.1);
809 h->Draw();
810
811 pad->GetPad(2)->cd(2);
812 gPad->SetBorderMode(0);
813 fHThetaProb.Draw();
814 AppendPad("thetandf");
815 fHThetaNDF.Draw("same");
816 DrawRightAxis("NDF");
817
818 pad->GetPad(2)->cd(3);
819 gPad->SetBorderMode(0);
820 fHThetaEffOn.Draw();
821 AppendPad("theta");
822 fHThetaLambda.Draw("same");
823 DrawRightAxis("\\lambda [s^{-1}]");
824}
Note: See TracBrowser for help on using the repository browser.