Context Navigation

source: tags/Mars-V2.2/mhflux/MHAlpha.cc

Visit:

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

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

Download in other formats: