<commit_before>// @(#)root/gl:$Id$
// Author:  <NAME>  01/04/2005
/*
  CSGLib - Software Library for Constructive Solid Geometry
  Copyright (C) 2003-2004  <NAME>

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Library General Public
  License as published by the Free Software Foundation; either
  version 2 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Library General Public License for more details.

  You should have received a copy of the GNU Library General Public
  License along with this library; if not, write to the Free
  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

  Please send remarks, questions and bug reports to <EMAIL>
*/

/*
  SOLID - Software Library for Interference Detection
  Copyright (C) 1997-1998  <NAME>

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Library General Public
  License as published by the Free Software Foundation; either
  version 2 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Library General Public License for more details.

  You should have received a copy of the GNU Library General Public
  License along with this library; if not, write to the Free
  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

  Please send remarks, questions and bug reports to <EMAIL>,
  or write to:
                  <NAME>
         Department of Mathematics and Computing Science
         Eindhoven University of Technology
         P.O. Box 513, 5600 MB Eindhoven, The Netherlands
*/

/*
   This file contains compressed version of CSGSolid and SOLID.
   All stuff is in ROOT::Csg namespace,
   I used ROOT's own typedefs and math functions, removed resulting triangulation
   (I use OpenGL's inner implementation), changed names from MT_xxx CSG_xxx
   to more natural etc.
   Interface to CSGSolid :
   ConvertToMesh(to convert TBuffer3D to mesh)
   BuildUnion
   BuildIntersection
   BuildDifference

   31.03.05 <NAME>.
*/

#include <ROOT/REveCsgOps.hxx>

#include <algorithm>
#include <vector>

#include "TBuffer3D.h"
#include "Rtypes.h"
#include "TMath.h"

#include <ROOT/REveGeoShape.hxx>
#include <ROOT/REveUtil.hxx>

// Includes for BuildFromCompositeShape - this could go somewhere else.
#include "TBuffer3DTypes.h"
#include "TGeoBoolNode.h"
#include "TGeoCompositeShape.h"
#include "TGeoMatrix.h"
#include "TPad.h"
#include "TVirtualViewer3D.h"


namespace ROOT {
namespace Experimental {
namespace EveCsg {

const Double_t epsilon = 1e-10;
const Double_t epsilon2 = 1e-20;
const Double_t infinity = 1e50;

/////////////////////////////////////////////////////////////////////////////

Int_t sign(Double_t x)
{
   return x < 0. ? -1 : x > 0. ? 1 : 0;
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_zero(Double_t x)
{
   return TMath::Abs(x) < epsilon;
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_zero2(Double_t x)
{
   return TMath::Abs(x) < epsilon2;
}

class Tuple2 {
protected:
   Double_t fCo[2];

public:
   Tuple2(){SetValue(0, 0);}
   Tuple2(const Double_t *vv){SetValue(vv);}
   Tuple2(Double_t xx, Double_t yy){SetValue(xx, yy);}

   Double_t       &operator [] (Int_t i){return fCo[i];}
   const Double_t &operator [] (Int_t i)const{return fCo[i];}

   Double_t       &X(){return fCo[0];}
   const Double_t &X()const{return fCo[0];}
   Double_t       &Y(){return fCo[1];}
   const Double_t &Y()const{return fCo[1];}
   Double_t       &U(){return fCo[0];}
   const Double_t &U()const{return fCo[0];}
   Double_t       &V(){return fCo[1];}
   const Double_t &V()const{return fCo[1];}

   Double_t       *GetValue(){return fCo;}
   const Double_t *GetValue()const{return fCo;}
   void            GetValue(Double_t *vv)const{vv[0] = fCo[0]; vv[1] = fCo[1];}

   void            SetValue(const Double_t *vv){fCo[0] = vv[0]; fCo[1] = vv[1];}
   void            SetValue(Double_t xx, Double_t yy){fCo[0] = xx; fCo[1] = yy;}
};

Bool_t operator == (const Tuple2 &t1, const Tuple2 &t2)
{
   return t1[0] == t2[0] && t1[1] == t2[1];
}

class TVector2 : public Tuple2 {
public:
   TVector2(){}
   TVector2(const Double_t *v) : Tuple2(v) {}
   TVector2(Double_t xx, Double_t yy) : Tuple2(xx, yy) {}

   TVector2 &operator += (const TVector2 &v);
   TVector2 &operator -= (const TVector2 &v);
   TVector2 &operator *= (Double_t s);
   TVector2 &operator /= (Double_t s);

   Double_t Dot(const TVector2 &v)const;
   Double_t Length2()const;
   Double_t Length()const;
   TVector2 Absolute()const;
   void     Normalize();
   TVector2 Normalized()const;
   void     Scale(Double_t x, Double_t y);
   TVector2 Scaled(Double_t x, Double_t y)const;
   Bool_t   FuzzyZero()const;
   Double_t Angle(const TVector2 &v)const;
   TVector2 Cross(const TVector2 &v)const;
   Double_t Triple(const TVector2 &v1, const TVector2 &v2)const;
};

/////////////////////////////////////////////////////////////////////////////
///

TVector2 &TVector2::operator+=(const TVector2 &vv)
{
   fCo[0] += vv[0]; fCo[1] += vv[1];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 &TVector2::operator-=(const TVector2 &vv)
{
   fCo[0] -= vv[0]; fCo[1] -= vv[1];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 &TVector2::operator *= (Double_t s)
{
   fCo[0] *= s; fCo[1] *= s; return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 &TVector2::operator /= (Double_t s)
{
   return *this *= 1. / s;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator + (const TVector2 &v1, const TVector2 &v2)
{
   return TVector2(v1[0] + v2[0], v1[1] + v2[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator - (const TVector2 &v1, const TVector2 &v2)
{
   return TVector2(v1[0] - v2[0], v1[1] - v2[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator - (const TVector2 &v)
{
   return TVector2(-v[0], -v[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator * (const TVector2 &v, Double_t s)
{
   return TVector2(v[0] * s, v[1] * s);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator * (Double_t s, const TVector2 &v)
{
   return v * s;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator / (const TVector2 & v, Double_t s)
{
   return v * (1.0 / s);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector2::Dot(const TVector2 &vv)const
{
   return fCo[0] * vv[0] + fCo[1] * vv[1];
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector2::Length2()const
{
   return Dot(*this);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector2::Length()const
{
   return TMath::Sqrt(Length2());
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 TVector2::Absolute()const
{
   return TVector2(TMath::Abs(fCo[0]), TMath::Abs(fCo[1]));
}

/////////////////////////////////////////////////////////////////////////////
///

Bool_t TVector2::FuzzyZero()const
{
   return fuzzy_zero2(Length2());
}

/////////////////////////////////////////////////////////////////////////////
///

void TVector2::Normalize()
{
   *this /= Length();
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 TVector2::Normalized()const
{
   return *this / Length();
}

/////////////////////////////////////////////////////////////////////////////
///

void TVector2::Scale(Double_t xx, Double_t yy)
{
   fCo[0] *= xx; fCo[1] *= yy;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 TVector2::Scaled(Double_t xx, Double_t yy)const
{
   return TVector2(fCo[0] * xx, fCo[1] * yy);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector2::Angle(const TVector2 &vv)const
{
   Double_t s = TMath::Sqrt(Length2() * vv.Length2());
   return TMath::ACos(Dot(vv) / s);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t  dot(const TVector2 &v1, const TVector2 &v2)
{
   return v1.Dot(v2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t length2(const TVector2 &v)
{
   return v.Length2();
}

/////////////////////////////////////////////////////////////////////////////

Double_t length(const TVector2 &v)
{
   return v.Length();
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_zero(const TVector2 &v)
{
   return v.FuzzyZero();
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_equal(const TVector2 &v1, const TVector2 &v2)
{
   return fuzzy_zero(v1 - v2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t Angle(const TVector2 &v1, const TVector2 &v2)
{
   return v1.Angle(v2);
}

class TPoint2 : public TVector2 {
public:
   TPoint2(){}
   TPoint2(const Double_t *v) : TVector2(v){}
   TPoint2(Double_t x, Double_t y) : TVector2(x, y) {}

   TPoint2 &operator += (const TVector2 &v);
   TPoint2 &operator -= (const TVector2 &v);
   TPoint2 &operator = (const TVector2 &v);

   Double_t Distance(const TPoint2 &p)const;
   Double_t Distance2(const TPoint2 &p)const;
   TPoint2   Lerp(const TPoint2 &p, Double_t t)const;
};

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 &TPoint2::operator += (const TVector2 &v)
{
   fCo[0] += v[0]; fCo[1] += v[1];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 &TPoint2::operator-=(const TVector2& v)
{
   fCo[0] -= v[0]; fCo[1] -= v[1];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 &TPoint2::operator=(const TVector2& v)
{
   fCo[0] = v[0]; fCo[1] = v[1];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPoint2::Distance(const TPoint2& p)const
{
   return (p - *this).Length();
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPoint2::Distance2(const TPoint2& p)const
{
   return (p - *this).Length2();
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 TPoint2::Lerp(const TPoint2 &p, Double_t t)const
{
   return TPoint2(fCo[0] + (p[0] - fCo[0]) * t,
                  fCo[1] + (p[1] - fCo[1]) * t);
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 operator + (const TPoint2 &p, const TVector2 &v)
{
   return TPoint2(p[0] + v[0], p[1] + v[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 operator - (const TPoint2 &p, const TVector2 &v)
{
   return TPoint2(p[0] - v[0], p[1] - v[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator - (const TPoint2 &p1, const TPoint2 &p2)
{
   return TVector2(p1[0] - p2[0], p1[1] - p2[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t distance(const TPoint2 &p1, const TPoint2 &p2)
{
   return p1.Distance(p2);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t distance2(const TPoint2 &p1, const TPoint2 &p2)
{
   return p1.Distance2(p2);
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 lerp(const TPoint2 &p1, const TPoint2 &p2, Double_t t)
{
   return p1.Lerp(p2, t);
}

class Tuple3 {
protected:
   Double_t fCo[3];
public:
   Tuple3(){SetValue(0, 0, 0);}
   Tuple3(const Double_t *v){SetValue(v);}
   Tuple3(Double_t xx, Double_t yy, Double_t zz){SetValue(xx, yy, zz);}

   Double_t       &operator [] (Int_t i){return fCo[i];}
   const Double_t &operator [] (Int_t i)const{return fCo[i];}

   Double_t       &X(){return fCo[0];}
   const Double_t &X()const{return fCo[0];}
   Double_t       &Y(){return fCo[1];}
   const Double_t &Y()const{return fCo[1];}
   Double_t       &Z(){return fCo[2];}
   const Double_t &Z()const{return fCo[2];}

   Double_t       *GetValue(){return fCo;}
   const Double_t *GetValue()const{ return fCo; }

   void GetValue(Double_t *v)const
   {
      v[0] = Double_t(fCo[0]), v[1] = Double_t(fCo[1]), v[2] = Double_t(fCo[2]);
   }

   void SetValue(const Double_t *v)
   {
      fCo[0] = Double_t(v[0]), fCo[1] = Double_t(v[1]), fCo[2] = Double_t(v[2]);
   }
   void SetValue(Double_t xx, Double_t yy, Double_t zz)
   {
      fCo[0] = xx; fCo[1] = yy; fCo[2] = zz;
   }
};

/////////////////////////////////////////////////////////////////////////////

Bool_t operator==(const Tuple3& t1, const Tuple3& t2)
{
   return t1[0] == t2[0] && t1[1] == t2[1] && t1[2] == t2[2];
}

class TVector3 : public Tuple3 {
public:
   TVector3(){}
   TVector3(const Double_t *v) : Tuple3(v){}
   TVector3(Double_t xx, Double_t yy, Double_t zz) : Tuple3(xx, yy, zz){}

   TVector3 &operator += (const TVector3& v);
   TVector3 &operator -= (const TVector3& v);
   TVector3 &operator *= (Double_t s);
   TVector3 &operator /= (Double_t s);

   Double_t Dot(const TVector3& v)const;
   Double_t Length2()const;
   Double_t Length()const;
   TVector3  Absolute()const;
   void     NoiseGate(Double_t threshold);
   void     Normalize();
   TVector3  Normalized()const;
   TVector3  SafeNormalized()const;
   void     Scale(Double_t x, Double_t y, Double_t z);
   TVector3  Scaled(Double_t x, Double_t y, Double_t z)const;
   Bool_t   FuzzyZero()const;
   Double_t Angle(const TVector3 &v)const;
   TVector3  Cross(const TVector3 &v)const;
   Double_t Triple(const TVector3 &v1, const TVector3 &v2)const;
   Int_t    ClosestAxis()const;

   static TVector3 Random();
};

/////////////////////////////////////////////////////////////////////////////
///

TVector3 &TVector3::operator += (const TVector3 &v)
{
   fCo[0] += v[0]; fCo[1] += v[1]; fCo[2] += v[2];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 &TVector3::operator -= (const TVector3 &v)
{
   fCo[0] -= v[0]; fCo[1] -= v[1]; fCo[2] -= v[2];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 &TVector3::operator *= (Double_t s)
{
   fCo[0] *= s; fCo[1] *= s; fCo[2] *= s;
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 &TVector3::operator /= (Double_t s)
{
   return *this *= Double_t(1.0) / s;
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector3::Dot(const TVector3 &v)const
{
   return fCo[0] * v[0] + fCo[1] * v[1] + fCo[2] * v[2];
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector3::Length2()const
{
   return Dot(*this);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector3::Length()const
{
   return TMath::Sqrt(Length2());
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TVector3::Absolute()const
{
   return TVector3(TMath::Abs(fCo[0]), TMath::Abs(fCo[1]), TMath::Abs(fCo[2]));
}

/////////////////////////////////////////////////////////////////////////////
///

Bool_t TVector3::FuzzyZero()const
{
   return fuzzy_zero(Length2());
}

/////////////////////////////////////////////////////////////////////////////
///

void TVector3::NoiseGate(Double_t threshold)
{
   if (Length2() < threshold) SetValue(0., 0., 0.);
}

/////////////////////////////////////////////////////////////////////////////
///

void TVector3::Normalize()
{
   *this /= Length();
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator * (const TVector3 &v, Double_t s)
{
   return TVector3(v[0] * s, v[1] * s, v[2] * s);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator / (const TVector3 &v, Double_t s)
{
   return v * (1. / s);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TVector3::Normalized()const
{
   return *this / Length();
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TVector3::SafeNormalized()const
{
   Double_t len = Length();
   return fuzzy_zero(len) ? TVector3(1., 0., 0.):*this / len;
}

/////////////////////////////////////////////////////////////////////////////
///

void TVector3::Scale(Double_t xx, Double_t yy, Double_t zz)
{
   fCo[0] *= xx; fCo[1] *= yy; fCo[2] *= zz;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TVector3::Scaled(Double_t xx, Double_t yy, Double_t zz)const
{
   return TVector3(fCo[0] * xx, fCo[1] * yy, fCo[2] * zz);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector3::Angle(const TVector3 &v)const
{
   Double_t s = TMath::Sqrt(Length2() * v.Length2());
   return TMath::ACos(Dot(v) / s);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TVector3::Cross(const TVector3 &v)const
{
   return TVector3(fCo[1] * v[2] - fCo[2] * v[1],
                   fCo[2] * v[0] - fCo[0] * v[2],
                   fCo[0] * v[1] - fCo[1] * v[0]);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector3::Triple(const TVector3 &v1, const TVector3 &v2)const
{
   return fCo[0] * (v1[1] * v2[2] - v1[2] * v2[1]) +
      fCo[1] * (v1[2] * v2[0] - v1[0] * v2[2]) +
      fCo[2] * (v1[0] * v2[1] - v1[1] * v2[0]);
}

/////////////////////////////////////////////////////////////////////////////
///

Int_t TVector3::ClosestAxis()const
{
   TVector3 a = Absolute();
   return a[0] < a[1] ? (a[1] < a[2] ? 2 : 1) : (a[0] < a[2] ? 2 : 0);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator + (const TVector3 &v1, const TVector3 &v2)
{
   return TVector3(v1[0] + v2[0], v1[1] + v2[1], v1[2] + v2[2]);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator - (const TVector3 &v1, const TVector3 &v2)
{
   return TVector3(v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2]);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator - (const TVector3 &v)
{
   return TVector3(-v[0], -v[1], -v[2]);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator * (Double_t s, const TVector3 &v)
{
   return v * s;
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator * (const TVector3 &v1, const TVector3 &v2)
{
   return TVector3(v1[0] * v2[0], v1[1] * v2[1], v1[2] * v2[2]);
}

/////////////////////////////////////////////////////////////////////////////

Double_t  dot(const TVector3 &v1, const TVector3 &v2)
{
   return v1.Dot(v2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t length2(const TVector3 &v)
{
   return v.Length2();
}

/////////////////////////////////////////////////////////////////////////////

Double_t length(const TVector3 &v)
{
   return v.Length();
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_zero(const TVector3 &v)
{
   return v.FuzzyZero();
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_equal(const TVector3 &v1, const TVector3 &v2)
{
   return fuzzy_zero(v1 - v2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t Angle(const TVector3 &v1, const TVector3 &v2)
{
   return v1.Angle(v2);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 cross(const TVector3 &v1, const TVector3 &v2)
{
   return v1.Cross(v2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t triple(const TVector3 &v1, const TVector3 &v2, const TVector3 &v3)
{
   return v1.Triple(v2, v3);
}

class TPoint3 : public TVector3 {
public:
   TPoint3(){}
   TPoint3(const Double_t *v) : TVector3(v) {}
   TPoint3(Double_t xx, Double_t yy, Double_t zz) : TVector3(xx, yy, zz) {}

   TPoint3 &operator += (const TVector3 &v);
   TPoint3 &operator -= (const TVector3 &v);
   TPoint3 &operator = (const TVector3 &v);

   Double_t Distance(const TPoint3 &p)const;
   Double_t Distance2(const TPoint3 &p)const;
   TPoint3   Lerp(const TPoint3 &p, Double_t t)const;
};

/////////////////////////////////////////////////////////////////////////////
///

TPoint3 &TPoint3::operator += (const TVector3 &v)
{
   fCo[0] += v[0]; fCo[1] += v[1]; fCo[2] += v[2];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint3 &TPoint3::operator-=(const TVector3 &v)
{
   fCo[0] -= v[0]; fCo[1] -= v[1]; fCo[2] -= v[2];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint3 &TPoint3::operator=(const TVector3 &v)
{
   fCo[0] = v[0]; fCo[1] = v[1]; fCo[2] = v[2];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPoint3::Distance(const TPoint3 &p)const
{
   return (p - *this).Length();
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPoint3::Distance2(const TPoint3 &p)const
{
   return (p - *this).Length2();
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint3 TPoint3::Lerp(const TPoint3 &p, Double_t t)const
{
   return TPoint3(fCo[0] + (p[0] - fCo[0]) * t,
                  fCo[1] + (p[1] - fCo[1]) * t,
                  fCo[2] + (p[2] - fCo[2]) * t);
}

/////////////////////////////////////////////////////////////////////////////

TPoint3 operator + (const TPoint3 &p, const TVector3 &v)
{
   return TPoint3(p[0] + v[0], p[1] + v[1], p[2] + v[2]);
}

/////////////////////////////////////////////////////////////////////////////

TPoint3 operator - (const TPoint3 &p, const TVector3 &v)
{
   return TPoint3(p[0] - v[0], p[1] - v[1], p[2] - v[2]);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator - (const TPoint3 &p1, const TPoint3 &p2)
{
   return TVector3(p1[0] - p2[0], p1[1] - p2[1], p1[2] - p2[2]);
}

/////////////////////////////////////////////////////////////////////////////

Double_t distance(const TPoint3 &p1, const TPoint3 &p2)
{
   return p1.Distance(p2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t distance2(const TPoint3 &p1, const TPoint3 &p2)
{
   return p1.Distance2(p2);
}

/////////////////////////////////////////////////////////////////////////////

TPoint3 lerp(const TPoint3 &p1, const TPoint3 &p2, Double_t t)
{
   return p1.Lerp(p2, t);
}

class Tuple4 {
protected:
   Double_t fCo[4];

public:
   Tuple4(){SetValue(0, 0, 0, 0);}
   Tuple4(const Double_t *v){SetValue(v);}
   Tuple4(Double_t xx, Double_t yy, Double_t zz, Double_t ww)
   {
      SetValue(xx, yy, zz, ww);
   }

   Double_t       &operator [] (Int_t i){return fCo[i];}
   const Double_t &operator [] (Int_t i)const{return fCo[i];}

   Double_t       &X(){return fCo[0];}
   const Double_t &X()const{return fCo[0];}
   Double_t       &Y(){return fCo[1];}
   const Double_t &Y()const{return fCo[1];}
   Double_t       &Z(){return fCo[2];}
   const Double_t &Z()const{return fCo[2];}
   Double_t       &W(){return fCo[3];}
   const Double_t &W()const{return fCo[3];}

   Double_t       *GetValue(){return fCo;}
   const Double_t *GetValue()const{return fCo;}

   void GetValue(Double_t *v)const
   {
      v[0] = fCo[0]; v[1] = fCo[1]; v[2] = fCo[2]; v[3] = fCo[3];
   }

   void SetValue(const Double_t *v)
   {
      fCo[0] = v[0]; fCo[1] = v[1]; fCo[2] = v[2]; fCo[3] = v[3];
   }
   void SetValue(Double_t xx, Double_t yy, Double_t zz, Double_t ww)
   {
      fCo[0] = xx; fCo[1] = yy; fCo[2] = zz; fCo[3] = ww;
   }
};

/////////////////////////////////////////////////////////////////////////////

Bool_t operator == (const Tuple4 &t1, const Tuple4 &t2)
{
   return t1[0] == t2[0] && t1[1] == t2[1] && t1[2] == t2[2] && t1[3] == t2[3];
}

class TMatrix3x3 {
protected:
   TVector3 fEl[3];

public:
   TMatrix3x3(){}
   TMatrix3x3(const Double_t *m){SetValue(m);}
   TMatrix3x3(const TVector3 &euler){SetEuler(euler);}
   TMatrix3x3(const TVector3 &euler, const TVector3 &s)
   {
      SetEuler(euler); Scale(s[0], s[1], s[2]);
   }
   TMatrix3x3(Double_t xx, Double_t xy, Double_t xz,
              Double_t yx, Double_t yy, Double_t yz,
              Double_t zx, Double_t zy, Double_t zz)
   {
      SetValue(xx, xy, xz, yx, yy, yz, zx, zy, zz);
   }

   TVector3 &operator       [] (Int_t i){return fEl[i];}
   const TVector3 &operator [] (Int_t i)const{return fEl[i];}

   void SetValue(const Double_t *m)
   {
      fEl[0][0] = *m++; fEl[1][0] = *m++; fEl[2][0] = *m++; m++;
      fEl[0][1] = *m++; fEl[1][1] = *m++; fEl[2][1] = *m++; m++;
      fEl[0][2] = *m++; fEl[1][2] = *m++; fEl[2][2] = *m;
   }
   void SetValue(Double_t xx, Double_t xy, Double_t xz,
                 Double_t yx, Double_t yy, Double_t yz,
                 Double_t zx, Double_t zy, Double_t zz)
   {
      fEl[0][0] = xx; fEl[0][1] = xy; fEl[0][2] = xz;
      fEl[1][0] = yx; fEl[1][1] = yy; fEl[1][2] = yz;
      fEl[2][0] = zx; fEl[2][1] = zy; fEl[2][2] = zz;
   }
   void SetEuler(const TVector3 &euler)
   {
      Double_t ci = TMath::Cos(euler[0]);
      Double_t cj = TMath::Cos(euler[1]);
      Double_t ch = TMath::Cos(euler[2]);
      Double_t si = TMath::Sin(euler[0]);
      Double_t sj = TMath::Sin(euler[1]);
      Double_t sh = TMath::Sin(euler[2]);
      Double_t cc = ci * ch;
      Double_t cs = ci * sh;
      Double_t sc = si * ch;
      Double_t ss = si * sh;
      SetValue(cj * ch, sj * sc - cs, sj * cc + ss,
               cj * sh, sj * ss + cc, sj * cs - sc,
               -sj, cj * si, cj * ci);
   }

   void Scale(Double_t x, Double_t y, Double_t z)
   {
      fEl[0][0] *= x; fEl[0][1] *= y; fEl[0][2] *= z;
      fEl[1][0] *= x; fEl[1][1] *= y; fEl[1][2] *= z;
      fEl[2][0] *= x; fEl[2][1] *= y; fEl[2][2] *= z;
   }
   TMatrix3x3 Scaled(Double_t x, Double_t y, Double_t z)const
   {
      return TMatrix3x3(fEl[0][0] * x, fEl[0][1] * y, fEl[0][2] * z,
                        fEl[1][0] * x, fEl[1][1] * y, fEl[1][2] * z,
                        fEl[2][0] * x, fEl[2][1] * y, fEl[2][2] * z);
   }

   void SetIdentity()
   {
      SetValue(1., 0., 0., 0., 1., 0., 0., 0., 1.);
   }
   void GetValue(Double_t *m)const
   {
      *m++ = fEl[0][0]; *m++ = fEl[1][0]; *m++ = fEl[2][0]; *m++ = 0.0;
      *m++ = fEl[0][1]; *m++ = fEl[1][1]; *m++ = fEl[2][1]; *m++ = 0.0;
      *m++ = fEl[0][2]; *m++ = fEl[1][2]; *m++ = fEl[2][2]; *m   = 0.0;
   }

   TMatrix3x3 &operator *= (const TMatrix3x3 &m);
   Double_t Tdot(Int_t c, const TVector3 &v)const
   {
      return fEl[0][c] * v[0] + fEl[1][c] * v[1] + fEl[2][c] * v[2];
   }
   Double_t Cofac(Int_t r1, Int_t c1, Int_t r2, Int_t c2)const
   {
      return fEl[r1][c1] * fEl[r2][c2] - fEl[r1][c2] * fEl[r2][c1];
   }

   Double_t  Determinant()const;
   TMatrix3x3 Adjoint()const;
   TMatrix3x3 Absolute()const;
   TMatrix3x3 Transposed()const;
   void      Transpose();
   TMatrix3x3 Inverse()const;
   void      Invert();
};

/////////////////////////////////////////////////////////////////////////////
///

TMatrix3x3 &TMatrix3x3::operator *= (const TMatrix3x3 &m)
{
   SetValue(m.Tdot(0, fEl[0]), m.Tdot(1, fEl[0]), m.Tdot(2, fEl[0]),
            m.Tdot(0, fEl[1]), m.Tdot(1, fEl[1]), m.Tdot(2, fEl[1]),
            m.Tdot(0, fEl[2]), m.Tdot(1, fEl[2]), m.Tdot(2, fEl[2]));
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TMatrix3x3::Determinant()const
{
   return triple((*this)[0], (*this)[1], (*this)[2]);
}

/////////////////////////////////////////////////////////////////////////////
///

TMatrix3x3 TMatrix3x3::Absolute()const
{
   return TMatrix3x3(TMath::Abs(fEl[0][0]), TMath::Abs(fEl[0][1]), TMath::Abs(fEl[0][2]),
                     TMath::Abs(fEl[1][0]), TMath::Abs(fEl[1][1]), TMath::Abs(fEl[1][2]),
                     TMath::Abs(fEl[2][0]), TMath::Abs(fEl[2][1]), TMath::Abs(fEl[2][2]));
}

/////////////////////////////////////////////////////////////////////////////
///

TMatrix3x3 TMatrix3x3::Transposed()const
{
   return TMatrix3x3(fEl[0][0], fEl[1][0], fEl[2][0],
                     fEl[0][1], fEl[1][1], fEl[2][1],
                     fEl[0][2], fEl[1][2], fEl[2][2]);
}

/////////////////////////////////////////////////////////////////////////////
///

void TMatrix3x3::Transpose()
{
   *this = Transposed();
}

/////////////////////////////////////////////////////////////////////////////
///

TMatrix3x3 TMatrix3x3::Adjoint()const
{
   return TMatrix3x3(Cofac(1, 1, 2, 2), Cofac(0, 2, 2, 1), Cofac(0, 1, 1, 2),
                     Cofac(1, 2, 2, 0), Cofac(0, 0, 2, 2), Cofac(0, 2, 1, 0),
                     Cofac(1, 0, 2, 1), Cofac(0, 1, 2, 0), Cofac(0, 0, 1, 1));
}

/////////////////////////////////////////////////////////////////////////////
///

TMatrix3x3 TMatrix3x3::Inverse()const
{
   TVector3 co(Cofac(1, 1, 2, 2), Cofac(1, 2, 2, 0), Cofac(1, 0, 2, 1));
   Double_t det = dot((*this)[0], co);
   Double_t s = 1. / det;
   return TMatrix3x3(co[0] * s, Cofac(0, 2, 2, 1) * s, Cofac(0, 1, 1, 2) * s,
                     co[1] * s, Cofac(0, 0, 2, 2) * s, Cofac(0, 2, 1, 0) * s,
                     co[2] * s, Cofac(0, 1, 2, 0) * s, Cofac(0, 0, 1, 1) * s);
}

/////////////////////////////////////////////////////////////////////////////
///

void TMatrix3x3::Invert()
{
   *this = Inverse();
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator * (const TMatrix3x3& m, const TVector3& v)
{
   return TVector3(dot(m[0], v), dot(m[1], v), dot(m[2], v));
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator * (const TVector3& v, const TMatrix3x3& m)
{
   return TVector3(m.Tdot(0, v), m.Tdot(1, v), m.Tdot(2, v));
}

/////////////////////////////////////////////////////////////////////////////

TMatrix3x3 operator * (const TMatrix3x3 &m1, const TMatrix3x3 &m2)
{
   return TMatrix3x3(m2.Tdot(0, m1[0]), m2.Tdot(1, m1[0]), m2.Tdot(2, m1[0]),
                     m2.Tdot(0, m1[1]), m2.Tdot(1, m1[1]), m2.Tdot(2, m1[1]),
                     m2.Tdot(0, m1[2]), m2.Tdot(1, m1[2]), m2.Tdot(2, m1[2]));
}

/////////////////////////////////////////////////////////////////////////////

TMatrix3x3 mmult_transpose_left(const TMatrix3x3 &m1, const TMatrix3x3 &m2)
{
   return TMatrix3x3(m1[0][0] * m2[0][0] + m1[1][0] * m2[1][0] + m1[2][0] * m2[2][0],
                     m1[0][0] * m2[0][1] + m1[1][0] * m2[1][1] + m1[2][0] * m2[2][1],
                     m1[0][0] * m2[0][2] + m1[1][0] * m2[1][2] + m1[2][0] * m2[2][2],
                     m1[0][1] * m2[0][0] + m1[1][1] * m2[1][0] + m1[2][1] * m2[2][0],
                     m1[0][1] * m2[0][1] + m1[1][1] * m2[1][1] + m1[2][1] * m2[2][1],
                     m1[0][1] * m2[0][2] + m1[1][1] * m2[1][2] + m1[2][1] * m2[2][2],
                     m1[0][2] * m2[0][0] + m1[1][2] * m2[1][0] + m1[2][2] * m2[2][0],
                     m1[0][2] * m2[0][1] + m1[1][2] * m2[1][1] + m1[2][2] * m2[2][1],
                     m1[0][2] * m2[0][2] + m1[1][2] * m2[1][2] + m1[2][2] * m2[2][2]);
}

/////////////////////////////////////////////////////////////////////////////

TMatrix3x3 mmult_transpose_right(const TMatrix3x3& m1, const TMatrix3x3& m2)
{
   return TMatrix3x3(m1[0].Dot(m2[0]), m1[0].Dot(m2[1]), m1[0].Dot(m2[2]),
                     m1[1].Dot(m2[0]), m1[1].Dot(m2[1]), m1[1].Dot(m2[2]),
                     m1[2].Dot(m2[0]), m1[2].Dot(m2[1]), m1[2].Dot(m2[2]));
}

class TLine3 {
private :
   Bool_t   fBounds[2];
   Double_t fParams[2];
   TPoint3   fOrigin;
   TVector3  fDir;

public :
   TLine3();
   TLine3(const TPoint3 &p1, const TPoint3 &p2);
   TLine3(const TPoint3 &p1, const TVector3 &v);
   TLine3(const TPoint3 &p1, const TVector3 &v, Bool_t bound1, Bool_t bound2);

   static TLine3   InfiniteRay(const TPoint3 &p1, const TVector3 &v);
   const TVector3 &Direction()const {return fDir;}
   const TPoint3  &Origin()const{ return fOrigin;}

   Bool_t Bounds(Int_t i)const
   {
      return (i == 0 ? fBounds[0] : fBounds[1]);
   }
   Bool_t &Bounds(Int_t i)
   {
      return (i == 0 ? fBounds[0] : fBounds[1]);
   }
   const Double_t &Param(Int_t i)const
   {
      return (i == 0 ? fParams[0] : fParams[1]);
   }
   Double_t &Param(Int_t i)
   {
      return (i == 0 ? fParams[0] : fParams[1]);
   }
   TVector3 UnboundSmallestVector(const TPoint3 &point)const
   {
      TVector3 diff(fOrigin - point);
      return diff - fDir * diff.Dot(fDir);
   }
   Double_t UnboundClosestParameter(const TPoint3 &point)const
   {
      TVector3 diff(fOrigin-point);
      return diff.Dot(fDir);
   }
   Double_t UnboundDistance(const TPoint3& point)const
   {
      return UnboundSmallestVector(point).Length();
   }
   Bool_t IsParameterOnLine(const Double_t &t) const
   {
      return ((fParams[0] - epsilon < t) || (!fBounds[0])) && ((fParams[1] > t + epsilon) || (!fBounds[1]));
   }
};

/////////////////////////////////////////////////////////////////////////////
///

TLine3::TLine3() : fOrigin(0,0,0), fDir(1,0,0)
{
   fBounds[0] = kFALSE; fBounds[1] = kFALSE;
   fParams[0] = 0; fParams[1] = 1;
}

/////////////////////////////////////////////////////////////////////////////
///

TLine3::TLine3(const TPoint3 &p1, const TPoint3 &p2) : fOrigin(p1), fDir(p2-p1)
{
   fBounds[0] = kTRUE; fBounds[1] = kTRUE;
   fParams[0] = 0; fParams[1] = 1;
}

/////////////////////////////////////////////////////////////////////////////
///

TLine3::TLine3(const TPoint3 &p1, const TVector3 &v): fOrigin(p1), fDir(v)
{
   fBounds[0] = kFALSE; fBounds[1] = kFALSE;
   fParams[0] = 0; fParams[1] = 1;
}

/////////////////////////////////////////////////////////////////////////////
///

TLine3::TLine3(const TPoint3 &p1, const TVector3 &v, Bool_t bound1, Bool_t bound2)
   : fOrigin(p1), fDir(v)
{
   fBounds[0] = bound1; fBounds[1] = bound2;
   fParams[0] = 0; fParams[1] = 1;
}

class TPlane3 : public Tuple4 {
public :
   TPlane3(const TVector3 &a, const TVector3 &b, const TVector3 &c);
   TPlane3(const TVector3 &n, const TVector3 &p);
   TPlane3();
   TPlane3(const TPlane3 & p):Tuple4(p){}

   TVector3  Normal()const;
   Double_t  Scalar()const;
   void     Invert();
   Double_t SignedDistance(const TVector3 &)const;

   TPlane3 &operator = (const TPlane3 & rhs);
};

/////////////////////////////////////////////////////////////////////////////
///

TPlane3::TPlane3(const TVector3 &a, const TVector3 &b, const TVector3 &c)
{
   TVector3 l1 = b-a;
   TVector3 l2 = c-b;
   TVector3 n = l1.Cross(l2);
   n = n.SafeNormalized();
   Double_t d = n.Dot(a);
   fCo[0] = n.X(); fCo[1] = n.Y(); fCo[2] = n.Z(); fCo[3] = -d;
}

/////////////////////////////////////////////////////////////////////////////
///

TPlane3::TPlane3(const TVector3 &n, const TVector3 &p)
{
   TVector3 mn = n.SafeNormalized();
   Double_t md = mn.Dot(p);
   fCo[0] = mn.X(); fCo[1] = mn.Y(); fCo[2] = mn.Z(); fCo[3] = -md;
}

/////////////////////////////////////////////////////////////////////////////
///

TPlane3::TPlane3() : Tuple4()
{
   fCo[0] = 1.; fCo[1] = 0.;
   fCo[2] = 0.; fCo[3] = 0.;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TPlane3::Normal()const
{
   return TVector3(fCo[0],fCo[1],fCo[2]);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPlane3::Scalar()const
{
   return fCo[3];
}

/////////////////////////////////////////////////////////////////////////////
///

void TPlane3::Invert()
{
   fCo[0] = -fCo[0]; fCo[1] = -fCo[1]; fCo[2] = -fCo[2]; fCo[3] = -fCo[3];
}

/////////////////////////////////////////////////////////////////////////////
///

TPlane3 &TPlane3::operator = (const TPlane3 &rhs)
{
   fCo[0] = rhs.fCo[0]; fCo[1] = rhs.fCo[1]; fCo[2] = rhs.fCo[2]; fCo[3] = rhs.fCo[3];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPlane3::SignedDistance(const TVector3 &v)const
{
   return Normal().Dot(v) + fCo[3];
}

class TBBox {
   friend Bool_t intersect(const TBBox &a, const TBBox &b);

private:
   TPoint3  fCenter;
   TVector3 fExtent;

public:
   TBBox(){}
   TBBox(const TPoint3 &mini, const TPoint3 &maxi){SetValue(mini,maxi);}

   const TPoint3   &Center()const{return fCenter;}
   const TVector3  &Extent()const{return fExtent;}
   TPoint3         &Center(){return fCenter;}
   TVector3        &Extent(){return fExtent;}

   void SetValue(const TPoint3 &mini,const TPoint3 &maxi)
   {
      fExtent = (maxi - mini) / 2;
      fCenter = mini + fExtent;
   }
   void Enclose(const TBBox &a, const TBBox &b)
   {
      TPoint3 lower(
                    TMath::Min(a.Lower(0), b.Lower(0)),
                    TMath::Min(a.Lower(1), b.Lower(1)),
                    TMath::Min(a.Lower(2), b.Lower(2))
                    );
      TPoint3 upper(
                    TMath::Max(a.Upper(0), b.Upper(0)),
                    TMath::Max(a.Upper(1), b.Upper(1)),
                    TMath::Max(a.Upper(2), b.Upper(2))
                    );
      SetValue(lower, upper);
   }

   void SetEmpty()
   {
      fCenter.SetValue(0., 0., 0.);
      fExtent.SetValue(-infinity, -infinity, -infinity);
   }
   void Include(const TPoint3 &p)
   {
      TPoint3 lower(
                    TMath::Min(Lower(0), p[0]),
                    TMath::Min(Lower(1), p[1]),
                    TMath::Min(Lower(2), p[2])
                    );
      TPoint3 upper(
                    TMath::Max(Upper(0), p[0]),
                    TMath::Max(Upper(1), p[1]),
                    TMath::Max(Upper(2), p[2])
                    );
      SetValue(lower, upper);
   }

   void Include(const TBBox &b)
   {
      Enclose(*this, b);
   }
   Double_t Lower(Int_t i)const
   {
      return fCenter[i] - fExtent[i];
   }
   Double_t Upper(Int_t i)const
   {
      return fCenter[i] + fExtent[i];
   }
   TPoint3 Lower()const
   {
      return fCenter - fExtent;
   }
   TPoint3 Upper()const
   {
      return fCenter + fExtent;
   }
   Double_t Size()const
   {
      return TMath::Max(TMath::Max(fExtent[0], fExtent[1]), fExtent[2]);
   }
   Int_t LongestAxis()const
   {
      return fExtent.ClosestAxis();
   }
   Bool_t IntersectXRay(const TPoint3 &xBase)const
   {
      if (xBase[0] <= Upper(0)) {
         if (xBase[1] <= Upper(1) && xBase[1] >= Lower(1)) {
            if (xBase[2] <= Upper(2) && xBase[2] >= Lower(2)) {return kTRUE;}
         }
      }
      return kFALSE;
   }
};

/////////////////////////////////////////////////////////////////////////////

Bool_t intersect(const TBBox &a, const TBBox &b)
{
   return TMath::Abs(a.fCenter[0] - b.fCenter[0]) <= a.fExtent[0] + b.fExtent[0] &&
      TMath::Abs(a.fCenter[1] - b.fCenter[1]) <= a.fExtent[1] + b.fExtent[1] &&
      TMath::Abs(a.fCenter[2] - b.fCenter[2]) <= a.fExtent[2] + b.fExtent[2];
}

class TBBoxNode {
public:
   enum ETagType {kLeaf, kInternal};
   TBBox     fBBox;
   ETagType fTag;
};

class TBBoxLeaf : public TBBoxNode {
public:
   Int_t fPolyIndex;

   TBBoxLeaf() : fPolyIndex(0) {}
   TBBoxLeaf(Int_t polyIndex, const TBBox &bbox) : fPolyIndex(polyIndex)
   {
      fBBox = bbox;
      fTag = kLeaf;
   }
};

typedef TBBoxLeaf *LeafPtr_t;
typedef TBBoxNode *NodePtr_t;

class TBBoxInternal : public TBBoxNode {
public:
   NodePtr_t fLeftSon;
   NodePtr_t fRightSon;
   TBBoxInternal() : fLeftSon(0) ,fRightSon(0) {}
   TBBoxInternal(Int_t n, LeafPtr_t leafIt);
};

typedef TBBoxInternal* InternalPtr_t;

class TBBoxTree {
private:
   Int_t         fBranch;
   LeafPtr_t     fLeaves;
   InternalPtr_t fInternals;
   Int_t         fNumLeaves;

public :
   TBBoxTree() : fBranch(0), fLeaves(0), fInternals(0), fNumLeaves(0) {}
   NodePtr_t RootNode()const{return fInternals;}
   ~TBBoxTree()
   {
      delete[] fLeaves;
      delete[] fInternals;
   }
   void BuildTree(LeafPtr_t leaves, Int_t numLeaves);

private :
   void RecursiveTreeBuild(Int_t n, LeafPtr_t leafIt);
};

/////////////////////////////////////////////////////////////////////////////
///

TBBoxInternal::TBBoxInternal(Int_t n, LeafPtr_t leafIt) :
   fLeftSon(0) ,fRightSon(0)
{
   fTag = kInternal;
   fBBox.SetEmpty();
   for (Int_t i=0;i<n;i++)
      fBBox.Include(leafIt[i].fBBox);
}

/////////////////////////////////////////////////////////////////////////////
///

void TBBoxTree::BuildTree(LeafPtr_t leaves, Int_t numLeaves)
{
   fBranch = 0;
   fLeaves = leaves;
   fNumLeaves = numLeaves;
   fInternals = new TBBoxInternal[numLeaves];
   RecursiveTreeBuild(fNumLeaves,fLeaves);
}

/////////////////////////////////////////////////////////////////////////////
///

void TBBoxTree::RecursiveTreeBuild(Int_t n, LeafPtr_t leafIt)
{
   fInternals[fBranch] = TBBoxInternal(n,leafIt);
   TBBoxInternal &aBBox = fInternals[fBranch];
   fBranch++;

   Int_t axis = aBBox.fBBox.LongestAxis();
   Int_t i = 0, mid = n;

   while (i < mid) {
      if (leafIt[i].fBBox.Center()[axis] < aBBox.fBBox.Center()[axis]) {
         ++i;
      } else {
         --mid;
         std::swap(leafIt[i], leafIt[mid]);
      }
   }

   if (mid == 0 || mid == n) {
      mid = n / 2;
   }
   if (mid >= 2) {
      aBBox.fRightSon = fInternals + fBranch;
      RecursiveTreeBuild(mid,leafIt);
   } else {
      aBBox.fRightSon = leafIt;
   }
   if (n - mid >= 2) {
      aBBox.fLeftSon = fInternals + fBranch;
      RecursiveTreeBuild(n - mid, leafIt + mid);
   } else {
      aBBox.fLeftSon = leafIt + mid;
   }
}

class TBlenderVProp {
private:
   Int_t fVertexIndex;

public:
   TBlenderVProp(Int_t vIndex) : fVertexIndex(vIndex){}
   TBlenderVProp(Int_t vIndex, const TBlenderVProp &,
                 const TBlenderVProp &, const Double_t &)
   {
      fVertexIndex = vIndex;
   }
   TBlenderVProp() : fVertexIndex(-1){}
   operator Int_t()const
   {
      return fVertexIndex;
   }
   TBlenderVProp &operator = (Int_t i)
   {
      fVertexIndex = i; return *this;
   }
};

template <class TMesh>
class TPolygonGeometry {
public:
   typedef typename TMesh::Polygon TPolygon;

private:
   const TMesh    &fMesh;
   const TPolygon &fPoly;
public:
   TPolygonGeometry(const TMesh &mesh, Int_t pIndex)
      : fMesh(mesh), fPoly(mesh.Polys()[pIndex])
   {}
   TPolygonGeometry(const TMesh &mesh, const TPolygon &poly)
      : fMesh(mesh), fPoly(poly)
   {}
   const TPoint3 &operator [] (Int_t i)const
   {
      return fMesh.Verts()[fPoly[i]].Pos();
   }
   Int_t Size()const
   {
      return fPoly.Size();
   }

private:
   TPolygonGeometry(const TPolygonGeometry &);
   TPolygonGeometry& operator = (TPolygonGeometry &);
};

template <class TPolygon, class TVertex>
class TMesh : public TBaseMesh {
public:
   typedef std::vector<TVertex> VLIST;
   typedef std::vector<TPolygon> PLIST;
   typedef TPolygon Polygon;
   typedef TVertex Vertex;
   typedef TPolygonGeometry<TMesh> TGBinder;

private:
   VLIST fVerts;
   PLIST fPolys;

public:
   VLIST       &Verts(){return fVerts;}
   const VLIST &Verts()const{return fVerts;}
   PLIST       &Polys(){return fPolys;}
   const PLIST &Polys()const{return fPolys;}

   //TBaseMesh's final-overriders
   Int_t          NumberOfPolys()const{return fPolys.size();}
   Int_t          NumberOfVertices()const{return fVerts.size();}
   Int_t          SizeOfPoly(Int_t polyIndex)const{return fPolys[polyIndex].Size();}
   const Double_t *GetVertex(Int_t vertexNum)const{return fVerts[vertexNum].GetValue();}

   Int_t GetVertexIndex(Int_t polyNum, Int_t vertexNum)const
   {
      return fPolys[polyNum][vertexNum];
   }
};

const Int_t cofacTable[3][2] = {{1,2}, {0,2}, {0,1}};
<commit_msg>eve7: identify methods which are overwritte in TBaseMatrix<commit_after>// @(#)root/gl:$Id$
// Author:  <NAME>  01/04/2005
/*
  CSGLib - Software Library for Constructive Solid Geometry
  Copyright (C) 2003-2004  <NAME>

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Library General Public
  License as published by the Free Software Foundation; either
  version 2 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Library General Public License for more details.

  You should have received a copy of the GNU Library General Public
  License along with this library; if not, write to the Free
  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

  Please send remarks, questions and bug reports to <EMAIL>
*/

/*
  SOLID - Software Library for Interference Detection
  Copyright (C) 1997-1998  <NAME>

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Library General Public
  License as published by the Free Software Foundation; either
  version 2 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Library General Public License for more details.

  You should have received a copy of the GNU Library General Public
  License along with this library; if not, write to the Free
  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

  Please send remarks, questions and bug reports to <EMAIL>,
  or write to:
                  <NAME>
         Department of Mathematics and Computing Science
         Eindhoven University of Technology
         P.O. Box 513, 5600 MB Eindhoven, The Netherlands
*/

/*
   This file contains compressed version of CSGSolid and SOLID.
   All stuff is in ROOT::Csg namespace,
   I used ROOT's own typedefs and math functions, removed resulting triangulation
   (I use OpenGL's inner implementation), changed names from MT_xxx CSG_xxx
   to more natural etc.
   Interface to CSGSolid :
   ConvertToMesh(to convert TBuffer3D to mesh)
   BuildUnion
   BuildIntersection
   BuildDifference

   31.03.05 <NAME>.
*/

#include <ROOT/REveCsgOps.hxx>

#include <algorithm>
#include <vector>

#include "TBuffer3D.h"
#include "Rtypes.h"
#include "TMath.h"

#include <ROOT/REveGeoShape.hxx>
#include <ROOT/REveUtil.hxx>

// Includes for BuildFromCompositeShape - this could go somewhere else.
#include "TBuffer3DTypes.h"
#include "TGeoBoolNode.h"
#include "TGeoCompositeShape.h"
#include "TGeoMatrix.h"
#include "TPad.h"
#include "TVirtualViewer3D.h"


namespace ROOT {
namespace Experimental {
namespace EveCsg {

const Double_t epsilon = 1e-10;
const Double_t epsilon2 = 1e-20;
const Double_t infinity = 1e50;

/////////////////////////////////////////////////////////////////////////////

Int_t sign(Double_t x)
{
   return x < 0. ? -1 : x > 0. ? 1 : 0;
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_zero(Double_t x)
{
   return TMath::Abs(x) < epsilon;
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_zero2(Double_t x)
{
   return TMath::Abs(x) < epsilon2;
}

class Tuple2 {
protected:
   Double_t fCo[2];

public:
   Tuple2(){SetValue(0, 0);}
   Tuple2(const Double_t *vv){SetValue(vv);}
   Tuple2(Double_t xx, Double_t yy){SetValue(xx, yy);}

   Double_t       &operator [] (Int_t i){return fCo[i];}
   const Double_t &operator [] (Int_t i)const{return fCo[i];}

   Double_t       &X(){return fCo[0];}
   const Double_t &X()const{return fCo[0];}
   Double_t       &Y(){return fCo[1];}
   const Double_t &Y()const{return fCo[1];}
   Double_t       &U(){return fCo[0];}
   const Double_t &U()const{return fCo[0];}
   Double_t       &V(){return fCo[1];}
   const Double_t &V()const{return fCo[1];}

   Double_t       *GetValue(){return fCo;}
   const Double_t *GetValue()const{return fCo;}
   void            GetValue(Double_t *vv)const{vv[0] = fCo[0]; vv[1] = fCo[1];}

   void            SetValue(const Double_t *vv){fCo[0] = vv[0]; fCo[1] = vv[1];}
   void            SetValue(Double_t xx, Double_t yy){fCo[0] = xx; fCo[1] = yy;}
};

Bool_t operator == (const Tuple2 &t1, const Tuple2 &t2)
{
   return t1[0] == t2[0] && t1[1] == t2[1];
}

class TVector2 : public Tuple2 {
public:
   TVector2(){}
   TVector2(const Double_t *v) : Tuple2(v) {}
   TVector2(Double_t xx, Double_t yy) : Tuple2(xx, yy) {}

   TVector2 &operator += (const TVector2 &v);
   TVector2 &operator -= (const TVector2 &v);
   TVector2 &operator *= (Double_t s);
   TVector2 &operator /= (Double_t s);

   Double_t Dot(const TVector2 &v)const;
   Double_t Length2()const;
   Double_t Length()const;
   TVector2 Absolute()const;
   void     Normalize();
   TVector2 Normalized()const;
   void     Scale(Double_t x, Double_t y);
   TVector2 Scaled(Double_t x, Double_t y)const;
   Bool_t   FuzzyZero()const;
   Double_t Angle(const TVector2 &v)const;
   TVector2 Cross(const TVector2 &v)const;
   Double_t Triple(const TVector2 &v1, const TVector2 &v2)const;
};

/////////////////////////////////////////////////////////////////////////////
///

TVector2 &TVector2::operator+=(const TVector2 &vv)
{
   fCo[0] += vv[0]; fCo[1] += vv[1];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 &TVector2::operator-=(const TVector2 &vv)
{
   fCo[0] -= vv[0]; fCo[1] -= vv[1];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 &TVector2::operator *= (Double_t s)
{
   fCo[0] *= s; fCo[1] *= s; return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 &TVector2::operator /= (Double_t s)
{
   return *this *= 1. / s;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator + (const TVector2 &v1, const TVector2 &v2)
{
   return TVector2(v1[0] + v2[0], v1[1] + v2[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator - (const TVector2 &v1, const TVector2 &v2)
{
   return TVector2(v1[0] - v2[0], v1[1] - v2[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator - (const TVector2 &v)
{
   return TVector2(-v[0], -v[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator * (const TVector2 &v, Double_t s)
{
   return TVector2(v[0] * s, v[1] * s);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator * (Double_t s, const TVector2 &v)
{
   return v * s;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator / (const TVector2 & v, Double_t s)
{
   return v * (1.0 / s);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector2::Dot(const TVector2 &vv)const
{
   return fCo[0] * vv[0] + fCo[1] * vv[1];
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector2::Length2()const
{
   return Dot(*this);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector2::Length()const
{
   return TMath::Sqrt(Length2());
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 TVector2::Absolute()const
{
   return TVector2(TMath::Abs(fCo[0]), TMath::Abs(fCo[1]));
}

/////////////////////////////////////////////////////////////////////////////
///

Bool_t TVector2::FuzzyZero()const
{
   return fuzzy_zero2(Length2());
}

/////////////////////////////////////////////////////////////////////////////
///

void TVector2::Normalize()
{
   *this /= Length();
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 TVector2::Normalized()const
{
   return *this / Length();
}

/////////////////////////////////////////////////////////////////////////////
///

void TVector2::Scale(Double_t xx, Double_t yy)
{
   fCo[0] *= xx; fCo[1] *= yy;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 TVector2::Scaled(Double_t xx, Double_t yy)const
{
   return TVector2(fCo[0] * xx, fCo[1] * yy);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector2::Angle(const TVector2 &vv)const
{
   Double_t s = TMath::Sqrt(Length2() * vv.Length2());
   return TMath::ACos(Dot(vv) / s);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t  dot(const TVector2 &v1, const TVector2 &v2)
{
   return v1.Dot(v2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t length2(const TVector2 &v)
{
   return v.Length2();
}

/////////////////////////////////////////////////////////////////////////////

Double_t length(const TVector2 &v)
{
   return v.Length();
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_zero(const TVector2 &v)
{
   return v.FuzzyZero();
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_equal(const TVector2 &v1, const TVector2 &v2)
{
   return fuzzy_zero(v1 - v2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t Angle(const TVector2 &v1, const TVector2 &v2)
{
   return v1.Angle(v2);
}

class TPoint2 : public TVector2 {
public:
   TPoint2(){}
   TPoint2(const Double_t *v) : TVector2(v){}
   TPoint2(Double_t x, Double_t y) : TVector2(x, y) {}

   TPoint2 &operator += (const TVector2 &v);
   TPoint2 &operator -= (const TVector2 &v);
   TPoint2 &operator = (const TVector2 &v);

   Double_t Distance(const TPoint2 &p)const;
   Double_t Distance2(const TPoint2 &p)const;
   TPoint2   Lerp(const TPoint2 &p, Double_t t)const;
};

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 &TPoint2::operator += (const TVector2 &v)
{
   fCo[0] += v[0]; fCo[1] += v[1];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 &TPoint2::operator-=(const TVector2& v)
{
   fCo[0] -= v[0]; fCo[1] -= v[1];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 &TPoint2::operator=(const TVector2& v)
{
   fCo[0] = v[0]; fCo[1] = v[1];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPoint2::Distance(const TPoint2& p)const
{
   return (p - *this).Length();
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPoint2::Distance2(const TPoint2& p)const
{
   return (p - *this).Length2();
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 TPoint2::Lerp(const TPoint2 &p, Double_t t)const
{
   return TPoint2(fCo[0] + (p[0] - fCo[0]) * t,
                  fCo[1] + (p[1] - fCo[1]) * t);
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 operator + (const TPoint2 &p, const TVector2 &v)
{
   return TPoint2(p[0] + v[0], p[1] + v[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 operator - (const TPoint2 &p, const TVector2 &v)
{
   return TPoint2(p[0] - v[0], p[1] - v[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector2 operator - (const TPoint2 &p1, const TPoint2 &p2)
{
   return TVector2(p1[0] - p2[0], p1[1] - p2[1]);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t distance(const TPoint2 &p1, const TPoint2 &p2)
{
   return p1.Distance(p2);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t distance2(const TPoint2 &p1, const TPoint2 &p2)
{
   return p1.Distance2(p2);
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint2 lerp(const TPoint2 &p1, const TPoint2 &p2, Double_t t)
{
   return p1.Lerp(p2, t);
}

class Tuple3 {
protected:
   Double_t fCo[3];
public:
   Tuple3(){SetValue(0, 0, 0);}
   Tuple3(const Double_t *v){SetValue(v);}
   Tuple3(Double_t xx, Double_t yy, Double_t zz){SetValue(xx, yy, zz);}

   Double_t       &operator [] (Int_t i){return fCo[i];}
   const Double_t &operator [] (Int_t i)const{return fCo[i];}

   Double_t       &X(){return fCo[0];}
   const Double_t &X()const{return fCo[0];}
   Double_t       &Y(){return fCo[1];}
   const Double_t &Y()const{return fCo[1];}
   Double_t       &Z(){return fCo[2];}
   const Double_t &Z()const{return fCo[2];}

   Double_t       *GetValue(){return fCo;}
   const Double_t *GetValue()const{ return fCo; }

   void GetValue(Double_t *v)const
   {
      v[0] = Double_t(fCo[0]), v[1] = Double_t(fCo[1]), v[2] = Double_t(fCo[2]);
   }

   void SetValue(const Double_t *v)
   {
      fCo[0] = Double_t(v[0]), fCo[1] = Double_t(v[1]), fCo[2] = Double_t(v[2]);
   }
   void SetValue(Double_t xx, Double_t yy, Double_t zz)
   {
      fCo[0] = xx; fCo[1] = yy; fCo[2] = zz;
   }
};

/////////////////////////////////////////////////////////////////////////////

Bool_t operator==(const Tuple3& t1, const Tuple3& t2)
{
   return t1[0] == t2[0] && t1[1] == t2[1] && t1[2] == t2[2];
}

class TVector3 : public Tuple3 {
public:
   TVector3(){}
   TVector3(const Double_t *v) : Tuple3(v){}
   TVector3(Double_t xx, Double_t yy, Double_t zz) : Tuple3(xx, yy, zz){}

   TVector3 &operator += (const TVector3& v);
   TVector3 &operator -= (const TVector3& v);
   TVector3 &operator *= (Double_t s);
   TVector3 &operator /= (Double_t s);

   Double_t Dot(const TVector3& v)const;
   Double_t Length2()const;
   Double_t Length()const;
   TVector3  Absolute()const;
   void     NoiseGate(Double_t threshold);
   void     Normalize();
   TVector3  Normalized()const;
   TVector3  SafeNormalized()const;
   void     Scale(Double_t x, Double_t y, Double_t z);
   TVector3  Scaled(Double_t x, Double_t y, Double_t z)const;
   Bool_t   FuzzyZero()const;
   Double_t Angle(const TVector3 &v)const;
   TVector3  Cross(const TVector3 &v)const;
   Double_t Triple(const TVector3 &v1, const TVector3 &v2)const;
   Int_t    ClosestAxis()const;

   static TVector3 Random();
};

/////////////////////////////////////////////////////////////////////////////
///

TVector3 &TVector3::operator += (const TVector3 &v)
{
   fCo[0] += v[0]; fCo[1] += v[1]; fCo[2] += v[2];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 &TVector3::operator -= (const TVector3 &v)
{
   fCo[0] -= v[0]; fCo[1] -= v[1]; fCo[2] -= v[2];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 &TVector3::operator *= (Double_t s)
{
   fCo[0] *= s; fCo[1] *= s; fCo[2] *= s;
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 &TVector3::operator /= (Double_t s)
{
   return *this *= Double_t(1.0) / s;
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector3::Dot(const TVector3 &v)const
{
   return fCo[0] * v[0] + fCo[1] * v[1] + fCo[2] * v[2];
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector3::Length2()const
{
   return Dot(*this);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector3::Length()const
{
   return TMath::Sqrt(Length2());
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TVector3::Absolute()const
{
   return TVector3(TMath::Abs(fCo[0]), TMath::Abs(fCo[1]), TMath::Abs(fCo[2]));
}

/////////////////////////////////////////////////////////////////////////////
///

Bool_t TVector3::FuzzyZero()const
{
   return fuzzy_zero(Length2());
}

/////////////////////////////////////////////////////////////////////////////
///

void TVector3::NoiseGate(Double_t threshold)
{
   if (Length2() < threshold) SetValue(0., 0., 0.);
}

/////////////////////////////////////////////////////////////////////////////
///

void TVector3::Normalize()
{
   *this /= Length();
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator * (const TVector3 &v, Double_t s)
{
   return TVector3(v[0] * s, v[1] * s, v[2] * s);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator / (const TVector3 &v, Double_t s)
{
   return v * (1. / s);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TVector3::Normalized()const
{
   return *this / Length();
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TVector3::SafeNormalized()const
{
   Double_t len = Length();
   return fuzzy_zero(len) ? TVector3(1., 0., 0.):*this / len;
}

/////////////////////////////////////////////////////////////////////////////
///

void TVector3::Scale(Double_t xx, Double_t yy, Double_t zz)
{
   fCo[0] *= xx; fCo[1] *= yy; fCo[2] *= zz;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TVector3::Scaled(Double_t xx, Double_t yy, Double_t zz)const
{
   return TVector3(fCo[0] * xx, fCo[1] * yy, fCo[2] * zz);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector3::Angle(const TVector3 &v)const
{
   Double_t s = TMath::Sqrt(Length2() * v.Length2());
   return TMath::ACos(Dot(v) / s);
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TVector3::Cross(const TVector3 &v)const
{
   return TVector3(fCo[1] * v[2] - fCo[2] * v[1],
                   fCo[2] * v[0] - fCo[0] * v[2],
                   fCo[0] * v[1] - fCo[1] * v[0]);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TVector3::Triple(const TVector3 &v1, const TVector3 &v2)const
{
   return fCo[0] * (v1[1] * v2[2] - v1[2] * v2[1]) +
      fCo[1] * (v1[2] * v2[0] - v1[0] * v2[2]) +
      fCo[2] * (v1[0] * v2[1] - v1[1] * v2[0]);
}

/////////////////////////////////////////////////////////////////////////////
///

Int_t TVector3::ClosestAxis()const
{
   TVector3 a = Absolute();
   return a[0] < a[1] ? (a[1] < a[2] ? 2 : 1) : (a[0] < a[2] ? 2 : 0);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator + (const TVector3 &v1, const TVector3 &v2)
{
   return TVector3(v1[0] + v2[0], v1[1] + v2[1], v1[2] + v2[2]);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator - (const TVector3 &v1, const TVector3 &v2)
{
   return TVector3(v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2]);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator - (const TVector3 &v)
{
   return TVector3(-v[0], -v[1], -v[2]);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator * (Double_t s, const TVector3 &v)
{
   return v * s;
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator * (const TVector3 &v1, const TVector3 &v2)
{
   return TVector3(v1[0] * v2[0], v1[1] * v2[1], v1[2] * v2[2]);
}

/////////////////////////////////////////////////////////////////////////////

Double_t  dot(const TVector3 &v1, const TVector3 &v2)
{
   return v1.Dot(v2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t length2(const TVector3 &v)
{
   return v.Length2();
}

/////////////////////////////////////////////////////////////////////////////

Double_t length(const TVector3 &v)
{
   return v.Length();
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_zero(const TVector3 &v)
{
   return v.FuzzyZero();
}

/////////////////////////////////////////////////////////////////////////////

Bool_t fuzzy_equal(const TVector3 &v1, const TVector3 &v2)
{
   return fuzzy_zero(v1 - v2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t Angle(const TVector3 &v1, const TVector3 &v2)
{
   return v1.Angle(v2);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 cross(const TVector3 &v1, const TVector3 &v2)
{
   return v1.Cross(v2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t triple(const TVector3 &v1, const TVector3 &v2, const TVector3 &v3)
{
   return v1.Triple(v2, v3);
}

class TPoint3 : public TVector3 {
public:
   TPoint3(){}
   TPoint3(const Double_t *v) : TVector3(v) {}
   TPoint3(Double_t xx, Double_t yy, Double_t zz) : TVector3(xx, yy, zz) {}

   TPoint3 &operator += (const TVector3 &v);
   TPoint3 &operator -= (const TVector3 &v);
   TPoint3 &operator = (const TVector3 &v);

   Double_t Distance(const TPoint3 &p)const;
   Double_t Distance2(const TPoint3 &p)const;
   TPoint3   Lerp(const TPoint3 &p, Double_t t)const;
};

/////////////////////////////////////////////////////////////////////////////
///

TPoint3 &TPoint3::operator += (const TVector3 &v)
{
   fCo[0] += v[0]; fCo[1] += v[1]; fCo[2] += v[2];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint3 &TPoint3::operator-=(const TVector3 &v)
{
   fCo[0] -= v[0]; fCo[1] -= v[1]; fCo[2] -= v[2];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint3 &TPoint3::operator=(const TVector3 &v)
{
   fCo[0] = v[0]; fCo[1] = v[1]; fCo[2] = v[2];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPoint3::Distance(const TPoint3 &p)const
{
   return (p - *this).Length();
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPoint3::Distance2(const TPoint3 &p)const
{
   return (p - *this).Length2();
}

/////////////////////////////////////////////////////////////////////////////
///

TPoint3 TPoint3::Lerp(const TPoint3 &p, Double_t t)const
{
   return TPoint3(fCo[0] + (p[0] - fCo[0]) * t,
                  fCo[1] + (p[1] - fCo[1]) * t,
                  fCo[2] + (p[2] - fCo[2]) * t);
}

/////////////////////////////////////////////////////////////////////////////

TPoint3 operator + (const TPoint3 &p, const TVector3 &v)
{
   return TPoint3(p[0] + v[0], p[1] + v[1], p[2] + v[2]);
}

/////////////////////////////////////////////////////////////////////////////

TPoint3 operator - (const TPoint3 &p, const TVector3 &v)
{
   return TPoint3(p[0] - v[0], p[1] - v[1], p[2] - v[2]);
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator - (const TPoint3 &p1, const TPoint3 &p2)
{
   return TVector3(p1[0] - p2[0], p1[1] - p2[1], p1[2] - p2[2]);
}

/////////////////////////////////////////////////////////////////////////////

Double_t distance(const TPoint3 &p1, const TPoint3 &p2)
{
   return p1.Distance(p2);
}

/////////////////////////////////////////////////////////////////////////////

Double_t distance2(const TPoint3 &p1, const TPoint3 &p2)
{
   return p1.Distance2(p2);
}

/////////////////////////////////////////////////////////////////////////////

TPoint3 lerp(const TPoint3 &p1, const TPoint3 &p2, Double_t t)
{
   return p1.Lerp(p2, t);
}

class Tuple4 {
protected:
   Double_t fCo[4];

public:
   Tuple4(){SetValue(0, 0, 0, 0);}
   Tuple4(const Double_t *v){SetValue(v);}
   Tuple4(Double_t xx, Double_t yy, Double_t zz, Double_t ww)
   {
      SetValue(xx, yy, zz, ww);
   }

   Double_t       &operator [] (Int_t i){return fCo[i];}
   const Double_t &operator [] (Int_t i)const{return fCo[i];}

   Double_t       &X(){return fCo[0];}
   const Double_t &X()const{return fCo[0];}
   Double_t       &Y(){return fCo[1];}
   const Double_t &Y()const{return fCo[1];}
   Double_t       &Z(){return fCo[2];}
   const Double_t &Z()const{return fCo[2];}
   Double_t       &W(){return fCo[3];}
   const Double_t &W()const{return fCo[3];}

   Double_t       *GetValue(){return fCo;}
   const Double_t *GetValue()const{return fCo;}

   void GetValue(Double_t *v)const
   {
      v[0] = fCo[0]; v[1] = fCo[1]; v[2] = fCo[2]; v[3] = fCo[3];
   }

   void SetValue(const Double_t *v)
   {
      fCo[0] = v[0]; fCo[1] = v[1]; fCo[2] = v[2]; fCo[3] = v[3];
   }
   void SetValue(Double_t xx, Double_t yy, Double_t zz, Double_t ww)
   {
      fCo[0] = xx; fCo[1] = yy; fCo[2] = zz; fCo[3] = ww;
   }
};

/////////////////////////////////////////////////////////////////////////////

Bool_t operator == (const Tuple4 &t1, const Tuple4 &t2)
{
   return t1[0] == t2[0] && t1[1] == t2[1] && t1[2] == t2[2] && t1[3] == t2[3];
}

class TMatrix3x3 {
protected:
   TVector3 fEl[3];

public:
   TMatrix3x3(){}
   TMatrix3x3(const Double_t *m){SetValue(m);}
   TMatrix3x3(const TVector3 &euler){SetEuler(euler);}
   TMatrix3x3(const TVector3 &euler, const TVector3 &s)
   {
      SetEuler(euler); Scale(s[0], s[1], s[2]);
   }
   TMatrix3x3(Double_t xx, Double_t xy, Double_t xz,
              Double_t yx, Double_t yy, Double_t yz,
              Double_t zx, Double_t zy, Double_t zz)
   {
      SetValue(xx, xy, xz, yx, yy, yz, zx, zy, zz);
   }

   TVector3 &operator       [] (Int_t i){return fEl[i];}
   const TVector3 &operator [] (Int_t i)const{return fEl[i];}

   void SetValue(const Double_t *m)
   {
      fEl[0][0] = *m++; fEl[1][0] = *m++; fEl[2][0] = *m++; m++;
      fEl[0][1] = *m++; fEl[1][1] = *m++; fEl[2][1] = *m++; m++;
      fEl[0][2] = *m++; fEl[1][2] = *m++; fEl[2][2] = *m;
   }
   void SetValue(Double_t xx, Double_t xy, Double_t xz,
                 Double_t yx, Double_t yy, Double_t yz,
                 Double_t zx, Double_t zy, Double_t zz)
   {
      fEl[0][0] = xx; fEl[0][1] = xy; fEl[0][2] = xz;
      fEl[1][0] = yx; fEl[1][1] = yy; fEl[1][2] = yz;
      fEl[2][0] = zx; fEl[2][1] = zy; fEl[2][2] = zz;
   }
   void SetEuler(const TVector3 &euler)
   {
      Double_t ci = TMath::Cos(euler[0]);
      Double_t cj = TMath::Cos(euler[1]);
      Double_t ch = TMath::Cos(euler[2]);
      Double_t si = TMath::Sin(euler[0]);
      Double_t sj = TMath::Sin(euler[1]);
      Double_t sh = TMath::Sin(euler[2]);
      Double_t cc = ci * ch;
      Double_t cs = ci * sh;
      Double_t sc = si * ch;
      Double_t ss = si * sh;
      SetValue(cj * ch, sj * sc - cs, sj * cc + ss,
               cj * sh, sj * ss + cc, sj * cs - sc,
               -sj, cj * si, cj * ci);
   }

   void Scale(Double_t x, Double_t y, Double_t z)
   {
      fEl[0][0] *= x; fEl[0][1] *= y; fEl[0][2] *= z;
      fEl[1][0] *= x; fEl[1][1] *= y; fEl[1][2] *= z;
      fEl[2][0] *= x; fEl[2][1] *= y; fEl[2][2] *= z;
   }
   TMatrix3x3 Scaled(Double_t x, Double_t y, Double_t z)const
   {
      return TMatrix3x3(fEl[0][0] * x, fEl[0][1] * y, fEl[0][2] * z,
                        fEl[1][0] * x, fEl[1][1] * y, fEl[1][2] * z,
                        fEl[2][0] * x, fEl[2][1] * y, fEl[2][2] * z);
   }

   void SetIdentity()
   {
      SetValue(1., 0., 0., 0., 1., 0., 0., 0., 1.);
   }
   void GetValue(Double_t *m)const
   {
      *m++ = fEl[0][0]; *m++ = fEl[1][0]; *m++ = fEl[2][0]; *m++ = 0.0;
      *m++ = fEl[0][1]; *m++ = fEl[1][1]; *m++ = fEl[2][1]; *m++ = 0.0;
      *m++ = fEl[0][2]; *m++ = fEl[1][2]; *m++ = fEl[2][2]; *m   = 0.0;
   }

   TMatrix3x3 &operator *= (const TMatrix3x3 &m);
   Double_t Tdot(Int_t c, const TVector3 &v)const
   {
      return fEl[0][c] * v[0] + fEl[1][c] * v[1] + fEl[2][c] * v[2];
   }
   Double_t Cofac(Int_t r1, Int_t c1, Int_t r2, Int_t c2)const
   {
      return fEl[r1][c1] * fEl[r2][c2] - fEl[r1][c2] * fEl[r2][c1];
   }

   Double_t  Determinant()const;
   TMatrix3x3 Adjoint()const;
   TMatrix3x3 Absolute()const;
   TMatrix3x3 Transposed()const;
   void      Transpose();
   TMatrix3x3 Inverse()const;
   void      Invert();
};

/////////////////////////////////////////////////////////////////////////////
///

TMatrix3x3 &TMatrix3x3::operator *= (const TMatrix3x3 &m)
{
   SetValue(m.Tdot(0, fEl[0]), m.Tdot(1, fEl[0]), m.Tdot(2, fEl[0]),
            m.Tdot(0, fEl[1]), m.Tdot(1, fEl[1]), m.Tdot(2, fEl[1]),
            m.Tdot(0, fEl[2]), m.Tdot(1, fEl[2]), m.Tdot(2, fEl[2]));
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TMatrix3x3::Determinant()const
{
   return triple((*this)[0], (*this)[1], (*this)[2]);
}

/////////////////////////////////////////////////////////////////////////////
///

TMatrix3x3 TMatrix3x3::Absolute()const
{
   return TMatrix3x3(TMath::Abs(fEl[0][0]), TMath::Abs(fEl[0][1]), TMath::Abs(fEl[0][2]),
                     TMath::Abs(fEl[1][0]), TMath::Abs(fEl[1][1]), TMath::Abs(fEl[1][2]),
                     TMath::Abs(fEl[2][0]), TMath::Abs(fEl[2][1]), TMath::Abs(fEl[2][2]));
}

/////////////////////////////////////////////////////////////////////////////
///

TMatrix3x3 TMatrix3x3::Transposed()const
{
   return TMatrix3x3(fEl[0][0], fEl[1][0], fEl[2][0],
                     fEl[0][1], fEl[1][1], fEl[2][1],
                     fEl[0][2], fEl[1][2], fEl[2][2]);
}

/////////////////////////////////////////////////////////////////////////////
///

void TMatrix3x3::Transpose()
{
   *this = Transposed();
}

/////////////////////////////////////////////////////////////////////////////
///

TMatrix3x3 TMatrix3x3::Adjoint()const
{
   return TMatrix3x3(Cofac(1, 1, 2, 2), Cofac(0, 2, 2, 1), Cofac(0, 1, 1, 2),
                     Cofac(1, 2, 2, 0), Cofac(0, 0, 2, 2), Cofac(0, 2, 1, 0),
                     Cofac(1, 0, 2, 1), Cofac(0, 1, 2, 0), Cofac(0, 0, 1, 1));
}

/////////////////////////////////////////////////////////////////////////////
///

TMatrix3x3 TMatrix3x3::Inverse()const
{
   TVector3 co(Cofac(1, 1, 2, 2), Cofac(1, 2, 2, 0), Cofac(1, 0, 2, 1));
   Double_t det = dot((*this)[0], co);
   Double_t s = 1. / det;
   return TMatrix3x3(co[0] * s, Cofac(0, 2, 2, 1) * s, Cofac(0, 1, 1, 2) * s,
                     co[1] * s, Cofac(0, 0, 2, 2) * s, Cofac(0, 2, 1, 0) * s,
                     co[2] * s, Cofac(0, 1, 2, 0) * s, Cofac(0, 0, 1, 1) * s);
}

/////////////////////////////////////////////////////////////////////////////
///

void TMatrix3x3::Invert()
{
   *this = Inverse();
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator * (const TMatrix3x3& m, const TVector3& v)
{
   return TVector3(dot(m[0], v), dot(m[1], v), dot(m[2], v));
}

/////////////////////////////////////////////////////////////////////////////

TVector3 operator * (const TVector3& v, const TMatrix3x3& m)
{
   return TVector3(m.Tdot(0, v), m.Tdot(1, v), m.Tdot(2, v));
}

/////////////////////////////////////////////////////////////////////////////

TMatrix3x3 operator * (const TMatrix3x3 &m1, const TMatrix3x3 &m2)
{
   return TMatrix3x3(m2.Tdot(0, m1[0]), m2.Tdot(1, m1[0]), m2.Tdot(2, m1[0]),
                     m2.Tdot(0, m1[1]), m2.Tdot(1, m1[1]), m2.Tdot(2, m1[1]),
                     m2.Tdot(0, m1[2]), m2.Tdot(1, m1[2]), m2.Tdot(2, m1[2]));
}

/////////////////////////////////////////////////////////////////////////////

TMatrix3x3 mmult_transpose_left(const TMatrix3x3 &m1, const TMatrix3x3 &m2)
{
   return TMatrix3x3(m1[0][0] * m2[0][0] + m1[1][0] * m2[1][0] + m1[2][0] * m2[2][0],
                     m1[0][0] * m2[0][1] + m1[1][0] * m2[1][1] + m1[2][0] * m2[2][1],
                     m1[0][0] * m2[0][2] + m1[1][0] * m2[1][2] + m1[2][0] * m2[2][2],
                     m1[0][1] * m2[0][0] + m1[1][1] * m2[1][0] + m1[2][1] * m2[2][0],
                     m1[0][1] * m2[0][1] + m1[1][1] * m2[1][1] + m1[2][1] * m2[2][1],
                     m1[0][1] * m2[0][2] + m1[1][1] * m2[1][2] + m1[2][1] * m2[2][2],
                     m1[0][2] * m2[0][0] + m1[1][2] * m2[1][0] + m1[2][2] * m2[2][0],
                     m1[0][2] * m2[0][1] + m1[1][2] * m2[1][1] + m1[2][2] * m2[2][1],
                     m1[0][2] * m2[0][2] + m1[1][2] * m2[1][2] + m1[2][2] * m2[2][2]);
}

/////////////////////////////////////////////////////////////////////////////

TMatrix3x3 mmult_transpose_right(const TMatrix3x3& m1, const TMatrix3x3& m2)
{
   return TMatrix3x3(m1[0].Dot(m2[0]), m1[0].Dot(m2[1]), m1[0].Dot(m2[2]),
                     m1[1].Dot(m2[0]), m1[1].Dot(m2[1]), m1[1].Dot(m2[2]),
                     m1[2].Dot(m2[0]), m1[2].Dot(m2[1]), m1[2].Dot(m2[2]));
}

class TLine3 {
private :
   Bool_t   fBounds[2];
   Double_t fParams[2];
   TPoint3   fOrigin;
   TVector3  fDir;

public :
   TLine3();
   TLine3(const TPoint3 &p1, const TPoint3 &p2);
   TLine3(const TPoint3 &p1, const TVector3 &v);
   TLine3(const TPoint3 &p1, const TVector3 &v, Bool_t bound1, Bool_t bound2);

   static TLine3   InfiniteRay(const TPoint3 &p1, const TVector3 &v);
   const TVector3 &Direction()const {return fDir;}
   const TPoint3  &Origin()const{ return fOrigin;}

   Bool_t Bounds(Int_t i)const
   {
      return (i == 0 ? fBounds[0] : fBounds[1]);
   }
   Bool_t &Bounds(Int_t i)
   {
      return (i == 0 ? fBounds[0] : fBounds[1]);
   }
   const Double_t &Param(Int_t i)const
   {
      return (i == 0 ? fParams[0] : fParams[1]);
   }
   Double_t &Param(Int_t i)
   {
      return (i == 0 ? fParams[0] : fParams[1]);
   }
   TVector3 UnboundSmallestVector(const TPoint3 &point)const
   {
      TVector3 diff(fOrigin - point);
      return diff - fDir * diff.Dot(fDir);
   }
   Double_t UnboundClosestParameter(const TPoint3 &point)const
   {
      TVector3 diff(fOrigin-point);
      return diff.Dot(fDir);
   }
   Double_t UnboundDistance(const TPoint3& point)const
   {
      return UnboundSmallestVector(point).Length();
   }
   Bool_t IsParameterOnLine(const Double_t &t) const
   {
      return ((fParams[0] - epsilon < t) || (!fBounds[0])) && ((fParams[1] > t + epsilon) || (!fBounds[1]));
   }
};

/////////////////////////////////////////////////////////////////////////////
///

TLine3::TLine3() : fOrigin(0,0,0), fDir(1,0,0)
{
   fBounds[0] = kFALSE; fBounds[1] = kFALSE;
   fParams[0] = 0; fParams[1] = 1;
}

/////////////////////////////////////////////////////////////////////////////
///

TLine3::TLine3(const TPoint3 &p1, const TPoint3 &p2) : fOrigin(p1), fDir(p2-p1)
{
   fBounds[0] = kTRUE; fBounds[1] = kTRUE;
   fParams[0] = 0; fParams[1] = 1;
}

/////////////////////////////////////////////////////////////////////////////
///

TLine3::TLine3(const TPoint3 &p1, const TVector3 &v): fOrigin(p1), fDir(v)
{
   fBounds[0] = kFALSE; fBounds[1] = kFALSE;
   fParams[0] = 0; fParams[1] = 1;
}

/////////////////////////////////////////////////////////////////////////////
///

TLine3::TLine3(const TPoint3 &p1, const TVector3 &v, Bool_t bound1, Bool_t bound2)
   : fOrigin(p1), fDir(v)
{
   fBounds[0] = bound1; fBounds[1] = bound2;
   fParams[0] = 0; fParams[1] = 1;
}

class TPlane3 : public Tuple4 {
public :
   TPlane3(const TVector3 &a, const TVector3 &b, const TVector3 &c);
   TPlane3(const TVector3 &n, const TVector3 &p);
   TPlane3();
   TPlane3(const TPlane3 & p):Tuple4(p){}

   TVector3  Normal()const;
   Double_t  Scalar()const;
   void     Invert();
   Double_t SignedDistance(const TVector3 &)const;

   TPlane3 &operator = (const TPlane3 & rhs);
};

/////////////////////////////////////////////////////////////////////////////
///

TPlane3::TPlane3(const TVector3 &a, const TVector3 &b, const TVector3 &c)
{
   TVector3 l1 = b-a;
   TVector3 l2 = c-b;
   TVector3 n = l1.Cross(l2);
   n = n.SafeNormalized();
   Double_t d = n.Dot(a);
   fCo[0] = n.X(); fCo[1] = n.Y(); fCo[2] = n.Z(); fCo[3] = -d;
}

/////////////////////////////////////////////////////////////////////////////
///

TPlane3::TPlane3(const TVector3 &n, const TVector3 &p)
{
   TVector3 mn = n.SafeNormalized();
   Double_t md = mn.Dot(p);
   fCo[0] = mn.X(); fCo[1] = mn.Y(); fCo[2] = mn.Z(); fCo[3] = -md;
}

/////////////////////////////////////////////////////////////////////////////
///

TPlane3::TPlane3() : Tuple4()
{
   fCo[0] = 1.; fCo[1] = 0.;
   fCo[2] = 0.; fCo[3] = 0.;
}

/////////////////////////////////////////////////////////////////////////////
///

TVector3 TPlane3::Normal()const
{
   return TVector3(fCo[0],fCo[1],fCo[2]);
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPlane3::Scalar()const
{
   return fCo[3];
}

/////////////////////////////////////////////////////////////////////////////
///

void TPlane3::Invert()
{
   fCo[0] = -fCo[0]; fCo[1] = -fCo[1]; fCo[2] = -fCo[2]; fCo[3] = -fCo[3];
}

/////////////////////////////////////////////////////////////////////////////
///

TPlane3 &TPlane3::operator = (const TPlane3 &rhs)
{
   fCo[0] = rhs.fCo[0]; fCo[1] = rhs.fCo[1]; fCo[2] = rhs.fCo[2]; fCo[3] = rhs.fCo[3];
   return *this;
}

/////////////////////////////////////////////////////////////////////////////
///

Double_t TPlane3::SignedDistance(const TVector3 &v)const
{
   return Normal().Dot(v) + fCo[3];
}

class TBBox {
   friend Bool_t intersect(const TBBox &a, const TBBox &b);

private:
   TPoint3  fCenter;
   TVector3 fExtent;

public:
   TBBox(){}
   TBBox(const TPoint3 &mini, const TPoint3 &maxi){SetValue(mini,maxi);}

   const TPoint3   &Center()const{return fCenter;}
   const TVector3  &Extent()const{return fExtent;}
   TPoint3         &Center(){return fCenter;}
   TVector3        &Extent(){return fExtent;}

   void SetValue(const TPoint3 &mini,const TPoint3 &maxi)
   {
      fExtent = (maxi - mini) / 2;
      fCenter = mini + fExtent;
   }
   void Enclose(const TBBox &a, const TBBox &b)
   {
      TPoint3 lower(
                    TMath::Min(a.Lower(0), b.Lower(0)),
                    TMath::Min(a.Lower(1), b.Lower(1)),
                    TMath::Min(a.Lower(2), b.Lower(2))
                    );
      TPoint3 upper(
                    TMath::Max(a.Upper(0), b.Upper(0)),
                    TMath::Max(a.Upper(1), b.Upper(1)),
                    TMath::Max(a.Upper(2), b.Upper(2))
                    );
      SetValue(lower, upper);
   }

   void SetEmpty()
   {
      fCenter.SetValue(0., 0., 0.);
      fExtent.SetValue(-infinity, -infinity, -infinity);
   }
   void Include(const TPoint3 &p)
   {
      TPoint3 lower(
                    TMath::Min(Lower(0), p[0]),
                    TMath::Min(Lower(1), p[1]),
                    TMath::Min(Lower(2), p[2])
                    );
      TPoint3 upper(
                    TMath::Max(Upper(0), p[0]),
                    TMath::Max(Upper(1), p[1]),
                    TMath::Max(Upper(2), p[2])
                    );
      SetValue(lower, upper);
   }

   void Include(const TBBox &b)
   {
      Enclose(*this, b);
   }
   Double_t Lower(Int_t i)const
   {
      return fCenter[i] - fExtent[i];
   }
   Double_t Upper(Int_t i)const
   {
      return fCenter[i] + fExtent[i];
   }
   TPoint3 Lower()const
   {
      return fCenter - fExtent;
   }
   TPoint3 Upper()const
   {
      return fCenter + fExtent;
   }
   Double_t Size()const
   {
      return TMath::Max(TMath::Max(fExtent[0], fExtent[1]), fExtent[2]);
   }
   Int_t LongestAxis()const
   {
      return fExtent.ClosestAxis();
   }
   Bool_t IntersectXRay(const TPoint3 &xBase)const
   {
      if (xBase[0] <= Upper(0)) {
         if (xBase[1] <= Upper(1) && xBase[1] >= Lower(1)) {
            if (xBase[2] <= Upper(2) && xBase[2] >= Lower(2)) {return kTRUE;}
         }
      }
      return kFALSE;
   }
};

/////////////////////////////////////////////////////////////////////////////

Bool_t intersect(const TBBox &a, const TBBox &b)
{
   return TMath::Abs(a.fCenter[0] - b.fCenter[0]) <= a.fExtent[0] + b.fExtent[0] &&
      TMath::Abs(a.fCenter[1] - b.fCenter[1]) <= a.fExtent[1] + b.fExtent[1] &&
      TMath::Abs(a.fCenter[2] - b.fCenter[2]) <= a.fExtent[2] + b.fExtent[2];
}

class TBBoxNode {
public:
   enum ETagType {kLeaf, kInternal};
   TBBox     fBBox;
   ETagType fTag;
};

class TBBoxLeaf : public TBBoxNode {
public:
   Int_t fPolyIndex;

   TBBoxLeaf() : fPolyIndex(0) {}
   TBBoxLeaf(Int_t polyIndex, const TBBox &bbox) : fPolyIndex(polyIndex)
   {
      fBBox = bbox;
      fTag = kLeaf;
   }
};

typedef TBBoxLeaf *LeafPtr_t;
typedef TBBoxNode *NodePtr_t;

class TBBoxInternal : public TBBoxNode {
public:
   NodePtr_t fLeftSon;
   NodePtr_t fRightSon;
   TBBoxInternal() : fLeftSon(0) ,fRightSon(0) {}
   TBBoxInternal(Int_t n, LeafPtr_t leafIt);
};

typedef TBBoxInternal* InternalPtr_t;

class TBBoxTree {
private:
   Int_t         fBranch;
   LeafPtr_t     fLeaves;
   InternalPtr_t fInternals;
   Int_t         fNumLeaves;

public :
   TBBoxTree() : fBranch(0), fLeaves(0), fInternals(0), fNumLeaves(0) {}
   NodePtr_t RootNode()const{return fInternals;}
   ~TBBoxTree()
   {
      delete[] fLeaves;
      delete[] fInternals;
   }
   void BuildTree(LeafPtr_t leaves, Int_t numLeaves);

private :
   void RecursiveTreeBuild(Int_t n, LeafPtr_t leafIt);
};

/////////////////////////////////////////////////////////////////////////////
///

TBBoxInternal::TBBoxInternal(Int_t n, LeafPtr_t leafIt) :
   fLeftSon(0) ,fRightSon(0)
{
   fTag = kInternal;
   fBBox.SetEmpty();
   for (Int_t i=0;i<n;i++)
      fBBox.Include(leafIt[i].fBBox);
}

/////////////////////////////////////////////////////////////////////////////
///

void TBBoxTree::BuildTree(LeafPtr_t leaves, Int_t numLeaves)
{
   fBranch = 0;
   fLeaves = leaves;
   fNumLeaves = numLeaves;
   fInternals = new TBBoxInternal[numLeaves];
   RecursiveTreeBuild(fNumLeaves,fLeaves);
}

/////////////////////////////////////////////////////////////////////////////
///

void TBBoxTree::RecursiveTreeBuild(Int_t n, LeafPtr_t leafIt)
{
   fInternals[fBranch] = TBBoxInternal(n,leafIt);
   TBBoxInternal &aBBox = fInternals[fBranch];
   fBranch++;

   Int_t axis = aBBox.fBBox.LongestAxis();
   Int_t i = 0, mid = n;

   while (i < mid) {
      if (leafIt[i].fBBox.Center()[axis] < aBBox.fBBox.Center()[axis]) {
         ++i;
      } else {
         --mid;
         std::swap(leafIt[i], leafIt[mid]);
      }
   }

   if (mid == 0 || mid == n) {
      mid = n / 2;
   }
   if (mid >= 2) {
      aBBox.fRightSon = fInternals + fBranch;
      RecursiveTreeBuild(mid,leafIt);
   } else {
      aBBox.fRightSon = leafIt;
   }
   if (n - mid >= 2) {
      aBBox.fLeftSon = fInternals + fBranch;
      RecursiveTreeBuild(n - mid, leafIt + mid);
   } else {
      aBBox.fLeftSon = leafIt + mid;
   }
}

class TBlenderVProp {
private:
   Int_t fVertexIndex;

public:
   TBlenderVProp(Int_t vIndex) : fVertexIndex(vIndex){}
   TBlenderVProp(Int_t vIndex, const TBlenderVProp &,
                 const TBlenderVProp &, const Double_t &)
   {
      fVertexIndex = vIndex;
   }
   TBlenderVProp() : fVertexIndex(-1){}
   operator Int_t()const
   {
      return fVertexIndex;
   }
   TBlenderVProp &operator = (Int_t i)
   {
      fVertexIndex = i; return *this;
   }
};

template <class TMesh>
class TPolygonGeometry {
public:
   typedef typename TMesh::Polygon TPolygon;

private:
   const TMesh    &fMesh;
   const TPolygon &fPoly;
public:
   TPolygonGeometry(const TMesh &mesh, Int_t pIndex)
      : fMesh(mesh), fPoly(mesh.Polys()[pIndex])
   {}
   TPolygonGeometry(const TMesh &mesh, const TPolygon &poly)
      : fMesh(mesh), fPoly(poly)
   {}
   const TPoint3 &operator [] (Int_t i)const
   {
      return fMesh.Verts()[fPoly[i]].Pos();
   }
   Int_t Size()const
   {
      return fPoly.Size();
   }

private:
   TPolygonGeometry(const TPolygonGeometry &);
   TPolygonGeometry& operator = (TPolygonGeometry &);
};

template <class TPolygon, class TVertex>
class TMesh : public TBaseMesh {
public:
   typedef std::vector<TVertex> VLIST;
   typedef std::vector<TPolygon> PLIST;
   typedef TPolygon Polygon;
   typedef TVertex Vertex;
   typedef TPolygonGeometry<TMesh> TGBinder;

private:
   VLIST fVerts;
   PLIST fPolys;

public:
   VLIST &Verts() { return fVerts; }
   const VLIST &Verts() const { return fVerts; }
   PLIST &Polys() { return fPolys; }
   const PLIST &Polys() const { return fPolys; }

   // TBaseMesh's final-overriders
   Int_t NumberOfPolys() const override { return fPolys.size(); }
   Int_t NumberOfVertices() const override { return fVerts.size(); }
   Int_t SizeOfPoly(Int_t polyIndex) const override { return fPolys[polyIndex].Size(); }
   const Double_t *GetVertex(Int_t vertexNum) const override { return fVerts[vertexNum].GetValue(); }

   Int_t GetVertexIndex(Int_t polyNum, Int_t vertexNum) const override { return fPolys[polyNum][vertexNum]; }
};

const Int_t cofacTable[3][2] = {{1,2}, {0,2}, {0,1}};

{'fraction_non_alphanumeric': 0.3428046948879123, 'fraction_numerical': 0.039372531568114814, 'mean_word_length': 1.4841366349767853, 'pattern_counts': {'":': 0, '<': 71, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1758824', 'n_tokens_mistral': 32068, 'n_tokens_neox': 31243, 'n_words': 6677}