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

Last change on this file since 1806 was 1806, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 14.7 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 CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
176 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
177 p += CRX;
178 p += CRY;
179
180 const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
181 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
182 p += NRX;
183 p += NRY;
184
185 const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
186 p += MAGIC;
187
188 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
189 const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
190 p += AW;
191 p += AN;
192
193// const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
194// p += MAGIC;
195
196
197 const AltAz CA(0, -fCa/cos(p.Alt()));
198 p += CA;
199
200 const AltAz NPAE(0, -fNpae*tan(p.Alt()));
201 p += NPAE;
202
203 const AltAz I(fIe, fIa);
204 p += I;
205
206 return p;
207}
208
209AltAz MBending::CorrectBack(const AltAz &aa) const
210{
211 // Correct [rad]
212 // zdaz [rad]
213 AltAz p = aa;
214
215 const AltAz I(fIe, fIa);
216 p -= I;
217
218 const AltAz NPAE(0, -fNpae*tan(p.Alt()));
219 p -= NPAE;
220
221 const AltAz CA(0, -fCa/cos(p.Alt()));
222 p -= CA;
223
224 const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
225 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
226 p -= AN;
227 p -= AW;
228
229 const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
230 p -= MAGIC;
231
232 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
233 const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
234 p -= NRY;
235 p -= NRX;
236
237 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
238 const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
239 p -= CRY;
240 p -= CRX;
241
242 const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
243 const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
244 p -= CEC;
245 p -= CES;
246
247 return p;
248}
249
250ZdAz MBending::Correct(const ZdAz &zdaz) const
251{
252 // Correct [rad]
253 // zdaz [rad]
254 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
255 AltAz c = Correct(p);
256 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
257}
258
259ZdAz MBending::CorrectBack(const ZdAz &zdaz) const
260{
261 // Correct [rad]
262 // zdaz [rad]
263 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
264 AltAz c = CorrectBack(p);
265 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
266}
267
268void MBending::SetParameters(const Double_t *par, Int_t n)
269{
270 Clear();
271
272 switch (n)
273 {
274 case 16:
275 fMagic2=par[15]/kRad2Deg; // Magic User Defined
276 case 15:
277 fMagic1=par[14]/kRad2Deg; // Magic User Defined
278 case 14:
279 fAcec =par[13]/kRad2Deg; // Azimuth Centering Error (cos)
280 case 13:
281 fEcec =par[12]/kRad2Deg; // Elevation Centering Error (cos)
282 case 12:
283 fAces =par[11]/kRad2Deg; // Azimuth Centering Error (sin)
284 case 11:
285 fEces =par[10]/kRad2Deg; // Elevation Centering Error (sin)
286 case 10:
287 fCry =par[9]/kRad2Deg; // Alt/Az Coude Displacement (E-W)
288 case 9:
289 fCrx =par[8]/kRad2Deg; // Alt/Az Coude Displacement (N-S)
290 case 8:
291 fNry =par[7]/kRad2Deg; // Nasmyth rotator displacement, vertical
292 case 7:
293 fNrx =par[6]/kRad2Deg; // Nasmyth rotator displacement, horizontan
294 case 6:
295 fAw =par[5]/kRad2Deg; // Azimuth Axis Misalignment (E-W)
296 case 5:
297 fAn =par[4]/kRad2Deg; // Azimuth Axis Misalignment (N-S)
298 case 4:
299 fCa =par[3]/kRad2Deg; // Left-Right Collimation Error
300 case 3:
301 fNpae =par[2]/kRad2Deg; // Az-El Nonperpendicularity
302 case 2:
303 fIe =par[1]/kRad2Deg; // Index Error in Elevation
304 case 1:
305 fIa =par[0]/kRad2Deg; // Index Error in Azimuth
306 }
307}
308
309void MBending::GetParameters(Double_t *par, Int_t n) const
310{
311 switch (n)
312 {
313 case 16:
314 par[15]=fMagic2*kRad2Deg; //
315 case 15:
316 par[14]=fMagic1*kRad2Deg; //
317 case 14:
318 par[13]=fAcec*kRad2Deg; // Azimuth Centering Error (cos)
319 case 13:
320 par[12]=fEcec*kRad2Deg; // Elevation Centering Error (cos)
321 case 12:
322 par[11]=fAces*kRad2Deg; // Azimuth Centering Error (sin)
323 case 11:
324 par[10]=fEces*kRad2Deg; // Elevation Centering Error (sin)
325 case 10:
326 par[9]=fCry*kRad2Deg; // Alt/Az Coude Displacement (E-W)
327 case 9:
328 par[8]=fCrx*kRad2Deg; // Alt/Az Coude Displacement (N-S)
329 case 8:
330 par[7]=fNry*kRad2Deg; // Nasmyth rotator displacement, vertical
331 case 7:
332 par[6]=fNrx*kRad2Deg; // Nasmyth rotator displacement, horizontan
333 case 6:
334 par[5]=fAw*kRad2Deg; // Azimuth Axis Misalignment (E-W)
335 case 5:
336 par[4]=fAn*kRad2Deg; // Azimuth Axis Misalignment (N-S)
337 case 4:
338 par[3]=fCa*kRad2Deg; // Left-Right Collimation Error
339 case 3:
340 par[2]=fNpae*kRad2Deg; // Az-El Nonperpendicularity
341 case 2:
342 par[1]=fIe*kRad2Deg; // Index Error in Elevation
343 case 1:
344 par[0]=fIa*kRad2Deg; // Index Error in Azimuth
345 }
346}
347
348void MBending::SetMinuitParameters(TMinuit &m, Int_t n=-1) const
349{
350 if (n<0)
351 n = m.GetNumPars();
352
353 Int_t ierflg = 0;
354
355 switch (n)
356 {
357 case 16:
358 m.mnparm(15,"MAGIC2", fMagic2*kRad2Deg, 1, -360, 360, ierflg);
359 // cout << "Init 3 CA: " << fCa << endl;
360 case 15:
361 m.mnparm(14,"MAGIC1", fMagic1*kRad2Deg, 1, -360, 360, ierflg);
362 // cout << "Init 3 CA: " << fCa << endl;
363 case 14:
364 m.mnparm(13,"ACEC", fAcec*kRad2Deg, 1, -360, 360, ierflg);
365 // cout << "Init 3 CA: " << fCa << endl;
366 case 13:
367 m.mnparm(12,"ECEC", fEcec*kRad2Deg, 1, -360, 360, ierflg);
368 // cout << "Init 3 CA: " << fCa << endl;
369 case 12:
370 m.mnparm(11,"ACES", fAcec*kRad2Deg, 1, -360, 360, ierflg);
371 // cout << "Init 3 CA: " << fCa << endl;
372 case 11:
373 m.mnparm(10,"ECES", fEcec*kRad2Deg, 1, -360, 360, ierflg);
374 // cout << "Init 3 CA: " << fCa << endl;
375 case 10:
376 m.mnparm(9, "CRY", fCry*kRad2Deg, 1, -360, 360, ierflg);
377 // cout << "Init 3 CA: " << fCa << endl;
378 case 9:
379 m.mnparm(8, "CRX", fCrx*kRad2Deg, 1, -360, 360, ierflg);
380 // cout << "Init 3 CA: " << fCa << endl;
381 case 8:
382 m.mnparm(7, "NRY", fNry*kRad2Deg, 1, -360, 360, ierflg);
383 // cout << "Init 3 CA: " << fCa << endl;
384 case 7:
385 m.mnparm(6, "NRX", fNrx*kRad2Deg, 1, -360, 360, ierflg);
386 // cout << "Init 3 CA: " << fCa << endl;
387 case 6:
388 m.mnparm(5, "AW", fAw*kRad2Deg, 1, -360, 360, ierflg);
389 // cout << "Init 3 CA: " << fCa << endl;
390 case 5:
391 m.mnparm(4, "AN", fAn*kRad2Deg, 1, -360, 360, ierflg);
392 // cout << "Init 3 CA: " << fCa << endl;
393 case 4:
394 m.mnparm(3, "CA", fCa*kRad2Deg, 1, -360, 360, ierflg);
395 // cout << "Init 3 CA: " << fCa << endl;
396 case 3:
397 m.mnparm(2, "NPAE", fNpae*kRad2Deg, 1, -360, 360, ierflg);
398 // cout << "Init 2 NPAE: " << fNpae << endl;
399 case 2:
400 m.mnparm(1, "IE", fIe*kRad2Deg, 1, -360, 360, ierflg);
401 // cout << "Init 1 IE: " << fIe << endl;
402 case 1:
403 m.mnparm(0, "IA", fIa*kRad2Deg, 1, -360, 360, ierflg);
404 // cout << "Init 0 IA: " << fIa << endl;
405 }
406}
407
408void MBending::GetMinuitParameters(TMinuit &m, Int_t n=-1)
409{
410 if (n<0)
411 n = m.GetNumPars();
412
413 Double_t err;
414
415 switch (n)
416 {
417 case 16:
418 m.GetParameter(15, fMagic2, err);
419 fMagic2 /= kRad2Deg;
420 case 15:
421 m.GetParameter(14, fMagic1, err);
422 fMagic1 /= kRad2Deg;
423 case 14:
424 m.GetParameter(13, fAcec, err);
425 fAcec /= kRad2Deg;
426 case 13:
427 m.GetParameter(12, fEcec, err);
428 fEcec /= kRad2Deg;
429 case 12:
430 m.GetParameter(11, fAces, err);
431 fAces /= kRad2Deg;
432 case 11:
433 m.GetParameter(10, fEces, err);
434 fEces /= kRad2Deg;
435 case 10:
436 m.GetParameter(9, fCry, err);
437 fCry /= kRad2Deg;
438 case 9:
439 m.GetParameter(8, fCrx, err);
440 fCrx /= kRad2Deg;
441 case 8:
442 m.GetParameter(7, fNry, err);
443 fNry /= kRad2Deg;
444 case 7:
445 m.GetParameter(6, fNrx, err);
446 fNrx /= kRad2Deg;
447 case 6:
448 m.GetParameter(5, fAw, err);
449 fAw /= kRad2Deg;
450 case 5:
451 m.GetParameter(4, fAn, err);
452 fAn /= kRad2Deg;
453 case 4:
454 m.GetParameter(3, fCa, err);
455 fCa /= kRad2Deg;
456 case 3:
457 m.GetParameter(2, fNpae, err);
458 fNpae /= kRad2Deg;
459 case 2:
460 m.GetParameter(1, fIe, err);
461 fIe /= kRad2Deg;
462 case 1:
463 m.GetParameter(0, fIa, err);
464 fIa /= kRad2Deg;
465 }
466}
467
468void MBending::PrintMinuitParameters(TMinuit &m, Int_t n=-1) const
469{
470 if (n<0)
471 n = m.GetNumPars();
472
473 Double_t par, err;
474
475 switch (n)
476 {
477 case 16:
478 m.GetParameter(15, par, err);
479 cout << " 15 MAGIC2: " << par << " +- " << err << endl;
480 case 15:
481 m.GetParameter(14, par, err);
482 cout << " 14 MAGIC1: " << par << " +- " << err << endl;
483 case 14:
484 m.GetParameter(13, par, err);
485 cout << " 13 ACEC: " << par << " +- " << err << endl;
486 case 13:
487 m.GetParameter(12, par, err);
488 cout << " 12 ECEC: " << par << " +- " << err << endl;
489 case 12:
490 m.GetParameter(11, par, err);
491 cout << " 11 ACES: " << par << " +- " << err << endl;
492 case 11:
493 m.GetParameter(10, par, err);
494 cout << " 10 ECES: " << par << " +- " << err << endl;
495 case 10:
496 m.GetParameter(9, par, err);
497 cout << " 9 CRY: " << par << " +- " << err << endl;
498 case 9:
499 m.GetParameter(8, par, err);
500 cout << " 8 CRX: " << par << " +- " << err << endl;
501 case 8:
502 m.GetParameter(7, par, err);
503 cout << " 7 NRY: " << par << " +- " << err << endl;
504 case 7:
505 m.GetParameter(6, par, err);
506 cout << " 6 NRX: " << par << " +- " << err << endl;
507 case 6:
508 m.GetParameter(5, par, err);
509 cout << " 5 AW: " << par << " +- " << err << endl;
510 case 5:
511 m.GetParameter(4, par, err);
512 cout << " 4 AN: " << par << " +- " << err << endl;
513 case 4:
514 m.GetParameter(3, par, err);
515 cout << " 3 CA: " << par << " +- " << err << endl;
516 case 3:
517 m.GetParameter(2, par, err);
518 cout << " 2 NPAE: " << par << " +- " << err << endl;
519 case 2:
520 m.GetParameter(1, par, err);
521 cout << " 1 IE: " << par << " +- " << err << endl;
522 case 1:
523 m.GetParameter(0, par, err);
524 cout << " 0 IA: " << par << " +- " << err << endl;
525 }
526}
Note: See TracBrowser for help on using the repository browser.