Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.cc (revision 5169) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.cc (revision 5170) @@ -36,4 +36,5 @@ #include "MControlPlots.h" #include "MIslands.h" +#include "MImgIsland.h" #include "MHCamera.h" #include "MGeomCamMagic.h" @@ -93,5 +94,4 @@ else fCameraHisto[fMode] = new MHCamera(*fGeomCam,"Survive","Pixels surviving Image Cleaning"); - return kTRUE; } @@ -101,8 +101,26 @@ // Int_t MControlPlots::Process() -{ +{ if(!fIslands) return kTRUE; - for (UInt_t i=0;iGetNumPixels();i++) + MImgIsland *imgIsl = new MImgIsland; + TIter Next(fIslands->GetList()); + + Int_t pixNum = 0; + Int_t idPix = -1; + + while ((imgIsl=(MImgIsland*)Next())) { + + pixNum = imgIsl->GetPixNum(); + + for(Int_t k = 0; kGetPixList(k); + fCameraHisto[fMode]->Fill(idPix,1); + fCameraHisto[fMode]->SetUsed(idPix); + } + } + + /* for (UInt_t i=0;iGetNumPixels();i++) { // cout << fIslands->GetIslId(i) << " "; @@ -112,5 +130,6 @@ fCameraHisto[fMode]->SetUsed(i); } - } + }*/ + // cout << endl; return kTRUE; Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.h (revision 5169) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.h (revision 5170) @@ -8,4 +8,5 @@ class TString; class MIslands; +class MImgIsland; class MGeomCam; class MHCamera; Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc (revision 5170) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc (revision 5170) @@ -0,0 +1,114 @@ +/* ======================================================================== *\ +! +! * +! * 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 expressed +! * or implied warranty. +! * +! +! +! Author(s): Ester Aliu , 9/2004 +! +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ +#include "MImgIsland.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MImgIsland); + +using namespace std; + +MImgIsland::MImgIsland(const char *name, const char *title) +{ + + fName = name ? name : "MImgIsland"; + fTitle = title ? title : "Storage container for one island information"; + + Reset(); +} + +void MImgIsland::Reset() +{ + fPixNum = 0; + fSigToNoise = 0; + fTimeSpread = 0; + fMeanX = 0; + fMeanY = 0; + fDist = 0; + fDistW = 0; + fDistL = 0; + + fPixList.Reset(-1); + fPeakPulse.Reset(-1); + +} + +// ------------------------------------------------------------------------- +// +// Initialize +// +void MImgIsland::InitSize(Int_t i) +{ + fPixList.Set(i); + fPeakPulse.Set(i); +} + + +// ------------------------------------------------------------------------- +// +// Set the pixels id in the island +// +void MImgIsland::SetPixList(const Int_t i, const Int_t idx) +{ + fPixList[i] = idx; +} + + + + +// ------------------------------------------------------------------------- +// +// Set the arrival time in one pixel of the island +// +void MImgIsland::SetPeakPulse(const Int_t i, const Float_t t) +{ + fPeakPulse[i] = t; +} + + +// ------------------------------------------------------------------------- +// +// Print the island information +// +void MImgIsland::Print(Option_t *opt) const +{ + + // TString o = opt; + + *fLog << inf << " Number of pixels = " << fPixNum << endl; + *fLog << inf << " Signal-To-Noise =" << fSigToNoise << endl; + *fLog << inf << " Time Spread (Core pixels) = " << fTimeSpread << endl; + *fLog << inf << " DistL = " << fDistL << endl; + *fLog << inf << " DistW = " << fDistW << endl; + *fLog << inf << " DistS = " << fDistS << endl; + +} + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.h (revision 5170) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.h (revision 5170) @@ -0,0 +1,83 @@ +#ifndef MARS_MImgIsland +#define MARS_MImgIsland + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayI +#include +#endif + +#ifndef ROOT_TArrayF +#include +#endif + + +class MImgIsland : public MParContainer +{ + private: + + Int_t fPixNum; + Float_t fSigToNoise; + Float_t fTimeSpread; + Float_t fMeanX; + Float_t fMeanY; + Float_t fWidth; + Float_t fLength; + Float_t fDist; + Float_t fDistL; + Float_t fDistW; + Float_t fDistS; + + TArrayI fPixList; + TArrayF fPeakPulse; + +public: + + MImgIsland(const char *name=NULL, const char *title=NULL); + //~MImgIsland(); + + Int_t GetPixNum() { return fPixNum; } + Float_t GetSigToNoise() { return fSigToNoise; } + Float_t GetTimeSpread() { return fTimeSpread; } + Float_t GetMeanX() { return fMeanX; } + Float_t GetMeanY() { return fMeanY; } + Float_t GetWidth() { return fWidth; } + Float_t GetLength() { return fLength; } + Float_t GetDist() { return fDist; } + Float_t GetDistL() { return fDistL; } + Float_t GetDistW() { return fDistW; } + Float_t GetDistS() { return fDistS; } + + void InitSize(Int_t i); + UInt_t GetSize() const { return fPixList.GetSize(); } + + Int_t GetPixList(const Int_t i = 0) const { return fPixList.At(i); }; + Float_t GetPeakPulse(const Int_t i = 0) const { return fPeakPulse.At(i); }; + + void Reset(); + + void SetPixNum (Int_t i) { fPixNum = i;} + void SetSigToNoise(Float_t val) { fSigToNoise = val;} + void SetTimeSpread(Float_t val) { fTimeSpread = val;} + void SetMeanX (Float_t val) { fMeanX = val;} + void SetMeanY (Float_t val) { fMeanY = val;} + void SetDist (Float_t val) { fDist = val;} + void SetWidth (Float_t val) { fWidth = val;} + void SetLength (Float_t val) { fLength = val;} + void SetDistL (Float_t val) { fDistL = val;} + void SetDistW (Float_t val) { fDistW = val;} + void SetDistS (Float_t val) { fDistS = val;} + + void SetPixList( const Int_t i,const Int_t id); + void SetPeakPulse( const Int_t i,const Float_t time); + + // void Paint(Option_t *opt=NULL); + void Print(Option_t *opt=NULL) const; + + ClassDef(MImgIsland, 1) // Container that holds the island information + +}; + +#endif Index: unk/MagicSoft/Mars/mtemp/mifae/library/MIslands.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslands.cc (revision 5169) +++ (revision ) @@ -1,75 +1,0 @@ -#include - -#include "MIslands.h" - -#include "MCerPhotPix.h" -#include "MCerPhotEvt.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MIslands); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MIslands::MIslands(const char *name, const char *title) -{ - fName = name ? name : "MIslands"; - fTitle = title ? title : "Storage container for the island information of one event"; - -} - -// -------------------------------------------------------------------------- -// -// Destructor. -// -MIslands::~MIslands() -{ - -} - - -// -------------------------------------------------------------------------- -// -// Getter functions -// - - -// -------------------------------------------------------------------------- -// -// Setter functions -// - - - - -// -------------------------------------------------------------------------- -// -// Print the island parameters to *fLog -// -void MIslands::Print(Option_t *opt) const -{ - *fLog << all; - *fLog << "Island Parameters (" << GetName() << ")" << endl; - *fLog << " - Island Number = " << fIslNum << endl; - for (Int_t i = 0; i< fIslNum ; i++) - { - *fLog << " Island Id " << i << endl; - *fLog << " + Pixel Number = " << fPixNum[i] << endl; - *fLog << " + SigToNoise = " << fSigToNoise[i] << endl; - *fLog << " + TimeSpread [time slices]= " << fTimeSpread[i] << endl; - *fLog << " + MeanX [mm]= " << fMeanX[i] << endl; - *fLog << " + MeanY [mm]= " << fMeanY[i] << endl; - *fLog << " + Dist [mm]= " << fDist[i] << endl; - *fLog << " + Length of the larger island [mm] = " << fLength < -#endif - -#ifndef ROOT_TArrayI -#include -#endif - -#ifndef ROOT_TObjArray -#include -#endif - -class MIslands : public MParContainer -{ -private: - // for description see MIslands.cc - Int_t fIslNum; // number of islands - // Int_t** fIslId; //[fIslNum] island Id - // TObjArray fIslId; // island Id - Int_t fIslId[577]; //island Id - Int_t* fPixNum; //[fIslNum] number of pixels in the island - Float_t* fSigToNoise; //[fIslNum] signal to noise of the island - // Float_t** fTime; //[fIslNum] mean of the arrival time - Float_t fTime[10][577]; // mean of the arrival time - Float_t* fTimeSpread; //[fIslNum] mean arrival time spread of the core pixels of the island - Float_t* fMeanX; //[fIslNum] mean X position of the island - Float_t* fMeanY; //[fIslNum] mean Y position of the island - Float_t* fDist; //[fIslNum] dist between islands and continent(larger island) - Float_t fLength; // major axis of the larger island ellipse - Float_t fWidth; // minor axis of the larger island ellipse - Float_t* fDistL; //[fIslNum] Dist of the island divided by Length of the larger island - Float_t* fDistW; //[fIslNum] Dist of the island divided by Width of the larger island - -public: - MIslands(const char *name=NULL, const char *title=NULL); - ~MIslands(); - - // void Clear(); - void Print(Option_t *opt=NULL) const; - - //getter methods - Int_t GetIslNum() const { return fIslNum; } - Int_t GetIslId(Int_t idx) { return fIslId[idx]; } - //Int_t GetIslId(Int_t isl, Int_t idx) { return fIslId[isl][idx]; } - // TObjArray GetIslId() {return fIslId;} - Float_t GetArrivalTime(Int_t isl, Int_t idx) { return fTime[isl][idx]; } - //TObjArray GetArrivalTime() { return fTime; } - Int_t GetPixNum(Int_t isl) { return fPixNum[isl]; } - Float_t GetSigToNoise(Int_t isl) { return fSigToNoise[isl]; } - Float_t GetTimeSpread(Int_t isl) { return fTimeSpread[isl];} - Float_t GetMeanX(Int_t isl) { return fMeanX[isl];} - Float_t GetMeanY(Int_t isl) { return fMeanY[isl];} - Float_t GetDist(Int_t isl) { return fDist[isl]; } - Float_t GetDistL(Int_t isl) { return fDistL[isl]; } - Float_t GetDistW(Int_t isl) { return fDistW[isl]; } - - Float_t GetLength() const { return fLength; } - Float_t GetWidth() const { return fWidth; } - - //setter functions - void SetIslNum(Int_t nisl) { fIslNum = nisl; } - void SetIslId(Int_t idx, Int_t isl) { fIslId[idx] = isl; } - - // void SetIslId(Int_t** vect) { fIslId = vect; } - void SetArrivalTime(Int_t isl, Int_t idx, Float_t val) { fTime[isl][idx] = val;} - // void SetArrivalTime(Float_t** vect) { fTime = vect;} - void SetPixNum(Int_t* npix) { fPixNum = npix; } - void SetSigToNoise(Float_t* val) { fSigToNoise = val; } - void SetTimeSpread(Float_t* val) { fTimeSpread = val; } - void SetMeanX( Float_t* val) { fMeanX = val; } - void SetMeanY(Float_t* val) { fMeanY = val; } - void SetDist(Float_t* val) { fDist = val; } - void SetDistL(Float_t* val) { fDistL=val; } - void SetDistW(Float_t* val) { fDistW=val; } - - void SetLength(Float_t val) { fLength=val; } - void SetWidth(Float_t val) { fWidth=val; } - - - - - ClassDef(MIslands, 2) // Storage Container for Island Parameters -}; - -#endif Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc (revision 5170) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc (revision 5170) @@ -0,0 +1,760 @@ +/* ======================================================================== *\ +! +! * +! * 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): Ester Aliu, 2/2004 +| +! Last Update: 7/2004 +! +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MIslandsCalc +// +// The Island Calc task calculates some islands parameters for each +// of the events such as: +// +// - fPixNum // number of pixels in the island +// - fSigToNoise // signal to noise of the island +// - fTimeSpread // mean arrival time spread of the island +// - fMeanX // mean X position of the island +// - fMeanY // mean Y position of the island +// - fDist // dist between an island and the continent +// - fLength // major axis of the larger island ellipse +// - fWidth // minor axis of the larger island ellipse +// - fDistL // dist divided by lenght of the larger island +// - fDistW // dist divided by width of the larger island +// - fDistS // dist divided by size of the larger island +// +// - fPixList // list of pixels in the island (TArrayI) +// - fPeakPulse // mean arrival time of the pixels in the island (TArrayF) +// +// Input Containers: +// MGeomCam +// MCerPhotEvt +// MPedestalCam +// MArrivalTimeCam +// +// Output Containers: +// MIslands +// MImgIsland +// +///////////////////////////////////////////////////////////////////////////// +#include "MIslandsCalc.h" + +#include // atof +#include // ofstream, SavePrimitive + +#include "MLog.h" +#include "MLogManip.h" + +#include "MIslands.h" +#include "MImgIsland.h" + +#include "MParList.h" + +#include "MGeomPix.h" +#include "MGeomCam.h" + +#include "MCerPhotPix.h" +#include "MCerPhotEvt.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" + +#include "MArrivalTimeCam.h" +#include "MArrivalTimePix.h" + +ClassImp(MIslandsCalc); + + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MIslandsCalc::MIslandsCalc(const char* name, const char* title) + : fIsl(NULL) +{ + fName = name ? name : "MIslandsCalc"; + fTitle = title ? title : "Calculate island parameters"; +} + + +// -------------------------------------------------------------------------- +Int_t MIslandsCalc::PreProcess (MParList *pList) +{ + fCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); + if (!fCam) + { + *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt")); + if (!fEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam")); + if (!fPed) + { + *fLog << dbginf << "MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + fTime = (MArrivalTimeCam*)pList->FindObject(AddSerialNumber("MArrivalTimeCam")); + if (!fTime) + { + *fLog << dbginf << "MArrivalTimeCam not found... aborting." << endl; + return kFALSE; + } + + if (strlen(fIslName) > 0) + { + fIsl = (MIslands*)pList->FindCreateObj("MIslands", AddSerialNumber(fIslName)); + } + else + { + fIsl = (MIslands*)pList->FindCreateObj(AddSerialNumber("MIslands")); + } + if (!fIsl) + return kFALSE; + + return kTRUE; +} + + +Int_t MIslandsCalc::Process(){ + + fIsl->GetList()->Delete(); + IslandPar(); + + return kTRUE; +} + + +Int_t MIslandsCalc::IslandPar(){ + + //calculates all the island parameters + + const Int_t nPix=fCam->GetNumPixels(); + const Int_t nVect=50; + Int_t numisl; + + Int_t** vect; + vect = new Int_t*[nVect]; + for(Int_t i=0;iSetIslNum(numisl); + + //examine each island... + + Float_t noise; + Float_t signal; + + Int_t PixelNumIsl[numisl]; + Float_t SigToNoise[numisl]; + Float_t time[nPix]; + Float_t timeVariance[numisl]; + Float_t meanX[numisl]; + Float_t meanY[numisl]; + Float_t dist[numisl]; + Float_t distL[numisl]; + Float_t distW[numisl]; + Float_t distS[numisl]; + + Float_t size[numisl], sizeLargeIsl, length, width; + + Float_t X = 0; + Float_t Y = 0; + sizeLargeIsl = 0; + + for(Int_t i = 1; i<=numisl ; i++) + { + + MImgIsland *imgIsl = new MImgIsland; + + imgIsl->InitSize(num[i]); + + Int_t n = 0; + + Float_t MIN = 10000.; + Float_t MAX = 0.; + + signal = 0; + noise = 0; + + size[i-1] = 0; + meanX[i-1] = 0; + meanY[i-1] = 0; + dist[i-1] = 0; + + PixelNumIsl[i-1] = 0; + timeVariance[i-1] = 0; + + + for(Int_t idx=0 ; idxGetPixById(idx); + if(!pix) continue; + const MGeomPix &gpix2 = (*fCam)[pix->GetPixId()]; + const MPedestalPix &ped = (*fPed)[idx]; + const MArrivalTimePix &timepix = (*fTime)[idx]; + const Float_t nphot = pix->GetNumPhotons(); + + if (vect[i][idx]==1){ + + PixelNumIsl[i-1]++; + signal += nphot * (fCam->GetPixRatio(idx)); + noise += pow(ped.GetPedestalRms(),2); + + size[i-1] += nphot; + if (i == 1) + sizeLargeIsl += nphot; + + meanX[i-1] += nphot * gpix2.GetX(); + meanY[i-1] += nphot * gpix2.GetY(); + + time[i-1] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain(); + + imgIsl->SetPixList(PixelNumIsl[i-1]-1,pix->GetPixId()); + imgIsl->SetPeakPulse(PixelNumIsl[i-1]-1,time[i-1]); + + //calculates the time spread only for core pixels + if (fEvt->IsPixelCore(idx)){ + + if (time[i-1] > MAX) + MAX = time[i-1]; + if (time[i-1] < MIN) + MIN = time[i-1]; + + } + + n++; + } + } + + meanX[i-1] /= size[i-1]; + meanY[i-1] /= size[i-1]; + + + if (i == 1){ + X = meanX[i-1]; + Y = meanY[i-1]; + } + + dist[i-1] = TMath::Power(meanX[i-1]-X,2) + TMath::Power(meanY[i-1]-Y,2); + dist[i-1] = TMath::Sqrt(dist[i-1]); + + timeVariance[i-1] = MAX-MIN; + + //noise = 0, in the case of MC w/o noise + if (noise == 0) noise = 1; + + SigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise); + + imgIsl->SetPixNum(PixelNumIsl[i-1]); + imgIsl->SetSigToNoise(SigToNoise[i-1]); + imgIsl->SetTimeSpread(timeVariance[i-1]); + imgIsl->SetMeanX(meanX[i-1]); + imgIsl->SetMeanY(meanY[i-1]); + imgIsl->SetDist(dist[i-1]); + + fIsl->GetList()->Add(imgIsl); + + } + + + //Length and Width of the larger island according the definition of the hillas parameters + + // calculate 2nd moments + // --------------------- + Double_t corrxx=0; // [m^2] + Double_t corrxy=0; // [m^2] + Double_t corryy=0; // [m^2] + + for(Int_t idx=0 ; idxGetPixById(idx); + if(!pix) continue; + const MGeomPix &gpix3 = (*fCam)[pix->GetPixId()]; + const Float_t nphot = pix->GetNumPhotons(); + + // if (pix == NULL) break; + + if (vect[1][idx]==1){ + + const Float_t dx = gpix3.GetX() - X; // [mm] + const Float_t dy = gpix3.GetY() - Y; // [mm] + + + corrxx += nphot * dx*dx; // [mm^2] + corrxy += nphot * dx*dy; // [mm^2] + corryy += nphot * dy*dy; // [mm^2] + + } + } + + + // calculate the hillas parameters Width and Length + + MImgIsland *imgIsl = new MImgIsland; + TIter Next(fIsl->GetList()); + + Int_t i = 1; + while ((imgIsl=(MImgIsland*)Next())) { + + const Double_t d0 = corryy - corrxx; + const Double_t d1 = corrxy*2; + const Double_t d2 = d0 + TMath::Sqrt(d0*d0 + d1*d1); + const Double_t tand = d2 / d1; + const Double_t tand2 = tand*tand; + + const Double_t s2 = tand2+1; + + const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/size[i-1]; + const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/size[i-1]; + + // + // fLength^2 is the second moment along the major axis of the ellipse + // fWidth^2 is the second moment along the minor axis of the ellipse + // + // From the algorithm we get: fWidth <= fLength is always true + // + // very small numbers can get negative by rounding + // + length = axis1<0 ? 0 : TMath::Sqrt(axis1); // [mm] + width = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm] + + + distL[i-1]=dist[i-1]/length; + distW[i-1]=dist[i-1]/width; + distS[i-1]= dist[i-1]/size[i-1]; + imgIsl->SetLength(length); + imgIsl->SetWidth(width); + imgIsl->SetDistL(distL[i-1]); + imgIsl->SetDistW(distW[i-1]); + imgIsl->SetDistS(distS[i-1]); + i++; + } + + //fIsl->SetReadyToSave(); + + for (Int_t i = 0; i< nVect; i++) + delete [] vect[i]; + + delete vect; + + return kTRUE; +} + +//------------------------------------------------------------------------------------------ +void MIslandsCalc::Calc1(Int_t& numisl, const Int_t nv, const Int_t npix, Int_t** vect, Int_t* num){ + + + ///////////////////////////// + // + // ALGORITHM # 1 + // counts the number of islands as you can see in + // the event display after doing the std image cleaning + // + ///////////////////////////// + + Int_t sflag; + Int_t control; + + Int_t nvect = 0; + + numisl = 0; + + Int_t zeros[nv]; + + for(Int_t m = 0; m < nv ; m++) + for(Int_t n = 0; n < npix ; n++) + vect[m][n] = 0; + + for(Int_t n = 0; n < nv ; n++) + zeros[n] = 0; + + //cout << "new event" <(Next()))) + { + const Int_t idx = pix->GetPixId(); + + const MGeomPix &gpix = (*fCam)[idx]; + const Int_t nnmax = gpix.GetNumNeighbors(); + + if( fEvt->IsPixelUsed(idx)) + { + //cout << idx <pixMAX) + { + pixMAX = numpixels; + idMAX = l; + } + l++; + } + } + + //the larger island will correspond to the 1st component of the vector + + num[nvect+1] = num[1]; + num[1] = num[idMAX]; + num[idMAX] = num[nvect+1]; + + + for(Int_t k = 0; k(Next0()))) + { + const Int_t idx = pix->GetPixId(); + + const MGeomPix &gpix = (*fCam)[idx]; + const Int_t nnmax = gpix.GetNumNeighbors(); + + if( fEvt->IsPixelUsed(idx)) + { + kk[idx] = 1 ; + for(Int_t j=0; j< nnmax; j++) + { + kk[gpix.GetNeighbor(j)] = 1; + } + } + + } + + //2nd loop over all pixels + MCerPhotEvtIter Next(fEvt, kFALSE); + + while ((pix=static_cast(Next()))) + { + const Int_t idx = pix->GetPixId(); + + const MGeomPix &gpix = (*fCam)[idx]; + const Int_t nnmax = gpix.GetNumNeighbors(); + + if ( kk[idx] > 0) + { + sflag = 0; + + for(Int_t j=0; j < nnmax ; j++) + { + const Int_t idx2 = gpix.GetNeighbor(j); + + if (idx2 < idx) + { + for(Int_t k = 1; k <= nvect; k++) + { + if (vect[k][idx2] == 1) + { + sflag = 1; + vect[k][idx] = 1; + } + } + } + } + + if (sflag == 0) + { + nvect++; + vect[nvect][idx] = 1; + } + + } + } + + numisl = nvect; + + + // Repeated Chain Corrections + + Int_t jmin = 0; + + for(Int_t i = 1; i <= nvect; i++) + { + for(Int_t j = i+1; j <= nvect; j++) + { + control = 0; + for(Int_t k = 0; k < npix; k++) + { + + if (vect[i][k] == 1 && vect[j][k] == 1) + { + control = 1; + break; + } + else if(vect[i][k] == 1 && vect[j][k] == 0){ + + for(Int_t jj=1 ;jjpixMAX) + { + pixMAX = numpixels; + idMAX = l; + } + l++; + } + } + + + //the larger island will correspond to the 1st component of the vector + + num[nvect +1] = num[1]; + num[1] = num[idMAX]; + num[idMAX]=num[nvect+1]; + + for(Int_t k = 0; k +#endif + +class MGeomCam; +class MSigmabar; +class MCerPhotPix; +class MCerPhotEvt; +class MPedestalCam; +class MArrivalTimeCam; +class MArrivalTimePix; +class MGeomCam; +class MIslands; +class MImgIsland; + +class MIslandsCalc : public MGTask +{ + private: + const MGeomCam *fCam; //! + MCerPhotEvt *fEvt; //! + MSigmabar *fSgb; //! + MPedestalCam *fPed; //! + MArrivalTimeCam *fTime; //! + MGeomCam *fGeomCam; + + MIslands *fIsl; //! output container to store result + + TString fIslName; // name of the 'MIslands' container + + Int_t fIslandAlgorithm; + + Int_t PreProcess(MParList *plist); + Int_t Process(); + Int_t IslandPar(); // + void Calc1(Int_t&,const Int_t,const Int_t,Int_t**,Int_t*); // algorithm of counting islands #1 + void Calc2(Int_t&,const Int_t,const Int_t,Int_t**,Int_t*); // algorithm of counting islands #2 + + + public: + MIslandsCalc(const char* name=NULL, const char* title=NULL); + void SetOutputName(TString outname) { fIslName = outname; } + + void SetAlgorithm(Int_t m) {fIslandAlgorithm = m;} + + ClassDef(MIslandsCalc, 0) // task doing the image cleaning +}; + +#endif Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsClean.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsClean.cc (revision 5170) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsClean.cc (revision 5170) @@ -0,0 +1,225 @@ +/* ======================================================================== *\ +! +! +! Author(s): Ester Aliu, 3/2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MIslandClean +// +// The Island Cleaning task selects the islands you use for the Hillas +// parameters calculation, after the normal image cleaning. +// +// There are two methods to make the selection: +// +// - No time method, as used is Whipple. It calculates an island parameter +// called "signal to noise" and adds a new threshold that eliminates some +// of the islands. The way the island is eliminated is seeting the pixels +// of the islands as UNUSED +// +// - Time method, taking profit of the time information in MAGIC. +// Calculates the maximum time difference (in time slices) for each island +// and corrects for the island size. With an other new threshold "noise" +// islands are supposed to be eliminated. +// +// Other cleanings that are allowed in this code are: +// +// - Resting only with the continent, i.e. the larger island +// +// Example: +// +// MIslands isl; +// isl.SetName("MIslands1"); + +// MImgCleanStd clean; +// MIslandClean islclean(0.2); +// islclean.SetInputName("MIslands1"); +// islclean.SetMethod(0); // for timing method +// +// tlist.AddToList(&clean); +// tlist.AddToList(&islclean); +// +// +// Input Containers: +// MGeomCam +// MCerPhotEvt +// MPedestalCam +// MArrivalTime +// MIslands +// +// Output Containers: +// MCerPhotEvt +// +///////////////////////////////////////////////////////////////////////////// +#include "MIslandsClean.h" + +#include // atof +#include // ofstream, SavePrimitive + +#include "MLog.h" +#include "MLogManip.h" + +#include "MIslands.h" +#include "MImgIsland.h" + +#include "MParList.h" + +#include "MGeomPix.h" +#include "MGeomCam.h" + +#include "MCerPhotPix.h" +#include "MCerPhotEvt.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" + +#include "MArrivalTimeCam.h" +#include "MArrivalTimePix.h" + +ClassImp(MIslandsClean); + + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MIslandsClean::MIslandsClean(const Float_t newThres, const char *name, const char *title) + : fIsl(NULL), fIslandCleaningMethod(kNoTiming), fIslCleanThres(newThres) +{ + fName = name ? name : "MIslandsClean"; + fTitle = title ? title : "Clean islands"; +} + + +Int_t MIslandsClean::PreProcess (MParList *pList) +{ + fCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); + if (!fCam) + { + *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt")); + if (!fEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam")); + if (!fPed) + { + *fLog << dbginf << "MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + fTime = (MArrivalTimeCam*)pList->FindObject(AddSerialNumber("MArrivalTimeCam")); + if (!fTime) + { + *fLog << dbginf << "MArrivalTimeCam not found... aborting." << endl; + return kFALSE; + } + + if (strlen(fIslName) > 0) + fIsl = (MIslands*)pList->FindObject(AddSerialNumber(fIslName)); + else + fIsl = (MIslands*)pList->FindObject(AddSerialNumber("MIslands")); + if (!fIsl) + { + *fLog << dbginf << "MIslands not found... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + + + +Int_t MIslandsClean::Process() +{ + // + //eliminates the island with a signal-to-noise + //lower than a given limit + // + //if ( fIslandCleaningMethod == kNoTiming ){ + + MImgIsland *imgIsl = new MImgIsland; + TIter Next(fIsl->GetList()); + + Int_t pixNum = 0; + Int_t idPix = -1; + + if ( fIslandCleaningMethod == 1 ) { + while ((imgIsl=(MImgIsland*)Next())) { + + if(imgIsl->GetSigToNoise() < fIslCleanThres) + { + pixNum = imgIsl->GetPixNum(); + + for(Int_t k = 0; kGetPixList(k); + MCerPhotPix &pix = (*fEvt)[idPix]; + pix.SetPixelUnused(); + } + } + } + } + + // + //eliminates the island with a time spread + //higher than a given limit + // + //else if( fIslandCleaningMethod == kTiming ){ + else if( fIslandCleaningMethod == 0 ){ + while ((imgIsl=(MImgIsland*)Next())) { + + //fIslCleanThreshold has different values, FIXME, put two variables + if(imgIsl->GetTimeSpread() > fIslCleanThres) + { + pixNum = imgIsl->GetPixNum(); + + for(Int_t k = 0; kGetPixList(k); + MCerPhotPix &pix = (*fEvt)[idPix]; + pix.SetPixelUnused(); + } + } + } + + } + + // + //eliminates all the islands except the continent, + //i.e. the larger island in the event + // + else if( fIslandCleaningMethod == 3 ){ + Int_t i = 0; + while ((imgIsl=(MImgIsland*)Next())) { + if (i != 0){ + + pixNum = imgIsl->GetPixNum(); + + for(Int_t k = 0; kGetPixList(k); + MCerPhotPix &pix = (*fEvt)[idPix]; + pix.SetPixelUnused(); + } + } + i++; + } + } + + fEvt->SetReadyToSave(); + + return kTRUE; + +} Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsClean.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsClean.h (revision 5170) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsClean.h (revision 5170) @@ -0,0 +1,60 @@ +#ifndef MARS_MIslandsClean +#define MARS_MIslandsClean + +#ifndef MARS_MGTask +#include "MGTask.h" +#endif + +#ifndef ROOT_TArrayF +#include +#endif + +class MGeomCam; +class MSigmabar; +class MCerPhotPix; +class MCerPhotEvt; +class MPedestalCam; +class MArrivalTimeCam; +class MArrivalTimePix; +class MIslands; +class MImgIsland; + +class MIslandsClean : public MGTask +{ + public: + + typedef enum { + kTiming, + kNoTiming + } IslandCleaningMethod_t; + + private: + const MGeomCam *fCam; //! + MCerPhotEvt *fEvt; //! + MSigmabar *fSgb; //! + MPedestalCam *fPed; //! + MArrivalTimeCam *fTime; //! + MIslands *fIsl; //! + + TString fIslName; + + // IslandCleaningMethod_t fIslandCleaningMethod; + Int_t fIslandCleaningMethod; + + Float_t fIslCleanThres; + + Int_t PreProcess(MParList *plist); + Int_t Process(); + + public: + MIslandsClean(const Float_t newThres=50, const char *name=NULL, const char *title=NULL); + + void SetInputName(TString inname) {fIslName = inname;} + + //void SetMethod(IslandCleaningMethod_t m) {fIslandCleaningMethod = m;} + void SetMethod(Int_t m) {fIslandCleaningMethod = m;} + + ClassDef(MIslandsClean, 0) // task doing the island cleaning +}; + +#endif