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