source: trunk/MagicSoft/Mars/mgeom/MGeomPix.cc@ 9350

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Harald Kornmayer 1/2001
!
!   Copyright: MAGIC Software Development, 2000-2008
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MGeomPix
//
// This container stores the geometry (position) information of
// a single pixel together with the information about next neighbors.
//
//
// Version 1:
// ----------
//  - first implementation
//
// Version 2:
// ----------
//  - added fA
//
// Version 3:
// ----------
//  - added fAidx
//
// Version 4:
// ----------
//  - added fUserBits
//
//
// FIXME: According to an agreement we have to change the name 'Id' to 'idx'
//
////////////////////////////////////////////////////////////////////////////
#include "MGeomPix.h"

#include <TMath.h>
#include <TVector2.h>

#include "MLog.h"
#include "MLogManip.h"

#include "MGeomCam.h"

ClassImp(MGeomPix);

using namespace std;

const Float_t MGeomPix::gsTan30 = TMath::Tan(30/kRad2Deg); // sqrt(3)/3
const Float_t MGeomPix::gsTan60 = TMath::Tan(60/kRad2Deg); // sqrt(3)

// --------------------------------------------------------------------------
//
// Initializes one pixel
//
MGeomPix::MGeomPix(Float_t x, Float_t y, Float_t r, UInt_t s, UInt_t a) : fUserBits(0)
{
    //  default constructor
    Set(x, y, r, s, a);
    SetNeighbors();
}

// --------------------------------------------------------------------------
//
// Return position as TVector2(fX, fY)
//
TVector2 MGeomPix::GetP() const
{
    return TVector2(fX, fY);
}

// --------------------------------------------------------------------------
//
// Initializes Next Neighbors.
//
// WARNING: This function is public, but it is not meant for user access.
// It should only be used from geometry classes (like MGeomCam)
//
void MGeomPix::SetNeighbors(Short_t i0, Short_t i1, Short_t i2,
                            Short_t i3, Short_t i4, Short_t i5)
{
    fNeighbors[0] = i0;
    fNeighbors[1] = i1;
    fNeighbors[2] = i2;
    fNeighbors[3] = i3;
    fNeighbors[4] = i4;
    fNeighbors[5] = i5;

    int i;
    for (i=0; i<6; i++)
        if (fNeighbors[i]<0)
            break;

    fNumNeighbors = i;

    fNumNeighbors<5 ? SETBIT(fUserBits, kIsInOutermostRing) : CLRBIT(fUserBits, kIsInOutermostRing);
}

// --------------------------------------------------------------------------
//
//  Set the kIsOuterRing flag if this pixel has a outermost pixel
//  as Next Neighbor and don't have the kIsOutermostRing flag itself.
//
void MGeomPix::CheckOuterRing(const MGeomCam &cam)
{
    if (IsInOutermostRing())
        return;

    CLRBIT(fUserBits, kIsInOuterRing);

    for (int i=0; i<fNumNeighbors; i++)
        if (cam[fNeighbors[i]].IsInOutermostRing())
        {
            SETBIT(fUserBits, kIsInOuterRing);
            return;
        }
}

// --------------------------------------------------------------------------
//
// Print the geometry information of one pixel.
//
void MGeomPix::Print(Option_t *opt) const
{ 
    //   information about a pixel
    *fLog << all << "MPixGeom:  x/y=" << fX << "/" << fY << "mm ";
    *fLog << "d= " << fD << "mm  A= " << fA << "mm² (";
    for (int i=0; i<fNumNeighbors; i++)
        *fLog << fNeighbors[i] << (i<fNumNeighbors-1?",":"");
    *fLog << ")" << endl;
}

// ------------------------------------------------------------------------
//
// Return distance of center to coordinate origin: hypot(fX,fY)
//
Float_t MGeomPix::GetDist() const
{
    return TMath::Hypot(fX, fY);
}

// ------------------------------------------------------------------------
//
// Return distance of center to center of pix: hypot(fX-pix.fX,fY-pix.fY)
//
Float_t MGeomPix::GetDist(const MGeomPix &pix) const
{
    return TMath::Hypot(fX-pix.fX, fY-pix.fY);
}

// ------------------------------------------------------------------------
//
// Return angle defined by the center and the center of pix:
//  atan2(fX-pix.fX,fY-pix.fY)
//
Float_t MGeomPix::GetAngle(const MGeomPix &pix) const
{
    return TMath::ATan2(fX - pix.GetX(), fY - pix.GetY());
}

// ------------------------------------------------------------------------
//
// compute the distance of a point (px,py) to the Hexagon center in
// MGeomPix coordinates. Return kTRUE if inside.
//
Bool_t MGeomPix::IsInside(Float_t px, Float_t py) const
{
    //
    //  compute the distance of the Point to the center of the Hexagon
    //
    const Double_t dx = px-fX;

    //
    // 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 kFALSE;

    const Double_t dy = py-fY;

    const static Double_t cos60 = TMath::Cos(60/kRad2Deg);
    const static Double_t sin60 = TMath::Sin(60/kRad2Deg);

    const Double_t dxc = dx*cos60;
    const Double_t dys = dy*sin60;

    if  (TMath::Abs(dxc + dys)>fD/2)
        return kFALSE;

    if (TMath::Abs(dxc - dys)>fD/2)
        return kFALSE;

    return kTRUE;
}

// ------------------------------------------------------------------------
//
// Return the direction of the pixel pix w.r.t. this pixel.
// Remark: This function assumes a simple geometry.
//
Int_t MGeomPix::GetDirection(const MGeomPix &pix) const
{
    const Double_t x1 = GetX();
    const Double_t y1 = GetY();

    const Double_t x2 = pix.GetX();
    const Double_t y2 = pix.GetY();

    if (x1<=x2 && y1<y2) return kRightTop;
    if (x1<=x2 && y1>y2) return kRightBottom;
    if (x1>=x2 && y1<y2) return kLeftTop;
    if (x1>=x2 && y1>y2) return kLeftBottom;
    if (x1<x2)           return kRight;
    if (x1>x2)           return kLeft;

    cout << -1 << endl;

    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 <TOrdCollection.h>
#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; 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;
}
*/