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

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