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

Last change on this file since 5741 was 5620, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 23.7 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.6, 0.99 };
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, 200); // 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<3*NDF;
338
339 if (paint)
340 {
341 func.SetLineWidth(2);
342 func.SetLineColor(ok ? kGreen : kRed);
343 func.Paint("same");
344 }
345
346 const Double_t lambda = func.GetParameter(0);
347 //const Double_t N0 = func.GetParameter(1);
348 const Double_t prob = func.GetProb();
349
350 /*
351 *fLog << all << "Nm/lambda=" << Nm/lambda << " chi2/NDF=";
352 *fLog << (NDF ? chi2/NDF : 0.0) << " lambda=";
353 *fLog << lambda << " N0=" << N0 << endl;
354 */
355
356 Double_t emat[2][2];
357 gMinuit->mnemat((Double_t*)emat, 2);
358
359 const Double_t dldl = emat[0][0];
360 //const Double_t dN0dN0 = emat[1][1];
361
362 const Double_t teff = Nm/lambda;
363 const Double_t dteff = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm);
364 const Double_t dl = TMath::Sqrt(dldl);
365
366 //const Double_t kappa = Nm/N0;
367 //const Double_t Rdead = 1.0 - kappa;
368 //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
369
370 // the effective on time is Nm/lambda
371 res[0] = teff;
372 res[1] = dteff;
373
374 // plot chi2-probability of fit
375 res[2] = prob*100;
376
377 // lambda of fit
378 res[3] = lambda;
379 res[4] = dl;
380
381 // NDoF of fit
382 res[5] = NDF;
383
384 // Chi2
385 res[6] = chi2;
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 ok;
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[7];
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[7];
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
554 // histogram - if it fails try again when 1.6% more events available
555 //
556 const Int_t n = (Int_t)fH1DeltaT.GetEntries();
557 if (n>=fNumEvents && n%(fNumEvents/60)==0)
558 FitTimeBin();
559
560 return kTRUE;
561}
562
563// --------------------------------------------------------------------------
564//
565// Fit the theta projections of the 2D histogram and the 1D Delta-T
566// distribution
567//
568Bool_t MHEffectiveOnTime::Finalize()
569{
570 FitThetaBins();
571 FitTimeBin();
572 MH::RemoveFirstBin(fHTimeEffOn);
573 MH::RemoveFirstBin(fHTimeProb);
574 MH::RemoveFirstBin(fHTimeLambda);
575 //MH::RemoveFirstBin(fHTimeNDF);
576
577 fIsFinalized = kTRUE;
578
579 return kTRUE;
580}
581
582// --------------------------------------------------------------------------
583//
584// Paint the integral and the error on top of the histogram
585//
586void MHEffectiveOnTime::PaintText(Double_t val, Double_t error) const
587{
588 TLatex text(0.45, 0.94, Form("T_{eff} = %.1fs \\pm %.1fs", val, error));
589 text.SetBit(TLatex::kTextNDC);
590 text.SetTextSize(0.04);
591 text.Paint();
592}
593
594void MHEffectiveOnTime::PaintText(Double_t *res) const
595{
596 TLatex text(0.27, 0.94, Form("T_{eff}=%.1fs\\pm%.1fs \\lambda=%.1f\\pm%.1fHz p=%.1f%% \\chi^{2}/%d=%.1f",
597 res[0], res[1], res[3], res[4], res[2], (int)res[5], res[6]/res[5]));
598 text.SetBit(TLatex::kTextNDC);
599 text.SetTextSize(0.04);
600 text.Paint();
601}
602
603void MHEffectiveOnTime::PaintProb(TH1 &h) const
604{
605 Double_t sum = 0;
606 Int_t n = 0;
607 for (int i=0; i<h.GetNbinsX(); i++)
608 if (h.GetBinContent(i+1)>0)
609 {
610 sum += h.GetBinContent(i+1);
611 n++;
612 }
613
614 if (n==0)
615 return;
616
617 TLatex text(0.47, 0.94, Form("\\bar{p} = %.1f%%", sum/n));
618 text.SetBit(TLatex::kTextNDC);
619 text.SetTextSize(0.04);
620 text.Paint();
621}
622
623void MHEffectiveOnTime::UpdateRightAxis(TH1 &h)
624{
625 const Double_t max = h.GetMaximum()*1.1;
626 if (max==0)
627 return;
628
629 h.SetNormFactor(h.Integral()*gPad->GetUymax()/max);
630
631 TGaxis *axis = (TGaxis*)gPad->FindObject("RightAxis");
632 if (!axis)
633 return;
634
635 axis->SetX1(gPad->GetUxmax());
636 axis->SetX2(gPad->GetUxmax());
637 axis->SetY1(gPad->GetUymin());
638 axis->SetY2(gPad->GetUymax());
639 axis->SetWmax(max);
640}
641
642// --------------------------------------------------------------------------
643//
644// Prepare painting the histograms
645//
646void MHEffectiveOnTime::Paint(Option_t *opt)
647{
648 *fLog << all << "Paint: '" << opt << "'" << endl;
649
650 TH1D *h=0;
651 TPaveStats *st=0;
652
653 TString o(opt);
654 if (o==(TString)"fit")
655 {
656 TVirtualPad *pad = gPad;
657
658 for (int x=0; x<2; x++)
659 for (int y=0; y<3; y++)
660 {
661 TVirtualPad *p=gPad->GetPad(x+1)->GetPad(y+1);
662 if (!(st = (TPaveStats*)p->GetPrimitive("stats")))
663 continue;
664
665 if (st->GetOptStat()==11)
666 continue;
667
668 const Double_t y1 = st->GetY1NDC();
669 const Double_t y2 = st->GetY2NDC();
670 const Double_t x1 = st->GetX1NDC();
671 const Double_t x2 = st->GetX2NDC();
672
673 st->SetY1NDC((y2-y1)/3+y1);
674 st->SetX1NDC((x2-x1)/3+x1);
675 st->SetOptStat(11);
676 }
677
678 pad->GetPad(1)->cd(1);
679 if ((h = (TH1D*)gPad->FindObject("ProjDeltaT"/*fNameProjDeltaT*/)))
680 {
681 h = fH2DeltaT.ProjectionX("ProjDeltaT"/*fNameProjDeltaT*/, -1, 9999, "E");
682 if (h->GetEntries()>0)
683 gPad->SetLogy();
684 }
685
686 pad->GetPad(2)->cd(1);
687 if ((h = (TH1D*)gPad->FindObject("ProjTheta"/*fNameProjTheta*/)))
688 fH2DeltaT.ProjectionY("ProjTheta"/*fNameProjTheta*/, -1, 9999, "E");
689
690 if (!fIsFinalized)
691 FitThetaBins();
692 return;
693 }
694 if (o==(TString)"paint")
695 {
696 if ((h = (TH1D*)gPad->FindObject("ProjDeltaT"/*fNameProjDeltaT*/)))
697 {
698 Double_t res[7];
699 FitH(h, res, kTRUE);
700 PaintText(res);
701 }
702 return;
703 }
704
705 if (o==(TString)"timendf")
706 {
707 // UpdateRightAxis(fHTimeNDF);
708 // FIXME: first bin?
709 PaintProb(fHTimeProb);
710 }
711
712 if (o==(TString)"thetandf")
713 {
714 UpdateRightAxis(fHThetaNDF);
715 // FIXME: first bin?
716 PaintProb(fHThetaProb);
717 }
718
719 h=0;
720 if (o==(TString)"theta")
721 {
722 h = &fHThetaEffOn;
723 UpdateRightAxis(fHThetaLambda);
724 }
725 if (o==(TString)"time")
726 {
727 h = &fHTimeEffOn;
728 UpdateRightAxis(fHTimeLambda);
729 }
730
731 if (!h)
732 return;
733
734 Double_t error = 0;
735 for (int i=0; i<h->GetXaxis()->GetNbins(); i++)
736 error += h->GetBinError(i);
737
738 PaintText(h->Integral(), error);
739}
740
741void MHEffectiveOnTime::DrawRightAxis(const char *title)
742{
743 TGaxis *axis = new TGaxis(gPad->GetUxmax(), gPad->GetUymin(),
744 gPad->GetUxmax(), gPad->GetUymax(),
745 0, 1, 510, "+L");
746 axis->SetName("RightAxis");
747 axis->SetTitle(title);
748 axis->SetTitleOffset(0.9);
749 axis->SetTextColor(kGreen);
750 axis->CenterTitle();
751 axis->SetBit(kCanDelete);
752 axis->Draw();
753}
754
755// --------------------------------------------------------------------------
756//
757// Draw the histogram
758//
759void MHEffectiveOnTime::Draw(Option_t *opt)
760{
761 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
762 pad->SetBorderMode(0);
763
764 AppendPad("fit");
765
766 pad->Divide(2, 1, 0, 0);
767
768 TH1 *h;
769
770 pad->cd(1);
771 gPad->SetBorderMode(0);
772 gPad->Divide(1, 3, 0, 0);
773 pad->GetPad(1)->cd(1);
774 gPad->SetBorderMode(0);
775 h = fH2DeltaT.ProjectionX("ProjDeltaT"/*fNameProjDeltaT*/, -1, 9999, "E");
776 h->SetTitle("Distribution of \\Delta t [s]");
777 h->SetXTitle("\\Delta t [s]");
778 h->SetYTitle("Counts");
779 h->SetDirectory(NULL);
780 h->SetMarkerStyle(kFullDotMedium);
781 h->SetBit(kCanDelete);
782 h->Draw();
783 AppendPad("paint");
784
785 pad->GetPad(1)->cd(2);
786 gPad->SetBorderMode(0);
787 fHTimeProb.Draw();
788 AppendPad("timendf");
789 //fHTimeNDF.Draw("same");
790 //DrawRightAxis("NDF");
791
792 pad->GetPad(1)->cd(3);
793 gPad->SetBorderMode(0);
794 fHTimeEffOn.Draw();
795 AppendPad("time");
796 fHTimeLambda.Draw("same");
797 DrawRightAxis("\\lambda [s^{-1}]");
798
799 pad->cd(2);
800 gPad->SetBorderMode(0);
801 gPad->Divide(1, 3, 0, 0);
802
803 pad->GetPad(2)->cd(1);
804 gPad->SetBorderMode(0);
805 h = fH2DeltaT.ProjectionY("ProjTheta"/*fNameProjTheta*/, -1, 9999, "E");
806 h->SetTitle("Distribution of \\Theta [\\circ]");
807 h->SetXTitle("\\Theta [\\circ]");
808 h->SetYTitle("Counts");
809 h->SetDirectory(NULL);
810 h->SetMarkerStyle(kFullDotMedium);
811 h->SetBit(kCanDelete);
812 h->GetYaxis()->SetTitleOffset(1.1);
813 h->Draw();
814
815 pad->GetPad(2)->cd(2);
816 gPad->SetBorderMode(0);
817 fHThetaProb.Draw();
818 AppendPad("thetandf");
819 fHThetaNDF.Draw("same");
820 DrawRightAxis("NDF");
821
822 pad->GetPad(2)->cd(3);
823 gPad->SetBorderMode(0);
824 fHThetaEffOn.Draw();
825 AppendPad("theta");
826 fHThetaLambda.Draw("same");
827 DrawRightAxis("\\lambda [s^{-1}]");
828}
Note: See TracBrowser for help on using the repository browser.