source: trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.cc@ 3242

/* ======================================================================== *\
!
! *
! * 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 02/2004 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MSourcePosfromStarPos
//
//  This is a task which
//  - calculates the position of the source in the camera
//    from the position of a known star in the camera
//  - and puts the source position into the container MSrcPosCam
//
//  Input :
//   ASCII file containing for each run 
//   - the run number
//   - the direction (theta, phi) the telescope is pointing to in [deg]
//   - the position  (xStar, yStar)   of a known star in the camera in [mm]
//   - the error     (dxStar, dyStar) of this position in [mm]
//
//  Output Containers :
//   MSrcPosCam
//
/////////////////////////////////////////////////////////////////////////////
#include <TList.h>
#include <TSystem.h>

#include <fstream>

#include "MSourcePosfromStarPos.h"

#include "MParList.h"
#include "MRawRunHeader.h"
#include "MGeomCam.h"
#include "MSrcPosCam.h"

#include "MLog.h"
#include "MLogManip.h"

ClassImp(MSourcePosfromStarPos);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MSourcePosfromStarPos::MSourcePosfromStarPos(
                                         const char *name, const char *title)
    : fIn(NULL)
{
    fName  = name  ? name  : "MSourcePosfromStarPos";
    fTitle = title ? title : "Calculate source position from star position";

    fFileNames = new TList;
    fFileNames->SetOwner();
}

// --------------------------------------------------------------------------
//
// Delete the filename list and the input stream if one exists.
//
MSourcePosfromStarPos::~MSourcePosfromStarPos()
{
    delete fFileNames;
    if (fIn)
        delete fIn;
}


// --------------------------------------------------------------------------
//
// Set the sky coordinates of the source and of the star
//
void MSourcePosfromStarPos::SetSourceAndStarPosition(
                            Double_t decSource, Double_t raSource,
                            Double_t decStar,   Double_t raStar)
{
  fDecSource = decSource;
  fRaSource  = raSource;

  fDecStar = decStar;
  fRaStar  = raStar;

  *fLog << all << "MSourcePosfromStarPos::SetSourceAndStarPosition; fDecSource, fRaSource, fDecStar, fRaStar were set to : "
        << fDecSource << ",  " << fRaSource << ",  "
        << fDecStar << ",  " << fRaStar << endl;
}

// --------------------------------------------------------------------------
//
//
Int_t MSourcePosfromStarPos::PreProcess(MParList *pList)
{
    MGeomCam *geom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
    if (!geom)
    {
        *fLog << err << "MSourcePosfromStarPos : MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
        return kFALSE;
    }
    fMm2Deg = geom->GetConvMm2Deg();


    fRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!fRun)
    {
        *fLog << err << "MRawRunHeader not found... aborting." << endl;
        return kFALSE;
    }


    fSrcPos = (MSrcPosCam*)pList->FindObject(AddSerialNumber("MSrcPosCam"));
    if (!fSrcPos)
    {
        *fLog << err << "MSourcePosfromStarPos : MSrcPosCam not found...  aborting" << endl;
        return kFALSE;
    }

    //---------------------------------------------------------------------
    // read all files and call ReadData() to read and store the information
    //

    fRuns = 0;
    fSize = 0;

    while(1)
    {
      if (!OpenNextFile()) break;

      if (fIn->eof())
        if (!OpenNextFile()) break;

      // FIXME! Set InputStreamID

      ReadData();
    }

    FixSize();
    //-------------------------------------------------------------

    return kTRUE;
}

//=========================================================================
//
// SourcefromStar
//
// this routine calculates the position of a source (for example Crab) in the camera
// from the position of a star (for example ZetaTauri) in the camera. The latter 
// position may have been determined by analysing the DC currents in the different
// pixels.
//
// Input  : thetaTel, phiTel          the direction the telescope is pointing to,
//                                    in local coordinates
//          f                         the distance between camera and reflector
//          decStar, raStar           the position of the star in sky coordinates
//          decSource, raSource       the position of the source in sky coordinates
//          xStar, yStar              the position of the star in the camera
//          dxStar, dyStar            error of the position of the star in the camera
//
// Output : xSource, ySource       the calculated position of the source in the camera
//          dxSource, dySource     error of the calculated position of the source in 
//                                 the camera
//
// Useful formulas can be found in TDAS 00-11 and TDAS 01-05
//

void MSourcePosfromStarPos::SourcefromStar(Double_t &f,
                    Double_t &decStar,     Double_t &raStar,
                    Double_t &decSource,   Double_t &raSource,
                    Double_t &thetaTel,    Double_t &phiTel,    
                    Double_t &xStar,       Double_t &yStar,
                    Double_t &dxStar,      Double_t &dyStar,
                    Double_t &xSource,     Double_t &ySource,
                    Double_t &dxSource,    Double_t &dySource)
{
  // the units are assumed to be radians for theta, phi, dec and ra
  //            and                   mm for f, x and y

  // calculate coordinates of star and source in system B (see TDAS 00-11, eqs. (2))
  Double_t xB  =  cos(decStar) * cos(raStar);
  Double_t yB  =  cos(decStar) * sin(raStar);
  Double_t zB  = -sin(decStar);
 
  Double_t xB0 =  cos(decSource) * cos(raSource);
  Double_t yB0 =  cos(decSource) * sin(raSource);
  Double_t zB0 = -sin(decSource);
  
  // calculate rotation angle alpha of sky image in camera 
  // (see TDAS 00-11, eqs. (18) and (20))
  // a1 = cos(Lat), a3 = -sin(Lat), where Lat is the geographical latitude of La Palma
  Double_t a1 =  0.876627;
  Double_t a3 = -0.481171;

  Double_t denom =  1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) +
                              ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) * 
                              ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) )   );
  Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom;
  Double_t sinal =    a1 * sin(phiTel) * denom;


  // calculate coordinates of star in a system with the basis vectors e1, e2, e3
  // where  e1 is in the direction (r0 x a)
  //        e2 is in the direction (e1 x r0)
  //   and  e3 is in the direction -r0;
  // r0 is the direction the telescope is pointing to
  // and a is the earth rotation axis (pointing to the celestial north pole)
  // 
  Double_t x = (-xB*yB0 + xB0*yB) / sqrt( xB0*xB0 + yB0*yB0 );
  Double_t y = ( xB*xB0*zB0 + yB*yB0*zB0 - zB*(xB0*xB0 + yB0*yB0) ) 
                                  / sqrt( xB0*xB0 + yB0*yB0 );
  Double_t z = -(xB*xB0 + yB*yB0 + zB*zB0);

  // calculate coordinates of star in camera
  Double_t xtilde = -f/z * (cosal*x - sinal*y);
  Double_t ytilde = -f/z * (sinal*x + cosal*y);

  // calculate coordinates of source in camera
  // note : in real camera signs are inverted (therefore s = -1.0)
  Double_t s = -1.0;
  xSource = s * (s*xStar - xtilde);
  ySource = s * (s*yStar - ytilde);

  *fLog << "thetaTel, phiTel, cosal, sinal = " << thetaTel << ",  "
        << phiTel << ",  " << cosal << ",  " << sinal << endl;
}


// --------------------------------------------------------------------------
//
// Get the source position and put it into MSrcPosCam
//
//
Bool_t MSourcePosfromStarPos::ReInit(MParList *pList)
{
  if (fDecStar   == 0.0 || fRaStar   == 0.0 || 
      fDecSource == 0.0 || fRaSource == 0.0)
  {
    *fLog << err << "MSourcePosfromStarPos::ReInit; sky coordinates of star and source are not defined ... aborting" 
          << endl;
    return kFALSE;
  }

  // f is the distance of the camera from the reflector center in [mm]
  Double_t f = kRad2Deg / fMm2Deg;   

  for (Int_t i=0; i<fSize; i++)
  {
    Int_t run = fRun->GetRunNumber();
    if (run == fRunNr[i])
    {
      MSourcePosfromStarPos::SourcefromStar( f,
                    fDecStar, fRaStar, fDecSource, fRaSource,
                    fThetaTel[i],   fPhiTel[i],   
                    fxStar[i],      fyStar[i],
                    fdxStar[i],     fdyStar[i],
                    fxSource,       fySource,
                    fdxSource,      fdySource);
      
      fSrcPos->SetXY(fxSource, fySource);

      *fLog << all << "MSourcePosfromStarPos::ReInit; f = " << f << endl;
      *fLog << all << "MSourcePosfromStarPos::ReInit; fRunNr, fxSource, fySource = "
            << fRunNr[i] << ",  " << fxSource << ",  " << fySource 
            << endl;
       
      fSrcPos->SetReadyToSave();       
    }
  }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//
Int_t MSourcePosfromStarPos::Process()
{

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//
Int_t MSourcePosfromStarPos::PostProcess()
{

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// read the data from the ASCII file and store them
//
void MSourcePosfromStarPos::FixSize()
{
    fSize = fRuns;

    fRunNr.Set(fSize);

    fThetaTel.Set(fSize);
    fPhiTel.Set(fSize);
    fdThetaTel.Set(fSize);
    fdPhiTel.Set(fSize);

    fxStar.Set(fSize);
    fyStar.Set(fSize);
    fdxStar.Set(fSize);
    fdyStar.Set(fSize);
}

// --------------------------------------------------------------------------
//
// read the data from the ASCII file and store them
//
void MSourcePosfromStarPos::ReadData()
{
  Float_t val;

  if (fRuns >= fSize)
  {
    fSize += 100;

    fRunNr.Set(fSize);

    fThetaTel.Set(fSize);
    fPhiTel.Set(fSize);
    fdThetaTel.Set(fSize);
    fdPhiTel.Set(fSize);

    fxStar.Set(fSize);
    fyStar.Set(fSize);
    fdxStar.Set(fSize);
    fdyStar.Set(fSize);
  }

  fRuns += 1;

  *fIn >> val;
  fRunNr.AddAt(val, fRuns-1);

  *fIn >> val;
  fThetaTel.AddAt(val, fRuns-1);
  *fIn >> val;
  fPhiTel.AddAt(val, fRuns-1);

  *fIn >> val;
  fdThetaTel.AddAt(val, fRuns-1);
  *fIn >> val;
  fdPhiTel.AddAt(val, fRuns-1);

  *fIn >> val;
  fxStar.AddAt(val, fRuns-1);
  *fIn >> val;
  fyStar.AddAt(val, fRuns-1);

  *fIn >> val;
  fdxStar.AddAt(val, fRuns-1);
  *fIn >> val;
  fdyStar.AddAt(val, fRuns-1);
}

// --------------------------------------------------------------------------
//
// Add this file as the last entry in the chain
//
Int_t MSourcePosfromStarPos::AddFile(const char *txt, Int_t)
{
    TNamed *name = new TNamed(txt, "");
    fFileNames->AddLast(name);
    return 1;
}

// --------------------------------------------------------------------------
//
// This opens the next file in the list and deletes its name from the list.
//
Bool_t MSourcePosfromStarPos::OpenNextFile()
{
    //
    // open the input stream and check if it is really open (file exists?)
    //
    if (fIn)
        delete fIn;
    fIn = NULL;

    //
    // Check for the existence of a next file to read
    //
    TNamed *file = (TNamed*)fFileNames->First();
    if (!file)
        return kFALSE;

    //
    // open the file which is the first one in the chain
    //
    const char *name = file->GetName();

    const char *expname = gSystem->ExpandPathName(name);
    fIn = new ifstream(expname);
    delete [] expname;

    const Bool_t noexist = !(*fIn);

    if (noexist)
        *fLog << dbginf << "Cannot open file '" << name << "'" << endl;
    else
        *fLog << "Open file: '" << name << "'" << endl;

    //
    // Remove this file from the list of pending files
    //
    fFileNames->Remove(file);

    return !noexist;
}

// --------------------------------------------------------------------------