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

Last change on this file since 4556 was 4545, checked in by wittek, 21 years ago
*** empty log message ***
File size: 32.2 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17! 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// MStarLocalCam[MStarLocalCam], MStarLocalCamSource[MStarLocalCam]
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#include "MSrcPosCam.h"
51
52#include "MLog.h"
53#include "MLogManip.h"
54
55#include "MStarLocalCam.h"
56#include "MStarLocalPos.h"
57#include "MSkyCamTrans.h"
58
59ClassImp(MTelAxisFromStars);
60
61using namespace std;
62
63// --------------------------------------------------------------------------
64//
65// Constructor
66//
67MTelAxisFromStars::MTelAxisFromStars(const char *name, const char *title)
68{
69 fName = name ? name : "MTelAxisFromStars";
70 fTitle = title ? title : "Calculate source position from star positions";
71
72 // if scale factor fLambda should NOT be fixed set fFixdScaleFactor to
73 // -1.0; otherwise set it to the value requested
74 fFixedScaleFactor = 1.0;
75
76 // if rotation angle fAlfa should NOT be fixed set fFixdRotationAngle to
77 // -1.0; otherwise set it to the requested value
78 fFixedRotationAngle = 0.0;
79
80 // default type of input is : the result of the correlated Gauss fit
81 // type 0 : result from the weighted average
82 // type 1 : result from the uncorrelated Gauss fit
83 fInputType = 2;
84}
85
86// --------------------------------------------------------------------------
87//
88// Destructor
89//
90MTelAxisFromStars::~MTelAxisFromStars()
91{
92}
93
94// --------------------------------------------------------------------------
95//
96// Set links to containers
97//
98
99Int_t MTelAxisFromStars::PreProcess(MParList *pList)
100{
101
102 fStarLocalCam = (MStarLocalCam*)pList->FindObject("MStarLocalCam", "MStarLocalCam");
103 if (!fStarLocalCam)
104 {
105 *fLog << err << "MTelAxisFromStars::PreProcess; container 'MStarLocalCam' not found... aborting." << endl;
106 return kFALSE;
107 }
108
109
110 fSourceLocalCam = (MStarLocalCam*)pList->FindObject("MSourceLocalCam", "MStarLocalCam");
111 if (!fSourceLocalCam)
112 {
113 *fLog << err << "MTelAxisFromStars::PreProcess; container 'MSourceLocalCam' not found... continue " << endl;
114 }
115
116
117 fSrcPos = (MSrcPosCam*)pList->FindCreateObj(AddSerialNumber("MSrcPosCam"));
118 if (!fSrcPos)
119 {
120 *fLog << err << "MTelAxisFromStars::PreProcess; MSrcPosCam not found... aborting" << endl;
121 return kFALSE;
122 }
123
124 fSkyCamTrans = (MSkyCamTrans*)pList->FindCreateObj(AddSerialNumber("MSkyCamTrans"));
125 if (!fSkyCamTrans)
126 {
127 *fLog << err << "MTelAxisFromStars::PreProcess; MSkyCamTrans not found... aborting" << endl;
128 return kFALSE;
129 }
130
131
132 return kTRUE;
133}
134
135// --------------------------------------------------------------------------
136//
137// Set the type of the input
138//
139// type = 0 calculated star positions (by averaging)
140// type = 1 fitted star positions (by uncorrelated Gauss fit)
141// type = 2 fitted star positions (by correlated Gauss fit)
142//
143void MTelAxisFromStars::SetInputType(Int_t type)
144{
145 *fLog << all << "MTelAxisFromStars::SetInputType; input type is set equal to : "
146 << type << endl;
147
148 fInputType = type;
149}
150
151// --------------------------------------------------------------------------
152//
153// Fix the scale factor fLambda
154//
155//
156void MTelAxisFromStars::FixScaleFactorAt(Double_t lambda)
157{
158 *fLog << all << "MTelAxisFromStars::FixScaleFactorAt; scale factor will be fixed at : "
159 << lambda << endl;
160
161 fFixedScaleFactor = lambda;
162}
163
164
165// --------------------------------------------------------------------------
166//
167// Fix rotation angle fAlfa
168//
169//
170void MTelAxisFromStars::FixRotationAngleAt(Double_t alfa)
171{
172 *fLog << all << "MTelAxisFromStars::FixRotationAngleAt; rotation angle will be fixed at : "
173 << alfa << endl;
174
175 fFixedRotationAngle = alfa; // [degrees]
176}
177
178
179// --------------------------------------------------------------------------
180//
181// Process
182//
183// call FindSkyCamTrans to find the Sky-Camera transformation
184// call TransSkyCam to transform some sky directions
185// into the camera system
186// put the estimated source position into MSrcPosCam
187//
188
189Int_t MTelAxisFromStars::Process()
190{
191 //Int_t run = fRun->GetRunNumber();
192 //*fLog << "MTelAxisFromStars::Process; run = " << run << endl;
193
194 //--------------------------------------
195 // Define the input for FindSkyCamTrans
196 //
197
198 // get the expected (axy[0], axy[1]) and the measured positions
199 // (bxy[0], bxy[1]) of stars in the camera from MStarLocalCam
200 Int_t fNumStars = fStarLocalCam->GetNumStars();
201
202 if (fNumStars <= 0)
203 return kTRUE;
204
205 TArrayD axy[2];
206 axy[0].Set(fNumStars);
207 axy[1].Set(fNumStars);
208
209 TArrayD bxy[2];
210 bxy[0].Set(fNumStars);
211 bxy[1].Set(fNumStars);
212
213 // error matrix of bxy
214 TArrayD exy[2][2];
215 exy[0][0].Set(fNumStars);
216 exy[0][1].Set(fNumStars);
217 exy[1][0].Set(fNumStars);
218 exy[1][1].Set(fNumStars);
219
220 // transformation parameters
221 Double_t fLambda;
222 Double_t fAlfa;
223 Double_t fA[2][2];
224 Double_t fD[2];
225 Double_t fErrD[2][2];
226 Int_t fNumIter;
227 Int_t fNdof;
228 Double_t fChi2;
229 Double_t fChi2Prob;
230
231 MStarLocalPos *star = 0;
232 TIter next(fStarLocalCam->GetList());
233 Int_t ix = 0;
234
235 // loop over all stars
236 while ( (star = (MStarLocalPos*)next()) )
237 {
238 axy[0][ix] = star->GetXExp();
239 axy[1][ix] = star->GetYExp();
240
241 if (fInputType == 0)
242 {
243 // values from averaging
244 bxy[0][ix] = star->GetMeanXCalc();
245 bxy[1][ix] = star->GetMeanYCalc();
246 exy[0][0][ix]= star->GetSigmaMinorAxisCalc()*star->GetSigmaMinorAxisCalc();
247 exy[0][1][ix] = 0.0;
248 exy[1][0][ix] = 0.0;
249 exy[1][1][ix]= star->GetSigmaMajorAxisCalc()*star->GetSigmaMajorAxisCalc();
250 }
251
252 else if (fInputType == 1)
253 {
254 // values from uncorrelatd Gauss fit
255 bxy[0][ix] = star->GetMeanXFit();
256 bxy[1][ix] = star->GetMeanYFit();
257 exy[0][0][ix]= star->GetSigmaMinorAxisFit()*star->GetSigmaMinorAxisFit();
258 exy[0][1][ix] = 0.0;
259 exy[1][0][ix] = 0.0;
260 exy[1][1][ix]= star->GetSigmaMajorAxisFit()*star->GetSigmaMajorAxisFit();
261 }
262
263 else if (fInputType == 2)
264 {
265 // values from correlatd Gauss fit
266 bxy[0][ix] = star->GetMeanXCGFit();
267 bxy[1][ix] = star->GetMeanYCGFit();
268
269 // this is the error matrix for (MeanXCGFit, MeanYCGFit);
270 // this is the error matrix which should be used
271 exy[0][0][ix] = star->GetXXErrCGFit();
272 exy[0][1][ix] = star->GetXYErrCGFit();
273 exy[1][0][ix] = star->GetXYErrCGFit();
274 exy[1][1][ix] = star->GetYYErrCGFit();
275
276 // this is the error matrix constructed from SigmaXCGFit and SigmaYCGFit;
277 // it is used because the errors above are too small, at present
278 //exy[0][0][ix] = star->GetSigmaXCGFit() * star->GetSigmaXCGFit();
279 //exy[0][1][ix] = star->GetCorrXYCGFit() *
280 // star->GetSigmaXCGFit() * star->GetSigmaYCGFit();
281 //exy[1][0][ix] = exy[0][1][ix];
282 //exy[1][1][ix] = star->GetSigmaYCGFit() * star->GetSigmaYCGFit();
283 }
284
285 else
286 {
287 *fLog << err << "MTelAxisFromStars::Process; type of input is not defined"
288 << endl;
289 return kFALSE;
290 }
291
292 // don't include stars with undefined error
293 Double_t deter = exy[0][0][ix]*exy[1][1][ix]
294 - exy[0][1][ix]*exy[1][0][ix];
295
296 //*fLog << "ix ,deter, xx, xy, yy = " << ix << ": "
297 // << deter << ", " << exy[0][0][ix] << ", "
298 // << exy[0][1][ix] << ", " << exy[1][1][ix] << endl;
299 if (deter <= 0.0)
300 continue;
301
302 //*fLog << "MTelAxisFromStars : " << endl;
303 //*fLog << " ix, XExp, YExp, XFit, YFit, SigmaX2, SigmaXY, SigmaY2 = "
304 // << ix << " : "
305 // << axy[0][ix] << ", " << axy[1][ix] << ", "
306 // << bxy[0][ix] << ", " << bxy[1][ix] << ", "
307 // << exy[0][0][ix] << ", " << exy[0][1][ix] << ", "
308 // << exy[1][1][ix] << endl;
309
310 ix++;
311 }
312
313 //--------------------------------------
314 // Find the transformation from expected positions (axy[1], axy[2])
315 // to measured positions (bxy[1], bxy[2]) in the camera
316
317 Int_t fNStars = ix;
318
319 if (ix < fNumStars)
320 {
321 // reset the sizes of the arrays
322 Int_t fNStars = ix;
323 axy[0].Set(fNStars);
324 axy[1].Set(fNStars);
325
326 bxy[0].Set(fNStars);
327 bxy[1].Set(fNStars);
328
329 exy[0][0].Set(fNStars);
330 exy[0][1].Set(fNStars);
331 exy[1][0].Set(fNStars);
332 exy[1][1].Set(fNStars);
333 }
334
335 Bool_t fitOK;
336 if (fNStars < 1)
337 {
338 *fLog << "MTelAxisFromStars::Process; no star for MTelAxisFromStars"
339 << endl;
340 fitOK = kFALSE;
341 }
342 else
343 {
344 fitOK = FindSkyCamTrans(axy, bxy, exy,
345 fFixedRotationAngle, fFixedScaleFactor, fLambda,
346 fAlfa , fA, fD, fErrD,
347 fNumIter, fNdof, fChi2, fChi2Prob);
348 }
349
350 if (!fitOK)
351 {
352 *fLog << err
353 << "MTelAxisFromStars::Process; Fit to find transformation from star to camera system failed"
354 << endl;
355
356 if (fNStars > 0)
357 {
358 *fLog << err
359 << " fNumIter, fNdof, fChi2, fChi2Prob = " << fNumIter
360 << ", " << fNdof << ", " << fChi2 << ", " << fChi2Prob << endl;
361 }
362
363 return kTRUE;
364 }
365
366
367 //--------------------------------------
368 // Put the transformation parameters into the MSkyCamTrans container
369
370 fSkyCamTrans->SetParameters(fLambda, fAlfa, fA, fD, fErrD,
371 fNumStars, fNumIter, fNdof, fChi2, fChi2Prob);
372 fSkyCamTrans->SetReadyToSave();
373
374
375 //--------------------------------------
376 // Put the estimated position, obtained by transforming the expected
377 // position (0, 0), into SrcPosCam
378
379 fSrcPos->SetXY(fD[0], fD[1]);
380 fSrcPos->SetReadyToSave();
381
382
383 //--------------------------------------
384 // Apply the transformation to some expected positions (asxy[1], asxy[2])
385 // to obtain estimated positions (bsxy[1], bsxy[2]) in the camera
386 // and their error matrices esxy[2][2]
387
388 // get the expected positions (asxy[1], asxy[2]) from another MStarLocalCam
389 // container (with the name "MSourceLocalCam")
390 Int_t fNumStarsSource = fSourceLocalCam->GetNumStars();
391
392 //*fLog << "MTelAxisFromStars::Process; fNumStarsSource = "
393 // << fNumStarsSource << endl;
394
395 if (fNumStarsSource > 0)
396 {
397 TArrayD asxy[2];
398 asxy[0].Set(fNumStarsSource);
399 asxy[1].Set(fNumStarsSource);
400
401 TArrayD bsxy[2];
402 bsxy[0].Set(fNumStarsSource);
403 bsxy[1].Set(fNumStarsSource);
404
405 TArrayD esxy[2][2];
406 esxy[0][0].Set(fNumStarsSource);
407 esxy[0][1].Set(fNumStarsSource);
408 esxy[1][0].Set(fNumStarsSource);
409 esxy[1][1].Set(fNumStarsSource);
410
411 MStarLocalPos *starSource = 0;
412 TIter nextSource(fSourceLocalCam->GetList());
413 ix = 0;
414 while ( (starSource = (MStarLocalPos*)nextSource()) )
415 {
416 asxy[0][ix] = starSource->GetXExp();
417 asxy[1][ix] = starSource->GetYExp();
418
419 ix++;
420 }
421
422 TransSkyCam(fLambda, fA, fD, fErrD, asxy, bsxy, esxy);
423
424 // put the estimated positions into the MStarLocalPos container
425 TIter setnextSource(fSourceLocalCam->GetList());
426 ix = 0;
427 while ( (starSource = (MStarLocalPos*)setnextSource()) )
428 {
429 if (fInputType == 1)
430 {
431 starSource->SetFitValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix],
432 sqrt(esxy[0][0][ix]), sqrt(esxy[1][1][ix]), fChi2, fNdof);
433 }
434 else if (fInputType == 2)
435 {
436 Double_t corr = esxy[0][1][ix]/ sqrt( esxy[0][0][ix] * esxy[1][1][ix] );
437 starSource->SetCGFitValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix],
438 sqrt(esxy[0][0][ix]), sqrt(esxy[1][1][ix]), corr,
439 esxy[0][0][ix], esxy[0][1][ix], esxy[1][1][ix],
440 fChi2, fNdof);
441 }
442
443 ix++;
444 }
445
446 }
447
448 //--------------------------------------
449
450 return kTRUE;
451}
452
453//---------------------------------------------------------------------------
454//
455// FindSkyCamTrans
456//
457// This routine determines the transformation
458//
459// ( cos(alfa) -sin(alfa) )
460// b = lambda * A * a + d A = ( )
461// ^ ^ ^ ( sin(alfa) cos(alfa) )
462// | | |
463// scale rotation shift
464// factor matrix
465//
466// from sky coordinates 'a' (projected onto the camera) to camera
467// coordinates 'b', using the positions of known stars in the camera.
468// The latter positions may have been determined by analysing the
469// DC currents in the different pixels.
470//
471// Input : a[2] x and y coordinates of stars projected onto the camera;
472// they were obtained from (RA, dec) of the stars and
473// (ThetaTel, PhiTel) and the time of observation;
474// these are the 'expected positions' of stars in the camera
475// b[2] 'measured positions' of these stars in the camera;
476// they may have been obtained from the DC currents
477// e[2][2] error matrix of b[2]
478// fixedrotationangle value [in degrees] at which rotation angle
479// alfa should be fixed; -1 means don't fix
480// fixedscalefactor value at which scale factor lambda
481// should be fixed; -1 means don't fix
482//
483// Output : lambda, alfadeg, A[2][2], d[2] fit results;
484// parameters describing the transformation
485// from 'expected positions' to the 'measured
486// positions' in the camera
487// errd[2][2] error matrix of d[2]
488// fNumIter number of iterations
489// fNdoF number of degrees of freedom
490// fChi2 chi-square value
491// fChi2Prob chi-square probability
492//
493// The units are assumed to be
494// [degrees] for alfadeg
495// [mm] for a, b, d
496// [1] for lambda
497
498Bool_t MTelAxisFromStars::FindSkyCamTrans(
499 TArrayD a[2], TArrayD b[2], TArrayD e[2][2],
500 Double_t &fixedrotationang, Double_t &fixedscalefac, Double_t &lambda,
501 Double_t &alfadeg, Double_t A[2][2], Double_t d[2], Double_t errd[2][2],
502 Int_t &fNumIter, Int_t &fNdof, Double_t &fChi2, Double_t &fChi2Prob)
503{
504 Int_t fNumStars = a[0].GetSize();
505
506 //*fLog << "MTelAxisFromStars::FindSkyCamTrans; expected and measured positions :"
507 // << endl;
508 for (Int_t ix=0; ix<fNumStars; ix++)
509 {
510 //*fLog << " ix, a[0], a[1], b[0], b[1], errxx, errxy, erryy = "
511 // << ix << " : "
512 // << a[0][ix] << ", " << a[1][ix] << "; "
513 // << b[0][ix] << ", " << b[1][ix] << "; "
514 // << e[0][0][ix] << ", " << e[0][1][ix] << ", "
515 // << e[1][1][ix] << endl;
516 }
517
518
519 //-------------------------------------------
520 // fix some parameters if the number of degrees of freedom is too low
521 // (<= 0.0)
522
523 Double_t fixedscalefactor = fixedscalefac;
524 Double_t fixedrotationangle = fixedrotationang;
525
526 // calculate number of degrees of freedom
527 fNdof = 2 * fNumStars - 4;
528 if (fixedscalefactor != -1.0)
529 fNdof += 1;
530 if (fixedrotationangle != -1.0)
531 fNdof += 1;
532
533 // if there is only 1 star fix both rotation angle and scale factor
534 if (fNumStars == 1)
535 {
536 if (fixedscalefactor == -1.0)
537 {
538 fixedscalefactor = 1.0;
539 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor is fixed at "
540 << fixedscalefactor << endl;
541 }
542 if (fixedrotationangle == -1.0)
543 {
544 fixedrotationangle = 0.0;
545 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; rotation angle is fixed at "
546 << fixedrotationangle << endl;
547 }
548 }
549 // otherwise fix only 1 parameter if possible
550 else if (fNdof < 0)
551 {
552 if (fNdof < -2)
553 {
554 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
555 << fNdof << "; fNumStars = " << fNumStars << endl;
556 return kFALSE;
557 }
558 else if (fNdof == -2)
559 {
560 if (fixedscalefactor == -1.0 && fixedrotationangle == -1.0)
561 {
562 fixedscalefactor = 1.0;
563 fixedrotationangle = 0.0;
564 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor and rotation angle are fixed at "
565 << fixedscalefactor << " and " << fixedrotationangle
566 << " respectively" << endl;
567 }
568 else
569 {
570 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
571 << fNdof << "; fNumStars = " << fNumStars << endl;
572 return kFALSE;
573 }
574 }
575 else if (fNdof == -1)
576 {
577 if (fixedrotationangle == -1.0)
578 {
579 fixedrotationangle = 0.0;
580 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; rotation angle is fixed at "
581 << fixedrotationangle << endl;
582 }
583 else if (fixedscalefactor == -1.0)
584 {
585 fixedscalefactor = 1.0;
586 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor is fixed at "
587 << fixedscalefactor << endl;
588 }
589 else
590 {
591 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
592 << fNdof << "; fNumStars = " << fNumStars<< endl;
593 return kFALSE;
594 }
595 }
596 }
597
598 // recalculate number of degrees of freedom
599 fNdof = 2 * fNumStars - 4;
600 if (fixedscalefactor != -1.0)
601 fNdof += 1;
602 if (fixedrotationangle != -1.0)
603 fNdof += 1;
604
605 if (fNdof < 0)
606 return kFALSE;
607 //-------------------------------------------
608
609
610 // get first approximation of scaling factor
611 if (fixedscalefactor != -1.0)
612 lambda = fixedscalefactor;
613 else
614 lambda = 1.0;
615
616 Double_t lambdaold = lambda;
617 Double_t dlambda = 0.0;
618
619 // get first approximation of rotation angle
620 Double_t alfa = 0.0;
621 if (fixedrotationangle != -1.0)
622 alfa = fixedrotationangle / kRad2Deg;
623
624
625
626 Double_t alfaold = alfa;
627 // maximum allowed change of alfa in 1 iteration step (5 degrees)
628 Double_t dalfamax = 5.0 / kRad2Deg;
629 Double_t dalfa = 0.0;
630
631 Double_t cosal = cos(alfa);
632 Double_t sinal = sin(alfa);
633
634 A[0][0] = cosal;
635 A[0][1] = -sinal;
636 A[1][0] = sinal;
637 A[1][1] = cosal;
638
639
640 Double_t absdold2 = 10000.0;
641 Double_t fChangeofd2 = 10000.0;
642
643
644 TArrayD Aa[2];
645 Aa[0].Set(fNumStars);
646 Aa[1].Set(fNumStars);
647
648
649 Double_t sumEbminlamAa[2];
650
651 TArrayD Ebminlambracd[2];
652 Ebminlambracd[0].Set(fNumStars);
653 Ebminlambracd[1].Set(fNumStars);
654
655 TArrayD EAa[2];
656 EAa[0].Set(fNumStars);
657 EAa[1].Set(fNumStars);
658
659 // invert the error matrices
660 TArrayD c[2][2];
661 c[0][0].Set(fNumStars);
662 c[0][1].Set(fNumStars);
663 c[1][0].Set(fNumStars);
664 c[1][1].Set(fNumStars);
665
666 for (Int_t ix=0; ix<fNumStars; ix++)
667 {
668 Double_t XX = e[0][0][ix];
669 Double_t XY = e[0][1][ix];
670 Double_t YY = e[1][1][ix];
671
672 // get inverse of error matrix
673 Double_t determ = XX*YY - XY*XY;
674 c[0][0][ix] = YY / determ;
675 c[0][1][ix] = -XY / determ;
676 c[1][0][ix] = -XY / determ;
677 c[1][1][ix] = XX / determ;
678 }
679
680
681
682 // calculate sum of inverted error matrices
683 Double_t determsumc;
684 Double_t sumc[2][2];
685 sumc[0][0] = 0.0;
686 sumc[0][1] = 0.0;
687 sumc[1][0] = 0.0;
688 sumc[1][1] = 0.0;
689
690 for (Int_t ix=0; ix<fNumStars; ix++)
691 {
692 sumc[0][0] += c[0][0][ix];
693 sumc[0][1] += c[0][1][ix];
694 sumc[1][0] += c[1][0][ix];
695 sumc[1][1] += c[1][1][ix];
696 }
697 determsumc = sumc[0][0]*sumc[1][1] - sumc[0][1]*sumc[1][0];
698
699 // calculate inverse of sum of inverted error matrices
700 Double_t sumcinv[2][2];
701 sumcinv[0][0] = sumc[1][1] / determsumc;
702 sumcinv[0][1] = -sumc[0][1] / determsumc;
703 sumcinv[1][0] = -sumc[1][0] / determsumc;
704 sumcinv[1][1] = sumc[0][0] / determsumc;
705
706 //*fLog << "sumcinv = " << sumcinv[0][0] << ", " << sumcinv[0][1]
707 // << ", " << sumcinv[1][1] << endl;
708
709
710 // minimize chi2 by iteration ***** start **********************
711
712 // stop iteration when change in |d|*|d| is less than 'told2'
713 // and change in alfa is less than 'toldalfa'
714 // and change in lambda is less than 'toldlambda'
715 // or chi2 is less than 'tolchi2'
716 Double_t told2 = 0.3*0.3; // [mm*mm]; 1/100 of an inner pixel diameter
717 Double_t toldalfa = 0.01 / kRad2Deg; // 0.01 degrees
718 Double_t toldlambda = 0.00006; // uncertainty of 1 mm of distance
719 // between camera and reflector
720 Double_t tolchi2 = 1.e-5;
721
722 Int_t fNumIterMax = 100;
723 fNumIter = 0;
724
725 for (Int_t i=0; i<fNumIterMax; i++)
726 {
727 fNumIter++;
728
729 // get next approximation of d ------------------
730 for (Int_t ix=0; ix<fNumStars; ix++)
731 {
732 Aa[0][ix] = A[0][0] * a[0][ix] + A[0][1]*a[1][ix];
733 Aa[1][ix] = A[1][0] * a[0][ix] + A[1][1]*a[1][ix];
734
735 //*fLog << "ix, Aa = " << ix << " : " << Aa[0][ix] << ", "
736 // << Aa[1][ix] << endl;
737 }
738
739 sumEbminlamAa[0] = 0.0;
740 sumEbminlamAa[1] = 0.0;
741
742 for (Int_t ix=0; ix<fNumStars; ix++)
743 {
744 sumEbminlamAa[0] += c[0][0][ix] * (b[0][ix] - lambda*Aa[0][ix])
745 + c[0][1][ix] * (b[1][ix] - lambda*Aa[1][ix]);
746
747 sumEbminlamAa[1] += c[1][0][ix] * (b[0][ix] - lambda*Aa[0][ix])
748 + c[1][1][ix] * (b[1][ix] - lambda*Aa[1][ix]);
749 }
750
751 //*fLog << "sumEbminlamAa = " << sumEbminlamAa[0] << ", "
752 // << sumEbminlamAa[1] << endl;
753
754 d[0] = sumcinv[0][0] * sumEbminlamAa[0]
755 + sumcinv[0][1] * sumEbminlamAa[1] ;
756
757 d[1] = sumcinv[1][0] * sumEbminlamAa[0]
758 + sumcinv[1][1] * sumEbminlamAa[1] ;
759
760 Double_t absdnew2 = d[0]*d[0] + d[1]*d[1];
761 fChangeofd2 = absdnew2 - absdold2;
762
763 //*fLog << "fNumIter : " << fNumIter
764 // << "; alfa, lambda, d[0], d[1], absdold2, absdnew2 = " << endl;
765 //*fLog << alfa << ", " << lambda << ", " << d[0] << ", " << d[1]
766 // << ", " << absdold2 << ", " << absdnew2 << endl;
767
768
769 if ( fabs(fChangeofd2) < told2 && fabs(dalfa) < toldalfa &&
770 fabs(dlambda) < toldlambda )
771 {
772 //*fLog << "Iteration stopped because of small changes : fChangeofd2, dalfa, dlambda = "
773 // << fChangeofd2 << ", " << dalfa << ", " << dlambda << endl;
774 break;
775 }
776 absdold2 = absdnew2;
777
778 // get next approximation of matrix A ----------------
779 if (fFixedRotationAngle == -1.0)
780 {
781 for (Int_t ix=0; ix<fNumStars; ix++)
782 {
783 Ebminlambracd[0][ix] =
784 c[0][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
785 + c[0][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
786
787 Ebminlambracd[1][ix] =
788 c[1][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
789 + c[1][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
790
791 //*fLog << "ix, Ebminlambracd = " << ix << " : "
792 // << Ebminlambracd[0][ix] << ", "
793 // << Ebminlambracd[1][ix] << endl;
794 }
795
796 // stop iteration if fChi2 is small enough
797 fChi2 = 0.0;
798 for (Int_t ix=0; ix<fNumStars; ix++)
799 {
800 fChi2 += (b[0][ix]-lambda*Aa[0][ix]-d[0] ) * Ebminlambracd[0][ix]
801 + (b[1][ix]-lambda*Aa[1][ix]-d[1] ) * Ebminlambracd[1][ix];
802 }
803 if ( fChi2 < tolchi2 )
804 {
805 //*fLog << "iteration stopped because of small fChi2 : "
806 // << fChi2 << endl;
807 break;
808 }
809
810
811 Double_t dchi2dA[2][2];
812 dchi2dA[0][0] = 0.0;
813 dchi2dA[0][1] = 0.0;
814 dchi2dA[1][0] = 0.0;
815 dchi2dA[1][1] = 0.0;
816
817 for (Int_t ix=0; ix<fNumStars; ix++)
818 {
819 dchi2dA[0][0] += Ebminlambracd[0][ix] * a[0][ix];
820 dchi2dA[0][1] += Ebminlambracd[0][ix] * a[1][ix];
821 dchi2dA[1][0] += Ebminlambracd[1][ix] * a[0][ix];
822 dchi2dA[1][1] += Ebminlambracd[1][ix] * a[1][ix];
823 }
824
825 //*fLog << "dchi2dA = " << dchi2dA[0][0] << ", " << dchi2dA[0][1]
826 // << ", " << dchi2dA[1][0] << ", " << dchi2dA[1][1] << endl;
827
828 // ********* 1st derivative (d chi2) / (d alfa) ************
829 Double_t dchi2dalfa = -2.0*lambda *
830 ( - sinal*(dchi2dA[0][0]+dchi2dA[1][1])
831 + cosal*(dchi2dA[1][0]-dchi2dA[0][1]) );
832
833
834 //Double_t dalfa1st = - fChi2 / dchi2dalfa;
835
836 //*fLog << "fChi2, dchi2dalfa = " << fChi2 << ", "
837 // << dchi2dalfa << endl;
838 //*fLog << "proposed change of alfa using 1st derivative = "
839 // << dalfa1st << endl;
840
841 // ********* 2nd derivative (d2 chi2) / (d alfa2) ******
842 Double_t term1 = 0.0;
843 Double_t term2 = 0.0;
844 Double_t term3 = 0.0;
845 Double_t term4 = 0.0;
846
847 for (Int_t ix=0; ix<fNumStars; ix++)
848 {
849 term1 += a[0][ix]*c[0][0][ix]*a[0][ix] + a[1][ix]*c[1][0][ix]*a[0][ix]
850 + a[0][ix]*c[0][1][ix]*a[1][ix] + a[1][ix]*c[1][1][ix]*a[1][ix];
851
852 term2 += a[0][ix]*c[1][0][ix]*a[0][ix] - a[1][ix]*c[0][0][ix]*a[0][ix]
853 + a[0][ix]*c[1][1][ix]*a[1][ix] - a[1][ix]*c[0][1][ix]*a[1][ix];
854
855 term3 = a[0][ix]*c[0][0][ix]*a[1][ix] + a[1][ix]*c[1][0][ix]*a[1][ix]
856 - a[0][ix]*c[0][1][ix]*a[0][ix] - a[1][ix]*c[1][1][ix]*a[0][ix];
857
858 term4 += a[0][ix]*c[1][0][ix]*a[1][ix] - a[1][ix]*c[0][0][ix]*a[1][ix]
859 - a[0][ix]*c[1][1][ix]*a[0][ix] + a[1][ix]*c[0][1][ix]*a[0][ix];
860 }
861
862 Double_t d2chi2dalfa2 =
863 - 2.0*lambda * ( - cosal*(dchi2dA[0][0]+dchi2dA[1][1])
864 - sinal*(dchi2dA[1][0]-dchi2dA[0][1]) )
865 + 2.0*lambda*lambda * ( sinal*sinal * term1 - sinal*cosal * term2
866 + sinal*cosal * term3 - cosal*cosal * term4);
867
868 // Gauss-Newton step
869 Double_t dalfa2nd = - dchi2dalfa / d2chi2dalfa2;
870
871 //*fLog << "proposed change of alfa using 2st derivative = "
872 // << dalfa2nd << endl;
873
874 //dalfa = dalfa1st;
875 dalfa = dalfa2nd;
876
877 // ******************************************
878
879
880 // restrict change of alfa
881 if ( fabs(dalfa) > dalfamax )
882 {
883 dalfa = TMath::Sign( dalfamax, dalfa );
884 }
885 alfa = alfaold + dalfa;
886
887 if ( alfa < -5.0/kRad2Deg )
888 alfa = -5.0/kRad2Deg;
889 else if ( alfa > 5.0/kRad2Deg )
890 alfa = 5.0/kRad2Deg;
891
892 dalfa = alfa - alfaold;
893
894 alfaold = alfa;
895
896 sinal = sin(alfa);
897 cosal = cos(alfa);
898
899 A[0][0] = cosal;
900 A[0][1] = -sinal;
901 A[1][0] = sinal;
902 A[1][1] = cosal;
903
904 //*fLog << "alfa-alfaold = " << dalfa << endl;
905 //*fLog << "new alfa = " << alfa << endl;
906 }
907
908
909 // get next approximation of lambda ----------------
910 if (fFixedScaleFactor == -1.0)
911 {
912 for (Int_t ix=0; ix<fNumStars; ix++)
913 {
914 Aa[0][ix] = A[0][0]*a[0][ix] + A[0][1]*a[1][ix];
915 Aa[1][ix] = A[1][0]*a[0][ix] + A[1][1]*a[1][ix];
916
917 EAa[0][ix] =
918 c[0][0][ix] * Aa[0][ix] + c[0][1][ix] * Aa[1][ix];
919 EAa[1][ix] =
920 c[1][0][ix] * Aa[0][ix] + c[1][1][ix] * Aa[1][ix];
921
922 //*fLog << "ix, Aa = " << ix << " : " << Aa[0][ix] << ", "
923 // << Aa[1][ix] << endl;
924
925 //*fLog << "ix, EAa = " << ix << " : " << EAa[0][ix] << ", "
926 // << EAa[1][ix] << endl;
927 }
928
929 Double_t num = 0.0;
930 Double_t denom = 0.0;
931 for (Int_t ix=0; ix<fNumStars; ix++)
932 {
933 num += (b[0][ix]-d[0]) * EAa[0][ix]
934 + (b[1][ix]-d[1]) * EAa[1][ix];
935
936 denom += Aa[0][ix] * EAa[0][ix]
937 + Aa[1][ix] * EAa[1][ix];
938
939 //*fLog << "ix : b-d = " << ix << " : " << b[0][ix]-d[0]
940 // << ", " << b[1][ix]-d[1] << endl;
941
942 //*fLog << "ix : Aa = " << ix << " : " << Aa[0][ix]
943 // << ", " << Aa[1][ix] << endl;
944 }
945
946 lambda = num / denom;
947
948 if ( lambda < 0.9 )
949 lambda = 0.9;
950 else if ( lambda > 1.1 )
951 lambda = 1.1;
952
953 dlambda = lambda - lambdaold;
954 lambdaold = lambda;
955
956 //*fLog << "num, denom, lambda, dlambda = " << num
957 // << ", " << denom << ", " << lambda << ", "
958 // << dlambda << endl;
959 }
960
961 }
962 //------- end of iteration *****************************************
963
964 alfadeg = alfa * kRad2Deg;
965
966 // calculate error matrix of d[2]
967 errd[0][0] = sumcinv[0][0];
968 errd[0][1] = sumcinv[0][1];
969 errd[1][0] = sumcinv[1][0];
970 errd[1][1] = sumcinv[1][1];
971
972 // evaluate quality of fit
973
974 // calculate chi2
975 for (Int_t ix=0; ix<fNumStars; ix++)
976 {
977 Ebminlambracd[0][ix] =
978 c[0][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
979 + c[0][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
980
981 Ebminlambracd[1][ix] =
982 c[1][0][ix] * (b[0][ix] - lambda*Aa[0][ix] - d[0] )
983 + c[1][1][ix] * (b[1][ix] - lambda*Aa[1][ix] - d[1] );
984 }
985
986 fChi2 = 0.0;
987 for (Int_t ix=0; ix<fNumStars; ix++)
988 {
989 fChi2 += (b[0][ix] - lambda*Aa[0][ix] - d[0] ) * Ebminlambracd[0][ix]
990 + (b[1][ix] - lambda*Aa[1][ix] - d[1] ) * Ebminlambracd[1][ix];
991 }
992
993 fChi2Prob = TMath::Prob(fChi2, fNdof);
994
995 *fLog << "MTelAxisFromStars::FindSkyCamTrans :" << endl;
996 *fLog << " fNumStars, fChi2, fNdof, fChi2Prob, fNumIter, fChangeofd2, dalfa, dlambda = "
997 << fNumStars << ", " << fChi2 << ", " << fNdof << ", "
998 << fChi2Prob << ", "
999 << fNumIter << ", " << fChangeofd2 << ", " << dalfa << ", "
1000 << dlambda << endl;
1001 *fLog << " lambda, alfadeg, d[0], d[1] = " << lambda << ", "
1002 << alfadeg << ", " << d[0] << ", " << d[1] << endl;
1003
1004 return kTRUE;
1005}
1006
1007// --------------------------------------------------------------------------
1008//
1009// Apply transformation (lambda, A, d)
1010// to the expected positions (a[1], a[2])
1011// to obtain the estimated positions (b[1], b[2])
1012//
1013// e[2][2] is the error matrix of b[2]
1014
1015void MTelAxisFromStars::TransSkyCam(
1016 Double_t &lambda, Double_t A[2][2], Double_t d[2], Double_t errd[2][2],
1017 TArrayD a[2], TArrayD b[2], TArrayD e[2][2])
1018{
1019 Int_t numpos = a[0].GetSize();
1020 if (numpos <= 0)
1021 return;
1022
1023 //*fLog << "MTelAxisFromStars::TransSkyCam; expected and estimated positions :"
1024 // << endl;
1025
1026 for (Int_t ix=0; ix<numpos; ix++)
1027 {
1028 //*fLog << "MTelAxisFromStars; ix = " << ix << endl;
1029
1030 b[0][ix] = lambda * (A[0][0]*a[0][ix] + A[0][1]*a[1][ix]) + d[0];
1031 b[1][ix] = lambda * (A[1][0]*a[0][ix] + A[1][1]*a[1][ix]) + d[1];
1032
1033 e[0][0][ix] = errd[0][0];
1034 e[0][1][ix] = errd[0][1];
1035 e[1][0][ix] = errd[1][0];
1036 e[1][1][ix] = errd[1][1];
1037
1038 //*fLog << " ix, a[0], a[1], b[0], b[1], errxx, errxy, erryy = "
1039 // << ix << " : "
1040 // << a[0][ix] << ", " << a[1][ix] << "; "
1041 // << b[0][ix] << ", " << b[1][ix] << "; "
1042 // << e[0][0][ix] << ", " << e[0][1][ix] << ", "
1043 // << e[1][1][ix] << endl;
1044 }
1045}
1046
1047// --------------------------------------------------------------------------
1048//
1049//
1050Int_t MTelAxisFromStars::PostProcess()
1051{
1052
1053 return kTRUE;
1054}
1055
1056
1057// --------------------------------------------------------------------------
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
Note: See TracBrowser for help on using the repository browser.