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

Last change on this file since 4716 was 4705, checked in by wittek, 22 years ago
*** empty log message ***
File size: 31.8 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#include "MSrcPosCam.h"
51
52#include "MLog.h"
53#include "MLogManip.h"
54
55#include "MStarCam.h"
56#include "MStarPos.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 Gauss fit
81 // type 0 : result from the weighted average
82 // type 1 : result from the Gauss fit
83 fInputType = 1;
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 fStarCam = (MStarCam*)pList->FindObject("MStarCam", "MStarCam");
103 if (!fStarCam)
104 {
105 *fLog << err << "MTelAxisFromStars::PreProcess; container 'MStarCam' not found... aborting." << endl;
106 return kFALSE;
107 }
108
109
110 fSourceCam = (MStarCam*)pList->FindObject("MSourceCam", "MStarCam");
111 if (!fSourceCam)
112 {
113 *fLog << err << "MTelAxisFromStars::PreProcess; container 'MSourceCam' 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 Gauss fit)
141//
142void MTelAxisFromStars::SetInputType(Int_t type)
143{
144 *fLog << all << "MTelAxisFromStars::SetInputType; input type is set equal to : "
145 << type ;
146 if (type == 0)
147 *fLog << " (calculated star positions)" << endl;
148 else
149 *fLog << " (fitted star positions)" << endl;
150
151 fInputType = type;
152}
153
154// --------------------------------------------------------------------------
155//
156// Fix the scale factor fLambda
157//
158//
159void MTelAxisFromStars::FixScaleFactorAt(Double_t lambda)
160{
161 *fLog << all << "MTelAxisFromStars::FixScaleFactorAt; scale factor will be fixed at : "
162 << lambda << endl;
163
164 fFixedScaleFactor = lambda;
165}
166
167
168// --------------------------------------------------------------------------
169//
170// Fix rotation angle fAlfa
171//
172//
173void MTelAxisFromStars::FixRotationAngleAt(Double_t alfa)
174{
175 *fLog << all << "MTelAxisFromStars::FixRotationAngleAt; rotation angle will be fixed at : "
176 << alfa << endl;
177
178 fFixedRotationAngle = alfa; // [degrees]
179}
180
181
182// --------------------------------------------------------------------------
183//
184// Process
185//
186// call FindSkyCamTrans to find the Sky-Camera transformation
187// call TransSkyCam to transform some sky directions
188// into the camera system
189// put the estimated source position into MSrcPosCam
190//
191
192Int_t MTelAxisFromStars::Process()
193{
194 //Int_t run = fRun->GetRunNumber();
195 //*fLog << "MTelAxisFromStars::Process; run = " << run << endl;
196
197 //--------------------------------------
198 // Define the input for FindSkyCamTrans
199 //
200
201 // get the expected (axy[0], axy[1]) and the measured positions
202 // (bxy[0], bxy[1]) of stars in the camera from MStarCam
203 Int_t fNumStars = fStarCam->GetNumStars();
204
205 if (fNumStars <= 0)
206 return kTRUE;
207
208 TArrayD axy[2];
209 axy[0].Set(fNumStars);
210 axy[1].Set(fNumStars);
211
212 TArrayD bxy[2];
213 bxy[0].Set(fNumStars);
214 bxy[1].Set(fNumStars);
215
216 // error matrix of bxy
217 TArrayD exy[2][2];
218 exy[0][0].Set(fNumStars);
219 exy[0][1].Set(fNumStars);
220 exy[1][0].Set(fNumStars);
221 exy[1][1].Set(fNumStars);
222
223 // transformation parameters
224 Double_t fLambda;
225 Double_t fAlfa;
226 Double_t fA[2][2];
227 Double_t fD[2];
228 Double_t fErrD[2][2];
229 Int_t fNumIter;
230 Int_t fNdof;
231 Double_t fChi2;
232 Double_t fChi2Prob;
233
234 MStarPos *star = 0;
235 TIter next(fStarCam->GetList());
236 Int_t ix = 0;
237
238 // loop over all stars
239 while ( (star = (MStarPos*)next()) )
240 {
241 axy[0][ix] = star->GetXExp();
242 axy[1][ix] = star->GetYExp();
243
244 if (fInputType == 0)
245 {
246 // values from averaging
247 bxy[0][ix] = star->GetMeanXCalc();
248 bxy[1][ix] = star->GetMeanYCalc();
249
250 // this is the error matrix for (MeanXCalc, MeanYCalc);
251 // this is the error matrix which should be used
252 exy[0][0][ix] = star->GetXXErrCalc();
253 exy[0][1][ix] = star->GetXYErrCalc();
254 exy[1][0][ix] = star->GetXYErrCalc();
255 exy[1][1][ix] = star->GetYYErrCalc();
256
257 //exy[0][0][ix] = star->GetSigmaXCalc()*star->GetSigmaXCalc();
258 //exy[0][1][ix] = 0.0;
259 //exy[1][0][ix] = 0.0;
260 //exy[1][1][ix] = star->GetSigmaYCalc()*star->GetSigmaYCalc();
261 }
262
263 else if (fInputType == 1)
264 {
265 // values from Gauss fit
266 bxy[0][ix] = star->GetMeanXFit();
267 bxy[1][ix] = star->GetMeanYFit();
268
269 // this is the error matrix for (MeanXFit, MeanYFit);
270 // this is the error matrix which should be used
271 exy[0][0][ix] = star->GetXXErrFit();
272 exy[0][1][ix] = star->GetXYErrFit();
273 exy[1][0][ix] = star->GetXYErrFit();
274 exy[1][1][ix] = star->GetYYErrFit();
275
276 // this is the error matrix constructed from SigmaXFit and SigmaYFit;
277 // it is used because the errors above are too small, at present
278 //exy[0][0][ix] = star->GetSigmaXFit() * star->GetSigmaXFit();
279 //exy[0][1][ix] = star->GetCorrXYFit() *
280 // star->GetSigmaXFit() * star->GetSigmaYFit();
281 //exy[1][0][ix] = exy[0][1][ix];
282 //exy[1][1][ix] = star->GetSigmaYFit() * star->GetSigmaYFit();
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 MStarCam
389 // container (with the name "MSourceCam")
390 Int_t fNumStarsSource = 0;
391
392 if (fSourceCam)
393 fNumStarsSource = fSourceCam->GetNumStars();
394
395 //*fLog << "MTelAxisFromStars::Process; fNumStarsSource = "
396 // << fNumStarsSource << endl;
397
398 if (fNumStarsSource > 0)
399 {
400 TArrayD asxy[2];
401 asxy[0].Set(fNumStarsSource);
402 asxy[1].Set(fNumStarsSource);
403
404 TArrayD bsxy[2];
405 bsxy[0].Set(fNumStarsSource);
406 bsxy[1].Set(fNumStarsSource);
407
408 TArrayD esxy[2][2];
409 esxy[0][0].Set(fNumStarsSource);
410 esxy[0][1].Set(fNumStarsSource);
411 esxy[1][0].Set(fNumStarsSource);
412 esxy[1][1].Set(fNumStarsSource);
413
414 MStarPos *starSource = 0;
415 TIter nextSource(fSourceCam->GetList());
416 ix = 0;
417 while ( (starSource = (MStarPos*)nextSource()) )
418 {
419 asxy[0][ix] = starSource->GetXExp();
420 asxy[1][ix] = starSource->GetYExp();
421
422 ix++;
423 }
424
425 TransSkyCam(fLambda, fA, fD, fErrD, asxy, bsxy, esxy);
426
427 // put the estimated positions into the MStarCam container
428 // with name "MSourceCam"
429 TIter setnextSource(fSourceCam->GetList());
430 ix = 0;
431 while ( (starSource = (MStarPos*)setnextSource()) )
432 {
433 Double_t corr = esxy[0][1][ix]/ sqrt( esxy[0][0][ix] * esxy[1][1][ix] );
434 starSource->SetFitValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix],
435 sqrt(esxy[0][0][ix]), sqrt(esxy[1][1][ix]), corr,
436 esxy[0][0][ix], esxy[0][1][ix], esxy[1][1][ix],
437 fChi2, fNdof);
438
439 ix++;
440 }
441
442 }
443
444 //--------------------------------------
445
446 return kTRUE;
447}
448
449//---------------------------------------------------------------------------
450//
451// FindSkyCamTrans
452//
453// This routine determines the transformation
454//
455// ( cos(alfa) -sin(alfa) )
456// b = lambda * A * a + d A = ( )
457// ^ ^ ^ ( sin(alfa) cos(alfa) )
458// | | |
459// scale rotation shift
460// factor matrix
461//
462// from sky coordinates 'a' (projected onto the camera) to camera
463// coordinates 'b', using the positions of known stars in the camera.
464// The latter positions may have been determined by analysing the
465// DC currents in the different pixels.
466//
467// Input : a[2] x and y coordinates of stars projected onto the camera;
468// they were obtained from (RA, dec) of the stars and
469// (ThetaTel, PhiTel) and the time of observation;
470// these are the 'expected positions' of stars in the camera
471// b[2] 'measured positions' of these stars in the camera;
472// they may have been obtained from the DC currents
473// e[2][2] error matrix of b[2]
474// fixedrotationangle value [in degrees] at which rotation angle
475// alfa should be fixed; -1 means don't fix
476// fixedscalefactor value at which scale factor lambda
477// should be fixed; -1 means don't fix
478//
479// Output : lambda, alfadeg, A[2][2], d[2] fit results;
480// parameters describing the transformation
481// from 'expected positions' to the 'measured
482// positions' in the camera
483// errd[2][2] error matrix of d[2]
484// fNumIter number of iterations
485// fNdoF number of degrees of freedom
486// fChi2 chi-square value
487// fChi2Prob chi-square probability
488//
489// The units are assumed to be
490// [degrees] for alfadeg
491// [mm] for a, b, d
492// [1] for lambda
493
494Bool_t MTelAxisFromStars::FindSkyCamTrans(
495 TArrayD a[2], TArrayD b[2], TArrayD e[2][2],
496 Double_t &fixedrotationang, Double_t &fixedscalefac, Double_t &lambda,
497 Double_t &alfadeg, Double_t A[2][2], Double_t d[2], Double_t errd[2][2],
498 Int_t &fNumIter, Int_t &fNdof, Double_t &fChi2, Double_t &fChi2Prob)
499{
500 Int_t fNumStars = a[0].GetSize();
501
502 //*fLog << "MTelAxisFromStars::FindSkyCamTrans; expected and measured positions :"
503 // << endl;
504 for (Int_t ix=0; ix<fNumStars; ix++)
505 {
506 //*fLog << " ix, a[0], a[1], b[0], b[1], errxx, errxy, erryy = "
507 // << ix << " : "
508 // << a[0][ix] << ", " << a[1][ix] << "; "
509 // << b[0][ix] << ", " << b[1][ix] << "; "
510 // << e[0][0][ix] << ", " << e[0][1][ix] << ", "
511 // << e[1][1][ix] << endl;
512 }
513
514
515 //-------------------------------------------
516 // fix some parameters if the number of degrees of freedom is too low
517 // (<= 0.0)
518
519 Double_t fixedscalefactor = fixedscalefac;
520 Double_t fixedrotationangle = fixedrotationang;
521
522 // calculate number of degrees of freedom
523 fNdof = 2 * fNumStars - 4;
524 if (fixedscalefactor != -1.0)
525 fNdof += 1;
526 if (fixedrotationangle != -1.0)
527 fNdof += 1;
528
529 // if there is only 1 star fix both rotation angle and scale factor
530 if (fNumStars == 1)
531 {
532 if (fixedscalefactor == -1.0)
533 {
534 fixedscalefactor = 1.0;
535 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor is fixed at "
536 << fixedscalefactor << endl;
537 }
538 if (fixedrotationangle == -1.0)
539 {
540 fixedrotationangle = 0.0;
541 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; rotation angle is fixed at "
542 << fixedrotationangle << endl;
543 }
544 }
545 // otherwise fix only 1 parameter if possible
546 else if (fNdof < 0)
547 {
548 if (fNdof < -2)
549 {
550 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
551 << fNdof << "; fNumStars = " << fNumStars << endl;
552 return kFALSE;
553 }
554 else if (fNdof == -2)
555 {
556 if (fixedscalefactor == -1.0 && fixedrotationangle == -1.0)
557 {
558 fixedscalefactor = 1.0;
559 fixedrotationangle = 0.0;
560 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor and rotation angle are fixed at "
561 << fixedscalefactor << " and " << fixedrotationangle
562 << " respectively" << endl;
563 }
564 else
565 {
566 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
567 << fNdof << "; fNumStars = " << fNumStars << endl;
568 return kFALSE;
569 }
570 }
571 else if (fNdof == -1)
572 {
573 if (fixedrotationangle == -1.0)
574 {
575 fixedrotationangle = 0.0;
576 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; rotation angle is fixed at "
577 << fixedrotationangle << endl;
578 }
579 else if (fixedscalefactor == -1.0)
580 {
581 fixedscalefactor = 1.0;
582 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor is fixed at "
583 << fixedscalefactor << endl;
584 }
585 else
586 {
587 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
588 << fNdof << "; fNumStars = " << fNumStars<< endl;
589 return kFALSE;
590 }
591 }
592 }
593
594 // recalculate number of degrees of freedom
595 fNdof = 2 * fNumStars - 4;
596 if (fixedscalefactor != -1.0)
597 fNdof += 1;
598 if (fixedrotationangle != -1.0)
599 fNdof += 1;
600
601 if (fNdof < 0)
602 return kFALSE;
603 //-------------------------------------------
604
605
606 // get first approximation of scaling factor
607 if (fixedscalefactor != -1.0)
608 lambda = fixedscalefactor;
609 else
610 lambda = 1.0;
611
612 Double_t lambdaold = lambda;
613 Double_t dlambda = 0.0;
614
615 // get first approximation of rotation angle
616 Double_t alfa = 0.0;
617 if (fixedrotationangle != -1.0)
618 alfa = fixedrotationangle / kRad2Deg;
619
620
621
622 Double_t alfaold = alfa;
623 // maximum allowed change of alfa in 1 iteration step (5 degrees)
624 Double_t dalfamax = 5.0 / kRad2Deg;
625 Double_t dalfa = 0.0;
626
627 Double_t cosal = cos(alfa);
628 Double_t sinal = sin(alfa);
629
630 A[0][0] = cosal;
631 A[0][1] = -sinal;
632 A[1][0] = sinal;
633 A[1][1] = cosal;
634
635
636 Double_t absdold2 = 10000.0;
637 Double_t fChangeofd2 = 10000.0;
638
639
640 TArrayD Aa[2];
641 Aa[0].Set(fNumStars);
642 Aa[1].Set(fNumStars);
643
644
645 Double_t sumEbminlamAa[2];
646
647 TArrayD Ebminlambracd[2];
648 Ebminlambracd[0].Set(fNumStars);
649 Ebminlambracd[1].Set(fNumStars);
650
651 TArrayD EAa[2];
652 EAa[0].Set(fNumStars);
653 EAa[1].Set(fNumStars);
654
655 // invert the error matrices
656 TArrayD c[2][2];
657 c[0][0].Set(fNumStars);
658 c[0][1].Set(fNumStars);
659 c[1][0].Set(fNumStars);
660 c[1][1].Set(fNumStars);
661
662 for (Int_t ix=0; ix<fNumStars; ix++)
663 {
664 Double_t XX = e[0][0][ix];
665 Double_t XY = e[0][1][ix];
666 Double_t YY = e[1][1][ix];
667
668 // get inverse of error matrix
669 Double_t determ = XX*YY - XY*XY;
670 c[0][0][ix] = YY / determ;
671 c[0][1][ix] = -XY / determ;
672 c[1][0][ix] = -XY / determ;
673 c[1][1][ix] = XX / determ;
674 }
675
676
677
678 // calculate sum of inverted error matrices
679 Double_t determsumc;
680 Double_t sumc[2][2];
681 sumc[0][0] = 0.0;
682 sumc[0][1] = 0.0;
683 sumc[1][0] = 0.0;
684 sumc[1][1] = 0.0;
685
686 for (Int_t ix=0; ix<fNumStars; ix++)
687 {
688 sumc[0][0] += c[0][0][ix];
689 sumc[0][1] += c[0][1][ix];
690 sumc[1][0] += c[1][0][ix];
691 sumc[1][1] += c[1][1][ix];
692 }
693 determsumc = sumc[0][0]*sumc[1][1] - sumc[0][1]*sumc[1][0];
694
695 // calculate inverse of sum of inverted error matrices
696 Double_t sumcinv[2][2];
697 sumcinv[0][0] = sumc[1][1] / determsumc;
698 sumcinv[0][1] = -sumc[0][1] / determsumc;
699 sumcinv[1][0] = -sumc[1][0] / determsumc;
700 sumcinv[1][1] = sumc[0][0] / determsumc;
701
702 //*fLog << "sumcinv = " << sumcinv[0][0] << ", " << sumcinv[0][1]
703 // << ", " << sumcinv[1][1] << endl;
704
705
706 // minimize chi2 by iteration ***** start **********************
707
708 // stop iteration when change in |d|*|d| is less than 'told2'
709 // and change in alfa is less than 'toldalfa'
710 // and change in lambda is less than 'toldlambda'
711 // or chi2 is less than 'tolchi2'
712 Double_t told2 = 0.3*0.3; // [mm*mm]; 1/100 of an inner pixel diameter
713 Double_t toldalfa = 0.01 / kRad2Deg; // 0.01 degrees
714 Double_t toldlambda = 0.00006; // uncertainty of 1 mm of distance
715 // between camera and reflector
716 Double_t tolchi2 = 1.e-5;
717
718 Int_t fNumIterMax = 100;
719 fNumIter = 0;
720
721 for (Int_t i=0; i<fNumIterMax; i++)
722 {
723 fNumIter++;
724
725 // get next approximation of d ------------------
726 for (Int_t ix=0; ix<fNumStars; ix++)
727 {
728 Aa[0][ix] = A[0][0] * a[0][ix] + A[0][1]*a[1][ix];
729 Aa[1][ix] = A[1][0] * a[0][ix] + A[1][1]*a[1][ix];
730
731 //*fLog << "ix, Aa = " << ix << " : " << Aa[0][ix] << ", "
732 // << Aa[1][ix] << endl;
733 }
734
735 sumEbminlamAa[0] = 0.0;
736 sumEbminlamAa[1] = 0.0;
737
738 for (Int_t ix=0; ix<fNumStars; ix++)
739 {
740 sumEbminlamAa[0] += c[0][0][ix] * (b[0][ix] - lambda*Aa[0][ix])
741 + c[0][1][ix] * (b[1][ix] - lambda*Aa[1][ix]);
742
743 sumEbminlamAa[1] += c[1][0][ix] * (b[0][ix] - lambda*Aa[0][ix])
744 + c[1][1][ix] * (b[1][ix] - lambda*Aa[1][ix]);
745 }
746
747 //*fLog << "sumEbminlamAa = " << sumEbminlamAa[0] << ", "
748 // << sumEbminlamAa[1] << endl;
749
750 d[0] = sumcinv[0][0] * sumEbminlamAa[0]
751 + sumcinv[0][1] * sumEbminlamAa[1] ;
752
753 d[1] = sumcinv[1][0] * sumEbminlamAa[0]
754 + sumcinv[1][1] * sumEbminlamAa[1] ;
755
756 Double_t absdnew2 = d[0]*d[0] + d[1]*d[1];
757 fChangeofd2 = absdnew2 - absdold2;
758
759 //*fLog << "fNumIter : " << fNumIter
760 // << "; alfa, lambda, d[0], d[1], absdold2, absdnew2 = " << endl;
761 //*fLog << alfa << ", " << lambda << ", " << d[0] << ", " << d[1]
762 // << ", " << absdold2 << ", " << absdnew2 << endl;
763
764
765 if ( fabs(fChangeofd2) < told2 && fabs(dalfa) < toldalfa &&
766 fabs(dlambda) < toldlambda )
767 {
768 //*fLog << "Iteration stopped because of small changes : fChangeofd2, dalfa, dlambda = "
769 // << fChangeofd2 << ", " << dalfa << ", " << dlambda << endl;
770 break;
771 }
772 absdold2 = absdnew2;
773
774 // get next approximation of matrix A ----------------
775 if (fFixedRotationAngle == -1.0)
776 {
777 for (Int_t ix=0; ix<fNumStars; ix++)
778 {
779 Ebminlambracd[0][ix] =
780 c[0][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
781 + c[0][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
782
783 Ebminlambracd[1][ix] =
784 c[1][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
785 + c[1][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
786
787 //*fLog << "ix, Ebminlambracd = " << ix << " : "
788 // << Ebminlambracd[0][ix] << ", "
789 // << Ebminlambracd[1][ix] << endl;
790 }
791
792 // stop iteration if fChi2 is small enough
793 fChi2 = 0.0;
794 for (Int_t ix=0; ix<fNumStars; ix++)
795 {
796 fChi2 += (b[0][ix]-lambda*Aa[0][ix]-d[0] ) * Ebminlambracd[0][ix]
797 + (b[1][ix]-lambda*Aa[1][ix]-d[1] ) * Ebminlambracd[1][ix];
798 }
799 if ( fChi2 < tolchi2 )
800 {
801 //*fLog << "iteration stopped because of small fChi2 : "
802 // << fChi2 << endl;
803 break;
804 }
805
806
807 Double_t dchi2dA[2][2];
808 dchi2dA[0][0] = 0.0;
809 dchi2dA[0][1] = 0.0;
810 dchi2dA[1][0] = 0.0;
811 dchi2dA[1][1] = 0.0;
812
813 for (Int_t ix=0; ix<fNumStars; ix++)
814 {
815 dchi2dA[0][0] += Ebminlambracd[0][ix] * a[0][ix];
816 dchi2dA[0][1] += Ebminlambracd[0][ix] * a[1][ix];
817 dchi2dA[1][0] += Ebminlambracd[1][ix] * a[0][ix];
818 dchi2dA[1][1] += Ebminlambracd[1][ix] * a[1][ix];
819 }
820
821 //*fLog << "dchi2dA = " << dchi2dA[0][0] << ", " << dchi2dA[0][1]
822 // << ", " << dchi2dA[1][0] << ", " << dchi2dA[1][1] << endl;
823
824 // ********* 1st derivative (d chi2) / (d alfa) ************
825 Double_t dchi2dalfa = -2.0*lambda *
826 ( - sinal*(dchi2dA[0][0]+dchi2dA[1][1])
827 + cosal*(dchi2dA[1][0]-dchi2dA[0][1]) );
828
829
830 //Double_t dalfa1st = - fChi2 / dchi2dalfa;
831
832 //*fLog << "fChi2, dchi2dalfa = " << fChi2 << ", "
833 // << dchi2dalfa << endl;
834 //*fLog << "proposed change of alfa using 1st derivative = "
835 // << dalfa1st << endl;
836
837 // ********* 2nd derivative (d2 chi2) / (d alfa2) ******
838 Double_t term1 = 0.0;
839 Double_t term2 = 0.0;
840 Double_t term3 = 0.0;
841 Double_t term4 = 0.0;
842
843 for (Int_t ix=0; ix<fNumStars; ix++)
844 {
845 term1 += a[0][ix]*c[0][0][ix]*a[0][ix] + a[1][ix]*c[1][0][ix]*a[0][ix]
846 + a[0][ix]*c[0][1][ix]*a[1][ix] + a[1][ix]*c[1][1][ix]*a[1][ix];
847
848 term2 += a[0][ix]*c[1][0][ix]*a[0][ix] - a[1][ix]*c[0][0][ix]*a[0][ix]
849 + a[0][ix]*c[1][1][ix]*a[1][ix] - a[1][ix]*c[0][1][ix]*a[1][ix];
850
851 term3 = a[0][ix]*c[0][0][ix]*a[1][ix] + a[1][ix]*c[1][0][ix]*a[1][ix]
852 - a[0][ix]*c[0][1][ix]*a[0][ix] - a[1][ix]*c[1][1][ix]*a[0][ix];
853
854 term4 += a[0][ix]*c[1][0][ix]*a[1][ix] - a[1][ix]*c[0][0][ix]*a[1][ix]
855 - a[0][ix]*c[1][1][ix]*a[0][ix] + a[1][ix]*c[0][1][ix]*a[0][ix];
856 }
857
858 Double_t d2chi2dalfa2 =
859 - 2.0*lambda * ( - cosal*(dchi2dA[0][0]+dchi2dA[1][1])
860 - sinal*(dchi2dA[1][0]-dchi2dA[0][1]) )
861 + 2.0*lambda*lambda * ( sinal*sinal * term1 - sinal*cosal * term2
862 + sinal*cosal * term3 - cosal*cosal * term4);
863
864 // Gauss-Newton step
865 Double_t dalfa2nd = - dchi2dalfa / d2chi2dalfa2;
866
867 //*fLog << "proposed change of alfa using 2st derivative = "
868 // << dalfa2nd << endl;
869
870 //dalfa = dalfa1st;
871 dalfa = dalfa2nd;
872
873 // ******************************************
874
875
876 // restrict change of alfa
877 if ( fabs(dalfa) > dalfamax )
878 {
879 dalfa = TMath::Sign( dalfamax, dalfa );
880 }
881 alfa = alfaold + dalfa;
882
883 if ( alfa < -5.0/kRad2Deg )
884 alfa = -5.0/kRad2Deg;
885 else if ( alfa > 5.0/kRad2Deg )
886 alfa = 5.0/kRad2Deg;
887
888 dalfa = alfa - alfaold;
889
890 alfaold = alfa;
891
892 sinal = sin(alfa);
893 cosal = cos(alfa);
894
895 A[0][0] = cosal;
896 A[0][1] = -sinal;
897 A[1][0] = sinal;
898 A[1][1] = cosal;
899
900 //*fLog << "alfa-alfaold = " << dalfa << endl;
901 //*fLog << "new alfa = " << alfa << endl;
902 }
903
904
905 // get next approximation of lambda ----------------
906 if (fFixedScaleFactor == -1.0)
907 {
908 for (Int_t ix=0; ix<fNumStars; ix++)
909 {
910 Aa[0][ix] = A[0][0]*a[0][ix] + A[0][1]*a[1][ix];
911 Aa[1][ix] = A[1][0]*a[0][ix] + A[1][1]*a[1][ix];
912
913 EAa[0][ix] =
914 c[0][0][ix] * Aa[0][ix] + c[0][1][ix] * Aa[1][ix];
915 EAa[1][ix] =
916 c[1][0][ix] * Aa[0][ix] + c[1][1][ix] * Aa[1][ix];
917
918 //*fLog << "ix, Aa = " << ix << " : " << Aa[0][ix] << ", "
919 // << Aa[1][ix] << endl;
920
921 //*fLog << "ix, EAa = " << ix << " : " << EAa[0][ix] << ", "
922 // << EAa[1][ix] << endl;
923 }
924
925 Double_t num = 0.0;
926 Double_t denom = 0.0;
927 for (Int_t ix=0; ix<fNumStars; ix++)
928 {
929 num += (b[0][ix]-d[0]) * EAa[0][ix]
930 + (b[1][ix]-d[1]) * EAa[1][ix];
931
932 denom += Aa[0][ix] * EAa[0][ix]
933 + Aa[1][ix] * EAa[1][ix];
934
935 //*fLog << "ix : b-d = " << ix << " : " << b[0][ix]-d[0]
936 // << ", " << b[1][ix]-d[1] << endl;
937
938 //*fLog << "ix : Aa = " << ix << " : " << Aa[0][ix]
939 // << ", " << Aa[1][ix] << endl;
940 }
941
942 lambda = num / denom;
943
944 if ( lambda < 0.9 )
945 lambda = 0.9;
946 else if ( lambda > 1.1 )
947 lambda = 1.1;
948
949 dlambda = lambda - lambdaold;
950 lambdaold = lambda;
951
952 //*fLog << "num, denom, lambda, dlambda = " << num
953 // << ", " << denom << ", " << lambda << ", "
954 // << dlambda << endl;
955 }
956
957 }
958 //------- end of iteration *****************************************
959
960 alfadeg = alfa * kRad2Deg;
961
962 // calculate error matrix of d[2]
963 errd[0][0] = sumcinv[0][0];
964 errd[0][1] = sumcinv[0][1];
965 errd[1][0] = sumcinv[1][0];
966 errd[1][1] = sumcinv[1][1];
967
968 // evaluate quality of fit
969
970 // calculate chi2
971 for (Int_t ix=0; ix<fNumStars; ix++)
972 {
973 Ebminlambracd[0][ix] =
974 c[0][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
975 + c[0][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
976
977 Ebminlambracd[1][ix] =
978 c[1][0][ix] * (b[0][ix] - lambda*Aa[0][ix] - d[0] )
979 + c[1][1][ix] * (b[1][ix] - lambda*Aa[1][ix] - d[1] );
980 }
981
982 fChi2 = 0.0;
983 for (Int_t ix=0; ix<fNumStars; ix++)
984 {
985 fChi2 += (b[0][ix] - lambda*Aa[0][ix] - d[0] ) * Ebminlambracd[0][ix]
986 + (b[1][ix] - lambda*Aa[1][ix] - d[1] ) * Ebminlambracd[1][ix];
987 }
988
989 fChi2Prob = TMath::Prob(fChi2, fNdof);
990
991 *fLog << "MTelAxisFromStars::FindSkyCamTrans :" << endl;
992 *fLog << " fNumStars, fChi2, fNdof, fChi2Prob, fNumIter, fChangeofd2, dalfa, dlambda = "
993 << fNumStars << ", " << fChi2 << ", " << fNdof << ", "
994 << fChi2Prob << ", "
995 << fNumIter << ", " << fChangeofd2 << ", " << dalfa << ", "
996 << dlambda << endl;
997 *fLog << " lambda, alfadeg, d[0], d[1] = " << lambda << ", "
998 << alfadeg << ", " << d[0] << ", " << d[1] << endl;
999
1000 return kTRUE;
1001}
1002
1003// --------------------------------------------------------------------------
1004//
1005// Apply transformation (lambda, A, d)
1006// to the expected positions (a[1], a[2])
1007// to obtain the estimated positions (b[1], b[2])
1008//
1009// e[2][2] is the error matrix of b[2]
1010
1011void MTelAxisFromStars::TransSkyCam(
1012 Double_t &lambda, Double_t A[2][2], Double_t d[2], Double_t errd[2][2],
1013 TArrayD a[2], TArrayD b[2], TArrayD e[2][2])
1014{
1015 Int_t numpos = a[0].GetSize();
1016 if (numpos <= 0)
1017 return;
1018
1019 //*fLog << "MTelAxisFromStars::TransSkyCam; expected and estimated positions :"
1020 // << endl;
1021
1022 for (Int_t ix=0; ix<numpos; ix++)
1023 {
1024 //*fLog << "MTelAxisFromStars; ix = " << ix << endl;
1025
1026 b[0][ix] = lambda * (A[0][0]*a[0][ix] + A[0][1]*a[1][ix]) + d[0];
1027 b[1][ix] = lambda * (A[1][0]*a[0][ix] + A[1][1]*a[1][ix]) + d[1];
1028
1029 e[0][0][ix] = errd[0][0];
1030 e[0][1][ix] = errd[0][1];
1031 e[1][0][ix] = errd[1][0];
1032 e[1][1][ix] = errd[1][1];
1033
1034 //*fLog << " ix, a[0], a[1], b[0], b[1], errxx, errxy, erryy = "
1035 // << ix << " : "
1036 // << a[0][ix] << ", " << a[1][ix] << "; "
1037 // << b[0][ix] << ", " << b[1][ix] << "; "
1038 // << e[0][0][ix] << ", " << e[0][1][ix] << ", "
1039 // << e[1][1][ix] << endl;
1040 }
1041}
1042
1043// --------------------------------------------------------------------------
1044//
1045//
1046Int_t MTelAxisFromStars::PostProcess()
1047{
1048
1049 return kTRUE;
1050}
1051
1052
1053// --------------------------------------------------------------------------
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
Note: See TracBrowser for help on using the repository browser.