Index: trunk/Mars/mimage/ImageLinkDef.h =================================================================== --- trunk/Mars/mimage/ImageLinkDef.h (revision 19968) +++ trunk/Mars/mimage/ImageLinkDef.h (revision 19969) @@ -8,4 +8,5 @@ #pragma link C++ class MImgCleanTime+; #pragma link C++ class MImgCleanSimple+; +#pragma link C++ class MImgCleanFAMOUS+; #pragma link C++ class MHillas+; Index: trunk/Mars/mimage/MImgCleanFAMOUS.cc =================================================================== --- trunk/Mars/mimage/MImgCleanFAMOUS.cc (revision 19969) +++ trunk/Mars/mimage/MImgCleanFAMOUS.cc (revision 19969) @@ -0,0 +1,283 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 06/2020 +! +! Copyright: MAGIC Software Development, 2000-2020 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MImgCleanFAMOUS +// +// This image cleaning algorithm was introduced by J. Audehm (Master +// thesis, RWTH Aachen, 2020) based on a cleaning by M. Schaufel +// (Master thesus, RWTH Aachen, 2017). It was slightly modified +// to fit into the framework. +// +///////////////////////////////////////////////////////////////////////////// +#include "MImgCleanFAMOUS.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MGeomPix.h" +#include "MGeomCam.h" + +#include "MSignalPix.h" +#include "MSignalCam.h" + +ClassImp(MImgCleanFAMOUS); + +using namespace std; + +static const TString gsDefName = "MImgCleanFAMOUS"; +static const TString gsDefTitle = "Task to perform image cleaning"; + +const TString MImgCleanFAMOUS::gsNameSignalCam ="MSignalCam"; // default name of the 'MSignalCam' container +const TString MImgCleanFAMOUS::gsNameGeomCam ="MGeomCam"; // default name of the 'MGeomCam' container + +// -------------------------------------------------------------------------- +// +// Default constructor. Here you can specify the cleaning method and levels. +// If you don't specify them the 'common standard' values 3.0 and 2.5 (sigma +// above mean) are used. +// Here you can also specify how many rings around the core pixels you want +// to analyze (with the fixed lvl2). The default value for "rings" is 1. +// +MImgCleanFAMOUS::MImgCleanFAMOUS(const char *name, const char *title) + : fSignalMin(40), fSignalMax(4000), fTimeMin(110), fTimeMax(150), + fSlopeMin(0), fSlopeMax(100), fTimeWindow(5), + fNameGeomCam(gsNameGeomCam), fNameSignalCam(gsNameSignalCam) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); +} + +void MImgCleanFAMOUS::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; idx=fTimeWindow) + continue; + + if (npix.GetNumPhotons()<0) + continue; + + // FIXME: Checks on minimum signal and slope parameter? + + if (npix.IsPixelUnmapped()) + continue; + + pix.SetPixelUsed(); + pix.SetIdxIsland(1); + + FindIsland(ipix); + } +} + + +void MImgCleanFAMOUS::DoCleaning() const +{ + UInt_t max_idx = -1; + + // Serach for the brightest none-saturating pixel + // that is consistent with originating form a shower + const UInt_t npixevt = fEvt->GetNumPixels(); + for (UInt_t idx=0; idxfSlopeMax) + continue; + + // Check that pixel aligns with the trigger position + if (pix.GetArrivalTime() <= fTimeMin || pix.GetArrivalTime()>fTimeMax) + continue; + + // Check if this is above the cleaning level and not saturating + if (pix.GetNumPhotons() <= fSignalMin || pix.GetNumPhotons()>fSignalMax) + continue; + + // Ignore unmapped pixels + if (pix.IsPixelUnmapped()) + continue; + + pix.SetPixelCore(); + + if (max_idx<0 || pix.GetNumPhotons()>(*fEvt)[max_idx].GetNumPhotons()) + max_idx = idx; + } + + if (max_idx<0) + return; + + (*fEvt)[max_idx].SetPixelUsed(); + FindIsland(max_idx); +} + +// -------------------------------------------------------------------------- +// +// Check if MEvtHeader exists in the Parameter list already. +// if not create one and add them to the list +// +Int_t MImgCleanFAMOUS::PreProcess(MParList *pList) +{ + fCam = (MGeomCam*)pList->FindObject(AddSerialNumber(fNameGeomCam), "MGeomCam"); + if (!fCam) + { + *fLog << err << fNameGeomCam << " [MGeomCam] not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + fEvt = (MSignalCam*)pList->FindObject(AddSerialNumber(fNameSignalCam), "MSignalCam"); + if (!fEvt) + { + *fLog << err << fNameSignalCam << " [MSignalCam] not found... aborting." << endl; + return kFALSE; + } + + Print(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Cleans the image. +// +Int_t MImgCleanFAMOUS::Process() +{ + ResetCleaning(); + DoCleaning(); + + // Takes roughly 10% of the time + fEvt->CalcIslands(*fCam); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Print descriptor and cleaning levels. +// +void MImgCleanFAMOUS::Print(Option_t *o) const +{ + *fLog << all << GetDescriptor() << " using" << endl; + *fLog << " * " << fSignalMin << " < Signal <= " << fSignalMax << endl; + *fLog << " * " << fTimeMin << " < Time <= " << fTimeMax << endl; + *fLog << " * " << fSlopeMin << " < Slope <= " << fSlopeMax << endl; + *fLog << " * " << "Delta T < " << fTimeWindow << endl; +} + +// -------------------------------------------------------------------------- +// +// Read the setup from a TEnv, eg: +// MImgCleanFAMOUS.SignalMin: 40 +// MImgCleanFAMOUS.SignalMax: 4000 +// MImgCleanFAMOUS.TimeMin: 110 +// MImgCleanFAMOUS.TimeMax: 150 +// MImgCleanFAMOUS.SlopeMin: 0 +// MImgCleanFAMOUS.SlopeMax: 100 +// MImgCleanFAMOUS.TimeWindow: 5 +// +Int_t MImgCleanFAMOUS::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = kFALSE; + if (IsEnvDefined(env, prefix, "SignalMin", print)) + { + rc = kTRUE; + fSignalMin = GetEnvValue(env, prefix, "SignalMin", fSignalMin); + } + if (IsEnvDefined(env, prefix, "SignalMax", print)) + { + rc = kTRUE; + fSignalMax = GetEnvValue(env, prefix, "SignalMax", fSignalMax); + } + if (IsEnvDefined(env, prefix, "TimeMin", print)) + { + rc = kTRUE; + fTimeMin = GetEnvValue(env, prefix, "TimeMin", fTimeMin); + } + if (IsEnvDefined(env, prefix, "TimeMax", print)) + { + rc = kTRUE; + fTimeMax = GetEnvValue(env, prefix, "TimeMax", fTimeMax); + } + if (IsEnvDefined(env, prefix, "SlopeMin", print)) + { + rc = kTRUE; + fSlopeMin = GetEnvValue(env, prefix, "SlopeMin", fSlopeMin); + } + if (IsEnvDefined(env, prefix, "SlopeMax", print)) + { + rc = kTRUE; + fSlopeMax = GetEnvValue(env, prefix, "SlopeMax", fSlopeMax); + } + if (IsEnvDefined(env, prefix, "TimeWindow", print)) + { + rc = kTRUE; + fTimeWindow = GetEnvValue(env, prefix, "TimeWindow", fTimeWindow); + } + + return rc; +} Index: trunk/Mars/mimage/MImgCleanFAMOUS.h =================================================================== --- trunk/Mars/mimage/MImgCleanFAMOUS.h (revision 19969) +++ trunk/Mars/mimage/MImgCleanFAMOUS.h (revision 19969) @@ -0,0 +1,61 @@ +#ifndef MARS_MImgCleanFAMOUS +#define MARS_MImgCleanFAMOUS + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MGeomCam; +class MGeom; +class MSignalCam; +class MArrivalTime; + +class MImgCleanFAMOUS : public MTask +{ +private: + static const TString gsNameGeomCam; // default name of the 'MGeomCam' container + static const TString gsNameSignalCam; // default name of the 'MSignalCam' container + + const MGeomCam *fCam; //! + MSignalCam *fEvt; //! + + Float_t fSignalMin; + Float_t fSignalMax; + Float_t fTimeMin; + Float_t fTimeMax; + Float_t fSlopeMin; + Float_t fSlopeMax; + Float_t fTimeWindow; + + TString fNameGeomCam; // name of the 'MGeomCam' container + TString fNameSignalCam; // name of the 'MSignalCam' container + + // MImgCleanFAMOUS + void ResetCleaning() const; + void FindIsland(Int_t) const; + void DoCleaning() const; + + // MGTask, MTask, MParContainer + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); + //void StreamPrimitive(ostream &out) const; + + Int_t PreProcess(MParList *pList); + Int_t Process(); + +public: + MImgCleanFAMOUS(const char *name=NULL, const char *title=NULL); + + void Print(Option_t *o="") const; + + void SetSignalMinMax(const Float_t &min, const Float_t &max) { fSignalMin=min; fSignalMax=max; } + void SetTimeMinMax(const Float_t &min, const Float_t &max) { fTimeMin=min; fTimeMax=max; } + void SetSlopeMinMax(const Float_t &min, const Float_t &max) { fSlopeMin=min; fSlopeMax=max; } + void SetTimeWindow(const Float_t &deltat) { fTimeWindow=deltat; } + + void SetNameSignalCam(const char *name) { fNameSignalCam = name; } + void SetNameGeomCam(const char *name) { fNameGeomCam = name; } + + ClassDef(MImgCleanFAMOUS, 1) // task doing the image cleaning +}; + +#endif Index: trunk/Mars/mimage/Makefile =================================================================== --- trunk/Mars/mimage/Makefile (revision 19968) +++ trunk/Mars/mimage/Makefile (revision 19969) @@ -29,4 +29,5 @@ MImgCleanTime.cc \ MImgCleanSimple.cc \ + MImgCleanFAMOUS.cc \ MHillas.cc \ MHillasSrc.cc \