Changeset 14581 for trunk/Cosy/tpoint/TPointStar.cc

Timestamp:

11/08/12 11:48:44 (12 years ago)

Author:

tbretz

Message:

Fixed the GetResidual.

File:

• trunk/Cosy/tpoint/TPointStar.cc (modified) (1 diff)

trunk/Cosy/tpoint/TPointStar.cc

@@ -61 +61 @@
 
 Double_t TPointStar::GetResidual(Double_t *err) const
+{
+    /*
+     TVector3 v1, v2;
+     v1.SetMagThetaPhi(1, TMath::Pi()/2-fRawEl, fRawAz);
+     v2.SetMagThetaPhi(1, TMath::Pi()/2-fStarEl, fStarAz);
+
+     return v1.Angle(v2)*TMath::RadToDeg();
+     */
+
+    const Double_t del = fRawEl-fStarEl;
+    const Double_t daz = fRawAz-fStarAz;
+
+    const double x = sin(fRawEl) * sin(fStarEl) * cos(fStarAz-fRawAz);
+    const double y = cos(fRawEl) * cos(fStarEl);
+
+    const Double_t d = acos(x + y);
+
+    if (err)
     {
-      /*
-       TVector3 v1, v2;
-       v1.SetMagThetaPhi(1, TMath::Pi()/2-fRawEl, fRawAz);
-       v2.SetMagThetaPhi(1, TMath::Pi()/2-fStarEl, fStarAz);
+        // Error of one pixel in the CCD
+        const Double_t e1 = 32./3600*TMath::DegToRad()   * 0.5;
 
-       return v1.Angle(v2)*TMath::RadToDeg();
-       */
+        // Error of one SE unit
+        const Double_t e2 = 360./16384*TMath::DegToRad() * 0.5;
 
-        const Double_t del = fRawEl-fStarEl;
-        const Double_t daz = fRawAz-fStarAz;
+        const Double_t e11 =  sin(del)+cos(fRawEl)*sin(fStarEl)*(1-cos(daz));
+        const Double_t e12 =  cos(fRawEl)*cos(fStarEl)*sin(daz);
 
-        /*
-        const Double_t dphi2 = daz/2.;
-        const Double_t cos2  = cos(dphi2)*cos(dphi2);
-        const Double_t sin2  = sin(dphi2)*sin(dphi2);
-        const Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
-        */
-        const Double_t d  = cos(del) - cos(fRawEl)*cos(fStarEl)*(1.-cos(daz));
+        const Double_t e21 = -sin(del)+sin(fRawEl)*cos(fStarEl)*(1-cos(daz));
+        const Double_t e22 = -cos(fRawEl)*cos(fStarEl)*sin(daz);
 
-        if (err)
-        {
-            // Error of one pixel in the CCD
-            const Double_t e1 = 32./3600*TMath::DegToRad()   * 0.5;
+        const Double_t err1  = sqrt(1-d*d);
+        const Double_t err2  = (e11*e11 + e12*e12)*e1*e1;
+        const Double_t err3  = (e21*e21 + e22*e22)*e2*e2;
 
-            // Error of one SE unit
-            const Double_t e2 = 360./16384*TMath::DegToRad() * 0.5;
+        *err = sqrt(err2+err3)/err1 * TMath::RadToDeg();
+    }
 
-            const Double_t e11 =  sin(del)+cos(fRawEl)*sin(fStarEl)*(1-cos(daz));
-            const Double_t e12 =  cos(fRawEl)*cos(fStarEl)*sin(daz);
-
-            const Double_t e21 = -sin(del)+sin(fRawEl)*cos(fStarEl)*(1-cos(daz));
-            const Double_t e22 = -cos(fRawEl)*cos(fStarEl)*sin(daz);
-
-            const Double_t err1  = sqrt(1-d*d);
-            const Double_t err2  = (e11*e11 + e12*e12)*e1*e1;
-            const Double_t err3  = (e21*e21 + e22*e22)*e2*e2;
-
-            *err = sqrt(err2+err3)/err1 * TMath::RadToDeg();
-        }
-
-        const Double_t dist = acos(d);
-        return dist * TMath::RadToDeg();
-    }
+    const Double_t dist = acos(d);
+    return dist * TMath::RadToDeg();
+}
 
 istream &operator>>(istream &fin, TPointStar &set)