Context Navigation

TPointGui.cc@ 9844

Visit:

Last change on this file since 9844 was 9551, checked in by tbretz, 15 years ago
* empty log message *
File size: 39.0 KB

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

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

Context Navigation

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

Download in other formats: