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