Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc	(revision 5394)
+++ 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 5394)
+++ 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 ; idx<nPix ; idx++)
 	{
@@ -247,10 +253,10 @@
 	    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]);     	  	
@@ -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 ; idx<nPix ; idx++){
+	
+	MCerPhotPix *pix = fEvt->GetPixById(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 ; idx<nPix ; idx++)
-    {
-      MCerPhotPix *pix = fEvt->GetPixById(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<MCerPhotPix*>(Next())))
     {
@@ -493,5 +511,5 @@
       if( fEvt->IsPixelUsed(idx)) 
      	{
-	  //cout << idx <<endl;
+	  //cout <<idx <<endl;
 	  sflag = 0;
 	  
@@ -650,6 +668,9 @@
   MCerPhotPix *pix;
 
-  //1st loop over all pixels
-  MCerPhotEvtIter Next0(fEvt, kFALSE);
+  //after the interpolation of the pixels, these can be disordered by index. This is important for this algorithm of calculating islands
+  fEvt->Sort();
+
+  // 1st loop over used pixels only
+  TIter Next0(*fEvt);
  
   while ((pix=static_cast<MCerPhotPix*>(Next0())))
@@ -671,7 +692,8 @@
     }
       
+ 
   //2nd loop over all pixels
-  MCerPhotEvtIter Next(fEvt, kFALSE);
-  
+  TIter Next(*fEvt);
+   
   while ((pix=static_cast<MCerPhotPix*>(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<j) {
-	  ii=i;
-	  jj=j;
-	}
-	else{
-	  ii=j;
-	  jj=i;
-      }
-	if (control == 1) {
-	  for(Int_t k = 0; k < npix; k++) {
-	    if(vect[jj][k] == 1) vect[ii][k] = 1;
-	    vect[jj][k] = 0;
-	    zeros[jj] = 1;
-	  }   
-	  numisl = numisl-1;       
-	}
-      }
-    }
-  }
+ // Repeated Chain Corrections
+
+  Int_t jmin = 0;
+  
+  for(Int_t i = 1; i <= nvect; i++){
+    control=0;
+    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;
+	  k=npix;
+	}
+      }
+      if (control == 1){
+	for(Int_t k = 0; k < npix; k++){
+	  if(vect[j][k] == 1) vect[i][k] = 1;
+	  vect[j][k] = 0;
+	  zeros[j] = 1;
+	}	
+	numisl = numisl-1;	    
+      }
+    }
+    
+    for(Int_t j = 1; j <= i-1; j++){
+      for(Int_t k = 0; k < npix; k++){
+	if (vect[i][k] == 1 && vect[j][k] == 1){
+	  control = 2; 
+	  jmin=j;
+	  k=npix;
+	  j=i;
+	}
+      }
+      
+      if (control == 2){
+	for (Int_t k = 0; k < npix; k++){
+	  if(vect[i][k]==1) vect[jmin][k]=1;
+	  vect[i][k] = 0;
+	  zeros[i] = 1;
+	}
+	numisl = numisl-1;	    
+      }	   
+    } 
+  } 
   
   Int_t l = 1;
