Context Navigation

source: trunk/Mars/mhflux/MHAlpha.cc@ 19948

Visit:

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

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

Download in other formats: