source: releases/Mars.2014.05.26/mimage/MImgCleanSimple.cc@ 18481

Last change on this file since 18481 was 14202, checked in by tbretz, 12 years ago
New simple image cleaning.
File size: 12.9 KB
Line 
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
45ClassImp(MImgCleanSimple);
46
47using namespace std;
48
49static const TString gsDefName = "MImgCleanSimple";
50static const TString gsDefTitle = "Task to perform image cleaning";
51
52const TString MImgCleanSimple::gsNameSignalCam ="MSignalCam"; // default name of the 'MSignalCam' container
53const 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//
63MImgCleanSimple::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
72void 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
90void 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
135void 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
175void 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//
229Int_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//
253Int_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//
280void 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//
319void 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//
373Int_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}
Note: See TracBrowser for help on using the repository browser.