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

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