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