Context Navigation

source: tags/Mars-V0.10.3/mhflux/MHFalseSource.cc

Visit:

Last change on this file was 8106, checked in by tbretz, 18 years ago
* empty log message *
File size: 38.3 KB

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

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

Download in other formats: