Context Navigation

source: trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc@ 5098

Visit:

Last change on this file since 5098 was 5080, checked in by tbretz, 20 years ago
* empty log message *
File size: 32.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-2004
21	!
22	!
23	\* ======================================================================== */
24
25	//////////////////////////////////////////////////////////////////////////////
26	//
27	// MHFalseSource
28	//
29	// Create a false source plot. For the Binning in x,y the object MBinning
30	// "BinningFalseSource" is taken from the paremeter list.
31	//
32	// The binning in alpha is currently fixed as 18bins from 0 to 90deg.
33	//
34	// If MTime, MObservatory and MPointingPos is available in the paremeter
35	// list each image is derotated.
36	//
37	// MHFalseSource fills a 3D histogram with alpha distribution for
38	// each (derotated) x and y position.
39	//
40	// Only a radius of 1.2deg around the camera center is taken into account.
41	//
42	// The displayed histogram is the projection of alpha<fAlphaCut into
43	// the x,y plain and the projection of alpha>90-fAlphaCut
44	//
45	// By using the context menu (find it in a small region between the single
46	// pads) you can change the AlphaCut 'online'
47	//
48	// Each Z-Projection (Alpha-histogram) is scaled such, that the number
49	// of entries fits the maximum number of entries in all Z-Projections.
50	// This should correct for the different acceptance in the center
51	// and at the vorder of the camera. This, however, produces more noise
52	// at the border.
53	//
54	// Here is a slightly simplified version of the algorithm:
55	// ------------------------------------------------------
56	// MHillas hil; // Taken as second argument in MFillH
57	//
58	// MSrcPosCam src;
59	// MHillasSrc hsrc;
60	// hsrc.SetSrcPos(&src);
61	//
62	// for (int ix=0; ix<nx; ix++)
63	// for (int iy=0; iy<ny; iy++)
64	// {
65	// TVector2 v(cx[ix], cy[iy]); //cx center of x-bin ix
66	// if (rho!=0) //cy center of y-bin iy
67	// v=v.Rotate(rho); //image rotation angle
68	//
69	// src.SetXY(v); //source position in the camera
70	//
71	// if (!hsrc.Calc(hil)) //calc source dependant hillas
72	// return;
73	//
74	// //fill absolute alpha into histogram
75	// const Double_t alpha = hsrc.GetAlpha();
76	// fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w);
77	// }
78	// }
79	//
80	// Use MHFalse Source like this:
81	// -----------------------------
82	// MFillH fill("MHFalseSource", "MHillas");
83	// or
84	// MHFalseSource hist;
85	// hist.SetAlphaCut(12.5); // The default value
86	// hist.SetBgmean(55); // The default value
87	// MFillH fill(&hist, "MHillas");
88	//
89	// To be implemented:
90	// ------------------
91	// - a switch to use alpha or \|alpha\|
92	// - taking the binning for alpha from the parlist (binning is
93	// currently fixed)
94	// - a possibility to change the fit-function (eg using a different TF1)
95	// - a possibility to change the fit algorithm (eg which paremeters
96	// are fixed at which time)
97	// - use a different varaible than alpha
98	// - a possiblity to change the algorithm which is used to calculate
99	// alpha (or another parameter) is missing (this could be done using
100	// a tasklist somewhere...)
101	// - a more clever (and faster) algorithm to fill the histogram, eg by
102	// calculating alpha once and fill the two coils around the mean
103	// - more drawing options...
104	// - Move Significance() to a more 'general' place and implement
105	// also different algorithms like (Li/Ma)
106	// - implement fit for best alpha distribution -- online (Paint)
107	// - currently a constant pointing position is assumed in Fill()
108	// - the center of rotation need not to be the camera center
109	// - ERRORS on each alpha bin to estimate the chi^2 correctly!
110	// (sqrt(N)/binwidth) needed for WOlfgangs proposed caluclation
111	// of alpha(Li/Ma)
112	//
113	//////////////////////////////////////////////////////////////////////////////
114	#include "MHFalseSource.h"
115
116	#include <TF1.h>
117	#include <TF2.h>
118	#include <TH2.h>
119	#include <TGraph.h>
120	#include <TStyle.h>
121	#include <TCanvas.h>
122	#include <TPaveText.h>
123	#include <TStopwatch.h>
124
125	#include "MGeomCam.h"
126	#include "MSrcPosCam.h"
127	#include "MHillasSrc.h"
128	#include "MTime.h"
129	#include "MObservatory.h"
130	#include "MPointingPos.h"
131	#include "MAstroCatalog.h"
132	#include "MAstroSky2Local.h"
133	#include "MStatusDisplay.h"
134
135	#include "MMath.h"
136	#include "MAlphaFitter.h"
137
138	#include "MBinning.h"
139	#include "MParList.h"
140
141	#include "MLog.h"
142	#include "MLogManip.h"
143
144	ClassImp(MHFalseSource);
145
146	using namespace std;
147
148	// --------------------------------------------------------------------------
149	//
150	// Default Constructor
151	//
152	MHFalseSource::MHFalseSource(const char name, const char title)
153	: fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), fAlphaCut(12.5),
154	fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1)
155	{
156	//
157	// set the name and title of this object
158	//
159	fName = name ? name : "MHFalseSource";
160	fTitle = title ? title : "3D-plot of Alpha vs x, y";
161
162	fHist.SetDirectory(NULL);
163
164	fHist.SetName("Alpha");
165	fHist.SetTitle("3D-plot of Alpha vs x, y");
166	fHist.SetXTitle("x [\\circ]");
167	fHist.SetYTitle("y [\\circ]");
168	fHist.SetZTitle("\\alpha [\\circ]");
169	}
170
171	void MHFalseSource::MakeSymmetric(TH1 *h)
172	{
173	const Float_t min = TMath::Abs(h->GetMinimum());
174	const Float_t max = TMath::Abs(h->GetMaximum());
175
176	const Float_t absmax = TMath::Max(min, max)*1.002;
177
178	h->SetMaximum( absmax);
179	h->SetMinimum(-absmax);
180	}
181
182	// --------------------------------------------------------------------------
183	//
184	// Set the alpha cut (\|alpha\|<fAlphaCut) which is use to estimate the
185	// number of excess events
186	//
187	void MHFalseSource::SetAlphaCut(Float_t alpha)
188	{
189	if (alpha<0)
190	*fLog << warn << "Alpha<0... taking absolute value." << endl;
191
192	fAlphaCut = TMath::Abs(alpha);
193
194	Modified();
195	}
196
197	// --------------------------------------------------------------------------
198	//
199	// Set mean alpha around which the off sample is determined
200	// (fBgMean-fAlphaCut/2<\|fAlpha\|<fBgMean+fAlphaCut/2) which is use
201	// to estimate the number of off events
202	//
203	void MHFalseSource::SetBgMean(Float_t alpha)
204	{
205	if (alpha<0)
206	*fLog << warn << "Alpha<0... taking absolute value." << endl;
207
208	fBgMean = TMath::Abs(alpha);
209
210	Modified();
211	}
212
213	// --------------------------------------------------------------------------
214	//
215	// Calculate Significance as
216	// significance = (s-b)/sqrt(s+kkb) mit k=s/b
217	//
218	// s: total number of events in signal region
219	// b: number of background events in signal region
220	//
221	Double_t MHFalseSource::Significance(Double_t s, Double_t b)
222	{
223	return MMath::SignificanceSym(s, b);
224	}
225
226	// --------------------------------------------------------------------------
227	//
228	// calculates the significance according to Li & Ma
229	// ApJ 272 (1983) 317, Formula 17
230	//
231	// s // s: number of on events
232	// b // b: number of off events
233	// alpha = t_on/t_off; // t: observation time
234	//
235	Double_t MHFalseSource::SignificanceLiMa(Double_t s, Double_t b, Double_t alpha)
236	{
237	return MMath::SignificanceLiMaSigned(s, b);
238	}
239
240	// --------------------------------------------------------------------------
241	//
242	// Set binnings (takes BinningFalseSource) and prepare filling of the
243	// histogram.
244	//
245	// Also search for MTime, MObservatory and MPointingPos
246	//
247	Bool_t MHFalseSource::SetupFill(const MParList *plist)
248	{
249	const MGeomCam geom = (MGeomCam)plist->FindObject("MGeomCam");
250	if (!geom)
251	{
252	*fLog << err << "MGeomCam not found... aborting." << endl;
253	return kFALSE;
254	}
255
256	fMm2Deg = geom->GetConvMm2Deg();
257
258	MBinning binsa;
259	binsa.SetEdges(18, 0, 90);
260
261	const MBinning bins = (MBinning)plist->FindObject("BinningFalseSource");
262	if (!bins)
263	{
264	const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg;
265
266	MBinning b;
267	b.SetEdges(20, -r, r);
268	SetBinning(&fHist, &b, &b, &binsa);
269	}
270	else
271	SetBinning(&fHist, bins, bins, &binsa);
272
273	fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos"));
274	if (!fPointPos)
275	*fLog << warn << "MPointingPos not found... no derotation." << endl;
276
277	fTime = (MTime*)plist->FindObject(AddSerialNumber("MTime"));
278	if (!fTime)
279	*fLog << warn << "MTime not found... no derotation." << endl;
280
281	fSrcPos = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam"));
282	if (!fSrcPos)
283	*fLog << warn << "MSrcPosCam not found... no translation." << endl;
284
285	fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory"));
286	if (!fObservatory)
287	*fLog << warn << "MObservatory not found... no derotation." << endl;
288
289	// FIXME: Because the pointing position could change we must check
290	// for the current pointing position and add a offset in the
291	// Fill function!
292	fRa = fPointPos->GetRa();
293	fDec = fPointPos->GetDec();
294
295	return kTRUE;
296	}
297
298	// --------------------------------------------------------------------------
299	//
300	// Fill the histogram. For details see the code or the class description
301	//
302	Bool_t MHFalseSource::Fill(const MParContainer *par, const Stat_t w)
303	{
304	const MHillas hil = dynamic_cast<const MHillas>(par);
305	if (!hil)
306	{
307	*fLog << err << "MHFalseSource::Fill: No container specified!" << endl;
308	return kFALSE;
309	}
310
311	// Get max radius...
312	const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1));
313
314	// Get camera rotation angle
315	Double_t rho = 0;
316	if (fTime && fObservatory && fPointPos)
317	rho = fPointPos->RotationAngle(fObservatory, fTime);
318	//if (fPointPos)
319	// rho = fPointPos->RotationAngle(*fObservatory);
320
321	// Create necessary containers for calculation
322	MSrcPosCam src;
323	MHillasSrc hsrc;
324	hsrc.SetSrcPos(&src);
325
326	// Get number of bins and bin-centers
327	const Int_t nx = fHist.GetNbinsX();
328	const Int_t ny = fHist.GetNbinsY();
329	Axis_t cx[nx];
330	Axis_t cy[ny];
331	fHist.GetXaxis()->GetCenter(cx);
332	fHist.GetYaxis()->GetCenter(cy);
333
334	for (int ix=0; ix<nx; ix++)
335	{
336	for (int iy=0; iy<ny; iy++)
337	{
338	// check distance... to get a circle plot
339	if (TMath::Hypot(cx[ix], cy[iy])>maxr)
340	continue;
341
342	// rotate center of bin
343	// precalculation of sin/cos doesn't accelerate
344	TVector2 v(cx[ix], cy[iy]);
345	if (rho!=0)
346	v=v.Rotate(rho);
347
348	// convert degrees to millimeters
349	v *= 1./fMm2Deg;
350
351	if (fSrcPos)
352	v += fSrcPos->GetXY();
353
354	src.SetXY(v);
355
356	// Source dependant hillas parameters
357	if (hsrc.Calc(*hil)>0)
358	{
359	*fLog << warn << "Calculation of MHillasSrc failed for x=" << cx[ix] << " y=" << cy[iy] << endl;
360	return kFALSE;
361	}
362
363	// FIXME: This should be replaced by an external MFilter
364	// and/or MTaskList
365	// Source dependant distance cut
366	if (fMinDist>0 && hsrc.GetDist()*fMm2Deg<fMinDist)
367	continue;
368	if (fMaxDist>0 && hsrc.GetDist()*fMm2Deg>fMaxDist)
369	continue;
370
371	if (fMaxLD>0 && hil->GetLength()>fMaxLD*hsrc.GetDist())
372	continue;
373	if (fMinLD>0 && hil->GetLength()<fMinLD*hsrc.GetDist())
374	continue;
375
376	// Fill histogram
377	const Double_t alpha = hsrc.GetAlpha();
378	fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w);
379	}
380	}
381
382	return kTRUE;
383	}
384
385	// --------------------------------------------------------------------------
386	//
387	// Create projection for off data, taking sum of bin contents of
388	// range (fBgMean-fAlphaCut/2, fBgMean+fAlphaCut) Making sure to take
389	// the same number of bins than for on-data
390	//
391	void MHFalseSource::ProjectOff(TH2D h2, TH2D all)
392	{
393	TAxis &axe = *fHist.GetZaxis();
394
395	// Find range to cut (left edge and width)
396	const Int_t f = axe.FindBin(fBgMean-fAlphaCut/2);
397	const Int_t l = axe.FindBin(fAlphaCut)+f-1;
398
399	axe.SetRange(f, l);
400	const Float_t cut1 = axe.GetBinLowEdge(f);
401	const Float_t cut2 = axe.GetBinUpEdge(l);
402	h2->SetTitle(Form("Distribution of %.1f\\circ<\|\\alpha\|<%.1f\\circ in x,y", cut1, cut2));
403
404	// Get projection for range
405	TH2D p = (TH2D)fHist.Project3D("yx_off");
406
407	// Reset range
408	axe.SetRange(0,9999);
409
410	// Move contents from projection to h2
411	h2->Reset();
412	h2->Add(p, all->GetMaximum());
413	h2->Divide(all);
414
415	// Delete p
416	delete p;
417
418	// Set Minimum as minimum value Greater Than 0
419	h2->SetMinimum(GetMinimumGT(*h2));
420	}
421
422	// --------------------------------------------------------------------------
423	//
424	// Create projection for on data, taking sum of bin contents of
425	// range (0, fAlphaCut)
426	//
427	void MHFalseSource::ProjectOn(TH2D h3, TH2D all)
428	{
429	TAxis &axe = *fHist.GetZaxis();
430
431	// Find range to cut
432	axe.SetRangeUser(0, fAlphaCut);
433	const Float_t cut = axe.GetBinUpEdge(axe.GetLast());
434	h3->SetTitle(Form("Distribution of \|\\alpha\|<%.1f\\circ in x,y", cut));
435
436	// Get projection for range
437	TH2D p = (TH2D)fHist.Project3D("yx_on");
438
439	// Reset range
440	axe.SetRange(0,9999);
441
442	// Move contents from projection to h3
443	h3->Reset();
444	h3->Add(p, all->GetMaximum());
445	h3->Divide(all);
446
447	// Delete p
448	delete p;
449
450	// Set Minimum as minimum value Greater Than 0
451	h3->SetMinimum(GetMinimumGT(*h3));
452	}
453
454	// --------------------------------------------------------------------------
455	//
456	// Create projection for all data, taking sum of bin contents of
457	// range (0, 90) - corresponding to the number of entries in this slice.
458	//
459	void MHFalseSource::ProjectAll(TH2D *h3)
460	{
461	h3->SetTitle("Number of entries");
462
463	// Get projection for range
464	TH2D p = (TH2D)fHist.Project3D("yx_all");
465
466	// Move contents from projection to h3
467	h3->Reset();
468	h3->Add(p);
469	delete p;
470
471	// Set Minimum as minimum value Greater Than 0
472	h3->SetMinimum(GetMinimumGT(*h3));
473	}
474
475	// --------------------------------------------------------------------------
476	//
477	// Update the projections and paint them
478	//
479	void MHFalseSource::Paint(Option_t *opt)
480	{
481	// Set pretty color palette
482	gStyle->SetPalette(1, 0);
483
484	TVirtualPad *padsave = gPad;
485
486	TH1D* h1;
487	TH2D* h0;
488	TH2D* h2;
489	TH2D* h3;
490	TH2D* h4;
491	TH2D* h5;
492
493	/*
494	fHistProjAll = Form("All_%p", this);
495	fHistProjOn = Form("On_%p", this);
496	fHistProjOff = Form("Off_%p", this);
497	fHistProjDiff = Form("Diff_%p", this);
498	fHistProjAll = Form("All_%p", this);
499	*/
500
501	// Update projection of all-events
502	padsave->GetPad(2)->cd(3);
503	if ((h0 = (TH2D*)gPad->FindObject("Alpha_yx_all")))
504	ProjectAll(h0);
505
506	// Update projection of on-events
507	padsave->GetPad(1)->cd(1);
508	if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on")))
509	ProjectOn(h3, h0);
510
511	// Update projection of off-events
512	padsave->GetPad(1)->cd(3);
513	if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off")))
514	ProjectOff(h2, h0);
515
516	padsave->GetPad(2)->cd(2);
517	if ((h5 = (TH2D*)gPad->FindObject("Alpha_yx_diff")))
518	{
519	h5->Add(h2, h3, -1);
520	MakeSymmetric(h5);
521	}
522
523	// Update projection of significance
524	padsave->GetPad(1)->cd(2);
525	if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_yx_sig")))
526	{
527	const Int_t nx = h4->GetXaxis()->GetNbins();
528	const Int_t ny = h4->GetYaxis()->GetNbins();
529	//const Int_t nr = nxnx + nyny;
530
531	Int_t maxx=nx/2;
532	Int_t maxy=ny/2;
533
534	Int_t max = h4->GetBin(nx, ny);
535
536	for (int ix=1; ix<=nx; ix++)
537	for (int iy=1; iy<=ny; iy++)
538	{
539	const Int_t n = h4->GetBin(ix, iy);
540
541	const Double_t s = h3->GetBinContent(n);
542	const Double_t b = h2->GetBinContent(n);
543
544	const Double_t sig = SignificanceLiMa(s, b);
545
546	h4->SetBinContent(n, sig);
547
548	if (sig>h4->GetBinContent(max) && sig>0/* && (ix-nx/2)(ix-nx/2)+(iy-ny/2)(iy-ny/2)<nrnr/9/)
549	{
550	max = n;
551	maxx=ix;
552	maxy=iy;
553	}
554	}
555
556	MakeSymmetric(h4);
557
558	// Update projection of 'the best alpha-plot'
559	padsave->GetPad(2)->cd(1);
560	if ((h1 = (TH1D*)gPad->FindObject("Alpha_z")) && max>0)
561	{
562	const Double_t x = h4->GetXaxis()->GetBinCenter(maxx);
563	const Double_t y = h4->GetYaxis()->GetBinCenter(maxy);
564	const Double_t s = h4->GetBinContent(max);
565
566	// This is because a 3D histogram x and y are vice versa
567	// Than for their projections
568	TH1 *h = fHist.ProjectionZ("Alpha_z", maxx, maxx, maxy, maxy);
569	h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f)", x, y, s));
570	}
571	}
572
573	gPad = padsave;
574	}
575
576	// --------------------------------------------------------------------------
577	//
578	// Get the MAstroCatalog corresponding to fRa, fDec. The limiting magnitude
579	// is set to 9, while the fov is equal to the current fov of the false
580	// source plot.
581	//
582	TObject *MHFalseSource::GetCatalog()
583	{
584	const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1));
585
586	// Create catalog...
587	MAstroCatalog *stars = new MAstroCatalog;
588	stars->SetLimMag(9);
589	stars->SetGuiActive(kFALSE);
590	stars->SetRadiusFOV(maxr);
591	stars->SetRaDec(fRaTMath::DegToRad()15, fDec*TMath::DegToRad());
592	stars->ReadBSC("bsc5.dat");
593
594	*fLog << err << "FIXME - The catalog will never be deleted, because this crashes!" << endl;
595
596	// stars->SetBit(kCanDelete);
597
598	return stars;
599	}
600
601	// --------------------------------------------------------------------------
602	//
603	// Draw the histogram
604	//
605	void MHFalseSource::Draw(Option_t *opt)
606	{
607	TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
608	pad->SetBorderMode(0);
609
610	AppendPad("");
611
612	pad->Divide(1, 2, 0, 0.03);
613
614	*fLog << err << "FIXME - Plotting the catalog is broken!" << endl;
615
616	TObject *catalog = GetCatalog();
617
618	// Initialize upper part
619	pad->cd(1);
620	gPad->SetBorderMode(0);
621	gPad->Divide(3, 1);
622
623	// PAD #1
624	pad->GetPad(1)->cd(1);
625	gPad->SetBorderMode(0);
626	fHist.GetZaxis()->SetRangeUser(0,fAlphaCut);
627	TH1 *h3 = fHist.Project3D("yx_on");
628	fHist.GetZaxis()->SetRange(0,9999);
629	h3->SetDirectory(NULL);
630	h3->SetXTitle(fHist.GetXaxis()->GetTitle());
631	h3->SetYTitle(fHist.GetYaxis()->GetTitle());
632	h3->Draw("colz");
633	h3->SetBit(kCanDelete);
634	catalog->Draw("mirror same *");
635
636	// PAD #2
637	pad->GetPad(1)->cd(2);
638	gPad->SetBorderMode(0);
639	fHist.GetZaxis()->SetRange(0,0);
640	TH1 *h4 = fHist.Project3D("yx_sig"); // Do this to get the correct binning....
641	fHist.GetZaxis()->SetRange(0,9999);
642	h4->SetTitle("Significance");
643	h4->SetDirectory(NULL);
644	h4->SetXTitle(fHist.GetXaxis()->GetTitle());
645	h4->SetYTitle(fHist.GetYaxis()->GetTitle());
646	h4->Draw("colz");
647	h4->SetBit(kCanDelete);
648	catalog->Draw("mirror same *");
649
650	// PAD #3
651	pad->GetPad(1)->cd(3);
652	gPad->SetBorderMode(0);
653	fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
654	TH1 *h2 = fHist.Project3D("yx_off");
655	h2->SetDirectory(NULL);
656	h2->SetXTitle(fHist.GetXaxis()->GetTitle());
657	h2->SetYTitle(fHist.GetYaxis()->GetTitle());
658	h2->Draw("colz");
659	h2->SetBit(kCanDelete);
660	catalog->Draw("mirror same *");
661
662	// Initialize lower part
663	pad->cd(2);
664	gPad->SetBorderMode(0);
665	gPad->Divide(3, 1);
666
667	// PAD #4
668	pad->GetPad(2)->cd(1);
669	gPad->SetBorderMode(0);
670	TH1 *h1 = fHist.ProjectionZ("Alpha_z");
671	h1->SetDirectory(NULL);
672	h1->SetTitle("Distribution of \\alpha");
673	h1->SetXTitle(fHist.GetZaxis()->GetTitle());
674	h1->SetYTitle("Counts");
675	h1->Draw(opt);
676	h1->SetBit(kCanDelete);
677
678	// PAD #5
679	pad->GetPad(2)->cd(2);
680	gPad->SetBorderMode(0);
681	TH1 h5 = (TH1)h3->Clone("Alpha_yx_diff");
682	h5->Add(h2, -1);
683	h5->SetTitle("Difference of on- and off-distribution");
684	h5->SetDirectory(NULL);
685	h5->Draw("colz");
686	h5->SetBit(kCanDelete);
687	catalog->Draw("mirror same *");
688
689	// PAD #6
690	pad->GetPad(2)->cd(3);
691	gPad->SetBorderMode(0);
692	TH1 *h0 = fHist.Project3D("yx_all");
693	h0->SetDirectory(NULL);
694	h0->SetXTitle(fHist.GetXaxis()->GetTitle());
695	h0->SetYTitle(fHist.GetYaxis()->GetTitle());
696	h0->Draw("colz");
697	h0->SetBit(kCanDelete);
698	catalog->Draw("mirror same *");
699	}
700
701	// --------------------------------------------------------------------------
702	//
703	// Everything which is in the main pad belongs to this class!
704	//
705	Int_t MHFalseSource::DistancetoPrimitive(Int_t px, Int_t py)
706	{
707	return 0;
708	}
709
710	// --------------------------------------------------------------------------
711	//
712	// Set all sub-pads to Modified()
713	//
714	void MHFalseSource::Modified()
715	{
716	if (!gPad)
717	return;
718
719	TVirtualPad *padsave = gPad;
720	padsave->Modified();
721	padsave->GetPad(1)->cd(1);
722	gPad->Modified();
723	padsave->GetPad(1)->cd(3);
724	gPad->Modified();
725	padsave->GetPad(2)->cd(1);
726	gPad->Modified();
727	padsave->GetPad(2)->cd(2);
728	gPad->Modified();
729	padsave->GetPad(2)->cd(3);
730	gPad->Modified();
731	gPad->cd();
732	}
733
734	Double_t FcnGauss2d(Double_t x, Double_t par)
735	{
736	TVector2 v = TVector2(x[0], x[1]).Rotate(par[5]*TMath::DegToRad());
737
738	const Double_t g0 = TMath::Gaus(v.X(), par[1], par[2]);
739	const Double_t g1 = TMath::Gaus(v.Y(), par[3], par[4]);
740
741	//Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE);
742	return par[0]*(g0+g1);
743	}
744
745	// --------------------------------------------------------------------------
746	//
747	// This is a preliminary implementation of a alpha-fit procedure for
748	// all possible source positions. It will be moved into its own
749	// more powerfull class soon.
750	//
751	// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
752	// [0]exp(-0.5((x-[1])/[2])^2) + [3] + [4]x + [5]x^2
753	// or
754	// Aexp(-0.5((x-mu)/sigma)^2) + a + bx + cx^2
755	//
756	// Parameter [1] is fixed to 0 while the alpha peak should be
757	// symmetric around alpha=0.
758	//
759	// Parameter [4] is fixed to 0 because the first derivative at
760	// alpha=0 should be 0, too.
761	//
762	// In a first step the background is fitted between bgmin and bgmax,
763	// while the parameters [0]=0 and [2]=1 are fixed.
764	//
765	// In a second step the signal region (alpha<sigmax) is fittet using
766	// the whole function with parameters [1], [3], [4] and [5] fixed.
767	//
768	// The number of excess and background events are calculated as
769	// s = int(0, sigint, gaus(0)+pol2(3))
770	// b = int(0, sigint, pol2(3))
771	//
772	// The Significance is calculated using the Significance() member
773	// function.
774	//
775	void MHFalseSource::FitSignificance(Float_t sigint, Float_t sigmax, Float_t bgmin, Float_t bgmax, Byte_t polynom)
776	{
777	TObject *catalog = GetCatalog();
778
779	TH1D h0a("A", "", 50, 0, 4000);
780	TH1D h4a("chisq1", "", 50, 0, 35);
781	//TH1D h5a("prob1", "", 50, 0, 1.1);
782	TH1D h6("signifcance", "", 50, -20, 20);
783
784	TH1D h1("mu", "Parameter \\mu", 50, -1, 1);
785	TH1D h2("sigma", "Parameter \\sigma", 50, 0, 90);
786	TH1D h3("b", "Parameter b", 50, -0.1, 0.1);
787
788	TH1D h0b("a", "Parameter a (red), A (blue)", 50, 0, 4000);
789	TH1D h4b("\\chi^{2}", "\\chi^{2} (red, green) / significance (black)", 50, 0, 35);
790	//TH1D h5b("prob", "Fit probability: Bg(red), F(blue)", 50, 0, 1.1);
791
792	h0a.SetLineColor(kBlue);
793	h4a.SetLineColor(kBlue);
794	//h5a.SetLineColor(kBlue);
795	h0b.SetLineColor(kRed);
796	h4b.SetLineColor(kRed);
797	//h5b.SetLineColor(kRed);
798
799	TH1 *hist = fHist.Project3D("yx_fit");
800	hist->SetDirectory(0);
801	TH1 *hists = fHist.Project3D("yx_fit");
802	hists->SetDirectory(0);
803	TH1 *histb = fHist.Project3D("yx_fit");
804	histb->SetDirectory(0);
805	hist->Reset();
806	hists->Reset();
807	histb->Reset();
808	hist->SetNameTitle("Significance",
809	Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
810	sigmax, bgmin, bgmax));
811	hists->SetName("Excess");
812	histb->SetName("Background");
813	hist->SetXTitle(fHist.GetXaxis()->GetTitle());
814	hists->SetXTitle(fHist.GetXaxis()->GetTitle());
815	histb->SetXTitle(fHist.GetXaxis()->GetTitle());
816	hist->SetYTitle(fHist.GetYaxis()->GetTitle());
817	hists->SetYTitle(fHist.GetYaxis()->GetTitle());
818	histb->SetYTitle(fHist.GetYaxis()->GetTitle());
819
820	const Double_t w = fHist.GetZaxis()->GetBinWidth(1);
821
822	TArrayD maxpar;
823
824	/* func.SetParName(0, "A");
825	* func.SetParName(1, "mu");
826	* func.SetParName(2, "sigma");
827	*/
828
829	const Int_t nx = hist->GetXaxis()->GetNbins();
830	const Int_t ny = hist->GetYaxis()->GetNbins();
831	//const Int_t nr = nxnx+nyny;
832
833	Double_t maxalpha0=0;
834	Double_t maxs=3;
835
836	Int_t maxx=0;
837	Int_t maxy=0;
838
839	TStopwatch clk;
840	clk.Start();
841
842	*fLog << inf;
843	*fLog << "Signal fit: alpha < " << sigmax << endl;
844	*fLog << "Integration: alpha < " << sigint << endl;
845	*fLog << "Background fit: " << bgmin << " < alpha < " << bgmax << endl;
846	*fLog << "Polynom order: " << (int)polynom << endl;
847	*fLog << "Fitting False Source Plot..." << flush;
848
849	TH1 *h0 = fHist.Project3D("yx_entries");
850	Float_t entries = h0->GetMaximum();
851	delete h0;
852
853	MAlphaFitter fit;
854	fit.SetSignalIntegralMax(sigint);
855	fit.SetSignalFitMax(sigmax);
856	fit.SetBackgroundFitMin(bgmin);
857	fit.SetBackgroundFitMax(bgmax);
858	fit.SetPolynomOrder(polynom);
859
860	TH1D *h=0;
861	for (int ix=1; ix<=nx; ix++)
862	for (int iy=1; iy<=ny; iy++)
863	{
864	// This is because a 3D histogram x and y are vice versa
865	// Than for their projections
866	h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy);
867
868	if (h->GetEntries()==0)
869	continue;
870
871
872	h->Scale(entries/h->GetEntries());
873
874	if (!fit.Fit(*h))
875	continue;
876
877	const Double_t alpha0 = h->GetBinContent(1);
878	if (alpha0>maxalpha0)
879	maxalpha0=alpha0;
880
881	// Fill results into some histograms
882	h0a.Fill(fit.GetGausA()); // gaus-A
883	h0b.Fill(fit.GetCoefficient(3)); // 3
884	h1.Fill(fit.GetGausMu()); // mu
885	h2.Fill(fit.GetGausSigma()); // sigma-gaus
886	if (polynom>1)
887	h3.Fill(fit.GetCoefficient(5));
888	h4b.Fill(fit.GetChiSqSignal());
889
890	const Double_t sig = fit.GetSignificance();
891	const Double_t b = fit.GetEventsBackground();
892	const Double_t s = fit.GetEventsSignal();
893
894	const Int_t n = hist->GetBin(ix, iy);
895	hists->SetBinContent(n, s-b);
896	histb->SetBinContent(n, b);
897
898	hist->SetBinContent(n, sig);
899	if (sig!=0)
900	h6.Fill(sig);
901
902	if (sig>maxs)
903	{
904	maxs = sig;
905	maxx = ix;
906	maxy = iy;
907	maxpar = fit.GetCoefficients();
908	}
909	}
910
911	*fLog << "Done." << endl;
912
913	h0a.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
914	h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
915
916	hists->SetTitle(Form("Excess events for \\alpha<%.0f\\circ (N_{max}=%d)", sigint, (int)hists->GetMaximum()));
917	histb->SetTitle(Form("Background events for \\alpha<%.0f\\circ", sigint));
918
919	//hists->SetMinimum(GetMinimumGT(*hists));
920	histb->SetMinimum(GetMinimumGT(*histb));
921
922	MakeSymmetric(hists);
923	MakeSymmetric(hist);
924
925	clk.Stop();
926	clk.Print("m");
927
928	TCanvas *c=new TCanvas;
929
930	gStyle->SetPalette(1, 0);
931
932	c->Divide(3,2, 0, 0);
933	c->cd(1);
934	gPad->SetBorderMode(0);
935	hists->Draw("colz");
936	hists->SetBit(kCanDelete);
937	catalog->Draw("mirror same *");
938	c->cd(2);
939	gPad->SetBorderMode(0);
940	hist->Draw("colz");
941	hist->SetBit(kCanDelete);
942
943
944	TF2 f2d("", FcnGauss2d, -1.5, 1.5, -1.5, 1.5, 6);
945	f2d.SetParName(0, "Max sigma");
946	f2d.SetParName(1, "Mean_1 deg");
947	f2d.SetParName(2, "Sigma_1 deg");
948	f2d.SetParName(3, "Mean_2 deg");
949	f2d.SetParName(4, "Sigma_2 deg");
950	f2d.SetParName(5, "Phi deg");
951	f2d.SetParLimits(1, -1, 1); // mu_1
952	f2d.SetParLimits(3, -1, 1); // mu_2
953	f2d.SetParLimits(2, 0, 1); // sigma_1
954	f2d.SetParLimits(4, 0, 1); // sigma_2
955	f2d.SetParLimits(5, 0, 90); // phi
956	f2d.SetParameter(0, maxs); // A
957	f2d.SetParameter(1, 0); // mu_1
958	f2d.SetParameter(2, 0.2); // sigma_1
959	f2d.SetParameter(3, 0); // mu_2
960	f2d.SetParameter(4, 0.2); // sigma_2
961	f2d.SetParameter(5, 0); // phi
962	hist->Fit(&f2d, "IN0");
963	//f2d.DrawClone("cont2same")->SetBit(kCanDelete);
964
965
966	catalog->Draw("mirror same *");
967	c->cd(3);
968	gPad->SetBorderMode(0);
969	histb->Draw("colz");
970	histb->SetBit(kCanDelete);
971	catalog->Draw("mirror same *");
972	c->cd(4);
973	gPad->Divide(1,3, 0, 0);
974	TVirtualPad *p=gPad;
975	p->SetBorderMode(0);
976	p->cd(1);
977	gPad->SetBorderMode(0);
978	h0b.DrawCopy();
979	h0a.DrawCopy("same");
980	p->cd(2);
981	gPad->SetBorderMode(0);
982	h3.DrawCopy();
983	p->cd(3);
984	gPad->SetBorderMode(0);
985	h2.DrawCopy();
986	c->cd(6);
987	gPad->Divide(1,2, 0, 0);
988	TVirtualPad *q=gPad;
989	q->SetBorderMode(0);
990	q->cd(1);
991	gPad->SetBorderMode(0);
992	gPad->SetBorderMode(0);
993	h4b.DrawCopy();
994	h4a.DrawCopy("same");
995	h6.DrawCopy("same");
996	q->cd(2);
997	gPad->SetBorderMode(0);
998	//h5b.DrawCopy();
999	//h5a.DrawCopy("same");
1000	c->cd(5);
1001	gPad->SetBorderMode(0);
1002	if (maxx>0 && maxy>0)
1003	{
1004	const char *title = Form(" \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f) ",
1005	hist->GetXaxis()->GetBinCenter(maxx),
1006	hist->GetYaxis()->GetBinCenter(maxy), maxs);
1007
1008	TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
1009	result->Scale(entries/h->GetEntries());
1010
1011	result->SetDirectory(NULL);
1012	result->SetNameTitle("Result \\alpha", title);
1013	result->SetBit(kCanDelete);
1014	result->SetXTitle("\\alpha [\\circ]");
1015	result->SetYTitle("Counts");
1016	result->Draw();
1017
1018	TF1 f1("f1", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
1019	TF1 f2("f2", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
1020	f1.SetParameters(maxpar.GetArray());
1021	f2.SetParameters(maxpar.GetArray());
1022	f2.FixParameter(0, 0);
1023	f2.FixParameter(1, 0);
1024	f2.FixParameter(2, 1);
1025	f1.SetLineColor(kGreen);
1026	f2.SetLineColor(kRed);
1027
1028	f2.DrawCopy("same");
1029	f1.DrawCopy("same");
1030
1031	TPaveText *leg = new TPaveText(0.35, 0.10, 0.90, 0.35, "brNDC");
1032	leg->SetBorderSize(1);
1033	leg->SetTextSize(0.04);
1034	leg->AddText(0.5, 0.82, Form("A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + pol%d", polynom))->SetTextAlign(22);
1035	//leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22);
1036	leg->AddLine(0, 0.65, 0, 0.65);
1037	leg->AddText(0.06, 0.54, Form("A=%.2f", maxpar[0]))->SetTextAlign(11);
1038	leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", maxpar[2]))->SetTextAlign(11);
1039	leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fix)", maxpar[1]))->SetTextAlign(11);
1040	leg->AddText(0.60, 0.54, Form("a=%.2f", maxpar[3]))->SetTextAlign(11);
1041	leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11);
1042	if (polynom>1)
1043	leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11);
1044	leg->SetBit(kCanDelete);
1045	leg->Draw();
1046
1047	q->cd(2);
1048
1049	TGraph *g = new TGraph;
1050	g->SetPoint(0, 0, 0);
1051
1052	Int_t max=0;
1053	Float_t maxsig=0;
1054	for (int i=1; i<89; i++)
1055	{
1056	const Double_t s = f1.Integral(0, (float)i)/w;
1057	const Double_t b = f2.Integral(0, (float)i)/w;
1058
1059	const Double_t sig = SignificanceLiMa(s, b);
1060
1061	g->SetPoint(g->GetN(), i, sig);
1062
1063	if (sig>maxsig)
1064	{
1065	max = i;
1066	maxsig = sig;
1067	}
1068	}
1069
1070	g->SetNameTitle("SigVs\\alpha", "Significance vs \\alpha");
1071	g->GetHistogram()->SetXTitle("\\alpha_{0} [\\circ]");
1072	g->GetHistogram()->SetYTitle("Significance");
1073	g->SetBit(kCanDelete);
1074	g->Draw("AP");
1075
1076	leg = new TPaveText(0.35, 0.10, 0.90, 0.25, "brNDC");
1077	leg->SetBorderSize(1);
1078	leg->SetTextSize(0.1);
1079	leg->AddText(Form("\\sigma_{max}=%.1f at \\alpha_{max}=%d\\circ", maxsig, max));
1080	leg->SetBit(kCanDelete);
1081	leg->Draw();
1082	}
1083	}

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

Download in other formats: