Changeset 8646 for trunk


Ignore:
Timestamp:
07/26/07 12:13:00 (17 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r8644 r8646  
    1818
    1919                                                 -*-*- END OF LINE -*-*-
     20
     21 2007/07/26 Thomas Bretz
     22
     23   * mfilter/MFMagicCuts.cc:
     24     - implemented new disp-parametrization
     25     - implemented new ghostbusting
     26     - implemented possibility to use an external ghostbuster
     27     - implemented possibility to use an external disp calculator
     28
     29   * mimage/MImgCleanStd.[h,cc]:
     30     - implemented the cleaning in a recursive way. It should be
     31       a little bit faster
     32     - implemented post-cleaning using timing information
     33     - implemented new options how to treat the removed single
     34       core pixels
     35
     36   * mjobs/MJCut.cc:
     37     - write an external ghostbuster to the output file if available
     38
     39   * mjtrain/MJTrainEnergy.cc, mjtrain/MJTrainSeparation.cc:
     40     - set display name of MRanForestCalc as title instead of name
     41
     42   * mranforest/MRanForestCalc.cc:
     43     - set fTitle as eventloop name instead of fName
     44
     45
     46
     47 2007/07/25 Thomas Bretz
     48
     49   * scripts/merppupdate:
     50     - small fixes
     51
     52
    2053
    2154 2007/07/24 Thomas Bretz
  • trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc

    r8633 r8646  
    5151//
    5252//    c[0],  c[5], c[6], c[7],  c[8], c[9]:
    53 //          xi          = c[0]+c[6]*pow(leakage1, c[7]);
     53//          xi          = c[0] + c[8]*slope + c[9]*leak +
     54//                        (lgsize>c[10])*c[11]*(lgsize-c[10])^2;
    5455//          p           = xi*(1-width/length);
    5556//          sign1       = m3long-c[5]
    56 //          sign2       = (dist-c[9])*c[8]-slope
     57//          sign2       = (dist-c[7])*c[6]-slope
    5758//          disp        = sign1<0 ||sign2<0 ? -disp : disp
    58 //          antisign1   = antim3long-c[5]
    59 //          antisign2   = (antidist-c[9])*c[8]-antislope
    60 //          antidisp    = antisign1<0 || antisign2<0 ? -antidisp : antidisp
    6159//          thetasq     = disp^2 + dist^2 - 2*disp*dist*alpha
    62 //          antithetasq = antidisp^2 + antidist^2 - 2*antidisp*antidist*antialpha
     60//
     61//    And the values with respect to the antisource position respectively.
     62//
    6363//
    6464//    c[1]:
     
    7171//
    7272//    HadronnessCut:
    73 //    c[10], c[11]:
    74 //          had       <= c[10]
    75 //          10^lgsize >= c[11]
     73//    c[13], c[14]:
     74//          had       <= c[13]
     75//          10^lgsize >= c[14]
    7676//
    7777//
     
    8181//   HadronnessCut:
    8282//    - none/nocut: No area cut
    83 //    - area:       Area cut
    84 //    - all:        same as area
     83//    - area:       Area cut <default>
     84//    - hadronness: Cut in size and hadronness (c[10], c[11])
     85//    - all:        area + hadronness
    8586//
    8687//   ThetaCut:
    87 //    - none/nocut: no theta cut
     88//    - none/nocut: no theta cut <default>
    8889//    - on:         cut in theta only
    8990//    - off:        anti-cut in anti-theta only
    9091//    - wobble:     same as on|off (both)
    9192//
     93//   CalcDisp:
     94//    - yes:        Disp is calculated as defined above <default>
     95//    - no:         abs(Disp.fVal) from the parameter list is used instead
     96//
     97//   CalcGhostbuster:
     98//    - yes:        The ghostbuster is calculated as defined above <default>
     99//    - no:         Ghostbuster.fVal<c[12] is used as ghostbusting condition
    92100//
    93101//  Input Container:
     
    99107//   MHillasSrc
    100108//   [MHillasSrcAnti [MHillasSrc]]
     109//   [Disp [MParameterD]]
     110//   [Ghostbuster [MParameterD]]
    101111//
    102112//  Output:
     
    138148MFMagicCuts::MFMagicCuts(const char *name, const char *title)
    139149    : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fNewImgPar(0),
    140     fThetaSq(0), fDisp(0), fHadronness(0), fMatrix(0), fVariables(30),
    141     fThetaCut(kNone), fHadronnessCut(kArea)
     150    fThetaSq(0), fDisp(0), fHadronness(0), fMatrix(0), fVariables(20),
     151    fThetaCut(kNone), fHadronnessCut(kArea), fCalcDisp(kTRUE),
     152    fCalcGhostbuster(kTRUE)
    142153{
    143154    fName  = name  ? name  : "MFMagicCuts";
     
    159170    AddToBranchList("MNewImagePar.fLeakage1");
    160171    AddToBranchList("Hadronness.fVal");
    161 /*
    162     fVariables[0] =  1.42547;      // Xi
    163     fVariables[1] =  0.233773;     // Theta^2
    164     fVariables[2] =  0.2539460;    // Area[0]
    165     fVariables[3] =  5.2149800;    // Area[1]
    166     fVariables[4] =  0.1139130;    // Area[2]
    167     fVariables[5] = -0.0889;       // M3long
    168     fVariables[6] =  0.18;         // Xi(theta)
    169     fVariables[7] =  0.18;         // Xi(theta)
    170     */
     172    AddToBranchList("Disp.fVal");
     173    AddToBranchList("Ghostbuster.fVal");
    171174}
    172175
     
    190193    if (!fThetaSq)
    191194        return kFALSE;
    192     fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
     195
     196    if (!fCalcDisp)
     197        fDisp = (MParameterD*)pList->FindObject("Disp", "MParameterD");
     198    else
     199        fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
    193200    if (!fDisp)
    194         return kFALSE;
     201    {
     202        *fLog << err << "Disp [MParameterD] not found... aborting." << endl;
     203        return kFALSE;
     204    }
     205
     206    if (!fCalcGhostbuster)
     207        fGhostbuster = (MParameterD*)pList->FindObject("Ghostbuster", "MParameterD");
     208    else
     209        fGhostbuster = (MParameterD*)pList->FindCreateObj("MParameterD", "Ghostbuster");
     210    if (!fGhostbuster)
     211    {
     212        *fLog << err << "Ghostbuster [MParameterD] not found... aborting." << endl;
     213        return kFALSE;
     214    }
    195215
    196216    // propagate Theta cut to the parameter list
     
    297317    fMap[kEM3Long]  = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
    298318
    299     fMap[kESrcSign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)");
     319    fMap[kESign]    = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)");
    300320
    301321    fMap[kESlope]   = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
     322
     323    if (!fCalcDisp)
     324        fMap[kEDisp] = fMatrix->AddColumn("abs(Disp.fVal)");
     325
     326    if (!fCalcGhostbuster)
     327        fMap[kEGhostbuster] = fMatrix->AddColumn("Ghostbuster.fVal");
    302328
    303329    if (fThetaCut&kOff)
     
    321347// If par!=NULL p is stored in this MParameterD
    322348//
    323 Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t wl, Double_t d, Double_t a) const
     349Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t d, Double_t a) const
    324350{
    325351    // wl = width/l [1]
     
    338364}
    339365
    340 // ---------------------------------------------------------------------------
    341 //
    342 // Evaluate dynamical magic-cuts
    343 //
    344 Int_t MFMagicCuts::Process()
    345 {
     366// --------------------------------------------------------------------------
     367//
     368// Return abs(Disp.fVal) if disp calculation is switched off.
     369// Otherwise return (c0+c6*leak^c7)*(1-width/length)
     370//
     371Double_t MFMagicCuts::GetDisp(Double_t slope, Double_t lgsize) const
     372{
     373    if (!fCalcDisp)
     374        return fMatrix ? GetVal(kEDisp) : TMath::Abs(fDisp->GetVal());
     375
    346376    // Get some variables
    347     const Double_t wdivl  = fMatrix ? GetVal(kEWdivL)   : fHil->GetWidth()/fHil->GetLength();
    348     const Double_t lgsize = fMatrix ? GetVal(kESize)    : TMath::Log10(fHil->GetSize());
    349     const Double_t leak   = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1);
     377    const Double_t wdivl = fMatrix ? GetVal(kEWdivL)   : fHil->GetWidth()/fHil->GetLength();
     378    const Double_t leak  = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1);
    350379
    351380    // For simplicity
    352381    const Double_t *c = fVariables.GetArray();
    353382
     383    // As rule for root or MDataPhrase:
     384    //   ((M[3]>[3])*[4]*(M[3]-[3])^2 + [2]*M[2] + [1]*M[1] + [0])*M[0]
     385    //
     386    Double_t xi = c[0] + c[8]*slope + c[9]*leak;
     387    if (lgsize>c[10])
     388        xi += (lgsize-c[10])*(lgsize-c[10])*c[11];
     389
     390    const Double_t disp = xi*(1-wdivl);
     391
     392    return disp;
     393}
     394
     395Bool_t MFMagicCuts::IsGhost(Double_t m3long, Double_t slope, Double_t dist) const
     396{
     397    // For simplicity
     398    const Double_t *c = fVariables.GetArray();
     399
     400    if (!fCalcGhostbuster)
     401        return (fMatrix ? GetVal(kEGhostbuster) : fGhostbuster->GetVal())<c[12];
     402
     403    // Do Ghostbusting with slope and m3l
     404    const Bool_t sign1 = m3long < c[5];
     405    const Bool_t sign2 = slope  > (dist-c[7])*c[6];
     406
     407    return sign1 || sign2;
     408}
     409
     410// ---------------------------------------------------------------------------
     411//
     412// Evaluate dynamical magic-cuts
     413//
     414Int_t MFMagicCuts::Process()
     415{
     416    // For simplicity
     417    const Double_t *c = fVariables.GetArray();
     418
     419    // Default if we return before the end
     420    fResult = kFALSE;
     421
    354422    // Value for Theta cut (Disp parametrization)
    355423    const Double_t cut  = GetThetaSqCut();
    356     const Double_t xi   = c[0]+c[6]*pow(leak, c[7]);
    357     const Double_t disp = xi*(1-wdivl);
    358 
    359     // Default if we return before the end
    360     fResult = kFALSE;
    361424
    362425    // ------------------- Most efficient cut -----------------------
    363426    // ---------------------- Theta cut ON --------------------------
    364     const Double_t dist   = fMatrix ? GetVal(kEDist)    : fHilSrc->GetDist()*fMm2Deg;
    365     const Double_t alpha  = fMatrix ? GetVal(kEAlpha)   : fHilSrc->GetAlpha();
    366     const Double_t sign   = fMatrix ? GetVal(kESrcSign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha());
    367     const Double_t m3long = fMatrix ? GetVal(kEM3Long)  : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
    368     const Double_t slope  = fMatrix ? GetVal(kESlope)   : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
    369 
    370     // Do Ghostbusting with slope and m3l
    371     const Double_t sign1 = (dist-c[9])*c[8]-slope;
    372     const Double_t sign2 = m3long-c[5];
     427    const Double_t dist   = fMatrix ? GetVal(kEDist)   : fHilSrc->GetDist()*fMm2Deg;
     428    const Double_t alpha  = fMatrix ? GetVal(kEAlpha)  : fHilSrc->GetAlpha();
     429    const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
     430    const Double_t slope  = fMatrix ? GetVal(kESlope)  : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
     431    const Double_t sign   = fMatrix ? GetVal(kESign)   : MMath::Sgn(fHilSrc->GetCosDeltaAlpha());
     432    const Double_t lgsize = fMatrix ? GetVal(kESize)   : TMath::Log10(fHil->GetSize());
    373433
    374434    // Calculate disp including sign
    375     const Double_t p = sign1<0 || sign2<0 ? -disp : disp;
    376 
    377     // Align disp along source direction depending on third moment
    378     //const Double_t p = TMath::Sign(disp, m3long-c[5]);
     435    const Double_t disp  = GetDisp(slope, lgsize);
     436    const Double_t ghost = IsGhost(m3long, slope, dist);
     437
     438    const Double_t p = ghost ? -disp : disp;
    379439
    380440    // Align sign of disp along shower axis like m3long
     
    382442
    383443    // Calculate corresponding Theta^2
    384     const Double_t thetasq = GetThetaSq(p, wdivl, dist, alpha);
     444    const Double_t thetasq = GetThetaSq(p, dist, alpha);
    385445    fThetaSq->SetVal(thetasq);
    386446
     
    405465    {
    406466        const Double_t had = fMatrix ? GetVal(kEHadronness) : fHadronness->GetVal();
    407         if (had>c[10])
     467        if (had>c[13])
    408468            return kTRUE;
    409469
    410         if (TMath::Power(10, lgsize)<c[11])
     470        if (TMath::Power(10, lgsize)<c[14])
    411471            return kTRUE;
    412472    }
     
    422482        const Double_t slopeanti   = fMatrix ? GetVal(kESlopeAnti)  : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
    423483
    424         // Do Ghostbusting with slope and m3l
    425         const Double_t sign3 = (distanti-c[9])*c[8]-slopeanti;
    426         const Double_t sign4 = m3lanti-c[5];
    427         const Double_t panti = sign3<0 || sign4<0 ? -disp : disp;
     484        const Double_t dispanti  = GetDisp(slopeanti, lgsize);
     485        const Double_t ghostanti = IsGhost(m3lanti, slopeanti, lgsize);
     486
     487        const Double_t panti = ghostanti ? -dispanti : dispanti;
    428488
    429489        // Align disp along source direction depending on third moment
    430         //const Double_t panti       = TMath::Sign(disp, m3lanti-c[5]);
    431         const Double_t thetasqanti = GetThetaSq(panti, wdivl, distanti, alphaanti);
     490        const Double_t thetasqanti = GetThetaSq(panti, distanti, alphaanti);
    432491
    433492        if (thetasqanti<cut)
     
    494553        break;
    495554    }
     555    if (fCalcDisp)
     556        *fLog << "Disp is calculated from c0, c7, c8." << endl;
     557    else
     558        *fLog << "Disp.fVal from the parameter list is used as disp." << endl;
     559    if (fCalcGhostbuster)
     560        *fLog << "Ghostbusting is calculated from c5, c6, c7." << endl;
     561    else
     562        *fLog << "Ghostbuster.fVal from the parameter list is used for ghostbusting." << endl;
     563
    496564    *fLog << "Filter is" << (IsInverted()?"":" not") << " inverted." << endl;
    497565
     
    568636        rc = kTRUE;
    569637    }
     638
    570639    if (IsEnvDefined(env, prefix, "HadronnessCut", print))
    571640    {
     
    586655    }
    587656
     657    if (IsEnvDefined(env, prefix, "CalcDisp", print))
     658    {
     659        fCalcDisp = GetEnvValue(env, prefix, "CalcDisp", fCalcDisp);
     660        rc = kTRUE;
     661    }
     662
     663    if (IsEnvDefined(env, prefix, "CalcGhostbuster", print))
     664    {
     665        fCalcGhostbuster = GetEnvValue(env, prefix, "CalcGhostbuster", fCalcGhostbuster);
     666        rc = kTRUE;
     667    }
     668
    588669    if (IsEnvDefined(env, prefix, "File", print))
    589670    {
     
    605686    }
    606687    return rc;
    607     //return kTRUE; // means: can use default values
    608     //return rc;  // means: require something in resource file
    609 }
     688}
  • trunk/MagicSoft/Mars/mfilter/MFMagicCuts.h

    r8619 r8646  
    4343    enum {
    4444        kESize, kEAlpha, kEAlphaAnti, kEArea, kEDist, kEDistAnti,
    45         kEM3Long, kEM3LongAnti, kEWdivL, kELeakage, kESrcSign,
    46         kESlope, kESlopeAnti, kEHadronness,
     45        kEM3Long, kEM3LongAnti, kEWdivL, kELeakage, kESlope,
     46        kESlopeAnti, kEHadronness, kESign, kEDisp, kEGhostbuster,
    4747        kLastElement
    4848    };
     
    5555    MParameterD    *fThetaSq;           //! Pointer to MParameterD container called ThetaSq
    5656    MParameterD    *fDisp;              //! Pointer to MParameterD container called Disp
     57    MParameterD    *fGhostbuster;       //! Pointer to MParameterD container called Ghostbuster
    5758    MParameterD    *fHadronness;        //! Pointer to MParameterD container called Hadronness
    5859
     
    6768    ThetaCut_t      fThetaCut;          // Which kind of theta cut should be evaluated
    6869    HadronnessCut_t fHadronnessCut;     // Which kind of hadronness cut should be evaluated
     70    Bool_t          fCalcDisp;          // Should we use Disp from the parameterlist?
     71    Bool_t          fCalcGhostbuster;   // Should we use Ghostbuster from the parameterlist?
    6972
    7073    // MTask
     
    7679
    7780    // MFMagicCuts
     81    Double_t GetDisp(Double_t slope, Double_t lgsize) const;
     82    Bool_t   IsGhost(Double_t m3long, Double_t slope, Double_t dist) const;
    7883    Double_t GetVal(Int_t i) const;
    7984    TString  GetParam(Int_t i) const;
    80     Double_t GetThetaSq(Double_t p, Double_t wl, Double_t d, Double_t a) const;
     85    Double_t GetThetaSq(Double_t p, Double_t d, Double_t a) const;
    8186
    8287public:
     
    8994    void   SetThetaCut(ThetaCut_t c) { fThetaCut=c; }
    9095    void   SetHadronnessCut(HadronnessCut_t c) { fHadronnessCut=c; }
     96    void   SetCalcDisp(Bool_t b=kTRUE) { fCalcDisp=b; }
     97    void   SetCalcGhostbuster(Bool_t b=kTRUE) { fCalcGhostbuster=b; }
    9198
    9299    // MFMagicCuts
  • trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc

    r8554 r8646  
    2121!   Author(s): Stefan Ruegamer, 03/2006 <mailto:snruegam@astro.uni-wuerzburg.de>
    2222!
    23 !   Copyright: MAGIC Software Development, 2000-2006
     23!   Copyright: MAGIC Software Development, 2000-2007
    2424!
    2525!
     
    279279//
    280280//
     281//  Class Version 4:
     282//  ----------------
     283//   + Float_t  fTimeLvl2;
     284//   - Bool_t   fKeepSinglePixels;
     285//   + Bool_t   fKeepIsolatedPixels;
     286//   + Int_t    fRecoverIsolatedPixels;
     287//
     288//
    281289//  Input Containers:
    282290//   MGeomCam
     
    302310#include "MLogManip.h"
    303311
     312#include "MArrayI.h"
    304313#include "MParList.h"
    305314#include "MCameraData.h"
     
    339348MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2,
    340349                           const char *name, const char *title)
    341     : fCleaningMethod(kStandard), fCleanLvl1(lvl1),
    342     fCleanLvl2(lvl2), fCleanRings(1), fKeepSinglePixels(kFALSE),
    343     fNamePedPhotCam(gsNamePedPhotCam), fNameGeomCam(gsNameGeomCam),
    344     fNameSignalCam(gsNameSignalCam)
     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)
    345355{
    346356    fName  = name  ? name  : gsDefName.Data();
    347357    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    }
    348376}
    349377
     
    367395//
    368396//
    369 void MImgCleanStd::CleanStep1()
    370 {
     397Bool_t MImgCleanStd::HasCoreNeighbors(const MGeomPix &gpix) const
     398{
     399//    if (fKeepIsolatedPixels)
     400//        return kTRUE;
     401
     402#ifdef DEBUG
    371403    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
    372511
    373512    //
     
    378517    for (UInt_t idx=0; idx<npixevt; idx++)
    379518    {
    380         // The default for pixels is "used" set by
    381         // MParContainer::Reset before processing
    382         if (data[idx]>fCleanLvl1)
    383             continue;
    384 
    385         // Setting a pixel to unused if it is unmapped would overwrite
    386         // the unmapped-status. Therefor this pixels are excluded.
    387519        MSignalPix &pix = (*fEvt)[idx];
    388         if (!pix.IsPixelUnmapped())
    389             pix.SetPixelUnused();
    390     }
    391 }
    392 
    393 // --------------------------------------------------------------------------
    394 //
    395 //  Check if the survived pixel have a neighbor, that also
    396 //  survived. Set all single pixels Unused if !fKeepSinglePixels. Now we
    397 //  declare all pixels that survived previous CleanSteps as CorePixels.
    398 //  Return number of single pixels, and there cumulative size in size.
    399 //
    400 Short_t MImgCleanStd::CleanStep2(Float_t &size)
    401 {
    402     Short_t n=0;
    403     size = 0;
    404 
     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    //
    405588    const UInt_t npixevt = fEvt->GetNumPixels();
    406589    for (UInt_t idx=0; idx<npixevt; idx++)
    407590    {
    408         // Exclude all unused (this includes all unmapped) pixels
    409591        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];
    410671        if (!pix.IsPixelUsed())
    411672            continue;
    412  
    413         // check for 'used' neighbors of this pixel
    414         const MGeomPix &gpix  = (*fCam)[idx];
    415         const Int_t     nnmax = gpix.GetNumNeighbors();
    416 
    417         Bool_t hasNeighbor = kFALSE;
    418 
    419         //loop on the neighbors to check if they are used
    420         for (Int_t j=0; j<nnmax; j++)
     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++)
    421682        {
    422             const Int_t idx2 = gpix.GetNeighbor(j);
    423 
    424             // when you find an used neighbor (this excludes unused
    425             // and unmapped pixels) break the loop
    426             if ((*fEvt)[idx2].IsPixelUsed())
    427             {
    428                 hasNeighbor = kTRUE;
     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)
    429698                break;
    430             }
    431699        }
    432700
    433         // If the pixel has at least one core-neighbor
    434         // go on with the next pixel
    435         if (hasNeighbor)
    436             continue;
    437 
    438         // If the pixel has no neighbors and the single pixels
    439         // should not be kept turn the used- into an unused-status
    440         if (!fKeepSinglePixels)
    441             pix.SetPixelUnused();
    442 
    443         // count size and number of single core-pixels
    444         size += pix.GetNumPhotons();
    445         n++;
    446     }
    447 
    448     // Now turn the used-status into the core-status
    449     // (FIXME: A more intelligent handling of used/core in clean step1/2
    450     //         would make this loop obsolete!)
    451     for (UInt_t idx=0; idx<npixevt; idx++)
    452     {
    453         MSignalPix &pix = (*fEvt)[idx];
    454         pix.SetPixelCore(pix.IsPixelUsed());
    455     }
    456 
    457     return n;
    458 }
    459 
    460 void MImgCleanStd::CleanStep3b(Int_t idx)
    461 {
    462     MSignalPix &pix = (*fEvt)[idx];
    463 
    464     //
    465     // check if the pixel's next neighbor is a core pixel.
    466     // if it is a core pixel set pixel state to: used.
    467     //
    468     MGeomPix   &gpix  = (*fCam)[idx];
    469     const Int_t nnmax = gpix.GetNumNeighbors();
    470 
    471     for (Int_t j=0; j<nnmax; j++)
    472     {
    473         const Int_t idx2 = gpix.GetNeighbor(j);
    474 
    475         // Check if the neighbor pixel is a core pixel. (Rem: Unampped
    476         // pixels are never assigned the core-pixel status)
    477         if (!(*fEvt)[idx2].IsPixelCore())
    478             continue;
    479 
    480         pix.SetPixelUsed();
    481         break;
    482     }
    483 }
    484 
    485 // --------------------------------------------------------------------------
    486 //
    487 //   NT: Add option "rings": default value = 1.
    488 //   Look n (n>1) times for the boundary pixels around the used pixels.
    489 //   If a pixel has more than 2.5 (clean level 2.5) sigma,
    490 //   it is declared as used.
    491 //
    492 //   If a value<2 for fCleanRings is used, no CleanStep4 is done.
    493 //
    494 void MImgCleanStd::CleanStep4(UShort_t r, Int_t idx)
    495 {
    496     MSignalPix &pix = (*fEvt)[idx];
    497 
    498     //
    499     // Skip events that have already a defined status;
    500     //
    501     if (pix.GetRing() != 0)
     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)
    502721        return;
    503722
    504     //
    505     // check if the pixel's next neighbor is a used pixel.
    506     // if it is a used pixel set pixel state to: used,
    507     // and tell to which ring it belongs to.
    508     //
    509     MGeomPix  &gpix = (*fCam)[idx];
    510 
    511     const Int_t nnmax = gpix.GetNumNeighbors();
    512 
    513     for (Int_t j=0; j<nnmax; j++)
    514     {
    515         const Int_t idx2 = gpix.GetNeighbor(j);
    516 
    517         const MSignalPix &npix = (*fEvt)[idx2];
    518         if (!npix.IsPixelUsed() || npix.GetRing()>r-1 )
    519             continue;
    520 
    521         pix.SetRing(r);
    522         break;
    523     }
    524 }
    525 
    526 // --------------------------------------------------------------------------
    527 //
    528 //  Look for the boundary pixels around the core pixels
    529 //  if a pixel has more than 2.5 (clean level 2.5) sigma, and
    530 //  a core neighbor, it is declared as used.
    531 //
    532 void MImgCleanStd::CleanStep3()
    533 {
    534     const TArrayD &data = fData->GetData();
    535 
    536     for (UShort_t r=1; r<fCleanRings+1; r++)
    537     {
    538         // Loop over all pixels
    539         const UInt_t npixevt = fEvt->GetNumPixels();
    540         for (UInt_t idx=0; idx<npixevt; idx++)
    541         {
    542             MSignalPix &pix = (*fEvt)[idx];
    543 
    544             //
    545             // if pixel is a core pixel or unmapped, go to the next pixel
    546             //
    547             if (pix.IsPixelCore() || pix.IsPixelUnmapped())
    548                 continue;
    549 
    550             if (data[idx] <= fCleanLvl2)
    551                 continue;
    552 
    553             if (r==1)
    554                 CleanStep3b(idx);
    555             else
    556                 CleanStep4(r, idx);
    557         }
    558     }
     723    if (fPostCleanType&2)
     724        CleanTime(2, fTimeLvl2);
     725
     726    if (fPostCleanType&1)
     727        CleanTime(1, fTimeLvl1);
    559728}
    560729
     
    572741        return kFALSE;
    573742    }
    574 
    575743    fEvt = (MSignalCam*)pList->FindObject(AddSerialNumber(fNameSignalCam), "MSignalCam");
    576744    if (!fEvt)
     
    595763        return kFALSE;
    596764
     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
    597794    Print();
    598795
     
    631828
    632829#ifdef DEBUG
    633     *fLog << all << "CleanStep 1" << endl;
     830    *fLog << all << "ResetCleaning" << endl;
    634831#endif
    635     CleanStep1();
    636 
     832    ResetCleaning();
    637833
    638834#ifdef DEBUG
    639     *fLog << all << "CleanStep 2" << endl;
     835    *fLog << all << "DoCleaning" << endl;
    640836#endif
    641837    Float_t size;
    642     const Short_t n = CleanStep2(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
    643858    fEvt->SetSinglePixels(n, size);
    644859
    645     // For speed reasons skip the rest of the cleaning if no
    646     // action will be taken!
    647     if (fCleanLvl1>fCleanLvl2)
    648     {
    649860#ifdef DEBUG
    650         *fLog << all << "CleanStep 3" << endl;
    651 #endif
    652         CleanStep3();
    653     }
    654 
    655 #ifdef DEBUG
    656     *fLog << all << "Calc Islands" << endl;
     861    *fLog << all << "CalcIslands" << endl;
    657862#endif
    658863    // Takes roughly 10% of the time
     
    698903    *fLog << "initialized with level " << fCleanLvl1 << " and " << fCleanLvl2;
    699904    *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    }
    700927
    701928    if (fCleaningMethod!=kAbsolute && fCleaningMethod!=kTime)
     
    8481075    if (gsNameSignalCam!=fNameSignalCam)
    8491076        out << "   " << GetUniqueName() << ".SetNameSignalCam(\"" << fNameSignalCam << "\");" << endl;
    850     if (fKeepSinglePixels)
    851         out << "   " << GetUniqueName() << ".SetKeepSinglePixels();" << endl;
     1077    if (fKeepIsolatedPixels)
     1078        out << "   " << GetUniqueName() << ".SetKeepIsolatedPixels();" << endl;
     1079    if (fRecoverIsolatedPixels)
     1080        out << "   " << GetUniqueName() << ".SetRecoverIsolatedPixels(" << fRecoverIsolatedPixels << ");" << endl;
    8521081
    8531082}
     
    8601089//   MImgCleanStd.CleanMethod: Standard, Scaled, Democratic, Probability, Absolute
    8611090//   MImgCleanStd.CleanRings:  1
    862 //   MImgCleanStd.KeepSinglePixels: yes, no
     1091//   MImgCleanStd.KeepIsolatedPixels: yes, no
    8631092//
    8641093Int_t MImgCleanStd::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
     
    8701099        SetCleanRings(GetEnvValue(env, prefix, "CleanRings", fCleanRings));
    8711100    }
     1101    if (IsEnvDefined(env, prefix, "CleanLevel0", print))
     1102    {
     1103        rc = kTRUE;
     1104        fCleanLvl0 = GetEnvValue(env, prefix, "CleanLevel0", fCleanLvl0);
     1105    }
    8721106    if (IsEnvDefined(env, prefix, "CleanLevel1", print))
    8731107    {
     
    8801114        fCleanLvl2 = GetEnvValue(env, prefix, "CleanLevel2", fCleanLvl2);
    8811115    }
    882     if (IsEnvDefined(env, prefix, "KeepSinglePixels", print))
     1116    if (IsEnvDefined(env, prefix, "TimeLevel1", print))
    8831117    {
    8841118        rc = kTRUE;
    885         fKeepSinglePixels = GetEnvValue(env, prefix, "KeepSinglePixels", fKeepSinglePixels);
     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);
    8861140    }
    8871141
  • trunk/MagicSoft/Mars/mimage/MImgCleanStd.h

    r7804 r8646  
    77
    88class MGeomCam;
    9 class MSigmabar;
     9class MGeomPix;
    1010class MSignalCam;
    1111class MPedPhotCam;
    1212class MArrivalTime;
    1313class MCameraData;
     14//class MRawRunHeader;
    1415
    1516class MGGroupFrame;
     
    3233    static const TString gsNameSignalCam;  // default name of the 'MSignalCam' container
    3334
    34     const MGeomCam     *fCam;  //!
    35           MSignalCam   *fEvt;  //!
    36           MPedPhotCam  *fPed;  //!
    37           MCameraData  *fData; //!
     35    const MGeomCam      *fCam;  //!
     36          MSignalCam    *fEvt;  //!
     37          MPedPhotCam   *fPed;  //!
     38          MCameraData   *fData; //!
     39//          MRawRunHeader *fHeader; //!
    3840
    3941    CleaningMethod_t fCleaningMethod;
    4042
     43    Float_t  fCleanLvl0;
    4144    Float_t  fCleanLvl1;
    4245    Float_t  fCleanLvl2;
    4346
     47    Float_t  fTimeLvl1;
     48    Float_t  fTimeLvl2;
     49
    4450    UShort_t fCleanRings;
    45     Bool_t   fKeepSinglePixels;
     51    Bool_t   fKeepIsolatedPixels;
     52    Int_t    fRecoverIsolatedPixels;
     53    Int_t    fPostCleanType;
    4654
    4755    TString  fNamePedPhotCam; // name of the 'MPedPhotCam' container
     
    5058
    5159    // MImgCleanStd
    52     void    CleanStep1();
    53     Short_t CleanStep2(Float_t &size);
    54     void    CleanStep3();
    55     void    CleanStep3b(Int_t idx);
    56     void    CleanStep4(UShort_t r, Int_t idx);
     60    Bool_t HasCoreNeighbors(const MGeomPix &gpix) const;
     61    Bool_t HasUsedNeighbors(const MGeomPix &gpix) const;
     62    void   SetUsedNeighbors(const MGeomPix &gpix, Int_t r=1) const;
     63    Int_t  DoCleaning(Float_t &size) const;
     64    void   ResetCleaning() const;
     65    Int_t  RecoverIsolatedPixels(Float_t &size) const;
     66    void   CleanTime(Int_t n, Double_t lvl) const;
     67
     68    void CleanStepTime() const;
    5769
    5870    // MGTask, MTask, MParContainer
     
    7082    void Print(Option_t *o="") const;
    7183
    72     Float_t  GetCleanLvl1() const { return fCleanLvl1; }
    73     Float_t  GetCleanLvl2() const { return fCleanLvl2; }
     84    Float_t GetCleanLvl0() const { return fCleanLvl0; }
     85    Float_t GetCleanLvl1() const { return fCleanLvl1; }
     86    Float_t GetCleanLvl2() const { return fCleanLvl2; }
     87
     88    Float_t GetTimeLvl1() const { return fTimeLvl1; }
     89    Float_t GetTimeLvl2() const { return fTimeLvl2; }
     90
     91    void SetCleanLvl0(Float_t lvl) { fCleanLvl0=lvl; }
     92    void SetCleanLvl1(Float_t lvl) { fCleanLvl1=lvl; }
     93    void SetCleanLvl2(Float_t lvl) { fCleanLvl2=lvl; }
     94
     95    void SetTimeLvl1(Float_t lvl) { fTimeLvl1=lvl; }
     96    void SetTimeLvl2(Float_t lvl) { fTimeLvl2=lvl; }
     97
     98    void SetCleanRings(UShort_t r) { fCleanRings=r; }
    7499    UShort_t GetCleanRings() const { return fCleanRings;}
    75100
    76     void SetCleanRings(UShort_t r) { if(r==0) r=1; fCleanRings=r; }
    77101    void SetMethod(CleaningMethod_t m) { fCleaningMethod = m; }
    78     void SetKeepSinglePixels(Bool_t b=kTRUE) { fKeepSinglePixels=b; }
     102    void SetKeepIsolatedPixels(Bool_t b=kTRUE) { fKeepIsolatedPixels=b; }
     103    void SetRecoverIsolatedPixels(Int_t n=-1) { fRecoverIsolatedPixels=n; }
    79104
    80105    Bool_t ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2);
     
    84109    void SetNameGeomCam(const char *name)     { fNameGeomCam = name; }
    85110
    86     ClassDef(MImgCleanStd, 3)    // task doing the image cleaning
     111    ClassDef(MImgCleanStd, 4)    // task doing the image cleaning
    87112};
    88113
  • trunk/MagicSoft/Mars/mjobs/MJCut.cc

    r8644 r8646  
    388388    write->AddContainer("OpticalAxis",    "Events", kFALSE);
    389389    write->AddContainer("Disp",           "Events", kFALSE);
     390    write->AddContainer("Ghostbuster",    "Events", kFALSE);
    390391    write->AddContainer("MEnergyEst",     "Events", kFALSE);
    391392    write->AddContainer("MTime",          "Events", kFALSE);
  • trunk/MagicSoft/Mars/mjtrain/MJTrainEnergy.cc

    r8644 r8646  
    145145
    146146    // ------------------------ Train RF --------------------------
    147     MRanForestCalc rf(fTitle);
     147    MRanForestCalc rf("TrainEnergy", fTitle);
    148148    rf.SetNumTrees(fNumTrees);
    149149    rf.SetNdSize(fNdSize);
  • trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc

    r8644 r8646  
    802802    // ------------------------ Train RF --------------------------
    803803
    804     MRanForestCalc rf;
     804    MRanForestCalc rf("TrainSeparation", fTitle);
    805805    rf.SetNumTrees(fNumTrees);
    806806    rf.SetNdSize(fNdSize);
  • trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc

    r8644 r8646  
    11/* ======================================================================== *\
    2 ! $Name: not supported by cvs2svn $:$Id: MRanForestCalc.cc,v 1.26 2007-07-24 13:35:39 tbretz Exp $
     2! $Name: not supported by cvs2svn $:$Id: MRanForestCalc.cc,v 1.27 2007-07-26 11:13:00 tbretz Exp $
    33! --------------------------------------------------------------------------
    44!
     
    229229        tlist.AddToList(&fillh);
    230230
    231         MEvtLoop evtloop(fName);
     231        MEvtLoop evtloop(fTitle);
    232232        evtloop.SetParList(&plist);
    233233        evtloop.SetDisplay(fDisplay);
Note: See TracChangeset for help on using the changeset viewer.