Ignore:
Timestamp:
10/20/03 15:54:10 (21 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Cosy/tpoint/gui.C

    r1811 r2407  
     1#include <fstream.h>
    12#include <iostream.h>
     3#include <iomanip.h>
    24
    35#include <TGFrame.h>
    … …  
    57#include <TGButton.h>
    68
     9#include <TF1.h>
     10#include <TH1.h>
     11#include <TH2.h>
     12#include <TGraph.h>
     13
     14#include <TList.h>
     15#include <TStyle.h>
     16#include <TMinuit.h>
     17
     18#include <TView.h>
     19#include <TLine.h>
     20#include <TMarker.h>
     21#include <TCanvas.h>
     22
     23#include "coord.h"
     24
    725#include "MGList.h"
    8 
     26#include "MBending.h"
     27
     28class Set : public TObject
     29{
     30    friend ifstream &operator>>(ifstream &fin, Set &set);
     31private:
     32    Double_t fStarAz;
     33    Double_t fStarEl;
     34
     35    Double_t fRawAz;
     36    Double_t fRawEl;
     37
     38public:
     39    Set(Double_t sel=0, Double_t saz=0, Double_t rel=0, Double_t raz=0) :
     40        fStarAz(saz*TMath::Pi()/180),
     41        fStarEl(sel*TMath::Pi()/180),
     42        fRawAz(raz*TMath::Pi()/180),
     43        fRawEl(rel*TMath::Pi()/180)
     44    {
     45    }
     46
     47    Double_t GetResidual() const
     48    {
     49         Double_t del = fRawEl-fStarEl;
     50         Double_t daz = fRawAz-fStarAz;
     51         Double_t dphi2 = daz/2.;
     52         Double_t cos2  = cos(dphi2)*cos(dphi2);
     53         Double_t sin2  = sin(dphi2)*sin(dphi2);
     54         Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
     55
     56         Double_t dist = acos(d);
     57
     58         return dist * 180 / TMath::Pi();
     59    }
     60
     61    void operator=(Set &set)
     62    {
     63        fStarAz = set.fStarAz;
     64        fStarEl = set.fStarEl;
     65        fRawAz  = set.fRawAz;
     66        fRawEl  = set.fRawEl;
     67    }
     68
     69    Double_t GetDEl() const     { return (fRawEl-fStarEl)*180/TMath::Pi(); }
     70    Double_t GetDZd() const     { return -GetDEl(); }
     71    Double_t GetDAz() const     { return (fRawAz-fStarAz)*180/TMath::Pi(); }
     72    Double_t GetStarEl() const  { return fStarEl*180/TMath::Pi(); }
     73    Double_t GetStarZd() const  { return 90.-fStarEl*180/TMath::Pi(); }
     74    Double_t GetStarAz() const  { return fStarAz*180/TMath::Pi(); }
     75    Double_t GetRawEl() const   { return fRawEl*180/TMath::Pi(); }
     76    Double_t GetRawAz() const   { return fRawAz*180/TMath::Pi(); }
     77    Double_t GetRawZd() const   { return 90.-fRawEl*180/TMath::Pi(); }
     78
     79    ZdAz  GetStarZdAz() const   { return ZdAz(TMath::Pi()/2-fStarEl, fStarAz); }
     80    AltAz GetStarAltAz() const  { return AltAz(fStarEl, fStarAz); }
     81
     82    ZdAz  GetRawZdAz() const    { return ZdAz(TMath::Pi()/2-fRawEl, fRawAz); }
     83    AltAz GetRawAltAz() const   { return AltAz(fRawEl, fRawAz); }
     84
     85    void AdjustEl(Double_t del) { fStarEl += del*TMath::Pi()/180; }
     86    void AdjustAz(Double_t daz) { fStarAz += daz*TMath::Pi()/180; }
     87
     88    void Adjust(const MBending &bend)
     89    {
     90        AltAz p = bend(GetStarAltAz());
     91        fStarEl = p.Alt();
     92        fStarAz = p.Az();
     93    }
     94    void AdjustBack(const MBending &bend)
     95    {
     96        AltAz p = bend.CorrectBack(GetRawAltAz());
     97        fRawEl = p.Alt();
     98        fRawAz = p.Az();
     99    }
     100    ClassDef(Set, 0)
     101};
     102
     103ClassImp(Set);
     104
     105ifstream &operator>>(ifstream &fin, Set &set)
     106{
     107    Double_t v[4];
     108    fin >> v[0];
     109    fin >> v[1];
     110    fin >> v[2];
     111    fin >> v[3];
     112
     113    Double_t dummy;
     114    fin>>dummy;
     115    fin>>dummy;
     116    fin>>dummy;
     117
     118    set.fStarAz = v[0]*TMath::Pi()/180;
     119    set.fStarEl = v[1]*TMath::Pi()/180;
     120
     121    set.fRawAz  = v[2]*TMath::Pi()/180;
     122    set.fRawEl  = v[3]*TMath::Pi()/180;
     123
     124    return fin;
     125}
    9126
    10127class MFit : public TGMainFrame
    … …  
    12129private:
    13130    enum {
    14         kCbIA,
    15         kCbIE,
    16         kCbNPAE,
    17         kCbCA,
    18         kCbAN,
    19         kCbAW,
    20         kCbCRX,
    21         kCbCRY,
    22         kCbNRX,
    23         kCbNRY
     131        kTbFit = 19, //MBending::GetNumPar(), // FIXME!!!
     132        kTbLoad,
     133        kTbSave,
     134        kTbLoadStars,
     135        kTbReset,
     136        kTbResetStars
    24137    };
    25138
    26     MGList *list;
     139    MGList *fList;
     140
     141    TList fOriginal;
     142    TList fCoordinates;
     143    TList fLabel;
     144
     145    MBending fBending;
     146
     147    FontStruct_t fFont;
     148
     149    void Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
     150    {
     151        f = 0;
     152
     153        MBending bend;
     154        bend.SetParameters(par); // Set Parameters [deg] to MBending
     155
     156        for (int i=0; i<fCoordinates.GetSize(); i++)
     157        {
     158            Set set = *(Set*)fCoordinates.At(i);
     159
     160            set.Adjust(bend);
     161
     162            Double_t err = 0.02; // [deg] = 1SE
     163            Double_t res = set.GetResidual()/err;
     164
     165            f += res*res;
     166        }
     167
     168        //f /= (fCoordinates.GetSize()-7)*(fCoordinates.GetSize()-7);
     169        //f /= fCoordinates.GetSize()*fCoordinates.GetSize();
     170        f /= fCoordinates.GetSize();
     171    }
     172
     173    static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
     174    {
     175        ((MFit*)gMinuit->GetObjectFit())->Fcn(npar, gin, f, par, iflag);
     176    }
     177
     178    void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
     179    {
     180        TView *view = pad->GetView();
     181
     182        if (!view)
     183        {
     184            cout << "No View!" << endl;
     185            return;
     186        }
     187
     188        TMarker mark0;
     189        TMarker mark1;
     190        mark0.SetMarkerStyle(kStar);
     191        mark1.SetMarkerStyle(kStar);
     192        mark1.SetMarkerColor(kRed);
     193   
     194        r0 /= 90;
     195        r1 /= 90;
     196        phi0 *= TMath::Pi()/180;
     197        phi1 *= TMath::Pi()/180;
     198   
     199        Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
     200        Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
     201   
     202        Double_t y0[3], y1[3];
     203   
     204        view->WCtoNDC(x0, y0);
     205        view->WCtoNDC(x1, y1);
     206   
     207        mark0.DrawMarker(y0[0], y0[1]);
     208        mark1.DrawMarker(y1[0], y1[1]);
     209    }
     210   
     211    void DrawPolLine(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
     212    {
     213        TView *view = pad->GetView();
     214   
     215        if (!view)
     216        {
     217            cout << "No View!" << endl;
     218            return;
     219        }
     220        /*
     221        if (r0<0)
     222        {
     223            r0 = -r0;
     224            phi0 += 180;
     225        }
     226        if (r1<0)
     227        {
     228            r1 = -r1;
     229            phi1 += 180;
     230        }
     231   
     232        phi0 = fmod(phi0+360, 360);
     233        phi1 = fmod(phi1+360, 360);
     234   
     235        if (phi1-phi0<-180)
     236            phi1+=360;
     237        */
     238        TLine line;
     239        line.SetLineWidth(2);
     240        line.SetLineColor(kBlue);
     241   
     242        Double_t p0 = phi0<phi1?phi0:phi1;
     243        Double_t p1 = phi0<phi1?phi1:phi0;
     244   
     245        if (phi0>phi1)
     246        {
     247            Double_t d = r1;
     248            r1 = r0;
     249            r0 = d;
     250        }
     251   
     252        r0 /= 90;
     253        r1 /= 90;
     254   
     255        Double_t dr = r1-r0;
     256        Double_t dp = p1-p0;
     257   
     258        Double_t x0[3] = { r0*cos(p0*TMath::Pi()/180), r0*sin(p0*TMath::Pi()/180), 0};
     259   
     260        for (double i=p0+10; i<p1+10; i+=10)
     261        {
     262            if (i>p1)
     263                i=p1;
     264   
     265            Double_t r = dr/dp*(i-p0)+r0;
     266            Double_t p = TMath::Pi()/180*i;
     267   
     268            Double_t x1[3] = { r*cos(p), r*sin(p), 0};
     269   
     270            Double_t y0[3], y1[3];
     271   
     272            view->WCtoNDC(x0, y0);
     273            view->WCtoNDC(x1, y1);
     274   
     275            line.DrawLine(y0[0], y0[1], y1[0], y1[1]);
     276   
     277            x0[0] = x1[0];
     278            x0[1] = x1[1];
     279        }
     280    }
     281   
     282    void DrawSet(TVirtualPad *pad, Set &set, Float_t scale=-1)
     283    {
     284        Double_t r0   = set.GetRawZd();
     285        Double_t phi0 = set.GetRawAz();
     286        Double_t r1   = set.GetStarZd();
     287        Double_t phi1 = set.GetStarAz();
     288   
     289        if (r0<0)
     290        {
     291            r0 = -r0;
     292            phi0 += 180;
     293        }
     294        if (r1<0)
     295        {
     296            r1 = -r1;
     297            phi1 += 180;
     298        }
     299   
     300        phi0 = fmod(phi0+360, 360);
     301        phi1 = fmod(phi1+360, 360);
     302   
     303        if (phi1-phi0<-180)
     304            phi1+=360;
     305   
     306        if (scale<0 || scale>1000)
     307            scale = -1;
     308   
     309        if (scale>0)
     310        {
     311            Double_t d = r1-r0;
     312            r0 += scale*d;
     313            r1 -= scale*d;
     314            d = phi1-phi0;
     315            phi0 += scale*d;
     316            phi1 -= scale*d;
     317   
     318            DrawPolLine(pad, r0, phi0, r1, phi1);
     319            DrawMarker(pad,  r0, phi0, r1, phi1);
     320        }
     321        else
     322            DrawMarker(pad,  r1, phi1, 0 ,0);
     323    }
     324
     325    void Fit()
     326    {
     327        if (fOriginal.GetSize()==0)
     328        {
     329            cout << "Sorry, no input data loaded..." << endl;
     330            return;
     331        }
     332
     333        fCoordinates.Delete();
     334        for (int i=0; i<fOriginal.GetSize(); i++)
     335            fCoordinates.Add(new Set(*(Set*)fOriginal.At(i)));
     336
     337        cout << "-----------------------------------------------------------------------" << endl;
     338   
     339        TH1F hres1("Res1", " Residuals before correction ", fOriginal.GetSize()/2, 0, 3);
     340        TH1F hres2("Res2", " Residuals after correction ",  fOriginal.GetSize(), 0, 3);
     341        TH1F hres3("Res3", " Residuals after backward correction ",  fOriginal.GetSize(), 0, 3);
     342   
     343        hres1.SetXTitle("\\Delta [\\circ]");
     344        hres1.SetYTitle("Counts");
     345   
     346        hres2.SetXTitle("\\Delta [\\circ]");
     347        hres2.SetYTitle("Counts");
     348   
     349        hres3.SetXTitle("\\Delta [\\circ]");
     350        hres3.SetYTitle("Counts");
     351   
     352        TGraph gdaz;
     353        TGraph gdzd;
     354        TGraph gaz;
     355        TGraph gzd;
     356        TGraph graz;
     357        TGraph grzd;
     358   
     359        gdaz.SetTitle(" \\Delta Az vs. Zd ");
     360        gdzd.SetTitle(" \\Delta Zd vs. Az ");
     361   
     362        gaz.SetTitle(" \\Delta Az vs. Az ");
     363        gzd.SetTitle(" \\Delta Zd vs. Zd ");
     364
     365        graz.SetTitle(" \\Delta vs. Az ");
     366        grzd.SetTitle(" \\Delta vs. Zd ");
     367   
     368        TMinuit minuit(MBending::GetNumPar());  //initialize TMinuit with a maximum of 5 params
     369        minuit.SetObjectFit(this);
     370        minuit.SetPrintLevel(-1);
     371        minuit.SetFCN(fcn);
     372   
     373        fBending.SetMinuitParameters(minuit, MBending::GetNumPar()); // Init Parameters [deg]
     374   
     375        for (int i=0; i<MBending::GetNumPar(); i++)
     376        {
     377            TGButton *l = (TGButton*)fList->FindWidget(i);
     378            minuit.FixParameter(i);
     379            if (l->GetState()==kButtonDown)
     380                minuit.Release(i);
     381        }
     382   
     383        //minuit.Command("SHOW PARAMETERS");
     384        //minuit.Command("SHOW LIMITS");
     385   
     386        cout << endl;
     387        cout << "Starting fit..." << endl;
     388        cout << "For the fit an measurement error in the residual of ";
     389        cout << "0.02deg (=1SE) is assumed." << endl;
     390        cout << endl;
     391   
     392        Int_t ierflg = 0;
     393        ierflg = minuit.Migrad();
     394        cout << "Migrad returns " << ierflg << endl;
     395        // minuit.Release(2);
     396        ierflg = minuit.Migrad();
     397        cout << "Migrad returns " << ierflg << endl << endl;
     398   
     399        //
     400        // Get Fit Results
     401        //
     402        fBending.GetMinuitParameters(minuit);
     403        fBending.PrintMinuitParameters(minuit);
     404        //fBending.Save("bending_magic.txt");
     405   
     406   
     407        //
     408        // Make a copy of all list entries
     409        //
     410        TList list;
     411        list.SetOwner();
     412        for (int i=0; i<fCoordinates.GetSize(); i++)
     413            list.Add(new Set(*(Set*)fCoordinates.At(i)));
     414   
     415        //
     416        // Correct for Offsets only
     417        //
     418        TArrayD par;
     419        fBending.GetParameters(par);
     420        for (int i=2; i<MBending::GetNumPar(); i++)
     421            par[i]=0;
     422   
     423        MBending b2;
     424        b2.SetParameters(par);
     425   
     426        //
     427        // Calculate correction and residuals
     428        //
     429        for (int i=0; i<fCoordinates.GetSize(); i++)
     430        {
     431            Set &set0 = *(Set*)fCoordinates.At(i);
     432   
     433            ZdAz za = set0.GetStarZdAz()*kRad2Deg;
     434   
     435            //
     436            // Correct for offsets only
     437            //
     438            Set set1(set0);
     439            set1.Adjust(b2);
     440   
     441            hres1.Fill(set1.GetResidual());
     442   
     443            set0.Adjust(fBending);
     444            hres2.Fill(set0.GetResidual());
     445   
     446            Double_t dz = fmod(set0.GetDAz()+720, 360);
     447            if (dz>180)
     448                dz -= 360;
     449   
     450            gdzd.SetPoint(i, za.Az(), set0.GetDZd());
     451            gdaz.SetPoint(i, za.Zd(), dz);
     452            graz.SetPoint(i, za.Az(), set0.GetResidual());
     453            grzd.SetPoint(i, za.Zd(), set0.GetResidual());
     454            gaz.SetPoint( i, za.Az(), dz);
     455            gzd.SetPoint( i, za.Zd(), set0.GetDZd());
     456
     457        }
     458   
     459        //
     460        // Check for overflows
     461        //
     462        const Stat_t ov = hres2.GetBinContent(hres2.GetNbinsX()+1);
     463        if (ov>0)
     464            cout << "WARNING: " << ov << " overflows in residuals." << endl;
     465   
     466   
     467   
     468        cout << dec << endl;
     469        cout << "              Number of calls to FCN: " << minuit.fNfcn << endl;
     470        cout << "Minimum value found for FCN (Chi^2?): " << minuit.fAmin << endl;
     471        cout << "                  Fit-Probability(?): " << TMath::Prob(minuit.fAmin/*fOriginal.GetSize()*/, fOriginal.GetSize()-minuit.GetNumFreePars())*100 << "%" << endl;
     472        cout << "                           Chi^2/NDF: " << minuit.fAmin/(fOriginal.GetSize()-minuit.GetNumFreePars()) << endl;
     473        //cout << "Prob(?): " << TMath::Prob(fChisquare,ndf);
     474   
     475   
     476   
     477        //
     478        // Print all data sets for which the backward correction is
     479        // twice times worse than the residual gotten from the
     480        // bending correction itself
     481        //
     482        cout << endl;
     483        cout << "Checking backward correction (raw-->star):" << endl;
     484        for (int i=0; i<fCoordinates.GetSize(); i++)
     485        {
     486            Set set0(*(Set*)list.At(i));
     487            Set &set1 = *(Set*)list.At(i);
     488   
     489            set0.AdjustBack(fBending);
     490            set1.Adjust(fBending);
     491   
     492            const Double_t res0 = set0.GetResidual();
     493            const Double_t res1 = set1.GetResidual();
     494            const Double_t diff = TMath::Abs(res0-res1);
     495   
     496            hres3.Fill(res0);
     497   
     498            if (diff<hres2.GetMean()*0.66)
     499                continue;
     500   
     501            cout << "DBack: " << setw(6) << set0.GetStarZd() << " " << setw(7) << set0.GetStarAz() << ":  ";
     502            cout << "ResB="<< setw(7) << res0*60 << "  ResF=" << setw(7) << res1*60 << "  |ResB-ResF|=" << setw(7) << diff*60 << " arcmin" << endl;
     503        }
     504        cout << "OK." << endl;
     505        cout << endl;
     506   
     507
     508
     509
     510        TCanvas *c1;
     511
     512        if ((c1 = (TCanvas*)gROOT->FindObject("CanvGraphs")))
     513            delete c1;
     514        if ((c1 = (TCanvas*)gROOT->FindObject("CanvResiduals")))
     515            delete c1;
     516
     517
     518
     519        c1=new TCanvas("CanvGraphs", "Graphs");
     520        c1->Divide(2,3,0,0);
     521
     522        c1->cd(1);
     523        gPad->SetBorderMode(0);
     524        TGraph *g=(TGraph*)gaz.DrawClone("A*");
     525        g->SetBit(kCanDelete);
     526        g->GetHistogram()->SetXTitle("Az [\\circ]");
     527        g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]");
     528   
     529        c1->cd(2);
     530        gPad->SetBorderMode(0);
     531        g=(TGraph*)gdaz.DrawClone("A*");
     532        g->SetBit(kCanDelete);
     533        g->GetHistogram()->SetXTitle("Zd [\\circ]");
     534        g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]");
     535        cout << "Mean dAz: " << g->GetMean(2) << " \xb1 " << g->GetRMS(2) <<  endl;
     536   
     537        c1->cd(3);
     538        gPad->SetBorderMode(0);
     539        g=(TGraph*)gdzd.DrawClone("A*");
     540        g->SetBit(kCanDelete);
     541        g->GetHistogram()->SetXTitle("Az [\\circ]");
     542        g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]");
     543        cout << "Mean dZd: " << g->GetMean(2) << " \xb1 " << g->GetRMS(2) <<  endl;
     544        cout << endl;
     545   
     546        c1->cd(4);
     547        gPad->SetBorderMode(0);
     548        g=(TGraph*)gzd.DrawClone("A*");
     549        g->SetBit(kCanDelete);
     550        g->GetHistogram()->SetXTitle("Zd [\\circ]");
     551        g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]");
     552
     553        c1->cd(5);
     554        gPad->SetBorderMode(0);
     555        g=(TGraph*)graz.DrawClone("A*");
     556        g->SetBit(kCanDelete);
     557        g->GetHistogram()->SetXTitle("Az [\\circ]");
     558        g->GetHistogram()->SetYTitle("\\Delta [\\circ]");
     559   
     560        c1->cd(6);
     561        gPad->SetBorderMode(0);
     562        g=(TGraph*)grzd.DrawClone("A*");
     563        g->SetBit(kCanDelete);
     564        g->GetHistogram()->SetXTitle("Zd [\\circ]");
     565        g->GetHistogram()->SetYTitle("\\Delta [\\circ]");
     566
     567
     568
     569        //
     570        // Print out the residual before and after correction in several
     571        // units
     572        //
     573        cout << fCoordinates.GetSize() << " data sets." << endl << endl;
     574        cout << "Total Spread of Residual:" << endl;
     575        cout << "-------------------------" << endl;
     576        cout << "before: " << hres1.GetMean() << " \xb1 " << hres1.GetRMS() << " deg " << endl;
     577        cout << "after:  " << hres2.GetMean() << " \xb1 " << hres2.GetRMS() << " deg " << endl;
     578        cout << "before: " << (int)(hres1.GetMean()*60+.5) << " \xb1 " << (int)(hres1.GetRMS()*60+.5) << " arcmin" << endl;
     579        cout << "after:  " << (int)(hres2.GetMean()*60+.5) << " \xb1 " << (int)(hres2.GetRMS()*60+.5) << " arcmin" << endl;
     580        cout << "before: " << (int)(hres1.GetMean()*60*60/23.4+.5) << " \xb1 " << (int)(hres1.GetRMS()*60*60/23.4+.5) << " pix" << endl;
     581        cout << "after:  " << (int)(hres2.GetMean()*60*60/23.4+.5) << " \xb1 " << (int)(hres2.GetRMS()*60*60/23.4+.5) << " pix" << endl;
     582        cout << "after:  " << (int)(hres2.GetMean()*16384/360+.5) << " \xb1 " << (int)(hres2.GetRMS()*16384/360+.5) << " SE" << endl;
     583        cout << endl;
     584        cout << "backw:  " << (int)(hres3.GetMean()*60+.5) << " \xb1 " << (int)(hres3.GetRMS()*60+.5) << " arcmin" << endl;
     585        cout << endl;                                            // ±
     586   
     587   
     588   
     589        gStyle->SetOptStat(1110);
     590        gStyle->SetStatFormat("6.2g");
     591
     592
     593        c1=new TCanvas("CanvResiduals", "Residuals", 800, 800);
     594        c1->Divide(2, 2, 0, 0);
     595   
     596        c1->cd(2);
     597        gPad->SetBorderMode(0);
     598        hres1.SetLineColor(kRed);
     599        hres1.DrawCopy();
     600   
     601        c1->cd(4);
     602        gPad->SetBorderMode(0);
     603        hres2.SetLineColor(kBlue);
     604        TH1 *h=hres2.DrawCopy();
     605        TF1 f("mygaus", "(gaus)", 0, 1);
     606        f.SetLineColor(kMagenta/*6*/);
     607        f.SetLineWidth(1);
     608        f.SetParameter(0, h->GetBinContent(1));
     609        f.FixParameter(1, 0);
     610        f.SetParameter(2, h->GetRMS());
     611        h->Fit("mygaus", "QR");
     612        hres3.SetLineColor(kGreen);
     613        hres3.SetLineStyle(kDashed);
     614        hres3.DrawCopy("same");
     615        cout << "Gaus-Fit  Sigma: " << f.GetParameter(2) << "\xb0" << endl;
     616        cout << "Fit-Probability: " << f.GetProb()*100 << "%" << endl;
     617        cout << "      Chi^2/NDF: " << f.GetChisquare()/f.GetNDF() << endl;
     618
     619        c1->cd(1);
     620        gPad->SetBorderMode(0);
     621        gPad->SetTheta(90);
     622        gPad->SetPhi(90);
     623        TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360,  10, 0, 90);
     624        h2res1.SetBit(TH1::kNoStats);
     625        h2res1.DrawCopy("surf1pol");
     626        gPad->Modified();
     627        gPad->Update();
     628        for (int i=0; i<fCoordinates.GetSize(); i++)
     629            DrawSet(gPad, *(Set*)list.At(i));//, 10./hres1.GetMean());
     630   
     631        c1->cd(3);
     632        gPad->SetBorderMode(0);
     633        gPad->SetTheta(90);
     634        gPad->SetPhi(90);
     635        h2res1.SetTitle(" Arb. Residuals after correction (scaled) ");
     636        h2res1.DrawCopy("surf1pol");
     637        gPad->Modified();
     638        gPad->Update();
     639        for (int i=0; i<fCoordinates.GetSize(); i++)
     640            DrawSet(gPad, *(Set*)fCoordinates.At(i), 10./hres2.GetMean());
     641
     642        RaiseWindow();
     643    }
     644
     645    void LoadStars()
     646    {
     647        const Int_t size = fOriginal.GetSize();
     648
     649        ifstream fin("../tpoint_new.txt");
     650
     651        while (fin && fin.get()!='\n');
     652        while (fin && fin.get()!='\n');
     653        while (fin && fin.get()!='\n');
     654        if (!fin)
     655        {
     656            cout << "File not found!" << endl;
     657            return;
     658        }
     659
     660        while (1)
     661        {
     662            Set set;
     663
     664            fin >> set;  // Read data from file [deg], it is stored in [rad]
     665
     666            if (!fin)
     667                break;
     668
     669            fOriginal.Add(new Set(set));
     670        }
     671
     672        cout << "Found " << fOriginal.GetSize()-size << " sets of coordinates ";
     673        cout << "(Total=" << fOriginal.GetSize() << ")" << endl;
     674    }
    27675
    28676    Bool_t ProcessMessage(Long_t msg, Long_t mp1, Long_t mp2)
    29677    {
    30         cout << "Msg: " << hex << GET_MSG(msg) << endl;
    31         cout << "SubMsg: " << hex << GET_SUBMSG(msg) << endl;
     678        // cout << "Msg: " << hex << GET_MSG(msg) << endl;
     679        // cout << "SubMsg: " << hex << GET_SUBMSG(msg) << dec << endl;
    32680        switch (GET_MSG(msg))
    33681        {
    … …  
    35683            switch (GET_SUBMSG(msg))
    36684            {
    37             case kCM_CHECKBUTTON:
    38                 /*
    39                  TGCheckButton *but = (TGCheckButton*)list->FindWidget(mp1);
    40                  cout << hex << but->GetName() << " " << mp2 << endl;
    41                  */
     685            case kCM_BUTTON:
     686                switch (mp1)
     687                {
     688                case kTbFit:
     689                    Fit();
     690                    DisplayBending();
     691                    return kTRUE;
     692                case kTbLoad:
     693                    fBending.Load("bending_magic.txt");
     694                    DisplayBending();
     695                    return kTRUE;
     696                case kTbSave:
     697                    fBending.Save("bending_magic.txt");
     698                    return kTRUE;
     699                case kTbLoadStars:
     700                    LoadStars();
     701                    DisplayData();
     702                    return kTRUE;
     703                case kTbReset:
     704                    fBending.Clear();
     705                    DisplayBending();
     706                    return kTRUE;
     707                case kTbResetStars:
     708                    fOriginal.Delete();
     709                    DisplayData();
     710                    return kTRUE;
     711                }
    42712                return kTRUE;
    43713            }
    … …  
    52722        but->Associate(this);
    53723        f->AddFrame(but, h);
    54         list->Add(but);
     724        fList->Add(but);
    55725
    56726    }
    … …  
    61731        but->Associate(this);
    62732        f->AddFrame(but, h);
    63         list->Add(but);
    64     }
    65 
    66     void AddLabel(TGCompositeFrame *f, TString txt, TGLayoutHints *h=0, TList *label=0)
    67     {
    68         TGLabel *l = new TGLabel(f, txt);
     733        fList->Add(but);
     734    }
     735
     736    TGLabel *AddLabel(TGCompositeFrame *f, TString txt, TGLayoutHints *h=0)
     737    {
     738        TGLabel *l = new TGLabel(f, txt/*, TGLabel::GetDefaultGC()(), fFont*/);
    69739        f->AddFrame(l, h);
    70         list->Add(l);
    71         if (label)
    72             label->Add(l);
     740        fList->Add(l);
     741        fLabel.Add(l);
     742        return l;
     743    }
     744
     745    void DisplayBending()
     746    {
     747        TArrayD par, err;
     748        fBending.GetParameters(par);
     749        fBending.GetError(err);
     750
     751        TGLabel *l;
     752
     753        for (int i=0; i<MBending::GetNumPar(); i++)
     754        {
     755            l = (TGLabel*)fLabel.At(i);
     756            l->SetText(Form("%.4f\xb0", par[i]));
     757
     758            l = (TGLabel*)fLabel.At(MBending::GetNumPar()+i);
     759            l->SetText(Form("\xb1 %8.4f\xb0", err[i]));
     760        }
     761    }
     762
     763    void DisplayData()
     764    {
     765        TGLabel *l = (TGLabel*)fLabel.At(3*MBending::GetNumPar());
     766        l->SetText(Form("%d data sets loaded.", fOriginal.GetSize()));
    73767    }
    74768
    75769public:
    76     MFit() : TGMainFrame(gClient->GetRoot(), 550, 350, kHorizontalFrame)
    77     {
    78         list = new MGList;
    79         list->SetOwner();
     770    ~MFit()
     771    {
     772        if (fFont)
     773            gVirtualX->DeleteFont(fFont);
     774    }
     775    MFit() : TGMainFrame(gClient->GetRoot(), 740, 370, kHorizontalFrame)
     776    {
     777        fCoordinates.SetOwner();
     778        fOriginal.SetOwner();
     779
     780        fList = new MGList;
     781        fList->SetOwner();
     782
     783        fFont = gVirtualX->LoadQueryFont("7x13bold");
    80784
    81785        TGLayoutHints *hints0 = new TGLayoutHints(kLHintsExpandY, 7, 5, 5, 6);
    82786        TGLayoutHints *hints1 = new TGLayoutHints(kLHintsExpandX|kLHintsExpandY, 5, 7, 5, 6);
    83787        TGLayoutHints *hints2 = new TGLayoutHints(kLHintsCenterX|kLHintsCenterY, 5, 5, 5, 5);
    84         list->Add(hints0);
    85         list->Add(hints1);
    86         list->Add(hints2);
     788        fList->Add(hints0);
     789        fList->Add(hints1);
     790        fList->Add(hints2);
    87791
    88792        TGGroupFrame *grp1 = new TGGroupFrame(this, "Control", kVerticalFrame);
    89793        AddFrame(grp1, hints0);
    90         list->Add(grp1);
     794        fList->Add(grp1);
    91795
    92796        TGGroupFrame *grp2 = new TGGroupFrame(this, "Parameters", kHorizontalFrame);
    93797        AddFrame(grp2, hints1);
    94         list->Add(grp2);
    95 
    96 
    97 
    98         TGLayoutHints *hints4 = new TGLayoutHints(kLHintsExpandX, 5, 5, 10);
    99         AddTextButton(grp1, "Load Bending Model", -1, hints4);
    100         AddTextButton(grp1, "Save Bending Model", -1, hints4);
    101         AddTextButton(grp1, "Load Observations",  -1, hints4);
    102         AddTextButton(grp1, "Fit Parameters",     -1, hints4);
    103         AddTextButton(grp1, "Reset Parameters",   -1, hints4);
     798        fList->Add(grp2);
     799
     800
     801
     802        TGLayoutHints *hints4 = new TGLayoutHints(kLHintsExpandX, 5, 5,  5);
     803        TGLayoutHints *hints5 = new TGLayoutHints(kLHintsExpandX, 5, 5, 15);
     804        AddTextButton(grp1, "Load Pointing Model", kTbLoad,       hints5);
     805        AddTextButton(grp1, "Save Pointing Model", kTbSave,       hints4);
     806        AddTextButton(grp1, "Fit Parameters",      kTbFit,        hints5);
     807        AddTextButton(grp1, "Reset Parameters",    kTbReset,      hints4);
     808        AddTextButton(grp1, "Load Stars",          kTbLoadStars,  hints5);
     809        AddTextButton(grp1, "Reset Stars",         kTbResetStars, hints4);
     810        fList->Add(hints4);
     811        fList->Add(hints5);
    104812
    105813
    106814        TGHorizontalFrame *comp = new TGHorizontalFrame(grp2, 1, 1);
    107815        grp2->AddFrame(comp);
    108         list->Add(comp);
     816        fList->Add(comp);
    109817
    110818        TGLayoutHints *hints3 = new TGLayoutHints(kLHintsLeft|kLHintsTop, 0, 20, 5, 0);
    111         list->Add(hints3);
     819        fList->Add(hints3);
    112820
    113821        TGLabel *l;
    … …  
    115823        TGVerticalFrame *vframe = new TGVerticalFrame(comp, 1, 1);
    116824        comp->AddFrame(vframe, hints3);
    117         list->Add(vframe);
    118 
    119         AddCheckButton(vframe, "IA",   kCbIA);
    120         AddCheckButton(vframe, "IE",   kCbIE);
    121         AddCheckButton(vframe, "NPAE", kCbNPAE);
    122         AddCheckButton(vframe, "CA",   kCbCA);
    123         AddCheckButton(vframe, "AN",   kCbAN);
    124         AddCheckButton(vframe, "AW",   kCbAW);
    125         AddCheckButton(vframe, "CRX",  kCbCRX);
    126         AddCheckButton(vframe, "CRY",  kCbCRY);
    127         AddCheckButton(vframe, "NRX",  kCbNRX);
    128         AddCheckButton(vframe, "NRY",  kCbNRY);
     825        fList->Add(vframe);
     826
     827        for (int i=0; i<MBending::GetNumPar(); i++)
     828            AddCheckButton(vframe, fBending.GetName(i), i);
    129829
    130830        vframe = new TGVerticalFrame(comp, 1, 1);
    131831        comp->AddFrame(vframe, hints3);
    132         list->Add(vframe);
    133 
    134 
    135         TList label;
    136 
    137         l = new TGLabel(vframe, "0");
    138         list->Add(l);
    139         label.Add(l);
    140 
    141         TGButton *but = (TGButton*)list->FindWidget(kCbIA);
     832        fList->Add(vframe);
     833
     834        l = new TGLabel(vframe, "+000.0000");
     835        l->SetTextJustify(kTextRight);
     836        fList->Add(l);
     837        fLabel.Add(l);
     838
     839        TGButton *but = (TGButton*)fList->FindWidget(0);
    142840
    143841        TGLayoutHints *h = new TGLayoutHints(kLHintsCenterY, 0, 0, but->GetHeight()-l->GetHeight());
    144         list->Add(h);
     842        fList->Add(h);
    145843
    146844        vframe->AddFrame(l,h);
    147845
    148         AddLabel(vframe, "0", h, &label);
    149         AddLabel(vframe, "0", h, &label);
    150         AddLabel(vframe, "0", h, &label);
    151         AddLabel(vframe, "0", h, &label);
    152         AddLabel(vframe, "0", h, &label);
    153         AddLabel(vframe, "0", h, &label);
    154         AddLabel(vframe, "0", h, &label);
    155         AddLabel(vframe, "0", h, &label);
    156         AddLabel(vframe, "0", h, &label);
     846        for (int i=1; i<MBending::GetNumPar(); i++)
     847            AddLabel(vframe, "+000.0000", h)->SetTextJustify(kTextRight);
    157848
    158849        vframe = new TGVerticalFrame(comp, 1, 1);
    159850        comp->AddFrame(vframe, hints3);
    160         list->Add(vframe);
    161 
    162         AddLabel(vframe, "Offset Azimuth",         h, &label);
    163         AddLabel(vframe, "Offset Zenith Distance", h, &label);
    164         AddLabel(vframe, "0", h, &label);
    165         AddLabel(vframe, "0", h, &label);
    166         AddLabel(vframe, "0", h, &label);
    167         AddLabel(vframe, "0", h, &label);
    168         AddLabel(vframe, "0", h, &label);
    169         AddLabel(vframe, "0", h, &label);
    170         AddLabel(vframe, "0", h, &label);
     851        fList->Add(vframe);
     852
     853        for (int i=0; i<MBending::GetNumPar(); i++)
     854            AddLabel(vframe, "\xb1 00.0000\xb0", h)->SetTextJustify(kTextRight);
     855
     856        vframe = new TGVerticalFrame(comp, 1, 1);
     857        comp->AddFrame(vframe, hints3);
     858        fList->Add(vframe);
     859
     860        for (int i=0; i<MBending::GetNumPar(); i++)
     861            AddLabel(vframe, fBending.GetDescription(i), h);
     862
     863        l = new TGLabel(grp1, "0000000 Data Sets loaded.");
     864        grp1->AddFrame(l, hints5);
     865        fList->Add(l);
     866        fLabel.Add(l);
     867
     868
     869        ((TGCheckButton*)fList->FindWidget(0))->SetState(kButtonDown);
     870        ((TGCheckButton*)fList->FindWidget(1))->SetState(kButtonDown);
     871        ((TGCheckButton*)fList->FindWidget(3))->SetState(kButtonDown);
     872        ((TGCheckButton*)fList->FindWidget(4))->SetState(kButtonDown);
     873        ((TGCheckButton*)fList->FindWidget(5))->SetState(kButtonDown);
     874        ((TGCheckButton*)fList->FindWidget(6))->SetState(kButtonDown);
     875        ((TGCheckButton*)fList->FindWidget(7))->SetState(kButtonDown);
     876
     877        SetWindowName("TPoint Fitting Window");
     878        SetIconName("TPoint++");
    171879
    172880        Layout();
    … …  
    175883        MapWindow();
    176884
    177 
    178 
    179 
    180         l = (TGLabel*)label.At(0);
    181         l->SetText("221.4567");
    182 
    183         l = (TGLabel*)label.At(1);
    184         l->SetText("-1.65423");
    185 
    186         Layout();
    187     }
     885        DisplayBending();
     886        DisplayData();
     887    }
     888    ClassDef(MFit, 0)
    188889};
     890
     891ClassImp(MFit);
    189892
    190893void gui()
    … …  
    192895    new MFit;
    193896}
    194 
Note: See TracChangeset for help on using the changeset viewer.