Changeset 6076 for trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc

Timestamp:

01/28/05 10:50:20 (20 years ago)

Author:

moralejo

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc (modified) (10 diffs)

trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc

@@ -17 +17 @@
 !
 !   Author(s): Thomas Bretz 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
+!              Abelardo Moralejo 1/2005 <mailto:moralejo@pd.infn.it>
 !
 !   Copyright: MAGIC Software Development, 2000-2004
@@ -27 +28 @@
 // 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 +46 @@
 //   MObservatory
 //   MGeomCam
+//   MTime
+//   [MSourcePos] (of type MPointingPos)
 //
 // Output Container:
@@ -45 +53 @@
 //
 // 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 +87 @@
 //
 MSrcPosCalc::MSrcPosCalc(const char *name, const char *title)
-    : fX(0), fY(0)
 {
     fName  = name  ? name  : "MSrcPosCalc";
@@ -82 +94 @@
 // --------------------------------------------------------------------------
 //
-// 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 +127 @@
     if (!fPointPos)
     {
-        *fLog << err << "MPointPos not found... aborting." << endl;
+        *fLog << err << "MPointingPos not found... aborting." << endl;
         return kFALSE;
     }
@@ -114 +145 @@
     }
 
-    fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam");
-    if (!fSrcPos)
-        return kFALSE;
-
     return kTRUE;
 }
@@ -123 +150 @@
 // --------------------------------------------------------------------------
 //
-// 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 +188 @@
 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;
+}