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

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