Changeset 3340 for trunk

Timestamp:

02/27/04 14:40:57 (21 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhist

Files:

• MHStarMap.cc (modified) (7 diffs)
• MHStarMap.h (modified) (4 diffs)

trunk/MagicSoft/Mars/mhist/MHStarMap.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
-!
-!   Copyright: MAGIC Software Development, 2000-2002
+!   Author(s): Thomas Bretz    12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
+!              Wolfgang Wittek 02/2004 <mailto:wittek@mppmu.mpg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2004
 !
 !
@@ -31 +32 @@
 // from the Hillas parameters (Fill) can be enhanced.
 //
+// For a given a shower, a series of points along its main axis are filled 
+// into the 2-dim histogram (x, y) of the camera plane.
+// 
+// Before filling a point (x, y) into the histogram it is
+//        - shifted by (xSrc, ySrc)   (the expected source position)
+//        - and rotated in order to compensate the rotation of the
+//          sky image in the camera
+//         
+//         
 /////////////////////////////////////////////////////////////////////////////
 #include "MHStarMap.h"
@@ -47 +57 @@
 #include "MHillas.h"
 #include "MBinning.h"
+#include "MMcEvt.hxx"
+#include "MSrcPosCam.h"
+#include "MObservatory.h"
 
 ClassImp(MHStarMap);
@@ -117 +130 @@
     }
 
+   fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
+   if (!fMcEvt)
+   {
+       *fLog << err << "MMcEvt not found... aborting." << endl;
+       return kFALSE;
+   }
+
+    fSrcPos = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam"));
+    if (!fSrcPos)
+    {
+        *fLog << err << "MSrcPosCam not found...  aborting" << endl;
+        return kFALSE;
+    }
+
+
     const MBinning *bins = (MBinning*)plist->FindObject("BinningStarMap");
     if (!bins)
@@ -134 +162 @@
     if (!geom)
     {
-        *fLog << warn << dbginf << "No Camera Geometry available. Using mm-scale for histograms." << endl;
+        *fLog << warn << "No Camera Geometry available. Using mm-scale for histograms." << endl;
         return kTRUE;
     }
@@ -141 +169 @@
 
     return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Get the observatory location "MObservatory" from parameter list
+// 
+Bool_t MHStarMap::ReInit(MParList *pList)
+{
+    MObservatory *obs = (MObservatory*)pList->FindObject(AddSerialNumber("MObservatory"));
+    if (!obs)
+    {
+        *fLog << err << "MObservatory not found...  aborting" << endl;
+        return kFALSE;
+    }
+
+    fSinLat = TMath::Sin(obs->GetPhi());
+    fCosLat = TMath::Cos(obs->GetPhi());
+
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// RotationAngle 
+// 
+// calculates the angle for the rotation of the sky image in the camera;
+// this angle is a function of the local coordinates (theta, phi) in rad.
+//
+// calculate rotation angle alpha of sky image in camera
+// (see TDAS 00-11, eqs. (18) and (20))
+//
+void MHStarMap::GetRotationAngle(Double_t &cosal, Double_t &sinal)
+{
+    const Double_t theta = fMcEvt->GetTelescopeTheta();
+    const Double_t phi   = fMcEvt->GetTelescopePhi();
+
+    const Double_t sint = TMath::Sin(theta);
+    const Double_t cost = TMath::Cos(theta);
+
+    const Double_t sinp = TMath::Sin(phi);
+    const Double_t cosp = TMath::Cos(phi);
+
+    const Double_t v1 = sint*sinp;
+    const Double_t v2 = fCosLat*cost - fSinLat*sint*cosp;
+
+    const Double_t denom =  1./ TMath::Sqrt(v1*v1 + v2*v2);
+
+    cosal = (fSinLat * sint) + fCosLat * cost * cosp * denom;
+    sinal = (fCosLat * sinp) * denom;
 }
 
@@ -152 +229 @@
     const MHillas &h = *(MHillas*)par;
 
-    const float delta = h.GetDelta();
-
-    const float m    = tan(delta);
-    const float cosd = 1.0/sqrt(1.0+m*m);
-    const float sind = sqrt(1.0-cosd*cosd);
-
-    float t = h.GetMeanY() - m*h.GetMeanX();
+    const Float_t delta = h.GetDelta();
+
+    const Float_t m    = TMath::Tan(delta);
+
+    const Float_t cosd = 1.0/TMath::Sqrt(m*m+1);
+    const Float_t sind = TMath::Sqrt(1-cosd*cosd);
+
+    Float_t t = h.GetMeanY() - m*h.GetMeanX();
+
+    Float_t xSource = fSrcPos->GetX();
+    Float_t ySource = fSrcPos->GetY();
 
     if (!fUseMmScale)
+    {
         t *= fMm2Deg;
+        xSource *= fMm2Deg;
+        ySource *= fMm2Deg;
+    }
 
     // get step size ds along the main axis of the ellipse
     const TAxis &axex = *fStarMap->GetXaxis();
-    const float xmin = axex.GetXmin();
-    const float xmax = axex.GetXmax();
+    const Float_t xmin = axex.GetXmin();
+    const Float_t xmax = axex.GetXmax();
 
     // FIXME: Fixed number?
-    const float ds = (xmax-xmin) / 200.0;
+    const Float_t ds = (xmax-xmin) / 200.0;
+
+    // Fill points along the main axis of the shower into the histogram;
+    // before filling 
+    //   - perform a translation by (xSource, ySource)
+    //   - and perform a rotation to compensate the rotation of the
+    //     sky image in the camera
+    Double_t cosal;
+    Double_t sinal;
+    GetRotationAngle(cosal, sinal);
 
     if (m>-1 && m<1)
     {
-        const float dx = ds * cosd;
+        const Float_t dx = ds * cosd;
 
         for (float x=xmin+dx/2; x<(xmax-xmin)+dx; x+=dx)
-            fStarMap->Fill(x, m*x+t, w);
+        {
+            const Float_t xorig = x     - xSource;
+            const Float_t yorig = m*x+t - ySource;
+
+            const Float_t xfill =  cosal*xorig + sinal*yorig;
+            const Float_t yfill = -sinal*xorig + cosal*yorig;
+            fStarMap->Fill(xfill, yfill, w);
+        }
     }
     else
     {
         const TAxis &axey = *fStarMap->GetYaxis();
-        const float ymin = axey.GetXmin();
-        const float ymax = axey.GetXmax();
+        const Float_t ymin = axey.GetXmin();
+        const Float_t ymax = axey.GetXmax();
 
         const float dy = ds * sind;
 
         for (float y=ymin+dy/2; y<(ymax-ymin)+dy; y+=dy)
-            fStarMap->Fill((y-t)/m, y, w);
+        {
+            const Float_t xorig = (y-t)/m - xSource;
+            const Float_t yorig = y       - ySource;
+
+            const Float_t xfill =  cosal*xorig + sinal*yorig;
+            const Float_t yfill = -sinal*xorig + cosal*yorig;
+            fStarMap->Fill(xfill, yfill, w);
+        }
     }
 

trunk/MagicSoft/Mars/mhist/MHStarMap.h

@@ -8 +8 @@
 class TH2F;
 class MHillas;
+class MSrcPosCam;
+class MMcEvt;
 
 class MHStarMap : public MH
 {
 private:
-    TH2F *fStarMap; //->
+    MSrcPosCam *fSrcPos; //!
+    MMcEvt     *fMcEvt;  //!
+
+    TH2F *fStarMap;      //->
 
     Float_t fMm2Deg;
@@ -18 +23 @@
     Bool_t fUseMmScale;
 
+    Float_t fCosLat; //!
+    Float_t fSinLat; //!
+
     void PrepareDrawing() const;
 
     void Paint(Option_t *opt="");
+
+    void GetRotationAngle(Double_t &cosangle, Double_t &sinangle);
+
+    Bool_t SetupFill(const MParList *pList);
+    Bool_t Fill(const MParContainer *par, const Stat_t w=1);
+    Bool_t ReInit(MParList *pList);
 
 public:
@@ -29 +43 @@
     void SetMm2Deg(Float_t mmdeg);
 
-    Bool_t SetupFill(const MParList *pList);
-    Bool_t Fill(const MParContainer *par, const Stat_t w=1);
-
-    TH1 *GetHistByName(const TString name) { return (TH1*)fStarMap; }
+     TH1 *GetHistByName(const TString name) { return (TH1*)fStarMap; }
 
     TH2F *GetHist() { return fStarMap; }
@@ -39 +50 @@
     TObject *DrawClone(Option_t *opt=NULL) const;
 
-    ClassDef(MHStarMap, 1) // Container to hold 2-dim histogram (starmap)
+    ClassDef(MHStarMap, 1) // Container to hold the Starmap
 };