Changeset 701 for trunk/MagicSoft

Timestamp:

03/22/01 16:28:06 (24 years ago)

Author:

tbretz

Message:

fixed some minor errors in Hillas calculation

File:

• trunk/MagicSoft/Mars/manalysis/MHillas.cc (modified) (3 diffs)

trunk/MagicSoft/Mars/manalysis/MHillas.cc

@@ -55 +55 @@
 
     fEllipse = new TEllipse(cos(fTheta)*fDist, sin(fTheta)*fDist,
-                            fWidth, fLength,
-                            0, 360, fAlpha);
+                            fLength, fWidth,
+                            0, 360, fTheta*kRad2Deg+fAlpha-180);
 
     fEllipse->SetLineWidth(2);
@@ -137 +137 @@
     float s = sin(theta); // [1]
 
-    const float direction = c*xmean + s*ymean;
-
-    if (direction<0)
-    {
-        c = -c;
-        s = -s;
-    }
-
+    //
+    // resolve four-fold ambiguity of solution
+    //
     float axis1 =  2.0*c*s*sigmaxy + c*c*sigmaxx + s*s*sigmayy; // [mm^2]
     float axis2 = -2.0*c*s*sigmaxy + s*s*sigmaxx + c*c*sigmayy; // [mm^2]
@@ -158 +153 @@
 
     //
-    // check the rotation of the axis
+    // check the rotation of the axis (maybe turn by 90ø)
     //
-    const int rotation = axis1<axis2;
+    const int   rotation = axis1<axis2;
 
-    fLength = rotation ? sqrt(axis2) : sqrt(axis1);      // [mm]
-    fWidth  = rotation ? sqrt(axis1) : sqrt(axis2);      // [mm]
+    fLength = rotation ? sqrt(axis2) : sqrt(axis1);     // [mm]
+    fWidth  = rotation ? sqrt(axis1) : sqrt(axis2);     // [mm]
 
-    fAlpha = 180/kPI*atan((-xmean*s+ymean*c)/direction); // [deg]
+    fAlpha  = rotation ?
+        fabs(atan((-xmean*c - ymean*s)/(c*ymean - s*xmean))) :
+        fabs(atan(( ymean*c - xmean*s)/(c*xmean + s*ymean))) ;
 
-    fDist  = sqrt(xmean*xmean + ymean*ymean);            // [mm]
+        // [deg]
 
-    fTheta = atan(ymean/xmean);                          // [rad]
-    if (xmean<0) fTheta += kPI;                          // [deg]
+    fAlpha *= kRad2Deg;                                  // [deg]
+
+    fDist   = sqrt(xmean*xmean + ymean*ymean);           // [mm]
+
+    fTheta  = atan(ymean/xmean);                         // [rad]
+    if (xmean<0) fTheta += kPI;                          // [rad]
 }