Changeset 1810 for trunk/MagicSoft

Timestamp:

03/11/03 13:52:48 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Cosy

Files:

• Changelog (modified) (1 diff)
• base/MStar.h (modified) (1 diff)
• base/MStarList.cc (modified) (3 diffs)
• base/MStarList.h (modified) (1 diff)
• base/timer.cc (modified) (1 diff)
• base/timer.h (modified) (1 diff)
• gui/MGAccuracy.cc (modified) (2 diffs)
• gui/MGPngReader.cc (modified) (12 diffs)
• gui/MGPngReader.h (modified) (3 diffs)
• main/MBending.cc (modified) (14 diffs)
• tpoint/tpointfit.C (modified) (16 diffs)
• videodev/CaosFilter.cc (modified) (2 diffs)
• videodev/CaosFilter.h (modified) (1 diff)

trunk/MagicSoft/Cosy/Changelog

@@ -1 +1 @@
                                                                   -*-*- END -*-*-
+ 2003/03/11 - Daniela Dorner, Thomas Bretz:
+ 
+   * base/MStar.h:
+     - added Compare
+     
+   * base/MStarList.cc:
+     - some small bugfixes
+   
+   * base/MStarList.h:
+     - added Sort
+     - added Expand
+     
+   * base/timer.[h,cc]:
+     - Added GetTimeval
+
+   * gui/MGAccuracy.cc:
+     - Exchanged zd/az in calculation of Residual!!!
+   
+   * gui/MGPngReader.[h,cc]:
+     - set default lim mag to 7.0
+     - added new ouput for the pointing position fPZdAz
+     - added/fixed TrackingError/CalcTrackingError
+     - changed Filter2 to CaosFilter
+     - reordered starguider stuff in Execute
+     - changed color of circles
+     
+   * main/MBending.[h,cc]:
+     - removed MAGIC1 and MAGIC2
+     - removed '-' from writing
+     - fixed some bugs in the enumerations of the coefficients
+     - added some formating option for output
+
+   * tpoint/tpointfit.C:
+     - removed usage of MyAdjust
+     - fixed the Calculation of the residuals
+     - fixed reading
+     - added some correction in case of an overflow (360deg/0deg)
+     - fixed drawing
+     - added second Migrad turn...
+     - changed the screen and graphical output
+     
+   * videodev/CaosFilter.[h,cc]:
+     - changed RemoveTwins to accept a radius
+
+
 
  2003/03/02 - Daniela Dorner, Thomas Bretz (LaPalma):

trunk/MagicSoft/Cosy/base/MStar.h

@@ -21 +21 @@
     void Set(Double_t mx, Double_t my) { fX=mx; fY=my; }
 
+    Int_t Compare(const TObject *obj) const
+    {
+        const MStar *const s = (MStar*)obj;
+
+        if (fMag<s->fMag)
+            return -1;
+
+        if (fMag>s->fMag)
+            return 1;
+
+        return 0;
+    }
+
+    Bool_t IsSortable() const { return kTRUE; }
+
+
     ClassDef(MStar, 1)
 };

trunk/MagicSoft/Cosy/base/MStarList.cc

@@ -1 +1 @@
 #include "MStarList.h"
+
+#include <iostream.h>
 
 void MStarList::RemoveTwins(Double_t radius)
@@ -18 +20 @@
             return;
 
+        fStars.RemoveAt(idx);
+
         MStarListIter Next(this, *first, radius);
-        Delete(idx);
 
         MStar *pos;
@@ -29 +32 @@
             mx += pos->GetX();
             my += pos->GetY();
-            Delete(pos);
+            fStars.Remove(pos);
             cnt++;
         }

trunk/MagicSoft/Cosy/base/MStarList.h

@@ -38 +38 @@
 
     void RemoveTwins(Double_t radius);
+
+    void Sort() { fStars.Sort(); }
+
+    void Expand(int n) { fStars.Expand(n); }
 };
 

trunk/MagicSoft/Cosy/base/timer.cc

@@ -26 +26 @@
 //    fDiv = fmod((fMs+fSecs)/(60*60*24), 1.0);
 }
