Changeset 9070 for trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc

Timestamp:

08/02/08 15:49:41 (16 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc (modified) (22 diffs)

trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc

@@ -17 +17 @@
 !
 !   Author(s): Thomas Bretz 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
-!   Author(s): Abelardo Moralejo 1/2005 <mailto:moralejo@pd.infn.it>
 !
-!   Copyright: MAGIC Software Development, 2000-2005
+!   Copyright: MAGIC Software Development, 2000-2008
 !
 !
@@ -43 +42 @@
 // coordinates of the observatory from MObservatory.
 //
+// If a reflector version >= 700 is detected the source position for
+// Gammas (and _only_ gammas) is dynamically calculated from the
+// telescope (e.g. MMcEvt::fTelescopePhi) and the shower (e.g. MMcEvt::fPhi)
+// orientation. In off-mode it is fixed at the camera center.
+//
 // FIXME: Add here some information about the more-cycle calculation
 //
@@ -50 +54 @@
 //   MGeomCam
 //   MTime
+//   [MMcRunHeader]
+//   [MMcCorsikaRunHeader]
 //   [MSourcePos] (of type MPointingPos)
 //
@@ -57 +63 @@
 // To be done:
 //   - a switch between using sky-coordinates and time or local-coordinates
-//     from MPointingPos for determine the rotation angle
-/////  NOTE, A. Moralejo: the above is not possible if MAstroSky2Local does not 
-/////  account for precession and nutation.
+//     from MPointingPos for determine the rotation angle. (Therefor a conversion
+//     needs to be implemented inlcuding asrometric corrections)
 //
 //   - the center of rotation need not to be the camera center
@@ -82 +87 @@
 #include "MSrcPosCam.h"
 #include "MRawRunHeader.h"
+
+#include "MMcEvt.hxx"
+#include "MMcRunHeader.hxx"
 #include "MMcCorsikaRunHeader.h"
 
@@ -97 +105 @@
 //
 MSrcPosCalc::MSrcPosCalc(const char *name, const char *title)
-    : fObservatory(NULL), fPointPos(NULL), fSourcePos(NULL), fDeviation(NULL),
-    fSrcPosCam(NULL), fSrcPosAnti(NULL), fGeom(NULL), fTime(NULL), fCallback(NULL),
-    fMode(kDefault), fNumRandomOffPositions(0)
+    : fObservatory(NULL), fPointPos(NULL), fDeviation(NULL), fMcEvt(NULL),
+    fMcHeader(NULL), fGeom(NULL), fTime(NULL), fCallback(NULL),
+    fSourcePos(NULL), fSrcPosCam(NULL), fSrcPosAnti(NULL), fMode(kDefault),
+    fNumRandomOffPositions(0)
 {
     fName  = name  ? name  : "MSrcPosCalc";
@@ -105 +114 @@
 
     AddToBranchList("MTime.*");
+    AddToBranchList("MMcEvt.*");
     AddToBranchList("MPointingPos.*");
 }
@@ -170 +180 @@
     if (!fSrcPosAnti)
         return kFALSE;
+
+    fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
+    if (!fGeom)
+    {
+        *fLog << err << "MGeomCam not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fPointPos = (MPointingPos*)pList->FindObject("MPointingPos");
+    if (!fPointPos)
+    {
+        *fLog << err << "MPointingPos not found... aborting." << endl;
+        return kFALSE;
+    }
 
     if (!fSourcePos)
@@ -184 +208 @@
     }
     // FIXME: Maybe we have to call FreeSourcePos in PostProcess()?
-
-    fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
-    if (!fGeom)
-    {
-        *fLog << err << "MGeomCam not found... aborting." << endl;
-        return kFALSE;
-    }
-
-    fPointPos = (MPointingPos*)pList->FindObject("MPointingPos");
-    if (!fPointPos)
-    {
-        *fLog << err << "MPointingPos not found... aborting." << endl;
-        return kFALSE;
-    }
 
     fDeviation = (MPointingDev*)pList->FindObject("MPointingDev");
@@ -262 +272 @@
 }
 
+void MSrcPosCalc::InitFixedPos() const
+{
+    // It is set here for the cases in which Process is not called at all
+    fSrcPosCam->SetXY(fFixedPos/fGeom->GetConvMm2Deg());
+
+    // Informal message output
+    *fLog << inf;
+    *fLog << "Source position set to x=" << fFixedPos.X() << "deg ";
+    *fLog << "y=" << fFixedPos.Y() << "deg" << endl;
+}
+
 // --------------------------------------------------------------------------
 //
@@ -269 +290 @@
 Bool_t MSrcPosCalc::ReInit(MParList *plist)
 {
+    // In OffMode set fixed position to camera center
     if (fMode==kOffData)
-    {
-        // Set fixed position to camera center
         fFixedPos = TVector2();
 
+    if (fMode==kOffData || fMode==kFixed)
+    {
         // It is set here for the cases in which Process is not called at all
-        fSrcPosCam->SetXY(fFixedPos);
+        InitFixedPos();
         return kTRUE;
     }
@@ -302 +324 @@
             return kFALSE;
         }
+
         return kTRUE;
     }
 
-    MMcCorsikaRunHeader *h = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
+    fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
+    if (!fMcEvt)
+    {
+        *fLog << err << "MMcEvt not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fMcHeader = (const MMcRunHeader*)plist->FindObject("MMcRunHeader");
+    if (!fMcHeader)
+    {
+        *fLog << err << "MMcRunHeader not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    const MMcCorsikaRunHeader *h = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
     if (!h)
     {
@@ -311 +348 @@
         return kFALSE;
     }
+
+    // We don't know here if these are gammas
+    //if (fMcHeader->GetReflVersion()>600)
+    //    *fLog << inf << "Source position will be calculated from shower and telescope orientation." << endl;
 
     // Determine Monte Carlo position from Monte Carlo header
     TVector2 v(0, 0);
     if (h->GetWobbleMode()>0.5)
-        v.Set(120., 0.);
+        v.Set(120.*fGeom->GetConvMm2Deg(), 0.);
     if (h->GetWobbleMode()<-0.5)
-        v.Set(-120., 0.);
+        v.Set(-120.*fGeom->GetConvMm2Deg(), 0.);
 
     // Set fixed position
     fFixedPos = v;
 
-    // It is set here for the cases in which Process is not called at all
-    fSrcPosCam->SetXY(fFixedPos);
-
-    // Informal message output
-    *fLog << inf;
-    *fLog << "Source Position set to x=" << fFixedPos.X() << "mm ";
-    *fLog << "y=" << fFixedPos.Y() << "mm" << endl;
+    InitFixedPos();
 
     return kTRUE;
 }
 
+void MSrcPosCalc::CalcResult(const MVector3 &pos0, const MVector3 &pos)
+{
+    // Calculate source position in camera, and convert to mm:
+    TVector2 v = MAstro::GetDistOnPlain(pos0, pos, -fGeom->GetCameraDist()*1000);
+
+    if (fDeviation)
+        v -= fDeviation->GetDevXY()/fGeom->GetConvMm2Deg();
+
+    SetSrcPos(v);
+}
+
 // --------------------------------------------------------------------------
 //
@@ -339 +385 @@
 Int_t MSrcPosCalc::Process()
 {
-    if (fRunType==MRawRunHeader::kRTMonteCarlo || !fSourcePos || !fTime || !fObservatory || fMode==kOffData)
-    {
-        SetSrcPos(fFixedPos);
+    // In off-data mode the position is set to a fixed value
+    // independent of the run type. In fixed-data mode the position
+    // is set to the predefined position converted to mm.
+    if (fMode==kOffData || fMode==kFixed)
+    {
+        SetSrcPos(fFixedPos/fGeom->GetConvMm2Deg());
+        return kTRUE;
+    }
+
+    // If it is a monte carlo run calculate source position from
+    // Monte Carlo headers
+    if (fRunType==MRawRunHeader::kRTMonteCarlo)
+    {
+        // If the reflector version is too old take source position
+        // from the WobbleMode in the header
+        if (fMcHeader->GetReflVersion()<=600)
+        {
+            SetSrcPos(fFixedPos/fGeom->GetConvMm2Deg());
+            return kTRUE;
+        }
+
+        // If the reflector version was new enough calculate the
+        // source position from shower and telescope orientation
+        // FIXME: Copy fMcEvt to fSourcePos!
+        MVector3 pos0, pos;
+        pos0.SetZdAz(fPointPos->GetZdRad(), fPointPos->GetAzRad());
+        pos.SetZdAz(fMcEvt->GetTheta(), fMcEvt->GetPhi());
+
+        CalcResult(pos0, pos);
+        return kTRUE;
+    }
+
+    // If the user requested a fixed source position use this position
+    if (!fSourcePos || !fTime || !fObservatory)
+    {
+        SetSrcPos(fFixedPos/fGeom->GetConvMm2Deg());
         return kTRUE;
     }
@@ -352 +431 @@
     pos.SetRaDec(fSourcePos->GetRaRad(), fSourcePos->GetDecRad());
 
+    // Set sky coordinates of camera center in pos0 (for details see below)
+    MVector3 pos0;  // pos0: camera center
+    pos0.SetRaDec(fPointPos->GetRaRad(), fPointPos->GetDecRad());
+
+    // Here we check that the source position which was given
+    // by the drive system is not off by more than 1deg from the
+    // source position given by the user. (Checking every individual
+    // event is not the fastest, but the safest)
+    if (pos.DeltaPhi(pos0)*TMath::RadToDeg()>1)
+    {
+        *fLog << err << "ERROR - Defined source position deviates from pointing-ra/dec by more than 1deg!" << endl;
+        *fLog << "User defined source pos:  ";
+        fSourcePos->Print();
+        *fLog << "Pointing position of tel: ";
+        fPointPos->Print();
+
+        return kFALSE;
+    }
+
     // Convert sky coordinates of source to local coordinates. Warning! These are not the "true" local
     // coordinates, since this transformation ignores precession and nutation effects.
@@ -357 +455 @@
     pos *= conv;
 
-    // Set sky coordinates of camera center in pos0, and then convert to
-    // local. Same comment as above. These coordinates differ from the true
+    // Convert sky coordinates of camera center convert to local.
+    // Same comment as above. These coordinates differ from the true
     // local coordinates of the camera center that one could get from
     // "MPointingPos", calculated by the Drive: the reason is that the Drive
@@ -369 +467 @@
     // would be wrong is to set in pos0 directly the local coordinates
     // found in MPointingPos!
-    MVector3 pos0;  // pos0: camera center
-    pos0.SetRaDec(fPointPos->GetRaRad(), fPointPos->GetDecRad());
     pos0 *= conv;
 
@@ -384 +480 @@
     }
 
-    // Calculate source position in camera, and convert to mm:
-    TVector2 v = MAstro::GetDistOnPlain(pos0, pos, -fGeom->GetCameraDist()*1000);
-
-    if (fDeviation)
-        v -= fDeviation->GetDevXY()/fGeom->GetConvMm2Deg();
-
-    SetSrcPos(v);
-
+    CalcResult(pos0, pos);
     return kTRUE;
 }
@@ -440 +529 @@
 Int_t MSrcPosCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
 {
+    Bool_t rc = kFALSE;
+
+    if (IsEnvDefined(env, prefix, "FixedPos", print))
+    {
+        TString str = GetEnvValue(env, prefix, "FixedPos", "");
+        str = str.Strip(TString::kBoth);
+
+        Double_t x, y;
+        const Int_t n=sscanf(str.Data(), "%lf %lf", &x, &y);
+        if (n!=2)
+        {
+            *fLog << warn << "WARNING - FixedPos requires two numbers." << endl;
+            return kERROR;
+        }
+
+        fFixedPos = TVector2(x, y);
+        fMode = kFixed;
+
+        rc = kTRUE;
+    }
+
     Double_t ra=0;
     Double_t dec=0;
@@ -461 +571 @@
 
         SetSourcePos(ra, dec);
-        return kTRUE;
-    }
-
-    Bool_t rc = kFALSE;
+        return rc;
+    }
+
+    Bool_t rcsrc = kFALSE;
     if (IsEnvDefined(env, prefix, "SourceRa", print))
     {
@@ -472 +582 @@
             return kERROR;
 
-        rc = kTRUE;
+        rcsrc = kTRUE;
     }
 
@@ -482 +592 @@
             return kERROR;
 
-        rc = kTRUE;
-    }
-
-    if (!rc)
-        return kFALSE;
-
-    SetSourcePos(ra, dec);
-    return kTRUE;
-}
+        rcsrc = kTRUE;
+    }
+
+    if (rcsrc)
+        SetSourcePos(ra, dec);
+
+    return rcsrc || rc;
+}