Index: /trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 5465) +++ /trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 5466) @@ -18,4 +18,12 @@ -*-*- END OF LINE -*-*- + + 2004/11/24 Ester Aliu + * library/MIslandsCalc.h + - Remove a memmory leak, and add a sort function to order the + pixels by index after their interpolation + - Other minor changes + * library/MImgIsland.[h,cc] + - Reset some variables 2004/11/11 Ester Aliu Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc (revision 5465) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc (revision 5466) @@ -37,6 +37,10 @@ fName = name ? name : "MImgIsland"; fTitle = title ? title : "Storage container for one island information"; - + // cout << "Creo 1 isla" << endl; Reset(); +} +MImgIsland::~MImgIsland() +{ + // cout << "Destruyo 1 isla" << endl; } @@ -44,14 +48,21 @@ { fPixNum = 0; - fSigToNoise = 0; - fTimeSpread = 0; - fMeanX = 0; - fMeanY = 0; - fDist = 0; - fDistW = 0; - fDistL = 0; + fSigToNoise = -1; + fTimeSpread = -1; + fMeanX = -1000; + fMeanY = -1000; + fDist = -1; + fDistW = -1; + fDistL = -1; + fDistS = -1; + + fWidth = -1; + fLength = -1; + fSizeIsl = -1; + fAlpha = -1000; fPixList.Reset(-1); fPeakPulse.Reset(-1); + } Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc (revision 5465) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc (revision 5466) @@ -35,10 +35,10 @@ // - fPixNum // number of pixels in the island // - fSigToNoise // signal to noise of the island -// - fTimeSpread // mean arrival time spread of the island +// - fTimeSpread // "time" 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 +// - fLength // major axis of the island ellipse +// - fWidth // minor axis of the island ellipse // - fDistL // dist divided by lenght of the larger island // - fDistW // dist divided by width of the larger island @@ -192,4 +192,6 @@ Float_t meanX[numisl]; Float_t meanY[numisl]; + Float_t length[numisl]; + Float_t width[numisl]; Float_t dist[numisl]; Float_t distL[numisl]; @@ -197,8 +199,7 @@ Float_t distS[numisl]; - Float_t size[numisl], sizeLargeIsl, length, width, distance, alpha, alphaW, sizetot; - - Float_t X = 0; - Float_t Y = 0; + + Float_t size[numisl], sizeLargeIsl, distance, alpha, alphaW, sizetot; + sizeLargeIsl = 0; alphaW = 0; @@ -211,22 +212,27 @@ imgIsl->InitSize(num[i]); - + Int_t n = 0; - - Float_t MIN = 10000.; - Float_t MAX = 0.; - + + Float_t minTime = 10000.; + Float_t maxTime = 0.; + + Float_t minX = 10000.; + Float_t maxX = 0.; + + Float_t minY = 10000.; + Float_t maxY = 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 ; idxSetPixList(PixelNumIsl[i-1]-1,pix->GetPixId()); imgIsl->SetPeakPulse(PixelNumIsl[i-1]-1,time[i-1]); @@ -259,34 +265,38 @@ if (fEvt->IsPixelCore(idx)){ - if (time[i-1] > MAX) - MAX = time[i-1]; - if (time[i-1] < MIN) - MIN = time[i-1]; - - } - + if (time[i-1] > maxTime){ + maxTime = time[i-1]; + maxX = gpix2.GetX(); + maxY = gpix2.GetY(); + } + if (time[i-1] < minTime){ + minTime = time[i-1]; + minX = gpix2.GetX(); + minY = gpix2.GetY(); + } + + } 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::Power(meanX[i-1]-meanX[0],2) + TMath::Power(meanY[i-1]-meanY[i-1],2); dist[i-1] = TMath::Sqrt(dist[i-1]); - - timeVariance[i-1] = MAX-MIN; + + //timeVariance[i-1] = (maxTime-minTime); + + if (maxX!=minX && maxY!=minY) + timeVariance[i-1] = (maxTime-minTime)/sqrt(TMath::Power(maxX-minX,2) + TMath::Power(maxY-minY,2)); + else + timeVariance[i-1] = -1; //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]); @@ -296,154 +306,159 @@ imgIsl->SetDist(dist[i-1]); imgIsl->SetSizeIsl(size[i-1]); + + + // sanity check: if one island has 2 or less pixels + if (num[i]>2){ + + + // calculate width and lenght of each island + + 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 (vect[i][idx]==1){ + + const Float_t dx = gpix3.GetX() - meanX[i-1]; // [mm] + const Float_t dy = gpix3.GetY() - meanY[i-1]; // [mm] + + corrxx += nphot * dx*dx; // [mm^2] + corrxy += nphot * dx*dy; // [mm^2] + corryy += nphot * dy*dy; // [mm^2] + + } + } + + 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 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]; + + // + // 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[i-1] = axis1<0 ? 0 : TMath::Sqrt(axis1); // [mm] + width[i-1] = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm] + + + // alpha calculation + + const Double_t mx = meanX[i-1]; // [mm] + const Double_t my = meanY[i-1]; // [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){ + + for (Int_t l = 0; l< nVect; l++) + delete [] vect[l]; + + delete [] vect; + delete [] num; + + 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){ + + for (Int_t l = 0; l< nVect; l++) + delete [] vect[l]; + + delete [] vect; + delete [] num; + + 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[0]; + distW[i-1]=dist[i-1]/width[0]; + distS[i-1]= dist[i-1]/size[0]; + + imgIsl->SetLength(length[i-1]); + imgIsl->SetWidth(width[i-1]); + + imgIsl->SetDistL(distL[i-1]); + imgIsl->SetDistW(distW[i-1]); + imgIsl->SetDistS(distS[i-1]); + + imgIsl->SetAlpha(alpha); + + } 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()); - //Next.Reset(); - - 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 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]; - - // - // 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] - - - // 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(); - - for (Int_t i = 0; i< nVect; i++) - delete [] vect[i]; - - delete vect; + + for (Int_t l = 0; l< nVect; l++) + delete [] vect[l]; + + delete [] vect; + delete [] num; return kTRUE; @@ -481,7 +496,10 @@ MCerPhotPix *pix; - //loop over all pixels - MCerPhotEvtIter Next(fEvt, kFALSE); - + // Loop over used pixels only + TIter Next(*fEvt); + + //after the interpolation of the pixels, these can be disordered by index. This is important for this algorithm of calculating islands + fEvt->Sort(); + while ((pix=static_cast(Next()))) { @@ -493,5 +511,5 @@ if( fEvt->IsPixelUsed(idx)) { - //cout << idx <Sort(); + + // 1st loop over used pixels only + TIter Next0(*fEvt); while ((pix=static_cast(Next0()))) @@ -671,7 +692,8 @@ } + //2nd loop over all pixels - MCerPhotEvtIter Next(fEvt, kFALSE); - + TIter Next(*fEvt); + while ((pix=static_cast(Next()))) { @@ -713,38 +735,48 @@ numisl = nvect; - - // Repeated Chain Corrections - - Int_t ii, jj; - - for(Int_t i = 1; i <= nvect; i++) { - for(Int_t j = 1; j <= nvect; j++) { - if (i!=j && zeros[j]!=1){ - control = 0; - for(Int_t k = 0; k < npix; k++) { - if (vect[i][k] == 1 && vect[j][k] == 1) { - control = 1; - break; - } - } - if(i