+
+void Timer::GetTimeval(struct timeval *tv) const
+{
+    tv->tv_sec  = fSecs;
+    tv->tv_usec = fMs;
+}
+
 /*
 void Timer::Set(const long mjd)

trunk/MagicSoft/Cosy/base/timer.h

@@ -33 +33 @@
     void SetTimer(const struct timeval *tv);
 
+    void GetTimeval(struct timeval *tv) const;
+
     int GetSecs() { return fSecs; }
 

trunk/MagicSoft/Cosy/gui/MGAccuracy.cc

@@ -195 +195 @@
     aaz *= d2r;
 
+    const double el = TMath::Pi()/2-pzd;
+
     const double dphi2 = aaz/2.;
     const double cos2  = cos(dphi2)*cos(dphi2);
     const double sin2  = sin(dphi2)*sin(dphi2);
-    const double d     = cos(azd)*cos2 - cos(2*pzd+azd)*sin2;
+    const double d     = cos(azd)*cos2 - cos(2*el+azd)*sin2;
 
     double dist = acos(d);
@@ -250 +252 @@
 
     UpdateCross(x, y);
-    UpdateText(pos.Zd(), acc.Az(), acc.Zd());
+    UpdateText(pos.Zd(), acc.Zd(), acc.Az());
 
     SetModified();

trunk/MagicSoft/Cosy/gui/MGPngReader.cc

@@ -3 +3 @@
 #include <fstream.h>    // ifstream
 #include <iostream.h>   // cout
+#include <iomanip.h>   // cout
 
 #include <TTimer.h>
 
 #include <TGMenu.h>
+#include <TGLabel.h>
 #include <TSystem.h>
 #include <TGSplitter.h>    // TGHorizontal3DLine
@@ -13 +15 @@
 #include "MGImage.h"
 #include "MGCoordinates.h"
+#include "MGDispCoordinates.h"
 
 #include "coord.h"
@@ -255 +258 @@
     fList->Add(fLimMag);
 
-    fSao->SetLimitMag(8.0);
+    fSao->SetLimitMag(7.0);
 
     fInterpol = new MGPopupMenu(p);
@@ -314 +317 @@
     fList->Add(fCRaDec);
 
-    fCZdAz = new MGCoordinates(this, kETypeZdAz);
+    fCZdAz = new MGCoordinates(this, kETypeZdAz, kFALSE);
     fCZdAz->Move(240+12+10, fMenu->GetDefaultHeight()+584);
     AddFrame(fCZdAz);
     fList->Add(fCZdAz);
 
+    fPZdAz = new MGCoordinates(this, kETypeZdAz, kFALSE);
+    fPZdAz->Move(240+12+10, fMenu->GetDefaultHeight()+630);
+    AddFrame(fPZdAz);
+    fList->Add(fPZdAz);
+
+    TGLabel *l = new TGLabel(this, "Arb.-Sky Pos");
+    l->SetTextJustify(kTextLeft);
+    l->Move(480+32, fMenu->GetDefaultHeight()+590);
+    AddFrame(l);
+    fList->Add(l);
+
+    l = new TGLabel(this, "arcsec/pix");
+    l->SetTextJustify(kTextLeft);
+    l->Move(605, fMenu->GetDefaultHeight()+619);
+    AddFrame(l);
+    fList->Add(l);
+
+    l = new TGLabel(this, "Pointing Pos");
+    l->SetTextJustify(kTextLeft);
+    l->Move(480+32, fMenu->GetDefaultHeight()+655);
+    AddFrame(l);
+    fList->Add(l);
+
     const Double_t pixsize = 23.4;
 
@@ -328 +354 @@
     fPixSize = new TGTextEntry(this, txt, IDM_kPixSize);
     fPixSize->SetAlignment(kTextCenterX);
-    fPixSize->Move(600, fMenu->GetDefaultHeight()+584);
+    fPixSize->Move(547, fMenu->GetDefaultHeight()+617);
     AddFrame(fPixSize);
     fList->Add(fPixSize);
@@ -357 +383 @@
     MapSubwindows();
     MapWindow();
+
+    //------------------------------------------------------------
+    //    XY xy(3.819444, 24.05333);
+    //    fCRaDec->SetCoordinates(xy);
+    //    fRaDec->Set(xy.X()*360/24, xy.Y());
+    //------------------------------------------------------------
 }
 
 MGPngReader::MGPngReader(MObservatory::LocationName_t obs)
-: TGMainFrame(gClient->GetRoot(), 768, 700), fFile(NULL), fDx((768-kZOOM)/2), fDy((512-kZOOM)/2)
+: TGMainFrame(gClient->GetRoot(), 768, 740), fFile(NULL), fDx((768-kZOOM)/2), fDy((512-kZOOM)/2)
 {
     fSao = new StarCatalog(obs);
@@ -811 +843 @@
 }
 
-void MGPngReader::CalcTrackingError(MStarList &spots, MStarList &stars)
-{
+ZdAz MGPngReader::TrackingError(TArrayF &x, TArrayF &y, TArrayF &mag) const
+{
+    //
+    // Viewable area (FIXME: AZ)
+    //
+    TH2F h("Hist", "dX/dY",  77, -768/2-.5,  768/2+.5, 58, -576/2-.5,  576/2+.5); // 3
+
+    /*
+    TH1F hmag("HistMag", "Mag", 19, 0, 100);
+    for (int i=0; i<mag.GetSize(); i++)
+        hmag.Fill(mag[i]);
+        */
+
+    //
+    // Search for matching Magnitudes
+    //
+    for (int i=0; i<mag.GetSize(); i++)
+    {
+        if (mag[i]>48-15 && mag[i]<48+15)
+            h.Fill(x[i], y[i]);
+    }
+
+    //
+    // Serach for an excess in the histogram
+    //
+    Int_t mx, my, dummy;
+    h.GetMaximumBin(mx, my, dummy);
+
+    const double xmax = h.GetXaxis()->GetBinCenter(mx);
+    const double dx   = h.GetXaxis()->GetBinWidth(mx);
+
+    const double ymax = h.GetYaxis()->GetBinCenter(my);
+    const double dy   = h.GetYaxis()->GetBinWidth(my);
+
+    cout << setprecision(3);
+    cout << "Cut-XY:       " << xmax << " +- " << dx << " / " << ymax << " +- " << dy << endl;
+
+    TGraph g;
+    for (int i=0; i<mag.GetSize(); i++)
+    {
+        if (!(x[i]>xmax-dx && x[i]<xmax+dx &&
+              y[i]>ymax-dy && y[i]<ymax+dy &&
+              mag[i]>48-15 && mag[i]<48+15))
+            continue;
+
+        g.SetPoint(g.GetN(), x[i], y[i]);
+    }
+
+    cout << "Offset-XY:    " << g.GetMean(1) << " +- " << g.GetRMS(1) << " / ";
+    cout << g.GetMean(2) << " +- " << g.GetRMS(2) << endl;
+
+    AltAz pos0 = fSao->CalcAltAzFromPix(768/2,              576/2)*kRad2Deg;
+    AltAz pos1 = fSao->CalcAltAzFromPix(768/2+g.GetMean(1), 576/2+g.GetMean(2))*kRad2Deg;
+
+    pos1 -= pos0;
+
+    ofstream fout("tracking_error.txt");
+    fout << setprecision(10) << fSao->GetMjd()-52000 << " " << -pos1.Alt() << " " << pos1.Az() << endl;
+    fout.close();
+
+    return ZdAz(-pos1.Alt(), pos1.Az());
+}
+
+void MGPngReader::CalcTrackingError(Leds &leds, MStarList &stars)
+{
+    const Int_t max = leds.GetEntries();
+
     if (stars.GetRealEntries() < 3)
     {
@@ -819 +916 @@
     }
 
-    if (spots.GetRealEntries() < 1)
+    if (max < 1)
     {
         cout << "Sorry, less than 1 detected spot in FOV!" << endl;
@@ -825 +922 @@
     }
 
-    Int_t idx = 0;
-
-    MStarList sortedspots;
+    stars.Sort(); // Sort by magnitude
+
+    TString str = "data/tracking_";
+    str += fSao->GetMjd()-52000;
+    str += ".txt";
+
+    ofstream fout(str);
+
+    TArrayF x, y, mag;
+
+    Int_t num = 0;
+
+    // FIXME: Is predifined value 3 a good idea?
 
     MStar *star;
-    MStar *spot;
     MStarListIter NextStar(&stars);
-    MStarListIter NextSpot(&spots);
-
-    while ((spot=NextSpot()))
-    {
-        AltAz aa = fSao->CalcAltAzFromPix(spot->GetX(), spot->GetY());
-        spot->Set(aa.Az(), aa.Alt());
-    }
-
-    while ((star=NextStar()))
-    {
-        AltAz aa = fSao->CalcAltAzFromPix(star->GetX(), star->GetY());
-        star->Set(aa.Az(), aa.Alt());
-
-        const double aaz   = star->GetX();
-        const double dphi2 = aaz/2.;
-        const double cos2  = cos(dphi2)*cos(dphi2);
-        const double sin2  = sin(dphi2)*sin(dphi2);
-
-        Double_t min = 800;
-
-        NextSpot.Reset();
-        while ((spot=NextSpot()))
+    while ((star=NextStar()) && num++<max+3)
+    {
+        TIter NextSp(&leds);
+        Led *spot=NULL;
+        while ((spot=(Led*)NextSp()))
         {
-            const double pzd = TMath::Pi()/2-spot->GetY();
-            const double azd = TMath::Pi()/2-star->GetY();
-
-            const double d = cos(azd)*cos2 - cos(2*pzd+azd)*sin2;
-
-            const Double_t dist = acos(d);
-
-            if (dist>=min)
-                continue;
-
-            min = dist;
-            sortedspots.AddAt(idx, spot->GetX(), spot->GetY(), spot->GetMag());
+            const XY dpos(spot->GetX()-star->GetX(), spot->GetY()-star->GetY());
+
+            const Int_t idx = x.GetSize();
+
+            x.Set(idx+1);
+            y.Set(idx+1);
+            mag.Set(idx+1);
+
+            x.AddAt(dpos.X(), idx);
+            y.AddAt(dpos.Y(), idx);
+            mag.AddAt(spot->GetMag()/star->GetMag(), idx);
+
+            if (fout)
+                fout << x[idx] << " " << y[idx] << " " << mag[idx] << endl;
         }
-        if (min>768)
-        {
-            cout << "ERROR!!!!!!!!" << endl;
-            return;
-        }
-        idx++;
-    }
-
-    //
-    // Now we have in sortedspots the entries with the shortest distances
-    // to the corresponding ones in stars.
-    // Now calculate the tracking error.
-    //
-    NextStar.Reset();
-    MStarListIter NextSpot2(&sortedspots);
-
-    Double_t meanx=0;
-    Double_t meany=0;
-
-    while ((star=NextStar()))
-    {
-        spot = NextSpot2();
-
-        meanx += star->GetX() - spot->GetX();
-        meany += star->GetY() - spot->GetY();
-    }
-
-    meanx /= idx;
-    meany /= idx;
-
-    cout << "Tracking Error:  dAlt=" << meany*180/TMath::Pi();
-    cout << "°  dAz=" << meanx*180/TMath::Pi() << "°    (calculated";
-    cout << " with " << idx << " stars/spots)" << endl;
-}
-
+    }
+
+    ZdAz d = TrackingError(x, y, mag);
+
+    //
+    // Calculated offsets
+    //
+
+    // round= floor(x+.5)
+    cout << "Offset-ZdAz: " << d.Zd()*60 << "' / " << d.Az()*60 << "'" << endl;
+    cout << "Offset-ZdAz: " << d.Zd()/360*16384 << " / " << d.Az()/360*16384 << " (SE) " << endl;
+
+    //
+    // Current Pointing position
+    //
+    ZdAz cpos = fSao->GetZdAz()-d;
+    fPZdAz->SetCoordinates(cpos);
+}
 
 
@@ -945 +1019 @@
 
     MStarList spots;
+
+    /*
     if (fDisplay->IsEntryChecked(IDM_kStarguider))
         Filter2::Execute(spots, c);
     else
-        if (fDisplay->IsEntryChecked(IDM_kFilter))
-            Filter::Execute(c);
+     */
+    if (fDisplay->IsEntryChecked(IDM_kFilter))
+        Filter::Execute(c);
 
     if (fDisplay->IsEntryChecked(IDM_kCaosFilter))
@@ -1031 +1108 @@
     if (fDisplay->IsEntryChecked(IDM_kCatalog))
     {
+        Timer time(tm);
+
+        GetCoordinates();
+
+        MStarList stars;
+        fSao->GetStars(stars, time.GetMjd(), *fRaDec);
+
+        if (fDisplay->IsEntryChecked(IDM_kStarguider))
+        {
+            Leds leds;
+            CaosFilter::Execute(c, leds, 1);
+
+            cout << "Found: " << leds.GetEntries() << " leds" << endl;
+
+            CaosFilter::RemoveTwins(leds, 3);
+            leds.Compress();
+
+            cout << "Rest: " << leds.GetEntries() << " leds" << endl;
+
+            CalcTrackingError(leds, stars);
+        }
+
+        byte cimg[768*576];
+        fSao->GetImg(c, cimg, stars);
+
         DrawCircle(c, 0.5);
         DrawCircle(c, 1.0);
         DrawCircle(c, 1.5);
 
-        byte cimg[768*576];
-
-        GetCoordinates();
-
-        Timer time(tm);
-
-        MStarList stars;
-        fSao->GetStars(stars, time.GetMjd(), *fRaDec);
-        fSao->GetImg(c, cimg, stars);
-        //fSao->GetImg(c, cimg, time.CalcMjd(), *fRaDec);
+        fCZdAz->SetCoordinates(fSao->GetZdAz());
 
         fImage->DrawColImg(c, cimg);
-
-        fCZdAz->SetCoordinates(fSao->GetZdAz());
-
-        if (fDisplay->IsEntryChecked(IDM_kStarguider))
-            CalcTrackingError(spots, stars);
     }
     else
@@ -1067 +1155 @@
     {
         const int dy = (int)sqrt(rpix*rpix-dx*dx);
-        img[cx+dx + (cy-dy)*768] = 0xd0;
-        img[cx+dx + (cy+dy)*768] = 0xd0;
-        img[cx-dy + (cy+dx)*768] = 0xd0;
-        img[cx+dy + (cy+dx)*768] = 0xd0;
+        img[cx+dx + (cy-dy)*768] = 0x40;
+        img[cx+dx + (cy+dy)*768] = 0x40;
+        img[cx-dy + (cy+dx)*768] = 0x40;
+        img[cx+dy + (cy+dx)*768] = 0x40;
     }
 }

trunk/MagicSoft/Cosy/gui/MGPngReader.h

@@ -13 +13 @@
 #include "MGImage.h"
 
-#include <TH1.h>
-#include <TH2.h>
-#include <TGraph.h>
-#include <TCanvas.h>
+#include "coord.h"
 
-class AltAz;
-class RaDec;
+class TArrayF;
+class TH1F;
+class TH2F;
+class TGraph;
 
 class TTimer;
@@ -87 +86 @@
     MGCoordinates *fCZdAz;
 
+    MGCoordinates *fPZdAz;
+
     TGTextEntry   *fPixSize;
 
@@ -104 +105 @@
     void Toggle(MGPopupMenu *p, UInt_t id);
     void GetCoordinates();
-    void CalcTrackingError(MStarList &, MStarList &);
+    void CalcTrackingError(Leds &, MStarList &);
+    ZdAz TrackingError(TArrayF &alt, TArrayF &az, TArrayF &mag) const;
 
     void InitHists();

trunk/MagicSoft/Cosy/main/MBending.cc

@@ -107 +107 @@
         if (str=="ECEC")   fEcec = val;
         if (str=="ACEC")   fAcec = val;
-        if (str=="MAGIC1") fMagic1 = val;
-        if (str=="MAGIC2") fMagic2 = val;
+        //if (str=="MAGIC1") fMagic1 = val;
+        //if (str=="MAGIC2") fMagic2 = val;
 
         fin >> val;
@@ -143 +143 @@
     fout << "S   00   000000   000000  0000000" << endl;
     fout << setprecision(8);
-    fout << " IA     " << -kRad2Deg*fIa   << " -1" << endl;
-    fout << " IE     " << -kRad2Deg*fIe   << " -1" << endl;
-    fout << " CA     " << -kRad2Deg*fNpae << " -1" << endl;
-    fout << " NPAE   " << -kRad2Deg*fCa   << " -1" << endl;
-    fout << " AN     " << -kRad2Deg*fAn   << " -1" << endl;
-    fout << " AW     " << -kRad2Deg*fAw   << " -1" << endl;
-    fout << " NRX    " << -kRad2Deg*fNrx  << " -1" << endl;
-    fout << " NRY    " << -kRad2Deg*fNry  << " -1" << endl;
-    fout << " CRX    " << -kRad2Deg*fCrx  << " -1" << endl;
-    fout << " CRY    " << -kRad2Deg*fCry  << " -1" << endl;
-    fout << " ECES   " << -kRad2Deg*fEces << " -1" << endl;
-    fout << " ACES   " << -kRad2Deg*fAces << " -1" << endl;
-    fout << " ECEC   " << -kRad2Deg*fEcec << " -1" << endl;
-    fout << " ACEC   " << -kRad2Deg*fAcec << " -1" << endl;
-    fout << " MAGIC1 " << -kRad2Deg*fMagic1 << " -1" << endl;
-    fout << " MAGIC2 " << -kRad2Deg*fMagic2 << " -1" << endl;
+    fout << " IA     " << kRad2Deg*fIa   << " -1" << endl;
+    fout << " IE     " << kRad2Deg*fIe   << " -1" << endl;
+    fout << " CA     " << kRad2Deg*fNpae << " -1" << endl;
+    fout << " NPAE   " << kRad2Deg*fCa   << " -1" << endl;
+    fout << " AN     " << kRad2Deg*fAn   << " -1" << endl;
+    fout << " AW     " << kRad2Deg*fAw   << " -1" << endl;
+    fout << " NRX    " << kRad2Deg*fNrx  << " -1" << endl;
+    fout << " NRY    " << kRad2Deg*fNry  << " -1" << endl;
+    fout << " CRX    " << kRad2Deg*fCrx  << " -1" << endl;
+    fout << " CRY    " << kRad2Deg*fCry  << " -1" << endl;
+    fout << " ECES   " << kRad2Deg*fEces << " -1" << endl;
+    fout << " ACES   " << kRad2Deg*fAces << " -1" << endl;
+    fout << " ECEC   " << kRad2Deg*fEcec << " -1" << endl;
+    fout << " ACEC   " << kRad2Deg*fAcec << " -1" << endl;
+    //fout << " MAGIC1 " << -kRad2Deg*fMagic1 << " -1" << endl;
+    //fout << " MAGIC2 " << -kRad2Deg*fMagic2 << " -1" << endl;
     fout << "END" << endl;
 }
@@ -173 +173 @@
     p += CEC;
 
+    //    const AltAz MAGIC(0, -fMagic1*tan(p.Alt())-fMagic2);
+    //    p += MAGIC;
+
     const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()),  fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
     const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
@@ -183 +186 @@
     p += NRY;
 
-    const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
-    p += MAGIC;
+    //    const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
+    //    p += MAGIC;
 
     const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
@@ -191 +194 @@
     p += AN;
 
-//    const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
-//    p += MAGIC;
-
+    //    const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
+    //    p += MAGIC;
 
     const AltAz CA(0, -fCa/cos(p.Alt()));
@@ -227 +229 @@
     p -= AW;
 
-    const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
-    p -= MAGIC;
+    //    const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
+    //    p -= MAGIC;
 
     const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
@@ -285 +287 @@
         fEces =par[10]/kRad2Deg; // Elevation Centering Error (sin)
     case 10:
-        fCry  =par[9]/kRad2Deg;  // Alt/Az Coude Displacement (E-W)
+        fCry  =par[9] /kRad2Deg; // Alt/Az Coude Displacement (E-W)
     case 9:
-        fCrx  =par[8]/kRad2Deg;  // Alt/Az Coude Displacement (N-S)
+        fCrx  =par[8] /kRad2Deg; // Alt/Az Coude Displacement (N-S)
     case 8:
-        fNry  =par[7]/kRad2Deg;  // Nasmyth rotator displacement, vertical
+        fNry  =par[7] /kRad2Deg; // Nasmyth rotator displacement, vertical
     case 7:
-        fNrx  =par[6]/kRad2Deg;  // Nasmyth rotator displacement, horizontan
+        fNrx  =par[6] /kRad2Deg; // Nasmyth rotator displacement, horizontan
     case 6:
-        fAw   =par[5]/kRad2Deg;  // Azimuth Axis Misalignment (E-W)
+        fAw   =par[5] /kRad2Deg; // Azimuth Axis Misalignment (E-W)
     case 5:
-        fAn   =par[4]/kRad2Deg;  // Azimuth Axis Misalignment (N-S)
+        fAn   =par[4] /kRad2Deg; // Azimuth Axis Misalignment (N-S)
     case 4:
-        fCa   =par[3]/kRad2Deg;  // Left-Right Collimation Error
+        fCa   =par[3] /kRad2Deg; // Left-Right Collimation Error
     case 3:
-        fNpae =par[2]/kRad2Deg;  // Az-El Nonperpendicularity
+        fNpae =par[2] /kRad2Deg; // Az-El Nonperpendicularity
     case 2:
-        fIe   =par[1]/kRad2Deg;  // Index Error in Elevation
+        fIe   =par[1] /kRad2Deg; // Index Error in Elevation
     case 1:
-        fIa   =par[0]/kRad2Deg;  // Index Error in Azimuth
+        fIa   =par[0] /kRad2Deg; // Index Error in Azimuth
     }
 }
@@ -312 +314 @@
     {
     case 16:
-        par[15]=fMagic2*kRad2Deg; //
+        par[15]=kRad2Deg*fMagic2; //
     case 15:
-        par[14]=fMagic1*kRad2Deg; //
+        par[14]=kRad2Deg*fMagic1; //
     case 14:
-        par[13]=fAcec*kRad2Deg; // Azimuth Centering Error (cos)
+        par[13]=kRad2Deg*fAcec; // Azimuth Centering Error (cos)
     case 13:
-        par[12]=fEcec*kRad2Deg; // Elevation Centering Error (cos)
+        par[12]=kRad2Deg*fEcec; // Elevation Centering Error (cos)
     case 12:
-        par[11]=fAces*kRad2Deg; // Azimuth Centering Error (sin)
+        par[11]=kRad2Deg*fAces; // Azimuth Centering Error (sin)
     case 11:
-        par[10]=fEces*kRad2Deg; // Elevation Centering Error (sin)
+        par[10]=kRad2Deg*fEces; // Elevation Centering Error (sin)
     case 10:
-        par[9]=fCry*kRad2Deg; // Alt/Az Coude Displacement (E-W)
+        par[9] =kRad2Deg*fCry;  // Alt/Az Coude Displacement (E-W)
     case 9:
-        par[8]=fCrx*kRad2Deg; // Alt/Az Coude Displacement (N-S)
+        par[8] =kRad2Deg*fCrx;  // Alt/Az Coude Displacement (N-S)
     case 8:
-        par[7]=fNry*kRad2Deg; // Nasmyth rotator displacement, vertical
+        par[7] =kRad2Deg*fNry;  // Nasmyth rotator displacement, vertical
     case 7:
-        par[6]=fNrx*kRad2Deg; // Nasmyth rotator displacement, horizontan
+        par[6] =kRad2Deg*fNrx;  // Nasmyth rotator displacement, horizontan
     case 6:
-        par[5]=fAw*kRad2Deg; // Azimuth Axis Misalignment (E-W)
+        par[5] =kRad2Deg*fAw;   // Azimuth Axis Misalignment (E-W)
     case 5:
-        par[4]=fAn*kRad2Deg; // Azimuth Axis Misalignment (N-S)
+        par[4] =kRad2Deg*fAn;   // Azimuth Axis Misalignment (N-S)
     case 4:
-        par[3]=fCa*kRad2Deg; // Left-Right Collimation Error
+        par[3] =kRad2Deg*fCa;   // Left-Right Collimation Error
     case 3:
-        par[2]=fNpae*kRad2Deg; // Az-El Nonperpendicularity
+        par[2] =kRad2Deg*fNpae; // Az-El Nonperpendicularity
     case 2:
-        par[1]=fIe*kRad2Deg; // Index Error in Elevation
+        par[1] =kRad2Deg*fIe;   // Index Error in Elevation
     case 1:
-        par[0]=fIa*kRad2Deg; // Index Error in Azimuth
+        par[0] =kRad2Deg*fIa;   // Index Error in Azimuth
     }
 }
