Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

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

Visit:

Last change on this file since 8823 was 8823, checked in by tbretz, 17 years ago
* empty log message *
File size: 33.7 KB

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

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

Download in other formats: