source: branches/FACT++_lidctrl_usb/drive/MPointing.cc@ 19364

Last change on this file since 19364 was 18618, checked in by tbretz, 10 years ago
Copied files from Mars and Cosy to compile toold which fit our current drivectrl.
File size: 22.6 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-2007
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MPointing
28// =========
29//
30// This is the class used for the pointing correction done in the MAGIC
31// drive software cosy. NEVER CHANGE IT WITHOUT CONTACTING THE AUTHOR FIRST!
32//
33// Variables/Coefficients
34// ----------------------
35//
36// Double_t fIe ; // [rad] Index Error in Elevation
37// Double_t fIa ; // [rad] Index Error in Azimuth
38//
39// Double_t fFlop ; // [rad] Vertical Sag
40// * do not use if not data: Zd<0
41//
42// Double_t fNpae ; // [rad] Az-El Nonperpendicularity
43//
44// Double_t fCa ; // [rad] Left-Right Collimation Error
45//
46// Double_t fAn ; // [rad] Azimuth Axis Misalignment (N-S)
47// Double_t fAw ; // [rad] Azimuth Axis Misalignment (E-W)
48//
49// Double_t fTf ; // [rad] Tube fluxture (sin)
50// * same as ecec if no data: Zd<0
51// Double_t fTx ; // [rad] Tube fluxture (tan)
52// * do not use with NPAE if no data: Zd<0
53//
54// Double_t fNrx ; // [rad] Nasmyth rotator displacement, horizontan
55// Double_t fNry ; // [rad] Nasmyth rotator displacement, vertical
56//
57// Double_t fCrx ; // [rad] Alt/Az Coude Displacement (N-S)
58// Double_t fCry ; // [rad] Alt/Az Coude Displacement (E-W)
59//
60// Double_t fEces ; // [rad] Elevation Centering Error (sin)
61// Double_t fAces ; // [rad] Azimuth Centering Error (sin)
62// Double_t fEcec ; // [rad] Elevation Centering Error (cos)
63// Double_t fAcec ; // [rad] Azimuth Centering Error (cos)
64//
65// Double_t fMagic1;// [rad] MAGIC culmination hysteresis
66// Double_t fMagic2;// [rad] undefined
67//
68// Double_t fDx; // [rad] X-offset in camera (for starguider calibration)
69// Double_t fDy; // [rad] Y-offset in camera (for starguider calibration)
70//
71//
72// Class Version 2:
73// ----------------
74// + fPx
75// + fPy
76// + fDx
77// + fDy
78//
79//
80////////////////////////////////////////////////////////////////////////////
81#include "MPointing.h"
82
83#include <fstream>
84
85#include <TMinuit.h>
86
87#include "MLog.h"
88#include "MLogManip.h"
89
90#include "MTime.h"
91
92ClassImp(AltAz);
93ClassImp(ZdAz);
94ClassImp(RaDec);
95ClassImp(MPointing);
96
97using namespace std;
98
99#undef DEBUG
100//#define DEBUG(txt) txt
101#define DEBUG(txt)
102
103void ZdAz::Round()
104{
105 fX = TMath::Nint(fX);
106 fY = TMath::Nint(fY);
107}
108
109void ZdAz::Abs()
110{
111 fX = TMath::Abs(fX);
112 fY = TMath::Abs(fY);
113}
114
115void MPointing::Init(const char *name, const char *title)
116{
117 fName = name ? name : "MPointing";
118 fTitle = title ? title : "Pointing correction model for the MAGIC telescope";
119
120 fCoeff = new Double_t*[kNumPar];
121 fNames = new TString[kNumPar];
122 fDescr = new TString[kNumPar];
123
124 fCoeff[kIA] = &fIa; fNames[kIA] = "IA";
125 fCoeff[kIE] = &fIe; fNames[kIE] = "IE";
126 fCoeff[kFLOP] = &fFlop; fNames[kFLOP] = "FLOP";
127 fCoeff[kAN] = &fAn; fNames[kAN] = "AN";
128 fCoeff[kAW] = &fAw; fNames[kAW] = "AW";
129 fCoeff[kNPAE] = &fNpae; fNames[kNPAE] = "NPAE";
130 fCoeff[kCA] = &fCa; fNames[kCA] = "CA";
131 fCoeff[kTF] = &fTf; fNames[kTF] = "TF";
132 fCoeff[kTX] = &fTx; fNames[kTX] = "TX";
133 fCoeff[kECES] = &fEces; fNames[kECES] = "ECES";
134 fCoeff[kACES] = &fAces; fNames[kACES] = "ACES";
135 fCoeff[kECEC] = &fEcec; fNames[kECEC] = "ECEC";
136 fCoeff[kACEC] = &fAcec; fNames[kACEC] = "ACEC";
137 fCoeff[kNRX] = &fNrx; fNames[kNRX] = "NRX";
138 fCoeff[kNRY] = &fNry; fNames[kNRY] = "NRY";
139 fCoeff[kCRX] = &fCrx; fNames[kCRX] = "CRX";
140 fCoeff[kCRY] = &fCry; fNames[kCRY] = "CRY";
141 fCoeff[kMAGIC1] = &fMagic1; fNames[kMAGIC1] = "MAGIC1";
142 fCoeff[kMAGIC2] = &fMagic2; fNames[kMAGIC2] = "MAGIC2";
143 fCoeff[kPX] = &fPx; fNames[kPX] = "PX";
144 fCoeff[kPY] = &fPy; fNames[kPY] = "PY";
145 fCoeff[kDX] = &fDx; fNames[kDX] = "DX";
146 fCoeff[kDY] = &fDy; fNames[kDY] = "DY";
147
148 fDescr[kIA] = "Index Error Azimuth";
149 fDescr[kIE] = "Index Error Zenith Distance";
150 fDescr[kFLOP] = "Vertical Sag";
151 fDescr[kAN] = "Azimuth Axis Misalignment (N-S)";
152 fDescr[kAW] = "Azimuth Axis Misalignment (E-W)";
153 fDescr[kNPAE] = "Az-El Nonperpendicularity";
154 fDescr[kCA] = "Left-Right Collimation Error";
155 fDescr[kTF] = "Tube fluxture (sin)";
156 fDescr[kTX] = "Tube fluxture (tan)";
157 fDescr[kECES] = "Elevation Centering Error (sin)";
158 fDescr[kACES] = "Azimuth Centering Error (sin)";
159 fDescr[kECEC] = "Elevation Centering Error (cos)";
160 fDescr[kACEC] = "Azimuth Centering Error (cos)";
161 fDescr[kNRX] = "Nasmyth rotator displacement (horizontal)";
162 fDescr[kNRY] = "Nasmyth rotator displacement (vertical)";
163 fDescr[kCRX] = "Alt/Az Coude Displacement (N-S)";
164 fDescr[kCRY] = "Alt/Az Coude Displacement (E-W)";
165 fDescr[kMAGIC1] = "MAGIC culmination hysteresis";
166 fDescr[kMAGIC2] = "n/a";
167 fDescr[kPX] = "Starguider calibration fixed offset x";
168 fDescr[kPY] = "Starguider calibration fixed offset y";
169 fDescr[kDX] = "Starguider calibration additional offset dx";
170 fDescr[kDY] = "Starguider calibration additional offset dy";
171}
172
173void MPointing::Reset()
174{
175 Clear();
176}
177
178Bool_t MPointing::Load(const char *name)
179{
180 /*
181 ! MMT 1987 July 8
182 ! T 36 7.3622 41.448 -0.0481
183 ! IA -37.5465 20.80602
184 ! IE -13.9180 1.25217
185 ! NPAE +7.0751 26.44763
186 ! CA -6.9149 32.05358
187 ! AN +0.5053 1.40956
188 ! AW -2.2016 1.37480
189 ! END
190 */
191
192 ifstream fin(name);
193 if (!fin)
194 {
195 *fLog << err << "ERROR - Cannot open file '" << name << "'" << endl;
196 return kFALSE;
197 }
198
199 char c;
200 while (fin && fin.get()!='\n');
201 fin >> c;
202
203 if (c!='S' && c!='s')
204 {
205 *fLog << err << "Error: This in not a model correcting the star position (" << c << ")" << endl;
206 return kFALSE;
207 }
208
209 Clear();
210
211 cout << endl;
212
213 Double_t val;
214 fin >> val;
215 *fLog << inf;
216 //*fLog << "Number of observed stars: " << val << endl;
217 fin >> val;
218 //*fLog << "Sky RMS: " << val << "\"" << endl;
219 fin >> val;
220 //*fLog << "Refraction Constant A: " << val << "\"" << endl;
221 fin >> val;
222 //*fLog << "Refraction Constant B: " << val << "\"" << endl;
223
224 *fLog << endl;
225
226 *fLog << " & = Name Value Sigma " << endl;
227 *fLog << "--------------------------------------------------" << endl;
228
229 while (fin)
230 {
231 TString str;
232 fin >> str;
233 if (!fin)
234 {
235 *fLog << err << "ERROR - Reading file " << name << endl;
236 return kFALSE;
237 }
238
239 str = str.Strip(TString::kBoth);
240
241 if (str=="END")
242 break;
243
244 TString sout;
245
246 if (str[0]=='#')
247 continue;
248
249 if (str[0]=='&')
250 {
251 sout += " & ";
252 str.Remove(0);
253 }
254 else
255 sout += " ";
256
257 if (str[1]=='=')
258 {
259 sout += "= ";
260 str.Remove(0);
261 }
262 else
263 sout += " ";
264
265 fin >> val;
266
267 sout += str;
268 sout += '\t';
269 sout += Form("%11f", val);
270 sout += UTF8::kDeg;
271 sout += " \t";
272 val *= TMath::DegToRad();
273
274 // Find parameter
275 Int_t n = -1;
276 for (int i=0; i<kNumPar; i++)
277 if (str==fNames[i])
278 {
279 n = i;
280 *fCoeff[i] = val;
281 break;
282 }
283
284 fin >> val;
285 sout += Form("%9f%s", val, UTF8::kDeg);
286
287 if (*fCoeff[n]!=0 || val>0)
288 *fLog << sout << endl;
289
290 if (!fin)
291 {
292 *fLog << err << "ERROR - Reading line " << str << endl;
293 return kFALSE;
294 }
295
296 if (n<0)
297 {
298 *fLog << warn << "WARNING - Parameter " << str << " unknown." << endl;
299 continue;
300 }
301
302 // corresponding error
303 fError[n] = val*TMath::DegToRad();
304 }
305 *fLog << endl;
306
307 fName = name;
308
309 return kTRUE;
310}
311
312Bool_t MPointing::Save(const char *name)
313{
314 /*
315 ! MMT 1987 July 8
316 ! T 36 7.3622 41.448 -0.0481
317 ! IA -37.5465 20.80602
318 ! IE -13.9180 1.25217
319 ! NPAE +7.0751 26.44763
320 ! CA -6.9149 32.05358
321 ! AN +0.5053 1.40956
322 ! AW -2.2016 1.37480
323 ! END
324 */
325
326 ofstream fout(name);
327 if (!fout)
328 {
329 cout << "Error: Cannot open file '" << name << "'" << endl;
330 return kFALSE;
331 }
332
333 MTime t;
334 t.Now();
335
336 fout << "MAGIC1 " << t << endl;
337 fout << "S 00 000000 000000 0000000" << endl;
338 fout << setprecision(8);
339 for (int i=0; i<kNumPar; i++)
340 {
341 fout << " " << setw(6) << GetVarName(i) << " ";
342 fout << setw(13) << *fCoeff[i]*kRad2Deg << " ";
343 fout << setw(11) << fError[i]*kRad2Deg << endl;
344 }
345 fout << "END" << endl;
346
347 fName = name;
348
349 return kTRUE;
350}
351
352Double_t MPointing::Sign(Double_t val, Double_t alt)
353{
354 // Some pointing corrections are defined as Delta ZA, which
355 // is (P. Wallace) defined [0,90]deg while Alt is defined
356 // [0,180]deg
357 return (TMath::Pi()/2-alt < 0 ? -val : val);
358}
359
360AltAz MPointing::AddOffsets(const AltAz &aa) const
361{
362 // Correct [rad]
363 // zdaz [rad]
364 AltAz p = aa;
365
366 const AltAz I(fIe, fIa);
367 p += I;
368
369 return p;
370}
371
372ZdAz MPointing::AddOffsets(const ZdAz &zdaz) const
373{
374 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
375
376 AltAz c = AddOffsets(p);
377
378 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
379}
380
381TVector3 MPointing::AddOffsets(const TVector3 &v) const
382{
383 AltAz p(TMath::Pi()/2-v.Theta(), v.Phi());
384 AltAz c = AddOffsets(p);
385
386 TVector3 rc;
387 rc.SetMagThetaPhi(1, TMath::Pi()/2-c.Alt(), c.Az());
388 return rc;
389}
390
391AltAz MPointing::SubtractOffsets(const AltAz &aa) const
392{
393 // Correct [rad]
394 // zdaz [rad]
395 AltAz p = aa;
396
397 const AltAz I(fIe, fIa);
398 p -= I;
399
400 return p;
401}
402
403ZdAz MPointing::SubtractOffsets(const ZdAz &zdaz) const
404{
405 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
406
407 AltAz c = SubtractOffsets(p);
408
409 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
410}
411
412TVector3 MPointing::SubtractOffsets(const TVector3 &v) const
413{
414 AltAz p(TMath::Pi()/2-v.Theta(), v.Phi());
415 AltAz c = SubtractOffsets(p);
416
417 TVector3 rc;
418 rc.SetMagThetaPhi(1, TMath::Pi()/2-c.Alt(), c.Az());
419 return rc;
420}
421
422AltAz MPointing::CalcAnAw(const AltAz &p, Int_t sign) const
423{
424 // Corrections for AN and AW without approximations
425 // as done by Patrick Wallace. The approximation cannot
426 // be used for MAGIC because the correctioon angle
427 // AW (~1.5deg) is not small enough.
428
429 // Vector in cartesian coordinates
430 TVector3 v1;
431
432 // Set Azimuth and Elevation
433 v1.SetMagThetaPhi(1, TMath::Pi()/2-p.Alt(), p.Az());
434
435
436 TVector3 v2(v1);
437// cout << sign << endl;
438
439// cout << "v1: " << v1.Theta()*TMath::RadToDeg() << " " << v1.Phi()*TMath::RadToDeg() << endl;
440
441 // Rotate around the x- and y-axis
442 v1.RotateY(sign*fAn);
443 v1.RotateX(sign*fAw);
444
445// cout << "v1: " << v1.Theta()*TMath::RadToDeg() << " " << v1.Phi()*TMath::RadToDeg() << endl;
446// cout << "v2: " << v2.Theta()*TMath::RadToDeg() << " " << v2.Theta()*TMath::RadToDeg() << endl;
447
448 // cout << "dv: " << (v2.Theta()-v1.Theta())*TMath::RadToDeg() << " " << (v2.Phi()-v1.Phi())*TMath::RadToDeg() << endl;
449
450 Double_t dalt = v1.Theta()-v2.Theta();
451 Double_t daz = v1.Phi() -v2.Phi();
452
453 //cout << dalt*TMath::RadToDeg() << " " << daz*TMath::RadToDeg() << endl;
454
455 if (daz>TMath::Pi())
456 daz -= TMath::TwoPi();
457 if (daz<-TMath::Pi())
458 daz += TMath::TwoPi();
459
460// if (daz>TMath::Pi()/2)
461// {
462// }
463
464 AltAz d(dalt, daz);
465 return d;
466
467 // Calculate Delta Azimuth and Delta Elevation
468 /*
469 AltAz d(TMath::Pi()/2-v1.Theta(), v1.Phi());
470
471 cout << "p : " << p.Alt()*TMath::RadToDeg() << " " << p.Az()*TMath::RadToDeg() << endl;
472 cout << "d : " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
473 d -= p;
474 cout << "d-p: " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
475 d *= sign;
476 cout << "d* : " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
477
478
479 cout << "p2: " << 90-p.Alt()*TMath::RadToDeg() << " " << p.Az()*TMath::RadToDeg() << endl;
480 cout << "d2: " << 90-d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
481
482 Int_t s1 = 90-d.Alt()*TMath::RadToDeg() < 0 ? -1 : 1;
483 Int_t s2 = 90-p.Alt()*TMath::RadToDeg() < 0 ? -1 : 1;
484
485
486 if (s1 != s2)
487 {
488 //90-d.Alt() <-- -90+d.Alt()
489
490 d.Alt(d.Alt()-TMath::Pi());
491 cout << "Alt-" << endl;
492 }
493 cout << "d': " << 90-d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;*/
494 /*
495 // Fix 'direction' of output depending on input vector
496 if (TMath::Pi()/2-sign*p.Alt()<0)
497 {
498 d.Alt(d.Alt()-TMath::Pi());
499 cout << "Alt-" << endl;
500 }
501 //if (TMath::Pi()/2-sign*p.Alt()>TMath::Pi())
502 //{
503 // d.Alt(TMath::Pi()-d.Alt());
504 // cout << "Alt+" << endl;
505 //}
506
507 // Align correction into [-180,180]
508 while (d.Az()>TMath::Pi())
509 {
510 d.Az(d.Az()-TMath::Pi()*2);
511 cout << "Az-" << endl;
512 }
513 while (d.Az()<-TMath::Pi())
514 {
515 d.Az(d.Az()+TMath::Pi()*2);
516 cout << "Az+" << endl;
517 }
518 */
519 return d;
520}
521
522
523AltAz MPointing::Correct(const AltAz &aa) const
524{
525 // Correct [rad]
526 // zdaz [rad]
527 AltAz p = aa;
528
529 DEBUG(cout << setprecision(16));
530 DEBUG(cout << "Bend7: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
531
532 const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
533 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
534 p += CRX;
535 p += CRY;
536
537 DEBUG(cout << "Bend6: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
538
539 const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
540 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
541 p += NRX;
542 p += NRY;
543
544 DEBUG(cout << "Bend5: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
545
546 const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
547 const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
548 p += CES;
549 p += CEC;
550
551 DEBUG(cout << "Bend4: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
552
553 const AltAz TX(Sign(fTx/tan(p.Alt()), p.Alt()), 0);
554 const AltAz TF(Sign(fTf*cos(p.Alt()), p.Alt()), 0);
555 //p += TX;
556 p += TF;
557
558
559 DEBUG(cout << "Bend3: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
560
561 /*
562 //New Corrections for NPAE and CA:
563 TVector3 v(1.,1.,1.); // Vector in cartesian coordinates
564
565 //Set Azimuth and Elevation
566 v.SetPhi(p.Az());
567 v.SetTheta(TMath::Pi()/2-p.Alt());
568 //Rotation Vectors:
569 TVector3 vNpae( cos(p.Az()), sin(p.Az()), 0);
570 TVector3 vCa( -cos(p.Az())*cos(p.Alt()), -sin(p.Az())*cos(p.Alt()), sin(p.Alt()));
571 //Rotate around the vectors vNpae and vCa
572 v.Rotate(fNpae, vNpae);
573 v.Rotate(fCa, vCa);
574
575 p.Az(v.Phi());
576 p.Alt(TMath::Pi()/2-v.Theta());
577 */
578
579 //Old correction terms for Npae and Ca:
580 const AltAz CA(0, -fCa/cos(p.Alt()));
581 p += CA;
582
583 const AltAz NPAE(0, -fNpae*tan(p.Alt()));
584 p += NPAE;
585
586 DEBUG(cout << "Bend2: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
587
588 const AltAz ANAW(CalcAnAw(p, -1));
589 p += ANAW;
590
591 /* Old correction terms for An and Aw:
592 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
593 const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
594 p += AW;
595 p += AN;
596 */
597
598 DEBUG(cout << "Bend1: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
599
600 const AltAz FLOP(Sign(fFlop, p.Alt()), 0);
601 p += FLOP;
602
603 const AltAz MAGIC1(fMagic1*TMath::Sign(1., sin(p.Az())), 0);
604 p += MAGIC1;
605
606 const AltAz I(fIe, fIa);
607 p += I;
608
609 DEBUG(cout << "Bend0: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
610
611 return p;
612}
613
614AltAz MPointing::CorrectBack(const AltAz &aa) const
615{
616 // Correct [rad]
617 // zdaz [rad]
618 AltAz p = aa;
619
620 DEBUG(cout << setprecision(16));
621 DEBUG(cout << "Back0: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
622
623 const AltAz I(fIe, fIa);
624 p -= I;
625
626 const AltAz MAGIC1(fMagic1*TMath::Sign(1., sin(p.Az())), 0);
627 p -= MAGIC1;
628
629 //const AltAz MAGIC1(fMagic1*sin(p.Az()), 0);
630 //p -= MAGIC1;
631
632 const AltAz FLOP(Sign(fFlop, p.Alt()), 0);
633 p -= FLOP;
634
635 DEBUG(cout << "Back1: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
636
637 /* Old correction terms for An and Aw:
638 const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
639 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
640 p -= AN;
641 p -= AW;
642 */
643
644 const AltAz ANAW(CalcAnAw(p, -1));
645 p -= ANAW;
646
647 DEBUG(cout << "Back2: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
648
649 //Old Correction terms for Npae and Ca:
650 const AltAz NPAE(0, -fNpae*tan(p.Alt()));
651 p -= NPAE;
652
653 const AltAz CA(0, -fCa/cos(p.Alt()));
654 p -= CA;
655
656 /*
657 //New Correction term for Npae and Ca:
658 TVector3 v2(1.,1.,1.); // Vector in cartesian coordinates
659 //Set Azimuth and Elevation
660 v2.SetPhi(p.Az());
661 v2.SetTheta(TMath::Pi()/2-p.Alt());
662 //Rotation Vectors:
663 TVector3 vNpae( cos(p.Az()), sin(p.Az()), 0);
664 TVector3 vCa( -cos(p.Az())*cos(p.Alt()), -sin(p.Az())*cos(p.Alt()), sin(p.Alt()));
665 //Rotate around the vectors vCa and vNpae
666 v2.Rotate(-fCa, vCa);
667 v2.Rotate(-fNpae, vNpae);
668
669 p.Az(v2.Phi());
670 p.Alt(TMath::Pi()/2-v2.Theta());
671 */
672
673 DEBUG(cout << "Back3: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
674
675 const AltAz TF(Sign(fTf*cos(p.Alt()), p.Alt()), 0);
676 const AltAz TX(Sign(fTx/tan(p.Alt()), p.Alt()), 0);
677 p -= TF;
678 //p -= TX;
679
680 DEBUG(cout << "Back4: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
681
682 const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
683 const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
684 p -= CEC;
685 p -= CES;
686
687 DEBUG(cout << "Back5: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
688
689 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
690 const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
691 p -= NRY;
692 p -= NRX;
693
694 DEBUG(cout << "Back6: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
695
696 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
697 const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
698 p -= CRY;
699 p -= CRX;
700
701 DEBUG(cout << "Back7: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
702
703 return p;
704}
705
706ZdAz MPointing::Correct(const ZdAz &zdaz) const
707{
708 // Correct [rad]
709 // zdaz [rad]
710 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
711 AltAz c = Correct(p);
712 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
713}
714
715TVector3 MPointing::Correct(const TVector3 &v) const
716{
717 // Correct [rad]
718 // zdaz [rad]
719 AltAz p(TMath::Pi()/2-v.Theta(), v.Phi());
720 AltAz c = Correct(p);
721 TVector3 rc;
722 rc.SetMagThetaPhi(1, TMath::Pi()/2-c.Alt(), c.Az());
723 return rc;
724}
725
726ZdAz MPointing::CorrectBack(const ZdAz &zdaz) const
727{
728 // Correct [rad]
729 // zdaz [rad]
730 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
731 AltAz c = CorrectBack(p);
732 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
733}
734
735TVector3 MPointing::CorrectBack(const TVector3 &v) const
736{
737 // Correct [rad]
738 // zdaz [rad]
739 AltAz p(TMath::Pi()/2-v.Theta(), v.Phi());
740 AltAz c = CorrectBack(p);
741 TVector3 rc;
742 rc.SetMagThetaPhi(1, TMath::Pi()/2-c.Alt(), c.Az());
743 return rc;
744}
745
746void MPointing::SetParameters(const Double_t *par, Int_t n)
747{
748 Clear();
749
750 while (n--)
751 *fCoeff[n] = par[n]/kRad2Deg;
752}
753
754void MPointing::GetParameters(Double_t *par, Int_t n) const
755{
756 while (n--)
757 par[n] = *fCoeff[n]*kRad2Deg;
758}
759
760void MPointing::GetError(TArrayD &par) const
761{
762 par = fError;
763 for (int i=0; i<kNumPar; i++)
764 par[i] *= TMath::RadToDeg();
765}
766
767TVector2 MPointing::GetDxy() const
768{
769 return TVector2(fDx, fDy)*TMath::RadToDeg();
770}
771
772Double_t MPointing::GetPx() const
773{
774 return fPx*TMath::RadToDeg();
775}
776
777Double_t MPointing::GetPy() const
778{
779 return fPy*TMath::RadToDeg();
780}
781
782void MPointing::SetMinuitParameters(TMinuit &m, Int_t n) const
783{
784 if (n<0)
785 n = kNumPar;
786
787 Int_t ierflg = 0;
788
789 while (n--)
790 m.mnparm(n, fNames[n], *fCoeff[n]*kRad2Deg, 1, -360, 360, ierflg);
791}
792
793void MPointing::GetMinuitParameters(TMinuit &m, Int_t n)
794{
795 if (n<0 || n>m.GetNumPars())
796 n = m.GetNumPars();
797
798 while (n--)
799 {
800 m.GetParameter(n, *fCoeff[n], fError[n]);
801 *fCoeff[n] /= kRad2Deg;
802 fError[n] /= kRad2Deg;
803 }
804}
805/*
806void FormatPar(TMinuit &m, Int_t n)
807{
808 Double_t par, err;
809 m.GetParameter(n, par, err);
810
811 int expp = (int)log10(par);
812 int expe = (int)log10(err);
813
814 if (err<2*pow(10, expe))
815 expe--;
816
817 Int_t exp = expe>expp ? expp : expe;
818
819 par = (int)(par/pow(10, exp)) * pow(10, exp);
820 err = (int)(err/pow(10, exp)) * pow(10, exp);
821
822 cout << par << " +- " << err << flush;
823}
824*/
825void MPointing::PrintMinuitParameters(TMinuit &m, Int_t n) const
826{
827 if (n<0)
828 n = m.GetNumPars();
829
830 cout << setprecision(3);
831
832 Double_t par, er;
833
834 while (n--)
835 {
836 m.GetParameter(n, par, er);
837 cout << Form(" %2d %6s: ", n, (const char*)fNames[n]);
838 cout << setw(8) << par << " \xb1 " << setw(6) << er << endl;
839 }
840}
Note: See TracBrowser for help on using the repository browser.