Changeset 2035


Ignore:
Timestamp:
04/28/03 19:02:00 (22 years ago)
Author:
tonello
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r2034 r2035  
    11                                                 -*-*- END OF LINE -*-*-
     2 2003/04/28: Nadia Tonello
     3
     4  * mimage/MImgCleanStd.[h,cc]
     5     - added the option kDemocratic using sigmabar of the
     6        inner pixels
     7     - added the option to select the number of rings of pixels
     8        to analyze around the core pixels
     9     - added documentation
     10
     11  * manalysis/MCerPhotPix.[h,cc]
     12     - added fRing and Get-Set functions
     13 
    214 2003/04/28: Oscar Blanch
    315
     
    2032       declared external.
    2133
    22  2003/04/27: Thomas Bretz
     34 2003/04/28: Thomas Bretz
    2335
    2436   * NEWS:
  • trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc

    r2014 r2035  
    2828//  MImgCleanStd                                                           //
    2929//                                                                         //
    30 //  This is the standard image cleaning. If you want to know how it works  //
    31 //  Please look at the three CleanSteps and Process                        //
     30// The Image Cleaning task selects the pixels you use for the Hillas       //
     31// parameters calculation.                                                 //
     32//                                                                         //
     33// There are two methods to make the selection: the standard one, as done  //
     34// in the analysis of CT1 data, and the democratic one, as suggested by    //
     35// W.Wittek. The number of photo-electrons of a pixel is compared with the //
     36// pedestal RMS of the pixel itself (standard method) or with the average  //
     37// RMS of the inner pixels (democratic method).                            //
     38// In both cases, the possibility to have a camera with pixels of          //
     39// different area is taken into account.                                   //
     40// The too noisy pixels can be recognized and eventally switched off       //
     41// (Unmap: set blinb pixels to UNUSED) separately, using the               //
     42// MBlindPixelCalc Class. In the MBlindPixelCalc class there is also the   //
     43// function to replace the value of the noisy pixels with the interpolation//
     44// of the content of the neighbors (SetUseInterpolation).                  //
     45//                                                                         //
     46// Example:                                                                //
     47// ...                                                                     //
     48// MBlindPixelCalc blind;                                                  //
     49// blind.SetUseInterpolation();                                            //
     50// blind.SetUseBlindPixels();                                              //
     51//                                                                         //
     52// MImgCleanStd clean;                                                     //
     53// ...                                                                     //
     54// tlist.AddToList(&blind);                                                //
     55// tlist.AddToList(&clean);                                                //
     56//                                                                         //
     57// Look at the MBlindPixelCalc Class for more details.                     //   
     58//                                                                         //
     59// Starting point: default values ---------------------------------------- //
     60//                                                                         //
     61// When an event is read, before the image cleaning, all the pixels that   //
     62// are in MCerPhotEvt are set as USED and NOT CORE. All the pixels belong  //
     63// to RING number 1 (like USED pixels).                                    //
     64// Look at MCerPhotPix.h to see how these informations of the pixel are    //
     65// stored.                                                                 //
     66// The default  cleaning METHOD is the STANDARD one and the number of the  //
     67// rings around the CORE pixel it analyzes is 1. Look at the Constructor   //
     68// of the class in MImgCleanStd.cc to see (or change) the default values.  //
     69//                                                                         //
     70// Example: To modify this setting, use the member functions               //
     71// SetMethod(MImgCleanStd::kDemocratic) and SetCleanRings(UShort_t n).     //
     72//                                                                         //
     73// MImgCleanStd:CleanStep1 -------------------------------------------- ---//
     74//                                                                         //
     75// The first step of cleaning defines the CORE pixels. The CORE pixels are //
     76// the ones which contain the informations about the core of the electro-  //
     77// magnetic shower.                                                        //
     78// The ratio  (A_0/A_i) is calculated from fCam->GetPixRatio(i). A_0 is    //
     79// the area of the central pixel of the camera, A_i is the area of the     //
     80// examined pixel. In this way, if we have a MAGIC-like camera, with the   //
     81// outer pixels bigger than the inner ones, the level of cleaning in the   //
     82// two different regions is weighted.                                      //
     83// This avoids problems of deformations of the shower images.              //
     84// The signal S_i and the pedestal RMS Prms_i of the pixel are called from //
     85// the object MCerPhotPix.                                                 //
     86// If (default method = kStandard)                                         //
     87//                                                                         //
     88//         S_i > L_1 * Prms_i / sqrt(A_0/A_i)                              //
     89//                                                                         //
     90// the pixel is set as CORE pixel. L_1 is called "first level of cleaning" //
     91// (default: fCleanLvl1 = 3).                                              //
     92// All the other pixels are set as UNUSED and belong to RING 0.            //
     93// After this point, only the CORE pixels are set as USED, with RING       //
     94// number 1.                                                               //
     95//                                                                         //
     96// MImgCleanStd:CleanStep2 ----------------------------------------------  //
     97//                                                                         //   
     98// The second step of cleaning looks at the isolated CORE pixels and sets  //
     99// them to UNUSED. An isolated pixel is a pixel without CORE neighbors.    //
     100// At the end of this point, we have set as USED only CORE pixels with at  //
     101// least one CORE neighbor.                                                //
     102//                                                                         //       
     103// MImgCleanStd:CleanStep3  ---------------------------------------------- //
     104//                                                                         //
     105// The third step of cleaning looks at all the pixels (USED or UNUSED) that//
     106// surround the USED pixels.                                               //
     107// If the content of the analyzed pixel survives at the second level of    //
     108// cleaning, i.e. if                                                       //
     109//                                                                         //
     110//          S_i > L_2 * Prms_i /  sqrt(A_0/A_i)                            //
     111//                                                                         //
     112// the pixel is set as USED. L_2 is called "second level of cleaning"      //
     113// (default:fCleanLvl2 = 2.5).                                             //
     114//                                                                         //
     115// When the number of RINGS to analyze is 1 (default value), only the      //
     116// pixels that have a neighbor CORE pixel are analyzed.                    //
     117//                                                                         //
     118// There is the option to decide the number of times you want to repeat    //
     119// this procedure (number of RINGS analyzed around the core pixels = n).   //
     120// Every time the level of cleaning is the same (fCleanLvl2) and the pixel //
     121// will belong to ring r+1, 1 < r < n+1. This is described in              //
     122// MImgCleanStd:CleanStep4 .                                               //
     123//                                                                         //
     124// Dictionary and member functions --------------------------------------- //
     125//                                                                         //
     126// Here there is the detailed description of the member functions and of   //
     127// the terms commonly used in the class.                                   // 
     128//                                                                         //
     129// STANDARD CLEANING:                                                      //
     130// =================                                                       //
     131// This is the method used for the CT1 data analysis. It is the default    //
     132// method of the class.                                                    //
     133// The number of photo-electrons of a pixel (S_i) is compared to the       //
     134// pedestal RMS of the pixel itself (Prms_i). To have the comparison to    //
     135// the same photon density for all the pixels, taking into account they    //
     136// can have different areas, we have to keep in mind that the number of    //
     137// photons that hit each pixel, goes linearly with the area of the pixel.  //
     138// The fluctuations of the LONS are proportional to sqrt(A_i), so when we  //
     139// compare S_i with Prms_i, only a factor  sqrt(A_0/A_i) is missing to     //
     140// have the same (N.photons/Area) threshold for all the pixels.            //
     141//                                                                         //
     142//              !!WARNING: if noise independent from the                   //
     143//                        pixel size is considered,                        //
     144//                  this weight can give wrong results!!                   //
     145// If                                                                      //
     146//                                                                         //
     147//           S_i > L_n * Prms_i / sqrt(A_0/A_i)                            //
     148//                                                                         //
     149// the pixel survives the cleaning and it is set as CORE (when L_n is the  //
     150// first level of cleaning, fCleanLvl1) or USED ((when L_n is the second   //
     151// level of cleaning, fCleanLvl2).                                         //
     152//                                                                         //
     153// Example:                                                                //
     154//                                                                         //
     155// MImgCleanStd clean;                                                     //
     156// //creates a default Cleaning object, with default setting               //
     157// ...                                                                     //
     158// tlist.AddToList(&clean);                                                //
     159// // add the image cleaning to the main task list                         //
     160//                                                                         //
     161// DEMOCRATIC CLEANING:                                                    //
     162// ===================                                                     //
     163// You use this cleaning method when you want to compare the number of     //
     164// photo-electons of each pixel with the average pedestal RMS              //
     165// (fInnerNoise = fSgb->GetSigmabarInner()) of the inner pixels (for the   //
     166// MAGIC camera they are the smaller ones):                                //
     167//                                                                         //
     168//  S_i > L_n * fInnerNoise / (A_0/A_i)                                    //
     169//                                                                         //
     170// In this case, the simple ratio (A_0/A_i) is used to weight the level of //
     171// cleaning, because both the inner and the outer pixels (that in MAGIC    //
     172// have a different area) are compared to the same pedestal RMS, coming    //
     173// from the inner pixels.                                                  //
     174// To calculate the average pedestal RMS of the inner pixels, you have to  //
     175// add to the main task list an object of type MSigmabarCalc before the    //
     176// MImgCleanStd object. To know how the calculation of fInnerNoise is done //
     177// look at the MSigmabarCalc Class.                                        //
     178//                                                                         //
     179// Example:                                                                //
     180//                                                                         //
     181// MSigmabarCalc   sbcalc;                                                 //
     182// //creates an object that calcutates the average pedestal RMS            //
     183// MImgCleanStd clean;                                                     //
     184// ...                                                                     //
     185// tlist.AddToList(&sbcalc);                                               //
     186// tlist.AddToList(&clean);                                                //
     187//                                                                         //
     188// Member Function:  SetMethod()                                           //
     189// ============================                                            //
     190// When you call the MImgCleanStd task, the default method is kStandard.   //
     191//                                                                         //
     192// If you want to switch to the kDemocratic method you have to             //
     193// call this member function.                                              //
     194//                                                                         //
     195// Example:                                                                //
     196//                                                                         //
     197// MImgCleanStd clean;                                                     //
     198// //creates a default Cleaning object, with default setting               //
     199//                                                                         //
     200// clean.SetMethod(MImgCleanStd::kDemocratic);                             //
     201// //now the method of cleaning is changed to Democratic                   //
     202//                                                                         //
     203// FIRST AND SECOND CLEANING LEVEL                                         //
     204// ===============================                                         //
     205// When you call the MImgCleanStd task, the default cleaning levels are    //
     206// fCleanLvl1 = 3, fCleanLvl2 = 2.5. You can change them easily when you   //
     207// create the MImgCleanStd object.                                         //
     208//                                                                         //
     209// Example:                                                                //
     210//                                                                         //
     211// MImgCleanStd clean(Float_t lvl1,Float_t lvl2);                          //
     212// //creates a default cleaning object, but the cleaning levels are now    //
     213// //lvl1 and lvl2.                                                        //
     214//                                                                         //
     215// RING NUMBER                                                             //
     216// ===========                                                             //
     217// The standard cleaning procedure is such that it looks for the           //
     218// informations of the boundary part of the shower only on the first       //
     219// neighbors of the CORE pixels.                                           //
     220// There is the possibility now to look not only at the firs neighbors     //
     221// (first ring),but also further away, around the CORE pixels. All the new //
     222// pixels you can find with this method, are tested with the second level  // 
     223// of cleaning and have to have at least an USED neighbor.                 //
     224//                                                                         //
     225// They will be also set as USED and will be taken into account during the //
     226// calculation of the image parameters.                                    //
     227// The only way to distinguish them from the other USED pixels, is the     //
     228// Ring number, that is bigger than 1.                                     //
     229//                                                                         //
     230// Example: You can decide how many rings you want to analyze using:       //
     231//                                                                         //
     232// MImgCleanStd clean;                                                     //
     233// //creates a default cleaning object (default number of rings =1)        //
     234// clean.SetCleanRings(UShort\_t r);                                       //
     235// //now it looks r times around the CORE pixels to find new pixels with   //
     236// //signal.                                                               //
     237//                                                                         //
    32238//                                                                         //
    33239//   FIXME: MImgCleanStd is not yet completely optimized for speed.        //
     
    35241//                                                                         //
    36242//  Input Containers:                                                      //
    37 //   MGeomCam, MCerPhotEvt                                                 //
     243//   MGeomCam, MCerPhotEvt, MSigmabar                                      //
    38244//                                                                         //
    39245//  Output Containers:                                                     //
    40246//   -/-                                                                   //
    41247//                                                                         //
    42 /////////////////////////////////////////////////////////////////////////////
    43 #include "MImgCleanStd.h"
    44 
    45 #include <stdlib.h>       // atof
     248////////////////////////////////////////////////////////////////////////// //
     249#include "MImgCleanStd.h"                                                         
     250#include <stdlib.h>       // atof                                         
    46251#include <fstream.h>      // ofstream, SavePrimitive
    47252
     
    58263#include "MCerPhotPix.h"
    59264#include "MCerPhotEvt.h"
     265#include "MSigmabar.h"
    60266
    61267#include "MGGroupFrame.h" // MGGroupFrame
     
    69275
    70276static const TString gsDefName  = "MImgCleanStd";
    71 static const TString gsDefTitle = "Perform standard image cleaning";
    72 
    73 // --------------------------------------------------------------------------
    74 //
    75 // Default constructor. Here you can specify the cleaning levels. If you
    76 // don't specify them the 'common standard' values 3.0 and 2.5 (sigma
    77 // above mean) are used
     277static const TString gsDefTitle = "Task to perform image cleaning";
     278
     279// --------------------------------------------------------------------------
     280//
     281// Default constructor. Here you can specify the cleaning method and levels.
     282// If you don't specify them the 'common standard' values 3.0 and 2.5 (sigma
     283// above mean) are used.
     284// Here you can also specify how many rings around the core pixels you want
     285// to analyze (with the fixed lvl2). The default value for "rings" is 1.
    78286//
    79287MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2,
    80                            const char *name, const char *title)
    81     : fCleanLvl1(lvl1), fCleanLvl2(lvl2)
     288                     const char *name, const char *title)
     289    : fSgb(NULL), fCleaningMethod(kStandard), fCleanLvl1(lvl1),
     290    fCleanLvl2(lvl2), fCleanRings(1)
     291
    82292{
    83293    fName  = name  ? name  : gsDefName.Data();
     
    89299// --------------------------------------------------------------------------
    90300//
    91 //  This method looks for all pixels with an entry (photons)
    92 //  that is three times bigger than the noise of the pixel
    93 //  (std: 3 sigma, clean level 1)
    94 //
    95 //
    96 //  AM 18/11/2002: now cut levels are proportional to the square root
    97 //  of the pixel area. In this way the cut corresponds to a fixed
     301//  NT 28/04/2003: now the option to use the standard method or the
     302//  democratic method is implemented:
     303//
     304//  KStandard: This method looks for all pixels with an entry (photons)
     305//             that is three times bigger than the noise of the pixel
     306//             (default: 3 sigma, clean level 1)
     307//  AM 18/11/2002: now cut levels are proportional to the pixel area.
     308//  In this way the cut corresponds to a fixed
    98309//  phe-density (otherwise, it would bias the images).
    99310//
    100311//  Returns the maximum Pixel Id (used for ispixused in CleanStep2)
    101312//
    102 Int_t MImgCleanStd::CleanStep1()
     313Int_t MImgCleanStd::CleanStep1Std()
    103314{
    104315    const Int_t entries = fEvt->GetNumPixels();
     
    109320    // check the number of all pixels against the noise level and
    110321    // set them to 'unused' state if necessary
    111     //
     322    // 
    112323    for (Int_t i=0; i<entries; i++ )
    113324    {
     
    118329
    119330        const Int_t id = pix.GetPixId();
    120 
    121         const Double_t ratio = fCam->GetPixRatio(id);
     331        const Double_t ratio = TMath::Sqrt(fCam->GetPixRatio(id));
    122332
    123333        // COBB: '<=' to skip entry=noise=0
    124         if (entry*ratio <= fCleanLvl1*noise)
     334        if (entry <= fCleanLvl1 * noise / ratio)
    125335            pix.SetPixelUnused();
    126336
     
    128338            max = id;
    129339    }
     340
    130341    return max;
    131342}
     
    133344// --------------------------------------------------------------------------
    134345//
    135 //  check if the survived pixel have a neighbor, that also
    136 //  survived, otherwise set pixel to unused. (removes pixels without
    137 //  neighbors)
    138 //
    139 //  takes the maximum pixel id from CleanStep1 as an argument
     346//  NT 28/04/2003: now the option to use the standard method or the
     347//  democratic method is implemented:
     348//
     349//  "KDemocratic": this method looks for all pixels with an entry (photons)
     350//                 that is n times bigger than the noise of the mean of the
     351//                 inner pixels (default: 3 sigmabar, clean level 1)
     352//
     353//  Returns the maximum Pixel Id (used for ispixused in CleanStep2)
     354//
     355Int_t MImgCleanStd::CleanStep1Dem()
     356{
     357    const Int_t entries = fEvt->GetNumPixels();
     358
     359    Int_t max = entries;
     360
     361    //
     362    // check the number of all pixels against the noise level and
     363    // set them to 'unused' state if necessary
     364    //
     365    for (Int_t i=0; i<entries; i++ )
     366    {
     367        MCerPhotPix &pix = (*fEvt)[i];
     368
     369        const Float_t entry = pix.GetNumPhotons();
     370
     371        const Int_t id = pix.GetPixId();
     372
     373        const Double_t ratio = fCam->GetPixRatio(id);
     374
     375            // COBB: '<=' to skip entry=noise=0
     376
     377        if (entry <= fCleanLvl1 * fInnerNoise / ratio)
     378            pix.SetPixelUnused();
     379
     380        if (id>max)
     381            max = id;
     382    }
     383    return max;
     384}
     385
     386// --------------------------------------------------------------------------
     387// The first step of cleaning defines the CORE pixels. All the other pixels
     388// are set as UNUSED and belong to RING 0.               
     389// After this point, only the CORE pixels are set as USED, with RING     
     390// number 1.                                                             
     391// Returns the maximum Pixel Id (used for ispixused in CleanStep2)
     392//
     393Int_t MImgCleanStd::CleanStep1()
     394{
     395   switch (fCleaningMethod)
     396    {
     397    case kStandard:
     398        return CleanStep1Std();
     399    case kDemocratic:
     400        return CleanStep1Dem();
     401    }
     402
     403    return 0;
     404}
     405
     406// --------------------------------------------------------------------------
     407//
     408//  Check if the survived pixel have a neighbor, that also
     409//  survived, otherwise set pixel to unused (removes pixels without
     410//  neighbors).
     411//
     412//  Takes the maximum pixel id from CleanStep1 as an argument
    140413//
    141414void MImgCleanStd::CleanStep2(Int_t max)
     
    150423    //
    151424    Byte_t *ispixused = new Byte_t[max+1];
    152     memset(ispixused, 0, max+1);
    153425
    154426    for (Int_t i=0; i<entries; i++)
    155427    {
    156         MCerPhotPix &pix = (*fEvt)[i];
    157         ispixused[pix.GetPixId()] = pix.IsPixelUsed();
     428        const MCerPhotPix &pix = (*fEvt)[i];
     429        ispixused[pix.GetPixId()] = pix.IsPixelUsed() ? 1 : 0 ;
    158430    }
    159431
    160432    for (Int_t i=0; i<entries; i++)
    161433    {
    162         //
    163434        // get entry i from list
    164435        //
    165436        MCerPhotPix &pix = (*fEvt)[i];
    166437
    167         //
    168         // check if pixel is in use, if not goto next pixel in list
    169         //
    170 #if 0
    171         if (!pix.IsPixelUsed())
    172             continue;
    173 #endif
    174 
    175         //
    176         // get pixel id of this entry
     438        // get pixel id of this entry
    177439        //
    178440        const Int_t id = pix.GetPixId();
    179441
     442        // check if pixel is in use, if not goto next pixel in list
    180443        //
    181         // check for 'used' neighbors this pixel which
     444        if (ispixused[id] == 0)
     445            continue;
     446 
     447        // check for 'used' neighbors of this pixel
    182448        //
    183449        const MGeomPix &gpix  = (*fCam)[id];
    184450        const Int_t     nnmax = gpix.GetNumNeighbors();
    185451
    186 #if 0
    187         Bool_t cnt = kFALSE;
     452        Bool_t hasNeighbor = kFALSE;
     453
     454        //loop on the neighbors to check if they are used
     455        //
    188456        for (Int_t j=0; j<nnmax; j++)
    189457        {
    190458            const Int_t id2 = gpix.GetNeighbor(j);
    191459
    192             if (id2>max || !ispixused[id2])
    193                 continue;
    194 
    195             cnt = kTRUE;
    196             break;
     460            // when you find an used  neighbor, break the loop
     461            if (ispixused[id2] == 1  )
     462              {
     463                hasNeighbor = kTRUE ;
     464                break;
     465              }
    197466        }
    198         if (cnt)
    199             continue;
    200 #else
    201         Int_t cnt = 0;
    202         for (Int_t j=0; j<nnmax; j++)
    203         {
    204             const Int_t id2 = gpix.GetNeighbor(j);
    205 
    206             if (id2>max || !ispixused[id2])
    207                 continue;
    208 
    209             if (cnt++>nnmax-4)
    210                 break;
    211         }
    212         if (cnt==nnmax-2 && nnmax>=4)
    213         {
    214             pix.SetPixelUsed();
    215             continue;
    216         }
    217         if (cnt>0)
    218             continue;
    219 #endif
    220 
    221         //
    222         // check if no next neighbor has the state 'used'
    223         // set this pixel to 'unused', too.
    224         //
    225         pix.SetPixelUnused();
     467       
     468        if (hasNeighbor == kFALSE)
     469          pix.SetPixelUnused();
    226470    }
    227471
     
    246490//   a core neigbor it is declared as used.
    247491//
     492Bool_t MImgCleanStd::CleanStep3Std(const MCerPhotPix &pix)
     493{
     494    //
     495    // get pixel id of this entry
     496    //
     497    const Int_t id = pix.GetPixId();
     498
     499    //
     500    // check the num of photons against the noise level
     501    //
     502    const Float_t entry = pix.GetNumPhotons();
     503    const Float_t noise = pix.GetErrorPhot();
     504
     505    const Double_t ratio = TMath::Sqrt(fCam->GetPixRatio(id));
     506
     507    return (entry <= fCleanLvl2 * noise / ratio);
     508}
     509
     510// --------------------------------------------------------------------------
     511//
     512//   Look for the boundary pixels around the core pixels
     513//   if a pixel has more than 2.5 (clean level 2.5) sigmabar and
     514//   a core neighbor, it is declared as used.
     515//
     516Bool_t MImgCleanStd::CleanStep3Dem(const MCerPhotPix &pix)
     517{
     518    //
     519    // get pixel id of this entry
     520    //
     521    const Int_t id = pix.GetPixId();
     522
     523    //
     524    // check the num of photons against the noise level
     525    //
     526    const Float_t entry = pix.GetNumPhotons();
     527
     528    const Double_t ratio = fCam->GetPixRatio(id);
     529
     530    return (entry <= fCleanLvl2 * fInnerNoise / ratio);
     531}
     532
     533void MImgCleanStd::CleanStep3b(MCerPhotPix &pix)
     534{
     535    const Int_t id = pix.GetPixId();
     536
     537    //
     538    // check if the pixel's next neighbor is a core pixel.
     539    // if it is a core pixel set pixel state to: used.
     540    //
     541    MGeomPix   &gpix  = (*fCam)[id];
     542    const Int_t nnmax = gpix.GetNumNeighbors();
     543
     544    for (Int_t j=0; j<nnmax; j++)
     545    {
     546        const Int_t id2 = gpix.GetNeighbor(j);
     547
     548        if (!fEvt->GetPixById(id2) || !fEvt->IsPixelCore(id2))
     549          continue;
     550
     551        pix.SetPixelUsed();
     552           break;
     553    }
     554}
     555
     556// --------------------------------------------------------------------------
     557//
     558//   NT: Add option "rings": default value = 1.
     559//   Look n (n>1) times for the boundary pixels around the used pixels.
     560//   If a pixel has more than 2.5 (clean level 2.5) sigma,
     561//   it is declared as used.
     562//
     563//   If a value<2 for fCleanRings is used, no CleanStep4 is done.
     564//
     565void MImgCleanStd::CleanStep4(UShort_t r, MCerPhotPix &pix)
     566{
     567    //
     568    // check if the pixel's next neighbor is a used pixel.
     569    // if it is a used pixel set pixel state to: used,
     570    // and tell to which ring it belongs to.
     571    //
     572    const Int_t id = pix.GetPixId();
     573    MGeomPix  &gpix  = (*fCam)[id];
     574
     575    const Int_t nnmax = gpix.GetNumNeighbors();
     576
     577    for (Int_t j=0; j<nnmax; j++)
     578    {
     579        const Int_t id2 = gpix.GetNeighbor(j);
     580
     581        MCerPhotPix &npix = *fEvt->GetPixById(id2);
     582
     583        // FIXME!
     584        // Needed check to read CT1 data without having a Segmentation fault
     585        if (!fEvt->GetPixById(id2))
     586            continue;
     587
     588        if (!npix.IsPixelUsed() || npix.GetRing()>r-1 )
     589            continue;
     590
     591        pix.SetRing(r);
     592            break;
     593    }
     594}
     595
     596// --------------------------------------------------------------------------
     597//
     598//   Look for the boundary pixels around the core pixels
     599//   if a pixel has more than 2.5 (clean level 2.5) sigma, and
     600//   a core neigbor, it is declared as used.
     601//
    248602void MImgCleanStd::CleanStep3()
    249603{
    250604    const Int_t entries = fEvt->GetNumPixels();
    251605
    252     for (Int_t i=0; i<entries; i++)
    253     {
    254         //
    255         // get pixel as entry il from list
    256         //
    257         MCerPhotPix &pix = (*fEvt)[i];
    258 
    259         //
    260         // if pixel is a core pixel go to the next pixel
    261         //
    262         if (pix.IsPixelCore())
    263             continue;
    264 
    265         //
    266         // check the num of photons against the noise level
    267         //
    268         const Float_t entry = pix.GetNumPhotons();
    269         const Float_t noise = pix.GetErrorPhot();
    270 
    271         //
    272         // get pixel id of this entry
    273         //
    274         const Int_t id = pix.GetPixId();
    275 
    276         const Double_t ratio = TMath::Sqrt(fCam->GetPixRatio(id));
    277 
    278         if (entry*ratio <= fCleanLvl2*noise)
    279             continue;
    280 
    281         //
    282         // check if the pixel's next neighbor is a core pixel.
    283         // if it is a core pixel set pixel state to: used.
    284         //
    285         MGeomPix   &gpix  = (*fCam)[id];
    286         const Int_t nnmax = gpix.GetNumNeighbors();
    287 
    288         for (Int_t j=0; j<nnmax; j++)
     606    for (UShort_t r=1; r<fCleanRings+1; r++)
     607    {
     608        for (Int_t i=0; i<entries; i++)
    289609        {
    290             const Int_t id2 = gpix.GetNeighbor(j);
    291 
    292             if (!fEvt->IsPixelCore(id2))
     610            //
     611            // get pixel as entry il from list
     612            //
     613            MCerPhotPix &pix = (*fEvt)[i];
     614
     615            //
     616            // if pixel is a core pixel go to the next pixel
     617            //
     618            if (pix.IsPixelCore())
    293619                continue;
    294620
    295             pix.SetPixelUsed();
    296             break;
     621            switch (fCleaningMethod)
     622            {
     623            case kStandard:
     624                if (CleanStep3Std(pix))
     625                    continue;
     626                break;
     627            case kDemocratic:
     628                if (CleanStep3Dem(pix))
     629                    continue;
     630                break;
     631            }
     632
     633            if (r==1)
     634                CleanStep3b(pix);
     635            else
     636                CleanStep4(r, pix);
    297637        }
    298638    }
     
    301641// --------------------------------------------------------------------------
    302642//
    303 //  check if MEvtHeader exists in the Parameter list already.
     643//  Check if MEvtHeader exists in the Parameter list already.
    304644//  if not create one and add them to the list
    305645//
     
    320660    }
    321661
     662    if (fCleaningMethod != kDemocratic)
     663        return kTRUE;
     664
     665    fSgb = (MSigmabar*)pList->FindObject("MSigmabar");
     666    if (!fSgb)
     667    {
     668        *fLog << dbginf << "MSigmabar not found... aborting." << endl;
     669        return kFALSE;
     670    }
     671
    322672    return kTRUE;
    323673}
    324674
    325 
    326675// --------------------------------------------------------------------------
    327676//
     
    330679Bool_t MImgCleanStd::Process()
    331680{
     681    if (fSgb)
     682        fInnerNoise = fSgb->GetSigmabarInner();
     683
    332684    const Int_t max = CleanStep1();
    333685    CleanStep2(max);
     
    343695void MImgCleanStd::Print(Option_t *o) const
    344696{
    345     *fLog << GetDescriptor() << " initialized with noise level ";
    346     *fLog << fCleanLvl1 << " and " << fCleanLvl2 << endl;
    347 }
    348 
    349 // --------------------------------------------------------------------------
    350 //
    351 //  Craete two text entry fields, one for each cleaning level and a
     697    *fLog << all << GetDescriptor() << " using ";
     698    switch (fCleaningMethod)
     699    {
     700    case kDemocratic:
     701        *fLog << "democratic";
     702        break;
     703    case kStandard:
     704        *fLog << "standard";
     705        break;
     706    }
     707    *fLog << " cleaning initialized with noise level " << fCleanLvl1 << " and " << fCleanLvl2;
     708    *fLog << " (CleanRings=" << fCleanRings << ")" << endl;
     709}
     710
     711// --------------------------------------------------------------------------
     712//
     713//  Create two text entry fields, one for each cleaning level and a
    352714//  describing text line.
    353715//
    354716void MImgCleanStd::CreateGuiElements(MGGroupFrame *f)
    355717{
    356     //
    357     // Create a frame for line 3 and 4 to be able
    358     // to align entry field and label in one line
    359     //
    360     TGHorizontalFrame *f1 = new TGHorizontalFrame(f, 0, 0);
    361     TGHorizontalFrame *f2 = new TGHorizontalFrame(f, 0, 0);
    362 
    363     /*
    364      * --> use with root >=3.02 <--
    365      *
    366 
    367      TGNumberEntry *fNumEntry1 = new TGNumberEntry(frame, 3.0, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
    368      TGNumberEntry *fNumEntry2 = new TGNumberEntry(frame, 2.5, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
    369 
    370      */
    371     TGTextEntry *entry1 = new TGTextEntry(f1, "****", kImgCleanLvl1);
    372     TGTextEntry *entry2 = new TGTextEntry(f2, "****", kImgCleanLvl2);
    373 
    374     // --- doesn't work like expected (until root 3.02?) --- fNumEntry1->SetAlignment(kTextRight);
    375     // --- doesn't work like expected (until root 3.02?) --- fNumEntry2->SetAlignment(kTextRight);
    376 
    377     entry1->SetText("3.0");
    378     entry2->SetText("2.5");
    379 
    380     entry1->Associate(f);
    381     entry2->Associate(f);
    382 
    383     TGLabel *l1 = new TGLabel(f1, "Cleaning Level 1");
    384     TGLabel *l2 = new TGLabel(f2, "Cleaning Level 2");
    385 
    386     l1->SetTextJustify(kTextLeft);
    387     l2->SetTextJustify(kTextLeft);
    388 
    389     //
    390     // Align the text of the label centered, left in the row
    391     // with a left padding of 10
    392     //
    393     TGLayoutHints *laylabel = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 10);
    394     TGLayoutHints *layframe = new TGLayoutHints(kLHintsCenterY|kLHintsLeft,  5, 0, 10);
    395 
    396     //
    397     // Add one entry field and the corresponding label to each line
    398     //
    399     f1->AddFrame(entry1);
    400     f2->AddFrame(entry2);
    401 
    402     f1->AddFrame(l1, laylabel);
    403     f2->AddFrame(l2, laylabel);
    404 
    405     f->AddFrame(f1, layframe);
    406     f->AddFrame(f2, layframe);
    407 
    408     f->AddToList(entry1);
    409     f->AddToList(entry2);
    410     f->AddToList(l1);
    411     f->AddToList(l2);
    412     f->AddToList(laylabel);
    413     f->AddToList(layframe);
    414 }
    415 
    416 void MImgCleanStd::SetLvl1(Float_t lvl)
    417 {
    418     fCleanLvl1 = lvl;
    419     *fLog << inf << "Cleaning level 1 set to " << lvl << " sigma." << endl;
    420 }
    421 
    422 void MImgCleanStd::SetLvl2(Float_t lvl)
    423 {
    424     fCleanLvl2 = lvl;
    425     *fLog << inf << "Cleaning level 2 set to " << lvl << " sigma." << endl;
     718  //
     719  // Create a frame for line 3 and 4 to be able
     720  // to align entry field and label in one line
     721  //
     722  TGHorizontalFrame *f1 = new TGHorizontalFrame(f, 0, 0);
     723  TGHorizontalFrame *f2 = new TGHorizontalFrame(f, 0, 0);
     724 
     725  /*
     726   * --> use with root >=3.02 <--
     727   *
     728   
     729   TGNumberEntry *fNumEntry1 = new TGNumberEntry(frame, 3.0, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
     730   TGNumberEntry *fNumEntry2 = new TGNumberEntry(frame, 2.5, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
     731
     732  */
     733  TGTextEntry *entry1 = new TGTextEntry(f1, "****", kImgCleanLvl1);
     734  TGTextEntry *entry2 = new TGTextEntry(f2, "****", kImgCleanLvl2);
     735 
     736  // --- doesn't work like expected (until root 3.02?) --- fNumEntry1->SetAlignment(kTextRight);
     737  // --- doesn't work like expected (until root 3.02?) --- fNumEntry2->SetAlignment(kTextRight);
     738 
     739  entry1->SetText("3.0");
     740  entry2->SetText("2.5");
     741 
     742  entry1->Associate(f);
     743  entry2->Associate(f);
     744 
     745  TGLabel *l1 = new TGLabel(f1, "Cleaning Level 1");
     746  TGLabel *l2 = new TGLabel(f2, "Cleaning Level 2");
     747 
     748  l1->SetTextJustify(kTextLeft);
     749  l2->SetTextJustify(kTextLeft);
     750 
     751  //
     752  // Align the text of the label centered, left in the row
     753  // with a left padding of 10
     754  //
     755  TGLayoutHints *laylabel = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 10);
     756  TGLayoutHints *layframe = new TGLayoutHints(kLHintsCenterY|kLHintsLeft,  5, 0, 10);
     757 
     758  //
     759  // Add one entry field and the corresponding label to each line
     760  //
     761  f1->AddFrame(entry1);
     762  f2->AddFrame(entry2);
     763 
     764  f1->AddFrame(l1, laylabel);
     765  f2->AddFrame(l2, laylabel);
     766 
     767  f->AddFrame(f1, layframe);
     768  f->AddFrame(f2, layframe);
     769 
     770  f->AddToList(entry1);
     771  f->AddToList(entry2);
     772  f->AddToList(l1);
     773  f->AddToList(l2);
     774  f->AddToList(laylabel);
     775  f->AddToList(layframe);
    426776}
    427777
     
    445795    {
    446796    case kImgCleanLvl1:
    447         SetLvl1(lvl);
     797        fCleanLvl1 = lvl;
     798        *fLog << "Cleaning level 1 set to " << lvl << " sigma." << endl;
    448799        return kTRUE;
    449800
    450801    case kImgCleanLvl2:
    451         SetLvl2(lvl);
     802        fCleanLvl2 = lvl;
     803        *fLog << "Cleaning level 2 set to " << lvl << " sigma." << endl;
    452804        return kTRUE;
    453805    }
     
    475827    out << ");" << endl;
    476828}
    477 
    478 // --------------------------------------------------------------------------
    479 //
    480 // Read the setup from a TEnv:
    481 //   Float_t fCleanLvl1: CleaningLevel1
    482 //   Float_t fCleanLvl2: CleaningLevel2
    483 //
    484 Bool_t MImgCleanStd::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
    485 {
    486     Bool_t rc = kTRUE;
    487     if (!IsEnvDefined(env, prefix, "CleaningLevel1", print))
    488         rc = kFALSE;
    489     else
    490     {
    491         SetLvl1(GetEnvValue(env, prefix, "CleaningLevel1", fCleanLvl1));
    492         if (fCleanLvl1<0)
    493         {
    494             *fLog << err << "ERROR - Negative values for Cleaning Level 1 forbidden." << endl;
    495             return kERROR;
    496         }
    497     }
    498 
    499     if (!IsEnvDefined(env, prefix, "CleaningLevel2", print))
    500         rc = kFALSE;
    501     else
    502     {
    503         SetLvl2(GetEnvValue(env, prefix, "CleaningLevel2", fCleanLvl2));
    504         if (fCleanLvl2<0)
    505         {
    506             *fLog << err << "ERROR - Negative values for Cleaning Level 2 forbidden." << endl;
    507             return kERROR;
    508         }
    509     }
    510 
    511     return rc;
    512 }
  • trunk/MagicSoft/Mars/mimage/MImgCleanStd.h

    r1940 r2035  
    77
    88class MGeomCam;
     9class MSigmabar;
     10class MCerPhotPix;
    911class MCerPhotEvt;
    1012
     
    1315class MImgCleanStd : public MGTask
    1416{
     17public:
     18    typedef enum {
     19        kStandard,
     20        kDemocratic
     21    } CleaningMethod_t;
     22
    1523private:
    1624    const MGeomCam    *fCam;  //!
    1725          MCerPhotEvt *fEvt;  //!
     26          MSigmabar   *fSgb;  //!
     27
     28    CleaningMethod_t fCleaningMethod;
    1829
    1930    Float_t fCleanLvl1;
    2031    Float_t fCleanLvl2;
    2132
    22     void SetLvl1(Float_t lvl);
    23     void SetLvl2(Float_t lvl);
     33    UShort_t fCleanRings;
     34
     35    Float_t fInnerNoise;      //!
    2436
    2537    void CreateGuiElements(MGGroupFrame *f);
    2638    void StreamPrimitive(ofstream &out) const;
    2739
     40    Int_t  CleanStep1Dem();
     41    Int_t  CleanStep1Std();
     42    Bool_t CleanStep3Dem(const MCerPhotPix &pix);
     43    Bool_t CleanStep3Std(const MCerPhotPix &pix);
     44    void   CleanStep3b(MCerPhotPix &pix);
     45    void   CleanStep4(UShort_t r, MCerPhotPix &pix);
     46
    2847public:
    2948    MImgCleanStd(const Float_t lvl1=3.0, const Float_t lvl2=2.5,
    30                  const char *name=NULL, const char *title=NULL);
     49              const char *name=NULL, const char *title=NULL);
    3150
    3251    Int_t CleanStep1();
     
    4160    Float_t GetCleanLvl1() const { return fCleanLvl1; }
    4261    Float_t GetCleanLvl2() const { return fCleanLvl2; }
     62    UShort_t  GetCleanRings() const { return fCleanRings;}
     63
     64    void SetCleanRings(UShort_t r) { fCleanRings=r; }
     65    void SetMethod(CleaningMethod_t m) { fCleaningMethod = m; }
    4366
    4467    Bool_t ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2);
    45     Bool_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
    4668
    47     ClassDef(MImgCleanStd, 0)    // task doing a standard image cleaning
     69    ClassDef(MImgCleanStd, 2)    // task doing the image cleaning
    4870};
    4971
    5072#endif
    5173
     74
     75
     76
     77
     78
     79
     80
Note: See TracChangeset for help on using the changeset viewer.