Changeset 4107 for trunk/MagicSoft/Cosy/tpoint

Timestamp:

05/20/04 05:40:59 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Cosy/tpoint/gui.C (modified) (23 diffs)

trunk/MagicSoft/Cosy/tpoint/gui.C

@@ +1 @@
+#include <fstream.h>
+#include <fstream.h>
 #include <fstream.h>
 #include <iostream.h>
 #include <iomanip.h>
+
+#include <TError.h>
 
 #include <TGFrame.h>
@@ -38 +42 @@
 public:
     Set(Double_t sel=0, Double_t saz=0, Double_t rel=0, Double_t raz=0) :
-        fStarAz(saz*TMath::Pi()/180),
-        fStarEl(sel*TMath::Pi()/180),
-        fRawAz(raz*TMath::Pi()/180),
-        fRawEl(rel*TMath::Pi()/180)
-    {
-    }
-
-    Double_t GetResidual() const
-    {
-         Double_t del = fRawEl-fStarEl;
-         Double_t daz = fRawAz-fStarAz;
-         Double_t dphi2 = daz/2.;
-         Double_t cos2  = cos(dphi2)*cos(dphi2);
-         Double_t sin2  = sin(dphi2)*sin(dphi2);
-         Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
-
-         Double_t dist = acos(d);
-
-         return dist * 180 / TMath::Pi();
+        fStarAz(saz*TMath::DegToRad()),
+        fStarEl(sel*TMath::DegToRad()),
+        fRawAz(raz*TMath::DegToRad()),
+        fRawEl(rel*TMath::DegToRad())
+    {
+    }
+
+    Double_t GetResidual(Double_t *err=0) 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_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));
+
+        if (err)
+        {
+            // Error of one pixel in the CCD
+            const Double_t e1 = 32./3600*TMath::DegToRad()   * 0.5;
+
+            // Error of one SE unit
+            const Double_t e2 = 360./16384*TMath::DegToRad() * 0.5;
+
+            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();
     }
 
@@ -67 +103 @@
     }
 
-    Double_t GetDEl() const     { return (fRawEl-fStarEl)*180/TMath::Pi(); }
+    Double_t GetDEl() const     { return (fRawEl-fStarEl)*TMath::RadToDeg(); }
     Double_t GetDZd() const     { return -GetDEl(); }
-    Double_t GetDAz() const     { return (fRawAz-fStarAz)*180/TMath::Pi(); }
-    Double_t GetStarEl() const  { return fStarEl*180/TMath::Pi(); }
-    Double_t GetStarZd() const  { return 90.-fStarEl*180/TMath::Pi(); }
-    Double_t GetStarAz() const  { return fStarAz*180/TMath::Pi(); }
-    Double_t GetRawEl() const   { return fRawEl*180/TMath::Pi(); }
-    Double_t GetRawAz() const   { return fRawAz*180/TMath::Pi(); }
-    Double_t GetRawZd() const   { return 90.-fRawEl*180/TMath::Pi(); }
+    Double_t GetDAz() const     { return (fRawAz-fStarAz)*TMath::RadToDeg(); }
+    Double_t GetStarEl() const  { return fStarEl*TMath::RadToDeg(); }
+    Double_t GetStarZd() const  { return 90.-fStarEl*TMath::RadToDeg(); }
+    Double_t GetStarAz() const  { return fStarAz*TMath::RadToDeg(); }
+    Double_t GetRawEl() const   { return fRawEl*TMath::RadToDeg(); }
+    Double_t GetRawAz() const   { return fRawAz*TMath::RadToDeg(); }
+    Double_t GetRawZd() const   { return 90.-fRawEl*TMath::RadToDeg(); }
 
     ZdAz  GetStarZdAz() const   { return ZdAz(TMath::Pi()/2-fStarEl, fStarAz); }
@@ -83 +119 @@
     AltAz GetRawAltAz() const   { return AltAz(fRawEl, fRawAz); }
 
-    void AdjustEl(Double_t del) { fStarEl += del*TMath::Pi()/180; }
-    void AdjustAz(Double_t daz) { fStarAz += daz*TMath::Pi()/180; }
+    void AdjustEl(Double_t del) { fStarEl += del*TMath::DegToRad(); }
+    void AdjustAz(Double_t daz) { fStarAz += daz*TMath::DegToRad(); }
 
     void Adjust(const MBending &bend)
