a Belle crystal in Geant4 More...

#include <BelleCrystal.h>

Inheritance diagram for BelleCrystal:

Public Member Functions
	BelleCrystal (const G4String &pName)
	Constructor for "nominal" BelleCrystal.

	BelleCrystal (const G4String &pName, int, const G4ThreeVector *)
	Constructor.

virtual	~BelleCrystal ()
	Destructor.

Plane_t	GetSidePlane (G4int n) const
	Get side plane.

G4double	GetCubicVolume ()
	get the Cubic volume

G4double	GetSurfaceArea ()
	get the surface area

void	ComputeDimensions (G4VPVParameterisation p, const G4int n, const G4VPhysicalVolume pRep)
	compute the dimensions

G4bool	CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
	calculate the extent of the volume

EInside	Inside (const G4ThreeVector &p) const
	Return whether point inside/outside/on surface, using tolerance.

G4ThreeVector	SurfaceNormal (const G4ThreeVector &p) const
	Calculate side nearest to p, and return normal.

G4double	DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const
	Calculate exact shortest distance to any boundary from outside.

G4double	DistanceToIn (const G4ThreeVector &p) const
	Calculate exact shortest distance to any boundary from outside.

G4double	DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool validNorm=0, G4ThreeVector n=0) const
	Calculate exact shortest distance to any boundary from inside.

G4double	DistanceToOut (const G4ThreeVector &p) const
	Calculate exact shortest distance to any boundary from inside.

G4GeometryType	GetEntityType () const
	Get entity type.

G4ThreeVector	GetPointOnSurface () const
	get point on surface

G4VSolid *	Clone () const
	Make a clone of the object.

std::ostream &	StreamInfo (std::ostream &os) const
	Stream object contents to an output stream.

void	DescribeYourselfTo (G4VGraphicsScene &scene) const
	Visualisation functions.

G4Polyhedron *	CreatePolyhedron () const
	create polyhedron

void	BoundingLimits (G4ThreeVector &pMin, G4ThreeVector &pMax) const
	Two vectors define an axis-parallel bounding box for the shape.

	BelleCrystal (__void__ &)
	Fake default constructor for usage restricted to direct object persistency for clients requiring preallocation of memory for persistifiable objects.

	BelleCrystal (const BelleCrystal &rhs)
	copy constructor

BelleCrystal &	operator= (const BelleCrystal &rhs)
	assignment operator

G4ThreeVector	vertex (unsigned int i) const
	return ith vertex

Protected Member Functions
G4bool	MakePlane (const G4ThreeVector &p1, const G4ThreeVector &p2, const G4ThreeVector &p3, const G4ThreeVector &p4, Plane_t &plane) const
	Calculate the coef's of the plane p1->p2->p3->p4->p1 where the ThreeVectors 1-4 are in anti-clockwise order when viewed from in front of the plane (i.e.

Private Member Functions
G4ThreeVector	GetPointOnTriangle (int) const
	Returns a random point on the surface of one of the faces.

double	area (int, double &) const
	triangle area

double	getvolarea () const
	get volume area

const unsigned int *	ivertx (unsigned int i) const
	get the ith vertex

Private Attributes
unsigned int	nsides
	the number of sides

double	fDz
	Dz.

std::vector< Plane_t >	fPlanes
	vector of planes

std::vector< Point_t >	fx
	vector of points

std::vector< double >	fareas
	vector of area values

Detailed Description

a Belle crystal in Geant4

Definition at line 37 of file BelleCrystal.h.

Constructor & Destructor Documentation

◆ BelleCrystal() [1/4]

BelleCrystal ( const G4String & pName )

explicit

Constructor for "nominal" BelleCrystal.

Definition at line 145 of file BelleCrystal.cc.

  : G4CSGSolid(pName), nsides(0), fDz(0)
{
}

◆ BelleCrystal() [2/4]

BelleCrystal	(	const G4String &	pName,
		int	n,
		const G4ThreeVector *	pt
	)

Constructor.

Definition at line 47 of file BelleCrystal.cc.

  : G4CSGSolid(pName), nsides(n), fPlanes(nsides), fx(2 * nsides)
{
#if PERFCOUNTER==1
  memset(fcounter, 0, sizeof(fcounter));
#endif
#ifdef MATCHNAME
  // cppcheck-suppress ConfigurationNotChecked
  fmatch = GetName() == QUOTE(MATCHNAME);
  cout.precision(17);
#endif
  //  ref = new G4Trap("dummy",pt);
  fDz = abs(pt[2 * nsides - 1].z());
  for (unsigned int i = 0; i < 2 * nsides; i++) {
    fx[i].x = pt[i].x();
    fx[i].y = pt[i].y();
  }
  //  cout<<pName<<" "<<fx.data()<<endl;
  //  for(int i=0; i<2*nsides + 2*(nsides-2); i++) ivertx(i);
 
  auto isConvex = [](std::vector<Point_t>::const_iterator begin, std::vector<Point_t>::const_iterator end) -> bool {
    bool sign = false;
    int np = end - begin;
    for (int i0 = 0; i0 < np; i0++)
    {
      int i1 = (i0 + 1) % np, i2 = (i0 + 2) % np;
      const Point_t& r2 = *(begin + i2), &r1 = *(begin + i1), &r0 = *(begin + i0);
      double dx1 = r2.x - r1.x, dy1 = r2.y - r1.y;
      double dx2 = r0.x - r1.x, dy2 = r0.y - r1.y;
      double x = dx1 * dy2 - dy1 * dx2;
      if (i0 == 0) sign = x > 0;
      else {
        if (sign != (x > 0)) return false;
      }
    }
    return true;
  };
 
  auto isClockwise = [](std::vector<Point_t>::const_iterator begin, std::vector<Point_t>::const_iterator end) -> bool {
    std::vector<Point_t>::const_iterator it = begin;
    Point_t r0 = *it++;
    double sum = 0;
    for (; it != end;)
    {
      Point_t r1 = *it++; sum += (r1.x - r0.x) * (r1.y + r0.y);
      r0 = r1;
    }
    Point_t r1 = *begin; sum += (r1.x - r0.x) * (r1.y + r0.y);
    return sum > 0;
  };
 
  if (!isConvex(fx.begin(), fx.begin() + nsides)) {
    std::ostringstream message;
    message << "At -z polygon is not convex: " << GetName();
    G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
  }
  if (!isConvex(fx.begin() + nsides, fx.end())) {
    std::ostringstream message;
    message << "At +z polygon is not convex: " << GetName();
    G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
  }
 
  if (isClockwise(fx.begin(), fx.begin() + nsides)) {
    std::ostringstream message;
    message << "At -z polygon is not in anti-clockwise order: " << GetName();
    G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
  }
  if (isClockwise(fx.begin() + nsides, fx.end())) {
    std::ostringstream message;
    message << "At +z polygon is not in anti-clockwise order: " << GetName();
    G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
  }
 
  // sanity check
  const G4ThreeVector* pm = pt, *pp = pt + nsides;
  bool zm = pm[0].z() < 0, zp = pp[0].z() > 0;
  for (unsigned int i = 1; i < nsides; i++) {
    zm = zm && (pm[i - 1].z() - pm[i].z()) < kCarTolerance;
    zp = zp && (pp[i - 1].z() - pp[i].z()) < kCarTolerance;
  }
  bool zpm = abs(pm[0].z() + pp[0].z()) < kCarTolerance;
  if (!(zm && zp && zpm)) {
    std::ostringstream message;
    message << "Invalid vertice coordinates for Solid: " << GetName() << " " << zp << " " << zm << " " << zpm;
    G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
  }
  if (nsides > 3) {
    for (unsigned int i = 0; i < nsides; i++)
      MakePlane(pt[i], pt[i + nsides], pt[((i + 1) % (nsides)) + nsides], pt[(i + 1) % nsides], fPlanes[i]);
  } else {
    std::ostringstream message;
    message << "Wrong number of sides for Belle Crystal: " << GetName() << " nsides = " << nsides;
    G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
  }
}

◆ ~BelleCrystal()

~BelleCrystal ( )

virtual

Destructor.

Definition at line 158 of file BelleCrystal.cc.

{
#if PERFCOUNTER==1
  cout << GetName() << " ";
  for (int i = 0; i < 6; i++) cout << counter[i] << " "; cout << endl;
  exit(0);
#endif
  //  delete ref;
}

◆ BelleCrystal() [3/4]

BelleCrystal ( __void__ & a )

explicit

Fake default constructor for usage restricted to direct object persistency for clients requiring preallocation of memory for persistifiable objects.

