source: trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc@ 8462

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