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

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