Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 8645) +++ trunk/MagicSoft/Mars/Changelog (revision 8646) @@ -18,4 +18,37 @@ -*-*- END OF LINE -*-*- + + 2007/07/26 Thomas Bretz + + * mfilter/MFMagicCuts.cc: + - implemented new disp-parametrization + - implemented new ghostbusting + - implemented possibility to use an external ghostbuster + - implemented possibility to use an external disp calculator + + * mimage/MImgCleanStd.[h,cc]: + - implemented the cleaning in a recursive way. It should be + a little bit faster + - implemented post-cleaning using timing information + - implemented new options how to treat the removed single + core pixels + + * mjobs/MJCut.cc: + - write an external ghostbuster to the output file if available + + * mjtrain/MJTrainEnergy.cc, mjtrain/MJTrainSeparation.cc: + - set display name of MRanForestCalc as title instead of name + + * mranforest/MRanForestCalc.cc: + - set fTitle as eventloop name instead of fName + + + + 2007/07/25 Thomas Bretz + + * scripts/merppupdate: + - small fixes + + 2007/07/24 Thomas Bretz Index: trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc (revision 8645) +++ trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc (revision 8646) @@ -51,14 +51,14 @@ // // c[0], c[5], c[6], c[7], c[8], c[9]: -// xi = c[0]+c[6]*pow(leakage1, c[7]); +// xi = c[0] + c[8]*slope + c[9]*leak + +// (lgsize>c[10])*c[11]*(lgsize-c[10])^2; // p = xi*(1-width/length); // sign1 = m3long-c[5] -// sign2 = (dist-c[9])*c[8]-slope +// sign2 = (dist-c[7])*c[6]-slope // disp = sign1<0 ||sign2<0 ? -disp : disp -// antisign1 = antim3long-c[5] -// antisign2 = (antidist-c[9])*c[8]-antislope -// antidisp = antisign1<0 || antisign2<0 ? -antidisp : antidisp // thetasq = disp^2 + dist^2 - 2*disp*dist*alpha -// antithetasq = antidisp^2 + antidist^2 - 2*antidisp*antidist*antialpha +// +// And the values with respect to the antisource position respectively. +// // // c[1]: @@ -71,7 +71,7 @@ // // HadronnessCut: -// c[10], c[11]: -// had <= c[10] -// 10^lgsize >= c[11] +// c[13], c[14]: +// had <= c[13] +// 10^lgsize >= c[14] // // @@ -81,13 +81,21 @@ // HadronnessCut: // - none/nocut: No area cut -// - area: Area cut -// - all: same as area +// - area: Area cut +// - hadronness: Cut in size and hadronness (c[10], c[11]) +// - all: area + hadronness // // ThetaCut: -// - none/nocut: no theta cut +// - none/nocut: no theta cut // - on: cut in theta only // - off: anti-cut in anti-theta only // - wobble: same as on|off (both) // +// CalcDisp: +// - yes: Disp is calculated as defined above +// - no: abs(Disp.fVal) from the parameter list is used instead +// +// CalcGhostbuster: +// - yes: The ghostbuster is calculated as defined above +// - no: Ghostbuster.fValFindCreateObj("MParameterD", "Disp"); + + if (!fCalcDisp) + fDisp = (MParameterD*)pList->FindObject("Disp", "MParameterD"); + else + fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp"); if (!fDisp) - return kFALSE; + { + *fLog << err << "Disp [MParameterD] not found... aborting." << endl; + return kFALSE; + } + + if (!fCalcGhostbuster) + fGhostbuster = (MParameterD*)pList->FindObject("Ghostbuster", "MParameterD"); + else + fGhostbuster = (MParameterD*)pList->FindCreateObj("MParameterD", "Ghostbuster"); + if (!fGhostbuster) + { + *fLog << err << "Ghostbuster [MParameterD] not found... aborting." << endl; + return kFALSE; + } // propagate Theta cut to the parameter list @@ -297,7 +317,13 @@ fMap[kEM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg"); - fMap[kESrcSign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)"); + fMap[kESign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)"); fMap[kESlope] = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg"); + + if (!fCalcDisp) + fMap[kEDisp] = fMatrix->AddColumn("abs(Disp.fVal)"); + + if (!fCalcGhostbuster) + fMap[kEGhostbuster] = fMatrix->AddColumn("Ghostbuster.fVal"); if (fThetaCut&kOff) @@ -321,5 +347,5 @@ // If par!=NULL p is stored in this MParameterD // -Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t wl, Double_t d, Double_t a) const +Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t d, Double_t a) const { // wl = width/l [1] @@ -338,43 +364,77 @@ } -// --------------------------------------------------------------------------- -// -// Evaluate dynamical magic-cuts -// -Int_t MFMagicCuts::Process() -{ +// -------------------------------------------------------------------------- +// +// Return abs(Disp.fVal) if disp calculation is switched off. +// Otherwise return (c0+c6*leak^c7)*(1-width/length) +// +Double_t MFMagicCuts::GetDisp(Double_t slope, Double_t lgsize) const +{ + if (!fCalcDisp) + return fMatrix ? GetVal(kEDisp) : TMath::Abs(fDisp->GetVal()); + // Get some variables - const Double_t wdivl = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength(); - const Double_t lgsize = fMatrix ? GetVal(kESize) : TMath::Log10(fHil->GetSize()); - const Double_t leak = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1); + const Double_t wdivl = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength(); + const Double_t leak = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1); // For simplicity const Double_t *c = fVariables.GetArray(); + // As rule for root or MDataPhrase: + // ((M[3]>[3])*[4]*(M[3]-[3])^2 + [2]*M[2] + [1]*M[1] + [0])*M[0] + // + Double_t xi = c[0] + c[8]*slope + c[9]*leak; + if (lgsize>c[10]) + xi += (lgsize-c[10])*(lgsize-c[10])*c[11]; + + const Double_t disp = xi*(1-wdivl); + + return disp; +} + +Bool_t MFMagicCuts::IsGhost(Double_t m3long, Double_t slope, Double_t dist) const +{ + // For simplicity + const Double_t *c = fVariables.GetArray(); + + if (!fCalcGhostbuster) + return (fMatrix ? GetVal(kEGhostbuster) : fGhostbuster->GetVal()) (dist-c[7])*c[6]; + + return sign1 || sign2; +} + +// --------------------------------------------------------------------------- +// +// Evaluate dynamical magic-cuts +// +Int_t MFMagicCuts::Process() +{ + // For simplicity + const Double_t *c = fVariables.GetArray(); + + // Default if we return before the end + fResult = kFALSE; + // Value for Theta cut (Disp parametrization) const Double_t cut = GetThetaSqCut(); - const Double_t xi = c[0]+c[6]*pow(leak, c[7]); - const Double_t disp = xi*(1-wdivl); - - // Default if we return before the end - fResult = kFALSE; // ------------------- Most efficient cut ----------------------- // ---------------------- Theta cut ON -------------------------- - const Double_t dist = fMatrix ? GetVal(kEDist) : fHilSrc->GetDist()*fMm2Deg; - const Double_t alpha = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha(); - const Double_t sign = fMatrix ? GetVal(kESrcSign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha()); - const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha()); - const Double_t slope = fMatrix ? GetVal(kESlope) : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha()); - - // Do Ghostbusting with slope and m3l - const Double_t sign1 = (dist-c[9])*c[8]-slope; - const Double_t sign2 = m3long-c[5]; + const Double_t dist = fMatrix ? GetVal(kEDist) : fHilSrc->GetDist()*fMm2Deg; + const Double_t alpha = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha(); + const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha()); + const Double_t slope = fMatrix ? GetVal(kESlope) : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha()); + const Double_t sign = fMatrix ? GetVal(kESign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha()); + const Double_t lgsize = fMatrix ? GetVal(kESize) : TMath::Log10(fHil->GetSize()); // Calculate disp including sign - const Double_t p = sign1<0 || sign2<0 ? -disp : disp; - - // Align disp along source direction depending on third moment - //const Double_t p = TMath::Sign(disp, m3long-c[5]); + const Double_t disp = GetDisp(slope, lgsize); + const Double_t ghost = IsGhost(m3long, slope, dist); + + const Double_t p = ghost ? -disp : disp; // Align sign of disp along shower axis like m3long @@ -382,5 +442,5 @@ // Calculate corresponding Theta^2 - const Double_t thetasq = GetThetaSq(p, wdivl, dist, alpha); + const Double_t thetasq = GetThetaSq(p, dist, alpha); fThetaSq->SetVal(thetasq); @@ -405,8 +465,8 @@ { const Double_t had = fMatrix ? GetVal(kEHadronness) : fHadronness->GetVal(); - if (had>c[10]) + if (had>c[13]) return kTRUE; - if (TMath::Power(10, lgsize)GetSlopeLong()/TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha()); - // Do Ghostbusting with slope and m3l - const Double_t sign3 = (distanti-c[9])*c[8]-slopeanti; - const Double_t sign4 = m3lanti-c[5]; - const Double_t panti = sign3<0 || sign4<0 ? -disp : disp; + const Double_t dispanti = GetDisp(slopeanti, lgsize); + const Double_t ghostanti = IsGhost(m3lanti, slopeanti, lgsize); + + const Double_t panti = ghostanti ? -dispanti : dispanti; // Align disp along source direction depending on third moment - //const Double_t panti = TMath::Sign(disp, m3lanti-c[5]); - const Double_t thetasqanti = GetThetaSq(panti, wdivl, distanti, alphaanti); + const Double_t thetasqanti = GetThetaSq(panti, distanti, alphaanti); if (thetasqanti ! -! Copyright: MAGIC Software Development, 2000-2006 +! Copyright: MAGIC Software Development, 2000-2007 ! ! @@ -279,4 +279,12 @@ // // +// Class Version 4: +// ---------------- +// + Float_t fTimeLvl2; +// - Bool_t fKeepSinglePixels; +// + Bool_t fKeepIsolatedPixels; +// + Int_t fRecoverIsolatedPixels; +// +// // Input Containers: // MGeomCam @@ -302,4 +310,5 @@ #include "MLogManip.h" +#include "MArrayI.h" #include "MParList.h" #include "MCameraData.h" @@ -339,11 +348,30 @@ MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2, const char *name, const char *title) - : fCleaningMethod(kStandard), fCleanLvl1(lvl1), - fCleanLvl2(lvl2), fCleanRings(1), fKeepSinglePixels(kFALSE), - fNamePedPhotCam(gsNamePedPhotCam), fNameGeomCam(gsNameGeomCam), - fNameSignalCam(gsNameSignalCam) + : fCleaningMethod(kStandard), fCleanLvl0(lvl1), fCleanLvl1(lvl1), + fCleanLvl2(lvl2), fTimeLvl1(1.5), fTimeLvl2(1.5), fCleanRings(1), + fKeepIsolatedPixels(kFALSE), fRecoverIsolatedPixels(0), + fPostCleanType(0), fNamePedPhotCam(gsNamePedPhotCam), + fNameGeomCam(gsNameGeomCam), fNameSignalCam(gsNameSignalCam) { fName = name ? name : gsDefName.Data(); fTitle = title ? title : gsDefTitle.Data(); +} + +void MImgCleanStd::ResetCleaning() const +{ + // + // check the number of all pixels against the noise level and + // set them to 'unused' state if necessary + // + const UInt_t npixevt = fEvt->GetNumPixels(); + for (UInt_t idx=0; idxGetData(); +#else + const Double_t *data = fData->GetData().GetArray(); +#endif + + //loop on the neighbors to check if they are used + const Int_t n = gpix.GetNumNeighbors(); + for (Int_t i=0; ifCleanRings) + return; + +#ifdef DEBUG + const TArrayD &data = fData->GetData(); +#else + const Double_t *data = fData->GetData().GetArray(); +#endif + + // At least one neighbor has been identified as core, + // that means we will keep the pixel + const Int_t n = gpix.GetNumNeighbors(); + for (Int_t i=0; i0) + // and unmapped status is exclusive + if (pix.IsPixelUnmapped()) + continue; + + // Set or reset the ring number + pix.SetRing(r); + + // Step forward to the surrounding pixels + SetUsedNeighbors((*fCam)[idx], r+1); + } +} + +// -------------------------------------------------------------------------- +// +// Here we do the cleaning. We search for all the possible core candidates +// and from them on we recursively search for used pixels with +// SetUsedNeighbors. To check the validity of a core pixel +// either fTimeLvl2 and/or HasCoreNeighbors is used. +// The size of all removed +Int_t MImgCleanStd::DoCleaning(Float_t &size) const +{ + Int_t n = 0; + size = 0; + +#ifdef DEBUG + const TArrayD &data = fData->GetData(); +#else + const Double_t *data = fData->GetData().GetArray(); +#endif // @@ -378,183 +517,213 @@ for (UInt_t idx=0; idxfCleanLvl1) - continue; - - // Setting a pixel to unused if it is unmapped would overwrite - // the unmapped-status. Therefor this pixels are excluded. MSignalPix &pix = (*fEvt)[idx]; - if (!pix.IsPixelUnmapped()) - pix.SetPixelUnused(); - } -} - -// -------------------------------------------------------------------------- -// -// Check if the survived pixel have a neighbor, that also -// survived. Set all single pixels Unused if !fKeepSinglePixels. Now we -// declare all pixels that survived previous CleanSteps as CorePixels. -// Return number of single pixels, and there cumulative size in size. -// -Short_t MImgCleanStd::CleanStep2(Float_t &size) -{ - Short_t n=0; - size = 0; - + + // Check if this pixel is a possible candidate for a core pixel + if (data[idx] <= fCleanLvl1) + continue; + + // Ignore unmapped pixels + if (pix.IsPixelUnmapped()) + continue; + + const MGeomPix &gpix = (*fCam)[idx]; + + // Check if the pixel is an isolated core pixel + if (!HasCoreNeighbors(gpix)) + { + // Count size and number of isolated core pixels + size += pix.GetNumPhotons(); + n++; + + // If isolated pixels should not be kept or the pixel + // is lower than the cleaning level for isolated core + // pixels. It is not treated as core pixel. + if (!fKeepIsolatedPixels || data[idx]<=fCleanLvl0) + continue; + } + + // Mark pixel as used and core + pix.SetPixelUsed(); + pix.SetPixelCore(); + + // Check if neighbor pixels should be marked as used + // This is done recursively depening on fCleanRings + SetUsedNeighbors(gpix); + } + + return n; +} + +/* +Float_t MImgCleanStd::GetArrivalTimeNeighbor(const MGeomPix &gpix) const +{ + Float_t min = FLT_MAX; + + const Int_t n = gpix.GetNumNeighbors(); + for (int i=0; iGetData(); + + // + // check the number of all pixels against the noise level and + // set them to 'unused' state if necessary + // const UInt_t npixevt = fEvt->GetNumPixels(); for (UInt_t idx=0; idxGetData(); +#else + const Double_t *data = fData->GetData().GetArray(); +#endif + + Int_t n = 0; + + // + // check the number of all pixels against the noise level and + // set them to 'unused' state if necessary + // + const UInt_t npixevt = fEvt->GetNumPixels(); + for (UInt_t idx=0; idxGetNumPixels(); + for (UInt_t idx=0; idxlvl) + continue; + + // Now count the pixel. If we did not found n pixels yet + // which fullfill the condition, go on searching + if (++cnt>=n) break; - } } - // If the pixel has at least one core-neighbor - // go on with the next pixel - if (hasNeighbor) - continue; - - // If the pixel has no neighbors and the single pixels - // should not be kept turn the used- into an unused-status - if (!fKeepSinglePixels) - pix.SetPixelUnused(); - - // count size and number of single core-pixels - size += pix.GetNumPhotons(); - n++; - } - - // Now turn the used-status into the core-status - // (FIXME: A more intelligent handling of used/core in clean step1/2 - // would make this loop obsolete!) - for (UInt_t idx=0; idx1) times for the boundary pixels around the used pixels. -// If a pixel has more than 2.5 (clean level 2.5) sigma, -// it is declared as used. -// -// If a value<2 for fCleanRings is used, no CleanStep4 is done. -// -void MImgCleanStd::CleanStep4(UShort_t r, Int_t idx) -{ - MSignalPix &pix = (*fEvt)[idx]; - - // - // Skip events that have already a defined status; - // - if (pix.GetRing() != 0) + // If we found at least n neighbors which are + // with a time difference of lvl keep the pixel + if (cnt>=n) + continue; + + indices.Set(indices.GetSize()+1); + indices[indices.GetSize()-1] = idx; + } + + // Now remove the pixels which didn't fullfill the requirement + for (UInt_t i=0; ir-1 ) - continue; - - pix.SetRing(r); - break; - } -} - -// -------------------------------------------------------------------------- -// -// Look for the boundary pixels around the core pixels -// if a pixel has more than 2.5 (clean level 2.5) sigma, and -// a core neighbor, it is declared as used. -// -void MImgCleanStd::CleanStep3() -{ - const TArrayD &data = fData->GetData(); - - for (UShort_t r=1; rGetNumPixels(); - for (UInt_t idx=0; idxFindObject(AddSerialNumber(fNameSignalCam), "MSignalCam"); if (!fEvt) @@ -595,4 +763,33 @@ return kFALSE; + if (fCleanLvl2>fCleanLvl1) + { + *fLog << warn << "WARNING - fCleanLvl2 (" << fCleanLvl2 << ") > fCleanLvl1 (" << fCleanLvl1 << ")... resetting fCleanLvl2." << endl; + fCleanLvl2 = fCleanLvl1; + } + + if (fCleanLvl2==fCleanLvl1 && fCleanRings>0) + { + *fLog << warn << "WARNING - fCleanLvl2 == fCleanLvl1 (" << fCleanLvl1 << ") but fCleanRings>0... resetting fCleanRings to 0." << endl; + fCleanRings=0; + } + + if (fKeepIsolatedPixels && fTimeLvl2fCleanLvl1) + { + *fLog << warn << "WARNING - fKeepIsolatedPixels not set but CleanLvl0 (" << fTimeLvl2 << ") > fCleanLvl1 (" << fCleanLvl1 << ")... setting fKeepIsolatedCorePixels." << endl; + fKeepIsolatedPixels=kTRUE; + } + + if (fKeepIsolatedPixels && fTimeLvl2==fCleanLvl1 && fRecoverIsolatedPixels!=0) + { + *fLog << warn << "WARNING - fTimeLvl2 == fCleanLvl1 (" << fTimeLvl2 << ") and fKeepSinglePixels and fRecoverIsolatedPixels!=0... setting fRecoverIsolatedPixels=0." << endl; + fRecoverIsolatedPixels = 0; + } + Print(); @@ -631,28 +828,36 @@ #ifdef DEBUG - *fLog << all << "CleanStep 1" << endl; + *fLog << all << "ResetCleaning" << endl; #endif - CleanStep1(); - + ResetCleaning(); #ifdef DEBUG - *fLog << all << "CleanStep 2" << endl; + *fLog << all << "DoCleaning" << endl; #endif Float_t size; - const Short_t n = CleanStep2(size); + Short_t n = DoCleaning(size); + +#ifdef DEBUG + *fLog << all << "RecoverIsolatedPixels" << endl; +#endif + for (UInt_t i=0; i<(UInt_t)fRecoverIsolatedPixels; i++) + { + const Int_t rc=RecoverIsolatedPixels(size); + if (rc==0) + break; + + n -= rc; + } + +#ifdef DEBUG + *fLog << all << "Time Cleaning" << endl; +#endif + // FIXME: Remove removed core piselx? + CleanStepTime(); + fEvt->SetSinglePixels(n, size); - // For speed reasons skip the rest of the cleaning if no - // action will be taken! - if (fCleanLvl1>fCleanLvl2) - { #ifdef DEBUG - *fLog << all << "CleanStep 3" << endl; -#endif - CleanStep3(); - } - -#ifdef DEBUG - *fLog << all << "Calc Islands" << endl; + *fLog << all << "CalcIslands" << endl; #endif // Takes roughly 10% of the time @@ -698,4 +903,26 @@ *fLog << "initialized with level " << fCleanLvl1 << " and " << fCleanLvl2; *fLog << " (CleanRings=" << fCleanRings << ")" << endl; + + if (fPostCleanType) + { + *fLog << "Time post cleaning is switched on with:" << endl; + if (fPostCleanType&2) + *fLog << " - Two pixel coincidence window: " << fTimeLvl2 << "ns" << endl; + if (fPostCleanType&1) + *fLog << " - One pixel coincidence window: " << fTimeLvl1 << "ns" << endl; + } + + if (fKeepIsolatedPixels) + *fLog << "isolated core pixels will be kept above " << fCleanLvl0 << endl; + + if (fRecoverIsolatedPixels) + { + if (fRecoverIsolatedPixels<0) + *fLog << "all "; + *fLog << "isolated core pixels with used pixels as neighbors will be recovered"; + if (fRecoverIsolatedPixels>0) + *fLog << " " << fRecoverIsolatedPixels << " times"; + *fLog << "." << endl;; + } if (fCleaningMethod!=kAbsolute && fCleaningMethod!=kTime) @@ -848,6 +1075,8 @@ if (gsNameSignalCam!=fNameSignalCam) out << " " << GetUniqueName() << ".SetNameSignalCam(\"" << fNameSignalCam << "\");" << endl; - if (fKeepSinglePixels) - out << " " << GetUniqueName() << ".SetKeepSinglePixels();" << endl; + if (fKeepIsolatedPixels) + out << " " << GetUniqueName() << ".SetKeepIsolatedPixels();" << endl; + if (fRecoverIsolatedPixels) + out << " " << GetUniqueName() << ".SetRecoverIsolatedPixels(" << fRecoverIsolatedPixels << ");" << endl; } @@ -860,5 +1089,5 @@ // MImgCleanStd.CleanMethod: Standard, Scaled, Democratic, Probability, Absolute // MImgCleanStd.CleanRings: 1 -// MImgCleanStd.KeepSinglePixels: yes, no +// MImgCleanStd.KeepIsolatedPixels: yes, no // Int_t MImgCleanStd::ReadEnv(const TEnv &env, TString prefix, Bool_t print) @@ -870,4 +1099,9 @@ SetCleanRings(GetEnvValue(env, prefix, "CleanRings", fCleanRings)); } + if (IsEnvDefined(env, prefix, "CleanLevel0", print)) + { + rc = kTRUE; + fCleanLvl0 = GetEnvValue(env, prefix, "CleanLevel0", fCleanLvl0); + } if (IsEnvDefined(env, prefix, "CleanLevel1", print)) { @@ -880,8 +1114,28 @@ fCleanLvl2 = GetEnvValue(env, prefix, "CleanLevel2", fCleanLvl2); } - if (IsEnvDefined(env, prefix, "KeepSinglePixels", print)) + if (IsEnvDefined(env, prefix, "TimeLevel1", print)) { rc = kTRUE; - fKeepSinglePixels = GetEnvValue(env, prefix, "KeepSinglePixels", fKeepSinglePixels); + fTimeLvl1 = GetEnvValue(env, prefix, "TimeLevel1", fTimeLvl1); + } + if (IsEnvDefined(env, prefix, "TimeLevel2", print)) + { + rc = kTRUE; + fTimeLvl2 = GetEnvValue(env, prefix, "TimeLevel2", fTimeLvl2); + } + if (IsEnvDefined(env, prefix, "KeepIsolatedPixels", print)) + { + rc = kTRUE; + fKeepIsolatedPixels = GetEnvValue(env, prefix, "KeepIsolatedPixels", fKeepIsolatedPixels); + } + if (IsEnvDefined(env, prefix, "RecoverIsolatedPixels", print)) + { + rc = kTRUE; + fRecoverIsolatedPixels = GetEnvValue(env, prefix, "RecoverIsolatedPixels", fKeepIsolatedPixels); + } + if (IsEnvDefined(env, prefix, "PostCleanType", print)) + { + rc = kTRUE; + fPostCleanType = GetEnvValue(env, prefix, "PostCleanType", fPostCleanType); } Index: trunk/MagicSoft/Mars/mimage/MImgCleanStd.h =================================================================== --- trunk/MagicSoft/Mars/mimage/MImgCleanStd.h (revision 8645) +++ trunk/MagicSoft/Mars/mimage/MImgCleanStd.h (revision 8646) @@ -7,9 +7,10 @@ class MGeomCam; -class MSigmabar; +class MGeomPix; class MSignalCam; class MPedPhotCam; class MArrivalTime; class MCameraData; +//class MRawRunHeader; class MGGroupFrame; @@ -32,16 +33,23 @@ static const TString gsNameSignalCam; // default name of the 'MSignalCam' container - const MGeomCam *fCam; //! - MSignalCam *fEvt; //! - MPedPhotCam *fPed; //! - MCameraData *fData; //! + const MGeomCam *fCam; //! + MSignalCam *fEvt; //! + MPedPhotCam *fPed; //! + MCameraData *fData; //! +// MRawRunHeader *fHeader; //! CleaningMethod_t fCleaningMethod; + Float_t fCleanLvl0; Float_t fCleanLvl1; Float_t fCleanLvl2; + Float_t fTimeLvl1; + Float_t fTimeLvl2; + UShort_t fCleanRings; - Bool_t fKeepSinglePixels; + Bool_t fKeepIsolatedPixels; + Int_t fRecoverIsolatedPixels; + Int_t fPostCleanType; TString fNamePedPhotCam; // name of the 'MPedPhotCam' container @@ -50,9 +58,13 @@ // MImgCleanStd - void CleanStep1(); - Short_t CleanStep2(Float_t &size); - void CleanStep3(); - void CleanStep3b(Int_t idx); - void CleanStep4(UShort_t r, Int_t idx); + Bool_t HasCoreNeighbors(const MGeomPix &gpix) const; + Bool_t HasUsedNeighbors(const MGeomPix &gpix) const; + void SetUsedNeighbors(const MGeomPix &gpix, Int_t r=1) const; + Int_t DoCleaning(Float_t &size) const; + void ResetCleaning() const; + Int_t RecoverIsolatedPixels(Float_t &size) const; + void CleanTime(Int_t n, Double_t lvl) const; + + void CleanStepTime() const; // MGTask, MTask, MParContainer @@ -70,11 +82,24 @@ void Print(Option_t *o="") const; - Float_t GetCleanLvl1() const { return fCleanLvl1; } - Float_t GetCleanLvl2() const { return fCleanLvl2; } + Float_t GetCleanLvl0() const { return fCleanLvl0; } + Float_t GetCleanLvl1() const { return fCleanLvl1; } + Float_t GetCleanLvl2() const { return fCleanLvl2; } + + Float_t GetTimeLvl1() const { return fTimeLvl1; } + Float_t GetTimeLvl2() const { return fTimeLvl2; } + + void SetCleanLvl0(Float_t lvl) { fCleanLvl0=lvl; } + void SetCleanLvl1(Float_t lvl) { fCleanLvl1=lvl; } + void SetCleanLvl2(Float_t lvl) { fCleanLvl2=lvl; } + + void SetTimeLvl1(Float_t lvl) { fTimeLvl1=lvl; } + void SetTimeLvl2(Float_t lvl) { fTimeLvl2=lvl; } + + void SetCleanRings(UShort_t r) { fCleanRings=r; } UShort_t GetCleanRings() const { return fCleanRings;} - void SetCleanRings(UShort_t r) { if(r==0) r=1; fCleanRings=r; } void SetMethod(CleaningMethod_t m) { fCleaningMethod = m; } - void SetKeepSinglePixels(Bool_t b=kTRUE) { fKeepSinglePixels=b; } + void SetKeepIsolatedPixels(Bool_t b=kTRUE) { fKeepIsolatedPixels=b; } + void SetRecoverIsolatedPixels(Int_t n=-1) { fRecoverIsolatedPixels=n; } Bool_t ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2); @@ -84,5 +109,5 @@ void SetNameGeomCam(const char *name) { fNameGeomCam = name; } - ClassDef(MImgCleanStd, 3) // task doing the image cleaning + ClassDef(MImgCleanStd, 4) // task doing the image cleaning }; Index: trunk/MagicSoft/Mars/mjobs/MJCut.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 8645) +++ trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 8646) @@ -388,4 +388,5 @@ write->AddContainer("OpticalAxis", "Events", kFALSE); write->AddContainer("Disp", "Events", kFALSE); + write->AddContainer("Ghostbuster", "Events", kFALSE); write->AddContainer("MEnergyEst", "Events", kFALSE); write->AddContainer("MTime", "Events", kFALSE); Index: trunk/MagicSoft/Mars/mjtrain/MJTrainEnergy.cc =================================================================== --- trunk/MagicSoft/Mars/mjtrain/MJTrainEnergy.cc (revision 8645) +++ trunk/MagicSoft/Mars/mjtrain/MJTrainEnergy.cc (revision 8646) @@ -145,5 +145,5 @@ // ------------------------ Train RF -------------------------- - MRanForestCalc rf(fTitle); + MRanForestCalc rf("TrainEnergy", fTitle); rf.SetNumTrees(fNumTrees); rf.SetNdSize(fNdSize); Index: trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc =================================================================== --- trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc (revision 8645) +++ trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc (revision 8646) @@ -802,5 +802,5 @@ // ------------------------ Train RF -------------------------- - MRanForestCalc rf; + MRanForestCalc rf("TrainSeparation", fTitle); rf.SetNumTrees(fNumTrees); rf.SetNdSize(fNdSize); Index: trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc (revision 8645) +++ trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc (revision 8646) @@ -1,4 +1,4 @@ /* ======================================================================== *\ -! $Name: not supported by cvs2svn $:$Id: MRanForestCalc.cc,v 1.26 2007-07-24 13:35:39 tbretz Exp $ +! $Name: not supported by cvs2svn $:$Id: MRanForestCalc.cc,v 1.27 2007-07-26 11:13:00 tbretz Exp $ ! -------------------------------------------------------------------------- ! @@ -229,5 +229,5 @@ tlist.AddToList(&fillh); - MEvtLoop evtloop(fName); + MEvtLoop evtloop(fTitle); evtloop.SetParList(&plist); evtloop.SetDisplay(fDisplay);