Index: trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 5349) +++ trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 5379) @@ -18,4 +18,14 @@ -*-*- END OF LINE -*-*- + + 2004/11/11 Ester Aliu + * library/MIslands.h + - Add the variable AlphaW + * library/MImgIsland.[h,cc] + - Add alpha for each island + * library/MIslandsCalc.cc + - Bug solved + * library/MIslandsClean.cc + - More different island cleanings added 2004/10/20 Javier Rico Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.cc (revision 5349) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.cc (revision 5379) @@ -16,5 +16,6 @@ ! ! Author(s): Javier Rico 04/2004 -! +! Ester Aliu 10/2004 (update according +! the new classes of the islands) ! Copyright: MAGIC Software Development, 2000-2004 ! Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc (revision 5349) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc (revision 5379) @@ -54,5 +54,4 @@ fPixList.Reset(-1); fPeakPulse.Reset(-1); - } @@ -105,4 +104,5 @@ *fLog << inf << " DistW = " << fDistW << endl; *fLog << inf << " DistS = " << fDistS << endl; + *fLog << inf << " Alpha = " << fAlpha << endl; } Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.h (revision 5349) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.h (revision 5379) @@ -22,12 +22,17 @@ 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; + + Float_t fWidth; + Float_t fLength; + Float_t fSizeIsl; + Float_t fMeanX; + Float_t fMeanY; + + Float_t fAlpha; + TArrayI fPixList; @@ -42,13 +47,19 @@ Float_t GetSigToNoise() { return fSigToNoise; } Float_t GetTimeSpread() { return fTimeSpread; } + Float_t GetDist() { return fDist; } + Float_t GetDistL() { return fDistL; } + Float_t GetDistW() { return fDistW; } + Float_t GetDistS() { return fDistS; } + + //hillas parameters + Float_t GetSizeIsl() { return fSizeIsl; } 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; } + // hillas src parameters + Float_t GetAlpha() { return fAlpha; } + void InitSize(Int_t i); UInt_t GetSize() const { return fPixList.GetSize(); } @@ -59,23 +70,29 @@ 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 SetPixNum (Int_t i) { fPixNum = i; } + void SetSigToNoise(Float_t val) { fSigToNoise = val; } + void SetTimeSpread(Float_t val) { fTimeSpread = val; } + void SetDist (Float_t val) { fDist = val; } + void SetDistL (Float_t val) { fDistL = val; } + void SetDistW (Float_t val) { fDistW = val; } + void SetDistS (Float_t val) { fDistS = val; } + + //hillas parameters + void SetSizeIsl (Float_t val) { fSizeIsl = val; } + void SetMeanX (Float_t val) { fMeanX = val; } + void SetMeanY (Float_t val) { fMeanY = val; } + void SetWidth (Float_t val) { fWidth = val; } + void SetLength (Float_t val) { fLength = val; } + + // hillas src parameters + void SetAlpha (Float_t val) { fAlpha = val; } - void SetPixList( const Int_t i,const Int_t id); - void SetPeakPulse( const Int_t i,const Float_t time); + 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 + ClassDef(MImgIsland, 2) // Container that holds the island information }; Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MIslands.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MIslands.h (revision 5349) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MIslands.h (revision 5379) @@ -23,4 +23,5 @@ Int_t fIslNum; + Float_t fAlphaW; public: @@ -34,7 +35,8 @@ TList* GetList() const { return fIslands; } UInt_t GetIslNum() const { return fIslands->GetSize(); } //number of islands + Float_t GetAlphaW() const { return fAlphaW; } //alpha weighted - void SetIslNum (Int_t i) { fIslNum = i; } - + void SetIslNum (Int_t i) { fIslNum = i; } + void SetAlphaW(Float_t val) { fAlphaW = val; } // MHCamera& GetDisplay() { return fDisplay; } Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc (revision 5349) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc (revision 5379) @@ -180,4 +180,5 @@ fIsl->SetIslNum(numisl); + //cout << "code numisl " << numisl << endl; //examine each island... @@ -196,9 +197,11 @@ Float_t distS[numisl]; - Float_t size[numisl], sizeLargeIsl, length, width; + Float_t size[numisl], sizeLargeIsl, length, width, distance, alpha, alphaW, sizetot; Float_t X = 0; Float_t Y = 0; sizeLargeIsl = 0; + alphaW = 0; + sizetot = 0; for(Int_t i = 1; i<=numisl ; i++) @@ -292,4 +295,5 @@ imgIsl->SetMeanY(meanY[i-1]); imgIsl->SetDist(dist[i-1]); + imgIsl->SetSizeIsl(size[i-1]); fIsl->GetList()->Add(imgIsl); @@ -333,5 +337,6 @@ MImgIsland *imgIsl = new MImgIsland; TIter Next(fIsl->GetList()); - + //Next.Reset(); + Int_t i = 1; while ((imgIsl=(MImgIsland*)Next())) { @@ -344,5 +349,9 @@ const Double_t s2 = tand2+1; - + const Double_t s = TMath::Sqrt(s2); + + const Double_t CosDelta = 1.0/s; // need these in derived classes + const Double_t SinDelta = tand/s; // like MHillasExt + const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/size[i-1]; const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/size[i-1]; @@ -359,16 +368,76 @@ width = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm] - + + // alpha calculation + + const Double_t mx = imgIsl->GetMeanX(); // [mm] + const Double_t my = imgIsl->GetMeanY(); // [mm] + + //FIXME: xpos, ypos from MSrcPos + const Double_t xpos = 0.; + const Double_t ypos = 0.; + + const Double_t sx = mx - xpos; // [mm] + const Double_t sy = my - ypos; // [mm] + + const Double_t sd = SinDelta; // [1] + const Double_t cd = CosDelta; // [1] + + // + // Distance from source position to center of ellipse. + // If the distance is 0 distance, Alpha is not specified. + // The calculation has failed and returnes kFALSE. + // + distance = TMath::Sqrt(sx*sx + sy*sy); // [mm] + if (distance==0) + return 1; + + // + // Calculate Alpha and Cosda = cos(d,a) + // The sign of Cosda will be used for quantities containing + // a head-tail information + // + // *OLD* const Double_t arg = (sy-tand*sx) / (dist*sqrt(tand*tand+1)); + // *OLD* fAlpha = asin(arg)*kRad2Deg; + // + const Double_t arg2 = cd*sx + sd*sy; // [mm] + if (arg2==0) + return 2; + + const Double_t arg1 = cd*sy - sd*sx; // [mm] + + // + // Due to numerical uncertanties in the calculation of the + // square root (dist) and arg1 it can happen (in less than 1e-5 cases) + // that the absolute value of arg exceeds 1. Because this uncertainty + // results in an Delta Alpha which is still less than 1e-3 we don't care + // about this uncertainty in general and simply set values which exceed + // to 1 saving its sign. + // + const Double_t arg = arg1/distance; + alpha = TMath::Abs(arg)>1 ? TMath::Sign(90., arg) : TMath::ASin(arg)*TMath::RadToDeg(); // [deg] + + alphaW += alpha*size[i-1]; + sizetot += size[i-1]; + + //////////////////////////////////////// + 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]); + + imgIsl->SetAlpha(alpha); i++; } - + + + fIsl->SetAlphaW(alphaW/sizetot); //fIsl->SetReadyToSave(); @@ -394,5 +463,5 @@ Int_t sflag; - Int_t control; + Int_t control = 0; Int_t nvect = 0; @@ -455,60 +524,49 @@ 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 ;jjGetTimeSpread() > fIslCleanThres) + { + pixNum = imgIsl->GetPixNum(); + + for(Int_t k = 0; kGetPixList(k); + MCerPhotPix &pix = (*fEvt)[idPix]; + pix.SetPixelUnused(); + } + } + } + + } + + + //////////////////////////////////////////////////// + // + // SIGNAL TO NOISE ISLAND CLEANING + // eliminates the island with a signal-to-noise + // lower than a given limit + /////////////////////////////////////////////////// + else if ( fIslandCleaningMethod == 1 ) { while ((imgIsl=(MImgIsland*)Next())) { @@ -173,14 +200,21 @@ } - // - //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) + + //////////////////////////////////////////////////// + // + // ISLAND SIZE OVER EVENT SIZE ISLAND CLEANING + // eliminates the island with an island size over event size + // lower than a given limit + /////////////////////////////////////////////////// + else if ( fIslandCleaningMethod == 2 ) { + Float_t size = 0; + while ((imgIsl=(MImgIsland*)Next())) { + size += imgIsl->GetSizeIsl(); + } + + Next.Reset(); + while ((imgIsl=(MImgIsland*)Next())) { + + if(imgIsl->GetSizeIsl()/size < fIslCleanThres) { pixNum = imgIsl->GetPixNum(); @@ -193,12 +227,15 @@ } } - } - - } - - // - //eliminates all the islands except the continent, - //i.e. the larger island in the event - // + } + } + + + //////////////////////////////////////////////////// + // + // CONTINENT ISLAND CLEANING + // eliminates all the islands except the continent + // i.e. the larger island in the event + // according the number of pixels + /////////////////////////////////////////////////// else if( fIslandCleaningMethod == 3 ){ Int_t i = 0; @@ -219,4 +256,115 @@ } + + //////////////////////////////////////////////////// + // + // LARGER and SECOND LARGER ISLAND CLEANING + // eliminates all the islands except the two biggest + // ones according size + // + /////////////////////////////////////////////////// + else if( fIslandCleaningMethod == 4 ){ + + Int_t i = 0; + Int_t islnum = fIsl->GetIslNum(); + Float_t islSize[islnum]; + Int_t islIdx[islnum]; + + for (Int_t j = 0; jGetSizeIsl(); + i++; + } + + + TMath::Sort(islnum, islSize, islIdx, kTRUE); + + i = 0; + Next.Reset(); + while ((imgIsl=(MImgIsland*)Next())) { + if (islnum > 1 && islIdx[0]!=i && islIdx[1]!=i){ + + //cout << "removed " << i << " isl 0" << islIdx[0] << " isl 1" << islIdx[1] << endl; + + pixNum = imgIsl->GetPixNum(); + + for(Int_t k = 0; kGetPixList(k); + MCerPhotPix &pix = (*fEvt)[idPix]; + pix.SetPixelUnused(); + } + } + i++; + } + } + + + + /////////////////////////////////////////////////////// + // + // LARGER and SECOND LARGER ISLAND CLEANING II + // eliminates all the islands except the two biggest + // ones according size, if the second one almost has + // the 80% of the size of the biggest one + // + // + ////////////////////////////////////////////////////// + else if( fIslandCleaningMethod == 5 ){ + + Int_t i = 0; + Int_t islnum = fIsl->GetIslNum(); + Float_t islSize[islnum]; + Int_t islIdx[islnum]; + + for (Int_t j = 0; jGetSizeIsl(); + i++; + } + + + TMath::Sort(islnum, islSize, islIdx, kTRUE); + + i = 0; + Next.Reset(); + while ((imgIsl=(MImgIsland*)Next())) { + + if (islnum > 1 && islIdx[0]!=i && islIdx[1]!=i){ + + pixNum = imgIsl->GetPixNum(); + + for(Int_t k = 0; kGetPixList(k); + MCerPhotPix &pix = (*fEvt)[idPix]; + pix.SetPixelUnused(); + } + } + else if(islnum>1 && islSize[islIdx[1]]<0.6*islSize[islIdx[0]]){ + + pixNum = imgIsl->GetPixNum(); + + for(Int_t k = 0; kGetPixList(k); + MCerPhotPix &pix = (*fEvt)[idPix]; + pix.SetPixelUnused(); + } + } + i++; + } + } + fEvt->SetReadyToSave();