Definition at line 152 of file BelleCrystal.cc.

  : G4CSGSolid(a), nsides(0), fDz(0)
{
}

◆ BelleCrystal() [4/4]

BelleCrystal ( const BelleCrystal & rhs )

copy constructor

Definition at line 169 of file BelleCrystal.cc.

  : G4CSGSolid(rhs), nsides(rhs.nsides), fDz(rhs.fDz), fPlanes(rhs.fPlanes), fx(rhs.fx)
{
}

Member Function Documentation

◆ area()

double area	(	int	it,
		double &	vol
	)		const

private

triangle area

Definition at line 648 of file BelleCrystal.cc.

{
  const unsigned int* iv = ivertx(it);
  G4ThreeVector p0 = vertex(iv[0]), p1 = vertex(iv[1]), p2 = vertex(iv[2]);
  //  std::cout<<it<<" "<<p0<<" "<<p1<<" "<<p2<<std::endl;
  G4ThreeVector n = (p1 - p0).cross(p2 - p0);
  vol = n * p0;
  return 0.5 * n.mag();
}

◆ BoundingLimits()

void BoundingLimits	(	G4ThreeVector &	pMin,
		G4ThreeVector &	pMax
	)		const

Two vectors define an axis-parallel bounding box for the shape.

Definition at line 733 of file BelleCrystal.cc.

{
  // All extremums will be at vertices
  int np = 2 * nsides;
  G4ThreeVector pt;
 
  // Placeholder vectors
  const double inf = std::numeric_limits<double>::infinity();
  G4ThreeVector minimum(inf, inf, inf), maximum(-inf, -inf, -inf);
 
  for (int i = 0; i < np; i++) {
    pt = vertex(i);
 
    // Assign new minimum and new maximum
    minimum = min(minimum, pt);
    maximum = max(maximum, pt);
  }
 
  // Assign
  pMin = minimum;
  pMax = maximum;
}

◆ CalculateExtent()

G4bool CalculateExtent	(	const EAxis	pAxis,
		const G4VoxelLimits &	pVoxelLimit,
		const G4AffineTransform &	pTransform,
		G4double &	pMin,
		G4double &	pMax
	)		const

calculate the extent of the volume

Definition at line 254 of file BelleCrystal.cc.

{
  const double inf = std::numeric_limits<double>::infinity();
  G4ThreeVector v0(inf, inf, inf), v1(-inf, -inf, -inf);
  for (unsigned int i = 0; i < 2 * nsides; i++) {
    G4ThreeVector v = pTransform.TransformPoint(vertex(i));
    v0 = min(v0, v);
    v1 = max(v1, v);
  }
  G4ThreeVector b0(bb.GetMinXExtent(), bb.GetMinYExtent(), bb.GetMinZExtent());
  G4ThreeVector b1(bb.GetMaxXExtent(), bb.GetMaxYExtent(), bb.GetMaxZExtent());
 
  v0 = max(v0, b0);
  v1 = min(v1, b1);
 
  switch (pAxis) {
    case kXAxis: pMin = v0.x(); pMax = v1.x(); break;
    case kYAxis: pMin = v0.y(); pMax = v1.y(); break;
    case kZAxis: pMin = v0.z(); pMax = v1.z(); break;
    default:                               break;
  }
 
  G4bool flag = false;
  if ((pMin < inf) || (pMax > -inf)) {
    flag = true;
    // Add tolerance to avoid precision troubles
    pMin -= kCarTolerance ;
    pMax += kCarTolerance ;
  }
  // do {
  //   double t_pMin, t_pMax;
  //   bool t_flag = ref->CalculateExtent(pAxis, bb, pTransform, t_pMin, t_pMax);
  //   if(flag!=t_flag||abs(t_pMin-pMin)>kCarTolerance||abs(t_pMax-pMax)>kCarTolerance){
  //     cout<<GetName()<<" "<<pAxis<<" "<<bb<<" "<<pTransform<<t_pMin-pMin<<" "<<t_pMax-pMax<<endl;
  //     exit(0);
  //   }
  // } while(0);
  return flag;
}

◆ Clone()

G4VSolid * Clone ( ) const

Make a clone of the object.

Definition at line 547 of file BelleCrystal.cc.

{
  return new BelleCrystal(*this);
}

◆ ComputeDimensions()

