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

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