source: trunk/Mars/mhflux/MHFalseSource.cc@ 11014

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