@@ -349 +351 @@
 {
     if (n<0)
-        n = m.GetNumPars();
+        n = 16;
 
     Int_t ierflg = 0;
@@ -357 +359 @@
     case 16:
         m.mnparm(15,"MAGIC2", fMagic2*kRad2Deg,  1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 15:
         m.mnparm(14,"MAGIC1", fMagic1*kRad2Deg,  1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 14:
         m.mnparm(13,"ACEC", fAcec*kRad2Deg,  1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 13:
         m.mnparm(12,"ECEC", fEcec*kRad2Deg,  1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 12:
         m.mnparm(11,"ACES", fAcec*kRad2Deg,  1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 11:
         m.mnparm(10,"ECES", fEcec*kRad2Deg,  1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 10:
         m.mnparm(9, "CRY",  fCry*kRad2Deg,  1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 9:
         m.mnparm(8, "CRX",  fCrx*kRad2Deg,  1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 8:
         m.mnparm(7, "NRY",  fNry*kRad2Deg,  1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 7:
         m.mnparm(6, "NRX",  fNrx*kRad2Deg,  1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 6:
         m.mnparm(5, "AW",   fAw*kRad2Deg,   1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 5:
         m.mnparm(4, "AN",   fAn*kRad2Deg,   1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 4:
         m.mnparm(3, "CA",   fCa*kRad2Deg,   1, -360, 360, ierflg);
-        // cout << "Init 3 CA:    " << fCa << endl;
     case 3:
         m.mnparm(2, "NPAE", fNpae*kRad2Deg, 1, -360, 360, ierflg);
-        // cout << "Init 2 NPAE:  " << fNpae << endl;
     case 2:
-        m.mnparm(1, "IE",   fIe*kRad2Deg,   1, -360, 360, ierflg);
-        // cout << "Init 1 IE:    " << fIe << endl;
+        m.mnparm(1, "IE",   fIe*kRad2Deg,   1,  -90,  90, ierflg);
     case 1:
         m.mnparm(0, "IA",   fIa*kRad2Deg,   1, -360, 360, ierflg);
-        // cout << "Init 0 IA:    " << fIa << endl;
     }
 }
@@ -408 +394 @@
 void MBending::GetMinuitParameters(TMinuit &m, Int_t n=-1)
 {
-    if (n<0)
+    if (n<0 || n>m.GetNumPars())
         n = m.GetNumPars();
 
@@ -465 +451 @@
     }
 }
-
+/*
+void FormatPar(TMinuit &m, Int_t n)
+{
+    Double_t par, err;
+    m.GetParameter(n, par, err);
+
+    int expp = (int)log10(par);
+    int expe = (int)log10(err);
+
+    if (err<2*pow(10, expe))
+        expe--;
+
+    Int_t exp = expe>expp ? expp : expe;
+
+    par = (int)(par/pow(10, exp)) * pow(10, exp);
+    err = (int)(err/pow(10, exp)) * pow(10, exp);
+
+    cout << par << " +- " << err << flush;
+}
+*/
 void MBending::PrintMinuitParameters(TMinuit &m, Int_t n=-1) const
 {
@@ -471 +476 @@
         n = m.GetNumPars();
 
+    cout << setprecision(3);
+
     Double_t par, err;
 
@@ -477 +484 @@
     case 16:
         m.GetParameter(15, par, err);
-        cout << " 15 MAGIC2: " << par << " +- " << err << endl;
+        cout << " 15 MAGIC2: " << setw(8) << par << " +- " << setw(4) <<  err << endl;
     case 15:
         m.GetParameter(14, par, err);
-        cout << " 14 MAGIC1: " << par << " +- " << err << endl;
+        cout << " 14 MAGIC1: " << setw(8) << par << " +- " << setw(4) <<  err << endl;
     case 14:
         m.GetParameter(13, par, err);
-        cout << " 13 ACEC: " << par << " +- " << err << endl;
+        cout << " 13 ACEC:   " << setw(8) << par << " +- " << setw(4) <<  err << "  Azimuth Centering Error (cos)" << endl;
     case 13:
         m.GetParameter(12, par, err);
-        cout << " 12 ECEC: " << par << " +- " << err << endl;
+        cout << " 12 ECEC:   " << setw(8) << par << " +- " << setw(4) <<  err << "  Elevation Centering Error (cos)" << endl;
     case 12:
         m.GetParameter(11, par, err);
-        cout << " 11 ACES: " << par << " +- " << err << endl;
+        cout << " 11 ACES:   " << setw(8) << par << " +- " << setw(4) <<  err << "  Azimuth Centering Error (sin)" << endl;
     case 11:
         m.GetParameter(10, par, err);
-        cout << " 10 ECES: " << par << " +- " << err << endl;
+        cout << " 10 ECES:   " << setw(8) << par << " +- " << setw(4) <<  err << "  Elevation Centering Error (sin)" << endl;
     case 10:
         m.GetParameter(9, par, err);
-        cout << "  9 CRY: " << par << " +- " << err << endl;
+        cout << "  9 CRY:    " << setw(8) << par << " +- " << setw(4) <<  err << "  Alt/Az Coude Displacement (E-W)" << endl;
     case 9:
         m.GetParameter(8, par, err);
-        cout << "  8 CRX: " << par << " +- " << err << endl;
+        cout << "  8 CRX:    " << setw(8) << par << " +- " << setw(4) <<  err << "  Alt/Az Coude Displacement (N-S)" << endl;
     case 8:
         m.GetParameter(7, par, err);
-        cout << "  7 NRY: " << par << " +- " << err << endl;
+        cout << "  7 NRY:    " << setw(8) << par << " +- " << setw(4) <<  err << "  Nasmyth rotator displacement, vertical" << endl;
     case 7:
         m.GetParameter(6, par, err);
-        cout << "  6 NRX: " << par << " +- " << err << endl;
+        cout << "  6 NRX:    " << setw(8) << par << " +- " << setw(4) <<  err << "  Nasmyth rotator displacement, horizontan" << endl;
     case 6:
         m.GetParameter(5, par, err);
-        cout << "  5 AW:  " << par << " +- " << err << endl;
+        cout << "  5 AW:     " << setw(8) << par << " +- " << setw(4) <<  err << "  Azimuth Axis Misalignment (E-W)" << endl;
     case 5:
         m.GetParameter(4, par, err);
-        cout << "  4 AN:  " << par << " +- " << err << endl;
+        cout << "  4 AN:     " << setw(8) << par << " +- " << setw(4) <<  err << "  Azimuth Axis Misalignment (N-S)" << endl;
     case 4:
         m.GetParameter(3, par, err);
-        cout << "  3 CA:  " << par << " +- " << err << endl;
+        cout << "  3 CA:     " << setw(8) << par << " +- " << setw(4) <<  err << "  Left-Right Collimation Error" << endl;
     case 3:
         m.GetParameter(2, par, err);
-        cout << "  2 NPAE:  " << par << " +- " << err << endl;
+        cout << "  2 NPAE:   " << setw(8) << par << " +- " << setw(4) <<  err << "  Az-El Nonperpendicularity" << endl;
     case 2:
         m.GetParameter(1, par, err);
-        cout << "  1 IE:  " << par << " +- " << err << endl;
+        cout << "  1 IE:     " << setw(8) << par << " +- " << setw(4) <<  err << "  Index Error Elevation (Offset)" << endl;
     case 1:
         m.GetParameter(0, par, err);
-        cout << "  0 IA:  " << par << " +- " << err << endl;
-    }
-}
+        cout << "  0 IA:     " << setw(8) << par << " +- " << setw(4) <<  err << "  Index Error Azimuth (Offset)" << endl;
+    }
+}

trunk/MagicSoft/Cosy/tpoint/tpointfit.C

@@ -1 +1 @@
 #include <fstream.h>
 #include <iostream.h>
+#include <iomanip.h>
 
 #include <TF1.h>
@@ -21 +22 @@
 // Sekans = 1/cos
 //
-
-void MyAdjust(AltAz &p, Double_t *par)
-{
-    // p[rad]
-    MBending bend;
-    bend.SetParameters(par); // [deg]
-    p=bend(p);
-}
 
 class Set : public TObject
@@ -54 +47 @@
     Double_t GetResidual() const
     {
-        Double_t daz = fRawEl-fStarEl;
-        Double_t del = fRawAz-fStarAz;
-
-        //Double_t rzd = TMath::Pi()/2-fRawEl;
-        Double_t rzd = TMath::Pi()/2-fOrigStarEl;
+        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(2*rzd+del)*sin2;
-
-        Double_t dist = acos(d);
+        Double_t d = cos(del)*cos2 - cos(2*fOrigStarEl+del)*sin2;
+
+        Double_t  dist = acos(d);
 
         return dist * 180 / TMath::Pi();
@@ -79 +69 @@
 
     Double_t GetDEl() const     { return (fRawEl-fStarEl)*180/TMath::Pi(); }
+    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(); }
@@ -102 +93 @@
         fStarAz = p.Az();
     }
-
-    void Adjust(const MBending &bend, void (*fcn)(AltAz &zdaz, Double_t *par))
-    {
-        AltAz star = GetStarAltAz();
-
-        AltAz p = bend(star, fcn);
-        fStarEl = p.Alt();
-        fStarAz = p.Az();
+    void AdjustBack(const MBending &bend)
+    {
+        AltAz p = bend.CorrectBack(GetRawAltAz());
+        fRawEl = p.Alt();
+        fRawAz = p.Az();
     }
 };
@@ -116 +104 @@
 {
     Double_t v[4];
+    fin >> v[0];
+    fin >> v[1];
     fin >> v[2];
     fin >> v[3];
-    fin >> v[0];
-    fin >> v[1];
 
     set.fStarAz = v[0]*TMath::Pi()/180;
@@ -140 +128 @@
 
     MBending bend;
-    bend.SetParameters(par/*, npar*/); // Set Parameters [deg] to MBending
+    bend.SetParameters(par); // Set Parameters [deg] to MBending
 
     for (int i=0; i<gCoordinates.GetSize(); i++)
@@ -146 +134 @@
         Set set = *(Set*)gCoordinates.At(i);
 
-        //set.Adjust(bend);
-        set.Adjust(bend, MyAdjust);
+        set.Adjust(bend);
 
         Double_t err = 1;
@@ -200 +187 @@
         return;
     }
-
+    /*
     if (r0<0)
     {
@@ -212 +199 @@
     }
 
-    /*
-     phi0 = fmod(phi0+360, 360);
-     phi1 = fmod(phi1+360, 360);
+    phi0 = fmod(phi0+360, 360);
+    phi1 = fmod(phi1+360, 360);
+
+    if (phi1-phi0<-180)
+        phi1+=360;
     */
-
     TLine line;
     line.SetLineWidth(2);
@@ -268 +256 @@
     Double_t phi1 = set.GetStarAz();
 
+    if (r0<0)
+    {
+        r0 = -r0;
+        phi0 += 180;
+    }
+    if (r1<0)
+    {
+        r1 = -r1;
+        phi1 += 180;
+    }
+
+    phi0 = fmod(phi0+360, 360);
+    phi1 = fmod(phi1+360, 360);
+
+    if (phi1-phi0<-180)
+        phi1+=360;
+
+    if (scale<0 || scale>1000)
+        scale = -1;
+
     if (scale>0)
     {
         Double_t d = r1-r0;
-        if (d<0) d=-d;
-        r0 -= scale*d;
-        r1 += scale*d;
+        r0 += scale*d;
+        r1 -= scale*d;
         d = phi1-phi0;
-        if (d<0) d=-d;
-        phi0 -= scale*d;
-        phi1 += scale*d;
+        phi0 += scale*d;
+        phi1 -= scale*d;
     }
 
@@ -290 +296 @@
     gCoordinates.SetOwner();
 
-    TH2F hdaz("dAz", "\\Delta Az vs. Alt",  32, 0,  90,  32, -1, 1);
-    TH2F hdzd("dZd", "\\Delta Alt vs. Az",  32, 0, 360,  32, -1, 1);
-
     TH1F hres1("Res1", "  Residuals before correction     ", 10, 0, 180);
     TH1F hres2("Res2", "  Residuals after correction     ",  10, 0,   1);
 
-    TH2F h2res1("Res2D1", "  Arb. Residuals before correction     ", 32, 0, 360,  10, 0, 90);
-    TH2F h2res2("Res2D2", "  Arb. Residuals after correction     ",  32, 0, 360,  10, 0, 90);
-
-    hdaz.SetXTitle("Zd [\\circ]");
-    hdaz.SetYTitle("\\Delta Az [\\circ]");
-
-    hdzd.SetXTitle("Az [\\circ]");
-    hdzd.SetYTitle("\\Delta Zd [\\circ]");
-
     hres1.SetXTitle("\\Delta [\\circ]");
     hres1.SetYTitle("Counts");
@@ -313 +307 @@
     TGraph gdaz;
     TGraph gdzd;
-
-    ifstream fin("tpoint/tpoint2.txt");
+    TGraph gaz;
+    TGraph gzd;
+
+    gdaz.SetTitle(" \\Delta Az vs. Zd ");
+    gdzd.SetTitle(" \\Delta Zd vs. Az ");
+
+    gaz.SetTitle(" \\Delta Az vs. Az ");
+    gzd.SetTitle(" \\Delta Zd vs. Zd ");
+
+    ifstream fin("tpoint/tpoint3.txt");
 
     while (fin && fin.get()!='\n');
@@ -346 +348 @@
     bending.SetMinuitParameters(minuit, 16); // Init Parameters [deg]
 
-    Int_t ierflg = 0;
-
-   //  minuit.FixParameter(2); // NPAE
-   //  minuit.FixParameter(3); // CA
-   //  minuit.FixParameter(4); // AN
-   //  minuit.FixParameter(5); // AW
-   // minuit.FixParameter(6);  // NRX
-   // minuit.FixParameter(7);  // NRY
+    //minuit.FixParameter(0);  //(1) IA
+    //minuit.FixParameter(1);  //(1) IE
+    minuit.FixParameter(2);  //(2) NPAE
+    //minuit.FixParameter(3);  //(1) CA
+    minuit.FixParameter(4);  //(2) AN
+    //minuit.FixParameter(5);  //(1) AW
+
+    minuit.FixParameter(6);  // NRX
+    minuit.FixParameter(7);  // NRY
     minuit.FixParameter(8);  // CRX
     minuit.FixParameter(9);  // CRY
-   //  minuit.FixParameter(10); // ECES
-   //  minuit.FixParameter(11); // ACES
-   //  minuit.FixParameter(12); // ECEC
-   //  minuit.FixParameter(13); // ACEC
-    // minuit.FixParameter(14); // MAGIC1
+    minuit.FixParameter(10); // ECES
+    minuit.FixParameter(11); //(2) ACES
+    minuit.FixParameter(12); // ECEC
+    minuit.FixParameter(13); //(2) ACEC
+    minuit.FixParameter(14); // MAGIC1
     minuit.FixParameter(15); // MAGIC2
 
+
     //minuit.Command("SHOW PARAMETERS");
+    //minuit.Command("SHOW LIMITS");
+    cout << endl;
+
+    Int_t ierflg = 0;
     ierflg = minuit.Migrad();
-    cout << endl << "Migrad returns " << ierflg << endl << endl;
-    minuit.Command("SHOW LIMITS");
-
-    cout << endl;
-
-    /*
-     minuit.FixParameter(0);
-     minuit.FixParameter(1);
-     minuit.Release(2);
-
-     ierflg = minuit.Migrad();
-     cout << endl << "Migrad returns " << ierflg << endl << endl;
-     */
-
+    cout << "Migrad returns " << ierflg << endl;
+    // minuit.Release(2);
+    ierflg = minuit.Migrad();
+    cout << "Migrad returns " << ierflg << endl << endl;
+
+    //
+    // Get Fit Results
+    //
     bending.GetMinuitParameters(minuit);
     bending.PrintMinuitParameters(minuit);
     bending.Save("bending_magic.txt");
 
+
+    //
+    // Make a copy of all list entries
+    //
     TList list;
     list.SetOwner();
-
     for (int i=0; i<gCoordinates.GetSize(); i++)
+        list.Add(new Set(*(Set*)gCoordinates.At(i)));
+
+    //
+    // Calculate correction and residuals
+    //
+    for (int i=0; i<gCoordinates.GetSize(); i++)
     {
         Set &set0 = *(Set*)gCoordinates.At(i);
 
-        list.Add(new Set(set0));
+        ZdAz za = set0.GetStarZdAz()*kRad2Deg;
 
         hres1.Fill(set0.GetResidual());
-        //set0.Adjust(bending);
-        set0.Adjust(bending, MyAdjust);
+        set0.Adjust(bending);
         hres2.Fill(set0.GetResidual());
 
-        hdzd.Fill(   set0.GetStarAz(), set0.GetDEl());
-        hdaz.Fill(90-set0.GetStarEl(), set0.GetDAz());
+        Double_t dz = fmod(set0.GetDAz()+720, 360);
+        if (dz>180)
+            dz -= 360;
 
         static int j=0;
-        gdzd.SetPoint(j,    set0.GetStarAz(), set0.GetDEl());
-        gdaz.SetPoint(j, 90-set0.GetStarEl(), set0.GetDAz());
+        gdzd.SetPoint(j, za.Az()/*   set0.GetStarAz()*/, set0.GetDZd());
+        gaz.SetPoint( j, za.Az()/*   set0.GetStarAz()*/, dz);
+        gdaz.SetPoint(j, za.Zd()/*90-set0.GetStarEl()*/, dz);
+        gzd.SetPoint( j, za.Zd()/*90-set0.GetStarEl()*/, set0.GetDZd());
         j++;
     }
 
+    //
+    // Check for overflows
+    //
+    const Stat_t ov = hres2.GetBinContent(hres2.GetNbinsX()+1);
+    if (ov>0)
+        cout << "WARNING: " << ov << " overflows in residuals." << endl;
+
+    //
+    // Print all data sets for which the backward correction is
+    // twice times worse than the residual gotten from the
+    // bending correction itself
+    //
     cout << endl;
+    cout << "Checking backward correction (raw-->star):" << endl;
+    for (int i=0; i<gCoordinates.GetSize(); i++)
+    {
+        Set set0(*(Set*)list.At(i));
+        Set set1(*(Set*)list.At(i));
+
+        set0.AdjustBack(bending);
+        set1.Adjust(bending);
+
+        Double_t res0 = set0.GetResidual();
+        Double_t res1 = set1.GetResidual();
+
+        Double_t diff = fabs(res0-res1);
+
+        if (diff<hres2.GetMean()*0.66)
+            continue;
+
+        cout << "DBack: " << setw(7) << set0.GetStarZd() << " " << setw(8) << set0.GetStarAz() << ":  ";
+        cout << "ResB="<< setw(7) << res0*60 << "  ResF=" << setw(7) << res1*60 << "  |ResB-ResF|=" << setw(7) << diff*60 << " arcmin" << endl;
+    }
+    cout << "OK." << endl;
+
+    //
+    // Print out the residual before and after correction in several
+    // units
+    //
+    cout << endl;
+    cout << gCoordinates.GetSize() << " data sets." << endl << endl;
     cout << "Spead before: " << hres1.GetMean() << "deg  +-  " << hres1.GetRMS() << "deg" << endl;
     cout << "Spead after:  " << hres2.GetMean() << "deg  +-  " << hres2.GetRMS() << "deg" << endl;
@@ -414 +467 @@
     cout << endl;
 
+
+
     TCanvas *c1=new TCanvas("c1", "c");
-    c1->Divide(3,2);
+    c1->Divide(2,2);
+
+    c1->cd(2);
+    gdaz.DrawClone("A*");
 
     c1->cd(1);
-    hdaz.DrawCopy("cont");
-    gdaz.DrawClone("*");
+    gaz.DrawClone("A*");
+
+    c1->cd(3);
+    gdzd.DrawClone("A*");
 
     c1->cd(4);
-    hdzd.DrawCopy("cont");
-    gdzd.DrawClone("*");
+    gzd.DrawClone("A*");
+
+
+
+    c1=new TCanvas("c2", "c", 800, 800);
+    c1->Divide(2,2);
 
     c1->cd(2);
@@ -429 +493 @@
     hres1.DrawCopy();
 
-    c1->cd(5);
+    c1->cd(4);
     hres2.SetLineColor(kBlue);
     hres2.DrawCopy();
 
-//     c1->cd(3);
-//     gPad->SetTheta(90);
-//     gPad->SetPhi(90);
-//     h2res1.DrawCopy("surf1pol");
-//     gPad->Modified();
-//     gPad->Update();
-//     for (int i=0; i<gCoordinates.GetSize(); i++)
-//         DrawSet(gPad, *(Set*)list.At(i), 10./hres1.GetMean());
-
-//     c1->cd(6);
-//     gPad->SetTheta(90);
-//     gPad->SetPhi(90);
-//     h2res1.DrawCopy("surf1pol");
-//     gPad->Modified();
-//     gPad->Update();
-//     for (int i=0; i<gCoordinates.GetSize(); i++)
-//         DrawSet(gPad, *(Set*)gCoordinates.At(i), 10./hres2.GetMean());
+    c1->cd(1);
+    gPad->SetTheta(90);
+    gPad->SetPhi(90);
+    TH2F h2res1("Res2D1", "   Arb. Residuals before correction (scaled)     ", 32, 0, 360,  10, 0, 90);
+    h2res1.DrawCopy("surf1pol");
+    gPad->Modified();
+    gPad->Update();
+    for (int i=0; i<gCoordinates.GetSize(); i++)
+        DrawSet(gPad, *(Set*)list.At(i), 10./hres1.GetMean());
+
+    c1->cd(3);
+    gPad->SetTheta(90);
+    gPad->SetPhi(90);
+    h2res1.SetTitle("   Arb. Residuals after correction (scaled)     ");
+    h2res1.DrawCopy("surf1pol");
+    gPad->Modified();
+    gPad->Update();
+    for (int i=0; i<gCoordinates.GetSize(); i++)
+        DrawSet(gPad, *(Set*)gCoordinates.At(i), 10./hres2.GetMean());
 
 }

trunk/MagicSoft/Cosy/videodev/CaosFilter.cc

@@ -345 +345 @@
     }
 
+    RemoveTwins(leds, 5);
+}
+
+void CaosFilter::RemoveTwins(Leds &leds, Double_t radius)
+{
     TIter Next1(&leds);
 
@@ -368 +373 @@
             const Double_t dy = y2-y1;
 
-            if (dx*dx+dy*dy<5*5)
+            if (dx*dx+dy*dy<radius*radius)
             {
                 // FIXME: Interpolation

trunk/MagicSoft/Cosy/videodev/CaosFilter.h

@@ -49 +49 @@
     static void Execute(byte *img, Leds &leds, Double_t conv);
 
-    static void  FilterLeds(Leds &leds);
+    static void FilterLeds(Leds &leds);
+    static void RemoveTwins(Leds &leds, Double_t radius);
 
     ClassDef(CaosFilter, 0)