Context Navigation

TPointGui.cc@ 17556

Visit:

Last change on this file since 17556 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
33	ClassImp(TPointGui);
34
35	using namespace std;
36
37	TPointGui::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+2MPointing::GetNumPar()))->SetState(kButtonDisabled);
269	((TGCheckButton)fList->FindWidget(20+2MPointing::GetNumPar()))->SetState(kButtonDisabled);
270	((TGCheckButton)fList->FindWidget(21+2MPointing::GetNumPar()))->SetState(kButtonDisabled);
271	((TGCheckButton)fList->FindWidget(22+2MPointing::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
290	TPointGui::~TPointGui()
291	{
292	if (fFont)
293	gVirtualX->DeleteFont(fFont);
294
295	delete fList;
296
297	if (fExitLoopOnClose)
298	gSystem->ExitLoop();
299	}
300
301	void 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
329	void 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
334	void 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
343	void 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
351	void 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
361	TGLabel 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
370	void 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
388	void TPointGui::DisplayData()
389	{
390	TGLabel l = (TGLabel)fLabel.At(3*MPointing::GetNumPar());
391	l->SetText(Form("%d data sets loaded.", fOriginal.GetSize()));
392	}
393
394	void 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
406	void 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] = { r0cos(phi0), r0sin(phi0), 0, 0, 0, 0};
424
425	view->WCtoNDC(x, x+3);
426
427	mark0.DrawMarker(-x[3], x[4]);
428	}
429
430	void 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] = { r0cos(p0TMath::DegToRad()), r0sin(p0TMath::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] = { rcos(p), rsin(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
501	void 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
544	void 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] = { zdcos(az), zdsin(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
590	TString 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
618	void 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
650	void 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
695	Float_t TPointGui::GetFloat(Int_t id) const
696	{
697	return atof(static_cast<TGTextEntry*>(fList->FindWidget(id))->GetText());
698	}
699
700	Bool_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
793	void 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) << res060 << " ResF=" << setw(7) << res160 << " \|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()6060/23.4) << " \xb1 " << Form("%.1f", hres1.GetRMS()6060/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()6060/23.4) << " \xb1 " << Form("%.1f", hres2.GetRMS()6060/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()6060/23.4) << " \xb1 " << Form("%.1f", hres3.GetRMS()6060/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.

Context Navigation

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

Download in other formats: