Changeset 3274

Timestamp:

02/24/04 15:49:17 (21 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• macros/ONOFFAnalysis.C (modified) (16 diffs)
• manalysis/MSourcePosfromStarPos.cc (modified) (19 diffs)
• manalysis/MSourcePosfromStarPos.h (modified) (5 diffs)
• mfilter/MFSelBasic.cc (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -5 +5 @@
                                                  -*-*- END OF LINE -*-*-
  
+
+
+
+ 2004/02/24: Wolfgang Wittek
+  * manalysis/MSourcPosfromStarPos.[h,cc]
+    - change member function SetSourceAndStarPosition() to expect sky
+      coordinates in the standard units
+    - generalize to more than 1 star
+    - the class is not yet fully tested
+
+  * mfilter/MFSelBasic.[h,cc]
+    - change default values of cuts
+
 
  2004/02/24: Markus Gaug

trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C

@@ -172 +172 @@
       //-----------------------------------------------
       //const char *offfile = "~magican/ct1test/wittek/offdata.preproc"; 
-      const char *offfile = "20040126_OffCrab_"; 
+      //const char *offfile = "20040126_OffCrab_"; 
+      //const char *offfile  = "*.OFF";
+      const char *offfile  = "12*.OFF";
 
       //const char *onfile  = "~magican/ct1test/wittek/mkn421_on.preproc"; 
       //      const char *onfile  = "~magican/ct1test/wittek/mkn421_00-01"; 
-      const char *onfile  = "20040126_Crab_";
+      //const char *onfile  = "20040126_Crab_";
+      //const char *onfile  = "MCerPhot_output";
+      //const char *onfile  = "*.ON";
+      const char *onfile  = "12*.ON";
 
      const char *mcfile  = "/data/MAGIC/mc_eth/magLQE_3/gh/0/0/G_M0_00_0_550*.root";
@@ -184 +189 @@
 
       // path for input for Mars
-      TString inPath = "/.magic/magicserv01/scratch/";
+      //TString inPath = "/.magic/magicserv01/scratch/";
+     TString inPath = "/mnt/data17a/hbartko/";
+     //TString inPath = "~wittek/datacrab_feb04/";
 
       // path for output from Mars
-      TString outPath = "~wittek/datacrab/";
+      TString outPath = "~wittek/datacrab_feb04/";
 
       //-----------------------------------------------
@@ -235 +242 @@
     Bool_t JobB_SC_UP  = kFALSE;
     Bool_t CMatrix     = kFALSE;  // create training and test matrices 
-    Bool_t RMatrix     = kTRUE;  // read training and test matrices from file
-    Bool_t WOptimize   = kTRUE;  // do optimization using the training sample
+    Bool_t RMatrix     = kFALSE;  // read training and test matrices from file
+    Bool_t WOptimize   = kFALSE;  // do optimization using the training sample
                                   // and write supercuts parameter values 
                                   // onto the file parSCfile
     Bool_t RTest       = kFALSE;  // test the supercuts using the test matrix
-    Bool_t WSC         = kFALSE;  // update input root file ?
+    Bool_t WSC         = kTRUE;  // update input root file ?
 
 
@@ -316 +323 @@
     TString fileON  = inPath;
     fileON += onfile;
-    fileON += "CalibratedEvts";
     fileON += ".root";
 
     TString fileOFF = inPath;
     fileOFF += offfile;
-    fileOFF += "CalibratedEvts";
     fileOFF += ".root";
 
     TString fileMC  = mcfile;
     fileMC += mcfile;
-    fileMC += "CalibratedEvts";
     fileMC += ".root";
 
@@ -339 +343 @@
     //--------------------------------------------------
     // type of data to be padded 
-    //TString typeInput = "ON";
-    TString typeInput = "OFF";
+    TString typeInput = "ON";
+    //TString typeInput = "OFF";
     //TString typeInput = "MC";
     gLog << "typeInput = " << typeInput << endl;
@@ -366 +370 @@
     outNameImage += "Hillas";
     outNameImage += typeInput;
-    outNameImage += "1.root";
+    outNameImage += "1a.root";
     gLog << "padded data to be written onto : " << outNameImage << endl;
 
@@ -680 +684 @@
     MGeomApply        apply;
 
-    MPedestalWorkaround waround;
+    //MPedestalWorkaround waround;
 
     // a way to find out whether one is dealing with MC :
@@ -693 +697 @@
     //                                             "MCT1PointingCorrCalc");
     //}
+    MSourcePosfromStarPos sourcefromstar;
+    sourcefromstar.SetSourceAndStarPosition("Crab", 22, 0, 52, 5, 34, 32,
+                                        "Zeta-Tau", 21, 8, 33, 5, 37, 38.7);
+    sourcefromstar.AddFile("~wittek/datacrab_feb04/positions.4.txt", 0);
 
     MBlindPixelCalc blindbeforepad;
@@ -706 +714 @@
     contbasic.SetName("SelBasic");
 
-    //MFillH fillblind("BlindPixels[MHBlindPixels]", "MBlindPixels");
-    //fillblind.SetName("HBlind");
+    MFillH fillblind("BlindPixels[MHBlindPixels]", "MBlindPixels");
+    fillblind.SetName("HBlind");
 
     MSigmabarCalc sigbarcalc;
@@ -750 +758 @@
     selstandard.SetHillasName(fHilName);
     selstandard.SetImgParName(fImgParName);
-    selstandard.SetCuts(100, 4, 20, 0.0, 1.0, 0.0, 0.0);
+    selstandard.SetCuts(200, 6, 600, 0.4, 1.1, 0.0, 0.0);
     MContinue contstandard(&selstandard);
     contstandard.SetName("SelStandard");
@@ -786 +794 @@
     //tliston.AddToList(&f2);
     tliston.AddToList(&apply);
-    tliston.AddToList(&waround);
+    //tliston.AddToList(&waround);
+    tliston.AddToList(&sourcefromstar);
 
     //tliston.AddToList(&blindbeforepad);
@@ -793 +802 @@
     //  tliston.AddToList(&pointcorr);
 
-    //tliston.AddToList(&blind);
-    tliston.AddToList(&contbasic);
-
-    //tliston.AddToList(&fillblind);
+    tliston.AddToList(&blind);
+    //tliston.AddToList(&contbasic);
+
+    tliston.AddToList(&fillblind);
     tliston.AddToList(&sigbarcalc);
     tliston.AddToList(&fillsigtheta);
@@ -1783 +1792 @@
     TString parSCinit = outPath;
     //parSCinit += "parSC_060204";
+    //parSCinit = "parSC_240204a";
     parSCinit = "";
 
@@ -1802 +1812 @@
 
     TString parSCfile = outPath;
-    parSCfile += "parSC_130204a";
+    parSCfile += "parSC_240204a";
 
     gLog << "parSCfile = " << parSCfile << endl;
@@ -1843 +1853 @@
     filenameTrain += typeInput;
     filenameTrain += "1.root";
-    Int_t howManyTrain = 7000;
+    Int_t howManyTrain = 20000;
     gLog << "filenameTrain = " << filenameTrain << ",   howManyTrain = "
          << howManyTrain  << endl; 
@@ -1853 +1863 @@
     filenameTest += typeInput;
     filenameTest += "1.root";
-    Int_t howManyTest = 7000;
+    Int_t howManyTest = 20000;
 
     gLog << "filenameTest = " << filenameTest << ",   howManyTest = "

trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.cc

@@ -44 +44 @@
 #include <TList.h>
 #include <TSystem.h>
+#include <TMatrixD.h>
 
 #include <fstream>
@@ -53 +54 @@
 #include "MGeomCam.h"
 #include "MSrcPosCam.h"
+#include "MMcEvt.hxx"
 
 #include "MLog.h"
@@ -67 +69 @@
 MSourcePosfromStarPos::MSourcePosfromStarPos(
                                          const char *name, const char *title)
-    : fIn(NULL)
+  : fIn(NULL)
 {
     fName  = name  ? name  : "MSourcePosfromStarPos";
@@ -74 +76 @@
     fFileNames = new TList;
     fFileNames->SetOwner();
+
+    fRuns  = 0;
+    fSize  = 100;
+    fStars = 1;
+
+    fRunNr.Set(fSize);
+
+    fThetaTel.Set(fSize);
+    fPhiTel.Set(fSize);
+    fdThetaTel.Set(fSize);
+    fdPhiTel.Set(fSize);
+
+    fDecStar.Set(fStars);
+    fRaStar.Set(fStars);
+    fxStar.ResizeTo(fStars,fSize);
+    fyStar.ResizeTo(fStars,fSize);
+    fdxStar.ResizeTo(fStars,fSize);
+    fdyStar.ResizeTo(fStars,fSize);
 }
 
@@ -92 +112 @@
 // Set the sky coordinates of the source and of the star
 //
+// Input :
+// declination in units of     (Deg,  Min, Sec)
+// right ascension in units of (Hour, Min, Sec)
+//
+
 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 : "
+         TString  nameSource,
+         Double_t decSourceDeg, Double_t decSourceMin, Double_t decSourceSec, 
+         Double_t raSourceHour, Double_t raSourceMin,  Double_t raSourceSec,
+         TString nameStar,
+         Double_t decStarDeg,   Double_t decStarMin,   Double_t decStarSec, 
+         Double_t raStarHour,   Double_t raStarMin,    Double_t raStarSec  )
+{
+  *fLog << "MSourcePosfromStarPos::SetSourceAndStarPosition :" << endl;
+  *fLog << "       Source : "  << nameSource << "   " << decSourceDeg << ":" 
+        << decSourceMin << ":" << decSourceSec << endl;
+  *fLog << "       Star   : "  << nameStar   << "   " << decStarDeg << ":" 
+        << decStarMin << ":"   << decStarSec << endl;
+
+  // convert into radians
+  fDecSource = (decSourceDeg + decSourceMin/60.0 + decSourceSec/3600.0)
+               / kRad2Deg;
+  fRaSource  = (raSourceHour + raSourceMin/60.0  + raSourceSec/3600.0)
+               * 360.0 / (24.0 * kRad2Deg);
+
+  fDecStar.Set(fStars);
+  fDecStar[fStars-1] = (decStarDeg + decStarMin/60.0 + decStarSec/3600.0)
+                       / kRad2Deg;
+  fRaStar[fStars-1]  = (raStarHour + raStarMin/60.0  + raStarSec/3600.0)
+                       * 360.0 / (24.0 * kRad2Deg);
+
+  *fLog << all << "MSourcePosfromStarPos::SetSourceAndStarPosition; fDecSource, fRaSource, fDecStar, fRaStar were set to : [radians]  "
         << fDecSource << ",  " << fRaSource << ",  "
-        << fDecStar << ",  " << fRaStar << endl;
+        << fDecStar[fStars-1] << ",  " << fRaStar[fStars-1] << endl;
+}
+
+// --------------------------------------------------------------------------
+//
+// Set the sky coordinates of another star
+//
+// Input :
+// declination in units of     (Deg,  Min, Sec)
+// right ascension in units of (Hour, Min, Sec)
+//
+
+void MSourcePosfromStarPos::AddStar(
+         TString nameStar,
+         Double_t decStarDeg,   Double_t decStarMin,   Double_t decStarSec, 
+         Double_t raStarHour,   Double_t raStarMin,    Double_t raStarSec  )
+{
+  *fLog << "MSourcePosfromStarPos::AddStar :" << endl;
+  *fLog << "       Star   : "  << nameStar   << "   " << decStarDeg << ":" 
+        << decStarMin << ":"   << decStarSec << endl;
+
+  // convert into radians
+  Int_t fStars = fDecStar.GetSize() + 1;
+  fDecStar.Set(fStars);
+  fDecStar[fStars-1] = (decStarDeg + decStarMin/60.0 + decStarSec/3600.0)
+                       / kRad2Deg;
+  fRaStar[fStars-1]  = (raStarHour + raStarMin/60.0  + raStarSec/3600.0)
+                       * 360.0 / (24.0 * kRad2Deg);
+
+  *fLog << all << "MSourcePosfromStarPos::AddStar; fDecStar, fRaStar were set to : [radians]  "
+        << fDecStar[fStars-1] << ",  " << fRaStar[fStars-1] << endl;
 }
 
@@ -119 +190 @@
     }
     fMm2Deg = geom->GetConvMm2Deg();
