Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

source: trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc@ 5138

Visit:

Last change on this file since 5138 was 4830, checked in by mazin, 20 years ago
* empty log message *
File size: 37.4 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	! Author(s): Daniel Mazin, 8/2004 <mailto:mazin@mppmu.mpg.de>
20	!
21	! Copyright: MAGIC Software Development, 2000-2004
22	!
23	!
24	\* ======================================================================== */
25
26	//////////////////////////////////////////////////////////////////////////////
27	//
28	// MSkyPlot
29	//
30	// Create a false source plot. For the Binning in x,y the object MBinning
31	// "BinningFalseSource" is taken from the paremeter list.
32	//
33	// The binning in alpha is currently fixed as 18bins from 0 to 90deg.
34	//
35	// If MTime, MObservatory and MPointingPos is available in the paremeter
36	// list each image is derotated.
37	//
38	// MSkyPlot fills a 3D histogram with alpha distribution for
39	// each (derotated) x and y position.
40	//
41	// Only a radius of 1.2deg around the camera center is taken into account.
42	//
43	// The displayed histogram is the projection of alpha<fAlphaCut into
44	// the x,y plain and the projection of alpha>90-fAlphaCut
45	//
46	// By using the context menu (find it in a small region between the single
47	// pads) you can change the AlphaCut 'online'
48	//
49	// Each Z-Projection (Alpha-histogram) is scaled such, that the number
50	// of entries fits the maximum number of entries in all Z-Projections.
51	// This should correct for the different acceptance in the center
52	// and at the border of the camera. This, however, produces more noise
53	// at the border.
54	//
55	// Here is a slightly simplified version of the algorithm:
56	// ------------------------------------------------------
57	// MHillas hil; // Taken as second argument in MFillH
58	//
59	// MSrcPosCam src;
60	// MHillasSrc hsrc;
61	// hsrc.SetSrcPos(&src);
62	//
63	// for (int ix=0; ix<nx; ix++)
64	// for (int iy=0; iy<ny; iy++)
65	// {
66	// TVector2 v(cx[ix], cy[iy]); //cx center of x-bin ix
67	// if (rho!=0) //cy center of y-bin iy
68	// v=v.Rotate(rho); //image rotation angle
69	//
70	// src.SetXY(v); //source position in the camera
71	//
72	// if (!hsrc.Calc(hil)) //calc source dependant hillas
73	// return;
74	//
75	// //fill absolute alpha into histogram
76	// const Double_t alpha = hsrc.GetAlpha();
77	// fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w);
78	// }
79	// }
80	//
81	// Use MHFalse Source like this:
82	// -----------------------------
83	// MFillH fill("MSkyPlot", "MHillas");
84	// or
85	// MSkyPlot hist;
86	// hist.SetAlphaCut(12.5); // The default value
87	// hist.SetBgmean(55); // The default value
88	// MFillH fill(&hist, "MHillas");
89	//
90	// To be implemented:
91	// ------------------
92	// - a switch to use alpha or \|alpha\|
93	// - taking the binning for alpha from the parlist (binning is
94	// currently fixed)
95	// - a possibility to change the fit-function (eg using a different TF1)
96	// - a possibility to change the fit algorithm (eg which paremeters
97	// are fixed at which time)
98	// - use a different varaible than alpha
99	// - a possiblity to change the algorithm which is used to calculate
100	// alpha (or another parameter) is missing (this could be done using
101	// a tasklist somewhere...)
102	// - a more clever (and faster) algorithm to fill the histogram, eg by
103	// calculating alpha once and fill the two coils around the mean
104	// - more drawing options...
105	// - Move Significance() to a more 'general' place and implement
106	// also different algorithms like (Li/Ma)
107	// - implement fit for best alpha distribution -- online (Paint)
108	// - currently a constant pointing position is assumed in Fill()
109	// - the center of rotation need not to be the camera center
110	// - ERRORS on each alpha bin to estimate the chi^2 correctly!
111	// (sqrt(N)/binwidth) needed for WOlfgangs proposed caluclation
112	// of alpha(Li/Ma)
113	//
114	//////////////////////////////////////////////////////////////////////////////
115	#include "MSkyPlot.h"
116
117	#include <TF1.h>
118	#include <TH2.h>
119	#include <TH1.h>
120	#include <TGraph.h>
121	#include <TStyle.h>
122	#include <TLatex.h>
123	#include <TCanvas.h>
124	#include <TPaveText.h>
125	#include <TStopwatch.h>
126	#include <TFile.h>
127
128	#include "MGeomCam.h"
129	#include "MSrcPosCam.h"
130	#include "MReportDrive.h"
131	#include "MHillasSrc.h"
132	#include "MHillas.h"
133	#include "MHillasExt.h"
134	#include "MNewImagePar.h"
135	#include "MTime.h"
136	#include "MObservatory.h"
137	#include "MPointingPos.h"
138	#include "MAstroCatalog.h"
139	#include "MAstroSky2Local.h"
140	#include "MStarCamTrans.h"
141	#include "MStatusDisplay.h"
142	#include "MMath.h"
143	#include "MSupercuts.h"
144
145	#include "MBinning.h"
146	#include "MParList.h"
147
148	#include "MLog.h"
149	#include "MLogManip.h"
150
151	ClassImp(MSkyPlot);
152
153	using namespace std;
154
155	// --------------------------------------------------------------------------
156	//
157	// Default Constructor
158	//
159	MSkyPlot::MSkyPlot(const char name, const char title)
160	: fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1)
161	// fAlphaCut(12.5), BgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1)
162	{
163
164	*fLog << warn << "entering default constructor in MSkyPlot" << endl;
165	//
166	// set the name and title of this object
167	//
168	fName = name ? name : "MSkyPlot";
169	fTitle = title ? title : "sky plot vs x, y";
170
171	fSetCenter=kTRUE;
172
173	fHistSignif.SetDirectory(NULL);
174	fHistNexcess.SetDirectory(NULL);
175	fHistOn.SetDirectory(NULL);
176
177	fHistSignif.SetName("SkyPlotSignif");
178	fHistSignif.SetTitle("Sky Plot of significance vs x, y");
179	fHistSignif.SetXTitle("x [\\circ]");
180	fHistSignif.SetYTitle("y [\\circ]");
181
182	fHistNexcess.SetName("SkyPlotNexcess");
183	fHistNexcess.SetTitle("Sky Plot of number of excess vs x, y");
184	fHistNexcess.SetXTitle("x [\\circ]");
185	fHistNexcess.SetYTitle("y [\\circ]");
186
187	fHistOn.SetName("SkyPlotOn");
188	fHistOn.SetTitle("Sky Plot of number of On events vs x, y");
189	fHistOn.SetXTitle("x [\\circ]");
190	fHistOn.SetYTitle("y [\\circ]");
191
192	// set some values for the sky plot geometry:
193	fMinXGrid =-1.; // [deg]
194	fMaxXGrid = 1.; // [deg] , right edge of the skyplot
195	fMinYGrid =-1.; // [deg] , upper edge of the skyplot
196	fMaxYGrid = 1.; // [deg] , lower edge of the skyplot
197	fBinStepGrid = 0.05; // [deg],
198	fAlphaONMax = 5.; // [deg] , upper cut for alpha ON region in the alpha plot
199	// [deg], ON region in the alpha plot, maybe 5 deg is better
200	// NOTE: up to now only values of 5, 10, 15, 20 degrees are possible
201	// ra,dec lines from wolfgang:
202	fGridBinning = 0.50; // degrees
203	fGridFineBin = 0.01; // degrees
204
205
206	// some filter cuts:
207	fSizeMin = 2000.; // min size in photons
208	fSizeMax = 100000.; // max size in photons
209	fLeakMax = 0.25; // leakmax in per cent
210	fMaxDist = 1.4; // dist max cut (ever possible)
211	fMinDist = 0.1; // dist max cut (ever possible)
212
213	fNumBinsAlpha = 36; // number of bins for alpha histograms
214	fAlphaLeftEdge = 0.; // [deg] left edge
215	fAlphaRightEdge = 90.; // [deg] left edge
216
217	fAlphaBgLow = 30.;
218	fAlphaBgUp = 90.;
219
220	{
221	// SET DEFAULT VALUES HERE
222	fLengthUp[0] = 0.2;
223	fLengthUp[1] = 0.0;
224	fLengthUp[2] = 0.0;
225	fLengthUp[3] = 0.0;
226	fLengthUp[4] = 0.0;
227	fLengthUp[5] = 0.0;
228	fLengthUp[6] = 0.0;
229	fLengthUp[7] = 0.0;
230
231	fLengthLo[0] = 0.;
232	fLengthLo[1] = 0.;
233	fLengthLo[2] = 0.;
234	fLengthLo[3] = 0.;
235	fLengthLo[4] = 0.;
236	fLengthLo[5] = 0.;
237	fLengthLo[6] = 0.;
238	fLengthLo[7] = 0.;
239
240	fWidthUp[0] = 0.1;
241	fWidthUp[1] = 0.0;
242	fWidthUp[2] = 0.0;
243	fWidthUp[3] = 0.0;
244	fWidthUp[4] = 0.0;
245	fWidthUp[5] = 0.0;
246	fWidthUp[6] = 0.0;
247	fWidthUp[7] = 0.0;
248
249	fWidthLo[0] = 0.;
250	fWidthLo[1] = 0.;
251	fWidthLo[2] = 0.;
252	fWidthLo[3] = 0.;
253	fWidthLo[4] = 0.;
254	fWidthLo[5] = 0.;
255	fWidthLo[6] = 0.;
256	fWidthLo[7] = 0.;
257
258	fDistUp[0] = 1.e10;
259	fDistUp[1] = 0.0;
260	fDistUp[2] = 0.0;
261	fDistUp[3] = 0.0;
262	fDistUp[4] = 0.0;
263	fDistUp[5] = 0.0;
264	fDistUp[6] = 0.0;
265	fDistUp[7] = 0.0;
266
267	fDistLo[0] = 0.0;
268	fDistLo[1] = 0.0;
269	fDistLo[2] = 0.0;
270	fDistLo[3] = 0.0;
271	fDistLo[4] = 0.0;
272	fDistLo[5] = 0.0;
273	fDistLo[6] = 0.0;
274	fDistLo[7] = 0.0;
275	}
276	}
277
278	void MSkyPlot::ReadCuts(const TString parSCinit="OptSCParametersONOFFThetaRange0_1570mRad.root")
279	{
280	//--------------------------------
281	// create container for the supercut parameters
282	// and set them to their initial values
283	MSupercuts super;
284
285	// read the cuts coefficients
286	TFile inparam(parSCinit);
287	super.Read("MSupercuts");
288	inparam.Close();
289	*fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from file "
290	<< parSCinit << endl;
291
292	TArrayD params = super.GetParameters();
293	TArrayD steps = super.GetStepsizes();
294	// TMP2
295	if (params.GetSize() == steps.GetSize())
296	{
297	*fLog << "SC parameters and Setps are: " << endl;
298	for (Int_t z = 0; z < params.GetSize(); z++)
299	{
300	cout << params[z] << setw(20) << steps[z] << endl;
301	}
302	}
303	// ENDTMP2
304	for (Int_t i=0; i<8; i++)
305	{
306	fLengthUp[i] = params[i];
307	fLengthLo[i] = params[i+8];
308	fWidthUp[i] = params[i+16];
309	fWidthLo[i] = params[i+24];
310	fDistUp[i] = params[i+32];
311	fDistLo[i] = params[i+40];
312	}
313	}
314
315	void MSkyPlot::SetSkyPlot(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t step)
316	{
317	Float_t temp;
318
319	if (xmax<xmin)
320	{
321	*fLog << warn << "SetSkyPlot: xmax is smaller xmin ... exchanging them." << endl;
322	temp = xmax;
323	xmax = xmin;
324	xmin = temp;
325	}
326
327	if (ymax<ymin)
328	{
329	*fLog << warn << "SetSkyPlot: ymax is smaller ymin ... exchanging them." << endl;
330	temp = ymax;
331	ymax = ymin;
332	ymin = temp;
333	}
334
335	if (step<0)
336	*fLog << warn << "SetSkyPlot: step<0... taking absolute value." << endl;
337
338	fBinStepGrid = TMath::Abs(step);
339	fMinXGrid = xmin;
340	fMaxXGrid = xmax;
341	fMinYGrid = ymin;
342	fMaxYGrid = ymax;
343	}
344
345	// --------------------------------------------------------------------------
346	//
347	// Set the alpha cut (\|alpha\|<fAlphaCut) which is use to estimate the
348	// number of excess events
349	//
350	void MSkyPlot::SetAlphaCut(Float_t alpha)
351	{
352	if (alpha<0)
353	*fLog << warn << "Alpha<0... taking absolute value." << endl;
354
355	fAlphaONMax = TMath::Abs(alpha);
356	}
357
358	// calculate the values for the cuts:
359	Double_t MSkyPlot::CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2)
360	{
361
362	Double_t limit = a[0] + a[1] * dd2 +
363	ls * (a[3] + a[4] * dd2) +
364	ls2 * (a[6] + a[7] * dd2);
365
366	return limit;
367	}
368
369
370
371	// --------------------------------------------------------------------------
372	//
373	// Set binnings (takes BinningFalseSource) and prepare filling of the
374	// histogram.
375	//
376	// Also search for MTime, MObservatory and MPointingPos
377	//
378	Int_t MSkyPlot::PreProcess(MParList *plist)
379	{
380	*fLog << warn << "entering PreProcess in MSkyPlot::PreProcess" << endl;
381	fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");
382	if (!fGeomCam)
383	{
384	*fLog << err << "MGeomCam not found... aborting." << endl;
385	return kFALSE;
386	}
387
388	fMm2Deg = fGeomCam->GetConvMm2Deg();
389
390	fRepDrive = (MReportDrive*)plist->FindObject(AddSerialNumber("MReportDrive"));
391	if (!fRepDrive)
392	*fLog << warn << "MReportDrive not found... no sky map." << endl;
393
394
395	fSrcPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam"));
396	if (!fSrcPosCam)
397	*fLog << warn << "MSrcPosCam not found... no sky map." << endl;
398
399	fPntPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MPntPosCam"));
400	if (!fPntPosCam)
401	*fLog << warn << "MPntPosCam not found... no sky map." << endl;
402
403	fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos"));
404	if (!fPointPos)
405	*fLog << warn << "MPointingPos not found... no sky map." << endl;
406
407	fTime = (MTime*)plist->FindObject(AddSerialNumber("MTime"));
408	if (!fTime)
409	*fLog << warn << "MTime not found... no sky map." << endl;
410
411	fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory"));
412	if (!fObservatory)
413	*fLog << warn << "MObservatory not found... no sky map." << endl;
414
415
416	fHillas = (MHillas*)plist->FindObject(AddSerialNumber("MHillas"));
417	if (!fHillas)
418	*fLog << err << "MHillas not found... no sky map." << endl;
419
420	fHillasExt = (MHillasExt*)plist->FindObject(AddSerialNumber("MHillasExt"));
421	if (!fHillasExt)
422	*fLog << err << "MHillasExt not found... no sky map." << endl;
423
424	fHillasSrc = (MHillasSrc*)plist->FindObject(AddSerialNumber("MHillasSrc"));
425	if (!fHillasSrc)
426	*fLog << err << "MHillasSrc not found... no sky map." << endl;
427
428	fNewImagePar = (MNewImagePar*)plist->FindObject(AddSerialNumber("MNewImagePar"));
429	if (!fNewImagePar)
430	*fLog << err << "MNewImagePar not found... no sky map." << endl;
431
432
433
434	// FIXME: Because the pointing position could change we must check
435	// for the current pointing position and add a offset in the
436	// Fill function!
437	kSaveAlphaPlots=kTRUE;
438	kSaveSkyPlots=kTRUE;
439
440	// prepare skyplot
441	fNumStepsX = (int) ((fMaxXGrid - fMinXGrid) / fBinStepGrid + 1.5);
442	fNumStepsY = (int) ((fMaxYGrid - fMinYGrid) / fBinStepGrid + 1.5);
443	fNumalphahist = (int) (fNumStepsX * fNumStepsY + 0.5);
444
445	fHistSignif.SetName("SPSignif2ndOrder");
446	fHistSignif.SetTitle("Sky Plot of significance (2nd order fit)");
447	fHistSignif.SetBins(fNumStepsX, fMinXGrid-0.5fBinStepGrid, fMaxXGrid+0.5fBinStepGrid,
448	fNumStepsY, fMinYGrid-0.5fBinStepGrid, fMaxYGrid+0.5fBinStepGrid);
449	fHistSignif.SetDirectory(NULL);
450	fHistSignif.SetFillStyle(4000);
451	fHistSignif.UseCurrentStyle();
452
453	fHistNexcess.SetName("SPNex2ndOrder");
454	fHistNexcess.SetTitle("Sky Plot of Number of excess events (2nd order fit)");
455	fHistNexcess.SetBins(fNumStepsX, fMinXGrid-0.5fBinStepGrid, fMaxXGrid+0.5fBinStepGrid,
456	fNumStepsY, fMinYGrid-0.5fBinStepGrid, fMaxYGrid+0.5fBinStepGrid);
457	fHistNexcess.SetDirectory(NULL);
458	fHistNexcess.SetFillStyle(4000);
459	fHistNexcess.UseCurrentStyle();
460
461	fHistOn.SetName("SPOnCounted");
462	fHistOn.SetTitle("Sky Plot of ON events (counted)");
463	fHistOn.SetBins(fNumStepsX, fMinXGrid-0.5fBinStepGrid, fMaxXGrid+0.5fBinStepGrid,
464	fNumStepsY, fMinYGrid-0.5fBinStepGrid, fMaxYGrid+0.5fBinStepGrid);
465	fHistOn.SetDirectory(NULL);
466	fHistOn.SetFillStyle(4000);
467	fHistOn.UseCurrentStyle();
468
469
470	// prepare alpha plots
471	// temp histogram for an alpha plot
472	TH1I *temp = new TH1I("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge);
473	temp->SetDirectory(NULL);
474	fHistAlphaBinWidth = temp->GetBinWidth(1);
475	// create alpha histograms
476	for (Int_t i=0; i< fNumalphahist; i++) {
477	fHistAlpha.push_back(*temp);
478	fHistAlpha[i].SetDirectory(NULL);
479	}
480
481	delete temp;
482	return kTRUE;
483	}
484
485	// --------------------------------------------------------------------------
486	//
487	// Fill the histogram. For details see the code or the class description
488	//
489	Int_t MSkyPlot::Process()
490	{
491	// check whether MPointingPos comtains something:
492	if (TMath::Abs(fPointPos->GetRa())<1e-3 && TMath::Abs(fPointPos->GetDec())<1e-3)
493	{
494	*fLog << warn << "MPointingPos is not filled ... event skipped" << endl;
495	return kTRUE;
496	}
497
498	// Get RA_0, DEC_0 for the camera center (Tracking MDrive?, can be set from outside)
499	if (fSetCenter==kTRUE)
500	{
501	fRa0 = fRepDrive->GetRa(); // [h]
502	fDec0 = fRepDrive->GetDec(); // [deg]
503	fSetCenter=kFALSE;
504	}
505
506
507	// some filter cuts:
508	if ( fHillas->GetSize() > fSizeMin && fHillas->GetSize() < fSizeMax && fNewImagePar->GetLeakage1() < fLeakMax)
509	{
510
511	Double_t xsource, ysource, cosgam, singam, x_0, y_0, xsourcam, ysourcam;
512	Double_t dx, dy, beta, tanbeta, alphapar, distpar;
513	Double_t logsize, lgsize, lgsize2, dist2;
514	const Double_t log3000 = log(3000.);
515	const Float_t fDistOffset = 0.9;
516
517	// Get camera rotation angle
518	Double_t rho = 0;
519	if (fTime && fObservatory && fPointPos)
520	{
521	rho = fPointPos->RotationAngle(fObservatory, fTime);
522	*fLog << inf << " rho = " << rho << ", Zd = " << fPointPos->GetZd() <<
523	", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl;
524	}
525	//if (fPointPos)
526	// rho = fPointPos->RotationAngle(*fObservatory);
527	/*
528	//TEMP
529	// theta = mcevt->GetTelescopeTheta();
530	Double_t theta, phi, sing, cosg;
531	theta = fPointPos->GetZd()*TMath::Pi()/180.;
532	phi = fPointPos->GetAz()*TMath::Pi()/180.;
533	// printf("theta: %5.3f, phi: %5.3f\n", theta180./4.1415, phi180./4.1415);
534	fObservatory->RotationAngle(theta, phi, sing, cosg);
535
536	// conclusion: diffference in rho = 7 deg
537	*fLog << "new thetaTel, phiTel, cosal, sinal, rho = " << theta << ", "
538	<< phi << ", " << cosg << ", " << sing << ", " << TMath::ATan2(sing,cosg)*180./TMath::Pi() << endl;
539
540	Double_t a1 = 0.876627;
541	Double_t a3 = -0.481171;
542	Double_t thetaTel=theta, phiTel=phi;
543
544	Double_t denom = 1./ sqrt( sin(thetaTel)sin(phiTel) sin(thetaTel)*sin(phiTel) +
545	( a1cos(thetaTel)+a3sin(thetaTel)cos(phiTel) )
546	( a1cos(thetaTel)+a3sin(thetaTel)*cos(phiTel) ) );
547	Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom;
548	Double_t sinal = a1 * sin(phiTel) * denom;
549
550	*fLog << "old thetaTel, phiTel, cosal, sinal, rho = " << thetaTel << ", "
551	<< phiTel << ", " << cosal << ", " << sinal << ", " << TMath::ATan2(sinal,cosal)*180./TMath::Pi() << endl;
552
553
554
555	// END TEMP
556	*/
557
558	// => coord. system: xtilde, ytilde with center in RA_0, DEC_0
559
560	// Get Rot. Angle:
561	cosgam = TMath::Cos(rho);
562	singam = TMath::Sin(rho);
563	// Get x_0, y_0 for RA_0,DEC_0 of the current event
564	x_0 = fPntPosCam->GetX()*fMm2Deg;
565	y_0 = fPntPosCam->GetY()*fMm2Deg;
566
567	Int_t index = 0; // index for alpha histograms
568	// loop over xtilde
569	for (Int_t gridy = 0; gridy < fNumStepsY; gridy++)
570	{
571	ysource = fMinYGrid + gridy * fBinStepGrid;
572	// loop over ytilde
573	for (Int_t gridx = 0; gridx < fNumStepsX; gridx++)
574	{
575	xsource = fMinXGrid + gridx * fBinStepGrid;
576	/* derotation : rotate into camera coordinates */
577	/* formel: (x_cam) (cos(gam) -sin(gam)) (xtilde) (x_0)
578	( ) = ( ) * ( ) + ( )
579	(y_cam) (sin(gam) sin(gam)) (ytilde) (y_0)
580	*/
581	xsourcam = cosgam * xsource - singam * ysource + x_0;
582	ysourcam = singam * xsource + cosgam * ysource + y_0;
583	/* end derotatiom */
584
585	/* calculate ALPHA und DIST according to the source position */
586	dx = fHillas->GetMeanX()*fMm2Deg - xsourcam;
587	dy = fHillas->GetMeanY()*fMm2Deg - ysourcam;
588	tanbeta = dy / dx ;
589	beta = TMath::ATan(tanbeta);
590	alphapar = (fHillas->GetDelta() - beta) * 180./ TMath::Pi();
591	distpar = sqrt( dydy + dxdx );
592	if(alphapar > 90.) alphapar -= 180.;
593	if(alphapar < -90.) alphapar += 180.;
594	alphapar = TMath::Abs(alphapar);
595
596	// apply cuts
597	logsize = log(fHillas->GetSize());
598	lgsize = logsize-log3000;
599	lgsize2 = lgsize*lgsize;
600	dist2 = (distpar-fDistOffset)*(distpar-fDistOffset);
601
602	if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) &&
603	(fHillas->GetLength()*fMm2Deg) > CalcLimit(fLengthLo, lgsize, lgsize2, dist2))
604	if ( (fHillas->GetWidth()*fMm2Deg) < CalcLimit(fWidthUp, lgsize, lgsize2, dist2) &&
605	(fHillas->GetWidth()*fMm2Deg) > CalcLimit(fWidthLo, lgsize, lgsize2, dist2))
606	if ( distpar < CalcLimit(fDistUp, lgsize, lgsize2, dist2) &&
607	distpar > CalcLimit(fDistLo, lgsize, lgsize2, dist2) &&
608	distpar < fMaxDist && distpar > fMinDist)
609	{
610	// gamma candidates!
611	fHistAlpha[index].Fill(alphapar);
612	index++;
613	}
614	}
615	}
616	}
617
618	return kTRUE;
619	}
620
621	// calculate number of events below alphacut, number of excess events, significance
622
623	Int_t MSkyPlot::PostProcess()
624	{
625
626	Int_t nrow, ncolumn;
627	Double_t Non, chisquarefit, numdegfreed, Noff_fit, Nex, Sign;
628	const Int_t onbinsalpha = TMath::Nint(fAlphaONMax/fHistAlphaBinWidth);
629
630	// fit function for the background
631	TF1 * fitbgpar = new TF1("fbgpar", "[0]xx + [1]", fAlphaBgLow, fAlphaBgUp);
632	fitbgpar->SetLineColor(2);
633
634	// number degrees of freedom:
635	numdegfreed = (fAlphaBgUp - fAlphaBgLow) / fHistAlphaBinWidth - 2.;
636
637	vector <TH1I>::iterator alpha_iterator;
638	int index2 = 0; // index of the TH2F histograms
639
640	alpha_iterator = fHistAlpha.begin();
641
642	while( alpha_iterator != fHistAlpha.end() ) {
643	// find the bin in the 2dim histogram
644	nrow = index2/fHistOn.GetNbinsX() + 1; // row of the histogram (y)
645	ncolumn = TMath::Nint(fmod(index2,fHistOn.GetNbinsX()))+1; // column of the histogram (x)
646
647	// number of ON events below fAlphaONMax
648	Non = 0.;
649	for(Int_t i=1; i<=onbinsalpha;i++) Non += (*alpha_iterator).GetBinContent(i);
650
651	fHistOn.SetBinContent(ncolumn, nrow, Non); // fill in the fHistOn
652
653	// fit parabola to a background region
654	(*alpha_iterator).Fit(fitbgpar,"NWR"); // NWR OK?????????????????????????
655	// get Chi2
656	chisquarefit = fitbgpar->GetChisquare();
657	if (chisquarefit/numdegfreed<2. && chisquarefit/numdegfreed>0.01);
658	else *fLog << warn << "Fit is bad, chi2/ndf = " << chisquarefit/numdegfreed << endl;
659
660	// estimate Noff from the fit:
661	Noff_fit = (1./3. * (fitbgpar->GetParameter(0)) * pow(fAlphaONMax,3.) +
662	(fitbgpar->GetParameter(1)) * fAlphaONMax) / fHistAlphaBinWidth;
663
664	Nex = Non - Noff_fit;
665
666	fHistNexcess.SetBinContent(ncolumn, nrow, Nex); // fill in the fHistNexcess
667
668	if (Noff_fit<0.) Sign = 0.;
669	// else Sign = LiMa17(Non,Noff_fit,1.);
670	else Sign = MMath::SignificanceLiMaSigned(Non, Noff_fit, 1.);
671
672	fHistSignif.SetBinContent(ncolumn, nrow, Sign); // fill in the fHistSignif
673
674	alpha_iterator++;
675	index2++;
676	}
677
678	//temp
679	Double_t decl, hang;
680	MStarCamTrans mstarc(fGeomCam, fObservatory);
681	mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),decl, hang);
682
683	*fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl;
684	*fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl;
685	//endtemp
686
687	// save alpha plots:
688	if(kSaveAlphaPlots==kTRUE) SaveAlphaPlots();
689
690	// save sky plots:
691	if(kSaveSkyPlots==kTRUE) SaveSkyPlots();
692
693	fHistAlpha.clear();
694	return kTRUE;
695	}
696
697
698	// --------------------------------------------------------------------------
699	//
700	// Get the MAstroCatalog corresponding to fRa, fDec. The limiting magnitude
701	// is set to 9, while the fov is equal to the current fov of the false
702	// source plot.
703	//
704	TObject *MSkyPlot::GetCatalog()
705	{
706	const Double_t maxr = 0.98*TMath::Abs(fHistSignif.GetBinCenter(1));
707
708	// Create catalog...
709	MAstroCatalog stars;
710	stars.SetLimMag(9);
711	stars.SetGuiActive(kFALSE);
712	stars.SetRadiusFOV(maxr);
713	stars.SetRaDec(fRaTMath::DegToRad()15, fDec*TMath::DegToRad());
714	stars.ReadBSC("bsc5.dat");
715
716	TObject o = (MAstroCatalog)stars.Clone();
717	o->SetBit(kCanDelete);
718
719	return o;
720	}
721
722
723	void MSkyPlot::SaveSkyPlots(const TString skyplotfilename="skyplots.root")
724	{
725	TFile rootfile(skyplotfilename, "RECREATE",
726	"sky plots in some variations");
727	fHistSignif.Write();
728	fHistNexcess.Write();
729	fHistOn.Write();
730
731	// from Wolfgang:
732	//--------------------------------------------------------------------
733	// Dec-Hour lines
734	Double_t rho, sinrho, cosrho;
735	// if (fTime && fObservatory && fPointPos)
736	rho = fPointPos->RotationAngle(fObservatory, fTime);
737	sinrho=TMath::Sin(rho);
738	cosrho=TMath::Cos(rho);
739
740	Double_t theta, phi, sing, cosg;
741	theta = fPointPos->GetZd()*TMath::Pi()/180.;
742	phi = fPointPos->GetAz()*TMath::Pi()/180.;
743	// printf("theta: %5.3f, phi: %5.3f\n", theta180./4.1415, phi180./4.1415);
744	fObservatory->RotationAngle(theta, phi, sing, cosg);
745
746	*fLog << "1: sinrho, cosrho = " << sinrho << ", " << cosrho << endl;
747	*fLog << "2: sinrho, cosrho = " << sing << ", " << cosg << endl;
748
749	Double_t fDistCam = fGeomCam->GetCameraDist() * 1000.0; // [mm]
750	Double_t gridbinning = fGridBinning;
751	Double_t gridfinebin = fGridFineBin;
752	Int_t numgridlines = (Int_t)(4.0/gridbinning);
753	Double_t aberr = 1.07;
754	Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ;
755
756	Double_t declin, hangle; // [deg], [h]
757	MStarCamTrans mstarc(fGeomCam, fObservatory);
758	mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle);
759
760	TGraph *graph1;
761	TLatex *pix;
762
763	Char_t tit[100];
764	Double_t xtxt;
765	Double_t ytxt;
766
767	Double_t xlo = -534.0 * mmtodeg;
768	Double_t xup = 534.0 * mmtodeg;
769
770	Double_t ylo = -534.0 * mmtodeg;
771	Double_t yup = 534.0 * mmtodeg;
772
773	Double_t xx, yy;
774
775
776	// direction for the center of the camera
777	Double_t dec0 = declin;
778	Double_t h0 = hangle*360./24.;
779	//trans.LocToCel(theta0, phi0, dec0, h0);
780
781	TCanvas *c1 = new TCanvas("SPlotsRaDecLines","SPlotsRaDecLines", 400, 0, 700, 600);
782	c1->cd(1);
783	fHistOn.Draw("colz");
784
785	//----- lines for fixed dec ------------------------------------
786
787	const Int_t Ndec = numgridlines;
788	Double_t dec[Ndec];
789	Double_t ddec = gridbinning;
790
791	Int_t nh = (Int_t)(4.0/gridfinebin);
792	const Int_t Nh = nh+1;
793	Double_t h[Nh];
794	Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926);
795	if ( dh > 360.0/((Double_t)(Nh-1)) )
796	dh = 360.0/((Double_t)(Nh-1));
797
798
799	for (Int_t j=0; j<Ndec; j++)
800	{
801	dec[j] = dec0 + ((Double_t)j)*ddec
802	- ((Double_t)(Ndec/2)-1.0)*ddec;
803	}
804
805	for (Int_t k=0; k<Nh; k++)
806	{
807	h[k] = h0 + ((Double_t)k)dh - ((Double_t)(Nh/2)-1.0)dh;
808	}
809
810	Double_t xh[Nh];
811	Double_t yh[Nh];
812
813
814
815	for (Int_t j=0; j<Ndec; j++)
816	{
817	if (fabs(dec[j]) > 90.0) continue;
818
819	for (Int_t k=0; k<Nh; k++)
820	{
821	Double_t hh0 = h0 *24.0/360.0;
822	Double_t hx = h[k]*24.0/360.0;
823	mstarc.Cel0CelToCam(dec0, hh0, dec[j], hx,
824	xx, yy);
825	xx = xx * mmtodeg * aberr;
826	yy = yy * mmtodeg * aberr;
827	// xh[k] = xx * mmtodeg * aberr;
828	// yh[k] = yy * mmtodeg * aberr;
829	xh[k] = cosg * xx + sing * yy;
830	yh[k] =-sing * xx + cosg * yy;
831	// xh[k] = cosrho * xx + sinrho * yy;
832	// yh[k] =-sinrho * xx + cosrho * yy;
833
834
835	//gLog << "dec0, h0 = " << dec0 << ", " << h0
836	// << " : dec, h, xh, yh = " << dec[j] << ", "
837	// << h[k] << "; " << xh[k] << ", " << yh[k] << endl;
838	}
839
840	// c1->cd(2);
841	graph1 = new TGraph(Nh,xh,yh);
842	graph1->SetLineColor(j+1);
843	graph1->SetLineStyle(1);
844	graph1->SetMarkerColor(j+1);
845	graph1->SetMarkerSize(.2);
846	graph1->SetMarkerStyle(20);
847	graph1->Draw("PL");
848	//delete graph1;
849
850	sprintf(tit,"Dec = %6.2f", dec[j]);
851	xtxt = xlo + (xup-xlo)*0.1;
852	ytxt = ylo + (yup-ylo)0.80 - ((Double_t)j) (yup-ylo)/20.0;
853	pix = new TLatex(xtxt, ytxt, tit);
854	pix->SetTextColor(j+1);
855	pix->SetTextSize(.03);
856	pix->Draw("");
857	//delete pix;
858
859	}
860
861	//----- lines for fixed hour angle ------------------------------------
862
863	Int_t ndec1 = (Int_t)(4.0/gridfinebin);
864	const Int_t Ndec1 = ndec1;
865	Double_t dec1[Ndec1];
866	Double_t ddec1 = gridfinebin;
867
868	const Int_t Nh1 = numgridlines;
869	Double_t h1[Nh1];
870	Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926);
871	if ( dh1 > 360.0/((Double_t)Nh1) )
872	dh1 = 360.0/((Double_t)Nh1);
873
874
875	for (Int_t j=0; j<Ndec1; j++)
876	{
877	dec1[j] = dec0 + ((Double_t)j)*ddec1
878	- ((Double_t)(Ndec1/2)-1.0)*ddec1;
879	}
880
881	for (Int_t k=0; k<Nh1; k++)
882	{
883	h1[k] = h0 + ((Double_t)k)dh1 - ((Double_t)(Nh1/2)-1.0)dh1;
884	}
885
886
887	Double_t xd[Ndec1];
888	Double_t yd[Ndec1];
889
890	for (Int_t k=0; k<Nh1; k++)
891	{
892	Int_t count = 0;
893	for (Int_t j=0; j<Ndec1; j++)
894	{
895	if (fabs(dec1[j]) > 90.0) continue;
896
897	Double_t hh0 = h0 *24.0/360.0;
898	Double_t hhx = h1[k]*24.0/360.0;
899	mstarc.Cel0CelToCam(dec0, hh0, dec1[j], hhx,
900	xx, yy);
901	// xd[count] = xx * mmtodeg * aberr;
902	// yd[count] = yy * mmtodeg * aberr;
903
904	xx = xx * mmtodeg * aberr;
905	yy = yy * mmtodeg * aberr;
906	xd[count] = cosg * xx + sing * yy;
907	yd[count] =-sing * xx + cosg * yy;
908
909	//gLog << "dec0, h0 = " << dec0 << ", " << h0
910	// << " : dec1, h1, xd, yd = " << dec1[j] << ", "
911	// << h1[k] << "; " << xd[count] << ", " << yd[count] << endl;
912
913	count++;
914	}
915
916	// c1->cd(2);
917	graph1 = new TGraph(count,xd,yd);
918	graph1->SetLineColor(k+1);
919	graph1->SetLineStyle(2);
920	graph1->SetMarkerColor(k+1);
921	graph1->SetMarkerSize(.2);
922	graph1->SetMarkerStyle(20);
923	graph1->Draw("PL");
924	//delete graph1;
925
926	sprintf(tit,"H = %6.2f", h1[k]);
927	Double_t xtxt = xlo + (xup-xlo)*0.75;
928	Double_t ytxt = ylo + (yup-ylo)0.80 - ((Double_t)k) (yup-ylo)/20.0;
929	pix = new TLatex(xtxt, ytxt, tit);
930	pix->SetTextColor(k+1);
931	pix->SetTextSize(.03);
932	pix->Draw("");
933	//delete pix;
934
935	}
936
937	// c1->cd(2);
938	sprintf(tit,"Dec0 = %6.2f H0 = %6.2f", dec0, h0);
939	xtxt = xlo + (xup-xlo)*0.05;
940	ytxt = ylo + (yup-ylo)*0.92;
941	pix = new TLatex(xtxt, ytxt, tit);
942	pix->SetTextColor(1);
943	pix->SetTextSize(.06);
944	pix->Draw("");
945	//delete pix;
946
947	sprintf(tit," [deg] [deg]");
948	xtxt = xlo + (xup-xlo)*0.05;
949	ytxt = ylo + (yup-ylo)*0.86;
950	pix = new TLatex(xtxt, ytxt, tit);
951	pix->SetTextColor(1);
952	pix->SetTextSize(.06);
953	pix->Draw("");
954
955	c1->Write();
956
957	TCanvas *c2 = new TCanvas("SkyPlotsWithRaDecLines","SkyPlotsWithRaDecLines", 0, 0, 300, 600);
958	c2->Divide(1,2);
959	c2->SetBorderMode(0);
960	c2->cd(1);
961	fHistSignif.Draw("colz");
962	c2->cd(2);
963	fHistNexcess.Draw("colz");
964	c2->Write();
965
966	rootfile.Close();
967	delete c1;
968	delete c2;
969	delete graph1;
970	delete pix;
971	}
972
973	void MSkyPlot::SaveAlphaPlots(const TString alphaplotfilename="alphaplots.root")
974	{
975	TFile rootfile(alphaplotfilename, "RECREATE",
976	"all the alpha plots");
977
978	vector <TH1I>::iterator alpha_iterator;
979	int index = 0; // index of the TH2F histograms
980	Char_t strtitle[100];
981	Char_t strname[100];
982	Float_t xpos, ypos, signif, nex;
983	Int_t nrow, ncolumn, non;
984
985
986	alpha_iterator = fHistAlpha.begin();
987
988	while( alpha_iterator != fHistAlpha.end() ) {
989
990	nrow = index/fHistOn.GetNbinsX() + 1; // row of the histogram (y)
991	ncolumn = TMath::Nint(fmod(index,fHistOn.GetNbinsX()))+1; // column of the histogram (x)
992	xpos = fMinXGrid + fBinStepGrid*(ncolumn-1);
993	ypos = fMinYGrid + fBinStepGrid*(nrow-1);
994	non = TMath::Nint(fHistOn.GetBinContent(ncolumn,nrow));
995	nex = fHistNexcess.GetBinContent(ncolumn,nrow);
996	signif = fHistSignif.GetBinContent(ncolumn,nrow);
997
998	sprintf(strname,"AlphaPlotX%.2fY%.2f", xpos, ypos);
999	sprintf(strtitle,"for x= %.2f deg, y= %.2f deg: S= %.2f sigma, Nex= %.2f, Non= %d",
1000	xpos, ypos, signif, nex, non);
1001
1002	(*alpha_iterator).SetName(strname);
1003	(*alpha_iterator).SetTitle(strtitle);
1004	(*alpha_iterator).Write();
1005	alpha_iterator++;
1006	index++;
1007	}
1008
1009	rootfile.Close();
1010	}
1011
1012	// --------------------------------------------------------------------------
1013	//
1014	// Draw the histogram
1015	//
1016	void MSkyPlot::Draw(Option_t *opt)
1017	{
1018	/*
1019	TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
1020	pad->SetBorderMode(0);
1021
1022	AppendPad("");
1023
1024	pad->Divide(1, 2, 0, 0.03);
1025
1026	TObject *catalog = GetCatalog();
1027
1028	// Initialize upper part
1029	pad->cd(1);
1030	gPad->SetBorderMode(0);
1031	gPad->Divide(3, 1);
1032	*/
1033	/*
1034	// PAD #1
1035	pad->GetPad(1)->cd(1);
1036	gPad->SetBorderMode(0);
1037	fHist.GetZaxis()->SetRangeUser(0,fAlphaONMax);
1038	TH1 *h3 = fHist.Project3D("yx_on");
1039	fHist.GetZaxis()->SetRange(0,9999);
1040	h3->SetDirectory(NULL);
1041	h3->SetXTitle(fHist.GetXaxis()->GetTitle());
1042	h3->SetYTitle(fHist.GetYaxis()->GetTitle());
1043	h3->Draw("colz");
1044	h3->SetBit(kCanDelete);
1045	catalog->Draw("mirror same");
1046
1047	// PAD #2
1048	pad->GetPad(1)->cd(2);
1049	gPad->SetBorderMode(0);
1050	fHist.GetZaxis()->SetRange(0,0);
1051	TH1 *h4 = fHist.Project3D("yx_sig"); // Do this to get the correct binning....
1052	fHist.GetZaxis()->SetRange(0,9999);
1053	h4->SetTitle("Significance");
1054	h4->SetDirectory(NULL);
1055	h4->SetXTitle(fHist.GetXaxis()->GetTitle());
1056	h4->SetYTitle(fHist.GetYaxis()->GetTitle());
1057	h4->Draw("colz");
1058	h4->SetBit(kCanDelete);
1059	catalog->Draw("mirror same");
1060
1061	// PAD #3
1062	pad->GetPad(1)->cd(3);
1063	gPad->SetBorderMode(0);
1064	fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
1065	TH1 *h2 = fHist.Project3D("yx_off");
1066	h2->SetDirectory(NULL);
1067	h2->SetXTitle(fHist.GetXaxis()->GetTitle());
1068	h2->SetYTitle(fHist.GetYaxis()->GetTitle());
1069	h2->Draw("colz");
1070	h2->SetBit(kCanDelete);
1071	catalog->Draw("mirror same");
1072
1073
1074	// Initialize lower part
1075	pad->cd(2);
1076	gPad->SetBorderMode(0);
1077	gPad->Divide(3, 1);
1078
1079	// PAD #4
1080	pad->GetPad(2)->cd(1);
1081	gPad->SetBorderMode(0);
1082	TH1 *h1 = fHist.ProjectionZ("Alpha");
1083	h1->SetDirectory(NULL);
1084	h1->SetTitle("Distribution of \\alpha");
1085	h1->SetXTitle(fHist.GetZaxis()->GetTitle());
1086	h1->SetYTitle("Counts");
1087	h1->Draw(opt);
1088	h1->SetBit(kCanDelete);
1089
1090	// PAD #5
1091	pad->GetPad(2)->cd(2);
1092	gPad->SetBorderMode(0);
1093	TH1 h5 = (TH1)h3->Clone("Alpha_yx_diff");
1094	h5->Add(h2, -1);
1095	h5->SetTitle("Difference of on- and off-distribution");
1096	h5->SetDirectory(NULL);
1097	h5->Draw("colz");
1098	h5->SetBit(kCanDelete);
1099	catalog->Draw("mirror same");
1100
1101	// PAD #6
1102	pad->GetPad(2)->cd(3);
1103	gPad->SetBorderMode(0);
1104	TH1 *h0 = fHist.Project3D("yx_all");
1105	h0->SetDirectory(NULL);
1106	h0->SetXTitle(fHist.GetXaxis()->GetTitle());
1107	h0->SetYTitle(fHist.GetYaxis()->GetTitle());
1108	h0->Draw("colz");
1109	h0->SetBit(kCanDelete);
1110	catalog->Draw("mirror same");
1111	*/
1112	}
1113
1114	// --------------------------------------------------------------------------
1115	//
1116	// Everything which is in the main pad belongs to this class!
1117	//
1118	Int_t MSkyPlot::DistancetoPrimitive(Int_t px, Int_t py)
1119	{
1120	return 0;
1121	}
1122
1123	// --------------------------------------------------------------------------
1124	//
1125	// Set all sub-pads to Modified()
1126	//
1127	/*
1128	void MSkyPlot::Modified()
1129	{
1130	if (!gPad)
1131	return;
1132
1133	TVirtualPad *padsave = gPad;
1134	padsave->Modified();
1135	padsave->GetPad(1)->cd(1);
1136	gPad->Modified();
1137	padsave->GetPad(1)->cd(3);
1138	gPad->Modified();
1139	padsave->GetPad(2)->cd(1);
1140	gPad->Modified();
1141	padsave->GetPad(2)->cd(2);
1142	gPad->Modified();
1143	padsave->GetPad(2)->cd(3);
1144	gPad->Modified();
1145	gPad->cd();
1146	}
1147	*/
1148	// --------------------------------------------------------------------------
1149	//
1150	// This is a preliminary implementation of a alpha-fit procedure for
1151	// all possible source positions. It will be moved into its own
1152	// more powerfull class soon.
1153	//
1154	// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
1155	// [0]exp(-0.5((x-[1])/[2])^2) + [3] + [4]x + [5]x^2
1156	// or
1157	// Aexp(-0.5((x-mu)/sigma)^2) + a + bx + cx^2
1158	//
1159	// Parameter [1] is fixed to 0 while the alpha peak should be
1160	// symmetric around alpha=0.
1161	//
1162	// Parameter [4] is fixed to 0 because the first derivative at
1163	// alpha=0 should be 0, too.
1164	//
1165	// In a first step the background is fitted between bgmin and bgmax,
1166	// while the parameters [0]=0 and [2]=1 are fixed.
1167	//
1168	// In a second step the signal region (alpha<sigmax) is fittet using
1169	// the whole function with parameters [1], [3], [4] and [5] fixed.
1170	//
1171	// The number of excess and background events are calculated as
1172	// s = int(0, sigint, gaus(0)+pol2(3))
1173	// b = int(0, sigint, pol2(3))
1174	//
1175	// The Significance is calculated using the Significance() member
1176	// function.
1177	//
1178
1179	// I might need this
1180	/*
1181	TCanvas *c=new TCanvas;
1182
1183	gStyle->SetPalette(1, 0);
1184
1185	c->Divide(3,2, 0, 0);
1186	c->cd(1);
1187	gPad->SetBorderMode(0);
1188	hists->Draw("colz");
1189	hists->SetBit(kCanDelete);
1190	catalog->Draw("mirror same");
1191	c->cd(2);
1192	gPad->SetBorderMode(0);
1193	hist->Draw("colz");
1194	hist->SetBit(kCanDelete);
1195	catalog->Draw("mirror same");
1196	c->cd(3);
1197	gPad->SetBorderMode(0);
1198	histb->Draw("colz");
1199	histb->SetBit(kCanDelete);
1200	catalog->Draw("mirror same");
1201	c->cd(4);
1202	gPad->Divide(1,3, 0, 0);
1203	TVirtualPad *p=gPad;
1204	p->SetBorderMode(0);
1205	p->cd(1);
1206	gPad->SetBorderMode(0);
1207	h0b.DrawCopy();
1208	h0a.DrawCopy("same");
1209	p->cd(2);
1210	gPad->SetBorderMode(0);
1211	h3.DrawCopy();
1212	p->cd(3);
1213	gPad->SetBorderMode(0);
1214	h2.DrawCopy();
1215
1216	*/

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

Download in other formats: