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