Geant4  10.01
UTrap.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 // UTrap
12 //
13 // 12.02.13 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include "UUtils.hh"
18 #include <string>
19 #include <cmath>
20 #include <sstream>
21 #include "UTrap.hh"
22 
23 using namespace std;
24 
26 //
27 // Accuracy of coplanarity
28 
29 const double kCoplanar_Tolerance = 1E-4 ;
30 
32 //
33 // Private enum: Not for external use
34 
35 enum Eside {kUndef, ks0, ks1, ks2, ks3, kPZ, kMZ};
36 
38 //
39 // Constructor - check and Set half-widths as well as angles:
40 // final check of coplanarity
41 
42 UTrap::UTrap(const std::string& pName,
43  double pDz,
44  double pTheta, double pPhi,
45  double pDy1, double pDx1, double pDx2,
46  double pAlp1,
47  double pDy2, double pDx3, double pDx4,
48  double pAlp2)
49  : VUSolid(pName)
50 {
51  if (pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 ||
52  pDx2 <= 0 || pDy2 <= 0 || pDx3 <= 0 || pDx4 <= 0)
53  {
54  std::ostringstream message;
55  message << "Invalid length parameters for Solid: " << GetName() << std::endl
56  << " X - "
57  << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << std::endl
58  << " Y - " << pDy1 << ", " << pDy2 << std::endl
59  << " Z - " << pDz;
60  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
61  FatalErrorInArguments, 1, message.str().c_str());
62  }
63 
64  fDz = pDz;
65  fTthetaCphi = std::tan(pTheta) * std::cos(pPhi);
66  fTthetaSphi = std::tan(pTheta) * std::sin(pPhi);
67 
68  fDy1 = pDy1;
69  fDx1 = pDx1;
70  fDx2 = pDx2;
71  fTalpha1 = std::tan(pAlp1);
72 
73  fDy2 = pDy2;
74  fDx3 = pDx3;
75  fDx4 = pDx4;
76  fTalpha2 = std::tan(pAlp2);
77 
78  MakePlanes();
79  fCubicVolume=0.;
80  fSurfaceArea=0.;
81 }
82 
84 //
85 // Constructor - Design of trapezoid based on 8 UVector3 parameters,
86 // which are its vertices. Checking of planarity with preparation of
87 // fPlanes[] and than calculation of other members
88 
89 UTrap::UTrap(const std::string& pName,
90  const UVector3 pt[8])
91  : VUSolid(pName)
92 {
93  SetPlanes(pt);
94  fCubicVolume=0.;
95  fSurfaceArea=0.;
96 }
97 
99 //
100 // Constructor for Right Angular Wedge from STEP
101 
102 UTrap::UTrap(const std::string& pName,
103  double pZ,
104  double pY,
105  double pX, double pLTX)
106  : VUSolid(pName)
107 {
108  bool good;
109 
110  if (pZ <= 0 || pY <= 0 || pX <= 0 || pLTX <= 0 || pLTX > pX)
111  {
112  std::ostringstream message;
113  message << "Invalid length parameters for Solid: " << GetName();
114  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
115  FatalErrorInArguments, 1, message.str().c_str());
116  }
117 
118  fDz = 0.5 * pZ ;
119  fTthetaCphi = 0 ;
120  fTthetaSphi = 0 ;
121 
122  fDy1 = 0.5 * pY;
123  fDx1 = 0.5 * pX ;
124  fDx2 = 0.5 * pLTX;
125  fTalpha1 = 0.5 * (pLTX - pX) / pY;
126 
127  fDy2 = fDy1 ;
128  fDx3 = fDx1;
129  fDx4 = fDx2 ;
130  fTalpha2 = fTalpha1 ;
131 
132  UVector3 pt[8] ;
133 
134  pt[0] = UVector3(-fDz * fTthetaCphi - fDy1 * fTalpha1 - fDx1,
135  -fDz * fTthetaSphi - fDy1, -fDz);
136  pt[1] = UVector3(-fDz * fTthetaCphi - fDy1 * fTalpha1 + fDx1,
137  -fDz * fTthetaSphi - fDy1, -fDz);
138  pt[2] = UVector3(-fDz * fTthetaCphi + fDy1 * fTalpha1 - fDx2,
139  -fDz * fTthetaSphi + fDy1, -fDz);
140  pt[3] = UVector3(-fDz * fTthetaCphi + fDy1 * fTalpha1 + fDx2,
141  -fDz * fTthetaSphi + fDy1, -fDz);
142  pt[4] = UVector3(+fDz * fTthetaCphi - fDy2 * fTalpha2 - fDx3,
143  +fDz * fTthetaSphi - fDy2, +fDz);
144  pt[5] = UVector3(+fDz * fTthetaCphi - fDy2 * fTalpha2 + fDx3,
145  +fDz * fTthetaSphi - fDy2, +fDz);
146  pt[6] = UVector3(+fDz * fTthetaCphi + fDy2 * fTalpha2 - fDx4,
147  +fDz * fTthetaSphi + fDy2, +fDz);
148  pt[7] = UVector3(+fDz * fTthetaCphi + fDy2 * fTalpha2 + fDx4,
149  +fDz * fTthetaSphi + fDy2, +fDz);
150 
151  // Bottom side with normal approx. -Y
152  //
153  good = MakePlane(pt[0], pt[4], pt[5], pt[1], fPlanes[0]);
154  if (!good)
155  {
156  std::ostringstream message;
157  message << "Face at ~-Y not planar for Solid: " << GetName();
158  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
159  FatalErrorInArguments, 1, message.str().c_str());
160  }
161 
162  // Top side with normal approx. +Y
163  //
164  good = MakePlane(pt[2], pt[3], pt[7], pt[6], fPlanes[1]);
165  if (!good)
166  {
167  std::ostringstream message;
168  message << "Face at ~+Y not planar for Solid: " << GetName();
169  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
170  FatalErrorInArguments, 1, message.str().c_str());
171  }
172 
173  // Front side with normal approx. -X
174  //
175  good = MakePlane(pt[0], pt[2], pt[6], pt[4], fPlanes[2]);
176  if (!good)
177  {
178  std::ostringstream message;
179  message << "Face at ~-X not planar for Solid: " << GetName();
180  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
181  FatalErrorInArguments, 1, message.str().c_str());
182  }
183 
184  // Back side iwth normal approx. +X
185  //
186  good = MakePlane(pt[1], pt[5], pt[7], pt[3], fPlanes[3]);
187  if (!good)
188  {
189  std::ostringstream message;
190  message << "Face at ~+X not planar for Solid: " << GetName();
191  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
192  FatalError, 1, message.str().c_str());
193  }
194  fCubicVolume=0.;
195  fSurfaceArea=0.;
196 }
197 
199 //
200 // Constructor for UTrd
201 
202 UTrap::UTrap(const std::string& pName,
203  double pDx1, double pDx2,
204  double pDy1, double pDy2,
205  double pDz)
206  : VUSolid(pName)
207 {
208  bool good;
209 
210  if (pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 || pDx2 <= 0 || pDy2 <= 0)
211  {
212  std::ostringstream message;
213  message << "Invalid length parameters for Solid: " << GetName();
214  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
215  FatalErrorInArguments, 1, message.str().c_str());
216  }
217 
218  fDz = pDz;
219  fTthetaCphi = 0 ;
220  fTthetaSphi = 0 ;
221 
222  fDy1 = pDy1 ;
223  fDx1 = pDx1 ;
224  fDx2 = pDx1 ;
225  fTalpha1 = 0 ;
226 
227  fDy2 = pDy2 ;
228  fDx3 = pDx2 ;
229  fDx4 = pDx2 ;
230  fTalpha2 = 0 ;
231 
232  UVector3 pt[8] ;
233 
234  pt[0] = UVector3(-fDz * fTthetaCphi - fDy1 * fTalpha1 - fDx1,
235  -fDz * fTthetaSphi - fDy1, -fDz);
236  pt[1] = UVector3(-fDz * fTthetaCphi - fDy1 * fTalpha1 + fDx1,
237  -fDz * fTthetaSphi - fDy1, -fDz);
238  pt[2] = UVector3(-fDz * fTthetaCphi + fDy1 * fTalpha1 - fDx2,
239  -fDz * fTthetaSphi + fDy1, -fDz);
240  pt[3] = UVector3(-fDz * fTthetaCphi + fDy1 * fTalpha1 + fDx2,
241  -fDz * fTthetaSphi + fDy1, -fDz);
242  pt[4] = UVector3(+fDz * fTthetaCphi - fDy2 * fTalpha2 - fDx3,
243  +fDz * fTthetaSphi - fDy2, +fDz);
244  pt[5] = UVector3(+fDz * fTthetaCphi - fDy2 * fTalpha2 + fDx3,
245  +fDz * fTthetaSphi - fDy2, +fDz);
246  pt[6] = UVector3(+fDz * fTthetaCphi + fDy2 * fTalpha2 - fDx4,
247  +fDz * fTthetaSphi + fDy2, +fDz);
248  pt[7] = UVector3(+fDz * fTthetaCphi + fDy2 * fTalpha2 + fDx4,
249  +fDz * fTthetaSphi + fDy2, +fDz);
250 
251  // Bottom side with normal approx. -Y
252  //
253  good = MakePlane(pt[0], pt[4], pt[5], pt[1], fPlanes[0]);
254  if (!good)
255  {
256  std::ostringstream message;
257  message << "Face at ~-Y not planar for Solid: " << GetName();
258  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
259  FatalErrorInArguments, 1, message.str().c_str());
260  }
261 
262  // Top side with normal approx. +Y
263  //
264  good = MakePlane(pt[2], pt[3], pt[7], pt[6], fPlanes[1]);
265  if (!good)
266  {
267  std::ostringstream message;
268  message << "Face at ~+Y not planar for Solid: " << GetName();
269  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002", FatalErrorInArguments, 1, message.str().c_str());
270  }
271 
272  // Front side with normal approx. -X
273  //
274  good = MakePlane(pt[0], pt[2], pt[6], pt[4], fPlanes[2]);
275  if (!good)
276  {
277  std::ostringstream message;
278  message << "Face at ~-X not planar for Solid: " << GetName();
279  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
280  FatalErrorInArguments, 1, message.str().c_str());
281  }
282 
283  // Back side iwth normal approx. +X
284  //
285  good = MakePlane(pt[1], pt[5], pt[7], pt[3], fPlanes[3]);
286  if (!good)
287  {
288  std::ostringstream message;
289  message << "Face at ~+X not planar for Solid: " << GetName();
290  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
291  FatalErrorInArguments, 1, message.str().c_str());
292  }
293  fCubicVolume=0.;
294  fSurfaceArea=0.;
295 }
296 
298 //
299 // Constructor for UPara
300 
301 UTrap::UTrap(const std::string& pName,
302  double pDx, double pDy,
303  double pDz,
304  double pAlpha,
305  double pTheta, double pPhi)
306  : VUSolid(pName)
307 {
308  bool good;
309 
310  if (pDz <= 0 || pDy <= 0 || pDx <= 0)
311  {
312  std::ostringstream message;
313  message << "Invalid length parameters for Solid: " << GetName();
314  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
315  FatalErrorInArguments, 1, message.str().c_str());
316  }
317 
318  fDz = pDz ;
319  fTthetaCphi = std::tan(pTheta) * std::cos(pPhi) ;
320  fTthetaSphi = std::tan(pTheta) * std::sin(pPhi) ;
321 
322  fDy1 = pDy ;
323  fDx1 = pDx ;
324  fDx2 = pDx ;
325  fTalpha1 = std::tan(pAlpha) ;
326 
327  fDy2 = pDy ;
328  fDx3 = pDx ;
329  fDx4 = pDx ;
330  fTalpha2 = fTalpha1 ;
331 
332  UVector3 pt[8] ;
333 
334  pt[0] = UVector3(-fDz * fTthetaCphi - fDy1 * fTalpha1 - fDx1,
335  -fDz * fTthetaSphi - fDy1, -fDz);
336  pt[1] = UVector3(-fDz * fTthetaCphi - fDy1 * fTalpha1 + fDx1,
337  -fDz * fTthetaSphi - fDy1, -fDz);
338  pt[2] = UVector3(-fDz * fTthetaCphi + fDy1 * fTalpha1 - fDx2,
339  -fDz * fTthetaSphi + fDy1, -fDz);
340  pt[3] = UVector3(-fDz * fTthetaCphi + fDy1 * fTalpha1 + fDx2,
341  -fDz * fTthetaSphi + fDy1, -fDz);
342  pt[4] = UVector3(+fDz * fTthetaCphi - fDy2 * fTalpha2 - fDx3,
343  +fDz * fTthetaSphi - fDy2, +fDz);
344  pt[5] = UVector3(+fDz * fTthetaCphi - fDy2 * fTalpha2 + fDx3,
345  +fDz * fTthetaSphi - fDy2, +fDz);
346  pt[6] = UVector3(+fDz * fTthetaCphi + fDy2 * fTalpha2 - fDx4,
347  +fDz * fTthetaSphi + fDy2, +fDz);
348  pt[7] = UVector3(+fDz * fTthetaCphi + fDy2 * fTalpha2 + fDx4,
349  +fDz * fTthetaSphi + fDy2, +fDz);
350 
351  // Bottom side with normal approx. -Y
352  //
353  good = MakePlane(pt[0], pt[4], pt[5], pt[1], fPlanes[0]);
354  if (!good)
355  {
356  std::ostringstream message;
357  message << "Face at ~-Y not planar for Solid: " << GetName();
358  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
359  FatalErrorInArguments, 1, message.str().c_str());
360  }
361 
362  // Top side with normal approx. +Y
363  //
364  good = MakePlane(pt[2], pt[3], pt[7], pt[6], fPlanes[1]);
365  if (!good)
366  {
367  std::ostringstream message;
368  message << "Face at ~+Y not planar for Solid: " << GetName();
369  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
370  FatalErrorInArguments, 1, message.str().c_str());
371  }
372 
373  // Front side with normal approx. -X
374  //
375  good = MakePlane(pt[0], pt[2], pt[6], pt[4], fPlanes[2]);
376  if (!good)
377  {
378  std::ostringstream message;
379  message << "Face at ~-X not planar for Solid: " << GetName();
380  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
381  FatalErrorInArguments, 1, message.str().c_str());
382  }
383 
384  // Back side iwth normal approx. +X
385  //
386  good = MakePlane(pt[1], pt[5], pt[7], pt[3], fPlanes[3]);
387  if (!good)
388  {
389  std::ostringstream message;
390  message << "Face at ~+X not planar for Solid: " << GetName();
391  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
392  FatalErrorInArguments, 1, message.str().c_str());
393  }
394  fCubicVolume=0.;
395  fSurfaceArea=0.;
396 }
397 
399 //
400 // Nominal constructor for UTrap whose parameters are to be Set by
401 // a UVParamaterisation later. Check and Set half-widths as well as
402 // angles: final check of coplanarity
403 
404 UTrap::UTrap(const std::string& pName)
405  : VUSolid(pName), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
406  fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
407  fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
408 {
409  MakePlanes();
410  fCubicVolume=0.;
411  fSurfaceArea=0.;
412 }
413 
414 
416 //
417 // Destructor
418 
419 UTrap::~UTrap()
420 {
421 }
422 
424 //
425 // Copy constructor
426 
427 UTrap::UTrap(const UTrap& rhs)
428  : VUSolid(rhs), fDz(rhs.fDz),
429  fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
430  fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
431  fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
432 {
433  for (size_t i = 0; i < 4; ++i)
434  {
435  fPlanes[i].a = rhs.fPlanes[i].a;
436  fPlanes[i].b = rhs.fPlanes[i].b;
437  fPlanes[i].c = rhs.fPlanes[i].c;
438  fPlanes[i].d = rhs.fPlanes[i].d;
439  }
440  fCubicVolume=rhs.fCubicVolume;
441  fSurfaceArea=rhs.fSurfaceArea;
442 }
443 
445 //
446 // Assignment operator
447 
448 UTrap& UTrap::operator = (const UTrap& rhs)
449 {
450  // Check assignment to self
451  //
452  if (this == &rhs)
453  {
454  return *this;
455  }
456 
457  // Copy base class data
458  //
459  VUSolid::operator=(rhs);
460 
461  // Copy data
462  //
463  fDz = rhs.fDz;
464  fTthetaCphi = rhs.fTthetaCphi;
465  fTthetaSphi = rhs.fTthetaSphi;
466  fDy1 = rhs.fDy1;
467  fDx1 = rhs.fDx1;
468  fDx2 = rhs.fDx2;
469  fTalpha1 = rhs.fTalpha1;
470  fDy2 = rhs.fDy2;
471  fDx3 = rhs.fDx3;
472  fDx4 = rhs.fDx4;
473  fTalpha2 = rhs.fTalpha2;
474  for (size_t i = 0; i < 4; ++i)
475  {
476  fPlanes[i].a = rhs.fPlanes[i].a;
477  fPlanes[i].b = rhs.fPlanes[i].b;
478  fPlanes[i].c = rhs.fPlanes[i].c;
479  fPlanes[i].d = rhs.fPlanes[i].d;
480  }
481  fCubicVolume=rhs.fCubicVolume;
482  fSurfaceArea=rhs.fSurfaceArea;
483  return *this;
484 }
485 
487 //
488 // Set all parameters, as for constructor - check and Set half-widths
489 // as well as angles: final check of coplanarity
490 
491 void UTrap::SetAllParameters(double pDz,
492  double pTheta,
493  double pPhi,
494  double pDy1,
495  double pDx1,
496  double pDx2,
497  double pAlp1,
498  double pDy2,
499  double pDx3,
500  double pDx4,
501  double pAlp2)
502 {
503  if (pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 || pDx2 <= 0 || pDy2 <= 0 || pDx3 <= 0 || pDx4 <= 0)
504  {
505  std::ostringstream message;
506  message << "Invalid Length Parameters for Solid: " << GetName() << std::endl
507  << " X - "
508  << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << std::endl
509  << " Y - " << pDy1 << ", " << pDy2 << std::endl
510  << " Z - " << pDz;
511  UUtils::Exception("UTrap::SetAllParameters()", "GeomSolids0002",
512  FatalErrorInArguments, 1, message.str().c_str());
513  }
514  fCubicVolume = 0.;
515  fSurfaceArea = 0.;
516  fDz = pDz;
517  fTthetaCphi = std::tan(pTheta) * std::cos(pPhi);
518  fTthetaSphi = std::tan(pTheta) * std::sin(pPhi);
519 
520  fDy1 = pDy1;
521  fDx1 = pDx1;
522  fDx2 = pDx2;
523  fTalpha1 = std::tan(pAlp1);
524 
525  fDy2 = pDy2;
526  fDx3 = pDx3;
527  fDx4 = pDx4;
528  fTalpha2 = std::tan(pAlp2);
529 
530  MakePlanes();
531 }
532 
533 void UTrap::SetPlanes(const UVector3 pt[8])
534 {
535  bool good;
536 
537  // Start with check of centering - the center of gravity trap line
538  // should Cross the origin of frame
539 
540  if (!(pt[0].z < 0
541  && pt[0].z == pt[1].z && pt[0].z == pt[2].z
542  && pt[0].z == pt[3].z
543  && pt[4].z > 0
544  && pt[4].z == pt[5].z && pt[4].z == pt[6].z
545  && pt[4].z == pt[7].z
546  && std::fabs(pt[0].z + pt[4].z) < VUSolid::Tolerance()
547  && pt[0].y == pt[1].y && pt[2].y == pt[3].y
548  && pt[4].y == pt[5].y && pt[6].y == pt[7].y
549  && std::fabs(pt[0].y + pt[2].y + pt[4].y + pt[6].y) < VUSolid::Tolerance()
550  && std::fabs(pt[0].x + pt[1].x + pt[4].x + pt[5].x +
551  pt[2].x + pt[3].x + pt[6].x + pt[7].x) < VUSolid::Tolerance()))
552  {
553  std::ostringstream message;
554  message << "Invalid vertice coordinates for Solid: " << GetName();
555  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
556  FatalErrorInArguments, 1, message.str().c_str());
557  }
558 
559  // Bottom side with normal approx. -Y
560 
561  good = MakePlane(pt[0], pt[4], pt[5], pt[1], fPlanes[0]);
562 
563  if (!good)
564  {
565 
566  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002", FatalError, 1,
567  "Face at ~-Y not planar.");
568  }
569 
570  // Top side with normal approx. +Y
571 
572  good = MakePlane(pt[2], pt[3], pt[7], pt[6], fPlanes[1]);
573 
574  if (!good)
575  {
576  std::ostringstream message;
577  message << "Face at ~+Y not planar for Solid: " << GetName();
578  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
579  FatalErrorInArguments, 1, message.str().c_str());
580  }
581 
582  // Front side with normal approx. -X
583 
584  good = MakePlane(pt[0], pt[2], pt[6], pt[4], fPlanes[2]);
585 
586  if (!good)
587  {
588  std::ostringstream message;
589  message << "Face at ~-X not planar for Solid: " << GetName();
590  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002", FatalErrorInArguments, 1, message.str().c_str());
591  }
592 
593  // Back side iwth normal approx. +X
594 
595  good = MakePlane(pt[1], pt[5], pt[7], pt[3], fPlanes[3]);
596  if (!good)
597  {
598  std::ostringstream message;
599  message << "Face at ~+X not planar for Solid: " << GetName();
600  UUtils::Exception("UTrap::UTrap()", "GeomSolids0002",
601  FatalErrorInArguments, 1, message.str().c_str());
602  }
603  fDz = (pt[7]).z ;
604 
605  fDy1 = ((pt[2]).y - (pt[1]).y) * 0.5;
606  fDx1 = ((pt[1]).x - (pt[0]).x) * 0.5;
607  fDx2 = ((pt[3]).x - (pt[2]).x) * 0.5;
608  fTalpha1 = ((pt[2]).x + (pt[3]).x - (pt[1]).x - (pt[0]).x) * 0.25 / fDy1;
609 
610  fDy2 = ((pt[6]).y - (pt[5]).y) * 0.5;
611  fDx3 = ((pt[5]).x - (pt[4]).x) * 0.5;
612  fDx4 = ((pt[7]).x - (pt[6]).x) * 0.5;
613  fTalpha2 = ((pt[6]).x + (pt[7]).x - (pt[5]).x - (pt[4]).x) * 0.25 / fDy2;
614 
615  fTthetaCphi = ((pt[4]).x + fDy2 * fTalpha2 + fDx3) / fDz;
616  fTthetaSphi = ((pt[4]).y + fDy2) / fDz;
617 }
618 
620 //
621 // Checking of coplanarity
622 
623 bool UTrap::MakePlanes()
624 {
625  bool good = true;
626 
627  UVector3 pt[8] ;
628 
629  pt[0] = UVector3(-fDz * fTthetaCphi - fDy1 * fTalpha1 - fDx1,
630  -fDz * fTthetaSphi - fDy1, -fDz);
631  pt[1] = UVector3(-fDz * fTthetaCphi - fDy1 * fTalpha1 + fDx1,
632  -fDz * fTthetaSphi - fDy1, -fDz);
633  pt[2] = UVector3(-fDz * fTthetaCphi + fDy1 * fTalpha1 - fDx2,
634  -fDz * fTthetaSphi + fDy1, -fDz);
635  pt[3] = UVector3(-fDz * fTthetaCphi + fDy1 * fTalpha1 + fDx2,
636  -fDz * fTthetaSphi + fDy1, -fDz);
637  pt[4] = UVector3(+fDz * fTthetaCphi - fDy2 * fTalpha2 - fDx3,
638  +fDz * fTthetaSphi - fDy2, +fDz);
639  pt[5] = UVector3(+fDz * fTthetaCphi - fDy2 * fTalpha2 + fDx3,
640  +fDz * fTthetaSphi - fDy2, +fDz);
641  pt[6] = UVector3(+fDz * fTthetaCphi + fDy2 * fTalpha2 - fDx4,
642  +fDz * fTthetaSphi + fDy2, +fDz);
643  pt[7] = UVector3(+fDz * fTthetaCphi + fDy2 * fTalpha2 + fDx4,
644  +fDz * fTthetaSphi + fDy2, +fDz);
645 
646  // Bottom side with normal approx. -Y
647  //
648  good = MakePlane(pt[0], pt[4], pt[5], pt[1], fPlanes[0]) ;
649  if (!good)
650  {
651  std::ostringstream message;
652  message << "Face at ~-Y not planar for Solid: " << GetName();
653  UUtils::Exception("UTrap::MakePlanes()", "GeomSolids0002",
654  FatalError, 1, message.str().c_str());
655  }
656 
657  // Top side with normal approx. +Y
658  //
659  good = MakePlane(pt[2], pt[3], pt[7], pt[6], fPlanes[1]);
660  if (!good)
661  {
662  std::ostringstream message;
663  message << "Face at ~+Y not planar for Solid: " << GetName();
664  UUtils::Exception("UTrap::MakePlanes()", "GeomSolids0002",
665  FatalError, 1, message.str().c_str());
666  }
667 
668  // Front side with normal approx. -X
669  //
670  good = MakePlane(pt[0], pt[2], pt[6], pt[4], fPlanes[2]);
671  if (!good)
672  {
673  std::ostringstream message;
674  message << "Face at ~-X not planar for Solid: " << GetName();
675  UUtils::Exception("UTrap::MakePlanes()", "GeomSolids0002",
676  FatalError, 1, message.str().c_str());
677  }
678 
679  // Back side iwth normal approx. +X
680  //
681  good = MakePlane(pt[1], pt[5], pt[7], pt[3], fPlanes[3]);
682  if (!good)
683  {
684  std::ostringstream message;
685  message << "Face at ~+X not planar for Solid: " << GetName();
686  UUtils::Exception("UTrap::MakePlanes()", "GeomSolids0002",
687  FatalError, 1, message.str().c_str());
688  }
689 
690  return good;
691 }
692 
694 //
695 // Calculate the coef's of the plane p1->p2->p3->p4->p1
696 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed from
697 // infront of the plane (i.e. from normal direction).
698 //
699 // Return true if the ThreeVectors are coplanar + Set coef;s
700 // false if ThreeVectors are not coplanar
701 
702 bool UTrap::MakePlane(const UVector3& p1,
703  const UVector3& p2,
704  const UVector3& p3,
705  const UVector3& p4,
706  UTrapSidePlane& plane)
707 {
708  double a, b, c, sd;
709  UVector3 v12, v13, v14, Vcross;
710 
711  bool good;
712 
713  v12 = p2 - p1;
714  v13 = p3 - p1;
715  v14 = p4 - p1;
716  Vcross = v12.Cross(v13);
717 
718  if (std::fabs(Vcross.Dot(v14) / (Vcross.Mag()*v14.Mag())) > kCoplanar_Tolerance)
719  {
720  good = false;
721  }
722  else
723  {
724  // a,b,c correspond to the x/y/z components of the
725  // normal vector to the plane
726 
727  // a = (p2.y-p1.y)*(p1.z+p2.z)+(p3.y-p2.y)*(p2.z+p3.z);
728  // a += (p4.y-p3.y)*(p3.z+p4.z)+(p1.y-p4.y)*(p4.z+p1.z); // ?
729  // b = (p2.z-p1.z)*(p1.x+p2.x)+(p3.z-p2.z)*(p2.x+p3.x);
730  // b += (p4.z-p3.z)*(p3.x+p4.x)+(p1.z-p4.z)*(p4.x+p1.x); // ?
731  // c = (p2.x-p1.x)*(p1.y+p2.y)+(p3.x-p2.x)*(p2.y+p3.y);
732  // c += (p4.x-p3.x)*(p3.y+p4.y)+(p1.x-p4.x)*(p4.y+p1.y); // ?
733 
734  // Let create diagonals 4-2 and 3-1 than (4-2)x(3-1) provides
735  // vector perpendicular to the plane directed to outside !!!
736  // and a,b,c, = f(1,2,3,4) external relative to trap normal
737 
738  a = +(p4.y - p2.y) * (p3.z - p1.z)
739  - (p3.y - p1.y) * (p4.z - p2.z);
740 
741  b = -(p4.x - p2.x) * (p3.z - p1.z)
742  + (p3.x - p1.x) * (p4.z - p2.z);
743 
744  c = +(p4.x - p2.x) * (p3.y - p1.y)
745  - (p3.x - p1.x) * (p4.y - p2.y);
746 
747  sd = std::sqrt(a * a + b * b + c * c); // so now vector plane.(a,b,c) is Unit
748 
749  if (sd > 0)
750  {
751  plane.a = a / sd;
752  plane.b = b / sd;
753  plane.c = c / sd;
754  }
755  else
756  {
757  std::ostringstream message;
758  message << "Invalid parameters: norm.mod() <= 0, for Solid: "
759  << GetName();
760  UUtils::Exception("UTrap::MakePlanes()", "GeomSolids0002",
761  FatalError, 1, message.str().c_str()) ;
762  }
763  // Calculate D: p1 in in plane so D=-n.p1.Vect()
764 
765  plane.d = -(plane.a * p1.x + plane.b * p1.y + plane.c * p1.z);
766 
767  good = true;
768  }
769  return good;
770 }
771 
772 
773 
774 
776 //
777 // Return whether point inside/outside/on surface, using tolerance
778 
779 VUSolid::EnumInside UTrap::Inside(const UVector3& p) const
780 {
781  VUSolid::EnumInside in;
782  double Dist;
783  int i;
784  if (std::fabs(p.z) <= fDz - VUSolid::Tolerance() * 0.5)
785  {
786  in = eInside;
787 
788  for (i = 0; i < 4; i++)
789  {
790  Dist = fPlanes[i].a * p.x + fPlanes[i].b * p.y
791  + fPlanes[i].c * p.z + fPlanes[i].d;
792 
793  if (Dist > VUSolid::Tolerance() * 0.5) return in = eOutside;
794  else if (Dist > -VUSolid::Tolerance() * 0.5) in = eSurface;
795 
796  }
797  }
798  else if (std::fabs(p.z) <= fDz + VUSolid::Tolerance() * 0.5)
799  {
800  in = eSurface;
801 
802  for (i = 0; i < 4; i++)
803  {
804  Dist = fPlanes[i].a * p.x + fPlanes[i].b * p.y
805  + fPlanes[i].c * p.z + fPlanes[i].d;
806 
807  if (Dist > VUSolid::Tolerance() * 0.5) return in = eOutside;
808  }
809  }
810  else in = eOutside;
811 
812  return in;
813 }
814 
816 //
817 // Calculate side nearest to p, and return normal
818 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
819 
820 bool UTrap::Normal(const UVector3& p, UVector3& aNormal) const
821 {
822  int i, noSurfaces = 0;
823  double dist, distz, distx, disty, distmx, distmy, safe = UUtils::kInfinity;
824  double delta = 0.5 * VUSolid::Tolerance();
825  UVector3 norm, sumnorm(0., 0., 0.);
826 
827  for (i = 0; i < 4; i++)
828  {
829  dist = std::fabs(fPlanes[i].a * p.x + fPlanes[i].b * p.y
830  + fPlanes[i].c * p.z + fPlanes[i].d);
831  if (dist < safe)
832  {
833  safe = dist;
834  }
835  }
836  distz = std::fabs(std::fabs(p.z) - fDz);
837 
838  distmy = std::fabs(fPlanes[0].a * p.x + fPlanes[0].b * p.y
839  + fPlanes[0].c * p.z + fPlanes[0].d);
840 
841  disty = std::fabs(fPlanes[1].a * p.x + fPlanes[1].b * p.y
842  + fPlanes[1].c * p.z + fPlanes[1].d);
843 
844  distmx = std::fabs(fPlanes[2].a * p.x + fPlanes[2].b * p.y
845  + fPlanes[2].c * p.z + fPlanes[2].d);
846 
847  distx = std::fabs(fPlanes[3].a * p.x + fPlanes[3].b * p.y
848  + fPlanes[3].c * p.z + fPlanes[3].d);
849 
850  UVector3 nX = UVector3(fPlanes[3].a, fPlanes[3].b, fPlanes[3].c);
851  UVector3 nmX = UVector3(fPlanes[2].a, fPlanes[2].b, fPlanes[2].c);
852  UVector3 nY = UVector3(fPlanes[1].a, fPlanes[1].b, fPlanes[1].c);
853  UVector3 nmY = UVector3(fPlanes[0].a, fPlanes[0].b, fPlanes[0].c);
854  UVector3 nZ = UVector3(0., 0., 1.0);
855 
856  if (distx <= delta)
857  {
858  noSurfaces ++;
859  sumnorm += nX;
860  }
861  if (distmx <= delta)
862  {
863  noSurfaces ++;
864  sumnorm += nmX;
865  }
866  if (disty <= delta)
867  {
868  noSurfaces ++;
869  sumnorm += nY;
870  }
871  if (distmy <= delta)
872  {
873  noSurfaces ++;
874  sumnorm += nmY;
875  }
876  if (distz <= delta)
877  {
878  noSurfaces ++;
879  if (p.z >= 0.) sumnorm += nZ;
880  else sumnorm -= nZ;
881  }
882  if (noSurfaces == 0)
883  {
884 #ifdef UDEBUG
885  UUtils::Exception("UTrap::SurfaceNormal(p)", "GeomSolids1002",
886  Warning, 1, "Point p is not on surface !?");
887 #endif
888  norm = ApproxSurfaceNormal(p);
889  }
890  else if (noSurfaces == 1) norm = sumnorm;
891  else norm = sumnorm.Unit();
892  aNormal = norm;
893  return noSurfaces != 0;
894 }
895 
897 //
898 // Algorithm for SurfaceNormal() following the original specification
899 // for points not on the surface
900 
901 UVector3 UTrap::ApproxSurfaceNormal(const UVector3& p) const
902 {
903  double safe = UUtils::kInfinity, Dist, safez;
904  int i, imin = 0;
905  for (i = 0; i < 4; i++)
906  {
907  Dist = std::fabs(fPlanes[i].a * p.x + fPlanes[i].b * p.y
908  + fPlanes[i].c * p.z + fPlanes[i].d);
909  if (Dist < safe)
910  {
911  safe = Dist;
912  imin = i;
913  }
914  }
915  safez = std::fabs(std::fabs(p.z) - fDz);
916  if (safe < safez)
917  {
918  return UVector3(fPlanes[imin].a, fPlanes[imin].b, fPlanes[imin].c);
919  }
920  else
921  {
922  if (p.z > 0)
923  {
924  return UVector3(0, 0, 1);
925  }
926  else
927  {
928  return UVector3(0, 0, -1);
929  }
930  }
931 }
932 
934 //
935 // Calculate distance to shape from outside - return UUtils::kInfinity if no intersection
936 //
937 // ALGORITHM:
938 // For each component, calculate pair of minimum and maximum intersection
939 // values for which the particle is in the extent of the shape
940 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
941 
942 double UTrap::DistanceToIn(const UVector3& p,
943  const UVector3& v, double) const
944 {
945 
946  double snxt; // snxt = default return value
947  double max, smax, smin;
948  double pdist, Comp, vdist;
949  int i;
950  //
951  // Z Intersection range
952  //
953  if (v.z > 0)
954  {
955  max = fDz - p.z ;
956  if (max > 0.5 * VUSolid::Tolerance())
957  {
958  smax = max / v.z;
959  smin = (-fDz - p.z) / v.z;
960  }
961  else
962  {
963  return snxt = UUtils::kInfinity;
964  }
965  }
966  else if (v.z < 0)
967  {
968  max = - fDz - p.z ;
969  if (max < -0.5 * VUSolid::Tolerance())
970  {
971  smax = max / v.z;
972  smin = (fDz - p.z) / v.z;
973  }
974  else
975  {
976  return snxt = UUtils::kInfinity;
977  }
978  }
979  else
980  {
981  if (std::fabs(p.z) < fDz - 0.5 * VUSolid::Tolerance()) // Inside was <=fDz
982  {
983  smin = 0;
984  smax = UUtils::kInfinity;
985  }
986  else
987  {
988  return snxt = UUtils::kInfinity;
989  }
990  }
991 
992  for (i = 0; i < 4; i++)
993  {
994  pdist = fPlanes[i].a * p.x + fPlanes[i].b * p.y
995  + fPlanes[i].c * p.z + fPlanes[i].d;
996  Comp = fPlanes[i].a * v.x + fPlanes[i].b * v.y + fPlanes[i].c * v.z;
997  if (pdist >= -0.5 * VUSolid::Tolerance()) // was >0
998  {
999  //
1000  // Outside the plane -> this is an extent entry distance
1001  //
1002  if (Comp >= 0) // was >0
1003  {
1004  return snxt = UUtils::kInfinity ;
1005  }
1006  else
1007  {
1008  vdist = -pdist / Comp;
1009  if (vdist > smin)
1010  {
1011  if (vdist < smax)
1012  {
1013  smin = vdist;
1014  }
1015  else
1016  {
1017  return snxt = UUtils::kInfinity;
1018  }
1019  }
1020  }
1021  }
1022  else
1023  {
1024  //
1025  // Inside the plane -> couble be an extent exit distance (smax)
1026  //
1027  if (Comp > 0) // Will leave extent
1028  {
1029  vdist = -pdist / Comp;
1030  if (vdist < smax)
1031  {
1032  if (vdist > smin)
1033  {
1034  smax = vdist;
1035  }
1036  else
1037  {
1038  return snxt = UUtils::kInfinity;
1039  }
1040  }
1041  }
1042  }
1043  }
1044  //
1045  // Checks in non z plane intersections ensure smin<smax
1046  //
1047  if (smin >= 0)
1048  {
1049  snxt = smin ;
1050  }
1051  else
1052  {
1053  snxt = 0 ;
1054  }
1055  return snxt;
1056 }
1057 
1059 //
1060 // Calculate exact shortest distance to any boundary from outside
1061 // This is the best fast estimation of the shortest distance to trap
1062 // - Returns 0 is ThreeVector inside
1063 
1064 double UTrap::SafetyFromInside(const UVector3& p, bool) const
1065 {
1066  double safe = 0.0, Dist;
1067  int i;
1068  safe = std::fabs(p.z) - fDz;
1069  for (i = 0; i < 4; i++)
1070  {
1071  Dist = fPlanes[i].a * p.x + fPlanes[i].b * p.y
1072  + fPlanes[i].c * p.z + fPlanes[i].d;
1073  if (Dist > safe) safe = Dist;
1074  }
1075  if (safe < 0) safe = 0;
1076  return safe;
1077 }
1078 
1080 //
1081 // Calculate distance to surface of shape from inside
1082 // Calculate distance to x/y/z planes - smallest is exiting distance
1083 
1084 double UTrap::DistanceToOut(const UVector3& p, const UVector3& v, UVector3& aNormalVector, bool& aConvex, double) const
1085 {
1086  Eside side = kUndef;
1087  double snxt; // snxt = return value
1088  double pdist, Comp, vdist, max;
1089  //
1090  // Z Intersections
1091  //
1092  if (v.z > 0)
1093  {
1094  max = fDz - p.z;
1095  if (max > VUSolid::Tolerance() / 2)
1096  {
1097  snxt = max / v.z;
1098  side = kPZ;
1099  }
1100  else
1101  {
1102  aConvex = true;
1103  aNormalVector = UVector3(0, 0, 1);
1104  return snxt = 0;
1105  }
1106  }
1107  else if (v.z < 0)
1108  {
1109  max = -fDz - p.z;
1110  if (max < -VUSolid::Tolerance() / 2)
1111  {
1112  snxt = max / v.z;
1113  side = kMZ;
1114  }
1115  else
1116  {
1117  aConvex = true;
1118  aNormalVector = UVector3(0, 0, -1);
1119  return snxt = 0;
1120  }
1121  }
1122  else
1123  {
1124  snxt = UUtils::kInfinity;
1125  }
1126 
1127  //
1128  // Intersections with planes[0] (expanded because of setting enum)
1129  //
1130  pdist = fPlanes[0].a * p.x + fPlanes[0].b * p.y + fPlanes[0].c * p.z + fPlanes[0].d;
1131  Comp = fPlanes[0].a * v.x + fPlanes[0].b * v.y + fPlanes[0].c * v.z;
1132  if (pdist > 0)
1133  {
1134  // Outside the plane
1135  if (Comp > 0)
1136  {
1137  // Leaving immediately
1138  aConvex = true;
1139  aNormalVector = UVector3(fPlanes[0].a, fPlanes[0].b, fPlanes[0].c);
1140  return snxt = 0;
1141  }
1142  }
1143  else if (pdist < -VUSolid::Tolerance() / 2)
1144  {
1145  // Inside the plane
1146  if (Comp > 0)
1147  {
1148  // Will leave extent
1149  vdist = -pdist / Comp;
1150  if (vdist < snxt)
1151  {
1152  snxt = vdist;
1153  side = ks0;
1154  }
1155  }
1156  }
1157  else
1158  {
1159  // On surface
1160  if (Comp > 0)
1161  {
1162  aConvex = true;
1163  aNormalVector = UVector3(fPlanes[0].a, fPlanes[0].b, fPlanes[0].c);
1164  return snxt = 0;
1165  }
1166  }
1167 
1168  //
1169  // Intersections with planes[1] (expanded because of setting enum)
1170  //
1171  pdist = fPlanes[1].a * p.x + fPlanes[1].b * p.y + fPlanes[1].c * p.z + fPlanes[1].d;
1172  Comp = fPlanes[1].a * v.x + fPlanes[1].b * v.y + fPlanes[1].c * v.z;
1173  if (pdist > 0)
1174  {
1175  // Outside the plane
1176  if (Comp > 0)
1177  {
1178  // Leaving immediately
1179  aConvex = true;
1180  aNormalVector = UVector3(fPlanes[1].a, fPlanes[1].b, fPlanes[1].c);
1181  return snxt = 0;
1182  }
1183  }
1184  else if (pdist < -VUSolid::Tolerance() / 2)
1185  {
1186  // Inside the plane
1187  if (Comp > 0)
1188  {
1189  // Will leave extent
1190  vdist = -pdist / Comp;
1191  if (vdist < snxt)
1192  {
1193  snxt = vdist;
1194  side = ks1;
1195  }
1196  }
1197  }
1198  else
1199  {
1200  // On surface
1201  if (Comp > 0)
1202  {
1203  aConvex = true;
1204  aNormalVector = UVector3(fPlanes[1].a, fPlanes[1].b, fPlanes[1].c);
1205  return snxt = 0;
1206  }
1207  }
1208 
1209  //
1210  // Intersections with planes[2] (expanded because of setting enum)
1211  //
1212  pdist = fPlanes[2].a * p.x + fPlanes[2].b * p.y + fPlanes[2].c * p.z + fPlanes[2].d;
1213  Comp = fPlanes[2].a * v.x + fPlanes[2].b * v.y + fPlanes[2].c * v.z;
1214  if (pdist > 0)
1215  {
1216  // Outside the plane
1217  if (Comp > 0)
1218  {
1219  // Leaving immediately
1220  aConvex = true;
1221  aNormalVector = UVector3(fPlanes[2].a, fPlanes[2].b, fPlanes[2].c);
1222  return snxt = 0;
1223  }
1224  }
1225  else if (pdist < -VUSolid::Tolerance() / 2)
1226  {
1227  // Inside the plane
1228  if (Comp > 0)
1229  {
1230  // Will leave extent
1231  vdist = -pdist / Comp;
1232  if (vdist < snxt)
1233  {
1234  snxt = vdist;
1235  side = ks2;
1236  }
1237  }
1238  }
1239  else
1240  {
1241  // On surface
1242  if (Comp > 0)
1243  {
1244  aConvex = true;
1245  aNormalVector = UVector3(fPlanes[2].a, fPlanes[2].b, fPlanes[2].c);
1246  return snxt = 0;
1247  }
1248  }
1249 
1250  //
1251  // Intersections with planes[3] (expanded because of setting enum)
1252  //
1253  pdist = fPlanes[3].a * p.x + fPlanes[3].b * p.y + fPlanes[3].c * p.z + fPlanes[3].d;
1254  Comp = fPlanes[3].a * v.x + fPlanes[3].b * v.y + fPlanes[3].c * v.z;
1255  if (pdist > 0)
1256  {
1257  // Outside the plane
1258  if (Comp > 0)
1259  {
1260  // Leaving immediately
1261  aConvex = true;
1262  aNormalVector = UVector3(fPlanes[3].a, fPlanes[3].b, fPlanes[3].c);
1263  return snxt = 0;
1264  }
1265  }
1266  else if (pdist < -VUSolid::Tolerance() / 2)
1267  {
1268  // Inside the plane
1269  if (Comp > 0)
1270  {
1271  // Will leave extent
1272  vdist = -pdist / Comp;
1273  if (vdist < snxt)
1274  {
1275  snxt = vdist;
1276  side = ks3;
1277  }
1278  }
1279  }
1280  else
1281  {
1282  // On surface
1283  if (Comp > 0)
1284  {
1285  aConvex = true;
1286  aNormalVector = UVector3(fPlanes[3].a, fPlanes[3].b, fPlanes[3].c);
1287  return snxt = 0;
1288  }
1289  }
1290 
1291  // Set normal
1292  aConvex = true;
1293  switch (side)
1294  {
1295  case ks0:
1296  aNormalVector = UVector3(fPlanes[0].a, fPlanes[0].b, fPlanes[0].c);
1297  break;
1298  case ks1:
1299  aNormalVector = UVector3(fPlanes[1].a, fPlanes[1].b, fPlanes[1].c);
1300  break;
1301  case ks2:
1302  aNormalVector = UVector3(fPlanes[2].a, fPlanes[2].b, fPlanes[2].c);
1303  break;
1304  case ks3:
1305  aNormalVector = UVector3(fPlanes[3].a, fPlanes[3].b, fPlanes[3].c);
1306  break;
1307  case kMZ:
1308  aNormalVector = UVector3(0, 0, -1);
1309  break;
1310  case kPZ:
1311  aNormalVector = UVector3(0, 0, 1);
1312  break;
1313  default:
1314  cout << std::endl;
1315 
1316  std::ostringstream message;
1317  int oldprc = message.precision(16);
1318  message << "Undefined side for valid surface normal to solid."
1319  << std::endl
1320  << "Position:" << std::endl << std::endl
1321  << "p.x = " << p.x << " mm" << std::endl
1322  << "p.y = " << p.y << " mm" << std::endl
1323  << "p.z = " << p.z << " mm" << std::endl << std::endl
1324  << "Direction:" << std::endl << std::endl
1325  << "v.x = " << v.x << std::endl
1326  << "v.y = " << v.y << std::endl
1327  << "v.z = " << v.z << std::endl << std::endl
1328  << "Proposed distance :" << std::endl << std::endl
1329  << "snxt = " << snxt << " mm" << std::endl;
1330  message.precision(oldprc);
1331  UUtils::Exception("UTrap::DistanceToOut(p,v,..)", "GeomSolids1002",
1332  Warning, 1, message.str().c_str());
1333  break;
1334  }
1335  return snxt;
1336 }
1337 
1339 //
1340 // Calculate exact shortest distance to any boundary from inside
1341 // - Returns 0 is ThreeVector outside
1342 
1343 double UTrap::SafetyFromOutside(const UVector3& p, bool /*precise*/) const
1344 {
1345  double safe = 0.0, Dist;
1346  int i;
1347 
1348 #ifdef UDEBUG
1349  if (Inside(p) == eOutside)
1350  {
1351  int oldprc = cout.precision(16) ;
1352  cout << std::endl ;
1353 
1354  cout << "Position:" << std::endl << std::endl ;
1355  cout << "p.x = " << p.x << " mm" << std::endl ;
1356  cout << "p.y = " << p.y << " mm" << std::endl ;
1357  cout << "p.z = " << p.z << " mm" << std::endl << std::endl ;
1358  cout.precision(oldprc) ;
1359  UUtils::Exception("UTrap::DistanceToOut(p)",
1360  "GeomSolids1002", Warning, 1, "Point p is outside !?");
1361  }
1362 #endif
1363 
1364  safe = fDz - std::fabs(p.z);
1365  if (safe < 0) safe = 0;
1366  else
1367  {
1368  for (i = 0; i < 4; i++)
1369  {
1370  Dist = -(fPlanes[i].a * p.x + fPlanes[i].b * p.y
1371  + fPlanes[i].c * p.z + fPlanes[i].d);
1372  if (Dist < safe) safe = Dist;
1373  }
1374  if (safe < 0) safe = 0;
1375  }
1376  return safe;
1377 }
1378 
1379 
1381 //
1382 // GetEntityType
1383 
1384 UGeometryType UTrap::GetEntityType() const
1385 {
1386  return std::string("UTrap");
1387 }
1388 
1390 //
1391 // Make a clone of the object
1392 //
1393 VUSolid* UTrap::Clone() const
1394 {
1395  return new UTrap(*this);
1396 }
1397 
1399 //
1400 // Stream object contents to an output stream
1401 
1402 std::ostream& UTrap::StreamInfo(std::ostream& os) const
1403 {
1404  int oldprc = os.precision(16);
1405  os << "-----------------------------------------------------------\n"
1406  << " *** Dump for solid - " << GetName() << " ***\n"
1407  << " ===================================================\n"
1408  << " Solid type: UTrap\n"
1409  << " Parameters: \n"
1410  << " half length Z: " << fDz << " mm \n"
1411  << " half length Y of face -fDz: " << fDy1 << " mm \n"
1412  << " half length X of side -fDy1, face -fDz: " << fDx1 << " mm \n"
1413  << " half length X of side +fDy1, face -fDz: " << fDx2 << " mm \n"
1414  << " half length Y of face +fDz: " << fDy2 << " mm \n"
1415  << " half length X of side -fDy2, face +fDz: " << fDx3 << " mm \n"
1416  << " half length X of side +fDy2, face +fDz: " << fDx4 << " mm \n"
1417  << " std::tan(theta)*std::cos(phi): " << fTthetaCphi / (UUtils::kPi / 180.0) << " degrees \n"
1418  << " std::tan(theta)*std::sin(phi): " << fTthetaSphi / (UUtils::kPi / 180.0) << " degrees \n"
1419  << " std::tan(alpha), -fDz: " << fTalpha1 / (UUtils::kPi / 180.0) << " degrees \n"
1420  << " std::tan(alpha), +fDz: " << fTalpha2 / (UUtils::kPi / 180.0) << " degrees \n"
1421  << " trap side plane equations:\n"
1422  << " " << fPlanes[0].a << " X + " << fPlanes[0].b << " Y + "
1423  << fPlanes[0].c << " Z + " << fPlanes[0].d << " = 0\n"
1424  << " " << fPlanes[1].a << " X + " << fPlanes[1].b << " Y + "
1425  << fPlanes[1].c << " Z + " << fPlanes[1].d << " = 0\n"
1426  << " " << fPlanes[2].a << " X + " << fPlanes[2].b << " Y + "
1427  << fPlanes[2].c << " Z + " << fPlanes[2].d << " = 0\n"
1428  << " " << fPlanes[3].a << " X + " << fPlanes[3].b << " Y + "
1429  << fPlanes[3].c << " Z + " << fPlanes[3].d << " = 0\n"
1430  << "-----------------------------------------------------------\n";
1431  os.precision(oldprc);
1432 
1433  return os;
1434 }
1435 
1437 //
1438 // GetPointOnPlane
1439 //
1440 // Auxiliary method for Get Point on Surface
1441 
1442 UVector3 UTrap::GetPointOnPlane(UVector3 p0, UVector3 p1,
1443  UVector3 p2, UVector3 p3,
1444  double& area) const
1445 {
1446  double lambda1, lambda2, chose, aOne, aTwo;
1447  UVector3 t, u, v, w, Area, normal;
1448 
1449  t = p1 - p0;
1450  u = p2 - p1;
1451  v = p3 - p2;
1452  w = p0 - p3;
1453 
1454  Area = UVector3(w.y * v.z - w.z * v.y,
1455  w.z * v.x - w.x * v.z,
1456  w.x * v.y - w.y * v.x);
1457 
1458  aOne = 0.5 * Area.Mag();
1459 
1460  Area = UVector3(t.y * u.z - t.z * u.y,
1461  t.z * u.x - t.x * u.z,
1462  t.x * u.y - t.y * u.x);
1463 
1464  aTwo = 0.5 * Area.Mag();
1465 
1466  area = aOne + aTwo;
1467 
1468  chose = UUtils::Random(0., aOne + aTwo);
1469 
1470  if ((chose >= 0.) && (chose < aOne))
1471  {
1472  lambda1 = UUtils::Random(0., 1.);
1473  lambda2 = UUtils::Random(0., lambda1);
1474  return (p2 + lambda1 * v + lambda2 * w);
1475  }
1476 
1477  // else
1478 
1479  lambda1 = UUtils::Random(0., 1.);
1480  lambda2 = UUtils::Random(0., lambda1);
1481 
1482  return (p0 + lambda1 * t + lambda2 * u);
1483 }
1484 
1486 //
1487 // GetPointOnSurface
1488 
1489 UVector3 UTrap::GetPointOnSurface() const
1490 {
1491  double aOne, aTwo, aThree, aFour, aFive, aSix, chose;
1492  UVector3 One, Two, Three, Four, Five, Six, test;
1493  UVector3 pt[8];
1494 
1495  pt[0] = UVector3(-fDz * fTthetaCphi - fDy1 * fTalpha1 - fDx1,
1496  -fDz * fTthetaSphi - fDy1, -fDz);
1497  pt[1] = UVector3(-fDz * fTthetaCphi - fDy1 * fTalpha1 + fDx1,
1498  -fDz * fTthetaSphi - fDy1, -fDz);
1499  pt[2] = UVector3(-fDz * fTthetaCphi + fDy1 * fTalpha1 - fDx2,
1500  -fDz * fTthetaSphi + fDy1, -fDz);
1501  pt[3] = UVector3(-fDz * fTthetaCphi + fDy1 * fTalpha1 + fDx2,
1502  -fDz * fTthetaSphi + fDy1, -fDz);
1503  pt[4] = UVector3(+fDz * fTthetaCphi - fDy2 * fTalpha2 - fDx3,
1504  +fDz * fTthetaSphi - fDy2, +fDz);
1505  pt[5] = UVector3(+fDz * fTthetaCphi - fDy2 * fTalpha2 + fDx3,
1506  +fDz * fTthetaSphi - fDy2, +fDz);
1507  pt[6] = UVector3(+fDz * fTthetaCphi + fDy2 * fTalpha2 - fDx4,
1508  +fDz * fTthetaSphi + fDy2, +fDz);
1509  pt[7] = UVector3(+fDz * fTthetaCphi + fDy2 * fTalpha2 + fDx4,
1510  +fDz * fTthetaSphi + fDy2, +fDz);
1511 
1512  // make sure we provide the points in a clockwise fashion
1513 
1514  One = GetPointOnPlane(pt[0], pt[1], pt[3], pt[2], aOne);
1515  Two = GetPointOnPlane(pt[4], pt[5], pt[7], pt[6], aTwo);
1516  Three = GetPointOnPlane(pt[6], pt[7], pt[3], pt[2], aThree);
1517  Four = GetPointOnPlane(pt[4], pt[5], pt[1], pt[0], aFour);
1518  Five = GetPointOnPlane(pt[0], pt[2], pt[6], pt[4], aFive);
1519  Six = GetPointOnPlane(pt[1], pt[3], pt[7], pt[5], aSix);
1520 
1521  chose = UUtils::Random(0., aOne + aTwo + aThree + aFour + aFive + aSix);
1522  if ((chose >= 0.) && (chose < aOne))
1523  {
1524  return One;
1525  }
1526  else if ((chose >= aOne) && (chose < aOne + aTwo))
1527  {
1528  return Two;
1529  }
1530  else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aThree))
1531  {
1532  return Three;
1533  }
1534  else if ((chose >= aOne + aTwo + aThree) && (chose < aOne + aTwo + aThree + aFour))
1535  {
1536  return Four;
1537  }
1538  else if ((chose >= aOne + aTwo + aThree + aFour)
1539  && (chose < aOne + aTwo + aThree + aFour + aFive))
1540  {
1541  return Five;
1542  }
1543  return Six;
1544 }
1545 
1546 void UTrap::Extent(UVector3& aMin, UVector3& aMax) const
1547 {
1548  //Z axis
1549  aMin.z = -fDz;
1550  aMax.z = fDz;
1551 
1552  double temp[8] ; // some points for intersection with zMin/zMax
1553  UVector3 pt[8]; // vertices after translation
1554 
1555  //X axis
1556  pt[0]=UVector3(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
1557  -fDz*fTthetaSphi-fDy1,-fDz);
1558  pt[1]=UVector3(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
1559  -fDz*fTthetaSphi-fDy1,-fDz);
1560  pt[2]=UVector3(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
1561  -fDz*fTthetaSphi+fDy1,-fDz);
1562  pt[3]=UVector3(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
1563  -fDz*fTthetaSphi+fDy1,-fDz);
1564  pt[4]=UVector3(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
1565  fDz*fTthetaSphi-fDy2,+fDz);
1566  pt[5]=UVector3(fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
1567  fDz*fTthetaSphi-fDy2,+fDz);
1568  pt[6]=UVector3(fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
1569  fDz*fTthetaSphi+fDy2,+fDz);
1570  pt[7]=UVector3(fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
1571  fDz*fTthetaSphi+fDy2,+fDz);
1572 
1573  temp[0] = pt[0].x+(pt[4].x-pt[0].x)
1574  *(aMin.z-pt[0].z)/(pt[4].z-pt[0].z) ;
1575  temp[1] = pt[0].x+(pt[4].x-pt[0].x)
1576  *(aMax.z-pt[0].z)/(pt[4].z-pt[0].z) ;
1577  temp[2] = pt[2].x+(pt[6].x-pt[2].x)
1578  *(aMin.z-pt[2].z)/(pt[6].z-pt[2].z) ;
1579  temp[3] = pt[2].x+(pt[6].x-pt[2].x)
1580  *(aMax.z-pt[2].z)/(pt[6].z-pt[2].z) ;
1581  temp[4] = pt[3].x+(pt[7].x-pt[3].x)
1582  *(aMin.z-pt[3].z)/(pt[7].z-pt[3].z) ;
1583  temp[5] = pt[3].x+(pt[7].x-pt[3].x)
1584  *(aMax.z-pt[3].z)/(pt[7].z-pt[3].z) ;
1585  temp[6] = pt[1].x+(pt[5].x-pt[1].x)
1586  *(aMin.z-pt[1].z)/(pt[5].z-pt[1].z) ;
1587  temp[7] = pt[1].x+(pt[5].x-pt[1].x)
1588  *(aMax.z-pt[1].z)/(pt[5].z-pt[1].z) ;
1589 
1590  aMax.x = - std::fabs(fDz*fTthetaCphi) - fDx1 - fDx2 -fDx3 - fDx4 ;
1591  aMin.x = -aMax.x ;
1592 
1593  for(int i = 0 ; i < 8 ; i++ )
1594  {
1595  if( temp[i] > aMax.x) aMax.x = temp[i] ;
1596  if( temp[i] < aMin.x) aMin.x = temp[i] ;
1597  }
1598  //Y axis
1599  temp[0] = pt[0].y+(pt[4].y-pt[0].y)*(aMin.z-pt[0].z)
1600  /(pt[4].z-pt[0].z) ;
1601  temp[1] = pt[0].y+(pt[4].y-pt[0].y)*(aMax.z-pt[0].z)
1602  /(pt[4].z-pt[0].z) ;
1603  temp[2] = pt[2].y+(pt[6].y-pt[2].y)*(aMin.z-pt[2].z)
1604  /(pt[6].z-pt[2].z) ;
1605  temp[3] = pt[2].y+(pt[6].y-pt[2].y)*(aMax.z-pt[2].z)
1606  /(pt[6].z-pt[2].z) ;
1607 
1608  aMax.y = - std::fabs(fDz*fTthetaSphi) - fDy1 - fDy2 ;
1609  aMin.y = -aMax.y ;
1610 
1611  for( int i = 0 ; i < 4 ; i++ )
1612  {
1613  if( temp[i] > aMax.y ) aMax.y = temp[i] ;
1614  if( temp[i] < aMin.y ) aMin.y = temp[i] ;
1615  }
1616 }
UTrap::MakePlane
bool MakePlane(const UVector3 &p1, const UVector3 &p2, const UVector3 &p3, const UVector3 &p4, UTrapSidePlane &plane)
Definition: UTrap.cc:702
UTrap::fDx1
double fDx1
Definition: UTrap.hh:246
ks1
Definition: UTrap.cc:35
ks3
Definition: UTrap.cc:35
UTrap::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const
Definition: UTrap.cc:1402
UTrapSidePlane::d
double d
Definition: UTrap.hh:76
UTrap::fTthetaSphi
double fTthetaSphi
Definition: UTrap.hh:245
UTrap::GetPointOnPlane
UVector3 GetPointOnPlane(UVector3 p0, UVector3 p1, UVector3 p2, UVector3 p3, double &area) const
Definition: UTrap.cc:1442
UTrap::SetAllParameters
void SetAllParameters(double pDz, double pTheta, double pPhi, double pDy1, double pDx1, double pDx2, double pAlp1, double pDy2, double pDx3, double pDx4, double pAlp2)
Definition: UTrap.cc:491
UTrap::Normal
bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const
Definition: UTrap.cc:820
Eside
Eside
Definition: UTrap.cc:35
UTrap.hh
UTrap::SafetyFromOutside
double SafetyFromOutside(const UVector3 &p, bool precise=false) const
Definition: UTrap.cc:1343
FatalErrorInArguments
Definition: UUtils.hh:35
z
G4double z
Definition: TRTMaterials.hh:39
VUSolid::GetName
const std::string & GetName() const
Definition: VUSolid.hh:103
UVector3::Cross
UVector3 Cross(const UVector3 &) const
Definition: UVector3.hh:262
ks0
Definition: UTrap.cc:35
UTrap::ApproxSurfaceNormal
UVector3 ApproxSurfaceNormal(const UVector3 &p) const
Definition: UTrap.cc:901
UTrap
Definition: UTrap.hh:80
Warning
Definition: UUtils.hh:35
UTrap::fDx2
double fDx2
Definition: UTrap.hh:246
a
G4double a
Definition: TRTMaterials.hh:39
UTrap::fTalpha2
double fTalpha2
Definition: UTrap.hh:247
UTrap::fDz
double fDz
Definition: UTrap.hh:245
VUSolid::Tolerance
static double Tolerance()
Definition: VUSolid.hh:127
UVector3::x
double x
Definition: UVector3.hh:136
smax
const G4int smax
Definition: G4UrQMD1_3Interface.hh:71
UUtils.hh
kUndef
Definition: UTrap.cc:35
CLHEP::normal
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
UTrap::operator=
UTrap & operator=(const UTrap &rhs)
Definition: UTrap.cc:448
UTrap::~UTrap
virtual ~UTrap()
Definition: UTrap.cc:419
UUtils::kInfinity
static const double kInfinity
Definition: UUtils.hh:53
UTrap::fSurfaceArea
double fSurfaceArea
Definition: UTrap.hh:251
kCoplanar_Tolerance
const double kCoplanar_Tolerance
Definition: UTrap.cc:29
UTrap::GetEntityType
UGeometryType GetEntityType() const
Definition: UTrap.cc:1384
UTrap::fDx4
double fDx4
Definition: UTrap.hh:247
UTrap::fTthetaCphi
double fTthetaCphi
Definition: UTrap.hh:245
VUSolid::EnumInside
EnumInside
Definition: VUSolid.hh:23
UVector3::Dot
double Dot(const UVector3 &) const
Definition: UVector3.hh:257
UTrap::fPlanes
UTrapSidePlane fPlanes[4]
Definition: UTrap.hh:248
UTrap::UTrap
UTrap(const std::string &pName, double pDz, double pTheta, double pPhi, double pDy1, double pDx1, double pDx2, double pAlp1, double pDy2, double pDx3, double pDx4, double pAlp2)
Definition: UTrap.cc:42
UVector3::Mag
double Mag() const
Definition: UVector3.cc:48
UTrap::Extent
virtual void Extent(UVector3 &aMin, UVector3 &aMax) const
Definition: UTrap.cc:1546
UTrap::SafetyFromInside
double SafetyFromInside(const UVector3 &p, bool precise=false) const
Definition: UTrap.cc:1064
UTrap::SetPlanes
void SetPlanes(const UVector3 pt[8])
Definition: UTrap.cc:533
UTrapSidePlane::c
double c
Definition: UTrap.hh:76
UTrap::GetPointOnSurface
UVector3 GetPointOnSurface() const
Definition: UTrap.cc:1489
UTrap::fDy1
double fDy1
Definition: UTrap.hh:246
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
UTrap::fDy2
double fDy2
Definition: UTrap.hh:247
UTrap::Clone
VUSolid * Clone() const
Definition: UTrap.cc:1393
UUtils::kPi
static const double kPi
Definition: UUtils.hh:48
UTrap::Inside
VUSolid::EnumInside Inside(const UVector3 &p) const
Definition: UTrap.cc:779
UVector3::Unit
UVector3 Unit() const
Definition: UVector3.cc:80
kPZ
Definition: UTrap.cc:35
UTrapSidePlane::b
double b
Definition: UTrap.hh:76
UTrap::fCubicVolume
double fCubicVolume
Definition: UTrap.hh:250
UGeometryType
std::string UGeometryType
Definition: UTypes.hh:39
UTrap::DistanceToIn
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
Definition: UTrap.cc:942
kMZ
Definition: UTrap.cc:35
UTrap::DistanceToOut
double DistanceToOut(const UVector3 &p, const UVector3 &v, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
Definition: UTrap.cc:1084
FatalError
Definition: UUtils.hh:35
UUtils::Exception
void Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
UVector3::z
double z
Definition: UVector3.hh:138
UUtils::Random
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
UTrap::fTalpha1
double fTalpha1
Definition: UTrap.hh:246
UTrapSidePlane::a
double a
Definition: UTrap.hh:76
UTrap::MakePlanes
bool MakePlanes()
Definition: UTrap.cc:623
UVector3::y
double y
Definition: UVector3.hh:137
ks2
Definition: UTrap.cc:35
UTrap::fDx3
double fDx3
Definition: UTrap.hh:247