source: trunk/MagicSoft/Cosy/main/MBending.cc@ 2164

Last change on this file since 2164 was 1810, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 15.4 KB
Line 
1#include "MBending.h"
2
3#include <fstream.h>
4#include <iomanip.h>
5
6#include <TMinuit.h>
7
8#include "timer.h"
9
10ClassImp(MBending);
11
12void MBending::Reset()
13{
14 Clear();
15}
16
17void MBending::Load(const char *name)
18{
19 /*
20 ! MMT 1987 July 8
21 ! T 36 7.3622 41.448 -0.0481
22 ! IA -37.5465 20.80602
23 ! IE -13.9180 1.25217
24 ! NPAE +7.0751 26.44763
25 ! CA -6.9149 32.05358
26 ! AN +0.5053 1.40956
27 ! AW -2.2016 1.37480
28 ! END
29 */
30
31 ifstream fin(name);
32 if (!fin)
33 {
34 cout << "Error: Cannot open file '" << name << "'" << endl;
35 return;
36 }
37
38 char c;
39 while (fin && fin.get()!='\n');
40 fin >> c;
41
42 if (c!='S' && c!='s')
43 {
44 cout << "Error: This in not a model correcting the star position (" << c << ")" << endl;
45 return;
46 }
47
48 Clear();
49
50 cout << endl;
51
52 Double_t val;
53 fin >> val;
54 cout << "Number of observed stars: " << val << endl;
55 fin >> val;
56 cout << "Sky RMS: " << val << "\"" << endl;
57 fin >> val;
58 cout << "Refraction Constant A: " << val << "\"" << endl;
59 fin >> val;
60 cout << "Refraction Constant B: " << val << "\"" << endl;
61
62 cout << endl;
63
64 cout << " & = Name Value Sigma" << endl;
65 cout << "------------------------------------------" << endl;
66
67 while (fin)
68 {
69 TString str;
70 fin >> str;
71
72 if (str=="END")
73 break;
74
75 if (str[0]=='&')
76 {
77 cout << " & ";
78 str.Remove(0);
79 }
80 else
81 cout << " ";
82
83 if (str[1]=='=')
84 {
85 cout << "= ";
86 str.Remove(0);
87 }
88 else
89 cout << " ";
90
91 fin >> val;
92 cout << str << "\t" << val << "° \t";
93 val *= kDeg2Rad;
94
95 if (str=="IA") fIa = val;
96 if (str=="IE") fIe = val;
97 if (str=="CA") fCa = val;
98 if (str=="AN") fAn = val;
99 if (str=="AW") fAw = val;
100 if (str=="NRX") fNrx = val;
101 if (str=="NRY") fNry = val;
102 if (str=="CRX") fCrx = val;
103 if (str=="CRY") fCry = val;
104 if (str=="NPAE") fNpae = val;
105 if (str=="ECES") fEces = val;
106 if (str=="ACES") fAces = val;
107 if (str=="ECEC") fEcec = val;
108 if (str=="ACEC") fAcec = val;
109 //if (str=="MAGIC1") fMagic1 = val;
110 //if (str=="MAGIC2") fMagic2 = val;
111
112 fin >> val;
113 cout << val*kDeg2Rad << endl;
114 }
115 cout << endl;
116}
117
118void MBending::Save(const char *name)
119{
120 /*
121 ! MMT 1987 July 8
122 ! T 36 7.3622 41.448 -0.0481
123 ! IA -37.5465 20.80602
124 ! IE -13.9180 1.25217
125 ! NPAE +7.0751 26.44763
126 ! CA -6.9149 32.05358
127 ! AN +0.5053 1.40956
128 ! AW -2.2016 1.37480
129 ! END
130 */
131
132 ofstream fout(name);
133 if (!fout)
134 {
135 cout << "Error: Cannot open file '" << name << "'" << endl;
136 return;
137 }
138
139 Timer t;
140 t.Now();
141
142 fout << "MAGIC1 " << t.GetTimeStr() << endl;
143 fout << "S 00 000000 000000 0000000" << endl;
144 fout << setprecision(8);
145 fout << " IA " << kRad2Deg*fIa << " -1" << endl;
146 fout << " IE " << kRad2Deg*fIe << " -1" << endl;
147 fout << " CA " << kRad2Deg*fNpae << " -1" << endl;
148 fout << " NPAE " << kRad2Deg*fCa << " -1" << endl;
149 fout << " AN " << kRad2Deg*fAn << " -1" << endl;
150 fout << " AW " << kRad2Deg*fAw << " -1" << endl;
151 fout << " NRX " << kRad2Deg*fNrx << " -1" << endl;
152 fout << " NRY " << kRad2Deg*fNry << " -1" << endl;
153 fout << " CRX " << kRad2Deg*fCrx << " -1" << endl;
154 fout << " CRY " << kRad2Deg*fCry << " -1" << endl;
155 fout << " ECES " << kRad2Deg*fEces << " -1" << endl;
156 fout << " ACES " << kRad2Deg*fAces << " -1" << endl;
157 fout << " ECEC " << kRad2Deg*fEcec << " -1" << endl;
158 fout << " ACEC " << kRad2Deg*fAcec << " -1" << endl;
159 //fout << " MAGIC1 " << -kRad2Deg*fMagic1 << " -1" << endl;
160 //fout << " MAGIC2 " << -kRad2Deg*fMagic2 << " -1" << endl;
161 fout << "END" << endl;
162}
163
164AltAz MBending::Correct(const AltAz &aa) const
165{
166 // Correct [rad]
167 // zdaz [rad]
168 AltAz p = aa;
169
170 const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
171 const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
172 p += CES;
173 p += CEC;
174
175 // const AltAz MAGIC(0, -fMagic1*tan(p.Alt())-fMagic2);
176 // p += MAGIC;
177
178 const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
179 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
180 p += CRX;
181 p += CRY;
182
183 const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
184 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
185 p += NRX;
186 p += NRY;
187
188 // const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
189 // p += MAGIC;
190
191 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
192 const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
193 p += AW;
194 p += AN;
195
196 // const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
197 // p += MAGIC;
198
199 const AltAz CA(0, -fCa/cos(p.Alt()));
200 p += CA;
201
202 const AltAz NPAE(0, -fNpae*tan(p.Alt()));
203 p += NPAE;
204
205 const AltAz I(fIe, fIa);
206 p += I;
207
208 return p;
209}
210
211AltAz MBending::CorrectBack(const AltAz &aa) const
212{
213 // Correct [rad]
214 // zdaz [rad]
215 AltAz p = aa;
216
217 const AltAz I(fIe, fIa);
218 p -= I;
219
220 const AltAz NPAE(0, -fNpae*tan(p.Alt()));
221 p -= NPAE;
222
223 const AltAz CA(0, -fCa/cos(p.Alt()));
224 p -= CA;
225
226 const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
227 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
228 p -= AN;
229 p -= AW;
230
231 // const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
232 // p -= MAGIC;
233
234 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
235 const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
236 p -= NRY;
237 p -= NRX;
238
239 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
240 const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
241 p -= CRY;
242 p -= CRX;
243
244 const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
245 const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
246 p -= CEC;
247 p -= CES;
248
249 return p;
250}
251
252ZdAz MBending::Correct(const ZdAz &zdaz) const
253{
254 // Correct [rad]
255 // zdaz [rad]
256 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
257 AltAz c = Correct(p);
258 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
259}
260
261ZdAz MBending::CorrectBack(const ZdAz &zdaz) const
262{
263 // Correct [rad]
264 // zdaz [rad]
265 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
266 AltAz c = CorrectBack(p);
267 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
268}
269
270void MBending::SetParameters(const Double_t *par, Int_t n)
271{
272 Clear();
273
274 switch (n)
275 {
276 case 16:
277 fMagic2=par[15]/kRad2Deg; // Magic User Defined
278 case 15:
279 fMagic1=par[14]/kRad2Deg; // Magic User Defined
280 case 14:
281 fAcec =par[13]/kRad2Deg; // Azimuth Centering Error (cos)
282 case 13:
283 fEcec =par[12]/kRad2Deg; // Elevation Centering Error (cos)
284 case 12:
285 fAces =par[11]/kRad2Deg; // Azimuth Centering Error (sin)
286 case 11:
287 fEces =par[10]/kRad2Deg; // Elevation Centering Error (sin)
288 case 10:
289 fCry =par[9] /kRad2Deg; // Alt/Az Coude Displacement (E-W)
290 case 9:
291 fCrx =par[8] /kRad2Deg; // Alt/Az Coude Displacement (N-S)
292 case 8:
293 fNry =par[7] /kRad2Deg; // Nasmyth rotator displacement, vertical
294 case 7:
295 fNrx =par[6] /kRad2Deg; // Nasmyth rotator displacement, horizontan
296 case 6:
297 fAw =par[5] /kRad2Deg; // Azimuth Axis Misalignment (E-W)
298 case 5:
299 fAn =par[4] /kRad2Deg; // Azimuth Axis Misalignment (N-S)
300 case 4:
301 fCa =par[3] /kRad2Deg; // Left-Right Collimation Error
302 case 3:
303 fNpae =par[2] /kRad2Deg; // Az-El Nonperpendicularity
304 case 2:
305 fIe =par[1] /kRad2Deg; // Index Error in Elevation
306 case 1:
307 fIa =par[0] /kRad2Deg; // Index Error in Azimuth
308 }
309}
310
311void MBending::GetParameters(Double_t *par, Int_t n) const
312{
313 switch (n)
314 {
315 case 16:
316 par[15]=kRad2Deg*fMagic2; //
317 case 15:
318 par[14]=kRad2Deg*fMagic1; //
319 case 14:
320 par[13]=kRad2Deg*fAcec; // Azimuth Centering Error (cos)
321 case 13:
322 par[12]=kRad2Deg*fEcec; // Elevation Centering Error (cos)
323 case 12:
324 par[11]=kRad2Deg*fAces; // Azimuth Centering Error (sin)
325 case 11:
326 par[10]=kRad2Deg*fEces; // Elevation Centering Error (sin)
327 case 10:
328 par[9] =kRad2Deg*fCry; // Alt/Az Coude Displacement (E-W)
329 case 9:
330 par[8] =kRad2Deg*fCrx; // Alt/Az Coude Displacement (N-S)
331 case 8:
332 par[7] =kRad2Deg*fNry; // Nasmyth rotator displacement, vertical
333 case 7:
334 par[6] =kRad2Deg*fNrx; // Nasmyth rotator displacement, horizontan
335 case 6:
336 par[5] =kRad2Deg*fAw; // Azimuth Axis Misalignment (E-W)
337 case 5:
338 par[4] =kRad2Deg*fAn; // Azimuth Axis Misalignment (N-S)
339 case 4:
340 par[3] =kRad2Deg*fCa; // Left-Right Collimation Error
341 case 3:
342 par[2] =kRad2Deg*fNpae; // Az-El Nonperpendicularity
343 case 2:
344 par[1] =kRad2Deg*fIe; // Index Error in Elevation
345 case 1:
346 par[0] =kRad2Deg*fIa; // Index Error in Azimuth
347 }
348}
349
350void MBending::SetMinuitParameters(TMinuit &m, Int_t n=-1) const
351{
352 if (n<0)
353 n = 16;
354
355 Int_t ierflg = 0;
356
357 switch (n)
358 {
359 case 16:
360 m.mnparm(15,"MAGIC2", fMagic2*kRad2Deg, 1, -360, 360, ierflg);
361 case 15:
362 m.mnparm(14,"MAGIC1", fMagic1*kRad2Deg, 1, -360, 360, ierflg);
363 case 14:
364 m.mnparm(13,"ACEC", fAcec*kRad2Deg, 1, -360, 360, ierflg);
365 case 13:
366 m.mnparm(12,"ECEC", fEcec*kRad2Deg, 1, -360, 360, ierflg);
367 case 12:
368 m.mnparm(11,"ACES", fAcec*kRad2Deg, 1, -360, 360, ierflg);
369 case 11:
370 m.mnparm(10,"ECES", fEcec*kRad2Deg, 1, -360, 360, ierflg);
371 case 10:
372 m.mnparm(9, "CRY", fCry*kRad2Deg, 1, -360, 360, ierflg);
373 case 9:
374 m.mnparm(8, "CRX", fCrx*kRad2Deg, 1, -360, 360, ierflg);
375 case 8:
376 m.mnparm(7, "NRY", fNry*kRad2Deg, 1, -360, 360, ierflg);
377 case 7:
378 m.mnparm(6, "NRX", fNrx*kRad2Deg, 1, -360, 360, ierflg);
379 case 6:
380 m.mnparm(5, "AW", fAw*kRad2Deg, 1, -360, 360, ierflg);
381 case 5:
382 m.mnparm(4, "AN", fAn*kRad2Deg, 1, -360, 360, ierflg);
383 case 4:
384 m.mnparm(3, "CA", fCa*kRad2Deg, 1, -360, 360, ierflg);
385 case 3:
386 m.mnparm(2, "NPAE", fNpae*kRad2Deg, 1, -360, 360, ierflg);
387 case 2:
388 m.mnparm(1, "IE", fIe*kRad2Deg, 1, -90, 90, ierflg);
389 case 1:
390 m.mnparm(0, "IA", fIa*kRad2Deg, 1, -360, 360, ierflg);
391 }
392}
393
394void MBending::GetMinuitParameters(TMinuit &m, Int_t n=-1)
395{
396 if (n<0 || n>m.GetNumPars())
397 n = m.GetNumPars();
398
399 Double_t err;
400
401 switch (n)
402 {
403 case 16:
404 m.GetParameter(15, fMagic2, err);
405 fMagic2 /= kRad2Deg;
406 case 15:
407 m.GetParameter(14, fMagic1, err);
408 fMagic1 /= kRad2Deg;
409 case 14:
410 m.GetParameter(13, fAcec, err);
411 fAcec /= kRad2Deg;
412 case 13:
413 m.GetParameter(12, fEcec, err);
414 fEcec /= kRad2Deg;
415 case 12:
416 m.GetParameter(11, fAces, err);
417 fAces /= kRad2Deg;
418 case 11:
419 m.GetParameter(10, fEces, err);
420 fEces /= kRad2Deg;
421 case 10:
422 m.GetParameter(9, fCry, err);
423 fCry /= kRad2Deg;
424 case 9:
425 m.GetParameter(8, fCrx, err);
426 fCrx /= kRad2Deg;
427 case 8:
428 m.GetParameter(7, fNry, err);
429 fNry /= kRad2Deg;
430 case 7:
431 m.GetParameter(6, fNrx, err);
432 fNrx /= kRad2Deg;
433 case 6:
434 m.GetParameter(5, fAw, err);
435 fAw /= kRad2Deg;
436 case 5:
437 m.GetParameter(4, fAn, err);
438 fAn /= kRad2Deg;
439 case 4:
440 m.GetParameter(3, fCa, err);
441 fCa /= kRad2Deg;
442 case 3:
443 m.GetParameter(2, fNpae, err);
444 fNpae /= kRad2Deg;
445 case 2:
446 m.GetParameter(1, fIe, err);
447 fIe /= kRad2Deg;
448 case 1:
449 m.GetParameter(0, fIa, err);
450 fIa /= kRad2Deg;
451 }
452}
453/*
454void FormatPar(TMinuit &m, Int_t n)
455{
456 Double_t par, err;
457 m.GetParameter(n, par, err);
458
459 int expp = (int)log10(par);
460 int expe = (int)log10(err);
461
462 if (err<2*pow(10, expe))
463 expe--;
464
465 Int_t exp = expe>expp ? expp : expe;
466
467 par = (int)(par/pow(10, exp)) * pow(10, exp);
468 err = (int)(err/pow(10, exp)) * pow(10, exp);
469
470 cout << par << " +- " << err << flush;
471}
472*/
473void MBending::PrintMinuitParameters(TMinuit &m, Int_t n=-1) const
474{
475 if (n<0)
476 n = m.GetNumPars();
477
478 cout << setprecision(3);
479
480 Double_t par, err;
481
482 switch (n)
483 {
484 case 16:
485 m.GetParameter(15, par, err);
486 cout << " 15 MAGIC2: " << setw(8) << par << " +- " << setw(4) << err << endl;
487 case 15:
488 m.GetParameter(14, par, err);
489 cout << " 14 MAGIC1: " << setw(8) << par << " +- " << setw(4) << err << endl;
490 case 14:
491 m.GetParameter(13, par, err);
492 cout << " 13 ACEC: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Centering Error (cos)" << endl;
493 case 13:
494 m.GetParameter(12, par, err);
495 cout << " 12 ECEC: " << setw(8) << par << " +- " << setw(4) << err << " Elevation Centering Error (cos)" << endl;
496 case 12:
497 m.GetParameter(11, par, err);
498 cout << " 11 ACES: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Centering Error (sin)" << endl;
499 case 11:
500 m.GetParameter(10, par, err);
501 cout << " 10 ECES: " << setw(8) << par << " +- " << setw(4) << err << " Elevation Centering Error (sin)" << endl;
502 case 10:
503 m.GetParameter(9, par, err);
504 cout << " 9 CRY: " << setw(8) << par << " +- " << setw(4) << err << " Alt/Az Coude Displacement (E-W)" << endl;
505 case 9:
506 m.GetParameter(8, par, err);
507 cout << " 8 CRX: " << setw(8) << par << " +- " << setw(4) << err << " Alt/Az Coude Displacement (N-S)" << endl;
508 case 8:
509 m.GetParameter(7, par, err);
510 cout << " 7 NRY: " << setw(8) << par << " +- " << setw(4) << err << " Nasmyth rotator displacement, vertical" << endl;
511 case 7:
512 m.GetParameter(6, par, err);
513 cout << " 6 NRX: " << setw(8) << par << " +- " << setw(4) << err << " Nasmyth rotator displacement, horizontan" << endl;
514 case 6:
515 m.GetParameter(5, par, err);
516 cout << " 5 AW: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Axis Misalignment (E-W)" << endl;
517 case 5:
518 m.GetParameter(4, par, err);
519 cout << " 4 AN: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Axis Misalignment (N-S)" << endl;
520 case 4:
521 m.GetParameter(3, par, err);
522 cout << " 3 CA: " << setw(8) << par << " +- " << setw(4) << err << " Left-Right Collimation Error" << endl;
523 case 3:
524 m.GetParameter(2, par, err);
525 cout << " 2 NPAE: " << setw(8) << par << " +- " << setw(4) << err << " Az-El Nonperpendicularity" << endl;
526 case 2:
527 m.GetParameter(1, par, err);
528 cout << " 1 IE: " << setw(8) << par << " +- " << setw(4) << err << " Index Error Elevation (Offset)" << endl;
529 case 1:
530 m.GetParameter(0, par, err);
531 cout << " 0 IA: " << setw(8) << par << " +- " << setw(4) << err << " Index Error Azimuth (Offset)" << endl;
532 }
533}
Note: See TracBrowser for help on using the repository browser.