Context Navigation

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

Visit:

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

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

Download in other formats: