source: trunk/MagicSoft/Cosy/tpoint/TPointGui.cc @ 9495

Last change on this file since 9495 was 9495, checked in by tbretz, 12 years ago
*** empty log message ***
File size: 39.0 KB
Line 
1#include "TPointGui.h"
2
3#include <fstream>
4
5#include <TGLabel.h>
6#include <TGButton.h>
7#include <TGTextEntry.h>
8
9#include <TView.h>
10#include <TStyle.h>
11#include <TCanvas.h>
12
13#include <TText.h>
14#include <TLine.h>
15#include <TMarker.h>
16#include <TPolyLine.h>
17
18#include <TF1.h>
19#include <TH2.h>
20#include <TMinuit.h>
21#include <TProfile.h>
22#include <TGraphErrors.h>
23
24#include "MLog.h"
25#include "MLogManip.h"
26
27#include "MGList.h"
28
29#include "TPointStar.h"
30
31ClassImp(TPointGui);
32
33using namespace std;
34
35TPointGui::TPointGui(const char *fname, const char *mod) : TGMainFrame(gClient->GetRoot(), 750, 435, kHorizontalFrame), fExitLoopOnClose(kFALSE),
36   fAzMin(-180), fAzMax(180), fZdMin(0), fZdMax(90), fMagMin(0), fLimit(0.05)
37{
38    fCoordinates.SetOwner();
39    fOriginal.SetOwner();
40
41    fList = new MGList;
42    fList->SetOwner();
43
44    fFont = gVirtualX->LoadQueryFont("7x13bold");
45
46    TGLayoutHints *hints0 = new TGLayoutHints(kLHintsExpandY, 7, 5, 5, 0);
47    TGLayoutHints *hints1 = new TGLayoutHints(kLHintsExpandX|kLHintsExpandY, 5, 7, 5, 6);
48    fList->Add(hints0);
49    fList->Add(hints1);
50
51    TGGroupFrame *grp1 = new TGGroupFrame(this, "Control", kVerticalFrame);
52    AddFrame(grp1, hints0);
53    fList->Add(grp1);
54
55    TGGroupFrame *grp2 = new TGGroupFrame(this, "Parameters", kHorizontalFrame);
56    AddFrame(grp2, hints1);
57    fList->Add(grp2);
58
59
60    TGLayoutHints *hints4 = new TGLayoutHints(kLHintsExpandX, 5, 5,  3);
61    TGLayoutHints *hints5 = new TGLayoutHints(kLHintsExpandX, 5, 5, 10);
62    AddTextButton(grp1, "Load Pointing Model", kTbLoad,        hints5);
63    AddTextButton(grp1, "Save Pointing Model", kTbSave,        hints4);
64    AddTextButton(grp1, "Fit Parameters",      kTbFit,         hints5);
65    AddTextButton(grp1, "Reset Parameters",    kTbReset,       hints4);
66    AddTextButton(grp1, "Load Stars",          kTbLoadStars,   hints5);
67    AddTextButton(grp1, "Reset Stars",         kTbResetStars,  hints4);
68    AddTextButton(grp1, "Reload Stars",        kTbReloadStars, hints4);
69    fList->Add(hints4);
70    fList->Add(hints5);
71
72
73
74
75
76
77
78    TGHorizontalFrame *comp = new TGHorizontalFrame(grp2, 1, 1);
79    grp2->AddFrame(comp);
80    fList->Add(comp);
81
82    TGLayoutHints *hints3 = new TGLayoutHints(kLHintsLeft|kLHintsTop, 0, 10, 5, 0);
83    fList->Add(hints3);
84
85    TGVerticalFrame *vframe = new TGVerticalFrame(comp, 1, 1);
86
87    for (int i=0; i<MPointing::GetNumPar(); i++)
88        AddCheckButton(vframe, fBending.GetVarName(i), i);
89
90    TGButton *but = (TGButton*)fList->FindWidget(0);
91
92    comp->AddFrame(vframe, hints3);
93    fList->Add(vframe);
94
95    vframe = new TGVerticalFrame(comp, 1, 1);
96    comp->AddFrame(vframe, hints3);
97    fList->Add(vframe);
98
99    hints3 = new TGLayoutHints(kLHintsLeft|kLHintsTop, 0, 10, 5, 0);
100    fList->Add(hints3);
101
102    TGLabel *l = new TGLabel(vframe, "+000.0000");
103    l->SetTextJustify(kTextRight);
104    fList->Add(l);
105    fLabel.Add(l);
106
107    TGLayoutHints *h = new TGLayoutHints(kLHintsCenterY, 0, 0, but->GetHeight()-l->GetHeight());
108    fList->Add(h);
109
110    vframe->AddFrame(l,h);
111
112    for (int i=1; i<MPointing::GetNumPar(); i++)
113        AddLabel(vframe, "+000.0000", h)->SetTextJustify(kTextRight);
114
115    vframe = new TGVerticalFrame(comp, 1, 1);
116    comp->AddFrame(vframe, hints3);
117    fList->Add(vframe);
118
119    for (int i=0; i<MPointing::GetNumPar(); i++)
120        AddLabel(vframe, "\xb1 00.0000\xb0", h)->SetTextJustify(kTextRight);
121
122    hints3 = new TGLayoutHints(kLHintsLeft|kLHintsTop, 0, 20, 5, 0);
123    fList->Add(hints3);
124
125    TGLayoutHints *hreset = new TGLayoutHints(kLHintsLeft|kLHintsTop, 0, 0, 3, 1);
126    fList->Add(hreset);
127
128    TGVerticalFrame *vframe2 = new TGVerticalFrame(comp, 1, 1);
129    comp->AddFrame(vframe2, hints3);
130    fList->Add(vframe2);
131    for (int i=0; i<MPointing::GetNumPar(); i++)
132        AddResetButton(vframe2, i+2*MPointing::GetNumPar(), hreset,
133                       but->GetHeight()-4);
134
135    vframe = new TGVerticalFrame(comp, 1, 1);
136    comp->AddFrame(vframe, hints3);
137    fList->Add(vframe);
138
139    for (int i=0; i<MPointing::GetNumPar(); i++)
140        AddLabel(vframe, fBending.GetDescription(i), h);
141
142    TGLayoutHints *hints6 = new TGLayoutHints(kLHintsExpandX, 5, 5, 4, 6);
143    fList->Add(hints6);
144
145    l = new TGLabel(grp1, "0000000 Data Sets loaded.");
146    grp1->AddFrame(l, hints6);
147    fList->Add(l);
148    fLabel.Add(l);
149
150    l = new TGLabel(grp1, "");
151    l->SetTextJustify(kTextLeft);
152    grp1->AddFrame(l, hints6);
153    fList->Add(l);
154    fLabel.Add(l);
155
156    l = new TGLabel(grp1, "");
157    l->SetTextJustify(kTextLeft);
158    grp1->AddFrame(l, hints6);
159    fList->Add(l);
160    fLabel.Add(l);
161
162    l = new TGLabel(grp1, "");
163    l->SetTextJustify(kTextLeft);
164    grp1->AddFrame(l, hints6);
165    fList->Add(l);
166    fLabel.Add(l);
167
168    // ------------------------------------------------------------------
169
170    TGLayoutHints *hintse1 = new TGLayoutHints(kLHintsExpandX|kLHintsBottom);
171    TGLayoutHints *hintse2 = new TGLayoutHints(kLHintsExpandX, 2, 2);
172    TGLayoutHints *hintse3 = new TGLayoutHints(kLHintsExpandX);
173    TGLayoutHints *hintsl  = new TGLayoutHints(kLHintsExpandX, 1, 0, 5);
174    fList->Add(hintse1);
175    fList->Add(hintse2);
176    fList->Add(hintse3);
177
178    TGHorizontalFrame *entries = new TGHorizontalFrame(grp1, 1, 1);
179    grp1->AddFrame(entries, hintse1);
180    fList->Add(entries);
181
182    TGVerticalFrame *v1 = new TGVerticalFrame(entries);
183    TGVerticalFrame *v2 = new TGVerticalFrame(entries);
184    TGVerticalFrame *v3 = new TGVerticalFrame(entries);
185    entries->AddFrame(v1, hintse2);
186    entries->AddFrame(v2, hintse2);
187    entries->AddFrame(v3, hintse2);
188    fList->Add(v1);
189    fList->Add(v2);
190    fList->Add(v3);
191
192    TGLabel *label1 = new TGLabel(v1, "Az min/°");
193    TGLabel *label2 = new TGLabel(v2, "Az max/°");
194    TGLabel *label3 = new TGLabel(v3, "Mag min");
195    TGLabel *label4 = new TGLabel(v1, "Zd min/°");
196    TGLabel *label5 = new TGLabel(v2, "Zd max/°");
197    TGLabel *label6 = new TGLabel(v3, "Limit/°");
198    label1->SetTextJustify(kTextLeft);
199    label2->SetTextJustify(kTextLeft);
200    label3->SetTextJustify(kTextLeft);
201    label4->SetTextJustify(kTextLeft);
202    label5->SetTextJustify(kTextLeft);
203    label6->SetTextJustify(kTextLeft);
204    fList->Add(label1);
205    fList->Add(label2);
206    fList->Add(label3);
207    fList->Add(label4);
208    fList->Add(label5);
209    fList->Add(label6);
210
211    TGTextEntry *entry1 = new TGTextEntry(v1, Form("%.1f", fAzMin),  kIdAzMin);
212    TGTextEntry *entry2 = new TGTextEntry(v2, Form("%.1f", fAzMax),  kIdAzMax);
213    TGTextEntry *entry3 = new TGTextEntry(v3, Form("%.1f", fMagMin), kIdMagMin);
214    TGTextEntry *entry4 = new TGTextEntry(v1, Form("%.1f", fZdMin),  kIdZdMin);
215    TGTextEntry *entry5 = new TGTextEntry(v2, Form("%.1f", fZdMax),  kIdZdMax);
216    TGTextEntry *entry6 = new TGTextEntry(v3, Form("%.3f", fLimit),  kIdLimit);
217    entry1->SetToolTipText("TPoints with a real star located at Az<Az min are ignored in the fit.");
218    entry2->SetToolTipText("TPoints with a real star located at Az>Az max are ignored in the fit.");
219    entry2->SetToolTipText("TPoints with a artifiical magnitude Mag<Mag min are ignored in the fit.");
220    entry4->SetToolTipText("TPoints with a real star located at Zd<Zd min are ignored in the fit.");
221    entry5->SetToolTipText("TPoints with a real star located at Zd>Zd max are ignored in the fit.");
222    entry6->SetToolTipText("TPoints with an residual after the fit > Limit are output.");
223    entry1->Associate(this);
224    entry2->Associate(this);
225    entry3->Associate(this);
226    entry4->Associate(this);
227    entry5->Associate(this);
228    entry6->Associate(this);
229    v1->AddFrame(label1, hintsl);
230    v1->AddFrame(entry1, hintse3);
231    v1->AddFrame(label4, hintsl);
232    v1->AddFrame(entry4, hintse3);
233    v2->AddFrame(label2, hintsl);
234    v2->AddFrame(entry2, hintse3);
235    v2->AddFrame(label5, hintsl);
236    v2->AddFrame(entry5, hintse3);
237    v3->AddFrame(label3, hintsl);
238    v3->AddFrame(entry3, hintse3);
239    v3->AddFrame(label6, hintsl);
240    v3->AddFrame(entry6, hintse3);
241    fList->Add(entry1);
242    fList->Add(entry2);
243    fList->Add(entry3);
244    fList->Add(entry4);
245    fList->Add(entry5);
246    fList->Add(entry6);
247
248    // ------------------------------------------------------------------
249
250    // FIXME: Move this to the rc-file.
251    ((TGCheckButton*)fList->FindWidget(0))->SetState(kButtonDown);
252    ((TGCheckButton*)fList->FindWidget(1))->SetState(kButtonDown);
253    ((TGCheckButton*)fList->FindWidget(5))->SetState(kButtonDown);
254    ((TGCheckButton*)fList->FindWidget(6))->SetState(kButtonDown);
255    ((TGCheckButton*)fList->FindWidget(9))->SetState(kButtonDown);
256    ((TGCheckButton*)fList->FindWidget(11))->SetState(kButtonDown);
257    ((TGCheckButton*)fList->FindWidget(15))->SetState(kButtonDown);
258    ((TGCheckButton*)fList->FindWidget(16))->SetState(kButtonDown);
259    ((TGCheckButton*)fList->FindWidget(17))->SetState(kButtonDown);
260
261    ((TGCheckButton*)fList->FindWidget(19))->SetState(kButtonDisabled);
262    ((TGCheckButton*)fList->FindWidget(20))->SetState(kButtonDisabled);
263    ((TGCheckButton*)fList->FindWidget(21))->SetState(kButtonDisabled);
264    ((TGCheckButton*)fList->FindWidget(22))->SetState(kButtonDisabled);
265
266    ((TGCheckButton*)fList->FindWidget(19+2*MPointing::GetNumPar()))->SetState(kButtonDisabled);
267    ((TGCheckButton*)fList->FindWidget(20+2*MPointing::GetNumPar()))->SetState(kButtonDisabled);
268    ((TGCheckButton*)fList->FindWidget(21+2*MPointing::GetNumPar()))->SetState(kButtonDisabled);
269    ((TGCheckButton*)fList->FindWidget(22+2*MPointing::GetNumPar()))->SetState(kButtonDisabled);
270
271    SetWindowName("Telesto");
272    SetIconName("Telesto");
273
274    Layout();
275
276    MapSubwindows();
277    MapWindow();
278
279    if (!TString(fname).IsNull())
280        LoadStars(fname);
281    if (!TString(mod).IsNull())
282        fBending.Load(mod);
283
284    DisplayBending();
285    DisplayData();
286}
287
288TPointGui::~TPointGui()
289{
290    if (fFont)
291        gVirtualX->DeleteFont(fFont);
292
293    delete fList;
294
295    if (fExitLoopOnClose)
296        gSystem->ExitLoop();
297}
298
299void TPointGui::Fcn(Int_t &/*npar*/, Double_t */*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
300{
301    f = 0;
302
303    MPointing bend;
304    bend.SetParameters(par); // Set Parameters [deg] to MPointing
305
306    for (int i=0; i<fCoordinates.GetSize(); i++)
307    {
308        TPointStar set = *(TPointStar*)fCoordinates.At(i);
309
310        if (set.GetStarZd()<fZdMin || set.GetStarZd()>fZdMax ||
311            set.GetStarAz()<fAzMin || set.GetStarAz()>fAzMax ||
312            set.GetMag()   <fMagMin)
313            continue;
314
315        set.Adjust(bend);
316
317        Double_t err = 0.01; // [deg] = 0.25SE
318        Double_t res = set.GetResidual();//(&err);
319        res /= err;
320
321        f += res*res;
322    }
323
324    f /= fCoordinates.GetSize();
325}
326
327void TPointGui::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
328{
329    ((TPointGui*)gMinuit->GetObjectFit())->Fcn(npar, gin, f, par, iflag);
330}
331
332void TPointGui::AddTextButton(TGCompositeFrame *f, TString txt, Int_t id, TGLayoutHints *h)
333{
334    TGButton *but = new TGTextButton(f, txt, id);
335    but->Associate(this);
336    f->AddFrame(but, h);
337    fList->Add(but);
338
339}
340
341void TPointGui::AddCheckButton(TGCompositeFrame *f, TString txt, Int_t id, TGLayoutHints *h)
342{
343    TGButton *but = new TGCheckButton(f, txt, id);
344    but->Associate(this);
345    f->AddFrame(but, h);
346    fList->Add(but);
347}
348
349void TPointGui::AddResetButton(TGCompositeFrame *f, Int_t id, TGLayoutHints *h, Int_t height)
350{
351    TGPictureButton *but = new TGPictureButton(f, "tpoint/skull.xpm", id);
352    but->SetHeight(height); // Offsets from TGLayout
353    but->SetWidth(height);
354    but->Associate(this);
355    f->AddFrame(but, h);
356    fList->Add(but);
357}
358
359TGLabel *TPointGui::AddLabel(TGCompositeFrame *f, TString txt, TGLayoutHints *h)
360{
361    TGLabel *l = new TGLabel(f, txt/*, TGLabel::GetDefaultGC()(), fFont*/);
362    f->AddFrame(l, h);
363    fList->Add(l);
364    fLabel.Add(l);
365    return l;
366}
367
368void TPointGui::DisplayBending()
369{
370    TArrayD par, err;
371    fBending.GetParameters(par);
372    fBending.GetError(err);
373
374    TGLabel *l;
375
376    for (int i=0; i<MPointing::GetNumPar(); i++)
377    {
378        l = (TGLabel*)fLabel.At(i);
379        l->SetText(Form("%.4f\xb0", par[i]));
380
381        l = (TGLabel*)fLabel.At(MPointing::GetNumPar()+i);
382        l->SetText(Form("\xb1 %8.4f\xb0", err[i]>0?err[i]:0));
383    }
384}
385
386void TPointGui::DisplayData()
387{
388    TGLabel *l = (TGLabel*)fLabel.At(3*MPointing::GetNumPar());
389    l->SetText(Form("%d data sets loaded.", fOriginal.GetSize()));
390}
391
392void TPointGui::DisplayResult(Double_t before, Double_t after, Double_t backw)
393{
394    TGLabel *l1 = (TGLabel*)fLabel.At(3*MPointing::GetNumPar()+1);
395    l1->SetText(Form("Before: %.1f +- %.1f SE", before, 0.));
396
397    TGLabel *l2 = (TGLabel*)fLabel.At(3*MPointing::GetNumPar()+2);
398    l2->SetText(Form("After:  %.1f +- %.1f SE", after, 0.));
399
400    TGLabel *l3 = (TGLabel*)fLabel.At(3*MPointing::GetNumPar()+3);
401    l3->SetText(Form("Backw:  %.1f +- %.1f SE", backw, 0.));
402}
403
404void TPointGui::DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0)
405{
406    TView *view = pad->GetView();
407
408    if (!view)
409    {
410        gLog << err << "No View!" << endl;
411        return;
412    }
413
414    TMarker mark0;
415    mark0.SetMarkerStyle(kFullDotLarge);
416    mark0.SetMarkerColor(kBlue);
417
418    r0 /= 90;
419    phi0 *= TMath::DegToRad();
420
421    Double_t x[6] = { r0*cos(phi0), r0*sin(phi0), 0, 0, 0, 0};
422
423    view->WCtoNDC(x, x+3);
424
425    mark0.DrawMarker(-x[3], x[4]);
426}
427
428void TPointGui::DrawPolLine(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
429{
430    TView *view = pad->GetView();
431
432    if (!view)
433    {
434        gLog << err << "No View!" << endl;
435        return;
436    }
437    /*
438    if (r0<0)
439    {
440        r0 = -r0;
441        phi0 += 180;
442    }
443    if (r1<0)
444    {
445        r1 = -r1;
446        phi1 += 180;
447    }
448
449    phi0 = fmod(phi0+360, 360);
450    phi1 = fmod(phi1+360, 360);
451
452    if (phi1-phi0<-180)
453        phi1+=360;
454    */
455    TLine line;
456    line.SetLineWidth(2);
457    line.SetLineColor(kBlue);
458
459    Double_t p0 = phi0<phi1?phi0:phi1;
460    Double_t p1 = phi0<phi1?phi1:phi0;
461
462    if (phi0>phi1)
463    {
464        Double_t d = r1;
465        r1 = r0;
466        r0 = d;
467    }
468
469    r0 /= 90;
470    r1 /= 90;
471
472    Double_t dr = r1-r0;
473    Double_t dp = p1-p0;
474
475    Double_t x0[3] = { r0*cos(p0*TMath::DegToRad()), r0*sin(p0*TMath::DegToRad()), 0};
476
477    for (double i=p0+10; i<p1+10; i+=10)
478    {
479        if (i>p1)
480            i=p1;
481
482        Double_t r = dr/dp*(i-p0)+r0;
483        Double_t p = TMath::DegToRad()*i;
484
485        Double_t x1[3] = { r*cos(p), r*sin(p), 0};
486
487        Double_t y0[3], y1[3];
488
489        view->WCtoNDC(x0, y0);
490        view->WCtoNDC(x1, y1);
491
492        line.DrawLine(y0[0], y0[1], y1[0], y1[1]);
493
494        x0[0] = x1[0];
495        x0[1] = x1[1];
496    }
497}
498
499void TPointGui::DrawSet(TVirtualPad *pad, TPointStar &set, Float_t scale, Float_t angle)
500{
501    Double_t r0   = set.GetRawZd();
502    Double_t phi0 = set.GetRawAz()-angle;
503    Double_t r1   = set.GetStarZd();
504    Double_t phi1 = set.GetStarAz()-angle;
505
506    if (r0<0)
507    {
508        r0 = -r0;
509        phi0 += 180;
510    }
511    if (r1<0)
512    {
513        r1 = -r1;
514        phi1 += 180;
515    }
516
517    phi0 = fmod(phi0+360, 360);
518    phi1 = fmod(phi1+360, 360);
519
520    if (phi1-phi0<-180)
521        phi1+=360;
522
523    if (scale<0 || scale>1000)
524        scale = -1;
525
526    if (scale>0)
527    {
528        Double_t d = r1-r0;
529        r0 += scale*d;
530        r1 -= scale*d;
531        d = phi1-phi0;
532        phi0 += scale*d;
533        phi1 -= scale*d;
534
535        DrawPolLine(pad, r0, phi0, r1, phi1);
536        DrawMarker(pad,  r0, phi0);
537    }
538    else
539        DrawMarker(pad,  r1, phi1);
540}
541
542void TPointGui::DrawHorizon(TVirtualPad *pad, const char *fname) const
543{
544    TView *view = pad->GetView();
545
546    if (!view)
547    {
548        gLog << err << "No View!" << endl;
549        return;
550    }
551
552    ifstream fin(fname);
553    if (!fin)
554    {
555        gLog << err << "ERROR - " << fname << " not found." << endl;
556        return;
557    }
558
559    TPolyLine poly;
560    poly.SetLineWidth(2);
561    poly.SetLineColor(12);
562    poly.SetLineStyle(8);
563
564    while (1)
565    {
566        TString line;
567        line.ReadLine(fin);
568        if (!fin)
569            break;
570
571        Float_t az, alt;
572        sscanf(line.Data(), "%f %f", &az, &alt);
573
574        Float_t zd = 90-alt;
575
576        az *= TMath::DegToRad();
577        zd /= 90;
578
579        Double_t x[6] = { zd*cos(az), zd*sin(az), 0, 0, 0, 0};
580        view->WCtoNDC(x, x+3);
581        poly.SetNextPoint(-x[3], x[4]);
582    }
583
584    poly.DrawClone()->SetBit(kCanDelete);
585
586}
587
588TString TPointGui::OpenDialog(TString &dir, EFileDialogMode mode)
589{
590    static const char *gOpenTypes[] =
591    {
592        "TPoint files",     "*.txt",
593        "Collection files", "*.col",
594        "Model files",      "*.mod",
595        "All files",        "*",
596        NULL,           NULL
597    };
598
599    //static TString dir("tpoint/");
600
601    TGFileInfo fi; // fFileName and fIniDir deleted in ~TGFileInfo
602
603    fi.fFileTypes = (const char**)gOpenTypes;
604    fi.fIniDir    = StrDup(dir);
605
606    new TGFileDialog(fClient->GetRoot(), this, mode, &fi);
607
608    if (!fi.fFilename)
609        return 0;
610
611    dir = fi.fIniDir;
612
613    return fi.fFilename;
614}
615
616void TPointGui::LoadCollection(TString fname)
617{
618    ifstream fin(fname);
619    if (!fin)
620    {
621        gLog << err << "Collection '" << fname << "' not found!" << endl;
622        return;
623    }
624
625    while (1)
626    {
627        TString line;
628        line.ReadLine(fin);
629        if (!fin)
630            break;
631
632        line = line.Strip(TString::kBoth);
633        if (line[0]=='#')
634            continue;
635        if (line.Length()==0)
636            continue;
637
638        if (!line.EndsWith(".txt"))
639        {
640            gLog << warn << "WARNING: " << line << endl;
641            continue;
642        }
643
644        LoadStars(line);
645    }
646}
647
648void TPointGui::LoadStars(TString fname)
649{
650    if (fname.EndsWith(".col"))
651    {
652        LoadCollection(fname);
653        fFileNameStars = fname;
654        SetWindowName(Form("Telesto (%s)", fFileNameStars.Data()));
655        return;
656    }
657
658    const Int_t size = fOriginal.GetSize();
659
660    ifstream fin(fname);
661
662    while (fin && fin.get()!='\n');
663    while (fin && fin.get()!='\n');
664    while (fin && fin.get()!='\n');
665    if (!fin)
666    {
667        gLog << err << "File '" << fname << "' not found!" << endl;
668        return;
669    }
670
671    while (1)
672    {
673        TPointStar set;
674
675        fin >> set;  // Read data from file [deg], it is stored in [rad]
676        if (!fin)
677            break;
678
679//        if (set.GetRawZd()>60)
680//            continue;
681
682        fOriginal.Add(new TPointStar(set));
683    }
684
685    gLog << all << "Found " << fOriginal.GetSize()-size;
686    gLog << " sets of coordinates in " << fname;
687    gLog << " (Total=" << fOriginal.GetSize() << ")" << endl;
688
689    fFileNameStars = fname;
690    SetWindowName(Form("Telesto (%s)", fFileNameStars.Data()));
691}
692
693Float_t TPointGui::GetFloat(Int_t id) const
694{
695    return atof(static_cast<TGTextEntry*>(fList->FindWidget(id))->GetText());
696}
697
698Bool_t TPointGui::ProcessMessage(Long_t msg, Long_t mp1, Long_t)
699{
700    // cout << "Msg: " << hex << GET_MSG(msg) << endl;
701    // cout << "SubMsg: " << hex << GET_SUBMSG(msg) << dec << endl;
702
703    static TString dirmod("tpoint/");
704    static TString dircol("tpoint/");
705
706    switch (GET_MSG(msg))
707    {
708    case kC_COMMAND:
709        switch (GET_SUBMSG(msg))
710        {
711        case kCM_BUTTON:
712            switch (mp1)
713            {
714            case kTbFit:
715                {
716                    Double_t before=0;
717                    Double_t after=0;
718                    Double_t backw=0;
719                    Fit(before, after, backw);
720                    DisplayBending();
721                    DisplayResult(before, after, backw);
722                }
723                return kTRUE;
724            case kTbLoad:
725                fBending.Load(OpenDialog(dirmod));
726                DisplayBending();
727                return kTRUE;
728            case kTbSave:
729                fBending.Save(OpenDialog(dirmod, kFDSave));
730                return kTRUE;
731            case kTbLoadStars:
732                LoadStars(OpenDialog(dircol));
733                DisplayData();
734                return kTRUE;
735            case kTbReset:
736                fBending.Reset();
737                DisplayBending();
738                return kTRUE;
739            case kTbReloadStars:
740                fOriginal.Delete();
741                LoadStars(fFileNameStars); // FIXME: Use TGLabel!
742                DisplayData();
743                return kTRUE;
744            case kTbResetStars:
745                fOriginal.Delete();
746                DisplayData();
747                return kTRUE;
748            }
749
750            // In the default cas a reset button must have been pressed
751            fBending[mp1-2*MPointing::GetNumPar()] = 0;
752            DisplayBending();
753            return kTRUE;
754        }
755        return kTRUE;
756
757    case kC_TEXTENTRY:
758        switch (GET_SUBMSG(msg))
759        {
760        case kTE_TEXTCHANGED:
761            switch (mp1)
762            {
763            case kIdAzMin:
764                fAzMin = GetFloat(kIdAzMin);
765                return kTRUE;
766            case kIdAzMax:
767                fAzMax = GetFloat(kIdAzMax);
768                return kTRUE;
769            case kIdZdMin:
770                fZdMin = GetFloat(kIdZdMin);
771                return kTRUE;
772            case kIdZdMax:
773                fZdMax = GetFloat(kIdZdMax);
774                return kTRUE;
775            case kIdMagMin:
776                fMagMin = GetFloat(kIdMagMin);
777                return kTRUE;
778            case kIdLimit:
779                fLimit = GetFloat(kIdLimit);
780                return kTRUE;
781            }
782            return kTRUE;
783
784        }
785        return kTRUE;
786
787    }
788    return kTRUE;
789}
790
791void TPointGui::Fit(Double_t &before, Double_t &after, Double_t &backw)
792{
793    if (fOriginal.GetSize()==0)
794    {
795        gLog << warn << "Sorry, no input data loaded..." << endl;
796        return;
797    }
798
799    fCoordinates.Delete();
800    for (int i=0; i<fOriginal.GetSize(); i++)
801        fCoordinates.Add(new TPointStar(*(TPointStar*)fOriginal.At(i)));
802
803    gLog << all << "-----------------------------------------------------------------------" << endl;
804
805    gStyle->SetOptStat("emro");
806
807    TH1F hres1("Res1", " Residuals before correction ", fOriginal.GetSize()/3, 0, 0.3);
808    TH1F hres2("Res2", " Residuals after correction ",  fOriginal.GetSize()/3, 0, 0.3);
809    TH1F hres3("Res3", " Residuals after backward correction ",  fOriginal.GetSize()/3, 0, 0.3);
810
811    TProfile proaz ("ProAz",  " \\Delta profile vs. Az",  48, -180, 180);
812    TProfile prozd ("ProZd",  " \\Delta profile vs. Zd",  60,    0, 90);
813    TProfile promag("ProMag", " \\Delta profile vs. Mag", 40,    1,  4);
814
815    hres1.SetXTitle("\\Delta [\\circ]");
816    hres1.SetYTitle("Counts");
817
818    hres2.SetXTitle("\\Delta [\\circ]");
819    hres2.SetYTitle("Counts");
820
821    hres3.SetXTitle("\\Delta [\\circ]");
822    hres3.SetYTitle("Counts");
823
824    TGraph gdaz;
825    TGraph gdzd;
826    TGraph gaz;
827    TGraph gzd;
828    TGraphErrors graz;
829    TGraphErrors grzd;
830    TGraphErrors grmag;
831    TGraph gmaz;
832    TGraph gmzd;
833
834    gdaz.SetTitle(" \\Delta Az vs. Zd ");
835    gdzd.SetTitle(" \\Delta Zd vs. Az ");
836
837    gaz.SetTitle(" \\Delta Az vs. Az ");
838    gzd.SetTitle(" \\Delta Zd vs. Zd ");
839
840    gmaz.SetTitle(" \\Delta Az vs. Mag ");
841    gmzd.SetTitle(" \\Delta Zd vs. Mag ");
842
843    graz.SetTitle(" \\Delta vs. Az ");
844    grzd.SetTitle(" \\Delta vs. Zd ");
845    grmag.SetTitle(" \\Delta vs. Mag ");
846
847    TMinuit minuit(MPointing::GetNumPar());  //initialize TMinuit with a maximum of 5 params
848    minuit.SetObjectFit(this);
849    minuit.SetPrintLevel(-1);
850    minuit.SetFCN(fcn);
851
852    fBending.SetMinuitParameters(minuit, MPointing::GetNumPar()); // Init Parameters [deg]
853
854    for (int i=0; i<MPointing::GetNumPar(); i++)
855    {
856        TGButton *l = (TGButton*)fList->FindWidget(i);
857        minuit.FixParameter(i);
858        if (l->GetState()==kButtonDown)
859            minuit.Release(i);
860    }
861
862    //minuit.Command("SHOW PARAMETERS");
863    //minuit.Command("SHOW LIMITS");
864
865    gLog << inf << endl;
866    gLog << "Starting fit..." << endl;
867    gLog << "For the fit an measurement error in the residual of ";
868    gLog << "0.02deg (=1SE) is assumed." << endl;
869    gLog << endl;
870
871    Int_t ierflg = 0;
872    ierflg = minuit.Migrad();
873    gLog << inf << "Migrad returns " << ierflg << endl;
874    // minuit.Release(2);
875    ierflg = minuit.Migrad();
876    gLog << inf << "Migrad returns " << ierflg << endl << endl;
877
878    //
879    // Get Fit Results
880    //
881    fBending.GetMinuitParameters(minuit);
882    fBending.PrintMinuitParameters(minuit);
883    gLog << endl;
884    //fBending.Save("bending_magic.txt");
885
886
887    //
888    // Make a copy of all list entries
889    //
890    TList list;
891    list.SetOwner();
892    for (int i=0; i<fCoordinates.GetSize(); i++)
893        list.Add(new TPointStar(*(TPointStar*)fCoordinates.At(i)));
894
895    //
896    // Correct for Offsets only
897    //
898    TArrayD par;
899    fBending.GetParameters(par);
900    for (int i=2; i<MPointing::GetNumPar(); i++)
901        par[i]=0;
902
903    MPointing b2;
904    b2.SetParameters(par);
905
906    gLog << all << endl << "Sets with Residual exceeding " << fLimit << "deg:" << endl;
907
908    //
909    // Calculate correction and residuals
910    //
911    for (int i=0; i<fCoordinates.GetSize(); i++)
912    {
913        TPointStar orig = *(TPointStar*)fCoordinates.At(i);
914
915        TPointStar &set0 = *(TPointStar*)fCoordinates.At(i);
916
917        ZdAz za(set0.GetStarZdAz());
918        za *=kRad2Deg;
919
920        //
921        // Correct for offsets only
922        //
923        TPointStar set1(set0);
924        set1.Adjust(b2);
925
926        hres1.Fill(set1.GetResidual());
927
928        set0.Adjust(fBending);
929        hres2.Fill(set0.GetResidual());
930
931        Double_t dz = fmod(set0.GetDAz()+720, 360);
932        if (dz>180)
933            dz -= 360;
934
935        Double_t err;
936        Double_t resi = set0.GetResidual(&err);
937
938        gdzd.SetPoint(i, za.Az(), set0.GetDZd());
939        gdaz.SetPoint(i, za.Zd(), dz);
940        graz.SetPoint(i, za.Az(), resi);
941        graz.SetPointError(i, 0, err);
942        grzd.SetPoint(i, za.Zd(), resi);
943        grzd.SetPointError(i, 0, err);
944
945        if (resi>fLimit) // 0.13
946            gLog << all << " " << orig << "  <" << resi << ">" << endl;
947
948        proaz.Fill(za.Az(), set0.GetResidual(&err));
949        prozd.Fill(za.Zd(), set0.GetResidual(&err));
950        promag.Fill(set0.GetMag(), set0.GetResidual(&err));
951
952        gaz.SetPoint( i, za.Az(), dz);
953        gzd.SetPoint( i, za.Zd(), set0.GetDZd());
954        if (set0.GetMag()>=-20)
955        {
956            grmag.SetPoint(i, set0.GetMag(), set0.GetResidual(&err));
957            grmag.SetPointError(i, 0, err);
958            gmaz.SetPoint( i, set0.GetMag(), dz);
959            gmzd.SetPoint( i, set0.GetMag(), set0.GetDZd());
960        }
961    }
962
963    gLog << "done." << endl << endl;
964
965    //
966    // Check for overflows
967    //
968    const Stat_t ov = hres2.GetBinContent(hres2.GetNbinsX()+1);
969    if (ov>0)
970        gLog << warn << "WARNING: " << ov << " overflows in residuals." << endl;
971
972
973
974    gLog << inf << dec << endl;
975    gLog << "              Number of calls to FCN: " << minuit.fNfcn << endl;
976    gLog << "Minimum value found for FCN (Chi^2?): " << minuit.fAmin << endl;
977    gLog << "                  Fit-Probability(?): " << TMath::Prob(minuit.fAmin/*fOriginal.GetSize()*/, fOriginal.GetSize()-minuit.GetNumFreePars())*100 << "%" << endl;
978    gLog << "                           Chi^2/NDF: " << minuit.fAmin/(fOriginal.GetSize()-minuit.GetNumFreePars()) << endl;
979    //cout << "Prob(?): " << TMath::Prob(fChisquare,ndf);
980
981
982
983    //
984    // Print all data sets for which the backward correction is
985    // twice times worse than the residual gotten from the
986    // bending correction itself
987    //
988    gLog << endl;
989    gLog << "Checking backward correction (raw-->star):" << endl;
990    for (int i=0; i<fCoordinates.GetSize(); i++)
991    {
992        TPointStar set0(*(TPointStar*)list.At(i));
993        TPointStar &set1 = *(TPointStar*)list.At(i);
994
995        set0.AdjustBack(fBending);
996        set1.Adjust(fBending);
997
998        const Double_t res0 = set0.GetResidual();
999        const Double_t res1 = set1.GetResidual();
1000        const Double_t diff = TMath::Abs(res0-res1);
1001
1002        hres3.Fill(res0);
1003
1004        if (diff<hres2.GetMean()*0.66)
1005            continue;
1006
1007        gLog << "DBack: " << setw(6) << set0.GetStarZd() << " " << setw(7) << set0.GetStarAz() << ":  ";
1008        gLog << "ResB="<< setw(7) << res0*60 << "  ResF=" << setw(7) << res1*60 << "  |ResB-ResF|=" << setw(7) << diff*60 << " arcmin" << endl;
1009    }
1010    gLog << "OK." << endl;
1011    gLog << endl;
1012
1013    const Double_t max1 = TMath::Max(gaz.GetHistogram()->GetMaximum(), gdaz.GetHistogram()->GetMaximum());
1014    const Double_t max2 = TMath::Max(gzd.GetHistogram()->GetMaximum(), gdzd.GetHistogram()->GetMaximum());
1015    const Double_t max3 = TMath::Max(grzd.GetHistogram()->GetMaximum(), graz.GetHistogram()->GetMaximum());
1016
1017    const Double_t min1 = TMath::Min(gaz.GetHistogram()->GetMinimum(), gdaz.GetHistogram()->GetMinimum());
1018    const Double_t min2 = TMath::Min(gzd.GetHistogram()->GetMinimum(), gdzd.GetHistogram()->GetMinimum());
1019    const Double_t min3 = TMath::Min(grzd.GetHistogram()->GetMinimum(), graz.GetHistogram()->GetMinimum());
1020
1021    const Double_t absmax1 = TMath::Max(max1, TMath::Abs(min1));
1022    const Double_t absmax2 = TMath::Max(max2, TMath::Abs(min2));
1023    const Double_t absmax3 = TMath::Max(max3, TMath::Abs(min3));
1024
1025    gaz.SetMaximum(absmax1);
1026    gzd.SetMaximum(absmax2);
1027    gdaz.SetMaximum(absmax1);
1028    gdzd.SetMaximum(absmax2);
1029    gmaz.SetMaximum(absmax1);
1030    gmzd.SetMaximum(absmax2);
1031    graz.SetMaximum(absmax3);
1032    grzd.SetMaximum(absmax3);
1033    grmag.SetMaximum(absmax3);
1034    gaz.SetMinimum(-absmax1);
1035    gzd.SetMinimum(-absmax2);
1036    gdaz.SetMinimum(-absmax1);
1037    gdzd.SetMinimum(-absmax2);
1038    gmaz.SetMinimum(-absmax1);
1039    gmzd.SetMinimum(-absmax2);
1040    graz.SetMinimum(0);
1041    grzd.SetMinimum(0);
1042    grmag.SetMinimum(0);
1043
1044    TCanvas *c1;
1045
1046    if (gROOT->FindObject("CanvGraphs"))
1047        c1 = dynamic_cast<TCanvas*>(gROOT->FindObject("CanvGraphs"));
1048    else
1049        c1=new TCanvas("CanvGraphs", "Graphs");
1050
1051    gROOT->SetSelectedPad(0);
1052    c1->SetSelectedPad(0);
1053    c1->SetBorderMode(0);
1054    c1->SetFrameBorderMode(0);
1055    c1->Clear();
1056
1057    c1->SetFillColor(kWhite);
1058#ifndef PRESENTATION
1059    c1->Divide(3,3,1e-10,1e-10);
1060#else
1061    c1->Divide(2,2,1e-10,1e-10);
1062#endif
1063    c1->SetFillColor(kWhite);
1064
1065    TGraph *g=0;
1066
1067    TLine line;
1068    line.SetLineColor(kGreen);
1069    line.SetLineWidth(2);
1070#ifndef PRESENTATION
1071    c1->cd(1);
1072    gPad->SetBorderMode(0);
1073    gPad->SetFrameBorderMode(0);
1074    gPad->SetGridx();
1075    gPad->SetGridy();
1076    g=(TGraph*)gaz.DrawClone("A*");
1077    g->SetBit(kCanDelete);
1078    g->GetHistogram()->SetXTitle("Az [\\circ]");
1079    g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]");
1080
1081    line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
1082    line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
1083
1084    c1->cd(2);
1085    gPad->SetBorderMode(0);
1086    gPad->SetFrameBorderMode(0);
1087    gPad->SetGridx();
1088    gPad->SetGridy();
1089    g=(TGraph*)gdaz.DrawClone("A*");
1090    g->SetBit(kCanDelete);
1091    g->GetHistogram()->SetXTitle("Zd [\\circ]");
1092    g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]");
1093    line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
1094    line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
1095    cout << "Mean dAz: " << g->GetMean(2) << " \xb1 " << g->GetRMS(2) <<  endl;
1096
1097    c1->cd(3);
1098    gPad->SetBorderMode(0);
1099    gPad->SetFrameBorderMode(0);
1100    gPad->SetGridx();
1101    gPad->SetGridy();
1102    if (gmaz.GetN()>0)
1103    {
1104        g=(TGraph*)gmaz.DrawClone("A*");
1105        g->SetBit(kCanDelete);
1106        g->GetHistogram()->SetXTitle("Mag");
1107        g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]");
1108        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
1109        line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
1110    }
1111#endif
1112
1113#ifndef PRESENTATION
1114    c1->cd(4);
1115#else
1116    c1->cd(1);
1117#endif
1118    gPad->SetBorderMode(0);
1119    gPad->SetFrameBorderMode(0);
1120    gPad->SetGridx();
1121    gPad->SetGridy();
1122    g=(TGraph*)gdzd.DrawClone("A*");
1123    g->SetBit(kCanDelete);
1124    g->GetHistogram()->SetXTitle("Az [\\circ]");
1125    g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]");
1126    line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
1127    line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
1128    cout << "Mean dZd: " << g->GetMean(2) << " \xb1 " << g->GetRMS(2) <<  endl;
1129    cout << endl;
1130
1131#ifndef PRESENTATION
1132    c1->cd(5);
1133#else
1134    c1->cd(2);
1135#endif
1136    gPad->SetBorderMode(0);
1137    gPad->SetFrameBorderMode(0);
1138    gPad->SetGridx();
1139    gPad->SetGridy();
1140    g=(TGraph*)gzd.DrawClone("A*");
1141    g->SetBit(kCanDelete);
1142    g->GetHistogram()->SetXTitle("Zd [\\circ]");
1143    g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]");
1144    line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
1145    line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
1146#ifndef PRESENTATION
1147    c1->cd(6);
1148    gPad->SetBorderMode(0);
1149    gPad->SetFrameBorderMode(0);
1150    gPad->SetGridx();
1151    gPad->SetGridy();
1152    if (gmzd.GetN()>0)
1153    {
1154        g=(TGraph*)gmzd.DrawClone("A*");
1155        g->SetBit(kCanDelete);
1156        g->GetHistogram()->SetXTitle("Mag");
1157        g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]");
1158        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
1159        line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
1160    }
1161#endif
1162
1163#ifndef PRESENTATION
1164    c1->cd(7);
1165#else
1166    c1->cd(3);
1167#endif
1168    gPad->SetBorderMode(0);
1169    gPad->SetFrameBorderMode(0);
1170    gPad->SetGridx();
1171    gPad->SetGridy();
1172    g=(TGraph*)graz.DrawClone("AP");
1173    g->SetBit(kCanDelete);
1174    g->GetHistogram()->SetXTitle("Az [\\circ]");
1175    g->GetHistogram()->SetYTitle("\\Delta [\\circ]");
1176    line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
1177
1178    proaz.SetLineWidth(2);
1179    proaz.SetLineColor(kBlue);
1180    proaz.SetMarkerColor(kBlue);
1181    proaz.DrawCopy("pc hist same");
1182
1183#ifndef PRESENTATION
1184    c1->cd(8);
1185#else
1186    c1->cd(4);
1187#endif
1188    gPad->SetBorderMode(0);
1189    gPad->SetFrameBorderMode(0);
1190    gPad->SetGridx();
1191    gPad->SetGridy();
1192    g=(TGraph*)grzd.DrawClone("AP");
1193    g->SetBit(kCanDelete);
1194    g->GetHistogram()->SetXTitle("Zd [\\circ]");
1195    g->GetHistogram()->SetYTitle("\\Delta [\\circ]");
1196    line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
1197
1198    prozd.SetLineWidth(2);
1199    prozd.SetLineColor(kBlue);
1200    prozd.SetMarkerColor(kBlue);
1201    prozd.DrawCopy("pc hist same");
1202
1203#ifndef PRESENTATION
1204    c1->cd(9);
1205    gPad->SetBorderMode(0);
1206    gPad->SetFrameBorderMode(0);
1207    gPad->SetGridx();
1208    gPad->SetGridy();
1209    if (grmag.GetN()>0)
1210    {
1211        g=(TGraph*)grmag.DrawClone("AP");
1212        g->SetBit(kCanDelete);
1213        g->GetHistogram()->SetXTitle("Mag");
1214        g->GetHistogram()->SetYTitle("\\Delta [\\circ]");
1215        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
1216    }
1217    promag.SetLineWidth(2);
1218    promag.SetLineColor(kBlue);
1219    promag.SetMarkerColor(kBlue);
1220    promag.DrawCopy("pc hist same");
1221#endif
1222
1223    //
1224    // Print out the residual before and after correction in several
1225    // units
1226    //
1227    cout << fCoordinates.GetSize() << " data sets." << endl << endl;
1228    cout << "Total Spread of Residual:" << endl;
1229    cout << "-------------------------" << endl;
1230    cout << "before: " << Form("%6.4f", hres1.GetMean()) << " \xb1 " << Form("%6.4f", hres1.GetRMS()) << " deg \t";
1231    cout << "before: " << Form("%4.1f", hres1.GetMean()*60) << " \xb1 " << Form("%.1f", hres1.GetRMS()*60) << " arcmin" << endl;
1232    cout << "after:  " << Form("%6.4f", hres2.GetMean()) << " \xb1 " << Form("%6.4f", hres2.GetRMS()) << " deg \t";
1233    cout << "after:  " << Form("%4.1f", hres2.GetMean()*60) << " \xb1 " << Form("%.1f", hres2.GetRMS()*60) << " arcmin" << endl;
1234    cout << "backw:  " << Form("%6.4f", hres3.GetMean()) << " \xb1 " << Form("%6.4f", hres3.GetRMS()) << " deg \t";
1235    cout << "backw:  " << Form("%4.1f", hres3.GetMean()*60) << " \xb1 " << Form("%.1f", hres3.GetRMS()*60) << " arcmin" << endl;
1236    cout << endl;
1237    cout << "before: " << Form("%4.1f", hres1.GetMean()*16348/360) << " \xb1 " << Form("%.1f", hres1.GetRMS()*16384/360) << " SE \t\t";
1238    cout << "before: " << Form("%4.1f", hres1.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres1.GetRMS()*60*60/23.4) << " pix" << endl;
1239    cout << "after:  " << Form("%4.1f", hres2.GetMean()*16384/360) << " \xb1 " << Form("%.1f", hres2.GetRMS()*16384/360) << " SE \t\t";
1240    cout << "after:  " << Form("%4.1f", hres2.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres2.GetRMS()*60*60/23.4) << " pix" << endl;
1241    cout << "backw:  " << Form("%4.1f", hres3.GetMean()*16384/360) << " \xb1 " << Form("%.1f", hres3.GetRMS()*16384/360) << " SE \t\t";
1242    cout << "backw:  " << Form("%4.1f", hres3.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres3.GetRMS()*60*60/23.4) << " pix" << endl;
1243    cout << endl;
1244    cout << endl; // ±
1245
1246
1247    before = hres1.GetMean()*16384/360;
1248    after  = hres2.GetMean()*16384/360;
1249    backw  = hres3.GetMean()*16384/360;
1250
1251
1252    gStyle->SetOptStat(1110);
1253    gStyle->SetStatFormat("6.2g");
1254
1255    if (gROOT->FindObject("CanvResiduals"))
1256        c1 = dynamic_cast<TCanvas*>(gROOT->FindObject("CanvResiduals"));
1257    else
1258        c1=new TCanvas("CanvResiduals", "Residuals", 800, 800);
1259
1260    gROOT->SetSelectedPad(0);
1261    c1->SetSelectedPad(0);
1262    c1->Clear();
1263    c1->SetFillColor(kWhite);
1264
1265    c1->Divide(2, 2, 1e-10, 1e-10);
1266
1267    c1->cd(2);
1268    gPad->SetBorderMode(0);
1269    gPad->SetFrameBorderMode(0);
1270    hres1.SetLineColor(kRed);
1271    hres1.DrawCopy();
1272
1273    gPad->Update();
1274
1275    line.DrawLine(360./16384, gPad->GetUymin(), 360./16384, gPad->GetUymax());
1276
1277    c1->cd(4);
1278    gPad->SetBorderMode(0);
1279    gPad->SetFrameBorderMode(0);
1280    hres2.SetLineColor(kBlue);
1281    TH1 *h=hres2.DrawCopy();
1282    TF1 f("mygaus", "(gaus)", 0, 1);
1283    f.SetLineColor(kMagenta/*6*/);
1284    f.SetLineWidth(1);
1285    f.SetParameter(0, h->GetBinContent(1));
1286    f.FixParameter(1, 0);
1287    f.SetParameter(2, h->GetRMS());
1288    h->Fit("mygaus", "QR");
1289    hres3.SetLineColor(kCyan);
1290    hres3.SetLineStyle(kDashed);
1291    hres3.DrawCopy("same");
1292    cout << "Gaus-Fit  Sigma: " << f.GetParameter(2) << "\xb0" << endl;
1293    cout << "Fit-Probability: " << f.GetProb()*100 << "%" << endl;
1294    cout << "      Chi^2/NDF: " << f.GetChisquare() << "/" << f.GetNDF() << " = " << f.GetChisquare()/f.GetNDF() << endl;
1295    gPad->Update();
1296    line.DrawLine(360./16384, gPad->GetUymin(), 360./16384, gPad->GetUymax());
1297
1298    c1->cd(1);
1299    gPad->SetBorderMode(0);
1300    gPad->SetFrameBorderMode(0);
1301    gPad->SetTheta(90);
1302    gPad->SetPhi(90);
1303    TH2F h2res1("Res2D1", " Dataset positions on the sky ", 36, 0, 360,  8, 0, 90);
1304    h2res1.SetBit(TH1::kNoStats);
1305    h2res1.DrawCopy("surf1pol");
1306    gPad->Modified();
1307    gPad->Update();
1308    DrawHorizon(gPad);
1309    for (int i=0; i<fOriginal.GetSize(); i++)
1310        DrawSet(gPad, *(TPointStar*)fOriginal.At(i));//, 10./hres1.GetMean());
1311
1312    TText text;
1313    text.SetTextAlign(22);
1314    text.DrawText( 0.00,  0.66, "N");
1315    text.DrawText( 0.66,  0.00, "E");
1316    text.DrawText( 0.00, -0.66, "S");
1317    text.DrawText(-0.66,  0.00, "W");
1318
1319    c1->cd(3);
1320    gPad->SetBorderMode(0);
1321    gPad->SetFrameBorderMode(0);
1322    gPad->SetTheta(90);
1323    gPad->SetPhi(90);
1324    h2res1.SetTitle(" Arb. Residuals after correction (scaled) ");
1325    h2res1.DrawCopy("surf1pol");
1326    gPad->Modified();
1327    gPad->Update();
1328//        for (int i=0; i<fCoordinates.GetSize(); i++)
1329//            DrawSet(gPad, *(Set*)fCoordinates.At(i), 10./hres2.GetMean(), par[0]);
1330
1331    RaiseWindow();
1332}
1333
Note: See TracBrowser for help on using the repository browser.