Context Navigation

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

Visit:

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

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

Download in other formats: