Geant4  10.00.p02
UTrd.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UTrd
12 //
13 // 19.10.12 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include <iostream>
18 #include <cmath>
19 
20 #include "UTrd.hh"
21 #include "UUtils.hh"
22 
23 //dx1 Half-length along x at the surface positioned at -dz
24 //dx2 Half-length along x at the surface positioned at +dz
25 //dy1 Half-length along y at the surface positioned at -dz
26 //dy2 Half-length along y at the surface positioned at +dz
27 //dz Half-length along z axis
28 //______________________________________________________________________________
29 UTrd::UTrd(const std::string& pName,
30  double pdx1, double pdx2,
31  double pdy1, double pdy2,
32  double pdz)
33  : VUSolid(pName), fCubicVolume(0), fSurfaceArea(0)
34 {
35  CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz);
36 }
37 
38 void UTrd::CheckAndSetAllParameters ( double pdx1, double pdx2,
39  double pdy1, double pdy2,
40  double pdz )
41 {
42  if (pdx1 > 0 && pdx2 > 0 && pdy1 > 0 && pdy2 > 0 && pdz > 0)
43  {
44  fDx1 = pdx1;
45  fDx2 = pdx2;
46  fDy1 = pdy1;
47  fDy2 = pdy2;
48  fDz = pdz;
49  }
50  else
51  {
52  if (pdx1 >= 0 && pdx2 >= 0 && pdy1 >= 0 && pdy2 >= 0 && pdz >= 0)
53  {
54  // double Minimum_length= (1+per_thousand) * VUSolid::fgTolerance/2.;
55  // FIX-ME : temporary solution for ZERO or very-small parameters
56  //
57  double Minimum_length = fgTolerance;
58  fDx1 = std::max(pdx1, Minimum_length);
59  fDx2 = std::max(pdx2, Minimum_length);
60  fDy1 = std::max(pdy1, Minimum_length);
61  fDy2 = std::max(pdy2, Minimum_length);
62  fDz = std::max(pdz, Minimum_length);
63  }
64  else
65  {
66  std::cout << "ERROR - UTrd()::CheckAndSetAllParameters(): " << GetName()
67  << std::endl
68  << " Invalid dimensions, some are < 0 !" << std::endl
69  << " X - " << pdx1 << ", " << pdx2 << std::endl
70  << " Y - " << pdy1 << ", " << pdy2 << std::endl
71  << " Z - " << pdz << std::endl;
72  UUtils::Exception("UTrd::CheckAndSetAllParameters()", "InvalidSetup", FatalErrorInArguments, 1, "Invalid parameters.");
73  }
74  }
75 }
76 
77 void UTrd::SetAllParameters ( double pdx1, double pdx2, double pdy1,
78  double pdy2, double pdz )
79 {
80  CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz);
81 }
82 
83 double UTrd::SafetyFromInside(const UVector3& p, bool aAccurate) const
84 {
85  if (aAccurate) return SafetyFromInsideAccurate(p);
86 
87  double safe, zbase, tanxz, xdist, saf1, tanyz, ydist, saf2;
88 
89 #ifdef UDEBUG
90  if (Inside(p) == kOutside)
91  {
92  int oldprc = std::cout.precision(16) ;
93  std::cout << std::endl ;
94  DumpInfo();
95  std::cout << "Position:" << std::endl << std::endl ;
96  std::cout << "p.x() = " << p.x() / mm << " mm" << std::endl ;
97  std::cout << "p.y() = " << p.y() / mm << " mm" << std::endl ;
98  std::cout << "p.z() = " << p.z() / mm << " mm" << std::endl << std::endl ;
99  std::cout.precision(oldprc) ;
100  UUtils::Exception("G4Trd::DistanceToOut(p)", "Notification", Warning, 1,
101  "Point p is outside !?");
102  }
103 #endif
104 
105  safe = fDz - std::fabs(p.z); // z perpendicular Dist
106 
107  zbase = fDz + p.z;
108 
109  // xdist = distance perpendicular to z axis to closest x plane from p
110  // = (x half width of shape at p.z) - std::abs(p.x)
111  //
112  tanxz = (fDx2 - fDx1) * 0.5 / fDz; // angle between two parts on trinngle, which reflects proportional part on triangle related to position of point, see it on picture
113  xdist = fDx1 + tanxz * zbase - std::fabs(p.x); // distance to closest point on x border on x-axis on which point is located to
114  saf1 = xdist / std::sqrt(1.0 + tanxz * tanxz); // x*std::cos(ang_xz) =
115  // shortest (perpendicular)
116  // distance to plane, see picture
117  tanyz = (fDy2 - fDy1) * 0.5 / fDz;
118  ydist = fDy1 + tanyz * zbase - std::fabs(p.y);
119  saf2 = ydist / std::sqrt(1.0 + tanyz * tanyz);
120 
121  // Return minimum x/y/z distance
122  //
123  if (safe > saf1) safe = saf1;
124  if (safe > saf2) safe = saf2;
125 
126  if (safe < 0) safe = 0;
127  return safe;
128 }
129 
130 
131 
133 {
134  // computes the closest distance from given point to this shape, according
135  // to option. The matching point on the shape is stored in spoint.
136 
137  double saf[3];
138  //--- Compute safety first
139  // check Z facettes
140  saf[0] = fDz - std::abs(p.z);
141  double fx = 0.5 * (fDx1 - fDx2) / fDz;
142  double calf = 1. / std::sqrt(1.0 + fx * fx);
143  // check X facettes
144  double distx = 0.5 * (fDx1 + fDx2) - fx * p.z;
145  if (distx < 0) saf[1] = UUtils::kInfinity;
146  else saf[1] = (distx - std::abs(p.x)) * calf;
147 
148  double fy = 0.5 * (fDy1 - fDy2) / fDz;
149  calf = 1. / std::sqrt(1.0 + fy * fy);
150  // check Y facettes
151  distx = 0.5 * (fDy1 + fDy2) - fy * p.z;
152  if (distx < 0) saf[2] = UUtils::kInfinity;
153  else saf[2] = (distx - std::abs(p.y)) * calf;
154 
155  return amin(3, saf);
156 }
157 
158 /*
159 * Estimates the isotropic safety from a point outside the current solid to any
160 * of its surfaces. The algorithm may be accurate or should provide a fast
161 * underestimate.
162 * Note: In geant4, this method is equivalent to DistanceToIn, without given direction
163 * ______________________________________________________________________________
164 *
165 * Calculate distance (<= actual) to closest surface of shape from outside
166 * - Calculate distance to radial plane
167 * - Return 0 if point inside
168 * OK
169 */
170 double UTrd::SafetyFromOutside(const UVector3& p, bool aAccurate) const
171 {
172  if (aAccurate) return SafetyFromOutsideAccurate(p);
173 
174  double safe = 0.0;
175  double tanxz, distx, safx;
176  double tanyz, disty, safy;
177  double zbase;
178 
179  safe = std::abs(p.z) - fDz;
180  if (safe < 0) safe = 0; // Also used to ensure x/y distances
181  // POSITIVE
182  zbase = fDz + p.z;
183 
184  // Find distance along x direction to closest x plane
185  //
186  tanxz = (fDx2 - fDx1) * 0.5 / fDz;
187  // widx=fDx1+tanxz*(fDz+p.z()); // x width at p.z
188  // distx=std::abs(p.x())-widx; // distance to plane
189  distx = std::abs(p.x) - (fDx1 + tanxz * zbase); // distance to point on border of trd, related to axis on which point p lies
190  if (distx > safe)
191  {
192  safx = distx / std::sqrt(1.0 + tanxz * tanxz); // perpendicular distance calculation; vector Dist=Dist*std::cos(ang), it can be probably negative, then comparing in next statement, we will get rid of such distance, because it directs away from the solid
193  if (safx > safe) safe = safx;
194  }
195 
196  // Find distance along y direction to slanted wall
197  tanyz = (fDy2 - fDy1) * 0.5 / fDz;
198  // widy=fDy1+tanyz*(fDz+p.z()); // y width at p.z
199  // disty=std::abs(p.y())-widy; // distance to plane
200  disty = std::abs(p.y) - (fDy1 + tanyz * zbase);
201  if (disty > safe)
202  {
203  safy = disty / std::sqrt(1.0 + tanyz * tanyz); // distance along vector
204  if (safy > safe) safe = safy;
205  }
206  return safe;
207 }
208 
210 {
211  // computes the closest distance from given point to this shape, according
212  // to option. The matching point on the shape is stored in spoint.
213 
214  double saf[3];
215  //--- Compute safety first
216  // check Z facettes
217  saf[0] = fDz - std::abs(p.z);
218  double fx = 0.5 * (fDx1 - fDx2) / fDz;
219  double calf = 1. / std::sqrt(1.0 + fx * fx);
220  // check X facettes
221  double distx = 0.5 * (fDx1 + fDx2) - fx * p.z;
222  if (distx < 0) saf[1] = UUtils::kInfinity;
223  else saf[1] = (distx - std::abs(p.x)) * calf;
224 
225  double fy = 0.5 * (fDy1 - fDy2) / fDz;
226  calf = 1. / std::sqrt(1.0 + fy * fy);
227  // check Y facettes
228  distx = 0.5 * (fDy1 + fDy2) - fy * p.z;
229  if (distx < 0) saf[2] = UUtils::kInfinity;
230  else saf[2] = (distx - std::abs(p.y)) * calf;
231 
232  for (int i = 0; i < 3; i++) saf[i] = -saf[i];
233  return amax(3, saf);
234 }
235 
236 //______________________________________________________________________________
247 // ok
248 // Geant4 routine was selected because:
249 // 1. it is able to detect if the point is on surface
250 // 2. it counts with tolerance
251 // 3. algorithm checks are similar:
252 // Geant4: x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz - fgTolerance/2;
253 // Root: double dx = 0.5*(fDx2*(point.z+fDz)+fDx1*(fDz-point.z))/fDz;
254 //
256 {
258  double x, y, zbase1, zbase2;
259 
260  if (std::abs(p.z) <= fDz - fgTolerance / 2)
261  {
262  zbase1 = p.z + fDz; // Dist from -ve z plane
263  zbase2 = fDz - p.z; // Dist from +ve z plane
264 
265  // Check whether inside x tolerance
266  //
267  x = 0.5 * (fDx2 * zbase1 + fDx1 * zbase2) / fDz - fgTolerance / 2; // calculate x coordinate of one corner point corresponding to p.z inside trd (on x axis), ... by using proportional calculation related to triangle
268  if (std::abs(p.x) <= x)
269  {
270  y = 0.5 * ((fDy2 * zbase1 + fDy1 * zbase2)) / fDz - fgTolerance / 2;
271  if (std::abs(p.y) <= y)
272  {
273  in = eInside;
274  }
275  else if (std::abs(p.y) <= y + fgTolerance)
276  {
277  in = eSurface;
278  }
279  }
280  else if (std::abs(p.x) <= x + fgTolerance)
281  {
282  // y = y half width of shape at z of point + tolerant boundary
283  //
284  y = 0.5 * ((fDy2 * zbase1 + fDy1 * zbase2)) / fDz + fgTolerance / 2;
285  if (std::abs(p.y) <= y)
286  {
287  in = eSurface;
288  }
289  }
290  }
291  else if (std::abs(p.z) <= fDz + fgTolerance / 2)
292  {
293  // Only need to check outer tolerant boundaries
294  //
295  zbase1 = p.z + fDz; // Dist from -ve z plane
296  zbase2 = fDz - p.z; // Dist from +ve z plane
297 
298  // x = x half width of shape at z of point plus tolerance
299  //
300  x = 0.5 * (fDx2 * zbase1 + fDx1 * zbase2) / fDz + fgTolerance / 2;
301  if (std::abs(p.x) <= x)
302  {
303  // y = y half width of shape at z of point
304  //
305  y = 0.5 * ((fDy2 * zbase1 + fDy1 * zbase2)) / fDz + fgTolerance / 2;
306  if (std::abs(p.y) <= y) in = eSurface;
307  }
308  }
309  return in;
310 }
311 
312 
313 /*
314 * Computes distance from a point presumably outside the solid to the solid
315 * surface. Ignores first surface if the point is actually inside. Early return
316 * infinity in case the safety to any surface is found greater than the proposed
317 * step aPstep.
318 * The normal vector to the crossed surface is filled only in case the Orb is
319 * crossed, otherwise aNormal.IsNull() is true.
320 */
321 double UTrd::DistanceToIn(const UVector3& p,
322  const UVector3& v,
323  // UVector3 &aNormal,
324  double) const
325 {
326  double snxt = UUtils::kInfinity; // snxt = default return value
327  double smin, smax;
328  double s1, s2, tanxz, tanyz, ds1, ds2;
329  double ss1, ss2, sn1 = 0., sn2 = 0., dist;
330 
331  if (v.z) // Calculate valid z intersect range
332  {
333  if (v.z > 0) // Calculate smax: must be +ve or no intersection.
334  {
335  dist = fDz - p.z ; // to plane at +dz
336 
337  if (dist >= 0.5 * VUSolid::fgTolerance)
338  {
339  smax = dist / v.z; // distance to intersection with +dz, maximum
340  smin = -(fDz + p.z) / v.z; // distance to intersection with +dz, minimum
341  }
342  else return snxt ;
343  }
344  else // v.z <0
345  {
346  dist = fDz + p.z; // plane at -dz
347 
348  if (dist >= 0.5 * VUSolid::fgTolerance)
349  {
350  smax = -dist / v.z;
351  smin = (fDz - p.z) / v.z;
352  }
353  else return snxt ;
354  }
355  if (smin < 0) smin = 0 ;
356  }
357  else // v.z=0
358  {
359  if (std::abs(p.z) >= fDz) return snxt ; // Outside & no intersect
360  else
361  {
362  smin = 0 ; // Always inside z range
363  smax = UUtils::kInfinity;
364  }
365  }
366 
367  // Calculate x intersection range
368  //
369  // Calc half width at p.z, and components towards planes
370 
371  tanxz = (fDx2 - fDx1) * 0.5 / fDz ;
372  s1 = 0.5 * (fDx1 + fDx2) + tanxz * p.z ; // x half width at p.z
373  ds1 = v.x - tanxz * v.z ; // Components of v towards faces at +-x
374  ds2 = v.x + tanxz * v.z ;
375  ss1 = s1 - p.x; // -delta x to +ve plane
376  // -ve when outside
377  ss2 = -s1 - p.x; // -delta x to -ve plane
378  // +ve when outside
379  if (ss1 < 0 && ss2 <= 0)
380  {
381  if (ds1 < 0) // In +ve coord Area
382  {
383  sn1 = ss1 / ds1 ;
384 
385  if (ds2 < 0) sn2 = ss2 / ds2 ;
386  else sn2 = UUtils::kInfinity ;
387  }
388  else return snxt ;
389  }
390  else if (ss1 >= 0 && ss2 > 0)
391  {
392  if (ds2 > 0) // In -ve coord Area
393  {
394  sn1 = ss2 / ds2 ;
395 
396  if (ds1 > 0) sn2 = ss1 / ds1 ;
397  else sn2 = UUtils::kInfinity;
398 
399  }
400  else return snxt ;
401  }
402  else if (ss1 >= 0 && ss2 <= 0)
403  {
404  // Inside Area - calculate leaving distance
405  // *Don't* use exact distance to side for tolerance
406  // = ss1*std::cos(ang xz)
407  // = ss1/std::sqrt(1.0+tanxz*tanxz)
408  sn1 = 0 ;
409 
410  if (ds1 > 0)
411  {
412  if (ss1 > 0.5 * VUSolid::fgTolerance) sn2 = ss1 / ds1 ; // Leave +ve side extent
413  else return snxt ; // Leave immediately by +ve
414  }
415  else sn2 = UUtils::kInfinity ;
416 
417  if (ds2 < 0)
418  {
419  if (ss2 < -0.5 * VUSolid::fgTolerance)
420  {
421  dist = ss2 / ds2 ; // Leave -ve side extent
422  if (dist < sn2) sn2 = dist ;
423  }
424  else return snxt ;
425  }
426  }
427  else if (ss1 < 0 && ss2 > 0)
428  {
429  // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0)
430 
431  if (ds1 >= 0 || ds2 <= 0)
432  {
433  return snxt ;
434  }
435  else // Will intersect & stay inside
436  {
437  sn1 = ss1 / ds1 ;
438  dist = ss2 / ds2 ;
439  if (dist > sn1) sn1 = dist ;
440  sn2 = UUtils::kInfinity ;
441  }
442  }
443 
444  // Reduce allowed range of distances as appropriate
445 
446  if (sn1 > smin) smin = sn1 ;
447  if (sn2 < smax) smax = sn2 ;
448 
449  // Check for incompatible ranges (eg z intersects between 50 ->100 and x
450  // only 10-40 -> no intersection)
451 
452  if (smax < smin) return snxt ;
453 
454  // Calculate valid y intersection range
455  // (repeat of x intersection code)
456 
457  tanyz = (fDy2 - fDy1) * 0.5 / fDz ;
458  s2 = 0.5 * (fDy1 + fDy2) + tanyz * p.z; // y half width at p.z
459  ds1 = v.y - tanyz * v.z; // Components of v towards faces at +-y
460  ds2 = v.y + tanyz * v.z;
461  ss1 = s2 - p.y; // -delta y to +ve plane
462  ss2 = -s2 - p.y; // -delta y to -ve plane
463 
464  if (ss1 < 0 && ss2 <= 0)
465  {
466  if (ds1 < 0) // In +ve coord Area
467  {
468  sn1 = ss1 / ds1 ;
469  if (ds2 < 0) sn2 = ss2 / ds2 ;
470  else sn2 = UUtils::kInfinity ;
471  }
472  else return snxt ;
473  }
474  else if (ss1 >= 0 && ss2 > 0)
475  {
476  if (ds2 > 0) // In -ve coord Area
477  {
478  sn1 = ss2 / ds2 ;
479  if (ds1 > 0) sn2 = ss1 / ds1 ;
480  else sn2 = UUtils::kInfinity ;
481  }
482  else return snxt ;
483  }
484  else if (ss1 >= 0 && ss2 <= 0)
485  {
486  // Inside Area - calculate leaving distance
487  // *Don't* use exact distance to side for tolerance
488  // = ss1*std::cos(ang yz)
489  // = ss1/std::sqrt(1.0+tanyz*tanyz)
490  sn1 = 0 ;
491 
492  if (ds1 > 0)
493  {
494  if (ss1 > 0.5 * VUSolid::fgTolerance) sn2 = ss1 / ds1 ; // Leave +ve side extent
495  else return snxt ; // Leave immediately by +ve
496  }
497  else sn2 = UUtils::kInfinity ;
498 
499  if (ds2 < 0)
500  {
501  if (ss2 < -0.5 * VUSolid::fgTolerance)
502  {
503  dist = ss2 / ds2 ; // Leave -ve side extent
504  if (dist < sn2) sn2 = dist;
505  }
506  else return snxt ;
507  }
508  }
509  else if (ss1 < 0 && ss2 > 0)
510  {
511  // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0)
512 
513  if (ds1 >= 0 || ds2 <= 0)
514  {
515  return snxt ;
516  }
517  else // Will intersect & stay inside
518  {
519  sn1 = ss1 / ds1 ;
520  dist = ss2 / ds2 ;
521  if (dist > sn1) sn1 = dist ;
522  sn2 = UUtils::kInfinity ;
523  }
524  }
525 
526  // Reduce allowed range of distances as appropriate
527 
528  if (sn1 > smin) smin = sn1 ;
529  if (sn2 < smax) smax = sn2 ;
530 
531  // Check for incompatible ranges (eg x intersects between 50 ->100 and y
532  // only 10-40 -> no intersection). Set snxt if ok
533 
534  if (smax > smin) snxt = smin ;
535  if (snxt < 0.5 * VUSolid::fgTolerance) snxt = 0.0 ;
536 
537  return snxt ;
538 }
539 
540 
541 
542 //_____________________________________________________________________________
543 
544 /*
545 * Computes distance from a point presumably intside the solid to the solid
546 * surface. Ignores first surface along each axis systematically (for points
547 * inside or outside. Early returns zero in case the second surface is behind
548 * the starting point.
549 * o The proposed step is ignored.
550 * o The normal vector to the crossed surface is always filled.
551 * ______________________________________________________________________________
552 */
553 
554 
555 double UTrd::DistanceToOut(const UVector3& p, const UVector3& v,
556  UVector3& n, bool& aConvex, double /*aPstep*/) const
557 {
558  ESide side = kUndefined, snside = kUndefined;
559  aConvex = true;
560  double snxt, pdist;
561  double central, ss1, ss2, ds1, ds2, sn = 0., sn2 = 0.;
562  double tanxz = 0., cosxz = 0., tanyz = 0., cosyz = 0.;
563 
564  // Calculate z plane intersection
565  if (v.z > 0)
566  {
567  pdist = fDz - p.z;
568  if (pdist > VUSolid::fgTolerance / 2)
569  {
570  snxt = pdist / v.z;
571  side = kPZ;
572  }
573  else
574  {
575  n = UVector3(0, 0, 1);
576  return snxt = 0;
577  }
578  }
579  else if (v.z < 0)
580  {
581  pdist = fDz + p.z;
582  if (pdist > VUSolid::fgTolerance / 2)
583  {
584  snxt = -pdist / v.z;
585  side = kMZ;
586  }
587  else
588  {
589  // if (calcNorm)
590  {
591  n = UVector3(0, 0, -1);
592  }
593  return snxt = 0;
594  }
595  }
596  else
597  {
598  snxt = UUtils::kInfinity;
599  }
600 
601  //
602  // Calculate x intersection
603  //
604  tanxz = (fDx2 - fDx1) * 0.5 / fDz;
605  central = 0.5 * (fDx1 + fDx2);
606 
607  // +ve plane (1)
608  //
609  ss1 = central + tanxz * p.z - p.x; // distance || x axis to plane
610  // (+ve if point inside)
611  ds1 = v.x - tanxz * v.z; // component towards plane at +x
612  // (-ve if +ve -> -ve direction)
613  // -ve plane (2)
614  //
615  ss2 = -tanxz * p.z - p.x - central; //distance || x axis to plane
616  // (-ve if point inside)
617  ds2 = tanxz * v.z + v.x; // component towards plane at -x
618 
619  if (ss1 > 0 && ss2 < 0)
620  {
621  // Normal case - entirely inside region
622  if (ds1 <= 0 && ds2 < 0)
623  {
624  if (ss2 < -VUSolid::fgTolerance / 2)
625  {
626  sn = ss2 / ds2; // Leave by -ve side
627  snside = kMX;
628  }
629  else
630  {
631  sn = 0; // Leave immediately by -ve side
632  snside = kMX;
633  }
634  }
635  else if (ds1 > 0 && ds2 >= 0)
636  {
637  if (ss1 > VUSolid::fgTolerance / 2)
638  {
639  sn = ss1 / ds1; // Leave by +ve side
640  snside = kPX;
641  }
642  else
643  {
644  sn = 0; // Leave immediately by +ve side
645  snside = kPX;
646  }
647  }
648  else if (ds1 > 0 && ds2 < 0)
649  {
650  if (ss1 > VUSolid::fgTolerance / 2)
651  {
652  // sn=ss1/ds1; // Leave by +ve side
653  if (ss2 < -VUSolid::fgTolerance / 2)
654  {
655  sn = ss1 / ds1; // Leave by +ve side
656  sn2 = ss2 / ds2;
657  if (sn2 < sn)
658  {
659  sn = sn2;
660  snside = kMX;
661  }
662  else
663  {
664  snside = kPX;
665  }
666  }
667  else
668  {
669  sn = 0; // Leave immediately by -ve
670  snside = kMX;
671  }
672  }
673  else
674  {
675  sn = 0; // Leave immediately by +ve side
676  snside = kPX;
677  }
678  }
679  else
680  {
681  // Must be || to both
682  //
683  sn = UUtils::kInfinity; // Don't leave by either side
684  }
685  }
686  else if (ss1 <= 0 && ss2 < 0)
687  {
688  // Outside, in +ve Area
689 
690  if (ds1 > 0)
691  {
692  sn = 0; // Away from shape
693  // Left by +ve side
694  snside = kPX;
695  }
696  else
697  {
698  if (ds2 < 0)
699  {
700  // Ignore +ve plane and use -ve plane intersect
701  //
702  sn = ss2 / ds2; // Leave by -ve side
703  snside = kMX;
704  }
705  else
706  {
707  // Must be || to both -> exit determined by other axes
708  //
709  sn = UUtils::kInfinity; // Don't leave by either side
710  }
711  }
712  }
713  else if (ss1 > 0 && ss2 >= 0)
714  {
715  // Outside, in -ve Area
716 
717  if (ds2 < 0)
718  {
719  sn = 0; // away from shape
720  // Left by -ve side
721  snside = kMX;
722  }
723  else
724  {
725  if (ds1 > 0)
726  {
727  // Ignore +ve plane and use -ve plane intersect
728  //
729  sn = ss1 / ds1; // Leave by +ve side
730  snside = kPX;
731  }
732  else
733  {
734  // Must be || to both -> exit determined by other axes
735  //
736  sn = UUtils::kInfinity; // Don't leave by either side
737  }
738  }
739  }
740 
741  // Update minimum exit distance
742 
743  if (sn < snxt)
744  {
745  snxt = sn;
746  side = snside;
747  }
748  if (snxt > 0)
749  {
750  // Calculate y intersection
751  tanyz = (fDy2 - fDy1) * 0.5 / fDz;
752  central = 0.5 * (fDy1 + fDy2);
753 
754  // +ve plane (1)
755  //
756  ss1 = central + tanyz * p.z - p.y; // distance || y axis to plane
757  // (+ve if point inside)
758  ds1 = v.y - tanyz * v.z; // component towards +ve plane
759  // (-ve if +ve -> -ve direction)
760  // -ve plane (2)
761  //
762  ss2 = -tanyz * p.z - p.y - central; // distance || y axis to plane
763  // (-ve if point inside)
764  ds2 = tanyz * v.z + v.y; // component towards -ve plane
765 
766  if (ss1 > 0 && ss2 < 0)
767  {
768  // Normal case - entirely inside region
769 
770  if (ds1 <= 0 && ds2 < 0)
771  {
772  if (ss2 < -VUSolid::fgTolerance / 2)
773  {
774  sn = ss2 / ds2; // Leave by -ve side
775  snside = kMY;
776  }
777  else
778  {
779  sn = 0; // Leave immediately by -ve side
780  snside = kMY;
781  }
782  }
783  else if (ds1 > 0 && ds2 >= 0)
784  {
785  if (ss1 > VUSolid::fgTolerance / 2)
786  {
787  sn = ss1 / ds1; // Leave by +ve side
788  snside = kPY;
789  }
790  else
791  {
792  sn = 0; // Leave immediately by +ve side
793  snside = kPY;
794  }
795  }
796  else if (ds1 > 0 && ds2 < 0)
797  {
798  if (ss1 > VUSolid::fgTolerance / 2)
799  {
800  // sn=ss1/ds1; // Leave by +ve side
801  if (ss2 < -VUSolid::fgTolerance / 2)
802  {
803  sn = ss1 / ds1; // Leave by +ve side
804  sn2 = ss2 / ds2;
805  if (sn2 < sn)
806  {
807  sn = sn2;
808  snside = kMY;
809  }
810  else
811  {
812  snside = kPY;
813  }
814  }
815  else
816  {
817  sn = 0; // Leave immediately by -ve
818  snside = kMY;
819  }
820  }
821  else
822  {
823  sn = 0; // Leave immediately by +ve side
824  snside = kPY;
825  }
826  }
827  else
828  {
829  // Must be || to both
830  //
831  sn = UUtils::kInfinity; // Don't leave by either side
832  }
833  }
834  else if (ss1 <= 0 && ss2 < 0)
835  {
836  // Outside, in +ve Area
837 
838  if (ds1 > 0)
839  {
840  sn = 0; // Away from shape
841  // Left by +ve side
842  snside = kPY;
843  }
844  else
845  {
846  if (ds2 < 0)
847  {
848  // Ignore +ve plane and use -ve plane intersect
849  //
850  sn = ss2 / ds2; // Leave by -ve side
851  snside = kMY;
852  }
853  else
854  {
855  // Must be || to both -> exit determined by other axes
856  //
857  sn = UUtils::kInfinity; // Don't leave by either side
858  }
859  }
860  }
861  else if (ss1 > 0 && ss2 >= 0)
862  {
863  // Outside, in -ve Area
864  if (ds2 < 0)
865  {
866  sn = 0; // away from shape
867  // Left by -ve side
868  snside = kMY;
869  }
870  else
871  {
872  if (ds1 > 0)
873  {
874  // Ignore +ve plane and use -ve plane intersect
875  //
876  sn = ss1 / ds1; // Leave by +ve side
877  snside = kPY;
878  }
879  else
880  {
881  // Must be || to both -> exit determined by other axes
882  //
883  sn = UUtils::kInfinity; // Don't leave by either side
884  }
885  }
886  }
887 
888  // Update minimum exit distance
889 
890  if (sn < snxt)
891  {
892  snxt = sn;
893  side = snside;
894  }
895  }
896 
897  {
898  switch (side)
899  {
900  case kPX:
901  cosxz = 1.0 / std::sqrt(1.0 + tanxz * tanxz);
902  n.Set(cosxz, 0, -tanxz * cosxz);
903  break;
904  case kMX:
905  cosxz = -1.0 / std::sqrt(1.0 + tanxz * tanxz);
906  n.Set(cosxz, 0, tanxz * cosxz);
907  break;
908  case kPY:
909  cosyz = 1.0 / std::sqrt(1.0 + tanyz * tanyz);
910  n.Set(0, cosyz, -tanyz * cosyz);
911  break;
912  case kMY:
913  cosyz = -1.0 / std::sqrt(1.0 + tanyz * tanyz);
914  n.Set(0, cosyz, tanyz * cosyz);
915  break;
916  case kPZ:
917  n.Set(0, 0, 1);
918  break;
919  case kMZ:
920  n.Set(0, 0, -1);
921  break;
922  default:
923  UUtils::Exception("G4Trd::DistanceToOut(p,v,..)", "Notification", Warning, 1, "Undefined side for valid surface normal to solid.");
924  break;
925  }
926  }
927 
928  return snxt;
929 }
930 
931 
932 bool UTrd::Normal(const UVector3& p, UVector3& norm) const
933 {
934  UVector3 sumnorm(0., 0., 0.);
935  int noSurfaces = 0;
936  double z = 2.0 * fDz;
937  double delta = 0.5 * fgTolerance;
938 
939  double tanx = (fDx2 - fDx1) / z;
940  double secx = std::sqrt(1.0 + tanx * tanx);
941  double newpx = std::abs(p.x) - p.z * tanx;
942  double widx = fDx2 - fDz * tanx;
943 
944  double tany = (fDy2 - fDy1) / z;
945  double secy = std::sqrt(1.0 + tany * tany);
946  double newpy = std::abs(p.y) - p.z * tany;
947  double widy = fDy2 - fDz * tany;
948 
949  double distx = std::abs(newpx - widx) / secx; // perp. distance to x side
950  double disty = std::abs(newpy - widy) / secy; // to y side
951  double distz = std::abs(std::abs(p.z) - fDz); // to z side
952 
953  if (distx <= delta) // we are near the x corner?
954  {
955  double fcos = 1.0 / secx;
956  noSurfaces ++;
957  if (p.x >= 0.) sumnorm.x += fcos;
958  else sumnorm.x -= fcos;
959  sumnorm.z -= tanx * fcos;
960  }
961  if (disty <= delta) // y corner ...
962  {
963  double fcos = 1.0 / secy;
964  noSurfaces++;
965  if (p.y >= 0.) sumnorm.y += fcos;
966  else sumnorm.y -= fcos;
967  sumnorm.z -= tany * fcos;
968  }
969  if (distz <= delta)
970  {
971  noSurfaces++;
972  sumnorm.z += (p.z >= 0.) ? 1.0 : -1.0;
973  }
974  if (noSurfaces == 0)
975  {
976 #ifdef UDEBUG
977  UUtils::Exception("G4Trd::SurfaceNormal(p)", "Notification", Warning, 1, "Point p is not on surface !?");
978 #endif
979  // point is not on surface, calculate normal closest to surface. this is not likely to be used too often..., normally the user gives point on surface...
980  // norm = ApproxSurfaceNormal(p);
981 
982  // find closest side
983  //
984  double fcos;
985 
986  if (distx <= disty)
987  {
988  if (distx <= distz)
989  {
990  // Closest to X
991  //
992  fcos = 1.0 / secx;
993  // normal=(+/-std::cos(ang),0,-std::sin(ang))
994  if (p.x >= 0)
995  norm = UVector3(fcos, 0, -tanx * fcos);
996  else
997  norm = UVector3(-fcos, 0, -tanx * fcos);
998  }
999  else
1000  {
1001  // Closest to Z
1002  //
1003  if (p.z >= 0)
1004  norm = UVector3(0, 0, 1);
1005  else
1006  norm = UVector3(0, 0, -1);
1007  }
1008  }
1009  else
1010  {
1011  if (disty <= distz)
1012  {
1013  // Closest to Y
1014  //
1015  fcos = 1.0 / secy;
1016  if (p.y >= 0)
1017  norm = UVector3(0, fcos, -tany * fcos);
1018  else
1019  norm = UVector3(0, -fcos, -tany * fcos);
1020  }
1021  else
1022  {
1023  // Closest to Z
1024  //
1025  if (p.z >= 0)
1026  norm = UVector3(0, 0, 1);
1027  else
1028  norm = UVector3(0, 0, -1);
1029  }
1030  }
1031  }
1032  else
1033  {
1034  norm = sumnorm;
1035  // if we added to the vector more than once, we have to normalize the vector
1036  if (noSurfaces > 1) norm.Normalize();
1037  }
1038  return noSurfaces > 0; // (Inside(p) == eSurface);
1039 }
1040 
1041 
1043 //
1044 // Algorithm for SurfaceNormal() following the original specification
1045 // for points not on the surface
1046 
1048 {
1049  UVector3 norm;
1050  double z, tanx, secx, newpx, widx;
1051  double tany, secy, newpy, widy;
1052  double distx, disty, distz, fcos;
1053 
1054  z = 2.0 * fDz;
1055 
1056  tanx = (fDx2 - fDx1) / z;
1057  secx = std::sqrt(1.0 + tanx * tanx);
1058  newpx = std::abs(p.x) - p.z * tanx;
1059  widx = fDx2 - fDz * tanx;
1060 
1061  tany = (fDy2 - fDy1) / z;
1062  secy = std::sqrt(1.0 + tany * tany);
1063  newpy = std::abs(p.y) - p.z * tany;
1064  widy = fDy2 - fDz * tany;
1065 
1066  distx = std::abs(newpx - widx) / secx; // perpendicular distance to x side
1067  disty = std::abs(newpy - widy) / secy; // to y side
1068  distz = std::abs(std::abs(p.z) - fDz); // to z side
1069 
1070  // find closest side
1071  //
1072  if (distx <= disty)
1073  {
1074  if (distx <= distz)
1075  {
1076  // Closest to X
1077  //
1078  fcos = 1.0 / secx;
1079  // normal=(+/-std::cos(ang),0,-std::sin(ang))
1080  if (p.x >= 0)
1081  norm = UVector3(fcos, 0, -tanx * fcos);
1082  else
1083  norm = UVector3(-fcos, 0, -tanx * fcos);
1084  }
1085  else
1086  {
1087  // Closest to Z
1088  //
1089  if (p.z >= 0)
1090  norm = UVector3(0, 0, 1);
1091  else
1092  norm = UVector3(0, 0, -1);
1093  }
1094  }
1095  else
1096  {
1097  if (disty <= distz)
1098  {
1099  // Closest to Y
1100  //
1101  fcos = 1.0 / secy;
1102  if (p.y >= 0)
1103  norm = UVector3(0, fcos, -tany * fcos);
1104  else
1105  norm = UVector3(0, -fcos, -tany * fcos);
1106  }
1107  else
1108  {
1109  // Closest to Z
1110  //
1111  if (p.z >= 0)
1112  norm = UVector3(0, 0, 1);
1113  else
1114  norm = UVector3(0, 0, -1);
1115  }
1116  }
1117  return norm;
1118 }
1119 
1120 
1121 
1126 void UTrd::Extent(UVector3& aMin, UVector3& aMax) const
1127 {
1128  aMin.Set(-std::max(fDx1, fDx2), -std::max(fDy1, fDy2), -fDz);
1129  aMax.Set(std::max(fDx1, fDx2), std::max(fDy1, fDy2), fDz);
1130 }
1131 
1132 
1134 //
1135 // Stream object contents to an output stream
1136 
1137 std::ostream& UTrd::StreamInfo(std::ostream& os) const
1138 {
1139  int oldprc = os.precision(16);
1140  os << "-----------------------------------------------------------\n"
1141  << " *** Dump for solid - " << GetName() << " ***\n"
1142  << " ===================================================\n"
1143  << " Solid type: UTrd\n"
1144  << " Parameters: \n"
1145  << " half length X, surface -dZ: " << fDx1 << " mm \n"
1146  << " half length X, surface +dZ: " << fDx2 << " mm \n"
1147  << " half length Y, surface -dZ: " << fDy1 << " mm \n"
1148  << " half length Y, surface +dZ: " << fDy2 << " mm \n"
1149  << " half length Z : " << fDz << " mm \n"
1150  << "-----------------------------------------------------------\n";
1151  os.precision(oldprc);
1152 
1153  return os;
1154 }
1155 
1157 //
1158 // GetPointOnSurface
1159 //
1160 // Return a point (UVector3) randomly and uniformly
1161 // selected on the solid surface
1162 
1164 {
1165  double px, py, pz, tgX, tgY, secX, secY, select, sumS, tmp;
1166  double Sxy1, Sxy2, Sxy, Sxz, Syz;
1167 
1168  tgX = 0.5 * (fDx2 - fDx1) / fDz;
1169  secX = std::sqrt(1 + tgX * tgX);
1170  tgY = 0.5 * (fDy2 - fDy1) / fDz;
1171  secY = std::sqrt(1 + tgY * tgY);
1172 
1173  // calculate 0.25 of side surfaces, sumS is 0.25 of total surface
1174 
1175  Sxy1 = fDx1 * fDy1;
1176  Sxy2 = fDx2 * fDy2;
1177  Sxy = Sxy1 + Sxy2;
1178  Sxz = (fDx1 + fDx2) * fDz * secY;
1179  Syz = (fDy1 + fDy2) * fDz * secX;
1180  sumS = Sxy + Sxz + Syz;
1181 
1182  select = sumS * UUtils::Random();
1183 
1184  if (select < Sxy) // Sxy1 or Sxy2
1185  {
1186  if (select < Sxy1)
1187  {
1188  pz = -fDz;
1189  px = -fDx1 + 2 * fDx1 * UUtils::Random();
1190  py = -fDy1 + 2 * fDy1 * UUtils::Random();
1191  }
1192  else
1193  {
1194  pz = fDz;
1195  px = -fDx2 + 2 * fDx2 * UUtils::Random();
1196  py = -fDy2 + 2 * fDy2 * UUtils::Random();
1197  }
1198  }
1199  else if ((select - Sxy) < Sxz) // Sxz
1200  {
1201  pz = -fDz + 2 * fDz * UUtils::Random();
1202  tmp = fDx1 + (pz + fDz) * tgX;
1203  px = -tmp + 2 * tmp * UUtils::Random();
1204  tmp = fDy1 + (pz + fDz) * tgY;
1205 
1206  if (UUtils::Random() > 0.5)
1207  {
1208  py = tmp;
1209  }
1210  else
1211  {
1212  py = -tmp;
1213  }
1214  }
1215  else // Syz
1216  {
1217  pz = -fDz + 2 * fDz * UUtils::Random();
1218  tmp = fDy1 + (pz + fDz) * tgY;
1219  py = -tmp + 2 * tmp * UUtils::Random();
1220  tmp = fDx1 + (pz + fDz) * tgX;
1221 
1222  if (UUtils::Random() > 0.5)
1223  {
1224  px = tmp;
1225  }
1226  else
1227  {
1228  px = -tmp;
1229  }
1230  }
1231  return UVector3(px, py, pz);
1232 }
1233 
1235 //
1236 // Copy constructor
1237 
1238 UTrd::UTrd(const UTrd& rhs)
1239  : VUSolid(rhs), fDx1(rhs.fDx1), fDx2(rhs.fDx2),
1240  fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fDz),
1241  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea)
1242 {
1243 }
1244 
1246 //
1247 // Assignment operator
1248 
1250 {
1251  // Check assignment to self
1252  //
1253  if (this == &rhs)
1254  {
1255  return *this;
1256  }
1257 
1258 
1259  // Copy data
1260  //
1261  fDx1 = rhs.fDx1;
1262  fDx2 = rhs.fDx2;
1263  fDy1 = rhs.fDy1;
1264  fDy2 = rhs.fDy2;
1265  fDz = rhs.fDz;
1266  fCubicVolume = rhs.fCubicVolume;
1267  fSurfaceArea = rhs.fSurfaceArea;
1268 
1269  return *this;
1270 
1271 }
1273 //
1274 // Get Parameters List for Visualisation
1275 
1276 void UTrd::GetParametersList(int, double* aArray) const
1277 {
1278  aArray[0] = GetXHalfLength1();
1279  aArray[1] = GetXHalfLength2();
1280  aArray[2] = GetYHalfLength1();
1281  aArray[3] = GetYHalfLength2();
1282  aArray[4] = GetZHalfLength();
1283 }
1285 //
1286 // Clone
1287 
1289 {
1290  return new UTrd(*this);
1291 }
1293 //
1294 // Get Entity Type
1295 
1297 {
1298  return "Trd";
1299 }
double fDy2
Definition: UTrd.hh:108
virtual double SafetyFromInside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UTrd.cc:83
double GetZHalfLength() const
VUSolid * Clone() const
Definition: UTrd.cc:1288
virtual bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const
Definition: UTrd.cc:932
G4double z
Definition: TRTMaterials.hh:39
const std::string & GetName() const
Definition: VUSolid.hh:103
double SafetyFromOutsideAccurate(const UVector3 &aPoint) const
Definition: UTrd.cc:209
static double fgTolerance
Definition: VUSolid.hh:30
double Normalize()
Definition: UVector3.cc:89
UVector3 GetPointOnSurface() const
Definition: UTrd.cc:1163
double GetXHalfLength2() const
double fDy1
Definition: UTrd.hh:108
UVector3 ApproxSurfaceNormal(const UVector3 &p) const
Definition: UTrd.cc:1047
virtual void GetParametersList(int, double *) const
Definition: UTrd.cc:1276
ESide
Definition: UTrd.hh:30
double fDx2
Definition: UTrd.hh:108
double fCubicVolume
Definition: UTrd.hh:109
G4double fcos(G4double arg)
double x
Definition: UVector3.hh:136
const G4int smax
double SafetyFromInsideAccurate(const UVector3 &aPoint) const
Definition: UTrd.cc:132
void CheckAndSetAllParameters(double pdx1, double pdx2, double pdy1, double pdy2, double pdz)
Definition: UTrd.cc:38
double amax(int n, const double *a) const
static const double kInfinity
Definition: UUtils.hh:53
double GetYHalfLength2() const
void SetAllParameters(double pdx1, double pdx2, double pdy1, double pdy2, double pdz)
Definition: UTrd.cc:77
double fDz
Definition: UTrd.hh:108
double amin(int n, const double *a) const
double fDx1
Definition: UTrd.hh:108
double GetYHalfLength1() const
double fSurfaceArea
Definition: UTrd.hh:110
EnumInside
Definition: VUSolid.hh:23
const G4int n
EnumInside Inside(const UVector3 &aPoint) const
Return whether point inside/outside/on surface Split into radius checks.
Definition: UTrd.cc:255
void Set(double xx, double yy, double zz)
Definition: UVector3.hh:245
void Extent(UVector3 &aMin, UVector3 &aMax) const
Returns the full 3D cartesian extent of the solid.
Definition: UTrd.cc:1126
UTrd & operator=(const UTrd &rhs)
Definition: UTrd.cc:1249
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: UTrd.hh:28
double GetXHalfLength1() const
std::string UGeometryType
Definition: UTypes.hh:70
UGeometryType GetEntityType() const
Definition: UTrd.cc:1296
virtual double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
Definition: UTrd.cc:555
void Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
double z
Definition: UVector3.hh:138
UTrd()
Definition: UTrd.hh:32
std::ostream & StreamInfo(std::ostream &os) const
Definition: UTrd.cc:1137
virtual double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UTrd.cc:170
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
static const double mm
Definition: G4SIunits.hh:102
double y
Definition: UVector3.hh:137
virtual double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const
Definition: UTrd.cc:321