@@ -105 +141 @@
 ifstream &operator>>(ifstream &fin, Set &set)
 {
-    Double_t v[4];
-    fin >> v[0];
-    fin >> v[1];
-    fin >> v[2];
-    fin >> v[3];
-
-    Double_t dummy;
-    fin>>dummy;
-    fin>>dummy;
-    fin>>dummy;
-
-    set.fStarAz = v[0]*TMath::Pi()/180;
-    set.fStarEl = v[1]*TMath::Pi()/180;
-
-    set.fRawAz  = v[2]*TMath::Pi()/180;
-    set.fRawEl  = v[3]*TMath::Pi()/180;
+    TString str;
+    str.ReadLine(fin);
+
+    Float_t v[4];
+    sscanf(str.Data(), "%f %f %f %f", v, v+1, v+2, v+3);
+
+    set.fStarAz = v[0]*TMath::DegToRad();
+    set.fStarEl = v[1]*TMath::DegToRad();
+
+    set.fRawAz  = v[2]*TMath::DegToRad();
+    set.fRawEl  = v[3]*TMath::DegToRad();
+
+    if (fin)
+    {
+        Double_t res, err;
+        res = set.GetResidual(&err);
+        cout << "Read: " << v[0] << " " << v[1] << "  :  " << v[2] << " " << v[3] << "  :  " << v[2]-v[0] << " " << v[3]-v[1] << "  :  " << res << " " << err << " " << err/res << endl;
+    }
 
     return fin;
@@ -160 +198 @@
             set.Adjust(bend);
 
-            Double_t err = 0.02; // [deg] = 1SE
-            Double_t res = set.GetResidual()/err;
+            Double_t err;// = 0.005; // [deg] = 0.25SE
+            Double_t res = set.GetResidual(&err);
+            res /= err;
 
             f += res*res;
@@ -168 +207 @@
         //f /= (fCoordinates.GetSize()-7)*(fCoordinates.GetSize()-7);
         //f /= fCoordinates.GetSize()*fCoordinates.GetSize();
+        //cout << f << ": " << fCoordinates.GetSize() << endl;
         f /= fCoordinates.GetSize();
     }
@@ -187 +227 @@
 
         TMarker mark0;
-        TMarker mark1;
+        //TMarker mark1;
         mark0.SetMarkerStyle(kStar);
-        mark1.SetMarkerStyle(kStar);
-        mark1.SetMarkerColor(kRed);
+        //mark1.SetMarkerStyle(kStar);
+        //mark1.SetMarkerColor(kRed);
     
         r0 /= 90;
-        r1 /= 90;
-        phi0 *= TMath::Pi()/180;
-        phi1 *= TMath::Pi()/180;
-    
+        //r1 /= 90;
+        phi0 *= TMath::DegToRad();
+        //phi1 *= TMath::DegToRad();
+
         Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
-        Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
-    
+        //Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
+
         Double_t y0[3], y1[3];
     
         view->WCtoNDC(x0, y0);
-        view->WCtoNDC(x1, y1);
+        //view->WCtoNDC(x1, y1);
     
         mark0.DrawMarker(y0[0], y0[1]);
-        mark1.DrawMarker(y1[0], y1[1]);
+        //mark1.DrawMarker(y1[0], y1[1]);
     }
     
@@ -256 +296 @@
         Double_t dp = p1-p0;
     
-        Double_t x0[3] = { r0*cos(p0*TMath::Pi()/180), r0*sin(p0*TMath::Pi()/180), 0};
+        Double_t x0[3] = { r0*cos(p0*TMath::DegToRad()), r0*sin(p0*TMath::DegToRad()), 0};
     
         for (double i=p0+10; i<p1+10; i+=10)
@@ -264 +304 @@
     
             Double_t r = dr/dp*(i-p0)+r0;
-            Double_t p = TMath::Pi()/180*i;
+            Double_t p = TMath::DegToRad()*i;
     
             Double_t x1[3] = { r*cos(p), r*sin(p), 0};
