Changeset 8178 for trunk/MagicSoft/Mars

Timestamp:

10/30/06 12:46:13 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (2 diffs)
• mbase/MMath.cc (modified) (3 diffs)
• mbase/MMath.h (modified) (2 diffs)
• mgui/MHexagon.cc (modified) (9 diffs)
• mgui/MHexagon.h (modified) (3 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -19 +19 @@
                                                  -*-*- END OF LINE -*-*-
 
+ 2006/10/30 Thomas Bretz
+
+   * mbase/MMath.[h,cc]:
+     - added a new function to calculate the intersection point of
+       two lines
+
+   * mgui/MHexagon.[h,cc]:
+     - added function to calculate the intersection of two hexagons
+
+
+
  2006/10/30 Daniela Dorner
 
@@ -61 +72 @@
      - do not print observation time statistics if observation time is
        zero
+
+   * datacenter/macros/fillstar.C:
+     - do not fill arbitrary negative value
+
+   * mfbase/MF.cc:
+     - fixed a compiler warning about a comment
+
+   * mranforest/MRanForestCalc.cc:
+     - removed the appearance of an obsolete error message
+     - the printing of weights is now done less often
+
+   * msignal/MExtractTimeAndCharge.cc:
+     - added a comment line
 
    * mjobs/MJSpectrum.[h,cc]:

trunk/MagicSoft/Mars/mbase/MMath.cc

@@ -1 +1 @@
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MMath.cc,v 1.30 2006-10-30 12:46:12 tbretz Exp $
+! --------------------------------------------------------------------------
 !
 ! *
@@ -32 +34 @@
 #include "MMath.h"
 
+#ifndef ROOT_TVector2
+#include <TVector2.h>
+#endif
+
 #ifndef ROOT_TVector3
 #include <TVector3.h>
@@ -505 +511 @@
 // --------------------------------------------------------------------------
 //
+// Calculate the intersection of two lines defined by (x1;y1) and (x2;x2)
+// Returns the intersection point.
+//
+// It is assumed that the lines intersect. If there is no intersection
+// TVector2() is returned (which is not destinguishable from
+// TVector2(0,0) if the intersection is at the coordinate source)
+//
+// Formula from: http://mathworld.wolfram.com/Line-LineIntersection.html
+//
+TVector2 MMath::GetIntersectionPoint(const TVector2 &x1, const TVector2 &y1, const TVector2 &x2, const TVector2 &y2)
+{
+    TMatrix d(2,2);
+    d[0][0] = x1.X()-y1.X();
+    d[0][1] = x2.X()-y2.X();
+    d[1][0] = x1.Y()-y1.Y();
+    d[1][1] = x2.Y()-y2.Y();
+
+    const Double_t denom = d.Determinant();
+    if (denom==0)
+        return TVector2();
+
+    TMatrix l1(2,2);
+    TMatrix l2(2,2);
+
+    l1[0][0] = x1.X();
+    l1[0][1] = y1.X();
+    l2[0][0] = x2.X();
+    l2[0][1] = y2.X();
+
+    l1[1][0] = x1.Y();
+    l1[1][1] = y1.Y();
+    l2[1][0] = x2.Y();
+    l2[1][1] = y2.Y();
+
+    TMatrix a(2,2);
+    a[0][0] = l1.Determinant();
+    a[0][1] = l2.Determinant();
+    a[1][0] = x1.X()-y1.X();
+    a[1][1] = x2.X()-y2.X();
+
+    const Double_t X = a.Determinant()/denom;
+
+    a[1][0] = x1.Y()-y1.Y();
+    a[1][1] = x2.Y()-y2.Y();
+
+    const Double_t Y = a.Determinant()/denom;
+
+    return TVector2(X, Y);
+}
+
+// --------------------------------------------------------------------------
+//
 // Solves: x^2 + ax + b = 0;
 // Return number of solutions returned as x1, x2

trunk/MagicSoft/Mars/mbase/MMath.h

@@ -6 +6 @@
 #endif
 
+class TVector2;
 class TVector3;
 class TArrayD;
@@ -47 +48 @@
     Double_t InterpolParabCos(const TVector3 &vx, const TVector3 &vy, Double_t x);
 
+    TVector2 GetIntersectionPoint(const TVector2 &x1, const TVector2 &y1, const TVector2 &x2, const TVector2 &y2);
+
     inline Int_t SolvePol1(Double_t c, Double_t d, Double_t &x1)
     {

trunk/MagicSoft/Mars/mgui/MHexagon.cc

@@ -1 +1 @@
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MHexagon.cc,v 1.30 2006-10-30 12:46:13 tbretz Exp $
+! --------------------------------------------------------------------------
 !
 ! *
@@ -19 +21 @@
 !   Author(s): Harald Kornmayer 1/2001
 !
-!   Copyright: MAGIC Software Development, 2000-2002
+!   Copyright: MAGIC Software Development, 2000-2006
 !
 !
@@ -25 +27 @@
 
 //////////////////////////////////////////////////////////////////////////////
-//                                                                          //
-// MHexagon                                                                 //
-//                                                                          //
+//
+// MHexagon
+//
 //////////////////////////////////////////////////////////////////////////////
 
@@ -35 +37 @@
 #include <iostream>
 
-#include <TVirtualPad.h>  // gPad
-
-#include "MGeomPix.h"     // GetX
+#include <TVector2.h>       // TVector2
+#include <TVirtualPad.h>    // gPad
+#include <TOrdCollection.h> // TOrdCollection
+
+#include "MMath.h"
+#include "MGeomPix.h"       // GetX
 
 ClassImp(MHexagon);
 
 using namespace std;
+
+const Double_t MHexagon::fgCos60 = 0.5;        // TMath::Cos(60/TMath::RadToDeg());
+const Double_t MHexagon::fgSin60 = sqrt(3.)/2; // TMath::Sin(60/TMath::RadToDeg());
+
+const Double_t MHexagon::fgDx[6] = { fgCos60,   0.,          -fgCos60,  -fgCos60,   0.,            fgCos60 };
+const Double_t MHexagon::fgDy[6] = { fgSin60/3, fgSin60*2/3, fgSin60/3, -fgSin60/3, -fgSin60*2/3, -fgSin60/3 };
 
 // ------------------------------------------------------------------------
@@ -157 +168 @@
 // coordinates. Return -1 if inside.
 //
-Float_t MHexagon::DistanceToPrimitive(Float_t px, Float_t py)
+Float_t MHexagon::DistanceToPrimitive(Float_t px, Float_t py) const
 {
     //
@@ -176 +187 @@
       return disthex;
 
-    const static Double_t cos60 = TMath::Cos(60/kRad2Deg);
-    const static Double_t sin60 = TMath::Sin(60/kRad2Deg);
-
-    const Double_t dx2 = dx*cos60 + dy*sin60;
+    const Double_t dx2 = dx*fgCos60 + dy*fgSin60;
 
     if  (TMath::Abs(dx2) > fD/2.)
       return disthex;
 
-    const Double_t dx3 = dx*cos60 - dy*sin60;
+    const Double_t dx3 = dx*fgCos60 - dy*fgSin60;
 
     if  (TMath::Abs(dx3) > fD/2.)
@@ -282 +290 @@
 void MHexagon::PaintHexagon(Float_t inX, Float_t inY, Float_t inD)
 { 
-    const Int_t np = 6;
-
-    const Double_t dx[np+1] = { .5   , 0.    , -.5   , -.5   , 0.    ,  .5   , .5    };
-    const Double_t dy[np+1] = { .2886,  .5772,  .2886, -.2886, -.5772, -.2886, .2886 };
-
     //
     //  calculate the positions of the pixel corners
     //
-    Double_t x[np+1], y[np+1];
-    for (Int_t i=0; i<np+1; i++)
-    {
-        x[i] = inX + dx[i]*inD;
-        y[i] = inY + dy[i]*inD;
+    Double_t x[7], y[7];
+    for (Int_t i=0; i<7; i++)
+    {
+        x[i] = inX + fgDx[i%6]*inD;
+        y[i] = inY + fgDy[i%6]*inD;
     }
 
@@ -304 +307 @@
     //
     if (GetFillColor())
-        gPad->PaintFillArea(np, x, y);
+        gPad->PaintFillArea(6, x, y);
 
     if (GetLineStyle())
-        gPad->PaintPolyLine(np+1, x, y);
+        gPad->PaintPolyLine(7, x, y);
 }
 
@@ -367 +370 @@
 #endif
 }
+
+// ------------------------------------------------------------------------
+//
+// Small helper class to allow fast sorting of TVector2 by angle
+//
+class HVector2 : public TVector2
+{
+    Double_t fAngle;
+public:
+    HVector2() : TVector2()  { }
+    HVector2(const TVector2 &v) : TVector2(v)  { }
+    HVector2(Double_t x, Double_t y) : TVector2 (x, y) { }
+    void CalcAngle() { fAngle = Phi(); }
+    Bool_t IsSortable() const { return kTRUE; }
+    Int_t Compare(const TObject *obj) const { return ((HVector2*)obj)->fAngle>fAngle ? 1 : -1; }
+    Double_t GetAngle() const { return fAngle; }
+};
+
+// ------------------------------------------------------------------------
+//
+// Calculate the edge points of the intersection area of two hexagons.
+// The points are added as TVector2 to the TOrdCollection.
+// The user is responsible of delete the objects.
+//
+void MHexagon::GetIntersectionBorder(TOrdCollection &col, const MHexagon &hex) const
+{
+    Bool_t inside0[6], inside1[6];
+
+    HVector2 v0[6];
+    HVector2 v1[6];
+
+    Int_t cnt0=0;
+    Int_t cnt1=0;
+
+    // Calculate teh edges of each hexagon and whether this edge
+    // is inside the other hexgon or not
+    for (int i=0; i<6; i++)
+    {
+        const Double_t x0 = fX+fgDx[i]*fD;
+        const Double_t y0 = fY+fgDy[i]*fD;
+
+        const Double_t x1 = hex.fX+fgDx[i]*hex.fD;
+        const Double_t y1 = hex.fY+fgDy[i]*hex.fD;
+
+        v0[i].Set(x0, y0);
+        v1[i].Set(x1, y1);
+
+        inside0[i] = hex.DistanceToPrimitive(x0, y0) < 0;
+        inside1[i] = DistanceToPrimitive(x1, y1)     < 0;
+
+        if (inside0[i])
+        {
+            col.Add(new HVector2(v0[i]));
+            cnt0++;
+        }
+        if (inside1[i])
+        {
+            col.Add(new HVector2(v1[i]));
+            cnt1++;
+        }
+    }
+
+    if (cnt0==0 || cnt1==0)
+        return;
+
+    // No calculate which vorder lines intersect
+    Bool_t iscross0[6], iscross1[6];
+    for (int i=0; i<6; i++)
+    {
+        iscross0[i] = (inside0[i] && !inside0[(i+1)%6]) || (!inside0[i] && inside0[(i+1)%6]);
+        iscross1[i] = (inside1[i] && !inside1[(i+1)%6]) || (!inside1[i] && inside1[(i+1)%6]);
+    }
+
+    // Calculate the border points of our intersection area
+    for (int i=0; i<6; i++)
+    {
+        // Skip non intersecting lines
+        if (!iscross0[i])
+            continue;
+
+        for (int j=0; j<6; j++)
+        {
+            // Skip non intersecting lines
+            if (!iscross1[j])
+                continue;
+
+            const TVector2 p = MMath::GetIntersectionPoint(v0[i], v0[(i+1)%6], v1[j], v1[(j+1)%6]);
+            if (hex.DistanceToPrimitive(p.X(), p.Y())<1e-9)
+                col.Add(new HVector2(p));
+        }
+    }
+}
+
+// ------------------------------------------------------------------------
+//
+// Calculate the overlapping area of the two hexagons.
+//
+Double_t MHexagon::CalcOverlapArea(const MHexagon &cam) const
+{
+    TOrdCollection col;
+    col.SetOwner();
+
+    GetIntersectionBorder(col, cam);
+
+    // Check if there is an intersection to proceed with
+    const Int_t n = col.GetEntries();
+    if (n==0)
+        return 0;
+
+    // Calculate the center of gravity
+    TVector2 cog;
+
+    TIter Next(&col);
+    HVector2 *v=0;
+    while ((v=(HVector2*)Next()))
+        cog += *v;
+    cog /= n;
+
+    // Shift the figure to its center-og-gravity and
+    // calculate the angle of the connection line between the
+    // border points of our intersesction area and its cog
+    Next.Reset();
+    while ((v=(HVector2*)Next()))
+    {
+        *v -= cog;
+        v->CalcAngle();
+    }
+
+    // Sort these points by this angle
+    col.Sort();
+
+    // Now sum up the area of all the triangles between two
+    // following points and the cog.
+    Double_t A = 0;
+    for (int i=0; i<n; i++)
+    {
+        // Vectors from cog to two nearby border points
+        const HVector2 *v1 = (HVector2*)col.At(i);
+        const HVector2 *v2 = (HVector2*)col.At((i+1)%n);
+
+        // Angle between both vectors
+        const Double_t a = fmod(v1->GetAngle()-v2->GetAngle()+TMath::TwoPi(), TMath::TwoPi());
+
+        // Length of both vectors
+        const Double_t d1 = v1->Mod();
+        const Double_t d2 = v2->Mod();
+
+        A += d1*d2/2*sin(a);
+    }
+    return A;
+}

trunk/MagicSoft/Mars/mgui/MHexagon.h

@@ -27 +27 @@
 
 class MGeomPix;
+class TOrdCollection;
 
 class MHexagon : public TObject, public TAttLine, public TAttFill
 {
+private:
+    static const Double_t fgCos60;
+    static const Double_t fgSin60;
+
+    static const Double_t fgDx[6];   // X coordinate of the six edges
+    static const Double_t fgDy[6];   // Y coordinate of the six edges
+
 protected:
 
@@ -54 +62 @@
         return DistancetoPrimitive(px, py, 1);
     }
-    virtual Float_t DistanceToPrimitive(Float_t px, Float_t py);
+    virtual Float_t DistanceToPrimitive(Float_t px, Float_t py) const;
     virtual void  DrawHexagon(Float_t x, Float_t y, Float_t d);
 
@@ -70 +78 @@
     Float_t GetD() const { return fD; }
 
+    void GetIntersectionBorder(TOrdCollection &col, const MHexagon &hex) const;
+    Double_t CalcOverlapArea(const MHexagon &cam) const;
+
     ClassDef(MHexagon, 1)    // A hexagon for MAGIC
 };