Context Navigation

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

Visit:

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

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

Download in other formats: