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 "MHillas.h"
|
---|
54 | #include "MSrcPosCam.h"
|
---|
55 | #include "MPointingPos.h"
|
---|
56 | #include "MPointingDev.h"
|
---|
57 |
|
---|
58 | ClassImp(MHDisp);
|
---|
59 |
|
---|
60 | using namespace std;
|
---|
61 |
|
---|
62 | // --------------------------------------------------------------------------
|
---|
63 | //
|
---|
64 | // Default Constructor
|
---|
65 | //
|
---|
66 | MHDisp::MHDisp(const char *name, const char *title)
|
---|
67 | : fDisp(0), fDeviation(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 | //
|
---|
96 | Bool_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 | fDeviation = (MPointingDev*)plist->FindObject("MPointingDev");
|
---|
109 | if (!fDeviation)
|
---|
110 | *fLog << warn << "MPointingDev not found... ignored." << endl;
|
---|
111 |
|
---|
112 | MBinning binsx, binsy;
|
---|
113 | binsx.SetEdges(fHist, 'x');
|
---|
114 | binsy.SetEdges(fHist, 'y');
|
---|
115 | MH::SetBinning(&fHistBg, &binsx, &binsy);
|
---|
116 |
|
---|
117 | if (!fHistOff)
|
---|
118 | {
|
---|
119 | MH::SetBinning(&fHistBg1, &binsx, &binsy);
|
---|
120 | MH::SetBinning(&fHistBg2, &binsx, &binsy);
|
---|
121 | }
|
---|
122 |
|
---|
123 | return kTRUE;
|
---|
124 | }
|
---|
125 |
|
---|
126 | // --------------------------------------------------------------------------
|
---|
127 | //
|
---|
128 | // Calculate the delta angle between fSrcPos->GetXY() and v.
|
---|
129 | // Return result in deg.
|
---|
130 | //
|
---|
131 | Double_t MHDisp::DeltaPhiSrc(const TVector2 &v) const
|
---|
132 | {
|
---|
133 | return TMath::Abs(fSrcPos->GetXY().DeltaPhi(v))*TMath::RadToDeg();
|
---|
134 | }
|
---|
135 |
|
---|
136 | // --------------------------------------------------------------------------
|
---|
137 | //
|
---|
138 | // Fill the histogram. For details see the code or the class description
|
---|
139 | //
|
---|
140 | Bool_t MHDisp::Fill(const MParContainer *par, const Stat_t w)
|
---|
141 | {
|
---|
142 | const MHillas *hil = dynamic_cast<const MHillas*>(par);
|
---|
143 | if (!hil)
|
---|
144 | {
|
---|
145 | *fLog << err << "MHDisp::Fill: No container specified!" << endl;
|
---|
146 | return kFALSE;
|
---|
147 | }
|
---|
148 |
|
---|
149 | // Get camera rotation angle
|
---|
150 | Double_t rho = 0;
|
---|
151 | if (fTime && fObservatory && fPointPos)
|
---|
152 | rho = fPointPos->RotationAngle(*fObservatory, *fTime, fDeviation);
|
---|
153 |
|
---|
154 | // Vector of length disp in direction of shower
|
---|
155 | // Move origin of vector to center-of-gravity of shower and derotate
|
---|
156 | TVector2 pos1 = hil->GetMean()*fMm2Deg + hil->GetNormAxis()*fDisp->GetVal();
|
---|
157 |
|
---|
158 | Double_t w0 = 1;
|
---|
159 | if (fWobble)
|
---|
160 | {
|
---|
161 | const Double_t delta = DeltaPhiSrc(pos1);
|
---|
162 |
|
---|
163 | // Skip off-data not in the same half than the source (here: anti-source)
|
---|
164 | // Could be replaced by some kind of theta cut...
|
---|
165 | if (!fHistOff)
|
---|
166 | {
|
---|
167 | if (delta>165)
|
---|
168 | return kTRUE;
|
---|
169 |
|
---|
170 | // Because afterwards the plots are normalized with the number
|
---|
171 | // of entries the non-overlap region corresponds to half the
|
---|
172 | // contents, the overlap region to full contents. By taking
|
---|
173 | // the mean of both distributions we get the same result than
|
---|
174 | // we would have gotten using full off-events with a slightly
|
---|
175 | // increased uncertainty
|
---|
176 | // FIXME: The delta stuff could be replaced by a 2*antitheta cut...
|
---|
177 | w0 = delta>15 ? 1 : 2;
|
---|
178 | }
|
---|
179 |
|
---|
180 | // Define by the source position which histogram to fill
|
---|
181 | if (DeltaPhiSrc(fFormerSrc)>90)
|
---|
182 | fHalf = !fHalf;
|
---|
183 | fFormerSrc = fSrcPos->GetXY();
|
---|
184 | }
|
---|
185 |
|
---|
186 | // If on-data: Derotate source and P. Translate P to center.
|
---|
187 | TVector2 m;
|
---|
188 | if (fHistOff)
|
---|
189 | {
|
---|
190 | // Derotate all position around camera center by -rho
|
---|
191 | // Move origin of vector to center-of-gravity of shower and derotate
|
---|
192 | pos1=pos1.Rotate(-rho);
|
---|
193 |
|
---|
194 | // Shift the source position to 0/0
|
---|
195 | if (fSrcPos)
|
---|
196 | {
|
---|
197 | // m: Position of the camera center in the FS plot
|
---|
198 | m = fSrcPos->GetXY().Rotate(-rho)*fMm2Deg;
|
---|
199 | pos1 -= m;
|
---|
200 | }
|
---|
201 | }
|
---|
202 |
|
---|
203 | // -------------------------------------------------
|
---|
204 | // The following algorithm may look complicated...
|
---|
205 | // It is optimized for speed
|
---|
206 | // -------------------------------------------------
|
---|
207 |
|
---|
208 | // Define a vector used to calculate rotated x and y component
|
---|
209 | const TVector2 rot(TMath::Sin(rho), TMath::Cos(rho));
|
---|
210 |
|
---|
211 | // Fold event position with the PSF and fill it into histogram
|
---|
212 | const TAxis &axex = *fHist.GetXaxis();
|
---|
213 | const TAxis &axey = *fHist.GetYaxis();
|
---|
214 |
|
---|
215 | const Int_t nx = axex.GetNbins();
|
---|
216 | const Int_t ny = axey.GetNbins();
|
---|
217 |
|
---|
218 | // normg: Normalization for Gauss [sqrt(2*pi)*sigma]^2
|
---|
219 | const Double_t weight = axex.GetBinWidth(1)*axey.GetBinWidth(1)*w*w0;
|
---|
220 | const Double_t psf = 2*fSmearing*fSmearing;
|
---|
221 | const Double_t pi2 = fSmearing * TMath::Pi()/2;
|
---|
222 | const Double_t normg = pi2*pi2 * TMath::Sqrt(TMath::TwoPi()) / weight;
|
---|
223 | const Int_t bz = fHist.GetZaxis()->FindFixBin(0);
|
---|
224 |
|
---|
225 | TH2 &bg = fHalf ? fHistBg1 : fHistBg2;
|
---|
226 |
|
---|
227 | const Bool_t smear = fSmearing>0;
|
---|
228 | if (!smear)
|
---|
229 | {
|
---|
230 | if (!fHistOff)
|
---|
231 | bg.Fill(pos1.X(), pos1.Y(), w*w0);
|
---|
232 | else
|
---|
233 | fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*w0);
|
---|
234 | }
|
---|
235 |
|
---|
236 | // To calculate significance map smear with 2*theta-cut and
|
---|
237 | // do not normalize gauss, for event map use theta-cut/2 instead
|
---|
238 | if (smear || fHistOff)
|
---|
239 | {
|
---|
240 | for (int x=1; x<=nx; x++)
|
---|
241 | {
|
---|
242 | const Double_t cx = axex.GetBinCenter(x);
|
---|
243 | const Double_t px = cx-pos1.X();
|
---|
244 |
|
---|
245 | for (int y=1; y<=ny; y++)
|
---|
246 | {
|
---|
247 | const Double_t cy = axey.GetBinCenter(y);
|
---|
248 | const Double_t sp = Sq(px, cy-pos1.Y());
|
---|
249 | if (smear)
|
---|
250 | {
|
---|
251 | const Double_t dp = sp/psf;
|
---|
252 |
|
---|
253 | // Values below 1e-3 (>3.5sigma) are not filled into the histogram
|
---|
254 | if (dp<4)
|
---|
255 | {
|
---|
256 | const Double_t rc = TMath::Exp(-dp)/normg;
|
---|
257 | if (!fHistOff)
|
---|
258 | bg.AddBinContent(bg.GetBin(x, y), rc);
|
---|
259 | else
|
---|
260 | fHist.AddBinContent(fHist.GetBin(x, y, bz), rc);
|
---|
261 | }
|
---|
262 | }
|
---|
263 |
|
---|
264 | if (!fHistOff)
|
---|
265 | continue;
|
---|
266 |
|
---|
267 | // If we are filling the signal already (fHistOff!=NULL)
|
---|
268 | // we also fill the background by projecting the
|
---|
269 | // background in the camera into the sky plot.
|
---|
270 | const TVector2 v = TVector2(cx+m.X(), cy+m.Y());
|
---|
271 |
|
---|
272 | // Speed up further: Xmax-fWobble
|
---|
273 | if (v.Mod()>axex.GetXmax()) // Gains 10% speed
|
---|
274 | continue;
|
---|
275 |
|
---|
276 | const Int_t bx = axex.FindFixBin(v^rot);
|
---|
277 | const Int_t by = axey.FindFixBin(v*rot);
|
---|
278 | const Double_t bg = fHistOff->GetBinContent(bx, by, bz);
|
---|
279 |
|
---|
280 | fHistBg.AddBinContent(fHistBg.GetBin(x, y), bg);
|
---|
281 | }
|
---|
282 | }
|
---|
283 | }
|
---|
284 |
|
---|
285 | if (fHistOff)
|
---|
286 | fHistBg.SetEntries(fHistBg.GetEntries()+1);
|
---|
287 |
|
---|
288 | if (!smear)
|
---|
289 | return kTRUE;
|
---|
290 |
|
---|
291 | if (!fHistOff)
|
---|
292 | bg.SetEntries(bg.GetEntries()+1);
|
---|
293 | else
|
---|
294 | fHist.SetEntries(fHist.GetEntries()+1);
|
---|
295 |
|
---|
296 | return kTRUE;
|
---|
297 | }
|
---|
298 |
|
---|
299 | // --------------------------------------------------------------------------
|
---|
300 | //
|
---|
301 | // Compile the background in camera coordinates from the two background
|
---|
302 | // histograms
|
---|
303 | //
|
---|
304 | Bool_t MHDisp::Finalize()
|
---|
305 | {
|
---|
306 | if (fHistOff)
|
---|
307 | return kTRUE;
|
---|
308 |
|
---|
309 | const Int_t bz = fHist.GetZaxis()->FindFixBin(0);
|
---|
310 |
|
---|
311 | const Double_t n1 = fHistBg1.GetEntries();
|
---|
312 | const Double_t n2 = fHistBg2.GetEntries();
|
---|
313 |
|
---|
314 | for (int x=1; x<=fHist.GetNbinsX(); x++)
|
---|
315 | for (int y=1; y<=fHist.GetNbinsY(); y++)
|
---|
316 | {
|
---|
317 | const Int_t bin1 = fHistBg1.GetBin(x, y);
|
---|
318 |
|
---|
319 | const Double_t rc =
|
---|
320 | (n1==0?0:fHistBg1.GetBinContent(bin1)/n1)+
|
---|
321 | (n2==0?0:fHistBg2.GetBinContent(bin1)/n2);
|
---|
322 |
|
---|
323 | fHist.SetBinContent(x, y, bz, rc/2);
|
---|
324 | }
|
---|
325 |
|
---|
326 | fHist.SetEntries(n1+n2);
|
---|
327 |
|
---|
328 | // Result corresponds to one smeared background event
|
---|
329 |
|
---|
330 | return kTRUE;
|
---|
331 | }
|
---|
332 |
|
---|
333 |
|
---|
334 | // --------------------------------------------------------------------------
|
---|
335 | //
|
---|
336 | // Return the mean signal in h around (0,0) in a distance between
|
---|
337 | // 0.325 and 0.475deg
|
---|
338 | //
|
---|
339 | Double_t MHDisp::GetOffSignal(TH1 &h) const
|
---|
340 | {
|
---|
341 | const TAxis &axex = *h.GetXaxis();
|
---|
342 | const TAxis &axey = *h.GetYaxis();
|
---|
343 |
|
---|
344 | Int_t cnt = 0;
|
---|
345 | Double_t sum = 0;
|
---|
346 | for (int x=0; x<h.GetNbinsX(); x++)
|
---|
347 | for (int y=0; y<h.GetNbinsY(); y++)
|
---|
348 | {
|
---|
349 | const Double_t d = TMath::Hypot(axex.GetBinCenter(x+1), axey.GetBinCenter(y+1));
|
---|
350 | if (d>0.325 && d<0.475)
|
---|
351 | {
|
---|
352 | sum += h.GetBinContent(x+1,y+1);
|
---|
353 | cnt++;
|
---|
354 | }
|
---|
355 | }
|
---|
356 | return sum/cnt;
|
---|
357 | }
|
---|
358 |
|
---|
359 | // --------------------------------------------------------------------------
|
---|
360 | //
|
---|
361 | // Fill h2 with the radial profile of h1 around (x0, y0)
|
---|
362 | //
|
---|
363 | void MHDisp::Profile(TH1 &h2, const TH2 &h1, Axis_t x0, Axis_t y0) const
|
---|
364 | {
|
---|
365 | const TAxis &axex = *h1.GetXaxis();
|
---|
366 | const TAxis &axey = *h1.GetYaxis();
|
---|
367 |
|
---|
368 | h2.Reset();
|
---|
369 |
|
---|
370 | for (int x=1; x<=axex.GetNbins(); x++)
|
---|
371 | for (int y=1; y<=axey.GetNbins(); y++)
|
---|
372 | {
|
---|
373 | const Double_t dx = axex.GetBinCenter(x)-x0;
|
---|
374 | const Double_t dy = axey.GetBinCenter(y)-y0;
|
---|
375 |
|
---|
376 | const Double_t r = TMath::Hypot(dx, dy);
|
---|
377 |
|
---|
378 | h2.Fill(r, h1.GetBinContent(x, y));
|
---|
379 | }
|
---|
380 | }
|
---|
381 |
|
---|
382 | // --------------------------------------------------------------------------
|
---|
383 | //
|
---|
384 | // Remove contents of histogram around a circle.
|
---|
385 | //
|
---|
386 | void MHDisp::MakeDot(TH2 &h2) const
|
---|
387 | {
|
---|
388 | const TAxis &axex = *h2.GetXaxis();
|
---|
389 | const TAxis &axey = *h2.GetYaxis();
|
---|
390 |
|
---|
391 | const Double_t rmax = (fWobble ? axex.GetXmax()-0.4 : axex.GetXmax()) - axex.GetBinWidth(1);
|
---|
392 |
|
---|
393 | for (int x=1; x<=axex.GetNbins(); x++)
|
---|
394 | for (int y=1; y<=axey.GetNbins(); y++)
|
---|
395 | {
|
---|
396 | const Int_t bin = h2.GetBin(x,y);
|
---|
397 |
|
---|
398 | const Double_t px = h2.GetBinCenter(x);
|
---|
399 | const Double_t py = h2.GetBinCenter(y);
|
---|
400 |
|
---|
401 | if (rmax<TMath::Hypot(px, py))
|
---|
402 | h2.SetBinContent(bin, 0);
|
---|
403 | //else
|
---|
404 | // if (h2.GetBinContent(bin)==0)
|
---|
405 | // h2.SetBinContent(bin, 1e-10);
|
---|
406 | }
|
---|
407 | }
|
---|
408 |
|
---|
409 | // --------------------------------------------------------------------------
|
---|
410 | //
|
---|
411 | // Calculate from signal and background the significance map
|
---|
412 | //
|
---|
413 | void MHDisp::MakeSignificance(TH2 &s, const TH2 &h1, const TH2 &h2, const Double_t scale) const
|
---|
414 | {
|
---|
415 | const TAxis &axex = *s.GetXaxis();
|
---|
416 | const TAxis &axey = *s.GetYaxis();
|
---|
417 |
|
---|
418 | const Int_t nx = axex.GetNbins();
|
---|
419 | const Int_t ny = axey.GetNbins();
|
---|
420 |
|
---|
421 | const Int_t n = TMath::Nint(0.2/axex.GetBinWidth(1));
|
---|
422 |
|
---|
423 | const Double_t sc = h2.GetEntries()/fHistOff->GetEntries();
|
---|
424 |
|
---|
425 | for (int x=1; x<=nx; x++)
|
---|
426 | for (int y=1; y<=ny; y++)
|
---|
427 | {
|
---|
428 | Double_t S=0;
|
---|
429 |
|
---|
430 | // Only calculate significances for pixels far enough
|
---|
431 | // from the border to get all expected pixels.
|
---|
432 | if (TMath::Hypot((float)x-0.5*nx, (float)y-0.5*ny)<0.5*axex.GetNbins()-n)
|
---|
433 | {
|
---|
434 | Double_t sig=0;
|
---|
435 | Double_t bg =0;
|
---|
436 |
|
---|
437 | // Integral a region of n pixels around x/y
|
---|
438 | for (int dx=-n; dx<=n; dx++)
|
---|
439 | for (int dy=-n; dy<=n; dy++)
|
---|
440 | {
|
---|
441 | if (TMath::Hypot((float)dx, (float)dy)>n)
|
---|
442 | continue;
|
---|
443 |
|
---|
444 | const Int_t bin = s.GetBin(x+dx,y+dy);
|
---|
445 |
|
---|
446 | sig += h1.GetBinContent(bin);
|
---|
447 | bg += h2.GetBinContent(bin);
|
---|
448 | }
|
---|
449 |
|
---|
450 | // Scale such, that the statistical error corresponds to
|
---|
451 | // the amount of off-data used to determin the background
|
---|
452 | // model, not to the background itself. This gives
|
---|
453 | // significances as calculated from the theta-sq plot.
|
---|
454 | S = sig>0 ? MMath::SignificanceLiMaSigned(sig, bg*scale/sc, sc) : 0;
|
---|
455 | }
|
---|
456 |
|
---|
457 | s.SetBinContent(x, y, S);
|
---|
458 | }
|
---|
459 | }
|
---|
460 |
|
---|
461 | void MHDisp::Profile1D(const char *name, const TH2 &h) const
|
---|
462 | {
|
---|
463 | if (!gPad)
|
---|
464 | return;
|
---|
465 |
|
---|
466 | TH1D *hrc = dynamic_cast<TH1D*>(gPad->FindObject(name));
|
---|
467 | if (!hrc)
|
---|
468 | return;
|
---|
469 |
|
---|
470 | hrc->Reset();
|
---|
471 |
|
---|
472 | //const_cast<TH2&>(h).SetMaximum();
|
---|
473 | const Double_t max = h.GetMaximum();
|
---|
474 |
|
---|
475 | MBinning(51, -max*1.1, max*1.1).Apply(*hrc);
|
---|
476 |
|
---|
477 | for (int x=1; x<=h.GetXaxis()->GetNbins(); x++)
|
---|
478 | for (int y=1; y<=h.GetXaxis()->GetNbins(); y++)
|
---|
479 | {
|
---|
480 | const Int_t bin = h.GetBin(x,y);
|
---|
481 | const Double_t sig = h.GetBinContent(bin);
|
---|
482 | if (sig!=0)
|
---|
483 | hrc->Fill(sig);
|
---|
484 | }
|
---|
485 |
|
---|
486 | gPad->SetLogy(hrc->GetMaximum()>0);
|
---|
487 |
|
---|
488 | if (!fHistOff || hrc->GetRMS()<0.1)
|
---|
489 | return;
|
---|
490 |
|
---|
491 | // ---------- Fix mean ----------
|
---|
492 | // TF1 *g = (TF1*)gROOT->GetFunction("gaus");
|
---|
493 | // g->FixParameter(1, 0);
|
---|
494 | // hrc->Fit("gaus", "BQI");
|
---|
495 |
|
---|
496 | hrc->Fit("gaus", "QI");
|
---|
497 |
|
---|
498 | TF1 *f = hrc->GetFunction("gaus");
|
---|
499 | if (f)
|
---|
500 | {
|
---|
501 | f->SetLineWidth(1);
|
---|
502 | f->SetLineColor(kBlue);
|
---|
503 | }
|
---|
504 | }
|
---|
505 |
|
---|
506 | void MHDisp::Paint(Option_t *o)
|
---|
507 | {
|
---|
508 | // Compile Background if necessary
|
---|
509 | Finalize();
|
---|
510 |
|
---|
511 | // Paint result
|
---|
512 | TVirtualPad *pad = gPad;
|
---|
513 |
|
---|
514 | pad->cd(1);
|
---|
515 |
|
---|
516 | // Project on data onto yx-plane
|
---|
517 | fHist.GetZaxis()->SetRange(0,9999);
|
---|
518 | TH2 *h1=(TH2*)fHist.Project3D("yx_on");
|
---|
519 |
|
---|
520 | // Set Glow-palette (PaintSurface doesn't allow more than 99 colors)
|
---|
521 | MH::SetPalette(fHistOff?"glowsym":"glow1", 99);
|
---|
522 | h1->SetContour(99);
|
---|
523 |
|
---|
524 | TH2 *hx=0;
|
---|
525 |
|
---|
526 | Double_t scale = 1;
|
---|
527 | if (fHistOff)
|
---|
528 | {
|
---|
529 | // Project off data onto yx-plane and subtract it from on-data
|
---|
530 | fHistOff->GetZaxis()->SetRange(0,9999);
|
---|
531 | TH2 *h=(TH2*)fHistOff->Project3D("yx_off");
|
---|
532 |
|
---|
533 | const Double_t h1off = GetOffSignal(*h1);
|
---|
534 | const Double_t hoff = GetOffSignal(fHistBg);
|
---|
535 |
|
---|
536 | scale = hoff==0 ? -1 : -h1off/hoff;
|
---|
537 |
|
---|
538 | hx = (TH2*)pad->GetPad(4)->FindObject("Alpha_yx_sig");
|
---|
539 | if (hx)
|
---|
540 | {
|
---|
541 | hx->SetContour(99);
|
---|
542 | MakeSignificance(*hx, *h1, fHistBg, TMath::Abs(scale));
|
---|
543 | MakeDot(*hx);
|
---|
544 | MakeSymmetric(hx);
|
---|
545 | }
|
---|
546 |
|
---|
547 | h1->Add(&fHistBg, scale);
|
---|
548 |
|
---|
549 | MakeDot(*h1);
|
---|
550 | MakeSymmetric(h1);
|
---|
551 |
|
---|
552 | delete h;
|
---|
553 | }
|
---|
554 |
|
---|
555 | pad->cd(3);
|
---|
556 | TH1 *h2 = (TH1*)gPad->FindObject("RadProf");
|
---|
557 |
|
---|
558 | TString opt(o);
|
---|
559 | opt.ToLower();
|
---|
560 |
|
---|
561 | if (h1 && h2)
|
---|
562 | {
|
---|
563 | Int_t ix, iy, iz;
|
---|
564 | h1->GetMaximumBin(ix, iy, iz);
|
---|
565 |
|
---|
566 | const Double_t x0 = h1->GetXaxis()->GetBinCenter(ix);
|
---|
567 | const Double_t y0 = h1->GetYaxis()->GetBinCenter(iy);
|
---|
568 | //const Double_t w0 = h1->GetXaxis()->GetBinWidth(1);
|
---|
569 |
|
---|
570 | Profile(*h2, *h1, 0, 0);
|
---|
571 | //Profile(*h2, *h1, x0, y0);
|
---|
572 |
|
---|
573 | if (h2->GetEntries()>0)
|
---|
574 | {
|
---|
575 | // Replace with MAlphaFitter?
|
---|
576 | TF1 func("fcn", "gaus + [3]*x*x + [4]");
|
---|
577 | func.SetLineWidth(1);
|
---|
578 | func.SetLineColor(kBlue);
|
---|
579 |
|
---|
580 | func.SetParLimits(2, h2->GetBinWidth(1), 1.0);
|
---|
581 |
|
---|
582 | func.SetParameter(0, h2->GetBinContent(1));
|
---|
583 | func.FixParameter(1, 0);
|
---|
584 | func.SetParameter(2, 0.12);
|
---|
585 | if (fHistOff)
|
---|
586 | func.FixParameter(3, 0);
|
---|
587 | func.SetParameter(4, 0);//h2->GetBinContent(20));
|
---|
588 | h2->Fit(&func, "IQ", "", 0, 1.0);
|
---|
589 |
|
---|
590 | h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ f=%.2f", x0, y0, func.GetParameter(2), TMath::Abs(scale)));
|
---|
591 | }
|
---|
592 | }
|
---|
593 |
|
---|
594 | pad->cd(5);
|
---|
595 | if (h1)
|
---|
596 | Profile1D("Exc", *h1);
|
---|
597 |
|
---|
598 | pad->cd(6);
|
---|
599 | if (hx)
|
---|
600 | Profile1D("Sig", *hx);
|
---|
601 | }
|
---|
602 |
|
---|
603 | void MHDisp::Draw(Option_t *o)
|
---|
604 | {
|
---|
605 | TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
|
---|
606 | const Int_t col = pad->GetFillColor();
|
---|
607 | pad->SetBorderMode(0);
|
---|
608 |
|
---|
609 | AppendPad(o);
|
---|
610 |
|
---|
611 | // ----- Pad number 4 -----
|
---|
612 | TString name = Form("%s_4", pad->GetName());
|
---|
613 | TPad *p = new TPad(name,name, 0.5025, 0.3355, 0.995, 0.995, col, 0, 0);
|
---|
614 | p->SetNumber(4);
|
---|
615 | p->Draw();
|
---|
616 | p->cd();
|
---|
617 |
|
---|
618 | TH1 *h3 = fHist.Project3D("yx_sig");
|
---|
619 | h3->SetTitle("Significance Map");
|
---|
620 | h3->SetDirectory(NULL);
|
---|
621 | h3->SetXTitle(fHist.GetXaxis()->GetTitle());
|
---|
622 | h3->SetYTitle(fHist.GetYaxis()->GetTitle());
|
---|
623 | h3->SetMinimum(0);
|
---|
624 | h3->Draw("colz");
|
---|
625 | h3->SetBit(kCanDelete);
|
---|
626 |
|
---|
627 | if (fHistOff)
|
---|
628 | GetCatalog()->Draw("mirror same *");
|
---|
629 |
|
---|
630 | // ----- Pad number 1 -----
|
---|
631 | pad->cd();
|
---|
632 | name = Form("%s_1", pad->GetName());
|
---|
633 | p = new TPad(name,name, 0.005, 0.3355, 0.4975, 0.995, col, 0, 0);
|
---|
634 | p->SetNumber(1);
|
---|
635 | p->Draw();
|
---|
636 | p->cd();
|
---|
637 |
|
---|
638 | h3 = fHist.Project3D("yx_on");
|
---|
639 | h3->SetTitle("Distribution of equivalent events");
|
---|
640 | h3->SetDirectory(NULL);
|
---|
641 | h3->SetXTitle(fHist.GetXaxis()->GetTitle());
|
---|
642 | h3->SetYTitle(fHist.GetYaxis()->GetTitle());
|
---|
643 | h3->SetMinimum(0);
|
---|
644 | h3->Draw("colz");
|
---|
645 | h3->SetBit(kCanDelete);
|
---|
646 |
|
---|
647 | if (fHistOff)
|
---|
648 | GetCatalog()->Draw("mirror same *");
|
---|
649 |
|
---|
650 | // ----- Pad number 2 -----
|
---|
651 | pad->cd();
|
---|
652 | name = Form("%s_2", pad->GetName());
|
---|
653 | p = new TPad(name,name, 0.005, 0.005, 0.2485, 0.3315, col, 0, 0);
|
---|
654 | p->SetNumber(2);
|
---|
655 | p->Draw();
|
---|
656 | p->cd();
|
---|
657 | h3->Draw("surf3");
|
---|
658 |
|
---|
659 | // ----- Pad number 3 -----
|
---|
660 | pad->cd();
|
---|
661 | name = Form("%s_3", pad->GetName());
|
---|
662 | p = new TPad(name,name, 0.2525, 0.005, 0.4985, 0.3315, col, 0, 0);
|
---|
663 | p->SetNumber(3);
|
---|
664 | p->SetGrid();
|
---|
665 | p->Draw();
|
---|
666 | p->cd();
|
---|
667 |
|
---|
668 | const Double_t maxr = TMath::Hypot(h3->GetXaxis()->GetXmax(), h3->GetYaxis()->GetXmax());
|
---|
669 | const Int_t nbin = (h3->GetNbinsX()+h3->GetNbinsY())/2;
|
---|
670 | TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr);
|
---|
671 | h->SetDirectory(0);
|
---|
672 | //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr);
|
---|
673 | //h->Sumw2();
|
---|
674 | h->SetXTitle("\\vartheta [\\circ]");
|
---|
675 | h->SetYTitle("<cts>/\\Delta R");
|
---|
676 | h->SetBit(kCanDelete);
|
---|
677 | h->Draw();
|
---|
678 |
|
---|
679 | // ----- Pad number 5 -----
|
---|
680 | pad->cd();
|
---|
681 | name = Form("%s_5", pad->GetName());
|
---|
682 | p = new TPad(name,name, 0.5025, 0.005, 0.7485, 0.3315, col, 0, 0);
|
---|
683 | p->SetNumber(5);
|
---|
684 | p->SetGrid();
|
---|
685 | p->Draw();
|
---|
686 | p->cd();
|
---|
687 |
|
---|
688 | h3 = new TH1D("Exc", "Distribution of excess events", 10, -1, 1);
|
---|
689 | h3->SetDirectory(NULL);
|
---|
690 | h3->SetXTitle("N");
|
---|
691 | h3->SetYTitle("Counts");
|
---|
692 | h3->Draw();
|
---|
693 | h3->SetBit(kCanDelete);
|
---|
694 |
|
---|
695 | // ----- Pad number 6 -----
|
---|
696 | pad->cd();
|
---|
697 | name = Form("%s_6", pad->GetName());
|
---|
698 | p = new TPad(name,name, 0.7525, 0.005, 0.9995, 0.3315, col, 0, 0);
|
---|
699 | p->SetNumber(6);
|
---|
700 | p->SetGrid();
|
---|
701 | p->Draw();
|
---|
702 | p->cd();
|
---|
703 |
|
---|
704 | h3 = new TH1D("Sig", "Distribution of significances", 10, -1, 1);
|
---|
705 | h3->SetDirectory(NULL);
|
---|
706 | h3->SetXTitle("N");
|
---|
707 | h3->SetYTitle("Counts");
|
---|
708 | h3->Draw();
|
---|
709 | h3->SetBit(kCanDelete);
|
---|
710 | }
|
---|
711 |
|
---|
712 | // --------------------------------------------------------------------------
|
---|
713 | //
|
---|
714 | // The following resources are available:
|
---|
715 | //
|
---|
716 | // MHDisp.Smearing: 0.1
|
---|
717 | // MHDisp.Wobble: on/off
|
---|
718 | //
|
---|
719 | Int_t MHDisp::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
|
---|
720 | {
|
---|
721 | Int_t rc = MHFalseSource::ReadEnv(env, prefix, print);
|
---|
722 | if (rc==kERROR)
|
---|
723 | return kERROR;
|
---|
724 |
|
---|
725 | if (IsEnvDefined(env, prefix, "Smearing", print))
|
---|
726 | {
|
---|
727 | rc = kTRUE;
|
---|
728 | SetSmearing(GetEnvValue(env, prefix, "Smearing", fSmearing));
|
---|
729 | }
|
---|
730 | if (IsEnvDefined(env, prefix, "Wobble", print))
|
---|
731 | {
|
---|
732 | rc = kTRUE;
|
---|
733 | SetWobble(GetEnvValue(env, prefix, "Wobble", fWobble));
|
---|
734 | }
|
---|
735 |
|
---|
736 | return rc;
|
---|
737 | }
|
---|
738 |
|
---|