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

Last change on this file since 2408 was 2407, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 12.3 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 2003 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MBending
28//
29// Double_t fIe ; // [rad] Index Error in Elevation
30// Double_t fIa ; // [rad] Index Error in Azimuth
31//
32// Double_t fFlop ; // [rad] Vertical Sag
33// * do not use if not data: Zd<0
34//
35// Double_t fNpae ; // [rad] Az-El Nonperpendicularity
36//
37// Double_t fCa ; // [rad] Left-Right Collimation Error
38//
39// Double_t fAn ; // [rad] Azimuth Axis Misalignment (N-S)
40// Double_t fAw ; // [rad] Azimuth Axis Misalignment (E-W)
41//
42// Double_t fTf ; // [rad] Tube fluxture (sin)
43// * same as ecec if no data: Zd<0
44// Double_t fTx ; // [rad] Tube fluxture (tan)
45// * do not use with NPAE if no data: Zd<0
46//
47// Double_t fNrx ; // [rad] Nasmyth rotator displacement, horizontan
48// Double_t fNry ; // [rad] Nasmyth rotator displacement, vertical
49//
50// Double_t fCrx ; // [rad] Alt/Az Coude Displacement (N-S)
51// Double_t fCry ; // [rad] Alt/Az Coude Displacement (E-W)
52//
53// Double_t fEces ; // [rad] Elevation Centering Error (sin)
54// Double_t fAces ; // [rad] Azimuth Centering Error (sin)
55// Double_t fEcec ; // [rad] Elevation Centering Error (cos)
56// Double_t fAcec ; // [rad] Azimuth Centering Error (cos)
57//
58////////////////////////////////////////////////////////////////////////////
59#include "MBending.h"
60
61#include <fstream.h>
62#include <iomanip.h>
63
64#include <TMinuit.h>
65
66#include "timer.h"
67
68ClassImp(MBending);
69
70const Int_t MBending::fNumPar=19;
71
72void MBending::Init()
73{
74 fCoeff = new Double_t*[fNumPar];
75 fName = new TString[fNumPar];
76 fDescr = new TString[fNumPar];
77
78 fCoeff[ 0] = &fIa; fName[ 0] = "IA";
79 fCoeff[ 1] = &fIe; fName[ 1] = "IE";
80 fCoeff[ 2] = &fFlop; fName[ 2] = "FLOP";
81 fCoeff[ 3] = &fNpae; fName[ 3] = "NPAE";
82 fCoeff[ 4] = &fCa; fName[ 4] = "CA";
83 fCoeff[ 5] = &fAn; fName[ 5] = "AN";
84 fCoeff[ 6] = &fAw; fName[ 6] = "AW";
85 fCoeff[ 7] = &fTf; fName[ 7] = "TF";
86 fCoeff[ 8] = &fTx; fName[ 8] = "TX";
87 fCoeff[ 9] = &fEces; fName[ 9] = "ECES";
88 fCoeff[10] = &fAces; fName[10] = "ACES";
89 fCoeff[11] = &fEcec; fName[11] = "ECEC";
90 fCoeff[12] = &fAcec; fName[12] = "ACEC";
91 fCoeff[13] = &fNrx; fName[13] = "NRX";
92 fCoeff[14] = &fNry; fName[14] = "NRY";
93 fCoeff[15] = &fCrx; fName[15] = "CRX";
94 fCoeff[16] = &fCry; fName[16] = "CRY";
95 fCoeff[17] = &fMagic1; fName[17] = "MAGIC1";
96 fCoeff[18] = &fMagic2; fName[18] = "MAGIC2";
97
98 fDescr[ 0] = "Index Error Azimuth";
99 fDescr[ 1] = "Index Error Zenith Distance";
100 fDescr[ 2] = "Vertical Sag";
101 fDescr[ 3] = "Az-El Nonperpendicularity";
102 fDescr[ 4] = "Left-Right Collimation Error";
103 fDescr[ 5] = "Azimuth Axis Misalignment (N-S)";
104 fDescr[ 6] = "Azimuth Axis Misalignment (E-W)";
105 fDescr[ 7] = "Tube fluxture (sin)";
106 fDescr[ 8] = "Tube fluxture (tan)";
107 fDescr[ 9] = "Elevation Centering Error (sin)";
108 fDescr[10] = "Azimuth Centering Error (sin)";
109 fDescr[11] = "Elevation Centering Error (cos)";
110 fDescr[12] = "Azimuth Centering Error (cos)";
111 fDescr[13] = "Nasmyth rotator displacement (horizontal)";
112 fDescr[14] = "Nasmyth rotator displacement (vertical)";
113 fDescr[15] = "Alt/Az Coude Displacement (N-S)";
114 fDescr[16] = "Alt/Az Coude Displacement (E-W)";
115 fDescr[17] = "n/a";
116 fDescr[18] = "n/a";
117}
118void MBending::Reset()
119{
120 Clear();
121}
122
123void MBending::Load(const char *name)
124{
125 /*
126 ! MMT 1987 July 8
127 ! T 36 7.3622 41.448 -0.0481
128 ! IA -37.5465 20.80602
129 ! IE -13.9180 1.25217
130 ! NPAE +7.0751 26.44763
131 ! CA -6.9149 32.05358
132 ! AN +0.5053 1.40956
133 ! AW -2.2016 1.37480
134 ! END
135 */
136
137 ifstream fin(name);
138 if (!fin)
139 {
140 cout << "Error: Cannot open file '" << name << "'" << endl;
141 return;
142 }
143
144 char c;
145 while (fin && fin.get()!='\n');
146 fin >> c;
147
148 if (c!='S' && c!='s')
149 {
150 cout << "Error: This in not a model correcting the star position (" << c << ")" << endl;
151 return;
152 }
153
154 Clear();
155
156 cout << endl;
157
158 Double_t val;
159 fin >> val;
160 cout << "Number of observed stars: " << val << endl;
161 fin >> val;
162 cout << "Sky RMS: " << val << "\"" << endl;
163 fin >> val;
164 cout << "Refraction Constant A: " << val << "\"" << endl;
165 fin >> val;
166 cout << "Refraction Constant B: " << val << "\"" << endl;
167
168 cout << endl;
169
170 cout << " & = Name Value Sigma" << endl;
171 cout << "--------------------------------------------------" << endl;
172
173 while (fin)
174 {
175 TString str;
176 fin >> str;
177
178 if (str=="END")
179 break;
180
181 if (str[0]=='&')
182 {
183 cout << " & ";
184 str.Remove(0);
185 }
186 else
187 cout << " ";
188
189 if (str[1]=='=')
190 {
191 cout << "= ";
192 str.Remove(0);
193 }
194 else
195 cout << " ";
196
197 fin >> val;
198 cout << str << "\t" << setw(11) << val << "° \t";
199 val *= kDeg2Rad;
200
201 Double_t *dest=NULL;
202
203 if (str=="IA") dest = &fIa;
204 if (str=="IE") dest = &fIe;
205 if (str=="FLOP") dest = &fFlop;
206 if (str=="NPAE") dest = &fNpae;
207 if (str=="CA") dest = &fCa;
208 if (str=="AN") dest = &fAn;
209 if (str=="AW") dest = &fAw;
210 if (str=="TF") dest = &fTf;
211 if (str=="TX") dest = &fTx;
212 if (str=="NRX") dest = &fNrx;
213 if (str=="NRY") dest = &fNry;
214 if (str=="CRX") dest = &fCrx;
215 if (str=="CRY") dest = &fCry;
216 if (str=="ECES") dest = &fEces;
217 if (str=="ACES") dest = &fAces;
218 if (str=="ECEC") dest = &fEcec;
219 if (str=="ACEC") dest = &fAcec;
220
221 if (dest)
222 *dest = val;
223
224 fin >> val;
225 cout << setw(9) << val << "°" << endl;
226
227 // Find corresponding error
228 for (int i=0; i<MBending::GetNumPar(); i++)
229 if (dest==fCoeff[i])
230 {
231 fError[i] = val*kDeg2Rad;
232 break;
233 }
234 }
235 cout << endl;
236}
237
238void MBending::Save(const char *name)
239{
240 /*
241 ! MMT 1987 July 8
242 ! T 36 7.3622 41.448 -0.0481
243 ! IA -37.5465 20.80602
244 ! IE -13.9180 1.25217
245 ! NPAE +7.0751 26.44763
246 ! CA -6.9149 32.05358
247 ! AN +0.5053 1.40956
248 ! AW -2.2016 1.37480
249 ! END
250 */
251
252 ofstream fout(name);
253 if (!fout)
254 {
255 cout << "Error: Cannot open file '" << name << "'" << endl;
256 return;
257 }
258
259 Timer t;
260 t.Now();
261
262 fout << "MAGIC1 " << t.GetTimeStr() << endl;
263 fout << "S 00 000000 000000 0000000" << endl;
264 fout << setprecision(8);
265 for (int i=0; i<fNumPar; i++)
266 {
267 fout << " " << setw(6) << GetName(i) << " ";
268 fout << setw(13) << *fCoeff[i]*kRad2Deg << " ";
269 fout << setw(11) << fError[i]*kRad2Deg << endl;
270 }
271 fout << "END" << endl;
272}
273
274Double_t MBending::Sign(Double_t val, Double_t alt)
275{
276 // Some pointing corrections are defined as Delta ZA, which
277 // is (P. Wallace) defined [0,90]deg while Alt is defined
278 // [0,180]deg
279 return (TMath::Pi()/2-alt < 0 ? -val : val);
280}
281
282AltAz MBending::Correct(const AltAz &aa) const
283{
284 // Correct [rad]
285 // zdaz [rad]
286 AltAz p = aa;
287
288 const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
289 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
290 p += CRX;
291 p += CRY;
292
293 const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
294 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
295 p += NRX;
296 p += NRY;
297
298 const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
299 const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
300 p += CES;
301 p += CEC;
302
303 const AltAz TX(Sign(fTx/tan(p.Alt()), p.Alt()), 0);
304 const AltAz TF(Sign(fTf*cos(p.Alt()), p.Alt()), 0);
305 p += TX;
306 p += TF;
307
308 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
309 const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
310 p += AW;
311 p += AN;
312
313 const AltAz CA(0, -fCa/cos(p.Alt()));
314 p += CA;
315
316 const AltAz NPAE(0, -fNpae*tan(p.Alt()));
317 p += NPAE;
318
319 const AltAz FLOP(Sign(fFlop, p.Alt()), 0);
320 p += FLOP;
321
322 const AltAz I(fIe, fIa);
323 p += I;
324
325 return p;
326}
327
328AltAz MBending::CorrectBack(const AltAz &aa) const
329{
330 // Correct [rad]
331 // zdaz [rad]
332 AltAz p = aa;
333
334 const AltAz I(fIe, fIa);
335 p -= I;
336
337 const AltAz FLOP(Sign(fFlop, p.Alt()), 0);
338 p -= FLOP;
339
340 const AltAz NPAE(0, -fNpae*tan(p.Alt()));
341 p -= NPAE;
342
343 const AltAz CA(0, -fCa/cos(p.Alt()));
344 p -= CA;
345
346 const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
347 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
348 p -= AN;
349 p -= AW;
350
351 const AltAz TF(Sign(fTf*cos(p.Alt()), p.Alt()), 0);
352 const AltAz TX(Sign(fTx/tan(p.Alt()), p.Alt()), 0);
353 p -= TF;
354 p -= TX;
355
356 const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
357 const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
358 p -= CEC;
359 p -= CES;
360
361 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
362 const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
363 p -= NRY;
364 p -= NRX;
365
366 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
367 const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
368 p -= CRY;
369 p -= CRX;
370
371 return p;
372}
373
374ZdAz MBending::Correct(const ZdAz &zdaz) const
375{
376 // Correct [rad]
377 // zdaz [rad]
378 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
379 AltAz c = Correct(p);
380 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
381}
382
383ZdAz MBending::CorrectBack(const ZdAz &zdaz) const
384{
385 // Correct [rad]
386 // zdaz [rad]
387 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
388 AltAz c = CorrectBack(p);
389 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
390}
391
392void MBending::SetParameters(const Double_t *par, Int_t n)
393{
394 Clear();
395
396 while (n--)
397 *fCoeff[n] = par[n]/kRad2Deg;
398}
399
400void MBending::GetParameters(Double_t *par, Int_t n) const
401{
402 while (n--)
403 par[n] = *fCoeff[n]*kRad2Deg;
404}
405
406void MBending::SetMinuitParameters(TMinuit &m, Int_t n=-1) const
407{
408 if (n<0)
409 n = fNumPar;
410
411 Int_t ierflg = 0;
412
413 while (n--)
414 m.mnparm(n, fName[n], *fCoeff[n]*kRad2Deg, 1, -360, 360, ierflg);
415}
416
417void MBending::GetMinuitParameters(TMinuit &m, Int_t n=-1)
418{
419 if (n<0 || n>m.GetNumPars())
420 n = m.GetNumPars();
421
422 while (n--)
423 {
424 m.GetParameter(n, *fCoeff[n], fError[n]);
425 *fCoeff[n] /= kRad2Deg;
426 fError[n] /= kRad2Deg;
427 }
428}
429/*
430void FormatPar(TMinuit &m, Int_t n)
431{
432 Double_t par, err;
433 m.GetParameter(n, par, err);
434
435 int expp = (int)log10(par);
436 int expe = (int)log10(err);
437
438 if (err<2*pow(10, expe))
439 expe--;
440
441 Int_t exp = expe>expp ? expp : expe;
442
443 par = (int)(par/pow(10, exp)) * pow(10, exp);
444 err = (int)(err/pow(10, exp)) * pow(10, exp);
445
446 cout << par << " +- " << err << flush;
447}
448*/
449void MBending::PrintMinuitParameters(TMinuit &m, Int_t n=-1) const
450{
451 if (n<0)
452 n = m.GetNumPars();
453
454 cout << setprecision(3);
455
456 Double_t par, err;
457
458 while (n--)
459 {
460 m.GetParameter(n, par, err);
461 cout << Form(" %2d %6s: ", n, (const char*)fName[n]);
462 cout << setw(8) << par << " \xb1 " << setw(6) << err << endl;
463 }
464}
Note: See TracBrowser for help on using the repository browser.