source: trunk/MagicSoft/Mars/mhflux/MHPhi.cc@ 7594

Last change on this file since 7594 was 7594, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 14.8 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, 07/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MHPhi
28//
29// Plot delta phi the angle between the reconstructed shower origin and
30// the source position.
31//
32// More detail can be found at:
33// http://www.astro.uni-wuerzburg.de/results/ringmethod/
34//
35// Class Version 2:
36// + TH1D fHPhi; // Phi plot of the signal w.r.t. source
37// + TH1D fHPhiOff; // Phi plot of the signal w.r.t. source+180deg
38// + TH1D fHTemplate; // Template how the background should look in the ideal case
39//
40// + TH1D fHInhom; // Phi plot off the signal w.r.t. source (out of the ring with the signal)
41// + TH1D fHInhomOff; // Phi plot off the signal w.r.t. source+180deg (out of the ring with the signal)
42//
43// + Int_t fNumBinsSignal; // Number of bins for signal region
44// + Float_t fThetaCut; // Theta cut
45// + Float_t fDistSrc; // Nominal distance of source from center in wobble
46//
47// + Bool_t fUseAntiPhiCut; // Whether to use a anti-phi cut or not
48//
49////////////////////////////////////////////////////////////////////////////
50#include "MHPhi.h"
51
52#include <TCanvas.h>
53#include <TVector2.h>
54#include <TLatex.h>
55#include <TLine.h>
56#include <TMarker.h>
57
58#include "MLog.h"
59#include "MLogManip.h"
60
61#include "MParList.h"
62
63#include "MHillas.h"
64#include "MSrcPosCam.h"
65#include "MParameters.h"
66#include "MGeomCam.h"
67#include "MBinning.h"
68#include "MMath.h"
69
70ClassImp(MHPhi);
71
72using namespace std;
73
74// --------------------------------------------------------------------------
75//
76// Setup histograms
77//
78MHPhi::MHPhi(const char *name, const char *title)
79 : fHillas(0), fSrcPos(0), fDisp(0),
80 fNumBinsSignal(2), fThetaCut(0.23), fDistSrc(0.4),
81 fUseAntiPhiCut(kTRUE)
82{
83 fName = name ? name : "MHPhi";
84 fTitle = title ? title : "Graphs for rate data";
85
86 // Init Graphs
87 fHPhi.SetNameTitle("Phi", "\\Delta\\Phi-Distribution");
88 fHPhi.SetXTitle("\\Delta\\Phi' [\\circ]");
89 fHPhi.SetYTitle("Counts");
90
91 fHPhi.SetMinimum(0);
92 fHPhi.SetDirectory(0);
93 fHPhi.Sumw2();
94 fHPhi.SetBit(TH1::kNoStats);
95 fHPhi.SetMarkerStyle(kFullDotMedium);
96 fHPhi.SetLineColor(kBlue);
97 fHPhi.SetMarkerColor(kBlue);
98 fHPhi.GetYaxis()->SetTitleOffset(1.3);
99
100 fHPhiOff.SetMinimum(0);
101 fHPhiOff.SetDirectory(0);
102 fHPhiOff.Sumw2();
103 fHPhiOff.SetBit(TH1::kNoStats);
104 fHPhiOff.SetLineColor(kRed);
105 fHPhiOff.SetMarkerColor(kRed);
106
107 fHTemplate.SetMinimum(0);
108 fHTemplate.SetDirectory(0);
109 fHTemplate.SetBit(TH1::kNoStats);
110 fHTemplate.SetLineColor(kGreen);
111
112 fHInhom.SetNameTitle("Inhomogeneity", "\\Delta\\Phi-Distribution");
113 fHInhom.SetXTitle("\\Delta\\Phi' [\\circ]");
114 fHInhom.SetYTitle("Counts");
115 fHInhom.Sumw2();
116 fHInhom.SetMinimum(0);
117 fHInhom.SetDirectory(0);
118 fHInhom.SetBit(TH1::kNoStats);
119 fHInhom.GetYaxis()->SetTitleOffset(1.3);
120
121 fHInhomOff.Sumw2();
122 fHInhomOff.SetMinimum(0);
123 fHInhomOff.SetDirectory(0);
124 fHInhomOff.SetBit(TH1::kNoStats);
125 fHInhomOff.SetLineColor(kRed);
126 fHInhomOff.SetMarkerColor(kRed);
127}
128
129// --------------------------------------------------------------------------
130//
131// Setup the Binning for the histograms automatically if the correct
132// instances of MBinning
133//
134Bool_t MHPhi::SetupFill(const MParList *plist)
135{
136 fHillas = (MHillas*)plist->FindObject("MHillas");
137 if (!fHillas)
138 {
139 *fLog << err << "MHillas not found... abort." << endl;
140 return kFALSE;
141 }
142 fSrcPos = (MSrcPosCam*)plist->FindObject("MSrcPosCam");
143 if (!fSrcPos)
144 {
145 *fLog << err << "MSrcPosCam not found... abort." << endl;
146 return kFALSE;
147 }
148 fDisp = (MParameterD*)plist->FindObject("Disp", "MParameterD");
149 if (!fDisp)
150 {
151 *fLog << err << "Disp [MParameterD] not found... abort." << endl;
152 return kFALSE;
153 }
154
155 MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
156 if (!geom)
157 {
158 *fLog << err << "MGeomCam not found... abort." << endl;
159 return kFALSE;
160 }
161
162 fConvMm2Deg = geom->GetConvMm2Deg();
163
164 MParameterD *cut = (MParameterD*)plist->FindObject("ThetaSquaredCut", "MParameterD");
165 if (!cut)
166 *fLog << warn << "ThetaSquareCut [MParameterD] not found... using default theta<" << fThetaCut << "." << endl;
167 else
168 fThetaCut = TMath::Sqrt(cut->GetVal());
169
170 const Double_t w = TMath::ATan(fThetaCut/fDistSrc);
171 const Float_t sz = TMath::RadToDeg()*w/fNumBinsSignal;
172 const Int_t n = TMath::Nint(TMath::Ceil(180/sz));
173
174 MBinning(n+3, 0, (n+3)*sz).Apply(fHPhi);
175 MBinning(n+3, 0, (n+3)*sz).Apply(fHPhiOff);
176 MBinning(n+3, 0, (n+3)*sz).Apply(fHTemplate);
177 MBinning(n+3, 0, (n+3)*sz).Apply(fHInhom);
178 MBinning(n+3, 0, (n+3)*sz).Apply(fHInhomOff);
179
180 return kTRUE;
181}
182
183// --------------------------------------------------------------------------
184//
185// Fill the histograms with data from a MMuonCalibPar and
186// MMuonSearchPar container.
187//
188Bool_t MHPhi::Fill(const MParContainer *par, const Stat_t weight)
189{
190 // Here we calculate an upper phi cut to take a
191 // possible anti-theta cut into account
192 const Double_t ulim = fUseAntiPhiCut ? 180-fHPhi.GetBinLowEdge(fNumBinsSignal+1)*1.1 : 180;
193
194 // Calculate the shower origin and the source position in units of deg
195 const TVector2 pos = fHillas->GetMean()*fConvMm2Deg + fHillas->GetNormAxis()*fDisp->GetVal();
196 const TVector2 src = fSrcPos->GetXY()*fConvMm2Deg;
197
198 // Calculate radial distance between shower origin and source
199 const Double_t d = pos.Mod() - src.Mod();
200
201 // define an upper and lower cut for the radial distance between both
202 const Double_t dR = fThetaCut;
203 const Double_t dr = fThetaCut*0.913;
204
205 // calculate the phi-angle of the shower origin w.r.t. the source position
206 const Double_t delta = src.DeltaPhi(pos)*TMath::RadToDeg();
207
208 // Define the radii of the upper and lower ring border
209 const Double_t R = src.Mod()+dR;
210 const Double_t r = src.Mod()-dr;
211
212 // Calculate a scale to scale all source positions to the
213 // nominal distance to center
214 const Double_t scale = src.Mod()/fDistSrc;
215
216 // Fill a phi-histograms with all events outside the ring
217 // Take the upper phi cut into account
218 if ((d<-dr || d>dR)/*TMath::Abs(d)>fThetaCut*1.2*/ && TMath::Abs(delta)<ulim)
219 {
220 fHInhom.Fill(TMath::Abs(delta)*scale, weight);
221 fHInhomOff.Fill((ulim-TMath::Abs(delta))*scale, weight);
222 }
223
224 // Now forget about all events which are not inside the ring
225 if (d<-dr || d>dR)
226 return kTRUE;
227
228 // Fill the histograms for on and off with the scaled phi
229 // only if we are below the upper phi cut
230 if (TMath::Abs(delta)<ulim)
231 {
232 fHPhi.Fill( TMath::Abs(delta)*scale, weight);
233 fHPhiOff.Fill((ulim-TMath::Abs(delta))*scale, weight);
234 }
235
236 // Calculate the maximum scaled phi taking the upper phi cut into account
237 const Double_t max = scale*ulim;
238
239 // Fill a template, this is how the phi-plot would look like
240 // without a signal and for an ideal camera.
241 const Int_t n = fHTemplate.GetNbinsX();
242 TArrayD arr(n);
243 for (int i=1; i<=n; i++)
244 {
245 const Double_t hi = fHTemplate.GetBinLowEdge(i+1);
246 const Double_t lo = fHTemplate.GetBinLowEdge(i);
247
248 // Decide whether the bin is fully contained in the upper phi-cut or
249 // the maximum possible phi is inside the bin
250 if (hi<max)
251 arr[i-1] = 1;
252 else
253 if (lo<max) // if its inside calculate the ratio
254 arr[i-1] = (max-lo)/fHTemplate.GetBinWidth(i+1);
255 else
256 break;
257 }
258
259 // The template is scaled with the current ring width. The upper phi-
260 // cut must not be taken into account because it is just a constant
261 // for all events.
262 const Double_t sum = arr.GetSum();
263 for (int i=1; i<=n; i++)
264 fHTemplate.AddBinContent(i, (R*R-r*r)*arr[i-1]/sum);
265
266 return kTRUE;
267}
268
269// --------------------------------------------------------------------------
270//
271// This displays the TGraph like you expect it to be (eg. time on the axis)
272// It should also make sure, that the displayed time always is UTC and
273// not local time...
274//
275void MHPhi::Draw(Option_t *opt)
276{
277 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
278 pad->SetBorderMode(0);
279
280 AppendPad("update");
281
282 pad->Divide(2,2);
283
284 // --------------------------
285
286 pad->cd(1);
287 gPad->SetBorderMode(0);
288 gPad->SetFrameBorderMode(0);
289
290 fHPhi.Draw();
291 fHPhiOff.Draw("same");
292
293 TH1D *h1 = new TH1D(fHTemplate);
294 h1->SetName("Template");
295 h1->SetBit(kCanDelete);
296 h1->SetDirectory(0);
297 h1->Draw("same");
298
299 // --------------------------
300
301 pad->cd(2);
302 gPad->SetBorderMode(0);
303 gPad->SetFrameBorderMode(0);
304
305 fHInhom.Draw();
306 fHInhomOff.Draw("same");
307
308 TH1D *h2 = new TH1D(fHTemplate);
309 h2->SetName("Template");
310 h2->SetBit(kCanDelete);
311 h2->SetDirectory(0);
312 h2->Draw("same");
313
314 // --------------------------
315
316 pad->cd(3);
317 gPad->SetBorderMode(0);
318 gPad->SetFrameBorderMode(0);
319
320 fHPhi.Draw();
321 TH1D *h4 = new TH1D(fHInhom);
322 h4->SetName("Inhomogeneity");
323 h4->SetBit(kCanDelete);
324 h4->SetDirectory(0);
325 h4->Draw("same");
326
327 h1->Draw("same");
328
329 // --------------------------
330
331 pad->cd(4);
332 gPad->SetBorderMode(0);
333 gPad->SetFrameBorderMode(0);
334
335 TH1D *h3 = new TH1D(fHPhi);
336 h3->SetName("Result");
337 h3->SetBit(kCanDelete);
338 h3->SetDirectory(0);
339 h3->Draw();
340
341 h1->Draw("same");
342
343 // --------------------------
344
345 pad->cd();
346 AppendPad("result");
347}
348
349void MHPhi::PaintUpdate() const
350{
351 TVirtualPad *pad1 = gPad->GetPad(1);
352 TVirtualPad *pad2 = gPad->GetPad(2);
353 TVirtualPad *pad3 = gPad->GetPad(3);
354 TVirtualPad *pad4 = gPad->GetPad(4);
355
356 Double_t sig2 = 0;
357 Double_t bg2 = 0;
358 Double_t f2 = 1;
359 TH1D *htemp = pad1 ? dynamic_cast<TH1D*>(pad1->FindObject("Template")) : NULL;
360 if (htemp)
361 {
362 fHTemplate.Copy(*htemp);
363 htemp->SetName("Template");
364 htemp->SetDirectory(0);
365
366 Double_t sc1 = 1;
367 Double_t sc2 = 1;
368
369 TH1D *res = pad4 ? dynamic_cast<TH1D*>(pad4->FindObject("Result")) : NULL;
370 if (res)
371 {
372 fHPhi.Copy(*res);
373
374 // Scale inhomogeneity to match the phi-plot in the off-region
375 sc1 = res->Integral(fNumBinsSignal+1, 9999)/fHInhom.Integral(fNumBinsSignal+1, 9999);
376 // Scale inhomogeneity to match the phi-plot in the off-region
377 sc2 = fHInhom.Integral()/htemp->Integral();
378
379 res->Add(&fHInhom, -sc1);
380 res->Add(htemp, sc1*sc2);
381 res->SetName("Result");
382 res->SetDirectory(0);
383
384 htemp->Scale(res->Integral(fNumBinsSignal+1, 9999)/htemp->Integral(fNumBinsSignal+1, 9999));
385
386 sig2 = res->Integral(1, fNumBinsSignal);
387 bg2 = fHPhi.Integral(fNumBinsSignal+1, 9999);
388 f2 = htemp->Integral(1, fNumBinsSignal)/htemp->Integral(fNumBinsSignal+1, 9999);
389 }
390
391 TH1D *hinhom = pad3 ? dynamic_cast<TH1D*>(pad3->FindObject("Inhomogeneity")) : NULL;
392 if (hinhom)
393 {
394 fHInhom.Copy(*hinhom);
395 hinhom->SetName("Inhomogeneity");
396 hinhom->SetDirectory(0);
397 hinhom->Scale(sc1);
398 }
399 }
400
401 htemp = pad2 ? dynamic_cast<TH1D*>(pad2->FindObject("Template")) : NULL;
402 if (htemp)
403 {
404 fHTemplate.Copy(*htemp);
405 htemp->Scale(fHInhom.Integral()/htemp->Integral());
406 htemp->SetName("Template");
407 htemp->SetDirectory(0);
408 }
409}
410
411void MHPhi::PaintText(const TH1D &res) const
412{
413 const Double_t cut = res.GetBinLowEdge(fNumBinsSignal+1);
414
415 const Double_t sig = res.Integral(1, fNumBinsSignal);
416 const Double_t bg = res.Integral(fNumBinsSignal+1, 9999);
417
418 const Double_t f = fHTemplate.Integral(1, fNumBinsSignal)/fHTemplate.Integral(fNumBinsSignal+1, 9999);
419
420 const Double_t S0 = MMath::SignificanceLiMaSigned(sig, bg*f);
421 const Double_t S = MMath::SignificanceLiMaSigned(sig, bg, f);
422
423 const TString fmt = Form("\\sigma_{L/M}=%.1f (\\sigma_{0}=%.1f) \\Delta\\Phi_{on}<%.1f\\circ E=%.0f B=%.0f f=%.2f",
424 S, S0, cut, sig-bg*f, bg*f, f);
425
426 const Double_t b = bg *f/fNumBinsSignal;
427 const Double_t e = TMath::Sqrt(bg)*f/fNumBinsSignal;
428
429 TLatex text(0.275, 0.95, fmt);
430 text.SetBit(TLatex::kTextNDC);
431 text.SetTextSize(0.042);
432 text.Paint();
433
434 TLine line;
435 line.SetLineColor(14);
436 line.PaintLine(cut, gPad->GetUymin(), cut, gPad->GetUymax());
437
438 // Here we calculate an upper phi cut to take a
439 // possible anti-theta cut into account
440 const Double_t ulim = fUseAntiPhiCut ? 180-fHPhi.GetBinLowEdge(fNumBinsSignal+1)*1.1 : 180;
441 line.SetLineStyle(kDotted);
442 line.PaintLine(ulim, gPad->GetUymin(), ulim, gPad->GetUymax());
443
444 line.SetLineStyle(kSolid);
445 line.SetLineColor(kBlack);
446 line.PaintLine(0, b, cut, b);
447 line.PaintLine(cut/2, b-e, cut/2, b+e);
448 line.SetLineStyle(kDashed);
449 line.PaintLine(cut, b, fHPhi.GetXaxis()->GetXmax(), b);
450
451 TMarker m;
452 m.SetMarkerStyle(kFullDotMedium);
453 m.PaintMarker(cut/2, b);
454}
455
456void MHPhi::PaintResult() const
457{
458 TVirtualPad *pad = gPad;
459
460 pad->cd(1);
461 PaintText(fHPhi);
462
463 pad->cd(4);
464 TH1D *res = gPad ? dynamic_cast<TH1D*>(gPad->FindObject("Result")) : NULL;
465 if (res)
466 PaintText(*res);
467}
468
469void MHPhi::Paint(Option_t *o)
470{
471 TString opt(o);
472 if (opt=="update")
473 PaintUpdate();
474 if (opt=="result")
475 PaintResult();
476}
477
478Int_t MHPhi::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
479{
480 Bool_t rc = kFALSE;
481 if (IsEnvDefined(env, prefix, "NumBinsSignal", print))
482 {
483 SetNumBinsSignal(GetEnvValue(env, prefix, "NumBinsSignal", (Int_t)fNumBinsSignal));
484 rc = kTRUE;
485 }
486 if (IsEnvDefined(env, prefix, "UseAntiPhiCut", print))
487 {
488 SetUseAntiPhiCut(GetEnvValue(env, prefix, "UseAntiPhiCut", (Int_t)fUseAntiPhiCut));
489 rc = kTRUE;
490 }
491
492 return rc;
493}
Note: See TracBrowser for help on using the repository browser.