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

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