source: trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.cc@ 4742

Last change on this file since 4742 was 4741, checked in by wittek, 20 years ago
*** empty log message ***
File size: 38.9 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! Author(s): Wolfgang Wittek 07/2004 <mailto:wittek@mppmu.mpg.de>
18!
19! Copyright: MAGIC Software Development, 2000-2004
20!
21!
22\* ======================================================================== */
23
24/////////////////////////////////////////////////////////////////////////////
25//
26// MTelAxisFromStars
27//
28// This task
29// - determines the transformation from expected positions of stars
30// in the camera to measured positions of these stars in the camera
31// - applies this transformation to expected positions of other objects
32// to obtain the estimated positions of these objects in the camera
33// - puts the estimated positions into the relevant containers
34//
35// Input Containers :
36// MStarCam[MStarCam], MStarCamSource[MStarCam]
37//
38// Output Containers :
39// MSkyCamTrans, MSrcPosCam
40//
41/////////////////////////////////////////////////////////////////////////////
42#include <TList.h>
43#include <TSystem.h>
44
45#include <fstream>
46
47#include "MTelAxisFromStars.h"
48
49#include "MParList.h"
50
51
52#include "MLog.h"
53#include "MLogManip.h"
54
55#include "MReportDrive.h"
56#include "MPointingPos.h"
57#include "MSrcPosCam.h"
58
59#include "MStarCam.h"
60#include "MStarPos.h"
61#include "MSkyCamTrans.h"
62#include "MStarCamTrans.h"
63
64ClassImp(MTelAxisFromStars);
65
66using namespace std;
67
68// --------------------------------------------------------------------------
69//
70// Constructor
71//
72MTelAxisFromStars::MTelAxisFromStars(const char *name, const char *title)
73{
74 fName = name ? name : "MTelAxisFromStars";
75 fTitle = title ? title : "Calculate source position from star positions";
76
77 // if scale factor fLambda should NOT be fixed set fFixdScaleFactor to
78 // -1.0; otherwise set it to the value requested
79 fFixedScaleFactor = 1.0;
80
81 // if rotation angle fAlfa should NOT be fixed set fFixdRotationAngle to
82 // -1.0; otherwise set it to the requested value
83 fFixedRotationAngle = 0.0;
84
85 // default type of input is : the result of the Gauss fit
86 // type 0 : result from the weighted average
87 // type 1 : result from the Gauss fit
88 fInputType = 1;
89
90 // default value of fAberr
91 // the value 1.07 is valid if the expected position (with aberration)
92 // in the camera is calculated as the average of the positions obtained
93 // from the reflections at the individual mirrors
94 fAberr = 1.07;
95
96}
97
98// --------------------------------------------------------------------------
99//
100// Destructor
101//
102MTelAxisFromStars::~MTelAxisFromStars()
103{
104 delete fStarCamTrans;
105}
106
107// --------------------------------------------------------------------------
108//
109// Set links to containers
110//
111
112Int_t MTelAxisFromStars::PreProcess(MParList *pList)
113{
114 fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));
115 if (!fDrive)
116 {
117 *fLog << err << AddSerialNumber("MReportDrive")
118 << " not found... aborting." << endl;
119 return kFALSE;
120 }
121
122
123 fStarCam = (MStarCam*)pList->FindObject("MStarCam", "MStarCam");
124 if (!fStarCam)
125 {
126 *fLog << err << "MTelAxisFromStars::PreProcess; container 'MStarCam' not found... aborting." << endl;
127 return kFALSE;
128 }
129
130
131 fSourceCam = (MStarCam*)pList->FindObject("MSourceCam", "MStarCam");
132 if (!fSourceCam)
133 {
134 *fLog << err << "MTelAxisFromStars::PreProcess; container 'MSourceCam' not found... continue " << endl;
135 }
136
137
138 //-----------------------------------------------------------------
139 fSrcPos = (MSrcPosCam*)pList->FindCreateObj(AddSerialNumber("MSrcPosCam"));
140 if (!fSrcPos)
141 {
142 *fLog << err << "MTelAxisFromStars::PreProcess; MSrcPosCam not found... aborting" << endl;
143 return kFALSE;
144 }
145
146 fPointPos = (MPointingPos*)pList->FindCreateObj(AddSerialNumber("MPointingPos"), "MPointingPos");
147 if (!fPointPos)
148 {
149 *fLog << err << "MTelAxisFromStars::PreProcess; MPointingPos not found... aborting" << endl;
150 return kFALSE;
151 }
152
153 fSourcePos = (MPointingPos*)pList->FindCreateObj(AddSerialNumber("MPointingPos"), "MSourcePos");
154 if (!fSourcePos)
155 {
156 *fLog << err << "MTelAxisFromStars::PreProcess; MSourcePos[MPointingPos] not found... aborting" << endl;
157 return kFALSE;
158 }
159
160 fSkyCamTrans = (MSkyCamTrans*)pList->FindCreateObj(AddSerialNumber("MSkyCamTrans"));
161 if (!fSkyCamTrans)
162 {
163 *fLog << err << "MTelAxisFromStars::PreProcess; MSkyCamTrans not found... aborting" << endl;
164 return kFALSE;
165 }
166
167 //-----------------------------------------------------------------
168 // book an MStarCamTrans object
169 // this is needed when calling one of the member functions of MStarCamTrans
170
171 MGeomCam *geom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
172 if (!geom)
173 {
174 *fLog << err << AddSerialNumber("MGeomCam")
175 << " not found... aborting." << endl;
176 return kFALSE;
177 }
178
179 MObservatory *obs = (MObservatory*)pList->FindObject(AddSerialNumber("MObservatory"));
180 if (!obs)
181 {
182 *fLog << err << AddSerialNumber("MObservatory")
183 << " not found... aborting." << endl;
184 return kFALSE;
185 }
186
187 fStarCamTrans = new MStarCamTrans(*geom, *obs);
188
189 return kTRUE;
190}
191
192// --------------------------------------------------------------------------
193//
194// Set optical aberration factor
195//
196// fAberr is the ratio between
197// the distance from the camera center with optical aberration and
198// the distance from the camera center with an ideal imaging
199//
200// fAberr = r_real/r_ideal
201//
202void MTelAxisFromStars::SetOpticalAberr(Double_t aberr)
203{
204 *fLog << all << "MTelAxisFromStars::SetOpticalAberr; the optical aberration factor is set equal to : "
205 << aberr ;
206
207 fAberr = aberr;
208}
209
210// --------------------------------------------------------------------------
211//
212// Set the type of the input
213//
214// type = 0 calculated star positions (by averaging)
215// type = 1 fitted star positions (by Gauss fit)
216//
217void MTelAxisFromStars::SetInputType(Int_t type)
218{
219 *fLog << all << "MTelAxisFromStars::SetInputType; input type is set equal to : "
220 << type ;
221 if (type == 0)
222 *fLog << " (calculated star positions)" << endl;
223 else
224 *fLog << " (fitted star positions)" << endl;
225
226 fInputType = type;
227}
228
229// --------------------------------------------------------------------------
230//
231// Fix the scale factor fLambda
232//
233//
234void MTelAxisFromStars::FixScaleFactorAt(Double_t lambda)
235{
236 *fLog << all << "MTelAxisFromStars::FixScaleFactorAt; scale factor will be fixed at : "
237 << lambda << endl;
238
239 fFixedScaleFactor = lambda;
240}
241
242
243// --------------------------------------------------------------------------
244//
245// Fix rotation angle fAlfa
246//
247//
248void MTelAxisFromStars::FixRotationAngleAt(Double_t alfa)
249{
250 *fLog << all << "MTelAxisFromStars::FixRotationAngleAt; rotation angle will be fixed at : "
251 << alfa << endl;
252
253 fFixedRotationAngle = alfa; // [degrees]
254}
255
256
257// --------------------------------------------------------------------------
258//
259// Process
260//
261// call FindSkyCamTrans to find the Sky-Camera transformation
262// call TransSkyCam to transform some sky directions
263// into the camera system
264// put the estimated source position into MSrcPosCam
265//
266
267Int_t MTelAxisFromStars::Process()
268{
269 //Int_t run = fRun->GetRunNumber();
270 //*fLog << "MTelAxisFromStars::Process; run = " << run << endl;
271
272 //--------------------------------------
273 // Define the input for FindSkyCamTrans
274 //
275
276 // get the expected (axy[0], axy[1]) and the measured positions
277 // (bxy[0], bxy[1]) of stars in the camera from MStarCam
278 Int_t fNumStars = fStarCam->GetNumStars();
279
280 if (fNumStars <= 0)
281 return kTRUE;
282
283 TArrayD axy[2];
284 axy[0].Set(fNumStars);
285 axy[1].Set(fNumStars);
286
287 TArrayD bxy[2];
288 bxy[0].Set(fNumStars);
289 bxy[1].Set(fNumStars);
290
291 // error matrix of bxy
292 TArrayD exy[2][2];
293 exy[0][0].Set(fNumStars);
294 exy[0][1].Set(fNumStars);
295 exy[1][0].Set(fNumStars);
296 exy[1][1].Set(fNumStars);
297
298 // transformation parameters
299 Double_t fLambda;
300 Double_t fAlfa;
301 Double_t fA[2][2];
302 Double_t fD[2];
303 Double_t fErrD[2][2];
304 Int_t fNumIter;
305 Int_t fNdof;
306 Double_t fChi2;
307 Double_t fChi2Prob;
308
309 MStarPos *star = 0;
310 TIter next(fStarCam->GetList());
311 Int_t ix = 0;
312
313 // loop over all stars
314 while ( (star = (MStarPos*)next()) )
315 {
316 axy[0][ix] = star->GetXExp();
317 axy[1][ix] = star->GetYExp();
318
319 if (fInputType == 0)
320 {
321 // values from averaging
322 bxy[0][ix] = star->GetMeanXCalc();
323 bxy[1][ix] = star->GetMeanYCalc();
324
325 // this is the error matrix for (MeanXCalc, MeanYCalc);
326 // this is the error matrix which should be used
327 exy[0][0][ix] = star->GetXXErrCalc();
328 exy[0][1][ix] = star->GetXYErrCalc();
329 exy[1][0][ix] = star->GetXYErrCalc();
330 exy[1][1][ix] = star->GetYYErrCalc();
331
332 //exy[0][0][ix] = star->GetSigmaXCalc()*star->GetSigmaXCalc();
333 //exy[0][1][ix] = 0.0;
334 //exy[1][0][ix] = 0.0;
335 //exy[1][1][ix] = star->GetSigmaYCalc()*star->GetSigmaYCalc();
336 }
337
338 else if (fInputType == 1)
339 {
340 // values from Gauss fit
341 bxy[0][ix] = star->GetMeanXFit();
342 bxy[1][ix] = star->GetMeanYFit();
343
344 // this is the error matrix for (MeanXFit, MeanYFit);
345 // this is the error matrix which should be used
346 exy[0][0][ix] = star->GetXXErrFit();
347 exy[0][1][ix] = star->GetXYErrFit();
348 exy[1][0][ix] = star->GetXYErrFit();
349 exy[1][1][ix] = star->GetYYErrFit();
350
351 // this is the error matrix constructed from SigmaXFit and SigmaYFit;
352 // it is used because the errors above are too small, at present
353 //exy[0][0][ix] = star->GetSigmaXFit() * star->GetSigmaXFit();
354 //exy[0][1][ix] = star->GetCorrXYFit() *
355 // star->GetSigmaXFit() * star->GetSigmaYFit();
356 //exy[1][0][ix] = exy[0][1][ix];
357 //exy[1][1][ix] = star->GetSigmaYFit() * star->GetSigmaYFit();
358 }
359
360 else
361 {
362 *fLog << err << "MTelAxisFromStars::Process; type of input is not defined"
363 << endl;
364 return kFALSE;
365 }
366
367 // don't include stars with undefined error
368 Double_t deter = exy[0][0][ix]*exy[1][1][ix]
369 - exy[0][1][ix]*exy[1][0][ix];
370
371 //*fLog << "ix ,deter, xx, xy, yy = " << ix << ": "
372 // << deter << ", " << exy[0][0][ix] << ", "
373 // << exy[0][1][ix] << ", " << exy[1][1][ix] << endl;
374 if (deter <= 0.0)
375 continue;
376
377 //*fLog << "MTelAxisFromStars : " << endl;
378 //*fLog << " ix, XExp, YExp, XFit, YFit, SigmaX2, SigmaXY, SigmaY2 = "
379 // << ix << " : "
380 // << axy[0][ix] << ", " << axy[1][ix] << ", "
381 // << bxy[0][ix] << ", " << bxy[1][ix] << ", "
382 // << exy[0][0][ix] << ", " << exy[0][1][ix] << ", "
383 // << exy[1][1][ix] << endl;
384
385 ix++;
386 }
387
388 //--------------------------------------
389 // Find the transformation from expected positions (axy[1], axy[2])
390 // to measured positions (bxy[1], bxy[2]) in the camera
391
392 Int_t fNStars = ix;
393
394 if (ix < fNumStars)
395 {
396 // reset the sizes of the arrays
397 Int_t fNStars = ix;
398 axy[0].Set(fNStars);
399 axy[1].Set(fNStars);
400
401 bxy[0].Set(fNStars);
402 bxy[1].Set(fNStars);
403
404 exy[0][0].Set(fNStars);
405 exy[0][1].Set(fNStars);
406 exy[1][0].Set(fNStars);
407 exy[1][1].Set(fNStars);
408 }
409
410 Bool_t fitOK;
411 if (fNStars < 1)
412 {
413 *fLog << "MTelAxisFromStars::Process; no star for MTelAxisFromStars"
414 << endl;
415 fitOK = kFALSE;
416 }
417 else
418 {
419 fitOK = FindSkyCamTrans(axy, bxy, exy,
420 fFixedRotationAngle, fFixedScaleFactor, fLambda,
421 fAlfa , fA, fD, fErrD,
422 fNumIter, fNdof, fChi2, fChi2Prob);
423 }
424
425 if (!fitOK)
426 {
427 *fLog << err
428 << "MTelAxisFromStars::Process; Fit to find transformation from star to camera system failed"
429 << endl;
430
431 if (fNStars > 0)
432 {
433 *fLog << err
434 << " fNumIter, fNdof, fChi2, fChi2Prob = " << fNumIter
435 << ", " << fNdof << ", " << fChi2 << ", " << fChi2Prob << endl;
436 }
437
438 return kTRUE;
439 }
440
441
442 //--------------------------------------
443 // Put the transformation parameters into the MSkyCamTrans container
444
445 fSkyCamTrans->SetParameters(fLambda, fAlfa, fA, fD, fErrD,
446 fNumStars, fNumIter, fNdof, fChi2, fChi2Prob);
447 fSkyCamTrans->SetReadyToSave();
448
449
450 //--------------------------------------
451 // Put the corrected pointing position into MPointingPos
452 //
453 SetPointingPosition(fStarCamTrans, fDrive, fSkyCamTrans, fPointPos);
454
455
456 //--------------------------------------
457 // Put the estimated position of the source into SrcPosCam
458 //
459 // get the source direction from MReportDrive
460 // Note : this has to be changed for the wobble mode, where the source
461 // isn't in the center of the camera
462 Double_t decsource = fDrive->GetDec();
463 Double_t rasource = fDrive->GetRa();
464 //
465 Double_t radrive = fDrive->GetRa();
466 Double_t hdrive = fDrive->GetHa();
467 Double_t hsource = hdrive + radrive - rasource;
468 fSourcePos->SetSkyPosition(rasource, decsource, hsource);
469
470 SetSourcePosition(fStarCamTrans, fPointPos, fSourcePos, fSrcPos);
471 //fSrcPos->SetXY(fD[0], fD[1]);
472 //fSrcPos->SetReadyToSave();
473 *fLog << "after calling SetSourcePosition : , X, Y ,fD = "
474 << fSrcPos->GetX() << ", " << fSrcPos->GetY() << ", "
475 << fD[0] << ", " << fD[1] << endl;
476
477 //--------------------------------------
478 // Apply the transformation to some expected positions (asxy[1], asxy[2])
479 // to obtain estimated positions (bsxy[1], bsxy[2]) in the camera
480 // and their error matrices esxy[2][2]
481
482 // get the expected positions (asxy[1], asxy[2]) from another MStarCam
483 // container (with the name "MSourceCam")
484 Int_t fNumStarsSource = 0;
485
486 if (fSourceCam)
487 fNumStarsSource = fSourceCam->GetNumStars();
488
489 //*fLog << "MTelAxisFromStars::Process; fNumStarsSource = "
490 // << fNumStarsSource << endl;
491
492 if (fNumStarsSource > 0)
493 {
494 TArrayD asxy[2];
495 asxy[0].Set(fNumStarsSource);
496 asxy[1].Set(fNumStarsSource);
497
498 TArrayD bsxy[2];
499 bsxy[0].Set(fNumStarsSource);
500 bsxy[1].Set(fNumStarsSource);
501
502 TArrayD esxy[2][2];
503 esxy[0][0].Set(fNumStarsSource);
504 esxy[0][1].Set(fNumStarsSource);
505 esxy[1][0].Set(fNumStarsSource);
506 esxy[1][1].Set(fNumStarsSource);
507
508 MStarPos *starSource = 0;
509 TIter nextSource(fSourceCam->GetList());
510 ix = 0;
511 while ( (starSource = (MStarPos*)nextSource()) )
512 {
513 asxy[0][ix] = starSource->GetXExp();
514 asxy[1][ix] = starSource->GetYExp();
515
516 ix++;
517 }
518
519 TransSkyCam(fLambda, fA, fD, fErrD, asxy, bsxy, esxy);
520
521 // put the estimated positions into the MStarCam container
522 // with name "MSourceCam"
523 TIter setnextSource(fSourceCam->GetList());
524 ix = 0;
525 while ( (starSource = (MStarPos*)setnextSource()) )
526 {
527 Double_t corr = esxy[0][1][ix]/ sqrt( esxy[0][0][ix] * esxy[1][1][ix] );
528 starSource->SetFitValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix],
529 sqrt(esxy[0][0][ix]), sqrt(esxy[1][1][ix]), corr,
530 esxy[0][0][ix], esxy[0][1][ix], esxy[1][1][ix],
531 fChi2, fNdof);
532
533 ix++;
534 }
535
536 }
537
538 //--------------------------------------
539
540 return kTRUE;
541}
542
543
544
545//---------------------------------------------------------------------------
546//
547// FindSkyCamTrans
548//
549// This routine determines the transformation
550//
551// ( cos(alfa) -sin(alfa) )
552// b = lambda * A * a + d A = ( )
553// ^ ^ ^ ( sin(alfa) cos(alfa) )
554// | | |
555// scale rotation shift
556// factor matrix
557//
558// from sky coordinates 'a' (projected onto the camera) to camera
559// coordinates 'b', using the positions of known stars in the camera.
560// The latter positions may have been determined by analysing the
561// DC currents in the different pixels.
562//
563// Input : a[2] x and y coordinates of stars projected onto the camera;
564// they were obtained from (RA, dec) of the stars and
565// (ThetaTel, PhiTel) and the time of observation;
566// these are the 'expected positions' of stars in the camera
567// b[2] 'measured positions' of these stars in the camera;
568// they may have been obtained from the DC currents
569// e[2][2] error matrix of b[2]
570// fixedrotationangle value [in degrees] at which rotation angle
571// alfa should be fixed; -1 means don't fix
572// fixedscalefactor value at which scale factor lambda
573// should be fixed; -1 means don't fix
574//
575// Output : lambda, alfadeg, A[2][2], d[2] fit results;
576// parameters describing the transformation
577// from 'expected positions' to the 'measured
578// positions' in the camera
579// errd[2][2] error matrix of d[2]
580// fNumIter number of iterations
581// fNdoF number of degrees of freedom
582// fChi2 chi-square value
583// fChi2Prob chi-square probability
584//
585// The units are assumed to be
586// [degrees] for alfadeg
587// [mm] for a, b, d
588// [1] for lambda
589
590Bool_t MTelAxisFromStars::FindSkyCamTrans(
591 TArrayD a[2], TArrayD b[2], TArrayD e[2][2],
592 Double_t &fixedrotationang, Double_t &fixedscalefac, Double_t &lambda,
593 Double_t &alfadeg, Double_t A[2][2], Double_t d[2], Double_t errd[2][2],
594 Int_t &fNumIter, Int_t &fNdof, Double_t &fChi2, Double_t &fChi2Prob)
595{
596 Int_t fNumStars = a[0].GetSize();
597
598 //*fLog << "MTelAxisFromStars::FindSkyCamTrans; expected and measured positions :"
599 // << endl;
600 for (Int_t ix=0; ix<fNumStars; ix++)
601 {
602 //*fLog << " ix, a[0], a[1], b[0], b[1], errxx, errxy, erryy = "
603 // << ix << " : "
604 // << a[0][ix] << ", " << a[1][ix] << "; "
605 // << b[0][ix] << ", " << b[1][ix] << "; "
606 // << e[0][0][ix] << ", " << e[0][1][ix] << ", "
607 // << e[1][1][ix] << endl;
608 }
609
610
611 //-------------------------------------------
612 // fix some parameters if the number of degrees of freedom is too low
613 // (<= 0.0)
614
615 Double_t fixedscalefactor = fixedscalefac;
616 Double_t fixedrotationangle = fixedrotationang;
617
618 // calculate number of degrees of freedom
619 fNdof = 2 * fNumStars - 4;
620 if (fixedscalefactor != -1.0)
621 fNdof += 1;
622 if (fixedrotationangle != -1.0)
623 fNdof += 1;
624
625 // if there is only 1 star fix both rotation angle and scale factor
626 if (fNumStars == 1)
627 {
628 if (fixedscalefactor == -1.0)
629 {
630 fixedscalefactor = 1.0;
631 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor is fixed at "
632 << fixedscalefactor << endl;
633 }
634 if (fixedrotationangle == -1.0)
635 {
636 fixedrotationangle = 0.0;
637 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; rotation angle is fixed at "
638 << fixedrotationangle << endl;
639 }
640 }
641 // otherwise fix only 1 parameter if possible
642 else if (fNdof < 0)
643 {
644 if (fNdof < -2)
645 {
646 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
647 << fNdof << "; fNumStars = " << fNumStars << endl;
648 return kFALSE;
649 }
650 else if (fNdof == -2)
651 {
652 if (fixedscalefactor == -1.0 && fixedrotationangle == -1.0)
653 {
654 fixedscalefactor = 1.0;
655 fixedrotationangle = 0.0;
656 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor and rotation angle are fixed at "
657 << fixedscalefactor << " and " << fixedrotationangle
658 << " respectively" << endl;
659 }
660 else
661 {
662 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
663 << fNdof << "; fNumStars = " << fNumStars << endl;
664 return kFALSE;
665 }
666 }
667 else if (fNdof == -1)
668 {
669 if (fixedrotationangle == -1.0)
670 {
671 fixedrotationangle = 0.0;
672 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; rotation angle is fixed at "
673 << fixedrotationangle << endl;
674 }
675 else if (fixedscalefactor == -1.0)
676 {
677 fixedscalefactor = 1.0;
678 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor is fixed at "
679 << fixedscalefactor << endl;
680 }
681 else
682 {
683 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
684 << fNdof << "; fNumStars = " << fNumStars<< endl;
685 return kFALSE;
686 }
687 }
688 }
689
690 // recalculate number of degrees of freedom
691 fNdof = 2 * fNumStars - 4;
692 if (fixedscalefactor != -1.0)
693 fNdof += 1;
694 if (fixedrotationangle != -1.0)
695 fNdof += 1;
696
697 if (fNdof < 0)
698 return kFALSE;
699 //-------------------------------------------
700
701
702 // get first approximation of scaling factor
703 if (fixedscalefactor != -1.0)
704 lambda = fixedscalefactor;
705 else
706 lambda = 1.0;
707
708 Double_t lambdaold = lambda;
709 Double_t dlambda = 0.0;
710
711 // get first approximation of rotation angle
712 Double_t alfa = 0.0;
713 if (fixedrotationangle != -1.0)
714 alfa = fixedrotationangle / kRad2Deg;
715
716
717
718 Double_t alfaold = alfa;
719 // maximum allowed change of alfa in 1 iteration step (5 degrees)
720 Double_t dalfamax = 5.0 / kRad2Deg;
721 Double_t dalfa = 0.0;
722
723 Double_t cosal = cos(alfa);
724 Double_t sinal = sin(alfa);
725
726 A[0][0] = cosal;
727 A[0][1] = -sinal;
728 A[1][0] = sinal;
729 A[1][1] = cosal;
730
731
732 Double_t absdold2 = 10000.0;
733 Double_t fChangeofd2 = 10000.0;
734
735
736 TArrayD Aa[2];
737 Aa[0].Set(fNumStars);
738 Aa[1].Set(fNumStars);
739
740
741 Double_t sumEbminlamAa[2];
742
743 TArrayD Ebminlambracd[2];
744 Ebminlambracd[0].Set(fNumStars);
745 Ebminlambracd[1].Set(fNumStars);
746
747 TArrayD EAa[2];
748 EAa[0].Set(fNumStars);
749 EAa[1].Set(fNumStars);
750
751 // invert the error matrices
752 TArrayD c[2][2];
753 c[0][0].Set(fNumStars);
754 c[0][1].Set(fNumStars);
755 c[1][0].Set(fNumStars);
756 c[1][1].Set(fNumStars);
757
758 for (Int_t ix=0; ix<fNumStars; ix++)
759 {
760 Double_t XX = e[0][0][ix];
761 Double_t XY = e[0][1][ix];
762 Double_t YY = e[1][1][ix];
763
764 // get inverse of error matrix
765 Double_t determ = XX*YY - XY*XY;
766 c[0][0][ix] = YY / determ;
767 c[0][1][ix] = -XY / determ;
768 c[1][0][ix] = -XY / determ;
769 c[1][1][ix] = XX / determ;
770 }
771
772
773
774 // calculate sum of inverted error matrices
775 Double_t determsumc;
776 Double_t sumc[2][2];
777 sumc[0][0] = 0.0;
778 sumc[0][1] = 0.0;
779 sumc[1][0] = 0.0;
780 sumc[1][1] = 0.0;
781
782 for (Int_t ix=0; ix<fNumStars; ix++)
783 {
784 sumc[0][0] += c[0][0][ix];
785 sumc[0][1] += c[0][1][ix];
786 sumc[1][0] += c[1][0][ix];
787 sumc[1][1] += c[1][1][ix];
788 }
789 determsumc = sumc[0][0]*sumc[1][1] - sumc[0][1]*sumc[1][0];
790
791 // calculate inverse of sum of inverted error matrices
792 Double_t sumcinv[2][2];
793 sumcinv[0][0] = sumc[1][1] / determsumc;
794 sumcinv[0][1] = -sumc[0][1] / determsumc;
795 sumcinv[1][0] = -sumc[1][0] / determsumc;
796 sumcinv[1][1] = sumc[0][0] / determsumc;
797
798 //*fLog << "sumcinv = " << sumcinv[0][0] << ", " << sumcinv[0][1]
799 // << ", " << sumcinv[1][1] << endl;
800
801
802 // minimize chi2 by iteration ***** start **********************
803
804 // stop iteration when change in |d|*|d| is less than 'told2'
805 // and change in alfa is less than 'toldalfa'
806 // and change in lambda is less than 'toldlambda'
807 // or chi2 is less than 'tolchi2'
808 Double_t told2 = 0.3*0.3; // [mm*mm]; 1/100 of an inner pixel diameter
809 Double_t toldalfa = 0.01 / kRad2Deg; // 0.01 degrees
810 Double_t toldlambda = 0.00006; // uncertainty of 1 mm of distance
811 // between camera and reflector
812 Double_t tolchi2 = 1.e-5;
813
814 Int_t fNumIterMax = 100;
815 fNumIter = 0;
816
817 for (Int_t i=0; i<fNumIterMax; i++)
818 {
819 fNumIter++;
820
821 // get next approximation of d ------------------
822 for (Int_t ix=0; ix<fNumStars; ix++)
823 {
824 Aa[0][ix] = A[0][0] * a[0][ix] + A[0][1]*a[1][ix];
825 Aa[1][ix] = A[1][0] * a[0][ix] + A[1][1]*a[1][ix];
826
827 //*fLog << "ix, Aa = " << ix << " : " << Aa[0][ix] << ", "
828 // << Aa[1][ix] << endl;
829 }
830
831 sumEbminlamAa[0] = 0.0;
832 sumEbminlamAa[1] = 0.0;
833
834 for (Int_t ix=0; ix<fNumStars; ix++)
835 {
836 sumEbminlamAa[0] += c[0][0][ix] * (b[0][ix] - lambda*Aa[0][ix])
837 + c[0][1][ix] * (b[1][ix] - lambda*Aa[1][ix]);
838
839 sumEbminlamAa[1] += c[1][0][ix] * (b[0][ix] - lambda*Aa[0][ix])
840 + c[1][1][ix] * (b[1][ix] - lambda*Aa[1][ix]);
841 }
842
843 //*fLog << "sumEbminlamAa = " << sumEbminlamAa[0] << ", "
844 // << sumEbminlamAa[1] << endl;
845
846 d[0] = sumcinv[0][0] * sumEbminlamAa[0]
847 + sumcinv[0][1] * sumEbminlamAa[1] ;
848
849 d[1] = sumcinv[1][0] * sumEbminlamAa[0]
850 + sumcinv[1][1] * sumEbminlamAa[1] ;
851
852 Double_t absdnew2 = d[0]*d[0] + d[1]*d[1];
853 fChangeofd2 = absdnew2 - absdold2;
854
855 //*fLog << "fNumIter : " << fNumIter
856 // << "; alfa, lambda, d[0], d[1], absdold2, absdnew2 = " << endl;
857 //*fLog << alfa << ", " << lambda << ", " << d[0] << ", " << d[1]
858 // << ", " << absdold2 << ", " << absdnew2 << endl;
859
860
861 if ( fabs(fChangeofd2) < told2 && fabs(dalfa) < toldalfa &&
862 fabs(dlambda) < toldlambda )
863 {
864 //*fLog << "Iteration stopped because of small changes : fChangeofd2, dalfa, dlambda = "
865 // << fChangeofd2 << ", " << dalfa << ", " << dlambda << endl;
866 break;
867 }
868 absdold2 = absdnew2;
869
870 // get next approximation of matrix A ----------------
871 if (fFixedRotationAngle == -1.0)
872 {
873 for (Int_t ix=0; ix<fNumStars; ix++)
874 {
875 Ebminlambracd[0][ix] =
876 c[0][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
877 + c[0][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
878
879 Ebminlambracd[1][ix] =
880 c[1][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
881 + c[1][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
882
883 //*fLog << "ix, Ebminlambracd = " << ix << " : "
884 // << Ebminlambracd[0][ix] << ", "
885 // << Ebminlambracd[1][ix] << endl;
886 }
887
888 // stop iteration if fChi2 is small enough
889 fChi2 = 0.0;
890 for (Int_t ix=0; ix<fNumStars; ix++)
891 {
892 fChi2 += (b[0][ix]-lambda*Aa[0][ix]-d[0] ) * Ebminlambracd[0][ix]
893 + (b[1][ix]-lambda*Aa[1][ix]-d[1] ) * Ebminlambracd[1][ix];
894 }
895 if ( fChi2 < tolchi2 )
896 {
897 //*fLog << "iteration stopped because of small fChi2 : "
898 // << fChi2 << endl;
899 break;
900 }
901
902
903 Double_t dchi2dA[2][2];
904 dchi2dA[0][0] = 0.0;
905 dchi2dA[0][1] = 0.0;
906 dchi2dA[1][0] = 0.0;
907 dchi2dA[1][1] = 0.0;
908
909 for (Int_t ix=0; ix<fNumStars; ix++)
910 {
911 dchi2dA[0][0] += Ebminlambracd[0][ix] * a[0][ix];
912 dchi2dA[0][1] += Ebminlambracd[0][ix] * a[1][ix];
913 dchi2dA[1][0] += Ebminlambracd[1][ix] * a[0][ix];
914 dchi2dA[1][1] += Ebminlambracd[1][ix] * a[1][ix];
915 }
916
917 //*fLog << "dchi2dA = " << dchi2dA[0][0] << ", " << dchi2dA[0][1]
918 // << ", " << dchi2dA[1][0] << ", " << dchi2dA[1][1] << endl;
919
920 // ********* 1st derivative (d chi2) / (d alfa) ************
921 Double_t dchi2dalfa = -2.0*lambda *
922 ( - sinal*(dchi2dA[0][0]+dchi2dA[1][1])
923 + cosal*(dchi2dA[1][0]-dchi2dA[0][1]) );
924
925
926 //Double_t dalfa1st = - fChi2 / dchi2dalfa;
927
928 //*fLog << "fChi2, dchi2dalfa = " << fChi2 << ", "
929 // << dchi2dalfa << endl;
930 //*fLog << "proposed change of alfa using 1st derivative = "
931 // << dalfa1st << endl;
932
933 // ********* 2nd derivative (d2 chi2) / (d alfa2) ******
934 Double_t term1 = 0.0;
935 Double_t term2 = 0.0;
936 Double_t term3 = 0.0;
937 Double_t term4 = 0.0;
938
939 for (Int_t ix=0; ix<fNumStars; ix++)
940 {
941 term1 += a[0][ix]*c[0][0][ix]*a[0][ix] + a[1][ix]*c[1][0][ix]*a[0][ix]
942 + a[0][ix]*c[0][1][ix]*a[1][ix] + a[1][ix]*c[1][1][ix]*a[1][ix];
943
944 term2 += a[0][ix]*c[1][0][ix]*a[0][ix] - a[1][ix]*c[0][0][ix]*a[0][ix]
945 + a[0][ix]*c[1][1][ix]*a[1][ix] - a[1][ix]*c[0][1][ix]*a[1][ix];
946
947 term3 = a[0][ix]*c[0][0][ix]*a[1][ix] + a[1][ix]*c[1][0][ix]*a[1][ix]
948 - a[0][ix]*c[0][1][ix]*a[0][ix] - a[1][ix]*c[1][1][ix]*a[0][ix];
949
950 term4 += a[0][ix]*c[1][0][ix]*a[1][ix] - a[1][ix]*c[0][0][ix]*a[1][ix]
951 - a[0][ix]*c[1][1][ix]*a[0][ix] + a[1][ix]*c[0][1][ix]*a[0][ix];
952 }
953
954 Double_t d2chi2dalfa2 =
955 - 2.0*lambda * ( - cosal*(dchi2dA[0][0]+dchi2dA[1][1])
956 - sinal*(dchi2dA[1][0]-dchi2dA[0][1]) )
957 + 2.0*lambda*lambda * ( sinal*sinal * term1 - sinal*cosal * term2
958 + sinal*cosal * term3 - cosal*cosal * term4);
959
960 // Gauss-Newton step
961 Double_t dalfa2nd = - dchi2dalfa / d2chi2dalfa2;
962
963 //*fLog << "proposed change of alfa using 2st derivative = "
964 // << dalfa2nd << endl;
965
966 //dalfa = dalfa1st;
967 dalfa = dalfa2nd;
968
969 // ******************************************
970
971
972 // restrict change of alfa
973 if ( fabs(dalfa) > dalfamax )
974 {
975 dalfa = TMath::Sign( dalfamax, dalfa );
976 }
977 alfa = alfaold + dalfa;
978
979 if ( alfa < -5.0/kRad2Deg )
980 alfa = -5.0/kRad2Deg;
981 else if ( alfa > 5.0/kRad2Deg )
982 alfa = 5.0/kRad2Deg;
983
984 dalfa = alfa - alfaold;
985
986 alfaold = alfa;
987
988 sinal = sin(alfa);
989 cosal = cos(alfa);
990
991 A[0][0] = cosal;
992 A[0][1] = -sinal;
993 A[1][0] = sinal;
994 A[1][1] = cosal;
995
996 //*fLog << "alfa-alfaold = " << dalfa << endl;
997 //*fLog << "new alfa = " << alfa << endl;
998 }
999
1000
1001 // get next approximation of lambda ----------------
1002 if (fFixedScaleFactor == -1.0)
1003 {
1004 for (Int_t ix=0; ix<fNumStars; ix++)
1005 {
1006 Aa[0][ix] = A[0][0]*a[0][ix] + A[0][1]*a[1][ix];
1007 Aa[1][ix] = A[1][0]*a[0][ix] + A[1][1]*a[1][ix];
1008
1009 EAa[0][ix] =
1010 c[0][0][ix] * Aa[0][ix] + c[0][1][ix] * Aa[1][ix];
1011 EAa[1][ix] =
1012 c[1][0][ix] * Aa[0][ix] + c[1][1][ix] * Aa[1][ix];
1013
1014 //*fLog << "ix, Aa = " << ix << " : " << Aa[0][ix] << ", "
1015 // << Aa[1][ix] << endl;
1016
1017 //*fLog << "ix, EAa = " << ix << " : " << EAa[0][ix] << ", "
1018 // << EAa[1][ix] << endl;
1019 }
1020
1021 Double_t num = 0.0;
1022 Double_t denom = 0.0;
1023 for (Int_t ix=0; ix<fNumStars; ix++)
1024 {
1025 num += (b[0][ix]-d[0]) * EAa[0][ix]
1026 + (b[1][ix]-d[1]) * EAa[1][ix];
1027
1028 denom += Aa[0][ix] * EAa[0][ix]
1029 + Aa[1][ix] * EAa[1][ix];
1030
1031 //*fLog << "ix : b-d = " << ix << " : " << b[0][ix]-d[0]
1032 // << ", " << b[1][ix]-d[1] << endl;
1033
1034 //*fLog << "ix : Aa = " << ix << " : " << Aa[0][ix]
1035 // << ", " << Aa[1][ix] << endl;
1036 }
1037
1038 lambda = num / denom;
1039
1040 if ( lambda < 0.9 )
1041 lambda = 0.9;
1042 else if ( lambda > 1.1 )
1043 lambda = 1.1;
1044
1045 dlambda = lambda - lambdaold;
1046 lambdaold = lambda;
1047
1048 //*fLog << "num, denom, lambda, dlambda = " << num
1049 // << ", " << denom << ", " << lambda << ", "
1050 // << dlambda << endl;
1051 }
1052
1053 }
1054 //------- end of iteration *****************************************
1055
1056 alfadeg = alfa * kRad2Deg;
1057
1058 // calculate error matrix of d[2]
1059 errd[0][0] = sumcinv[0][0];
1060 errd[0][1] = sumcinv[0][1];
1061 errd[1][0] = sumcinv[1][0];
1062 errd[1][1] = sumcinv[1][1];
1063
1064 // evaluate quality of fit
1065
1066 // calculate chi2
1067 for (Int_t ix=0; ix<fNumStars; ix++)
1068 {
1069 Ebminlambracd[0][ix] =
1070 c[0][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
1071 + c[0][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
1072
1073 Ebminlambracd[1][ix] =
1074 c[1][0][ix] * (b[0][ix] - lambda*Aa[0][ix] - d[0] )
1075 + c[1][1][ix] * (b[1][ix] - lambda*Aa[1][ix] - d[1] );
1076 }
1077
1078 fChi2 = 0.0;
1079 for (Int_t ix=0; ix<fNumStars; ix++)
1080 {
1081 fChi2 += (b[0][ix] - lambda*Aa[0][ix] - d[0] ) * Ebminlambracd[0][ix]
1082 + (b[1][ix] - lambda*Aa[1][ix] - d[1] ) * Ebminlambracd[1][ix];
1083 }
1084
1085 fChi2Prob = TMath::Prob(fChi2, fNdof);
1086
1087 *fLog << "MTelAxisFromStars::FindSkyCamTrans :" << endl;
1088 *fLog << " fNumStars, fChi2, fNdof, fChi2Prob, fNumIter, fChangeofd2, dalfa, dlambda = "
1089 << fNumStars << ", " << fChi2 << ", " << fNdof << ", "
1090 << fChi2Prob << ", "
1091 << fNumIter << ", " << fChangeofd2 << ", " << dalfa << ", "
1092 << dlambda << endl;
1093 *fLog << " lambda, alfadeg, d[0], d[1] = " << lambda << ", "
1094 << alfadeg << ", " << d[0] << ", " << d[1] << endl;
1095
1096 return kTRUE;
1097}
1098
1099// --------------------------------------------------------------------------
1100//
1101// Apply transformation (lambda, A, d)
1102// to the expected positions (a[1], a[2])
1103// to obtain the estimated positions (b[1], b[2])
1104//
1105// e[2][2] is the error matrix of b[2]
1106
1107void MTelAxisFromStars::TransSkyCam(
1108 Double_t &lambda, Double_t A[2][2], Double_t d[2], Double_t errd[2][2],
1109 TArrayD a[2], TArrayD b[2], TArrayD e[2][2])
1110{
1111 Int_t numpos = a[0].GetSize();
1112 if (numpos <= 0)
1113 return;
1114
1115 //*fLog << "MTelAxisFromStars::TransSkyCam; expected and estimated positions :"
1116 // << endl;
1117
1118 for (Int_t ix=0; ix<numpos; ix++)
1119 {
1120 //*fLog << "MTelAxisFromStars; ix = " << ix << endl;
1121
1122 b[0][ix] = lambda * (A[0][0]*a[0][ix] + A[0][1]*a[1][ix]) + d[0];
1123 b[1][ix] = lambda * (A[1][0]*a[0][ix] + A[1][1]*a[1][ix]) + d[1];
1124
1125 e[0][0][ix] = errd[0][0];
1126 e[0][1][ix] = errd[0][1];
1127 e[1][0][ix] = errd[1][0];
1128 e[1][1][ix] = errd[1][1];
1129
1130 //*fLog << " ix, a[0], a[1], b[0], b[1], errxx, errxy, erryy = "
1131 // << ix << " : "
1132 // << a[0][ix] << ", " << a[1][ix] << "; "
1133 // << b[0][ix] << ", " << b[1][ix] << "; "
1134 // << e[0][0][ix] << ", " << e[0][1][ix] << ", "
1135 // << e[1][1][ix] << endl;
1136 }
1137}
1138
1139// --------------------------------------------------------------------------
1140//
1141// SetPointingPosition
1142//
1143// put the corrected pointing direction into MPointingPos[MPointingPos];
1144// this direction corresponds to the position (0,0) in the camera
1145//
1146
1147Bool_t MTelAxisFromStars::SetPointingPosition(MStarCamTrans *fstarcamtrans,
1148 MReportDrive *fdrive, MSkyCamTrans *ftrans, MPointingPos *fpointpos)
1149{
1150 Double_t decdrive = fdrive->GetDec();
1151 Double_t hdrive = fdrive->GetHa();
1152 Double_t radrive = fdrive->GetRa();
1153
1154 // this is the estimated position (with optical aberration) in the camera
1155 // corresponding to the direction in MReportDrive
1156 Double_t Xpoint = (ftrans->GetShiftD())[0];
1157 Double_t Ypoint = (ftrans->GetShiftD())[1];
1158
1159 // get the sky direction corresponding to the position (0,0) in the camera
1160 Double_t decpoint = 0.0;
1161 Double_t hpoint = 0.0;
1162 fstarcamtrans->CelCamToCel0(decdrive, hdrive,
1163 Xpoint/fAberr, Ypoint/fAberr, decpoint, hpoint);
1164 Double_t rapoint = radrive - hpoint + hdrive;
1165 fpointpos->SetSkyPosition(rapoint, decpoint, hpoint);
1166
1167 // get the local direction corresponding to the position (0,0) in the camera
1168 Double_t thetadrive = fdrive->GetNominalZd();
1169 Double_t phidrive = fdrive->GetNominalAz();
1170 Double_t thetapoint = 0.0;
1171 Double_t phipoint = 0.0;
1172 fstarcamtrans->LocCamToLoc0(thetadrive, phidrive,
1173 Xpoint/fAberr, Ypoint/fAberr, thetapoint, phipoint);
1174 fpointpos->SetLocalPosition(thetapoint, phipoint);
1175 fpointpos->SetReadyToSave();
1176
1177 *fLog << "SetPointingPosition : decdrive, hdrive, radrive Xpoint, Ypoint = "
1178 << decdrive << ", " << hdrive << ", " << radrive << ", "
1179 << Xpoint << ", " << Ypoint << endl;
1180
1181 *fLog << "SetPointingPosition : thetadrive, phidrive, thetapoint, phipoint = "
1182 << thetadrive << ", " << phidrive << ", " << thetapoint << ", "
1183 << phipoint << endl;
1184
1185 return kTRUE;
1186}
1187
1188// --------------------------------------------------------------------------
1189//
1190// SetSourcePosition
1191//
1192// put the estimated position of the source in the camera into
1193// MSrcPosCam[MSrcPosCam]
1194//
1195// and the estimated local direction of the source into
1196// MSourcePos[MPointingPos]
1197//
1198
1199Bool_t MTelAxisFromStars::SetSourcePosition(MStarCamTrans *fstarcamtrans,
1200 MPointingPos *fpointpos, MPointingPos *fsourcepos, MSrcPosCam *fsrcpos)
1201{
1202 // get the corrected pointing direction
1203 // corresponding to the position (0,0) in the camera
1204 Double_t decpoint = fpointpos->GetDec();
1205 Double_t hpoint = fpointpos->GetHa();
1206
1207 // get the sky direction of the source
1208 Double_t decsource = fsourcepos->GetDec();
1209 Double_t hsource = fsourcepos->GetHa();
1210
1211 // get the estimated position (Xsource, Ysource) of the source in the camera;
1212 // this is a position for an ideal imaging, without optical aberration
1213 Double_t Xsource = 0.0;
1214 Double_t Ysource = 0.0;
1215 fstarcamtrans->Cel0CelToCam(decpoint, hpoint,
1216 decsource, hsource, Xsource, Ysource);
1217 fsrcpos->SetXY(Xsource*fAberr, Ysource*fAberr);
1218 fsrcpos->SetReadyToSave();
1219
1220 // get the estimated local direction of the source
1221 Double_t thetapoint = fpointpos->GetZd();
1222 Double_t phipoint = fpointpos->GetAz();
1223 Double_t thetasource = 0.0;
1224 Double_t phisource = 0.0;
1225 fstarcamtrans->Loc0CamToLoc(thetapoint, phipoint,
1226 Xsource, Ysource, thetasource, phisource);
1227 fsourcepos->SetLocalPosition(thetasource, phisource);
1228 fsourcepos->SetReadyToSave();
1229
1230 *fLog << "SetSourcePosition : decpoint, hpoint, decsource, hsource, Xsource, Ysource = "
1231 << decpoint << ", " << hpoint << ", " << decsource << ", "
1232 << hsource << ", " << Xsource << ", " << Ysource << endl;
1233 *fLog << "SetSourcePosition : thetapoint, phipoint, thetasource, phisource = "
1234 << thetapoint << ", " << phipoint << ", " << thetasource << ", "
1235 << phisource << endl;
1236
1237 return kTRUE;
1238}
1239
1240// --------------------------------------------------------------------------
1241//
1242//
1243Int_t MTelAxisFromStars::PostProcess()
1244{
1245
1246 return kTRUE;
1247}
1248
1249
1250// --------------------------------------------------------------------------
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
Note: See TracBrowser for help on using the repository browser.