Context Navigation

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

Visit:

Last change on this file since 5700 was 5692, checked in by tbretz, 20 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, 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	: fOffData(0), fResult(0), /fExcess(0),/ fEnergy(0), fHillas(0),
78	fPointPos(0), fTimeEffOn(0), fTime(0),
79	fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE), fSkipHistEnergy(kFALSE),
80	//fEnergyMin(-1), fEnergyMax(-1), fSizeMin(-1), fSizeMax(-1),
81	fMatrix(0)
82	{
83	//
84	// set the name and title of this object
85	//
86	fName = name ? name : "MHAlpha";
87	fTitle = title ? title : "Alpha plot";
88
89	fHAlpha.SetName("Alpha");
90	fHAlpha.SetTitle("Alpha");
91	fHAlpha.SetXTitle("\\Theta [deg]");
92	fHAlpha.SetYTitle("E_{est} [GeV]");
93	fHAlpha.SetZTitle("\|\\alpha\| [\\circ]");
94	fHAlpha.SetDirectory(NULL);
95	fHAlpha.UseCurrentStyle();
96
97	// Main histogram
98	fHAlphaTime.SetName("Alpha");
99	fHAlphaTime.SetXTitle("\|\\alpha\| [\\circ]");
100	fHAlphaTime.SetYTitle("Counts");
101	fHAlphaTime.UseCurrentStyle();
102	fHAlphaTime.SetDirectory(NULL);
103
104
105	//fHEnergy.SetName("Energy");
106	//fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
107	//fHEnergy.SetXTitle("E_{est} [GeV]");
108	fHEnergy.SetYTitle("N_{exc}");
109	fHEnergy.SetDirectory(NULL);
110	fHEnergy.UseCurrentStyle();
111
112	fHTheta.SetName("Theta");
113	fHTheta.SetTitle(" N_{exc} vs. \\Theta ");
114	fHTheta.SetXTitle("\\Theta [\\circ]");
115	fHTheta.SetYTitle("N_{exc}");
116	fHTheta.SetDirectory(NULL);
117	fHTheta.UseCurrentStyle();
118
119	// effective on time versus time
120	fHTime.SetName("Time");
121	fHTime.SetTitle(" N_{exc} vs. Time ");
122	fHTime.SetXTitle("Time");
123	fHTime.SetYTitle("N_{exc} [s]");
124	fHTime.UseCurrentStyle();
125	fHTime.SetDirectory(NULL);
126	fHTime.GetYaxis()->SetTitleOffset(1.2);
127	fHTime.GetXaxis()->SetLabelSize(0.033);
128	fHTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
129	fHTime.GetXaxis()->SetTimeDisplay(1);
130	fHTime.Sumw2();
131
132	MBinning binsa, binse, binst;
133	binsa.SetEdges(18, 0, 90);
134	binse.SetEdgesLog(25, 10, 100000);
135	binst.SetEdgesCos(50, 0, 60);
136	binse.Apply(fHEnergy);
137	binst.Apply(fHTheta);
138	binsa.Apply(fHAlphaTime);
139
140	MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
141	}
142
143	void MHAlpha::FitEnergySpec(Bool_t paint)
144	{
145	/*
146	TF1 f("Spectrum", "pow(x, [0])exp(-x/[1])[2]*x");//
147	f.SetParameter(0, -2.1);
148	f.SetParameter(1, 1400); // 50
149	f.SetParameter(2, fHEnergy.Integral()); // 50
150	f.SetLineWidth(2);
151
152	fHEnergy.Fit(&f, "0QI", "", 100, 2400); // Integral? 90, 2800
153
154	const Double_t chi2 = f.GetChisquare();
155	const Int_t NDF = f.GetNDF();
156
157	// was fit successful ?
158	const Bool_t ok = NDF>0 && chi2<2.5*NDF;
159	f.SetLineColor(ok ? kGreen : kRed);
160
161	if (paint)
162	{
163	f.Paint("same");
164
165	if (ok)
166	{
167	TString s;
168	s += Form("\\alpha=%.1f ", f.GetParameter(0));
169	s += Form("E_{c}=%.1f TeV ", f.GetParameter(1)/1000);
170	s += Form("N=%.1e", f.GetParameter(2));
171
172	TLatex text(0.25, 0.94, s.Data());
173	text.SetBit(TLatex::kTextNDC);
174	text.SetTextSize(0.04);
175	text.Paint();
176	}
177	}*/
178	}
179
180	void MHAlpha::FitEnergyBins(Bool_t paint)
181	{
182	// Do not store the 'final' result in fFit
183	MAlphaFitter fit(fFit);
184
185	const Int_t n = fHAlpha.GetNbinsY();
186
187	for (int i=1; i<=n; i++)
188	{
189	if (fit.FitEnergy(fHAlpha, fOffData, i))
190	{
191	fHEnergy.SetBinContent(i, fit.GetEventsExcess());
192	fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
193	}
194	}
195	}
196
197	void MHAlpha::FitThetaBins(Bool_t paint)
198	{
199	// Do not store the 'final' result in fFit
200	MAlphaFitter fit(fFit);
201
202	const Int_t n = fHAlpha.GetNbinsX();
203
204	for (int i=1; i<=n; i++)
205	{
206	if (fit.FitTheta(fHAlpha, fOffData, i))
207	{
208	fHTheta.SetBinContent(i, fit.GetEventsExcess());
209	fHTheta.SetBinError(i, fit.GetEventsExcess()*0.2);
210	}
211	}
212	}
213
214	Bool_t MHAlpha::SetupFill(const MParList *pl)
215	{
216	fHAlpha.Reset();
217	fHEnergy.Reset();
218	fHTheta.Reset();
219
220	if (fName!=(TString)"MHAlphaOff")
221	{
222	MHAlpha hoff = (MHAlpha)pl->FindObject("MHAlphaOff");
223	if (!hoff)
224	*fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl;
225	else
226	{
227	*fLog << inf << "MHAlphaOff [MHAlpha] found... using on-off mode!" << endl;
228	SetOffData(*hoff);
229	}
230	}
231
232	fHillas = 0;
233	/*
234	if (fSizeMin>=0 \|\| fSizeMax>=0)
235	{
236	fHillas = (MHillas*)pl->FindObject("MHillas");
237	if (!fHillas)
238	{
239	*fLog << warn << "Size cut set, but MHillas not found... abort." << endl;
240	return kFALSE;
241	}
242	}
243	*/
244	fEnergy = (MEnergyEst*)pl->FindObject("MEnergyEst");
245	if (!fEnergy)
246	{ /*
247	if (fEnergyMin>=0 \|\| fEnergyMax>=0)
248	{
249	*fLog << warn << "Energy cut set, but MEnergyEst not found... abort." << endl;
250	return kFALSE;
251	} */
252
253	*fLog << warn << "MEnergyEst not found... " << flush;
254
255	if (!fHillas)
256	fHillas = (MHillas*)pl->FindObject("MHillas");
257	if (fHillas)
258	*fLog << "using SIZE instead!" << endl;
259	else
260	*fLog << "ignored." << endl;
261
262	fHEnergy.SetName("Size");
263	fHEnergy.SetTitle(" N_{exc} vs. Size ");
264	fHEnergy.SetXTitle("Size [\\gamma]");
265	}
266	else
267	{
268	fHEnergy.SetName("Energy");
269	fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
270	fHEnergy.SetXTitle("E_{est} [GeV]");
271	}
272
273	fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
274	if (!fPointPos)
275	*fLog << warn << "MPointingPos not found... ignored." << endl;
276
277	fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime");
278	if (!fTimeEffOn)
279	*fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;
280	else
281	*fTimeEffOn = MTime(); // FIXME: How to do it different?
282
283	fTime = (MTime*)pl->FindObject("MTime");
284	if (!fTime)
285	*fLog << warn << "MTime not found... ignored." << endl;
286
287	fResult = (MParameterD)const_cast<MParList>(pl)->FindCreateObj("MParameterD", "Significance");
288	if (!fResult)
289	return kFALSE;
290
291	//fExcess = (MParameterD)const_cast<MParList>(pl)->FindCreateObj("MParameterD", "MExcess");
292	//if (!fExcess)
293	// return kFALSE;
294
295	fLastTime = MTime();
296
297	MBinning binst, binse, binsa;
298	binst.SetEdges(fOffData ? *fOffData : fHAlpha, 'x');
299	binse.SetEdges(fOffData ? *fOffData : fHAlpha, 'y');
300	binsa.SetEdges(fOffData ? *fOffData : fHAlpha, 'z');
301
302	if (!fOffData)
303	{
304	MBinning bins = (MBinning)pl->FindObject("BinningTheta", "MBinning");
305	if (fPointPos && bins)
306	binst.SetEdges(*bins);
307	if (!fPointPos)
308	binst.SetEdges(1, 0, 60);
309
310	bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning");
311	if ((fEnergy\|\|fHillas) && bins)
312	binse.SetEdges(*bins);
313	if (!fEnergy && !fHillas)
314	binse.SetEdges(1, 10, 100000);
315
316	bins = (MBinning*)pl->FindObject("BinningAlpha", "MBinning");
317	if (bins)
318	binsa.SetEdges(*bins);
319	}
320
321	binse.Apply(fHEnergy);
322	binst.Apply(fHTheta);
323	binsa.Apply(fHAlphaTime);
324	MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
325
326	MAlphaFitter fit = (MAlphaFitter)pl->FindObject("MAlphaFitter");
327	if (!fit)
328	*fLog << warn << "MAlphaFitter not found... ignored." << endl;
329	else
330	fit->Copy(fFit);
331
332	fFit.Print();
333
334	return kTRUE;
335	}
336
337	void MHAlpha::InitAlphaTime(const MTime &t)
338	{
339	//
340	// If this is the first call we have to initialize the time-histogram
341	//
342	MBinning bins;
343	bins.SetEdges(1, t.GetAxisTime()-60, t.GetAxisTime());
344	bins.Apply(fHTime);
345
346	fLastTime=t;
347	}
348
349	void MHAlpha::UpdateAlphaTime(Bool_t final)
350	{
351	if (!fTimeEffOn)
352	return;
353
354	const Int_t steps = 10;
355
356	static int rebin = steps;
357
358	if (!final)
359	{
360	if (fLastTime==MTime() && *fTimeEffOn!=MTime())
361	{
362	InitAlphaTime(*fTimeEffOn);
363	return;
364	}
365
366	if (fLastTime!=*fTimeEffOn)
367	{
368	fLastTime=*fTimeEffOn;
369	rebin--;
370	}
371
372	if (rebin>0)
373	return;
374	}
375
376	// Check new 'last time'
377	MTime *time = final ? fTime : fTimeEffOn;
378
379	if (time->GetAxisTime()<=fHTime.GetXaxis()->GetXmax())
380	{
381	fLog << warn << "WARNING - New time-stamp " << time << " lower" << endl;
382	*fLog << "than upper edge of histogram... skipped." << endl;
383	*fLog << "This should not happen. Maybe you started you eventloop with" << endl;
384	*fLog << "an already initialized time stamp MTimeEffectiveOnTime?" << endl;
385	rebin++;
386	return;
387	}
388
389	// Fit time histogram
390	MAlphaFitter fit(fFit);
391
392	TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", -1, 9999, -1, 9999, "E") : 0;
393	const Bool_t rc = fit.Fit(fHAlphaTime, h);
394	if (h)
395	delete h;
396
397	if (!rc)
398	return;
399
400	// Reset Histogram
401	fHAlphaTime.Reset();
402
403	// Get number of bins
404	const Int_t n = fHTime.GetNbinsX();
405
406	//
407	// Prepare Histogram
408	//
409	if (final)
410	time->Plus1ns();
411
412	// Enhance binning
413	MBinning bins;
414	bins.SetEdges(fHTime, 'x');
415	bins.AddEdge(time->GetAxisTime());
416	bins.Apply(fHTime);
417
418	//
419	// Fill histogram
420	//
421	fHTime.SetBinContent(n+1, fit.GetEventsExcess());
422	fHTime.SetBinError(n+1, fit.GetEventsExcess()*0.1);
423
424	fLog << all << fTimeEffOn << ": " << fit.GetEventsExcess() << endl;
425
426	rebin = steps;
427	}
428
429	// --------------------------------------------------------------------------
430	//
431	// Fill the histogram. For details see the code or the class description
432	//
433	Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
434	{
435	Double_t alpha, energy, theta;
436	Double_t size=-1;
437
438	if (fMatrix)
439	{
440	alpha = (*fMatrix)[fMap[0]];
441	energy = fMap[1]<0 ? -1 : (*fMatrix)[fMap[1]];
442	size = fMap[2]<0 ? -1 : (*fMatrix)[fMap[2]];
443	//<0 ? 1000 : (*fMatrix)[fMap[1]];
444	theta = 0;
445
446	if (energy<0)
447	energy=size;
448	if (size<0)
449	size=energy;
450
451	if (energy<0 && size<0)
452	energy = size = 1000;
453	}
454	else
455	{
456	const MHillasSrc hil = dynamic_cast<const MHillasSrc>(par);
457	if (!par)
458	{
459	*fLog << err << dbginf << "MHillasSrc not found... abort." << endl;
460	return kFALSE;
461	}
462
463	alpha = hil->GetAlpha();
464	if (fHillas)
465	size = fHillas->GetSize();
466	energy = fEnergy ? fEnergy->GetEnergy() : (fHillas?fHillas->GetSize():1000);
467	theta = fPointPos ? fPointPos->GetZd() : 0;
468	}
469
470	// enhance histogram if necessary
471	UpdateAlphaTime();
472
473	// Fill histograms
474	fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w);
475
476	// Check cuts
477	/*
478	if ( (fEnergyMin>=0 && energy<fEnergyMin) \|\|
479	(fEnergyMax>=0 && energy>fEnergyMax) \|\|
480	(fSizeMin>=0 && size <fSizeMin) \|\|
481	(fSizeMax>=0 && size >fSizeMin) )
482	return kTRUE;
483	*/
484
485	if (!fSkipHistTime)
486	fHAlphaTime.Fill(TMath::Abs(alpha), w);
487
488	return kTRUE;
489	}
490
491	// --------------------------------------------------------------------------
492	//
493	// Paint the integral and the error on top of the histogram
494	//
495	void MHAlpha::PaintText(Double_t val, Double_t error) const
496	{
497	TLatex text(0.45, 0.94, Form("N_{exc} = %.1f \\pm %.1f", val, error));
498	text.SetBit(TLatex::kTextNDC);
499	text.SetTextSize(0.04);
500	text.Paint();
501	}
502
503	// --------------------------------------------------------------------------
504	//
505	// Update the projections
506	//
507	void MHAlpha::Paint(Option_t *opt)
508	{
509	// Note: Any object cannot be added to a pad twice!
510	// The object is searched by gROOT->FindObject only in
511	// gPad itself!
512	TVirtualPad *padsave = gPad;
513
514	TH1D *h0=0;
515
516	TString o(opt);
517	if (o==(TString)"proj")
518	{
519	TPaveStats *st=0;
520	for (int x=0; x<4; x++)
521	{
522	TVirtualPad *p=padsave->GetPad(x+1);
523	if (!p \|\| !(st = (TPaveStats*)p->GetPrimitive("stats")))
524	continue;
525
526	if (st->GetOptStat()==11)
527	continue;
528
529	const Double_t y1 = st->GetY1NDC();
530	const Double_t y2 = st->GetY2NDC();
531	const Double_t x1 = st->GetX1NDC();
532	const Double_t x2 = st->GetX2NDC();
533
534	st->SetY1NDC((y2-y1)/3+y1);
535	st->SetX1NDC((x2-x1)/3+x1);
536	st->SetOptStat(11);
537	}
538
539	padsave->cd(1);
540	TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha");
541	if (fOffData)
542	{
543	TH1D *hoff = fOffData->ProjectionZ("ProjAlphaOff");
544	fFit.Scale(hoff, hon);
545
546	hon->SetMaximum();
547	hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05);
548
549	if ((h0=(TH1D*)gPad->FindObject("ProjAlphaOnOff")))
550	{
551	h0->Reset();
552	h0->Add(hoff, hon, -1);
553	const Float_t min = h0->GetMinimum()*1.05;
554	hon->SetMinimum(min<0 ? min : 0);
555	}
556	}
557	FitEnergyBins();
558	FitThetaBins();
559	}
560
561	if (o==(TString)"alpha")
562	if ((h0 = (TH1D*)gPad->FindObject("ProjAlpha")))
563	{
564	// Do not store the 'final' result in fFit
565	MAlphaFitter fit(fFit);
566	TH1D hoff = (TH1D)gPad->FindObject("ProjAlphaOff");
567	fit.Fit(*h0, hoff, kTRUE);
568	fit.PaintResult();
569	}
570
571	if (o==(TString)"time")
572	PaintText(fHTime.Integral(), 0);
573
574	if (o==(TString)"theta")
575	{
576	TH1 h = (TH1)gPad->FindObject("Alpha_x");
577	if (h)
578	{
579	TH1D h2 = (TH1D)fHAlpha.Project3D("dum_x");
580	h2->SetDirectory(0);
581	h2->Scale(fHTheta.Integral()/h2->Integral());
582	h->Reset();
583	h->Add(h2);
584	delete h2;
585	}
586	PaintText(fHTheta.Integral(), 0);
587	}
588
589	if (o==(TString)"energy")
590	{
591	TH1 h = (TH1)gPad->FindObject("Alpha_y");
592	if (h)
593	{
594	TH1D h2 = (TH1D)fHAlpha.Project3D("dum_y");
595	h2->SetDirectory(0);
596	h2->Scale(fHEnergy.Integral()/h2->Integral());
597	h->Reset();
598	h->Add(h2);
599	delete h2;
600	}
601
602	if (fHEnergy.GetMaximum()>1)
603	{
604	gPad->SetLogx();
605	gPad->SetLogy();
606	}
607	FitEnergySpec(kTRUE);
608	}
609
610	gPad = padsave;
611	}
612
613	// --------------------------------------------------------------------------
614	//
615	// Draw the histogram
616	//
617	void MHAlpha::Draw(Option_t *opt)
618	{
619	TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
620
621	// Do the projection before painting the histograms into
622	// the individual pads
623	AppendPad("proj");
624
625	pad->SetBorderMode(0);
626	pad->Divide(2,2);
627
628	TH1D *h=0;
629
630	pad->cd(1);
631	gPad->SetBorderMode(0);
632
633	h = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, -1, 9999, "E");
634	h->SetBit(TH1::kNoTitle);
635	h->SetStats(kTRUE);
636	h->SetXTitle("\\alpha [\\circ]");
637	h->SetYTitle("Counts");
638	h->SetDirectory(NULL);
639	h->SetMarkerStyle(kFullDotMedium);
640	h->SetBit(kCanDelete);
641	h->Draw("");
642
643	if (fOffData)
644	{
645	h->SetMarkerColor(kGreen);
646
647	h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E");
648	h->SetBit(TH1::kNoTitle\|TH1::kNoStats);
649	h->SetXTitle("\\alpha [\\circ]");
650	h->SetYTitle("Counts");
651	h->SetDirectory(NULL);
652	h->SetMarkerStyle(kFullDotMedium);
653	h->SetBit(kCanDelete);
654	h->SetMarkerColor(kRed);
655	h->Draw("same");
656
657	h = (TH1D*)h->Clone("ProjAlphaOnOff");
658	h->SetBit(TH1::kNoTitle);
659	h->SetXTitle("\\alpha [\\circ]");
660	h->SetYTitle("Counts");
661	h->SetDirectory(NULL);
662	h->SetMarkerStyle(kFullDotMedium);
663	h->SetBit(kCanDelete);
664	h->SetMarkerColor(kBlue);
665	h->Draw("same");
666
667	TLine lin;
668	lin.SetLineStyle(kDashed);
669	lin.DrawLine(0, 0, 90, 0);
670	}
671
672	// After the Alpha-projection has been drawn. Fit the histogram
673	// and paint the result into this pad
674	AppendPad("alpha");
675
676	if (fHEnergy.GetNbinsX()>1)
677	{
678	pad->cd(2);
679	gPad->SetBorderMode(0);
680	fHEnergy.Draw();
681
682	AppendPad("energy");
683
684	h = (TH1D*)fHAlpha.Project3D("y");
685	h->SetBit(TH1::kNoTitle\|TH1::kNoStats);
686	h->SetXTitle("E [GeV]");
687	h->SetYTitle("Counts");
688	h->SetDirectory(NULL);
689	h->SetMarkerStyle(kFullDotMedium);
690	h->SetBit(kCanDelete);
691	h->SetMarkerColor(kCyan);
692	h->SetLineColor(kCyan);
693	h->Draw("Psame");
694	}
695	else
696	delete pad->GetPad(2);
697
698	if (fTimeEffOn && fTime \|\| fHTime.GetNbinsX()>1)
699	{
700	pad->cd(3);
701	gPad->SetBorderMode(0);
702	fHTime.Draw();
703	AppendPad("time");
704	}
705	else
706	delete pad->GetPad(3);
707
708	if (fHTheta.GetNbinsX()>1)
709	{
710	pad->cd(4);
711	gPad->SetBorderMode(0);
712	fHTheta.Draw();
713
714	AppendPad("theta");
715
716	h = (TH1D*)fHAlpha.Project3D("x");
717	h->SetBit(TH1::kNoTitle\|TH1::kNoStats);
718	h->SetXTitle("\\theta [\\circ]");
719	h->SetYTitle("Counts");
720	h->SetDirectory(NULL);
721	h->SetMarkerStyle(kFullDotMedium);
722	h->SetBit(kCanDelete);
723	h->SetMarkerColor(kCyan);
724	h->SetLineColor(kCyan);
725	h->Draw("Psame");
726	}
727	else
728	delete pad->GetPad(4);
729	}
730
731	void MHAlpha::DrawAll()
732	{
733	// FIXME: Do in Paint if special option given!
734	TCanvas *c = new TCanvas;
735	Int_t n = fHAlpha.GetNbinsY();
736	Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
737	c->Divide(nc, nc, 0, 0);
738
739	// Do not store the 'final' result in fFit
740	MAlphaFitter fit(fFit);
741
742	for (int i=1; i<=fHAlpha.GetNbinsY(); i++)
743	{
744	c->cd(i);
745
746	TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, i, i, "E");
747	hon->SetBit(TH1::kNoTitle);
748	hon->SetStats(kTRUE);
749	hon->SetXTitle("\\alpha [\\circ]");
750	hon->SetYTitle("Counts");
751	hon->SetDirectory(NULL);
752	hon->SetMarkerStyle(kFullDotMedium);
753	hon->SetBit(kCanDelete);
754	hon->Draw("");
755
756	TH1D *hof = 0;
757
758	if (fOffData)
759	{
760	hon->SetMarkerColor(kGreen);
761
762	hof = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, i, i, "E");
763	hof->SetBit(TH1::kNoTitle\|TH1::kNoStats);
764	hof->SetXTitle("\\alpha [\\circ]");
765	hof->SetYTitle("Counts");
766	hof->SetDirectory(NULL);
767	hof->SetMarkerStyle(kFullDotMedium);
768	hof->SetBit(kCanDelete);
769	hof->SetMarkerColor(kRed);
770	hof->Draw("same");
771
772	fit.Scale(hof, hon);
773
774	hon->SetMaximum();
775	hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05);
776
777	TH1D diff = new TH1D(hon);
778	diff->Add(hof, -1);
779	diff->SetBit(TH1::kNoTitle);
780	diff->SetXTitle("\\alpha [\\circ]");
781	diff->SetYTitle("Counts");
782	diff->SetDirectory(NULL);
783	diff->SetMarkerStyle(kFullDotMedium);
784	diff->SetBit(kCanDelete);
785	diff->SetMarkerColor(kBlue);
786	diff->Draw("same");
787
788	TLine lin;
789	lin.SetLineStyle(kDashed);
790	lin.DrawLine(0, 0, 90, 0);
791
792	const Float_t min = diff->GetMinimum()*1.05;
793	hon->SetMinimum(min<0 ? min : 0);
794	}
795
796	if (hof ? fit.Fit(hon, hof) : fit.Fit(*hon))
797	*fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl;
798	/*
799	if (fit.FitEnergy(fHAlpha, fOffData, i, kTRUE))
800	{
801	fHEnergy.SetBinContent(i, fit.GetEventsExcess());
802	fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
803	}*/
804	}
805
806	}
807
808
809	Bool_t MHAlpha::Finalize()
810	{
811	//TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
812	//h->SetDirectory(0);
813	//Bool_t rc = fFit.Fit(*h);
814	//delete h;
815	if (!fFit.FitAlpha(fHAlpha, fOffData))
816	{
817	*fLog << warn << "Histogram empty." << endl;
818	return kTRUE;
819	}
820
821	// Store the final result in fFit
822	fFit.Print("result");
823
824	if (fResult)
825	fResult->SetVal(fFit.GetSignificance());
826
827	if (!fSkipHistEnergy)
828	{
829	*fLog << inf << "Processing energy bins..." << endl;
830	FitEnergyBins();
831	}
832	if (!fSkipHistTheta)
833	{
834	*fLog << inf << "Processing theta bins..." << endl;
835	FitThetaBins();
836	}
837	if (!fSkipHistTime)
838	{
839	*fLog << inf << "Processing time bins..." << endl;
840	UpdateAlphaTime(kTRUE);
841	MH::RemoveFirstBin(fHTime);
842	}
843
844	return kTRUE;
845	}
846
847	// --------------------------------------------------------------------------
848	//
849	// You can use this function if you want to use a MHMatrix instead of
850	// MMcEvt. This function adds all necessary columns to the
851	// given matrix. Afterward you should fill the matrix with the corresponding
852	// data (eg from a file by using MHMatrix::Fill). If you now loop
853	// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
854	// will take the values from the matrix instead of the containers.
855	//
856	// It takes fSkipHist* into account!
857	//
858	void MHAlpha::InitMapping(MHMatrix *mat, Int_t type)
859	{
860	if (fMatrix)
861	return;
862
863	fMatrix = mat;
864
865	fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha");
866	fMap[1] = -1;
867	fMap[2] = -1;
868	fMap[3] = -1;
869	fMap[4] = -1;
870
871	if (!fSkipHistEnergy)
872	if (type==0)
873	{
874	fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy");
875	fMap[2] = -1;
876	}
877	else
878	{
879	fMap[1] = -1;
880	fMap[2] = fMatrix->AddColumn("MHillas.fSize");
881	}
882
883	if (!fSkipHistTheta)
884	fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
885
886	// if (!fSkipHistTime)
887	// fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
888	}
889
890	void MHAlpha::StopMapping()
891	{
892	fMatrix = NULL;
893	}

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

Download in other formats: