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

Last change on this file since 4308 was 4256, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 31.7 KB
Line 
1#include <fstream.h>
2#include <fstream.h>
3#include <fstream.h>
4#include <iostream.h>
5#include <iomanip.h>
6
7#include <TError.h>
8
9#include <TGFrame.h>
10#include <TGLabel.h>
11#include <TGButton.h>
12
13#include <TF1.h>
14#include <TH1.h>
15#include <TH2.h>
16#include <TGraphErrors.h>
17
18#include <TList.h>
19#include <TStyle.h>
20#include <TMinuit.h>
21
22#include <TView.h>
23#include <TLine.h>
24#include <TMarker.h>
25#include <TCanvas.h>
26
27#include "coord.h"
28
29#include "MGList.h"
30#include "MBending.h"
31
32class Set : public TObject
33{
34 friend ifstream &operator>>(ifstream &fin, Set &set);
35private:
36 Double_t fStarAz;
37 Double_t fStarEl;
38
39 Double_t fRawAz;
40 Double_t fRawEl;
41
42public:
43 Set(Double_t sel=0, Double_t saz=0, Double_t rel=0, Double_t raz=0) :
44 fStarAz(saz*TMath::DegToRad()),
45 fStarEl(sel*TMath::DegToRad()),
46 fRawAz(raz*TMath::DegToRad()),
47 fRawEl(rel*TMath::DegToRad())
48 {
49 }
50
51 Double_t GetResidual(Double_t *err=0) const
52 {
53 /*
54 TVector3 v1, v2;
55 v1.SetMagThetaPhi(1, TMath::Pi()/2-fRawEl, fRawAz);
56 v2.SetMagThetaPhi(1, TMath::Pi()/2-fStarEl, fStarAz);
57
58 return v1.Angle(v2)*TMath::RadToDeg();
59 */
60
61 const Double_t del = fRawEl-fStarEl;
62 const Double_t daz = fRawAz-fStarAz;
63 /*
64 const Double_t dphi2 = daz/2.;
65 const Double_t cos2 = cos(dphi2)*cos(dphi2);
66 const Double_t sin2 = sin(dphi2)*sin(dphi2);
67 const Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
68 */
69
70 const Double_t d = cos(del) - cos(fRawEl)*cos(fStarEl)*(1.-cos(daz));
71
72 if (err)
73 {
74 // Error of one pixel in the CCD
75 const Double_t e1 = 32./3600*TMath::DegToRad() * 0.5;
76
77 // Error of one SE unit
78 const Double_t e2 = 360./16384*TMath::DegToRad() * 0.5;
79
80 const Double_t e11 = sin(del)+cos(fRawEl)*sin(fStarEl)*(1-cos(daz));
81 const Double_t e12 = cos(fRawEl)*cos(fStarEl)*sin(daz);
82
83 const Double_t e21 = -sin(del)+sin(fRawEl)*cos(fStarEl)*(1-cos(daz));
84 const Double_t e22 = -cos(fRawEl)*cos(fStarEl)*sin(daz);
85
86 const Double_t err1 = sqrt(1-d*d);
87 const Double_t err2 = (e11*e11 + e12*e12)*e1*e1;
88 const Double_t err3 = (e21*e21 + e22*e22)*e2*e2;
89
90 *err = sqrt(err2+err3)/err1 * TMath::RadToDeg();
91 }
92
93 const Double_t dist = acos(d);
94 return dist * TMath::RadToDeg();
95 }
96
97 void operator=(Set &set)
98 {
99 fStarAz = set.fStarAz;
100 fStarEl = set.fStarEl;
101 fRawAz = set.fRawAz;
102 fRawEl = set.fRawEl;
103 }
104
105 Double_t GetDEl() const { return (fRawEl-fStarEl)*TMath::RadToDeg(); }
106 Double_t GetDZd() const { return -GetDEl(); }
107 Double_t GetDAz() const { return (fRawAz-fStarAz)*TMath::RadToDeg(); }
108 Double_t GetStarEl() const { return fStarEl*TMath::RadToDeg(); }
109 Double_t GetStarZd() const { return 90.-fStarEl*TMath::RadToDeg(); }
110 Double_t GetStarAz() const { return fStarAz*TMath::RadToDeg(); }
111 Double_t GetRawEl() const { return fRawEl*TMath::RadToDeg(); }
112 Double_t GetRawAz() const { return fRawAz*TMath::RadToDeg(); }
113 Double_t GetRawZd() const { return 90.-fRawEl*TMath::RadToDeg(); }
114
115 ZdAz GetStarZdAz() const { return ZdAz(TMath::Pi()/2-fStarEl, fStarAz); }
116 AltAz GetStarAltAz() const { return AltAz(fStarEl, fStarAz); }
117
118 ZdAz GetRawZdAz() const { return ZdAz(TMath::Pi()/2-fRawEl, fRawAz); }
119 AltAz GetRawAltAz() const { return AltAz(fRawEl, fRawAz); }
120
121 void AdjustEl(Double_t del) { fStarEl += del*TMath::DegToRad(); }
122 void AdjustAz(Double_t daz) { fStarAz += daz*TMath::DegToRad(); }
123
124 void Adjust(const MBending &bend)
125 {
126 AltAz p = bend(GetStarAltAz());
127 fStarEl = p.Alt();
128 fStarAz = p.Az();
129 }
130 void AdjustBack(const MBending &bend)
131 {
132 AltAz p = bend.CorrectBack(GetRawAltAz());
133 fRawEl = p.Alt();
134 fRawAz = p.Az();
135 }
136 ClassDef(Set, 0)
137};
138
139ClassImp(Set);
140
141ifstream &operator>>(ifstream &fin, Set &set)
142{
143 TString str;
144 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 // Check for overflows
502 //
503 const Stat_t ov = hres2.GetBinContent(hres2.GetNbinsX()+1);
504 if (ov>0)
505 cout << "WARNING: " << ov << " overflows in residuals." << endl;
506
507
508
509 cout << dec << endl;
510 cout << " Number of calls to FCN: " << minuit.fNfcn << endl;
511 cout << "Minimum value found for FCN (Chi^2?): " << minuit.fAmin << endl;
512 cout << " Fit-Probability(?): " << TMath::Prob(minuit.fAmin/*fOriginal.GetSize()*/, fOriginal.GetSize()-minuit.GetNumFreePars())*100 << "%" << endl;
513 cout << " Chi^2/NDF: " << minuit.fAmin/(fOriginal.GetSize()-minuit.GetNumFreePars()) << endl;
514 //cout << "Prob(?): " << TMath::Prob(fChisquare,ndf);
515
516
517
518 //
519 // Print all data sets for which the backward correction is
520 // twice times worse than the residual gotten from the
521 // bending correction itself
522 //
523 cout << endl;
524 cout << "Checking backward correction (raw-->star):" << endl;
525 for (int i=0; i<fCoordinates.GetSize(); i++)
526 {
527 Set set0(*(Set*)list.At(i));
528 Set &set1 = *(Set*)list.At(i);
529
530 set0.AdjustBack(fBending);
531 set1.Adjust(fBending);
532
533 const Double_t res0 = set0.GetResidual();
534 const Double_t res1 = set1.GetResidual();
535 const Double_t diff = TMath::Abs(res0-res1);
536
537 hres3.Fill(res0);
538
539 if (diff<hres2.GetMean()*0.66)
540 continue;
541
542 cout << "DBack: " << setw(6) << set0.GetStarZd() << " " << setw(7) << set0.GetStarAz() << ": ";
543 cout << "ResB="<< setw(7) << res0*60 << " ResF=" << setw(7) << res1*60 << " |ResB-ResF|=" << setw(7) << diff*60 << " arcmin" << endl;
544 }
545 cout << "OK." << endl;
546 cout << endl;
547
548
549 TCanvas *c1;
550
551 if ((c1 = (TCanvas*)gROOT->FindObject("CanvGraphs")))
552 delete c1;
553 if ((c1 = (TCanvas*)gROOT->FindObject("CanvResiduals")))
554 delete c1;
555
556
557 const Double_t max1 = TMath::Max(gaz.GetHistogram()->GetMaximum(), gdaz.GetHistogram()->GetMaximum());
558 const Double_t max2 = TMath::Max(gzd.GetHistogram()->GetMaximum(), gdzd.GetHistogram()->GetMaximum());
559 const Double_t max3 = TMath::Max(grzd.GetHistogram()->GetMaximum(), graz.GetHistogram()->GetMaximum());
560
561 const Double_t min1 = TMath::Min(gaz.GetHistogram()->GetMinimum(), gdaz.GetHistogram()->GetMinimum());
562 const Double_t min2 = TMath::Min(gzd.GetHistogram()->GetMinimum(), gdzd.GetHistogram()->GetMinimum());
563 const Double_t min3 = TMath::Min(grzd.GetHistogram()->GetMinimum(), graz.GetHistogram()->GetMinimum());
564
565 const Double_t absmax1 = TMath::Max(max1, TMath::Abs(min1));
566 const Double_t absmax2 = TMath::Max(max2, TMath::Abs(min2));
567 const Double_t absmax3 = TMath::Max(max3, TMath::Abs(min3));
568
569 gaz.SetMaximum(absmax1);
570 gzd.SetMaximum(absmax2);
571 gdaz.SetMaximum(absmax1);
572 gdzd.SetMaximum(absmax2);
573 graz.SetMaximum(absmax3);
574 grzd.SetMaximum(absmax3);
575 gaz.SetMinimum(-absmax1);
576 gzd.SetMinimum(-absmax2);
577 gdaz.SetMinimum(-absmax1);
578 gdzd.SetMinimum(-absmax2);
579
580 c1=new TCanvas("CanvGraphs", "Graphs");
581 c1->Divide(2,3,0,0);
582
583 TLine line;
584 line.SetLineColor(kGreen);
585
586 c1->cd(1);
587 gPad->SetBorderMode(0);
588 gPad->SetGridx();
589 gPad->SetGridy();
590 TGraph *g=(TGraph*)gaz.DrawClone("A*");
591 g->SetBit(kCanDelete);
592 g->GetHistogram()->SetXTitle("Az [\\circ]");
593 g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]");
594 line.DrawLine(g->GetXaxis()->GetXmin(), 360./16384, g->GetXaxis()->GetXmax(), 360./16384);
595 line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
596
597 c1->cd(2);
598 gPad->SetBorderMode(0);
599 gPad->SetGridx();
600 gPad->SetGridy();
601 g=(TGraph*)gdaz.DrawClone("A*");
602 g->SetBit(kCanDelete);
603 g->GetHistogram()->SetXTitle("Zd [\\circ]");
604 g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]");
605 line.DrawLine(g->GetXaxis()->GetXmin(), 360./16384, g->GetXaxis()->GetXmax(), 360./16384);
606 line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
607 cout << "Mean dAz: " << g->GetMean(2) << " \xb1 " << g->GetRMS(2) << endl;
608
609 c1->cd(3);
610 gPad->SetBorderMode(0);
611 gPad->SetGridx();
612 gPad->SetGridy();
613 g=(TGraph*)gdzd.DrawClone("A*");
614 g->SetBit(kCanDelete);
615 g->GetHistogram()->SetXTitle("Az [\\circ]");
616 g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]");
617 line.DrawLine(g->GetXaxis()->GetXmin(), 360./16384, g->GetXaxis()->GetXmax(), 360./16384);
618 line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
619 cout << "Mean dZd: " << g->GetMean(2) << " \xb1 " << g->GetRMS(2) << endl;
620 cout << endl;
621
622 c1->cd(4);
623 gPad->SetBorderMode(0);
624 gPad->SetGridx();
625 gPad->SetGridy();
626 g=(TGraph*)gzd.DrawClone("A*");
627 g->SetBit(kCanDelete);
628 g->GetHistogram()->SetXTitle("Zd [\\circ]");
629 g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]");
630 line.DrawLine(g->GetXaxis()->GetXmin(), 360./16384, g->GetXaxis()->GetXmax(), 360./16384);
631 line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
632
633 c1->cd(5);
634 gPad->SetBorderMode(0);
635 gPad->SetGridx();
636 gPad->SetGridy();
637 g=(TGraph*)graz.DrawClone("A*");
638 g->SetBit(kCanDelete);
639 g->GetHistogram()->SetXTitle("Az [\\circ]");
640 g->GetHistogram()->SetYTitle("\\Delta [\\circ]");
641 line.DrawLine(g->GetXaxis()->GetXmin(), 360./16384, g->GetXaxis()->GetXmax(), 360./16384);
642
643 c1->cd(6);
644 gPad->SetBorderMode(0);
645 gPad->SetGridx();
646 gPad->SetGridy();
647 g=(TGraph*)grzd.DrawClone("A*");
648 g->SetBit(kCanDelete);
649 g->GetHistogram()->SetXTitle("Zd [\\circ]");
650 g->GetHistogram()->SetYTitle("\\Delta [\\circ]");
651 line.DrawLine(g->GetXaxis()->GetXmin(), 360./16384, g->GetXaxis()->GetXmax(), 360./16384);
652
653
654
655 //
656 // Print out the residual before and after correction in several
657 // units
658 //
659 cout << fCoordinates.GetSize() << " data sets." << endl << endl;
660 cout << "Total Spread of Residual:" << endl;
661 cout << "-------------------------" << endl;
662 cout << "before: " << hres1.GetMean() << " \xb1 " << hres1.GetRMS() << " deg " << endl;
663 cout << "after: " << hres2.GetMean() << " \xb1 " << hres2.GetRMS() << " deg " << endl;
664 cout << endl;
665 cout << "before: " << Form("%.1f", hres1.GetMean()*60) << " \xb1 " << Form("%.1f", hres1.GetRMS()*60) << " arcmin" << endl;
666 cout << "after: " << Form("%.1f", hres2.GetMean()*60) << " \xb1 " << Form("%.1f", hres2.GetRMS()*60) << " arcmin" << endl;
667 cout << endl;
668 cout << "before: " << Form("%.1f", hres1.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres1.GetRMS()*60*60/23.4) << " pix" << endl;
669 cout << "after: " << Form("%.1f", hres2.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres2.GetRMS()*60*60/23.4) << " pix" << endl;
670 cout << "after: " << Form("%.1f", hres2.GetMean()*16384/360) << " \xb1 " << Form("%.1f", hres2.GetRMS()*16384/360) << " SE" << endl;
671 cout << endl;
672 cout << "backw: " << Form("%.1f", hres3.GetMean()*60) << " \xb1 " << Form("%.1f", hres3.GetRMS()*60) << " arcmin" << endl;
673 cout << endl; // ±
674
675
676 before = hres1.GetMean()*16384/360;
677 after = hres2.GetMean()*16384/360;
678 backw = hres3.GetMean()*16384/360;
679
680
681 gStyle->SetOptStat(1110);
682 gStyle->SetStatFormat("6.2g");
683
684
685 c1=new TCanvas("CanvResiduals", "Residuals", 800, 800);
686 c1->Divide(2, 2, 0, 0);
687
688 c1->cd(2);
689 gPad->SetBorderMode(0);
690 hres1.SetLineColor(kRed);
691 hres1.DrawCopy();
692
693 line.DrawLine(360./16384, 0, 360./16384, hres1.GetMaximum());
694
695 c1->cd(4);
696 gPad->SetBorderMode(0);
697 hres2.SetLineColor(kBlue);
698 TH1 *h=hres2.DrawCopy();
699 TF1 f("mygaus", "(gaus)", 0, 1);
700 f.SetLineColor(kMagenta/*6*/);
701 f.SetLineWidth(1);
702 f.SetParameter(0, h->GetBinContent(1));
703 f.FixParameter(1, 0);
704 f.SetParameter(2, h->GetRMS());
705 h->Fit("mygaus", "QR");
706 hres3.SetLineColor(kCyan);
707 hres3.SetLineStyle(kDashed);
708 hres3.DrawCopy("same");
709 cout << "Gaus-Fit Sigma: " << f.GetParameter(2) << "\xb0" << endl;
710 cout << "Fit-Probability: " << f.GetProb()*100 << "%" << endl;
711 cout << " Chi^2/NDF: " << f.GetChisquare() << "/" << f.GetNDF() << " = " << f.GetChisquare()/f.GetNDF() << endl;
712 line.DrawLine(360./16384, 0, 360./16384, h->GetMaximum());
713
714 c1->cd(1);
715 gPad->SetBorderMode(0);
716 gPad->SetTheta(90);
717 gPad->SetPhi(90);
718 TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360, 10, 0, 90);
719 h2res1.SetBit(TH1::kNoStats);
720 h2res1.DrawCopy("surf1pol");
721 gPad->Modified();
722 gPad->Update();
723 for (int i=0; i<fOriginal.GetSize(); i++)
724 DrawSet(gPad, *(Set*)fOriginal.At(i));//, 10./hres1.GetMean());
725
726 c1->cd(3);
727 gPad->SetBorderMode(0);
728 gPad->SetTheta(90);
729 gPad->SetPhi(90);
730 h2res1.SetTitle(" Arb. Residuals after correction (scaled) ");
731 h2res1.DrawCopy("surf1pol");
732 gPad->Modified();
733 gPad->Update();
734// for (int i=0; i<fCoordinates.GetSize(); i++)
735// DrawSet(gPad, *(Set*)fCoordinates.At(i), 10./hres2.GetMean(), par[0]);
736
737 RaiseWindow();
738 }
739
740 void LoadStars()
741 {
742 const Int_t size = fOriginal.GetSize();
743
744 ifstream fin("tpoint.txt");
745
746 while (fin && fin.get()!='\n');
747 while (fin && fin.get()!='\n');
748 while (fin && fin.get()!='\n');
749 if (!fin)
750 {
751 cout << "File 'tpoint.txt' not found!" << endl;
752 return;
753 }
754
755 while (1)
756 {
757 Set set;
758
759 fin >> set; // Read data from file [deg], it is stored in [rad]
760 if (!fin)
761 break;
762
763 fOriginal.Add(new Set(set));
764 }
765
766 cout << "Found " << fOriginal.GetSize()-size << " sets of coordinates ";
767 cout << "(Total=" << fOriginal.GetSize() << ")" << endl;
768 }
769
770 Bool_t ProcessMessage(Long_t msg, Long_t mp1, Long_t mp2)
771 {
772 // cout << "Msg: " << hex << GET_MSG(msg) << endl;
773 // cout << "SubMsg: " << hex << GET_SUBMSG(msg) << dec << endl;
774 switch (GET_MSG(msg))
775 {
776 case kC_COMMAND:
777 switch (GET_SUBMSG(msg))
778 {
779 case kCM_BUTTON:
780 switch (mp1)
781 {
782 case kTbFit:
783 {
784 Double_t before=0;
785 Double_t after=0;
786 Double_t backw=0;
787 Fit(before, after, backw);
788 DisplayBending();
789 DisplayResult(before, after, backw);
790 }
791 return kTRUE;
792 case kTbLoad:
793 fBending.Load("bending_new.txt");
794 DisplayBending();
795 return kTRUE;
796 case kTbSave:
797 fBending.Save("bending_new.txt");
798 return kTRUE;
799 case kTbLoadStars:
800 LoadStars();
801 DisplayData();
802 return kTRUE;
803 case kTbReset:
804 fBending.Clear();
805 DisplayBending();
806 return kTRUE;
807 case kTbResetStars:
808 fOriginal.Delete();
809 DisplayData();
810 return kTRUE;
811 }
812 return kTRUE;
813 }
814 return kTRUE;
815 }
816 return kTRUE;
817 }
818
819 void AddTextButton(TGCompositeFrame *f, TString txt, Int_t id=-1, TGLayoutHints *h=0)
820 {
821 TGButton *but = new TGTextButton(f, txt, id);
822 but->Associate(this);
823 f->AddFrame(but, h);
824 fList->Add(but);
825
826 }
827
828 void AddCheckButton(TGCompositeFrame *f, TString txt, Int_t id=-1, TGLayoutHints *h=0)
829 {
830 TGButton *but = new TGCheckButton(f, txt, id);
831 but->Associate(this);
832 f->AddFrame(but, h);
833 fList->Add(but);
834 }
835
836 TGLabel *AddLabel(TGCompositeFrame *f, TString txt, TGLayoutHints *h=0)
837 {
838 TGLabel *l = new TGLabel(f, txt/*, TGLabel::GetDefaultGC()(), fFont*/);
839 f->AddFrame(l, h);
840 fList->Add(l);
841 fLabel.Add(l);
842 return l;
843 }
844
845 void DisplayBending()
846 {
847 TArrayD par, err;
848 fBending.GetParameters(par);
849 fBending.GetError(err);
850
851 TGLabel *l;
852
853 for (int i=0; i<MBending::GetNumPar(); i++)
854 {
855 l = (TGLabel*)fLabel.At(i);
856 l->SetText(Form("%.4f\xb0", par[i]));
857
858 l = (TGLabel*)fLabel.At(MBending::GetNumPar()+i);
859 l->SetText(Form("\xb1 %8.4f\xb0", err[i]));
860 }
861 }
862
863 void DisplayData()
864 {
865 TGLabel *l = (TGLabel*)fLabel.At(3*MBending::GetNumPar());
866 l->SetText(Form("%d data sets loaded.", fOriginal.GetSize()));
867 }
868
869 void DisplayResult(Double_t before, Double_t after, Double_t backw)
870 {
871 TGLabel *l1 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+1);
872 l1->SetText(Form("Before: %.1f +- %.1f SE", before));
873
874 TGLabel *l2 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+2);
875 l2->SetText(Form("After: %.1f +- %.1f SE", after));
876
877 TGLabel *l3 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+3);
878 l3->SetText(Form("Backw: %.1f +- %.1f SE", backw));
879 }
880
881public:
882 ~MFit()
883 {
884 if (fFont)
885 gVirtualX->DeleteFont(fFont);
886 }
887 MFit() : TGMainFrame(gClient->GetRoot(), 740, 370, kHorizontalFrame)
888 {
889 fCoordinates.SetOwner();
890 fOriginal.SetOwner();
891
892 fList = new MGList;
893 fList->SetOwner();
894
895 fFont = gVirtualX->LoadQueryFont("7x13bold");
896
897 TGLayoutHints *hints0 = new TGLayoutHints(kLHintsExpandY, 7, 5, 5, 6);
898 TGLayoutHints *hints1 = new TGLayoutHints(kLHintsExpandX|kLHintsExpandY, 5, 7, 5, 6);
899 TGLayoutHints *hints2 = new TGLayoutHints(kLHintsCenterX|kLHintsCenterY, 5, 5, 5, 5);
900 fList->Add(hints0);
901 fList->Add(hints1);
902 fList->Add(hints2);
903
904 TGGroupFrame *grp1 = new TGGroupFrame(this, "Control", kVerticalFrame);
905 AddFrame(grp1, hints0);
906 fList->Add(grp1);
907
908 TGGroupFrame *grp2 = new TGGroupFrame(this, "Parameters", kHorizontalFrame);
909 AddFrame(grp2, hints1);
910 fList->Add(grp2);
911
912
913
914 TGLayoutHints *hints4 = new TGLayoutHints(kLHintsExpandX, 5, 5, 5);
915 TGLayoutHints *hints5 = new TGLayoutHints(kLHintsExpandX, 5, 5, 15);
916 AddTextButton(grp1, "Load Pointing Model", kTbLoad, hints5);
917 AddTextButton(grp1, "Save Pointing Model", kTbSave, hints4);
918 AddTextButton(grp1, "Fit Parameters", kTbFit, hints5);
919 AddTextButton(grp1, "Reset Parameters", kTbReset, hints4);
920 AddTextButton(grp1, "Load Stars", kTbLoadStars, hints5);
921 AddTextButton(grp1, "Reset Stars", kTbResetStars, hints4);
922 fList->Add(hints4);
923 fList->Add(hints5);
924
925
926 TGHorizontalFrame *comp = new TGHorizontalFrame(grp2, 1, 1);
927 grp2->AddFrame(comp);
928 fList->Add(comp);
929
930 TGLayoutHints *hints3 = new TGLayoutHints(kLHintsLeft|kLHintsTop, 0, 20, 5, 0);
931 fList->Add(hints3);
932
933 TGLabel *l;
934
935 TGVerticalFrame *vframe = new TGVerticalFrame(comp, 1, 1);
936 comp->AddFrame(vframe, hints3);
937 fList->Add(vframe);
938
939 for (int i=0; i<MBending::GetNumPar(); i++)
940 AddCheckButton(vframe, fBending.GetName(i), i);
941
942 vframe = new TGVerticalFrame(comp, 1, 1);
943 comp->AddFrame(vframe, hints3);
944 fList->Add(vframe);
945
946 l = new TGLabel(vframe, "+000.0000");
947 l->SetTextJustify(kTextRight);
948 fList->Add(l);
949 fLabel.Add(l);
950
951 TGButton *but = (TGButton*)fList->FindWidget(0);
952
953 TGLayoutHints *h = new TGLayoutHints(kLHintsCenterY, 0, 0, but->GetHeight()-l->GetHeight());
954 fList->Add(h);
955
956 vframe->AddFrame(l,h);
957
958 for (int i=1; i<MBending::GetNumPar(); i++)
959 AddLabel(vframe, "+000.0000", h)->SetTextJustify(kTextRight);
960
961 vframe = new TGVerticalFrame(comp, 1, 1);
962 comp->AddFrame(vframe, hints3);
963 fList->Add(vframe);
964
965 for (int i=0; i<MBending::GetNumPar(); i++)
966 AddLabel(vframe, "\xb1 00.0000\xb0", h)->SetTextJustify(kTextRight);
967
968 vframe = new TGVerticalFrame(comp, 1, 1);
969 comp->AddFrame(vframe, hints3);
970 fList->Add(vframe);
971
972 for (int i=0; i<MBending::GetNumPar(); i++)
973 AddLabel(vframe, fBending.GetDescription(i), h);
974
975 l = new TGLabel(grp1, "0000000 Data Sets loaded.");
976 grp1->AddFrame(l, hints5);
977 fList->Add(l);
978 fLabel.Add(l);
979
980 l = new TGLabel(grp1, "");
981 l->SetTextJustify(kTextLeft);
982 grp1->AddFrame(l, hints5);
983 fList->Add(l);
984 fLabel.Add(l);
985
986 l = new TGLabel(grp1, "");
987 l->SetTextJustify(kTextLeft);
988 grp1->AddFrame(l, hints5);
989 fList->Add(l);
990 fLabel.Add(l);
991
992 l = new TGLabel(grp1, "");
993 l->SetTextJustify(kTextLeft);
994 grp1->AddFrame(l, hints5);
995 fList->Add(l);
996 fLabel.Add(l);
997
998
999 ((TGCheckButton*)fList->FindWidget(0))->SetState(kButtonDown);
1000 ((TGCheckButton*)fList->FindWidget(1))->SetState(kButtonDown);
1001 ((TGCheckButton*)fList->FindWidget(3))->SetState(kButtonDown);
1002 ((TGCheckButton*)fList->FindWidget(4))->SetState(kButtonDown);
1003 ((TGCheckButton*)fList->FindWidget(5))->SetState(kButtonDown);
1004 ((TGCheckButton*)fList->FindWidget(6))->SetState(kButtonDown);
1005 ((TGCheckButton*)fList->FindWidget(7))->SetState(kButtonDown);
1006
1007 SetWindowName("TPoint Fitting Window");
1008 SetIconName("TPoint++");
1009
1010 Layout();
1011
1012 MapSubwindows();
1013 MapWindow();
1014
1015 DisplayBending();
1016 DisplayData();
1017 }
1018 ClassDef(MFit, 0)
1019};
1020
1021ClassImp(MFit);
1022
1023void gui()
1024{
1025 gErrorIgnoreLevel = kError;
1026 new MFit;
1027 // TF1 f1("f1", "[0]/cos((90-x)*3.1415/180)", 0, 90)
1028}
Note: See TracBrowser for help on using the repository browser.