source: trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MTelAxisFromStars.cc@ 5138

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