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

Last change on this file since 7202 was 7202, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 17.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#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 // -------------------------------------------------
216 // The following algorithm may look complicated...
217 // It is optimized for speed
218 // -------------------------------------------------
219
220 // Define a vector used to calculate rotated x and y component
221 const TVector2 rot(TMath::Sin(rho), TMath::Cos(rho));
222
223 // Fold event position with the PSF and fill it into histogram
224 const TAxis &axex = *fHist.GetXaxis();
225 const TAxis &axey = *fHist.GetYaxis();
226
227 const Int_t nx = axex.GetNbins();
228 const Int_t ny = axey.GetNbins();
229
230 // normg: Normalization for Gauss [sqrt(2*pi)*sigma]^2
231 const Double_t weight = axex.GetBinWidth(1)*axey.GetBinWidth(1)*w*w0;
232 const Double_t psf = 2*fSmearing*fSmearing;
233 const Double_t pi2 = fSmearing * TMath::Pi()/2;
234 const Double_t normg = pi2*pi2 * TMath::Sqrt(TMath::TwoPi()) / weight;
235 const Int_t bz = fHist.GetZaxis()->FindFixBin(0);
236
237 TH2 &bg = fHalf ? fHistBg1 : fHistBg2;
238
239 const Bool_t smear = fSmearing>0;
240 if (!smear)
241 {
242 if (!fHistOff)
243 bg.Fill(pos1.X(), pos1.Y(), w*w0);
244 else
245 fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*w0);
246 }
247
248 // To calculate significance map smear with 2*theta-cut and
249 // do not normalize gauss, for event map use theta-cut/2 instead
250 for (int x=1; x<=nx; x++)
251 {
252 const Double_t cx = axex.GetBinCenter(x);
253 const Double_t px = cx-pos1.X();
254
255 for (int y=1; y<=ny; y++)
256 {
257 const Double_t cy = axey.GetBinCenter(y);
258 const Double_t sp = Sq(px, cy-pos1.Y());
259 if (smear)
260 {
261 const Double_t dp = sp/psf;
262
263 // Values below 1e-3 (>3.5sigma) are not filled into the histogram
264 if (dp<4)
265 {
266 const Double_t rc = TMath::Exp(-dp)/normg;
267 if (!fHistOff)
268 bg.AddBinContent(bg.GetBin(x, y), rc);
269 else
270 fHist.AddBinContent(fHist.GetBin(x, y, bz), rc);
271 }
272 }
273
274 if (!fHistOff)
275 continue;
276
277 // If we are filling the signal already (fHistOff!=NULL)
278 // we also fill the background by projecting the
279 // background in the camera into the sky plot.
280 const TVector2 v = TVector2(cx+m.X(), cy+m.Y());
281
282 // Speed up further: Xmax-fWobble
283 if (v.Mod()>axex.GetXmax()) // Gains 10% speed
284 continue;
285
286 const Int_t bx = axex.FindFixBin(v^rot);
287 const Int_t by = axey.FindFixBin(v*rot);
288 const Double_t bg = fHistOff->GetBinContent(bx, by, bz);
289
290 fHistBg.AddBinContent(fHistBg.GetBin(x, y), bg);
291 }
292 }
293
294 if (fHistOff)
295 fHistBg.SetEntries(fHistBg.GetEntries()+1);
296
297 if (!smear)
298 return kTRUE;
299
300 if (!fHistOff)
301 bg.SetEntries(bg.GetEntries()+1);
302 else
303 fHist.SetEntries(fHist.GetEntries()+1);
304
305 return kTRUE;
306}
307
308// --------------------------------------------------------------------------
309//
310// Compile the background in camera coordinates from the two background
311// histograms
312//
313Bool_t MHDisp::Finalize()
314{
315 if (fHistOff)
316 return kTRUE;
317
318 const Int_t bz = fHist.GetZaxis()->FindFixBin(0);
319
320 const Double_t n1 = fHistBg1.GetEntries();
321 const Double_t n2 = fHistBg2.GetEntries();
322
323 for (int x=1; x<=fHist.GetNbinsX(); x++)
324 for (int y=1; y<=fHist.GetNbinsY(); y++)
325 {
326 const Int_t bin1 = fHistBg1.GetBin(x, y);
327
328 const Double_t rc =
329 (n1==0?0:fHistBg1.GetBinContent(bin1)/n1)+
330 (n2==0?0:fHistBg2.GetBinContent(bin1)/n2);
331
332 fHist.SetBinContent(x, y, bz, rc/2);
333 }
334
335 fHist.SetEntries(n1+n2);
336
337 // Result corresponds to one smeared background event
338
339 return kTRUE;
340}
341
342
343// --------------------------------------------------------------------------
344//
345// Return the mean signal in h around (0,0) in a distance between
346// 0.33 and 0.47deg
347//
348Double_t MHDisp::GetOffSignal(TH1 &h) const
349{
350 const TAxis &axex = *h.GetXaxis();
351 const TAxis &axey = *h.GetYaxis();
352
353 Int_t cnt = 0;
354 Double_t sum = 0;
355 for (int x=0; x<h.GetNbinsX(); x++)
356 for (int y=0; y<h.GetNbinsY(); y++)
357 {
358 const Double_t d = TMath::Hypot(axex.GetBinCenter(x+1), axey.GetBinCenter(y+1));
359 if (d>0.33 && d<0.47)
360 //if (d>0.4 && d<0.8)
361 {
362 sum += h.GetBinContent(x+1,y+1);
363 cnt++;
364 }
365 }
366
367 return sum/cnt;
368}
369
370// --------------------------------------------------------------------------
371//
372// Fill h2 with the radial profile of h1 around (x0, y0)
373//
374void MHDisp::Profile(TH1 &h2, const TH2 &h1, Axis_t x0, Axis_t y0) const
375{
376 const TAxis &axex = *h1.GetXaxis();
377 const TAxis &axey = *h1.GetYaxis();
378
379 h2.Reset();
380
381 for (int x=1; x<=axex.GetNbins(); x++)
382 for (int y=1; y<=axey.GetNbins(); y++)
383 {
384 const Double_t dx = axex.GetBinCenter(x)-x0;
385 const Double_t dy = axey.GetBinCenter(y)-y0;
386
387 const Double_t r = TMath::Hypot(dx, dy);
388
389 h2.Fill(r, h1.GetBinContent(x, y));
390 }
391}
392
393// --------------------------------------------------------------------------
394//
395// Remove contents of histogram around a circle.
396//
397void MHDisp::MakeDot(TH2 &h2) const
398{
399 const TAxis &axex = *h2.GetXaxis();
400 const TAxis &axey = *h2.GetYaxis();
401
402 const Double_t rmax = fWobble ? axex.GetXmax()-0.4 : axex.GetXmax();
403
404 for (int x=1; x<=axex.GetNbins(); x++)
405 for (int y=1; y<=axey.GetNbins(); y++)
406 {
407 const Int_t bin = h2.GetBin(x,y);
408
409 const Double_t px = h2.GetBinCenter(x);
410 const Double_t py = h2.GetBinCenter(y);
411
412 if (rmax<TMath::Hypot(px, py))
413 h2.SetBinContent(bin, 0);
414 }
415}
416
417// --------------------------------------------------------------------------
418//
419// Calculate from signal and background the significance map
420//
421void MHDisp::MakeSignificance(TH2 &s, const TH2 &h1, const TH2 &h2, const Double_t scale) const
422{
423 const TAxis &axex = *s.GetXaxis();
424 const TAxis &axey = *s.GetYaxis();
425
426 for (int x=1; x<=axex.GetNbins(); x++)
427 for (int y=1; y<=axey.GetNbins(); y++)
428 {
429 const Int_t bin = s.GetBin(x,y);
430
431 const Double_t sig = h1.GetBinContent(bin);
432 const Double_t bg = h2.GetBinContent(bin);
433
434 const Double_t S = MMath::SignificanceLiMaSigned(sig, bg, TMath::Abs(scale));
435
436 s.SetBinContent(bin, S);
437 }
438}
439
440
441void MHDisp::Paint(Option_t *o)
442{
443 // Compile Background if necessary
444 Finalize();
445
446 // Paint result
447 TVirtualPad *pad = gPad;
448
449 pad->cd(1);
450
451 // Project on data onto yx-plane
452 fHist.GetZaxis()->SetRange(0,9999);
453 TH2 *h1=(TH2*)fHist.Project3D("yx_on");
454
455 // Set Glow-palette (PaintSurface doesn't allow more than 99 colors)
456 MH::SetPalette(fHistOff?"glowsym":"glow1", 99);
457 h1->SetContour(99);
458
459 Double_t scale = 1;
460 if (fHistOff)
461 {
462 // Project off data onto yx-plane and subtract it from on-data
463 fHistOff->GetZaxis()->SetRange(0,9999);
464 TH2 *h=(TH2*)fHistOff->Project3D("yx_off");
465
466 scale = -1;
467
468 const Double_t h1off = GetOffSignal(*h1);
469 const Double_t hoff = GetOffSignal(fHistBg);
470
471 scale = hoff==0 ? -1 : -h1off/hoff;
472
473 const Bool_t sig = kFALSE;
474
475 if (!sig)
476 h1->Add(&fHistBg, scale);
477 else
478 MakeSignificance(*h1, *h1, fHistBg, scale);
479
480 MakeDot(*h1);
481
482 delete h;
483
484 MakeSymmetric(h1);
485 }
486
487 pad->cd(3);
488 TH1 *h2 = (TH1*)gPad->FindObject("RadProf");
489
490 TString opt(o);
491 opt.ToLower();
492
493 if (h1 && h2)
494 {
495 Int_t ix, iy, iz;
496 h1->GetMaximumBin(ix, iy, iz);
497
498 const Double_t x0 = h1->GetXaxis()->GetBinCenter(ix);
499 const Double_t y0 = h1->GetYaxis()->GetBinCenter(iy);
500 //const Double_t w0 = h1->GetXaxis()->GetBinWidth(1);
501
502 Profile(*h2, *h1, 0, 0);
503 //Profile(*h2, *h1, x0, y0);
504
505 // Replace with MAlphaFitter?
506 TF1 func("fcn", "gaus + [3]*x*x + [4]");
507 func.SetLineWidth(1);
508 func.SetLineColor(kBlue);
509
510 func.SetParLimits(2, h2->GetBinWidth(1), 1.0);
511
512 func.SetParameter(0, h2->GetBinContent(1));
513 func.FixParameter(1, 0);
514 func.SetParameter(2, 0.15);
515 if (fHistOff)
516 func.FixParameter(3, 0);
517 func.SetParameter(4, h2->GetBinContent(15));
518 h2->Fit(&func, "IQ", "", 0, 1.0);
519
520 h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ f=%.2f", x0, y0, func.GetParameter(2), TMath::Abs(scale)));
521 }
522}
523
524void MHDisp::Draw(Option_t *o)
525{
526 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
527 const Int_t col = pad->GetFillColor();
528 pad->SetBorderMode(0);
529
530 AppendPad(o);
531
532 TString name = Form("%s_1", pad->GetName());
533 TPad *p = new TPad(name,name, 0.005, 0.005, 0.65, 0.995, col, 0, 0);
534 p->SetNumber(1);
535 p->Draw();
536 p->cd();
537
538 TH1 *h3 = fHist.Project3D("yx_on");
539 h3->SetTitle("Distribution of equivalent events");
540 h3->SetDirectory(NULL);
541 h3->SetXTitle(fHist.GetXaxis()->GetTitle());
542 h3->SetYTitle(fHist.GetYaxis()->GetTitle());
543 h3->SetMinimum(0);
544 h3->Draw("colz");
545 h3->SetBit(kCanDelete);
546
547 if (fHistOff)
548 GetCatalog()->Draw("white mirror same *");
549
550 pad->cd();
551 name = Form("%s_2", pad->GetName());
552 p = new TPad(name,name, 0.66, 0.005, 0.995, 0.5, col, 0, 0);
553 p->SetNumber(2);
554 p->Draw();
555 p->cd();
556 h3->Draw("surf3");
557
558 pad->cd();
559 name = Form("%s_3", pad->GetName());
560 p = new TPad(name,name, 0.66, 0.5, 0.995, 0.995, col, 0, 0);
561 p->SetNumber(3);
562 p->SetGrid();
563 p->Draw();
564 p->cd();
565
566 const Double_t maxr = TMath::Hypot(h3->GetXaxis()->GetXmax(), h3->GetYaxis()->GetXmax());
567 const Int_t nbin = (h3->GetNbinsX()+h3->GetNbinsY())/2;
568 TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr);
569 h->SetDirectory(0);
570 //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr);
571 //h->Sumw2();
572 h->SetXTitle("\\vartheta [\\circ]");
573 h->SetYTitle("<cts>/\\Delta R");
574 h->SetBit(kCanDelete);
575 h->Draw();
576}
577
578// --------------------------------------------------------------------------
579//
580// The following resources are available:
581//
582// MHDisp.Smearing: 0.1
583// MHDisp.Wobble: on/off
584//
585Int_t MHDisp::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
586{
587 Int_t rc = MHFalseSource::ReadEnv(env, prefix, print);
588 if (rc==kERROR)
589 return kERROR;
590
591 if (IsEnvDefined(env, prefix, "Smearing", print))
592 {
593 rc = kTRUE;
594 SetSmearing(GetEnvValue(env, prefix, "Smearing", fSmearing));
595 }
596 if (IsEnvDefined(env, prefix, "Wobble", print))
597 {
598 rc = kTRUE;
599 SetWobble(GetEnvValue(env, prefix, "Wobble", fWobble));
600 }
601
602 return rc;
603}
604
Note: See TracBrowser for help on using the repository browser.