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

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