| 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, 12/2000 <mailto:tbretz@uni-sw.gwdg.de>
|
|---|
| 19 | ! Author(s): Harald Kornmayer, 1/2001
|
|---|
| 20 | ! Author(s): Nadia Tonello, 4/2003 <mailto:tonello@mppmu.mpg.de>
|
|---|
| 21 | !
|
|---|
| 22 | ! Copyright: MAGIC Software Development, 2000-2003
|
|---|
| 23 | !
|
|---|
| 24 | !
|
|---|
| 25 | \* ======================================================================== */
|
|---|
| 26 |
|
|---|
| 27 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 28 | //
|
|---|
| 29 | // MImgCleanStd
|
|---|
| 30 | //
|
|---|
| 31 | // The Image Cleaning task selects the pixels you use for the Hillas
|
|---|
| 32 | // parameters calculation.
|
|---|
| 33 | //
|
|---|
| 34 | // There are two methods to make the selection: the standard one, as done
|
|---|
| 35 | // in the analysis of CT1 data, and the democratic one, as suggested by
|
|---|
| 36 | // W.Wittek. The number of photo-electrons of a pixel is compared with the
|
|---|
| 37 | // pedestal RMS of the pixel itself (standard method) or with the average
|
|---|
| 38 | // RMS of the inner pixels (democratic method).
|
|---|
| 39 | // In both cases, the possibility to have a camera with pixels of
|
|---|
| 40 | // different area is taken into account.
|
|---|
| 41 | // The too noisy pixels can be recognized and eventally switched off
|
|---|
| 42 | // (Unmap: set blind pixels to UNUSED) separately, using the
|
|---|
| 43 | // MBlindPixelCalc Class. In the MBlindPixelCalc class there is also the
|
|---|
| 44 | // function to replace the value of the noisy pixels with the interpolation
|
|---|
| 45 | // of the content of the neighbors (SetUseInterpolation).
|
|---|
| 46 | //
|
|---|
| 47 | // Example:
|
|---|
| 48 | // ...
|
|---|
| 49 | // MBlindPixelCalc blind;
|
|---|
| 50 | // blind.SetUseInterpolation();
|
|---|
| 51 | // blind.SetUseBlindPixels();
|
|---|
| 52 | //
|
|---|
| 53 | // MImgCleanStd clean;
|
|---|
| 54 | // ...
|
|---|
| 55 | // tlist.AddToList(&blind);
|
|---|
| 56 | // tlist.AddToList(&clean);
|
|---|
| 57 | //
|
|---|
| 58 | // Look at the MBlindPixelCalc Class for more details.
|
|---|
| 59 | //
|
|---|
| 60 | // Starting point: default values ----------------------------------------
|
|---|
| 61 | //
|
|---|
| 62 | // When an event is read, before the image cleaning, all the pixels that
|
|---|
| 63 | // are in MCerPhotEvt are set as USED and NOT CORE. All the pixels belong
|
|---|
| 64 | // to RING number 1 (like USED pixels).
|
|---|
| 65 | // Look at MCerPhotPix.h to see how these informations of the pixel are
|
|---|
| 66 | // stored.
|
|---|
| 67 | // The default cleaning METHOD is the STANDARD one and the number of the
|
|---|
| 68 | // rings around the CORE pixel it analyzes is 1. Look at the Constructor
|
|---|
| 69 | // of the class in MImgCleanStd.cc to see (or change) the default values.
|
|---|
| 70 | //
|
|---|
| 71 | // Example: To modify this setting, use the member functions
|
|---|
| 72 | // SetMethod(MImgCleanStd::kDemocratic) and SetCleanRings(UShort_t n).
|
|---|
| 73 | //
|
|---|
| 74 | // MImgCleanStd:CleanStep1 -----------------------------------------------
|
|---|
| 75 | //
|
|---|
| 76 | // The first step of cleaning defines the CORE pixels. The CORE pixels are
|
|---|
| 77 | // the ones which contain the informations about the core of the electro-
|
|---|
| 78 | // magnetic shower.
|
|---|
| 79 | // The ratio (A_0/A_i) is calculated from fCam->GetPixRatio(i). A_0 is
|
|---|
| 80 | // the area of the central pixel of the camera, A_i is the area of the
|
|---|
| 81 | // examined pixel. In this way, if we have a MAGIC-like camera, with the
|
|---|
| 82 | // outer pixels bigger than the inner ones, the level of cleaning in the
|
|---|
| 83 | // two different regions is weighted.
|
|---|
| 84 | // This avoids problems of deformations of the shower images.
|
|---|
| 85 | // The signal S_i and the pedestal RMS Prms_i of the pixel are called from
|
|---|
| 86 | // the object MCerPhotPix.
|
|---|
| 87 | // If (default method = kStandard)
|
|---|
| 88 | //Begin_Html
|
|---|
| 89 | // <img src="MImgCleanStd-f1.png">
|
|---|
| 90 | //End_Html
|
|---|
| 91 | // the pixel is set as CORE pixel. L_1 (n=1) is called "first level of
|
|---|
| 92 | // cleaning" (default: fCleanLvl1 = 3).
|
|---|
| 93 | // All the other pixels are set as UNUSED and belong to RING 0.
|
|---|
| 94 | // After this point, only the CORE pixels are set as USED, with RING
|
|---|
| 95 | // number 1.
|
|---|
| 96 | //
|
|---|
| 97 | // MImgCleanStd:CleanStep2 ----------------------------------------------
|
|---|
| 98 | //
|
|---|
| 99 | // The second step of cleaning looks at the isolated CORE pixels and sets
|
|---|
| 100 | // them to UNUSED. An isolated pixel is a pixel without CORE neighbors.
|
|---|
| 101 | // At the end of this point, we have set as USED only CORE pixels with at
|
|---|
| 102 | // least one CORE neighbor.
|
|---|
| 103 | //
|
|---|
| 104 | // MImgCleanStd:CleanStep3 ----------------------------------------------
|
|---|
| 105 | //
|
|---|
| 106 | // The third step of cleaning looks at all the pixels (USED or UNUSED) that
|
|---|
| 107 | // surround the USED pixels.
|
|---|
| 108 | // If the content of the analyzed pixel survives at the second level of
|
|---|
| 109 | // cleaning, i.e. if
|
|---|
| 110 | //Begin_Html
|
|---|
| 111 | // <img src="MImgCleanStd-f1.png">
|
|---|
| 112 | //End_Html
|
|---|
| 113 | // the pixel is set as USED. L_2 (n=2) is called "second level of cleaning"
|
|---|
| 114 | // (default:fCleanLvl2 = 2.5).
|
|---|
| 115 | //
|
|---|
| 116 | // When the number of RINGS to analyze is 1 (default value), only the
|
|---|
| 117 | // pixels that have a neighbor CORE pixel are analyzed.
|
|---|
| 118 | //
|
|---|
| 119 | // There is the option to decide the number of times you want to repeat
|
|---|
| 120 | // this procedure (number of RINGS analyzed around the core pixels = n).
|
|---|
| 121 | // Every time the level of cleaning is the same (fCleanLvl2) and the pixel
|
|---|
| 122 | // will belong to ring r+1, 1 < r < n+1. This is described in
|
|---|
| 123 | // MImgCleanStd:CleanStep4 .
|
|---|
| 124 | //
|
|---|
| 125 | // Dictionary and member functions ---------------------------------------
|
|---|
| 126 | //
|
|---|
| 127 | // Here there is the detailed description of the member functions and of
|
|---|
| 128 | // the terms commonly used in the class.
|
|---|
| 129 | //
|
|---|
| 130 | // STANDARD CLEANING:
|
|---|
| 131 | // =================
|
|---|
| 132 | // This is the method used for the CT1 data analysis. It is the default
|
|---|
| 133 | // method of the class.
|
|---|
| 134 | // The number of photo-electrons of a pixel (S_i) is compared to the
|
|---|
| 135 | // pedestal RMS of the pixel itself (Prms_i). To have the comparison to
|
|---|
| 136 | // the same photon density for all the pixels, taking into account they
|
|---|
| 137 | // can have different areas, we have to keep in mind that the number of
|
|---|
| 138 | // photons that hit each pixel, goes linearly with the area of the pixel.
|
|---|
| 139 | // The fluctuations of the LONS are proportional to sqrt(A_i), so when we
|
|---|
| 140 | // compare S_i with Prms_i, only a factor sqrt(A_0/A_i) is missing to
|
|---|
| 141 | // have the same (N.photons/Area) threshold for all the pixels.
|
|---|
| 142 | //
|
|---|
| 143 | // !!WARNING: if noise independent from the
|
|---|
| 144 | // pixel size is considered,
|
|---|
| 145 | // this weight can give wrong results!!
|
|---|
| 146 | // If
|
|---|
| 147 | //Begin_Html
|
|---|
| 148 | // <img src="MImgCleanStd-f1.png">
|
|---|
| 149 | //End_Html
|
|---|
| 150 | // the pixel survives the cleaning and it is set as CORE (when L_n is the
|
|---|
| 151 | // first level of cleaning, fCleanLvl1) or USED (when L_n is the second
|
|---|
| 152 | // level of cleaning, fCleanLvl2).
|
|---|
| 153 | //
|
|---|
| 154 | // Example:
|
|---|
| 155 | //
|
|---|
| 156 | // MImgCleanStd clean;
|
|---|
| 157 | // //creates a default Cleaning object, with default setting
|
|---|
| 158 | // ...
|
|---|
| 159 | // tlist.AddToList(&clean);
|
|---|
| 160 | // // add the image cleaning to the main task list
|
|---|
| 161 | //
|
|---|
| 162 | // DEMOCRATIC CLEANING:
|
|---|
| 163 | // ===================
|
|---|
| 164 | // You use this cleaning method when you want to compare the number of
|
|---|
| 165 | // photo-electons of each pixel with the average pedestal RMS
|
|---|
| 166 | // (fInnerNoise = fSgb->GetSigmabarInner()) of the inner pixels (for the
|
|---|
| 167 | // MAGIC camera they are the smaller ones):
|
|---|
| 168 | //Begin_Html
|
|---|
| 169 | // <img src="MImgCleanStd-f2.png">
|
|---|
| 170 | //End_Html
|
|---|
| 171 | // In this case, the simple ratio (A_0/A_i) is used to weight the level of
|
|---|
| 172 | // cleaning, because both the inner and the outer pixels (that in MAGIC
|
|---|
| 173 | // have a different area) are compared to the same pedestal RMS, coming
|
|---|
| 174 | // from the inner pixels.
|
|---|
| 175 | // To calculate the average pedestal RMS of the inner pixels, you have to
|
|---|
| 176 | // add to the main task list an object of type MSigmabarCalc before the
|
|---|
| 177 | // MImgCleanStd object. To know how the calculation of fInnerNoise is done
|
|---|
| 178 | // look at the MSigmabarCalc Class.
|
|---|
| 179 | //
|
|---|
| 180 | // Example:
|
|---|
| 181 | //
|
|---|
| 182 | // MSigmabarCalc sbcalc;
|
|---|
| 183 | // //creates an object that calcutates the average pedestal RMS
|
|---|
| 184 | // MImgCleanStd clean;
|
|---|
| 185 | // ...
|
|---|
| 186 | // tlist.AddToList(&sbcalc);
|
|---|
| 187 | // tlist.AddToList(&clean);
|
|---|
| 188 | //
|
|---|
| 189 | // Member Function: SetMethod()
|
|---|
| 190 | // ============================
|
|---|
| 191 | // When you call the MImgCleanStd task, the default method is kStandard.
|
|---|
| 192 | //
|
|---|
| 193 | // If you want to switch to the kDemocratic method you have to
|
|---|
| 194 | // call this member function.
|
|---|
| 195 | //
|
|---|
| 196 | // Example:
|
|---|
| 197 | //
|
|---|
| 198 | // MImgCleanStd clean;
|
|---|
| 199 | // //creates a default Cleaning object, with default setting
|
|---|
| 200 | //
|
|---|
| 201 | // clean.SetMethod(MImgCleanStd::kDemocratic);
|
|---|
| 202 | // //now the method of cleaning is changed to Democratic
|
|---|
| 203 | //
|
|---|
| 204 | // FIRST AND SECOND CLEANING LEVEL
|
|---|
| 205 | // ===============================
|
|---|
| 206 | // When you call the MImgCleanStd task, the default cleaning levels are
|
|---|
| 207 | // fCleanLvl1 = 3, fCleanLvl2 = 2.5. You can change them easily when you
|
|---|
| 208 | // create the MImgCleanStd object.
|
|---|
| 209 | //
|
|---|
| 210 | // Example:
|
|---|
| 211 | //
|
|---|
| 212 | // MImgCleanStd clean(Float_t lvl1,Float_t lvl2);
|
|---|
| 213 | // //creates a default cleaning object, but the cleaning levels are now
|
|---|
| 214 | // //lvl1 and lvl2.
|
|---|
| 215 | //
|
|---|
| 216 | // RING NUMBER
|
|---|
| 217 | // ===========
|
|---|
| 218 | // The standard cleaning procedure is such that it looks for the
|
|---|
| 219 | // informations of the boundary part of the shower only on the first
|
|---|
| 220 | // neighbors of the CORE pixels.
|
|---|
| 221 | // There is the possibility now to look not only at the firs neighbors
|
|---|
| 222 | // (first ring),but also further away, around the CORE pixels. All the new
|
|---|
| 223 | // pixels you can find with this method, are tested with the second level
|
|---|
| 224 | // of cleaning and have to have at least an USED neighbor.
|
|---|
| 225 | //
|
|---|
| 226 | // They will be also set as USED and will be taken into account during the
|
|---|
| 227 | // calculation of the image parameters.
|
|---|
| 228 | // The only way to distinguish them from the other USED pixels, is the
|
|---|
| 229 | // Ring number, that is bigger than 1.
|
|---|
| 230 | //
|
|---|
| 231 | // Example: You can decide how many rings you want to analyze using:
|
|---|
| 232 | //
|
|---|
| 233 | // MImgCleanStd clean;
|
|---|
| 234 | // //creates a default cleaning object (default number of rings =1)
|
|---|
| 235 | // clean.SetCleanRings(UShort_t r);
|
|---|
| 236 | // //now it looks r times around the CORE pixels to find new pixels with
|
|---|
| 237 | // //signal.
|
|---|
| 238 | //
|
|---|
| 239 | //
|
|---|
| 240 | // Input Containers:
|
|---|
| 241 | // MGeomCam, MCerPhotEvt, MSigmabar
|
|---|
| 242 | //
|
|---|
| 243 | // Output Containers:
|
|---|
| 244 | // MCerPhotEvt
|
|---|
| 245 | //
|
|---|
| 246 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 247 | #include "MImgCleanStd.h"
|
|---|
| 248 |
|
|---|
| 249 | #include <stdlib.h> // atof
|
|---|
| 250 | #include <fstream.h> // ofstream, SavePrimitive
|
|---|
| 251 |
|
|---|
| 252 | #include <TGFrame.h> // TGFrame
|
|---|
| 253 | #include <TGLabel.h> // TGLabel
|
|---|
| 254 | #include <TGTextEntry.h> // TGTextEntry
|
|---|
| 255 |
|
|---|
| 256 | #include "MLog.h"
|
|---|
| 257 | #include "MLogManip.h"
|
|---|
| 258 |
|
|---|
| 259 | #include "MParList.h"
|
|---|
| 260 | #include "MGeomPix.h"
|
|---|
| 261 | #include "MGeomCam.h"
|
|---|
| 262 | #include "MCerPhotPix.h"
|
|---|
| 263 | #include "MCerPhotEvt.h"
|
|---|
| 264 | #include "MSigmabar.h"
|
|---|
| 265 |
|
|---|
| 266 | #include "MGGroupFrame.h" // MGGroupFrame
|
|---|
| 267 |
|
|---|
| 268 | ClassImp(MImgCleanStd);
|
|---|
| 269 |
|
|---|
| 270 | enum {
|
|---|
| 271 | kImgCleanLvl1,
|
|---|
| 272 | kImgCleanLvl2
|
|---|
| 273 | };
|
|---|
| 274 |
|
|---|
| 275 | static const TString gsDefName = "MImgCleanStd";
|
|---|
| 276 | static const TString gsDefTitle = "Task to perform image cleaning";
|
|---|
| 277 |
|
|---|
| 278 | // --------------------------------------------------------------------------
|
|---|
| 279 | //
|
|---|
| 280 | // Default constructor. Here you can specify the cleaning method and levels.
|
|---|
| 281 | // If you don't specify them the 'common standard' values 3.0 and 2.5 (sigma
|
|---|
| 282 | // above mean) are used.
|
|---|
| 283 | // Here you can also specify how many rings around the core pixels you want
|
|---|
| 284 | // to analyze (with the fixed lvl2). The default value for "rings" is 1.
|
|---|
| 285 | //
|
|---|
| 286 | MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2,
|
|---|
| 287 | const char *name, const char *title)
|
|---|
| 288 | : fSgb(NULL), fCleaningMethod(kStandard), fCleanLvl1(lvl1),
|
|---|
| 289 | fCleanLvl2(lvl2), fCleanRings(1)
|
|---|
| 290 |
|
|---|
| 291 | {
|
|---|
| 292 | fName = name ? name : gsDefName.Data();
|
|---|
| 293 | fTitle = title ? title : gsDefTitle.Data();
|
|---|
| 294 |
|
|---|
| 295 | Print();
|
|---|
| 296 | }
|
|---|
| 297 |
|
|---|
| 298 | // --------------------------------------------------------------------------
|
|---|
| 299 | //
|
|---|
| 300 | // NT 28/04/2003: now the option to use the standard method or the
|
|---|
| 301 | // democratic method is implemented:
|
|---|
| 302 | //
|
|---|
| 303 | // KStandard: This method looks for all pixels with an entry (photons)
|
|---|
| 304 | // that is three times bigger than the noise of the pixel
|
|---|
| 305 | // (default: 3 sigma, clean level 1)
|
|---|
| 306 | // AM 18/11/2002: now cut levels are proportional to the pixel area.
|
|---|
| 307 | // In this way the cut corresponds to a fixed
|
|---|
| 308 | // phe-density (otherwise, it would bias the images).
|
|---|
| 309 | //
|
|---|
| 310 | // Returns the maximum Pixel Id (used for ispixused in CleanStep2)
|
|---|
| 311 | //
|
|---|
| 312 | Int_t MImgCleanStd::CleanStep1Std()
|
|---|
| 313 | {
|
|---|
| 314 | const Int_t entries = fEvt->GetNumPixels();
|
|---|
| 315 |
|
|---|
| 316 | Int_t max = entries;
|
|---|
| 317 |
|
|---|
| 318 | //
|
|---|
| 319 | // check the number of all pixels against the noise level and
|
|---|
| 320 | // set them to 'unused' state if necessary
|
|---|
| 321 | //
|
|---|
| 322 | for (Int_t i=0; i<entries; i++ )
|
|---|
| 323 | {
|
|---|
| 324 | MCerPhotPix &pix = (*fEvt)[i];
|
|---|
| 325 |
|
|---|
| 326 | const Int_t id = pix.GetPixId();
|
|---|
| 327 |
|
|---|
| 328 | const Float_t entry = pix.GetNumPhotons();
|
|---|
| 329 | const Float_t noise = pix.GetErrorPhot();
|
|---|
| 330 | const Double_t ratio = TMath::Sqrt(fCam->GetPixRatio(id));
|
|---|
| 331 |
|
|---|
| 332 | // COBB: '<=' to skip entry=noise=0
|
|---|
| 333 | if (entry * ratio <= fCleanLvl1 * noise)
|
|---|
| 334 | pix.SetPixelUnused();
|
|---|
| 335 |
|
|---|
| 336 | if (id>max)
|
|---|
| 337 | max = id;
|
|---|
| 338 | }
|
|---|
| 339 |
|
|---|
| 340 | return max;
|
|---|
| 341 | }
|
|---|
| 342 |
|
|---|
| 343 | // --------------------------------------------------------------------------
|
|---|
| 344 | //
|
|---|
| 345 | // NT 28/04/2003: now the option to use the standard method or the
|
|---|
| 346 | // democratic method is implemented:
|
|---|
| 347 | //
|
|---|
| 348 | // "KDemocratic": this method looks for all pixels with an entry (photons)
|
|---|
| 349 | // that is n times bigger than the noise of the mean of the
|
|---|
| 350 | // inner pixels (default: 3 sigmabar, clean level 1)
|
|---|
| 351 | //
|
|---|
| 352 | // Returns the maximum Pixel Id (used for ispixused in CleanStep2)
|
|---|
| 353 | //
|
|---|
| 354 | Int_t MImgCleanStd::CleanStep1Dem()
|
|---|
| 355 | {
|
|---|
| 356 | const Int_t entries = fEvt->GetNumPixels();
|
|---|
| 357 |
|
|---|
| 358 | Int_t max = entries;
|
|---|
| 359 |
|
|---|
| 360 | //
|
|---|
| 361 | // check the number of all pixels against the noise level and
|
|---|
| 362 | // set them to 'unused' state if necessary
|
|---|
| 363 | //
|
|---|
| 364 | for (Int_t i=0; i<entries; i++ )
|
|---|
| 365 | {
|
|---|
| 366 | MCerPhotPix &pix = (*fEvt)[i];
|
|---|
| 367 |
|
|---|
| 368 | const Int_t id = pix.GetPixId();
|
|---|
| 369 |
|
|---|
| 370 | const Float_t entry = pix.GetNumPhotons();
|
|---|
| 371 | const Double_t ratio = fCam->GetPixRatio(id);
|
|---|
| 372 |
|
|---|
| 373 | // COBB: '<=' to skip entry=noise=0
|
|---|
| 374 | if (entry * ratio <= fCleanLvl1 * fInnerNoise)
|
|---|
| 375 | pix.SetPixelUnused();
|
|---|
| 376 |
|
|---|
| 377 | if (id>max)
|
|---|
| 378 | max = id;
|
|---|
| 379 | }
|
|---|
| 380 | return max;
|
|---|
| 381 | }
|
|---|
| 382 |
|
|---|
| 383 | // --------------------------------------------------------------------------
|
|---|
| 384 | // The first step of cleaning defines the CORE pixels. All the other pixels
|
|---|
| 385 | // are set as UNUSED and belong to RING 0.
|
|---|
| 386 | // After this point, only the CORE pixels are set as USED, with RING
|
|---|
| 387 | // number 1.
|
|---|
| 388 | // Returns the maximum Pixel Id (used for ispixused in CleanStep2)
|
|---|
| 389 | //
|
|---|
| 390 | Int_t MImgCleanStd::CleanStep1()
|
|---|
| 391 | {
|
|---|
| 392 | switch (fCleaningMethod)
|
|---|
| 393 | {
|
|---|
| 394 | case kStandard:
|
|---|
| 395 | return CleanStep1Std();
|
|---|
| 396 | case kDemocratic:
|
|---|
| 397 | return CleanStep1Dem();
|
|---|
| 398 | }
|
|---|
| 399 |
|
|---|
| 400 | return 0;
|
|---|
| 401 | }
|
|---|
| 402 |
|
|---|
| 403 | // --------------------------------------------------------------------------
|
|---|
| 404 | //
|
|---|
| 405 | // Check if the survived pixel have a neighbor, that also
|
|---|
| 406 | // survived, otherwise set pixel to unused (removes pixels without
|
|---|
| 407 | // neighbors).
|
|---|
| 408 | //
|
|---|
| 409 | // Takes the maximum pixel id from CleanStep1 as an argument
|
|---|
| 410 | //
|
|---|
| 411 | void MImgCleanStd::CleanStep2(Int_t max)
|
|---|
| 412 | {
|
|---|
| 413 | const Int_t entries = fEvt->GetNumPixels();
|
|---|
| 414 |
|
|---|
| 415 | //
|
|---|
| 416 | // In the worst case we have to loop 6 times 577 times, to
|
|---|
| 417 | // catch the behaviour of all next neighbors. Here we can gain
|
|---|
| 418 | // much by using an array instead of checking through all pixels
|
|---|
| 419 | // (MCerPhotEvt::IsPixelUsed) all the time.
|
|---|
| 420 | //
|
|---|
| 421 | Byte_t *ispixused = new Byte_t[max+1];
|
|---|
| 422 |
|
|---|
| 423 | for (Int_t i=0; i<entries; i++)
|
|---|
| 424 | {
|
|---|
| 425 | const MCerPhotPix &pix = (*fEvt)[i];
|
|---|
| 426 | ispixused[pix.GetPixId()] = pix.IsPixelUsed() ? 1 : 0 ;
|
|---|
| 427 | }
|
|---|
| 428 |
|
|---|
| 429 | for (Int_t i=0; i<entries; i++)
|
|---|
| 430 | {
|
|---|
| 431 | // get entry i from list
|
|---|
| 432 | MCerPhotPix &pix = (*fEvt)[i];
|
|---|
| 433 |
|
|---|
| 434 | // get pixel id of this entry
|
|---|
| 435 | const Int_t id = pix.GetPixId();
|
|---|
| 436 |
|
|---|
| 437 | // check if pixel is in use, if not goto next pixel in list
|
|---|
| 438 | if (ispixused[id] == 0)
|
|---|
| 439 | continue;
|
|---|
| 440 |
|
|---|
| 441 | // check for 'used' neighbors of this pixel
|
|---|
| 442 | const MGeomPix &gpix = (*fCam)[id];
|
|---|
| 443 | const Int_t nnmax = gpix.GetNumNeighbors();
|
|---|
| 444 |
|
|---|
| 445 | Bool_t hasNeighbor = kFALSE;
|
|---|
| 446 |
|
|---|
| 447 | //loop on the neighbors to check if they are used
|
|---|
| 448 | for (Int_t j=0; j<nnmax; j++)
|
|---|
| 449 | {
|
|---|
| 450 | const Int_t id2 = gpix.GetNeighbor(j);
|
|---|
| 451 |
|
|---|
| 452 | // when you find an used neighbor, break the loop
|
|---|
| 453 | if (ispixused[id2] == 1)
|
|---|
| 454 | {
|
|---|
| 455 | hasNeighbor = kTRUE;
|
|---|
| 456 | break;
|
|---|
| 457 | }
|
|---|
| 458 | }
|
|---|
| 459 |
|
|---|
| 460 | if (hasNeighbor == kFALSE)
|
|---|
| 461 | pix.SetPixelUnused();
|
|---|
| 462 | }
|
|---|
| 463 |
|
|---|
| 464 | delete ispixused;
|
|---|
| 465 |
|
|---|
| 466 | //
|
|---|
| 467 | // now we declare all pixels that survive as CorePixels
|
|---|
| 468 | //
|
|---|
| 469 | for (Int_t i=0; i<entries; i++)
|
|---|
| 470 | {
|
|---|
| 471 | MCerPhotPix &pix = (*fEvt)[i];
|
|---|
| 472 |
|
|---|
| 473 | if (pix.IsPixelUsed())
|
|---|
| 474 | pix.SetPixelCore();
|
|---|
| 475 | }
|
|---|
| 476 | }
|
|---|
| 477 |
|
|---|
| 478 | // --------------------------------------------------------------------------
|
|---|
| 479 | //
|
|---|
| 480 | // Look for the boundary pixels around the core pixels
|
|---|
| 481 | // if a pixel has more than 2.5 (clean level 2.5) sigma, and
|
|---|
| 482 | // a core neigbor it is declared as used.
|
|---|
| 483 | //
|
|---|
| 484 | Bool_t MImgCleanStd::CleanStep3Std(const MCerPhotPix &pix)
|
|---|
| 485 | {
|
|---|
| 486 | //
|
|---|
| 487 | // get pixel id of this entry
|
|---|
| 488 | //
|
|---|
| 489 | const Int_t id = pix.GetPixId();
|
|---|
| 490 |
|
|---|
| 491 | //
|
|---|
| 492 | // check the num of photons against the noise level
|
|---|
| 493 | //
|
|---|
| 494 | const Float_t entry = pix.GetNumPhotons();
|
|---|
| 495 | const Float_t noise = pix.GetErrorPhot();
|
|---|
| 496 | const Double_t ratio = TMath::Sqrt(fCam->GetPixRatio(id));
|
|---|
| 497 |
|
|---|
| 498 | return (entry * ratio <= fCleanLvl2 * noise);
|
|---|
| 499 | }
|
|---|
| 500 |
|
|---|
| 501 | // --------------------------------------------------------------------------
|
|---|
| 502 | //
|
|---|
| 503 | // Look for the boundary pixels around the core pixels
|
|---|
| 504 | // if a pixel has more than 2.5 (clean level 2.5) sigmabar and
|
|---|
| 505 | // a core neighbor, it is declared as used.
|
|---|
| 506 | //
|
|---|
| 507 | Bool_t MImgCleanStd::CleanStep3Dem(const MCerPhotPix &pix)
|
|---|
| 508 | {
|
|---|
| 509 | //
|
|---|
| 510 | // get pixel id of this entry
|
|---|
| 511 | //
|
|---|
| 512 | const Int_t id = pix.GetPixId();
|
|---|
| 513 |
|
|---|
| 514 | //
|
|---|
| 515 | // check the num of photons against the noise level
|
|---|
| 516 | //
|
|---|
| 517 | const Float_t entry = pix.GetNumPhotons();
|
|---|
| 518 | const Double_t ratio = fCam->GetPixRatio(id);
|
|---|
| 519 |
|
|---|
| 520 | return (entry * ratio <= fCleanLvl2 * fInnerNoise);
|
|---|
| 521 | }
|
|---|
| 522 |
|
|---|
| 523 | void MImgCleanStd::CleanStep3b(MCerPhotPix &pix)
|
|---|
| 524 | {
|
|---|
| 525 | const Int_t id = pix.GetPixId();
|
|---|
| 526 |
|
|---|
| 527 | //
|
|---|
| 528 | // check if the pixel's next neighbor is a core pixel.
|
|---|
| 529 | // if it is a core pixel set pixel state to: used.
|
|---|
| 530 | //
|
|---|
| 531 | MGeomPix &gpix = (*fCam)[id];
|
|---|
| 532 | const Int_t nnmax = gpix.GetNumNeighbors();
|
|---|
| 533 |
|
|---|
| 534 | for (Int_t j=0; j<nnmax; j++)
|
|---|
| 535 | {
|
|---|
| 536 | const Int_t id2 = gpix.GetNeighbor(j);
|
|---|
| 537 |
|
|---|
| 538 | if (!fEvt->GetPixById(id2) || !fEvt->IsPixelCore(id2))
|
|---|
| 539 | continue;
|
|---|
| 540 |
|
|---|
| 541 | pix.SetPixelUsed();
|
|---|
| 542 | break;
|
|---|
| 543 | }
|
|---|
| 544 | }
|
|---|
| 545 |
|
|---|
| 546 | // --------------------------------------------------------------------------
|
|---|
| 547 | //
|
|---|
| 548 | // NT: Add option "rings": default value = 1.
|
|---|
| 549 | // Look n (n>1) times for the boundary pixels around the used pixels.
|
|---|
| 550 | // If a pixel has more than 2.5 (clean level 2.5) sigma,
|
|---|
| 551 | // it is declared as used.
|
|---|
| 552 | //
|
|---|
| 553 | // If a value<2 for fCleanRings is used, no CleanStep4 is done.
|
|---|
| 554 | //
|
|---|
| 555 | void MImgCleanStd::CleanStep4(UShort_t r, MCerPhotPix &pix)
|
|---|
| 556 | {
|
|---|
| 557 | //
|
|---|
| 558 | // check if the pixel's next neighbor is a used pixel.
|
|---|
| 559 | // if it is a used pixel set pixel state to: used,
|
|---|
| 560 | // and tell to which ring it belongs to.
|
|---|
| 561 | //
|
|---|
| 562 | const Int_t id = pix.GetPixId();
|
|---|
| 563 | MGeomPix &gpix = (*fCam)[id];
|
|---|
| 564 |
|
|---|
| 565 | const Int_t nnmax = gpix.GetNumNeighbors();
|
|---|
| 566 |
|
|---|
| 567 | for (Int_t j=0; j<nnmax; j++)
|
|---|
| 568 | {
|
|---|
| 569 | const Int_t id2 = gpix.GetNeighbor(j);
|
|---|
| 570 |
|
|---|
| 571 | MCerPhotPix &npix = *fEvt->GetPixById(id2);
|
|---|
| 572 |
|
|---|
| 573 | // FIXME!
|
|---|
| 574 | // Needed check to read CT1 data without having a Segmentation fault
|
|---|
| 575 | if (!fEvt->GetPixById(id2))
|
|---|
| 576 | continue;
|
|---|
| 577 |
|
|---|
| 578 | if (!npix.IsPixelUsed() || npix.GetRing()>r-1 )
|
|---|
| 579 | continue;
|
|---|
| 580 |
|
|---|
| 581 | pix.SetRing(r);
|
|---|
| 582 | break;
|
|---|
| 583 | }
|
|---|
| 584 | }
|
|---|
| 585 |
|
|---|
| 586 | // --------------------------------------------------------------------------
|
|---|
| 587 | //
|
|---|
| 588 | // Look for the boundary pixels around the core pixels
|
|---|
| 589 | // if a pixel has more than 2.5 (clean level 2.5) sigma, and
|
|---|
| 590 | // a core neigbor, it is declared as used.
|
|---|
| 591 | //
|
|---|
| 592 | void MImgCleanStd::CleanStep3()
|
|---|
| 593 | {
|
|---|
| 594 | const Int_t entries = fEvt->GetNumPixels();
|
|---|
| 595 |
|
|---|
| 596 | for (UShort_t r=1; r<fCleanRings+1; r++)
|
|---|
| 597 | {
|
|---|
| 598 | for (Int_t i=0; i<entries; i++)
|
|---|
| 599 | {
|
|---|
| 600 | //
|
|---|
| 601 | // get pixel as entry il from list
|
|---|
| 602 | //
|
|---|
| 603 | MCerPhotPix &pix = (*fEvt)[i];
|
|---|
| 604 |
|
|---|
| 605 | //
|
|---|
| 606 | // if pixel is a core pixel go to the next pixel
|
|---|
| 607 | //
|
|---|
| 608 | if (pix.IsPixelCore())
|
|---|
| 609 | continue;
|
|---|
| 610 |
|
|---|
| 611 | switch (fCleaningMethod)
|
|---|
| 612 | {
|
|---|
| 613 | case kStandard:
|
|---|
| 614 | if (CleanStep3Std(pix))
|
|---|
| 615 | continue;
|
|---|
| 616 | break;
|
|---|
| 617 | case kDemocratic:
|
|---|
| 618 | if (CleanStep3Dem(pix))
|
|---|
| 619 | continue;
|
|---|
| 620 | break;
|
|---|
| 621 | }
|
|---|
| 622 |
|
|---|
| 623 | if (r==1)
|
|---|
| 624 | CleanStep3b(pix);
|
|---|
| 625 | else
|
|---|
| 626 | CleanStep4(r, pix);
|
|---|
| 627 | }
|
|---|
| 628 | }
|
|---|
| 629 | }
|
|---|
| 630 |
|
|---|
| 631 | // --------------------------------------------------------------------------
|
|---|
| 632 | //
|
|---|
| 633 | // Check if MEvtHeader exists in the Parameter list already.
|
|---|
| 634 | // if not create one and add them to the list
|
|---|
| 635 | //
|
|---|
| 636 | Bool_t MImgCleanStd::PreProcess (MParList *pList)
|
|---|
| 637 | {
|
|---|
| 638 | fCam = (MGeomCam*)pList->FindObject("MGeomCam");
|
|---|
| 639 | if (!fCam)
|
|---|
| 640 | {
|
|---|
| 641 | *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
|
|---|
| 642 | return kFALSE;
|
|---|
| 643 | }
|
|---|
| 644 |
|
|---|
| 645 | fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
|
|---|
| 646 | if (!fEvt)
|
|---|
| 647 | {
|
|---|
| 648 | *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
|
|---|
| 649 | return kFALSE;
|
|---|
| 650 | }
|
|---|
| 651 |
|
|---|
| 652 | if (fCleaningMethod != kDemocratic)
|
|---|
| 653 | return kTRUE;
|
|---|
| 654 |
|
|---|
| 655 | fSgb = (MSigmabar*)pList->FindObject("MSigmabar");
|
|---|
| 656 | if (!fSgb)
|
|---|
| 657 | {
|
|---|
| 658 | *fLog << dbginf << "MSigmabar not found... aborting." << endl;
|
|---|
| 659 | return kFALSE;
|
|---|
| 660 | }
|
|---|
| 661 |
|
|---|
| 662 | return kTRUE;
|
|---|
| 663 | }
|
|---|
| 664 |
|
|---|
| 665 | // --------------------------------------------------------------------------
|
|---|
| 666 | //
|
|---|
| 667 | // Cleans the image.
|
|---|
| 668 | //
|
|---|
| 669 | Bool_t MImgCleanStd::Process()
|
|---|
| 670 | {
|
|---|
| 671 | if (fSgb)
|
|---|
| 672 | fInnerNoise = fSgb->GetSigmabarInner();
|
|---|
| 673 |
|
|---|
| 674 | const Int_t max = CleanStep1();
|
|---|
| 675 | CleanStep2(max);
|
|---|
| 676 | CleanStep3();
|
|---|
| 677 |
|
|---|
| 678 | return kTRUE;
|
|---|
| 679 | }
|
|---|
| 680 |
|
|---|
| 681 | // --------------------------------------------------------------------------
|
|---|
| 682 | //
|
|---|
| 683 | // Print descriptor and cleaning levels.
|
|---|
| 684 | //
|
|---|
| 685 | void MImgCleanStd::Print(Option_t *o) const
|
|---|
| 686 | {
|
|---|
| 687 | *fLog << all << GetDescriptor() << " using ";
|
|---|
| 688 | switch (fCleaningMethod)
|
|---|
| 689 | {
|
|---|
| 690 | case kDemocratic:
|
|---|
| 691 | *fLog << "democratic";
|
|---|
| 692 | break;
|
|---|
| 693 | case kStandard:
|
|---|
| 694 | *fLog << "standard";
|
|---|
| 695 | break;
|
|---|
| 696 | }
|
|---|
| 697 | *fLog << " cleaning initialized with noise level " << fCleanLvl1 << " and " << fCleanLvl2;
|
|---|
| 698 | *fLog << " (CleanRings=" << fCleanRings << ")" << endl;
|
|---|
| 699 | }
|
|---|
| 700 |
|
|---|
| 701 | // --------------------------------------------------------------------------
|
|---|
| 702 | //
|
|---|
| 703 | // Create two text entry fields, one for each cleaning level and a
|
|---|
| 704 | // describing text line.
|
|---|
| 705 | //
|
|---|
| 706 | void MImgCleanStd::CreateGuiElements(MGGroupFrame *f)
|
|---|
| 707 | {
|
|---|
| 708 | //
|
|---|
| 709 | // Create a frame for line 3 and 4 to be able
|
|---|
| 710 | // to align entry field and label in one line
|
|---|
| 711 | //
|
|---|
| 712 | TGHorizontalFrame *f1 = new TGHorizontalFrame(f, 0, 0);
|
|---|
| 713 | TGHorizontalFrame *f2 = new TGHorizontalFrame(f, 0, 0);
|
|---|
| 714 |
|
|---|
| 715 | /*
|
|---|
| 716 | * --> use with root >=3.02 <--
|
|---|
| 717 | *
|
|---|
| 718 |
|
|---|
| 719 | TGNumberEntry *fNumEntry1 = new TGNumberEntry(frame, 3.0, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
|
|---|
| 720 | TGNumberEntry *fNumEntry2 = new TGNumberEntry(frame, 2.5, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
|
|---|
| 721 |
|
|---|
| 722 | */
|
|---|
| 723 | TGTextEntry *entry1 = new TGTextEntry(f1, "****", kImgCleanLvl1);
|
|---|
| 724 | TGTextEntry *entry2 = new TGTextEntry(f2, "****", kImgCleanLvl2);
|
|---|
| 725 |
|
|---|
| 726 | // --- doesn't work like expected (until root 3.02?) --- fNumEntry1->SetAlignment(kTextRight);
|
|---|
| 727 | // --- doesn't work like expected (until root 3.02?) --- fNumEntry2->SetAlignment(kTextRight);
|
|---|
| 728 |
|
|---|
| 729 | entry1->SetText("3.0");
|
|---|
| 730 | entry2->SetText("2.5");
|
|---|
| 731 |
|
|---|
| 732 | entry1->Associate(f);
|
|---|
| 733 | entry2->Associate(f);
|
|---|
| 734 |
|
|---|
| 735 | TGLabel *l1 = new TGLabel(f1, "Cleaning Level 1");
|
|---|
| 736 | TGLabel *l2 = new TGLabel(f2, "Cleaning Level 2");
|
|---|
| 737 |
|
|---|
| 738 | l1->SetTextJustify(kTextLeft);
|
|---|
| 739 | l2->SetTextJustify(kTextLeft);
|
|---|
| 740 |
|
|---|
| 741 | //
|
|---|
| 742 | // Align the text of the label centered, left in the row
|
|---|
| 743 | // with a left padding of 10
|
|---|
| 744 | //
|
|---|
| 745 | TGLayoutHints *laylabel = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 10);
|
|---|
| 746 | TGLayoutHints *layframe = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 5, 0, 10);
|
|---|
| 747 |
|
|---|
| 748 | //
|
|---|
| 749 | // Add one entry field and the corresponding label to each line
|
|---|
| 750 | //
|
|---|
| 751 | f1->AddFrame(entry1);
|
|---|
| 752 | f2->AddFrame(entry2);
|
|---|
| 753 |
|
|---|
| 754 | f1->AddFrame(l1, laylabel);
|
|---|
| 755 | f2->AddFrame(l2, laylabel);
|
|---|
| 756 |
|
|---|
| 757 | f->AddFrame(f1, layframe);
|
|---|
| 758 | f->AddFrame(f2, layframe);
|
|---|
| 759 |
|
|---|
| 760 | f->AddToList(entry1);
|
|---|
| 761 | f->AddToList(entry2);
|
|---|
| 762 | f->AddToList(l1);
|
|---|
| 763 | f->AddToList(l2);
|
|---|
| 764 | f->AddToList(laylabel);
|
|---|
| 765 | f->AddToList(layframe);
|
|---|
| 766 | }
|
|---|
| 767 |
|
|---|
| 768 | // --------------------------------------------------------------------------
|
|---|
| 769 | //
|
|---|
| 770 | // Process the GUI Events comming from the two text entry fields.
|
|---|
| 771 | //
|
|---|
| 772 | Bool_t MImgCleanStd::ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2)
|
|---|
| 773 | {
|
|---|
| 774 | if (msg!=kC_TEXTENTRY || submsg!=kTE_ENTER)
|
|---|
| 775 | return kTRUE;
|
|---|
| 776 |
|
|---|
| 777 | TGTextEntry *txt = (TGTextEntry*)FindWidget(param1);
|
|---|
| 778 |
|
|---|
| 779 | if (!txt)
|
|---|
| 780 | return kTRUE;
|
|---|
| 781 |
|
|---|
| 782 | Float_t lvl = atof(txt->GetText());
|
|---|
| 783 |
|
|---|
| 784 | switch (param1)
|
|---|
| 785 | {
|
|---|
| 786 | case kImgCleanLvl1:
|
|---|
| 787 | fCleanLvl1 = lvl;
|
|---|
| 788 | *fLog << "Cleaning level 1 set to " << lvl << " sigma." << endl;
|
|---|
| 789 | return kTRUE;
|
|---|
| 790 |
|
|---|
| 791 | case kImgCleanLvl2:
|
|---|
| 792 | fCleanLvl2 = lvl;
|
|---|
| 793 | *fLog << "Cleaning level 2 set to " << lvl << " sigma." << endl;
|
|---|
| 794 | return kTRUE;
|
|---|
| 795 | }
|
|---|
| 796 |
|
|---|
| 797 | return kTRUE;
|
|---|
| 798 | }
|
|---|
| 799 |
|
|---|
| 800 | // --------------------------------------------------------------------------
|
|---|
| 801 | //
|
|---|
| 802 | // Implementation of SavePrimitive. Used to write the call to a constructor
|
|---|
| 803 | // to a macro. In the original root implementation it is used to write
|
|---|
| 804 | // gui elements to a macro-file.
|
|---|
| 805 | //
|
|---|
| 806 | void MImgCleanStd::StreamPrimitive(ofstream &out) const
|
|---|
| 807 | {
|
|---|
| 808 | out << " MImgCleanStd " << GetUniqueName() << "(";
|
|---|
| 809 | out << fCleanLvl1 << ", " << fCleanLvl2;
|
|---|
| 810 |
|
|---|
| 811 | if (fName!=gsDefName || fTitle!=gsDefTitle)
|
|---|
| 812 | {
|
|---|
| 813 | out << ", \"" << fName << "\"";
|
|---|
| 814 | if (fTitle!=gsDefTitle)
|
|---|
| 815 | out << ", \"" << fTitle << "\"";
|
|---|
| 816 | }
|
|---|
| 817 | out << ");" << endl;
|
|---|
| 818 |
|
|---|
| 819 | if (fCleaningMethod!=kDemocratic)
|
|---|
| 820 | return;
|
|---|
| 821 |
|
|---|
| 822 | out << " " << GetUniqueName() << ".SetMethod(MImgCleanStd::kDemocratic);" << endl;
|
|---|
| 823 |
|
|---|
| 824 | if (fCleanRings==1)
|
|---|
| 825 | return;
|
|---|
| 826 |
|
|---|
| 827 | out << " " << GetUniqueName() << ".SetCleanRings(" << fCleanRings << ");" << endl;
|
|---|
| 828 | }
|
|---|