source: trunk/MagicSoft/Cosy/tpoint/gui.C@ 4468

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