Index: trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc (revision 6048) +++ trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc (revision 6076) @@ -17,4 +17,5 @@ ! ! Author(s): Thomas Bretz 3/2004 +! Abelardo Moralejo 1/2005 ! ! Copyright: MAGIC Software Development, 2000-2004 @@ -27,12 +28,17 @@ // MSrcPosCalc // -// Calculate the current source position in the camera from the (already -// corrected) local pointing position of the telescope (MPointingPos) by -// derotating the position (given in x and y coordinates in the camera -// plane) at culmination time by the current rotation angle of the -// starfield dtermined from MObservatory and MPointingPos +// Calculate the current source position in the camera from the (possibly already +// corrected, by starguider) J2000 sky coordinates of the camera center (contained +// in MPointingPos), and the source J2000 sky coordinates contained in MSourcePos +// (of type MPointingPos as well). If no MSourcePos is found in the parameter list, +// source position is assumed to be the same for all events, that specified in +// MSrcPosCam (if it already existed in the parameter list), or (0,0), the center +// of the camera, if no MSrcPosCam was present in the parameter list. In both cases, +// no calculation is necessary and then the PreProcess returns kSKIP so that the task +// is removed from the task list. // // The conversion factor between the camera plain (mm) and the -// sky (deg) is taken from MGeomCam. +// sky (deg) is taken from MGeomCam. The time is taken from MTime, and the coordinates +// of the observatory from MObservatory. // // Input Container: @@ -40,4 +46,6 @@ // MObservatory // MGeomCam +// MTime +// [MSourcePos] (of type MPointingPos) // // Output Container: @@ -45,8 +53,13 @@ // // To be done: -// - wobble mode missing +// - wobble mode missing /// NOTE, A. Moralejo: I see no need for special code +// // - a switch between using sky-coordinates and time or local-coordinates -// from MPointingPos for determin the rotation angle +// from MPointingPos for determine the rotation angle +///// NOTE, A. Moralejo: the above is not possible if MAstroSky2Local does not +///// account for precession and nutation. +// // - the center of rotation need not to be the camera center +///// NOTE, A. Moralejo: I don't see the need for special code either. // ////////////////////////////////////////////////////////////////////////////// @@ -74,5 +87,4 @@ // MSrcPosCalc::MSrcPosCalc(const char *name, const char *title) - : fX(0), fY(0) { fName = name ? name : "MSrcPosCalc"; @@ -82,8 +94,27 @@ // -------------------------------------------------------------------------- // -// Search and if necessary create MSrcPosCam in the parameter list +// Search and if necessary create MSrcPosCam in the parameter list. Search MSourcePos. +// If not found, do nothing else, and skip the task. If MSrcPosCam did not exist +// before and has been created here, it will contain as source position the camera +// center (0,0). +// In the case that MSourcePos is found, go ahead in searching the rest of necessary +// containers. The source position will be calculated for each event in Process. // Int_t MSrcPosCalc::PreProcess(MParList *pList) { + fSrcPosCam = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam"); + if (!fSrcPosCam) + return kFALSE; + + + fSourcePos = (MPointingPos*)pList->FindObject("MSourcePos", "MPointingPos"); + if (!fSourcePos) + { + *fLog << warn << "MSourcePos [MPointPos] not found... The source position" << endl; + *fLog << warn << "set in MSrcPosCam (camera center if not set explicitely) will" << endl; + *fLog << warn << "be left unchanged, same for all events." << endl; + return kSKIP; + } + fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); if (!fGeom) @@ -96,5 +127,5 @@ if (!fPointPos) { - *fLog << err << "MPointPos not found... aborting." << endl; + *fLog << err << "MPointingPos not found... aborting." << endl; return kFALSE; } @@ -114,8 +145,4 @@ } - fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam"); - if (!fSrcPos) - return kFALSE; - return kTRUE; } @@ -123,46 +150,34 @@ // -------------------------------------------------------------------------- // -// Search the parameter list for MObservatory and MTime -// -Bool_t MSrcPosCalc::ReInit(MParList *pList) -{ - if (fX!=0 || fY!=0) - return kTRUE; - - *fLog << warn << "fX==0 && fY==0: Using arbitrary source position!" << endl; - - // - // This is a workaround for current crab misspointing - DO NOT USE! - // It will be removed soon! - // - const MRawRunHeader *h = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!h) - { - *fLog << err << "MRawRunHeader not found... aborting." << endl; - return kFALSE; - } - - const MTime &t = h->GetRunStart(); - - const Double_t rho = fPointPos->RotationAngle(*fObservatory, t); - - // Calculate x, y of Zeta Tauri - - Double_t tm = t.GetAxisTime(); - - Double_t x = 1.7267e6-6.03285e-3*tm; // [mm] - Double_t y = -189.823+974.908*exp(-52.3083*(tm*1e-5-2861.5)); // [mm] - - const Double_t cs = TMath::Cos(rho-fDrho); - const Double_t sn = TMath::Sin(rho-fDrho); - - const Double_t dx = fR*cs; - const Double_t dy = fR*sn; - - fSrcPos->SetXY(x-dx, y-dy); - - *fLog << dbg << t << " - Position: x=" << x-dx << "mm y=" << y-dy << "mm" << 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) +// +TVector2 MSrcPosCalc::CalcXYinCamera(const MVector3 &pos0, const MVector3 &pos) const +{ + const Double_t theta0 = pos0.Theta(); + const Double_t phi0 = pos0.Phi(); + + const Double_t theta = pos.Theta(); + const Double_t phi = pos.Phi(); + + //-------------------------------------------- + + // pos0[3] = TMath::Cos(theta0) + + const Double_t YC0 = TMath::Cos(theta0)*TMath::Tan(theta)*TMath::Cos(phi-phi0) - TMath::Sin(theta0); + const Double_t YC1 = TMath::Cos(theta0) + TMath::Sin(theta0)*TMath::Tan(theta); + const Double_t YC = YC0 / YC1; + + //-------------------------------------------- + + const Double_t XC0 = TMath::Cos(theta0) - YC*TMath::Sin(theta0); + const Double_t XC = -TMath::Sin(phi-phi0) * TMath::Tan(theta) * XC0; + + //-------------------------------------------- + return TVector2(XC, YC); } @@ -173,20 +188,55 @@ Int_t MSrcPosCalc::Process() { - if (fX==0 && fY==0) - return kTRUE; - - // rotate the source position by the current rotation angle - const Double_t rho = fPointPos->RotationAngle(*fObservatory, *fTime); - - // Define source position in the camera plain - TVector2 v(fX, fY); - v=v.Rotate(rho); - - // Convert coordinates into camera plain (mm) - v *= 1./fGeom->GetConvMm2Deg(); - - // Set current source position - fSrcPos->SetXY(v); - - return kTRUE; -} + // *fLog << dbg << "Camera center : Zd=" << fPointPos->GetZd() << " Az=" << fPointPos->GetAz() << endl; + // *fLog << dbg << "Camera center : RA=" << fPointPos->GetRa() << " Dec=" << fPointPos->GetDec() << endl; + // *fLog << dbg << "MJD: " << fTime->GetMjd() << ", time [ms]: " << fTime->GetTime() << endl; + + MVector3 pos, pos0; // pos: source position; pos0: camera center + + if (fSourcePos) + { + // Set Sky coordinates of source, taken from container "MSourcePos" of type MPointingPos. The sky + // coordinates must be J2000, as the sky coordinates of the camera center that we get from the container + // "MPointingPos" filled by the Drive. + + pos.SetRaDec(fSourcePos->GetRaRad(), fSourcePos->GetDecRad()); + + // *fLog << dbg << "Sky position of source: Ra=" << fSourcePos->GetRa() << + // " Dec=" << fSourcePos->GetDec() << endl; + } + + // Convert sky coordinates of source to local coordinates. Warning! These are not the "true" local + // coordinates, since this transformation ignores precession and nutation effects. + pos *= MAstroSky2Local(*fTime, *fObservatory); + + + // Set sky coordinates of camera center in pos0, and then convert to local. Same comment as above. These + // coordinates differ from the true local coordinates of the camera center that one could get from + // "MPointingPos", calculated by the Drive: the reason is that the Drive takes into account precession + // and nutation corrections, while MAstroSky2Local (as of Jan 27 2005 at least) does not. Since we just + // want to get the source position on the camera from the local coordinates of the center and the source, + // it does not matter that the coordinates contained in pos and pos0 ignore precession and nutation... + // since the shift would be the same in both cases. What would be wrong is to set in pos0 directly the + // local coordinates found in MPointingPos! + + pos0.SetRaDec(fPointPos->GetRaRad(), fPointPos->GetDecRad()); + pos0 *= MAstroSky2Local(*fTime, *fObservatory); + + // *fLog << dbg << "From MAstroSky2Local, without precession and nutation corrections:" << endl; + // *fLog << dbg << "- Camera center (2) from RA,dec and time : Zd=" << pos0.Theta()*180./TMath::Pi() + // << " Az=" << pos0.Phi()*180./TMath::Pi() << endl; + // *fLog << dbg << "- Local position of source: Zd=" << pos.Theta()*TMath::RadToDeg() << " Az=" << pos.Phi()*TMath::RadToDeg() << endl; + + + // Calculate source position in camera, and convert to mm: + TVector2 v = CalcXYinCamera(pos0, pos)*fGeom->GetCameraDist()*1000; + + fSrcPosCam->SetXY(v); + + // v *= fGeom->GetConvMm2Deg(); + // *fLog << dbg << "X=" << v.X() << " deg, Y=" << v.Y() << " deg" << endl; + // *fLog << *fTime << endl; + // *fLog << endl; + + return kTRUE; +}