source: trunk/Mars/mimage/MImgCleanStd.cc@ 19344

Last change on this file since 19344 was 9615, checked in by tbretz, 14 years ago
*** empty log message ***
File size: 36.8 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 MGeom &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 MGeom &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 MGeom &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 MGeom &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 MGeom &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 MGeom &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 MGeom &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 // Add conversion factor for dx here
667 const Double_t dtmax = lvl*fCam->GetConvMm2Deg();
668
669 const UInt_t npixevt = fEvt->GetNumPixels();
670 for (UInt_t idx=0; idx<npixevt; idx++)
671 {
672 // check if pixel is used or not
673 const MSignalPix &pix = (*fEvt)[idx];
674 if (!pix.IsPixelUsed())
675 continue;
676
677 // get arrival time
678 const Double_t tm0 = pix.GetArrivalTime();
679
680 // loop over its neighbpors
681 const MGeom &gpix = (*fCam)[idx];
682
683 Int_t cnt = 0;
684 for (Int_t i=0; i<gpix.GetNumNeighbors(); i++)
685 {
686 // Get index of neighbor
687 const Int_t idx2 = gpix.GetNeighbor(i);
688
689 // check if neighbor is used or not
690 const MSignalPix &npix = (*fEvt)[idx2];
691 if (!npix.IsPixelUsed())
692 continue;
693
694 const Double_t dt = TMath::Abs(npix.GetArrivalTime()-tm0);
695 const Double_t dx = gpix.GetDist((*fCam)[idx2]);
696
697 // If this pixel is to far away (in arrival time) don't count
698 if (dt>dtmax*dx)
699 continue;
700
701 // Now count the pixel. If we did not found n pixels yet
702 // which fullfill the condition, go on searching
703 if (++cnt>=n)
704 break;
705 }
706
707 // If we found at least n neighbors which are
708 // with a time difference of lvl keep the pixel
709 if (cnt>=n)
710 continue;
711
712 indices.Set(indices.GetSize()+1);
713 indices[indices.GetSize()-1] = idx;
714 }
715
716 // Now remove the pixels which didn't fullfill the requirement
717 for (UInt_t i=0; i<indices.GetSize(); i++)
718 {
719 (*fEvt)[indices[i]].SetPixelUnused();
720 (*fEvt)[indices[i]].SetPixelCore(kFALSE);
721 }
722}
723
724void MImgCleanStd::CleanStepTime() const
725{
726 if (fPostCleanType<=0)
727 return;
728
729 if (fPostCleanType&2)
730 CleanTime(2, fTimeLvl2);
731
732 if (fPostCleanType&1)
733 CleanTime(1, fTimeLvl1);
734}
735
736// --------------------------------------------------------------------------
737//
738// Check if MEvtHeader exists in the Parameter list already.
739// if not create one and add them to the list
740//
741Int_t MImgCleanStd::PreProcess (MParList *pList)
742{
743 fCam = (MGeomCam*)pList->FindObject(AddSerialNumber(fNameGeomCam), "MGeomCam");
744 if (!fCam)
745 {
746 *fLog << err << fNameGeomCam << " [MGeomCam] not found (no geometry information available)... aborting." << endl;
747 return kFALSE;
748 }
749 fEvt = (MSignalCam*)pList->FindObject(AddSerialNumber(fNameSignalCam), "MSignalCam");
750 if (!fEvt)
751 {
752 *fLog << err << fNameSignalCam << " [MSignalCam] not found... aborting." << endl;
753 return kFALSE;
754 }
755
756 fPed=0;
757 if (fCleaningMethod!=kAbsolute && fCleaningMethod!=kTime)
758 {
759 fPed = (MPedPhotCam*)pList->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam");
760 if (!fPed)
761 {
762 *fLog << err << fNamePedPhotCam << " [MPedPhotCam] not found... aborting." << endl;
763 return kFALSE;
764 }
765 }
766
767 fData = (MCameraData*)pList->FindCreateObj(AddSerialNumber("MCameraData"));
768 if (!fData)
769 return kFALSE;
770
771 if (fCleanLvl2>fCleanLvl1)
772 {
773 *fLog << warn << "WARNING - fCleanLvl2 (" << fCleanLvl2 << ") > fCleanLvl1 (" << fCleanLvl1 << ")... resetting fCleanLvl2." << endl;
774 fCleanLvl2 = fCleanLvl1;
775 }
776
777 if (fCleanLvl2==fCleanLvl1 && fCleanRings>0)
778 {
779 *fLog << warn << "WARNING - fCleanLvl2 == fCleanLvl1 (" << fCleanLvl1 << ") but fCleanRings>0... resetting fCleanRings to 0." << endl;
780 fCleanRings=0;
781 }
782
783 if (fKeepIsolatedPixels && fCleanLvl0<fCleanLvl1)
784 {
785 *fLog << warn << "WARNING - fKeepIsolatedPixels set but CleanLvl0 (" << fTimeLvl2 << ") < fCleanLvl1 (" << fCleanLvl1 << ")... resetting fTimeLvl2." << endl;
786 fTimeLvl2 = fCleanLvl1;
787 }
788 if (!fKeepIsolatedPixels && fCleanLvl0>fCleanLvl1)
789 {
790 *fLog << warn << "WARNING - fKeepIsolatedPixels not set but CleanLvl0 (" << fTimeLvl2 << ") > fCleanLvl1 (" << fCleanLvl1 << ")... setting fKeepIsolatedCorePixels." << endl;
791 fKeepIsolatedPixels=kTRUE;
792 }
793
794 if (fKeepIsolatedPixels && fTimeLvl2==fCleanLvl1 && fRecoverIsolatedPixels!=0)
795 {
796 *fLog << warn << "WARNING - fTimeLvl2 == fCleanLvl1 (" << fTimeLvl2 << ") and fKeepSinglePixels and fRecoverIsolatedPixels!=0... setting fRecoverIsolatedPixels=0." << endl;
797 fRecoverIsolatedPixels = 0;
798 }
799
800 Print();
801
802 return kTRUE;
803}
804
805// --------------------------------------------------------------------------
806//
807// Cleans the image.
808//
809Int_t MImgCleanStd::Process()
810{
811 switch (fCleaningMethod)
812 {
813 case kStandard:
814 fData->CalcCleaningLevel(*fEvt, *fPed, *fCam);
815 break;
816 case kScaled:
817 fData->CalcCleaningLevel2(*fEvt, *fPed, *fCam);
818 break;
819 case kDemocratic:
820 fData->CalcCleaningLevelDemocratic(*fEvt, *fPed, *fCam);
821 break;
822 case kProbability:
823 fData->CalcCleaningProbability(*fEvt, *fPed, *fCam);
824 break;
825 case kAbsolute:
826 fData->CalcCleaningAbsolute(*fEvt, *fCam);
827 break;
828 case kTime:
829 fData->CalcCleaningArrivalTime(*fEvt, *fCam);
830 break;
831 default:
832 break;
833 }
834
835#ifdef DEBUG
836 *fLog << all << "ResetCleaning" << endl;
837#endif
838 ResetCleaning();
839
840#ifdef DEBUG
841 *fLog << all << "DoCleaning" << endl;
842#endif
843 Float_t size;
844 Short_t n = DoCleaning(size);
845
846#ifdef DEBUG
847 *fLog << all << "RecoverIsolatedPixels" << endl;
848#endif
849 for (UInt_t i=0; i<(UInt_t)fRecoverIsolatedPixels; i++)
850 {
851 const Int_t rc=RecoverIsolatedPixels(size);
852 if (rc==0)
853 break;
854
855 n -= rc;
856 }
857
858#ifdef DEBUG
859 *fLog << all << "Time Cleaning" << endl;
860#endif
861 // FIXME: Remove removed core piselx?
862 CleanStepTime();
863
864 fEvt->SetSinglePixels(n, size);
865
866#ifdef DEBUG
867 *fLog << all << "CalcIslands" << endl;
868#endif
869 // Takes roughly 10% of the time
870 fEvt->CalcIslands(*fCam);
871
872#ifdef DEBUG
873 *fLog << all << "Done." << endl;
874#endif
875
876 return kTRUE;
877}
878
879// --------------------------------------------------------------------------
880//
881// Print descriptor and cleaning levels.
882//
883void MImgCleanStd::Print(Option_t *o) const
884{
885 *fLog << all << GetDescriptor() << " using ";
886 switch (fCleaningMethod)
887 {
888 case kDemocratic:
889 *fLog << "democratic";
890 break;
891 case kStandard:
892 *fLog << "standard";
893 break;
894 case kScaled:
895 *fLog << "scaled";
896 break;
897 case kProbability:
898 *fLog << "probability";
899 break;
900 case kAbsolute:
901 *fLog << "absolute";
902 break;
903 case kTime:
904 *fLog << "time";
905 break;
906
907 }
908 *fLog << " cleaning" << endl;
909 *fLog << "initialized with level " << fCleanLvl1 << " and " << fCleanLvl2;
910 *fLog << " (CleanRings=" << fCleanRings << ")" << endl;
911
912 if (fPostCleanType)
913 {
914 *fLog << "Time post cleaning is switched on with:" << endl;
915 if (fPostCleanType&2)
916 *fLog << " - Two pixel coincidence window: " << fTimeLvl2 << "ns" << endl;
917 if (fPostCleanType&1)
918 *fLog << " - One pixel coincidence window: " << fTimeLvl1 << "ns" << endl;
919 }
920
921 if (fKeepIsolatedPixels)
922 *fLog << "isolated core pixels will be kept above " << fCleanLvl0 << endl;
923
924 if (fRecoverIsolatedPixels)
925 {
926 if (fRecoverIsolatedPixels<0)
927 *fLog << "all ";
928 *fLog << "isolated core pixels with used pixels as neighbors will be recovered";
929 if (fRecoverIsolatedPixels>0)
930 *fLog << " " << fRecoverIsolatedPixels << " times";
931 *fLog << "." << endl;;
932 }
933
934 if (fCleaningMethod!=kAbsolute && fCleaningMethod!=kTime)
935 {
936 *fLog << "Name of MPedPhotCam container used: ";
937 *fLog << (fPed?((MParContainer*)fPed)->GetName():(const char*)fNamePedPhotCam) << endl;
938 }
939}
940
941// --------------------------------------------------------------------------
942//
943// Create two text entry fields, one for each cleaning level and a
944// describing text line.
945//
946void MImgCleanStd::CreateGuiElements(MGGroupFrame *f)
947{
948 //
949 // Create a frame for line 3 and 4 to be able
950 // to align entry field and label in one line
951 //
952 TGHorizontalFrame *f1 = new TGHorizontalFrame(f, 0, 0);
953 TGHorizontalFrame *f2 = new TGHorizontalFrame(f, 0, 0);
954
955 /*
956 * --> use with root >=3.02 <--
957 *
958
959 TGNumberEntry *fNumEntry1 = new TGNumberEntry(frame, 3.0, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
960 TGNumberEntry *fNumEntry2 = new TGNumberEntry(frame, 2.5, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
961
962 */
963 TGTextEntry *entry1 = new TGTextEntry(f1, "****", kImgCleanLvl1);
964 TGTextEntry *entry2 = new TGTextEntry(f2, "****", kImgCleanLvl2);
965
966 // --- doesn't work like expected (until root 3.02?) --- fNumEntry1->SetAlignment(kTextRight);
967 // --- doesn't work like expected (until root 3.02?) --- fNumEntry2->SetAlignment(kTextRight);
968
969 entry1->SetText("3.0");
970 entry2->SetText("2.5");
971
972 entry1->Associate(f);
973 entry2->Associate(f);
974
975 TGLabel *l1 = new TGLabel(f1, "Cleaning Level 1");
976 TGLabel *l2 = new TGLabel(f2, "Cleaning Level 2");
977
978 l1->SetTextJustify(kTextLeft);
979 l2->SetTextJustify(kTextLeft);
980
981 //
982 // Align the text of the label centered, left in the row
983 // with a left padding of 10
984 //
985 TGLayoutHints *laylabel = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 10);
986 TGLayoutHints *layframe = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 5, 0, 10);
987
988 //
989 // Add one entry field and the corresponding label to each line
990 //
991 f1->AddFrame(entry1);
992 f2->AddFrame(entry2);
993
994 f1->AddFrame(l1, laylabel);
995 f2->AddFrame(l2, laylabel);
996
997 f->AddFrame(f1, layframe);
998 f->AddFrame(f2, layframe);
999
1000 f->AddToList(entry1);
1001 f->AddToList(entry2);
1002 f->AddToList(l1);
1003 f->AddToList(l2);
1004 f->AddToList(laylabel);
1005 f->AddToList(layframe);
1006}
1007
1008// --------------------------------------------------------------------------
1009//
1010// Process the GUI Events coming from the two text entry fields.
1011//
1012Bool_t MImgCleanStd::ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2)
1013{
1014 if (msg!=kC_TEXTENTRY || submsg!=kTE_ENTER)
1015 return kTRUE;
1016
1017 TGTextEntry *txt = (TGTextEntry*)FindWidget(param1);
1018
1019 if (!txt)
1020 return kTRUE;
1021
1022 Float_t lvl = atof(txt->GetText());
1023
1024 switch (param1)
1025 {
1026 case kImgCleanLvl1:
1027 fCleanLvl1 = lvl;
1028 *fLog << "Cleaning level 1 set to " << lvl << endl;
1029 return kTRUE;
1030
1031 case kImgCleanLvl2:
1032 fCleanLvl2 = lvl;
1033 *fLog << "Cleaning level 2 set to " << lvl << endl;
1034 return kTRUE;
1035 }
1036
1037 return kTRUE;
1038}
1039
1040// --------------------------------------------------------------------------
1041//
1042// Implementation of SavePrimitive. Used to write the call to a constructor
1043// to a macro. In the original root implementation it is used to write
1044// gui elements to a macro-file.
1045//
1046void MImgCleanStd::StreamPrimitive(ostream &out) const
1047{
1048 out << " MImgCleanStd " << GetUniqueName() << "(";
1049 out << fCleanLvl1 << ", " << fCleanLvl2;
1050
1051 if (fName!=gsDefName || fTitle!=gsDefTitle)
1052 {
1053 out << ", \"" << fName << "\"";
1054 if (fTitle!=gsDefTitle)
1055 out << ", \"" << fTitle << "\"";
1056 }
1057 out << ");" << endl;
1058
1059 if (fCleaningMethod!=kStandard)
1060 {
1061 out << " " << GetUniqueName() << ".SetMethod(MImgCleanStd::k";
1062 switch (fCleaningMethod)
1063 {
1064 case kScaled: out << "Scaled"; break;
1065 case kDemocratic: out << "Democratic"; break;
1066 case kProbability: out << "Probability"; break;
1067 case kAbsolute: out << "Absolute"; break;
1068 case kTime : out << "Time"; break;
1069 default:
1070 break;
1071 }
1072 out << ");" << endl;
1073 }
1074 if (fCleanRings!=1)
1075 out << " " << GetUniqueName() << ".SetCleanRings(" << fCleanRings << ");" << endl;
1076
1077 if (gsNamePedPhotCam!=fNamePedPhotCam)
1078 out << " " << GetUniqueName() << ".SetNamePedPhotCam(\"" << fNamePedPhotCam << "\");" << endl;
1079 if (gsNameGeomCam!=fNameGeomCam)
1080 out << " " << GetUniqueName() << ".SetNameGeomCam(\"" << fNameGeomCam << "\");" << endl;
1081 if (gsNameSignalCam!=fNameSignalCam)
1082 out << " " << GetUniqueName() << ".SetNameSignalCam(\"" << fNameSignalCam << "\");" << endl;
1083 if (fKeepIsolatedPixels)
1084 out << " " << GetUniqueName() << ".SetKeepIsolatedPixels();" << endl;
1085 if (fRecoverIsolatedPixels)
1086 out << " " << GetUniqueName() << ".SetRecoverIsolatedPixels(" << fRecoverIsolatedPixels << ");" << endl;
1087
1088}
1089
1090// --------------------------------------------------------------------------
1091//
1092// Read the setup from a TEnv, eg:
1093// MImgCleanStd.CleanLevel1: 3.0
1094// MImgCleanStd.CleanLevel2: 2.5
1095// MImgCleanStd.CleanMethod: Standard, Scaled, Democratic, Probability, Absolute
1096// MImgCleanStd.CleanRings: 1
1097// MImgCleanStd.KeepIsolatedPixels: yes, no
1098//
1099Int_t MImgCleanStd::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
1100{
1101 Bool_t rc = kFALSE;
1102 if (IsEnvDefined(env, prefix, "CleanRings", print))
1103 {
1104 rc = kTRUE;
1105 SetCleanRings(GetEnvValue(env, prefix, "CleanRings", fCleanRings));
1106 }
1107 if (IsEnvDefined(env, prefix, "CleanLevel0", print))
1108 {
1109 rc = kTRUE;
1110 fCleanLvl0 = GetEnvValue(env, prefix, "CleanLevel0", fCleanLvl0);
1111 }
1112 if (IsEnvDefined(env, prefix, "CleanLevel1", print))
1113 {
1114 rc = kTRUE;
1115 fCleanLvl1 = GetEnvValue(env, prefix, "CleanLevel1", fCleanLvl1);
1116 }
1117 if (IsEnvDefined(env, prefix, "CleanLevel2", print))
1118 {
1119 rc = kTRUE;
1120 fCleanLvl2 = GetEnvValue(env, prefix, "CleanLevel2", fCleanLvl2);
1121 }
1122 if (IsEnvDefined(env, prefix, "TimeLevel1", print))
1123 {
1124 rc = kTRUE;
1125 fTimeLvl1 = GetEnvValue(env, prefix, "TimeLevel1", fTimeLvl1);
1126 }
1127 if (IsEnvDefined(env, prefix, "TimeLevel2", print))
1128 {
1129 rc = kTRUE;
1130 fTimeLvl2 = GetEnvValue(env, prefix, "TimeLevel2", fTimeLvl2);
1131 }
1132 if (IsEnvDefined(env, prefix, "KeepIsolatedPixels", print))
1133 {
1134 rc = kTRUE;
1135 fKeepIsolatedPixels = GetEnvValue(env, prefix, "KeepIsolatedPixels", fKeepIsolatedPixels);
1136 }
1137 if (IsEnvDefined(env, prefix, "RecoverIsolatedPixels", print))
1138 {
1139 rc = kTRUE;
1140 fRecoverIsolatedPixels = GetEnvValue(env, prefix, "RecoverIsolatedPixels", fKeepIsolatedPixels);
1141 }
1142 if (IsEnvDefined(env, prefix, "PostCleanType", print))
1143 {
1144 rc = kTRUE;
1145 fPostCleanType = GetEnvValue(env, prefix, "PostCleanType", fPostCleanType);
1146 }
1147
1148 if (IsEnvDefined(env, prefix, "CleanMethod", print))
1149 {
1150 rc = kTRUE;
1151 TString s = GetEnvValue(env, prefix, "CleanMethod", "");
1152 s.ToLower();
1153 if (s.BeginsWith("standard"))
1154 SetMethod(kStandard);
1155 if (s.BeginsWith("scaled"))
1156 SetMethod(kScaled);
1157 if (s.BeginsWith("democratic"))
1158 SetMethod(kDemocratic);
1159 if (s.BeginsWith("probability"))
1160 SetMethod(kProbability);
1161 if (s.BeginsWith("absolute"))
1162 SetMethod(kAbsolute);
1163 if (s.BeginsWith("time"))
1164 SetMethod(kTime);
1165 }
1166
1167 return rc;
1168}
Note: See TracBrowser for help on using the repository browser.