Context Navigation

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

Visit:

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

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

Download in other formats: