Geant4  10.02
G4Trd.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 result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4Trd.cc 83572 2014-09-01 15:23:27Z gcosmo $
28 //
29 //
30 // Implementation for G4Trd class
31 //
32 // History:
33 //
34 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal
35 // 26.04.05, V.Grichine, new SurfaceNoramal is default
36 // 07.12.04, V.Grichine, SurfaceNoramal with edges/vertices.
37 // 07.05.00, V.Grichine, in d = DistanceToIn(p,v), if d<0.5*kCarTolerance, d=0
38 // ~1996, V.Grichine, 1st implementation based on old code of P.Kent
39 //
41 
42 #include "G4Trd.hh"
43 
44 #if !defined(G4GEOM_USE_UTRD)
45 
46 #include "G4VPVParameterisation.hh"
47 #include "G4VoxelLimits.hh"
48 #include "G4AffineTransform.hh"
49 #include "Randomize.hh"
50 
51 #include "G4VGraphicsScene.hh"
52 
53 using namespace CLHEP;
54 
56 //
57 // Constructor - check & set half widths
58 
59 G4Trd::G4Trd( const G4String& pName,
60  G4double pdx1, G4double pdx2,
61  G4double pdy1, G4double pdy2,
62  G4double pdz )
63  : G4CSGSolid(pName)
64 {
65  CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz);
66 }
67 
69 //
70 // Set and check (coplanarity) of trd parameters
71 
73  G4double pdy1, G4double pdy2,
74  G4double pdz )
75 {
76  if ( pdx1>0&&pdx2>0&&pdy1>0&&pdy2>0&&pdz>0 )
77  {
78  fDx1=pdx1; fDx2=pdx2;
79  fDy1=pdy1; fDy2=pdy2;
80  fDz=pdz;
81  }
82  else
83  {
84  if ( pdx1>=0 && pdx2>=0 && pdy1>=0 && pdy2>=0 && pdz>=0 )
85  {
86  // G4double Minimum_length= (1+per_thousand) * kCarTolerance/2.;
87  // FIX-ME : temporary solution for ZERO or very-small parameters
88  //
89  G4double Minimum_length= kCarTolerance/2.;
90  fDx1=std::max(pdx1,Minimum_length);
91  fDx2=std::max(pdx2,Minimum_length);
92  fDy1=std::max(pdy1,Minimum_length);
93  fDy2=std::max(pdy2,Minimum_length);
94  fDz=std::max(pdz,Minimum_length);
95  }
96  else
97  {
98  std::ostringstream message;
99  message << "Invalid negative dimensions for Solid: " << GetName()
100  << G4endl
101  << " X - " << pdx1 << ", " << pdx2 << G4endl
102  << " Y - " << pdy1 << ", " << pdy2 << G4endl
103  << " Z - " << pdz;
104  G4Exception("G4Trd::CheckAndSetAllParameters()",
105  "GeomSolids0002", FatalException, message);
106  }
107  }
108  fCubicVolume= 0.;
109  fSurfaceArea= 0.;
110  fRebuildPolyhedron = true;
111 }
112 
114 //
115 // Fake default constructor - sets only member data and allocates memory
116 // for usage restricted to object persistency.
117 //
118 G4Trd::G4Trd( __void__& a )
119  : G4CSGSolid(a), fDx1(0.), fDx2(0.), fDy1(0.), fDy2(0.), fDz(0.)
120 {
121 }
122 
124 //
125 // Destructor
126 
128 {
129 }
130 
132 //
133 // Copy constructor
134 
135 G4Trd::G4Trd(const G4Trd& rhs)
136  : G4CSGSolid(rhs), fDx1(rhs.fDx1), fDx2(rhs.fDx2),
137  fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fDz)
138 {
139 }
140 
142 //
143 // Assignment operator
144 
146 {
147  // Check assignment to self
148  //
149  if (this == &rhs) { return *this; }
150 
151  // Copy base class data
152  //
154 
155  // Copy data
156  //
157  fDx1 = rhs.fDx1; fDx2 = rhs.fDx2;
158  fDy1 = rhs.fDy1; fDy2 = rhs.fDy2;
159  fDz = rhs.fDz;
160 
161  return *this;
162 }
163 
165 //
166 //
167 
169  G4double pdy2, G4double pdz )
170 {
171  CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz);
172 }
173 
174 
176 //
177 // Dispatch to parameterisation for replication mechanism dimension
178 // computation & modification.
179 
181  const G4int n,
182  const G4VPhysicalVolume* pRep )
183 {
184  p->ComputeDimensions(*this,n,pRep);
185 }
186 
187 
189 //
190 // Calculate extent under transform and specified limit
191 
193  const G4VoxelLimits& pVoxelLimit,
194  const G4AffineTransform& pTransform,
195  G4double& pMin, G4double& pMax ) const
196 {
197  if (!pTransform.IsRotated())
198  {
199  // Special case handling for unrotated solids
200  // Compute x/y/z mins and maxs respecting limits, with early returns
201  // if outside limits. Then switch() on pAxis
202 
203  G4double xoffset,xMin,xMax;
204  G4double yoffset,yMin,yMax;
205  G4double zoffset,zMin,zMax;
206 
207  zoffset=pTransform.NetTranslation().z();
208  zMin=zoffset-fDz;
209  zMax=zoffset+fDz;
210  if (pVoxelLimit.IsZLimited())
211  {
212  if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
213  || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
214  {
215  return false;
216  }
217  else
218  {
219  if (zMin<pVoxelLimit.GetMinZExtent())
220  {
221  zMin=pVoxelLimit.GetMinZExtent();
222  }
223  if (zMax>pVoxelLimit.GetMaxZExtent())
224  {
225  zMax=pVoxelLimit.GetMaxZExtent();
226  }
227  }
228  }
229  xoffset=pTransform.NetTranslation().x();
230  if (fDx2 >= fDx1)
231  {
232  xMax = xoffset+(fDx1+fDx2)/2+(zMax-zoffset)*(fDx2-fDx1)/(2*fDz) ;
233  xMin = 2*xoffset - xMax ;
234  }
235  else
236  {
237  xMax = xoffset+(fDx1+fDx2)/2+(zMin-zoffset)*(fDx2-fDx1)/(2*fDz) ;
238  xMin = 2*xoffset - xMax ;
239  }
240  if (pVoxelLimit.IsXLimited())
241  {
242  if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
243  || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
244  {
245  return false;
246  }
247  else
248  {
249  if (xMin<pVoxelLimit.GetMinXExtent())
250  {
251  xMin=pVoxelLimit.GetMinXExtent();
252  }
253  if (xMax>pVoxelLimit.GetMaxXExtent())
254  {
255  xMax=pVoxelLimit.GetMaxXExtent();
256  }
257  }
258  }
259  yoffset= pTransform.NetTranslation().y() ;
260  if(fDy2 >= fDy1)
261  {
262  yMax = yoffset+(fDy2+fDy1)/2+(zMax-zoffset)*(fDy2-fDy1)/(2*fDz) ;
263  yMin = 2*yoffset - yMax ;
264  }
265  else
266  {
267  yMax = yoffset+(fDy2+fDy1)/2+(zMin-zoffset)*(fDy2-fDy1)/(2*fDz) ;
268  yMin = 2*yoffset - yMax ;
269  }
270  if (pVoxelLimit.IsYLimited())
271  {
272  if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
273  || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
274  {
275  return false;
276  }
277  else
278  {
279  if (yMin<pVoxelLimit.GetMinYExtent())
280  {
281  yMin=pVoxelLimit.GetMinYExtent();
282  }
283  if (yMax>pVoxelLimit.GetMaxYExtent())
284  {
285  yMax=pVoxelLimit.GetMaxYExtent();
286  }
287  }
288  }
289 
290  switch (pAxis)
291  {
292  case kXAxis:
293  pMin=xMin;
294  pMax=xMax;
295  break;
296  case kYAxis:
297  pMin=yMin;
298  pMax=yMax;
299  break;
300  case kZAxis:
301  pMin=zMin;
302  pMax=zMax;
303  break;
304  default:
305  break;
306  }
307 
308  // Add 2*Tolerance to avoid precision troubles ?
309  //
310  pMin-=kCarTolerance;
311  pMax+=kCarTolerance;
312 
313  return true;
314  }
315  else
316  {
317  // General rotated case - create and clip mesh to boundaries
318 
319  G4bool existsAfterClip=false;
320  G4ThreeVectorList *vertices;
321 
322  pMin=+kInfinity;
323  pMax=-kInfinity;
324 
325  // Calculate rotated vertex coordinates
326  //
327  vertices=CreateRotatedVertices(pTransform);
328  ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
329  ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
330  ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
331 
332  if (pMin!=kInfinity||pMax!=-kInfinity)
333  {
334  existsAfterClip=true;
335 
336  // Add 2*tolerance to avoid precision troubles
337  //
338  pMin-=kCarTolerance;
339  pMax+=kCarTolerance;
340 
341  }
342  else
343  {
344  // Check for case where completely enveloping clipping volume
345  // If point inside then we are confident that the solid completely
346  // envelopes the clipping volume. Hence set min/max extents according
347  // to clipping volume extents along the specified axis.
348 
349  G4ThreeVector clipCentre(
350  (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
351  (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
352  (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
353 
354  if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
355  {
356  existsAfterClip=true;
357  pMin=pVoxelLimit.GetMinExtent(pAxis);
358  pMax=pVoxelLimit.GetMaxExtent(pAxis);
359  }
360  }
361  delete vertices;
362  return existsAfterClip;
363  }
364 }
365 
367 //
368 // Return whether point inside/outside/on surface, using tolerance
369 
371 {
372  EInside in=kOutside;
373  G4double x,y,zbase1,zbase2;
374 
375  if (std::fabs(p.z())<=fDz-kCarTolerance/2)
376  {
377  zbase1=p.z()+fDz; // Dist from -ve z plane
378  zbase2=fDz-p.z(); // Dist from +ve z plane
379 
380  // Check whether inside x tolerance
381  //
382  x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz - kCarTolerance/2;
383  if (std::fabs(p.x())<=x)
384  {
385  y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz - kCarTolerance/2;
386  if (std::fabs(p.y())<=y)
387  {
388  in=kInside;
389  }
390  else if (std::fabs(p.y())<=y+kCarTolerance)
391  {
392  in=kSurface;
393  }
394  }
395  else if (std::fabs(p.x())<=x+kCarTolerance)
396  {
397  // y = y half width of shape at z of point + tolerant boundary
398  //
399  y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz + kCarTolerance/2;
400  if (std::fabs(p.y())<=y)
401  {
402  in=kSurface;
403  }
404  }
405  }
406  else if (std::fabs(p.z())<=fDz+kCarTolerance/2)
407  {
408  // Only need to check outer tolerant boundaries
409  //
410  zbase1=p.z()+fDz; // Dist from -ve z plane
411  zbase2=fDz-p.z(); // Dist from +ve z plane
412 
413  // x = x half width of shape at z of point plus tolerance
414  //
415  x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz + kCarTolerance/2;
416  if (std::fabs(p.x())<=x)
417  {
418  // y = y half width of shape at z of point
419  //
420  y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz + kCarTolerance/2;
421  if (std::fabs(p.y())<=y) in=kSurface;
422  }
423  }
424  return in;
425 }
426 
428 //
429 // Calculate side nearest to p, and return normal
430 // If two sides are equidistant, normal of first side (x/y/z)
431 // encountered returned
432 
434 {
435  G4ThreeVector norm, sumnorm(0.,0.,0.);
436  G4int noSurfaces = 0;
437  G4double z = 2.0*fDz, tanx, secx, newpx, widx;
438  G4double tany, secy, newpy, widy;
439  G4double distx, disty, distz, fcos;
440  G4double delta = 0.5*kCarTolerance;
441 
442  tanx = (fDx2 - fDx1)/z;
443  secx = std::sqrt(1.0+tanx*tanx);
444  newpx = std::fabs(p.x())-p.z()*tanx;
445  widx = fDx2 - fDz*tanx;
446 
447  tany = (fDy2 - fDy1)/z;
448  secy = std::sqrt(1.0+tany*tany);
449  newpy = std::fabs(p.y())-p.z()*tany;
450  widy = fDy2 - fDz*tany;
451 
452  distx = std::fabs(newpx-widx)/secx; // perp. distance to x side
453  disty = std::fabs(newpy-widy)/secy; // to y side
454  distz = std::fabs(std::fabs(p.z())-fDz); // to z side
455 
456  fcos = 1.0/secx;
457  G4ThreeVector nX = G4ThreeVector( fcos,0,-tanx*fcos);
458  G4ThreeVector nmX = G4ThreeVector(-fcos,0,-tanx*fcos);
459 
460  fcos = 1.0/secy;
461  G4ThreeVector nY = G4ThreeVector(0, fcos,-tany*fcos);
462  G4ThreeVector nmY = G4ThreeVector(0,-fcos,-tany*fcos);
463  G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0);
464 
465  if (distx <= delta)
466  {
467  noSurfaces ++;
468  if ( p.x() >= 0.) sumnorm += nX;
469  else sumnorm += nmX;
470  }
471  if (disty <= delta)
472  {
473  noSurfaces ++;
474  if ( p.y() >= 0.) sumnorm += nY;
475  else sumnorm += nmY;
476  }
477  if (distz <= delta)
478  {
479  noSurfaces ++;
480  if ( p.z() >= 0.) sumnorm += nZ;
481  else sumnorm -= nZ;
482  }
483  if ( noSurfaces == 0 )
484  {
485 #ifdef G4CSGDEBUG
486  G4Exception("G4Trd::SurfaceNormal(p)", "GeomSolids1002", JustWarning,
487  "Point p is not on surface !?" );
488 #endif
489  norm = ApproxSurfaceNormal(p);
490  }
491  else if ( noSurfaces == 1 ) norm = sumnorm;
492  else norm = sumnorm.unit();
493  return norm;
494 }
495 
496 
498 //
499 // Algorithm for SurfaceNormal() following the original specification
500 // for points not on the surface
501 
503 {
504  G4ThreeVector norm;
505  G4double z,tanx,secx,newpx,widx;
506  G4double tany,secy,newpy,widy;
507  G4double distx,disty,distz,fcos;
508 
509  z=2.0*fDz;
510 
511  tanx=(fDx2-fDx1)/z;
512  secx=std::sqrt(1.0+tanx*tanx);
513  newpx=std::fabs(p.x())-p.z()*tanx;
514  widx=fDx2-fDz*tanx;
515 
516  tany=(fDy2-fDy1)/z;
517  secy=std::sqrt(1.0+tany*tany);
518  newpy=std::fabs(p.y())-p.z()*tany;
519  widy=fDy2-fDz*tany;
520 
521  distx=std::fabs(newpx-widx)/secx; // perpendicular distance to x side
522  disty=std::fabs(newpy-widy)/secy; // to y side
523  distz=std::fabs(std::fabs(p.z())-fDz); // to z side
524 
525  // find closest side
526  //
527  if (distx<=disty)
528  {
529  if (distx<=distz)
530  {
531  // Closest to X
532  //
533  fcos=1.0/secx;
534  // normal=(+/-std::cos(ang),0,-std::sin(ang))
535  if (p.x()>=0)
536  norm=G4ThreeVector(fcos,0,-tanx*fcos);
537  else
538  norm=G4ThreeVector(-fcos,0,-tanx*fcos);
539  }
540  else
541  {
542  // Closest to Z
543  //
544  if (p.z()>=0)
545  norm=G4ThreeVector(0,0,1);
546  else
547  norm=G4ThreeVector(0,0,-1);
548  }
549  }
550  else
551  {
552  if (disty<=distz)
553  {
554  // Closest to Y
555  //
556  fcos=1.0/secy;
557  if (p.y()>=0)
558  norm=G4ThreeVector(0,fcos,-tany*fcos);
559  else
560  norm=G4ThreeVector(0,-fcos,-tany*fcos);
561  }
562  else
563  {
564  // Closest to Z
565  //
566  if (p.z()>=0)
567  norm=G4ThreeVector(0,0,1);
568  else
569  norm=G4ThreeVector(0,0,-1);
570  }
571  }
572  return norm;
573 }
574 
576 //
577 // Calculate distance to shape from outside
578 // - return kInfinity if no intersection
579 //
580 // ALGORITHM:
581 // For each component, calculate pair of minimum and maximum intersection
582 // values for which the particle is in the extent of the shape
583 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
584 // - Z plane intersectin uses tolerance
585 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
586 // (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
587 // cases)
588 // - Note: XZ and YZ planes each divide space into four regions,
589 // characterised by ss1 ss2
590 // NOTE:
591 //
592 // `Inside' safe - meaningful answers given if point is inside the exact
593 // shape.
594 
596  const G4ThreeVector& v ) const
597 {
598  G4double snxt = kInfinity ; // snxt = default return value
599  G4double smin,smax;
600  G4double s1,s2,tanxz,tanyz,ds1,ds2;
601  G4double ss1,ss2,sn1=0.,sn2=0.,Dist;
602 
603  if ( v.z() ) // Calculate valid z intersect range
604  {
605  if ( v.z() > 0 ) // Calculate smax: must be +ve or no intersection.
606  {
607  Dist = fDz - p.z() ; // to plane at +dz
608 
609  if (Dist >= 0.5*kCarTolerance)
610  {
611  smax = Dist/v.z() ;
612  smin = -(fDz + p.z())/v.z() ;
613  }
614  else return snxt ;
615  }
616  else // v.z <0
617  {
618  Dist=fDz+p.z(); // plane at -dz
619 
620  if ( Dist >= 0.5*kCarTolerance )
621  {
622  smax = -Dist/v.z() ;
623  smin = (fDz - p.z())/v.z() ;
624  }
625  else return snxt ;
626  }
627  if (smin < 0 ) smin = 0 ;
628  }
629  else // v.z=0
630  {
631  if (std::fabs(p.z()) >= fDz ) return snxt ; // Outside & no intersect
632  else
633  {
634  smin = 0 ; // Always inside z range
635  smax = kInfinity;
636  }
637  }
638 
639  // Calculate x intersection range
640  //
641  // Calc half width at p.z, and components towards planes
642 
643  tanxz = (fDx2 - fDx1)*0.5/fDz ;
644  s1 = 0.5*(fDx1+fDx2) + tanxz*p.z() ; // x half width at p.z
645  ds1 = v.x() - tanxz*v.z() ; // Components of v towards faces at +-x
646  ds2 = v.x() + tanxz*v.z() ;
647  ss1 = s1 - p.x() ; // -delta x to +ve plane
648  // -ve when outside
649  ss2 = -s1 - p.x() ; // -delta x to -ve plane
650  // +ve when outside
651 
652  if (ss1 < 0 && ss2 <= 0 )
653  {
654  if (ds1 < 0) // In +ve coord Area
655  {
656  sn1 = ss1/ds1 ;
657 
658  if ( ds2 < 0 ) sn2 = ss2/ds2 ;
659  else sn2 = kInfinity ;
660  }
661  else return snxt ;
662  }
663  else if ( ss1 >= 0 && ss2 > 0 )
664  {
665  if ( ds2 > 0 ) // In -ve coord Area
666  {
667  sn1 = ss2/ds2 ;
668 
669  if (ds1 > 0) sn2 = ss1/ds1 ;
670  else sn2 = kInfinity;
671 
672  }
673  else return snxt ;
674  }
675  else if (ss1 >= 0 && ss2 <= 0 )
676  {
677  // Inside Area - calculate leaving distance
678  // *Don't* use exact distance to side for tolerance
679  // = ss1*std::cos(ang xz)
680  // = ss1/std::sqrt(1.0+tanxz*tanxz)
681  sn1 = 0 ;
682 
683  if ( ds1 > 0 )
684  {
685  if (ss1 > 0.5*kCarTolerance) sn2 = ss1/ds1 ; // Leave +ve side extent
686  else return snxt ; // Leave immediately by +ve
687  }
688  else sn2 = kInfinity ;
689 
690  if ( ds2 < 0 )
691  {
692  if ( ss2 < -0.5*kCarTolerance )
693  {
694  Dist = ss2/ds2 ; // Leave -ve side extent
695  if ( Dist < sn2 ) sn2 = Dist ;
696  }
697  else return snxt ;
698  }
699  }
700  else if (ss1 < 0 && ss2 > 0 )
701  {
702  // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0)
703 
704  if ( ds1 >= 0 || ds2 <= 0 )
705  {
706  return snxt ;
707  }
708  else // Will intersect & stay inside
709  {
710  sn1 = ss1/ds1 ;
711  Dist = ss2/ds2 ;
712  if (Dist > sn1 ) sn1 = Dist ;
713  sn2 = kInfinity ;
714  }
715  }
716 
717  // Reduce allowed range of distances as appropriate
718 
719  if ( sn1 > smin ) smin = sn1 ;
720  if ( sn2 < smax ) smax = sn2 ;
721 
722  // Check for incompatible ranges (eg z intersects between 50 ->100 and x
723  // only 10-40 -> no intersection)
724 
725  if ( smax < smin ) return snxt ;
726 
727  // Calculate valid y intersection range
728  // (repeat of x intersection code)
729 
730  tanyz = (fDy2-fDy1)*0.5/fDz ;
731  s2 = 0.5*(fDy1+fDy2) + tanyz*p.z() ; // y half width at p.z
732  ds1 = v.y() - tanyz*v.z() ; // Components of v towards faces at +-y
733  ds2 = v.y() + tanyz*v.z() ;
734  ss1 = s2 - p.y() ; // -delta y to +ve plane
735  ss2 = -s2 - p.y() ; // -delta y to -ve plane
736 
737  if ( ss1 < 0 && ss2 <= 0 )
738  {
739  if (ds1 < 0 ) // In +ve coord Area
740  {
741  sn1 = ss1/ds1 ;
742  if ( ds2 < 0 ) sn2 = ss2/ds2 ;
743  else sn2 = kInfinity ;
744  }
745  else return snxt ;
746  }
747  else if ( ss1 >= 0 && ss2 > 0 )
748  {
749  if ( ds2 > 0 ) // In -ve coord Area
750  {
751  sn1 = ss2/ds2 ;
752  if ( ds1 > 0 ) sn2 = ss1/ds1 ;
753  else sn2 = kInfinity ;
754  }
755  else return snxt ;
756  }
757  else if (ss1 >= 0 && ss2 <= 0 )
758  {
759  // Inside Area - calculate leaving distance
760  // *Don't* use exact distance to side for tolerance
761  // = ss1*std::cos(ang yz)
762  // = ss1/std::sqrt(1.0+tanyz*tanyz)
763  sn1 = 0 ;
764 
765  if ( ds1 > 0 )
766  {
767  if (ss1 > 0.5*kCarTolerance) sn2 = ss1/ds1 ; // Leave +ve side extent
768  else return snxt ; // Leave immediately by +ve
769  }
770  else sn2 = kInfinity ;
771 
772  if ( ds2 < 0 )
773  {
774  if ( ss2 < -0.5*kCarTolerance )
775  {
776  Dist = ss2/ds2 ; // Leave -ve side extent
777  if (Dist < sn2) sn2=Dist;
778  }
779  else return snxt ;
780  }
781  }
782  else if (ss1 < 0 && ss2 > 0 )
783  {
784  // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0)
785 
786  if (ds1 >= 0 || ds2 <= 0 )
787  {
788  return snxt ;
789  }
790  else // Will intersect & stay inside
791  {
792  sn1 = ss1/ds1 ;
793  Dist = ss2/ds2 ;
794  if (Dist > sn1 ) sn1 = Dist ;
795  sn2 = kInfinity ;
796  }
797  }
798 
799  // Reduce allowed range of distances as appropriate
800 
801  if ( sn1 > smin) smin = sn1 ;
802  if ( sn2 < smax) smax = sn2 ;
803 
804  // Check for incompatible ranges (eg x intersects between 50 ->100 and y
805  // only 10-40 -> no intersection). Set snxt if ok
806 
807  if ( smax > smin ) snxt = smin ;
808  if (snxt < 0.5*kCarTolerance ) snxt = 0.0 ;
809 
810  return snxt ;
811 }
812 
814 //
815 // Approximate distance to shape
816 // Calculate perpendicular distances to z/x/y surfaces, return largest
817 // which is the most fast estimation of shortest distance to Trd
818 // - Safe underestimate
819 // - If point within exact shape, return 0
820 
822 {
823  G4double safe=0.0;
824  G4double tanxz,distx,safx;
825  G4double tanyz,disty,safy;
826  G4double zbase;
827 
828  safe=std::fabs(p.z())-fDz;
829  if (safe<0) safe=0; // Also used to ensure x/y distances
830  // POSITIVE
831 
832  zbase=fDz+p.z();
833 
834  // Find distance along x direction to closest x plane
835  //
836  tanxz=(fDx2-fDx1)*0.5/fDz;
837  // widx=fDx1+tanxz*(fDz+p.z()); // x width at p.z
838  // distx=std::fabs(p.x())-widx; // distance to plane
839  distx=std::fabs(p.x())-(fDx1+tanxz*zbase);
840  if (distx>safe)
841  {
842  safx=distx/std::sqrt(1.0+tanxz*tanxz); // vector Dist=Dist*std::cos(ang)
843  if (safx>safe) safe=safx;
844  }
845 
846  // Find distance along y direction to slanted wall
847  tanyz=(fDy2-fDy1)*0.5/fDz;
848  // widy=fDy1+tanyz*(fDz+p.z()); // y width at p.z
849  // disty=std::fabs(p.y())-widy; // distance to plane
850  disty=std::fabs(p.y())-(fDy1+tanyz*zbase);
851  if (disty>safe)
852  {
853  safy=disty/std::sqrt(1.0+tanyz*tanyz); // distance along vector
854  if (safy>safe) safe=safy;
855  }
856  return safe;
857 }
858 
860 //
861 // Calcluate distance to surface of shape from inside
862 // Calculate distance to x/y/z planes - smallest is exiting distance
863 // - z planes have std. check for tolerance
864 // - xz yz planes have check based on distance || to x or y axis
865 // (not corrected for slope of planes)
866 // ?BUG? If v.z==0 are there cases when snside not set????
867 
869  const G4ThreeVector& v,
870  const G4bool calcNorm,
871  G4bool *validNorm,
872  G4ThreeVector *n ) const
873 {
874  ESide side = kUndefined, snside = kUndefined;
875  G4double snxt,pdist;
876  G4double central,ss1,ss2,ds1,ds2,sn=0.,sn2=0.;
877  G4double tanxz=0.,cosxz=0.,tanyz=0.,cosyz=0.;
878 
879  if (calcNorm) *validNorm=true; // All normals are valid
880 
881  // Calculate z plane intersection
882  if (v.z()>0)
883  {
884  pdist=fDz-p.z();
885  if (pdist>kCarTolerance/2)
886  {
887  snxt=pdist/v.z();
888  side=kPZ;
889  }
890  else
891  {
892  if (calcNorm)
893  {
894  *n=G4ThreeVector(0,0,1);
895  }
896  return snxt=0;
897  }
898  }
899  else if (v.z()<0)
900  {
901  pdist=fDz+p.z();
902  if (pdist>kCarTolerance/2)
903  {
904  snxt=-pdist/v.z();
905  side=kMZ;
906  }
907  else
908  {
909  if (calcNorm)
910  {
911  *n=G4ThreeVector(0,0,-1);
912  }
913  return snxt=0;
914  }
915  }
916  else
917  {
918  snxt=kInfinity;
919  }
920 
921  //
922  // Calculate x intersection
923  //
924  tanxz=(fDx2-fDx1)*0.5/fDz;
925  central=0.5*(fDx1+fDx2);
926 
927  // +ve plane (1)
928  //
929  ss1=central+tanxz*p.z()-p.x(); // distance || x axis to plane
930  // (+ve if point inside)
931  ds1=v.x()-tanxz*v.z(); // component towards plane at +x
932  // (-ve if +ve -> -ve direction)
933  // -ve plane (2)
934  //
935  ss2=-tanxz*p.z()-p.x()-central; //distance || x axis to plane
936  // (-ve if point inside)
937  ds2=tanxz*v.z()+v.x(); // component towards plane at -x
938 
939  if (ss1>0&&ss2<0)
940  {
941  // Normal case - entirely inside region
942  if (ds1<=0&&ds2<0)
943  {
944  if (ss2<-kCarTolerance/2)
945  {
946  sn=ss2/ds2; // Leave by -ve side
947  snside=kMX;
948  }
949  else
950  {
951  sn=0; // Leave immediately by -ve side
952  snside=kMX;
953  }
954  }
955  else if (ds1>0&&ds2>=0)
956  {
957  if (ss1>kCarTolerance/2)
958  {
959  sn=ss1/ds1; // Leave by +ve side
960  snside=kPX;
961  }
962  else
963  {
964  sn=0; // Leave immediately by +ve side
965  snside=kPX;
966  }
967  }
968  else if (ds1>0&&ds2<0)
969  {
970  if (ss1>kCarTolerance/2)
971  {
972  // sn=ss1/ds1; // Leave by +ve side
973  if (ss2<-kCarTolerance/2)
974  {
975  sn=ss1/ds1; // Leave by +ve side
976  sn2=ss2/ds2;
977  if (sn2<sn)
978  {
979  sn=sn2;
980  snside=kMX;
981  }
982  else
983  {
984  snside=kPX;
985  }
986  }
987  else
988  {
989  sn=0; // Leave immediately by -ve
990  snside=kMX;
991  }
992  }
993  else
994  {
995  sn=0; // Leave immediately by +ve side
996  snside=kPX;
997  }
998  }
999  else
1000  {
1001  // Must be || to both
1002  //
1003  sn=kInfinity; // Don't leave by either side
1004  }
1005  }
1006  else if (ss1<=0&&ss2<0)
1007  {
1008  // Outside, in +ve Area
1009 
1010  if (ds1>0)
1011  {
1012  sn=0; // Away from shape
1013  // Left by +ve side
1014  snside=kPX;
1015  }
1016  else
1017  {
1018  if (ds2<0)
1019  {
1020  // Ignore +ve plane and use -ve plane intersect
1021  //
1022  sn=ss2/ds2; // Leave by -ve side
1023  snside=kMX;
1024  }
1025  else
1026  {
1027  // Must be || to both -> exit determined by other axes
1028  //
1029  sn=kInfinity; // Don't leave by either side
1030  }
1031  }
1032  }
1033  else if (ss1>0&&ss2>=0)
1034  {
1035  // Outside, in -ve Area
1036 
1037  if (ds2<0)
1038  {
1039  sn=0; // away from shape
1040  // Left by -ve side
1041  snside=kMX;
1042  }
1043  else
1044  {
1045  if (ds1>0)
1046  {
1047  // Ignore +ve plane and use -ve plane intersect
1048  //
1049  sn=ss1/ds1; // Leave by +ve side
1050  snside=kPX;
1051  }
1052  else
1053  {
1054  // Must be || to both -> exit determined by other axes
1055  //
1056  sn=kInfinity; // Don't leave by either side
1057  }
1058  }
1059  }
1060 
1061  // Update minimum exit distance
1062 
1063  if (sn<snxt)
1064  {
1065  snxt=sn;
1066  side=snside;
1067  }
1068  if (snxt>0)
1069  {
1070  // Calculate y intersection
1071 
1072  tanyz=(fDy2-fDy1)*0.5/fDz;
1073  central=0.5*(fDy1+fDy2);
1074 
1075  // +ve plane (1)
1076  //
1077  ss1=central+tanyz*p.z()-p.y(); // distance || y axis to plane
1078  // (+ve if point inside)
1079  ds1=v.y()-tanyz*v.z(); // component towards +ve plane
1080  // (-ve if +ve -> -ve direction)
1081  // -ve plane (2)
1082  //
1083  ss2=-tanyz*p.z()-p.y()-central; // distance || y axis to plane
1084  // (-ve if point inside)
1085  ds2=tanyz*v.z()+v.y(); // component towards -ve plane
1086 
1087  if (ss1>0&&ss2<0)
1088  {
1089  // Normal case - entirely inside region
1090 
1091  if (ds1<=0&&ds2<0)
1092  {
1093  if (ss2<-kCarTolerance/2)
1094  {
1095  sn=ss2/ds2; // Leave by -ve side
1096  snside=kMY;
1097  }
1098  else
1099  {
1100  sn=0; // Leave immediately by -ve side
1101  snside=kMY;
1102  }
1103  }
1104  else if (ds1>0&&ds2>=0)
1105  {
1106  if (ss1>kCarTolerance/2)
1107  {
1108  sn=ss1/ds1; // Leave by +ve side
1109  snside=kPY;
1110  }
1111  else
1112  {
1113  sn=0; // Leave immediately by +ve side
1114  snside=kPY;
1115  }
1116  }
1117  else if (ds1>0&&ds2<0)
1118  {
1119  if (ss1>kCarTolerance/2)
1120  {
1121  // sn=ss1/ds1; // Leave by +ve side
1122  if (ss2<-kCarTolerance/2)
1123  {
1124  sn=ss1/ds1; // Leave by +ve side
1125  sn2=ss2/ds2;
1126  if (sn2<sn)
1127  {
1128  sn=sn2;
1129  snside=kMY;
1130  }
1131  else
1132  {
1133  snside=kPY;
1134  }
1135  }
1136  else
1137  {
1138  sn=0; // Leave immediately by -ve
1139  snside=kMY;
1140  }
1141  }
1142  else
1143  {
1144  sn=0; // Leave immediately by +ve side
1145  snside=kPY;
1146  }
1147  }
1148  else
1149  {
1150  // Must be || to both
1151  //
1152  sn=kInfinity; // Don't leave by either side
1153  }
1154  }
1155  else if (ss1<=0&&ss2<0)
1156  {
1157  // Outside, in +ve Area
1158 
1159  if (ds1>0)
1160  {
1161  sn=0; // Away from shape
1162  // Left by +ve side
1163  snside=kPY;
1164  }
1165  else
1166  {
1167  if (ds2<0)
1168  {
1169  // Ignore +ve plane and use -ve plane intersect
1170  //
1171  sn=ss2/ds2; // Leave by -ve side
1172  snside=kMY;
1173  }
1174  else
1175  {
1176  // Must be || to both -> exit determined by other axes
1177  //
1178  sn=kInfinity; // Don't leave by either side
1179  }
1180  }
1181  }
1182  else if (ss1>0&&ss2>=0)
1183  {
1184  // Outside, in -ve Area
1185  if (ds2<0)
1186  {
1187  sn=0; // away from shape
1188  // Left by -ve side
1189  snside=kMY;
1190  }
1191  else
1192  {
1193  if (ds1>0)
1194  {
1195  // Ignore +ve plane and use -ve plane intersect
1196  //
1197  sn=ss1/ds1; // Leave by +ve side
1198  snside=kPY;
1199  }
1200  else
1201  {
1202  // Must be || to both -> exit determined by other axes
1203  //
1204  sn=kInfinity; // Don't leave by either side
1205  }
1206  }
1207  }
1208 
1209  // Update minimum exit distance
1210 
1211  if (sn<snxt)
1212  {
1213  snxt=sn;
1214  side=snside;
1215  }
1216  }
1217 
1218  if (calcNorm)
1219  {
1220  switch (side)
1221  {
1222  case kPX:
1223  cosxz=1.0/std::sqrt(1.0+tanxz*tanxz);
1224  *n=G4ThreeVector(cosxz,0,-tanxz*cosxz);
1225  break;
1226  case kMX:
1227  cosxz=-1.0/std::sqrt(1.0+tanxz*tanxz);
1228  *n=G4ThreeVector(cosxz,0,tanxz*cosxz);
1229  break;
1230  case kPY:
1231  cosyz=1.0/std::sqrt(1.0+tanyz*tanyz);
1232  *n=G4ThreeVector(0,cosyz,-tanyz*cosyz);
1233  break;
1234  case kMY:
1235  cosyz=-1.0/std::sqrt(1.0+tanyz*tanyz);
1236  *n=G4ThreeVector(0,cosyz,tanyz*cosyz);
1237  break;
1238  case kPZ:
1239  *n=G4ThreeVector(0,0,1);
1240  break;
1241  case kMZ:
1242  *n=G4ThreeVector(0,0,-1);
1243  break;
1244  default:
1245  DumpInfo();
1246  G4Exception("G4Trd::DistanceToOut(p,v,..)",
1247  "GeomSolids1002", JustWarning,
1248  "Undefined side for valid surface normal to solid.");
1249  break;
1250  }
1251  }
1252  return snxt;
1253 }
1254 
1256 //
1257 // Calculate exact shortest distance to any boundary from inside
1258 // - Returns 0 is point outside
1259 
1261 {
1262  G4double safe=0.0;
1263  G4double tanxz,xdist,saf1;
1264  G4double tanyz,ydist,saf2;
1265  G4double zbase;
1266 
1267 #ifdef G4CSGDEBUG
1268  if( Inside(p) == kOutside )
1269  {
1270  G4int oldprc = G4cout.precision(16) ;
1271  G4cout << G4endl ;
1272  DumpInfo();
1273  G4cout << "Position:" << G4endl << G4endl ;
1274  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1275  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1276  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1277  G4cout.precision(oldprc) ;
1278  G4Exception("G4Trd::DistanceToOut(p)", "GeomSolids1002", JustWarning,
1279  "Point p is outside !?" );
1280  }
1281 #endif
1282 
1283  safe=fDz-std::fabs(p.z()); // z perpendicular Dist
1284 
1285  zbase=fDz+p.z();
1286 
1287  // xdist = distance perpendicular to z axis to closest x plane from p
1288  // = (x half width of shape at p.z) - std::fabs(p.x)
1289  //
1290  tanxz=(fDx2-fDx1)*0.5/fDz;
1291  xdist=fDx1+tanxz*zbase-std::fabs(p.x());
1292  saf1=xdist/std::sqrt(1.0+tanxz*tanxz); // x*std::cos(ang_xz) =
1293  // shortest (perpendicular)
1294  // distance to plane
1295  tanyz=(fDy2-fDy1)*0.5/fDz;
1296  ydist=fDy1+tanyz*zbase-std::fabs(p.y());
1297  saf2=ydist/std::sqrt(1.0+tanyz*tanyz);
1298 
1299  // Return minimum x/y/z distance
1300  //
1301  if (safe>saf1) safe=saf1;
1302  if (safe>saf2) safe=saf2;
1303 
1304  if (safe<0) safe=0;
1305  return safe;
1306 }
1307 
1309 //
1310 // Create a List containing the transformed vertices
1311 // Ordering [0-3] -fDz cross section
1312 // [4-7] +fDz cross section such that [0] is below [4],
1313 // [1] below [5] etc.
1314 // Note:
1315 // Caller has deletion resposibility
1316 
1319 {
1320  G4ThreeVectorList *vertices;
1321  vertices=new G4ThreeVectorList();
1322  if (vertices)
1323  {
1324  vertices->reserve(8);
1325  G4ThreeVector vertex0(-fDx1,-fDy1,-fDz);
1326  G4ThreeVector vertex1(fDx1,-fDy1,-fDz);
1327  G4ThreeVector vertex2(fDx1,fDy1,-fDz);
1328  G4ThreeVector vertex3(-fDx1,fDy1,-fDz);
1329  G4ThreeVector vertex4(-fDx2,-fDy2,fDz);
1330  G4ThreeVector vertex5(fDx2,-fDy2,fDz);
1331  G4ThreeVector vertex6(fDx2,fDy2,fDz);
1332  G4ThreeVector vertex7(-fDx2,fDy2,fDz);
1333 
1334  vertices->push_back(pTransform.TransformPoint(vertex0));
1335  vertices->push_back(pTransform.TransformPoint(vertex1));
1336  vertices->push_back(pTransform.TransformPoint(vertex2));
1337  vertices->push_back(pTransform.TransformPoint(vertex3));
1338  vertices->push_back(pTransform.TransformPoint(vertex4));
1339  vertices->push_back(pTransform.TransformPoint(vertex5));
1340  vertices->push_back(pTransform.TransformPoint(vertex6));
1341  vertices->push_back(pTransform.TransformPoint(vertex7));
1342  }
1343  else
1344  {
1345  DumpInfo();
1346  G4Exception("G4Trd::CreateRotatedVertices()",
1347  "GeomSolids0003", FatalException,
1348  "Error in allocation of vertices. Out of memory !");
1349  }
1350  return vertices;
1351 }
1352 
1354 //
1355 // GetEntityType
1356 
1358 {
1359  return G4String("G4Trd");
1360 }
1361 
1363 //
1364 // Make a clone of the object
1365 //
1367 {
1368  return new G4Trd(*this);
1369 }
1370 
1372 //
1373 // Stream object contents to an output stream
1374 
1375 std::ostream& G4Trd::StreamInfo( std::ostream& os ) const
1376 {
1377  G4int oldprc = os.precision(16);
1378  os << "-----------------------------------------------------------\n"
1379  << " *** Dump for solid - " << GetName() << " ***\n"
1380  << " ===================================================\n"
1381  << " Solid type: G4Trd\n"
1382  << " Parameters: \n"
1383  << " half length X, surface -dZ: " << fDx1/mm << " mm \n"
1384  << " half length X, surface +dZ: " << fDx2/mm << " mm \n"
1385  << " half length Y, surface -dZ: " << fDy1/mm << " mm \n"
1386  << " half length Y, surface +dZ: " << fDy2/mm << " mm \n"
1387  << " half length Z : " << fDz/mm << " mm \n"
1388  << "-----------------------------------------------------------\n";
1389  os.precision(oldprc);
1390 
1391  return os;
1392 }
1393 
1394 
1396 //
1397 // GetPointOnSurface
1398 //
1399 // Return a point (G4ThreeVector) randomly and uniformly
1400 // selected on the solid surface
1401 
1403 {
1404  G4double px, py, pz, tgX, tgY, secX, secY, select, sumS, tmp;
1405  G4double Sxy1, Sxy2, Sxy, Sxz, Syz;
1406 
1407  tgX = 0.5*(fDx2-fDx1)/fDz;
1408  secX = std::sqrt(1+tgX*tgX);
1409  tgY = 0.5*(fDy2-fDy1)/fDz;
1410  secY = std::sqrt(1+tgY*tgY);
1411 
1412  // calculate 0.25 of side surfaces, sumS is 0.25 of total surface
1413 
1414  Sxy1 = fDx1*fDy1;
1415  Sxy2 = fDx2*fDy2;
1416  Sxy = Sxy1 + Sxy2;
1417  Sxz = (fDx1 + fDx2)*fDz*secY;
1418  Syz = (fDy1 + fDy2)*fDz*secX;
1419  sumS = Sxy + Sxz + Syz;
1420 
1421  select = sumS*G4UniformRand();
1422 
1423  if( select < Sxy ) // Sxy1 or Sxy2
1424  {
1425  if( select < Sxy1 )
1426  {
1427  pz = -fDz;
1428  px = -fDx1 + 2*fDx1*G4UniformRand();
1429  py = -fDy1 + 2*fDy1*G4UniformRand();
1430  }
1431  else
1432  {
1433  pz = fDz;
1434  px = -fDx2 + 2*fDx2*G4UniformRand();
1435  py = -fDy2 + 2*fDy2*G4UniformRand();
1436  }
1437  }
1438  else if ( ( select - Sxy ) < Sxz ) // Sxz
1439  {
1440  pz = -fDz + 2*fDz*G4UniformRand();
1441  tmp = fDx1 + (pz + fDz)*tgX;
1442  px = -tmp + 2*tmp*G4UniformRand();
1443  tmp = fDy1 + (pz + fDz)*tgY;
1444 
1445  if(G4UniformRand() > 0.5) { py = tmp; }
1446  else { py = -tmp; }
1447  }
1448  else // Syz
1449  {
1450  pz = -fDz + 2*fDz*G4UniformRand();
1451  tmp = fDy1 + (pz + fDz)*tgY;
1452  py = -tmp + 2*tmp*G4UniformRand();
1453  tmp = fDx1 + (pz + fDz)*tgX;
1454 
1455  if(G4UniformRand() > 0.5) { px = tmp; }
1456  else { px = -tmp; }
1457  }
1458  return G4ThreeVector(px,py,pz);
1459 }
1460 
1462 //
1463 // Methods for visualisation
1464 
1466 {
1467  scene.AddSolid (*this);
1468 }
1469 
1471 {
1472  return new G4PolyhedronTrd2 (fDx1, fDx2, fDy1, fDy2, fDz);
1473 }
1474 
1475 #endif
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Trd.cc:1465
void SetAllParameters(G4double pdx1, G4double pdx2, G4double pdy1, G4double pdy2, G4double pdz)
Definition: G4Trd.cc:168
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Trd.cc:1318
G4double fDy2
Definition: G4Trd.hh:181
G4String GetName() const
~G4Trd()
Definition: G4Trd.cc:127
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:347
ESide
Definition: G4Trd.hh:157
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4bool fRebuildPolyhedron
Definition: G4CSGSolid.hh:80
G4AffineTransform Inverse() const
G4double z
Definition: TRTMaterials.hh:39
G4bool IsYLimited() const
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Trd.cc:1375
G4bool IsRotated() const
G4double fCubicVolume
Definition: G4CSGSolid.hh:78
G4ThreeVector NetTranslation() const
G4GeometryType GetEntityType() const
Definition: G4Trd.cc:1357
G4bool IsXLimited() const
G4double a
Definition: TRTMaterials.hh:39
Definition: G4Trd.hh:72
G4double fcos(G4double arg)
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
const G4int smax
void DumpInfo() const
G4VSolid * Clone() const
Definition: G4Trd.cc:1366
G4ThreeVector GetPointOnSurface() const
Definition: G4Trd.cc:1402
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Trd.cc:502
G4double fDy1
Definition: G4Trd.hh:181
G4double GetMaxXExtent() const
G4double GetMinZExtent() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79
G4GLOB_DLL std::ostream G4cout
G4Polyhedron * CreatePolyhedron() const
Definition: G4Trd.cc:1470
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Trd.cc:595
bool G4bool
Definition: G4Types.hh:79
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
EInside Inside(const G4ThreeVector &p) const
Definition: G4Trd.cc:370
const G4int n
G4double fDx2
Definition: G4Trd.hh:181
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4double GetMinXExtent() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double GetMaxZExtent() const
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
const G4double x[NPOINTSGL]
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Trd.cc:868
void CheckAndSetAllParameters(G4double pdx1, G4double pdx2, G4double pdy1, G4double pdy2, G4double pdz)
Definition: G4Trd.cc:72
G4Trd & operator=(const G4Trd &rhs)
Definition: G4Trd.cc:145
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Trd.cc:180
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Trd.cc:433
#define G4endl
Definition: G4ios.hh:61
G4double GetMaxYExtent() const
G4double kCarTolerance
Definition: G4VSolid.hh:304
G4Trd(const G4String &pName, G4double pdx1, G4double pdx2, G4double pdy1, G4double pdy2, G4double pdz)
Definition: G4Trd.cc:59
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
Definition: G4Trd.cc:192
double G4double
Definition: G4Types.hh:76
G4double GetMaxExtent(const EAxis pAxis) const
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:91
G4double fDz
Definition: G4Trd.hh:181
G4double fDx1
Definition: G4Trd.hh:181
G4bool IsZLimited() const
static const double mm
Definition: G4SIunits.hh:114
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:378
G4double GetMinExtent(const EAxis pAxis) const