source: trunk/Cosy/tpoint/TPointGui.cc@ 12958

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