| 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): Wolfgang Wittek 1/2003 <mailto:wittek@mppmu.mpg.de>
|
|---|
| 19 | ! Author(s): Thomas Bretz 4/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
|
|---|
| 20 | !
|
|---|
| 21 | ! Copyright: MAGIC Software Development, 2000-2003
|
|---|
| 22 | !
|
|---|
| 23 | !
|
|---|
| 24 | \* ======================================================================== */
|
|---|
| 25 |
|
|---|
| 26 | //////////////////////////////////////////////////////////////////////////////
|
|---|
| 27 | // //
|
|---|
| 28 | // MHSigmaTheta //
|
|---|
| 29 | // //
|
|---|
| 30 | // calculates - the 2D-histogram sigmabar vs. Theta (Inner), and //
|
|---|
| 31 | // - the 2D-histogram sigmabar vs. Theta (Outer) //
|
|---|
| 32 | // - the 3D-histogram sigma, pixel no., Theta //
|
|---|
| 33 | // - the 3D-histogram (sigma^2-sigmabar^2), pixel no., Theta //
|
|---|
| 34 | // //
|
|---|
| 35 | //////////////////////////////////////////////////////////////////////////////
|
|---|
| 36 |
|
|---|
| 37 | #include "MHSigmaTheta.h"
|
|---|
| 38 |
|
|---|
| 39 | #include <TCanvas.h>
|
|---|
| 40 |
|
|---|
| 41 | #include "MTime.h"
|
|---|
| 42 | #include "MPointingPos.h"
|
|---|
| 43 |
|
|---|
| 44 | #include "MBinning.h"
|
|---|
| 45 | #include "MParList.h"
|
|---|
| 46 |
|
|---|
| 47 | #include "MGeomCam.h"
|
|---|
| 48 | #include "MGeomPix.h"
|
|---|
| 49 | #include "MBadPixelsCam.h"
|
|---|
| 50 |
|
|---|
| 51 | #include "MPedPhotCam.h"
|
|---|
| 52 | #include "MPedPhotPix.h"
|
|---|
| 53 |
|
|---|
| 54 | #include "MCerPhotEvt.h"
|
|---|
| 55 | #include "MCerPhotPix.h"
|
|---|
| 56 |
|
|---|
| 57 | #include "MLog.h"
|
|---|
| 58 | #include "MLogManip.h"
|
|---|
| 59 |
|
|---|
| 60 | ClassImp(MHSigmaTheta);
|
|---|
| 61 |
|
|---|
| 62 | using namespace std;
|
|---|
| 63 |
|
|---|
| 64 | // --------------------------------------------------------------------------
|
|---|
| 65 | //
|
|---|
| 66 | // Default Constructor. It sets name and title of the histogram.
|
|---|
| 67 | //
|
|---|
| 68 | MHSigmaTheta::MHSigmaTheta(const char *name, const char *title)
|
|---|
| 69 | {
|
|---|
| 70 | fName = name ? name : "MHSigmaTheta";
|
|---|
| 71 | fTitle = title ? title : "2D histogram sigmabar vs. Theta";
|
|---|
| 72 |
|
|---|
| 73 | fSigmaTheta.SetDirectory(NULL);
|
|---|
| 74 | fSigmaTheta.SetName("2D-ThetaSigmabar(Inner)");
|
|---|
| 75 | fSigmaTheta.SetTitle("2D: \\bar{\\sigma}, \\Theta");
|
|---|
| 76 | fSigmaTheta.SetXTitle("\\Theta [\\circ]");
|
|---|
| 77 | fSigmaTheta.SetYTitle("Sigmabar(Inner) / SQRT(Area)");
|
|---|
| 78 |
|
|---|
| 79 | fSigmaThetaOuter.SetDirectory(NULL);
|
|---|
| 80 | fSigmaThetaOuter.SetName("2D-ThetaSigmabar(Outer)");
|
|---|
| 81 | fSigmaThetaOuter.SetTitle("2D: \\bar{\\sigma}, \\Theta");
|
|---|
| 82 | fSigmaThetaOuter.SetXTitle("\\Theta [\\circ]");
|
|---|
| 83 | fSigmaThetaOuter.SetYTitle("Sigmabar(Outer) / SQRT(Area)");
|
|---|
| 84 |
|
|---|
| 85 | fSigmaPixTheta.SetDirectory(NULL);
|
|---|
| 86 | fSigmaPixTheta.SetName("3D-ThetaPixSigma");
|
|---|
| 87 | fSigmaPixTheta.SetTitle("3D: \\Theta, Pixel Id, \\sigma");
|
|---|
| 88 | fSigmaPixTheta.SetXTitle("\\Theta [\\circ]");
|
|---|
| 89 | fSigmaPixTheta.SetYTitle("Pixel Id");
|
|---|
| 90 | fSigmaPixTheta.SetZTitle("Sigma");
|
|---|
| 91 |
|
|---|
| 92 | fDiffPixTheta.SetDirectory(NULL);
|
|---|
| 93 | fDiffPixTheta.SetName("3D-ThetaPixDiff");
|
|---|
| 94 | fDiffPixTheta.SetTitle("3D: \\Theta, Pixel Id, {\\sigma}^{2}-\\bar{\\sigma}^{2}");
|
|---|
| 95 | fDiffPixTheta.SetXTitle("\\Theta [\\circ]");
|
|---|
| 96 | fDiffPixTheta.SetYTitle("Pixel Id");
|
|---|
| 97 | fDiffPixTheta.SetZTitle("(Sigma2 - Sigmabar2)/Area");
|
|---|
| 98 |
|
|---|
| 99 | // Set default binning
|
|---|
| 100 | // FIXME: Maybe ist's necessary to adapt the value to the
|
|---|
| 101 | // Magic default values
|
|---|
| 102 | MBinning binsb;
|
|---|
| 103 | MBinning binst;
|
|---|
| 104 | MBinning binsd;
|
|---|
| 105 | binsd.SetEdges(100, -10, 20);
|
|---|
| 106 | binsb.SetEdges(100, 0, 10);
|
|---|
| 107 | binst.SetEdgesCos(10, 0, 90);
|
|---|
| 108 |
|
|---|
| 109 | MBinning binspix("BinningPixel");
|
|---|
| 110 | binspix.SetEdges(578, -0.5, 577.5);
|
|---|
| 111 |
|
|---|
| 112 | SetBinning(&fSigmaTheta, &binst, &binsb);
|
|---|
| 113 | SetBinning(&fSigmaThetaOuter, &binst, &binsb);
|
|---|
| 114 | SetBinning(&fSigmaPixTheta, &binst, &binspix, &binsb);
|
|---|
| 115 | SetBinning(&fDiffPixTheta, &binst, &binspix, &binsd);
|
|---|
| 116 |
|
|---|
| 117 | fNamePedPhotCam = "MPedPhotCamFromData";
|
|---|
| 118 | }
|
|---|
| 119 |
|
|---|
| 120 | // --------------------------------------------------------------------------
|
|---|
| 121 | //
|
|---|
| 122 | // Set the binnings and prepare the filling of the histogram
|
|---|
| 123 | //
|
|---|
| 124 | Bool_t MHSigmaTheta::SetupFill(const MParList *plist)
|
|---|
| 125 | {
|
|---|
| 126 | fCam = (MGeomCam*)plist->FindObject("MGeomCam");
|
|---|
| 127 | if (!fCam)
|
|---|
| 128 | {
|
|---|
| 129 | *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
|
|---|
| 130 | return kFALSE;
|
|---|
| 131 | }
|
|---|
| 132 |
|
|---|
| 133 | fPointPos = (MPointingPos*)plist->FindObject("MPointingPos");
|
|---|
| 134 | if (!fPointPos)
|
|---|
| 135 | *fLog << warn << "MPointingPos not found... aborting." << endl;
|
|---|
| 136 |
|
|---|
| 137 |
|
|---|
| 138 |
|
|---|
| 139 | fPed = (MPedPhotCam*)plist->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam");
|
|---|
| 140 | if (!fPed)
|
|---|
| 141 | {
|
|---|
| 142 | *fLog << err << AddSerialNumber(fNamePedPhotCam)
|
|---|
| 143 | << "[MPedPhotCam] not found... aborting." << endl;
|
|---|
| 144 | return kFALSE;
|
|---|
| 145 | }
|
|---|
| 146 | //fPed->InitSize(fCam->GetNumPixels());
|
|---|
| 147 |
|
|---|
| 148 |
|
|---|
| 149 | fBad = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam");
|
|---|
| 150 | if (!fBad)
|
|---|
| 151 | {
|
|---|
| 152 | *fLog << err << "MBadPixelsCam not found... continue. " << endl;
|
|---|
| 153 | }
|
|---|
| 154 |
|
|---|
| 155 |
|
|---|
| 156 | fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt");
|
|---|
| 157 | if (!fEvt)
|
|---|
| 158 | {
|
|---|
| 159 | *fLog << err << "MCerPhotEvt not found... aborting." << endl;
|
|---|
| 160 | return kFALSE;
|
|---|
| 161 | }
|
|---|
| 162 |
|
|---|
| 163 |
|
|---|
| 164 | // Get Theta Binning
|
|---|
| 165 | MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta", "MBinning");
|
|---|
| 166 | if (!binstheta)
|
|---|
| 167 | {
|
|---|
| 168 | *fLog << warn << "Object 'BinningTheta' [MBinning] not found... no binning applied." << endl;
|
|---|
| 169 | return kTRUE;
|
|---|
| 170 | }
|
|---|
| 171 |
|
|---|
| 172 | // Get Sigmabar binning
|
|---|
| 173 | MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigmabar", "MBinning");
|
|---|
| 174 | if (!binssigma)
|
|---|
| 175 | *fLog << warn << "Object 'BinningSigmabar' [MBinning] not found... no binning applied." << endl;
|
|---|
| 176 |
|
|---|
| 177 | // Get binning for (sigma^2-sigmabar^2)
|
|---|
| 178 | MBinning* binsdiff = (MBinning*)plist->FindObject("BinningDiffsigma2", "MBinning");
|
|---|
| 179 | if (!binsdiff)
|
|---|
| 180 | *fLog << warn << "Object 'BinningDiffsigma2' [MBinning] not found... no binning applied." << endl;
|
|---|
| 181 |
|
|---|
| 182 | //FIXME: Missing: Apply new binning to only one axis!
|
|---|
| 183 |
|
|---|
| 184 | // Get binning for pixel number
|
|---|
| 185 | const UInt_t npix1 = fPed->GetSize()+1;
|
|---|
| 186 | *fLog << "npix1 = " << npix1 << endl;
|
|---|
| 187 | MBinning binspix("BinningPixel");
|
|---|
| 188 | binspix.SetEdges(npix1, -0.5, npix1-0.5);
|
|---|
| 189 |
|
|---|
| 190 | // Set binnings in histograms
|
|---|
| 191 | if (binssigma)
|
|---|
| 192 | {
|
|---|
| 193 | SetBinning(&fSigmaTheta, binstheta, binssigma);
|
|---|
| 194 | SetBinning(&fSigmaThetaOuter, binstheta, binssigma);
|
|---|
| 195 | SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma);
|
|---|
| 196 | }
|
|---|
| 197 |
|
|---|
| 198 | if (binsdiff)
|
|---|
| 199 | SetBinning(&fDiffPixTheta, binstheta, &binspix, binsdiff);
|
|---|
| 200 |
|
|---|
| 201 | return kTRUE;
|
|---|
| 202 | }
|
|---|
| 203 |
|
|---|
| 204 | // --------------------------------------------------------------------------
|
|---|
| 205 | //
|
|---|
| 206 | // Fill the histograms
|
|---|
| 207 | //
|
|---|
| 208 | // ignore pixels if they are unused or blind
|
|---|
| 209 | //
|
|---|
| 210 | Bool_t MHSigmaTheta::Fill(const MParContainer *par, const Stat_t w)
|
|---|
| 211 | {
|
|---|
| 212 | Double_t theta = fPointPos->GetZd();
|
|---|
| 213 |
|
|---|
| 214 | Double_t mysig = (fPed->GetArea(0)).GetRms();
|
|---|
| 215 | Double_t mysigouter = (fPed->GetArea(1)).GetRms();
|
|---|
| 216 |
|
|---|
| 217 | *fLog << "theta, mysig, mysigouter = " << theta << ", " << mysig
|
|---|
| 218 | << ", " << mysigouter << endl;
|
|---|
| 219 |
|
|---|
| 220 | fSigmaTheta.Fill(theta, mysig);
|
|---|
| 221 | fSigmaThetaOuter.Fill(theta, mysigouter);
|
|---|
| 222 |
|
|---|
| 223 | const UInt_t npix = fEvt->GetNumPixels();
|
|---|
| 224 |
|
|---|
| 225 | for (UInt_t i=0; i<npix; i++)
|
|---|
| 226 | {
|
|---|
| 227 | MCerPhotPix &cerpix = (*fEvt)[i];
|
|---|
| 228 | const Int_t id = cerpix.GetPixId();
|
|---|
| 229 |
|
|---|
| 230 | if (!cerpix.IsPixelUsed())
|
|---|
| 231 | {
|
|---|
| 232 | //*fLog << all << "MHSigmaTheta::Fill; unused pixel found, id = "
|
|---|
| 233 | // << id << endl;
|
|---|
| 234 | continue;
|
|---|
| 235 | }
|
|---|
| 236 |
|
|---|
| 237 | const MPedPhotPix &pix = (*fPed)[id];
|
|---|
| 238 |
|
|---|
| 239 | if ( fBad != NULL && ((*fBad)[id]).IsUnsuitable() )
|
|---|
| 240 | {
|
|---|
| 241 | // this should never occur, because bad pixels should have
|
|---|
| 242 | // been set unused
|
|---|
| 243 | *fLog << all << "MHSigmaTheta::Fill; bad pixel found which is used, id = "
|
|---|
| 244 | << id << "... go to next pixel." << endl;
|
|---|
| 245 | continue;
|
|---|
| 246 | }
|
|---|
| 247 |
|
|---|
| 248 | // ratio is the area of pixel 0
|
|---|
| 249 | // divided by the area of the current pixel
|
|---|
| 250 | const Double_t ratio = fCam->GetPixRatio(id);
|
|---|
| 251 | const Double_t sigma = pix.GetRms();
|
|---|
| 252 |
|
|---|
| 253 | fSigmaPixTheta.Fill(theta, (Double_t)id, sigma*sqrt(ratio));
|
|---|
| 254 |
|
|---|
| 255 | Double_t diff;
|
|---|
| 256 | const Byte_t aidx = (*fCam)[id].GetAidx();
|
|---|
| 257 | if (aidx == 0)
|
|---|
| 258 | {
|
|---|
| 259 | // inner pixel
|
|---|
| 260 | diff = sigma*sigma*ratio - mysig*mysig;
|
|---|
| 261 | }
|
|---|
| 262 | else
|
|---|
| 263 | {
|
|---|
| 264 | // outer pixel
|
|---|
| 265 | diff = sigma*sigma*ratio - mysigouter*mysigouter;
|
|---|
| 266 | }
|
|---|
| 267 |
|
|---|
| 268 | fDiffPixTheta.Fill(theta, (Double_t)id, diff);
|
|---|
| 269 | }
|
|---|
| 270 |
|
|---|
| 271 | return kTRUE;
|
|---|
| 272 | }
|
|---|
| 273 |
|
|---|
| 274 | // --------------------------------------------------------------------------
|
|---|
| 275 | //
|
|---|
| 276 | // Update the projections and (if possible) set log scales before painting
|
|---|
| 277 | //
|
|---|
| 278 | void MHSigmaTheta::Paint(Option_t *opt)
|
|---|
| 279 | {
|
|---|
| 280 | TVirtualPad *padsave = gPad;
|
|---|
| 281 |
|
|---|
| 282 | TH1D* h;
|
|---|
| 283 |
|
|---|
| 284 | padsave->cd(1);
|
|---|
| 285 | if ((h = (TH1D*)gPad->FindObject("ProjX-Theta")))
|
|---|
| 286 | {
|
|---|
| 287 | ProjectionX(*h, fSigmaTheta);
|
|---|
| 288 | if (h->GetEntries()!=0)
|
|---|
| 289 | gPad->SetLogy();
|
|---|
| 290 | }
|
|---|
| 291 |
|
|---|
| 292 | padsave->cd(4);
|
|---|
| 293 | if ((h = (TH1D*)gPad->FindObject("ProjY-sigma")))
|
|---|
| 294 | ProjectionY(*h, fSigmaTheta);
|
|---|
| 295 |
|
|---|
| 296 | gPad = padsave;
|
|---|
| 297 | }
|
|---|
| 298 |
|
|---|
| 299 | // --------------------------------------------------------------------------
|
|---|
| 300 | //
|
|---|
| 301 | // Draw the histogram
|
|---|
| 302 | //
|
|---|
| 303 | void MHSigmaTheta::Draw(Option_t *opt)
|
|---|
| 304 | {
|
|---|
| 305 | TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
|
|---|
| 306 | pad->SetBorderMode(0);
|
|---|
| 307 |
|
|---|
| 308 | AppendPad("");
|
|---|
| 309 |
|
|---|
| 310 | pad->Divide(3, 3);
|
|---|
| 311 |
|
|---|
| 312 | // draw the 2D histogram Sigmabar versus Theta
|
|---|
| 313 | TH1 *h;
|
|---|
| 314 |
|
|---|
| 315 | pad->cd(1);
|
|---|
| 316 | gPad->SetBorderMode(0);
|
|---|
| 317 | h = fSigmaTheta.ProjectionX("ProjX-Theta", -1, 9999, "E");
|
|---|
| 318 | h->SetDirectory(NULL);
|
|---|
| 319 | h->SetTitle("Distribution of \\Theta");
|
|---|
| 320 | h->SetXTitle("\\Theta [\\circ]");
|
|---|
| 321 | h->SetYTitle("No.of events");
|
|---|
| 322 | h->Draw(opt);
|
|---|
| 323 | h->SetBit(kCanDelete);
|
|---|
| 324 |
|
|---|
| 325 | pad->cd(2);
|
|---|
| 326 | gPad->SetBorderMode(0);
|
|---|
| 327 | h = fDiffPixTheta.Project3D("zx");
|
|---|
| 328 | h->SetDirectory(NULL);
|
|---|
| 329 | h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)");
|
|---|
| 330 | h->SetXTitle("\\Theta [\\circ]");
|
|---|
| 331 | h->SetYTitle("(Sigma2 - Sigmabar2) / Area");
|
|---|
| 332 | h->Draw("box");
|
|---|
| 333 | h->SetBit(kCanDelete);
|
|---|
| 334 |
|
|---|
| 335 | pad->cd(3);
|
|---|
| 336 | gPad->SetBorderMode(0);
|
|---|
| 337 | h = fSigmaPixTheta.Project3D("zx");
|
|---|
| 338 | h->SetDirectory(NULL);
|
|---|
| 339 | h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)");
|
|---|
| 340 | h->SetXTitle("\\Theta [\\circ]");
|
|---|
| 341 | h->SetYTitle("Sigma");
|
|---|
| 342 | h->Draw("box");
|
|---|
| 343 | h->SetBit(kCanDelete);
|
|---|
| 344 |
|
|---|
| 345 | //pad->cd(7);
|
|---|
| 346 | //gPad->SetBorderMode(0);
|
|---|
| 347 | //h = fSigmaTheta.ProjectionY("ProjY-sigma", -1, 9999, "E");
|
|---|
| 348 | //h->SetDirectory(NULL);
|
|---|
| 349 | //h->SetTitle("Distribution of \\bar{\\sigma}_{ped}");
|
|---|
| 350 | //h->SetXTitle("\\bar{\\sigma}_{ped}");
|
|---|
| 351 | //h->SetYTitle("No.of events");
|
|---|
| 352 | //h->Draw(opt);
|
|---|
| 353 | //h->SetBit(kCanDelete);
|
|---|
| 354 |
|
|---|
| 355 | pad->cd(5);
|
|---|
| 356 | gPad->SetBorderMode(0);
|
|---|
| 357 | h = fDiffPixTheta.Project3D("zy");
|
|---|
| 358 | h->SetDirectory(NULL);
|
|---|
| 359 | h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all \\Theta)");
|
|---|
| 360 | h->SetXTitle("Pixel Id");
|
|---|
| 361 | h->SetYTitle("Sigma2 - sigmabar2) / Area");
|
|---|
| 362 | h->Draw("box");
|
|---|
| 363 | h->SetBit(kCanDelete);
|
|---|
| 364 |
|
|---|
| 365 | pad->cd(6);
|
|---|
| 366 | gPad->SetBorderMode(0);
|
|---|
| 367 | h = fSigmaPixTheta.Project3D("zy");
|
|---|
| 368 | h->SetDirectory(NULL);
|
|---|
| 369 | h->SetTitle("\\sigma_{ped} vs. pixel Id (all \\Theta)");
|
|---|
| 370 | h->SetXTitle("Pixel Id");
|
|---|
| 371 | h->SetYTitle("Sigma");
|
|---|
| 372 | h->Draw("box");
|
|---|
| 373 | h->SetBit(kCanDelete);
|
|---|
| 374 |
|
|---|
| 375 | pad->cd(4);
|
|---|
| 376 | fSigmaTheta.Draw(opt);
|
|---|
| 377 |
|
|---|
| 378 | pad->cd(7);
|
|---|
| 379 | fSigmaThetaOuter.Draw(opt);
|
|---|
| 380 |
|
|---|
| 381 | //pad->cd(8);
|
|---|
| 382 | //fDiffPixTheta.Draw(opt);
|
|---|
| 383 |
|
|---|
| 384 | //pad->cd(9);
|
|---|
| 385 | //fSigmaPixTheta.Draw(opt);
|
|---|
| 386 | }
|
|---|
| 387 |
|
|---|
| 388 |
|
|---|
| 389 |
|
|---|
| 390 |
|
|---|
| 391 |
|
|---|
| 392 |
|
|---|
| 393 |
|
|---|
| 394 |
|
|---|
| 395 |
|
|---|
| 396 |
|
|---|
| 397 |
|
|---|
| 398 |
|
|---|
| 399 |
|
|---|
| 400 |
|
|---|
