Index: /trunk/MagicSoft/Cosy/Changelog
===================================================================
--- /trunk/MagicSoft/Cosy/Changelog	(revision 7621)
+++ /trunk/MagicSoft/Cosy/Changelog	(revision 7622)
@@ -1,3 +1,10 @@
                                                                   -*-*- END -*-*-
+ 2006/03/23 - Daniela Dorner, Thomas Bretz
+
+   * main/MStarguider.[h,cc]:
+     - changed starguider algorithm
+
+
+
  2006/03/19 - Daniela Dorner
 
@@ -28,4 +35,6 @@
      - removed starguider analyis (writing root-files) from starguider mode to 
        stabilize the code 
+     - added 'fGetter->ExitLoop()' before each 'delete fGetter'
+     - added SetDirectory(0) for histogram in starguider
 
 
Index: /trunk/MagicSoft/Cosy/main/MStarguider.cc
===================================================================
--- /trunk/MagicSoft/Cosy/main/MStarguider.cc	(revision 7621)
+++ /trunk/MagicSoft/Cosy/main/MStarguider.cc	(revision 7622)
@@ -643,4 +643,5 @@
         if (ch0!=ch1)
         {
+//            fGetter->ExitLoop();
             delete fGetter;
             usleep(150000); // FIX: Device or resource busy.
@@ -952,4 +953,5 @@
                     if (ch0!=ch1)
                     {
+//                        fGetter->ExitLoop();
                         delete fGetter;
                         usleep(150000); // FIX: Device or resource busy.
@@ -1159,4 +1161,5 @@
                     fChannel->CheckEntry  (ch1==0?IDM_kChannel1:IDM_kChannel2);
                     fChannel->UnCheckEntry(ch1==1?IDM_kChannel1:IDM_kChannel2);
+//                    fGetter->ExitLoop();
                     delete fGetter;
                     usleep(150000); // FIX: Device or resource busy.
@@ -1243,65 +1246,129 @@
 }
 
-ZdAz MStarguider::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
-    // chose a bit coarser binning to enhance excess 
-    // important: chose binning symmetrical around (0|0)!
-    //    TH2F h("Hist", "dX/dY", 49, -576/2-8, 576/2+8, 37, -576/2-8, 576/2+8); // 3
-    // reduced range in which histogram is filled, made bins a bit smaller
-    // 72pix equivalent to 0.6deg
-    TH2F h("Hist", "dX/dY", 7, -72, 72, 7, -72, 72);
-
-//     TH1F hmag("HistMag", "Mag", 19, 0, 100);
-//     for (int i=0; i<mag.GetSize(); i++)
-//         hmag.Fill(mag[i]);
-    
+ZdAz MStarguider::TrackingError(TArrayF &x, TArrayF &y, TArrayF &mag,
+                                Int_t &num) const
+{
+    num = -1;
+
+    // Width of the final 2D gaussian
+    const Double_t width = 0.8;
+
+    // The binning is made 1.7 sigma (which, in the 1D case, gives the 
+    // highest significance of a gaussian signal)  (1bin equiv 3x3 sigma)
+    const Double_t max =  50;
+    const Int_t    n   = TMath::Nint(2*max/(1.77*1.7*width));
+    //1.77=sqrt(pi) comes from pir=bin
+
+    // Fill a histogram with the dx/dy values from the file
+    TH2F hist("Hist", "dX/dY", n, -max, max, n, -max, max);
+//    hist.SetDirectory(0);
+
     //
     // 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]);
+    //}
+
+    // fill dx and dy into histogram
     for (int i=0; i<mag.GetSize(); i++)
-    {
-        if (mag[i]>48-15 && mag[i]<48+15)
-            h.Fill(x[i], y[i]);
-    }
-
-    //
-    // Search 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)/2;
-
-    const double ymax = h.GetYaxis()->GetBinCenter(my);
-    const double dy   = h.GetYaxis()->GetBinWidth(my)/2;
+        hist.Fill(x[i], y[i]);
+
+    // Remove all bins which have content of only a single event
+    // including under- and overflow bins
+    for (int i=0; i<hist.GetSize(); i++)
+        if (hist.GetBinContent(i)<1.5)
+            hist.SetBinContent(i, 0);
+
+    // Find the bin containing the maximum
+    Int_t bx, by, bz;
+    hist.GetMaximumBin(bx, by, bz);
+
+    // In the worst case the events are spread through the
+    // bins which are the neighbors of the maximum
+    // Normally neighbors which do not belong to the signal are empty!
+    hist.GetXaxis()->SetRange(bx-1, bx+1); // 3x3bins  ~8x8pix  ~9x9sigma
+    hist.GetYaxis()->SetRange(by-1, by+1); // 3x3bins  ~8x8pix  ~9x9sigma
+
+    // Check whether this region contains events
+    num = TMath::Nint(hist.Integral());
+    if (num<1)
+        return ZdAz();
+
+    // Get Mean for the given axis range
+    Double_t mx = hist.GetMean(1);
+    Double_t my = hist.GetMean(2);
+
+    const Int_t max1 = TMath::Nint(hist.GetBinContent(bx, by));
+
+    // Check if the maximum is unique
+    Int_t bx2, by2, bz2;
+    hist.GetXaxis()->SetRange(-1, 9999);
+    hist.GetYaxis()->SetRange(-1, 9999);
+    hist.SetBinContent(bx, by, 0);
+    hist.GetMaximumBin(bx2, by2, bz2);
+
+    const Int_t max2 = TMath::Nint(hist.GetBinContent(bx2, by2));
+    if (max1==max2 && TMath::Abs(bx2-bx)+TMath::Abs(by2-by)>1)
+        return ZdAz();
+
+    // loop again over the file and fill the histogram again.
+    // to get closer to the correct value
+    Double_t sumx = 0;
+    Double_t sumy = 0;
+             num  = 0;
+    for (int i=0; i<mag.GetSize(); i++)
+    {
+        // only fill values into the histogram which
+        // are inside 2*1.7=3.4 sigma (99.7%)
+        if (TMath::Hypot(x[i]-mx, y[i]-my)>width*3.4)
+            continue;
+
+        sumx += x[i];
+        sumy += y[i];
+        num++;
+    }
+
+    if (num<1)
+        return ZdAz();
+
+    // calc MeanX and MeanY
+    mx = sumx/num;
+    my = sumy/num;
+
+    // loop again to fill the mispointing corrected histograms
+    // and fill another histogram to check the quality of the calculated
+    // mispointing.
+    sumx = 0;
+    sumy = 0;
+    num  = 0;
+    for (int i=0; i<mag.GetSize(); i++)
+    {
+        // only fill values into the histogram which
+        // are inside 1.7=3.4 sigma (92%)
+        // Cut at 1.7 sigma which gives the best background supression
+        if (TMath::Hypot(x[i]-mx, y[i]-my)>width*1.7)
+            continue;
+
+        sumx += x[i];
+        sumy += y[i];
+        num++;
+    }
+
+    if (num<1)
+        return ZdAz(); // FIXME!!!!!!!!!!!!
+
+    // Mispointing
+    mx = sumx/num;
+    my = sumy/num;
 
 #ifdef EXPERT
-    cout << setprecision(3);
-    cout << "Cut-XY:       " << xmax << " +- " << dx << " / " << ymax << " +- " << dy << endl;
+    cout << "Offset-XY:    " << mx << " +- " << my << endl;
 #endif
 
-    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]);	
-    }
-  
-#ifdef EXPERT
-    cout << "Offset-XY:    " << g.GetMean(1) << " +- " << g.GetRMS(1) << " / ";
-    cout << g.GetMean(2) << " +- " << g.GetRMS(2) << endl;
-#endif
-
-    AltAz pos0 = fSao->CalcAltAzFromPix(768/2,              576/2)*kRad2Deg;
-    AltAz pos1 = fSao->CalcAltAzFromPix(768/2+g.GetMean(1), 576/2+g.GetMean(2))*kRad2Deg;
+    AltAz pos0 = fSao->CalcAltAzFromPix(768/2,    576/2)*kRad2Deg;
+    AltAz pos1 = fSao->CalcAltAzFromPix(768/2+mx, 576/2+my)*kRad2Deg;
 
     ofstream fout1("pointingpos.txt");
