Context Navigation

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

Visit:

Last change on this file since 5672 was 5301, checked in by tbretz, 20 years ago
* empty log message *
File size: 33.2 KB

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

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

Download in other formats: