Changeset 6076 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
01/28/05 10:50:20 (20 years ago)
Author:
moralejo
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r6075 r6076  
    2020
    2121                                                 -*-*- END OF LINE -*-*-
     22 2005/01/28 Abelardo Moralejo
     23
     24   * mpointing/MSrcPosCalc.cc
     25     - updated. Make it work as desired: obtain for each event the x,y
     26       position on the camera of a source whose celestial coordinates
     27       (J2000) have been set by the user in the container "MSourcePos"
     28       of type MPointingPos, added to the parameter list. If the
     29       container MSourcePos is not found, the center of the camera (or
     30       other fixed position x,y set in the parameter list by the user -
     31       in MSrcPosCam) is used as source position for all the events in
     32       the loop.
     33
    2234 2005/01/28 Thomas Bretz
    2335
  • trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc

    r3666 r6076  
    1717!
    1818!   Author(s): Thomas Bretz 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
     19!              Abelardo Moralejo 1/2005 <mailto:moralejo@pd.infn.it>
    1920!
    2021!   Copyright: MAGIC Software Development, 2000-2004
     
    2728// MSrcPosCalc
    2829//
    29 // Calculate the current source position in the camera from the (already
    30 // corrected) local pointing position of the telescope (MPointingPos) by
    31 // derotating the position (given in x and y coordinates in the camera
    32 // plane) at culmination time by the current rotation angle of the
    33 // starfield dtermined from MObservatory and MPointingPos
     30// Calculate the current source position in the camera from the (possibly already
     31// corrected, by starguider) J2000 sky coordinates of the camera center (contained
     32// in MPointingPos), and the source J2000 sky coordinates contained in MSourcePos
     33// (of type MPointingPos as well). If no MSourcePos is found in the parameter list,
     34// source position is assumed to be the same for all events, that specified in
     35// MSrcPosCam (if it already existed in the parameter list), or (0,0), the center
     36// of the camera, if no MSrcPosCam was present in the parameter list. In both cases,
     37// no calculation is necessary and then the PreProcess returns kSKIP so that the task
     38// is removed from the task list.
    3439//
    3540// The conversion factor between the camera plain (mm) and the
    36 // sky (deg) is taken from MGeomCam.
     41// sky (deg) is taken from MGeomCam. The time is taken from MTime, and the coordinates
     42// of the observatory from MObservatory.
    3743//
    3844// Input Container:
     
    4046//   MObservatory
    4147//   MGeomCam
     48//   MTime
     49//   [MSourcePos] (of type MPointingPos)
    4250//
    4351// Output Container:
     
    4553//
    4654// To be done:
    47 //   - wobble mode missing
     55//   - wobble mode missing   /// NOTE, A. Moralejo: I see no need for special code
     56//
    4857//   - a switch between using sky-coordinates and time or local-coordinates
    49 //     from MPointingPos for determin the rotation angle
     58//     from MPointingPos for determine the rotation angle
     59/////  NOTE, A. Moralejo: the above is not possible if MAstroSky2Local does not
     60/////  account for precession and nutation.
     61//
    5062//   - the center of rotation need not to be the camera center
     63/////  NOTE, A. Moralejo: I don't see the need for special code either.
    5164//
    5265//////////////////////////////////////////////////////////////////////////////
     
    7487//
    7588MSrcPosCalc::MSrcPosCalc(const char *name, const char *title)
    76     : fX(0), fY(0)
    7789{
    7890    fName  = name  ? name  : "MSrcPosCalc";
     
    8294// --------------------------------------------------------------------------
    8395//
    84 // Search and if necessary create MSrcPosCam in the parameter list
     96// Search and if necessary create MSrcPosCam in the parameter list. Search MSourcePos.
     97// If not found, do nothing else, and skip the task. If MSrcPosCam did not exist
     98// before and has been created here, it will contain as source position the camera
     99// center (0,0).
     100// In the case that MSourcePos is found, go ahead in searching the rest of necessary
     101// containers. The source position will be calculated for each event in Process.
    85102//
    86103Int_t MSrcPosCalc::PreProcess(MParList *pList)
    87104{
     105    fSrcPosCam = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam");
     106    if (!fSrcPosCam)
     107        return kFALSE;
     108
     109
     110    fSourcePos = (MPointingPos*)pList->FindObject("MSourcePos", "MPointingPos");
     111    if (!fSourcePos)
     112    {
     113        *fLog << warn << "MSourcePos [MPointPos] not found... The source position" << endl;
     114        *fLog << warn << "set in MSrcPosCam (camera center if not set explicitely) will" << endl;
     115        *fLog << warn << "be left unchanged, same for all events." << endl;
     116        return kSKIP;
     117    }
     118
    88119    fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
    89120    if (!fGeom)
     
    96127    if (!fPointPos)
    97128    {
    98         *fLog << err << "MPointPos not found... aborting." << endl;
     129        *fLog << err << "MPointingPos not found... aborting." << endl;
    99130        return kFALSE;
    100131    }
     
    114145    }
    115146
    116     fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam");
    117     if (!fSrcPos)
    118         return kFALSE;
    119 
    120147    return kTRUE;
    121148}
     
    123150// --------------------------------------------------------------------------
    124151//
    125 // Search the parameter list for MObservatory and MTime
    126 //
    127 Bool_t MSrcPosCalc::ReInit(MParList *pList)
    128 {
    129     if (fX!=0 || fY!=0)
    130         return kTRUE;
    131 
    132     *fLog << warn << "fX==0 && fY==0: Using arbitrary source position!" << endl;
    133 
    134     //
    135     // This is a workaround for current crab misspointing - DO NOT USE!
    136     // It will be removed soon!
    137     //
    138     const MRawRunHeader *h = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    139     if (!h)
    140     {
    141         *fLog << err << "MRawRunHeader not found... aborting." << endl;
    142         return kFALSE;
    143     }
    144 
    145     const MTime &t = h->GetRunStart();
    146 
    147     const Double_t rho = fPointPos->RotationAngle(*fObservatory, t);
    148 
    149     // Calculate x, y of Zeta Tauri
    150 
    151     Double_t tm = t.GetAxisTime();
    152 
    153     Double_t x = 1.7267e6-6.03285e-3*tm; // [mm]
    154     Double_t y = -189.823+974.908*exp(-52.3083*(tm*1e-5-2861.5)); // [mm]
    155 
    156     const Double_t cs = TMath::Cos(rho-fDrho);
    157     const Double_t sn = TMath::Sin(rho-fDrho);
    158 
    159     const Double_t dx = fR*cs;
    160     const Double_t dy = fR*sn;
    161 
    162     fSrcPos->SetXY(x-dx, y-dy);
    163 
    164     *fLog << dbg << t << " - Position: x=" << x-dx << "mm y=" << y-dy << "mm" << endl;
    165 
    166     return kTRUE;
     152// Loc0LocToCam
     153//
     154// Input :   (theta0, phi0)   direction for the position (0,0) in the camera 
     155//           ( theta,  phi)   some other direction
     156//
     157// Output :  (X, Y)      position in the camera corresponding to (theta, phi)
     158//
     159TVector2 MSrcPosCalc::CalcXYinCamera(const MVector3 &pos0, const MVector3 &pos) const
     160{
     161    const Double_t theta0 = pos0.Theta();
     162    const Double_t phi0   = pos0.Phi();
     163
     164    const Double_t theta  = pos.Theta();
     165    const Double_t phi    = pos.Phi();
     166
     167    //--------------------------------------------
     168
     169    // pos0[3] = TMath::Cos(theta0)
     170
     171    const Double_t YC0 = TMath::Cos(theta0)*TMath::Tan(theta)*TMath::Cos(phi-phi0) - TMath::Sin(theta0);
     172    const Double_t YC1 = TMath::Cos(theta0) + TMath::Sin(theta0)*TMath::Tan(theta);
     173    const Double_t YC  = YC0 / YC1;
     174
     175    //--------------------------------------------
     176
     177    const Double_t XC0 =  TMath::Cos(theta0) - YC*TMath::Sin(theta0);
     178    const Double_t XC  = -TMath::Sin(phi-phi0) * TMath::Tan(theta) * XC0;
     179
     180    //--------------------------------------------
     181    return TVector2(XC, YC);
    167182}
    168183
     
    173188Int_t MSrcPosCalc::Process()
    174189{
    175     if (fX==0 && fY==0)
    176         return kTRUE;
    177 
    178     // rotate the source position by the current rotation angle
    179     const Double_t rho = fPointPos->RotationAngle(*fObservatory, *fTime);
    180 
    181     // Define source position in the camera plain
    182     TVector2 v(fX, fY);
    183     v=v.Rotate(rho);
    184 
    185     // Convert coordinates into camera plain (mm)
    186     v *= 1./fGeom->GetConvMm2Deg();
    187 
    188     // Set current source position
    189     fSrcPos->SetXY(v);
    190 
    191     return kTRUE;
    192 }
     190  //     *fLog << dbg << "Camera center : Zd=" << fPointPos->GetZd() << " Az=" << fPointPos->GetAz() << endl;
     191  //     *fLog << dbg << "Camera center : RA=" << fPointPos->GetRa() << " Dec=" << fPointPos->GetDec() << endl;
     192  //     *fLog << dbg << "MJD: " << fTime->GetMjd() << ", time [ms]: " << fTime->GetTime() << endl;
     193
     194  MVector3 pos, pos0;  // pos: source position;  pos0: camera center
     195
     196  if (fSourcePos)
     197    {
     198      // Set Sky coordinates of source, taken from container "MSourcePos" of type MPointingPos. The sky
     199      // coordinates must be J2000, as the sky coordinates of the camera center that we get from the container
     200      // "MPointingPos" filled by the Drive.
     201
     202      pos.SetRaDec(fSourcePos->GetRaRad(), fSourcePos->GetDecRad());
     203
     204      //        *fLog << dbg << "Sky position of source: Ra=" << fSourcePos->GetRa() <<
     205      //          " Dec=" << fSourcePos->GetDec() << endl;
     206    }
     207
     208  // Convert sky coordinates of source to local coordinates. Warning! These are not the "true" local
     209  // coordinates, since this transformation ignores precession and nutation effects.
     210  pos *= MAstroSky2Local(*fTime, *fObservatory);
     211
     212
     213  // Set sky coordinates of camera center in pos0, and then convert to local. Same comment as above. These
     214  // coordinates differ from the true local coordinates of the camera center that one could get from
     215  // "MPointingPos", calculated by the Drive: the reason is that the Drive takes into account precession
     216  // and nutation corrections, while MAstroSky2Local (as of Jan 27 2005 at least) does not. Since we just
     217  // want to get the source position on the camera from the local coordinates of the center and the source,
     218  // it does not matter that the coordinates contained in pos and pos0 ignore precession and nutation...
     219  // since the shift would be the same in both cases. What would be wrong is to set in pos0 directly the
     220  // local coordinates found in MPointingPos!
     221
     222  pos0.SetRaDec(fPointPos->GetRaRad(), fPointPos->GetDecRad());
     223  pos0 *= MAstroSky2Local(*fTime, *fObservatory);
     224
     225  //     *fLog << dbg << "From MAstroSky2Local, without precession and nutation corrections:" << endl;
     226  //     *fLog << dbg << "- Camera center (2) from RA,dec and time : Zd=" << pos0.Theta()*180./TMath::Pi()
     227  //      << " Az=" << pos0.Phi()*180./TMath::Pi() << endl;
     228  //     *fLog << dbg << "- Local position of source: Zd=" << pos.Theta()*TMath::RadToDeg() << " Az=" << pos.Phi()*TMath::RadToDeg() << endl;
     229
     230
     231  // Calculate source position in camera, and convert to mm:
     232  TVector2 v = CalcXYinCamera(pos0, pos)*fGeom->GetCameraDist()*1000;
     233
     234  fSrcPosCam->SetXY(v);
     235
     236  //     v *= fGeom->GetConvMm2Deg();
     237  //     *fLog << dbg << "X=" << v.X() << " deg,  Y=" << v.Y() << " deg" << endl;
     238  //     *fLog << *fTime << endl;
     239  //     *fLog << endl;
     240
     241  return kTRUE;
     242}
Note: See TracChangeset for help on using the changeset viewer.