@@ -1319,27 +1386,9 @@
     fout2 << -pos1.Alt() << " " << pos1.Az() << endl;
 
-//    if (g.GetMean(1)>9 || g.GetMean(2)>9) {
-//      TFile f1("sguider-highoffset.root","UPDATE");
-//      h.Write();
-//      hmag.Write();
-//      g.Write();
-//      f1.Close(); 
-//    } else {
-//      TFile f1("sguider-niceoffset.root","UPDATE");
-//     h.Write();
-//      hmag.Write();
-//      g.Write();
- //     f1.Close();
-	
-
-//	}
-
-    //deleting histogram and graph
-    h.Delete();
-    g.Delete();
     return ZdAz(-pos1.Alt(), pos1.Az());
 }
 
-bool MStarguider::CalcTrackingError(Leds &leds, MStarList &stars, ZdAz &d, MTime &t, double &bright)
+Int_t MStarguider::CalcTrackingError(Leds &leds, MStarList &stars,
+                                     ZdAz &d, MTime &t, double &bright)
 {
     const Int_t max = leds.GetEntries();
@@ -1349,5 +1398,5 @@
         if (fStargTPoint->IsDown())
             fStargTPoint->SetDown(kFALSE);
-        return kFALSE;
+        return 0;
     }
     if (max < 3) //was 1
