Context Navigation

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

Visit:

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

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

Download in other formats: