source: trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc@ 9029

Last change on this file since 9029 was 8646, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 36.7 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, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Harald Kornmayer, 1/2001
20! Author(s): Nadia Tonello, 4/2003 <mailto:tonello@mppmu.mpg.de>
21! Author(s): Stefan Ruegamer, 03/2006 <mailto:snruegam@astro.uni-wuerzburg.de>
22!
23! Copyright: MAGIC Software Development, 2000-2007
24!
25!
26\* ======================================================================== */
27
28/////////////////////////////////////////////////////////////////////////////
29//
30// MImgCleanStd
31//
32// The Image Cleaning task selects the pixels you use for the Hillas
33// parameters calculation.
34//
35// There are two methods to make the selection: the standard one, as done
36// in the analysis of CT1 data, and the democratic one, as suggested by
37// W.Wittek. The number of photo-electrons of a pixel is compared with the
38// pedestal RMS of the pixel itself (standard method) or with the average
39// RMS of the inner pixels (democratic method).
40// In both cases, the possibility to have a camera with pixels of
41// different area is taken into account.
42// The too noisy pixels can be recognized and eventually switched off
43// (Unmap: set blind pixels to UNUSED) separately, using the
44// MBlindPixelCalc Class. In the MBlindPixelCalc class there is also the
45// function to replace the value of the noisy pixels with the interpolation
46// of the content of the neighbors (SetUseInterpolation).
47//
48// Example:
49// ...
50// MBlindPixelCalc blind;
51// blind.SetUseInterpolation();
52// blind.SetUseBlindPixels();
53//
54// MImgCleanStd clean;
55// ...
56// tlist.AddToList(&blind);
57// tlist.AddToList(&clean);
58//
59// Look at the MBlindPixelCalc Class for more details.
60//
61// Starting point: default values ----------------------------------------
62//
63// When an event is read, before the image cleaning, all the pixels that
64// are in MSignalCam are set as USED and NOT CORE. All the pixels belong
65// to RING number 1 (like USED pixels).
66// Look at MSignalPix.h to see how these informations of the pixel are
67// stored.
68// The default cleaning METHOD is the STANDARD one and the number of the
69// rings around the CORE pixel it analyzes is 1. Look at the Constructor
70// of the class in MImgCleanStd.cc to see (or change) the default values.
71//
72// Example: To modify this setting, use the member functions
73// SetMethod(MImgCleanStd::kDemocratic) and SetCleanRings(UShort_t n).
74//
75// MImgCleanStd:CleanStep1 -----------------------------------------------
76//
77// The first step of cleaning defines the CORE pixels. The CORE pixels are
78// the ones which contain the informations about the core of the electro-
79// magnetic shower.
80// The ratio (A_0/A_i) is calculated from fCam->GetPixRatio(i). A_0 is
81// the area of the central pixel of the camera, A_i is the area of the
82// examined pixel. In this way, if we have a MAGIC-like camera, with the
83// outer pixels bigger than the inner ones, the level of cleaning in the
84// two different regions is weighted.
85// This avoids problems of deformations of the shower images.
86// The signal S_i and the pedestal RMS Prms_i of the pixel are called from
87// the object MSignalPix.
88// If (default method = kStandard)
89//Begin_Html
90//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<img src="images/MImgCleanStd-f1.png">
91//End_Html
92// the pixel is set as CORE pixel. L_1 (n=1) is called "first level of
93// cleaning" (default: fCleanLvl1 = 3).
94// All the other pixels are set as UNUSED and belong to RING 0.
95// After this point, only the CORE pixels are set as USED, with RING
96// number 1.
97//
98// MImgCleanStd:CleanStep2 ----------------------------------------------
99//
100// The second step of cleaning looks at the isolated CORE pixels and sets
101// them to UNUSED. An isolated pixel is a pixel without CORE neighbors.
102// At the end of this point, we have set as USED only CORE pixels with at
103// least one CORE neighbor.
104//
105// MImgCleanStd:CleanStep3 ----------------------------------------------
106//
107// The third step of cleaning looks at all the pixels (USED or UNUSED) that
108// surround the USED pixels.
109// If the content of the analyzed pixel survives at the second level of
110// cleaning, i.e. if
111//Begin_Html
112//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<img src="images/MImgCleanStd-f1.png">
113//End_Html
114// the pixel is set as USED. L_2 (n=2) is called "second level of cleaning"
115// (default:fCleanLvl2 = 2.5).
116//
117// When the number of RINGS to analyze is 1 (default value), only the
118// pixels that have a neighbor CORE pixel are analyzed.
119//
120// There is the option to decide the number of times you want to repeat
121// this procedure (number of RINGS analyzed around the core pixels = n).
122// Every time the level of cleaning is the same (fCleanLvl2) and the pixel
123// will belong to ring r+1, 1 < r < n+1. This is described in
124// MImgCleanStd:CleanStep4 .
125//
126// Dictionary and member functions ---------------------------------------
127//
128// Here there is the detailed description of the member functions and of
129// the terms commonly used in the class.
130//
131//
132// STANDARD CLEANING:
133// =================
134// This is the method used for the CT1 data analysis. It is the default
135// method of the class.
136// The number of photo-electrons of a pixel (S_i) is compared to the
137// pedestal RMS of the pixel itself (Prms_i). To have the comparison to
138// the same photon density for all the pixels, taking into account they
139// can have different areas, we have to keep in mind that the number of
140// photons that hit each pixel, goes linearly with the area of the pixel.
141// The fluctuations of the LONS are proportional to sqrt(A_i), so when we
142// compare S_i with Prms_i, only a factor sqrt(A_0/A_i) is missing to
143// have the same (N.photons/Area) threshold for all the pixels.
144//
145// !!WARNING: if noise independent from the
146// pixel size (example: electronic noise) is introduced,
147// then the noise fluctuations are no longer proportional
148// to sqrt(A_i), and then the cut value (for a camera with
149// pixels of different sizes) resulting from the above
150// procedure would not be proportional to pixel size as we
151// intend. In that case, democratic cleaning is preferred.
152//
153// If
154//Begin_Html
155//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<img src="images/MImgCleanStd-f1.png">
156//End_Html
157// the pixel survives the cleaning and it is set as CORE (when L_n is the
158// first level of cleaning, fCleanLvl1) or USED (when L_n is the second
159// level of cleaning, fCleanLvl2).
160//
161// Example:
162//
163// MImgCleanStd clean;
164// //creates a default Cleaning object, with default setting
165// ...
166// tlist.AddToList(&clean);
167// // add the image cleaning to the main task list
168//
169//
170// DEMOCRATIC CLEANING:
171// ===================
172// You use this cleaning method when you want to compare the number of
173// photo-electrons of each pixel with the average pedestal RMS of the
174// inner pixels (for the MAGIC camera they are the smaller ones):
175//Begin_Html
176//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<img src="images/MImgCleanStd-f2.png">
177//End_Html
178// In this case, the simple ratio (A_0/A_i) is used to weight the level of
179// cleaning, because both the inner and the outer pixels (that in MAGIC
180// have a different area) are compared to the same pedestal RMS, coming
181// from the inner pixels.
182//
183// Make sure that you used a class calculating the MPedPhotCam which also
184// updated the contents of the mean values (Recalc) correctly.
185//
186//
187// PROBABILITY CLEANING
188// ====================
189// This method takes signal height (over signal noise) and arrival time
190// into account. Instead of comparing signal/Noise with cleaning level
191// one and two, we calculate
192// - P_ped: The probability that a signal is a pedestal (using
193// the signal height and the pedestal) For this probability the
194// same algorithm like in kScaled is used (which is a standard
195// cleaning which scales the noise with the mean noise of pixels
196// with the same size)
197// - P_sig: The probability that the signal corresponds to the pixel
198// with the highest signal. For this probability we use the
199// arrival time only.
200//
201// The cleaning now is done in levels of Probability (eg. 0.2, 0.05)
202// The meaning of the cleaning levels is essentially the same (the same cleaning
203// algorithm is used) but the cleaning is not done in levels of signal/noise
204// but in level of this probability.
205//
206// This probability is calculated as (1-P_ped)*P_sig
207//
208// Example:
209//
210// MImgCleanStd clean(0.2, 0.05);
211// clean.SetMethod(MImgCleanStd::kProbability);
212//
213//
214// ABSOLUTE CLEANING
215// =================
216// This method takes signal height (photons) times area ratio
217// ad the cleaning levels.
218//
219// The cleaning now is done in these levels (eg. 16, 20)
220// The meaning of the cleaning levels is essentially the same (the same cleaning
221// algorithm is used) but the cleaning is not done in different 'units'
222//
223// Example:
224//
225// MImgCleanStd clean(20, 16);
226// clean.SetMethod(MImgCleanStd::kAbsolulte);
227//
228//
229// Member Function: SetMethod()
230// ============================
231// When you call the MImgCleanStd task, the default method is kStandard.
232//
233// If you want to switch to the kDemocratic method you have to
234// call this member function.
235//
236// Example:
237//
238// MImgCleanStd clean;
239// //creates a default Cleaning object, with default setting
240//
241// clean.SetMethod(MImgCleanStd::kDemocratic);
242// //now the method of cleaning is changed to Democratic
243//
244//
245// FIRST AND SECOND CLEANING LEVEL
246// ===============================
247// When you call the MImgCleanStd task, the default cleaning levels are
248// fCleanLvl1 = 3, fCleanLvl2 = 2.5. You can change them easily when you
249// create the MImgCleanStd object.
250//
251// Example:
252//
253// MImgCleanStd clean(Float_t lvl1,Float_t lvl2);
254// //creates a default cleaning object, but the cleaning levels are now
255// //lvl1 and lvl2.
256//
257// RING NUMBER
258// ===========
259// The standard cleaning procedure is such that it looks for the
260// informations of the boundary part of the shower only on the first
261// neighbors of the CORE pixels.
262// There is the possibility now to look not only at the first neighbors
263// (first ring),but also further away, around the CORE pixels. All the new
264// pixels you can find with this method, are tested with the second level
265// of cleaning and have to have at least an USED neighbor.
266//
267// They will be also set as USED and will be taken into account during the
268// calculation of the image parameters.
269// The only way to distinguish them from the other USED pixels, is the
270// Ring number, that is bigger than 1.
271//
272// Example: You can decide how many rings you want to analyze using:
273//
274// MImgCleanStd clean;
275// //creates a default cleaning object (default number of rings =1)
276// clean.SetCleanRings(UShort_t r);
277// //now it looks r times around the CORE pixels to find new pixels with
278// //signal.
279//
280//
281// Class Version 4:
282// ----------------
283// + Float_t fTimeLvl2;
284// - Bool_t fKeepSinglePixels;
285// + Bool_t fKeepIsolatedPixels;
286// + Int_t fRecoverIsolatedPixels;
287//
288//
289// Input Containers:
290// MGeomCam
291// MPedPhotCam
292// MSignalCam
293//
294// Output Containers:
295// MSignalCam
296//
297/////////////////////////////////////////////////////////////////////////////
298#include "MImgCleanStd.h"
299
300#include <stdlib.h> // atof
301#include <fstream> // ofstream, SavePrimitive
302
303#include <TEnv.h>
304
305#include <TGFrame.h> // TGFrame
306#include <TGLabel.h> // TGLabel
307#include <TGTextEntry.h> // TGTextEntry
308
309#include "MLog.h"
310#include "MLogManip.h"
311
312#include "MArrayI.h"
313#include "MParList.h"
314#include "MCameraData.h"
315
316#include "MGeomPix.h"
317#include "MGeomCam.h"
318
319#include "MSignalPix.h"
320#include "MSignalCam.h"
321
322#include "MGGroupFrame.h" // MGGroupFrame
323
324ClassImp(MImgCleanStd);
325
326using namespace std;
327
328enum {
329 kImgCleanLvl1,
330 kImgCleanLvl2
331};
332
333static const TString gsDefName = "MImgCleanStd";
334static const TString gsDefTitle = "Task to perform image cleaning";
335
336const TString MImgCleanStd::gsNamePedPhotCam="MPedPhotCam"; // default name of the 'MPedPhotCam' container
337const TString MImgCleanStd::gsNameSignalCam ="MSignalCam"; // default name of the 'MSignalCam' container
338const TString MImgCleanStd::gsNameGeomCam ="MGeomCam"; // default name of the 'MGeomCam' container
339
340// --------------------------------------------------------------------------
341//
342// Default constructor. Here you can specify the cleaning method and levels.
343// If you don't specify them the 'common standard' values 3.0 and 2.5 (sigma
344// above mean) are used.
345// Here you can also specify how many rings around the core pixels you want
346// to analyze (with the fixed lvl2). The default value for "rings" is 1.
347//
348MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2,
349 const char *name, const char *title)
350 : fCleaningMethod(kStandard), fCleanLvl0(lvl1), fCleanLvl1(lvl1),
351 fCleanLvl2(lvl2), fTimeLvl1(1.5), fTimeLvl2(1.5), fCleanRings(1),
352 fKeepIsolatedPixels(kFALSE), fRecoverIsolatedPixels(0),
353 fPostCleanType(0), fNamePedPhotCam(gsNamePedPhotCam),
354 fNameGeomCam(gsNameGeomCam), fNameSignalCam(gsNameSignalCam)
355{
356 fName = name ? name : gsDefName.Data();
357 fTitle = title ? title : gsDefTitle.Data();
358}
359
360void MImgCleanStd::ResetCleaning() const
361{
362 //
363 // check the number of all pixels against the noise level and
364 // set them to 'unused' state if necessary
365 //
366 const UInt_t npixevt = fEvt->GetNumPixels();
367 for (UInt_t idx=0; idx<npixevt; idx++)
368 {
369 MSignalPix &pix = (*fEvt)[idx];
370 if (pix.IsPixelUnmapped())
371 continue;
372
373 pix.SetPixelUnused();
374 pix.SetPixelCore(kFALSE);
375 }
376}
377
378// --------------------------------------------------------------------------
379//
380// The first step of cleaning defines the CORE pixels. All the other pixels
381// are set as UNUSED and belong to RING 0.
382// After this point, only the CORE pixels are set as USED, with RING
383// number 1.
384//
385// NT 28/04/2003: now the option to use the standard method or the
386// democratic method is implemented:
387//
388// kStandard: This method looks for all pixels with an entry (photons)
389// that is three times bigger than the noise of the pixel
390// (default: 3 sigma, clean level 1)
391//
392// kDemocratic: this method looks for all pixels with an entry (photons)
393// that is n times bigger than the noise of the mean of the
394// inner pixels (default: 3 sigmabar, clean level 1)
395//
396//
397Bool_t MImgCleanStd::HasCoreNeighbors(const MGeomPix &gpix) const
398{
399// if (fKeepIsolatedPixels)
400// return kTRUE;
401
402#ifdef DEBUG
403 const TArrayD &data = fData->GetData();
404#else
405 const Double_t *data = fData->GetData().GetArray();
406#endif
407
408 //loop on the neighbors to check if they are used
409 const Int_t n = gpix.GetNumNeighbors();
410 for (Int_t i=0; i<n; i++)
411 {
412 const Int_t idx = gpix.GetNeighbor(i);
413
414 // Check if a neighborpixel of our core pixel is
415 // also a core pixel
416 if (data[idx]<=fCleanLvl1)
417 continue;
418
419 // Ignore unmapped pixels
420 MSignalPix &pix = (*fEvt)[idx];
421 if (pix.IsPixelUnmapped())
422 continue;
423
424 return kTRUE;
425 }
426
427 return kFALSE;
428}
429
430Bool_t MImgCleanStd::HasUsedNeighbors(const MGeomPix &gpix) const
431{
432 //loop on the neighbors to check if they are used
433 const Int_t n = gpix.GetNumNeighbors();
434 for (Int_t i=0; i<n; i++)
435 {
436 const Int_t idx = gpix.GetNeighbor(i);
437
438 MSignalPix &pix = (*fEvt)[idx];
439
440 // Check if a neighborpixel of our core pixel is
441 // also a core pixel
442 if (pix.IsPixelUsed() && !pix.IsPixelUnmapped())
443 return kTRUE;
444 }
445
446 return kFALSE;
447}
448
449
450void MImgCleanStd::SetUsedNeighbors(const MGeomPix &gpix, Int_t r) const
451{
452 if (r>fCleanRings)
453 return;
454
455#ifdef DEBUG
456 const TArrayD &data = fData->GetData();
457#else
458 const Double_t *data = fData->GetData().GetArray();
459#endif
460
461 // At least one neighbor has been identified as core,
462 // that means we will keep the pixel
463 const Int_t n = gpix.GetNumNeighbors();
464 for (Int_t i=0; i<n; i++)
465 {
466 const Int_t idx = gpix.GetNeighbor(i);
467
468 MSignalPix &pix = (*fEvt)[idx];
469
470 // If the pixel has been assigned to the same or a previous
471 // ring we don't have to proceed. We have to try to set the
472 // ring number of each pixel as low as possible. This means
473 // we are only allowed to increase the ring number.
474 if (pix.IsPixelUsed() && pix.GetRing()<=r)
475 continue;
476
477 // All pixels below the second cleaning level should be ignored
478 if (data[idx] <= fCleanLvl2)
479 continue;
480
481 // Ignore unmapped pixels (remark: used (aka. ring>0)
482 // and unmapped status is exclusive
483 if (pix.IsPixelUnmapped())
484 continue;
485
486 // Set or reset the ring number
487 pix.SetRing(r);
488
489 // Step forward to the surrounding pixels
490 SetUsedNeighbors((*fCam)[idx], r+1);
491 }
492}
493
494// --------------------------------------------------------------------------
495//
496// Here we do the cleaning. We search for all the possible core candidates
497// and from them on we recursively search for used pixels with
498// SetUsedNeighbors. To check the validity of a core pixel
499// either fTimeLvl2 and/or HasCoreNeighbors is used.
500// The size of all removed
501Int_t MImgCleanStd::DoCleaning(Float_t &size) const
502{
503 Int_t n = 0;
504 size = 0;
505
506#ifdef DEBUG
507 const TArrayD &data = fData->GetData();
508#else
509 const Double_t *data = fData->GetData().GetArray();
510#endif
511
512 //
513 // check the number of all pixels against the noise level and
514 // set them to 'unused' state if necessary
515 //
516 const UInt_t npixevt = fEvt->GetNumPixels();
517 for (UInt_t idx=0; idx<npixevt; idx++)
518 {
519 MSignalPix &pix = (*fEvt)[idx];
520
521 // Check if this pixel is a possible candidate for a core pixel
522 if (data[idx] <= fCleanLvl1)
523 continue;
524
525 // Ignore unmapped pixels
526 if (pix.IsPixelUnmapped())
527 continue;
528
529 const MGeomPix &gpix = (*fCam)[idx];
530
531 // Check if the pixel is an isolated core pixel
532 if (!HasCoreNeighbors(gpix))
533 {
534 // Count size and number of isolated core pixels
535 size += pix.GetNumPhotons();
536 n++;
537
538 // If isolated pixels should not be kept or the pixel
539 // is lower than the cleaning level for isolated core
540 // pixels. It is not treated as core pixel.
541 if (!fKeepIsolatedPixels || data[idx]<=fCleanLvl0)
542 continue;
543 }
544
545 // Mark pixel as used and core
546 pix.SetPixelUsed();
547 pix.SetPixelCore();
548
549 // Check if neighbor pixels should be marked as used
550 // This is done recursively depening on fCleanRings
551 SetUsedNeighbors(gpix);
552 }
553
554 return n;
555}
556
557/*
558Float_t MImgCleanStd::GetArrivalTimeNeighbor(const MGeomPix &gpix) const
559{
560 Float_t min = FLT_MAX;
561
562 const Int_t n = gpix.GetNumNeighbors();
563 for (int i=0; i<n; i++)
564 {
565 const Int_t idx = gpix.GetNeighbor(i);
566
567 const MSignalPix &pix = (*fEvt)[idx];
568 // FIXME: Check also used pixels?
569 if (!pix.IsCorePixel())
570 continue;
571
572 const Float_t tm = pix.GetArrivalTime();
573 if (tm<min)
574 min = tm;
575 }
576
577 return tm;
578}
579
580void MImgCleanStd::CheckUsedPixelsForArrivalTime(Float_t timediff) const
581{
582 const MArrayD &data = fData->GetData();
583
584 //
585 // check the number of all pixels against the noise level and
586 // set them to 'unused' state if necessary
587 //
588 const UInt_t npixevt = fEvt->GetNumPixels();
589 for (UInt_t idx=0; idx<npixevt; idx++)
590 {
591 MSignalPix &pix = (*fEvt)[idx];
592
593 // If pixel has previously been marked used, ignore
594 if (!pix.IsPixelUsed() || pix.IsPixelCore())
595 continue;
596
597 const MGeomPix &gpix = (*fCam)[idx];
598
599 // If isolated possible-corepixel doesn't have used
600 // neighbors, ignore it
601 const Float_t tm0 = pix.GetArrivalTime();
602 const Float_t tm1 = GetArrivalTimeCoreNeighbor(gpix);
603
604 if (TMath::Abs(tm0-tm1)<timediff)
605 continue;
606
607 // Mark pixel as used and core
608 pix.SetPixelUnused();
609 }
610}
611*/
612
613Int_t MImgCleanStd::RecoverIsolatedPixels(Float_t &size) const
614{
615#ifdef DEBUG
616 const TArrayD &data = fData->GetData();
617#else
618 const Double_t *data = fData->GetData().GetArray();
619#endif
620
621 Int_t n = 0;
622
623 //
624 // check the number of all pixels against the noise level and
625 // set them to 'unused' state if necessary
626 //
627 const UInt_t npixevt = fEvt->GetNumPixels();
628 for (UInt_t idx=0; idx<npixevt; idx++)
629 {
630 MSignalPix &pix = (*fEvt)[idx];
631
632 // If pixel has previously been marked used, ignore
633 if (pix.IsPixelUsed())
634 continue;
635
636 // If pixel is not a candidate for a core pixel, ignore
637 if (data[idx] <= fCleanLvl1)
638 continue;
639
640 const MGeomPix &gpix = (*fCam)[idx];
641
642 // If isolated possible-corepixel doesn't have used
643 // neighbors, ignore it
644 if (!HasUsedNeighbors(gpix))
645 continue;
646
647 // Mark pixel as used and core
648 pix.SetPixelUsed();
649 pix.SetPixelCore();
650
651 // Check if neighbor pixels should be marked as used
652 // This is done recursively depening on fCleanRings
653 SetUsedNeighbors(gpix);
654
655 size -= pix.GetNumPhotons();
656 n++;
657 }
658
659 return n;
660}
661
662void MImgCleanStd::CleanTime(Int_t n, Double_t lvl) const
663{
664 MArrayI indices;
665
666 const UInt_t npixevt = fEvt->GetNumPixels();
667 for (UInt_t idx=0; idx<npixevt; idx++)
668 {
669 // check if pixel is used or not
670 const MSignalPix &pix = (*fEvt)[idx];
671 if (!pix.IsPixelUsed())
672 continue;
673
674 // get arrival time
675 const Double_t tm0 = pix.GetArrivalTime();
676
677 // loop over its neighbpors
678 const MGeomPix &gpix = (*fCam)[idx];
679
680 Int_t cnt = 0;
681 for (Int_t i=0; i<gpix.GetNumNeighbors(); i++)
682 {
683 // Get index of neighbor
684 const Int_t idx2 = gpix.GetNeighbor(i);
685
686 // check if neighbor is used or not
687 const MSignalPix &npix = (*fEvt)[idx2];
688 if (!npix.IsPixelUsed())
689 continue;
690
691 // If this pixel is to far away (in arrival time) don't count
692 if (TMath::Abs(npix.GetArrivalTime()-tm0)>lvl)
693 continue;
694
695 // Now count the pixel. If we did not found n pixels yet
696 // which fullfill the condition, go on searching
697 if (++cnt>=n)
698 break;
699 }
700
701 // If we found at least n neighbors which are
702 // with a time difference of lvl keep the pixel
703 if (cnt>=n)
704 continue;
705
706 indices.Set(indices.GetSize()+1);
707 indices[indices.GetSize()-1] = idx;
708 }
709
710 // Now remove the pixels which didn't fullfill the requirement
711 for (UInt_t i=0; i<indices.GetSize(); i++)
712 {
713 (*fEvt)[indices[i]].SetPixelUnused();
714 (*fEvt)[indices[i]].SetPixelCore(kFALSE);
715 }
716}
717
718void MImgCleanStd::CleanStepTime() const
719{
720 if (fPostCleanType<=0)
721 return;
722
723 if (fPostCleanType&2)
724 CleanTime(2, fTimeLvl2);
725
726 if (fPostCleanType&1)
727 CleanTime(1, fTimeLvl1);
728}
729
730// --------------------------------------------------------------------------
731//
732// Check if MEvtHeader exists in the Parameter list already.
733// if not create one and add them to the list
734//
735Int_t MImgCleanStd::PreProcess (MParList *pList)
736{
737 fCam = (MGeomCam*)pList->FindObject(AddSerialNumber(fNameGeomCam), "MGeomCam");
738 if (!fCam)
739 {
740 *fLog << err << fNameGeomCam << " [MGeomCam] not found (no geometry information available)... aborting." << endl;
741 return kFALSE;
742 }
743 fEvt = (MSignalCam*)pList->FindObject(AddSerialNumber(fNameSignalCam), "MSignalCam");
744 if (!fEvt)
745 {
746 *fLog << err << fNameSignalCam << " [MSignalCam] not found... aborting." << endl;
747 return kFALSE;
748 }
749
750 fPed=0;
751 if (fCleaningMethod!=kAbsolute && fCleaningMethod!=kTime)
752 {
753 fPed = (MPedPhotCam*)pList->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam");
754 if (!fPed)
755 {
756 *fLog << err << fNamePedPhotCam << " [MPedPhotCam] not found... aborting." << endl;
757 return kFALSE;
758 }
759 }
760
761 fData = (MCameraData*)pList->FindCreateObj(AddSerialNumber("MCameraData"));
762 if (!fData)
763 return kFALSE;
764
765 if (fCleanLvl2>fCleanLvl1)
766 {
767 *fLog << warn << "WARNING - fCleanLvl2 (" << fCleanLvl2 << ") > fCleanLvl1 (" << fCleanLvl1 << ")... resetting fCleanLvl2." << endl;
768 fCleanLvl2 = fCleanLvl1;
769 }
770
771 if (fCleanLvl2==fCleanLvl1 && fCleanRings>0)
772 {
773 *fLog << warn << "WARNING - fCleanLvl2 == fCleanLvl1 (" << fCleanLvl1 << ") but fCleanRings>0... resetting fCleanRings to 0." << endl;
774 fCleanRings=0;
775 }
776
777 if (fKeepIsolatedPixels && fTimeLvl2<fCleanLvl1)
778 {
779 *fLog << warn << "WARNING - fKeepIsolatedPixels set but CleanLvl0 (" << fTimeLvl2 << ") < fCleanLvl1 (" << fCleanLvl1 << ")... resetting fTimeLvl2." << endl;
780 fTimeLvl2 = fCleanLvl1;
781 }
782 if (!fKeepIsolatedPixels && fTimeLvl2>fCleanLvl1)
783 {
784 *fLog << warn << "WARNING - fKeepIsolatedPixels not set but CleanLvl0 (" << fTimeLvl2 << ") > fCleanLvl1 (" << fCleanLvl1 << ")... setting fKeepIsolatedCorePixels." << endl;
785 fKeepIsolatedPixels=kTRUE;
786 }
787
788 if (fKeepIsolatedPixels && fTimeLvl2==fCleanLvl1 && fRecoverIsolatedPixels!=0)
789 {
790 *fLog << warn << "WARNING - fTimeLvl2 == fCleanLvl1 (" << fTimeLvl2 << ") and fKeepSinglePixels and fRecoverIsolatedPixels!=0... setting fRecoverIsolatedPixels=0." << endl;
791 fRecoverIsolatedPixels = 0;
792 }
793
794 Print();
795
796 return kTRUE;
797}
798
799// --------------------------------------------------------------------------
800//
801// Cleans the image.
802//
803Int_t MImgCleanStd::Process()
804{
805 switch (fCleaningMethod)
806 {
807 case kStandard:
808 fData->CalcCleaningLevel(*fEvt, *fPed, *fCam);
809 break;
810 case kScaled:
811 fData->CalcCleaningLevel2(*fEvt, *fPed, *fCam);
812 break;
813 case kDemocratic:
814 fData->CalcCleaningLevelDemocratic(*fEvt, *fPed, *fCam);
815 break;
816 case kProbability:
817 fData->CalcCleaningProbability(*fEvt, *fPed, *fCam);
818 break;
819 case kAbsolute:
820 fData->CalcCleaningAbsolute(*fEvt, *fCam);
821 break;
822 case kTime:
823 fData->CalcCleaningArrivalTime(*fEvt, *fCam);
824 break;
825 default:
826 break;
827 }
828
829#ifdef DEBUG
830 *fLog << all << "ResetCleaning" << endl;
831#endif
832 ResetCleaning();
833
834#ifdef DEBUG
835 *fLog << all << "DoCleaning" << endl;
836#endif
837 Float_t size;
838 Short_t n = DoCleaning(size);
839
840#ifdef DEBUG
841 *fLog << all << "RecoverIsolatedPixels" << endl;
842#endif
843 for (UInt_t i=0; i<(UInt_t)fRecoverIsolatedPixels; i++)
844 {
845 const Int_t rc=RecoverIsolatedPixels(size);
846 if (rc==0)
847 break;
848
849 n -= rc;
850 }
851
852#ifdef DEBUG
853 *fLog << all << "Time Cleaning" << endl;
854#endif
855 // FIXME: Remove removed core piselx?
856 CleanStepTime();
857
858 fEvt->SetSinglePixels(n, size);
859
860#ifdef DEBUG
861 *fLog << all << "CalcIslands" << endl;
862#endif
863 // Takes roughly 10% of the time
864 fEvt->CalcIslands(*fCam);
865
866#ifdef DEBUG
867 *fLog << all << "Done." << endl;
868#endif
869
870 return kTRUE;
871}
872
873// --------------------------------------------------------------------------
874//
875// Print descriptor and cleaning levels.
876//
877void MImgCleanStd::Print(Option_t *o) const
878{
879 *fLog << all << GetDescriptor() << " using ";
880 switch (fCleaningMethod)
881 {
882 case kDemocratic:
883 *fLog << "democratic";
884 break;
885 case kStandard:
886 *fLog << "standard";
887 break;
888 case kScaled:
889 *fLog << "scaled";
890 break;
891 case kProbability:
892 *fLog << "probability";
893 break;
894 case kAbsolute:
895 *fLog << "absolute";
896 break;
897 case kTime:
898 *fLog << "time";
899 break;
900
901 }
902 *fLog << " cleaning" << endl;
903 *fLog << "initialized with level " << fCleanLvl1 << " and " << fCleanLvl2;
904 *fLog << " (CleanRings=" << fCleanRings << ")" << endl;
905
906 if (fPostCleanType)
907 {
908 *fLog << "Time post cleaning is switched on with:" << endl;
909 if (fPostCleanType&2)
910 *fLog << " - Two pixel coincidence window: " << fTimeLvl2 << "ns" << endl;
911 if (fPostCleanType&1)
912 *fLog << " - One pixel coincidence window: " << fTimeLvl1 << "ns" << endl;
913 }
914
915 if (fKeepIsolatedPixels)
916 *fLog << "isolated core pixels will be kept above " << fCleanLvl0 << endl;
917
918 if (fRecoverIsolatedPixels)
919 {
920 if (fRecoverIsolatedPixels<0)
921 *fLog << "all ";
922 *fLog << "isolated core pixels with used pixels as neighbors will be recovered";
923 if (fRecoverIsolatedPixels>0)
924 *fLog << " " << fRecoverIsolatedPixels << " times";
925 *fLog << "." << endl;;
926 }
927
928 if (fCleaningMethod!=kAbsolute && fCleaningMethod!=kTime)
929 {
930 *fLog << "Name of MPedPhotCam container used: ";
931 *fLog << (fPed?((MParContainer*)fPed)->GetName():(const char*)fNamePedPhotCam) << endl;
932 }
933}
934
935// --------------------------------------------------------------------------
936//
937// Create two text entry fields, one for each cleaning level and a
938// describing text line.
939//
940void MImgCleanStd::CreateGuiElements(MGGroupFrame *f)
941{
942 //
943 // Create a frame for line 3 and 4 to be able
944 // to align entry field and label in one line
945 //
946 TGHorizontalFrame *f1 = new TGHorizontalFrame(f, 0, 0);
947 TGHorizontalFrame *f2 = new TGHorizontalFrame(f, 0, 0);
948
949 /*
950 * --> use with root >=3.02 <--
951 *
952
953 TGNumberEntry *fNumEntry1 = new TGNumberEntry(frame, 3.0, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
954 TGNumberEntry *fNumEntry2 = new TGNumberEntry(frame, 2.5, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
955
956 */
957 TGTextEntry *entry1 = new TGTextEntry(f1, "****", kImgCleanLvl1);
958 TGTextEntry *entry2 = new TGTextEntry(f2, "****", kImgCleanLvl2);
959
960 // --- doesn't work like expected (until root 3.02?) --- fNumEntry1->SetAlignment(kTextRight);
961 // --- doesn't work like expected (until root 3.02?) --- fNumEntry2->SetAlignment(kTextRight);
962
963 entry1->SetText("3.0");
964 entry2->SetText("2.5");
965
966 entry1->Associate(f);
967 entry2->Associate(f);
968
969 TGLabel *l1 = new TGLabel(f1, "Cleaning Level 1");
970 TGLabel *l2 = new TGLabel(f2, "Cleaning Level 2");
971
972 l1->SetTextJustify(kTextLeft);
973 l2->SetTextJustify(kTextLeft);
974
975 //
976 // Align the text of the label centered, left in the row
977 // with a left padding of 10
978 //
979 TGLayoutHints *laylabel = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 10);
980 TGLayoutHints *layframe = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 5, 0, 10);
981
982 //
983 // Add one entry field and the corresponding label to each line
984 //
985 f1->AddFrame(entry1);
986 f2->AddFrame(entry2);
987
988 f1->AddFrame(l1, laylabel);
989 f2->AddFrame(l2, laylabel);
990
991 f->AddFrame(f1, layframe);
992 f->AddFrame(f2, layframe);
993
994 f->AddToList(entry1);
995 f->AddToList(entry2);
996 f->AddToList(l1);
997 f->AddToList(l2);
998 f->AddToList(laylabel);
999 f->AddToList(layframe);
1000}
1001
1002// --------------------------------------------------------------------------
1003//
1004// Process the GUI Events coming from the two text entry fields.
1005//
1006Bool_t MImgCleanStd::ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2)
1007{
1008 if (msg!=kC_TEXTENTRY || submsg!=kTE_ENTER)
1009 return kTRUE;
1010
1011 TGTextEntry *txt = (TGTextEntry*)FindWidget(param1);
1012
1013 if (!txt)
1014 return kTRUE;
1015
1016 Float_t lvl = atof(txt->GetText());
1017
1018 switch (param1)
1019 {
1020 case kImgCleanLvl1:
1021 fCleanLvl1 = lvl;
1022 *fLog << "Cleaning level 1 set to " << lvl << endl;
1023 return kTRUE;
1024
1025 case kImgCleanLvl2:
1026 fCleanLvl2 = lvl;
1027 *fLog << "Cleaning level 2 set to " << lvl << endl;
1028 return kTRUE;
1029 }
1030
1031 return kTRUE;
1032}
1033
1034// --------------------------------------------------------------------------
1035//
1036// Implementation of SavePrimitive. Used to write the call to a constructor
1037// to a macro. In the original root implementation it is used to write
1038// gui elements to a macro-file.
1039//
1040void MImgCleanStd::StreamPrimitive(ostream &out) const
1041{
1042 out << " MImgCleanStd " << GetUniqueName() << "(";
1043 out << fCleanLvl1 << ", " << fCleanLvl2;
1044
1045 if (fName!=gsDefName || fTitle!=gsDefTitle)
1046 {
1047 out << ", \"" << fName << "\"";
1048 if (fTitle!=gsDefTitle)
1049 out << ", \"" << fTitle << "\"";
1050 }
1051 out << ");" << endl;
1052
1053 if (fCleaningMethod!=kStandard)
1054 {
1055 out << " " << GetUniqueName() << ".SetMethod(MImgCleanStd::k";
1056 switch (fCleaningMethod)
1057 {
1058 case kScaled: out << "Scaled"; break;
1059 case kDemocratic: out << "Democratic"; break;
1060 case kProbability: out << "Probability"; break;
1061 case kAbsolute: out << "Absolute"; break;
1062 case kTime : out << "Time"; break;
1063 default:
1064 break;
1065 }
1066 out << ");" << endl;
1067 }
1068 if (fCleanRings!=1)
1069 out << " " << GetUniqueName() << ".SetCleanRings(" << fCleanRings << ");" << endl;
1070
1071 if (gsNamePedPhotCam!=fNamePedPhotCam)
1072 out << " " << GetUniqueName() << ".SetNamePedPhotCam(\"" << fNamePedPhotCam << "\");" << endl;
1073 if (gsNameGeomCam!=fNameGeomCam)
1074 out << " " << GetUniqueName() << ".SetNameGeomCam(\"" << fNameGeomCam << "\");" << endl;
1075 if (gsNameSignalCam!=fNameSignalCam)
1076 out << " " << GetUniqueName() << ".SetNameSignalCam(\"" << fNameSignalCam << "\");" << endl;
1077 if (fKeepIsolatedPixels)
1078 out << " " << GetUniqueName() << ".SetKeepIsolatedPixels();" << endl;
1079 if (fRecoverIsolatedPixels)
1080 out << " " << GetUniqueName() << ".SetRecoverIsolatedPixels(" << fRecoverIsolatedPixels << ");" << endl;
1081
1082}
1083
1084// --------------------------------------------------------------------------
1085//
1086// Read the setup from a TEnv, eg:
1087// MImgCleanStd.CleanLevel1: 3.0
1088// MImgCleanStd.CleanLevel2: 2.5
1089// MImgCleanStd.CleanMethod: Standard, Scaled, Democratic, Probability, Absolute
1090// MImgCleanStd.CleanRings: 1
1091// MImgCleanStd.KeepIsolatedPixels: yes, no
1092//
1093Int_t MImgCleanStd::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
1094{
1095 Bool_t rc = kFALSE;
1096 if (IsEnvDefined(env, prefix, "CleanRings", print))
1097 {
1098 rc = kTRUE;
1099 SetCleanRings(GetEnvValue(env, prefix, "CleanRings", fCleanRings));
1100 }
1101 if (IsEnvDefined(env, prefix, "CleanLevel0", print))
1102 {
1103 rc = kTRUE;
1104 fCleanLvl0 = GetEnvValue(env, prefix, "CleanLevel0", fCleanLvl0);
1105 }
1106 if (IsEnvDefined(env, prefix, "CleanLevel1", print))
1107 {
1108 rc = kTRUE;
1109 fCleanLvl1 = GetEnvValue(env, prefix, "CleanLevel1", fCleanLvl1);
1110 }
1111 if (IsEnvDefined(env, prefix, "CleanLevel2", print))
1112 {
1113 rc = kTRUE;
1114 fCleanLvl2 = GetEnvValue(env, prefix, "CleanLevel2", fCleanLvl2);
1115 }
1116 if (IsEnvDefined(env, prefix, "TimeLevel1", print))
1117 {
1118 rc = kTRUE;
1119 fTimeLvl1 = GetEnvValue(env, prefix, "TimeLevel1", fTimeLvl1);
1120 }
1121 if (IsEnvDefined(env, prefix, "TimeLevel2", print))
1122 {
1123 rc = kTRUE;
1124 fTimeLvl2 = GetEnvValue(env, prefix, "TimeLevel2", fTimeLvl2);
1125 }
1126 if (IsEnvDefined(env, prefix, "KeepIsolatedPixels", print))
1127 {
1128 rc = kTRUE;
1129 fKeepIsolatedPixels = GetEnvValue(env, prefix, "KeepIsolatedPixels", fKeepIsolatedPixels);
1130 }
1131 if (IsEnvDefined(env, prefix, "RecoverIsolatedPixels", print))
1132 {
1133 rc = kTRUE;
1134 fRecoverIsolatedPixels = GetEnvValue(env, prefix, "RecoverIsolatedPixels", fKeepIsolatedPixels);
1135 }
1136 if (IsEnvDefined(env, prefix, "PostCleanType", print))
1137 {
1138 rc = kTRUE;
1139 fPostCleanType = GetEnvValue(env, prefix, "PostCleanType", fPostCleanType);
1140 }
1141
1142 if (IsEnvDefined(env, prefix, "CleanMethod", print))
1143 {
1144 rc = kTRUE;
1145 TString s = GetEnvValue(env, prefix, "CleanMethod", "");
1146 s.ToLower();
1147 if (s.BeginsWith("standard"))
1148 SetMethod(kStandard);
1149 if (s.BeginsWith("scaled"))
1150 SetMethod(kScaled);
1151 if (s.BeginsWith("democratic"))
1152 SetMethod(kDemocratic);
1153 if (s.BeginsWith("probability"))
1154 SetMethod(kProbability);
1155 if (s.BeginsWith("absolute"))
1156 SetMethod(kAbsolute);
1157 if (s.BeginsWith("time"))
1158 SetMethod(kTime);
1159 }
1160
1161 return rc;
1162}
Note: See TracBrowser for help on using the repository browser.