Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h
===================================================================
--- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h	(revision 3204)
+++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h	(revision 3209)
@@ -80,4 +80,5 @@
 #pragma link C++ class MArrivalTimeCalc2+;
 
+#pragma link C++ class MSourcePosfromStarPos+;
 #pragma link C++ class MMcCalibrationUpdate+;
 
@@ -86,2 +87,4 @@
 
 #endif
+
+
Index: trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.cc
===================================================================
--- trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.cc	(revision 3209)
+++ trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.cc	(revision 3209)
@@ -0,0 +1,448 @@
+/* ======================================================================== *\
+!
+! *
+! * 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;
+}
+
+// --------------------------------------------------------------------------
+
+
+
+
+
+
+
+
+
+
Index: trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.h
===================================================================
--- trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.h	(revision 3209)
+++ trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.h	(revision 3209)
@@ -0,0 +1,99 @@
+#ifndef MARS_MSourcePosfromStarPos
+#define MARS_MSourcePosfromStarPos
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// MSourcePosfromStarPos                                                   //
+//                                                                         //
+// Task to calculate the position of the source in the camera from
+//      the position of a known star in the camera
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+
+#ifndef MARS_MTask
+#include "MTask.h"
+#endif
+
+#ifndef ROOT_TArrayI
+#include <TArrayI.h>
+#endif
+
+#ifndef ROOT_TArrayD
+#include <TArrayD.h>
+#endif
+
+class TList;
+class MRawRunHeader;
+class MGeomCam;
+class MSrcPosCam;
+
+
+class MSourcePosfromStarPos : public MTask
+{
+ private:
+    const MRawRunHeader *fRun;      //!
+    const MGeomCam      *fGeomCam;  //! Camera Geometry used to calculate Hillas
+    MSrcPosCam    *fSrcPos;   //!
+
+    ifstream    *fIn;             // input file
+    TList       *fFileNames;      // array which contains the \0-terminated file names
+
+    Float_t fMm2Deg;
+
+    Int_t   fRuns;                // current number of entries in TArray
+    Int_t   fSize;                // final   number of entries in TArray
+
+    Double_t fDecStar;
+    Double_t fRaStar;
+
+    Double_t fDecSource;
+    Double_t fRaSource;
+
+    Double_t fxSource;
+    Double_t fySource;
+    Double_t fdxSource;
+    Double_t fdySource;
+
+    TArrayI fRunNr;
+    TArrayD fThetaTel;
+    TArrayD fPhiTel;
+    TArrayD fdThetaTel;
+    TArrayD fdPhiTel;
+    TArrayD fxStar;
+    TArrayD fyStar;
+    TArrayD fdxStar;
+    TArrayD fdyStar;
+
+    Int_t  AddFile(const char *fname, Int_t dummy=-1);
+    Bool_t OpenNextFile();
+    void   ReadData();
+    void   FixSize();
+
+    Int_t  PreProcess(MParList *pList);
+    Bool_t ReInit(MParList *pList);
+    Int_t  Process();
+    Int_t  PostProcess();
+
+public:
+    MSourcePosfromStarPos(const char *name=NULL, const char *title=NULL);
+    ~MSourcePosfromStarPos();
+
+    void SetSourceAndStarPosition(Double_t decSource, Double_t raSource,
+				  Double_t decStar,   Double_t raStar);
+
+    void SourcefromStar(Double_t &, Double_t &, Double_t &,
+      Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &,
+      Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t & );
+
+    ClassDef(MSourcePosfromStarPos, 0) // Task to calculate the source position from a star position
+};
+
+#endif
+
+
+
+
+
+
+
+
Index: trunk/MagicSoft/Mars/manalysis/Makefile
===================================================================
--- trunk/MagicSoft/Mars/manalysis/Makefile	(revision 3204)
+++ trunk/MagicSoft/Mars/manalysis/Makefile	(revision 3209)
@@ -88,4 +88,5 @@
            MArrivalTimeCalc.cc \
            MArrivalTimeCalc2.cc \
+           MSourcePosfromStarPos.cc \
            MMcCalibrationUpdate.cc
 
