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

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