void ComputeDimensions	(	G4VPVParameterisation *	p,
		const G4int	n,
		const G4VPhysicalVolume *	pRep
	)

compute the dimensions

Definition at line 230 of file BelleCrystal.cc.

{
  G4Exception("BelleCrystal::ComputeDimensions()",
              "GeomSolids0001", FatalException,
              "BelleCrystal does not support Parameterisation.");
  // std::cout<<"ComputeDimensions"<<std::endl;
  // p->ComputeDimensions(*this,n,pRep);
}

◆ CreatePolyhedron()

G4Polyhedron * CreatePolyhedron ( ) const

create polyhedron

Definition at line 724 of file BelleCrystal.cc.

{
  int np = 2 * nsides;
  vector<G4ThreeVector> pt(np);
  for (int i = 0; i < np; i++) pt[i] = vertex(i);
  return new PolyhedronBelleCrystal(np, pt.data());
}

◆ DescribeYourselfTo()

void DescribeYourselfTo ( G4VGraphicsScene & scene ) const

Visualisation functions.

Definition at line 703 of file BelleCrystal.cc.

{
  scene.AddSolid(*this);
}

◆ DistanceToIn() [1/2]

G4double DistanceToIn ( const G4ThreeVector & p ) const

Calculate exact shortest distance to any boundary from outside.

Definition at line 382 of file BelleCrystal.cc.

{
  //  return ref->DistanceToIn(p);
  COUNTER(2);
  MATCHOUT("DistanceToIn(p) " << p.x() << " " << p.y() << " " << p.z());
  G4double d = std::fabs(p.z()) - fDz;
  unsigned int i = 0;
  do {
    const Plane_t& t = fPlanes[i++];
    d = max(t.n * p + t.d, d);
  } while (i < nsides);
  d = max(d, 0.0);
 
  // if(abs(d - ref->DistanceToIn(p))>kCarTolerance){
  //   cout<<GetName()<<" "<<p<<" "<<d-ref->DistanceToIn(p)<<endl;
  //   exit(0);
  // }
  //  if(fmatch&&d==0.0) cout<<"inside\n";
  return d;
}

◆ DistanceToIn() [2/2]

G4double DistanceToIn	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v
	)		const

Calculate exact shortest distance to any boundary from outside.

Definition at line 430 of file BelleCrystal.cc.

{
  //  return ref->DistanceToIn(p,v);
  COUNTER(4);
  MATCHOUT("DistanceToIn(p,v) " << p.x() << " " << p.y() << " " << p.z() << " " << v.x() << " " << v.y() << " " << v.z());
  G4double delta = 0.5 * kCarTolerance;
  double tin = 0, tout = kInfinity;
  auto outside = [delta, &tin, &tout](double d, double nv) -> bool {
    double t = -d / nv;
    if (nv < 0)
      tin  = max(tin, t);
    else
    {
      if (d >= -delta) return true;
      tout = min(tout, t);
    }
    return false;
  };
 
  unsigned int i = 0;
  do {
    const Plane_t& t = fPlanes[i++];
    if (outside(t.n * p + t.d, t.n * v)) {tin = kInfinity; goto exit;}
  } while (i < nsides);
  do {
    if (outside(p.z() - fDz, v.z())) {tin = kInfinity; goto exit;}
    if (outside(-p.z() - fDz, -v.z())) {tin = kInfinity; goto exit;}
  } while (0);
exit:
  double res = (tin > tout) ? kInfinity : tin;
  // G4int oldprc = cout.precision(16);
  // if(nsides==5) cout<<GetName()<<" "<<p.x()<<" "<<p.y()<<" "<<p.z()<<" "<<v.x()<<" "<<v.y()<<" "<<v.z()<<" "<<res<<endl;
  // cout.precision(oldprc);
 
  // if(abs(res - ref->DistanceToIn(p,v))>kCarTolerance){
  //   cout<<GetName()<<" "<<p<<" "<<v<<" "<<res-ref->DistanceToIn(p,v)<<endl;
  //   exit(0);
  // }
  return res;
}

◆ DistanceToOut() [1/2]

G4double DistanceToOut ( const G4ThreeVector & p ) const

Calculate exact shortest distance to any boundary from inside.

Definition at line 405 of file BelleCrystal.cc.

{
  //  return ref->DistanceToOut(p);
  COUNTER(3);
  MATCHOUT("DistanceToOut(p) " << p.x() << " " << p.y() << " " << p.z());
  G4double d = fDz - std::fabs(p.z());
  unsigned int i = 0;
  do {
    const Plane_t& t = fPlanes[i++];
    d = min(-(t.n * p + t.d), d);
  } while (i < nsides);
  d = max(d, 0.0);
  // if(abs(d - ref->DistanceToOut(p))>0){
  //   cout<<GetName()<<" "<<p<<" "<<d<<" "<<ref->DistanceToOut(p)<<" "<<fDz-ref->GetZHalfLength()<<endl;
  //   exit(0);
  // }
  return d;
}

◆ DistanceToOut() [2/2]

G4double DistanceToOut	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v,
		const G4bool	calcNorm = `false`,
		G4bool *	validNorm = `0`,
		G4ThreeVector *	n = `0`
	)		const

Calculate exact shortest distance to any boundary from inside.

Definition at line 474 of file BelleCrystal.cc.

{
  //  return ref->DistanceToOut(p,v,calcNorm,IsValid,n);
  COUNTER(5);
  MATCHOUT("DistanceToOut(p,v) " << p.x() << " " << p.y() << " " << p.z() << " " << v.x() << " " << v.y() << " " << v.z() << " " <<
           calcNorm);
  const G4double delta = 0.5 * kCarTolerance;
  unsigned int iside = 10;
  G4double lmin = kInfinity;
 
  auto outside = [delta, &lmin, &iside](double d, double nv, unsigned int i) -> bool {
    if (d > delta) // definitely outside
    {
      lmin = 0; iside = i; return true;
    } else
    {
      bool c = nv > 0;
      if (d < -delta) { // definitely inside
        if (c) { // has to point to the same direction
          d = -d;
          if (d < nv * lmin) {lmin = d / nv; iside = i;}
        }
      } else { // surface
        if (c) { // points outside
          lmin = 0; iside = i; return true;
        }
      }
    }
    return false;
  };
 
  unsigned int i = 0;
  do {
    const Plane_t& t = fPlanes[i];
    if (outside(t.n * p + t.d, t.n * v, i++)) goto exit;
  } while (i < nsides);
  do {
    if (outside(p.z() - fDz, v.z(), i++)) goto exit;
    if (outside(-p.z() - fDz, -v.z(), i++)) goto exit;
  } while (0);
exit:
  if (calcNorm) {
    *IsValid = true;
    if (iside < nsides) {
      *n = fPlanes[iside].n;
    } else if (iside == nsides) {
      *n = G4ThreeVector(0, 0, 1);
    } else {
      *n = G4ThreeVector(0, 0, -1);
    }
  }
  // G4int oldprc = cout.precision(16);
  // if(nsides==5) cout<<GetName()<<" "<<p.x()<<" "<<p.y()<<" "<<p.z()<<" "<<v.x()<<" "<<v.y()<<" "<<v.z()<<" "<<calcNorm<<" "<<lmin<<endl;
  // cout.precision(oldprc);
  // do {
  //   G4bool t_IsValid; G4ThreeVector t_n;
  //   double t_lmin = ref->DistanceToOut(p,v,calcNorm,&t_IsValid,&t_n);
  //   if(abs(lmin - t_lmin)>kCarTolerance||*IsValid!=t_IsValid||(*n - t_n).mag()>kCarTolerance){
  //     cout<<GetName()<<" "<<p<<" "<<v<<" "<<lmin-t_lmin<<(*n-t_n)<<endl;
  //     exit(0);
  //   }
  // } while(0);
  return lmin;
}

◆ GetCubicVolume()

G4double GetCubicVolume ( )

get the Cubic volume

Definition at line 674 of file BelleCrystal.cc.

{
  //  cout<<GetName()<<" GetCubicVolume "<<fCubicVolume<<endl;
  if (fCubicVolume == 0.0) {
    fCubicVolume = getvolarea();
    fSurfaceArea = fareas.back();
  }
  return fCubicVolume;
}

◆ GetEntityType()

G4GeometryType GetEntityType ( ) const

Get entity type.

Definition at line 541 of file BelleCrystal.cc.

{
  return G4String("BelleCrystal");
}

◆ GetPointOnSurface()

G4ThreeVector GetPointOnSurface ( ) const

get point on surface

Definition at line 694 of file BelleCrystal.cc.

{
  if (fareas.size() == 0) getvolarea();
  double r = CLHEP::RandFlat::shoot(0., fareas.back());
  std::vector<double>::const_iterator it = std::lower_bound(fareas.begin(), fareas.end(), r);
  return (it != fareas.end()) ? GetPointOnTriangle(it - fareas.begin()) : GetPointOnSurface();
}

◆ GetPointOnTriangle()

G4ThreeVector GetPointOnTriangle ( int it ) const

private

Returns a random point on the surface of one of the faces.

Definition at line 634 of file BelleCrystal.cc.

{
  // barycentric coordinates
  double a1 = CLHEP::RandFlat::shoot(0., 1.), a2 = CLHEP::RandFlat::shoot(0., 1.);
  if (a1 + a2 > 1) { a1 = 1 - a1; a2 = 1 - a2;}
  double a0 = 1 - (a1 + a2);
  const unsigned int* iv = ivertx(it);
  G4ThreeVector p0 = vertex(iv[0]), p1 = vertex(iv[1]), p2 = vertex(iv[2]);
  G4ThreeVector r1(p0.x()*a0 + p1.x()*a1 + p2.x()*a2,
                   p0.y()*a0 + p1.y()*a1 + p2.y()*a2,
                   p0.z()*a0 + p1.z()*a1 + p2.z()*a2);
  return r1;
}

◆ GetSidePlane()

Plane_t GetSidePlane ( G4int n ) const

inline

Get side plane.

Definition at line 48 of file BelleCrystal.h.

48{return fPlanes[n];}

◆ GetSurfaceArea()

G4double GetSurfaceArea ( )

get the surface area

Definition at line 684 of file BelleCrystal.cc.

{
  if (fSurfaceArea == 0.0) {
    fCubicVolume = getvolarea();
    fSurfaceArea = fareas.back();
  }
  return fSurfaceArea;
}

◆ getvolarea()

double getvolarea ( ) const

private

get volume area

Definition at line 658 of file BelleCrystal.cc.

{
  fareas.clear();
  int nt = 2 * nsides + 2 * (nsides - 2); // total number of triangles
  fareas.reserve(nt);
  double totarea = 0, totvol = 0;
  for (int i = 0; i < nt; i++) {
    double v, s = area(i, v);
    totarea += s;
    totvol  += v;
    fareas.push_back(totarea);
    //    cout<<i<<" "<<s<<" "<<v<<endl;
  }
  return totvol / 6;
}

◆ Inside()

EInside Inside ( const G4ThreeVector & p ) const

Return whether point inside/outside/on surface, using tolerance.

Definition at line 298 of file BelleCrystal.cc.

{
  COUNTER(0);
  MATCHOUT("Inside(p) " << p.x() << " " << p.y() << " " << p.z());
  //  return ref->Inside(p);
  double d = std::fabs(p.z()) - fDz;
  unsigned int i = 0;
  do {
    const Plane_t& t = fPlanes[i++];
    d = max(t.n * p + t.d, d);
  } while (i < nsides);
  const G4double delta = 0.5 * kCarTolerance;
  int in = 0;
  in += d <= delta;
  in += d <= -delta;
 
  EInside res = EInside(in);
  // if(res != ref->Inside(p)){
  //   cout<<GetName()<<" "<<p<<endl;
  //   exit(0);
  // }
  return res;
}

◆ ivertx()

const unsigned int * ivertx ( unsigned int i ) const

private

get the ith vertex

Definition at line 572 of file BelleCrystal.cc.

{
  static unsigned int vert[3];
  //  if(nsides==4) return ivertx4[it];
  if (nsides == 0) return ivertx4[it];
  //  else if(nsides==5) { vert[0] = ivertx5[it][0], vert[1] = ivertx5[it][1], vert[2] = ivertx5[it][2];}
  else {
    if (it < 2 * nsides) {
      int j = it / 2;
      if ((it & 1) == 0) {
        vert[0] = j, vert[2] = j + nsides, vert[1] = ((j + 1) % nsides) + nsides;
      } else {
        vert[0] = j, vert[2] = ((j + 1) % nsides) + nsides, vert[1] = (j + 1) % nsides;
      }
    } else if (it < 2 * nsides + (nsides - 2)) {
      int j = it - 2 * nsides;
      vert[0] = 0, vert[2] = 1 + j, vert[1] = 2 + j;
    } else {
      int j = it - (2 * nsides + (nsides - 2));
      vert[0] = nsides, vert[2] = nsides + 2 + j, vert[1] = nsides + 1 + j;
    }
  }
  //  cout<<it<<" {"<<vert[0]<<", "<<vert[1]<<", "<<vert[2]<<"}"<<endl;
  return vert;
}

◆ MakePlane()

G4bool MakePlane	(	const G4ThreeVector &	p1,
		const G4ThreeVector &	p2,
		const G4ThreeVector &	p3,
		const G4ThreeVector &	p4,
		Plane_t &	plane
	)		const

protected

Calculate the coef's of the plane p1->p2->p3->p4->p1 where the ThreeVectors 1-4 are in anti-clockwise order when viewed from in front of the plane (i.e.

from normal direction).

Return true if the ThreeVectors are coplanar + set coef;s false if ThreeVectors are not coplanar

Definition at line 198 of file BelleCrystal.cc.

{
  G4ThreeVector v12    = p2 - p1;
  G4ThreeVector v13    = p3 - p1;
  G4ThreeVector v14    = p4 - p1;
  G4ThreeVector Vcross = v12.cross(v13);
 
  if (std::fabs(Vcross.dot(v14) / (Vcross.mag()*v14.mag())) > kCoplanar_Tolerance) {
    std::ostringstream message;
    message << "Vertices are not in the same plane: " << GetName() << " volume =" << Vcross.dot(v14);
    G4Exception("BelleCrystal::MakePlane()", "BelleCrystal", FatalException, message);
  }
 
  double a = +(p4.y() - p2.y()) * (p3.z() - p1.z()) - (p3.y() - p1.y()) * (p4.z() - p2.z());
  double b = -(p4.x() - p2.x()) * (p3.z() - p1.z()) + (p3.x() - p1.x()) * (p4.z() - p2.z());
  double c = +(p4.x() - p2.x()) * (p3.y() - p1.y()) - (p3.x() - p1.x()) * (p4.y() - p2.y());
  double sd = sqrt(a * a + b * b + c * c); // so now vector plane.(a,b,c) is unit
  a /= sd;
  b /= sd;
  c /= sd;
  plane.n = G4ThreeVector(a, b, c);
  plane.d = -(a * p1.x() + b * p1.y() + c * p1.z());
 
  // plane.n = (p4-p2).cross(p3-p1).unit();
  // plane.d =-(plane.n*p1);
  return true;
}

◆ operator=()

BelleCrystal & operator= ( const BelleCrystal & rhs )

assignment operator

Definition at line 175 of file BelleCrystal.cc.

{
  // Check assignment to self
  if (this == &rhs)  { return *this; }
 
  // Copy base class data
  G4CSGSolid::operator=(rhs);
 
  // Copy data
  nsides = rhs.nsides;
  fDz = rhs.fDz;
  fx = rhs.fx;
  fPlanes = rhs.fPlanes;
 
  return *this;
}

◆ StreamInfo()

std::ostream & StreamInfo ( std::ostream & os ) const

Stream object contents to an output stream.

Definition at line 553 of file BelleCrystal.cc.

{
  G4int oldprc = os.precision(16);
  os << "-----------------------------------------------------------\n"
     << "    *** Dump for solid - " << GetName() << " ***\n"
     << "    ===================================================\n"
     << " Solid type: BelleCrystal\n"
     << " Parameters: \n"
     << "    crystal side plane equations:\n"
     << "-----------------------------------------------------------\n";
  os.precision(oldprc);
 
  return os;
}

◆ SurfaceNormal()

G4ThreeVector SurfaceNormal ( const G4ThreeVector & p ) const

Calculate side nearest to p, and return normal.

Definition at line 324 of file BelleCrystal.cc.

{
  //  return ref->SurfaceNormal(p);
  COUNTER(1);
  MATCHOUT("SurfaceNormal(p) " << p.x() << " " << p.y() << " " << p.z());
  const G4double delta = 0.5 * kCarTolerance;
  G4double safe = kInfinity;
  unsigned int iside = 0, kside = 0;
  vector<double> adist(nsides + 2);
  auto dist = [delta, &safe, &iside, &kside, &adist](double d, unsigned int i) -> void {
    d = std::abs(d);
    adist[i] = d;
    iside = (d < safe) ? i : iside;
    safe = min(d, safe);
    kside += d < delta;
  };
 
  unsigned int i = 0;
  do {
    const Plane_t& t = fPlanes[i];
    dist(t.n * p + t.d, i++);
  } while (i < nsides);
  do {
    dist(p.z() - fDz, i++);
    dist(p.z() + fDz, i++);
  } while (0);
 
  G4ThreeVector res;
  if (kside > 1) {
    i = 0;
    do {
      if (adist[i++] < delta) res += fPlanes[i].n;
    } while (i < nsides);
    do {
      if (adist[i++] < delta) res.setZ(res.z() + 1.0);
      if (adist[i++] < delta) res.setZ(res.z() - 1.0);
    } while (0);
    res = res.unit();
  } else {
    if (iside < nsides) {
      res = fPlanes[iside].n;
    } else if (iside == nsides) {
      res = G4ThreeVector(0, 0, 1);
    } else {
      res = G4ThreeVector(0, 0, -1);
    }
  }
  // if((res- ref->SurfaceNormal(p)).mag()>kCarTolerance){
  //   cout<<GetName()<<" "<<p<<" "<<res-ref->SurfaceNormal(p)<<endl;
  //   exit(0);
  // }
 
  return res;
}

◆ vertex()

G4ThreeVector vertex ( unsigned int i ) const

return ith vertex

Definition at line 612 of file BelleCrystal.cc.

{
  double Dz = (i < nsides) ? -fDz : fDz;
 
  // if(nsides==4) {
  //   const unsigned char w = 2, mask = ((1<<w)-1);
  //   unsigned int map = PACK(0,2,w,0)+PACK(0,3,w,1)+PACK(1,2,w,2)+PACK(1,3,w,3); // 20 bits
  //   map >>= (2*w)*(i%4);
  //   int k0 = map&mask, k1 = (map>>2)&mask;
  //   return crossplanes(fPlanes[k0],fPlanes[k1],Dz);
  // } else if(nsides==5){
  //   const unsigned char w = 3, mask = ((1<<w)-1);
  //   unsigned int map = PACK(0,4,w,0)+PACK(1,0,w,1)+PACK(2,1,w,2)+PACK(3,2,w,3)+PACK(4,3,w,4); // 30 bits
  //   map >>= (2*w)*(i%5);
  //   int k0 = map&mask, k1 = (map>>w)&mask;
  //   return crossplanes(fPlanes[k0],fPlanes[k1],Dz);
  // }
  return G4ThreeVector(fx[i].x, fx[i].y, Dz);
}

Member Data Documentation

◆ fareas

std::vector<double> fareas

mutableprivate

vector of area values

Definition at line 149 of file BelleCrystal.h.

◆ fDz

double fDz

private

Dz.

Definition at line 145 of file BelleCrystal.h.

◆ fPlanes

std::vector<Plane_t> fPlanes

private

vector of planes

Definition at line 146 of file BelleCrystal.h.

◆ fx

std::vector<Point_t> fx

private

vector of points

Definition at line 147 of file BelleCrystal.h.

◆ nsides

unsigned int nsides

private

the number of sides

Definition at line 144 of file BelleCrystal.h.

The documentation for this class was generated from the following files:

ecl/geometry/include/BelleCrystal.h
ecl/geometry/src/BelleCrystal.cc

Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ BelleCrystal() [1/4]

◆ BelleCrystal() [2/4]

◆ ~BelleCrystal()

◆ BelleCrystal() [3/4]

◆ BelleCrystal() [4/4]

Member Function Documentation

◆ area()

◆ BoundingLimits()

◆ CalculateExtent()

◆ Clone()

◆ ComputeDimensions()

◆ CreatePolyhedron()

◆ DescribeYourselfTo()

◆ DistanceToIn() [1/2]

◆ DistanceToIn() [2/2]

◆ DistanceToOut() [1/2]

◆ DistanceToOut() [2/2]

◆ GetCubicVolume()

◆ GetEntityType()

◆ GetPointOnSurface()

◆ GetPointOnTriangle()

◆ GetSidePlane()

◆ GetSurfaceArea()

◆ getvolarea()

◆ Inside()

◆ ivertx()

◆ MakePlane()

◆ operator=()

◆ StreamInfo()

◆ SurfaceNormal()

◆ vertex()

Member Data Documentation

◆ fareas

◆ fDz

◆ fPlanes

◆ fx

◆ nsides