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

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