Context Navigation

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

Visit:

Last change on this file since 4983 was 4975, 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
99	ClassImp(MHEffectiveOnTime);
100
101	using namespace std;
102
103	// --------------------------------------------------------------------------
104	//
105	// Default Constructor. It initializes all histograms.
106	//
107	MHEffectiveOnTime::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("S");
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("S");
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	//
253	Bool_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	//
292	Bool_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	static 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	// plot chi2/NDF of fit
381	//res[3] = NDF ? chi2/NDF : 0.0;
382
383	// lambda of fit
384	res[3] = lambda;
385	res[4] = dl;
386
387	// NDoF of fit
388	res[5] = NDF;
389
390	// Rdead (from fit) is the fraction from real time lost by the dead time
391	//fHRdead.SetBinContent(i, Rdead);
392	//fHRdead.SetBinError (i,dRdead);
393
394	return kTRUE;
395	}
396
397	// --------------------------------------------------------------------------
398	//
399	// Fit a all bins of the distribution in theta. Store the result in the
400	// Theta-Histograms
401	//
402	void MHEffectiveOnTime::FitThetaBins()
403	{
404	fHThetaEffOn.Reset();
405	fHThetaProb.Reset();
406	fHThetaLambda.Reset();
407	fHThetaNDF.Reset();
408
409	const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999));
410
411	// nbins = number of Theta bins
412	const Int_t nbins = fH2DeltaT.GetNbinsY();
413
414	TH1D *h=0;
415	for (int i=1; i<=nbins; i++)
416	{
417	// TH1D &h = *hist->ProjectionX("Calc-theta", i, i);
418	h = fH2DeltaT.ProjectionX(name, i, i, "E");
419
420	Double_t res[6];
421	if (!FitH(h, res))
422	continue;
423
424	// the effective on time is Nm/lambda
425	fHThetaEffOn.SetBinContent(i, res[0]);
426	fHThetaEffOn.SetBinError (i, res[1]);
427
428	// plot chi2-probability of fit
429	fHThetaProb.SetBinContent(i, res[2]);
430
431	// plot chi2/NDF of fit
432	//fHChi2.SetBinContent(i, res[3]);
433
434	// lambda of fit
435	fHThetaLambda.SetBinContent(i, res[3]);
436	fHThetaLambda.SetBinError (i, res[4]);
437
438	// NDoF of fit
439	fHThetaNDF.SetBinContent(i, res[5]);
440
441	// Rdead (from fit) is the fraction from real time lost by the dead time
442	//fHRdead.SetBinContent(i, Rdead);
443	//fHRdead.SetBinError (i,dRdead);
444	}
445
446	// Histogram is reused via gROOT->FindObject()
447	// Need to be deleted only once
448	if (h)
449	delete h;
450	}
451
452	// --------------------------------------------------------------------------
453	//
454	// Fit the single-time-bin histogram. Store the result in the
455	// Time-Histograms
456	//
457	void MHEffectiveOnTime::FitTimeBin()
458	{
459	//
460	// Fit histogram
461	//
462	Double_t res[6];
463	if (!FitH(&fH1DeltaT, res))
464	return;
465
466	// Reset Histogram
467	fH1DeltaT.Reset();
468
469	//
470	// Prepare Histogram
471	//
472
473	// Get number of bins
474	const Int_t n = fHTimeEffOn.GetNbinsX();
475
476	// Enhance binning
477	MBinning bins;
478	bins.SetEdges(fHTimeEffOn, 'x');
479	bins.AddEdge(fLastTime.GetAxisTime());
480	bins.Apply(fHTimeEffOn);
481	bins.Apply(fHTimeProb);
482	bins.Apply(fHTimeLambda);
483	//bins.Apply(fHTimeNDF);
484
485	//
486	// Fill histogram
487	//
488	fHTimeEffOn.SetBinContent(n, res[0]);
489	fHTimeEffOn.SetBinError(n, res[1]);
490
491	fHTimeProb.SetBinContent(n, res[2]);
492
493	fHTimeLambda.SetBinContent(n, res[3]);
494	fHTimeLambda.SetBinError(n, res[4]);
495
496	//fHTimeNDF.SetBinContent(n, res[5]);
497
498	//
499	// Now prepare output
500	//
501	fParam->SetVal(res[0], res[1]);
502	fParam->SetReadyToSave();
503
504	*fTime = fLastTime;
505
506	// Include the current event
507	fTime->Plus1ns();
508
509	*fLog << fLastTime << ": Val=" << res[0] << " Err=" << res[1] << endl;
510	}
511
512	// --------------------------------------------------------------------------
513	//
514	// Fill the histogram
515	//
516	Bool_t MHEffectiveOnTime::Fill(const MParContainer *par, const Stat_t w)
517	{
518	const MTime time = dynamic_cast<const MTime>(par);
519	if (!time)
520	{
521	*fLog << err << "ERROR - MHEffectiveOnTime::Fill without argument or container doesn't inherit from MTime... abort." << endl;
522	return kFALSE;
523	}
524
525	//
526	// If this is the first call we have to initialize the time-histogram
527	//
528	if (fLastTime==MTime())
529	{
530	MBinning bins;
531	bins.SetEdges(1, time->GetAxisTime()-fNumEvents/200, time->GetAxisTime());
532	bins.Apply(fHTimeEffOn);
533	bins.Apply(fHTimeProb);
534	bins.Apply(fHTimeLambda);
535	//bins.Apply(fHTimeNDF);
536
537	fParam->SetVal(0, 0);
538	fParam->SetReadyToSave();
539
540	fTime = time;
541
542	// Make this 1ns before the first event!
543	fTime->Minus1ns();
544	}
545
546	//
547	// Fill time difference into the histograms
548	//
549	const Double_t dt = *time-fLastTime;
550
551	fH2DeltaT.Fill(dt, fPointPos->GetZd(), w);
552	fH1DeltaT.Fill(dt, w);
553
554	fLastTime = *time;
555
556	//
557	// If we reached the event number limit for the time-bins fit the histogram
558	//
559	if (fH1DeltaT.GetEntries()>=fNumEvents)
560	FitTimeBin();
561
562	return kTRUE;
563	}
564
565	// --------------------------------------------------------------------------
566	//
567	// Fit the theta projections of the 2D histogram and the 1D Delta-T
568	// distribution
569	//
570	Bool_t MHEffectiveOnTime::Finalize()
571	{
572	FitThetaBins();
573	FitTimeBin();
574	MH::RemoveFirstBin(fHTimeEffOn);
575	MH::RemoveFirstBin(fHTimeProb);
576	MH::RemoveFirstBin(fHTimeLambda);
577	//MH::RemoveFirstBin(fHTimeNDF);
578
579	fIsFinalized = kTRUE;
580
581	return kTRUE;
582	}
583
584	// --------------------------------------------------------------------------
585	//
586	// Paint the integral and the error on top of the histogram
587	//
588	void MHEffectiveOnTime::PaintText(Double_t val, Double_t error) const
589	{
590	TLatex text(0.45, 0.94, Form("T_{eff} = %.1fs \\pm %.1fs", val, error));
591	text.SetBit(TLatex::kTextNDC);
592	text.SetTextSize(0.04);
593	text.Paint();
594	}
595
596	void MHEffectiveOnTime::PaintProb(TH1 &h) const
597	{
598	Double_t sum = 0;
599	Int_t n = 0;
600	for (int i=0; i<h.GetNbinsX(); i++)
601	if (h.GetBinContent(i+1)>0)
602	{
603	sum += h.GetBinContent(i+1);
604	n++;
605	}
606
607	if (n==0)
608	return;
609
610	TLatex text(0.45, 0.94, Form("\\bar{p} = %.1f%% (n=%d)", sum/n, n));
611	text.SetBit(TLatex::kTextNDC);
612	text.SetTextSize(0.04);
613	text.Paint();
614	}
615
616	void MHEffectiveOnTime::UpdateRightAxis(TH1 &h)
617	{
618	cout << "Update... " << flush;
619	const Double_t max = h.GetMaximum()*1.1;
620	if (max==0)
621	return;
622
623	h.SetNormFactor(h.Integral()*gPad->GetUymax()/max);
624
625	TGaxis axis = (TGaxis)gPad->FindObject("RightAxis");
626	if (!axis)
627	return;
628
629	cout << axis << " ";
630	cout << "X: " << gPad->GetUxmax() << " ";
631	cout << "Y: " << gPad->GetUymin() << "/" << gPad->GetUymax() << " ";
632	cout << "max=" << max << endl;
633
634	axis->SetX1(gPad->GetUxmax());
635	axis->SetX2(gPad->GetUxmax());
636	axis->SetY1(gPad->GetUymin());
637	axis->SetY2(gPad->GetUymax());
638	axis->SetWmax(max);
639	}
640
641	// --------------------------------------------------------------------------
642	//
643	// Prepare painting the histograms
644	//
645	void MHEffectiveOnTime::Paint(Option_t *opt)
646	{
647	TH1D *h=0;
648	TPaveStats *st=0;
649
650	TString o(opt);
651	if (o==(TString)"fit")
652	{
653	TVirtualPad *pad = gPad;
654
655	for (int x=0; x<2; x++)
656	for (int y=0; y<3; y++)
657	{
658	TVirtualPad *p=gPad->GetPad(x+1)->GetPad(y+1);
659	if (!(st = (TPaveStats*)p->GetPrimitive("stats")))
660	continue;
661
662	if (st->GetOptStat()==11)
663	continue;
664
665	const Double_t y1 = st->GetY1NDC();
666	const Double_t y2 = st->GetY2NDC();
667	const Double_t x1 = st->GetX1NDC();
668	const Double_t x2 = st->GetX2NDC();
669
670	st->SetY1NDC((y2-y1)/3+y1);
671	st->SetX1NDC((x2-x1)/3+x1);
672	st->SetOptStat(11);
673	}
674
675	pad->GetPad(1)->cd(1);
676	if ((h = (TH1D*)gPad->FindObject(fNameProjDeltaT)))
677	{
678	h = fH2DeltaT.ProjectionX(fNameProjDeltaT, -1, 9999, "E");
679	if (h->GetEntries()>0)
680	gPad->SetLogy();
681	}
682
683	pad->GetPad(2)->cd(1);
684	if ((h = (TH1D*)gPad->FindObject(fNameProjTheta)))
685	fH2DeltaT.ProjectionY(fNameProjTheta, -1, 9999, "E");
686
687	if (!fIsFinalized)
688	FitThetaBins();
689	return;
690	}
691	if (o==(TString)"paint")
692	{
693	if ((h = (TH1D*)gPad->FindObject(fNameProjDeltaT)))
694	{
695	Double_t res[6];
696	if (FitH(h, res, kTRUE))
697	PaintText(res[0], res[1]);
698	}
699	return;
700	}
701
702	if (o==(TString)"timendf")
703	{
704	// UpdateRightAxis(fHTimeNDF);
705	// FIXME: first bin?
706	PaintProb(fHTimeProb);
707	}
708
709	if (o==(TString)"thetandf")
710	{
711	UpdateRightAxis(fHThetaNDF);
712	// FIXME: first bin?
713	PaintProb(fHThetaProb);
714	}
715
716	h=0;
717	if (o==(TString)"theta")
718	{
719	h = &fHThetaEffOn;
720	UpdateRightAxis(fHThetaLambda);
721	}
722	if (o==(TString)"time")
723	{
724	h = &fHTimeEffOn;
725	UpdateRightAxis(fHTimeLambda);
726	}
727
728	if (!h)
729	return;
730
731	Double_t error = 0;
732	for (int i=0; i<h->GetXaxis()->GetNbins(); i++)
733	error += h->GetBinError(i);
734
735	PaintText(h->Integral(), error);
736	}
737
738	void MHEffectiveOnTime::DrawRightAxis(const char *title)
739	{
740	TGaxis *axis = new TGaxis(gPad->GetUxmax(), gPad->GetUymin(),
741	gPad->GetUxmax(), gPad->GetUymax(),
742	0, 1, 510, "+L");
743	axis->SetName("RightAxis");
744	axis->SetTitle(title);
745	axis->SetTitleOffset(0.9);
746	axis->SetTextColor(kGreen);
747	axis->CenterTitle();
748	axis->SetBit(kCanDelete);
749	axis->Draw();
750	}
751
752	// --------------------------------------------------------------------------
753	//
754	// Draw the histogram
755	//
756	void MHEffectiveOnTime::Draw(Option_t *opt)
757	{
758	TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
759	pad->SetBorderMode(0);
760
761	AppendPad("fit");
762
763	pad->Divide(2, 1, 0, 0);
764
765	TH1 *h;
766
767	pad->cd(1);
768	gPad->SetBorderMode(0);
769	gPad->Divide(1, 3, 0, 0);
770	pad->GetPad(1)->cd(1);
771	gPad->SetBorderMode(0);
772	h = fH2DeltaT.ProjectionX(fNameProjDeltaT, -1, 9999, "E");
773	h->SetTitle("Distribution of \\Delta t [s]");
774	h->SetXTitle("\\Delta t [s]");
775	h->SetYTitle("Counts");
776	h->SetDirectory(NULL);
777	h->SetMarkerStyle(kFullDotMedium);
778	h->SetBit(kCanDelete);
779	h->Draw();
780	AppendPad("paint");
781
782	pad->GetPad(1)->cd(2);
783	gPad->SetBorderMode(0);
784	fHTimeProb.Draw();
785	AppendPad("timendf");
786	//fHTimeNDF.Draw("same");
787	//DrawRightAxis("NDF");
788
789	pad->GetPad(1)->cd(3);
790	gPad->SetBorderMode(0);
791	fHTimeEffOn.Draw();
792	AppendPad("time");
793	fHTimeLambda.Draw("same");
794	DrawRightAxis("Slope");
795
796	pad->cd(2);
797	gPad->SetBorderMode(0);
798	gPad->Divide(1, 3, 0, 0);
799
800	pad->GetPad(2)->cd(1);
801	gPad->SetBorderMode(0);
802	h = fH2DeltaT.ProjectionY(fNameProjTheta, -1, 9999, "E");
803	h->SetTitle("Distribution of \\Theta [\\circ]");
804	h->SetXTitle("\\Theta [\\circ]");
805	h->SetYTitle("Counts");
806	h->SetDirectory(NULL);
807	h->SetMarkerStyle(kFullDotMedium);
808	h->SetBit(kCanDelete);
809	h->GetYaxis()->SetTitleOffset(1.1);
810	h->Draw();
811
812	pad->GetPad(2)->cd(2);
813	gPad->SetBorderMode(0);
814	fHThetaProb.Draw();
815	AppendPad("thetandf");
816	fHThetaNDF.Draw("same");
817	DrawRightAxis("NDF");
818
819	pad->GetPad(2)->cd(3);
820	gPad->SetBorderMode(0);
821	fHThetaEffOn.Draw();
822	AppendPad("theta");
823	fHThetaLambda.Draw("same");
824	DrawRightAxis("Slope");
825	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: