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

Last change on this file since 5971 was 5971, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 37.5 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//
113//////////////////////////////////////////////////////////////////////////////
114#include "MHFalseSource.h"
115
116#include <TF1.h>
117#include <TF2.h>
118#include <TH2.h>
119#include <TGraph.h>
120#include <TStyle.h>
121#include <TCanvas.h>
122#include <TRandom.h>
123#include <TPaveText.h>
124#include <TStopwatch.h>
125
126#include "MGeomCam.h"
127#include "MSrcPosCam.h"
128#include "MHillasSrc.h"
129#include "MTime.h"
130#include "MObservatory.h"
131#include "MPointingPos.h"
132#include "MAstroCatalog.h"
133#include "MAstroSky2Local.h"
134#include "MStatusDisplay.h"
135
136#include "MMath.h"
137#include "MAlphaFitter.h"
138
139#include "MBinning.h"
140#include "MParList.h"
141
142#include "MLog.h"
143#include "MLogManip.h"
144
145ClassImp(MHFalseSource);
146
147using namespace std;
148
149//class MHillasExt;
150
151// --------------------------------------------------------------------------
152//
153// Default Constructor
154//
155MHFalseSource::MHFalseSource(const char *name, const char *title)
156 : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), fAlphaCut(12.5),
157 fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinDW(-1), fMaxDW(-1),
158 fHistOff(0)
159{
160 //
161 // set the name and title of this object
162 //
163 fName = name ? name : "MHFalseSource";
164 fTitle = title ? title : "3D-plot of Alpha vs x, y";
165
166 fHist.SetDirectory(NULL);
167
168 fHist.SetName("Alpha");
169 fHist.SetTitle("3D-plot of Alpha vs x, y");
170 fHist.SetXTitle("x [\\circ]");
171 fHist.SetYTitle("y [\\circ]");
172 fHist.SetZTitle("\\alpha [\\circ]");
173}
174
175void MHFalseSource::MakeSymmetric(TH1 *h)
176{
177 h->SetMinimum();
178 h->SetMaximum();
179
180 const Float_t min = TMath::Abs(h->GetMinimum());
181 const Float_t max = TMath::Abs(h->GetMaximum());
182
183 const Float_t absmax = TMath::Max(min, max)*1.002;
184
185 h->SetMaximum( absmax);
186 h->SetMinimum(-absmax);
187}
188
189// --------------------------------------------------------------------------
190//
191// Set the alpha cut (|alpha|<fAlphaCut) which is use to estimate the
192// number of excess events
193//
194void MHFalseSource::SetAlphaCut(Float_t alpha)
195{
196 if (alpha<0)
197 *fLog << warn << "Alpha<0... taking absolute value." << endl;
198
199 fAlphaCut = TMath::Abs(alpha);
200
201 Modified();
202}
203
204// --------------------------------------------------------------------------
205//
206// Set mean alpha around which the off sample is determined
207// (fBgMean-fAlphaCut/2<|fAlpha|<fBgMean+fAlphaCut/2) which is use
208// to estimate the number of off events
209//
210void MHFalseSource::SetBgMean(Float_t alpha)
211{
212 if (alpha<0)
213 *fLog << warn << "Alpha<0... taking absolute value." << endl;
214
215 fBgMean = TMath::Abs(alpha);
216
217 Modified();
218}
219
220// --------------------------------------------------------------------------
221//
222// Set binnings (takes BinningFalseSource) and prepare filling of the
223// histogram.
224//
225// Also search for MTime, MObservatory and MPointingPos
226//
227//MHillasExt *ext=0;
228Bool_t MHFalseSource::SetupFill(const MParList *plist)
229{
230 const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
231 if (!geom)
232 {
233 *fLog << err << "MGeomCam not found... aborting." << endl;
234 return kFALSE;
235 }
236 /*
237 ext = (MHillasExt*)plist->FindObject("MHillasExt");
238 if (!ext)
239 {
240 *fLog << err << "MHillasExt not found... aborting." << endl;
241 return kFALSE;
242 }
243 */
244 fMm2Deg = geom->GetConvMm2Deg();
245
246 if (fName!=(TString)"MHFalseSourceOff" && fHistOff==NULL)
247 {
248 MHFalseSource *hoff = (MHFalseSource*)plist->FindObject("MHFalseSourceOff", "MHFalseSource");
249 if (!hoff)
250 *fLog << inf << "No MHFalseSourceOff [MHFalseSource] found... using current data only!" << endl;
251 else
252 {
253 *fLog << inf << "MHFalseSource [MHFalseSource] found... using on-off mode!" << endl;
254 SetOffData(*hoff);
255 }
256 }
257
258 if (fHistOff)
259 MH::SetBinning(&fHist, fHistOff);
260 else
261 {
262 MBinning binsa(18, 0, 90);
263 binsa.SetEdges(*plist, "BinningAlpha");
264
265 const MBinning *bins = (MBinning*)plist->FindObject("BinningFalseSource");
266 if (!bins || bins->IsDefault())
267 {
268 const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg;
269
270 MBinning b;
271 b.SetEdges(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 // FIXME: Because the pointing position could change we must check
295 // for the current pointing position and add a offset in the
296 // Fill function!
297 fRa = fPointPos ? fPointPos->GetRa() : 0;
298 fDec = fPointPos ? fPointPos->GetDec() : 90;
299
300 return kTRUE;
301}
302
303// --------------------------------------------------------------------------
304//
305// Fill the histogram. For details see the code or the class description
306//
307Bool_t MHFalseSource::Fill(const MParContainer *par, const Stat_t w)
308{
309 const MHillas *hil = dynamic_cast<const MHillas*>(par);
310 if (!hil)
311 {
312 *fLog << err << "MHFalseSource::Fill: No container specified!" << endl;
313 return kFALSE;
314 }
315
316 // Get max radius...
317 const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1));
318
319 // Get camera rotation angle
320 Double_t rho = 0;
321 if (fTime && fObservatory && fPointPos)
322 rho = fPointPos->RotationAngle(*fObservatory, *fTime);
323 //if (fPointPos)
324 // rho = fPointPos->RotationAngle(*fObservatory);
325
326 // Create necessary containers for calculation
327 MSrcPosCam src;
328 MHillasSrc hsrc;
329 hsrc.SetSrcPos(&src);
330
331 // Get number of bins and bin-centers
332 const Int_t nx = fHist.GetNbinsX();
333 const Int_t ny = fHist.GetNbinsY();
334 Axis_t cx[nx];
335 Axis_t cy[ny];
336 fHist.GetXaxis()->GetCenter(cx);
337 fHist.GetYaxis()->GetCenter(cy);
338
339 for (int ix=0; ix<nx; ix++)
340 {
341 for (int iy=0; iy<ny; iy++)
342 {
343 // check distance... to get a circle plot
344 if (TMath::Hypot(cx[ix], cy[iy])>maxr)
345 continue;
346
347 // rotate center of bin
348 // precalculation of sin/cos doesn't accelerate
349 TVector2 v(cx[ix], cy[iy]);
350 if (rho!=0)
351 v=v.Rotate(rho);
352
353 // convert degrees to millimeters
354 v *= 1./fMm2Deg;
355
356 if (fSrcPos)
357 v += fSrcPos->GetXY();
358
359 src.SetXY(v);
360
361 // Source dependant hillas parameters
362 if (hsrc.Calc(*hil/*, ext*/)>0)
363 {
364 *fLog << warn << "Calculation of MHillasSrc failed for x=" << cx[ix] << " y=" << cy[iy] << endl;
365 return kFALSE;
366 }
367
368 // FIXME: This should be replaced by an external MFilter
369 // and/or MTaskList
370 // Source dependant distance cut
371 if (fMinDist>0 && hsrc.GetDist()*fMm2Deg<fMinDist)
372 continue;
373 if (fMaxDist>0 && hsrc.GetDist()*fMm2Deg>fMaxDist)
374 continue;
375
376 if (fMaxDW>0 && hsrc.GetDist()>fMaxDW*hil->GetWidth())
377 continue;
378 if (fMinDW<0 && hsrc.GetDist()<fMinDW*hil->GetWidth())
379 continue;
380
381 // Fill histogram
382 const Double_t alpha = hsrc.GetAlpha();
383 fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w);
384 }
385 }
386
387 return kTRUE;
388}
389
390// --------------------------------------------------------------------------
391//
392// Create projection for off data, taking sum of bin contents of
393// range (fBgMean-fAlphaCut/2, fBgMean+fAlphaCut) Making sure to take
394// the same number of bins than for on-data
395//
396void MHFalseSource::ProjectOff(const TH3D &src, TH2D *h2, TH2D *all)
397{
398 TAxis &axe = *src.GetZaxis();
399
400 // Find range to cut (left edge and width)
401 const Int_t f = axe.FindBin(fBgMean-fAlphaCut/2);
402 const Int_t l = axe.FindBin(fAlphaCut)+f-1;
403
404 axe.SetRange(f, l);
405 const Float_t cut1 = axe.GetBinLowEdge(f);
406 const Float_t cut2 = axe.GetBinUpEdge(l);
407 h2->SetTitle(Form("Distribution of %.1f\\circ<|\\alpha|<%.1f\\circ in x,y", cut1, cut2));
408
409 // Get projection for range
410 TH2D *p = (TH2D*)src.Project3D("yx_off");
411
412 // Reset range
413 axe.SetRange(0,9999);
414
415 // Move contents from projection to h2
416 h2->Reset();
417 h2->Add(p, all->GetMaximum());
418 h2->Divide(all);
419
420 // Delete p
421 delete p;
422
423 // Set Minimum as minimum value Greater Than 0
424 h2->SetMinimum(GetMinimumGT(*h2));
425}
426
427// --------------------------------------------------------------------------
428//
429// Create projection for on data, taking sum of bin contents of
430// range (0, fAlphaCut)
431//
432void MHFalseSource::ProjectOn(const TH3D &src, TH2D *h3, TH2D *all)
433{
434 TAxis &axe = *src.GetZaxis();
435
436 // Find range to cut
437 axe.SetRangeUser(0, fAlphaCut);
438 const Float_t cut = axe.GetBinUpEdge(axe.GetLast());
439 h3->SetTitle(Form("Distribution of |\\alpha|<%.1f\\circ in x,y", cut));
440
441 // Get projection for range
442 TH2D *p = (TH2D*)src.Project3D("yx_on");
443
444 // Reset range
445 axe.SetRange(0,9999);
446
447 // Move contents from projection to h3
448 h3->Reset();
449 h3->Add(p, all->GetMaximum());
450 h3->Divide(all);
451
452 // Delete p
453 delete p;
454
455 // Set Minimum as minimum value Greater Than 0
456 h3->SetMinimum(GetMinimumGT(*h3));
457}
458
459// --------------------------------------------------------------------------
460//
461// Create projection for all data, taking sum of bin contents of
462// range (0, 90) - corresponding to the number of entries in this slice.
463//
464void MHFalseSource::ProjectAll(TH2D *h3)
465{
466 h3->SetTitle("Number of entries");
467
468 // Get projection for range
469 TH2D *p = (TH2D*)fHist.Project3D(Form("yx_%d", gRandom->Uniform(999999)));
470 p->SetDirectory(0);
471
472 // Move contents from projection to h3
473 h3->Reset();
474 h3->Add(p);
475 delete p;
476
477 // Set Minimum as minimum value Greater Than 0
478 h3->SetMinimum(GetMinimumGT(*h3));
479}
480
481void MHFalseSource::ProjectOnOff(TH2D *h2, TH2D *h0)
482{
483 ProjectOn(*fHistOff, h2, h0);
484
485 TH2D h;
486 MH::SetBinning(&h, h2);
487
488 // Divide by number of entries in off region (of off-data)
489 ProjectOff(*fHistOff, &h, h0);
490 h2->Divide(&h);
491
492 // Multiply by number of entries in off region (of on-data)
493 ProjectOff(fHist, &h, h0);
494 h2->Multiply(&h);
495
496 // Recalculate Minimum
497 h2->SetMinimum(GetMinimumGT(*h2));
498}
499
500// --------------------------------------------------------------------------
501//
502// Update the projections and paint them
503//
504void MHFalseSource::Paint(Option_t *opt)
505{
506 // Set pretty color palette
507 gStyle->SetPalette(1, 0);
508
509 TVirtualPad *padsave = gPad;
510
511 TH1D* h1;
512 TH2D* h0;
513 TH2D* h2;
514 TH2D* h3;
515 TH2D* h4;
516 TH2D* h5;
517
518 /*
519 fHistProjAll = Form("All_%p", this);
520 fHistProjOn = Form("On_%p", this);
521 fHistProjOff = Form("Off_%p", this);
522 fHistProjDiff = Form("Diff_%p", this);
523 fHistProjAll = Form("All_%p", this);
524 */
525
526 // Update projection of all-events
527 padsave->GetPad(2)->cd(3);
528 if ((h0 = (TH2D*)gPad->FindObject("Alpha_yx_all")))
529 ProjectAll(h0);
530
531 // Update projection of on-events
532 padsave->GetPad(1)->cd(1);
533 if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on")))
534 ProjectOn(fHist, h3, h0);
535
536 // Update projection of off-events
537 padsave->GetPad(1)->cd(3);
538 if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off")))
539 {
540 if (!fHistOff)
541 ProjectOff(fHist, h2, h0);
542 else
543 ProjectOnOff(h2, h0);
544 }
545
546 padsave->GetPad(2)->cd(2);
547 if ((h5 = (TH2D*)gPad->FindObject("Alpha_yx_diff")))
548 {
549 h5->Add(h2, h3, -1);
550 MakeSymmetric(h5);
551 }
552
553 // Update projection of significance
554 padsave->GetPad(1)->cd(2);
555 if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_yx_sig")))
556 {
557 const Int_t nx = h4->GetXaxis()->GetNbins();
558 const Int_t ny = h4->GetYaxis()->GetNbins();
559 //const Int_t nr = nx*nx + ny*ny;
560
561 Int_t maxx=nx/2;
562 Int_t maxy=ny/2;
563
564 Int_t max = h4->GetBin(nx, ny);
565
566 for (int ix=1; ix<=nx; ix++)
567 for (int iy=1; iy<=ny; iy++)
568 {
569 const Int_t n = h4->GetBin(ix, iy);
570
571 const Double_t s = h3->GetBinContent(n);
572 const Double_t b = h2->GetBinContent(n);
573
574 const Double_t sig = MMath::SignificanceLiMaSigned(s, b);
575
576 h4->SetBinContent(n, sig);
577
578 if (sig>h4->GetBinContent(max) && sig>0/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)<nr*nr/9*/)
579 {
580 max = n;
581 maxx=ix;
582 maxy=iy;
583 }
584 }
585
586 MakeSymmetric(h4);
587
588 // Update projection of 'the best alpha-plot'
589 padsave->GetPad(2)->cd(1);
590 if ((h1 = (TH1D*)gPad->FindObject("Alpha_z")) && max>0)
591 {
592 const Double_t x = h4->GetXaxis()->GetBinCenter(maxx);
593 const Double_t y = h4->GetYaxis()->GetBinCenter(maxy);
594 const Double_t s = h4->GetBinContent(max);
595
596 // This is because a 3D histogram x and y are vice versa
597 // Than for their projections
598 TH1 *h = fHist.ProjectionZ("Alpha_z", maxx, maxx, maxy, maxy);
599 h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (S_{max}=%.1f\\sigma)", x, y, s));
600
601 TH1D *h0=0;
602 if ((h0 = (TH1D*)gPad->FindObject("AlphaOff_z")))
603 {
604 fHistOff->ProjectionZ("AlphaOff_z", maxx, maxx, maxy, maxy);
605
606 const Int_t f = h0->GetXaxis()->FindFixBin(fBgMean-fAlphaCut/2);
607 const Int_t l = h0->GetXaxis()->FindFixBin(fAlphaCut)+f-1;
608 h0->Scale(h1->Integral(f, l)/h0->Integral(f, l));
609 //h0->Scale(h1->GetEntries()/h0->GetEntries());
610
611 }
612 }
613 }
614
615 gPad = padsave;
616}
617
618// --------------------------------------------------------------------------
619//
620// Get the MAstroCatalog corresponding to fRa, fDec. The limiting magnitude
621// is set to 9, while the fov is equal to the current fov of the false
622// source plot.
623//
624TObject *MHFalseSource::GetCatalog()
625{
626 const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1));
627
628 // Create catalog...
629 MAstroCatalog *stars = new MAstroCatalog;
630 stars->SetLimMag(9);
631 stars->SetGuiActive(kFALSE);
632 stars->SetRadiusFOV(maxr);
633 stars->SetRaDec(fRa*TMath::DegToRad()*15, fDec*TMath::DegToRad());
634 stars->ReadBSC("bsc5.dat");
635
636 stars->SetBit(kCanDelete);
637 return stars;
638}
639
640// --------------------------------------------------------------------------
641//
642// Draw the histogram
643//
644void MHFalseSource::Draw(Option_t *opt)
645{
646 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
647 pad->SetBorderMode(0);
648
649 AppendPad("");
650
651 pad->Divide(1, 2, 0, 0.03);
652
653// TObject *catalog = GetCatalog();
654
655 // Initialize upper part
656 pad->cd(1);
657 // Make sure that the catalog is deleted only once
658 // Normally this is not done by root, because it is not necessary...
659 // but in this case it is necessary, because the catalog is
660 // deleted by the first pad and the second one tries to do the same!
661// gROOT->GetListOfCleanups()->Add(gPad);
662 gPad->SetBorderMode(0);
663 gPad->Divide(3, 1);
664
665 // PAD #1
666 pad->GetPad(1)->cd(1);
667 gPad->SetBorderMode(0);
668 fHist.GetZaxis()->SetRangeUser(0,fAlphaCut);
669 TH1 *h3 = fHist.Project3D("yx_on");
670 fHist.GetZaxis()->SetRange(0,9999);
671 h3->SetDirectory(NULL);
672 h3->SetXTitle(fHist.GetXaxis()->GetTitle());
673 h3->SetYTitle(fHist.GetYaxis()->GetTitle());
674 h3->Draw("colz");
675 h3->SetBit(kCanDelete);
676// catalog->Draw("mirror same *");
677
678 // PAD #2
679 pad->GetPad(1)->cd(2);
680 gPad->SetBorderMode(0);
681 fHist.GetZaxis()->SetRange(0,0);
682 TH1 *h4 = fHist.Project3D("yx_sig"); // Do this to get the correct binning....
683 fHist.GetZaxis()->SetRange(0,9999);
684 h4->SetTitle("Significance");
685 h4->SetDirectory(NULL);
686 h4->SetXTitle(fHist.GetXaxis()->GetTitle());
687 h4->SetYTitle(fHist.GetYaxis()->GetTitle());
688 h4->Draw("colz");
689 h4->SetBit(kCanDelete);
690// catalog->Draw("mirror same *");
691
692 // PAD #3
693 pad->GetPad(1)->cd(3);
694 gPad->SetBorderMode(0);
695 TH1 *h2 = 0;
696 if (fHistOff)
697 {
698 fHistOff->GetZaxis()->SetRangeUser(0,fAlphaCut);
699 h2 = fHistOff->Project3D("yx_off");
700 fHistOff->GetZaxis()->SetRange(0,9999);
701 }
702 else
703 {
704 fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
705 h2 = fHist.Project3D("yx_off");
706 fHist.GetZaxis()->SetRange(0,9999);
707 }
708 h2->SetDirectory(NULL);
709 h2->SetXTitle(fHist.GetXaxis()->GetTitle());
710 h2->SetYTitle(fHist.GetYaxis()->GetTitle());
711 h2->Draw("colz");
712 h2->SetBit(kCanDelete);
713// catalog->Draw("mirror same *");
714
715 // Initialize lower part
716 pad->cd(2);
717 // Make sure that the catalog is deleted only once
718// gROOT->GetListOfCleanups()->Add(gPad);
719 gPad->SetBorderMode(0);
720 gPad->Divide(3, 1);
721
722 // PAD #4
723 pad->GetPad(2)->cd(1);
724 gPad->SetBorderMode(0);
725 TH1 *h1 = fHist.ProjectionZ("Alpha_z");
726 h1->SetDirectory(NULL);
727 h1->SetTitle("Distribution of \\alpha");
728 h1->SetXTitle(fHist.GetZaxis()->GetTitle());
729 h1->SetYTitle("Counts");
730 h1->Draw();
731 h1->SetBit(kCanDelete);
732 if (fHistOff)
733 {
734 h1->SetLineColor(kGreen);
735
736 h1 = fHistOff->ProjectionZ("AlphaOff_z");
737 h1->SetDirectory(NULL);
738 h1->SetTitle("Distribution of \\alpha");
739 h1->SetXTitle(fHistOff->GetZaxis()->GetTitle());
740 h1->SetYTitle("Counts");
741 h1->Draw("same");
742 h1->SetBit(kCanDelete);
743 h1->SetLineColor(kRed);
744 }
745
746 // PAD #5
747 pad->GetPad(2)->cd(2);
748 gPad->SetBorderMode(0);
749 TH1 *h5 = (TH1*)h3->Clone("Alpha_yx_diff");
750 h5->Add(h2, -1);
751 h5->SetTitle("Difference of on- and off-distribution");
752 h5->SetDirectory(NULL);
753 h5->Draw("colz");
754 h5->SetBit(kCanDelete);
755// catalog->Draw("mirror same *");
756
757 // PAD #6
758 pad->GetPad(2)->cd(3);
759 gPad->SetBorderMode(0);
760 TH1 *h0 = fHist.Project3D("yx_all");
761 h0->SetDirectory(NULL);
762 h0->SetXTitle(fHist.GetXaxis()->GetTitle());
763 h0->SetYTitle(fHist.GetYaxis()->GetTitle());
764 h0->Draw("colz");
765 h0->SetBit(kCanDelete);
766// catalog->Draw("mirror same *");
767}
768
769// --------------------------------------------------------------------------
770//
771// Everything which is in the main pad belongs to this class!
772//
773Int_t MHFalseSource::DistancetoPrimitive(Int_t px, Int_t py)
774{
775 return 0;
776}
777
778// --------------------------------------------------------------------------
779//
780// Set all sub-pads to Modified()
781//
782void MHFalseSource::Modified()
783{
784 if (!gPad)
785 return;
786
787 TVirtualPad *padsave = gPad;
788 padsave->Modified();
789 padsave->GetPad(1)->cd(1);
790 gPad->Modified();
791 padsave->GetPad(1)->cd(3);
792 gPad->Modified();
793 padsave->GetPad(2)->cd(1);
794 gPad->Modified();
795 padsave->GetPad(2)->cd(2);
796 gPad->Modified();
797 padsave->GetPad(2)->cd(3);
798 gPad->Modified();
799 gPad->cd();
800}
801
802// --------------------------------------------------------------------------
803//
804// The following resources are available:
805//
806// MHFalseSource.DistMin: 0.5
807// MHFalseSource.DistMax: 1.4
808// MHFalseSource.DWMin: 0.1
809// MHFalseSource.DWMax: 0.3
810//
811Int_t MHFalseSource::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
812{
813 Bool_t rc = kFALSE;
814 if (IsEnvDefined(env, prefix, "DistMin", print))
815 {
816 rc = kTRUE;
817 SetMinDist(GetEnvValue(env, prefix, "DistMin", fMinDist));
818 }
819 if (IsEnvDefined(env, prefix, "DistMax", print))
820 {
821 rc = kTRUE;
822 SetMaxDist(GetEnvValue(env, prefix, "DistMax", fMaxDist));
823 }
824
825 if (IsEnvDefined(env, prefix, "DWMin", print))
826 {
827 rc = kTRUE;
828 SetMinDW(GetEnvValue(env, prefix, "DWMin", fMinDist));
829 }
830 if (IsEnvDefined(env, prefix, "DWMax", print))
831 {
832 rc = kTRUE;
833 SetMaxDW(GetEnvValue(env, prefix, "DWMax", fMaxDist));
834 }
835
836 return rc;
837}
838
839Double_t FcnGauss2d(Double_t *x, Double_t *par)
840{
841 TVector2 v = TVector2(x[0], x[1]).Rotate(par[5]*TMath::DegToRad());
842
843 const Double_t g0 = TMath::Gaus(v.X(), par[1], par[2]);
844 const Double_t g1 = TMath::Gaus(v.Y(), par[3], par[4]);
845
846 //Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE);
847 return par[0]*g0*g1;
848}
849
850// --------------------------------------------------------------------------
851//
852// This is a preliminary implementation of a alpha-fit procedure for
853// all possible source positions. It will be moved into its own
854// more powerfull class soon.
855//
856// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
857// [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
858// or
859// A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
860//
861// Parameter [1] is fixed to 0 while the alpha peak should be
862// symmetric around alpha=0.
863//
864// Parameter [4] is fixed to 0 because the first derivative at
865// alpha=0 should be 0, too.
866//
867// In a first step the background is fitted between bgmin and bgmax,
868// while the parameters [0]=0 and [2]=1 are fixed.
869//
870// In a second step the signal region (alpha<sigmax) is fittet using
871// the whole function with parameters [1], [3], [4] and [5] fixed.
872//
873// The number of excess and background events are calculated as
874// s = int(0, sigint, gaus(0)+pol2(3))
875// b = int(0, sigint, pol2(3))
876//
877// The Significance is calculated using the Significance() member
878// function.
879//
880void MHFalseSource::FitSignificance(Float_t sigint, Float_t sigmax, Float_t bgmin, Float_t bgmax, Byte_t polynom)
881{
882// TObject *catalog = GetCatalog();
883
884 TH1D h0a("A", "", 50, 0, 4000);
885 TH1D h4a("chisq1", "", 50, 0, 35);
886 //TH1D h5a("prob1", "", 50, 0, 1.1);
887 TH1D h6("signifcance", "", 50, -20, 20);
888
889 TH1D h1("mu", "Parameter \\mu", 50, -1, 1);
890 TH1D h2("sigma", "Parameter \\sigma", 50, 0, 90);
891 TH1D h3("b", "Parameter b", 50, -0.1, 0.1);
892
893 TH1D h0b("a", "Parameter a (red), A (blue)", 50, 0, 4000);
894 TH1D h4b("\\chi^{2}", "\\chi^{2} (red, green) / significance (black)", 50, 0, 35);
895 //TH1D h5b("prob", "Fit probability: Bg(red), F(blue)", 50, 0, 1.1);
896
897 h0a.SetLineColor(kBlue);
898 h4a.SetLineColor(kBlue);
899 //h5a.SetLineColor(kBlue);
900 h0b.SetLineColor(kRed);
901 h4b.SetLineColor(kRed);
902 //h5b.SetLineColor(kRed);
903
904 TH1 *hist = fHist.Project3D("yx_fit");
905 hist->SetDirectory(0);
906 TH1 *hists = fHist.Project3D("yx_fit");
907 hists->SetDirectory(0);
908 TH1 *histb = fHist.Project3D("yx_fit");
909 histb->SetDirectory(0);
910 hist->Reset();
911 hists->Reset();
912 histb->Reset();
913 hist->SetNameTitle("Significance",
914 Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
915 sigmax, bgmin, bgmax));
916 hists->SetName("Excess");
917 histb->SetName("Background");
918 hist->SetXTitle(fHist.GetXaxis()->GetTitle());
919 hists->SetXTitle(fHist.GetXaxis()->GetTitle());
920 histb->SetXTitle(fHist.GetXaxis()->GetTitle());
921 hist->SetYTitle(fHist.GetYaxis()->GetTitle());
922 hists->SetYTitle(fHist.GetYaxis()->GetTitle());
923 histb->SetYTitle(fHist.GetYaxis()->GetTitle());
924
925 const Double_t w = fHist.GetZaxis()->GetBinWidth(1);
926
927 TArrayD maxpar;
928
929 /* func.SetParName(0, "A");
930 * func.SetParName(1, "mu");
931 * func.SetParName(2, "sigma");
932 */
933
934 const Int_t nx = hist->GetXaxis()->GetNbins();
935 const Int_t ny = hist->GetYaxis()->GetNbins();
936 //const Int_t nr = nx*nx+ny*ny;
937
938 Double_t maxalpha0=0;
939 Double_t maxs=3;
940
941 Int_t maxx=0;
942 Int_t maxy=0;
943
944 TStopwatch clk;
945 clk.Start();
946
947 *fLog << inf;
948 *fLog << "Signal fit: alpha < " << sigmax << endl;
949 *fLog << "Integration: alpha < " << sigint << endl;
950 *fLog << "Background fit: " << bgmin << " < alpha < " << bgmax << endl;
951 *fLog << "Polynom order: " << (int)polynom << endl;
952 *fLog << "Fitting False Source Plot..." << flush;
953
954 TH1 *h0 = fHist.Project3D("yx_entries");
955 Float_t entries = h0->GetMaximum();
956 delete h0;
957
958 MAlphaFitter fit;
959 fit.SetSignalIntegralMax(sigint);
960 fit.SetSignalFitMax(sigmax);
961 fit.SetBackgroundFitMin(bgmin);
962 fit.SetBackgroundFitMax(bgmax);
963 fit.SetPolynomOrder(polynom);
964
965 TH1D *h=0, *hoff=0;
966 Double_t scale = 1;
967 for (int ix=1; ix<=nx; ix++)
968 for (int iy=1; iy<=ny; iy++)
969 {
970 // This is because a 3D histogram x and y are vice versa
971 // Than for their projections
972 h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy);
973
974 if (h->GetEntries()==0)
975 continue;
976
977 // This is a quick hack to correct for the lower acceptance at
978 // the edge of the camera
979 h->Scale(entries/h->GetEntries());
980
981 if (fHistOff)
982 {
983 hoff = fHistOff->ProjectionZ("AlphaFitOff", ix, ix, iy, iy);
984 // This is a quick hack to correct for the lower acceptance at
985 // the edge of the camera
986 hoff->Scale(entries/h->GetEntries());
987 scale = fit.Scale(*hoff, *h);
988 }
989
990 if (!fit.Fit(*h, hoff, scale))
991 continue;
992
993 const Double_t alpha0 = h->GetBinContent(1);
994 if (alpha0>maxalpha0)
995 maxalpha0=alpha0;
996
997 // Fill results into some histograms
998 h0a.Fill(fit.GetGausA()); // gaus-A
999 h0b.Fill(fit.GetCoefficient(3)); // 3
1000 h1.Fill(fit.GetGausMu()); // mu
1001 h2.Fill(fit.GetGausSigma()); // sigma-gaus
1002 if (polynom>1 && !fHistOff)
1003 h3.Fill(fit.GetCoefficient(5));
1004 h4b.Fill(fit.GetChiSqSignal());
1005
1006 const Double_t sig = fit.GetSignificance();
1007 const Double_t b = fit.GetEventsBackground();
1008 const Double_t s = fit.GetEventsSignal();
1009
1010 const Int_t n = hist->GetBin(ix, iy);
1011 hists->SetBinContent(n, s-b);
1012 histb->SetBinContent(n, b);
1013
1014 hist->SetBinContent(n, sig);
1015 if (sig!=0)
1016 h6.Fill(sig);
1017
1018 if (sig>maxs)
1019 {
1020 maxs = sig;
1021 maxx = ix;
1022 maxy = iy;
1023 maxpar = fit.GetCoefficients();
1024 }
1025 }
1026
1027 *fLog << "Done." << endl;
1028
1029 h0a.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
1030 h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
1031
1032 hists->SetTitle(Form("Excess events for \\alpha<%.0f\\circ (N_{max}=%d)", sigint, (int)hists->GetMaximum()));
1033 histb->SetTitle(Form("Background events for \\alpha<%.0f\\circ", sigint));
1034
1035 //hists->SetMinimum(GetMinimumGT(*hists));
1036 histb->SetMinimum(GetMinimumGT(*histb));
1037
1038 MakeSymmetric(hists);
1039 MakeSymmetric(hist);
1040
1041 clk.Stop();
1042 clk.Print("m");
1043
1044 TCanvas *c=new TCanvas;
1045
1046 gStyle->SetPalette(1, 0);
1047
1048 c->Divide(3,2, 0, 0);
1049 c->cd(1);
1050 gPad->SetBorderMode(0);
1051 hists->Draw("colz");
1052 hists->SetBit(kCanDelete);
1053// catalog->Draw("mirror same *");
1054 c->cd(2);
1055 gPad->SetBorderMode(0);
1056 hist->Draw("colz");
1057 hist->SetBit(kCanDelete);
1058
1059
1060 const Double_t maxr = 0.9*TMath::Abs(fHist.GetBinCenter(1));
1061 TF2 f2d("Gaus-2D", FcnGauss2d, -maxr, maxr, -maxr, maxr, 6);
1062 f2d.SetLineWidth(1);
1063 f2d.SetParName(0, "Max sigma");
1064 f2d.SetParName(1, "Mean_1 deg");
1065 f2d.SetParName(2, "Sigma_1 deg");
1066 f2d.SetParName(3, "Mean_2 deg");
1067 f2d.SetParName(4, "Sigma_2 deg");
1068 f2d.SetParName(5, "Phi deg");
1069 f2d.SetParLimits(1, -maxr/2, maxr/2); // mu_1
1070 f2d.SetParLimits(3, -maxr/2, maxr/2); // mu_2
1071 f2d.SetParLimits(2, 0, maxr); // sigma_1
1072 f2d.SetParLimits(4, 0, maxr); // sigma_2
1073 f2d.SetParLimits(5, 0, 45); // phi
1074 f2d.SetParameter(0, maxs); // A
1075 f2d.SetParameter(1, hist->GetXaxis()->GetBinCenter(maxx)); // mu_1
1076 f2d.SetParameter(2, 0.1); // sigma_1
1077 f2d.SetParameter(3, hist->GetYaxis()->GetBinCenter(maxy)); // mu_2
1078 f2d.SetParameter(4, 0.1); // sigma_2
1079 f2d.FixParameter(5, 0); // phi
1080 hist->Fit(&f2d, "NI0R");
1081 f2d.DrawCopy("cont2same");
1082
1083
1084// catalog->Draw("mirror same *");
1085 c->cd(3);
1086 gPad->SetBorderMode(0);
1087 histb->Draw("colz");
1088 histb->SetBit(kCanDelete);
1089// catalog->Draw("mirror same *");
1090 c->cd(4);
1091 gPad->Divide(1,3, 0, 0);
1092 TVirtualPad *p=gPad;
1093 p->SetBorderMode(0);
1094 p->cd(1);
1095 gPad->SetBorderMode(0);
1096 h0b.DrawCopy();
1097 h0a.DrawCopy("same");
1098 p->cd(2);
1099 gPad->SetBorderMode(0);
1100 h3.DrawCopy();
1101 p->cd(3);
1102 gPad->SetBorderMode(0);
1103 h2.DrawCopy();
1104 c->cd(6);
1105 gPad->Divide(1,2, 0, 0);
1106 TVirtualPad *q=gPad;
1107 q->SetBorderMode(0);
1108 q->cd(1);
1109 gPad->SetBorderMode(0);
1110 gPad->SetBorderMode(0);
1111 h4b.DrawCopy();
1112 h4a.DrawCopy("same");
1113 h6.DrawCopy("same");
1114 q->cd(2);
1115 gPad->SetBorderMode(0);
1116 //h5b.DrawCopy();
1117 //h5a.DrawCopy("same");
1118 c->cd(5);
1119 gPad->SetBorderMode(0);
1120 if (maxx>0 && maxy>0)
1121 {
1122 const char *title = Form(" \\alpha for x=%.2f y=%.2f (S_{max}=%.1f\\sigma) ",
1123 hist->GetXaxis()->GetBinCenter(maxx),
1124 hist->GetYaxis()->GetBinCenter(maxy), maxs);
1125
1126 h = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
1127 h->Scale(entries/h->GetEntries());
1128
1129 h->SetDirectory(NULL);
1130 h->SetNameTitle("Result \\alpha", title);
1131 h->SetBit(kCanDelete);
1132 h->SetXTitle("\\alpha [\\circ]");
1133 h->SetYTitle("Counts");
1134 h->Draw();
1135
1136 if (fHistOff)
1137 {
1138 h->SetLineColor(kGreen);
1139
1140 TH1D *hof=fHistOff->ProjectionZ("AlphaFitOff", maxx, maxx, maxy, maxy);
1141 hof->Scale(entries/hof->GetEntries());
1142
1143 fit.Scale(*(TH1D*)hof, *(TH1D*)h);
1144
1145 hof->SetLineColor(kRed);
1146 hof->SetDirectory(NULL);
1147 hof->SetNameTitle("Result \\alpha", title);
1148 hof->SetBit(kCanDelete);
1149 hof->SetXTitle("\\alpha [\\circ]");
1150 hof->SetYTitle("Counts");
1151 hof->Draw("same");
1152
1153 TH1D *diff = (TH1D*)h->Clone("AlphaFitOnOff");
1154 diff->Add(hof, -1);
1155 diff->SetLineColor(kBlue);
1156 diff->SetNameTitle("Result On-Off \\alpha", title);
1157 diff->SetBit(kCanDelete);
1158 diff->SetXTitle("\\alpha [\\circ]");
1159 diff->SetYTitle("Counts");
1160 diff->Draw("same");
1161
1162 h->SetMinimum(diff->GetMinimum()<0 ? diff->GetMinimum()*1.2 : 0);
1163
1164 TLine l;
1165 l.DrawLine(0, 0, 90, 0);
1166 }
1167
1168 TF1 f1("f1", Form("gaus(0) + pol%d(3)", fHistOff ? 0 : polynom), 0, 90);
1169 TF1 f2("f2", Form("gaus(0) + pol%d(3)", fHistOff ? 0 : polynom), 0, 90);
1170 f1.SetParameters(maxpar.GetArray());
1171 f2.SetParameters(maxpar.GetArray());
1172 f2.FixParameter(0, 0);
1173 f2.FixParameter(1, 0);
1174 f2.FixParameter(2, 1);
1175 f1.SetLineColor(kGreen);
1176 f2.SetLineColor(kRed);
1177
1178 f2.DrawCopy("same");
1179 f1.DrawCopy("same");
1180
1181 TPaveText *leg = new TPaveText(0.35, 0.10, 0.90, 0.35, "brNDC");
1182 leg->SetBorderSize(1);
1183 leg->SetTextSize(0.04);
1184 leg->AddText(0.5, 0.82, Form("A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + pol%d", polynom))->SetTextAlign(22);
1185 //leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22);
1186 leg->AddLine(0, 0.65, 0, 0.65);
1187 leg->AddText(0.06, 0.54, Form("A=%.2f", maxpar[0]))->SetTextAlign(11);
1188 leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", maxpar[2]))->SetTextAlign(11);
1189 leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fix)", maxpar[1]))->SetTextAlign(11);
1190 leg->AddText(0.60, 0.54, Form("a=%.2f", maxpar[3]))->SetTextAlign(11);
1191 leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11);
1192 if (polynom>1)
1193 leg->AddText(0.60, 0.14, Form("c=%.2f", !fHistOff?maxpar[5]:0))->SetTextAlign(11);
1194 leg->SetBit(kCanDelete);
1195 leg->Draw();
1196
1197 q->cd(2);
1198
1199 TGraph *g = new TGraph;
1200 g->SetPoint(0, 0, 0);
1201
1202 Int_t max=0;
1203 Float_t maxsig=0;
1204 for (int i=1; i<89; i++)
1205 {
1206 const Double_t s = f1.Integral(0, (float)i)/w;
1207 const Double_t b = f2.Integral(0, (float)i)/w;
1208
1209 const Double_t sig = MMath::SignificanceLiMaSigned(s, b);
1210
1211 g->SetPoint(g->GetN(), i, sig);
1212
1213 if (sig>maxsig)
1214 {
1215 max = i;
1216 maxsig = sig;
1217 }
1218 }
1219
1220 g->SetNameTitle("SigVs\\alpha", "Significance vs \\alpha");
1221 g->GetHistogram()->SetXTitle("\\alpha_{0} [\\circ]");
1222 g->GetHistogram()->SetYTitle("Significance");
1223 g->SetBit(kCanDelete);
1224 g->Draw("AP");
1225
1226 leg = new TPaveText(0.35, 0.10, 0.90, 0.25, "brNDC");
1227 leg->SetBorderSize(1);
1228 leg->SetTextSize(0.1);
1229 leg->AddText(Form("S_{max}=%.1f\\sigma at \\alpha_{max}=%d\\circ", maxsig, max));
1230 leg->SetBit(kCanDelete);
1231 leg->Draw();
1232 }
1233}
Note: See TracBrowser for help on using the repository browser.