source: trunk/MagicSoft/Mars/mhflux/MHDisp.cc@ 7115

Last change on this file since 7115 was 7115, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 12.0 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, 5/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MHDisp
28//
29// Create a false source plot using disp.
30//
31// Currently the use of this class requires to be after MFMagicCuts
32// in the tasklist. Switching of the M3L cut in MFMagicCuts is recommended.
33//
34//////////////////////////////////////////////////////////////////////////////
35#include "MHDisp.h"
36
37#include <TStyle.h>
38#include <TCanvas.h>
39
40#include <TF1.h>
41#include <TF2.h>
42#include <TProfile.h>
43
44#include "MLog.h"
45#include "MLogManip.h"
46
47#include "MMath.h"
48
49#include "MParList.h"
50#include "MParameters.h"
51
52#include "MHillasExt.h"
53#include "MHillasSrc.h"
54#include "MSrcPosCam.h"
55#include "MPointingPos.h"
56
57ClassImp(MHDisp);
58
59using namespace std;
60
61// --------------------------------------------------------------------------
62//
63// Default Constructor
64//
65MHDisp::MHDisp(const char *name, const char *title)
66 : fHilExt(0), fDisp(0)
67{
68 //
69 // set the name and title of this object
70 //
71 fName = name ? name : "MHDisp";
72 fTitle = title ? title : "3D-plot using Disp vs x, y";
73
74 fHist.SetDirectory(NULL);
75
76 fHist.SetName("Alpha");
77 fHist.SetTitle("3D-plot of ThetaSq vs x, y");
78 fHist.SetXTitle("x [\\circ]");
79 fHist.SetYTitle("y [\\circ]");
80 fHist.SetZTitle("\\vartheta [deg^{2}]");
81}
82
83// --------------------------------------------------------------------------
84//
85// Set binnings (takes BinningFalseSource) and prepare filling of the
86// histogram.
87//
88// Also search for MTime, MObservatory and MPointingPos
89//
90Bool_t MHDisp::SetupFill(const MParList *plist)
91{
92 if (!MHFalseSource::SetupFill(plist))
93 return kFALSE;
94
95 fDisp = (MParameterD*)plist->FindObject("Disp", "MParameterD");
96 if (!fDisp)
97 {
98 *fLog << err << "Disp [MParameterD] not found... abort." << endl;
99 return kFALSE;
100 }
101
102 MParameterD *p = (MParameterD*)plist->FindObject("M3lCut", "MParameterD");
103 if (!p)
104 {
105 *fLog << err << "M3lCut [MParameterD] not found... abort." << endl;
106 return kFALSE;
107 }
108 fM3lCut = TMath::Abs(p->GetVal());
109
110 p = (MParameterD*)plist->FindObject("DispXi", "MParameterD");
111 if (!p)
112 {
113 *fLog << err << "DispXi [MParameterD] not found... abort." << endl;
114 return kFALSE;
115 }
116 fXi = p->GetVal();
117
118 p = (MParameterD*)plist->FindObject("DispXiTheta", "MParameterD");
119 if (!p)
120 {
121 *fLog << err << "DispXiTheta [MParameterD] not found... abort." << endl;
122 return kFALSE;
123 }
124 fXiTheta = p->GetVal();
125
126 fHilExt = (MHillasExt*)plist->FindObject("MHillasExt");
127 if (!fHilExt)
128 {
129 *fLog << err << "MHillasExt not found... aborting." << endl;
130 return kFALSE;
131 }
132
133 // Initialize all bins with a small (=0) value otherwise
134 // these bins are not displayed
135 for (int x=0; x<fHist.GetNbinsX(); x++)
136 for (int y=0; y<fHist.GetNbinsX(); y++)
137 fHist.Fill(fHist.GetXaxis()->GetBinCenter(x+1),
138 fHist.GetYaxis()->GetBinCenter(y+1),
139 0.0, 1e-10);
140
141 return kTRUE;
142}
143
144// --------------------------------------------------------------------------
145//
146// Fill the histogram. For details see the code or the class description
147//
148Bool_t MHDisp::Fill(const MParContainer *par, const Stat_t w)
149{
150 const MHillas *hil = dynamic_cast<const MHillas*>(par);
151 if (!hil)
152 {
153 *fLog << err << "MHDisp::Fill: No container specified!" << endl;
154 return kFALSE;
155 }
156
157 // Get camera rotation angle
158 Double_t rho = 0;
159 if (fTime && fObservatory && fPointPos)
160 rho = fPointPos->RotationAngle(*fObservatory, *fTime);
161
162 // Get Disp from Parlist
163 const Double_t disp = fDisp->GetVal();
164
165 // Calculate where disp is pointing
166 TVector2 pos1(hil->GetCosDelta(), hil->GetSinDelta());
167 TVector2 pos2(hil->GetCosDelta(), hil->GetSinDelta());
168 pos1 *= disp;
169 pos2 *= -disp;
170
171 pos1 += hil->GetMean()*fMm2Deg;
172 pos2 += hil->GetMean()*fMm2Deg;
173
174 // gammaweight: If we couldn't decide which position makes the
175 // event a gamma, both position are assigned 'half of a gamma'
176 Double_t gweight = 0.5;
177
178 // Check whether our g/h-separation allows to asign definitly
179 // to one unique position. Therefor we requeire that the event
180 // would be considered a gamma for one, but not for the other
181 // position. This can only be the case if the third moment
182 // has a value higher than the absolute cut value.
183 if (TMath::Abs(fHilExt->GetM3Long()) > fM3lCut)
184 {
185 // Because at one position the event is considered a gamma
186 // we have to find out which position it is...
187 MSrcPosCam src;
188 MHillasSrc hsrc;
189 hsrc.SetSrcPos(&src);
190
191 // Calculate the sign for one of the desired source positions
192 // The other position must have the opposite sign
193 src.SetXY(pos1/fMm2Deg);
194 if (hsrc.Calc(*hil)>0)
195 {
196 *fLog << warn << "Calculation of MHillasSrc failed." << endl;
197 return kFALSE;
198 }
199 const Double_t m3l = fHilExt->GetM3Long()*TMath::Sign(1.0f, hsrc.GetCosDeltaAlpha())*fMm2Deg;
200
201 gweight = m3l>fM3lCut ? 1 : 0;
202 }
203
204 // Now we can safly derotate both position...
205 TVector2 srcp(fSrcPos->GetXY());
206 if (rho!=0)
207 {
208 pos1=pos1.Rotate(-rho);
209 pos2=pos2.Rotate(-rho);
210 srcp=srcp.Rotate(-rho);
211 }
212
213 if (fSrcPos)
214 {
215 pos1 -= srcp*fMm2Deg;
216 pos2 -= srcp*fMm2Deg;
217 }
218
219 // Workaround: Number-of-entries
220 if (gweight>0.25)
221 fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*gweight);
222 if (gweight<0.75)
223 fHist.Fill(pos2.X(), pos2.Y(), 0.0, w*(1-gweight));
224
225 return kTRUE;
226}
227/*
228static Double_t FcnGauss2d(Double_t *x, Double_t *par)
229{
230 TVector2 v = TVector2(x[0], x[1]).Rotate(par[5]*TMath::DegToRad());
231
232 const Double_t g0 = TMath::Gaus(v.X(), par[1], par[2]);
233 const Double_t g1 = TMath::Gaus(v.Y(), par[3], par[4]);
234
235 //Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE);
236 return par[0]*g0*g1 + par[6]*(v.X()*v.X() + v.Y()*v.Y()) + par[7];
237}*/
238
239void MHDisp::Paint(Option_t *o)
240{
241 TVirtualPad *pad = gPad;
242
243 pad->cd(1);
244
245 fHist.GetZaxis()->SetRange(0,9999);
246 TH1 *h1=fHist.Project3D("yx_on");
247 gStyle->SetPalette(1, 0);
248
249 Int_t ix, iy, iz;
250 TF1 func("fcn", "gaus + [3]*x*x + [4]");
251
252 pad->cd(3);
253 TH1 *h2 = (TH1*)gPad->FindObject("RadProf");
254 /*
255 if (!h2)
256 {
257 const Double_t maxr = TMath::Hypot(h1->GetXaxis()->GetXmax(), h1->GetYaxis()->GetXmax());
258 const Int_t nbin = (h1->GetNbinsX()+h1->GetNbinsY())/2;
259 TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr);
260 //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr);
261 h->Sumw2();
262 h->SetXTitle("\\vartheta [\\circ]");
263 h->SetYTitle("<cts>/\\Delta R");
264 h->SetBit(kCanDelete);
265 h->Draw();
266 }*/
267 if (h1 && h2)
268 {
269 h2->Reset();
270
271 h1->GetMaximumBin(ix, iy, iz);
272
273 const Double_t x0 = h1->GetXaxis()->GetBinCenter(ix);
274 const Double_t y0 = h1->GetYaxis()->GetBinCenter(iy);
275
276 for (int x=0; x<h1->GetNbinsX(); x++)
277 for (int y=0; y<h1->GetNbinsY(); y++)
278 {
279 const Double_t r = TMath::Hypot(h1->GetXaxis()->GetBinCenter(x+1)-x0,
280 h1->GetYaxis()->GetBinCenter(y+1)-y0);
281 h2->Fill(r, h1->GetBinContent(x+1, y+1));
282 }
283
284 // Replace with MAlphaFitter?
285 func.SetLineWidth(1);
286 func.SetLineColor(kBlue);
287
288 func.SetParLimits(2, h2->GetBinWidth(1), 1.0);
289
290 func.SetParameter(0, h2->GetBinContent(1));
291 func.FixParameter(1, 0);
292 func.SetParameter(2, 0.15);
293 func.SetParameter(4, h2->GetBinContent(10));
294 h2->Fit(&func, "IMQ", "", 0, 1.0);
295
296 const Double_t r0 = 2*func.GetParameter(2);
297 const Double_t e = func.Integral(0, r0)/h1->GetBinWidth(1);
298 func.SetParameter(0, 0);
299 const Double_t b = func.Integral(0, r0)/h1->GetBinWidth(1);
300 const Double_t s = MMath::SignificanceLiMa(e, b);
301
302 h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ \\sigma=%.1f E=%.0f B=%.0f", x0, y0, func.GetParameter(2), s, e-b, b));
303 }
304 /*
305 if (h1)
306 {
307 const Double_t maxr = 0.9*TMath::Abs(fHist.GetXaxis()->GetXmax());
308
309 TF2 f2d("Gaus2D", FcnGauss2d, -maxr, maxr, -maxr, maxr, 8);
310 f2d.SetLineWidth(1);
311
312 f2d.SetParameter(0, h1->GetMaximum()*5); // A
313 f2d.SetParLimits(1, h1->GetXaxis()->GetBinCenter(ix)-h1->GetXaxis()->GetBinWidth(ix)*5, h1->GetXaxis()->GetBinCenter(ix)+h1->GetXaxis()->GetBinWidth(ix)); // mu_1
314 f2d.SetParLimits(3, h1->GetYaxis()->GetBinCenter(iy)-h1->GetYaxis()->GetBinWidth(iy)*5, h1->GetYaxis()->GetBinCenter(iy)+h1->GetYaxis()->GetBinWidth(iy)); // mu_2
315 f2d.SetParLimits(2, 0, func.GetParameter(2)*5); // sigma_1
316 f2d.SetParLimits(4, 0, func.GetParameter(2)*5); // sigma_2
317 f2d.SetParLimits(5, 0, 45); // phi
318 f2d.SetParLimits(6, 0, func.GetParameter(3)*5);
319 f2d.SetParLimits(7, 0, func.GetParameter(4)*5);
320
321 f2d.SetParameter(0, h1->GetMaximum()); // A
322 f2d.SetParameter(1, h1->GetXaxis()->GetBinCenter(ix)); // mu_1
323 f2d.SetParameter(2, func.GetParameter(2)); // sigma_1
324 f2d.SetParameter(3, h1->GetYaxis()->GetBinCenter(iy)); // mu_2
325 f2d.SetParameter(4, func.GetParameter(2)); // sigma_2
326 f2d.FixParameter(5, 0); // phi
327 f2d.SetParameter(6, func.GetParameter(3));
328 f2d.SetParameter(7, func.GetParameter(4));
329 h1->Fit(&f2d, "Q", "cont2");
330 //f2d.DrawCopy("cont2same");
331 }*/
332}
333
334void MHDisp::Draw(Option_t *o)
335{
336 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
337 const Int_t col = pad->GetFillColor();
338 pad->SetBorderMode(0);
339
340 AppendPad("");
341
342 TString name = Form("%s_1", pad->GetName());
343 TPad *p = new TPad(name,name, 0.005, 0.005, 0.65, 0.995, col, 0, 0);
344 p->SetNumber(1);
345 p->Draw();
346 p->cd();
347
348 TH1 *h3 = fHist.Project3D("yx_on");
349 h3->SetTitle("Distribution of equivalent events");
350 h3->SetDirectory(NULL);
351 h3->SetXTitle(fHist.GetXaxis()->GetTitle());
352 h3->SetYTitle(fHist.GetYaxis()->GetTitle());
353 h3->SetMinimum(0);
354 h3->Draw("colz");
355 h3->SetBit(kCanDelete);
356// catalog->Draw("mirror same *");
357
358 pad->cd();
359 name = Form("%s_2", pad->GetName());
360 p = new TPad(name,name, 0.66, 0.005, 0.995, 0.5, col, 0, 0);
361 p->SetNumber(2);
362 p->Draw();
363 p->cd();
364 h3->Draw("surf3");
365
366 pad->cd();
367 name = Form("%s_3", pad->GetName());
368 p = new TPad(name,name, 0.66, 0.5, 0.995, 0.995, col, 0, 0);
369 p->SetNumber(3);
370 p->Draw();
371 p->cd();
372
373 const Double_t maxr = TMath::Hypot(h3->GetXaxis()->GetXmax(), h3->GetYaxis()->GetXmax());
374 const Int_t nbin = (h3->GetNbinsX()+h3->GetNbinsY())/2;
375 TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr);
376 h->SetDirectory(0);
377 //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr);
378 //h->Sumw2();
379 h->SetXTitle("\\vartheta [\\circ]");
380 h->SetYTitle("<cts>/\\Delta R");
381 h->SetBit(kCanDelete);
382 h->Draw();
383}
Note: See TracBrowser for help on using the repository browser.