@@ -1356,5 +1405,5 @@
         if (fStargTPoint->IsDown())
             fStargTPoint->SetDown(kFALSE);
-        return kFALSE;
+        return 0;
     }
 
@@ -1383,5 +1432,6 @@
         while ((spot=(Led*)NextSp())) 
 	{
-            const XY dpos(spot->GetX()-(768-star->GetX()), spot->GetY()-star->GetY());
+            const XY dpos(spot->GetX()-(768-star->GetX()),
+                          spot->GetY()-star->GetY());
 	    
             const Int_t idx = x.GetSize();
@@ -1409,7 +1459,11 @@
     }
 
-    d = TrackingError(x, y, mag);
+    d = TrackingError(x, y, mag, num);
+
+    if (num<1)
+        return 0;
+
     fDZdAz->SetCoordinates(d);
-    
+
     //
     // Calculated offsets
@@ -1468,6 +1522,8 @@
     if ((t.GetMjd()-t2.GetMjd())>0.001) //1min20sec
     {
-        cout << "time difference between tpoint and starguider-tpoint > " <<
-            t.GetMjd()-t2.GetMjd() << "s => starguider tpoint hasn't been stored. Please repeat whole procedure. " << endl;
+        cout << "     time difference between tpoint and starguider-tpoint > 1 min *" <<
+            t.GetMjd()-t2.GetMjd() << "s) " << endl;
+        cout << "     => starguider tpoint hasn't been stored. " << endl;
+        cout << "     Please repeat whole procedure. " << endl;
         return kTRUE;
     }
@@ -1489,13 +1545,7 @@
     *fOutStargTp << endl;
 
-
-
-
-
-
     fTimeFromTp.Set(1970,1,1);
 
-    return kTRUE;
-
+    return num;
 }
 
@@ -1980,27 +2030,28 @@
             fSkyBright->SetBackgroundColor(color);
 
-  	    bool rc = CalcTrackingError(spots, stars, fD, t, bright);
-
-            if (rc && (bright <= 1.75* fLastBright) && (bright < 110))
+            const Int_t rc = CalcTrackingError(spots, stars, fD, t, bright);
+
+            if ((bright <= 1.75* fLastBright) && (bright < 110))
                 fStatus = MDriveCom::kMonitoring;
             else
                 fStatus = MDriveCom::kError;
 
-	    if (fCosy)
-		fPos = fCosy->GetPointingPos();   
-	 
-  	    if (fOperations->IsEntryChecked(IDM_kStargAnalysis)) 
+            if (fCosy)
+                fPos = fCosy->GetPointingPos();
+
+            if (fOperations->IsEntryChecked(IDM_kStargAnalysis))
                 fStargHistograms->Fill(spots, stars, fD,
                                        fSao->GetZdAz(), sgcenter, sgcenterzdaz,
                                        star, bright, fPos, t);
 
-	    fLastBright = bright;
+            fLastBright = bright;
 
             if (fCosy)
             {
                 MDriveCom &com = *fCosy->GetDriveCom();
-                com.SendStargReport(fStatus, fD, fSao->GetZdAz(), sgcenter, spots.GetEntries(), bright, time.GetMjd(), 0, 0);    // Report
+                com.SendStargReport(fStatus, fD, fSao->GetZdAz(),
+                                    sgcenter, rc, bright,
+                                    time.GetMjd(), 0, 0);    // Report
             }
-
         } //kStarguider
 
Index: /trunk/MagicSoft/Cosy/main/MStarguider.h
===================================================================
--- /trunk/MagicSoft/Cosy/main/MStarguider.h	(revision 7621)
+++ /trunk/MagicSoft/Cosy/main/MStarguider.h	(revision 7622)
@@ -129,7 +129,7 @@
     void ToggleCaosFilter();
     //void GetCoordinates();
-    bool CalcTrackingError(Leds &, MStarList &, ZdAz &, MTime &, double &bright);
+    Int_t CalcTrackingError(Leds &, MStarList &, ZdAz &, MTime &, double &bright);
     //void CalcTrackingError(Leds &, MStarList &);
-    ZdAz TrackingError(TArrayF &alt, TArrayF &az, TArrayF &mag) const;
+    ZdAz TrackingError(TArrayF &alt, TArrayF &az, TArrayF &mag, Int_t &num) const;
     bool Interpolate(const unsigned long n, byte *img) const;