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

Last change on this file since 8634 was 7618, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 14.7 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 AppendPad("result1");
300
301 // --------------------------
302
303 pad->cd(2);
304 gPad->SetBorderMode(0);
305 gPad->SetFrameBorderMode(0);
306
307 fHInhom.Draw();
308 fHInhomOff.Draw("same");
309
310 TH1D *h2 = new TH1D(fHTemplate);
311 h2->SetName("Template");
312 h2->SetBit(kCanDelete);
313 h2->SetDirectory(0);
314 h2->Draw("same");
315
316 // --------------------------
317
318 pad->cd(3);
319 gPad->SetBorderMode(0);
320 gPad->SetFrameBorderMode(0);
321
322 fHPhi.Draw();
323 TH1D *h4 = new TH1D(fHInhom);
324 h4->SetName("Inhomogeneity");
325 h4->SetBit(kCanDelete);
326 h4->SetDirectory(0);
327 h4->Draw("same");
328
329 h1->Draw("same");
330
331 // --------------------------
332
333 pad->cd(4);
334 gPad->SetBorderMode(0);
335 gPad->SetFrameBorderMode(0);
336
337 TH1D *h3 = new TH1D(fHPhi);
338 h3->SetName("Result");
339 h3->SetBit(kCanDelete);
340 h3->SetDirectory(0);
341 h3->Draw();
342
343 h1->Draw("same");
344
345 AppendPad("result4");
346
347 // --------------------------
348}
349
350void MHPhi::PaintUpdate() const
351{
352 TVirtualPad *pad1 = gPad->GetPad(1);
353 TVirtualPad *pad2 = gPad->GetPad(2);
354 TVirtualPad *pad3 = gPad->GetPad(3);
355 TVirtualPad *pad4 = gPad->GetPad(4);
356
357 Double_t sig2 = 0;
358 Double_t bg2 = 0;
359 Double_t f2 = 1;
360 TH1D *htemp = pad1 ? dynamic_cast<TH1D*>(pad1->FindObject("Template")) : NULL;
361 if (htemp)
362 {
363 fHTemplate.Copy(*htemp);
364 htemp->SetName("Template");
365 htemp->SetDirectory(0);
366
367 Double_t sc1 = 1;
368 Double_t sc2 = 1;
369
370 TH1D *res = pad4 ? dynamic_cast<TH1D*>(pad4->FindObject("Result")) : NULL;
371 if (res)
372 {
373 fHPhi.Copy(*res);
374
375 // Scale inhomogeneity to match the phi-plot in the off-region
376 sc1 = res->Integral(fNumBinsSignal+1, 9999)/fHInhom.Integral(fNumBinsSignal+1, 9999);
377 // Scale inhomogeneity to match the phi-plot in the off-region
378 sc2 = fHInhom.Integral()/htemp->Integral();
379
380 res->Add(&fHInhom, -sc1);
381 res->Add(htemp, sc1*sc2);
382 res->SetName("Result");
383 res->SetDirectory(0);
384
385 htemp->Scale(res->Integral(fNumBinsSignal+1, 9999)/htemp->Integral(fNumBinsSignal+1, 9999));
386
387 sig2 = res->Integral(1, fNumBinsSignal);
388 bg2 = fHPhi.Integral(fNumBinsSignal+1, 9999);
389 f2 = htemp->Integral(1, fNumBinsSignal)/htemp->Integral(fNumBinsSignal+1, 9999);
390 }
391
392 TH1D *hinhom = pad3 ? dynamic_cast<TH1D*>(pad3->FindObject("Inhomogeneity")) : NULL;
393 if (hinhom)
394 {
395 fHInhom.Copy(*hinhom);
396 hinhom->SetName("Inhomogeneity");
397 hinhom->SetDirectory(0);
398 hinhom->Scale(sc1);
399 }
400 }
401
402 htemp = pad2 ? dynamic_cast<TH1D*>(pad2->FindObject("Template")) : NULL;
403 if (htemp)
404 {
405 fHTemplate.Copy(*htemp);
406 htemp->Scale(fHInhom.Integral()/htemp->Integral());
407 htemp->SetName("Template");
408 htemp->SetDirectory(0);
409 }
410}
411
412void MHPhi::PaintText(const TH1D &res) const
413{
414 const Double_t cut = res.GetBinLowEdge(fNumBinsSignal+1);
415
416 const Double_t sig = res.Integral(1, fNumBinsSignal);
417 const Double_t bg = res.Integral(fNumBinsSignal+1, 9999);
418
419 const Double_t f = fHTemplate.Integral(1, fNumBinsSignal)/fHTemplate.Integral(fNumBinsSignal+1, 9999);
420
421 const Double_t S0 = MMath::SignificanceLiMaSigned(sig, bg*f);
422 const Double_t S = MMath::SignificanceLiMaSigned(sig, bg, f);
423
424 const TString fmt = Form("\\sigma_{L/M}=%.1f (\\sigma_{0}=%.1f) \\Delta\\Phi_{on}<%.1f\\circ E=%.0f B=%.0f f=%.2f",
425 S, S0, cut, sig-bg*f, bg*f, f);
426
427 const Double_t b = bg *f/fNumBinsSignal;
428 const Double_t e = TMath::Sqrt(bg)*f/fNumBinsSignal;
429
430 TLatex text(0.275, 0.95, fmt);
431 text.SetBit(TLatex::kTextNDC);
432 text.SetTextSize(0.042);
433 text.Paint();
434
435 TLine line;
436 line.SetLineColor(14);
437 line.PaintLine(cut, gPad->GetUymin(), cut, gPad->GetUymax());
438
439 // Here we calculate an upper phi cut to take a
440 // possible anti-theta cut into account
441 const Double_t ulim = fUseAntiPhiCut ? 180-fHPhi.GetBinLowEdge(fNumBinsSignal+1)*1.1 : 180;
442 line.SetLineStyle(kDotted);
443 line.PaintLine(ulim, gPad->GetUymin(), ulim, gPad->GetUymax());
444
445 line.SetLineStyle(kSolid);
446 line.SetLineColor(kBlack);
447 line.PaintLine(0, b, cut, b);
448 line.PaintLine(cut/2, b-e, cut/2, b+e);
449 line.SetLineStyle(kDashed);
450 line.PaintLine(cut, b, fHPhi.GetXaxis()->GetXmax(), b);
451
452 TMarker m;
453 m.SetMarkerStyle(kFullDotMedium);
454 m.PaintMarker(cut/2, b);
455}
456
457void MHPhi::Paint(Option_t *o)
458{
459 TString opt(o);
460 if (opt=="update")
461 PaintUpdate();
462 if (opt=="result1")
463 PaintText(fHPhi);
464 if (opt=="result4")
465 {
466 TH1D *res = gPad ? dynamic_cast<TH1D*>(gPad->FindObject("Result")) : NULL;
467 if (res)
468 PaintText(*res);
469 }
470}
471
472Int_t MHPhi::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
473{
474 Bool_t rc = kFALSE;
475 if (IsEnvDefined(env, prefix, "NumBinsSignal", print))
476 {
477 SetNumBinsSignal(GetEnvValue(env, prefix, "NumBinsSignal", (Int_t)fNumBinsSignal));
478 rc = kTRUE;
479 }
480 if (IsEnvDefined(env, prefix, "UseAntiPhiCut", print))
481 {
482 SetUseAntiPhiCut(GetEnvValue(env, prefix, "UseAntiPhiCut", (Int_t)fUseAntiPhiCut));
483 rc = kTRUE;
484 }
485
486 return rc;
487}
Note: See TracBrowser for help on using the repository browser.