Context Navigation

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

Visit:

Last change on this file since 5814 was 5807, checked in by tbretz, 20 years ago
* empty log message *
File size: 23.6 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	fHEnergy.SetEntries(0);
188
189	for (int i=1; i<=n; i++)
190	{
191	if (fit.FitEnergy(fHAlpha, fOffData, i))
192	{
193	fHEnergy.SetBinContent(i, fit.GetEventsExcess());
194	fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
195	}
196	}
197	}
198
199	void MHAlpha::FitThetaBins(Bool_t paint)
200	{
201	// Do not store the 'final' result in fFit
202	MAlphaFitter fit(fFit);
203
204	const Int_t n = fHAlpha.GetNbinsX();
205
206	fHTheta.SetEntries(0);
207
208	for (int i=1; i<=n; i++)
209	{
210	if (fit.FitTheta(fHAlpha, fOffData, i))
211	{
212	fHTheta.SetBinContent(i, fit.GetEventsExcess());
213	fHTheta.SetBinError(i, fit.GetEventsExcess()*0.2);
214	}
215	}
216	}
217
218	Bool_t MHAlpha::SetupFill(const MParList *pl)
219	{
220	fHAlpha.Reset();
221	fHEnergy.Reset();
222	fHTheta.Reset();
223
224	if (fName!=(TString)"MHAlphaOff" && fOffData==NULL)
225	{
226	MHAlpha hoff = (MHAlpha)pl->FindObject("MHAlphaOff", "MHAlpha");
227	if (!hoff)
228	*fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl;
229	else
230	{
231	*fLog << inf << "MHAlphaOff [MHAlpha] found... using on-off mode!" << endl;
232	SetOffData(*hoff);
233	}
234	}
235
236	fHillas = 0;
237	/*
238	if (fSizeMin>=0 \|\| fSizeMax>=0)
239	{
240	fHillas = (MHillas*)pl->FindObject("MHillas");
241	if (!fHillas)
242	{
243	*fLog << warn << "Size cut set, but MHillas not found... abort." << endl;
244	return kFALSE;
245	}
246	}
247	*/
248	fEnergy = (MEnergyEst*)pl->FindObject("MEnergyEst");
249	if (!fEnergy)
250	{ /*
251	if (fEnergyMin>=0 \|\| fEnergyMax>=0)
252	{
253	*fLog << warn << "Energy cut set, but MEnergyEst not found... abort." << endl;
254	return kFALSE;
255	} */
256
257	*fLog << warn << "MEnergyEst not found... " << flush;
258
259	if (!fHillas)
260	fHillas = (MHillas*)pl->FindObject("MHillas");
261	if (fHillas)
262	*fLog << "using SIZE instead!" << endl;
263	else
264	*fLog << "ignored." << endl;
265
266	fHEnergy.SetName("Size");
267	fHEnergy.SetTitle(" N_{exc} vs. Size ");
268	fHEnergy.SetXTitle("Size [\\gamma]");
269	}
270	else
271	{
272	fHEnergy.SetName("Energy");
273	fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
274	fHEnergy.SetXTitle("E_{est} [GeV]");
275	}
276
277	fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
278	if (!fPointPos)
279	*fLog << warn << "MPointingPos not found... ignored." << endl;
280
281	fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime");
282	if (!fTimeEffOn)
283	*fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;
284	else
285	*fTimeEffOn = MTime(); // FIXME: How to do it different?
286
287	fTime = (MTime*)pl->FindObject("MTime");
288	if (!fTime)
289	*fLog << warn << "MTime not found... ignored." << endl;
290
291	fResult = (MParameterD)const_cast<MParList>(pl)->FindCreateObj("MParameterD", "MinimizationValue");
292	if (!fResult)
293	return kFALSE;
294
295	//fExcess = (MParameterD)const_cast<MParList>(pl)->FindCreateObj("MParameterD", "MExcess");
296	//if (!fExcess)
297	// return kFALSE;
298
299	fLastTime = MTime();
300
301	MBinning binst, binse, binsa;
302	binst.SetEdges(fOffData ? *fOffData : fHAlpha, 'x');
303	binse.SetEdges(fOffData ? *fOffData : fHAlpha, 'y');
304	binsa.SetEdges(fOffData ? *fOffData : fHAlpha, 'z');
305
306	if (!fOffData)
307	{
308	if (!fPointPos)
309	binst.SetEdges(1, 0, 60);
310	else
311	binst.SetEdges(*pl, "BinningTheta");
312
313	if (!fEnergy && !fHillas)
314	binse.SetEdges(1, 10, 100000);
315	else
316	binse.SetEdges(*pl, "BinningEnergyEst");
317
318	binsa.SetEdges(*pl, "BinningAlpha");
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	fHAlphaTime.SetEntries(0);
403
404	// Get number of bins
405	const Int_t n = fHTime.GetNbinsX();
406
407	//
408	// Prepare Histogram
409	//
410	if (final)
411	time->Plus1ns();
412
413	// Enhance binning
414	MBinning bins;
415	bins.SetEdges(fHTime, 'x');
416	bins.AddEdge(time->GetAxisTime());
417	bins.Apply(fHTime);
418
419	//
420	// Fill histogram
421	//
422	fHTime.SetBinContent(n+1, fit.GetEventsExcess());
423	fHTime.SetBinError(n+1, fit.GetEventsExcess()*0.1);
424
425	fLog << all << fTimeEffOn << ": " << fit.GetEventsExcess() << endl;
426
427	rebin = steps;
428	}
429
430	// --------------------------------------------------------------------------
431	//
432	// Fill the histogram. For details see the code or the class description
433	//
434	Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
435	{
436	Double_t alpha, energy, theta;
437	Double_t size=-1;
438
439	if (fMatrix)
440	{
441	alpha = (*fMatrix)[fMap[0]];
442	energy = fMap[1]<0 ? -1 : (*fMatrix)[fMap[1]];
443	size = fMap[2]<0 ? -1 : (*fMatrix)[fMap[2]];
444	//<0 ? 1000 : (*fMatrix)[fMap[1]];
445	theta = 0;
446
447	if (energy<0)
448	energy=size;
449	if (size<0)
450	size=energy;
451
452	if (energy<0 && size<0)
453	energy = size = 1000;
454	}
455	else
456	{
457	const MHillasSrc hil = dynamic_cast<const MHillasSrc>(par);
458	if (!par)
459	{
460	*fLog << err << dbginf << "MHillasSrc not found... abort." << endl;
461	return kFALSE;
462	}
463
464	alpha = hil->GetAlpha();
465	if (fHillas)
466	size = fHillas->GetSize();
467	energy = fEnergy ? fEnergy->GetEnergy() : (fHillas?fHillas->GetSize():1000);
468	theta = fPointPos ? fPointPos->GetZd() : 0;
469	}
470
471	// enhance histogram if necessary
472	UpdateAlphaTime();
473
474	// Fill histograms
475	fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w);
476
477	// Check cuts
478	/*
479	if ( (fEnergyMin>=0 && energy<fEnergyMin) \|\|
480	(fEnergyMax>=0 && energy>fEnergyMax) \|\|
481	(fSizeMin>=0 && size <fSizeMin) \|\|
482	(fSizeMax>=0 && size >fSizeMin) )
483	return kTRUE;
484	*/
485
486	if (!fSkipHistTime)
487	fHAlphaTime.Fill(TMath::Abs(alpha), w);
488
489	return kTRUE;
490	}
491
492	// --------------------------------------------------------------------------
493	//
494	// Paint the integral and the error on top of the histogram
495	//
496	void MHAlpha::PaintText(Double_t val, Double_t error) const
497	{
498	TLatex text(0.45, 0.94, Form("N_{exc} = %.1f \\pm %.1f", val, error));
499	text.SetBit(TLatex::kTextNDC);
500	text.SetTextSize(0.04);
501	text.Paint();
502	}
503
504	// --------------------------------------------------------------------------
505	//
506	// Update the projections
507	//
508	void MHAlpha::Paint(Option_t *opt)
509	{
510	// Note: Any object cannot be added to a pad twice!
511	// The object is searched by gROOT->FindObject only in
512	// gPad itself!
513	TVirtualPad *padsave = gPad;
514
515	TH1D *h0=0;
516
517	TString o(opt);
518	if (o==(TString)"proj")
519	{
520	TPaveStats *st=0;
521	for (int x=0; x<4; x++)
522	{
523	TVirtualPad *p=padsave->GetPad(x+1);
524	if (!p \|\| !(st = (TPaveStats*)p->GetPrimitive("stats")))
525	continue;
526
527	if (st->GetOptStat()==11)
528	continue;
529
530	const Double_t y1 = st->GetY1NDC();
531	const Double_t y2 = st->GetY2NDC();
532	const Double_t x1 = st->GetX1NDC();
533	const Double_t x2 = st->GetX2NDC();
534
535	st->SetY1NDC((y2-y1)/3+y1);
536	st->SetX1NDC((x2-x1)/3+x1);
537	st->SetOptStat(11);
538	}
539
540	padsave->cd(1);
541	TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha");
542	if (fOffData)
543	{
544	TH1D *hoff = fOffData->ProjectionZ("ProjAlphaOff");
545	fFit.Scale(hoff, hon);
546
547	hon->SetMaximum();
548	hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05);
549
550	if ((h0=(TH1D*)gPad->FindObject("ProjAlphaOnOff")))
551	{
552	h0->Reset();
553	h0->Add(hoff, hon, -1);
554	const Float_t min = h0->GetMinimum()*1.05;
555	hon->SetMinimum(min<0 ? min : 0);
556	}
557	}
558	FitEnergyBins();
559	FitThetaBins();
560	}
561
562	if (o==(TString)"alpha")
563	if ((h0 = (TH1D*)gPad->FindObject("ProjAlpha")))
564	{
565	// Do not store the 'final' result in fFit
566	MAlphaFitter fit(fFit);
567	TH1D hoff = (TH1D)gPad->FindObject("ProjAlphaOff");
568	fit.Fit(*h0, hoff, kTRUE);
569	fit.PaintResult();
570	}
571
572	if (o==(TString)"time")
573	PaintText(fHTime.Integral(), 0);
574
575	if (o==(TString)"theta")
576	{
577	TH1 h = (TH1)gPad->FindObject("Alpha_x");
578	if (h)
579	{
580	TH1D h2 = (TH1D)fHAlpha.Project3D("dum_x");
581	h2->SetDirectory(0);
582	h2->Scale(fHTheta.Integral()/h2->Integral());
583	h->Reset();
584	h->Add(h2);
585	delete h2;
586	}
587	PaintText(fHTheta.Integral(), 0);
588	}
589
590	if (o==(TString)"energy")
591	{
592	TH1 h = (TH1)gPad->FindObject("Alpha_y");
593	if (h)
594	{
595	TH1D h2 = (TH1D)fHAlpha.Project3D("dum_y");
596	h2->SetDirectory(0);
597	h2->Scale(fHEnergy.Integral()/h2->Integral());
598	h->Reset();
599	h->Add(h2);
600	delete h2;
601	}
602
603	if (fHEnergy.GetMaximum()>1)
604	{
605	gPad->SetLogx();
606	gPad->SetLogy();
607	}
608	FitEnergySpec(kTRUE);
609	}
610
611	gPad = padsave;
612	}
613
614	// --------------------------------------------------------------------------
615	//
616	// Draw the histogram
617	//
618	void MHAlpha::Draw(Option_t *opt)
619	{
620	TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
621
622	// Do the projection before painting the histograms into
623	// the individual pads
624	AppendPad("proj");
625
626	pad->SetBorderMode(0);
627	pad->Divide(2,2);
628
629	TH1D *h=0;
630
631	pad->cd(1);
632	gPad->SetBorderMode(0);
633
634	h = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, -1, 9999, "E");
635	h->SetBit(TH1::kNoTitle);
636	h->SetStats(kTRUE);
637	h->SetXTitle("\\alpha [\\circ]");
638	h->SetYTitle("Counts");
639	h->SetDirectory(NULL);
640	h->SetMarkerStyle(kFullDotMedium);
641	h->SetBit(kCanDelete);
642	h->Draw("");
643
644	if (fOffData)
645	{
646	h->SetMarkerColor(kGreen);
647
648	h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E");
649	h->SetBit(TH1::kNoTitle);
650	h->SetXTitle("\\alpha [\\circ]");
651	h->SetYTitle("Counts");
652	h->SetDirectory(NULL);
653	h->SetMarkerStyle(kFullDotMedium);
654	h->SetBit(kCanDelete);
655	h->SetMarkerColor(kRed);
656	h->Draw("same");
657
658	h = (TH1D*)h->Clone("ProjAlphaOnOff");
659	h->SetBit(TH1::kNoTitle);
660	h->SetXTitle("\\alpha [\\circ]");
661	h->SetYTitle("Counts");
662	h->SetDirectory(NULL);
663	h->SetMarkerStyle(kFullDotMedium);
664	h->SetBit(kCanDelete);
665	h->SetMarkerColor(kBlue);
666	h->Draw("same");
667
668	TLine lin;
669	lin.SetLineStyle(kDashed);
670	lin.DrawLine(0, 0, 90, 0);
671	}
672
673	// After the Alpha-projection has been drawn. Fit the histogram
674	// and paint the result into this pad
675	AppendPad("alpha");
676
677	if (fHEnergy.GetNbinsX()>1)
678	{
679	pad->cd(2);
680	gPad->SetBorderMode(0);
681	fHEnergy.Draw();
682
683	AppendPad("energy");
684
685	h = (TH1D*)fHAlpha.Project3D("y");
686	h->SetBit(TH1::kNoTitle\|TH1::kNoStats);
687	h->SetXTitle("E [GeV]");
688	h->SetYTitle("Counts");
689	h->SetDirectory(NULL);
690	h->SetMarkerStyle(kFullDotMedium);
691	h->SetBit(kCanDelete);
692	h->SetMarkerColor(kCyan);
693	h->SetLineColor(kCyan);
694	h->Draw("Psame");
695	}
696	else
697	delete pad->GetPad(2);
698
699	if (fTimeEffOn && fTime \|\| fHTime.GetNbinsX()>1)
700	{
701	pad->cd(3);
702	gPad->SetBorderMode(0);
703	fHTime.Draw();
704	AppendPad("time");
705	}
706	else
707	delete pad->GetPad(3);
708
709	if (fHTheta.GetNbinsX()>1)
710	{
711	pad->cd(4);
712	gPad->SetBorderMode(0);
713	fHTheta.Draw();
714
715	AppendPad("theta");
716
717	h = (TH1D*)fHAlpha.Project3D("x");
718	h->SetBit(TH1::kNoTitle\|TH1::kNoStats);
719	h->SetXTitle("\\theta [\\circ]");
720	h->SetYTitle("Counts");
721	h->SetDirectory(NULL);
722	h->SetMarkerStyle(kFullDotMedium);
723	h->SetBit(kCanDelete);
724	h->SetMarkerColor(kCyan);
725	h->SetLineColor(kCyan);
726	h->Draw("Psame");
727	}
728	else
729	delete pad->GetPad(4);
730	}
731
732	void MHAlpha::DrawAll()
733	{
734	// FIXME: Do in Paint if special option given!
735	TCanvas *c = new TCanvas;
736	Int_t n = fHAlpha.GetNbinsY();
737	Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
738	c->Divide(nc, nc, 0, 0);
739
740	// Do not store the 'final' result in fFit
741	MAlphaFitter fit(fFit);
742
743	for (int i=1; i<=fHAlpha.GetNbinsY(); i++)
744	{
745	c->cd(i);
746
747	TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, i, i, "E");
748	hon->SetBit(TH1::kNoTitle);
749	hon->SetStats(kTRUE);
750	hon->SetXTitle("\\alpha [\\circ]");
751	hon->SetYTitle("Counts");
752	hon->SetDirectory(NULL);
753	hon->SetMarkerStyle(kFullDotMedium);
754	hon->SetBit(kCanDelete);
755	hon->Draw("");
756
757	TH1D *hof = 0;
758
759	if (fOffData)
760	{
761	hon->SetMarkerColor(kGreen);
762
763	hof = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, i, i, "E");
764	hof->SetBit(TH1::kNoTitle\|TH1::kNoStats);
765	hof->SetXTitle("\\alpha [\\circ]");
766	hof->SetYTitle("Counts");
767	hof->SetDirectory(NULL);
768	hof->SetMarkerStyle(kFullDotMedium);
769	hof->SetBit(kCanDelete);
770	hof->SetMarkerColor(kRed);
771	hof->Draw("same");
772
773	fit.Scale(hof, hon);
774
775	hon->SetMaximum();
776	hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05);
777
778	TH1D diff = new TH1D(hon);
779	diff->Add(hof, -1);
780	diff->SetBit(TH1::kNoTitle);
781	diff->SetXTitle("\\alpha [\\circ]");
782	diff->SetYTitle("Counts");
783	diff->SetDirectory(NULL);
784	diff->SetMarkerStyle(kFullDotMedium);
785	diff->SetBit(kCanDelete);
786	diff->SetMarkerColor(kBlue);
787	diff->Draw("same");
788
789	TLine lin;
790	lin.SetLineStyle(kDashed);
791	lin.DrawLine(0, 0, 90, 0);
792
793	const Float_t min = diff->GetMinimum()*1.05;
794	hon->SetMinimum(min<0 ? min : 0);
795	}
796
797	if (hof ? fit.Fit(hon, hof) : fit.Fit(*hon))
798	*fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl;
799	/*
800	if (fit.FitEnergy(fHAlpha, fOffData, i, kTRUE))
801	{
802	fHEnergy.SetBinContent(i, fit.GetEventsExcess());
803	fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
804	}*/
805	}
806
807	}
808
809
810	Bool_t MHAlpha::Finalize()
811	{
812	//TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
813	//h->SetDirectory(0);
814	//Bool_t rc = fFit.Fit(*h);
815	//delete h;
816	if (!fFit.FitAlpha(fHAlpha, fOffData))
817	{
818	*fLog << warn << "Histogram empty." << endl;
819	return kTRUE;
820	}
821
822	// Store the final result in fFit
823	fFit.Print("result");
824
825	if (fResult)
826	fResult->SetVal(-fFit.GetSignificance());
827
828	if (!fSkipHistEnergy)
829	{
830	*fLog << inf << "Processing energy bins..." << endl;
831	FitEnergyBins();
832	}
833	if (!fSkipHistTheta)
834	{
835	*fLog << inf << "Processing theta bins..." << endl;
836	FitThetaBins();
837	}
838	if (!fSkipHistTime)
839	{
840	*fLog << inf << "Processing time bins..." << endl;
841	UpdateAlphaTime(kTRUE);
842	MH::RemoveFirstBin(fHTime);
843	}
844
845	return kTRUE;
846	}
847
848	// --------------------------------------------------------------------------
849	//
850	// You can use this function if you want to use a MHMatrix instead of
851	// MMcEvt. This function adds all necessary columns to the
852	// given matrix. Afterward you should fill the matrix with the corresponding
853	// data (eg from a file by using MHMatrix::Fill). If you now loop
854	// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
855	// will take the values from the matrix instead of the containers.
856	//
857	// It takes fSkipHist* into account!
858	//
859	void MHAlpha::InitMapping(MHMatrix *mat, Int_t type)
860	{
861	if (fMatrix)
862	return;
863
864	fMatrix = mat;
865
866	fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha");
867	fMap[1] = -1;
868	fMap[2] = -1;
869	fMap[3] = -1;
870	fMap[4] = -1;
871
872	if (!fSkipHistEnergy)
873	if (type==0)
874	{
875	fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy");
876	fMap[2] = -1;
877	}
878	else
879	{
880	fMap[1] = -1;
881	fMap[2] = fMatrix->AddColumn("MHillas.fSize");
882	}
883
884	if (!fSkipHistTheta)
885	fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
886
887	// if (!fSkipHistTime)
888	// fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
889	}
890
891	void MHAlpha::StopMapping()
892	{
893	fMatrix = NULL;
894	}

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

Download in other formats: