Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 9234) +++ /trunk/MagicSoft/Mars/Changelog (revision 9235) @@ -19,4 +19,44 @@ -*-*- END OF LINE -*-*- + 2009/01/21 Thomas Bretz + + * mbase/BaseLinkDef.h, mbase/Makefile: + - added MQuaternion + + * mbase/MQuaternion.[h,cc]: + - added + + * mbase/MArrayI.[h,cc]: + - added a few new functions (ReSort, Add, Find and AddUniq) + + * mbase/MMath.cc: + - replaced delete by delete[] where appropriate + + * mgeom/MGeomCam.h: + - added new virtual functions HitFrame and HitDetector + - included MQuaternion + + * mgeom/MGeomCamDwarf.[h,cc], mgeom/MGeomCamMagic.[h,cc]: + - added new functions HitFrame and HitDetector: + + * msim/MHPhotonEvent.cc: + - removed the reference to MCamera + + * mgeom/MGeomPix.[h,cc], mgui/MHexagon.[h,cc]: + - moved DistanceToPrimitive from MHexagon to MGeomPix + - moved the funtions to calculate overlapping areas + also, but they are still in a comment + + * mgui/MHexagon.[h,cc]: + - started implementing a rotation angle + - added data meber fPhi + - increased class version number + + * mhist/MHCamera.cc, mtools/MagicJam.cc: + - we don't have to use a MHexagon anymore caluclating + DistanceToPrimitive + + + 2009/01/21 Stefan Ruegamer @@ -31,4 +71,9 @@ - added Esc member function + * msim/MPhotonData.[h,cc], msim/MPhotonEvent.[h,cc], + msim/MSimAbsorption.[h,cc], msim/MSimPointingPos.[h,cc], + MHPhotonEvent.[h,cc], msim/Makefile, msim/SimIncl.h + msim/SimLinkDef.h: + - added @@ -95,4 +140,5 @@ * readcorsika.cc: - added + Index: /trunk/MagicSoft/Mars/mgeom/MGeomPix.cc =================================================================== --- /trunk/MagicSoft/Mars/mgeom/MGeomPix.cc (revision 9234) +++ /trunk/MagicSoft/Mars/mgeom/MGeomPix.cc (revision 9235) @@ -240,2 +240,223 @@ return -1; } + +// ------------------------------------------------------------------------ +// +// compute the distance of a point (px,py) to the Hexagon center in world +// coordinates. Return -1 if inside. +// +Float_t MGeomPix::DistanceToPrimitive(Float_t px, Float_t py) const +{ + //FIXME: Rotation phi missing! + + static Double_t fgCos60 = 0.5; //TMath::Cos(TMath::DegToRad()*60); + static Double_t fgSin60 = TMath::Sqrt(3.)/2; //TMath::Sin(TMath::DegToRad()*60); + + // + // compute the distance of the Point to the center of the Hexagon + // + //const Double_t dx = px-fX; + //const Double_t dy = py-fY; + + TVector2 v(px-fX, py-fY); + // FIXME: fPhi or -fPhi? + // v = v.Rotate(-fPhi); // FIXME: Replace with a precalculates sin/cos vector + + const Double_t dx = v.X(); + const Double_t dy = v.Y(); + + const Double_t disthex = TMath::Sqrt(dx*dx + dy*dy); + + // + // Now check if point is outside of hexagon; just check x coordinate + // in three coordinate systems: the default one, in which two sides of + // the hexagon are paralel to the y axis (see camera displays) and two + // more, rotated with respect to that one by +- 60 degrees. + // + + if (TMath::Abs(dx) > fD/2.) + return disthex; + + const Double_t dx2 = dx*fgCos60 + dy*fgSin60; + + if (TMath::Abs(dx2) > fD/2.) + return disthex; + + const Double_t dx3 = dx*fgCos60 - dy*fgSin60; + + if (TMath::Abs(dx3) > fD/2.) + return disthex; + + return -1; +} + +// ============================================================================== +/* +#include +#include "MMath.h" +#include "../mgui/MHexagon.h" // MHexagon::fgDx + +// ------------------------------------------------------------------------ +// +// 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 MGeomPix::GetIntersectionBorder(TOrdCollection &col, const MGeomPix &hex) const +{ +// if (fPhi!=0) +// { +// gLog << warn << "MGeomPix::GetIntersectionBorder: Only phi=0 supported yet." << endl; +// return; +// } + + 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+MHexagon::fgDx[i]*fD; + const Double_t y0 = fY+MHexagon::fgDy[i]*fD; + + const Double_t x1 = hex.fX+MHexagon::fgDx[i]*hex.fD; + const Double_t y1 = hex.fY+MHexagon::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 MGeomPix::CalcOverlapArea(const MGeomPix &cam) const +{ +// if (fPhi!=0) +// { +// gLog << warn << "MGeomPix::CalcOverlapArea: Only phi=0 supported yet." << endl; +// return -1; +// } + + 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/mgeom/MGeomPix.h =================================================================== --- /trunk/MagicSoft/Mars/mgeom/MGeomPix.h (revision 9234) +++ /trunk/MagicSoft/Mars/mgeom/MGeomPix.h (revision 9235) @@ -8,4 +8,5 @@ class MGeomCam; class TVector2; +class TOrdCollection; class MGeomPix : public MParContainer @@ -96,4 +97,8 @@ Int_t GetDirection(const MGeomPix &pix) const; + virtual Float_t DistanceToPrimitive(Float_t px, Float_t py) const; + //void GetIntersectionBorder(TOrdCollection &col, const MGeomPix &hex) const; + //Double_t CalcOverlapArea(const MGeomPix &cam) const; + ClassDef(MGeomPix, 4) // Geometry class describing the geometry of one pixel }; Index: /trunk/MagicSoft/Mars/mgui/MHexagon.cc =================================================================== --- /trunk/MagicSoft/Mars/mgui/MHexagon.cc (revision 9234) +++ /trunk/MagicSoft/Mars/mgui/MHexagon.cc (revision 9235) @@ -1,4 +1,4 @@ /* ======================================================================== *\ -! $Name: not supported by cvs2svn $:$Id: MHexagon.cc,v 1.30 2006-10-30 12:46:13 tbretz Exp $ +! $Name: not supported by cvs2svn $:$Id: MHexagon.cc,v 1.31 2009-01-21 14:34:48 tbretz Exp $ ! -------------------------------------------------------------------------- ! @@ -21,5 +21,5 @@ ! Author(s): Harald Kornmayer 1/2001 ! -! Copyright: MAGIC Software Development, 2000-2006 +! Copyright: Software Development, 2000-2009 ! ! @@ -29,4 +29,7 @@ // // MHexagon +// +// Class Version 2: +// - added fPhi // ////////////////////////////////////////////////////////////////////////////// @@ -41,4 +44,7 @@ #include // TOrdCollection +#include "MLog.h" +#include "MLogManip.h" + #include "MMath.h" #include "MGeomPix.h" // GetX @@ -66,6 +72,6 @@ // normal constructor for MHexagon // -MHexagon::MHexagon(Float_t x, Float_t y, Float_t d) -: TAttLine(1, 1, 1), TAttFill(0, 1001), fX(x), fY(y), fD(d) +MHexagon::MHexagon(Float_t x, Float_t y, Float_t d, Float_t phi) +: TAttLine(1, 1, 1), TAttFill(0, 1001), fX(x), fY(y), fD(d), fPhi(phi) { } @@ -76,9 +82,11 @@ // MHexagon::MHexagon(const MGeomPix &pix) -: TAttLine(1, 1, 1), TAttFill(0, 1001) +: TAttLine(1, 1, 1), TAttFill(0, 1001), fPhi(0) { fX = pix.GetX(); fY = pix.GetY(); fD = pix.GetD(); + + // fPhi = pix.GetPhi(); } @@ -89,7 +97,8 @@ MHexagon::MHexagon(const MHexagon &hexagon) : TObject(hexagon), TAttLine(hexagon), TAttFill(hexagon) { - fX = hexagon.fX; - fY = hexagon.fY; - fD = hexagon.fD; + fX = hexagon.fX; + fY = hexagon.fY; + fD = hexagon.fD; + fPhi = hexagon.fPhi; } @@ -109,7 +118,8 @@ TAttFill::Copy(hex); - hex.fX = fX; - hex.fY = fY; - hex.fD = fD; + hex.fX = fX; + hex.fY = fY; + hex.fD = fD; + hex.fPhi = fPhi; } @@ -125,4 +135,6 @@ Int_t MHexagon::DistancetoPrimitive(Int_t px, Int_t py, Float_t conv) { + //FIXME: Rotation phi missing! + // // compute the distance of the Point to the center of the Hexagon @@ -131,6 +143,13 @@ const Int_t pyhex = gPad->YtoAbsPixel(fY*conv); - const Double_t x = TMath::Abs(px-pxhex); - const Double_t y = TMath::Abs(py-pyhex); + //const Double_t x = TMath::Abs(px-pxhex); + //const Double_t y = TMath::Abs(py-pyhex); + + TVector2 v(TMath::Abs(px-pxhex), TMath::Abs(py-pyhex)); + // FIXME: fPhi or -fPhi? + v = v.Rotate(-fPhi); // FIXME: Replace with a precalculates sin/cos vector + + const Double_t x = TMath::Abs(v.X()); + const Double_t y = TMath::Abs(v.Y()); const Double_t disthex = TMath::Sqrt(x*x + y*y); @@ -162,5 +181,5 @@ return distborder < disthex ? (int)((disthex-distborder)/conv+1) : 0; } - +/* // ------------------------------------------------------------------------ // @@ -170,9 +189,19 @@ Float_t MHexagon::DistanceToPrimitive(Float_t px, Float_t py) const { + //FIXME: Rotation phi missing! + // // compute the distance of the Point to the center of the Hexagon // - const Double_t dx = px-fX; - const Double_t dy = py-fY; + //const Double_t dx = px-fX; + //const Double_t dy = py-fY; + + TVector2 v(px-fX, py-fY); + // FIXME: fPhi or -fPhi? + v = v.Rotate(-fPhi); // FIXME: Replace with a precalculates sin/cos vector + + const Double_t dx = v.X(); + const Double_t dy = v.Y(); + const Double_t disthex = TMath::Sqrt(dx*dx + dy*dy); @@ -199,5 +228,5 @@ return -1; } - +*/ /* Float_t MHexagon::DistanceToPrimitive(Float_t px, Float_t py) @@ -227,7 +256,7 @@ // Draw this ellipse with new coordinate // -void MHexagon::DrawHexagon(Float_t x, Float_t y, Float_t d) +void MHexagon::DrawHexagon(Float_t x, Float_t y, Float_t d, Float_t phi) { - MHexagon *newhexagon = new MHexagon(x, y, d); + MHexagon *newhexagon = new MHexagon(x, y, d, phi); TAttLine::Copy(*newhexagon); @@ -281,5 +310,5 @@ void MHexagon::Paint(Option_t *) { - PaintHexagon(fX, fY, fD); + PaintHexagon(fX, fY, fD, fPhi); } @@ -288,5 +317,5 @@ // draw this hexagon with the coordinates // -void MHexagon::PaintHexagon(Float_t inX, Float_t inY, Float_t inD) +void MHexagon::PaintHexagon(Float_t inX, Float_t inY, Float_t inD, Float_t phi) { // @@ -296,6 +325,10 @@ for (Int_t i=0; i<7; i++) { - x[i] = inX + fgDx[i%6]*inD; - y[i] = inY + fgDy[i%6]*inD; + TVector2 v(fgDx[i%6], fgDy[i%6]); + + v = v.Rotate(phi); // FIXME: Replace with a precalculates sin/cos vector + + x[i] = inX + v.X()*inD; + y[i] = inY + v.Y()*inD; } @@ -320,5 +353,5 @@ { cout << "MHexagon: "; - cout << "x=" << fX << "mm y=" << fY << "mm r=" << fD << "mm" << endl; + cout << "x=" << fX << "mm y=" << fY << "mm r=" << fD << "mm phi=" << TMath::RadToDeg() << "deg" << endl; cout << " Line:"; @@ -343,5 +376,5 @@ out << " MHexagon *"; - out << "hexagon = new MHexagon(" << fX << "," << fY << "," << fD << ");" << endl; + out << "hexagon = new MHexagon(" << fX << "," << fY << "," << fD << "," << fPhi << ");" << endl; SaveFillAttributes(out, "hexagon"); @@ -360,5 +393,5 @@ out << " MHexagon *"; - out << "hexagon = new MHexagon(" << fX << "," << fY << "," << fD << ");" << endl; + out << "hexagon = new MHexagon(" << fX << "," << fY << "," << fD << "," << fPhi << ");" << endl; SaveFillAttributes(out, "hexagon"); @@ -370,153 +403,2 @@ #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 9234) +++ /trunk/MagicSoft/Mars/mgui/MHexagon.h (revision 9235) @@ -35,4 +35,5 @@ static const Double_t fgSin60; +public: static const Double_t fgDx[6]; // X coordinate of the six edges static const Double_t fgDy[6]; // Y coordinate of the six edges @@ -44,8 +45,10 @@ Float_t fD; // diameter D or better distance between opposite sides + Float_t fPhi; // Rotation angle + public: MHexagon(); - MHexagon(Float_t x, Float_t y, Float_t d); + MHexagon(Float_t x, Float_t y, Float_t d, Float_t phi=0); MHexagon(const MGeomPix &pix); MHexagon(const MHexagon &hexagon); @@ -62,6 +65,6 @@ return DistancetoPrimitive(px, py, 1); } - virtual Float_t DistanceToPrimitive(Float_t px, Float_t py) const; - virtual void DrawHexagon(Float_t x, Float_t y, Float_t d); + + virtual void DrawHexagon(Float_t x, Float_t y, Float_t d, Float_t phi=0); //virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py); @@ -69,5 +72,5 @@ virtual void ls(const Option_t *Option="") const; virtual void Paint(Option_t *Option=""); - virtual void PaintHexagon(Float_t x, Float_t y, Float_t d); + virtual void PaintHexagon(Float_t x, Float_t y, Float_t d, Float_t phi=0); virtual void Print(Option_t *Option="") const; // *MENU* virtual void SavePrimitive(ostream &out, Option_t *); @@ -77,9 +80,7 @@ Float_t GetY() const { return fY; } Float_t GetD() const { return fD; } + Float_t GetPhi() const { return fPhi; } - void GetIntersectionBorder(TOrdCollection &col, const MHexagon &hex) const; - Double_t CalcOverlapArea(const MHexagon &cam) const; - - ClassDef(MHexagon, 1) // A hexagon for MAGIC + ClassDef(MHexagon, 2) // A hexagon for MAGIC }; Index: /trunk/MagicSoft/Mars/mhist/MHCamera.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 9234) +++ /trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 9235) @@ -1,4 +1,4 @@ /* ======================================================================== *\ -! $Name: not supported by cvs2svn $:$Id: MHCamera.cc,v 1.111 2008-12-21 18:09:49 tbretz Exp $ +! $Name: not supported by cvs2svn $:$Id: MHCamera.cc,v 1.112 2009-01-21 14:37:07 tbretz Exp $ ! -------------------------------------------------------------------------- ! @@ -307,6 +307,5 @@ for (Int_t idx=0; idx0) + if ((*fGeomCam)[idx].DistanceToPrimitive(x, y)>0) continue; Index: /trunk/MagicSoft/Mars/mtools/MagicJam.cc =================================================================== --- /trunk/MagicSoft/Mars/mtools/MagicJam.cc (revision 9234) +++ /trunk/MagicSoft/Mars/mtools/MagicJam.cc (revision 9235) @@ -158,6 +158,5 @@ for (i=0; i0) + if ((*fGeomCam)[i].DistanceToPrimitive(px, py)>0) continue;