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

Last change on this file since 7178 was 7173, checked in by tbretz, 20 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, 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)//, fIsWobble(kFALSE)
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("Eq.cts");
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 //const Double_t xmax = fHist.GetXaxis()->GetXmax();
134
135 // Initialize all bins with a small (=0) value otherwise
136 // these bins are not displayed
137 for (int x=0; x<fHist.GetNbinsX(); x++)
138 for (int y=0; y<fHist.GetNbinsY(); y++)
139 {
140 const Double_t cx = fHist.GetXaxis()->GetBinCenter(x+1);
141 const Double_t cy = fHist.GetYaxis()->GetBinCenter(y+1);
142 //if (TMath::Hypot(cx, cy)<xmax)
143 fHist.Fill(cx, cy, 0.0, 1e-10);
144 }
145
146 return kTRUE;
147}
148
149// --------------------------------------------------------------------------
150//
151// Fill the histogram. For details see the code or the class description
152//
153Bool_t MHDisp::Fill(const MParContainer *par, const Stat_t w)
154{
155 const MHillas *hil = dynamic_cast<const MHillas*>(par);
156 if (!hil)
157 {
158 *fLog << err << "MHDisp::Fill: No container specified!" << endl;
159 return kFALSE;
160 }
161
162 // Get camera rotation angle
163 Double_t rho = 0;
164 if (fTime && fObservatory && fPointPos)
165 rho = fPointPos->RotationAngle(*fObservatory, *fTime);
166
167 // FIXME: Do wobble-rotation when subtracting?
168 if (!fHistOff/* && fIsWobble*/)
169 rho += TMath::Pi();
170
171 // Get Disp from Parlist
172 const Double_t disp = fDisp->GetVal();
173
174 // Calculate both postiions where disp could point
175 TVector2 pos1(hil->GetCosDelta(), hil->GetSinDelta());
176 TVector2 pos2(hil->GetCosDelta(), hil->GetSinDelta());
177 pos1 *= disp; // Vector of length disp in direction of shower
178 pos2 *= -disp; // Vector of length -disp in direction of shower
179
180 // Move origin of vector to center-of-gravity of shower
181 pos1 += hil->GetMean()*fMm2Deg;
182 pos2 += hil->GetMean()*fMm2Deg;
183
184 // gammaweight: If we couldn't decide which position makes the
185 // event a gamma, both position are assigned 'half of a gamma'
186 Double_t gweight = 0.5;
187
188 // Check whether our g/h-separation allows to asign definitly
189 // to one unique position. Therefor we requeire that the event
190 // would be considered a gamma for one, but not for the other
191 // position. This can only be the case if the third moment
192 // has a value higher than the absolute cut value.
193 if (TMath::Abs(fHilExt->GetM3Long()*fMm2Deg) > fM3lCut)
194 {
195 // Because at one position the event is considered a gamma
196 // we have to find out which position it is...
197 MSrcPosCam src;
198 MHillasSrc hsrc;
199 hsrc.SetSrcPos(&src);
200
201 // Calculate the sign for one of the desired source positions
202 // The other position must have the opposite sign
203 src.SetXY(pos1/fMm2Deg);
204 if (hsrc.Calc(*hil)>0)
205 {
206 *fLog << warn << "Calculation of MHillasSrc failed." << endl;
207 return kFALSE;
208 }
209 const Double_t m3l = fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, hsrc.GetCosDeltaAlpha());
210
211 gweight = m3l>fM3lCut ? 1 : 0;
212 }
213
214 // Now we can safly derotate both position...
215 TVector2 srcp;
216 if (fSrcPos)
217 srcp = fSrcPos->GetXY();
218
219 // Derotate all position around camera center by -rho
220 if (rho!=0)
221 {
222 pos1=pos1.Rotate(-rho);
223 pos2=pos2.Rotate(-rho);
224 srcp=srcp.Rotate(-rho);
225 }
226
227 // Shift the source position to 0/0
228 pos1 -= srcp*fMm2Deg;
229 pos2 -= srcp*fMm2Deg;
230
231 //const Double_t xmax = fHist.GetXaxis()->GetXmax();
232
233 // Fill histograms
234 //if (pos1.Mod()<xmax)
235 fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*gweight);
236 //if (pos2.Mod()<xmax)
237 fHist.Fill(pos2.X(), pos2.Y(), 0.0, w*(1-gweight));
238
239 return kTRUE;
240}
241
242Double_t MHDisp::GetOffSignal(TH1 &h) const
243{
244 const TAxis &axex = *h.GetXaxis();
245 const TAxis &axey = *h.GetYaxis();
246
247 Double_t sum = 0;
248 for (int x=0; x<h.GetNbinsX(); x++)
249 for (int y=0; y<h.GetNbinsY(); y++)
250 {
251 if (TMath::Hypot(axex.GetBinCenter(x+1), axey.GetBinCenter(y+1))>0.35)
252 sum += h.GetBinContent(x+1,y+1);
253 }
254
255 return sum;
256}
257
258void MHDisp::Paint(Option_t *o)
259{
260 TVirtualPad *pad = gPad;
261
262 pad->cd(1);
263
264 // Project on data onto yx-plane
265 fHist.GetZaxis()->SetRange(0,9999);
266 TH1 *h1=fHist.Project3D("yx_on");
267
268 // Set Glow-palette (PaintSurface doesn't allow more than 99 colors)
269 MH::SetPalette(fHistOff?"glowsym":"glow1", 99);
270 h1->SetContour(99);
271
272 Double_t scale = 1;
273 if (fHistOff)
274 {
275 // Project off data onto yx-plane and subtract it from on-data
276 fHistOff->GetZaxis()->SetRange(0,9999);
277 TH1 *h=fHistOff->Project3D("yx_off");
278
279 scale = -1;
280
281 //if (!fIsWobble)
282 {
283 const Double_t h1off = GetOffSignal(*h1);
284 const Double_t hoff = GetOffSignal(*h);
285 scale = hoff==0?-1:-h1off/hoff;
286 }
287
288 h1->Add(h, scale);
289 delete h;
290
291 // Force calculation of minimum, maximum
292 h1->SetMinimum();
293 h1->SetMaximum();
294
295 // Find maximum
296 const Double_t max = TMath::Max(TMath::Abs(h1->GetMinimum()),
297 TMath::Abs(h1->GetMaximum()));
298
299 // Set new minimum, maximum
300 h1->SetMinimum(-max);
301 h1->SetMaximum( max);
302 }
303
304 Int_t ix, iy, iz;
305 TF1 func("fcn", "gaus + [3]*x*x + [4]");
306
307 pad->cd(3);
308 TH1 *h2 = (TH1*)gPad->FindObject("RadProf");
309 /*
310 if (!h2)
311 {
312 const Double_t maxr = TMath::Hypot(h1->GetXaxis()->GetXmax(), h1->GetYaxis()->GetXmax());
313 const Int_t nbin = (h1->GetNbinsX()+h1->GetNbinsY())/2;
314 TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr);
315 //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr);
316 h->Sumw2();
317 h->SetXTitle("\\vartheta [\\circ]");
318 h->SetYTitle("<cts>/\\Delta R");
319 h->SetBit(kCanDelete);
320 h->Draw();
321 }*/
322 if (h1 && h2)
323 {
324 h2->Reset();
325
326 h1->GetMaximumBin(ix, iy, iz);
327
328 const Double_t x0 = h1->GetXaxis()->GetBinCenter(ix);
329 const Double_t y0 = h1->GetYaxis()->GetBinCenter(iy);
330 const Double_t w0 = h1->GetXaxis()->GetBinWidth(1);
331
332 for (int x=0; x<h1->GetNbinsX(); x++)
333 for (int y=0; y<h1->GetNbinsY(); y++)
334 {
335 const Double_t r = TMath::Hypot(h1->GetXaxis()->GetBinCenter(x+1)-x0,
336 h1->GetYaxis()->GetBinCenter(y+1)-y0);
337 h2->Fill(r, h1->GetBinContent(x+1, y+1));
338 }
339
340 // Replace with MAlphaFitter?
341 func.SetLineWidth(1);
342 func.SetLineColor(kBlue);
343
344 func.SetParLimits(2, h2->GetBinWidth(1), 1.0);
345
346 func.SetParameter(0, h2->GetBinContent(1));
347 func.FixParameter(1, 0);
348 func.SetParameter(2, 0.15);
349 func.SetParameter(4, h2->GetBinContent(10));
350 h2->Fit(&func, "IMQ", "", 0, 1.0);
351
352 // No wintegrate the function f(x) per Delta Area
353 // which is f(x)/(pi*delta r*(2*r+delta r))
354 TF1 func2("fcn2", Form("(gaus + [3]*x*x + [4])/(2*x+%.5f)", w0));
355 for (int i=0; i<5; i++)
356 func2.SetParameter(i, func.GetParameter(i));
357
358 const Double_t r0 = 2*func2.GetParameter(2);
359 const Double_t e = func2.Integral(0, r0)/(w0*TMath::Pi());
360 func2.SetParameter(0, 0);
361 const Double_t b = func2.Integral(0, r0)/(w0*TMath::Pi());
362 const Double_t s = MMath::SignificanceLiMa(e, b);
363
364 h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ \\sigma=%.1f E=%.0f B=%.0f f=%.2f", x0, y0, func.GetParameter(2), s, e-b, b, scale));
365 }
366 /*
367 if (h1)
368 {
369 const Double_t maxr = 0.9*TMath::Abs(fHist.GetXaxis()->GetXmax());
370
371 TF2 f2d("Gaus2D", FcnGauss2d, -maxr, maxr, -maxr, maxr, 8);
372 f2d.SetLineWidth(1);
373
374 f2d.SetParameter(0, h1->GetMaximum()*5); // A
375 f2d.SetParLimits(1, h1->GetXaxis()->GetBinCenter(ix)-h1->GetXaxis()->GetBinWidth(ix)*5, h1->GetXaxis()->GetBinCenter(ix)+h1->GetXaxis()->GetBinWidth(ix)); // mu_1
376 f2d.SetParLimits(3, h1->GetYaxis()->GetBinCenter(iy)-h1->GetYaxis()->GetBinWidth(iy)*5, h1->GetYaxis()->GetBinCenter(iy)+h1->GetYaxis()->GetBinWidth(iy)); // mu_2
377 f2d.SetParLimits(2, 0, func.GetParameter(2)*5); // sigma_1
378 f2d.SetParLimits(4, 0, func.GetParameter(2)*5); // sigma_2
379 f2d.SetParLimits(5, 0, 45); // phi
380 f2d.SetParLimits(6, 0, func.GetParameter(3)*5);
381 f2d.SetParLimits(7, 0, func.GetParameter(4)*5);
382
383 f2d.SetParameter(0, h1->GetMaximum()); // A
384 f2d.SetParameter(1, h1->GetXaxis()->GetBinCenter(ix)); // mu_1
385 f2d.SetParameter(2, func.GetParameter(2)); // sigma_1
386 f2d.SetParameter(3, h1->GetYaxis()->GetBinCenter(iy)); // mu_2
387 f2d.SetParameter(4, func.GetParameter(2)); // sigma_2
388 f2d.FixParameter(5, 0); // phi
389 f2d.SetParameter(6, func.GetParameter(3));
390 f2d.SetParameter(7, func.GetParameter(4));
391 h1->Fit(&f2d, "Q", "cont2");
392 //f2d.DrawCopy("cont2same");
393 }*/
394}
395
396void MHDisp::Draw(Option_t *o)
397{
398 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
399 const Int_t col = pad->GetFillColor();
400 pad->SetBorderMode(0);
401
402 AppendPad("");
403
404 TString name = Form("%s_1", pad->GetName());
405 TPad *p = new TPad(name,name, 0.005, 0.005, 0.65, 0.995, col, 0, 0);
406 p->SetNumber(1);
407 p->Draw();
408 p->cd();
409
410 TH1 *h3 = fHist.Project3D("yx_on");
411 h3->SetTitle("Distribution of equivalent events");
412 h3->SetDirectory(NULL);
413 h3->SetXTitle(fHist.GetXaxis()->GetTitle());
414 h3->SetYTitle(fHist.GetYaxis()->GetTitle());
415 h3->SetMinimum(0);
416 h3->Draw("colz");
417 h3->SetBit(kCanDelete);
418// catalog->Draw("mirror same *");
419
420 pad->cd();
421 name = Form("%s_2", pad->GetName());
422 p = new TPad(name,name, 0.66, 0.005, 0.995, 0.5, col, 0, 0);
423 p->SetNumber(2);
424 p->Draw();
425 p->cd();
426 h3->Draw("surf3");
427
428 pad->cd();
429 name = Form("%s_3", pad->GetName());
430 p = new TPad(name,name, 0.66, 0.5, 0.995, 0.995, col, 0, 0);
431 p->SetNumber(3);
432 p->Draw();
433 p->cd();
434
435 const Double_t maxr = TMath::Hypot(h3->GetXaxis()->GetXmax(), h3->GetYaxis()->GetXmax());
436 const Int_t nbin = (h3->GetNbinsX()+h3->GetNbinsY())/2;
437 TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr);
438 h->SetDirectory(0);
439 //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr);
440 //h->Sumw2();
441 h->SetXTitle("\\vartheta [\\circ]");
442 h->SetYTitle("<cts>/\\Delta R");
443 h->SetBit(kCanDelete);
444 h->Draw();
445}
446
447// --------------------------------------------------------------------------
448//
449// The following resources are available:
450//
451// MHDisp.Wobble: on/off
452//
453/*
454Int_t MHFalseSource::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
455{
456 Bool_t rc = kFALSE;
457 if (IsEnvDefined(env, prefix, "DistMin", print))
458 {
459 rc = kTRUE;
460 SetMinDist(GetEnvValue(env, prefix, "DistMin", fMinDist));
461 }
462 if (IsEnvDefined(env, prefix, "DistMax", print))
463 {
464 rc = kTRUE;
465 SetMaxDist(GetEnvValue(env, prefix, "DistMax", fMaxDist));
466 }
467
468 if (IsEnvDefined(env, prefix, "DWMin", print))
469 {
470 rc = kTRUE;
471 SetMinDW(GetEnvValue(env, prefix, "DWMin", fMinDist));
472 }
473 if (IsEnvDefined(env, prefix, "DWMax", print))
474 {
475 rc = kTRUE;
476 SetMaxDW(GetEnvValue(env, prefix, "DWMax", fMaxDist));
477 }
478
479 return rc;
480}
481*/
Note: See TracBrowser for help on using the repository browser.