source: trunk/MagicSoft/Mars/mpointing/MPointing.cc@ 8915

Last change on this file since 8915 was 8915, checked in by tbretz, 16 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 += "° \t";
271 val *= TMath::DegToRad();
272
273 // Find parameter
274 Int_t n = -1;
275 for (int i=0; i<kNumPar; i++)
276 if (str==fNames[i])
277 {
278 n = i;
279 *fCoeff[i] = val;
280 break;
281 }
282
283 fin >> val;
284 sout += Form("%9f°", val);
285
286 if (*fCoeff[n]!=0 || val>0)
287 *fLog << sout << endl;
288
289 if (!fin)
290 {
291 *fLog << err << "ERROR - Reading line " << str << endl;
292 return kFALSE;
293 }
294
295 if (n<0)
296 {
297 *fLog << warn << "WARNING - Parameter " << str << " unknown." << endl;
298 continue;
299 }
300
301 // corresponding error
302 fError[n] = val*TMath::DegToRad();
303 }
304 *fLog << endl;
305
306 fName = name;
307
308 return kTRUE;
309}
310
311Bool_t MPointing::Save(const char *name)
312{
313 /*
314 ! MMT 1987 July 8
315 ! T 36 7.3622 41.448 -0.0481
316 ! IA -37.5465 20.80602
317 ! IE -13.9180 1.25217
318 ! NPAE +7.0751 26.44763
319 ! CA -6.9149 32.05358
320 ! AN +0.5053 1.40956
321 ! AW -2.2016 1.37480
322 ! END
323 */
324
325 ofstream fout(name);
326 if (!fout)
327 {
328 cout << "Error: Cannot open file '" << name << "'" << endl;
329 return kFALSE;
330 }
331
332 MTime t;
333 t.Now();
334
335 fout << "MAGIC1 " << t << endl;
336 fout << "S 00 000000 000000 0000000" << endl;
337 fout << setprecision(8);
338 for (int i=0; i<kNumPar; i++)
339 {
340 fout << " " << setw(6) << GetVarName(i) << " ";
341 fout << setw(13) << *fCoeff[i]*kRad2Deg << " ";
342 fout << setw(11) << fError[i]*kRad2Deg << endl;
343 }
344 fout << "END" << endl;
345
346 fName = name;
347
348 return kTRUE;
349}
350
351Double_t MPointing::Sign(Double_t val, Double_t alt)
352{
353 // Some pointing corrections are defined as Delta ZA, which
354 // is (P. Wallace) defined [0,90]deg while Alt is defined
355 // [0,180]deg
356 return (TMath::Pi()/2-alt < 0 ? -val : val);
357}
358
359AltAz MPointing::AddOffsets(const AltAz &aa) const
360{
361 // Correct [rad]
362 // zdaz [rad]
363 AltAz p = aa;
364
365 const AltAz I(fIe, fIa);
366 p += I;
367
368 return p;
369}
370
371ZdAz MPointing::AddOffsets(const ZdAz &zdaz) const
372{
373 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
374
375 AltAz c = AddOffsets(p);
376
377 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
378}
379
380TVector3 MPointing::AddOffsets(const TVector3 &v) const
381{
382 AltAz p(TMath::Pi()/2-v.Theta(), v.Phi());
383 AltAz c = AddOffsets(p);
384
385 TVector3 rc;
386 rc.SetMagThetaPhi(1, TMath::Pi()/2-c.Alt(), c.Az());
387 return rc;
388}
389
390AltAz MPointing::SubtractOffsets(const AltAz &aa) const
391{
392 // Correct [rad]
393 // zdaz [rad]
394 AltAz p = aa;
395
396 const AltAz I(fIe, fIa);
397 p -= I;
398
399 return p;
400}
401
402ZdAz MPointing::SubtractOffsets(const ZdAz &zdaz) const
403{
404 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
405
406 AltAz c = SubtractOffsets(p);
407
408 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
409}
410
411TVector3 MPointing::SubtractOffsets(const TVector3 &v) const
412{
413 AltAz p(TMath::Pi()/2-v.Theta(), v.Phi());
414 AltAz c = SubtractOffsets(p);
415
416 TVector3 rc;
417 rc.SetMagThetaPhi(1, TMath::Pi()/2-c.Alt(), c.Az());
418 return rc;
419}
420
421AltAz MPointing::CalcAnAw(const AltAz &p, Int_t sign) const
422{
423 // Corrections for AN and AW without approximations
424 // as done by Patrick Wallace. The approximation cannot
425 // be used for MAGIC because the correctioon angle
426 // AW (~1.5deg) is not small enough.
427
428 // Vector in cartesian coordinates
429 TVector3 v1;
430
431 // Set Azimuth and Elevation
432 v1.SetMagThetaPhi(1, TMath::Pi()/2-p.Alt(), p.Az());
433
434
435 TVector3 v2(v1);
436// cout << sign << endl;
437
438// cout << "v1: " << v1.Theta()*TMath::RadToDeg() << " " << v1.Phi()*TMath::RadToDeg() << endl;
439
440 // Rotate around the x- and y-axis
441 v1.RotateY(sign*fAn);
442 v1.RotateX(sign*fAw);
443
444// cout << "v1: " << v1.Theta()*TMath::RadToDeg() << " " << v1.Phi()*TMath::RadToDeg() << endl;
445// cout << "v2: " << v2.Theta()*TMath::RadToDeg() << " " << v2.Theta()*TMath::RadToDeg() << endl;
446
447 // cout << "dv: " << (v2.Theta()-v1.Theta())*TMath::RadToDeg() << " " << (v2.Phi()-v1.Phi())*TMath::RadToDeg() << endl;
448
449 Double_t dalt = v1.Theta()-v2.Theta();
450 Double_t daz = v1.Phi() -v2.Phi();
451
452 //cout << dalt*TMath::RadToDeg() << " " << daz*TMath::RadToDeg() << endl;
453
454 if (daz>TMath::Pi())
455 daz -= TMath::TwoPi();
456 if (daz<-TMath::Pi())
457 daz += TMath::TwoPi();
458
459// if (daz>TMath::Pi()/2)
460// {
461// }
462
463 AltAz d(dalt, daz);
464 return d;
465
466 // Calculate Delta Azimuth and Delta Elevation
467 /*
468 AltAz d(TMath::Pi()/2-v1.Theta(), v1.Phi());
469
470 cout << "p : " << p.Alt()*TMath::RadToDeg() << " " << p.Az()*TMath::RadToDeg() << endl;
471 cout << "d : " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
472 d -= p;
473 cout << "d-p: " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
474 d *= sign;
475 cout << "d* : " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
476
477
478 cout << "p2: " << 90-p.Alt()*TMath::RadToDeg() << " " << p.Az()*TMath::RadToDeg() << endl;
479 cout << "d2: " << 90-d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
480
481 Int_t s1 = 90-d.Alt()*TMath::RadToDeg() < 0 ? -1 : 1;
482 Int_t s2 = 90-p.Alt()*TMath::RadToDeg() < 0 ? -1 : 1;
483
484
485 if (s1 != s2)
486 {
487 //90-d.Alt() <-- -90+d.Alt()
488
489 d.Alt(d.Alt()-TMath::Pi());
490 cout << "Alt-" << endl;
491 }
492 cout << "d': " << 90-d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;*/
493 /*
494 // Fix 'direction' of output depending on input vector
495 if (TMath::Pi()/2-sign*p.Alt()<0)
496 {
497 d.Alt(d.Alt()-TMath::Pi());
498 cout << "Alt-" << endl;
499 }
500 //if (TMath::Pi()/2-sign*p.Alt()>TMath::Pi())
501 //{
502 // d.Alt(TMath::Pi()-d.Alt());
503 // cout << "Alt+" << endl;
504 //}
505
506 // Align correction into [-180,180]
507 while (d.Az()>TMath::Pi())
508 {
509 d.Az(d.Az()-TMath::Pi()*2);
510 cout << "Az-" << endl;
511 }
512 while (d.Az()<-TMath::Pi())
513 {
514 d.Az(d.Az()+TMath::Pi()*2);
515 cout << "Az+" << endl;
516 }
517 */
518 return d;
519}
520
521
522AltAz MPointing::Correct(const AltAz &aa) const
523{
524 // Correct [rad]
525 // zdaz [rad]
526 AltAz p = aa;
527
528 DEBUG(cout << setprecision(16));
529 DEBUG(cout << "Bend7: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
530
531 const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
532 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
533 p += CRX;
534 p += CRY;
535
536 DEBUG(cout << "Bend6: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
537
538 const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
539 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
540 p += NRX;
541 p += NRY;
542
543 DEBUG(cout << "Bend5: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
544
545 const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
546 const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
547 p += CES;
548 p += CEC;
549
550 DEBUG(cout << "Bend4: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
551
552 const AltAz TX(Sign(fTx/tan(p.Alt()), p.Alt()), 0);
553 const AltAz TF(Sign(fTf*cos(p.Alt()), p.Alt()), 0);
554 //p += TX;
555 p += TF;
556
557
558 DEBUG(cout << "Bend3: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
559
560 /*
561 //New Corrections for NPAE and CA:
562 TVector3 v(1.,1.,1.); // Vector in cartesian coordinates
563
564 //Set Azimuth and Elevation
565 v.SetPhi(p.Az());
566 v.SetTheta(TMath::Pi()/2-p.Alt());
567 //Rotation Vectors:
568 TVector3 vNpae( cos(p.Az()), sin(p.Az()), 0);
569 TVector3 vCa( -cos(p.Az())*cos(p.Alt()), -sin(p.Az())*cos(p.Alt()), sin(p.Alt()));
570 //Rotate around the vectors vNpae and vCa
571 v.Rotate(fNpae, vNpae);
572 v.Rotate(fCa, vCa);
573
574 p.Az(v.Phi());
575 p.Alt(TMath::Pi()/2-v.Theta());
576 */
577
578 //Old correction terms for Npae and Ca:
579 const AltAz CA(0, -fCa/cos(p.Alt()));
580 p += CA;
581
582 const AltAz NPAE(0, -fNpae*tan(p.Alt()));
583 p += NPAE;
584
585 DEBUG(cout << "Bend2: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
586
587 const AltAz ANAW(CalcAnAw(p, -1));
588 p += ANAW;
589
590 /* Old correction terms for An and Aw:
591 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
592 const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
593 p += AW;
594 p += AN;
595 */
596
597 DEBUG(cout << "Bend1: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
598
599 const AltAz FLOP(Sign(fFlop, p.Alt()), 0);
600 p += FLOP;
601
602 const AltAz MAGIC1(fMagic1*TMath::Sign(1., sin(p.Az())), 0);
603 p += MAGIC1;
604
605 const AltAz I(fIe, fIa);
606 p += I;
607
608 DEBUG(cout << "Bend0: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
609
610 return p;
611}
612
613AltAz MPointing::CorrectBack(const AltAz &aa) const
614{
615 // Correct [rad]
616 // zdaz [rad]
617 AltAz p = aa;
618
619 DEBUG(cout << setprecision(16));
620 DEBUG(cout << "Back0: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
621
622 const AltAz I(fIe, fIa);
623 p -= I;
624
625 const AltAz MAGIC1(fMagic1*TMath::Sign(1., sin(p.Az())), 0);
626 p -= MAGIC1;
627
628 //const AltAz MAGIC1(fMagic1*sin(p.Az()), 0);
629 //p -= MAGIC1;
630
631 const AltAz FLOP(Sign(fFlop, p.Alt()), 0);
632 p -= FLOP;
633
634 DEBUG(cout << "Back1: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
635
636 /* Old correction terms for An and Aw:
637 const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
638 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
639 p -= AN;
640 p -= AW;
641 */
642
643 const AltAz ANAW(CalcAnAw(p, -1));
644 p -= ANAW;
645
646 DEBUG(cout << "Back2: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
647
648 //Old Correction terms for Npae and Ca:
649 const AltAz NPAE(0, -fNpae*tan(p.Alt()));
650 p -= NPAE;
651
652 const AltAz CA(0, -fCa/cos(p.Alt()));
653 p -= CA;
654
655 /*
656 //New Correction term for Npae and Ca:
657 TVector3 v2(1.,1.,1.); // Vector in cartesian coordinates
658 //Set Azimuth and Elevation
659 v2.SetPhi(p.Az());
660 v2.SetTheta(TMath::Pi()/2-p.Alt());
661 //Rotation Vectors:
662 TVector3 vNpae( cos(p.Az()), sin(p.Az()), 0);
663 TVector3 vCa( -cos(p.Az())*cos(p.Alt()), -sin(p.Az())*cos(p.Alt()), sin(p.Alt()));
664 //Rotate around the vectors vCa and vNpae
665 v2.Rotate(-fCa, vCa);
666 v2.Rotate(-fNpae, vNpae);
667
668 p.Az(v2.Phi());
669 p.Alt(TMath::Pi()/2-v2.Theta());
670 */
671
672 DEBUG(cout << "Back3: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
673
674 const AltAz TF(Sign(fTf*cos(p.Alt()), p.Alt()), 0);
675 const AltAz TX(Sign(fTx/tan(p.Alt()), p.Alt()), 0);
676 p -= TF;
677 //p -= TX;
678
679 DEBUG(cout << "Back4: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
680
681 const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
682 const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
683 p -= CEC;
684 p -= CES;
685
686 DEBUG(cout << "Back5: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
687
688 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
689 const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
690 p -= NRY;
691 p -= NRX;
692
693 DEBUG(cout << "Back6: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
694
695 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
696 const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
697 p -= CRY;
698 p -= CRX;
699
700 DEBUG(cout << "Back7: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
701
702 return p;
703}
704
705ZdAz MPointing::Correct(const ZdAz &zdaz) const
706{
707 // Correct [rad]
708 // zdaz [rad]
709 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
710 AltAz c = Correct(p);
711 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
712}
713
714TVector3 MPointing::Correct(const TVector3 &v) const
715{
716 // Correct [rad]
717 // zdaz [rad]
718 AltAz p(TMath::Pi()/2-v.Theta(), v.Phi());
719 AltAz c = Correct(p);
720 TVector3 rc;
721 rc.SetMagThetaPhi(1, TMath::Pi()/2-c.Alt(), c.Az());
722 return rc;
723}
724
725ZdAz MPointing::CorrectBack(const ZdAz &zdaz) const
726{
727 // Correct [rad]
728 // zdaz [rad]
729 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
730 AltAz c = CorrectBack(p);
731 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
732}
733
734TVector3 MPointing::CorrectBack(const TVector3 &v) const
735{
736 // Correct [rad]
737 // zdaz [rad]
738 AltAz p(TMath::Pi()/2-v.Theta(), v.Phi());
739 AltAz c = CorrectBack(p);
740 TVector3 rc;
741 rc.SetMagThetaPhi(1, TMath::Pi()/2-c.Alt(), c.Az());
742 return rc;
743}
744
745void MPointing::SetParameters(const Double_t *par, Int_t n)
746{
747 Clear();
748
749 while (n--)
750 *fCoeff[n] = par[n]/kRad2Deg;
751}
752
753void MPointing::GetParameters(Double_t *par, Int_t n) const
754{
755 while (n--)
756 par[n] = *fCoeff[n]*kRad2Deg;
757}
758
759void MPointing::GetError(TArrayD &par) const
760{
761 par = fError;
762 for (int i=0; i<kNumPar; i++)
763 par[i] *= TMath::RadToDeg();
764}
765
766TVector2 MPointing::GetDxy() const
767{
768 return TVector2(fDx, fDy)*TMath::RadToDeg();
769}
770
771Double_t MPointing::GetPx() const
772{
773 return fPx*TMath::RadToDeg();
774}
775
776Double_t MPointing::GetPy() const
777{
778 return fPy*TMath::RadToDeg();
779}
780
781void MPointing::SetMinuitParameters(TMinuit &m, Int_t n) const
782{
783 if (n<0)
784 n = kNumPar;
785
786 Int_t ierflg = 0;
787
788 while (n--)
789 m.mnparm(n, fNames[n], *fCoeff[n]*kRad2Deg, 1, -360, 360, ierflg);
790}
791
792void MPointing::GetMinuitParameters(TMinuit &m, Int_t n)
793{
794 if (n<0 || n>m.GetNumPars())
795 n = m.GetNumPars();
796
797 while (n--)
798 {
799 m.GetParameter(n, *fCoeff[n], fError[n]);
800 *fCoeff[n] /= kRad2Deg;
801 fError[n] /= kRad2Deg;
802 }
803}
804/*
805void FormatPar(TMinuit &m, Int_t n)
806{
807 Double_t par, err;
808 m.GetParameter(n, par, err);
809
810 int expp = (int)log10(par);
811 int expe = (int)log10(err);
812
813 if (err<2*pow(10, expe))
814 expe--;
815
816 Int_t exp = expe>expp ? expp : expe;
817
818 par = (int)(par/pow(10, exp)) * pow(10, exp);
819 err = (int)(err/pow(10, exp)) * pow(10, exp);
820
821 cout << par << " +- " << err << flush;
822}
823*/
824void MPointing::PrintMinuitParameters(TMinuit &m, Int_t n) const
825{
826 if (n<0)
827 n = m.GetNumPars();
828
829 cout << setprecision(3);
830
831 Double_t par, err;
832
833 while (n--)
834 {
835 m.GetParameter(n, par, err);
836 cout << Form(" %2d %6s: ", n, (const char*)fNames[n]);
837 cout << setw(8) << par << " \xb1 " << setw(6) << err << endl;
838 }
839}
Note: See TracBrowser for help on using the repository browser.