+    // fDistCameraReflector is the distance of the camera from the reflector center in [mm]
+    fDistCameraReflector = kRad2Deg / fMm2Deg;   
+        *fLog << all << "MSourcePosfromStarPos::PreProcess; fMm2Deg, fDistCameraReflector = " 
+              << fMm2Deg << ",  " << fDistCameraReflector << endl;
 
 
@@ -124 +199 @@
     if (!fRun)
     {
-        *fLog << err << "MRawRunHeader not found... aborting." << endl;
+        *fLog << err << "MSourcePosfromStarPos::PreProcess; MRawRunHeader not found... aborting." << endl;
         return kFALSE;
     }
+
+
+   fMcEvt = (MMcEvt*)pList->FindCreateObj("MMcEvt");
+   if (!fMcEvt)
+   {
+       *fLog << err << "MSourcePosfromStarPos::PreProcess; MMcEvt not found... aborting." << endl;
+       return kFALSE;
+   }
 
 
@@ -132 +215 @@
     if (!fSrcPos)
     {
-        *fLog << err << "MSourcePosfromStarPos : MSrcPosCam not found...  aborting" << endl;
+        *fLog << err << "MSourcePosfromStarPos::PreProcess; MSrcPosCam not found...  aborting" << endl;
         return kFALSE;
     }
