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