Changeset 5379 for trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc

Timestamp:

11/11/04 09:55:25 (20 years ago)

Author:

aliu

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc (modified) (10 diffs)

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

@@ -180 +180 @@
   fIsl->SetIslNum(numisl);
 
+  //cout << "code numisl " << numisl << endl;
   //examine each island...
 
@@ -196 +197 @@
   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 +295 @@
       imgIsl->SetMeanY(meanY[i-1]);
       imgIsl->SetDist(dist[i-1]);
+      imgIsl->SetSizeIsl(size[i-1]);
 
       fIsl->GetList()->Add(imgIsl);
@@ -333 +337 @@
   MImgIsland *imgIsl = new MImgIsland;
   TIter Next(fIsl->GetList());
-  
+  //Next.Reset();
+
   Int_t i = 1;
   while ((imgIsl=(MImgIsland*)Next())) {
@@ -344 +349 @@
     
     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 +368 @@
     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 +463 @@
   
   Int_t    sflag;
-  Int_t    control;
+  Int_t    control = 0;
   
   Int_t    nvect = 0;
@@ -455 +524 @@
   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 ;jj<i ; jj++){
-                 
-                  if(vect[jj][k]==1){ 
-                    jmin = jj;
-                    control = 2;
-                    break;
-                  }
-                }
-              }
-            }
-          
-          
-          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;            
-            }
-
-          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;            
-            }
-        }
-    }
+  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 pixMAX = 0;
@@ -540 +598 @@
         }
     }
-  
+
+
   //the larger island will correspond to the 1st component of the vector
   
@@ -657 +716 @@
   // 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 ;jj<i ; jj++){
-                  
-                  if(vect[jj][k]==1){ 
-                    jmin = jj;
-                    control = 2;
-                    break;
-                  }
-                }
-              }
-            }
-
-       
-          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;            
-            }
-
-          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 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;       
+        }
+      }
+    }
+  }
   
   Int_t l = 1;