@@ -140 +223 @@
     //
 
-    fRuns = 0;
-    fSize = 0;
-
     while(1)
     {
-      if (!OpenNextFile()) break;
-
-      if (fIn->eof())
-        if (!OpenNextFile()) break;
-
-      // FIXME! Set InputStreamID
-
-      ReadData();
+      if (!OpenNextFile()) 
+      {
+        *fLog << "there is no more file to open" << endl;
+        break;
+      }
+
+      *fLog << "read data" << endl;
+      while (1)
+      {
+        if (fIn->eof())
+        {
+          *fLog << "eof encountered; open next file" << endl;
+
+          if (!OpenNextFile()) break;
+        }
+
+        // FIXME! Set InputStreamID
+      
+        ReadData();
+      }
     }
 
+    *fLog << "all data were read" << endl;
     FixSize();
     //-------------------------------------------------------------
@@ -186 +279 @@
 
 void MSourcePosfromStarPos::SourcefromStar(Double_t &f,
-                    Double_t &decStar,     Double_t &raStar,
+                    TArrayD  &decStar,     TArrayD  &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,
+                    TArrayD  &xStar,       TArrayD  &yStar,
+                    TArrayD  &dxStar,      TArrayD  &dyStar,
                     Double_t &xSource,     Double_t &ySource,
                     Double_t &dxSource,    Double_t &dySource)
 {
+  *fLog << "MSourcePosfromStarPos::SourcefromStar :  printout in degrees" << endl;
+  *fLog << "       decStar, raStar = " << decStar[0]*kRad2Deg << ",  " 
+        << raStar[0]*kRad2Deg << endl;
+  *fLog << "       decSource, raSource = " << decSource*kRad2Deg << ",  " 
+        << raSource*kRad2Deg << endl;
+  *fLog << "       thetaTel, phiTel = " << thetaTel*kRad2Deg << ",  " 
+        << phiTel*kRad2Deg << endl;
+  *fLog << "       xStar, yStar = " << xStar[0]*fMm2Deg << ",  " 
+        << yStar[0]*fMm2Deg << endl;
+
+  *fLog << "MSourcePosfromStarPos::SourcefromStar :  printout in radians and mm" << endl;
+  *fLog << "       decStar, raStar = " << decStar[0] << ",  " 
+        << raStar[0] << endl;
+  *fLog << "       decSource, raSource = " << decSource << ",  " 
+        << raSource << endl;
+  *fLog << "       thetaTel, phiTel = " << thetaTel << ",  " 
+        << phiTel << endl;
+  *fLog << "       xStar, yStar = " << xStar[0] << ",  " 
+        << yStar[0] << endl;
+
+
   // 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))
@@ -218 +324 @@
   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;
-}
-
+
+
+  // calculate coordinates of source in system B (see TDAS 00-11, eqs. (2))
+  Double_t xB0 =  cos(decSource) * cos(raSource);
+  Double_t yB0 =  cos(decSource) * sin(raSource);
+  Double_t zB0 = -sin(decSource);
+
+  *fLog << "xB0, yB0, zB0 = " << xB0 << ",  " << yB0 << ",  "
+        << zB0 << endl;
+
+  //-----------------------------------------------------
+  // loop over stars
+  Double_t sumx  = 0.0;
+  Double_t sumy  = 0.0;
+  Double_t sumwx = 0.0;
+  Double_t sumwy = 0.0;
+
+  for (Int_t i=0; i<decStar.GetSize(); i++)
+  {
+    // calculate weights
+    Double_t weightx = 1.0 / (dxStar[i]*dxStar[i]);
+    Double_t weighty = 1.0 / (dyStar[i]*dyStar[i]);
+    sumwx += weightx;
+    sumwy += weighty;
+
+  *fLog << "weightx, weighty = " << weightx << ",  " << weighty << endl;
+
+    // calculate coordinates of star in system B (see TDAS 00-11, eqs. (2))
+    Double_t xB  =  cos(decStar[i]) * cos(raStar[i]);
+    Double_t yB  =  cos(decStar[i]) * sin(raStar[i]);
+    Double_t zB  = -sin(decStar[i]);
+
+
+  *fLog << "xB, yB, zB = " << xB << ",  " << yB << ",  "
+        << zB << endl;
+
+ 
+    // 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 to the source
+    // 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);
+
+  *fLog << "x, y, z = " << x << ",  " << y << ",  "
+        << z << endl;
+
+
+    // calculate coordinates of star in camera
+    Double_t xtilde = -f/z * (cosal*x - sinal*y);
+    Double_t ytilde = -f/z * (sinal*x + cosal*y);
+
+  *fLog << "xtilde, ytile = " << xtilde << ",  " << ytilde << endl;
+
+
+    // calculate coordinates of source in camera
+    // note : in real camera signs are inverted (therefore s = -1.0)
+    Double_t s = -1.0;
+    sumx += s * (s*xStar[i] - xtilde) * weightx;
+    sumy += s * (s*yStar[i] - ytilde) * weighty;
+  }
+  //-----------------------------------------------------
+
+  xSource  = sumx / sumwx;
+  ySource  = sumy / sumwy;
+  dxSource = 1.0 / sqrt(sumwx);
+  dySource = 1.0 / sqrt(sumwy);
+    
+  *fLog << all << "MSourcePosfromStarPos::SourcefromStar; xSource, ySource = "
+        << xSource << " +- " << dxSource << ",   " 
+        << ySource << " +- " << dySource << endl; 
+}
 
 // --------------------------------------------------------------------------
@@ -253 +411 @@
 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;   
-
+  //if (fDecSource == 0.0  ||  fRaSource == 0.0  ||  fStars == 0) 
+  //{
+  //  *fLog << err << "MSourcePosfromStarPos::ReInit; sky coordinates of star and source are not defined ... aborting" 
+  //        << endl;
+    //return kFALSE;
+  //}
+
+  Int_t run = fRun->GetRunNumber();
+
+  *fLog << "MSourcePosfromStarPos::ReInit; run = " << run << endl;
+
+  //-------------------------------------------------------------------
+  // search this run in the list 
   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);
+      //-----------------------------------------
+      // put the zenith angle into MMcEvt
+
+      Double_t thetarad = fThetaTel[i];
+      Double_t phirad   = fPhiTel[i];
+      fMcEvt->SetTelescopeTheta(thetarad);
+      fMcEvt->SetTelescopePhi(phirad);
+      fMcEvt->SetReadyToSave();
+
+      *fLog << "theta, phi = " << thetarad*kRad2Deg << ",  "
+            << phirad*kRad2Deg << endl;
+       
+      //-----------------------------------------
+      // Get source position and put it into MSrcPosCam
+
+      /*
+      if (fStars > 0)
+      {
+        TArrayD xStar(fxStar.GetNrows());
+        TArrayD dxStar(fdxStar.GetNrows());
+        TArrayD yStar(fyStar.GetNrows());
+        TArrayD dyStar(fdyStar.GetNrows());
+        for (Int_t j=0; j<fxStar.GetNrows(); j++)
+        {
+          xStar[j]  = fxStar(j, i);
+          dxStar[j] = fdxStar(j, i);
+          yStar[j]  = fyStar(j, i);
+          dyStar[j] = fdyStar(j, i);
+        }
+
+        MSourcePosfromStarPos::SourcefromStar( fDistCameraReflector,
+                      fDecStar, fRaStar, fDecSource, fRaSource,
+                      fThetaTel[i],   fPhiTel[i],   
+                      xStar,          yStar,
+                      dxStar,         dyStar,
+                      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->SetXY(fxSource, fySource);
+
+        *fLog << all << "MSourcePosfromStarPos::ReInit; fRunNr, fxSource, fySource = "
+              << fRunNr[i] << ",  " << fxSource << ",  " << fySource 
+              << endl;
        
-      fSrcPos->SetReadyToSave();       
+        fSrcPos->SetReadyToSave();       
+      }
+      */
+
+      return kTRUE;
     }
   }
+  //-------------------------------------------------------------------  
+    *fLog << warn << "MSourcePosfromStarPos::ReInit;  no information for run number = "
+          << run << endl;
+
+    Double_t thetadeg = 90.0;
+    Double_t thetarad = thetadeg / kRad2Deg;
+    fMcEvt->SetTelescopeTheta(thetarad);
+
+    Double_t phideg = 0.0;
+    Double_t phirad = phideg / kRad2Deg;
+    fMcEvt->SetTelescopePhi(phirad);
+    fMcEvt->SetReadyToSave();
 
     return kTRUE;
@@ -296 +500 @@
 Int_t MSourcePosfromStarPos::Process()
 {
+  Int_t run = fRun->GetRunNumber();
+
+  //*fLog << "MSourcePosfromStarPos::Process; run = " << run << endl;
+    
 
     return kTRUE;
@@ -315 +523 @@
 void MSourcePosfromStarPos::FixSize()
 {
+  *fLog << "MSourcePosfromStarPos::FixSize; fix size of arrays : fRuns = "
+        << fRuns << endl;
+
     fSize = fRuns;
 
@@ -324 +535 @@
     fdPhiTel.Set(fSize);
 
-    fxStar.Set(fSize);
-    fyStar.Set(fSize);
-    fdxStar.Set(fSize);
-    fdyStar.Set(fSize);
+    fStars = fxStar.GetNrows();
+    fxStar.ResizeTo(fStars, fSize);
+    fyStar.ResizeTo(fStars, fSize);
+    fdxStar.ResizeTo(fStars, fSize);
+    fdyStar.ResizeTo(fStars, fSize);
 }
 
@@ -337 +549 @@
 {
   Float_t val;
-
-  if (fRuns >= fSize)
+  Int_t   ival;
+
+  // extend size of arrays if necessary
+  if ( fRuns >= (fSize-1) )
   {
     fSize += 100;
@@ -349 +563 @@
     fdPhiTel.Set(fSize);
 
-    fxStar.Set(fSize);
-    fyStar.Set(fSize);
-    fdxStar.Set(fSize);
-    fdyStar.Set(fSize);
+    fStars = fxStar.GetNrows();
+    fxStar.ResizeTo(fStars, fSize);
+    fyStar.ResizeTo(fStars, fSize);
+    fdxStar.ResizeTo(fStars, fSize);
+    fdyStar.ResizeTo(fStars, fSize);
   }
 
+  //-------------------
+  // read header line
+  //*fIn >> val;
+
+  //*fLog << "val =" << val << endl;
+
+  //*fIn >> val;
+  //*fIn >> val;
+  //*fIn >> val;
+  //*fIn >> val;
+  //*fIn >> val;
+  //*fIn >> val;
+  //*fIn >> val;
+  //*fIn >> val;
+  //*fIn >> val;
+  //*fIn >> val;
+  //-------------------
+
+
   fRuns += 1;
 
+  *fIn >> ival;
+  
+  *fLog << fRuns <<"th run : " << ival << endl;  
+
+  fRunNr.AddAt(ival, fRuns-1);
+
+  *fLog << "check : fRuns, fRunNr[fRuns-1], fRunNr[fRuns] = " << fRuns << ",  "
+        << fRunNr[fRuns-1] << ",  " << fRunNr[fRuns] << endl;
+
+  // read mjdS, hmsS, mjdE, hmsE
+  *fIn >> val; 
+  *fIn >> val; 
+  *fIn >> val; 
+  *fIn >> val; 
+
+  *fIn >> val; 
+  *fIn >> val; 
+  *fIn >> val; 
+  *fIn >> val; 
+
+
   *fIn >> val;
-  fRunNr.AddAt(val, fRuns-1);
+  fThetaTel.AddAt(val/kRad2Deg, fRuns-1);
+  *fLog << "val, fThetaTel[fRuns-1] = " << val << ",  "
+        << fThetaTel[fRuns-1] << endl;
+
 
   *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);
+  fPhiTel.AddAt(val/kRad2Deg, fRuns-1);
+  *fLog << "val, fPhiTel[fRuns-1] = " << val << ",  "
+        << fPhiTel[fRuns-1] << endl;
+
+
+  //*fIn >> val;
+  //fdThetaTel.AddAt(val/kRad2Deg, fRuns-1);
+  //*fIn >> val;
+  //fdPhiTel.AddAt(val/kRad2Deg, fRuns-1);
+
+  // input is in [deg], convert to [mm]
+  fStars = fxStar.GetNrows();
+  for (Int_t i=0; i<fStars; i++)
+  {
+    *fIn >> val;
+    fxStar(i, fRuns-1) = val / fMm2Deg;;
+    *fLog << "val, fxStar(i, fRuns-1) = " << val << ",  "
+          << fxStar(i, fRuns-1) << endl;
+
+    *fIn >> val;
+    fyStar(i, fRuns-1) = val / fMm2Deg;
+    *fLog << "val, fyStar(i, fRuns-1) = " << val << ",  "
+          << fyStar(i, fRuns-1) << endl;
+
+
+    *fIn >> val;
+    //*fLog << "y=dxStar = " << val << endl;
+
+    fdxStar(i, fRuns-1) = val / fMm2Deg;
+    *fIn >> val;
+    //*fLog << "y=dyStar = " << val << endl;
+
+    fdyStar(i, fRuns-1) = val / fMm2Deg;
+  }
+
 }
 
@@ -389 +664 @@
     TNamed *name = new TNamed(txt, "");
     fFileNames->AddLast(name);
+
+    *fLog << "MSourcePosfromStarPos::AddFile; add file '" << txt << "'"
+          << endl;
+
     return 1;
 }
@@ -447 +726 @@
 
 
+

trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.h

@@ -23 +23 @@
 #endif
 
+#ifndef ROOT_TMatrixD
+#include <TMatrixD.h>
+#endif
+
 class TList;
 class MRawRunHeader;
+class MMcEvt;
 class MGeomCam;
 class MSrcPosCam;
@@ -34 +39 @@
     const MRawRunHeader *fRun;      //!
     const MGeomCam      *fGeomCam;  //! Camera Geometry used to calculate Hillas
+    MMcEvt              *fMcEvt;
     MSrcPosCam    *fSrcPos;   //!
 
@@ -39 +45 @@
     TList       *fFileNames;      // array which contains the \0-terminated file names
 
-    Float_t fMm2Deg;
+    Double_t fMm2Deg;
+    Double_t fDistCameraReflector;
 
-    Int_t   fRuns;                // current number of entries in TArray
-    Int_t   fSize;                // final   number of entries in TArray
+    Int_t   fRuns;                // current number of runs
+    Int_t   fSize;                // final   number of runs
+    Int_t   fStars;               // number of stars
 
-    Double_t fDecStar;
-    Double_t fRaStar;
 
     Double_t fDecSource;
     Double_t fRaSource;
-
     Double_t fxSource;
     Double_t fySource;
@@ -60 +65 @@
     TArrayD fdThetaTel;
     TArrayD fdPhiTel;
-    TArrayD fxStar;
-    TArrayD fyStar;
-    TArrayD fdxStar;
-    TArrayD fdyStar;
 
-    Int_t  AddFile(const char *fname, Int_t dummy=-1);
+    TArrayD  fDecStar;
+    TArrayD  fRaStar;
+    TMatrixD fxStar;
+    TMatrixD fyStar;
+    TMatrixD fdxStar;
+    TMatrixD fdyStar;
+
+
     Bool_t OpenNextFile();
     void   ReadData();
@@ -79 +87 @@
     ~MSourcePosfromStarPos();
 
-    void SetSourceAndStarPosition(Double_t decSource, Double_t raSource,
-                                  Double_t decStar,   Double_t raStar);
+    void SetSourceAndStarPosition(
+         TString  nameSource,
+         Double_t decSourceDeg, Double_t decSourceMin, Double_t decSourceSec, 
+         Double_t raSourceHour, Double_t raSourceMin,  Double_t raSourceSec,
+         TString  nameStar,
+         Double_t decStarDeg,   Double_t decStarMin,   Double_t decStarSec, 
+         Double_t raStarHour,   Double_t raStarMin,    Double_t raStarSec  );
 
-    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 & );
+    Int_t  AddFile(const char *fname, Int_t dummy=-1);
+
+    void AddStar(
+         TString  nameStar,
+         Double_t decStarDeg,   Double_t decStarMin,   Double_t decStarSec, 
+         Double_t raStarHour,   Double_t raStarMin,    Double_t raStarSec  );
+
+    void SourcefromStar(Double_t &, 
+      TArrayD  &, TArrayD  &, Double_t &, Double_t &, Double_t &, Double_t &, 
+      TArrayD  &, TArrayD  &, TArrayD  &, TArrayD  &, Double_t &, Double_t &, 
+      Double_t &, Double_t & );
 
     ClassDef(MSourcePosfromStarPos, 0) // Task to calculate the source position from a star position

trunk/MagicSoft/Mars/mfilter/MFSelBasic.cc

@@ -71 +71 @@
 
     // default values of cuts
-    SetCuts(13.0, 0.0, 60.0);
+    SetCuts(20.0, 0.0, 60.0);
 }
 
@@ -151 +151 @@
 
     // remove bad runs for MC gammas
-    if (fMcEvt->GetEnergy() == 0.0  &&  fMcEvt->GetImpact() == 0.0)
-    {
-      if (fRawRun->GetRunNumber() == 601  ||
-          fRawRun->GetRunNumber() == 613  ||
-          fRawRun->GetRunNumber() == 614    )
-        return Set(1);
-    }
+    //if (fMcEvt->GetEnergy() == 0.0  &&  fMcEvt->GetImpact() == 0.0)
+    //{
+    //  if (fRawRun->GetRunNumber() == 601  ||
+    //      fRawRun->GetRunNumber() == 613  ||
+    //      fRawRun->GetRunNumber() == 614    )
+    //  return Set(1);
+    //}
 
     if (theta<fThetaMin)