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

Last change on this file since 4737 was 4726, checked in by wittek, 21 years ago
*** empty log message ***
File size: 31.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#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 // Note : this has to be changed for the wobble mode, where the source
379 // won't be in the center of the camera
380
381 fSrcPos->SetXY(fD[0], fD[1]);
382 fSrcPos->SetReadyToSave();
383
384
385 //--------------------------------------
386 // Apply the transformation to some expected positions (asxy[1], asxy[2])
387 // to obtain estimated positions (bsxy[1], bsxy[2]) in the camera
388 // and their error matrices esxy[2][2]
389
390 // get the expected positions (asxy[1], asxy[2]) from another MStarCam
391 // container (with the name "MSourceCam")
392 Int_t fNumStarsSource = 0;
393
394 if (fSourceCam)
395 fNumStarsSource = fSourceCam->GetNumStars();
396
397 //*fLog << "MTelAxisFromStars::Process; fNumStarsSource = "
398 // << fNumStarsSource << endl;
399
400 if (fNumStarsSource > 0)
401 {
402 TArrayD asxy[2];
403 asxy[0].Set(fNumStarsSource);
404 asxy[1].Set(fNumStarsSource);
405
406 TArrayD bsxy[2];
407 bsxy[0].Set(fNumStarsSource);
408 bsxy[1].Set(fNumStarsSource);
409
410 TArrayD esxy[2][2];
411 esxy[0][0].Set(fNumStarsSource);
412 esxy[0][1].Set(fNumStarsSource);
413 esxy[1][0].Set(fNumStarsSource);
414 esxy[1][1].Set(fNumStarsSource);
415
416 MStarPos *starSource = 0;
417 TIter nextSource(fSourceCam->GetList());
418 ix = 0;
419 while ( (starSource = (MStarPos*)nextSource()) )
420 {
421 asxy[0][ix] = starSource->GetXExp();
422 asxy[1][ix] = starSource->GetYExp();
423
424 ix++;
425 }
426
427 TransSkyCam(fLambda, fA, fD, fErrD, asxy, bsxy, esxy);
428
429 // put the estimated positions into the MStarCam container
430 // with name "MSourceCam"
431 TIter setnextSource(fSourceCam->GetList());
432 ix = 0;
433 while ( (starSource = (MStarPos*)setnextSource()) )
434 {
435 Double_t corr = esxy[0][1][ix]/ sqrt( esxy[0][0][ix] * esxy[1][1][ix] );
436 starSource->SetFitValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix],
437 sqrt(esxy[0][0][ix]), sqrt(esxy[1][1][ix]), corr,
438 esxy[0][0][ix], esxy[0][1][ix], esxy[1][1][ix],
439 fChi2, fNdof);
440
441 ix++;
442 }
443
444 }
445
446 //--------------------------------------
447
448 return kTRUE;
449}
450
451//---------------------------------------------------------------------------
452//
453// FindSkyCamTrans
454//
455// This routine determines the transformation
456//
457// ( cos(alfa) -sin(alfa) )
458// b = lambda * A * a + d A = ( )
459// ^ ^ ^ ( sin(alfa) cos(alfa) )
460// | | |
461// scale rotation shift
462// factor matrix
463//
464// from sky coordinates 'a' (projected onto the camera) to camera
465// coordinates 'b', using the positions of known stars in the camera.
466// The latter positions may have been determined by analysing the
467// DC currents in the different pixels.
468//
469// Input : a[2] x and y coordinates of stars projected onto the camera;
470// they were obtained from (RA, dec) of the stars and
471// (ThetaTel, PhiTel) and the time of observation;
472// these are the 'expected positions' of stars in the camera
473// b[2] 'measured positions' of these stars in the camera;
474// they may have been obtained from the DC currents
475// e[2][2] error matrix of b[2]
476// fixedrotationangle value [in degrees] at which rotation angle
477// alfa should be fixed; -1 means don't fix
478// fixedscalefactor value at which scale factor lambda
479// should be fixed; -1 means don't fix
480//
481// Output : lambda, alfadeg, A[2][2], d[2] fit results;
482// parameters describing the transformation
483// from 'expected positions' to the 'measured
484// positions' in the camera
485// errd[2][2] error matrix of d[2]
486// fNumIter number of iterations
487// fNdoF number of degrees of freedom
488// fChi2 chi-square value
489// fChi2Prob chi-square probability
490//
491// The units are assumed to be
492// [degrees] for alfadeg
493// [mm] for a, b, d
494// [1] for lambda
495
496Bool_t MTelAxisFromStars::FindSkyCamTrans(
497 TArrayD a[2], TArrayD b[2], TArrayD e[2][2],
498 Double_t &fixedrotationang, Double_t &fixedscalefac, Double_t &lambda,
499 Double_t &alfadeg, Double_t A[2][2], Double_t d[2], Double_t errd[2][2],
500 Int_t &fNumIter, Int_t &fNdof, Double_t &fChi2, Double_t &fChi2Prob)
501{
502 Int_t fNumStars = a[0].GetSize();
503
504 //*fLog << "MTelAxisFromStars::FindSkyCamTrans; expected and measured positions :"
505 // << endl;
506 for (Int_t ix=0; ix<fNumStars; ix++)
507 {
508 //*fLog << " ix, a[0], a[1], b[0], b[1], errxx, errxy, erryy = "
509 // << ix << " : "
510 // << a[0][ix] << ", " << a[1][ix] << "; "
511 // << b[0][ix] << ", " << b[1][ix] << "; "
512 // << e[0][0][ix] << ", " << e[0][1][ix] << ", "
513 // << e[1][1][ix] << endl;
514 }
515
516
517 //-------------------------------------------
518 // fix some parameters if the number of degrees of freedom is too low
519 // (<= 0.0)
520
521 Double_t fixedscalefactor = fixedscalefac;
522 Double_t fixedrotationangle = fixedrotationang;
523
524 // calculate number of degrees of freedom
525 fNdof = 2 * fNumStars - 4;
526 if (fixedscalefactor != -1.0)
527 fNdof += 1;
528 if (fixedrotationangle != -1.0)
529 fNdof += 1;
530
531 // if there is only 1 star fix both rotation angle and scale factor
532 if (fNumStars == 1)
533 {
534 if (fixedscalefactor == -1.0)
535 {
536 fixedscalefactor = 1.0;
537 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor is fixed at "
538 << fixedscalefactor << endl;
539 }
540 if (fixedrotationangle == -1.0)
541 {
542 fixedrotationangle = 0.0;
543 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; rotation angle is fixed at "
544 << fixedrotationangle << endl;
545 }
546 }
547 // otherwise fix only 1 parameter if possible
548 else if (fNdof < 0)
549 {
550 if (fNdof < -2)
551 {
552 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
553 << fNdof << "; fNumStars = " << fNumStars << endl;
554 return kFALSE;
555 }
556 else if (fNdof == -2)
557 {
558 if (fixedscalefactor == -1.0 && fixedrotationangle == -1.0)
559 {
560 fixedscalefactor = 1.0;
561 fixedrotationangle = 0.0;
562 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor and rotation angle are fixed at "
563 << fixedscalefactor << " and " << fixedrotationangle
564 << " respectively" << endl;
565 }
566 else
567 {
568 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
569 << fNdof << "; fNumStars = " << fNumStars << endl;
570 return kFALSE;
571 }
572 }
573 else if (fNdof == -1)
574 {
575 if (fixedrotationangle == -1.0)
576 {
577 fixedrotationangle = 0.0;
578 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; rotation angle is fixed at "
579 << fixedrotationangle << endl;
580 }
581 else if (fixedscalefactor == -1.0)
582 {
583 fixedscalefactor = 1.0;
584 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; scale factor is fixed at "
585 << fixedscalefactor << endl;
586 }
587 else
588 {
589 *fLog << warn << "MTelAxisFromStars::FindSkyCamTrans; number of degrees of freedom is too low : "
590 << fNdof << "; fNumStars = " << fNumStars<< endl;
591 return kFALSE;
592 }
593 }
594 }
595
596 // recalculate number of degrees of freedom
597 fNdof = 2 * fNumStars - 4;
598 if (fixedscalefactor != -1.0)
599 fNdof += 1;
600 if (fixedrotationangle != -1.0)
601 fNdof += 1;
602
603 if (fNdof < 0)
604 return kFALSE;
605 //-------------------------------------------
606
607
608 // get first approximation of scaling factor
609 if (fixedscalefactor != -1.0)
610 lambda = fixedscalefactor;
611 else
612 lambda = 1.0;
613
614 Double_t lambdaold = lambda;
615 Double_t dlambda = 0.0;
616
617 // get first approximation of rotation angle
618 Double_t alfa = 0.0;
619 if (fixedrotationangle != -1.0)
620 alfa = fixedrotationangle / kRad2Deg;
621
622
623
624 Double_t alfaold = alfa;
625 // maximum allowed change of alfa in 1 iteration step (5 degrees)
626 Double_t dalfamax = 5.0 / kRad2Deg;
627 Double_t dalfa = 0.0;
628
629 Double_t cosal = cos(alfa);
630 Double_t sinal = sin(alfa);
631
632 A[0][0] = cosal;
633 A[0][1] = -sinal;
634 A[1][0] = sinal;
635 A[1][1] = cosal;
636
637
638 Double_t absdold2 = 10000.0;
639 Double_t fChangeofd2 = 10000.0;
640
641
642 TArrayD Aa[2];
643 Aa[0].Set(fNumStars);
644 Aa[1].Set(fNumStars);
645
646
647 Double_t sumEbminlamAa[2];
648
649 TArrayD Ebminlambracd[2];
650 Ebminlambracd[0].Set(fNumStars);
651 Ebminlambracd[1].Set(fNumStars);
652
653 TArrayD EAa[2];
654 EAa[0].Set(fNumStars);
655 EAa[1].Set(fNumStars);
656
657 // invert the error matrices
658 TArrayD c[2][2];
659 c[0][0].Set(fNumStars);
660 c[0][1].Set(fNumStars);
661 c[1][0].Set(fNumStars);
662 c[1][1].Set(fNumStars);
663
664 for (Int_t ix=0; ix<fNumStars; ix++)
665 {
666 Double_t XX = e[0][0][ix];
667 Double_t XY = e[0][1][ix];
668 Double_t YY = e[1][1][ix];
669
670 // get inverse of error matrix
671 Double_t determ = XX*YY - XY*XY;
672 c[0][0][ix] = YY / determ;
673 c[0][1][ix] = -XY / determ;
674 c[1][0][ix] = -XY / determ;
675 c[1][1][ix] = XX / determ;
676 }
677
678
679
680 // calculate sum of inverted error matrices
681 Double_t determsumc;
682 Double_t sumc[2][2];
683 sumc[0][0] = 0.0;
684 sumc[0][1] = 0.0;
685 sumc[1][0] = 0.0;
686 sumc[1][1] = 0.0;
687
688 for (Int_t ix=0; ix<fNumStars; ix++)
689 {
690 sumc[0][0] += c[0][0][ix];
691 sumc[0][1] += c[0][1][ix];
692 sumc[1][0] += c[1][0][ix];
693 sumc[1][1] += c[1][1][ix];
694 }
695 determsumc = sumc[0][0]*sumc[1][1] - sumc[0][1]*sumc[1][0];
696
697 // calculate inverse of sum of inverted error matrices
698 Double_t sumcinv[2][2];
699 sumcinv[0][0] = sumc[1][1] / determsumc;
700 sumcinv[0][1] = -sumc[0][1] / determsumc;
701 sumcinv[1][0] = -sumc[1][0] / determsumc;
702 sumcinv[1][1] = sumc[0][0] / determsumc;
703
704 //*fLog << "sumcinv = " << sumcinv[0][0] << ", " << sumcinv[0][1]
705 // << ", " << sumcinv[1][1] << endl;
706
707
708 // minimize chi2 by iteration ***** start **********************
709
710 // stop iteration when change in |d|*|d| is less than 'told2'
711 // and change in alfa is less than 'toldalfa'
712 // and change in lambda is less than 'toldlambda'
713 // or chi2 is less than 'tolchi2'
714 Double_t told2 = 0.3*0.3; // [mm*mm]; 1/100 of an inner pixel diameter
715 Double_t toldalfa = 0.01 / kRad2Deg; // 0.01 degrees
716 Double_t toldlambda = 0.00006; // uncertainty of 1 mm of distance
717 // between camera and reflector
718 Double_t tolchi2 = 1.e-5;
719
720 Int_t fNumIterMax = 100;
721 fNumIter = 0;
722
723 for (Int_t i=0; i<fNumIterMax; i++)
724 {
725 fNumIter++;
726
727 // get next approximation of d ------------------
728 for (Int_t ix=0; ix<fNumStars; ix++)
729 {
730 Aa[0][ix] = A[0][0] * a[0][ix] + A[0][1]*a[1][ix];
731 Aa[1][ix] = A[1][0] * a[0][ix] + A[1][1]*a[1][ix];
732
733 //*fLog << "ix, Aa = " << ix << " : " << Aa[0][ix] << ", "
734 // << Aa[1][ix] << endl;
735 }
736
737 sumEbminlamAa[0] = 0.0;
738 sumEbminlamAa[1] = 0.0;
739
740 for (Int_t ix=0; ix<fNumStars; ix++)
741 {
742 sumEbminlamAa[0] += c[0][0][ix] * (b[0][ix] - lambda*Aa[0][ix])
743 + c[0][1][ix] * (b[1][ix] - lambda*Aa[1][ix]);
744
745 sumEbminlamAa[1] += c[1][0][ix] * (b[0][ix] - lambda*Aa[0][ix])
746 + c[1][1][ix] * (b[1][ix] - lambda*Aa[1][ix]);
747 }
748
749 //*fLog << "sumEbminlamAa = " << sumEbminlamAa[0] << ", "
750 // << sumEbminlamAa[1] << endl;
751
752 d[0] = sumcinv[0][0] * sumEbminlamAa[0]
753 + sumcinv[0][1] * sumEbminlamAa[1] ;
754
755 d[1] = sumcinv[1][0] * sumEbminlamAa[0]
756 + sumcinv[1][1] * sumEbminlamAa[1] ;
757
758 Double_t absdnew2 = d[0]*d[0] + d[1]*d[1];
759 fChangeofd2 = absdnew2 - absdold2;
760
761 //*fLog << "fNumIter : " << fNumIter
762 // << "; alfa, lambda, d[0], d[1], absdold2, absdnew2 = " << endl;
763 //*fLog << alfa << ", " << lambda << ", " << d[0] << ", " << d[1]
764 // << ", " << absdold2 << ", " << absdnew2 << endl;
765
766
767 if ( fabs(fChangeofd2) < told2 && fabs(dalfa) < toldalfa &&
768 fabs(dlambda) < toldlambda )
769 {
770 //*fLog << "Iteration stopped because of small changes : fChangeofd2, dalfa, dlambda = "
771 // << fChangeofd2 << ", " << dalfa << ", " << dlambda << endl;
772 break;
773 }
774 absdold2 = absdnew2;
775
776 // get next approximation of matrix A ----------------
777 if (fFixedRotationAngle == -1.0)
778 {
779 for (Int_t ix=0; ix<fNumStars; ix++)
780 {
781 Ebminlambracd[0][ix] =
782 c[0][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
783 + c[0][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
784
785 Ebminlambracd[1][ix] =
786 c[1][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
787 + c[1][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
788
789 //*fLog << "ix, Ebminlambracd = " << ix << " : "
790 // << Ebminlambracd[0][ix] << ", "
791 // << Ebminlambracd[1][ix] << endl;
792 }
793
794 // stop iteration if fChi2 is small enough
795 fChi2 = 0.0;
796 for (Int_t ix=0; ix<fNumStars; ix++)
797 {
798 fChi2 += (b[0][ix]-lambda*Aa[0][ix]-d[0] ) * Ebminlambracd[0][ix]
799 + (b[1][ix]-lambda*Aa[1][ix]-d[1] ) * Ebminlambracd[1][ix];
800 }
801 if ( fChi2 < tolchi2 )
802 {
803 //*fLog << "iteration stopped because of small fChi2 : "
804 // << fChi2 << endl;
805 break;
806 }
807
808
809 Double_t dchi2dA[2][2];
810 dchi2dA[0][0] = 0.0;
811 dchi2dA[0][1] = 0.0;
812 dchi2dA[1][0] = 0.0;
813 dchi2dA[1][1] = 0.0;
814
815 for (Int_t ix=0; ix<fNumStars; ix++)
816 {
817 dchi2dA[0][0] += Ebminlambracd[0][ix] * a[0][ix];
818 dchi2dA[0][1] += Ebminlambracd[0][ix] * a[1][ix];
819 dchi2dA[1][0] += Ebminlambracd[1][ix] * a[0][ix];
820 dchi2dA[1][1] += Ebminlambracd[1][ix] * a[1][ix];
821 }
822
823 //*fLog << "dchi2dA = " << dchi2dA[0][0] << ", " << dchi2dA[0][1]
824 // << ", " << dchi2dA[1][0] << ", " << dchi2dA[1][1] << endl;
825
826 // ********* 1st derivative (d chi2) / (d alfa) ************
827 Double_t dchi2dalfa = -2.0*lambda *
828 ( - sinal*(dchi2dA[0][0]+dchi2dA[1][1])
829 + cosal*(dchi2dA[1][0]-dchi2dA[0][1]) );
830
831
832 //Double_t dalfa1st = - fChi2 / dchi2dalfa;
833
834 //*fLog << "fChi2, dchi2dalfa = " << fChi2 << ", "
835 // << dchi2dalfa << endl;
836 //*fLog << "proposed change of alfa using 1st derivative = "
837 // << dalfa1st << endl;
838
839 // ********* 2nd derivative (d2 chi2) / (d alfa2) ******
840 Double_t term1 = 0.0;
841 Double_t term2 = 0.0;
842 Double_t term3 = 0.0;
843 Double_t term4 = 0.0;
844
845 for (Int_t ix=0; ix<fNumStars; ix++)
846 {
847 term1 += a[0][ix]*c[0][0][ix]*a[0][ix] + a[1][ix]*c[1][0][ix]*a[0][ix]
848 + a[0][ix]*c[0][1][ix]*a[1][ix] + a[1][ix]*c[1][1][ix]*a[1][ix];
849
850 term2 += a[0][ix]*c[1][0][ix]*a[0][ix] - a[1][ix]*c[0][0][ix]*a[0][ix]
851 + a[0][ix]*c[1][1][ix]*a[1][ix] - a[1][ix]*c[0][1][ix]*a[1][ix];
852
853 term3 = a[0][ix]*c[0][0][ix]*a[1][ix] + a[1][ix]*c[1][0][ix]*a[1][ix]
854 - a[0][ix]*c[0][1][ix]*a[0][ix] - a[1][ix]*c[1][1][ix]*a[0][ix];
855
856 term4 += a[0][ix]*c[1][0][ix]*a[1][ix] - a[1][ix]*c[0][0][ix]*a[1][ix]
857 - a[0][ix]*c[1][1][ix]*a[0][ix] + a[1][ix]*c[0][1][ix]*a[0][ix];
858 }
859
860 Double_t d2chi2dalfa2 =
861 - 2.0*lambda * ( - cosal*(dchi2dA[0][0]+dchi2dA[1][1])
862 - sinal*(dchi2dA[1][0]-dchi2dA[0][1]) )
863 + 2.0*lambda*lambda * ( sinal*sinal * term1 - sinal*cosal * term2
864 + sinal*cosal * term3 - cosal*cosal * term4);
865
866 // Gauss-Newton step
867 Double_t dalfa2nd = - dchi2dalfa / d2chi2dalfa2;
868
869 //*fLog << "proposed change of alfa using 2st derivative = "
870 // << dalfa2nd << endl;
871
872 //dalfa = dalfa1st;
873 dalfa = dalfa2nd;
874
875 // ******************************************
876
877
878 // restrict change of alfa
879 if ( fabs(dalfa) > dalfamax )
880 {
881 dalfa = TMath::Sign( dalfamax, dalfa );
882 }
883 alfa = alfaold + dalfa;
884
885 if ( alfa < -5.0/kRad2Deg )
886 alfa = -5.0/kRad2Deg;
887 else if ( alfa > 5.0/kRad2Deg )
888 alfa = 5.0/kRad2Deg;
889
890 dalfa = alfa - alfaold;
891
892 alfaold = alfa;
893
894 sinal = sin(alfa);
895 cosal = cos(alfa);
896
897 A[0][0] = cosal;
898 A[0][1] = -sinal;
899 A[1][0] = sinal;
900 A[1][1] = cosal;
901
902 //*fLog << "alfa-alfaold = " << dalfa << endl;
903 //*fLog << "new alfa = " << alfa << endl;
904 }
905
906
907 // get next approximation of lambda ----------------
908 if (fFixedScaleFactor == -1.0)
909 {
910 for (Int_t ix=0; ix<fNumStars; ix++)
911 {
912 Aa[0][ix] = A[0][0]*a[0][ix] + A[0][1]*a[1][ix];
913 Aa[1][ix] = A[1][0]*a[0][ix] + A[1][1]*a[1][ix];
914
915 EAa[0][ix] =
916 c[0][0][ix] * Aa[0][ix] + c[0][1][ix] * Aa[1][ix];
917 EAa[1][ix] =
918 c[1][0][ix] * Aa[0][ix] + c[1][1][ix] * Aa[1][ix];
919
920 //*fLog << "ix, Aa = " << ix << " : " << Aa[0][ix] << ", "
921 // << Aa[1][ix] << endl;
922
923 //*fLog << "ix, EAa = " << ix << " : " << EAa[0][ix] << ", "
924 // << EAa[1][ix] << endl;
925 }
926
927 Double_t num = 0.0;
928 Double_t denom = 0.0;
929 for (Int_t ix=0; ix<fNumStars; ix++)
930 {
931 num += (b[0][ix]-d[0]) * EAa[0][ix]
932 + (b[1][ix]-d[1]) * EAa[1][ix];
933
934 denom += Aa[0][ix] * EAa[0][ix]
935 + Aa[1][ix] * EAa[1][ix];
936
937 //*fLog << "ix : b-d = " << ix << " : " << b[0][ix]-d[0]
938 // << ", " << b[1][ix]-d[1] << endl;
939
940 //*fLog << "ix : Aa = " << ix << " : " << Aa[0][ix]
941 // << ", " << Aa[1][ix] << endl;
942 }
943
944 lambda = num / denom;
945
946 if ( lambda < 0.9 )
947 lambda = 0.9;
948 else if ( lambda > 1.1 )
949 lambda = 1.1;
950
951 dlambda = lambda - lambdaold;
952 lambdaold = lambda;
953
954 //*fLog << "num, denom, lambda, dlambda = " << num
955 // << ", " << denom << ", " << lambda << ", "
956 // << dlambda << endl;
957 }
958
959 }
960 //------- end of iteration *****************************************
961
962 alfadeg = alfa * kRad2Deg;
963
964 // calculate error matrix of d[2]
965 errd[0][0] = sumcinv[0][0];
966 errd[0][1] = sumcinv[0][1];
967 errd[1][0] = sumcinv[1][0];
968 errd[1][1] = sumcinv[1][1];
969
970 // evaluate quality of fit
971
972 // calculate chi2
973 for (Int_t ix=0; ix<fNumStars; ix++)
974 {
975 Ebminlambracd[0][ix] =
976 c[0][0][ix] * ( b[0][ix] - lambda*Aa[0][ix] - d[0] )
977 + c[0][1][ix] * ( b[1][ix] - lambda*Aa[1][ix] - d[1] );
978
979 Ebminlambracd[1][ix] =
980 c[1][0][ix] * (b[0][ix] - lambda*Aa[0][ix] - d[0] )
981 + c[1][1][ix] * (b[1][ix] - lambda*Aa[1][ix] - d[1] );
982 }
983
984 fChi2 = 0.0;
985 for (Int_t ix=0; ix<fNumStars; ix++)
986 {
987 fChi2 += (b[0][ix] - lambda*Aa[0][ix] - d[0] ) * Ebminlambracd[0][ix]
988 + (b[1][ix] - lambda*Aa[1][ix] - d[1] ) * Ebminlambracd[1][ix];
989 }
990
991 fChi2Prob = TMath::Prob(fChi2, fNdof);
992
993 *fLog << "MTelAxisFromStars::FindSkyCamTrans :" << endl;
994 *fLog << " fNumStars, fChi2, fNdof, fChi2Prob, fNumIter, fChangeofd2, dalfa, dlambda = "
995 << fNumStars << ", " << fChi2 << ", " << fNdof << ", "
996 << fChi2Prob << ", "
997 << fNumIter << ", " << fChangeofd2 << ", " << dalfa << ", "
998 << dlambda << endl;
999 *fLog << " lambda, alfadeg, d[0], d[1] = " << lambda << ", "
1000 << alfadeg << ", " << d[0] << ", " << d[1] << endl;
1001
1002 return kTRUE;
1003}
1004
1005// --------------------------------------------------------------------------
1006//
1007// Apply transformation (lambda, A, d)
1008// to the expected positions (a[1], a[2])
1009// to obtain the estimated positions (b[1], b[2])
1010//
1011// e[2][2] is the error matrix of b[2]
1012
1013void MTelAxisFromStars::TransSkyCam(
1014 Double_t &lambda, Double_t A[2][2], Double_t d[2], Double_t errd[2][2],
1015 TArrayD a[2], TArrayD b[2], TArrayD e[2][2])
1016{
1017 Int_t numpos = a[0].GetSize();
1018 if (numpos <= 0)
1019 return;
1020
1021 //*fLog << "MTelAxisFromStars::TransSkyCam; expected and estimated positions :"
1022 // << endl;
1023
1024 for (Int_t ix=0; ix<numpos; ix++)
1025 {
1026 //*fLog << "MTelAxisFromStars; ix = " << ix << endl;
1027
1028 b[0][ix] = lambda * (A[0][0]*a[0][ix] + A[0][1]*a[1][ix]) + d[0];
1029 b[1][ix] = lambda * (A[1][0]*a[0][ix] + A[1][1]*a[1][ix]) + d[1];
1030
1031 e[0][0][ix] = errd[0][0];
1032 e[0][1][ix] = errd[0][1];
1033 e[1][0][ix] = errd[1][0];
1034 e[1][1][ix] = errd[1][1];
1035
1036 //*fLog << " ix, a[0], a[1], b[0], b[1], errxx, errxy, erryy = "
1037 // << ix << " : "
1038 // << a[0][ix] << ", " << a[1][ix] << "; "
1039 // << b[0][ix] << ", " << b[1][ix] << "; "
1040 // << e[0][0][ix] << ", " << e[0][1][ix] << ", "
1041 // << e[1][1][ix] << endl;
1042 }
1043}
1044
1045// --------------------------------------------------------------------------
1046//
1047//
1048Int_t MTelAxisFromStars::PostProcess()
1049{
1050
1051 return kTRUE;
1052}
1053
1054
1055// --------------------------------------------------------------------------
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
Note: See TracBrowser for help on using the repository browser.