TGeoHalfSpace.cxx

Go to the documentation of this file.

// @(#):$Id$
// Author: Mihaela Gheata   03/08/04
 
/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/
 
/** \class TGeoHalfSpace
\ingroup Shapes_classes
 
A half space is limited just by a plane, defined by a point and the
normal direction. The point lies on the plane and the normal vector
points outside the half space. The half space is the only shape
which is infinite and can be used only in Boolean operations that
result in non-infinite composite shapes (see also TGeoCompositeShape).
A half space has to be defined using the constructor:
 
~~~{.cpp}
TGeoHalfSpace (const char *name, Double_t *point[3],
Double_t *norm[3]);
~~~
 
*/
 
#include <iostream>
#include "TGeoHalfSpace.h"
#include "TMath.h"
 
ClassImp(TGeoHalfSpace);
 
////////////////////////////////////////////////////////////////////////////////
/// Dummy constructor
 
TGeoHalfSpace::TGeoHalfSpace()
{
   SetShapeBit(TGeoShape::kGeoHalfSpace);
   SetShapeBit(TGeoShape::kGeoInvalidShape);
   memset(fP, 0, 3*sizeof(Double_t));
   memset(fN, 0, 3*sizeof(Double_t));
}
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor with name, point on the plane and normal
 
TGeoHalfSpace::TGeoHalfSpace(const char *name, Double_t *p, Double_t *n)
              :TGeoBBox(name, 0,0,0)
{
   SetShapeBit(TGeoShape::kGeoHalfSpace);
   SetShapeBit(TGeoShape::kGeoInvalidShape);
   Double_t param[6];
   memcpy(param, p, 3*sizeof(Double_t));
   memcpy(&param[3], n, 3*sizeof(Double_t));
   SetDimensions(param);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Default constructor specifying minimum and maximum radius
 
TGeoHalfSpace::TGeoHalfSpace(Double_t *param)
              :TGeoBBox(0,0,0)
{
   SetShapeBit(TGeoShape::kGeoHalfSpace);
   SetShapeBit(TGeoShape::kGeoInvalidShape);
   SetDimensions(param);
}
 
////////////////////////////////////////////////////////////////////////////////
/// destructor
 
TGeoHalfSpace::~TGeoHalfSpace()
{
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute normal to closest surface from POINT.
 
void TGeoHalfSpace::ComputeNormal(const Double_t * /*point*/, const Double_t *dir, Double_t *norm)
{
   memcpy(norm, fN, 3*sizeof(Double_t));
   if (norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2]<0) {
      norm[0] = -norm[0];
      norm[1] = -norm[1];
      norm[2] = -norm[2];
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// test if point is inside the half-space
 
Bool_t TGeoHalfSpace::Contains(const Double_t *point) const
{
   Double_t r[3];
   r[0] = fP[0]-point[0];
   r[1] = fP[1]-point[1];
   r[2] = fP[2]-point[2];
   Double_t rdotn = r[0]*fN[0]+r[1]*fN[1]+r[2]*fN[2];
   if (rdotn < 0) return kFALSE;
   return kTRUE;
}
 
////////////////////////////////////////////////////////////////////////////////
/// A half-space does not have a mesh primitive
 
Int_t TGeoHalfSpace::DistancetoPrimitive(Int_t /*px*/, Int_t /*py*/)
{
   return 999;
}
 
////////////////////////////////////////////////////////////////////////////////
/// compute distance from inside point to the plane
 
Double_t TGeoHalfSpace::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
{
   Double_t r[3];
   r[0] = fP[0]-point[0];
   r[1] = fP[1]-point[1];
   r[2] = fP[2]-point[2];
   Double_t rdotn = r[0]*fN[0]+r[1]*fN[1]+r[2]*fN[2];
   if (iact<3 && safe) {
      *safe = rdotn;
      if (iact==0) return TGeoShape::Big();
      if ((iact==1) && (*safe>step)) return TGeoShape::Big();
   }
   // compute distance to plane
   Double_t snxt = TGeoShape::Big();
   Double_t ddotn = dir[0]*fN[0]+dir[1]*fN[1]+dir[2]*fN[2];
   if (TMath::Abs(ddotn)<TGeoShape::Tolerance()) return snxt;
   snxt = rdotn/ddotn;
   if (snxt<0) return TGeoShape::Big();
   return snxt;
}
 
////////////////////////////////////////////////////////////////////////////////
/// compute distance from inside point to the plane
 
Double_t TGeoHalfSpace::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
{
   Double_t r[3];
   r[0] = fP[0]-point[0];
   r[1] = fP[1]-point[1];
   r[2] = fP[2]-point[2];
   Double_t rdotn = r[0]*fN[0]+r[1]*fN[1]+r[2]*fN[2];
   if (iact<3 && safe) {
      *safe = -rdotn;
      if (iact==0) return TGeoShape::Big();
      if ((iact==1) && (step<*safe)) return TGeoShape::Big();
   }
   // compute distance to plane
   Double_t snxt = TGeoShape::Big();
   Double_t ddotn = dir[0]*fN[0]+dir[1]*fN[1]+dir[2]*fN[2];
   if (TMath::Abs(ddotn)<TGeoShape::Tolerance()) return snxt;
   snxt = rdotn/ddotn;
   if (snxt<0) return TGeoShape::Big();
   return snxt;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Divide the shape along one axis.
 
TGeoVolume *TGeoHalfSpace::Divide(TGeoVolume * /*voldiv*/, const char * /*divname*/, Int_t /*iaxis*/, Int_t /*ndiv*/,
                             Double_t /*start*/, Double_t /*step*/)
{
   Error("Divide", "Half-spaces cannot be divided");
   return 0;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Returns numbers of vertices, segments and polygons composing the shape mesh.
 
void TGeoHalfSpace::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
{
   nvert = 0;
   nsegs = 0;
   npols = 0;
}
 
////////////////////////////////////////////////////////////////////////////////
/// print shape parameters
 
void TGeoHalfSpace::InspectShape() const
{
   printf("*** Shape %s: TGeoHalfSpace ***\n", GetName());
   printf("    Point    : %11.5f, %11.5f, %11.5f\n", fP[0], fP[1], fP[2]);
   printf("    Normal   : %11.5f, %11.5f, %11.5f\n", fN[0], fN[1], fN[2]);
}
 
////////////////////////////////////////////////////////////////////////////////
/// computes the closest distance from given point to this shape, according
/// to option. The matching point on the shape is stored in spoint.
 
Double_t TGeoHalfSpace::Safety(const Double_t *point, Bool_t /*in*/) const
{
   Double_t r[3];
   r[0] = fP[0]-point[0];
   r[1] = fP[1]-point[1];
   r[2] = fP[2]-point[2];
   Double_t rdotn = r[0]*fN[0]+r[1]*fN[1]+r[2]*fN[2];
   return TMath::Abs(rdotn);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Save a primitive as a C++ statement(s) on output stream "out".
 
void TGeoHalfSpace::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
{
   if (TObject::TestBit(kGeoSavePrimitive)) return;
   out << "   // Shape: " << GetName() << " type: " << ClassName() << std::endl;
   out << "   point[0] = " << fP[0] << ";" << std::endl;
   out << "   point[1] = " << fP[1] << ";" << std::endl;
   out << "   point[2] = " << fP[2] << ";" << std::endl;
   out << "   norm[0]  = " << fN[0] << ";" << std::endl;
   out << "   norm[1]  = " << fN[1] << ";" << std::endl;
   out << "   norm[2]  = " << fN[2] << ";" << std::endl;
   out << "   TGeoShape *" << GetPointerName() << " = new TGeoHalfSpace(\"" << GetName() << "\", point,norm);" << std::endl;
   TObject::SetBit(TGeoShape::kGeoSavePrimitive);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Set half-space parameters as stored in an array.
 
void TGeoHalfSpace::SetDimensions(Double_t *param)
{
   memcpy(fP, param, 3*sizeof(Double_t));
   memcpy(fN, &param[3], 3*sizeof(Double_t));
   Double_t nsq = TMath::Sqrt(fN[0]*fN[0]+fN[1]*fN[1]+fN[2]*fN[2]);
   fN[0] /= nsq;
   fN[1] /= nsq;
   fN[2] /= nsq;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check the inside status for each of the points in the array.
/// Input: Array of point coordinates + vector size
/// Output: Array of Booleans for the inside of each point
 
void TGeoHalfSpace::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
{
   for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute the normal for an array o points so that norm.dot.dir is positive
/// Input: Arrays of point coordinates and directions + vector size
/// Output: Array of normal directions
 
void TGeoHalfSpace::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
{
   for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute distance from array of input points having directions specified by dirs. Store output in dists
 
void TGeoHalfSpace::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
{
   for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute distance from array of input points having directions specified by dirs. Store output in dists
 
void TGeoHalfSpace::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
{
   for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute safe distance from each of the points in the input array.
/// Input: Array of point coordinates, array of statuses for these points, size of the arrays
/// Output: Safety values
 
void TGeoHalfSpace::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
{
   for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
}