Context Navigation

source: trunk/MagicSoft/Mars/mhflux/MHAlpha.cc@ 5002

Visit:

Last change on this file since 5002 was 5002, checked in by tbretz, 20 years ago
* empty log message *
File size: 20.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, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19	!
20	! Copyright: MAGIC Software Development, 2000-2004
21	!
22	!
23	\* ======================================================================== */
24
25	//////////////////////////////////////////////////////////////////////////////
26	//
27	// MHAlpha
28	//
29	// Create a single Alpha-Plot. The alpha-plot is fitted online. You can
30	// check the result when it is filles in the MStatusDisplay
31	// For more information see MHFalseSource::FitSignificance
32	//
33	// For convinience (fit) the output significance is stored in a
34	// container in the parlisrt
35	//
36	// PRELIMINARY!
37	//
38	//////////////////////////////////////////////////////////////////////////////
39	#include "MHAlpha.h"
40
41	#include <TF1.h>
42	#include <TGraph.h>
43	#include <TStyle.h>
44	#include <TLatex.h>
45	#include <TCanvas.h>
46	#include <TPaveStats.h>
47	#include <TStopwatch.h>
48
49	#include "MGeomCam.h"
50	#include "MSrcPosCam.h"
51	#include "MHillasSrc.h"
52	#include "MEnergyEst.h"
53	#include "MTime.h"
54	#include "MObservatory.h"
55	#include "MPointingPos.h"
56	#include "MAstroSky2Local.h"
57	#include "MStatusDisplay.h"
58	#include "MParameters.h"
59	#include "MHMatrix.h"
60
61	#include "MBinning.h"
62	#include "MParList.h"
63
64	#include "MLog.h"
65	#include "MLogManip.h"
66
67	ClassImp(MHAlpha);
68	ClassImp(MAlphaFitter);
69
70	using namespace std;
71
72	// --------------------------------------------------------------------------
73	//
74	// Default Constructor
75	//
76	MHAlpha::MHAlpha(const char name, const char title)
77	: fResult(0), /fExcess(0),/ fEnergy(0), fPointPos(0), fTimeEffOn(0),
78	fTime(0), fNameProjAlpha(Form("Alpha_%p", this)), fMatrix(0)
79	{
80	//
81	// set the name and title of this object
82	//
83	fName = name ? name : "MHAlpha";
84	fTitle = title ? title : "Alpha plot";
85
86	fHAlpha.SetName("Alpha");
87	fHAlpha.SetTitle("Alpha");
88	fHAlpha.SetXTitle("\\Theta [deg]");
89	fHAlpha.SetYTitle("E_{est} [GeV]");
90	fHAlpha.SetZTitle("\|\\alpha\| [\\circ]");
91	fHAlpha.SetDirectory(NULL);
92	fHAlpha.UseCurrentStyle();
93
94	// Main histogram
95	fHAlphaTime.SetName("Alpha");
96	fHAlphaTime.SetXTitle("\|\\alpha\| [\\circ]");
97	fHAlphaTime.SetYTitle("Counts");
98	fHAlphaTime.UseCurrentStyle();
99	fHAlphaTime.SetDirectory(NULL);
100
101
102	fHEnergy.SetName("Energy");
103	fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
104	fHEnergy.SetXTitle("E_{est} [GeV]");
105	fHEnergy.SetYTitle("N_{exc}");
106	fHEnergy.SetDirectory(NULL);
107	fHEnergy.UseCurrentStyle();
108
109	fHTheta.SetName("Theta");
110	fHTheta.SetTitle(" N_{exc} vs. \\Theta ");
111	fHTheta.SetXTitle("\\Theta [\\circ]");
112	fHTheta.SetYTitle("N_{exc}");
113	fHTheta.SetDirectory(NULL);
114	fHTheta.UseCurrentStyle();
115
116	// effective on time versus time
117	fHTime.SetName("Time");
118	fHTime.SetTitle(" N_{exc} vs. Time ");
119	fHTime.SetXTitle("Time");
120	fHTime.SetYTitle("N_{exc} [s]");
121	fHTime.UseCurrentStyle();
122	fHTime.SetDirectory(NULL);
123	fHTime.GetYaxis()->SetTitleOffset(1.2);
124	fHTime.GetXaxis()->SetLabelSize(0.033);
125	fHTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
126	fHTime.GetXaxis()->SetTimeDisplay(1);
127	fHTime.Sumw2();
128
129	MBinning binsa, binse, binst;
130	binsa.SetEdges(18, 0, 90);
131	binse.SetEdgesLog(25, 10, 100000);
132	binst.SetEdgesCos(50, 0, 60);
133	binse.Apply(fHEnergy);
134	binst.Apply(fHTheta);
135	binsa.Apply(fHAlphaTime);
136
137	MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
138	}
139
140	void MHAlpha::FitEnergySpec(Bool_t paint)
141	{
142	TF1 f("Spectrum", "pow(x, [0])exp(-x/[1])[2]");
143	f.SetParameter(0, -2.0);
144	f.SetParameter(1, 500); // 50
145	f.SetParameter(2, 1); // 50
146	f.SetLineWidth(2);
147	f.SetLineColor(kGreen);
148
149	fHEnergy.Fit(&f, "0QI", "", 90, 1900); // Integral?
150
151	if (paint)
152	{
153	f.Paint("same");
154
155	TLatex text(0.2, 0.94, Form("\\alpha=%.1f E_{0}=%.1fGeV N=%.1f",
156	f.GetParameter(0),
157	f.GetParameter(1),
158	0/f.GetParameter(2)/));
159	text.SetBit(TLatex::kTextNDC);
160	text.SetTextSize(0.04);
161	text.Paint();
162	}
163	}
164
165	void MHAlpha::FitEnergyBins(Bool_t paint)
166	{
167	// Do not store the 'final' result in fFit
168	MAlphaFitter fit(fFit);
169
170	const Int_t n = fHAlpha.GetNbinsY();
171
172	for (int i=1; i<=n; i++)
173	{
174	TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", -1, 9999, i, i, i==1?"E":"");
175	if (fit.Fit(*h))
176	{
177	fHEnergy.SetBinContent(i, fit.GetExcessEvents());
178	fHEnergy.SetBinError(i, fit.GetExcessEvents()*0.2);
179	}
180	delete h;
181	}
182	}
183
184	void MHAlpha::FitThetaBins(Bool_t paint)
185	{
186	// Do not store the 'final' result in fFit
187	MAlphaFitter fit(fFit);
188
189	const Int_t n = fHAlpha.GetNbinsX();
190
191	for (int i=1; i<=n; i++)
192	{
193	TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", i, i, -1, 9999, i==1?"E":"");
194	if (fit.Fit(*h))
195	{
196	fHTheta.SetBinContent(i, fit.GetExcessEvents());
197	fHTheta.SetBinError(i, fit.GetExcessEvents()*0.2);
198	}
199	delete h;
200	}
201	}
202
203	Bool_t MHAlpha::SetupFill(const MParList *pl)
204	{
205	fHAlpha.Reset();
206	fHEnergy.Reset();
207	fHTheta.Reset();
208
209	fEnergy = (MEnergyEst*)pl->FindObject("MEnergyEst");
210	if (!fEnergy)
211	*fLog << warn << "MEnergyEst not found... ignored." << endl;
212
213	fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
214	if (!fPointPos)
215	*fLog << warn << "MPointingPos not found... ignored." << endl;
216
217	fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime");
218	if (!fTimeEffOn)
219	*fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;
220
221	fTime = (MTime*)pl->FindObject("MTime");
222	if (!fTime)
223	*fLog << warn << "MTime not found... ignored." << endl;
224
225	fResult = (MParameterD)const_cast<MParList>(pl)->FindCreateObj("MParameterD", "Significance");
226	if (!fResult)
227	return kFALSE;
228
229	//fExcess = (MParameterD)const_cast<MParList>(pl)->FindCreateObj("MParameterD", "MExcess");
230	//if (!fExcess)
231	// return kFALSE;
232
233	fLastTime = MTime();
234
235	MBinning binst, binse, binsa;
236	binst.SetEdges(fHAlpha, 'x');
237	binse.SetEdges(fHAlpha, 'y');
238	binsa.SetEdges(fHAlpha, 'z');
239
240	MBinning bins = (MBinning)pl->FindObject("BinningTheta", "MBinning");
241	if (fPointPos && bins)
242	binst.SetEdges(*bins);
243	if (!fPointPos)
244	binst.SetEdges(1, 0, 90);
245
246	bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning");
247	if (fEnergy && bins)
248	binse.SetEdges(*bins);
249	if (!fEnergy)
250	binse.SetEdges(1, 10, 100000);
251
252	bins = (MBinning*)pl->FindObject("BinningAlpha", "MBinning");
253	if (bins)
254	binsa.SetEdges(*bins);
255
256	binse.Apply(fHEnergy);
257	binst.Apply(fHTheta);
258	binsa.Apply(fHAlphaTime);
259	MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
260
261	MAlphaFitter fit = (MAlphaFitter)pl->FindObject("MAlphaFitter");
262	if (!fit)
263	*fLog << warn << "MAlphaFitter not found... ignored." << endl;
264	else
265	fit->Copy(fFit);
266
267	return kTRUE;
268	}
269
270	void MHAlpha::InitAlphaTime(const MTime &t)
271	{
272	//
273	// If this is the first call we have to initialize the time-histogram
274	//
275	MBinning bins;
276	bins.SetEdges(1, t.GetAxisTime()-60, t.GetAxisTime());
277	bins.Apply(fHTime);
278
279	fLastTime=t;
280	}
281
282	void MHAlpha::UpdateAlphaTime(Bool_t final)
283	{
284	if (!fTimeEffOn)
285	return;
286
287	const Int_t steps = 6;
288
289	static int rebin = steps;
290
291	if (!final)
292	{
293	if (fLastTime==MTime() && *fTimeEffOn!=MTime())
294	{
295	InitAlphaTime(*fTimeEffOn);
296	return;
297	}
298
299	if (fLastTime!=*fTimeEffOn)
300	{
301	fLastTime=*fTimeEffOn;
302	rebin--;
303	}
304
305	if (rebin>0)
306	return;
307	}
308
309	MAlphaFitter fit(fFit);
310	if (!fit.Fit(fHAlphaTime))
311	return;
312
313	// Reset Histogram
314	fHAlphaTime.Reset();
315
316	// Get number of bins
317	const Int_t n = fHTime.GetNbinsX();
318
319	//
320	// Prepare Histogram
321	//
322
323	MTime *time = final ? fTime : fTimeEffOn;
324	if (final)
325	time->Plus1ns();
326
327	// Enhance binning
328	MBinning bins;
329	bins.SetEdges(fHTime, 'x');
330	bins.AddEdge(time->GetAxisTime());
331	bins.Apply(fHTime);
332
333	//
334	// Fill histogram
335	//
336	fHTime.SetBinContent(n+1, fit.GetExcessEvents());
337	fHTime.SetBinError(n+1, fit.GetExcessEvents()*0.1);
338
339	fLog << all << fTimeEffOn << ": " << fit.GetExcessEvents() << endl;
340
341	rebin = steps;
342	}
343
344	// --------------------------------------------------------------------------
345	//
346	// Fill the histogram. For details see the code or the class description
347	//
348	Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
349	{
350	Double_t alpha, energy, theta;
351
352	if (fMatrix)
353	{
354	alpha = (*fMatrix)[fMap[0]];
355	energy = 1000; //(*fMatrix)[fMap[1]];
356	theta = 0; //(*fMatrix)[fMap[2]];
357	}
358	else
359	{
360	const MHillasSrc hil = dynamic_cast<const MHillasSrc>(par);
361	if (!par)
362	{
363	*fLog << err << dbginf << "MHillasSrc not found... abort." << endl;
364	return kFALSE;
365	}
366
367	alpha = hil->GetAlpha();
368	energy = fEnergy ? fEnergy->GetEnergy() : 1000;
369	theta = fPointPos ? fPointPos->GetZd() : 0;
370	}
371
372	UpdateAlphaTime();
373
374	fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w);
375	fHAlphaTime.Fill(TMath::Abs(alpha), w);
376
377	/*
378	// N_gamma vs Energy and Theta
379	const Double_t Ngam = fHUnfold.GetBinContent(m,n);
380	// A_eff vs Energy and Theta
381	const Double_t Aeff = fHAeff.GetBinContent(m,n);
382	// T_eff vs Theta
383	const Double_t Effon = teff.GetBinContent(n);
384
385	const Double_t c1 = fHUnfold.GetBinError(m,n)/Ngam;
386	const Double_t c2 = teff.GetBinError(n) /Effon;
387	const Double_t c3 = fHAeff.GetBinError(m,n) /Aeff;
388
389	const Double_t cont = Ngam / (DeltaE * Effon * Aeff);
390	const Double_t dcont = sqrt(c1c1 + c2c2 + c3*c3);
391
392	//
393	// Out of Range
394	//
395	const Bool_t oor = Ngam<=0 \|\| DeltaE<=0 \|\| Effon<=0 \|\| Aeff<=0;
396
397	if (oor)
398	*fLog << warn << "MHFlux::CalcAbsGammaFlux(" << m << "," << n << ") ";
399
400	if (Ngam<=0)
401	*fLog << " Ngam=0";
402	if (DeltaE<=0)
403	*fLog << " DeltaE=0";
404	if (Effon<=0)
405	*fLog << " Effon=0";
406	if (Aeff<=0)
407	*fLog << " Aeff=0";
408
409	if (oor)
410	*fLog << endl;
411
412	fHFlux.SetBinContent(m,n, oor ? 1e-20 : cont);
413	fHFlux.SetBinError(m,n, oor ? 1e-20 : dcont*cont);
414	*/
415	return kTRUE;
416	}
417
418	// --------------------------------------------------------------------------
419	//
420	// Paint the integral and the error on top of the histogram
421	//
422	void MHAlpha::PaintText(Double_t val, Double_t error) const
423	{
424	TLatex text(0.45, 0.94, Form("N_{exc} = %.1fs \\pm %.1fs", val, error));
425	text.SetBit(TLatex::kTextNDC);
426	text.SetTextSize(0.04);
427	text.Paint();
428	}
429
430	// --------------------------------------------------------------------------
431	//
432	// Update the projections
433	//
434	void MHAlpha::Paint(Option_t *opt)
435	{
436	// Note: Any object cannot be added to a pad twice!
437	// The object is searched by gROOT->FindObject only in
438	// gPad itself!
439	TVirtualPad *padsave = gPad;
440
441	TH1D *h0=0;
442
443	TString o(opt);
444	if (o==(TString)"proj")
445	{
446	TPaveStats *st=0;
447	for (int x=0; x<4; x++)
448	{
449	TVirtualPad *p=padsave->GetPad(x+1);
450	if (!p \|\| !(st = (TPaveStats*)p->GetPrimitive("stats")))
451	continue;
452
453	if (st->GetOptStat()==11)
454	continue;
455
456	const Double_t y1 = st->GetY1NDC();
457	const Double_t y2 = st->GetY2NDC();
458	const Double_t x1 = st->GetX1NDC();
459	const Double_t x2 = st->GetX2NDC();
460
461	st->SetY1NDC((y2-y1)/3+y1);
462	st->SetX1NDC((x2-x1)/3+x1);
463	st->SetOptStat(11);
464	}
465
466	padsave->cd(1);
467	fHAlpha.ProjectionZ(fNameProjAlpha);
468	FitEnergyBins();
469	FitThetaBins();
470	}
471
472	if (o==(TString)"alpha")
473	if ((h0 = (TH1D*)gPad->FindObject(fNameProjAlpha)))
474	{
475	// Do not store the 'final' result in fFit
476	MAlphaFitter fit(fFit);
477	fit.Fit(*h0, kTRUE);
478	fit.PaintResult();
479	}
480
481	if (o==(TString)"time")
482	PaintText(fHTime.Integral(), 0);
483
484	if (o==(TString)"theta")
485	PaintText(fHTheta.Integral(), 0);
486
487	if (o==(TString)"energy")
488	{
489	if (fHEnergy.GetEntries()>0)
490	{
491	gPad->SetLogx();
492	gPad->SetLogy();
493	}
494	FitEnergySpec(kTRUE);
495	}
496
497	gPad = padsave;
498	}
499
500	// --------------------------------------------------------------------------
501	//
502	// Draw the histogram
503	//
504	void MHAlpha::Draw(Option_t *opt)
505	{
506	TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
507
508	// Do the projection before painting the histograms into
509	// the individual pads
510	AppendPad("proj");
511
512	pad->SetBorderMode(0);
513	pad->Divide(2,2);
514
515	TH1D *h=0;
516
517	pad->cd(1);
518	gPad->SetBorderMode(0);
519	h = fHAlpha.ProjectionZ(fNameProjAlpha, -1, 9999, -1, 9999, "E");
520	h->SetBit(TH1::kNoTitle);
521	h->SetXTitle("\\alpha [\\circ]");
522	h->SetYTitle("Counts");
523	h->SetDirectory(NULL);
524	h->SetMarkerStyle(kFullDotMedium);
525	h->SetBit(kCanDelete);
526	h->Draw();
527	// After the Alpha-projection has been drawn. Fit the histogram
528	// and paint the result into this pad
529	AppendPad("alpha");
530
531	if (fHEnergy.GetNbinsX()>1)
532	{
533	pad->cd(2);
534	gPad->SetBorderMode(0);
535	fHEnergy.Draw();
536	AppendPad("energy");
537	}
538	else
539	delete pad->GetPad(2);
540
541	if (fTimeEffOn && fTime \|\| fHTime.GetNbinsX()>1)
542	{
543	pad->cd(3);
544	gPad->SetBorderMode(0);
545	fHTime.Draw();
546	AppendPad("time");
547	}
548	else
549	delete pad->GetPad(3);
550
551	if (fHTheta.GetNbinsX()>1)
552	{
553	pad->cd(4);
554	gPad->SetBorderMode(0);
555	fHTheta.Draw();
556	AppendPad("theta");
557	}
558	else
559	delete pad->GetPad(4);
560
561	}
562
563	// --------------------------------------------------------------------------
564	//
565	// This is a preliminary implementation of a alpha-fit procedure for
566	// all possible source positions. It will be moved into its own
567	// more powerfull class soon.
568	//
569	// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
570	// [0]exp(-0.5((x-[1])/[2])^2) + [3] + [4]x + [5]x^2
571	// or
572	// Aexp(-0.5((x-mu)/sigma)^2) + a + bx + cx^2
573	//
574	// Parameter [1] is fixed to 0 while the alpha peak should be
575	// symmetric around alpha=0.
576	//
577	// Parameter [4] is fixed to 0 because the first derivative at
578	// alpha=0 should be 0, too.
579	//
580	// In a first step the background is fitted between bgmin and bgmax,
581	// while the parameters [0]=0 and [2]=1 are fixed.
582	//
583	// In a second step the signal region (alpha<sigmax) is fittet using
584	// the whole function with parameters [1], [3], [4] and [5] fixed.
585	//
586	// The number of excess and background events are calculated as
587	// s = int(hist, 0, 1.25*sigint)
588	// b = int(pol2(3), 0, 1.25*sigint)
589	//
590	// The Significance is calculated using the Significance() member
591	// function.
592	//
593	/*
594	Bool_t MAlphaFitter::Fit(TH1D &h, Bool_t paint)
595	{
596	Double_t sigint=fSigInt;
597	Double_t sigmax=fSigMax;
598	Double_t bgmin=fBgMin;
599	Double_t bgmax=fBgMax;
600	Byte_t polynom=fPolynom;
601
602	// Implementing the function yourself is only about 5% faster
603	TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
604	//TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90);
605	TArrayD maxpar(func.GetNpar());
606
607	func.FixParameter(1, 0);
608	func.FixParameter(4, 0);
609	func.SetParLimits(2, 0, 90);
610	func.SetParLimits(3, -1, 1);
611
612	const Double_t alpha0 = h.GetBinContent(1);
613	const Double_t alphaw = h.GetXaxis()->GetBinWidth(1);
614
615	// Check for the regios which is not filled...
616	if (alpha0==0)
617	return kFALSE; //*fLog << warn << "Histogram empty." << endl;
618
619	// First fit a polynom in the off region
620	func.FixParameter(0, 0);
621	func.FixParameter(2, 1);
622	func.ReleaseParameter(3);
623
624	for (int i=5; i<func.GetNpar(); i++)
625	func.ReleaseParameter(i);
626
627	// options : N do not store the function, do not draw
628	// I use integral of function in bin rather than value at bin center
629	// R use the range specified in the function range
630	// Q quiet mode
631	h.Fit(&func, "NQI", "", bgmin, bgmax);
632
633	fChiSqBg = func.GetChisquare()/func.GetNDF();
634
635	// ------------------------------------
636	if (paint)
637	{
638	func.SetRange(0, 90);
639	func.SetLineColor(kRed);
640	func.SetLineWidth(2);
641	func.Paint("same");
642	}
643	// ------------------------------------
644
645	func.ReleaseParameter(0);
646	//func.ReleaseParameter(1);
647	func.ReleaseParameter(2);
648	func.FixParameter(3, func.GetParameter(3));
649	for (int i=5; i<func.GetNpar(); i++)
650	func.FixParameter(i, func.GetParameter(i));
651
652	// Do not allow signals smaller than the background
653	const Double_t A = alpha0-func.GetParameter(3);
654	const Double_t dA = TMath::Abs(A);
655	func.SetParLimits(0, -dA4, dA4);
656	func.SetParLimits(2, 0, 90);
657
658	// Now fit a gaus in the on region on top of the polynom
659	func.SetParameter(0, A);
660	func.SetParameter(2, sigmax*0.75);
661
662	// options : N do not store the function, do not draw
663	// I use integral of function in bin rather than value at bin center
664	// R use the range specified in the function range
665	// Q quiet mode
666	h.Fit(&func, "NQI", "", 0, sigmax);
667
668	fChiSqSignal = func.GetChisquare()/func.GetNDF();
669	fSigmaGaus = func.GetParameter(2);
670
671	//const Bool_t ok = NDF>0 && chi2<2.5*NDF;
672
673	// ------------------------------------
674	if (paint)
675	{
676	func.SetLineColor(kGreen);
677	func.SetLineWidth(2);
678	func.Paint("same");
679	}
680	// ------------------------------------
681	const Double_t s = func.Integral(0, sigint)/alphaw;
682	func.SetParameter(0, 0);
683	func.SetParameter(2, 1);
684	const Double_t b = func.Integral(0, sigint)/alphaw;
685
686	fSignificance = MMath::SignificanceLiMaSigned(s, b);
687	//exc = s-b;
688
689	const Double_t uin = 1.25*sigint;
690	const Int_t bin = h.GetXaxis()->FindFixBin(uin);
691	fIntegralMax = h.GetBinLowEdge(bin+1);
692	fExcessEvents = h.Integral(0, bin)-func.Integral(0, fIntegralMax)/alphaw;
693
694	return kTRUE;
695	}
696
697	void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size) const
698	{
699	TLatex text(x, y, Form("\\sigma_{Li/Ma}=%.1f (\\alpha<%.1f\\circ) \\omega=%.1f\\circ E=%d (\\alpha<%.1f\\circ) (\\chi_{b}^{2}=%.1f \\chi_{s}^{2}=%.1f)",
700	fSignificance, fSigInt, fSigmaGaus,
701	(int)fExcessEvents, fIntegralMax,
702	fChiSqBg, fChiSqSignal));
703
704	text.SetBit(TLatex::kTextNDC);
705	text.SetTextSize(size);
706	text.Paint();
707	}
708	*/
709	Bool_t MHAlpha::Finalize()
710	{
711	// Store the final result in fFit
712	TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
713	Bool_t rc = fFit.Fit(*h);
714	delete h;
715	if (!rc)
716	{
717	*fLog << warn << "Histogram empty." << endl;
718	return kTRUE;
719	}
720
721	if (fResult)
722	fResult->SetVal(fFit.GetSignificance());
723
724	FitEnergyBins();
725	FitThetaBins();
726	UpdateAlphaTime(kTRUE);
727	MH::RemoveFirstBin(fHTime);
728
729	return kTRUE;
730	}
731
732	// --------------------------------------------------------------------------
733	//
734	// You can use this function if you want to use a MHMatrix instead of
735	// MMcEvt. This function adds all necessary columns to the
736	// given matrix. Afterward you should fill the matrix with the corresponding
737	// data (eg from a file by using MHMatrix::Fill). If you now loop
738	// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
739	// will take the values from the matrix instead of the containers.
740	//
741	void MHAlpha::InitMapping(MHMatrix *mat)
742	{
743	if (fMatrix)
744	return;
745
746	fMatrix = mat;
747
748	fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha");
749	//fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy");
750	}
751
752	void MHAlpha::StopMapping()
753	{
754	fMatrix = NULL;
755	}

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

Download in other formats: