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