Context Navigation

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

Visit:

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

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

Download in other formats: