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

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