Context Navigation

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

Visit:

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

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

Download in other formats: