source: trunk/MagicSoft/Cosy/tpoint/tpointfit.C@ 2273

Last change on this file since 2273 was 1819, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 12.9 KB
Line 
1#include <fstream.h>
2#include <iostream.h>
3#include <iomanip.h>
4
5#include <TF1.h>
6#include <TH1.h>
7#include <TH2.h>
8#include <TGraph.h>
9
10#include <TList.h>
11#include <TMinuit.h>
12
13#include <TView.h>
14#include <TLine.h>
15#include <TMarker.h>
16#include <TCanvas.h>
17
18#include "coord.h"
19#include "MBending.h"
20
21//
22// Sekans = 1/cos
23//
24
25class Set : public TObject
26{
27 friend ifstream &operator>>(ifstream &fin, Set &set);
28private:
29 Double_t fStarAz;
30 Double_t fStarEl;
31
32 Double_t fRawAz;
33 Double_t fRawEl;
34
35public:
36 Set(Double_t sel=0, Double_t saz=0, Double_t rel=0, Double_t raz=0) :
37 fStarAz(saz*TMath::Pi()/180),
38 fStarEl(sel*TMath::Pi()/180),
39 fRawAz(raz*TMath::Pi()/180),
40 fRawEl(rel*TMath::Pi()/180)
41 {
42 }
43
44 Double_t GetResidual() const
45 {
46 Double_t del = fRawEl-fStarEl;
47 Double_t daz = fRawAz-fStarAz;
48 Double_t dphi2 = daz/2.;
49 Double_t cos2 = cos(dphi2)*cos(dphi2);
50 Double_t sin2 = sin(dphi2)*sin(dphi2);
51 Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
52
53 Double_t dist = acos(d);
54
55 return dist * 180 / TMath::Pi();
56 }
57
58 void operator=(Set &set)
59 {
60 fStarAz = set.fStarAz;
61 fStarEl = set.fStarEl;
62 fRawAz = set.fRawAz;
63 fRawEl = set.fRawEl;
64 }
65
66 Double_t GetDEl() const { return (fRawEl-fStarEl)*180/TMath::Pi(); }
67 Double_t GetDZd() const { return -GetDEl(); }
68 Double_t GetDAz() const { return (fRawAz-fStarAz)*180/TMath::Pi(); }
69 Double_t GetStarEl() const { return fStarEl*180/TMath::Pi(); }
70 Double_t GetStarZd() const { return 90.-fStarEl*180/TMath::Pi(); }
71 Double_t GetStarAz() const { return fStarAz*180/TMath::Pi(); }
72 Double_t GetRawEl() const { return fRawEl*180/TMath::Pi(); }
73 Double_t GetRawAz() const { return fRawAz*180/TMath::Pi(); }
74 Double_t GetRawZd() const { return 90.-fRawEl*180/TMath::Pi(); }
75
76 ZdAz GetStarZdAz() const { return ZdAz(TMath::Pi()/2-fStarEl, fStarAz); }
77 AltAz GetStarAltAz() const { return AltAz(fStarEl, fStarAz); }
78
79 ZdAz GetRawZdAz() const { return ZdAz(TMath::Pi()/2-fRawEl, fRawAz); }
80 AltAz GetRawAltAz() const { return AltAz(fRawEl, fRawAz); }
81
82 void AdjustEl(Double_t del) { fStarEl += del*TMath::Pi()/180; }
83 void AdjustAz(Double_t daz) { fStarAz += daz*TMath::Pi()/180; }
84
85 void Adjust(const MBending &bend)
86 {
87 AltAz p = bend(GetStarAltAz());
88 fStarEl = p.Alt();
89 fStarAz = p.Az();
90 }
91 void AdjustBack(const MBending &bend)
92 {
93 AltAz p = bend.CorrectBack(GetRawAltAz());
94 fRawEl = p.Alt();
95 fRawAz = p.Az();
96 }
97};
98
99ifstream &operator>>(ifstream &fin, Set &set)
100{
101 Double_t v[4];
102 fin >> v[0];
103 fin >> v[1];
104 fin >> v[2];
105 fin >> v[3];
106
107 set.fStarAz = v[0]*TMath::Pi()/180;
108 set.fStarEl = v[1]*TMath::Pi()/180;
109
110 set.fRawAz = v[2]*TMath::Pi()/180;
111 set.fRawEl = v[3]*TMath::Pi()/180;
112
113 return fin;
114}
115
116TList gCoordinates;
117
118void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
119{
120 f = 0;
121
122 MBending bend;
123 bend.SetParameters(par); // Set Parameters [deg] to MBending
124
125 for (int i=0; i<gCoordinates.GetSize(); i++)
126 {
127 Set set = *(Set*)gCoordinates.At(i);
128
129 set.Adjust(bend);
130
131 Double_t err = 1;
132 Double_t res = set.GetResidual()/err;
133
134 f += res*res;
135 }
136
137 f /= gCoordinates.GetSize()*gCoordinates.GetSize();
138}
139
140void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
141{
142 TView *view = pad->GetView();
143
144 if (!view)
145 {
146 cout << "No View!" << endl;
147 return;
148 }
149
150 TMarker mark0;
151 TMarker mark1;
152 mark0.SetMarkerStyle(kStar);
153 mark1.SetMarkerStyle(kStar);
154 mark1.SetMarkerColor(kRed);
155
156 r0 /= 90;
157 r1 /= 90;
158 phi0 *= TMath::Pi()/180;
159 phi1 *= TMath::Pi()/180;
160
161 Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
162 Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
163
164 Double_t y0[3], y1[3];
165
166 view->WCtoNDC(x0, y0);
167 view->WCtoNDC(x1, y1);
168
169 mark0.DrawMarker(y0[0], y0[1]);
170 mark1.DrawMarker(y1[0], y1[1]);
171}
172
173void DrawPolLine(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
174{
175 TView *view = pad->GetView();
176
177 if (!view)
178 {
179 cout << "No View!" << endl;
180 return;
181 }
182 /*
183 if (r0<0)
184 {
185 r0 = -r0;
186 phi0 += 180;
187 }
188 if (r1<0)
189 {
190 r1 = -r1;
191 phi1 += 180;
192 }
193
194 phi0 = fmod(phi0+360, 360);
195 phi1 = fmod(phi1+360, 360);
196
197 if (phi1-phi0<-180)
198 phi1+=360;
199 */
200 TLine line;
201 line.SetLineWidth(2);
202 line.SetLineColor(kBlue);
203
204 Double_t p0 = phi0<phi1?phi0:phi1;
205 Double_t p1 = phi0<phi1?phi1:phi0;
206
207 if (phi0>phi1)
208 {
209 Double_t d = r1;
210 r1 = r0;
211 r0 = d;
212 }
213
214 r0 /= 90;
215 r1 /= 90;
216
217 Double_t dr = r1-r0;
218 Double_t dp = p1-p0;
219
220 Double_t x0[3] = { r0*cos(p0*TMath::Pi()/180), r0*sin(p0*TMath::Pi()/180), 0};
221
222 for (double i=p0+10; i<p1+10; i+=10)
223 {
224 if (i>p1)
225 i=p1;
226
227 Double_t r = dr/dp*(i-p0)+r0;
228 Double_t p = TMath::Pi()/180*i;
229
230 Double_t x1[3] = { r*cos(p), r*sin(p), 0};
231
232 Double_t y0[3], y1[3];
233
234 view->WCtoNDC(x0, y0);
235 view->WCtoNDC(x1, y1);
236
237 line.DrawLine(y0[0], y0[1], y1[0], y1[1]);
238
239 x0[0] = x1[0];
240 x0[1] = x1[1];
241 }
242}
243
244void DrawSet(TVirtualPad *pad, Set &set, Float_t scale=-1)
245{
246 Double_t r0 = set.GetRawZd();
247 Double_t phi0 = set.GetRawAz();
248 Double_t r1 = set.GetStarZd();
249 Double_t phi1 = set.GetStarAz();
250
251 if (r0<0)
252 {
253 r0 = -r0;
254 phi0 += 180;
255 }
256 if (r1<0)
257 {
258 r1 = -r1;
259 phi1 += 180;
260 }
261
262 phi0 = fmod(phi0+360, 360);
263 phi1 = fmod(phi1+360, 360);
264
265 if (phi1-phi0<-180)
266 phi1+=360;
267
268 if (scale<0 || scale>1000)
269 scale = -1;
270
271 if (scale>0)
272 {
273 Double_t d = r1-r0;
274 r0 += scale*d;
275 r1 -= scale*d;
276 d = phi1-phi0;
277 phi0 += scale*d;
278 phi1 -= scale*d;
279
280 DrawPolLine(pad, r0, phi0, r1, phi1);
281 DrawMarker(pad, r0, phi0, r1, phi1);
282 }
283 else
284 DrawMarker(pad, r1, phi1, 0 ,0);
285}
286
287void tpointfit()
288{
289 cout << "Starting..." << endl;
290
291 gCoordinates.SetOwner();
292
293 TH1F hres1("Res1", " Residuals before correction ", 10, 0, 4);
294 TH1F hres2("Res2", " Residuals after correction ", 40, 0, 4);
295
296 hres1.SetXTitle("\\Delta [\\circ]");
297 hres1.SetYTitle("Counts");
298
299 hres2.SetXTitle("\\Delta [\\circ]");
300 hres2.SetYTitle("Counts");
301
302 TGraph gdaz;
303 TGraph gdzd;
304 TGraph gaz;
305 TGraph gzd;
306
307 gdaz.SetTitle(" \\Delta Az vs. Zd ");
308 gdzd.SetTitle(" \\Delta Zd vs. Az ");
309
310 gaz.SetTitle(" \\Delta Az vs. Az ");
311 gzd.SetTitle(" \\Delta Zd vs. Zd ");
312
313 ifstream fin("/home/tbretz/pc4/results/tpoint3.txt");
314
315 while (fin && fin.get()!='\n');
316 while (fin && fin.get()!='\n');
317 while (fin && fin.get()!='\n');
318 if (!fin)
319 {
320 cout << "File not found!" << endl;
321 return;
322 }
323
324 while (1)
325 {
326 Set set;
327
328 fin >> set; // Read data from file [deg], it is stored in [rad]
329
330 if (!fin)
331 break;
332
333 gCoordinates.Add(new Set(set));
334 }
335
336 cout << "Found " << gCoordinates.GetSize() << " sets of coordinates." << endl;
337
338 TMinuit minuit(16); //initialize TMinuit with a maximum of 5 params
339 minuit.SetPrintLevel(-1);
340 minuit.SetFCN(fcn);
341
342 MBending bending;
343 bending.SetMinuitParameters(minuit, 16); // Init Parameters [deg]
344
345 //minuit.FixParameter(0); //(1) IA
346 //minuit.FixParameter(1); //(1) IE
347 minuit.FixParameter(2); //(2) NPAE
348 //minuit.FixParameter(3); //(1) CA
349 minuit.FixParameter(4); //(2) AN
350 //minuit.FixParameter(5); //(1) AW
351
352 minuit.FixParameter(6); // NRX
353 minuit.FixParameter(7); // NRY
354 minuit.FixParameter(8); // CRX
355 minuit.FixParameter(9); // CRY
356 minuit.FixParameter(10); // ECES
357 minuit.FixParameter(11); //(2) ACES
358 minuit.FixParameter(12); // ECEC
359 minuit.FixParameter(13); //(2) ACEC
360 minuit.FixParameter(14); // MAGIC1
361 minuit.FixParameter(15); // MAGIC2
362
363
364 //minuit.Command("SHOW PARAMETERS");
365 //minuit.Command("SHOW LIMITS");
366 cout << endl;
367
368 Int_t ierflg = 0;
369 ierflg = minuit.Migrad();
370 cout << "Migrad returns " << ierflg << endl;
371 // minuit.Release(2);
372 ierflg = minuit.Migrad();
373 cout << "Migrad returns " << ierflg << endl << endl;
374
375 //
376 // Get Fit Results
377 //
378 bending.GetMinuitParameters(minuit);
379 bending.PrintMinuitParameters(minuit);
380 bending.Save("bending_magic.txt");
381
382
383 //
384 // Make a copy of all list entries
385 //
386 TList list;
387 list.SetOwner();
388 for (int i=0; i<gCoordinates.GetSize(); i++)
389 list.Add(new Set(*(Set*)gCoordinates.At(i)));
390
391 //
392 // Correct for Offsets only
393 //
394 double par[16];
395 bending.GetParameters(par);
396 for (int i=2; i<16; i++) par[i]=0;
397 MBending b2;
398 b2.SetParameters(par);
399
400 //
401 // Calculate correction and residuals
402 //
403 for (int i=0; i<gCoordinates.GetSize(); i++)
404 {
405 Set &set0 = *(Set*)gCoordinates.At(i);
406
407 ZdAz za = set0.GetStarZdAz()*kRad2Deg;
408
409 //
410 // Correct for offsets only
411 //
412 Set set1(set0);
413 set1.Adjust(b2);
414
415 hres1.Fill(set1.GetResidual());
416
417 set0.Adjust(bending);
418 hres2.Fill(set0.GetResidual());
419
420 Double_t dz = fmod(set0.GetDAz()+720, 360);
421 if (dz>180)
422 dz -= 360;
423
424 static int j=0;
425 gdzd.SetPoint(j, za.Az(), set0.GetDZd());
426 gaz.SetPoint( j, za.Az(), dz);
427 gdaz.SetPoint(j, za.Zd(), dz);
428 gzd.SetPoint( j, za.Zd(), set0.GetDZd());
429 j++;
430 }
431
432 //
433 // Check for overflows
434 //
435 const Stat_t ov = hres2.GetBinContent(hres2.GetNbinsX()+1);
436 if (ov>0)
437 cout << "WARNING: " << ov << " overflows in residuals." << endl;
438
439 //
440 // Print all data sets for which the backward correction is
441 // twice times worse than the residual gotten from the
442 // bending correction itself
443 //
444 cout << endl;
445 cout << "Checking backward correction (raw-->star):" << endl;
446 for (int i=0; i<gCoordinates.GetSize(); i++)
447 {
448 Set set0(*(Set*)list.At(i));
449 Set set1(*(Set*)list.At(i));
450
451 set0.AdjustBack(bending);
452 set1.Adjust(bending);
453
454 Double_t res0 = set0.GetResidual();
455 Double_t res1 = set1.GetResidual();
456
457 Double_t diff = fabs(res0-res1);
458
459 if (diff<hres2.GetMean()*0.66)
460 continue;
461
462 cout << "DBack: " << setw(7) << set0.GetStarZd() << " " << setw(8) << set0.GetStarAz() << ": ";
463 cout << "ResB="<< setw(7) << res0*60 << " ResF=" << setw(7) << res1*60 << " |ResB-ResF|=" << setw(7) << diff*60 << " arcmin" << endl;
464 }
465 cout << "OK." << endl;
466
467 //
468 // Print out the residual before and after correction in several
469 // units
470 //
471 cout << endl;
472 cout << gCoordinates.GetSize() << " data sets." << endl << endl;
473 cout << "Spead before: " << hres1.GetMean() << "deg +- " << hres1.GetRMS() << "deg" << endl;
474 cout << "Spead after: " << hres2.GetMean() << "deg +- " << hres2.GetRMS() << "deg" << endl;
475 cout << "Spead before: " << (int)(hres1.GetMean()*60+.5) << "' +- " << (int)(hres1.GetRMS()*60+.5) << "'" << endl;
476 cout << "Spead before: " << (int)(hres1.GetMean()*60*60/23.4+.5) << "pix +- " << (int)(hres1.GetRMS()*60*60/23.4+.5) << "pix" << endl;
477 cout << "Spead after: " << (int)(hres2.GetMean()*60+.5) << "' +- " << (int)(hres2.GetRMS()*60+.5) << "'" << endl;
478 cout << "Spead after: " << (int)(hres2.GetMean()*60*60/23.4+.5) << "pix +- " << (int)(hres2.GetRMS()*60*60/23.4+.5) << "pix" << endl;
479 cout << "Spead after: " << (int)(hres2.GetMean()*16384/360+.5) << "SE +- " << (int)(hres2.GetRMS()*16384/360+.5) << "SE" << endl;
480 cout << endl;
481
482
483
484 TCanvas *c1=new TCanvas("c1", "c");
485 c1->Divide(2,2);
486
487 c1->cd(2);
488 gdaz.DrawClone("A*");
489
490 c1->cd(1);
491 gaz.DrawClone("A*");
492
493 c1->cd(3);
494 gdzd.DrawClone("A*");
495
496 c1->cd(4);
497 gzd.DrawClone("A*");
498
499
500
501 c1=new TCanvas("c2", "c", 800, 800);
502 c1->Divide(2,2);
503
504 c1->cd(2);
505 hres1.SetLineColor(kRed);
506 hres1.DrawCopy();
507
508 c1->cd(4);
509 hres2.SetLineColor(kBlue);
510 hres2.DrawCopy();
511
512 c1->cd(1);
513 gPad->SetTheta(90);
514 gPad->SetPhi(90);
515 TH2F h2res1("Res2D1", " Arb. Residuals before correction (scaled) ", 32, 0, 360, 10, 0, 90);
516 h2res1.DrawCopy("surf1pol");
517 gPad->Modified();
518 gPad->Update();
519 for (int i=0; i<gCoordinates.GetSize(); i++)
520 DrawSet(gPad, *(Set*)list.At(i));//, 10./hres1.GetMean());
521
522 c1->cd(3);
523 gPad->SetTheta(90);
524 gPad->SetPhi(90);
525 h2res1.SetTitle(" Arb. Residuals after correction (scaled) ");
526 h2res1.DrawCopy("surf1pol");
527 gPad->Modified();
528 gPad->Update();
529 for (int i=0; i<gCoordinates.GetSize(); i++)
530 DrawSet(gPad, *(Set*)gCoordinates.At(i), 10./hres2.GetMean());
531
532}
Note: See TracBrowser for help on using the repository browser.