Index: trunk/MagicSoft/Cosy/tpoint/tpointfit.C
===================================================================
--- trunk/MagicSoft/Cosy/tpoint/tpointfit.C	(revision 1806)
+++ trunk/MagicSoft/Cosy/tpoint/tpointfit.C	(revision 1810)
@@ -1,4 +1,5 @@
 #include <fstream.h>
 #include <iostream.h>
+#include <iomanip.h>
 
 #include <TF1.h>
@@ -21,12 +22,4 @@
 // 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,16 +47,13 @@
     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,4 +69,5 @@
 
     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,12 +93,9 @@
         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,8 +104,8 @@
 {
     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,5 +128,5 @@
 
     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,6 +134,5 @@
         Set set = *(Set*)gCoordinates.At(i);
 
-        //set.Adjust(bend);
-        set.Adjust(bend, MyAdjust);
+        set.Adjust(bend);
 
         Double_t err = 1;
@@ -200,5 +187,5 @@
         return;
     }
-
+    /*
     if (r0<0)
     {
@@ -212,9 +199,10 @@
     }
 
-    /*
-     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,14 +256,32 @@
     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,19 +296,7 @@
     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,6 +307,14 @@
     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,65 +348,116 @@
     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,14 +467,25 @@
     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,25 +493,27 @@
     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());
 
 }