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

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