Context Navigation

TPointGui.cc@ 18678

Visit:

Last change on this file since 18678 was 18625, checked in by tbretz, 8 years ago
Removed reference to root dictionary, some minor update to disentengle from root and Mars as much as possible, changed path to horizon.dat
File size: 39.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/FACT++/drive/TPointGui.cc@ 18678

Download in other formats: