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

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