source: trunk/MagicSoft/Mars/mstarcam/MStarCamTrans.cc@ 4751

Last change on this file since 4751 was 4741, checked in by wittek, 21 years ago
*** empty log message ***
File size: 35.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!
18! Author(s): Wolfgang Wittek, 08/2004 <mailto:wittek@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MStarCamTrans
28// ---------------
29//
30// this is a collection of transformations between
31//
32// a) celestial (declination, hourangle) = ( dec, h); units (deg,hour)
33// b) local (zenith angle, azimuth angle) = (theta, phi); units (deg, deg)
34// c) camera (x-coordinate, y-coordinate) = ( X, Y); units ( mm, mm)
35//
36// coordinates. As to the definition of the coordinates, the conventions
37// of TDAS 00-11 are used.
38//
39// The transformations use
40// - cos(Lat) and sin(Lat) from an MObservatory container,
41// where 'Lat' is the geographical latitude of the telescope.
42// - fDistCam from an MGeomCam container,
43// which is the distance of the camera from the reflector center
44//
45// in order to apply one of the transformations the
46// 'MGeomCam' and 'MObservatory' containers have to be present;
47// the corresponding objects 'mgeom' and 'mobserv' have to be used as
48// arguments of the constructor
49//
50// MStarCamTrans(mgeom, mobserv);
51//
52// The tranformations ignore effects like nutation, precession, refraction, ...
53// This is not a real problem as long as relative positions are considered, as
54// in all the transformations, except in CelToLoc and LocToCel.
55//
56// The camera coordinates (X, Y) are assumed to be the coordinates for an ideal
57// imaging, without imaging errors. The x-axis is assumed to be horizonthal,
58// the y-axis is pointing upwards.
59//
60// the azimuthal angle is defined as : 0 degrees = north
61// 90 degrees = east
62// 180 degrees = south
63// 270 degrees = west
64//
65// Examples for the use of the transformations :
66//
67// - if the local coordinates (theta1, phi1) of a point (X1, Y1) in the camera
68// are known one may caculate the local coordinates (theta2, phi2) of any
69// other point (X2, Y2) in the camera :
70//
71// LocCamCamToLoc(theta1, phi1, X1, Y1, X2, Y2, theta2, phi2);
72//
73// if (X1, Y1) = (0, 0) one may use
74// Loc0CamToLoc(theta1, phi1, X2, Y2, theta2, phi2);
75//
76
77// - if the local coordinates (theta1, phi1) of a point (X1, Y1) in the camera
78// are known one may caculate the position (X2, Y2) for another direction
79// (theta2, phi2) :
80//
81// LocCamLocToCam(theta1, phi1, X1, Y1, theta2, phi2, X2, Y2);
82//
83// if (X1, Y1) = (0, 0) one may use
84// Loc0LocToCam(theta1, phi1, X2, Y2, theta2, phi2);
85//
86// - if the celestial coordinates (dec1, h1) of a point (X1, Y1) in the camera
87// are known one may caculate the celestial coordinates (dec2, h2) of any
88// other point (X2, Y2) in the camera :
89//
90// CelCamCamToCel(dec1, h1, X1, Y1, X2, Y2, dec2, h2);
91//
92// if (X1, Y1) = (0, 0) one may use
93// Cel0CamToCel(dec1, h1, X2, Y2, dec2, h2);
94//
95//
96// - if the celestial coordinates (dec1, h1) of a point (X1, Y1) in the camera
97// are known one may caculate the position (X2, Y2) for any other direction :
98//
99// CelCamCelToCam(dec1, h1, X1, Y1, dec2, h2, X2, Y2);
100//
101// if (X1, Y1) = (0, 0) one may use
102// Cel0CelToCam(dec1, h1, dec2, h2, X2, Y2);
103//
104//
105//
106////////////////////////////////////////////////////////////////////////////
107#include "MStarCamTrans.h"
108
109#include <math.h>
110
111#include <TH1.h>
112#include <TH2.h>
113#include <TPad.h>
114#include <TStyle.h>
115#include <TGraph.h>
116#include <TLatex.h>
117#include <TCanvas.h>
118
119//#include "MH.h"
120#include "MLog.h"
121#include "MLogManip.h"
122#include "MGeomCam.h"
123#include "MObservatory.h"
124
125using namespace std;
126
127ClassImp(MStarCamTrans);
128
129
130// --------------------------------------------------------------------------
131//
132// Constructor
133//
134// get distance between camera and reflector from an MGeomCam container
135// get cos(Lat) and sin(Lat) from an MObservatory container
136//
137MStarCamTrans::MStarCamTrans(const MGeomCam &cam, const MObservatory &obs)
138{
139 if (&cam == NULL)
140 {
141 gLog << err << "MStarCamTrans::MStarCamTrans; MGeomCam container is not available"
142 << endl;
143 return;
144 }
145
146 if (&obs == NULL)
147 {
148 gLog << err << "MStarCamTrans::MStarCamTrans; MObservatory container is not available"
149 << endl;
150 return;
151 }
152
153 fDistCam = cam.GetCameraDist() * 1000.0; // [mm]
154 fCLat = obs.GetCosPhi();
155 fSLat = obs.GetSinPhi();
156
157 gLog << "MStarCamTrans::TransCelLocCam; fDistCam, fCLat, fSLat = "
158 << fDistCam << ", " << fCLat << ", " << fSLat << endl;
159
160 if (fDistCam*fCLat*fSLat == 0.0)
161 {
162 gLog << "MStarCamTrans::TransCelLocCam; one of 'fDistCam, fCLat, fSLat' is zero !!!$$$$$$$$$$$$$$$$$ : "
163 << fDistCam << ", " << fCLat << ", " << fSLat << endl;
164 }
165
166
167 fGridBinning = 0.50; // degrees
168 fGridFineBin = 0.01; // degrees
169
170}
171
172// --------------------------------------------------------------------------
173//
174// Loc0CamToLoc
175//
176// Input : (theta0, phi0) direction for the position (0,0) in the camera
177// ( X, Y) a position in the camera
178//
179// Output : ( theta, phi) direction for the position (X,Y) in the camera
180//
181Bool_t MStarCamTrans::Loc0CamToLoc(Double_t theta0deg, Double_t phi0deg,
182 Double_t X, Double_t Y,
183 Double_t &thetadeg, Double_t &phideg)
184{
185 Double_t theta0 = theta0deg / kRad2Deg;
186 Double_t phi0 = phi0deg / kRad2Deg;
187
188 Double_t XC = X / fDistCam;
189 Double_t YC = Y / fDistCam;
190
191 Double_t theta;
192 Double_t phiminphi0;
193
194 //--------------------------------------------
195 Double_t sip = -XC;
196 Double_t cop = sin(theta0) + YC*cos(theta0);
197
198 Double_t sit = sqrt(cop*cop + XC*XC);
199 Double_t cot = cos(theta0) - YC*sin(theta0);
200
201 // there is an ambiguity in the sign of cos(theta)
202 // - if theta0 is close to 0 or Pi
203 // choose that theta which is closer to theta0,
204 // i.e. require the same sign for cos(theta) and cos(theta0)
205 // - otherwise choose that theta which is compatible with a small
206 // difference |phi-phi0|, i.e. cos(phi-phi0) > 0
207
208 if( fabs(theta0deg - 90.0) > 45.0 )
209 theta = TMath::ATan2( sit, TMath::Sign(cot, cos(theta0)) );
210 else
211 theta = TMath::ATan2( sit, TMath::Sign(cot, cop/cot) );
212
213 Double_t sig = TMath::Sign(1.0, cot*tan(theta));
214 phiminphi0 = TMath::ATan2(sig*sip, sig*cop);
215
216 //--------------------------------------------
217 phideg = (phi0 + phiminphi0) * kRad2Deg;
218 thetadeg = theta * kRad2Deg;
219
220 //gLog << "MStarCamTrans::Loc0Cam2Log; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
221 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
222 // << thetadeg << ", " << phideg << endl;
223
224 return kTRUE;
225}
226
227// --------------------------------------------------------------------------
228//
229// Loc0LocToCam
230//
231// Input : (theta0, phi0) direction for the position (0,0) in the camera
232// ( theta, phi) some other direction
233//
234// Output : (X, Y) position in the camera corresponding to (theta, phi)
235//
236Bool_t MStarCamTrans::Loc0LocToCam(Double_t theta0deg, Double_t phi0deg,
237 Double_t thetadeg, Double_t phideg,
238 Double_t &X, Double_t &Y)
239{
240 Double_t theta0 = theta0deg / kRad2Deg;
241 Double_t phi0 = phi0deg / kRad2Deg;
242
243 Double_t theta = thetadeg / kRad2Deg;
244 Double_t phi = phideg / kRad2Deg;
245
246 //--------------------------------------------
247 Double_t YC = (cos(theta0)*tan(theta)*cos(phi-phi0) - sin(theta0))
248 / (cos(theta0) + sin(theta0)*tan(theta));
249 Double_t XC = -sin(phi-phi0) * (cos(theta0) - YC*sin(theta0)) * tan(theta);
250
251 //--------------------------------------------
252 X = XC * fDistCam;
253 Y = YC * fDistCam;
254
255
256 //gLog << "MStarCamTrans::Loc0LocToCam; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
257 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
258 // << thetadeg << ", " << phideg << endl;
259
260 return kTRUE;
261}
262
263// --------------------------------------------------------------------------
264//
265// LocCamToLoc0
266//
267// Input : ( X, Y) a position in the camera
268// (theta, phi) direction for the position (X,Y) in the camera
269//
270// Output : ( theta0, phi0) direction for the position (0,0) in the camera
271//
272Bool_t MStarCamTrans::LocCamToLoc0(Double_t thetadeg, Double_t phideg,
273 Double_t X, Double_t Y,
274 Double_t &theta0deg, Double_t &phi0deg)
275{
276 Double_t theta = thetadeg / kRad2Deg;
277 Double_t phi = phideg / kRad2Deg;
278
279 Double_t XC = X / fDistCam;
280 Double_t YC = Y / fDistCam;
281
282 //--------------------------------------------
283
284 Double_t theta0;
285 Double_t phiminphi0;
286
287 Double_t tant1 = 0.0;
288 Double_t tant2 = 0.0;
289
290 Double_t ip;
291
292 // calculate tan(theta0)
293 // dummy loop to avoid a 'go to'
294 for (Int_t i=0; i<1; i++)
295 {
296 if ( tan(theta) == 0.0 )
297 {
298 if (XC != 0.0)
299 {
300 // this can never occur
301 gLog << "MStarCamTrans::LocCam2Loc0; thetadeg, XC = " << thetadeg
302 << ", " << XC << endl;
303 return kFALSE;
304 }
305 tant1 = -YC;
306
307 theta0 = TMath::Pi() *
308 modf( (atan(tant1)+TMath::Pi()) / TMath::Pi(), &ip );
309 phiminphi0 = 0.0;
310
311 break;
312 }
313
314 Double_t a = 1.0 + XC*XC - YC*YC*tan(theta)*tan(theta);
315 Double_t b = 2.0 * YC * (1.0 + tan(theta)*tan(theta));
316 Double_t c = XC*XC + YC*YC - tan(theta)*tan(theta);
317
318 if (a == 0.0)
319 {
320 if (b == 0.0)
321 {
322 // this can never occur
323 gLog << "MStarCamTrans::LocCam2Loc0; a, b = " << a << ", "
324 << b << endl;
325 return kFALSE;
326 }
327 tant1 = -c/b;
328 }
329 else
330 {
331 Double_t discr = b*b - 4.0*a*c;
332 if (discr < 0.0)
333 {
334 gLog << "MStarCamTrans::LocCam2Loc0; discr = " << discr << endl;
335 return kFALSE;
336 }
337
338 // two solutions for tan(theta0)
339 tant1 = (-b + sqrt(discr)) / (2.0*a);
340 tant2 = (-b - sqrt(discr)) / (2.0*a);
341
342 if (fabs(tant1-tant2) < 1.e-5)
343 tant2 = 0.0;
344 }
345
346 // define the sign of tan(theta0) and
347 // reject the solution with the wrong sign
348
349 if ( fabs(thetadeg - 90.0) > 45.0 )
350 {
351 Double_t sig = TMath::Sign(1.0, cos(theta));
352 if ( tant1 != 0.0 && TMath::Sign(1.0, tant1) != sig )
353 tant1 = 0.0;
354
355 if ( tant2 != 0.0 && TMath::Sign(1.0, tant2) != sig )
356 tant2 = 0.0;
357 }
358 else
359 {
360 // require sign of cos(phi-phi0) to be > 0
361 if ( tant1 != 0.0 &&
362 TMath::Sign( 1.0, (tant1+YC)/ ((1.0-YC*tant1)*tan(theta)) ) != 1.0 )
363 tant1 = 0.0;
364
365 if ( tant2 != 0.0 &&
366 TMath::Sign( 1.0, (tant2+YC)/ ((1.0-YC*tant2)*tan(theta)) ) != 1.0 )
367 tant2 = 0.0;
368 }
369
370 if (tant1 != 0.0 && tant2 != 0.0)
371 {
372 gLog << "MStarCamTrans::LocCam2Loc0; there are 2 solutions for tan(theta0), tant1, tant2 = "
373 << tant1 << ", " << tant2 << endl;
374 }
375
376 if (tant1 != 0.0)
377 theta0 = TMath::Pi() *
378 modf( (atan(tant1)+TMath::Pi()) / TMath::Pi(), &ip );
379 else if (tant2 != 0.0)
380 theta0 = TMath::Pi() *
381 modf( (atan(tant2)+TMath::Pi()) / TMath::Pi(), &ip );
382 else
383 {
384 theta0 = theta;
385 gLog << "MStarCamTrans::LocCam2Loc0; there is no solution for tan(theta0)"
386 << endl;
387 }
388
389 Double_t sip = -XC;
390 Double_t cop = sin(theta0) + YC*cos(theta0);
391 Double_t cot = cos(theta0) - YC*sin(theta0);
392 Double_t sig = TMath::Sign(1.0, cot*tan(theta));
393
394 phiminphi0 = TMath::ATan2(sig*sip, sig*cop);
395 } // end of dummy loop
396
397 //--------------------------------------------
398 phi0deg = (phi - phiminphi0) * kRad2Deg;
399 theta0deg = theta0 * kRad2Deg;
400
401 //gLog << "MStarCamTrans::LocCamToLoc0; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
402 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
403 // << thetadeg << ", " << phideg << endl;
404
405
406 return kTRUE;
407}
408
409// --------------------------------------------------------------------------
410//
411// Cel0CamToCel
412//
413// Input : (dec0, h0) direction for the position (0,0) in the camera
414// ( X, Y) a position in the camera
415//
416// Output : (dec, h) direction for the position (X,Y) in the camera
417//
418Bool_t MStarCamTrans::Cel0CamToCel(Double_t dec0deg, Double_t h0hour,
419 Double_t X, Double_t Y,
420 Double_t &decdeg, Double_t &hhour)
421{
422 //--------------------------------------------
423
424 // transform celestial coordinates ( dec0, h0)
425 // into local coordinates (theta0, phi0)
426 Double_t theta0deg;
427 Double_t phi0deg;
428 CelToLoc(dec0deg, h0hour, theta0deg, phi0deg);
429
430 // get the local coordinates (theta, phi)
431 // for the position (X, Y) in the camera
432 Double_t thetadeg;
433 Double_t phideg;
434 Loc0CamToLoc(theta0deg, phi0deg, X, Y, thetadeg, phideg);
435
436 // transform local coordinates (theta, phi)
437 // into celestial coordinates ( dec, h)
438 LocToCel(thetadeg, phideg, decdeg, hhour);
439
440 //--------------------------------------------
441
442
443 //gLog << "MStarCamTrans::Cel0CamToCel; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
444 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
445 // << thetadeg << ", " << phideg << endl;
446
447 return kTRUE;
448}
449
450// --------------------------------------------------------------------------
451//
452// Cel0CelToCam
453//
454// Input : (dec0, h0) direction for the position (0,0) in the camera
455// (dec, h) some other direction
456//
457// Output : (X, Y) position in the camera corresponding to (dec, h)
458//
459Bool_t MStarCamTrans::Cel0CelToCam(Double_t dec0deg, Double_t h0hour,
460 Double_t decdeg, Double_t hhour,
461 Double_t &X, Double_t &Y)
462{
463 //--------------------------------------------
464
465 // transform celestial coordinates ( dec0, h0)
466 // into local coordinates (theta0, phi0)
467 Double_t theta0deg;
468 Double_t phi0deg;
469 CelToLoc(dec0deg, h0hour, theta0deg, phi0deg);
470
471 // transform celestial coordinates ( dec, h)
472 // into local coordinates (theta, phi)
473 Double_t thetadeg;
474 Double_t phideg;
475 CelToLoc(decdeg, hhour, thetadeg, phideg);
476
477 // get the position (X, Y) in the camera
478 // from the local coordinates (theta, phi)
479 Loc0LocToCam(theta0deg, phi0deg, thetadeg, phideg, X, Y);
480
481 //--------------------------------------------
482
483 //gLog << "MStarCamTrans::Cel0CelToCam; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
484 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
485 // << thetadeg << ", " << phideg << endl;
486
487 return kTRUE;
488}
489
490// --------------------------------------------------------------------------
491//
492// CelCamToCel0
493//
494// Input : ( X, Y) a position in the camera
495// (dec, h) direction for the position (X,Y) in the camera
496//
497// Output : (dec0, h0) direction for the position (0,0) in the camera
498//
499
500Bool_t MStarCamTrans::CelCamToCel0(Double_t decdeg, Double_t hhour,
501 Double_t X, Double_t Y,
502 Double_t &dec0deg, Double_t &h0hour)
503{
504 //--------------------------------------------
505
506 // transform celestial coordinates ( dec, h)
507 // into local coordinates (theta, phi)
508 Double_t thetadeg;
509 Double_t phideg;
510 CelToLoc(decdeg, hhour, thetadeg, phideg);
511
512 // get the local coordinates (theta, phi)
513 // for the position (0, 0) in the camera
514 Double_t theta0deg;
515 Double_t phi0deg;
516 LocCamToLoc0(thetadeg, phideg, X, Y, theta0deg, phi0deg);
517
518 // transform local coordinates (theta0, phi0)
519 // into celestial coordinates ( dec0, h0)
520 LocToCel(theta0deg, phi0deg, dec0deg, h0hour);
521
522 //--------------------------------------------
523
524
525 //gLog << "MStarCamTrans::CelCamToCel0; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
526 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
527 // << thetadeg << ", " << phideg << endl;
528
529 return kTRUE;
530}
531
532// --------------------------------------------------------------------------
533//
534// LocCamCamToLoc
535//
536// Input : (theta1, phi1) direction for the position (X1,Y1) in the camera
537// ( X1, Y1) a position in the camera
538// ( X2, Y2) another position in the camera
539//
540// Output : (theta2, phi2) direction for the position (X2,Y2) in the camera
541//
542Bool_t MStarCamTrans::LocCamCamToLoc(Double_t theta1deg, Double_t phi1deg,
543 Double_t X1, Double_t Y1,
544 Double_t X2, Double_t Y2,
545 Double_t &theta2deg, Double_t &phi2deg)
546{
547 if (X1 == 0.0 && Y1 == 0.0)
548 {
549 Loc0CamToLoc(theta1deg, phi1deg, X2, Y2, theta2deg, phi2deg);
550 return kTRUE;
551 }
552
553 if (X2 == 0.0 && Y2 == 0.0)
554 {
555 LocCamToLoc0(theta1deg, phi1deg, X1, Y1, theta2deg, phi2deg);
556 return kTRUE;
557 }
558
559 Double_t theta0deg;
560 Double_t phi0deg;
561 LocCamToLoc0(theta1deg, phi1deg, X1, Y1, theta0deg, phi0deg);
562
563 Loc0CamToLoc(theta0deg, phi0deg, X2, Y2, theta2deg, phi2deg);
564
565 return kTRUE;
566}
567
568// --------------------------------------------------------------------------
569//
570// LocCamLocToCam
571//
572// Input : (theta1, phi1) direction for the position (X1,Y1) in the camera
573// ( X1, Y1) a position in the camera
574// (theta2, phi2) another direction
575//
576// Output : ( X2, Y2) position corresponding to (theta2, phi2)
577//
578Bool_t MStarCamTrans::LocCamLocToCam(Double_t theta1deg, Double_t phi1deg,
579 Double_t X1, Double_t Y1,
580 Double_t theta2deg, Double_t phi2deg,
581 Double_t &X2, Double_t &Y2)
582
583{
584 if (X1 == 0.0 && Y1 == 0.0)
585 {
586 Loc0LocToCam(theta1deg, phi1deg, theta2deg, phi2deg, X2, Y2);
587
588 return kTRUE;
589 }
590
591
592 Double_t theta0deg;
593 Double_t phi0deg;
594 LocCamToLoc0(theta1deg, phi1deg, X1, Y1, theta0deg, phi0deg);
595
596 Loc0LocToCam(theta0deg, phi0deg, theta2deg, phi2deg, X2, Y2);
597
598 return kTRUE;
599}
600
601// --------------------------------------------------------------------------
602//
603// CelCamCamToCel
604//
605// Input : (dec1, h1) direction for the position (X1,Y1) in the camera
606// ( X1, Y1) a position in the camera
607// ( X2, Y2) another position in the camera
608//
609// Output : (dec2, h2) direction for the position (X2,Y2) in the camera
610//
611Bool_t MStarCamTrans::CelCamCamToCel(Double_t dec1deg, Double_t h1deg,
612 Double_t X1, Double_t Y1,
613 Double_t X2, Double_t Y2,
614 Double_t &dec2deg, Double_t &h2deg)
615{
616 if (X1 == 0.0 && Y1 == 0.0)
617 {
618 Cel0CamToCel(dec1deg, h1deg, X2, Y2, dec2deg, h2deg);
619 return kTRUE;
620 }
621
622 if (X2 == 0.0 && Y2 == 0.0)
623 {
624 CelCamToCel0(dec1deg, h1deg, X1, Y1, dec2deg, h2deg);
625 return kTRUE;
626 }
627
628 Double_t dec0deg;
629 Double_t h0deg;
630 CelCamToCel0(dec1deg, h1deg, X1, Y1, dec0deg, h0deg);
631
632 Cel0CamToCel(dec0deg, h0deg, X2, Y2, dec2deg, h2deg);
633
634 return kTRUE;
635}
636
637// --------------------------------------------------------------------------
638//
639// CelCamCelToCam
640//
641// Input : (dec1, h1) direction for the position (X1,Y1) in the camera
642// ( X1, Y1) a position in the camera
643// (dec2, h2) another direction
644//
645// Output : ( X2, Y2) position corresponding to (dec2, h2)
646//
647Bool_t MStarCamTrans::CelCamCelToCam(Double_t dec1deg, Double_t h1deg,
648 Double_t X1, Double_t Y1,
649 Double_t dec2deg, Double_t h2deg,
650 Double_t &X2, Double_t &Y2)
651
652{
653 if (X1 == 0.0 && Y1 == 0.0)
654 {
655 Cel0CelToCam(dec1deg, h1deg, dec2deg, h2deg, X2, Y2);
656
657 return kTRUE;
658 }
659
660
661 Double_t dec0deg;
662 Double_t h0deg;
663 CelCamToCel0(dec1deg, h1deg, X1, Y1, dec0deg, h0deg);
664
665 Cel0CelToCam(dec0deg, h0deg, dec2deg, h2deg, X2, Y2);
666
667 return kTRUE;
668}
669
670// --------------------------------------------------------------------------
671//
672// CelToLoc
673//
674// Input : (dec, h) celestial coordinates
675//
676// Output : (theta, phi) corresponding local coordinates
677//
678// (see also MAstroCatalog and MAstroSky2Local)
679//
680Bool_t MStarCamTrans::CelToLoc(Double_t decdeg, Double_t hhour,
681 Double_t &thetadeg, Double_t &phideg)
682{
683 Double_t a1 = fCLat;
684 Double_t a3 = -fSLat;
685
686 Double_t dec = decdeg / kRad2Deg;
687 Double_t h = hhour / 24.0 * TMath::TwoPi();
688
689 // transform celestial coordinates ( dec, h)
690 // into local coordinates (theta, phi)
691 Double_t xB = cos(dec) * cos(h);
692 Double_t yB = cos(dec) * sin(h);
693 Double_t zB = - sin(dec);
694
695 Double_t xA = a3*xB - a1*zB;
696 Double_t yA = - yB;
697 Double_t zA = -a1*xB - a3*zB;
698
699 thetadeg = acos(-zA) * kRad2Deg;
700 phideg = TMath::ATan2(yA, xA) * kRad2Deg;
701
702 return kTRUE;
703}
704
705// --------------------------------------------------------------------------
706//
707// LocToCel
708//
709// Input : (theta, phi) local coordinates
710//
711// Output : (dec, h) corresponding celestial coordinates
712//
713// (see also MAstroCatalog and MAstroSky2Local)
714//
715Bool_t MStarCamTrans::LocToCel(Double_t thetadeg, Double_t phideg,
716 Double_t &decdeg, Double_t &hhour)
717{
718 Double_t a1 = fCLat;
719 Double_t a3 = -fSLat;
720
721 Double_t theta = thetadeg / kRad2Deg;
722 Double_t phi = phideg / kRad2Deg;
723
724 //--------------------------------------------
725
726 // transform local coordinates (theta, phi)
727 // into celestial coordinates ( dec, h)
728 Double_t xA = sin(theta) * cos(phi);
729 Double_t yA = sin(theta) * sin(phi);
730 Double_t zA = - cos(theta);
731
732 Double_t xB = a3*xA - a1*zA;
733 Double_t yB = - yA;
734 Double_t zB = -a1*xA - a3*zA;
735
736 Double_t dec = asin(-zB);
737 Double_t h = TMath::ATan2(yB, xB);
738
739 //--------------------------------------------
740 decdeg = dec * kRad2Deg;
741 hhour = h * 24.0 / TMath::TwoPi();
742
743
744 return kTRUE;
745}
746
747// --------------------------------------------------------------------------
748//
749// PlotGridAtDec0H0
750//
751// set the celestial coordinates of the camera center
752// and plot the lines of constant (Theta, Phi)
753// and the lines of constant (Dec, H )
754//
755// (see also MAstroCatalog::Draw and MAstroCamera::Draw)
756//
757// Warning: Leaks Memory!
758//
759Bool_t MStarCamTrans::PlotGridAtDec0H0(TString name,
760 Double_t dec0deg, Double_t h0hour)
761{
762 Double_t theta0deg;
763 Double_t phi0deg;
764 CelToLoc(dec0deg, h0hour, theta0deg, phi0deg);
765
766 PlotGrid(name, theta0deg, phi0deg, dec0deg, h0hour);
767
768 return kTRUE;
769}
770
771// --------------------------------------------------------------------------
772//
773// PlotGridAtTheta0Phi0
774//
775// set the local coordinates of the camera center
776// and plot the lines of constant (Theta, Phi)
777// and the lines of constant (Dec, H )
778//
779// (see also MAstroCatalog::Draw and MAstroCamera::Draw)
780//
781// Warning: Leaks Memory!
782//
783Bool_t MStarCamTrans::PlotGridAtTheta0Phi0(TString name,
784 Double_t theta0deg, Double_t phi0deg)
785{
786 Double_t dec0deg;
787 Double_t h0hour;
788 LocToCel(theta0deg, phi0deg, dec0deg, h0hour);
789
790 PlotGrid(name, theta0deg, phi0deg, dec0deg, h0hour);
791
792 return kTRUE;
793}
794
795// --------------------------------------------------------------------------
796//
797// SetGridParameters
798//
799// set the binning for the grid (fGridBinning)
800// set the binning along the lines (fGridFineBin)
801//
802Bool_t MStarCamTrans::SetGridParameters(
803 Double_t gridbinning, Double_t gridfinebin)
804{
805 fGridBinning = gridbinning;
806 fGridFineBin = gridfinebin;
807
808 return kTRUE;
809}
810
811// --------------------------------------------------------------------------
812//
813// PlotGrid
814//
815// - plot the lines of constant (Theta, Phi)
816// with the camera center at (Theta0, Phi0)
817// - and the lines of constant (Dec, H )
818// with the camera center at (Dec0, H0 )
819//
820// (see also MAstroCatalog::Draw and MAstroCamera::Draw)
821//
822// Warning: Leaks Memory!
823//
824Bool_t MStarCamTrans::PlotGrid(TString name,
825 Double_t theta0deg, Double_t phi0deg,
826 Double_t dec0deg, Double_t h0hour)
827{
828 Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ;
829 //gLog << "mmtodeg = " << mmtodeg << endl;
830
831 // aberr is the ratio r_optaberr/r_ideal between
832 // the distance from the camera center with optical aberration and
833 // the distance from the camera center with an ideal imaging
834 // the value 1.07 is valid if the expected position (with aberration)
835 // in the camera is calculated as the average of the positions obtained
836 // from the reflections at the individual mirrors
837 Double_t aberr = 1.07;
838
839 //--------------------------------------------------------------------
840
841 TCanvas *c1 = new TCanvas(name, name, 0, 0, 300, 600);
842
843 //gROOT->Reset();
844 gStyle->SetCanvasColor(0);
845 gStyle->SetCanvasBorderMode(0);
846 gStyle->SetPadBorderMode(0);
847 gStyle->SetFrameBorderMode(0);
848 gStyle->SetOptStat(0000000);
849 gStyle->SetPalette(1);
850
851 c1->Divide(1,2);
852 c1->SetBorderMode(0);
853
854 Double_t xlo = -534.0 * mmtodeg;
855 Double_t xup = 534.0 * mmtodeg;
856
857 Double_t ylo = -534.0 * mmtodeg;
858 Double_t yup = 534.0 * mmtodeg;
859
860 TString nam1 = name;
861 nam1 += "_:_Theta-Phi";
862 TString tit1 = name;
863 tit1 += "_:_Theta-Phi-lines";
864 TH2D *plot1 = new TH2D(nam1, tit1, 1, xlo, xup, 1, ylo, yup);
865 plot1->GetXaxis()->SetTitle("x [deg]");
866 plot1->GetYaxis()->SetTitle("y [deg]");
867
868
869 TString nam2 = name;
870 nam2 += "_:_Dec-Hour";
871 TString tit2 = name;
872 tit2 += "_:_Dec-Hour-lines";
873 TH2D *plot2 = new TH2D(nam2, tit2, 1, xlo, xup, 1, ylo, yup);
874 plot2->GetXaxis()->SetTitle("x [deg]");
875 plot2->GetYaxis()->SetTitle("y [deg]");
876
877
878 c1->cd(1);
879 plot1->Draw();
880 //delete plot1;
881
882 c1->cd(2);
883 plot2->Draw();
884 //delete plot2;
885
886 TGraph *graph1;
887 TLatex *pix;
888
889 Char_t tit[100];
890 Double_t xtxt;
891 Double_t ytxt;
892
893 Double_t xx;
894 Double_t yy;
895
896 Double_t gridbinning = fGridBinning;
897 Double_t gridfinebin = fGridFineBin;
898 Int_t numgridlines = (Int_t)(4.0/gridbinning);
899
900
901 //--------------------------------------------------------------------
902 // Theta-Phi lines
903
904 // direction for the center of the camera
905 Double_t theta0 = theta0deg;
906 Double_t phi0 = phi0deg;
907
908
909 //----- lines for fixed theta ------------------------------------
910
911 const Int_t Ntheta = numgridlines;
912 Double_t theta[Ntheta];
913 Double_t dtheta = gridbinning;
914
915 Int_t nphi = (Int_t)(4.0/gridfinebin);
916 const Int_t Nphi = nphi+1;
917 Double_t phi[Nphi];
918 Double_t dphi = gridfinebin/sin(theta0/180.0*3.1415926);
919 if ( dphi > 360.0/((Double_t)(Nphi-1)) )
920 dphi = 360.0/((Double_t)(Nphi-1));
921
922 for (Int_t j=0; j<Ntheta; j++)
923 {
924 theta[j] = theta0 + ((Double_t)j)*dtheta
925 - ((Double_t)(Ntheta/2)-1.0)*dtheta;
926 }
927
928 for (Int_t k=0; k<Nphi; k++)
929 {
930 phi[k] = phi0 + ((Double_t)k)*dphi - ((Double_t)(Nphi/2)-1.0)*dphi;
931 }
932
933
934 Double_t x[Nphi];
935 Double_t y[Nphi];
936
937
938
939 for (Int_t j=0; j<Ntheta; j++)
940 {
941 if (theta[j] < 0.0) continue;
942
943 for (Int_t k=0; k<Nphi; k++)
944 {
945 Loc0LocToCam(theta0, phi0, theta[j], phi[k],
946 xx, yy);
947 x[k] = xx * mmtodeg * aberr;
948 y[k] = yy * mmtodeg * aberr;
949
950 //gLog << "theta0, phi0 = " << theta0 << ", " << phi0
951 // << " : theta, phi, x, y = " << theta[j] << ", "
952 // << phi[k] << "; " << x[k] << ", " << y[k] << endl;
953 }
954
955 c1->cd(1);
956 graph1 = new TGraph(Nphi,x,y);
957 graph1->SetLineColor(j+1);
958 graph1->SetLineStyle(1);
959 graph1->SetMarkerColor(j+1);
960 graph1->SetMarkerSize(.2);
961 graph1->SetMarkerStyle(20);
962 graph1->Draw("PL");
963 //delete graph1;
964
965 sprintf(tit,"Theta = %6.2f", theta[j]);
966 xtxt = xlo + (xup-xlo)*0.1;
967 ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)j) *(yup-ylo)/20.0;
968 pix = new TLatex(xtxt, ytxt, tit);
969 pix->SetTextColor(j+1);
970 pix->SetTextSize(.03);
971 pix->Draw("");
972 //delete pix;
973
974 }
975
976 //----- lines for fixed phi ------------------------------------
977
978 Int_t ntheta1 = (Int_t)(4.0/gridfinebin);
979 const Int_t Ntheta1 = ntheta1;
980 Double_t theta1[Ntheta1];
981 Double_t dtheta1 = gridfinebin;
982
983 const Int_t Nphi1 = numgridlines;
984 Double_t phi1[Nphi1];
985 Double_t dphi1 = gridbinning/sin(theta0/180.0*3.1415926);
986 if ( dphi1 > 360.0/((Double_t)Nphi1) )
987 dphi1 = 360.0/((Double_t)Nphi1);
988
989 for (Int_t j=0; j<Ntheta1; j++)
990 {
991 theta1[j] = theta0 + ((Double_t)j)*dtheta1
992 - ((Double_t)(Ntheta1/2)-1.0)*dtheta1;
993 }
994
995 for (Int_t k=0; k<Nphi1; k++)
996 {
997 phi1[k] = phi0 + ((Double_t)k)*dphi1 - ((Double_t)(Nphi1/2)-1.0)*dphi1;
998 }
999
1000
1001 Double_t x1[Ntheta1];
1002 Double_t y1[Ntheta1];
1003
1004 for (Int_t k=0; k<Nphi1; k++)
1005 {
1006 Int_t count = 0;
1007 for (Int_t j=0; j<Ntheta1; j++)
1008 {
1009 if (theta1[j] < 0.0) continue;
1010
1011 Loc0LocToCam(theta0, phi0, theta1[j], phi1[k],
1012 xx, yy);
1013 x1[count] = xx * mmtodeg * aberr;
1014 y1[count] = yy * mmtodeg * aberr;
1015
1016 //gLog << "theta0, phi0 = " << theta0 << ", " << phi0
1017 // << " : theta1, phi1, x1, y1 = " << theta1[j] << ", "
1018 // << phi1[k] << "; " << x1[count] << ", " << y1[count] << endl;
1019 count++;
1020 }
1021
1022 c1->cd(1);
1023 graph1 = new TGraph(count,x1,y1);
1024 graph1->SetLineColor(k+1);
1025 graph1->SetLineStyle(2);
1026 graph1->SetMarkerColor(k+1);
1027 graph1->SetMarkerSize(.2);
1028 graph1->SetMarkerStyle(20);
1029 graph1->Draw("PL");
1030 //delete graph1;
1031
1032 sprintf(tit,"Phi = %6.2f", phi1[k]);
1033 Double_t xtxt = xlo + (xup-xlo)*0.75;
1034 Double_t ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)k) *(yup-ylo)/20.0;
1035 pix = new TLatex(xtxt, ytxt, tit);
1036 pix->SetTextColor(k+1);
1037 pix->SetTextSize(.03);
1038 pix->Draw("");
1039 //delete pix;
1040
1041 }
1042
1043 c1->cd(1);
1044 sprintf(tit,"Theta0 = %6.2f Phi0 = %6.2f", theta0, phi0);
1045 xtxt = xlo + (xup-xlo)*0.05;
1046 ytxt = ylo + (yup-ylo)*0.92;
1047 pix = new TLatex(xtxt, ytxt, tit);
1048 pix->SetTextColor(1);
1049 pix->SetTextSize(.06);
1050 pix->Draw("");
1051 //delete pix;
1052
1053 sprintf(tit," [deg] [deg]");
1054 xtxt = xlo + (xup-xlo)*0.05;
1055 ytxt = ylo + (yup-ylo)*0.86;
1056 pix = new TLatex(xtxt, ytxt, tit);
1057 pix->SetTextColor(1);
1058 pix->SetTextSize(.06);
1059 pix->Draw("");
1060 //delete pix;
1061
1062
1063 //--------------------------------------------------------------------
1064 // Dec-Hour lines
1065
1066 // direction for the center of the camera
1067 Double_t dec0 = dec0deg;
1068 Double_t h0 = h0hour;
1069 //trans.LocToCel(theta0, phi0, dec0, h0);
1070
1071
1072 //----- lines for fixed dec ------------------------------------
1073
1074 const Int_t Ndec = numgridlines;
1075 Double_t dec[Ndec];
1076 Double_t ddec = gridbinning;
1077
1078 Int_t nh = (Int_t)(4.0/gridfinebin);
1079 const Int_t Nh = nh+1;
1080 Double_t h[Nh];
1081 Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926);
1082 if ( dh > 360.0/((Double_t)(Nh-1)) )
1083 dh = 360.0/((Double_t)(Nh-1));
1084
1085
1086 for (Int_t j=0; j<Ndec; j++)
1087 {
1088 dec[j] = dec0 + ((Double_t)j)*ddec
1089 - ((Double_t)(Ndec/2)-1.0)*ddec;
1090 }
1091
1092 for (Int_t k=0; k<Nh; k++)
1093 {
1094 h[k] = h0 + ((Double_t)k)*dh - ((Double_t)(Nh/2)-1.0)*dh;
1095 }
1096
1097
1098 Double_t xh[Nh];
1099 Double_t yh[Nh];
1100
1101
1102
1103 for (Int_t j=0; j<Ndec; j++)
1104 {
1105 if (fabs(dec[j]) > 90.0) continue;
1106
1107 for (Int_t k=0; k<Nh; k++)
1108 {
1109 Double_t hh0 = h0 *24.0/360.0;
1110 Double_t hx = h[k]*24.0/360.0;
1111 Cel0CelToCam(dec0, hh0, dec[j], hx,
1112 xx, yy);
1113 xh[k] = xx * mmtodeg * aberr;
1114 yh[k] = yy * mmtodeg * aberr;
1115
1116 //gLog << "dec0, h0 = " << dec0 << ", " << h0
1117 // << " : dec, h, xh, yh = " << dec[j] << ", "
1118 // << h[k] << "; " << xh[k] << ", " << yh[k] << endl;
1119 }
1120
1121 c1->cd(2);
1122 graph1 = new TGraph(Nh,xh,yh);
1123 graph1->SetLineColor(j+1);
1124 graph1->SetLineStyle(1);
1125 graph1->SetMarkerColor(j+1);
1126 graph1->SetMarkerSize(.2);
1127 graph1->SetMarkerStyle(20);
1128 graph1->Draw("PL");
1129 //delete graph1;
1130
1131 sprintf(tit,"Dec = %6.2f", dec[j]);
1132 xtxt = xlo + (xup-xlo)*0.1;
1133 ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)j) *(yup-ylo)/20.0;
1134 pix = new TLatex(xtxt, ytxt, tit);
1135 pix->SetTextColor(j+1);
1136 pix->SetTextSize(.03);
1137 pix->Draw("");
1138 //delete pix;
1139
1140 }
1141
1142 //----- lines for fixed hour angle ------------------------------------
1143
1144 Int_t ndec1 = (Int_t)(4.0/gridfinebin);
1145 const Int_t Ndec1 = ndec1;
1146 Double_t dec1[Ndec1];
1147 Double_t ddec1 = gridfinebin;
1148
1149 const Int_t Nh1 = numgridlines;
1150 Double_t h1[Nh1];
1151 Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926);
1152 if ( dh1 > 360.0/((Double_t)Nh1) )
1153 dh1 = 360.0/((Double_t)Nh1);
1154
1155
1156 for (Int_t j=0; j<Ndec1; j++)
1157 {
1158 dec1[j] = dec0 + ((Double_t)j)*ddec1
1159 - ((Double_t)(Ndec1/2)-1.0)*ddec1;
1160 }
1161
1162 for (Int_t k=0; k<Nh1; k++)
1163 {
1164 h1[k] = h0 + ((Double_t)k)*dh1 - ((Double_t)(Nh1/2)-1.0)*dh1;
1165 }
1166
1167
1168 Double_t xd[Ndec1];
1169 Double_t yd[Ndec1];
1170
1171 for (Int_t k=0; k<Nh1; k++)
1172 {
1173 Int_t count = 0;
1174 for (Int_t j=0; j<Ndec1; j++)
1175 {
1176 if (fabs(dec1[j]) > 90.0) continue;
1177
1178 Double_t hh0 = h0 *24.0/360.0;
1179 Double_t hhx = h1[k]*24.0/360.0;
1180 Cel0CelToCam(dec0, hh0, dec1[j], hhx,
1181 xx, yy);
1182 xd[count] = xx * mmtodeg * aberr;
1183 yd[count] = yy * mmtodeg * aberr;
1184
1185 //gLog << "dec0, h0 = " << dec0 << ", " << h0
1186 // << " : dec1, h1, xd, yd = " << dec1[j] << ", "
1187 // << h1[k] << "; " << xd[count] << ", " << yd[count] << endl;
1188
1189 count++;
1190 }
1191
1192 c1->cd(2);
1193 graph1 = new TGraph(count,xd,yd);
1194 graph1->SetLineColor(k+1);
1195 graph1->SetLineStyle(2);
1196 graph1->SetMarkerColor(k+1);
1197 graph1->SetMarkerSize(.2);
1198 graph1->SetMarkerStyle(20);
1199 graph1->Draw("PL");
1200 //delete graph1;
1201
1202 sprintf(tit,"H = %6.2f", h1[k]);
1203 Double_t xtxt = xlo + (xup-xlo)*0.75;
1204 Double_t ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)k) *(yup-ylo)/20.0;
1205 pix = new TLatex(xtxt, ytxt, tit);
1206 pix->SetTextColor(k+1);
1207 pix->SetTextSize(.03);
1208 pix->Draw("");
1209 //delete pix;
1210
1211 }
1212
1213 c1->cd(2);
1214 sprintf(tit,"Dec0 = %6.2f H0 = %6.2f", dec0, h0);
1215 xtxt = xlo + (xup-xlo)*0.05;
1216 ytxt = ylo + (yup-ylo)*0.92;
1217 pix = new TLatex(xtxt, ytxt, tit);
1218 pix->SetTextColor(1);
1219 pix->SetTextSize(.06);
1220 pix->Draw("");
1221 //delete pix;
1222
1223 sprintf(tit," [deg] [deg]");
1224 xtxt = xlo + (xup-xlo)*0.05;
1225 ytxt = ylo + (yup-ylo)*0.86;
1226 pix = new TLatex(xtxt, ytxt, tit);
1227 pix->SetTextColor(1);
1228 pix->SetTextSize(.06);
1229 pix->Draw("");
1230 //delete pix;
1231
1232 return kTRUE;
1233}
1234
Note: See TracBrowser for help on using the repository browser.