Context Navigation

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

Visit:

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

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

Download in other formats: