| 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, 06/2012 <mailto:thomas.bretz@epfl.ch>
|
|---|
| 19 | !
|
|---|
| 20 | ! Copyright: MAGIC Software Development, 2000-2012
|
|---|
| 21 | !
|
|---|
| 22 | !
|
|---|
| 23 | \* ======================================================================== */
|
|---|
| 24 |
|
|---|
| 25 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 26 | //
|
|---|
| 27 | // MImgCleanSimple
|
|---|
| 28 | //
|
|---|
| 29 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 30 | #include "MImgCleanSimple.h"
|
|---|
| 31 |
|
|---|
| 32 | #include <TEnv.h>
|
|---|
| 33 |
|
|---|
| 34 | #include "MLog.h"
|
|---|
| 35 | #include "MLogManip.h"
|
|---|
| 36 |
|
|---|
| 37 | #include "MParList.h"
|
|---|
| 38 |
|
|---|
| 39 | #include "MGeomPix.h"
|
|---|
| 40 | #include "MGeomCam.h"
|
|---|
| 41 |
|
|---|
| 42 | #include "MSignalPix.h"
|
|---|
| 43 | #include "MSignalCam.h"
|
|---|
| 44 |
|
|---|
| 45 | ClassImp(MImgCleanSimple);
|
|---|
| 46 |
|
|---|
| 47 | using namespace std;
|
|---|
| 48 |
|
|---|
| 49 | static const TString gsDefName = "MImgCleanSimple";
|
|---|
| 50 | static const TString gsDefTitle = "Task to perform image cleaning";
|
|---|
| 51 |
|
|---|
| 52 | const TString MImgCleanSimple::gsNameSignalCam ="MSignalCam"; // default name of the 'MSignalCam' container
|
|---|
| 53 | const TString MImgCleanSimple::gsNameGeomCam ="MGeomCam"; // default name of the 'MGeomCam' container
|
|---|
| 54 |
|
|---|
| 55 | // --------------------------------------------------------------------------
|
|---|
| 56 | //
|
|---|
| 57 | // Default constructor. Here you can specify the cleaning method and levels.
|
|---|
| 58 | // If you don't specify them the 'common standard' values 3.0 and 2.5 (sigma
|
|---|
| 59 | // above mean) are used.
|
|---|
| 60 | // Here you can also specify how many rings around the core pixels you want
|
|---|
| 61 | // to analyze (with the fixed lvl2). The default value for "rings" is 1.
|
|---|
| 62 | //
|
|---|
| 63 | MImgCleanSimple::MImgCleanSimple(const Float_t lvl1, const Float_t lvl2,
|
|---|
| 64 | const char *name, const char *title) : fCleaningMethod(kStandard),
|
|---|
| 65 | fCleanLvl1(lvl1), fCleanLvl2(-1), fTimeLvl1(lvl2), fTimeLvl2(-1),
|
|---|
| 66 | fNameGeomCam(gsNameGeomCam), fNameSignalCam(gsNameSignalCam)
|
|---|
| 67 | {
|
|---|
| 68 | fName = name ? name : gsDefName.Data();
|
|---|
| 69 | fTitle = title ? title : gsDefTitle.Data();
|
|---|
| 70 | }
|
|---|
| 71 |
|
|---|
| 72 | void MImgCleanSimple::ResetCleaning() const
|
|---|
| 73 | {
|
|---|
| 74 | //
|
|---|
| 75 | // check the number of all pixels against the noise level and
|
|---|
| 76 | // set them to 'unused' state if necessary
|
|---|
| 77 | //
|
|---|
| 78 | const UInt_t npixevt = fEvt->GetNumPixels();
|
|---|
| 79 | for (UInt_t idx=0; idx<npixevt; idx++)
|
|---|
| 80 | {
|
|---|
| 81 | MSignalPix &pix = (*fEvt)[idx];
|
|---|
| 82 | if (pix.IsPixelUnmapped())
|
|---|
| 83 | continue;
|
|---|
| 84 |
|
|---|
| 85 | pix.SetPixelUnused();
|
|---|
| 86 | pix.SetPixelCore(kFALSE);
|
|---|
| 87 | }
|
|---|
| 88 | }
|
|---|
| 89 |
|
|---|
| 90 | void MImgCleanSimple::DoCleanStd() const
|
|---|
| 91 | {
|
|---|
| 92 | const UInt_t npixevt = fEvt->GetNumPixels();
|
|---|
| 93 | for (UInt_t idx=0; idx<npixevt; idx++)
|
|---|
| 94 | {
|
|---|
| 95 | MSignalPix &pix = (*fEvt)[idx];
|
|---|
| 96 |
|
|---|
| 97 | // Check if this is above the cleaning level
|
|---|
| 98 | if (pix.GetNumPhotons() <= fCleanLvl1)
|
|---|
| 99 | continue;
|
|---|
| 100 |
|
|---|
| 101 | // Ignore unmapped pixels
|
|---|
| 102 | if (pix.IsPixelUnmapped())
|
|---|
| 103 | continue;
|
|---|
| 104 |
|
|---|
| 105 | const MGeom &gpix = (*fCam)[idx];
|
|---|
| 106 |
|
|---|
| 107 | const Int_t n = gpix.GetNumNeighbors();
|
|---|
| 108 | for (Int_t i=0; i<n; i++)
|
|---|
| 109 | {
|
|---|
| 110 | const MSignalPix &npix = (*fEvt)[gpix.GetNeighbor(i)];
|
|---|
| 111 |
|
|---|
| 112 | // Check if one neighbor pixel is also above the cleaning level
|
|---|
| 113 | if (npix.GetNumPhotons() <= fCleanLvl1)
|
|---|
| 114 | continue;
|
|---|
| 115 |
|
|---|
| 116 | // Check if both pixels are within the given time window
|
|---|
| 117 | if (TMath::Abs(npix.GetArrivalTime()-pix.GetArrivalTime())>=fTimeLvl1)
|
|---|
| 118 | continue;
|
|---|
| 119 |
|
|---|
| 120 | // Ignore unmapped pixels
|
|---|
| 121 | if (npix.IsPixelUnmapped())
|
|---|
| 122 | continue;
|
|---|
| 123 |
|
|---|
| 124 | // The pixel is now marked to survive the cleaning. It's neighbor will
|
|---|
| 125 | // be marked later, because the condition is symmetric and the neighbor
|
|---|
| 126 | // will fullfil the condition, too.
|
|---|
| 127 | pix.SetPixelUsed();
|
|---|
| 128 | pix.SetPixelCore();
|
|---|
| 129 |
|
|---|
| 130 | break;
|
|---|
| 131 | }
|
|---|
| 132 | }
|
|---|
| 133 | }
|
|---|
| 134 |
|
|---|
| 135 | void MImgCleanSimple::DoCleanSum2() const
|
|---|
| 136 | {
|
|---|
| 137 | const UInt_t npixevt = fEvt->GetNumPixels();
|
|---|
| 138 | for (UInt_t idx=0; idx<npixevt; idx++)
|
|---|
| 139 | {
|
|---|
| 140 | MSignalPix &pix = (*fEvt)[idx];
|
|---|
| 141 |
|
|---|
| 142 | // Ignore unmapped pixels
|
|---|
| 143 | if (pix.IsPixelUnmapped())
|
|---|
| 144 | continue;
|
|---|
| 145 |
|
|---|
| 146 | const MGeom &gpix = (*fCam)[idx];
|
|---|
| 147 |
|
|---|
| 148 | const Int_t n = gpix.GetNumNeighbors();
|
|---|
| 149 | for (Int_t i=0; i<n; i++)
|
|---|
| 150 | {
|
|---|
| 151 | const MSignalPix &npix = (*fEvt)[gpix.GetNeighbor(i)];
|
|---|
| 152 |
|
|---|
| 153 | // Check if the sum (or average) charge of the two neighbors is below threshold
|
|---|
| 154 | if (pix.GetNumPhotons()+npix.GetNumPhotons()<=fCleanLvl1*2)
|
|---|
| 155 | continue;
|
|---|
| 156 |
|
|---|
| 157 | // Check if both pixels are within the given time window
|
|---|
| 158 | if (TMath::Abs(npix.GetArrivalTime()-pix.GetArrivalTime())>=fTimeLvl1)
|
|---|
| 159 | continue;
|
|---|
| 160 |
|
|---|
| 161 | // Ignore unmapped pixels
|
|---|
| 162 | if (npix.IsPixelUnmapped())
|
|---|
| 163 | continue;
|
|---|
| 164 |
|
|---|
| 165 | // The pixel is now marked to survive the cleaning. It's neighbor will
|
|---|
| 166 | // be marked later, because the condition is symmetric and the neighbor
|
|---|
| 167 | // will fullfil the condition, too.
|
|---|
| 168 | pix.SetPixelUsed();
|
|---|
| 169 | pix.SetPixelCore();
|
|---|
| 170 | break;
|
|---|
| 171 | }
|
|---|
| 172 | }
|
|---|
| 173 | }
|
|---|
| 174 |
|
|---|
| 175 | void MImgCleanSimple::DoCleanSum3() const
|
|---|
| 176 | {
|
|---|
| 177 | const Double_t clvl = fCleanLvl2>0 ? fCleanLvl2*3 : fCleanLvl1*3;
|
|---|
| 178 | const Double_t tlvl = fTimeLvl2>0 ? fTimeLvl2 : fTimeLvl1;
|
|---|
| 179 |
|
|---|
| 180 | const UInt_t npixevt = fEvt->GetNumPixels();
|
|---|
| 181 | for (UInt_t idx=0; idx<npixevt; idx++)
|
|---|
| 182 | {
|
|---|
| 183 | MSignalPix &pix = (*fEvt)[idx];
|
|---|
| 184 |
|
|---|
| 185 | // Ignore unmapped pixels
|
|---|
| 186 | if (pix.IsPixelUnmapped())
|
|---|
| 187 | continue;
|
|---|
| 188 |
|
|---|
| 189 | const MGeom &gpix = (*fCam)[idx];
|
|---|
| 190 |
|
|---|
| 191 | const Int_t n = gpix.GetNumNeighbors();
|
|---|
| 192 | for (Int_t i=0; i<n; i++)
|
|---|
| 193 | {
|
|---|
| 194 | const MSignalPix &npix1 = (*fEvt)[gpix.GetNeighbor(i)];
|
|---|
| 195 | const MSignalPix &npix2 = (*fEvt)[gpix.GetNeighbor((i+1)%n)];
|
|---|
| 196 |
|
|---|
| 197 | // Check if the sum (or average) charge of the three neighbors (closed package)
|
|---|
| 198 | // is below threshold
|
|---|
| 199 | if (pix.GetNumPhotons()+npix1.GetNumPhotons()+npix2.GetNumPhotons()<=clvl)
|
|---|
| 200 | continue;
|
|---|
| 201 |
|
|---|
| 202 | // Check if all thre pixels fullfil the time window condition
|
|---|
| 203 | if (TMath::Abs(pix.GetArrivalTime()-npix1.GetArrivalTime())>=tlvl)
|
|---|
| 204 | continue;
|
|---|
| 205 | if (TMath::Abs(npix1.GetArrivalTime()-npix2.GetArrivalTime())>=tlvl)
|
|---|
| 206 | continue;
|
|---|
| 207 | if (TMath::Abs(npix2.GetArrivalTime()-pix.GetArrivalTime())>=tlvl)
|
|---|
| 208 | continue;
|
|---|
| 209 |
|
|---|
| 210 | // Ignore unmapped pixels
|
|---|
| 211 | if (npix1.IsPixelUnmapped() || npix2.IsPixelUnmapped())
|
|---|
| 212 | continue;
|
|---|
| 213 |
|
|---|
| 214 | // The pixel is now marked to survive the cleaning. It's neighbors will
|
|---|
| 215 | // be marked later, because the condition is symmetric and the neighbors
|
|---|
| 216 | // will fullfil the condition, too.
|
|---|
| 217 | pix.SetPixelUsed();
|
|---|
| 218 | pix.SetPixelCore();
|
|---|
| 219 | break;
|
|---|
| 220 | }
|
|---|
| 221 | }
|
|---|
| 222 | }
|
|---|
| 223 |
|
|---|
| 224 | // --------------------------------------------------------------------------
|
|---|
| 225 | //
|
|---|
| 226 | // Check if MEvtHeader exists in the Parameter list already.
|
|---|
| 227 | // if not create one and add them to the list
|
|---|
| 228 | //
|
|---|
| 229 | Int_t MImgCleanSimple::PreProcess(MParList *pList)
|
|---|
| 230 | {
|
|---|
| 231 | fCam = (MGeomCam*)pList->FindObject(AddSerialNumber(fNameGeomCam), "MGeomCam");
|
|---|
| 232 | if (!fCam)
|
|---|
| 233 | {
|
|---|
| 234 | *fLog << err << fNameGeomCam << " [MGeomCam] not found (no geometry information available)... aborting." << endl;
|
|---|
| 235 | return kFALSE;
|
|---|
| 236 | }
|
|---|
| 237 | fEvt = (MSignalCam*)pList->FindObject(AddSerialNumber(fNameSignalCam), "MSignalCam");
|
|---|
| 238 | if (!fEvt)
|
|---|
| 239 | {
|
|---|
| 240 | *fLog << err << fNameSignalCam << " [MSignalCam] not found... aborting." << endl;
|
|---|
| 241 | return kFALSE;
|
|---|
| 242 | }
|
|---|
| 243 |
|
|---|
| 244 | Print();
|
|---|
| 245 |
|
|---|
| 246 | return kTRUE;
|
|---|
| 247 | }
|
|---|
| 248 |
|
|---|
| 249 | // --------------------------------------------------------------------------
|
|---|
| 250 | //
|
|---|
| 251 | // Cleans the image.
|
|---|
| 252 | //
|
|---|
| 253 | Int_t MImgCleanSimple::Process()
|
|---|
| 254 | {
|
|---|
| 255 | ResetCleaning();
|
|---|
| 256 |
|
|---|
| 257 | switch (fCleaningMethod)
|
|---|
| 258 | {
|
|---|
| 259 | case kStandard:
|
|---|
| 260 | DoCleanStd(); // Check each of the neighbors charge individually
|
|---|
| 261 | break;
|
|---|
| 262 | case kSum2:
|
|---|
| 263 | DoCleanSum2(); // Check the sum of two neighbors
|
|---|
| 264 | case kSum3:
|
|---|
| 265 | DoCleanSum2(); // Check the sum of two neighbors
|
|---|
| 266 | DoCleanSum3(); // Check the sum of three neighbors
|
|---|
| 267 | break;
|
|---|
| 268 | }
|
|---|
| 269 |
|
|---|
| 270 | // Takes roughly 10% of the time
|
|---|
| 271 | fEvt->CalcIslands(*fCam);
|
|---|
| 272 |
|
|---|
| 273 | return kTRUE;
|
|---|
| 274 | }
|
|---|
| 275 |
|
|---|
| 276 | // --------------------------------------------------------------------------
|
|---|
| 277 | //
|
|---|
| 278 | // Print descriptor and cleaning levels.
|
|---|
| 279 | //
|
|---|
| 280 | void MImgCleanSimple::Print(Option_t *o) const
|
|---|
| 281 | {
|
|---|
| 282 | *fLog << all << GetDescriptor() << " using ";
|
|---|
| 283 | switch (fCleaningMethod)
|
|---|
| 284 | {
|
|---|
| 285 | case kStandard:
|
|---|
| 286 | *fLog << "standard";
|
|---|
| 287 | break;
|
|---|
| 288 | case kSum2:
|
|---|
| 289 | *fLog << "sum2";
|
|---|
| 290 | break;
|
|---|
| 291 | case kSum3:
|
|---|
| 292 | *fLog << "sum3";
|
|---|
| 293 | break;
|
|---|
| 294 |
|
|---|
| 295 | }
|
|---|
| 296 | /*
|
|---|
| 297 | *fLog << " cleaning" << endl;
|
|---|
| 298 | *fLog << "initialized with level " << fCleanLvl1 << " and " << fCleanLvl2;
|
|---|
| 299 | *fLog << " (CleanRings=" << fCleanRings << ")" << endl;
|
|---|
| 300 |
|
|---|
| 301 | if (fPostCleanType)
|
|---|
| 302 | {
|
|---|
| 303 | *fLog << "Time post cleaning is switched on with:" << endl;
|
|---|
| 304 | if (fPostCleanType&2)
|
|---|
| 305 | *fLog << " - Two pixel coincidence window: " << fTimeLvl2 << "ns" << endl;
|
|---|
| 306 | if (fPostCleanType&1)
|
|---|
| 307 | *fLog << " - One pixel coincidence window: " << fTimeLvl1 << "ns" << endl;
|
|---|
| 308 | }
|
|---|
| 309 | */
|
|---|
| 310 | }
|
|---|
| 311 |
|
|---|
| 312 | /*
|
|---|
| 313 | // --------------------------------------------------------------------------
|
|---|
| 314 | //
|
|---|
| 315 | // Implementation of SavePrimitive. Used to write the call to a constructor
|
|---|
| 316 | // to a macro. In the original root implementation it is used to write
|
|---|
| 317 | // gui elements to a macro-file.
|
|---|
| 318 | //
|
|---|
| 319 | void MImgCleanSimple::StreamPrimitive(ostream &out) const
|
|---|
| 320 | {
|
|---|
| 321 | out << " MImgCleanSimple " << GetUniqueName() << "(";
|
|---|
| 322 | out << fCleanLvl1 << ", " << fCleanLvl2;
|
|---|
| 323 |
|
|---|
| 324 | if (fName!=gsDefName || fTitle!=gsDefTitle)
|
|---|
| 325 | {
|
|---|
| 326 | out << ", \"" << fName << "\"";
|
|---|
| 327 | if (fTitle!=gsDefTitle)
|
|---|
| 328 | out << ", \"" << fTitle << "\"";
|
|---|
| 329 | }
|
|---|
| 330 | out << ");" << endl;
|
|---|
| 331 |
|
|---|
| 332 | if (fCleaningMethod!=kStandard)
|
|---|
| 333 | {
|
|---|
| 334 | out << " " << GetUniqueName() << ".SetMethod(MImgCleanSimple::k";
|
|---|
| 335 | switch (fCleaningMethod)
|
|---|
| 336 | {
|
|---|
| 337 | case kScaled: out << "Scaled"; break;
|
|---|
| 338 | case kDemocratic: out << "Democratic"; break;
|
|---|
| 339 | case kProbability: out << "Probability"; break;
|
|---|
| 340 | case kAbsolute: out << "Absolute"; break;
|
|---|
| 341 | case kTime : out << "Time"; break;
|
|---|
| 342 | default:
|
|---|
| 343 | break;
|
|---|
| 344 | }
|
|---|
| 345 | out << ");" << endl;
|
|---|
| 346 | }
|
|---|
| 347 | if (fCleanRings!=1)
|
|---|
| 348 | out << " " << GetUniqueName() << ".SetCleanRings(" << fCleanRings << ");" << endl;
|
|---|
| 349 |
|
|---|
| 350 | if (gsNamePedPhotCam!=fNamePedPhotCam)
|
|---|
| 351 | out << " " << GetUniqueName() << ".SetNamePedPhotCam(\"" << fNamePedPhotCam << "\");" << endl;
|
|---|
| 352 | if (gsNameGeomCam!=fNameGeomCam)
|
|---|
| 353 | out << " " << GetUniqueName() << ".SetNameGeomCam(\"" << fNameGeomCam << "\");" << endl;
|
|---|
| 354 | if (gsNameSignalCam!=fNameSignalCam)
|
|---|
| 355 | out << " " << GetUniqueName() << ".SetNameSignalCam(\"" << fNameSignalCam << "\");" << endl;
|
|---|
| 356 | if (fKeepIsolatedPixels)
|
|---|
| 357 | out << " " << GetUniqueName() << ".SetKeepIsolatedPixels();" << endl;
|
|---|
| 358 | if (fRecoverIsolatedPixels)
|
|---|
| 359 | out << " " << GetUniqueName() << ".SetRecoverIsolatedPixels(" << fRecoverIsolatedPixels << ");" << endl;
|
|---|
| 360 |
|
|---|
| 361 | }
|
|---|
| 362 | */
|
|---|
| 363 |
|
|---|
| 364 | // --------------------------------------------------------------------------
|
|---|
| 365 | //
|
|---|
| 366 | // Read the setup from a TEnv, eg:
|
|---|
| 367 | // MImgCleanSimple.CleanLevel1: 3.0
|
|---|
| 368 | // MImgCleanSimple.CleanLevel2: 2.5
|
|---|
| 369 | // MImgCleanSimple.CleanMethod: Standard, Scaled, Democratic, Probability, Absolute
|
|---|
| 370 | // MImgCleanSimple.CleanRings: 1
|
|---|
| 371 | // MImgCleanSimple.KeepIsolatedPixels: yes, no
|
|---|
| 372 | //
|
|---|
| 373 | Int_t MImgCleanSimple::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
|
|---|
| 374 | {
|
|---|
| 375 | Bool_t rc = kFALSE;
|
|---|
| 376 | if (IsEnvDefined(env, prefix, "CleanLevel1", print))
|
|---|
| 377 | {
|
|---|
| 378 | rc = kTRUE;
|
|---|
| 379 | fCleanLvl1 = GetEnvValue(env, prefix, "CleanLevel1", fCleanLvl1);
|
|---|
| 380 | }
|
|---|
| 381 | if (IsEnvDefined(env, prefix, "CleanLevel2", print))
|
|---|
| 382 | {
|
|---|
| 383 | rc = kTRUE;
|
|---|
| 384 | fCleanLvl2 = GetEnvValue(env, prefix, "CleanLevel2", fCleanLvl2);
|
|---|
| 385 | }
|
|---|
| 386 | if (IsEnvDefined(env, prefix, "TimeLevel1", print))
|
|---|
| 387 | {
|
|---|
| 388 | rc = kTRUE;
|
|---|
| 389 | fTimeLvl1 = GetEnvValue(env, prefix, "TimeLevel1", fTimeLvl1);
|
|---|
| 390 | }
|
|---|
| 391 | if (IsEnvDefined(env, prefix, "TimeLevel2", print))
|
|---|
| 392 | {
|
|---|
| 393 | rc = kTRUE;
|
|---|
| 394 | fTimeLvl2 = GetEnvValue(env, prefix, "TimeLevel2", fTimeLvl2);
|
|---|
| 395 | }
|
|---|
| 396 |
|
|---|
| 397 | if (IsEnvDefined(env, prefix, "CleanMethod", print))
|
|---|
| 398 | {
|
|---|
| 399 | rc = kTRUE;
|
|---|
| 400 | TString s = GetEnvValue(env, prefix, "CleanMethod", "");
|
|---|
| 401 | s.ToLower();
|
|---|
| 402 | if (s.BeginsWith("standard"))
|
|---|
| 403 | SetMethod(kStandard);
|
|---|
| 404 | if (s.BeginsWith("sum2"))
|
|---|
| 405 | SetMethod(kSum2);
|
|---|
| 406 | if (s.BeginsWith("sum3"))
|
|---|
| 407 | SetMethod(kSum3);
|
|---|
| 408 | }
|
|---|
| 409 |
|
|---|
| 410 | return rc;
|
|---|
| 411 | }
|
|---|
