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

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