source: branches/removing_cpp11_features/mpointing/MPointing.cc@ 18481

Last change on this file since 18481 was 9573, checked in by tbretz, 15 years ago
*** empty log message ***
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.