Context Navigation

source: branches/Mars_MC/mhflux/MHAlpha.cc@ 17063

Visit:

Last change on this file since 17063 was 9851, checked in by tbretz, 14 years ago
Changed MH::SetBinning and similar functions to take references instead of pointers as arguments. For convenience wrappers for the old style functions are provided.
File size: 31.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 <TStyle.h>
49	#include <TLatex.h>
50	#include <TCanvas.h>
51	#include <TPaveStats.h>
52
53	#include "MSrcPosCam.h"
54	#include "MHillas.h"
55	#include "MHillasSrc.h"
56	#include "MTime.h"
57	#include "MObservatory.h"
58	#include "MPointingPos.h"
59	#include "MAstroSky2Local.h"
60	#include "MStatusDisplay.h"
61	#include "MParameters.h"
62	#include "MHMatrix.h"
63
64	#include "MString.h"
65	#include "MBinning.h"
66	#include "MParList.h"
67
68	#include "MLog.h"
69	#include "MLogManip.h"
70
71	ClassImp(MHAlpha);
72
73	using namespace std;
74
75	// --------------------------------------------------------------------------
76	//
77	// Default Constructor
78	//
79	MHAlpha::MHAlpha(const char name, const char title)
80	: fNameParameter("MHillasSrc"), fParameter(0),
81	fOffData(0), fResult(0), fSigma(0), fEnergy(0), fBin(0),
82	fPointPos(0), fTimeEffOn(0), fTime(0), fNumTimeBins(10),
83	fHillas(0), fMatrix(0), fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE),
84	fSkipHistEnergy(kFALSE), fForceUsingSize(kFALSE)
85	{
86	//
87	// set the name and title of this object
88	//
89	fName = name ? name : "MHAlpha";
90	fTitle = title ? title : "Alpha plot";
91
92	fHist.SetName("Alpha");
93	fHist.SetTitle("Alpha");
94	fHist.SetXTitle("\\Theta [deg]");
95	//fHist.SetYTitle("E_{est} [GeV]");
96	fHist.SetZTitle("\|\\alpha\| [\\circ]");
97	fHist.SetDirectory(NULL);
98	fHist.UseCurrentStyle();
99
100	// Main histogram
101	fHistTime.SetName("Alpha");
102	fHistTime.SetXTitle("\|\\alpha\| [\\circ]");
103	fHistTime.SetYTitle("Counts");
104	fHistTime.UseCurrentStyle();
105	fHistTime.SetDirectory(NULL);
106
107	fHEnergy.SetName("Excess");
108	//fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
109	//fHEnergy.SetXTitle("E_{est} [GeV]");
110	fHEnergy.SetYTitle("N_{exc}");
111	fHEnergy.SetDirectory(NULL);
112	fHEnergy.UseCurrentStyle();
113
114	fHTheta.SetName("ExcessTheta");
115	fHTheta.SetTitle(" N_{exc} vs. \\Theta ");
116	fHTheta.SetXTitle("\\Theta [\\circ]");
117	fHTheta.SetYTitle("N_{exc}");
118	fHTheta.SetDirectory(NULL);
119	fHTheta.UseCurrentStyle();
120	fHTheta.SetMinimum(0);
121
122	// effective on time versus time
123	fHTime.SetName("ExcessTime");
124	fHTime.SetTitle(" N_{exc} vs. Time ");
125	fHTime.SetXTitle("Time");
126	fHTime.SetYTitle("N_{exc} [s]");
127	fHTime.UseCurrentStyle();
128	fHTime.SetDirectory(NULL);
129	fHTime.GetYaxis()->SetTitleOffset(1.2);
130	fHTime.GetXaxis()->SetLabelSize(0.033);
131	fHTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
132	fHTime.GetXaxis()->SetTimeDisplay(1);
133	fHTime.SetMinimum(0);
134	fHTime.Sumw2();
135
136	MBinning binsa, binse, binst;
137	binsa.SetEdges(18, 0, 90);
138	binse.SetEdgesLog(15, 10, 100000);
139	binst.SetEdgesASin(67, -0.005, 0.665);
140	binse.Apply(fHEnergy);
141	binst.Apply(fHTheta);
142	binsa.Apply(fHistTime);
143
144	MH::SetBinning(fHist, binst, binse, binsa);
145	}
146
147	Float_t MHAlpha::FitEnergyBins(Bool_t paint)
148	{
149	// Do not store the 'final' result in fFit
150	MAlphaFitter fit(fFit);
151
152	const Int_t n = fHist.GetNbinsY();
153
154	fHEnergy.SetEntries(0);
155
156	Float_t mean = 0;
157	Int_t num = 0;
158	for (int i=1; i<=n; i++)
159	{
160	if (!fit.FitEnergy(fHist, fOffData, i))
161	continue;
162
163	// FIXME: Calculate UL!
164	if (fit.GetEventsExcess()<=0)
165	continue;
166
167	fHEnergy.SetBinContent(i, fit.GetEventsExcess());
168	fHEnergy.SetBinError(i, fit.GetErrorExcess());
169
170	mean += fit.GetEventsExcess()*fit.GetEventsExcess()/fit.GetErrorExcess()/fit.GetErrorExcess();
171	num++;
172	}
173	return TMath::Sqrt(mean)/num;
174	}
175
176	void MHAlpha::FitThetaBins(Bool_t paint)
177	{
178	// Do not store the 'final' result in fFit
179	MAlphaFitter fit(fFit);
180
181	const Int_t n = fHist.GetNbinsX();
182
183	fHTheta.SetEntries(0);
184
185	for (int i=1; i<=n; i++)
186	{
187	if (!fit.FitTheta(fHist, fOffData, i))
188	continue;
189
190	// FIXME: Calculate UL!
191	if (fit.GetEventsExcess()<=0)
192	continue;
193
194	fHTheta.SetBinContent(i, fit.GetEventsExcess());
195	fHTheta.SetBinError(i, fit.GetErrorExcess());
196	}
197	}
198
199	// --------------------------------------------------------------------------
200	//
201	// Return the value from fParemeter which should be filled into the plots
202	//
203	Double_t MHAlpha::GetVal() const
204	{
205	return static_cast<const MHillasSrc*>(fParameter)->GetAlpha();
206	}
207
208
209	// --------------------------------------------------------------------------
210	//
211	// Store the pointer to the parameter container storing the plotted value
212	// (here MHillasSrc) in fParameter.
213	//
214	// return whether it was found or not.
215	//
216	Bool_t MHAlpha::GetParameter(const MParList &pl)
217	{
218	fParameter = (MParContainer*)pl.FindObject(fNameParameter, "MHillasSrc");
219	if (fParameter)
220	return kTRUE;
221
222	*fLog << err << fNameParameter << " [MHillasSrc] not found... abort." << endl;
223	return kFALSE;
224	}
225
226	Bool_t MHAlpha::SetupFill(const MParList *pl)
227	{
228	fHist.Reset();
229	fHEnergy.Reset();
230	fHTheta.Reset();
231	fHTime.Reset();
232
233	const TString off(MString::Format("%sOff", fName.Data()));
234	if (fName!=off && fOffData==NULL)
235	{
236	const TString desc(MString::Format("%s [%s] found... using ", off.Data(), ClassName()));
237	MHAlpha hoff = (MHAlpha)pl->FindObject(off, ClassName());
238	if (!hoff)
239	*fLog << inf << "No " << desc << "current data only!" << endl;
240	else
241	{
242	*fLog << inf << desc << "on-off mode!" << endl;
243	SetOffData(*hoff);
244	}
245	}
246
247	if (!GetParameter(*pl))
248	return kFALSE;
249
250	fHillas = 0;
251	fEnergy = fForceUsingSize ? 0 : (MParameterD*)pl->FindObject("MEnergyEst", "MParameterD");
252	if (!fEnergy)
253	{
254	*fLog << warn << "MEnergyEst [MParameterD] not found... " << flush;
255
256	if (!fHillas)
257	fHillas = (MHillas*)pl->FindObject("MHillas");
258	if (fHillas)
259	*fLog << "using SIZE instead!" << endl;
260	else
261	*fLog << "ignored." << endl;
262
263	fHEnergy.SetTitle(" N_{exc} vs. Size ");
264	fHEnergy.SetXTitle("Size [phe]");
265	fHist.SetYTitle("Size [phe]");
266	}
267	else
268	{
269	fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
270	fHEnergy.SetXTitle("E_{est} [GeV]");
271	fHist.SetYTitle("E_{est} [GeV]");
272	}
273
274	fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
275	if (!fPointPos)
276	*fLog << warn << "MPointingPos not found... ignored." << endl;
277
278	fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime");
279	if (!fTimeEffOn)
280	*fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;
281	else
282	*fTimeEffOn = MTime(); // FIXME: How to do it different?
283
284	fTime = (MTime*)pl->FindObject("MTime");
285	if (!fTime)
286	*fLog << warn << "MTime not found... ignored." << endl;
287
288	fResult = (MParameterD)const_cast<MParList>(pl)->FindCreateObj("MParameterD", "MinimizationValue");
289	if (!fResult)
290	return kFALSE;
291	fSigma = (MParameterD)const_cast<MParList>(pl)->FindCreateObj("MParameterD", "GaussSigma");
292	if (!fSigma)
293	return kFALSE;
294	fBin = (MParameterI)const_cast<MParList>(pl)->FindCreateObj("MParameterI", "Bin");
295	if (!fBin)
296	return kFALSE;
297
298	//fExcess = (MParameterD)const_cast<MParList>(pl)->FindCreateObj("MParameterD", "MExcess");
299	//if (!fExcess)
300	// return kFALSE;
301
302	fLastTime = MTime();
303	fNumRebin = fNumTimeBins;
304
305	MBinning binst, binse, binsa;
306	binst.SetEdges(fOffData ? *fOffData : fHist, 'x');
307	binse.SetEdges(fOffData ? *fOffData : fHist, 'y');
308	binsa.SetEdges(fOffData ? *fOffData : fHist, 'z');
309	if (!fOffData)
310	{
311	if (!fPointPos)
312	binst.SetEdges(1, 0, 60);
313	else
314	binst.SetEdges(*pl, "BinningTheta");
315
316	if (!fEnergy && !fHillas)
317	binse.SetEdges(1, 10, 100000);
318	else
319	if (fEnergy)
320	binse.SetEdges(*pl, "BinningEnergyEst");
321	else
322	binse.SetEdges(*pl, "BinningSize");
323
324	binsa.SetEdges(*pl, MString::Format("Binning%s", ClassName()+2));
325	}
326	else
327	{
328	fHEnergy.SetTitle(fOffData->GetTitle());
329	fHEnergy.SetXTitle(fOffData->GetYaxis()->GetTitle());
330	fHist.SetYTitle(fOffData->GetYaxis()->GetTitle());
331	}
332
333	binse.Apply(fHEnergy);
334	binst.Apply(fHTheta);
335	binsa.Apply(fHistTime);
336	MH::SetBinning(fHist, binst, binse, binsa);
337
338	MAlphaFitter fit = (MAlphaFitter)pl->FindObject("MAlphaFitter");
339	if (!fit)
340	*fLog << warn << "MAlphaFitter not found... ignored." << endl;
341	else
342	fit->Copy(fFit);
343
344	*fLog << inf;
345	fFit.Print();
346
347	return kTRUE;
348	}
349
350	void MHAlpha::InitAlphaTime(const MTime &t)
351	{
352	//
353	// If this is the first call we have to initialize the time-histogram
354	//
355	const MBinning bins(1, t.GetAxisTime()-60, t.GetAxisTime());
356	bins.Apply(fHTime);
357
358	fLastTime=t;
359	}
360
361	void MHAlpha::UpdateAlphaTime(Bool_t final)
362	{
363	if (!fTimeEffOn)
364	return;
365
366	if (!final)
367	{
368	if (fLastTime==MTime() && *fTimeEffOn!=MTime())
369	{
370	InitAlphaTime(*fTimeEffOn);
371	return;
372	}
373
374	if (fLastTime!=*fTimeEffOn)
375	{
376	fLastTime=*fTimeEffOn;
377	fNumRebin--;
378	}
379
380	if (fNumRebin>0)
381	return;
382	}
383
384	// Check new 'last time'
385	MTime *time = final ? fTime : fTimeEffOn;
386
387	if (time->GetAxisTime()<=fHTime.GetXaxis()->GetXmax())
388	{
389	fLog << warn << "WARNING - New time-stamp " << time << " lower" << endl;
390	*fLog << "than upper edge of histogram... skipped." << endl;
391	*fLog << "This should not happen. Maybe you started you eventloop with" << endl;
392	*fLog << "an already initialized time stamp MTimeEffectiveOnTime?" << endl;
393	fNumRebin++;
394	return;
395	}
396
397	// Fit time histogram
398	MAlphaFitter fit(fFit);
399
400	TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1, "E") : 0;
401	const Bool_t rc = fit.ScaleAndFit(fHistTime, h);
402
403	if (h)
404	delete h;
405
406	if (!rc)
407	return;
408
409	// Reset Histogram
410	fHistTime.Reset();
411	fHistTime.SetEntries(0);
412
413	//
414	// Prepare Histogram
415	//
416	if (final)
417	time->Plus1ns();
418
419	// Enhance binning
420	MBinning bins;
421	bins.SetEdges(fHTime, 'x');
422	bins.AddEdge(time->GetAxisTime());
423	bins.Apply(fHTime);
424
425	//
426	// Fill histogram
427	//
428	// Get number of bins
429	const Int_t n = fHTime.GetNbinsX();
430
431	fHTime.SetBinContent(n, fit.GetEventsExcess());
432	fHTime.SetBinError(n, fit.GetErrorExcess());
433
434	//fLog << inf3 << fTimeEffOn << " (" << n << "): " << fit.GetEventsExcess() << endl;
435
436	fNumRebin = fNumTimeBins;
437	}
438
439	void MHAlpha::SetBin(Int_t ibin)
440	{
441	// Is this necessary?
442	// Could be speed up up searching for it only once.
443	const Float_t max = fFit.GetSignalIntegralMax();
444	const Int_t bin0 = fHist.GetZaxis()->FindFixBin(max);
445
446	const Int_t nbinsx = fHist.GetNbinsX();
447	const Int_t nbinsy = fHist.GetNbinsY();
448	const Int_t nxy = (nbinsx+2)*(nbinsy+2);
449
450	const Int_t binz = ibin/nxy;
451
452	const Bool_t issignal = binz>0 && binz<bin0;
453
454	fBin->SetVal(issignal ? binz : -binz);
455	}
456
457	// --------------------------------------------------------------------------
458	//
459	// Fill the histogram. For details see the code or the class description
460	//
461	Int_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
462	{
463	Double_t alpha, energy, theta;
464	Double_t size=-1;
465
466	if (fMatrix)
467	{
468	alpha = fMap[0]<0 ? GetVal() : (*fMatrix)[fMap[0]];
469	energy = fMap[1]<0 ? -1 : (*fMatrix)[fMap[1]];
470	size = fMap[2]<0 ? -1 : (*fMatrix)[fMap[2]];
471	//<0 ? 1000 : (*fMatrix)[fMap[1]];
472	theta = 0;
473
474	if (energy<0)
475	energy=size;
476	if (size<0)
477	size=energy;
478
479	if (energy<0 && size<0)
480	energy = size = 1000;
481	}
482	else
483	{
484	alpha = GetVal();
485
486	if (fHillas)
487	size = fHillas->GetSize();
488	energy = fEnergy ? fEnergy->GetVal() : (fHillas?fHillas->GetSize():1000);
489	theta = fPointPos ? fPointPos->GetZd() : 0;
490	}
491
492	// enhance histogram if necessary
493	UpdateAlphaTime();
494
495	// Fill histograms
496	const Int_t ibin = fHist.Fill(theta, energy, TMath::Abs(alpha), w);
497	SetBin(ibin);
498
499	if (!fSkipHistTime)
500	fHistTime.Fill(TMath::Abs(alpha), w);
501
502	return kTRUE;
503	}
504
505	// --------------------------------------------------------------------------
506	//
507	// Paint the integral and the error on top of the histogram
508	//
509	void MHAlpha::PaintText(Double_t val, Double_t error) const
510	{
511	TLatex text(0.45, 0.94, MString::Format("N_{exc} = %.1f \\pm %.1f", val, error));
512	text.SetBit(TLatex::kTextNDC);
513	text.SetTextSize(0.04);
514	text.Paint();
515	}
516
517	// --------------------------------------------------------------------------
518	//
519	// Paint the integral and the error on top of the histogram
520	//
521	void MHAlpha::PaintText(const TH1D &h) const
522	{
523	Double_t sumv = 0;
524	Double_t sume = 0;
525
526	for (int i=0; i<h.GetNbinsX(); i++)
527	{
528	sumv += h.GetBinContent(i+1);
529	sume += h.GetBinError(i+1);
530	}
531	PaintText(sumv, sume);
532	}
533	// --------------------------------------------------------------------------
534	//
535	// Update the projections
536	//
537	void MHAlpha::Paint(Option_t *opt)
538	{
539	// Note: Any object cannot be added to a pad twice!
540	// The object is searched by gROOT->FindObject only in
541	// gPad itself!
542	TVirtualPad *padsave = gPad;
543
544	TH1D *h0=0;
545
546	TString o(opt);
547
548	if (o==(TString)"proj")
549	{
550	TPaveStats *st=0;
551	for (int x=0; x<4; x++)
552	{
553	TVirtualPad *p=padsave->GetPad(x+1);
554	if (!p \|\| !(st = (TPaveStats*)p->GetPrimitive("stats")))
555	continue;
556
557	if (st->GetOptStat()==11)
558	continue;
559
560	const Double_t y1 = st->GetY1NDC();
561	const Double_t y2 = st->GetY2NDC();
562	const Double_t x1 = st->GetX1NDC();
563	const Double_t x2 = st->GetX2NDC();
564
565	st->SetY1NDC((y2-y1)/3+y1);
566	st->SetX1NDC((x2-x1)/3+x1);
567	st->SetOptStat(11);
568	}
569
570	padsave->cd(1);
571
572	TH1D hon = (TH1D)gPad->FindObject("Proj");
573	if (hon)
574	{
575	TH1D *dum = fHist.ProjectionZ("dumab", 0, fHist.GetNbinsX()+1, 0, fHist.GetNbinsY()+1);
576	dum->SetDirectory(0);
577	hon->Reset();
578	hon->Add(dum);
579	delete dum;
580
581	if (fOffData)
582	{
583	TH1D hoff = (TH1D)gPad->FindObject("ProjOff");
584	if (hoff)
585	{
586	TH1D *dum = fOffData->ProjectionZ("dumxy", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1);
587	dum->SetDirectory(0);
588	hoff->Reset();
589	hoff->Add(dum);
590	delete dum;
591
592	const Double_t alpha = fFit.Scale(hoff, hon);
593
594	hon->SetMaximum();
595	hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05);
596
597	// BE CARFEULL: This is a really weird workaround!
598	hoff->SetMaximum(alpha);
599
600	// For some reason the line-color is resetted
601	hoff->SetLineColor(kRed);
602
603	if ((h0=(TH1D*)gPad->FindObject("ProjOnOff")))
604	{
605	h0->Reset();
606	h0->Add(hoff, hon, -1);
607	const Float_t min = h0->GetMinimum()*1.05;
608	hon->SetMinimum(min<0 ? min : 0);
609	}
610
611	}
612	}
613	else
614	hon->SetMinimum(0);
615	}
616	FitEnergyBins();
617	FitThetaBins();
618	}
619
620	if (o==(TString)"variable")
621	if ((h0 = (TH1D*)gPad->FindObject("Proj")))
622	{
623	// Check whether we are dealing with on-off analysis
624	TH1D hoff = (TH1D)gPad->FindObject("ProjOff");
625
626	// BE CARFEULL: This is a really weird workaround!
627	const Double_t scale = !hoff \|\| hoff->GetMaximum()<0 ? 1 : hoff->GetMaximum();
628
629	// Do not store the 'final' result in fFit
630	MAlphaFitter fit(fFit);
631	fit.Fit(*h0, hoff, scale, kTRUE);
632	fit.PaintResult();
633	}
634
635	if (o==(TString)"time")
636	PaintText(fHTime);
637
638	if (o==(TString)"theta")
639	{
640	TH1 h = (TH1)gPad->FindObject(MString::Format("%s_x", fHist.GetName()));
641	if (h)
642	{
643	TH1D h2 = (TH1D)fHist.Project3D("dum_x");
644	h2->SetDirectory(0);
645	h2->Scale(fHTheta.Integral()/h2->Integral());
646	h->Reset();
647	h->Add(h2);
648	delete h2;
649	}
650	PaintText(fHTheta);
651	}
652
653	if (o==(TString)"energy")
654	{
655	TH1 h = (TH1)gPad->FindObject(MString::Format("%s_y", fHist.GetName()));
656	if (h)
657	{
658	TH1D h2 = (TH1D)fHist.Project3D("dum_y");
659	h2->SetDirectory(0);
660	h2->Scale(fHEnergy.Integral()/h2->Integral());
661	h->Reset();
662	h->Add(h2);
663	delete h2;
664	}
665	PaintText(fHEnergy);
666
667	if (fHEnergy.GetMaximum()>1)
668	{
669	gPad->SetLogx();
670	gPad->SetLogy();
671	}
672	}
673
674	gPad = padsave;
675	}
676
677	// --------------------------------------------------------------------------
678	//
679	// Draw the histogram
680	//
681	void MHAlpha::Draw(Option_t *opt)
682	{
683	TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
684
685	/*
686	if (TString(opt).Contains("sizebins", TString::kIgnoreCase))
687	{
688	AppendPad("sizebins");
689	return;
690	}
691	*/
692
693	// Do the projection before painting the histograms into
694	// the individual pads
695	AppendPad("proj");
696
697	pad->SetBorderMode(0);
698	pad->Divide(2,2);
699
700	TH1D *h=0;
701
702	pad->cd(1);
703	gPad->SetBorderMode(0);
704
705	h = fHist.ProjectionZ("Proj", 0, fHist.GetNbinsX()+1, 0, fHist.GetNbinsY()+1, "E");
706	h->SetBit(TH1::kNoTitle);
707	h->SetStats(kTRUE);
708	h->SetXTitle(fHist.GetZaxis()->GetTitle());
709	h->SetYTitle("Counts");
710	h->SetDirectory(NULL);
711	h->SetMarkerStyle(0);
712	h->SetBit(kCanDelete);
713	h->Draw("");
714
715	if (fOffData)
716	{
717	// To get green on-data
718	//h->SetMarkerColor(kGreen);
719	//h->SetLineColor(kGreen);
720
721	h = fOffData->ProjectionZ("ProjOff", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1, "E");
722	h->SetBit(TH1::kNoTitle);
723	h->SetXTitle(fHist.GetZaxis()->GetTitle());
724	h->SetYTitle("Counts");
725	h->SetDirectory(NULL);
726	h->SetMarkerStyle(0);
727	h->SetBit(kCanDelete);
728	h->SetMarkerColor(kRed);
729	h->SetLineColor(kRed);
730	//h->SetFillColor(18);
731	h->Draw("same"/"bar same"/);
732
733	// This is the only way to make it work...
734	// Clone and copy constructor give strange crashes :(
735	h = fOffData->ProjectionZ("ProjOnOff", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1, "E");
736	h->SetBit(TH1::kNoTitle);
737	h->SetXTitle(fHist.GetZaxis()->GetTitle());
738	h->SetYTitle("Counts");
739	h->SetDirectory(NULL);
740	h->SetMarkerStyle(0);
741	h->SetBit(kCanDelete);
742	h->SetMarkerColor(kBlue);
743	h->SetLineColor(kBlue);
744	h->Draw("same");
745
746	TLine lin;
747	lin.SetLineStyle(kDashed);
748	lin.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0);
749	}
750
751	// After the Alpha-projection has been drawn. Fit the histogram
752	// and paint the result into this pad
753	AppendPad("variable");
754
755	if (fHEnergy.GetNbinsX()>1 \|\| fHEnergy.GetBinContent(1)>0)
756	{
757	pad->cd(2);
758	gPad->SetBorderMode(0);
759	gPad->SetGridx();
760	gPad->SetGridy();
761	fHEnergy.Draw();
762
763	AppendPad("energy");
764
765	h = (TH1D*)fHist.Project3D("y");
766	h->SetBit(TH1::kNoTitle\|TH1::kNoStats);
767	h->SetXTitle("E [GeV]");
768	h->SetYTitle("Counts");
769	h->SetDirectory(NULL);
770	h->SetMarkerStyle(kFullDotSmall);
771	h->SetBit(kCanDelete);
772	h->SetMarkerColor(kCyan);
773	h->SetLineColor(kCyan);
774	h->Draw("Psame");
775	}
776	else
777	delete pad->GetPad(2);
778
779	if ((fTimeEffOn && fTime) \|\| fHTime.GetNbinsX()>1 \|\| fHTime.GetBinError(1)>0)
780	{
781	pad->cd(3);
782	gPad->SetBorderMode(0);
783	gPad->SetGridx();
784	gPad->SetGridy();
785	fHTime.Draw();
786	AppendPad("time");
787	}
788	else
789	delete pad->GetPad(3);
790
791	if (fHTheta.GetNbinsX()>1 \|\| fHTheta.GetBinContent(1)>0)
792	{
793	pad->cd(4);
794	gPad->SetGridx();
795	gPad->SetGridy();
796	gPad->SetBorderMode(0);
797	fHTheta.Draw();
798
799	AppendPad("theta");
800
801	h = (TH1D*)fHist.Project3D("x");
802	h->SetBit(TH1::kNoTitle\|TH1::kNoStats);
803	h->SetXTitle("\\theta [\\circ]");
804	h->SetYTitle("Counts");
805	h->SetDirectory(NULL);
806	h->SetMarkerStyle(kFullDotSmall);
807	h->SetBit(kCanDelete);
808	h->SetMarkerColor(kCyan);
809	h->SetLineColor(kCyan);
810	h->Draw("Psame");
811	}
812	else
813	delete pad->GetPad(4);
814	}
815
816	void MHAlpha::DrawAll(Bool_t newc)
817	{
818	if (!newc && !fDisplay)
819	return;
820
821	// FIXME: Do in Paint if special option given!
822	TCanvas &c = newc ? *new TCanvas : fDisplay->AddTab("SizeBins");
823	Int_t n = fHist.GetNbinsY();
824	Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
825	c.Divide(nc, nc, 1e-10, 1e-10);
826	gPad->SetBorderMode(0);
827	gPad->SetFrameBorderMode(0);
828
829	// Do not store the 'final' result in fFit
830	MAlphaFitter fit(fFit);
831
832	for (int i=1; i<=fHist.GetNbinsY(); i++)
833	{
834	c.cd(i);
835	gPad->SetBorderMode(0);
836	gPad->SetFrameBorderMode(0);
837
838	TH1D *hon = fHist.ProjectionZ("Proj", 0, fHist.GetNbinsX()+1, i, i, "E");
839	hon->SetBit(TH1::kNoTitle);
840	hon->SetStats(kTRUE);
841	hon->SetXTitle(fHist.GetZaxis()->GetTitle());
842	hon->SetYTitle("Counts");
843	hon->SetDirectory(NULL);
844	hon->SetMarkerStyle(0);
845	hon->SetBit(kCanDelete);
846	hon->Draw("");
847
848	TH1D *hof = 0;
849	Double_t alpha = 1;
850
851	if (fOffData)
852	{
853	hon->SetMarkerColor(kGreen);
854
855	hof = fOffData->ProjectionZ("ProjOff", 0, fOffData->GetNbinsX()+1, i, i, "E");
856	hof->SetBit(TH1::kNoTitle\|TH1::kNoStats);
857	hof->SetXTitle(fHist.GetZaxis()->GetTitle());
858	hof->SetYTitle("Counts");
859	hof->SetDirectory(NULL);
860	hof->SetMarkerStyle(0);
861	hof->SetBit(kCanDelete);
862	hof->SetMarkerColor(kRed);
863	hof->SetLineColor(kRed);
864	hof->Draw("same");
865
866	alpha = fit.Scale(hof, hon);
867
868	hon->SetMaximum();
869	hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05);
870
871	TH1D diff = new TH1D(hon);
872	diff->Add(hof, -1);
873	diff->SetBit(TH1::kNoTitle);
874	diff->SetXTitle(fHist.GetZaxis()->GetTitle());
875	diff->SetYTitle("Counts");
876	diff->SetDirectory(NULL);
877	diff->SetMarkerStyle(0);
878	diff->SetBit(kCanDelete);
879	diff->SetMarkerColor(kBlue);
880	diff->SetLineColor(kBlue);
881	diff->Draw("same");
882
883	TLine lin;
884	lin.SetLineStyle(kDashed);
885	lin.DrawLine(diff->GetXaxis()->GetXmin(), 0, diff->GetXaxis()->GetXmax(), 0);
886
887	const Float_t min = diff->GetMinimum()*1.05;
888	hon->SetMinimum(min<0 ? min : 0);
889	}
890
891	if (hof ? fit.Fit(hon, hof, alpha) : fit.Fit(*hon))
892	{
893	*fLog << dbg << "Bin " << i << ": sigmaexc=" << fit.GetEventsExcess()/fit.GetErrorExcess() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << " scale=" << fit.GetScaleFactor() << endl;
894	fit.DrawResult();
895	}
896	/*
897	if (fit.FitEnergy(fHist, fOffData, i, kTRUE))
898	{
899	fHEnergy.SetBinContent(i, fit.GetEventsExcess());
900	fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
901	}*/
902	}
903
904	}
905
906	void MHAlpha::DrawNicePlot(Bool_t newc, const char title, const char watermark, Int_t binlo, Int_t binhi)
907	{
908	if (!newc && !fDisplay)
909	return;
910
911	if (!fOffData)
912	return;
913
914	// Open and setup canvas/pad
915	TCanvas &c = newc ? *new TCanvas : fDisplay->AddTab("ThetsSq");
916
917	c.SetBorderMode(0);
918	c.SetFrameBorderMode(0);
919	c.SetFillColor(kWhite);
920
921	c.SetLeftMargin(0.12);
922	c.SetRightMargin(0.01);
923	c.SetBottomMargin(0.16);
924	c.SetTopMargin(0.18);
925
926	c.SetGridy();
927
928	gStyle->SetOptStat(0);
929
930	// Get on-data
931	TH1D hon = (TH1D)fHist.ProjectionZ("Proj", 0, fHist.GetNbinsX()+1, binlo, binhi, "E");
932	hon->SetDirectory(NULL);
933	hon->SetBit(kCanDelete);
934	hon->SetMarkerSize(0);
935	hon->SetLineWidth(2);
936	hon->SetLineColor(kBlack);
937	hon->SetMarkerColor(kBlack);
938
939	// Get off-data
940	TH1D *hoff = 0;
941	if (fOffData)
942	{
943	hoff = (TH1D*)fOffData->ProjectionZ("ProjOff", 0, fHist.GetNbinsX()+1, binlo, binhi, "E");
944	hoff->SetDirectory(NULL);
945	hoff->SetBit(kCanDelete);
946	hoff->SetFillColor(17);
947	hoff->SetMarkerSize(0);
948	hoff->SetLineColor(kBlack);
949	hoff->SetMarkerColor(kBlack);
950	}
951
952	// Setup plot which is drawn first
953	TH1D *h = hoff ? hoff : hon;
954	h->GetXaxis()->SetLabelSize(0.06);
955	h->GetXaxis()->SetTitleSize(0.06);
956	h->GetXaxis()->SetTitleOffset(0.95);
957	h->GetYaxis()->SetLabelSize(0.06);
958	h->GetYaxis()->SetTitleSize(0.06);
959	h->GetYaxis()->CenterTitle();
960	h->SetYTitle("Counts");
961	h->SetTitleSize(0.07);
962	h->SetTitle("");
963
964	const Double_t imax = fFit.GetSignalIntegralMax();
965	if (imax<1)
966	h->GetXaxis()->SetRangeUser(0, 0.6*0.6);
967
968	// scale off-data
969
970	MAlphaFitter fit(fFit);
971	fit.ScaleAndFit(*hon, hoff);
972
973	hon->SetMinimum(0);
974	hoff->SetMinimum(0);
975
976	// draw data
977	if (hoff)
978	{
979	hoff->SetMaximum(TMath::Max(hon->GetMaximum(),hoff->GetMaximum())*1.1);
980	hoff->Draw("bar");
981	hon->Draw("same");
982	}
983	else
984	{
985	hon->SetMaximum();
986	hon->Draw();
987	}
988
989	// draw a preliminary tag
990	TLatex text;
991	text.SetTextColor(kWhite);
992	text.SetBit(TLatex::kTextNDC);
993	text.SetTextSize(0.07);
994	text.SetTextAngle(2.5);
995
996	TString wm(watermark);
997	if (binlo>=1 \|\| binhi<hon->GetNbinsX())
998	{
999	wm += wm.IsNull() ? "(" : " (";
1000	if (binlo>=1)
1001	wm += MString::Format("%.1fGeV", fHist.GetYaxis()->GetBinLowEdge(binlo));
1002	wm += "-";
1003	if (binhi<hon->GetNbinsX())
1004	wm += MString::Format("%.1fGeV", fHist.GetYaxis()->GetBinLowEdge(binhi+1));
1005	wm += ")";
1006	}
1007	if (!wm.IsNull())
1008	text.DrawLatex(0.45, 0.2, wm);
1009	//enum { kTextNDC = BIT(14) };
1010
1011	// draw line showing cut
1012	TLine line;
1013	line.SetLineColor(14);
1014	line.SetLineStyle(7);
1015	line.DrawLine(imax, 0, imax, h->GetMaximum());
1016
1017	// Add a title above the plot
1018	TPaveText *pave=new TPaveText(0.12, 0.83, 0.99, 0.975, "blNDC");
1019	pave->SetBorderSize(1);
1020	pave->SetLabel(title);
1021
1022	TText *ptxt = pave->AddText(" ");
1023	ptxt->SetTextAlign(23);
1024
1025	ptxt = pave->AddText(MString::Format("Significance %.1f\\sigma, off-scale %.2f",
1026	fit.GetSignificance(), fit.GetScaleFactor()));
1027	ptxt->SetTextAlign(23);
1028
1029	ptxt = pave->AddText(MString::Format("%.1f excess events, %.1f background events",
1030	fit.GetEventsExcess(), fit.GetEventsBackground()));
1031	ptxt->SetTextAlign(23);
1032	pave->SetBit(kCanDelete);
1033	pave->Draw();
1034	}
1035
1036	Bool_t MHAlpha::Finalize()
1037	{
1038	if (!FitAlpha())
1039	{
1040	*fLog << warn << "MAlphaFitter - Fit failed..." << endl;
1041	return kTRUE;
1042	}
1043
1044	// Store the final result in fFit
1045	*fLog << all;
1046	fFit.Print("result");
1047
1048	fResult->SetVal(fFit.GetMinimizationValue());
1049	fSigma->SetVal(fFit.GetGausSigma());
1050
1051	if (!fSkipHistEnergy)
1052	{
1053	*fLog << inf3 << "Processing energy bins..." << endl;
1054	FitEnergyBins();
1055	}
1056	if (!fSkipHistTheta)
1057	{
1058	*fLog << inf3 << "Processing theta bins..." << endl;
1059	FitThetaBins();
1060	}
1061	if (!fSkipHistTime)
1062	{
1063	*fLog << inf3 << "Processing time bins..." << endl;
1064	UpdateAlphaTime(kTRUE);
1065	MH::RemoveFirstBin(fHTime);
1066	}
1067
1068	if (fOffData)
1069	DrawAll(kFALSE);
1070
1071	return kTRUE;
1072	}
1073
1074	// --------------------------------------------------------------------------
1075	//
1076	// You can use this function if you want to use a MHMatrix instead of
1077	// MMcEvt. This function adds all necessary columns to the
1078	// given matrix. Afterward you should fill the matrix with the corresponding
1079	// data (eg from a file by using MHMatrix::Fill). If you now loop
1080	// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
1081	// will take the values from the matrix instead of the containers.
1082	//
1083	// It takes fSkipHist* into account!
1084	//
1085	void MHAlpha::InitMapping(MHMatrix *mat, Int_t type)
1086	{
1087	if (fMatrix)
1088	return;
1089
1090	fMatrix = mat;
1091
1092	fMap[0] = fMatrix->AddColumn(GetParameterRule());
1093	fMap[1] = -1;
1094	fMap[2] = -1;
1095	fMap[3] = -1;
1096	fMap[4] = -1;
1097
1098	if (!fSkipHistEnergy)
1099	{
1100	fMap[1] = type==0 ? fMatrix->AddColumn("MEnergyEst.fVal") : -1;
1101	fMap[2] = type==0 ? -1 : fMatrix->AddColumn("MHillas.fSize");
1102	}
1103
1104	if (!fSkipHistTheta)
1105	fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
1106
1107	// if (!fSkipHistTime)
1108	// fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
1109	}
1110
1111	void MHAlpha::StopMapping()
1112	{
1113	fMatrix = NULL;
1114	}
1115
1116	void MHAlpha::ApplyScaling()
1117	{
1118	if (!fOffData)
1119	return;
1120
1121	fFit.ApplyScaling(fHist, const_cast<TH3D>(fOffData));
1122	}
1123
1124	Int_t MHAlpha::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
1125	{
1126	Bool_t rc = kFALSE;
1127	if (IsEnvDefined(env, prefix, "NumTimeBins", print))
1128	{
1129	SetNumTimeBins(GetEnvValue(env, prefix, "NumTimeBins", fNumTimeBins));
1130	rc = kTRUE;
1131	}
1132	if (IsEnvDefined(env, prefix, "ForceUsingSize", print))
1133	{
1134	fForceUsingSize = GetEnvValue(env, prefix, "ForceUsingSize", fForceUsingSize);
1135	rc = kTRUE;
1136	}
1137	return rc;
1138	}
1139
1140	Int_t MHAlpha::DistancetoPrimitive(Int_t px, Int_t py)
1141	{
1142	// If pad has no subpad return (we are in one of the subpads)
1143	if (gPad->GetPad(1)==NULL)
1144	return 9999;
1145
1146	// If pad has a subpad we are in the main pad. Check its value.
1147	return gPad->GetPad(1)->DistancetoPrimitive(px,py)==0 ? 0 : 9999;
1148	}
1149
1150	void MHAlpha::RecursiveRemove(TObject *obj)
1151	{
1152	if (obj==fOffData)
1153	fOffData = 0;
1154
1155	MH::RecursiveRemove(obj);
1156	}

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

Download in other formats: