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

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