Geant4_10
G4Tet.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the intellectual property of the *
19 // * Vanderbilt University Free Electron Laser Center *
20 // * Vanderbilt University, Nashville, TN, USA *
21 // * Development supported by: *
22 // * United States MFEL program under grant FA9550-04-1-0045 *
23 // * and NASA under contract number NNG04CT05P *
24 // * Written by Marcus H. Mendenhall and Robert A. Weller. *
25 // * *
26 // * Contributed to the Geant4 Core, January, 2005. *
27 // * *
28 // ********************************************************************
29 //
30 // $Id: G4Tet.cc 76263 2013-11-08 11:41:52Z gcosmo $
31 //
32 // class G4Tet
33 //
34 // Implementation for G4Tet class
35 //
36 // History:
37 //
38 // 20040903 - Marcus Mendenhall, created G4Tet
39 // 20041101 - Marcus Mendenhall, optimized constant dot products with
40 // fCdotNijk values
41 // 20041101 - MHM removed tracking error by clipping DistanceToOut to 0
42 // for surface cases
43 // 20041101 - MHM many speed optimizations in if statements
44 // 20041101 - MHM changed vdotn comparisons to 1e-12 instead of 0.0 to
45 // avoid nearly-parallel problems
46 // 20041102 - MHM Added extra distance into solid to DistanceToIn(p,v)
47 // hit testing
48 // 20041102 - MHM added ability to check for degeneracy without throwing
49 // G4Exception
50 // 20041103 - MHM removed many unused variables from class
51 // 20040803 - Dionysios Anninos, added GetPointOnSurface() method
52 // 20061112 - MHM added code for G4VSolid GetSurfaceArea()
53 // 20100920 - Gabriele Cosmo added copy-ctor and operator=()
54 //
55 // --------------------------------------------------------------------
56 
57 #include "G4Tet.hh"
58 
59 #if !defined(G4GEOM_USE_UTET)
60 
61 const char G4Tet::CVSVers[]="$Id: G4Tet.cc 76263 2013-11-08 11:41:52Z gcosmo $";
62 
63 #include "G4VoxelLimits.hh"
64 #include "G4AffineTransform.hh"
65 
66 #include "G4VPVParameterisation.hh"
67 
68 #include "Randomize.hh"
69 
70 #include "G4VGraphicsScene.hh"
71 #include "G4Polyhedron.hh"
72 #include "G4VisExtent.hh"
73 
74 #include "G4ThreeVector.hh"
75 
76 #include <cmath>
77 
78 using namespace CLHEP;
79 
81 //
82 // Constructor - create a tetrahedron
83 // This class is implemented separately from general polyhedra,
84 // because the simplex geometry can be computed very quickly,
85 // which may become important in situations imported from mesh generators,
86 // in which a very large number of G4Tets are created.
87 // A Tet has all of its geometrical information precomputed
88 
89 G4Tet::G4Tet(const G4String& pName,
90  G4ThreeVector anchor,
91  G4ThreeVector p2,
92  G4ThreeVector p3,
93  G4ThreeVector p4, G4bool *degeneracyFlag)
94  : G4VSolid(pName), fpPolyhedron(0), warningFlag(0)
95 {
96  // fV<x><y> is vector from vertex <y> to vertex <x>
97  //
98  G4ThreeVector fV21=p2-anchor;
99  G4ThreeVector fV31=p3-anchor;
100  G4ThreeVector fV41=p4-anchor;
101 
102  // make sure this is a correctly oriented set of points for the tetrahedron
103  //
104  G4double signed_vol=fV21.cross(fV31).dot(fV41);
105  if(signed_vol<0.0)
106  {
107  G4ThreeVector temp(p4);
108  p4=p3;
109  p3=temp;
110  temp=fV41;
111  fV41=fV31;
112  fV31=temp;
113  }
114  fCubicVolume = std::fabs(signed_vol) / 6.;
115 
116  G4ThreeVector fV24=p2-p4;
117  G4ThreeVector fV43=p4-p3;
118  G4ThreeVector fV32=p3-p2;
119 
120  fXMin=std::min(std::min(std::min(anchor.x(), p2.x()),p3.x()),p4.x());
121  fXMax=std::max(std::max(std::max(anchor.x(), p2.x()),p3.x()),p4.x());
122  fYMin=std::min(std::min(std::min(anchor.y(), p2.y()),p3.y()),p4.y());
123  fYMax=std::max(std::max(std::max(anchor.y(), p2.y()),p3.y()),p4.y());
124  fZMin=std::min(std::min(std::min(anchor.z(), p2.z()),p3.z()),p4.z());
125  fZMax=std::max(std::max(std::max(anchor.z(), p2.z()),p3.z()),p4.z());
126 
127  fDx=(fXMax-fXMin)*0.5; fDy=(fYMax-fYMin)*0.5; fDz=(fZMax-fZMin)*0.5;
128 
129  fMiddle=G4ThreeVector(fXMax+fXMin, fYMax+fYMin, fZMax+fZMin)*0.5;
130  fMaxSize=std::max(std::max(std::max((anchor-fMiddle).mag(),
131  (p2-fMiddle).mag()),
132  (p3-fMiddle).mag()),
133  (p4-fMiddle).mag());
134 
135  G4bool degenerate=std::fabs(signed_vol) < 1e-9*fMaxSize*fMaxSize*fMaxSize;
136 
137  if(degeneracyFlag) *degeneracyFlag=degenerate;
138  else if (degenerate)
139  {
140  G4Exception("G4Tet::G4Tet()", "GeomSolids0002", FatalException,
141  "Degenerate tetrahedron not allowed.");
142  }
143 
144  fTol=1e-9*(std::fabs(fXMin)+std::fabs(fXMax)+std::fabs(fYMin)
145  +std::fabs(fYMax)+std::fabs(fZMin)+std::fabs(fZMax));
146  //fTol=kCarTolerance;
147 
148  fAnchor=anchor;
149  fP2=p2;
150  fP3=p3;
151  fP4=p4;
152 
153  G4ThreeVector fCenter123=(anchor+p2+p3)*(1.0/3.0); // face center
154  G4ThreeVector fCenter134=(anchor+p4+p3)*(1.0/3.0);
155  G4ThreeVector fCenter142=(anchor+p4+p2)*(1.0/3.0);
156  G4ThreeVector fCenter234=(p2+p3+p4)*(1.0/3.0);
157 
158  // compute area of each triangular face by cross product
159  // and sum for total surface area
160 
161  G4ThreeVector normal123=fV31.cross(fV21);
162  G4ThreeVector normal134=fV41.cross(fV31);
163  G4ThreeVector normal142=fV21.cross(fV41);
164  G4ThreeVector normal234=fV32.cross(fV43);
165 
166  fSurfaceArea=(
167  normal123.mag()+
168  normal134.mag()+
169  normal142.mag()+
170  normal234.mag()
171  )/2.0;
172 
173  fNormal123=normal123.unit();
174  fNormal134=normal134.unit();
175  fNormal142=normal142.unit();
176  fNormal234=normal234.unit();
177 
178  fCdotN123=fCenter123.dot(fNormal123);
179  fCdotN134=fCenter134.dot(fNormal134);
180  fCdotN142=fCenter142.dot(fNormal142);
181  fCdotN234=fCenter234.dot(fNormal234);
182 }
183 
185 //
186 // Fake default constructor - sets only member data and allocates memory
187 // for usage restricted to object persistency.
188 //
189 G4Tet::G4Tet( __void__& a )
190  : G4VSolid(a), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0),
191  fAnchor(0,0,0), fP2(0,0,0), fP3(0,0,0), fP4(0,0,0), fMiddle(0,0,0),
192  fNormal123(0,0,0), fNormal142(0,0,0), fNormal134(0,0,0),
193  fNormal234(0,0,0), warningFlag(0),
194  fCdotN123(0.), fCdotN142(0.), fCdotN134(0.), fCdotN234(0.),
195  fXMin(0.), fXMax(0.), fYMin(0.), fYMax(0.), fZMin(0.), fZMax(0.),
196  fDx(0.), fDy(0.), fDz(0.), fTol(0.), fMaxSize(0.)
197 {
198 }
199 
201 //
202 // Destructor
203 
205 {
206  delete fpPolyhedron;
207 }
208 
210 //
211 // Copy constructor
212 
213 G4Tet::G4Tet(const G4Tet& rhs)
214  : G4VSolid(rhs),
215  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
216  fpPolyhedron(0), fAnchor(rhs.fAnchor),
217  fP2(rhs.fP2), fP3(rhs.fP3), fP4(rhs.fP4), fMiddle(rhs.fMiddle),
218  fNormal123(rhs.fNormal123), fNormal142(rhs.fNormal142),
219  fNormal134(rhs.fNormal134), fNormal234(rhs.fNormal234),
220  warningFlag(rhs.warningFlag), fCdotN123(rhs.fCdotN123),
221  fCdotN142(rhs.fCdotN142), fCdotN134(rhs.fCdotN134),
222  fCdotN234(rhs.fCdotN234), fXMin(rhs.fXMin), fXMax(rhs.fXMax),
223  fYMin(rhs.fYMin), fYMax(rhs.fYMax), fZMin(rhs.fZMin), fZMax(rhs.fZMax),
224  fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), fTol(rhs.fTol),
225  fMaxSize(rhs.fMaxSize)
226 {
227 }
228 
229 
231 //
232 // Assignment operator
233 
235 {
236  // Check assignment to self
237  //
238  if (this == &rhs) { return *this; }
239 
240  // Copy base class data
241  //
242  G4VSolid::operator=(rhs);
243 
244  // Copy data
245  //
246  fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
247  fpPolyhedron = 0; fAnchor = rhs.fAnchor;
248  fP2 = rhs.fP2; fP3 = rhs.fP3; fP4 = rhs.fP4; fMiddle = rhs.fMiddle;
249  fNormal123 = rhs.fNormal123; fNormal142 = rhs.fNormal142;
250  fNormal134 = rhs.fNormal134; fNormal234 = rhs.fNormal234;
251  warningFlag = rhs.warningFlag; fCdotN123 = rhs.fCdotN123;
252  fCdotN142 = rhs.fCdotN142; fCdotN134 = rhs.fCdotN134;
253  fCdotN234 = rhs.fCdotN234; fXMin = rhs.fXMin; fXMax = rhs.fXMax;
254  fYMin = rhs.fYMin; fYMax = rhs.fYMax; fZMin = rhs.fZMin; fZMax = rhs.fZMax;
255  fDx = rhs.fDx; fDy = rhs.fDy; fDz = rhs.fDz; fTol = rhs.fTol;
256  fMaxSize = rhs.fMaxSize;
257 
258  return *this;
259 }
260 
262 //
263 // CheckDegeneracy
264 
266  G4ThreeVector p2,
267  G4ThreeVector p3,
268  G4ThreeVector p4 )
269 {
270  G4bool result;
271  G4Tet *object=new G4Tet("temp",anchor,p2,p3,p4,&result);
272  delete object;
273  return result;
274 }
275 
277 //
278 // Dispatch to parameterisation for replication mechanism dimension
279 // computation & modification.
280 
282  const G4int ,
283  const G4VPhysicalVolume* )
284 {
285 }
286 
288 //
289 // Calculate extent under transform and specified limit
290 
292  const G4VoxelLimits& pVoxelLimit,
293  const G4AffineTransform& pTransform,
294  G4double& pMin, G4double& pMax) const
295 {
296  G4double xMin,xMax;
297  G4double yMin,yMax;
298  G4double zMin,zMax;
299 
300  if (pTransform.IsRotated())
301  {
302  G4ThreeVector pp0=pTransform.TransformPoint(fAnchor);
303  G4ThreeVector pp1=pTransform.TransformPoint(fP2);
304  G4ThreeVector pp2=pTransform.TransformPoint(fP3);
305  G4ThreeVector pp3=pTransform.TransformPoint(fP4);
306 
307  xMin = std::min(std::min(std::min(pp0.x(), pp1.x()),pp2.x()),pp3.x());
308  xMax = std::max(std::max(std::max(pp0.x(), pp1.x()),pp2.x()),pp3.x());
309  yMin = std::min(std::min(std::min(pp0.y(), pp1.y()),pp2.y()),pp3.y());
310  yMax = std::max(std::max(std::max(pp0.y(), pp1.y()),pp2.y()),pp3.y());
311  zMin = std::min(std::min(std::min(pp0.z(), pp1.z()),pp2.z()),pp3.z());
312  zMax = std::max(std::max(std::max(pp0.z(), pp1.z()),pp2.z()),pp3.z());
313 
314  }
315  else
316  {
317  G4double xoffset = pTransform.NetTranslation().x() ;
318  xMin = xoffset + fXMin;
319  xMax = xoffset + fXMax;
320  G4double yoffset = pTransform.NetTranslation().y() ;
321  yMin = yoffset + fYMin;
322  yMax = yoffset + fYMax;
323  G4double zoffset = pTransform.NetTranslation().z() ;
324  zMin = zoffset + fZMin;
325  zMax = zoffset + fZMax;
326  }
327 
328  if (pVoxelLimit.IsXLimited())
329  {
330  if ( (xMin > pVoxelLimit.GetMaxXExtent()+fTol) ||
331  (xMax < pVoxelLimit.GetMinXExtent()-fTol) ) { return false; }
332  else
333  {
334  xMin = std::max(xMin, pVoxelLimit.GetMinXExtent());
335  xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent());
336  }
337  }
338 
339  if (pVoxelLimit.IsYLimited())
340  {
341  if ( (yMin > pVoxelLimit.GetMaxYExtent()+fTol) ||
342  (yMax < pVoxelLimit.GetMinYExtent()-fTol) ) { return false; }
343  else
344  {
345  yMin = std::max(yMin, pVoxelLimit.GetMinYExtent());
346  yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent());
347  }
348  }
349 
350  if (pVoxelLimit.IsZLimited())
351  {
352  if ( (zMin > pVoxelLimit.GetMaxZExtent()+fTol) ||
353  (zMax < pVoxelLimit.GetMinZExtent()-fTol) ) { return false; }
354  else
355  {
356  zMin = std::max(zMin, pVoxelLimit.GetMinZExtent());
357  zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent());
358  }
359  }
360 
361  switch (pAxis)
362  {
363  case kXAxis:
364  pMin=xMin;
365  pMax=xMax;
366  break;
367  case kYAxis:
368  pMin=yMin;
369  pMax=yMax;
370  break;
371  case kZAxis:
372  pMin=zMin;
373  pMax=zMax;
374  break;
375  default:
376  break;
377  }
378 
379  return true;
380 }
381 
383 //
384 // Return whether point inside/outside/on surface, using tolerance
385 
387 {
388  G4double r123, r134, r142, r234;
389 
390  // this is written to allow if-statement truncation so the outside test
391  // (where most of the world is) can fail very quickly and efficiently
392 
393  if ( (r123=p.dot(fNormal123)-fCdotN123) > fTol ||
394  (r134=p.dot(fNormal134)-fCdotN134) > fTol ||
395  (r142=p.dot(fNormal142)-fCdotN142) > fTol ||
396  (r234=p.dot(fNormal234)-fCdotN234) > fTol )
397  {
398  return kOutside; // at least one is out!
399  }
400  else if( (r123 < -fTol)&&(r134 < -fTol)&&(r142 < -fTol)&&(r234 < -fTol) )
401  {
402  return kInside; // all are definitively inside
403  }
404  else
405  {
406  return kSurface; // too close to tell
407  }
408 }
409 
411 //
412 // Calculate side nearest to p, and return normal
413 // If two sides are equidistant, normal of first side (x/y/z)
414 // encountered returned.
415 // This assumes that we are looking from the inside!
416 
418 {
419  G4double r123=std::fabs(p.dot(fNormal123)-fCdotN123);
420  G4double r134=std::fabs(p.dot(fNormal134)-fCdotN134);
421  G4double r142=std::fabs(p.dot(fNormal142)-fCdotN142);
422  G4double r234=std::fabs(p.dot(fNormal234)-fCdotN234);
423 
424  const G4double delta = 0.5*kCarTolerance;
425  G4ThreeVector sumnorm(0., 0., 0.);
426  G4int noSurfaces=0;
427 
428  if (r123 <= delta)
429  {
430  noSurfaces ++;
431  sumnorm= fNormal123;
432  }
433 
434  if (r134 <= delta)
435  {
436  noSurfaces ++;
437  sumnorm += fNormal134;
438  }
439 
440  if (r142 <= delta)
441  {
442  noSurfaces ++;
443  sumnorm += fNormal142;
444  }
445  if (r234 <= delta)
446  {
447  noSurfaces ++;
448  sumnorm += fNormal234;
449  }
450 
451  if( noSurfaces > 0 )
452  {
453  if( noSurfaces == 1 )
454  {
455  return sumnorm;
456  }
457  else
458  {
459  return sumnorm.unit();
460  }
461  }
462  else // Approximative Surface Normal
463  {
464 
465  if( (r123<=r134) && (r123<=r142) && (r123<=r234) ) { return fNormal123; }
466  else if ( (r134<=r142) && (r134<=r234) ) { return fNormal134; }
467  else if (r142 <= r234) { return fNormal142; }
468  return fNormal234;
469  }
470 }
472 //
473 // Calculate distance to box from an outside point
474 // - return kInfinity if no intersection.
475 // All this is very unrolled, for speed.
476 
478  const G4ThreeVector& v) const
479 {
480  G4ThreeVector vu(v.unit()), hp;
481  G4double vdotn, t, tmin=kInfinity;
482 
483  G4double extraDistance=10.0*fTol; // a little ways into the solid
484 
485  vdotn=-vu.dot(fNormal123);
486  if(vdotn > 1e-12)
487  { // this is a candidate face, since it is pointing at us
488  t=(p.dot(fNormal123)-fCdotN123)/vdotn; // # distance to intersection
489  if( (t>=-fTol) && (t<tmin) )
490  { // if not true, we're going away from this face or it's not close
491  hp=p+vu*(t+extraDistance); // a little beyond point of intersection
492  if ( ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
493  ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) &&
494  ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
495  {
496  tmin=t;
497  }
498  }
499  }
500 
501  vdotn=-vu.dot(fNormal134);
502  if(vdotn > 1e-12)
503  { // # this is a candidate face, since it is pointing at us
504  t=(p.dot(fNormal134)-fCdotN134)/vdotn; // # distance to intersection
505  if( (t>=-fTol) && (t<tmin) )
506  { // if not true, we're going away from this face
507  hp=p+vu*(t+extraDistance); // a little beyond point of intersection
508  if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
509  ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) &&
510  ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
511  {
512  tmin=t;
513  }
514  }
515  }
516 
517  vdotn=-vu.dot(fNormal142);
518  if(vdotn > 1e-12)
519  { // # this is a candidate face, since it is pointing at us
520  t=(p.dot(fNormal142)-fCdotN142)/vdotn; // # distance to intersection
521  if( (t>=-fTol) && (t<tmin) )
522  { // if not true, we're going away from this face
523  hp=p+vu*(t+extraDistance); // a little beyond point of intersection
524  if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
525  ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
526  ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
527  {
528  tmin=t;
529  }
530  }
531  }
532 
533  vdotn=-vu.dot(fNormal234);
534  if(vdotn > 1e-12)
535  { // # this is a candidate face, since it is pointing at us
536  t=(p.dot(fNormal234)-fCdotN234)/vdotn; // # distance to intersection
537  if( (t>=-fTol) && (t<tmin) )
538  { // if not true, we're going away from this face
539  hp=p+vu*(t+extraDistance); // a little beyond point of intersection
540  if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
541  ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
542  ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) )
543  {
544  tmin=t;
545  }
546  }
547  }
548 
549  return std::max(0.0,tmin);
550 }
551 
553 //
554 // Approximate distance to tet.
555 // returns distance to sphere centered on bounding box
556 // - If inside return 0
557 
559 {
560  G4double dd=(p-fMiddle).mag() - fMaxSize - fTol;
561  return std::max(0.0, dd);
562 }
563 
565 //
566 // Calcluate distance to surface of box from inside
567 // by calculating distances to box's x/y/z planes.
568 // Smallest distance is exact distance to exiting.
569 
571  const G4bool calcNorm,
572  G4bool *validNorm, G4ThreeVector *n) const
573 {
574  G4ThreeVector vu(v.unit());
575  G4double t1=kInfinity,t2=kInfinity,t3=kInfinity,t4=kInfinity, vdotn, tt;
576 
577  vdotn=vu.dot(fNormal123);
578  if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
579  {
580  t1=(fCdotN123-p.dot(fNormal123))/vdotn; // # distance to intersection
581  }
582 
583  vdotn=vu.dot(fNormal134);
584  if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
585  {
586  t2=(fCdotN134-p.dot(fNormal134))/vdotn; // # distance to intersection
587  }
588 
589  vdotn=vu.dot(fNormal142);
590  if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
591  {
592  t3=(fCdotN142-p.dot(fNormal142))/vdotn; // # distance to intersection
593  }
594 
595  vdotn=vu.dot(fNormal234);
596  if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
597  {
598  t4=(fCdotN234-p.dot(fNormal234))/vdotn; // # distance to intersection
599  }
600 
601  tt=std::min(std::min(std::min(t1,t2),t3),t4);
602 
603  if (warningFlag && (tt == kInfinity || tt < -fTol))
604  {
605  DumpInfo();
606  std::ostringstream message;
607  message << "No good intersection found or already outside!?" << G4endl
608  << "p = " << p / mm << "mm" << G4endl
609  << "v = " << v << G4endl
610  << "t1, t2, t3, t4 (mm) "
611  << t1/mm << ", " << t2/mm << ", " << t3/mm << ", " << t4/mm;
612  G4Exception("G4Tet::DistanceToOut(p,v,...)", "GeomSolids1002",
613  JustWarning, message);
614  if(validNorm)
615  {
616  *validNorm=false; // flag normal as meaningless
617  }
618  }
619  else if(calcNorm && n)
620  {
621  G4ThreeVector normal;
622  if(tt==t1) { normal=fNormal123; }
623  else if (tt==t2) { normal=fNormal134; }
624  else if (tt==t3) { normal=fNormal142; }
625  else if (tt==t4) { normal=fNormal234; }
626  *n=normal;
627  if(validNorm) { *validNorm=true; }
628  }
629 
630  return std::max(tt,0.0); // avoid tt<0.0 by a tiny bit
631  // if we are right on a face
632 }
633 
635 //
636 // Calculate exact shortest distance to any boundary from inside
637 // - If outside return 0
638 
640 {
641  G4double t1,t2,t3,t4;
642  t1=fCdotN123-p.dot(fNormal123); // distance to plane, positive if inside
643  t2=fCdotN134-p.dot(fNormal134); // distance to plane
644  t3=fCdotN142-p.dot(fNormal142); // distance to plane
645  t4=fCdotN234-p.dot(fNormal234); // distance to plane
646 
647  // if any one of these is negative, we are outside,
648  // so return zero in that case
649 
650  G4double tmin=std::min(std::min(std::min(t1,t2),t3),t4);
651  return (tmin < fTol)? 0:tmin;
652 }
653 
655 //
656 // Create a List containing the transformed vertices
657 // Note: Caller has deletion responsibility
658 
661 {
662  G4ThreeVectorList* vertices = new G4ThreeVectorList();
663 
664  if (vertices)
665  {
666  vertices->reserve(4);
667  G4ThreeVector vertex0(fAnchor);
668  G4ThreeVector vertex1(fP2);
669  G4ThreeVector vertex2(fP3);
670  G4ThreeVector vertex3(fP4);
671 
672  vertices->push_back(pTransform.TransformPoint(vertex0));
673  vertices->push_back(pTransform.TransformPoint(vertex1));
674  vertices->push_back(pTransform.TransformPoint(vertex2));
675  vertices->push_back(pTransform.TransformPoint(vertex3));
676  }
677  else
678  {
679  DumpInfo();
680  G4Exception("G4Tet::CreateRotatedVertices()",
681  "GeomSolids0003", FatalException,
682  "Error in allocation of vertices. Out of memory !");
683  }
684  return vertices;
685 }
686 
688 //
689 // GetEntityType
690 
692 {
693  return G4String("G4Tet");
694 }
695 
697 //
698 // Make a clone of the object
699 
701 {
702  return new G4Tet(*this);
703 }
704 
706 //
707 // Stream object contents to an output stream
708 
709 std::ostream& G4Tet::StreamInfo(std::ostream& os) const
710 {
711  G4int oldprc = os.precision(16);
712  os << "-----------------------------------------------------------\n"
713  << " *** Dump for solid - " << GetName() << " ***\n"
714  << " ===================================================\n"
715  << " Solid type: G4Tet\n"
716  << " Parameters: \n"
717  << " anchor: " << fAnchor/mm << " mm \n"
718  << " p2: " << fP2/mm << " mm \n"
719  << " p3: " << fP3/mm << " mm \n"
720  << " p4: " << fP4/mm << " mm \n"
721  << " normal123: " << fNormal123 << " \n"
722  << " normal134: " << fNormal134 << " \n"
723  << " normal142: " << fNormal142 << " \n"
724  << " normal234: " << fNormal234 << " \n"
725  << "-----------------------------------------------------------\n";
726  os.precision(oldprc);
727 
728  return os;
729 }
730 
731 
733 //
734 // GetPointOnFace
735 //
736 // Auxiliary method for get point on surface
737 
738 G4ThreeVector G4Tet::GetPointOnFace(G4ThreeVector p1, G4ThreeVector p2,
739  G4ThreeVector p3, G4double& area) const
740 {
741  G4double lambda1,lambda2;
742  G4ThreeVector v, w;
743 
744  v = p3 - p1;
745  w = p1 - p2;
746 
747  lambda1 = RandFlat::shoot(0.,1.);
748  lambda2 = RandFlat::shoot(0.,lambda1);
749 
750  area = 0.5*(v.cross(w)).mag();
751 
752  return (p2 + lambda1*w + lambda2*v);
753 }
754 
756 //
757 // GetPointOnSurface
758 
760 {
761  G4double chose,aOne,aTwo,aThree,aFour;
762  G4ThreeVector p1, p2, p3, p4;
763 
764  p1 = GetPointOnFace(fAnchor,fP2,fP3,aOne);
765  p2 = GetPointOnFace(fAnchor,fP4,fP3,aTwo);
766  p3 = GetPointOnFace(fAnchor,fP4,fP2,aThree);
767  p4 = GetPointOnFace(fP4,fP3,fP2,aFour);
768 
769  chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour);
770  if( (chose>=0.) && (chose <aOne) ) {return p1;}
771  else if( (chose>=aOne) && (chose < aOne+aTwo) ) {return p2;}
772  else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) ) {return p3;}
773  return p4;
774 }
775 
777 //
778 // GetVertices
779 
780 std::vector<G4ThreeVector> G4Tet::GetVertices() const
781 {
782  std::vector<G4ThreeVector> vertices(4);
783  vertices[0] = fAnchor;
784  vertices[1] = fP2;
785  vertices[2] = fP3;
786  vertices[3] = fP4;
787 
788  return vertices;
789 }
790 
792 //
793 // GetCubicVolume
794 
796 {
797  return fCubicVolume;
798 }
799 
801 //
802 // GetSurfaceArea
803 
805 {
806  return fSurfaceArea;
807 }
808 
809 // Methods for visualisation
810 
812 //
813 // DescribeYourselfTo
814 
816 {
817  scene.AddSolid (*this);
818 }
819 
821 //
822 // GetExtent
823 
825 {
826  return G4VisExtent (fXMin, fXMax, fYMin, fYMax, fZMin, fZMax);
827 }
828 
830 //
831 // CreatePolyhedron
832 
834 {
835  G4Polyhedron *ph=new G4Polyhedron;
836  G4double xyz[4][3];
837  const G4int faces[4][4]={{1,3,2,0},{1,4,3,0},{1,2,4,0},{2,3,4,0}};
838  xyz[0][0]=fAnchor.x(); xyz[0][1]=fAnchor.y(); xyz[0][2]=fAnchor.z();
839  xyz[1][0]=fP2.x(); xyz[1][1]=fP2.y(); xyz[1][2]=fP2.z();
840  xyz[2][0]=fP3.x(); xyz[2][1]=fP3.y(); xyz[2][2]=fP3.z();
841  xyz[3][0]=fP4.x(); xyz[3][1]=fP4.y(); xyz[3][2]=fP4.z();
842 
843  ph->createPolyhedron(4,4,xyz,faces);
844 
845  return ph;
846 }
847 
849 //
850 // GetPolyhedron
851 
853 {
854  if (!fpPolyhedron ||
856  fpPolyhedron->GetNumberOfRotationSteps())
857  {
858  delete fpPolyhedron;
859  fpPolyhedron = CreatePolyhedron();
860  }
861  return fpPolyhedron;
862 }
863 
864 #endif
G4String GetName() const
ThreeVector shoot(const G4int Ap, const G4int Af)
TTree * t1
Definition: plottest35.C:26
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Tet.cc:570
tuple a
Definition: test.py:11
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector GetPointOnSurface() const
Definition: G4Tet.cc:759
double x() const
double dot(const Hep3Vector &) const
G4bool IsYLimited() const
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Tet.cc:477
const char * p
Definition: xmltok.h:285
G4GeometryType GetEntityType() const
Definition: G4Tet.cc:691
G4bool IsRotated() const
G4double G4NeutronHPJENDLHEData::G4double result
G4ThreeVector NetTranslation() const
G4bool IsXLimited() const
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
double z() const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Tet.cc:291
void DumpInfo() const
std::vector< G4ThreeVector > GetVertices() const
Definition: G4Tet.cc:780
G4double GetMaxXExtent() const
Char_t n[5]
G4double GetMinZExtent() const
Definition: G4Tet.hh:65
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Tet.cc:815
bool G4bool
Definition: G4Types.hh:79
static G4bool CheckDegeneracy(G4ThreeVector anchor, G4ThreeVector p2, G4ThreeVector p3, G4ThreeVector p4)
Definition: G4Tet.cc:265
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
G4Polyhedron * CreatePolyhedron() const
Definition: G4Tet.cc:833
tuple v
Definition: test.py:18
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetMinXExtent() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double GetMaxZExtent() const
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Tet.cc:709
static G4int GetNumberOfRotationSteps()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
Hep3Vector unit() const
G4double GetCubicVolume()
Definition: G4Tet.cc:795
double y() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetSurfaceArea()
Definition: G4Tet.cc:804
EInside Inside(const G4ThreeVector &p) const
Definition: G4Tet.cc:386
G4VisExtent GetExtent() const
Definition: G4Tet.cc:824
TTree * t2
Definition: plottest35.C:36
#define G4endl
Definition: G4ios.hh:61
G4double GetMaxYExtent() const
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:110
G4double kCarTolerance
Definition: G4VSolid.hh:305
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
Hep3Vector cross(const Hep3Vector &) const
G4VSolid * Clone() const
Definition: G4Tet.cc:700
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Tet.cc:660
G4Tet & operator=(const G4Tet &rhs)
Definition: G4Tet.cc:234
double G4double
Definition: G4Types.hh:76
virtual ~G4Tet()
Definition: G4Tet.cc:204
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Tet.cc:281
double mag() const
G4bool IsZLimited() const
G4Polyhedron * GetPolyhedron() const
Definition: G4Tet.cc:852
G4Tet(const G4String &pName, G4ThreeVector anchor, G4ThreeVector p2, G4ThreeVector p3, G4ThreeVector p4, G4bool *degeneracyFlag=0)
Definition: G4Tet.cc:89
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Tet.cc:417