Context Navigation

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

Visit:

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

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

Download in other formats: