Index: trunk/MagicSoft/Mars/mastro/MTransCelLocCam.cc =================================================================== --- trunk/MagicSoft/Mars/mastro/MTransCelLocCam.cc (revision 4674) +++ trunk/MagicSoft/Mars/mastro/MTransCelLocCam.cc (revision 4674) @@ -0,0 +1,1236 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek, 08/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MTransCelLocCam +// --------------- +// +// this is a collection of transformations between +// +// a) celestial (declination, hourangle) = ( dec, h); units (deg,hour) +// b) local (zenith angle, azimuth angle) = (theta, phi); units (deg, deg) +// c) camera (x-coordinate, y-coordinate) = ( X, Y); units ( mm, mm) +// +// coordinates. As to the definition of the coordinates, the conventions +// of TDAS 00-11 are used. +// +// The transformations use +// - cos(Lat) and sin(Lat) from an MObservatory container, +// where 'Lat' is the geographical latitude of the telescope. +// - fDistCam from an MGeomCam container, +// which is the distance of the camera from the reflector center +// +// in order to apply one of the transformations the +// 'MGeomCam' and 'MObservatory' containers have to be present; +// the corresponding objects 'mgeom' and 'mobserv' have to be used as +// arguments of the constructor +// +// MTransCelLocCam(mgeom, mobserv); +// +// The tranformations ignore effects like nutation, precession, refraction, ... +// This is not a real problem as long as relative positions are considered, as +// in all the transformations, except in CelToLoc and LocToCel. +// +// The camera coordinates (X, Y) are assumed to be the coordinates for an ideal +// imaging, without imaging errors. The x-axis is assumed to be horizonthal, +// the y-axis is pointing upwards. +// +// the azimuthal angle is defined as : 0 degrees = north +// 90 degrees = east +// 180 degrees = south +// 270 degrees = west +// +// Examples for the use of the transformations : +// +// - if the local coordinates (theta1, phi1) of a point (X1, Y1) in the camera +// are known one may caculate the local coordinates (theta2, phi2) of any +// other point (X2, Y2) in the camera : +// +// LocCamCamToLoc(theta1, phi1, X1, Y1, X2, Y2, theta2, phi2); +// +// if (X1, Y1) = (0, 0) one may use +// Loc0CamToLoc(theta1, phi1, X2, Y2, theta2, phi2); +// + +// - if the local coordinates (theta1, phi1) of a point (X1, Y1) in the camera +// are known one may caculate the position (X2, Y2) for another direction +// (theta2, phi2) : +// +// LocCamLocToCam(theta1, phi1, X1, Y1, theta2, phi2, X2, Y2); +// +// if (X1, Y1) = (0, 0) one may use +// Loc0LocToCam(theta1, phi1, X2, Y2, theta2, phi2); +// +// - if the celestial coordinates (dec1, h1) of a point (X1, Y1) in the camera +// are known one may caculate the celestial coordinates (dec2, h2) of any +// other point (X2, Y2) in the camera : +// +// CelCamCamToCel(dec1, h1, X1, Y1, X2, Y2, dec2, h2); +// +// if (X1, Y1) = (0, 0) one may use +// Cel0CamToCel(dec1, h1, X2, Y2, dec2, h2); +// +// +// - if the celestial coordinates (dec1, h1) of a point (X1, Y1) in the camera +// are known one may caculate the position (X2, Y2) for any other direction : +// +// CelCamCelToCam(dec1, h1, X1, Y1, dec2, h2, X2, Y2); +// +// if (X1, Y1) = (0, 0) one may use +// Cel0CelToCam(dec1, h1, dec2, h2, X2, Y2); +// +// +// +//////////////////////////////////////////////////////////////////////////// +#include "MTransCelLocCam.h" + +#include + +#include +#include +#include +#include +#include +#include +#include + +//#include "MH.h" +#include "MLog.h" +#include "MLogManip.h" +#include "MGeomCam.h" +#include "MObservatory.h" + +using namespace std; + +ClassImp(MTransCelLocCam); + + +// -------------------------------------------------------------------------- +// +// Constructor +// +// get distance between camera and reflector from an MGeomCam container +// get cos(Lat) and sin(Lat) from an MObservatory container +// +MTransCelLocCam::MTransCelLocCam(const MGeomCam &cam, const MObservatory &obs) +{ + if (&cam == NULL) + { + gLog << err << "MTransCelLocCam::MTransCelLocCam; MGeomCam container is not available" + << endl; + return; + } + + if (&obs == NULL) + { + gLog << err << "MTransCelLocCam::MTransCelLocCam; MObservatory container is not available" + << endl; + return; + } + + fDistCam = cam.GetCameraDist() * 1000.0; // [mm] + fCLat = obs.GetCosPhi(); + fSLat = obs.GetSinPhi(); + + gLog << "MTransCelLocCam::TransCelLocCam; fDistCam, fCLat, fSLat = " + << fDistCam << ", " << fCLat << ", " << fSLat << endl; + + if (fDistCam*fCLat*fSLat == 0.0) + { + gLog << "MTransCelLocCam::TransCelLocCam; one of 'fDistCam, fCLat, fSLat' is zero !!!$$$$$$$$$$$$$$$$$ : " + << fDistCam << ", " << fCLat << ", " << fSLat << endl; + } + + + fGridBinning = 0.50; // degrees + fGridFineBin = 0.01; // degrees + +} + +// -------------------------------------------------------------------------- +// +// Loc0CamToLoc +// +// Input : (theta0, phi0) direction for the position (0,0) in the camera +// ( X, Y) a position in the camera +// +// Output : ( theta, phi) direction for the position (X,Y) in the camera +// + +Bool_t MTransCelLocCam::Loc0CamToLoc(Double_t theta0deg, Double_t phi0deg, + Double_t X, Double_t Y, + Double_t &thetadeg, Double_t &phideg) +{ + Double_t theta0 = theta0deg / kRad2Deg; + Double_t phi0 = phi0deg / kRad2Deg; + + Double_t XC = X / fDistCam; + Double_t YC = Y / fDistCam; + + Double_t theta; + Double_t phiminphi0; + + //-------------------------------------------- + Double_t sip = -XC; + Double_t cop = sin(theta0) + YC*cos(theta0); + + Double_t sit = sqrt(cop*cop + XC*XC); + Double_t cot = cos(theta0) - YC*sin(theta0); + + // there is an ambiguity in the sign of cos(theta) + // - if theta0 is close to 0 or Pi + // choose that theta which is closer to theta0, + // i.e. require the same sign for cos(theta) and cos(theta0) + // - otherwise choose that theta which is compatible with a small + // difference |phi-phi0|, i.e. cos(phi-phi0) > 0 + + if( fabs(theta0deg - 90.0) > 45.0 ) + theta = TMath::ATan2( sit, TMath::Sign(cot, cos(theta0)) ); + else + theta = TMath::ATan2( sit, TMath::Sign(cot, cop/cot) ); + + Double_t sig = TMath::Sign(1.0, cot*tan(theta)); + phiminphi0 = TMath::ATan2(sig*sip, sig*cop); + + //-------------------------------------------- + phideg = (phi0 + phiminphi0) * kRad2Deg; + thetadeg = theta * kRad2Deg; + + //gLog << "MTransCelLocCam::Loc0Cam2Log; theta0deg, phi0deg, X, Y, thetadeg, phideg = " + // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", " + // << thetadeg << ", " << phideg << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Loc0LocToCam +// +// Input : (theta0, phi0) direction for the position (0,0) in the camera +// ( theta, phi) some other direction +// +// Output : (X, Y) position in the camera corresponding to (theta, phi) +// + +Bool_t MTransCelLocCam::Loc0LocToCam(Double_t theta0deg, Double_t phi0deg, + Double_t thetadeg, Double_t phideg, + Double_t &X, Double_t &Y) +{ + Double_t theta0 = theta0deg / kRad2Deg; + Double_t phi0 = phi0deg / kRad2Deg; + + Double_t theta = thetadeg / kRad2Deg; + Double_t phi = phideg / kRad2Deg; + + //-------------------------------------------- + Double_t YC = (cos(theta0)*tan(theta)*cos(phi-phi0) - sin(theta0)) + / (cos(theta0) + sin(theta0)*tan(theta)); + Double_t XC = -sin(phi-phi0) * (cos(theta0) - YC*sin(theta0)) * tan(theta); + + //-------------------------------------------- + X = XC * fDistCam; + Y = YC * fDistCam; + + + //gLog << "MTransCelLocCam::Loc0LocToCam; theta0deg, phi0deg, X, Y, thetadeg, phideg = " + // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", " + // << thetadeg << ", " << phideg << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// LocCamToLoc0 +// +// Input : ( X, Y) a position in the camera +// (theta, phi) direction for the position (X,Y) in the camera +// +// Output : ( theta0, phi0) direction for the position (0,0) in the camera +// + +Bool_t MTransCelLocCam::LocCamToLoc0(Double_t thetadeg, Double_t phideg, + Double_t X, Double_t Y, + Double_t &theta0deg, Double_t &phi0deg) +{ + Double_t theta = thetadeg / kRad2Deg; + Double_t phi = phideg / kRad2Deg; + + Double_t XC = X / fDistCam; + Double_t YC = Y / fDistCam; + + //-------------------------------------------- + + Double_t theta0; + Double_t phiminphi0; + + Double_t tant1 = 0.0; + Double_t tant2 = 0.0; + + Double_t ip; + + // calculate tan(theta0) + // dummy loop to avoid a 'go to' + for (Int_t i=0; i<1; i++) + { + if ( tan(theta) == 0.0 ) + { + if (XC != 0.0) + { + // this can never occur + gLog << "MTransCelLocCam::LocCam2Loc0; thetadeg, XC = " << thetadeg + << ", " << XC << endl; + return kFALSE; + } + tant1 = -YC; + + theta0 = TMath::Pi() * + modf( (atan(tant1)+TMath::Pi()) / TMath::Pi(), &ip ); + phiminphi0 = 0.0; + + break; + } + + Double_t a = 1.0 + XC*XC - YC*YC*tan(theta)*tan(theta); + Double_t b = 2.0 * YC * (1.0 + tan(theta)*tan(theta)); + Double_t c = XC*XC + YC*YC - tan(theta)*tan(theta); + + if (a == 0.0) + { + if (b == 0.0) + { + // this can never occur + gLog << "MTransCelLocCam::LocCam2Loc0; a, b = " << a << ", " + << b << endl; + return kFALSE; + } + tant1 = -c/b; + } + else + { + Double_t discr = b*b - 4.0*a*c; + if (discr < 0.0) + { + gLog << "MTransCelLocCam::LocCam2Loc0; discr = " << discr << endl; + return kFALSE; + } + + // two solutions for tan(theta0) + tant1 = (-b + sqrt(discr)) / (2.0*a); + tant2 = (-b - sqrt(discr)) / (2.0*a); + + if (fabs(tant1-tant2) < 1.e-5) + tant2 = 0.0; + } + + // define the sign of tan(theta0) and + // reject the solution with the wrong sign + + if ( fabs(thetadeg - 90.0) > 45.0 ) + { + Double_t sig = TMath::Sign(1.0, cos(theta)); + if ( tant1 != 0.0 && TMath::Sign(1.0, tant1) != sig ) + tant1 = 0.0; + + if ( tant2 != 0.0 && TMath::Sign(1.0, tant2) != sig ) + tant2 = 0.0; + } + else + { + // require sign of cos(phi-phi0) to be > 0 + if ( tant1 != 0.0 && + TMath::Sign( 1.0, (tant1+YC)/ ((1.0-YC*tant1)*tan(theta)) ) != 1.0 ) + tant1 = 0.0; + + if ( tant2 != 0.0 && + TMath::Sign( 1.0, (tant2+YC)/ ((1.0-YC*tant2)*tan(theta)) ) != 1.0 ) + tant2 = 0.0; + } + + if (tant1 != 0.0 && tant2 != 0.0) + { + gLog << "MTransCelLocCam::LocCam2Loc0; there are 2 solutions for tan(theta0), tant1, tant2 = " + << tant1 << ", " << tant2 << endl; + } + + if (tant1 != 0.0) + theta0 = TMath::Pi() * + modf( (atan(tant1)+TMath::Pi()) / TMath::Pi(), &ip ); + else if (tant2 != 0.0) + theta0 = TMath::Pi() * + modf( (atan(tant2)+TMath::Pi()) / TMath::Pi(), &ip ); + else + { + theta0 = theta; + gLog << "MTransCelLocCam::LocCam2Loc0; there is no solution for tan(theta0)" + << endl; + } + + Double_t sip = -XC; + Double_t cop = sin(theta0) + YC*cos(theta0); + Double_t cot = cos(theta0) - YC*sin(theta0); + Double_t sig = TMath::Sign(1.0, cot*tan(theta)); + + phiminphi0 = TMath::ATan2(sig*sip, sig*cop); + } // end of dummy loop + + //-------------------------------------------- + phi0deg = (phi - phiminphi0) * kRad2Deg; + theta0deg = theta0 * kRad2Deg; + + //gLog << "MTransCelLocCam::LocCamToLoc0; theta0deg, phi0deg, X, Y, thetadeg, phideg = " + // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", " + // << thetadeg << ", " << phideg << endl; + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Cel0CamToCel +// +// Input : (dec0, h0) direction for the position (0,0) in the camera +// ( X, Y) a position in the camera +// +// Output : (dec, h) direction for the position (X,Y) in the camera +// + +Bool_t MTransCelLocCam::Cel0CamToCel(Double_t dec0deg, Double_t h0hour, + Double_t X, Double_t Y, + Double_t &decdeg, Double_t &hhour) +{ + //-------------------------------------------- + + // transform celestial coordinates ( dec0, h0) + // into local coordinates (theta0, phi0) + Double_t theta0deg; + Double_t phi0deg; + CelToLoc(dec0deg, h0hour, theta0deg, phi0deg); + + // get the local coordinates (theta, phi) + // for the position (X, Y) in the camera + Double_t thetadeg; + Double_t phideg; + Loc0CamToLoc(theta0deg, phi0deg, X, Y, thetadeg, phideg); + + // transform local coordinates (theta, phi) + // into celestial coordinates ( dec, h) + LocToCel(thetadeg, phideg, decdeg, hhour); + + //-------------------------------------------- + + + //gLog << "MTransCelLocCam::Cel0CamToCel; theta0deg, phi0deg, X, Y, thetadeg, phideg = " + // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", " + // << thetadeg << ", " << phideg << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Cel0CelToCam +// +// Input : (dec0, h0) direction for the position (0,0) in the camera +// (dec, h) some other direction +// +// Output : (X, Y) position in the camera corresponding to (dec, h) +// + +Bool_t MTransCelLocCam::Cel0CelToCam(Double_t dec0deg, Double_t h0hour, + Double_t decdeg, Double_t hhour, + Double_t &X, Double_t &Y) +{ + //-------------------------------------------- + + // transform celestial coordinates ( dec0, h0) + // into local coordinates (theta0, phi0) + Double_t theta0deg; + Double_t phi0deg; + CelToLoc(dec0deg, h0hour, theta0deg, phi0deg); + + // transform celestial coordinates ( dec, h) + // into local coordinates (theta, phi) + Double_t thetadeg; + Double_t phideg; + CelToLoc(decdeg, hhour, thetadeg, phideg); + + // get the position (X, Y) in the camera + // from the local coordinates (theta, phi) + Loc0LocToCam(theta0deg, phi0deg, thetadeg, phideg, X, Y); + + //-------------------------------------------- + + //gLog << "MTransCelLocCam::Cel0CelToCam; theta0deg, phi0deg, X, Y, thetadeg, phideg = " + // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", " + // << thetadeg << ", " << phideg << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// CelCamToCel0 +// +// Input : ( X, Y) a position in the camera +// (dec, h) direction for the position (X,Y) in the camera +// +// Output : (dec0, h0) direction for the position (0,0) in the camera +// + +Bool_t MTransCelLocCam::CelCamToCel0(Double_t decdeg, Double_t hhour, + Double_t X, Double_t Y, + Double_t &dec0deg, Double_t &h0hour) +{ + //-------------------------------------------- + + // transform celestial coordinates ( dec, h) + // into local coordinates (theta, phi) + Double_t thetadeg; + Double_t phideg; + CelToLoc(decdeg, hhour, thetadeg, phideg); + + // get the local coordinates (theta, phi) + // for the position (0, 0) in the camera + Double_t theta0deg; + Double_t phi0deg; + LocCamToLoc0(thetadeg, phideg, X, Y, theta0deg, phi0deg); + + // transform local coordinates (theta0, phi0) + // into celestial coordinates ( dec0, h0) + LocToCel(theta0deg, phi0deg, dec0deg, h0hour); + + //-------------------------------------------- + + + //gLog << "MTransCelLocCam::CelCamToCel0; theta0deg, phi0deg, X, Y, thetadeg, phideg = " + // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", " + // << thetadeg << ", " << phideg << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// LocCamCamToLoc +// +// Input : (theta1, phi1) direction for the position (X1,Y1) in the camera +// ( X1, Y1) a position in the camera +// ( X2, Y2) another position in the camera +// +// Output : (theta2, phi2) direction for the position (X2,Y2) in the camera +// + +Bool_t MTransCelLocCam::LocCamCamToLoc(Double_t theta1deg, Double_t phi1deg, + Double_t X1, Double_t Y1, + Double_t X2, Double_t Y2, + Double_t &theta2deg, Double_t &phi2deg) +{ + if (X1 == 0.0 && Y1 == 0.0) + { + Loc0CamToLoc(theta1deg, phi1deg, X2, Y2, theta2deg, phi2deg); + return kTRUE; + } + + if (X2 == 0.0 && Y2 == 0.0) + { + LocCamToLoc0(theta1deg, phi1deg, X1, Y1, theta2deg, phi2deg); + return kTRUE; + } + + Double_t theta0deg; + Double_t phi0deg; + LocCamToLoc0(theta1deg, phi1deg, X1, Y1, theta0deg, phi0deg); + + Loc0CamToLoc(theta0deg, phi0deg, X2, Y2, theta2deg, phi2deg); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// LocCamLocToCam +// +// Input : (theta1, phi1) direction for the position (X1,Y1) in the camera +// ( X1, Y1) a position in the camera +// (theta2, phi2) another direction +// +// Output : ( X2, Y2) position corresponding to (theta2, phi2) +// + +Bool_t MTransCelLocCam::LocCamLocToCam(Double_t theta1deg, Double_t phi1deg, + Double_t X1, Double_t Y1, + Double_t theta2deg, Double_t phi2deg, + Double_t &X2, Double_t &Y2) + +{ + if (X1 == 0.0 && Y1 == 0.0) + { + Loc0LocToCam(theta1deg, phi1deg, theta2deg, phi2deg, X2, Y2); + + return kTRUE; + } + + + Double_t theta0deg; + Double_t phi0deg; + LocCamToLoc0(theta1deg, phi1deg, X1, Y1, theta0deg, phi0deg); + + Loc0LocToCam(theta0deg, phi0deg, theta2deg, phi2deg, X2, Y2); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// CelCamCamToCel +// +// Input : (dec1, h1) direction for the position (X1,Y1) in the camera +// ( X1, Y1) a position in the camera +// ( X2, Y2) another position in the camera +// +// Output : (dec2, h2) direction for the position (X2,Y2) in the camera +// + +Bool_t MTransCelLocCam::CelCamCamToCel(Double_t dec1deg, Double_t h1deg, + Double_t X1, Double_t Y1, + Double_t X2, Double_t Y2, + Double_t &dec2deg, Double_t &h2deg) +{ + if (X1 == 0.0 && Y1 == 0.0) + { + Cel0CamToCel(dec1deg, h1deg, X2, Y2, dec2deg, h2deg); + return kTRUE; + } + + if (X2 == 0.0 && Y2 == 0.0) + { + CelCamToCel0(dec1deg, h1deg, X1, Y1, dec2deg, h2deg); + return kTRUE; + } + + Double_t dec0deg; + Double_t h0deg; + CelCamToCel0(dec1deg, h1deg, X1, Y1, dec0deg, h0deg); + + Cel0CamToCel(dec0deg, h0deg, X2, Y2, dec2deg, h2deg); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// CelCamCelToCam +// +// Input : (dec1, h1) direction for the position (X1,Y1) in the camera +// ( X1, Y1) a position in the camera +// (dec2, h2) another direction +// +// Output : ( X2, Y2) position corresponding to (dec2, h2) +// + +Bool_t MTransCelLocCam::CelCamCelToCam(Double_t dec1deg, Double_t h1deg, + Double_t X1, Double_t Y1, + Double_t dec2deg, Double_t h2deg, + Double_t &X2, Double_t &Y2) + +{ + if (X1 == 0.0 && Y1 == 0.0) + { + Cel0CelToCam(dec1deg, h1deg, dec2deg, h2deg, X2, Y2); + + return kTRUE; + } + + + Double_t dec0deg; + Double_t h0deg; + CelCamToCel0(dec1deg, h1deg, X1, Y1, dec0deg, h0deg); + + Cel0CelToCam(dec0deg, h0deg, dec2deg, h2deg, X2, Y2); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// CelToLoc +// +// Input : (dec, h) celestial coordinates +// +// Output : (theta, phi) corresponding local coordinates +// + +Bool_t MTransCelLocCam::CelToLoc(Double_t decdeg, Double_t hhour, + Double_t &thetadeg, Double_t &phideg) +{ + Double_t a1 = fCLat; + Double_t a3 = -fSLat; + + Double_t dec = decdeg / kRad2Deg; + Double_t h = hhour / 24.0 * TMath::TwoPi(); + + // transform celestial coordinates ( dec, h) + // into local coordinates (theta, phi) + Double_t xB = cos(dec) * cos(h); + Double_t yB = cos(dec) * sin(h); + Double_t zB = - sin(dec); + + Double_t xA = a3*xB - a1*zB; + Double_t yA = - yB; + Double_t zA = -a1*xB - a3*zB; + + thetadeg = acos(-zA) * kRad2Deg; + phideg = TMath::ATan2(yA, xA) * kRad2Deg; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// LocToCel +// +// Input : (theta, phi) local coordinates +// +// Output : (dec, h) corresponding celestial coordinates +// + +Bool_t MTransCelLocCam::LocToCel(Double_t thetadeg, Double_t phideg, + Double_t &decdeg, Double_t &hhour) +{ + Double_t a1 = fCLat; + Double_t a3 = -fSLat; + + Double_t theta = thetadeg / kRad2Deg; + Double_t phi = phideg / kRad2Deg; + + //-------------------------------------------- + + // transform local coordinates (theta, phi) + // into celestial coordinates ( dec, h) + Double_t xA = sin(theta) * cos(phi); + Double_t yA = sin(theta) * sin(phi); + Double_t zA = - cos(theta); + + Double_t xB = a3*xA - a1*zA; + Double_t yB = - yA; + Double_t zB = -a1*xA - a3*zA; + + Double_t dec = asin(-zB); + Double_t h = TMath::ATan2(yB, xB); + + //-------------------------------------------- + decdeg = dec * kRad2Deg; + hhour = h * 24.0 / TMath::TwoPi(); + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// PlotGridAtDec0H0 +// +// set the celestial coordinates of the camera center +// and plot the lines of constant (Theta, Phi) +// and the lines of constant (Dec, H ) +// + +Bool_t MTransCelLocCam::PlotGridAtDec0H0(TString name, + Double_t dec0deg, Double_t h0hour) +{ + Double_t theta0deg; + Double_t phi0deg; + CelToLoc(dec0deg, h0hour, theta0deg, phi0deg); + + PlotGrid(name, theta0deg, phi0deg, dec0deg, h0hour); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// PlotGridAtTheta0Phi0 +// +// set the local coordinates of the camera center +// and plot the lines of constant (Theta, Phi) +// and the lines of constant (Dec, H ) +// + +Bool_t MTransCelLocCam::PlotGridAtTheta0Phi0(TString name, + Double_t theta0deg, Double_t phi0deg) +{ + Double_t dec0deg; + Double_t h0hour; + LocToCel(theta0deg, phi0deg, dec0deg, h0hour); + + PlotGrid(name, theta0deg, phi0deg, dec0deg, h0hour); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// SetGridParameters +// +// set the binning for the grid (fGridBinning) +// set the binning along the lines (fGridFineBin) +// + +Bool_t MTransCelLocCam::SetGridParameters( + Double_t gridbinning, Double_t gridfinebin) +{ + fGridBinning = gridbinning; + fGridFineBin = gridfinebin; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// PlotGrid +// +// - plot the lines of constant (Theta, Phi) +// with the camera center at (Theta0, Phi0) +// - and the lines of constant (Dec, H ) +// with the camera center at (Dec0, H0 ) +// + +Bool_t MTransCelLocCam::PlotGrid(TString name, + Double_t theta0deg, Double_t phi0deg, + Double_t dec0deg, Double_t h0hour) +{ + Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ; + //gLog << "mmtodeg = " << mmtodeg << endl; + + //-------------------------------------------------------------------- + + TCanvas *c1 = new TCanvas(name, name, 0, 0, 300, 600); + + //gROOT->Reset(); + gStyle->SetCanvasColor(0); + gStyle->SetCanvasBorderMode(0); + gStyle->SetPadBorderMode(0); + gStyle->SetFrameBorderMode(0); + gStyle->SetOptStat(0000000); + gStyle->SetPalette(1); + + c1->Divide(1,2); + c1->SetBorderMode(0); + + Double_t xlo = -534.0 * mmtodeg; + Double_t xup = 534.0 * mmtodeg; + + Double_t ylo = -534.0 * mmtodeg; + Double_t yup = 534.0 * mmtodeg; + + TString nam1 = name; + nam1 += "_:_Theta-Phi"; + TString tit1 = name; + tit1 += "_:_Theta-Phi-lines"; + TH2D *plot1 = new TH2D(nam1, tit1, 1, xlo, xup, 1, ylo, yup); + plot1->GetXaxis()->SetTitle("x [deg]"); + plot1->GetYaxis()->SetTitle("y [deg]"); + + + TString nam2 = name; + nam2 += "_:_Dec-Hour"; + TString tit2 = name; + tit2 += "_:_Dec-Hour-lines"; + TH2D *plot2 = new TH2D(nam2, tit2, 1, xlo, xup, 1, ylo, yup); + plot2->GetXaxis()->SetTitle("x [deg]"); + plot2->GetYaxis()->SetTitle("y [deg]"); + + + c1->cd(1); + plot1->Draw(); + //delete plot1; + + c1->cd(2); + plot2->Draw(); + //delete plot2; + + TGraph *graph1; + TLatex *pix; + + Char_t tit[100]; + Double_t xtxt; + Double_t ytxt; + + Double_t xx; + Double_t yy; + + Double_t gridbinning = fGridBinning; + Double_t gridfinebin = fGridFineBin; + Int_t numgridlines = (Int_t)(4.0/gridbinning); + + + //-------------------------------------------------------------------- + // Theta-Phi lines + + // direction for the center of the camera + Double_t theta0 = theta0deg; + Double_t phi0 = phi0deg; + + + //----- lines for fixed theta ------------------------------------ + + const Int_t Ntheta = numgridlines; + Double_t theta[Ntheta]; + Double_t dtheta = gridbinning; + + Int_t nphi = (Int_t)(4.0/gridfinebin); + const Int_t Nphi = nphi+1; + Double_t phi[Nphi]; + Double_t dphi = gridfinebin/sin(theta0/180.0*3.1415926); + if ( dphi > 360.0/((Double_t)(Nphi-1)) ) + dphi = 360.0/((Double_t)(Nphi-1)); + + for (Int_t j=0; jcd(1); + graph1 = new TGraph(Nphi,x,y); + graph1->SetLineColor(j+1); + graph1->SetLineStyle(1); + graph1->SetMarkerColor(j+1); + graph1->SetMarkerSize(.2); + graph1->SetMarkerStyle(20); + graph1->Draw("PL"); + //delete graph1; + + sprintf(tit,"Theta = %6.2f", theta[j]); + xtxt = xlo + (xup-xlo)*0.1; + ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)j) *(yup-ylo)/20.0; + pix = new TLatex(xtxt, ytxt, tit); + pix->SetTextColor(j+1); + pix->SetTextSize(.03); + pix->Draw(""); + //delete pix; + + } + + //----- lines for fixed phi ------------------------------------ + + Int_t ntheta1 = (Int_t)(4.0/gridfinebin); + const Int_t Ntheta1 = ntheta1; + Double_t theta1[Ntheta1]; + Double_t dtheta1 = gridfinebin; + + const Int_t Nphi1 = numgridlines; + Double_t phi1[Nphi1]; + Double_t dphi1 = gridbinning/sin(theta0/180.0*3.1415926); + if ( dphi1 > 360.0/((Double_t)Nphi1) ) + dphi1 = 360.0/((Double_t)Nphi1); + + for (Int_t j=0; jcd(1); + graph1 = new TGraph(count,x1,y1); + graph1->SetLineColor(k+1); + graph1->SetLineStyle(2); + graph1->SetMarkerColor(k+1); + graph1->SetMarkerSize(.2); + graph1->SetMarkerStyle(20); + graph1->Draw("PL"); + //delete graph1; + + sprintf(tit,"Phi = %6.2f", phi1[k]); + Double_t xtxt = xlo + (xup-xlo)*0.75; + Double_t ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)k) *(yup-ylo)/20.0; + pix = new TLatex(xtxt, ytxt, tit); + pix->SetTextColor(k+1); + pix->SetTextSize(.03); + pix->Draw(""); + //delete pix; + + } + + c1->cd(1); + sprintf(tit,"Theta0 = %6.2f Phi0 = %6.2f", theta0, phi0); + xtxt = xlo + (xup-xlo)*0.05; + ytxt = ylo + (yup-ylo)*0.92; + pix = new TLatex(xtxt, ytxt, tit); + pix->SetTextColor(1); + pix->SetTextSize(.06); + pix->Draw(""); + //delete pix; + + sprintf(tit," [deg] [deg]"); + xtxt = xlo + (xup-xlo)*0.05; + ytxt = ylo + (yup-ylo)*0.86; + pix = new TLatex(xtxt, ytxt, tit); + pix->SetTextColor(1); + pix->SetTextSize(.06); + pix->Draw(""); + //delete pix; + + + //-------------------------------------------------------------------- + // Dec-Hour lines + + // direction for the center of the camera + Double_t dec0 = dec0deg; + Double_t h0 = h0hour; + //trans.LocToCel(theta0, phi0, dec0, h0); + + + //----- lines for fixed dec ------------------------------------ + + const Int_t Ndec = numgridlines; + Double_t dec[Ndec]; + Double_t ddec = gridbinning; + + Int_t nh = (Int_t)(4.0/gridfinebin); + const Int_t Nh = nh+1; + Double_t h[Nh]; + Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926); + if ( dh > 360.0/((Double_t)(Nh-1)) ) + dh = 360.0/((Double_t)(Nh-1)); + + + for (Int_t j=0; j 90.0) continue; + + for (Int_t k=0; kcd(2); + graph1 = new TGraph(Nh,xh,yh); + graph1->SetLineColor(j+1); + graph1->SetLineStyle(1); + graph1->SetMarkerColor(j+1); + graph1->SetMarkerSize(.2); + graph1->SetMarkerStyle(20); + graph1->Draw("PL"); + //delete graph1; + + sprintf(tit,"Dec = %6.2f", dec[j]); + xtxt = xlo + (xup-xlo)*0.1; + ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)j) *(yup-ylo)/20.0; + pix = new TLatex(xtxt, ytxt, tit); + pix->SetTextColor(j+1); + pix->SetTextSize(.03); + pix->Draw(""); + //delete pix; + + } + + //----- lines for fixed hour angle ------------------------------------ + + Int_t ndec1 = (Int_t)(4.0/gridfinebin); + const Int_t Ndec1 = ndec1; + Double_t dec1[Ndec1]; + Double_t ddec1 = gridfinebin; + + const Int_t Nh1 = numgridlines; + Double_t h1[Nh1]; + Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926); + if ( dh1 > 360.0/((Double_t)Nh1) ) + dh1 = 360.0/((Double_t)Nh1); + + + for (Int_t j=0; j 90.0) continue; + + Double_t hh0 = h0 *24.0/360.0; + Double_t hhx = h1[k]*24.0/360.0; + Cel0CelToCam(dec0, hh0, dec1[j], hhx, + xx, yy); + xd[count] = xx * mmtodeg; + yd[count] = yy * mmtodeg; + + //gLog << "dec0, h0 = " << dec0 << ", " << h0 + // << " : dec1, h1, xd, yd = " << dec1[j] << ", " + // << h1[k] << "; " << xd[count] << ", " << yd[count] << endl; + + count++; + } + + c1->cd(2); + graph1 = new TGraph(count,xd,yd); + graph1->SetLineColor(k+1); + graph1->SetLineStyle(2); + graph1->SetMarkerColor(k+1); + graph1->SetMarkerSize(.2); + graph1->SetMarkerStyle(20); + graph1->Draw("PL"); + //delete graph1; + + sprintf(tit,"H = %6.2f", h1[k]); + Double_t xtxt = xlo + (xup-xlo)*0.75; + Double_t ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)k) *(yup-ylo)/20.0; + pix = new TLatex(xtxt, ytxt, tit); + pix->SetTextColor(k+1); + pix->SetTextSize(.03); + pix->Draw(""); + //delete pix; + + } + + c1->cd(2); + sprintf(tit,"Dec0 = %6.2f H0 = %6.2f", dec0, h0); + xtxt = xlo + (xup-xlo)*0.05; + ytxt = ylo + (yup-ylo)*0.92; + pix = new TLatex(xtxt, ytxt, tit); + pix->SetTextColor(1); + pix->SetTextSize(.06); + pix->Draw(""); + //delete pix; + + sprintf(tit," [deg] [deg]"); + xtxt = xlo + (xup-xlo)*0.05; + ytxt = ylo + (yup-ylo)*0.86; + pix = new TLatex(xtxt, ytxt, tit); + pix->SetTextColor(1); + pix->SetTextSize(.06); + pix->Draw(""); + //delete pix; + + return kTRUE; +} + +//=========================================================================== + + + + + + + + + + Index: trunk/MagicSoft/Mars/mastro/MTransCelLocCam.h =================================================================== --- trunk/MagicSoft/Mars/mastro/MTransCelLocCam.h (revision 4674) +++ trunk/MagicSoft/Mars/mastro/MTransCelLocCam.h (revision 4674) @@ -0,0 +1,103 @@ +#ifndef MARS_MTransCelLocCam +#define MARS_MTransCelLocCam + +#ifndef ROOT_TROOT +#include +#endif + +#ifndef MARS_MAGIC +#include "MAGIC.h" +#endif + +class MGeomCam; +class MObservatory; + +class MTransCelLocCam +{ +private: + Double_t fDistCam; // distance camera reflector [mm] + Double_t fCLat; // cos(Lat) + Double_t fSLat; // sin(Lat) + + Double_t fGridBinning; // degrees + Double_t fGridFineBin; // degrees + +public: + MTransCelLocCam(const MGeomCam &cam, const MObservatory &obs); + + Bool_t Loc0CamToLoc(Double_t theta0deg, Double_t phi0deg, + Double_t X, Double_t Y, + Double_t &thetadeg, Double_t &phideg); + + Bool_t Loc0LocToCam(Double_t theta0deg, Double_t phi0deg, + Double_t thetadeg, Double_t phideg, + Double_t &X, Double_t &Y); + + Bool_t LocCamToLoc0(Double_t thetadeg, Double_t phideg, + Double_t X, Double_t Y, + Double_t &theta0deg, Double_t &phi0deg); + + Bool_t Cel0CamToCel(Double_t dec0deg, Double_t h0hour, + Double_t X, Double_t Y, + Double_t &decdeg, Double_t &hhour); + + Bool_t Cel0CelToCam(Double_t dec0deg, Double_t h0hour, + Double_t decdeg, Double_t hhour, + Double_t &X, Double_t &Y); + + Bool_t CelCamToCel0(Double_t decdeg, Double_t hhour, + Double_t X, Double_t Y, + Double_t &dec0deg, Double_t &h0hour); + + Bool_t LocCamCamToLoc(Double_t theta1deg, Double_t phi1deg, + Double_t X1, Double_t Y1, + Double_t X2, Double_t Y2, + Double_t &theta2deg, Double_t &phi2deg); + + Bool_t LocCamLocToCam(Double_t theta1deg, Double_t phi1deg, + Double_t X1, Double_t Y1, + Double_t theta2deg, Double_t phi2deg, + Double_t &X2, Double_t &Y2); + + Bool_t CelCamCamToCel(Double_t dec1deg, Double_t h1hour, + Double_t X1, Double_t Y1, + Double_t X2, Double_t Y2, + Double_t &dec2deg, Double_t &h2hour); + + Bool_t CelCamCelToCam(Double_t dec1deg, Double_t h1hour, + Double_t X1, Double_t Y1, + Double_t dec2deg, Double_t h2hour, + Double_t &X2, Double_t &Y2); + + Bool_t CelToLoc(Double_t decdeg, Double_t hhour, + Double_t &thetadeg, Double_t &phideg); + + Bool_t LocToCel(Double_t thetadeg, Double_t phideg, + Double_t &decdeg, Double_t &hhour); + + Bool_t PlotGridAtDec0H0(TString name, + Double_t dec0deg, Double_t h0hour); + Bool_t PlotGridAtTheta0Phi0(TString name, + Double_t theta0deg, Double_t phi0deg); + + Bool_t SetGridParameters(Double_t gridbinning, Double_t gridfinebin); + + Bool_t PlotGrid(TString name, + Double_t theta0deg, Double_t phi0deg, + Double_t dec0deg, Double_t h0hour); + + ClassDef(MTransCelLocCam, 0) // Class for transformations between Celestial, Local and Camera coordinates +}; + +#endif + + + + + + + + + + +