Changeset 5466 for trunk/MagicSoft/Mars/mtemp/mifae/library

Timestamp:

11/24/04 16:33:39 (20 years ago)

Author:

aliu

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mtemp/mifae/library

Files:

• MImgIsland.cc (modified) (2 diffs)
• MIslandsCalc.cc (modified) (12 diffs)

trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc

@@ -37 +37 @@
     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 +48 @@
 {
   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);
+   
 }
 

trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc

@@ -35 +35 @@
 //   - 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 +192 @@
   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 +199 @@
   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 +212 @@
       
       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 +253 @@
             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 +265 @@
             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 +306 @@
       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 +496 @@
   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 +511 @@
       if( fEvt->IsPixelUsed(idx)) 
         {
-          //cout << idx <<endl;
+          //cout <<idx <<endl;
           sflag = 0;
           
@@ -650 +668 @@
   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 +692 @@
     }
       
+ 
   //2nd loop over all pixels
-  MCerPhotEvtIter Next(fEvt, kFALSE);
-  
+  TIter Next(*fEvt);
+   
   while ((pix=static_cast<MCerPhotPix*>(Next())))
     {
@@ -713 +735 @@
   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;