Context Navigation

source: tags/Mars-V0.9.2/mhflux/MHAlpha.cc

Visit:

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

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

Download in other formats: