source: trunk/MagicSoft/Mars/mastro/MTransCelLocCam.cc@ 4674

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