Changeset 4741 for trunk/MagicSoft/Mars/mtemp

Timestamp:

08/26/04 14:13:36 (20 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mtemp

Files:

• MTelAxisFromStars.cc (modified) (10 diffs)
• MTelAxisFromStars.h (modified) (4 diffs)
• findTelAxisFromStars.C (modified) (7 diffs)

trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.cc

@@ -48 +48 @@
 
 #include "MParList.h"
-#include "MSrcPosCam.h"
+
 
 #include "MLog.h"
 #include "MLogManip.h"
+
+#include "MReportDrive.h"
+#include "MPointingPos.h"
+#include "MSrcPosCam.h"
 
 #include "MStarCam.h"
 #include "MStarPos.h"
 #include "MSkyCamTrans.h"
+#include "MStarCamTrans.h"
 
 ClassImp(MTelAxisFromStars);
@@ -82 +87 @@
     // type 1 : result from the Gauss fit
     fInputType = 1;
+
+    // default value of fAberr
+    // the value 1.07 is valid if the expected position (with aberration)
+    // in the camera is calculated as the average of the positions obtained
+    // from the reflections at the individual mirrors
+    fAberr = 1.07;
+
 }
 
@@ -90 +102 @@
 MTelAxisFromStars::~MTelAxisFromStars()
 {
+  delete fStarCamTrans;
 }
 
@@ -99 +112 @@
 Int_t MTelAxisFromStars::PreProcess(MParList *pList)
 {
+   fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));
+   if (!fDrive)
+   {
+       *fLog << err << AddSerialNumber("MReportDrive") 
+             << " not found... aborting." << endl;
+       return kFALSE;
+   }
+ 
   
    fStarCam = (MStarCam*)pList->FindObject("MStarCam", "MStarCam");
@@ -115 +136 @@
 
 
+   //-----------------------------------------------------------------
     fSrcPos = (MSrcPosCam*)pList->FindCreateObj(AddSerialNumber("MSrcPosCam"));
     if (!fSrcPos)
@@ -122 +144 @@
     }
 
+    fPointPos = (MPointingPos*)pList->FindCreateObj(AddSerialNumber("MPointingPos"), "MPointingPos");
+    if (!fPointPos)
+    {
+        *fLog << err << "MTelAxisFromStars::PreProcess;  MPointingPos not found...  aborting" << endl;
+        return kFALSE;
+    }
+
+    fSourcePos = (MPointingPos*)pList->FindCreateObj(AddSerialNumber("MPointingPos"), "MSourcePos");
+    if (!fSourcePos)
+    {
+        *fLog << err << "MTelAxisFromStars::PreProcess;  MSourcePos[MPointingPos] not found...  aborting" << endl;
+        return kFALSE;
+    }
+
     fSkyCamTrans = (MSkyCamTrans*)pList->FindCreateObj(AddSerialNumber("MSkyCamTrans"));
     if (!fSkyCamTrans)
@@ -129 +165 @@
     }
 
+   //-----------------------------------------------------------------
+   // book an MStarCamTrans object
+   // this is needed when calling one of the member functions of MStarCamTrans
+
+   MGeomCam *geom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
+   if (!geom)
+   {
+       *fLog << err << AddSerialNumber("MGeomCam") 
+             << " not found... aborting." << endl;
+       return kFALSE;
+   }
+
+   MObservatory *obs = (MObservatory*)pList->FindObject(AddSerialNumber("MObservatory"));
+   if (!obs)
+   {
+       *fLog << err << AddSerialNumber("MObservatory") 
+             << " not found... aborting." << endl;
+       return kFALSE;
+   }
+
+    fStarCamTrans = new MStarCamTrans(*geom, *obs);
 
     return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Set optical aberration factor
+//
+//   fAberr   is the ratio between 
+//            the distance from the camera center with optical aberration and
+//            the distance from the camera center with an ideal imaging
+//
+//   fAberr = r_real/r_ideal 
+//
+void MTelAxisFromStars::SetOpticalAberr(Double_t aberr)
+{
+  *fLog << all   << "MTelAxisFromStars::SetOpticalAberr; the optical aberration factor is set equal to : " 
+        << aberr ;
+
+  fAberr = aberr;
 }
 
@@ -374 +449 @@
 
   //--------------------------------------
-  // Put the estimated position, obtained by transforming the expected
-  // position (0,0), into SrcPosCam
+  // Put the corrected pointing position into MPointingPos 
+  //
+  SetPointingPosition(fStarCamTrans, fDrive, fSkyCamTrans, fPointPos);
+
+
+  //--------------------------------------
+  // Put the estimated position of the source into SrcPosCam
+  //
+  // get the source direction from MReportDrive
   // Note : this has to be changed for the wobble mode, where the source 
-  //        won't be in the center of the camera
-
-  fSrcPos->SetXY(fD[0], fD[1]);
-  fSrcPos->SetReadyToSave();       
-
+  //        isn't in the center of the camera
+  Double_t decsource = fDrive->GetDec();
+  Double_t rasource  = fDrive->GetRa();
+  //
+  Double_t radrive = fDrive->GetRa();
+  Double_t hdrive  = fDrive->GetHa();
+  Double_t hsource = hdrive + radrive - rasource;
+  fSourcePos->SetSkyPosition(rasource, decsource, hsource);
+
+  SetSourcePosition(fStarCamTrans, fPointPos, fSourcePos, fSrcPos);
+  //fSrcPos->SetXY(fD[0], fD[1]);
+  //fSrcPos->SetReadyToSave();       
+  *fLog << "after calling SetSourcePosition : , X, Y ,fD = " 
+        << fSrcPos->GetX() << ",  " << fSrcPos->GetY() << ",  "
+        << fD[0] << ",  " << fD[1] << endl;
 
   //--------------------------------------
@@ -448 +540 @@
   return kTRUE;
 }
+
+
 
 //---------------------------------------------------------------------------
@@ -1045 +1139 @@
 // --------------------------------------------------------------------------
 //
+// SetPointingPosition
+//
+//    put the corrected pointing direction into MPointingPos[MPointingPos];
+//    this direction corresponds to the position (0,0) in the camera
+//
+
+Bool_t MTelAxisFromStars::SetPointingPosition(MStarCamTrans *fstarcamtrans,
+       MReportDrive *fdrive, MSkyCamTrans *ftrans, MPointingPos *fpointpos)
+{
+  Double_t decdrive = fdrive->GetDec();
+  Double_t hdrive   = fdrive->GetHa();
+  Double_t radrive  = fdrive->GetRa();
+
+  // this is the estimated position (with optical aberration) in the camera 
+  // corresponding to the direction in MReportDrive
+  Double_t Xpoint = (ftrans->GetShiftD())[0];    
+  Double_t Ypoint = (ftrans->GetShiftD())[1];    
+
+  // get the sky direction corresponding to the position (0,0) in the camera
+  Double_t decpoint = 0.0;
+  Double_t hpoint   = 0.0;
+  fstarcamtrans->CelCamToCel0(decdrive, hdrive, 
+                         Xpoint/fAberr, Ypoint/fAberr, decpoint, hpoint);
+  Double_t rapoint = radrive - hpoint + hdrive;
+  fpointpos->SetSkyPosition(rapoint, decpoint, hpoint);
+
+  // get the local direction corresponding to the position (0,0) in the camera
+  Double_t thetadrive = fdrive->GetNominalZd();
+  Double_t phidrive   = fdrive->GetNominalAz();
+  Double_t thetapoint = 0.0;
+  Double_t phipoint   = 0.0;
+  fstarcamtrans->LocCamToLoc0(thetadrive, phidrive, 
+               Xpoint/fAberr, Ypoint/fAberr, thetapoint, phipoint); 
+  fpointpos->SetLocalPosition(thetapoint, phipoint);
+  fpointpos->SetReadyToSave();
+
+  *fLog << "SetPointingPosition : decdrive, hdrive, radrive Xpoint, Ypoint = "
+        << decdrive << ",  " << hdrive << ",  " << radrive << ",  "
+        << Xpoint << ",  " << Ypoint << endl;
+
+  *fLog << "SetPointingPosition : thetadrive, phidrive, thetapoint, phipoint = "       
+        << thetadrive << ",  " << phidrive << ",  " << thetapoint << ",  "
+        << phipoint << endl;
+
+  return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// SetSourcePosition
+//
+//    put the estimated position of the source in the camera into 
+//    MSrcPosCam[MSrcPosCam]
+//
+//    and the estimated local direction of the source into
+//    MSourcePos[MPointingPos]
+//
+
+Bool_t MTelAxisFromStars::SetSourcePosition(MStarCamTrans *fstarcamtrans,
+       MPointingPos *fpointpos, MPointingPos *fsourcepos, MSrcPosCam *fsrcpos)
+{
+  // get the corrected pointing direction
+  // corresponding to the position (0,0) in the camera
+  Double_t decpoint = fpointpos->GetDec();
+  Double_t hpoint   = fpointpos->GetHa();
+
+  // get the sky direction of the source
+  Double_t decsource = fsourcepos->GetDec();
+  Double_t hsource   = fsourcepos->GetHa();
+
+  // get the estimated position (Xsource, Ysource) of the source in the camera;
+  // this is a position for an ideal imaging, without optical aberration
+  Double_t Xsource = 0.0;
+  Double_t Ysource = 0.0;
+  fstarcamtrans->Cel0CelToCam(decpoint,  hpoint, 
+                              decsource, hsource, Xsource, Ysource);
+  fsrcpos->SetXY(Xsource*fAberr, Ysource*fAberr);
+  fsrcpos->SetReadyToSave();
+
+  // get the estimated local direction of the source
+  Double_t thetapoint = fpointpos->GetZd();
+  Double_t phipoint   = fpointpos->GetAz();
+  Double_t thetasource = 0.0;
+  Double_t phisource   = 0.0;
+  fstarcamtrans->Loc0CamToLoc(thetapoint, phipoint, 
+                              Xsource, Ysource, thetasource, phisource);
+  fsourcepos->SetLocalPosition(thetasource, phisource);
+  fsourcepos->SetReadyToSave();
+
+  *fLog << "SetSourcePosition : decpoint, hpoint, decsource, hsource, Xsource, Ysource = "
+        << decpoint << ",  " << hpoint << ",  " << decsource << ",  "
+        << hsource  << ",  " << Xsource << ",  " << Ysource << endl;
+  *fLog << "SetSourcePosition : thetapoint, phipoint, thetasource, phisource = "
+        << thetapoint << ",  " << phipoint << ",  " << thetasource << ",  "
+        << phisource << endl;
+
+  return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
 //
 Int_t MTelAxisFromStars::PostProcess()

trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.h

@@ -20 +20 @@
 #endif
 
+class MReportDrive;
+class MPointingPos;
 class MSrcPosCam;
 class MStarCam;
 class MSkyCamTrans;
+class MStarCamTrans;
 
 class MTelAxisFromStars : public MTask
@@ -31 +34 @@
     const MStarCam *fSourceCam;        //!
 
+    MReportDrive        *fDrive;          //!
+    MPointingPos        *fPointPos;       //!
+    MPointingPos        *fSourcePos;      //!
     MSrcPosCam          *fSrcPos;         //!
     MSkyCamTrans        *fSkyCamTrans;    //!
+    MStarCamTrans       *fStarCamTrans;   //!
 
     Double_t fFixedScaleFactor;           //!
     Double_t fFixedRotationAngle;         //! [degrees]
+    Double_t fAberr;                      //! optical aberration factor
 
     Int_t fInputType;                     //! type of input
@@ -50 +58 @@
     void FixScaleFactorAt(Double_t lambda = 1.0);
     void FixRotationAngleAt(Double_t alfa = 0.0);  // alfa in [degrees]
-   
-    void SetInputType(Int_t type = 2);
+
+    void SetOpticalAberr(Double_t aberr);   
+    void SetInputType(Int_t type = 1);
 
     Bool_t FindSkyCamTrans(TArrayD[2],      TArrayD[2],  TArrayD[2][2], 
@@ -60 +69 @@
     void TransSkyCam(Double_t &,  Double_t[2][2], Double_t[2], Double_t[2][2],
                      TArrayD[2],  TArrayD[2],     TArrayD[2][2]);   
+
+    Bool_t SetPointingPosition(MStarCamTrans *fstarcamtrans,
+       MReportDrive *fdrive, MSkyCamTrans *ftrans, MPointingPos *fpointpos);
+    Bool_t SetSourcePosition(MStarCamTrans *fstarcamtrans,
+       MPointingPos *fpointpos, MPointingPos *fsourcepos, MSrcPosCam *fsrcpos);
 
     ClassDef(MTelAxisFromStars, 0) // Task to calculate the source position from star positions

trunk/MagicSoft/Mars/mtemp/findTelAxisFromStars.C

@@ -15 +15 @@
 }
 
+//--------------------------------------------------------------------------
+//
+//  findTelAxisFromStars
+//
+//  macro to 
+//    - determine the positions of stars in the camera from the DC currents 
+//    - and to find the expected position of the source in the camera
+//
+
 void findTelAxisFromStars(const TString filename="20040422_23213_D_Mrk421_E.root", const TString directory="/.magic/magicserv01/MAGIC/rootdata/2004_04_22/", const UInt_t numEvents = 0)
 {
   gLog.SetNoColors();
 
+  //-------------------------------------------------------------------
+  // book some histograms which are to be filled in the event loop
+  //
+
+  TH2D *aberr = new TH2D("OpticalAberr", "Opt_Aberr",
+                         100, 0.0, 650.0,  100, 0.0, 650.0);
+  aberr->SetXTitle("Ideal distance from center [mm]");
+  aberr->SetYTitle("Expected distance from center [mm]");
+
+  TProfile *reldiff = new TProfile("RelativeDiff", "Rel_Diff",
+                                   100, 0.0, 650.0, -0.2, 0.2, "S");
+  reldiff->SetXTitle("Ideal distance from center [mm]");
+  reldiff->SetYTitle("(R - R0)/R0");
+
+  //-------------------------------------------------------------------
+
+
   MParList  plist;
   MTaskList tlist;
@@ -28 +54 @@
 
   //$$$$$$$$$$$$$$$$ ww
+  MObservatory  obs;
   MStarCam      sourcecam;
   sourcecam.SetName("MSourceCam");
 
   MHTelAxisFromStars htelaxis;
+
+  plist.AddToList(&obs);
   plist.AddToList(&htelaxis);
   //$$$$$$$$$$$$$$$$ ww
@@ -99 +128 @@
   const Double_t dec = MAstro::Dms2Rad(38, 14, 29);
   sdirs.AddDirection(ra,dec,1,"My_UMa 51");
+
+  // dummy source dirtections
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(38, 20,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy1");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(38, 20,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy2");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(38, 30,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy3");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(38, 40,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy4");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(38, 50,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy5");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(39,  0,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy6");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(39, 10,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy7");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(39, 20,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy8");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(39, 30,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy9");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(39, 40,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy10");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(39, 50,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy11");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(40,  0,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy12");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(40, 10,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy13");
+
+  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
+  const Double_t dec = MAstro::Dms2Rad(40, 20,  0);
+  sdirs.AddDirection(ra,dec,1,"Dummy14");
 
 
@@ -195 +281 @@
         // This line prints the results:
         //        starcam.Print();
+
         // This is how to access the TList of stars:
-        //        TList* starlist = starcam.GetList();
+        //        TList* starlist = starcam.GetList();
 
         // This is how to iterate over stars found:
-        //        TIter Next(starlist);
+        //        TIter Next(starlist);
         //        MStarPos* star;
         //        UInt_t starnum = 0;
-        //        cout << filename << " ";
-        //        cout << "Iterating over list" << endl;
-        //        while ((star=(MStarPos*)Next())) 
-        //          {
-        //            cout << "star[" << starnum << "] ";
-        //            cout << star->GetMeanX() << " " 
-        //                 << star->GetMeanY() << " ";
-        //            starnum++;
-        //          }
-        //        cout << endl;
+        //        cout << filename << " ";
+        //        cout << "Iterating over list" << endl;
+        //        while ( (star=(MStarPos*)Next()) ) 
+        //        {
+        //          cout << "star[" << starnum << "] ";
+        //          cout << star->GetMeanX() << " " 
+        //               << star->GetMeanY() << " ";
+        //          starnum++;
+        //        }
+        //        cout << endl;
+        }
+
+      // fill the histograms
+      // loop over stars in "MSourceCam"
+      TList* sourcelist = sourcecam.GetList();
+      TIter Nextsource(sourcelist);
+      MStarPos* source;
+      while ( (source=(MStarPos*)Nextsource()) ) 
+      {
+        Double_t R = sqrt( source->GetXExp()*source->GetXExp()
+                   + source->GetYExp()*source->GetYExp() );
+
+        Double_t R0= sqrt( source->GetXIdeal()*source->GetXIdeal()
+                   + source->GetYIdeal()*source->GetYIdeal() );
+
+        cout << "R0, R = " << R0 << ",  " << R << endl;
+
+        aberr->Fill(R0, R, 1.0);
+        reldiff->Fill(R0, (R-R0)/R0), 1.0);
+      }
         
-        }//integratedevents
-
+      //--------------------------------------
       MTime time;
       time.Set(2004, 4, 22, 21, 51, 15);
@@ -225 +331 @@
       
       // wait after each event
-      if (!HandleInput())
-        break;
+      //if (!HandleInput())
+      //  break;
 
     }
@@ -256 +362 @@
   //$$$$$$$$$$$$$$$$ ww
 
-
+  // plot the histograms
+  TCanvas *canv = new TCanvas("Optical_Aberration", "Optical aberration",
+                              0, 0, 600, 300);
+  canv->Divide(2,1);
+  canv->SetBorderMode(0); 
+  
+  canv->cd(1);
+  aberr->Draw();
+
+  canv->cd(2);
+  reldiff->Draw();
 
 }
 
+//----------------------------------------------------------------------
+//
+//
 
 Bool_t HandleInput()
@@ -283 +402 @@
     return kFALSE;
 }
+//=========================================================================
+
+
+
+
+
+
+
+
+