Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 8177) +++ /trunk/MagicSoft/Mars/Changelog (revision 8178) @@ -19,4 +19,15 @@ -*-*- 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,4 +72,17 @@ - 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]: Index: /trunk/MagicSoft/Mars/mbase/MMath.cc =================================================================== --- /trunk/MagicSoft/Mars/mbase/MMath.cc (revision 8177) +++ /trunk/MagicSoft/Mars/mbase/MMath.cc (revision 8178) @@ -1,3 +1,5 @@ /* ======================================================================== *\ +! $Name: not supported by cvs2svn $:$Id: MMath.cc,v 1.30 2006-10-30 12:46:12 tbretz Exp $ +! -------------------------------------------------------------------------- ! ! * @@ -32,4 +34,8 @@ #include "MMath.h" +#ifndef ROOT_TVector2 +#include +#endif + #ifndef ROOT_TVector3 #include @@ -505,4 +511,56 @@ // -------------------------------------------------------------------------- // +// 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 Index: /trunk/MagicSoft/Mars/mbase/MMath.h =================================================================== --- /trunk/MagicSoft/Mars/mbase/MMath.h (revision 8177) +++ /trunk/MagicSoft/Mars/mbase/MMath.h (revision 8178) @@ -6,4 +6,5 @@ #endif +class TVector2; class TVector3; class TArrayD; @@ -47,4 +48,6 @@ 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) { Index: /trunk/MagicSoft/Mars/mgui/MHexagon.cc =================================================================== --- /trunk/MagicSoft/Mars/mgui/MHexagon.cc (revision 8177) +++ /trunk/MagicSoft/Mars/mgui/MHexagon.cc (revision 8178) @@ -1,3 +1,5 @@ /* ======================================================================== *\ +! $Name: not supported by cvs2svn $:$Id: MHexagon.cc,v 1.30 2006-10-30 12:46:13 tbretz Exp $ +! -------------------------------------------------------------------------- ! ! * @@ -19,5 +21,5 @@ ! Author(s): Harald Kornmayer 1/2001 ! -! Copyright: MAGIC Software Development, 2000-2002 +! Copyright: MAGIC Software Development, 2000-2006 ! ! @@ -25,7 +27,7 @@ ////////////////////////////////////////////////////////////////////////////// -// // -// MHexagon // -// // +// +// MHexagon +// ////////////////////////////////////////////////////////////////////////////// @@ -35,11 +37,20 @@ #include -#include // gPad - -#include "MGeomPix.h" // GetX +#include // TVector2 +#include // gPad +#include // 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,5 +168,5 @@ // 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,13 +187,10 @@ 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,17 +290,12 @@ 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; iPaintFillArea(np, x, y); + gPad->PaintFillArea(6, x, y); if (GetLineStyle()) - gPad->PaintPolyLine(np+1, x, y); + gPad->PaintPolyLine(7, x, y); } @@ -367,2 +370,153 @@ #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; iGetAngle()-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; +} Index: /trunk/MagicSoft/Mars/mgui/MHexagon.h =================================================================== --- /trunk/MagicSoft/Mars/mgui/MHexagon.h (revision 8177) +++ /trunk/MagicSoft/Mars/mgui/MHexagon.h (revision 8178) @@ -27,7 +27,15 @@ 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,5 +62,5 @@ 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,4 +78,7 @@ 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 };