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

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