Context Navigation

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

Visit:

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

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

Download in other formats: