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

Last change on this file since 7193 was 7193, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 16.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, 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#include "MBinning.h"
49
50#include "MParList.h"
51#include "MParameters.h"
52
53#include "MHillasExt.h"
54#include "MHillasSrc.h"
55#include "MSrcPosCam.h"
56#include "MPointingPos.h"
57
58ClassImp(MHDisp);
59
60using namespace std;
61
62// --------------------------------------------------------------------------
63//
64// Default Constructor
65//
66MHDisp::MHDisp(const char *name, const char *title)
67 : fHilExt(0), fDisp(0), fSmearing(-1), fWobble(kFALSE)
68{
69 //
70 // set the name and title of this object
71 //
72 fName = name ? name : "MHDisp";
73 fTitle = title ? title : "3D-plot using Disp vs x, y";
74
75 fHist.SetName("Alpha");
76 fHist.SetTitle("3D-plot of ThetaSq vs x, y");
77 fHist.SetXTitle("x [\\circ]");
78 fHist.SetYTitle("y [\\circ]");
79 fHist.SetZTitle("Eq.cts");
80
81 fHist.SetDirectory(NULL);
82 fHistBg.SetDirectory(NULL);
83 fHistBg1.SetDirectory(NULL);
84 fHistBg2.SetDirectory(NULL);
85
86 fHist.SetBit(TH1::kNoStats);
87}
88
89// --------------------------------------------------------------------------
90//
91// Set binnings (takes BinningFalseSource) and prepare filling of the
92// histogram.
93//
94// Also search for MTime, MObservatory and MPointingPos
95//
96Bool_t MHDisp::SetupFill(const MParList *plist)
97{
98 if (!MHFalseSource::SetupFill(plist))
99 return kFALSE;
100
101 fDisp = (MParameterD*)plist->FindObject("Disp", "MParameterD");
102 if (!fDisp)
103 {
104 *fLog << err << "Disp [MParameterD] not found... abort." << endl;
105 return kFALSE;
106 }
107
108 fHilExt = (MHillasExt*)plist->FindObject("MHillasExt");
109 if (!fHilExt)
110 {
111 *fLog << err << "MHillasExt not found... aborting." << endl;
112 return kFALSE;
113 }
114
115 fHilSrc = (MHillasSrc*)plist->FindObject("MHillasSrc");
116 if (!fHilSrc)
117 {
118 *fLog << err << "MHillasSrc not found... aborting." << endl;
119 return kFALSE;
120 }
121
122 MBinning binsx, binsy;
123 binsx.SetEdges(fHist, 'x');
124 binsy.SetEdges(fHist, 'y');
125 MH::SetBinning(&fHistBg, &binsx, &binsy);
126
127 if (!fHistOff)
128 {
129 MH::SetBinning(&fHistBg1, &binsx, &binsy);
130 MH::SetBinning(&fHistBg2, &binsx, &binsy);
131 }
132
133 return kTRUE;
134}
135
136// --------------------------------------------------------------------------
137//
138// Calculate the delta angle between fSrcPos->GetXY() and v.
139// Return result in deg.
140//
141Double_t MHDisp::DeltaPhiSrc(const TVector2 &v) const
142{
143 return TMath::Abs(fSrcPos->GetXY().DeltaPhi(v))*TMath::RadToDeg();
144}
145
146// --------------------------------------------------------------------------
147//
148// Fill the histogram. For details see the code or the class description
149//
150Bool_t MHDisp::Fill(const MParContainer *par, const Stat_t w)
151{
152 const MHillas *hil = dynamic_cast<const MHillas*>(par);
153 if (!hil)
154 {
155 *fLog << err << "MHDisp::Fill: No container specified!" << endl;
156 return kFALSE;
157 }
158
159 // Get camera rotation angle
160 Double_t rho = 0;
161 if (fTime && fObservatory && fPointPos)
162 rho = fPointPos->RotationAngle(*fObservatory, *fTime);
163
164 // Vector of length disp in direction of shower
165 const TVector2 p(hil->GetCosDelta(), hil->GetSinDelta());
166
167 // Move origin of vector to center-of-gravity of shower and derotate
168 TVector2 pos1 = hil->GetMean()*fMm2Deg + p*fDisp->GetVal();
169
170 Double_t w0 = 1;
171 if (fWobble)
172 {
173 const Double_t delta = DeltaPhiSrc(pos1);
174
175 // Skip off-data not in the same half than the source (here: anti-source)
176 // Could be replaced by some kind of theta cut...
177 if (!fHistOff)
178 {
179 if (delta>165)
180 return kTRUE;
181
182 // Because afterwards the plots are normalized with the number
183 // of entries the non-overlap region corresponds to half the
184 // contents, the overlap region to full contents. By taking
185 // the mean of both distributions we get the same result than
186 // we would have gotten using full off-events with a slightly
187 // increased uncertainty
188 // FIXME: The delta stuff could be replaced by a 2*antitheta cut...
189 w0 = delta>15 ? 1 : 2;
190 }
191
192 // Define by the source position which histogram to fill
193 if (DeltaPhiSrc(fFormerSrc)>90)
194 fHalf = !fHalf;
195 fFormerSrc = fSrcPos->GetXY();
196 }
197
198 // If on-data: Derotate source and P. Translate P to center.
199 TVector2 m;
200 if (fHistOff)
201 {
202 // Derotate all position around camera center by -rho
203 // Move origin of vector to center-of-gravity of shower and derotate
204 pos1=pos1.Rotate(-rho);
205
206 // Shift the source position to 0/0
207 if (fSrcPos)
208 {
209 // m: Position of the camera center in the FS plot
210 m = fSrcPos->GetXY().Rotate(-rho)*fMm2Deg;
211 pos1 -= m;
212 }
213 }
214
215 if (fSmearing<=0)
216 {
217 fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*w0);
218 return kTRUE;
219 }
220
221 // -------------------------------------------------
222 // The following algorithm may look complicated...
223 // It is optimized for speed
224 // -------------------------------------------------
225
226 // Define a vector used to calculate rotated x and y component
227 const TVector2 rot(TMath::Sin(rho), TMath::Cos(rho));
228
229 // Fold event position with the PSF and fill it into histogram
230 const TAxis &axex = *fHist.GetXaxis();
231 const TAxis &axey = *fHist.GetYaxis();
232
233 const Int_t nx = axex.GetNbins();
234 const Int_t ny = axey.GetNbins();
235
236 // normg: Normalization for Gauss [sqrt(2*pi)*sigma]^2
237 const Double_t weight = axex.GetBinWidth(1)*axey.GetBinWidth(1)*w*w0;
238 const Double_t psf = 2*fSmearing*fSmearing;
239 const Double_t pi2 = fSmearing * TMath::Pi()/2;
240 const Double_t normg = pi2*pi2 * TMath::Sqrt(TMath::TwoPi()) / weight;
241 const Int_t bz = fHist.GetZaxis()->FindFixBin(0);
242
243 TH1 &bg = fHalf ? fHistBg1 : fHistBg2;
244
245 // To calculate significance map smear with 2*theta-cut and
246 // do not normalize gauss, for event map use theta-cut/2 instead
247 for (int x=1; x<=nx; x++)
248 {
249 const Double_t cx = axex.GetBinCenter(x);
250 const Double_t px = cx-pos1.X();
251
252 for (int y=1; y<=ny; y++)
253 {
254 const Double_t cy = axey.GetBinCenter(y);
255 const Double_t sp = Sq(px, cy-pos1.Y());
256 const Double_t dp = sp/psf;
257
258 // Values below 1e-3 (>3.5sigma) are not filled into the histogram
259 if (dp<4)
260 {
261 const Double_t rc = TMath::Exp(-dp)/normg;
262 if (!fHistOff)
263 bg.AddBinContent(bg.GetBin(x, y), rc);
264 else
265 fHist.AddBinContent(fHist.GetBin(x, y, bz), rc);
266 }
267
268 if (!fHistOff)
269 continue;
270
271 // If we are filling the signal already (fHistOff!=NULL)
272 // we also fill the background by projecting the
273 // background in the camera into the sky plot.
274 const TVector2 v = TVector2(cx+m.X(), cy+m.Y());
275
276 // Speed up further: Xmax-fWobble
277 if (v.Mod()>axex.GetXmax()) // Gains 10% speed
278 continue;
279
280 const Int_t bx = axex.FindFixBin(v^rot);
281 const Int_t by = axey.FindFixBin(v*rot);
282 const Double_t bg = fHistOff->GetBinContent(bx, by, bz);
283
284 fHistBg.AddBinContent(fHistBg.GetBin(x, y), bg);
285 }
286 }
287 fHist.SetEntries(fHist.GetEntries()+1);
288
289 if (!fHistOff)
290 bg.SetEntries(fHistBg1.GetEntries()+1);
291
292 return kTRUE;
293}
294
295// --------------------------------------------------------------------------
296//
297// Compile the background in camera coordinates from the two background
298// histograms
299//
300Bool_t MHDisp::Finalize()
301{
302 if (fHistOff)
303 return kTRUE;
304
305 const Int_t bz = fHist.GetZaxis()->FindFixBin(0);
306
307 const Double_t n1 = fHistBg1.GetEntries();
308 const Double_t n2 = fHistBg2.GetEntries();
309
310 for (int x=1; x<=fHist.GetNbinsX(); x++)
311 for (int y=1; y<=fHist.GetNbinsY(); y++)
312 {
313 const Int_t bin1 = fHistBg1.GetBin(x, y);
314
315 const Double_t rc =
316 (n1==0?0:fHistBg1.GetBinContent(bin1)/n1)+
317 (n2==0?0:fHistBg2.GetBinContent(bin1)/n2);
318
319 fHist.SetBinContent(x, y, bz, rc/2);
320 }
321
322 fHist.SetEntries(n1+n2);
323
324 // Result corresponds to one smeared background event
325
326 return kTRUE;
327}
328
329
330// --------------------------------------------------------------------------
331//
332// Return the mean signal in h around (0,0) in a distance between
333// 0.33 and 0.47deg
334//
335Double_t MHDisp::GetOffSignal(TH1 &h) const
336{
337 const TAxis &axex = *h.GetXaxis();
338 const TAxis &axey = *h.GetYaxis();
339
340 Int_t cnt = 0;
341 Double_t sum = 0;
342 for (int x=0; x<h.GetNbinsX(); x++)
343 for (int y=0; y<h.GetNbinsY(); y++)
344 {
345 const Double_t d = TMath::Hypot(axex.GetBinCenter(x+1), axey.GetBinCenter(y+1));
346 if (d>0.33 && d<0.47)
347 //if (d>0.4 && d<0.8)
348 {
349 sum += h.GetBinContent(x+1,y+1);
350 cnt++;
351 }
352 }
353
354 return sum/cnt;
355}
356
357// --------------------------------------------------------------------------
358//
359// Fill h2 with the radial profile of h1 around (x0, y0)
360//
361void MHDisp::Profile(TH1 &h2, const TH2 &h1, Axis_t x0, Axis_t y0) const
362{
363 const TAxis &axex = *h1.GetXaxis();
364 const TAxis &axey = *h1.GetYaxis();
365
366 h2.Reset();
367
368 for (int x=1; x<=axex.GetNbins(); x++)
369 for (int y=1; y<=axey.GetNbins(); y++)
370 {
371 const Double_t dx = axex.GetBinCenter(x)-x0;
372 const Double_t dy = axey.GetBinCenter(y)-y0;
373
374 const Double_t r = TMath::Hypot(dx, dy);
375
376 h2.Fill(r, h1.GetBinContent(x, y));
377 }
378}
379
380// --------------------------------------------------------------------------
381//
382// Remove contents of histogram around a circle.
383//
384void MHDisp::MakeDot(TH2 &h2) const
385{
386 const TAxis &axex = *h2.GetXaxis();
387 const TAxis &axey = *h2.GetYaxis();
388
389 const Double_t rmax = fWobble ? axex.GetXmax()-0.4 : axex.GetXmax();
390
391 for (int x=1; x<=axex.GetNbins(); x++)
392 for (int y=1; y<=axey.GetNbins(); y++)
393 {
394 const Int_t bin = h2.GetBin(x,y);
395
396 const Double_t px = h2.GetBinCenter(x);
397 const Double_t py = h2.GetBinCenter(y);
398
399 if (rmax<TMath::Hypot(px, py))
400 h2.SetBinContent(bin, 0);
401 }
402}
403
404// --------------------------------------------------------------------------
405//
406// Calculate from signal and background the significance map
407//
408void MHDisp::MakeSignificance(TH2 &s, const TH2 &h1, const TH2 &h2, const Double_t scale) const
409{
410 const TAxis &axex = *s.GetXaxis();
411 const TAxis &axey = *s.GetYaxis();
412
413 for (int x=1; x<=axex.GetNbins(); x++)
414 for (int y=1; y<=axey.GetNbins(); y++)
415 {
416 const Int_t bin = s.GetBin(x,y);
417
418 const Double_t sig = h1.GetBinContent(bin);
419 const Double_t bg = h2.GetBinContent(bin);
420
421 const Double_t S = MMath::SignificanceLiMaSigned(sig, bg, TMath::Abs(scale));
422
423 s.SetBinContent(bin, S);
424 }
425}
426
427
428void MHDisp::Paint(Option_t *o)
429{
430 // Compile Background if necessary
431 Finalize();
432
433 // Paint result
434 TVirtualPad *pad = gPad;
435
436 pad->cd(1);
437
438 // Project on data onto yx-plane
439 fHist.GetZaxis()->SetRange(0,9999);
440 TH2 *h1=(TH2*)fHist.Project3D("yx_on");
441
442 // Set Glow-palette (PaintSurface doesn't allow more than 99 colors)
443 MH::SetPalette(fHistOff?"glowsym":"glow1", 99);
444 h1->SetContour(99);
445
446 Double_t scale = 1;
447 if (fHistOff)
448 {
449 // Project off data onto yx-plane and subtract it from on-data
450 fHistOff->GetZaxis()->SetRange(0,9999);
451 TH2 *h=(TH2*)fHistOff->Project3D("yx_off");
452
453 scale = -1;
454
455 const Double_t h1off = GetOffSignal(*h1);
456 const Double_t hoff = GetOffSignal(fHistBg);
457
458 scale = hoff==0 ? -1 : -h1off/hoff;
459
460 const Bool_t sig = kFALSE;
461
462 if (!sig)
463 h1->Add(&fHistBg, scale);
464 else
465 MakeSignificance(*h1, *h1, fHistBg, scale);
466
467 MakeDot(*h1);
468
469 delete h;
470
471 MakeSymmetric(h1);
472 }
473
474 pad->cd(3);
475 TH1 *h2 = (TH1*)gPad->FindObject("RadProf");
476
477 TString opt(o);
478 opt.ToLower();
479
480 if (h1 && h2)
481 {
482 Int_t ix, iy, iz;
483 h1->GetMaximumBin(ix, iy, iz);
484
485 const Double_t x0 = h1->GetXaxis()->GetBinCenter(ix);
486 const Double_t y0 = h1->GetYaxis()->GetBinCenter(iy);
487 //const Double_t w0 = h1->GetXaxis()->GetBinWidth(1);
488
489 Profile(*h2, *h1, 0, 0);
490 //Profile(*h2, *h1, x0, y0);
491
492 // Replace with MAlphaFitter?
493 TF1 func("fcn", "gaus + [3]*x*x + [4]");
494 func.SetLineWidth(1);
495 func.SetLineColor(kBlue);
496
497 func.SetParLimits(2, h2->GetBinWidth(1), 1.0);
498
499 func.SetParameter(0, h2->GetBinContent(1));
500 func.FixParameter(1, 0);
501 func.SetParameter(2, 0.15);
502 if (fHistOff)
503 func.FixParameter(3, 0);
504 func.SetParameter(4, h2->GetBinContent(15));
505 h2->Fit(&func, "IQ", "", 0, 1.0);
506
507 h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ f=%.2f", x0, y0, func.GetParameter(2), TMath::Abs(scale)));
508 }
509}
510
511void MHDisp::Draw(Option_t *o)
512{
513 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
514 const Int_t col = pad->GetFillColor();
515 pad->SetBorderMode(0);
516
517 AppendPad(o);
518
519 TString name = Form("%s_1", pad->GetName());
520 TPad *p = new TPad(name,name, 0.005, 0.005, 0.65, 0.995, col, 0, 0);
521 p->SetNumber(1);
522 p->Draw();
523 p->cd();
524
525 TH1 *h3 = fHist.Project3D("yx_on");
526 h3->SetTitle("Distribution of equivalent events");
527 h3->SetDirectory(NULL);
528 h3->SetXTitle(fHist.GetXaxis()->GetTitle());
529 h3->SetYTitle(fHist.GetYaxis()->GetTitle());
530 h3->SetMinimum(0);
531 h3->Draw("colz");
532 h3->SetBit(kCanDelete);
533
534 if (fHistOff)
535 GetCatalog()->Draw("white mirror same *");
536
537 pad->cd();
538 name = Form("%s_2", pad->GetName());
539 p = new TPad(name,name, 0.66, 0.005, 0.995, 0.5, col, 0, 0);
540 p->SetNumber(2);
541 p->Draw();
542 p->cd();
543 h3->Draw("surf3");
544
545 pad->cd();
546 name = Form("%s_3", pad->GetName());
547 p = new TPad(name,name, 0.66, 0.5, 0.995, 0.995, col, 0, 0);
548 p->SetNumber(3);
549 p->SetGrid();
550 p->Draw();
551 p->cd();
552
553 const Double_t maxr = TMath::Hypot(h3->GetXaxis()->GetXmax(), h3->GetYaxis()->GetXmax());
554 const Int_t nbin = (h3->GetNbinsX()+h3->GetNbinsY())/2;
555 TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr);
556 h->SetDirectory(0);
557 //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr);
558 //h->Sumw2();
559 h->SetXTitle("\\vartheta [\\circ]");
560 h->SetYTitle("<cts>/\\Delta R");
561 h->SetBit(kCanDelete);
562 h->Draw();
563}
564
565// --------------------------------------------------------------------------
566//
567// The following resources are available:
568//
569// MHDisp.Smearing: 0.1
570// MHDisp.Wobble: on/off
571//
572Int_t MHDisp::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
573{
574 Int_t rc = MHFalseSource::ReadEnv(env, prefix, print);
575 if (rc==kERROR)
576 return kERROR;
577
578 if (IsEnvDefined(env, prefix, "Smearing", print))
579 {
580 rc = kTRUE;
581 SetSmearing(GetEnvValue(env, prefix, "Smearing", fSmearing));
582 }
583 if (IsEnvDefined(env, prefix, "Wobble", print))
584 {
585 rc = kTRUE;
586 SetWobble(GetEnvValue(env, prefix, "Wobble", fWobble));
587 }
588
589 return rc;
590}
591
Note: See TracBrowser for help on using the repository browser.