@@ -280 +320 @@
     }
     
-    void DrawSet(TVirtualPad *pad, Set &set, Float_t scale=-1)
+    void DrawSet(TVirtualPad *pad, Set &set, Float_t scale=-1, Float_t angle=0)
     {
         Double_t r0   = set.GetRawZd();
-        Double_t phi0 = set.GetRawAz();
+        Double_t phi0 = set.GetRawAz()-angle;
         Double_t r1   = set.GetStarZd();
-        Double_t phi1 = set.GetStarAz();
+        Double_t phi1 = set.GetStarAz()-angle;
     
         if (r0<0)
@@ -323 +363 @@
     }
 
-    void Fit()
+    void Fit(Double_t &before, Double_t &after, Double_t &backw)
     {
         if (fOriginal.GetSize()==0)
@@ -336 +376 @@
 
         cout << "-----------------------------------------------------------------------" << endl;
-    
-        TH1F hres1("Res1", " Residuals before correction ", fOriginal.GetSize()/2, 0, 3);
-        TH1F hres2("Res2", " Residuals after correction ",  fOriginal.GetSize(), 0, 3);
-        TH1F hres3("Res3", " Residuals after backward correction ",  fOriginal.GetSize(), 0, 3);
+
+        gStyle->SetOptStat("emro");
+
+        TH1F hres1("Res1", " Residuals before correction ", fOriginal.GetSize()/3, 0, 0.3);
+        TH1F hres2("Res2", " Residuals after correction ",  fOriginal.GetSize()/3, 0, 0.3);
+        TH1F hres3("Res3", " Residuals after backward correction ",  fOriginal.GetSize()/3, 0, 0.3);
     
         hres1.SetXTitle("\\Delta [\\circ]");
@@ -505 +547 @@
         cout << endl;
     
-
-
 
         TCanvas *c1;
@@ -576 +616 @@
         cout << "before: " << hres1.GetMean() << " \xb1 " << hres1.GetRMS() << " deg " << endl;
         cout << "after:  " << hres2.GetMean() << " \xb1 " << hres2.GetRMS() << " deg " << endl;
-        cout << "before: " << (int)(hres1.GetMean()*60+.5) << " \xb1 " << (int)(hres1.GetRMS()*60+.5) << " arcmin" << endl;
-        cout << "after:  " << (int)(hres2.GetMean()*60+.5) << " \xb1 " << (int)(hres2.GetRMS()*60+.5) << " arcmin" << endl;
-        cout << "before: " << (int)(hres1.GetMean()*60*60/23.4+.5) << " \xb1 " << (int)(hres1.GetRMS()*60*60/23.4+.5) << " pix" << endl;
-        cout << "after:  " << (int)(hres2.GetMean()*60*60/23.4+.5) << " \xb1 " << (int)(hres2.GetRMS()*60*60/23.4+.5) << " pix" << endl;
-        cout << "after:  " << (int)(hres2.GetMean()*16384/360+.5) << " \xb1 " << (int)(hres2.GetRMS()*16384/360+.5) << " SE" << endl;
         cout << endl;
-        cout << "backw:  " << (int)(hres3.GetMean()*60+.5) << " \xb1 " << (int)(hres3.GetRMS()*60+.5) << " arcmin" << endl;
+        cout << "before: " << Form("%.1f", hres1.GetMean()*60) << " \xb1 " << Form("%.1f", hres1.GetRMS()*60) << " arcmin" << endl;
+        cout << "after:  " << Form("%.1f", hres2.GetMean()*60) << " \xb1 " << Form("%.1f", hres2.GetRMS()*60) << " arcmin" << endl;
+        cout << endl;
+        cout << "before: " << Form("%.1f", hres1.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres1.GetRMS()*60*60/23.4) << " pix" << endl;
+        cout << "after:  " << Form("%.1f", hres2.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres2.GetRMS()*60*60/23.4) << " pix" << endl;
+        cout << "after:  " << Form("%.1f", hres2.GetMean()*16384/360) << " \xb1 " << Form("%.1f", hres2.GetRMS()*16384/360) << " SE" << endl;
+        cout << endl;
+        cout << "backw:  " << Form("%.1f", hres3.GetMean()*60) << " \xb1 " << Form("%.1f", hres3.GetRMS()*60) << " arcmin" << endl;
         cout << endl;                                            // ±
     
-    
+
+        before = hres1.GetMean()*16384/360;
+        after  = hres2.GetMean()*16384/360;
+        backw  = hres3.GetMean()*16384/360;
+
     
         gStyle->SetOptStat(1110);
@@ -615 +661 @@
         cout << "Gaus-Fit  Sigma: " << f.GetParameter(2) << "\xb0" << endl;
         cout << "Fit-Probability: " << f.GetProb()*100 << "%" << endl;
-        cout << "      Chi^2/NDF: " << f.GetChisquare()/f.GetNDF() << endl;
+        cout << "      Chi^2/NDF: " << f.GetChisquare() << "/" << f.GetNDF() << " = " << f.GetChisquare()/f.GetNDF() << endl;
 
         c1->cd(1);
@@ -626 +672 @@
         gPad->Modified();
         gPad->Update();
-        for (int i=0; i<fCoordinates.GetSize(); i++)
-            DrawSet(gPad, *(Set*)list.At(i));//, 10./hres1.GetMean());
+        for (int i=0; i<fOriginal.GetSize(); i++)
+            DrawSet(gPad, *(Set*)fOriginal.At(i));//, 10./hres1.GetMean());
     
         c1->cd(3);
@@ -638 +684 @@
         gPad->Update();
         for (int i=0; i<fCoordinates.GetSize(); i++)
-            DrawSet(gPad, *(Set*)fCoordinates.At(i), 10./hres2.GetMean());
+            DrawSet(gPad, *(Set*)fCoordinates.At(i), 10./hres2.GetMean(), par[0]);
 
         RaiseWindow();
@@ -663 +709 @@
 
             fin >> set;  // Read data from file [deg], it is stored in [rad]
-
             if (!fin)
                 break;
@@ -687 +732 @@
                 {
                 case kTbFit:
-                    Fit();
-                    DisplayBending();
+                    {
+                        Double_t before=0;
+                        Double_t after=0;
+                        Double_t backw=0;
+                        Fit(before, after, backw);
+                        DisplayBending();
+                        DisplayResult(before, after, backw);
+                    }
                     return kTRUE;
                 case kTbLoad:
-                    fBending.Load("bending_magic.txt");
+                    fBending.Load("bending_new.txt");
                     DisplayBending();
                     return kTRUE;
                 case kTbSave:
-                    fBending.Save("bending_magic.txt");
+                    fBending.Save("bending_new.txt");
                     return kTRUE;
                 case kTbLoadStars:
@@ -767 +818 @@
     }
 
+    void DisplayResult(Double_t before, Double_t after, Double_t backw)
+    {
+        TGLabel *l1 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+1);
+        l1->SetText(Form("Before: %.1f +- %.1f SE", before));
+
+        TGLabel *l2 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+2);
+        l2->SetText(Form("After:  %.1f +- %.1f SE", after));
+
+        TGLabel *l3 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+3);
+        l3->SetText(Form("Backw:  %.1f +- %.1f SE", backw));
+    }
+
 public:
     ~MFit()
@@ -866 +929 @@
         fLabel.Add(l);
 
+        l = new TGLabel(grp1, "");
+        l->SetTextJustify(kTextLeft);
+        grp1->AddFrame(l, hints5);
+        fList->Add(l);
+        fLabel.Add(l);
+
+        l = new TGLabel(grp1, "");
+        l->SetTextJustify(kTextLeft);
+        grp1->AddFrame(l, hints5);
+        fList->Add(l);
+        fLabel.Add(l);
+
+        l = new TGLabel(grp1, "");
+        l->SetTextJustify(kTextLeft);
+        grp1->AddFrame(l, hints5);
+        fList->Add(l);
+        fLabel.Add(l);
+
 
         ((TGCheckButton*)fList->FindWidget(0))->SetState(kButtonDown);
@@ -893 +974 @@
 void gui()
 {
+    gErrorIgnoreLevel = kError;
     new MFit;
 }