Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Trap.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4Trap.cc 69788 2013-05-15 12:06:57Z gcosmo $
28 //
29 // class G4Trap
30 //
31 // Implementation for G4Trap class
32 //
33 // History:
34 //
35 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal
36 // 26.04.05 V.Grichine: new SurfaceNormal is default
37 // 19.04.05 V.Grichine: bug fixed in G4Trap("name",G4ThreeVector[8] vp)
38 // 12.12.04 V.Grichine: SurfaceNormal with edges/vertices
39 // 15.11.04 V.Grichine: bug fixed in G4Trap("name",G4ThreeVector[8] vp)
40 // 13.12.99 V.Grichine: bug fixed in DistanceToIn(p,v)
41 // 19.11.99 V.Grichine: kUndef was added to Eside enum
42 // 04.06.99 S.Giani: Fixed CalculateExtent in rotated case.
43 // 08.12.97 J.Allison: Added "nominal" constructor and method SetAllParameters.
44 // 01.11.96 V.Grichine: Costructor for Right Angular Wedge from STEP, G4Trd/Para
45 // 09.09.96 V.Grichine: Final modifications before to commit
46 // 21.03.95 P.Kent: Modified for `tolerant' geometry
47 //
49 
50 #include "G4Trap.hh"
51 #include "globals.hh"
52 
53 #include "G4VoxelLimits.hh"
54 #include "G4AffineTransform.hh"
55 
56 #include "G4VPVParameterisation.hh"
57 
58 #include "Randomize.hh"
59 
60 #include "G4VGraphicsScene.hh"
61 #include "G4Polyhedron.hh"
62 #include "G4NURBS.hh"
63 #include "G4NURBSbox.hh"
64 
65 using namespace CLHEP;
66 
68 //
69 // Accuracy of coplanarity
70 
71 const G4double kCoplanar_Tolerance = 1E-4 ;
72 
74 //
75 // Private enum: Not for external use
76 
77 enum Eside {kUndef,ks0,ks1,ks2,ks3,kPZ,kMZ};
78 
80 //
81 // Constructor - check and set half-widths as well as angles:
82 // final check of coplanarity
83 
84 G4Trap::G4Trap( const G4String& pName,
85  G4double pDz,
86  G4double pTheta, G4double pPhi,
87  G4double pDy1, G4double pDx1, G4double pDx2,
88  G4double pAlp1,
89  G4double pDy2, G4double pDx3, G4double pDx4,
90  G4double pAlp2)
91  : G4CSGSolid(pName)
92 {
93  if ( pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 ||
94  pDx2 <= 0 || pDy2 <= 0 || pDx3 <= 0 || pDx4 <= 0 )
95  {
96  std::ostringstream message;
97  message << "Invalid length parameters for Solid: " << GetName() << G4endl
98  << " X - "
99  << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
100  << " Y - " << pDy1 << ", " << pDy2 << G4endl
101  << " Z - " << pDz;
102  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
103  FatalException, message);
104  }
105 
106  fDz=pDz;
107  fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
108  fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
109 
110  fDy1=pDy1;
111  fDx1=pDx1;
112  fDx2=pDx2;
113  fTalpha1=std::tan(pAlp1);
114 
115  fDy2=pDy2;
116  fDx3=pDx3;
117  fDx4=pDx4;
118  fTalpha2=std::tan(pAlp2);
119 
120  MakePlanes();
121 }
122 
124 //
125 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
126 // which are its vertices. Checking of planarity with preparation of
127 // fPlanes[] and than calculation of other members
128 
129 G4Trap::G4Trap( const G4String& pName,
130  const G4ThreeVector pt[8] )
131  : G4CSGSolid(pName)
132 {
133  G4bool good;
134 
135  // Start with check of centering - the center of gravity trap line
136  // should cross the origin of frame
137 
138  if (!( pt[0].z() < 0
139  && pt[0].z() == pt[1].z() && pt[0].z() == pt[2].z()
140  && pt[0].z() == pt[3].z()
141  && pt[4].z() > 0
142  && pt[4].z() == pt[5].z() && pt[4].z() == pt[6].z()
143  && pt[4].z() == pt[7].z()
144  && std::fabs( pt[0].z() + pt[4].z() ) < kCarTolerance
145  && pt[0].y() == pt[1].y() && pt[2].y() == pt[3].y()
146  && pt[4].y() == pt[5].y() && pt[6].y() == pt[7].y()
147  && std::fabs( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) < kCarTolerance
148  && std::fabs( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() +
149  pt[2].x() + pt[3].x() + pt[6].x() + pt[7].x() ) < kCarTolerance ) )
150  {
151  std::ostringstream message;
152  message << "Invalid vertice coordinates for Solid: " << GetName();
153  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
154  FatalException, message);
155  }
156 
157  // Bottom side with normal approx. -Y
158 
159  good = MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
160 
161  if (!good)
162  {
163  DumpInfo();
164  G4Exception("G4Trap::G4Trap()", "GeomSolids0002", FatalException,
165  "Face at ~-Y not planar.");
166  }
167 
168  // Top side with normal approx. +Y
169 
170  good = MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
171 
172  if (!good)
173  {
174  std::ostringstream message;
175  message << "Face at ~+Y not planar for Solid: " << GetName();
176  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
177  FatalException, message);
178  }
179 
180  // Front side with normal approx. -X
181 
182  good = MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
183 
184  if (!good)
185  {
186  std::ostringstream message;
187  message << "Face at ~-X not planar for Solid: " << GetName();
188  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
189  FatalException, message);
190  }
191 
192  // Back side iwth normal approx. +X
193 
194  good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
195  if (!good)
196  {
197  std::ostringstream message;
198  message << "Face at ~+X not planar for Solid: " << GetName();
199  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
200  FatalException, message);
201  }
202  fDz = (pt[7]).z() ;
203 
204  fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5;
205  fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5;
206  fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5;
207  fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1;
208 
209  fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5;
210  fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5;
211  fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5;
212  fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2;
213 
214  fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz;
215  fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;
216 }
217 
219 //
220 // Constructor for Right Angular Wedge from STEP
221 
222 G4Trap::G4Trap( const G4String& pName,
223  G4double pZ,
224  G4double pY,
225  G4double pX, G4double pLTX )
226  : G4CSGSolid(pName)
227 {
228  G4bool good;
229 
230  if ( pZ<=0 || pY<=0 || pX<=0 || pLTX<=0 || pLTX>pX )
231  {
232  std::ostringstream message;
233  message << "Invalid length parameters for Solid: " << GetName();
234  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
235  FatalException, message);
236  }
237 
238  fDz = 0.5*pZ ;
239  fTthetaCphi = 0 ;
240  fTthetaSphi = 0 ;
241 
242  fDy1 = 0.5*pY;
243  fDx1 = 0.5*pX ;
244  fDx2 = 0.5*pLTX;
245  fTalpha1 = 0.5*(pLTX - pX)/pY;
246 
247  fDy2 = fDy1 ;
248  fDx3 = fDx1;
249  fDx4 = fDx2 ;
250  fTalpha2 = fTalpha1 ;
251 
252  G4ThreeVector pt[8] ;
253 
254  pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
255  -fDz*fTthetaSphi-fDy1,-fDz);
256  pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
257  -fDz*fTthetaSphi-fDy1,-fDz);
258  pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
259  -fDz*fTthetaSphi+fDy1,-fDz);
260  pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
261  -fDz*fTthetaSphi+fDy1,-fDz);
262  pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
263  +fDz*fTthetaSphi-fDy2,+fDz);
264  pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
265  +fDz*fTthetaSphi-fDy2,+fDz);
266  pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
267  +fDz*fTthetaSphi+fDy2,+fDz);
268  pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
269  +fDz*fTthetaSphi+fDy2,+fDz);
270 
271  // Bottom side with normal approx. -Y
272  //
273  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
274  if (!good)
275  {
276  std::ostringstream message;
277  message << "Face at ~-Y not planar for Solid: " << GetName();
278  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
279  FatalException, message);
280  }
281 
282  // Top side with normal approx. +Y
283  //
284  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
285  if (!good)
286  {
287  std::ostringstream message;
288  message << "Face at ~+Y not planar for Solid: " << GetName();
289  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
290  FatalException, message);
291  }
292 
293  // Front side with normal approx. -X
294  //
295  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
296  if (!good)
297  {
298  std::ostringstream message;
299  message << "Face at ~-X not planar for Solid: " << GetName();
300  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
301  FatalException, message);
302  }
303 
304  // Back side iwth normal approx. +X
305  //
306  good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
307  if (!good)
308  {
309  std::ostringstream message;
310  message << "Face at ~+X not planar for Solid: " << GetName();
311  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
312  FatalException, message);
313  }
314 }
315 
317 //
318 // Constructor for G4Trd
319 
320 G4Trap::G4Trap( const G4String& pName,
321  G4double pDx1, G4double pDx2,
322  G4double pDy1, G4double pDy2,
323  G4double pDz )
324  : G4CSGSolid(pName)
325 {
326  G4bool good;
327 
328  if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 )
329  {
330  std::ostringstream message;
331  message << "Invalid length parameters for Solid: " << GetName();
332  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
333  FatalException, message);
334  }
335 
336  fDz = pDz;
337  fTthetaCphi = 0 ;
338  fTthetaSphi = 0 ;
339 
340  fDy1 = pDy1 ;
341  fDx1 = pDx1 ;
342  fDx2 = pDx1 ;
343  fTalpha1 = 0 ;
344 
345  fDy2 = pDy2 ;
346  fDx3 = pDx2 ;
347  fDx4 = pDx2 ;
348  fTalpha2 = 0 ;
349 
350  G4ThreeVector pt[8] ;
351 
352  pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
353  -fDz*fTthetaSphi-fDy1,-fDz);
354  pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
355  -fDz*fTthetaSphi-fDy1,-fDz);
356  pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
357  -fDz*fTthetaSphi+fDy1,-fDz);
358  pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
359  -fDz*fTthetaSphi+fDy1,-fDz);
360  pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
361  +fDz*fTthetaSphi-fDy2,+fDz);
362  pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
363  +fDz*fTthetaSphi-fDy2,+fDz);
364  pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
365  +fDz*fTthetaSphi+fDy2,+fDz);
366  pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
367  +fDz*fTthetaSphi+fDy2,+fDz);
368 
369  // Bottom side with normal approx. -Y
370  //
371  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
372  if (!good)
373  {
374  std::ostringstream message;
375  message << "Face at ~-Y not planar for Solid: " << GetName();
376  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
377  FatalException, message);
378  }
379 
380  // Top side with normal approx. +Y
381  //
382  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
383  if (!good)
384  {
385  std::ostringstream message;
386  message << "Face at ~+Y not planar for Solid: " << GetName();
387  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
388  FatalException, message);
389  }
390 
391  // Front side with normal approx. -X
392  //
393  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
394  if (!good)
395  {
396  std::ostringstream message;
397  message << "Face at ~-X not planar for Solid: " << GetName();
398  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
399  FatalException, message);
400  }
401 
402  // Back side iwth normal approx. +X
403  //
404  good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
405  if (!good)
406  {
407  std::ostringstream message;
408  message << "Face at ~+X not planar for Solid: " << GetName();
409  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
410  FatalException, message);
411  }
412 }
413 
415 //
416 // Constructor for G4Para
417 
418 G4Trap::G4Trap( const G4String& pName,
419  G4double pDx, G4double pDy,
420  G4double pDz,
421  G4double pAlpha,
422  G4double pTheta, G4double pPhi )
423  : G4CSGSolid(pName)
424 {
425  G4bool good;
426 
427  if ( pDz<=0 || pDy<=0 || pDx<=0 )
428  {
429  std::ostringstream message;
430  message << "Invalid length parameters for Solid: " << GetName();
431  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
432  FatalException, message);
433  }
434 
435  fDz = pDz ;
436  fTthetaCphi = std::tan(pTheta)*std::cos(pPhi) ;
437  fTthetaSphi = std::tan(pTheta)*std::sin(pPhi) ;
438 
439  fDy1 = pDy ;
440  fDx1 = pDx ;
441  fDx2 = pDx ;
442  fTalpha1 = std::tan(pAlpha) ;
443 
444  fDy2 = pDy ;
445  fDx3 = pDx ;
446  fDx4 = pDx ;
447  fTalpha2 = fTalpha1 ;
448 
449  G4ThreeVector pt[8] ;
450 
451  pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
452  -fDz*fTthetaSphi-fDy1,-fDz);
453  pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
454  -fDz*fTthetaSphi-fDy1,-fDz);
455  pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
456  -fDz*fTthetaSphi+fDy1,-fDz);
457  pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
458  -fDz*fTthetaSphi+fDy1,-fDz);
459  pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
460  +fDz*fTthetaSphi-fDy2,+fDz);
461  pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
462  +fDz*fTthetaSphi-fDy2,+fDz);
463  pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
464  +fDz*fTthetaSphi+fDy2,+fDz);
465  pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
466  +fDz*fTthetaSphi+fDy2,+fDz);
467 
468  // Bottom side with normal approx. -Y
469  //
470  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
471  if (!good)
472  {
473  std::ostringstream message;
474  message << "Face at ~-Y not planar for Solid: " << GetName();
475  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
476  FatalException, message);
477  }
478 
479  // Top side with normal approx. +Y
480  //
481  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
482  if (!good)
483  {
484  std::ostringstream message;
485  message << "Face at ~+Y not planar for Solid: " << GetName();
486  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
487  FatalException, message);
488  }
489 
490  // Front side with normal approx. -X
491  //
492  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
493  if (!good)
494  {
495  std::ostringstream message;
496  message << "Face at ~-X not planar for Solid: " << GetName();
497  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
498  FatalException, message);
499  }
500 
501  // Back side iwth normal approx. +X
502  //
503  good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
504  if (!good)
505  {
506  std::ostringstream message;
507  message << "Face at ~+X not planar for Solid: " << GetName();
508  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
509  FatalException, message);
510  }
511 }
512 
514 //
515 // Nominal constructor for G4Trap whose parameters are to be set by
516 // a G4VParamaterisation later. Check and set half-widths as well as
517 // angles: final check of coplanarity
518 
519 G4Trap::G4Trap( const G4String& pName )
520  : G4CSGSolid (pName), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
521  fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
522  fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
523 {
524  MakePlanes();
525 }
526 
528 //
529 // Fake default constructor - sets only member data and allocates memory
530 // for usage restricted to object persistency.
531 //
532 G4Trap::G4Trap( __void__& a )
533  : G4CSGSolid(a), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
534  fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
535  fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
536 {
537  MakePlanes();
538 }
539 
541 //
542 // Destructor
543 
544 G4Trap::~G4Trap()
545 {
546 }
547 
549 //
550 // Copy constructor
551 
552 G4Trap::G4Trap(const G4Trap& rhs)
553  : G4CSGSolid(rhs), fDz(rhs.fDz),
554  fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
555  fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
556  fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
557 {
558  for (size_t i=0; i<4; ++i)
559  {
560  fPlanes[i].a = rhs.fPlanes[i].a;
561  fPlanes[i].b = rhs.fPlanes[i].b;
562  fPlanes[i].c = rhs.fPlanes[i].c;
563  fPlanes[i].d = rhs.fPlanes[i].d;
564  }
565 }
566 
568 //
569 // Assignment operator
570 
571 G4Trap& G4Trap::operator = (const G4Trap& rhs)
572 {
573  // Check assignment to self
574  //
575  if (this == &rhs) { return *this; }
576 
577  // Copy base class data
578  //
579  G4CSGSolid::operator=(rhs);
580 
581  // Copy data
582  //
583  fDz = rhs.fDz;
584  fTthetaCphi = rhs.fTthetaCphi; fTthetaSphi = rhs.fTthetaSphi;
585  fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
586  fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
587  for (size_t i=0; i<4; ++i)
588  {
589  fPlanes[i].a = rhs.fPlanes[i].a;
590  fPlanes[i].b = rhs.fPlanes[i].b;
591  fPlanes[i].c = rhs.fPlanes[i].c;
592  fPlanes[i].d = rhs.fPlanes[i].d;
593  }
594 
595  return *this;
596 }
597 
599 //
600 // Set all parameters, as for constructor - check and set half-widths
601 // as well as angles: final check of coplanarity
602 
603 void G4Trap::SetAllParameters ( G4double pDz,
604  G4double pTheta,
605  G4double pPhi,
606  G4double pDy1,
607  G4double pDx1,
608  G4double pDx2,
609  G4double pAlp1,
610  G4double pDy2,
611  G4double pDx3,
612  G4double pDx4,
613  G4double pAlp2 )
614 {
615  if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 || pDx3<=0 || pDx4<=0 )
616  {
617  std::ostringstream message;
618  message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
619  << " X - "
620  << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
621  << " Y - " << pDy1 << ", " << pDy2 << G4endl
622  << " Z - " << pDz;
623  G4Exception("G4Trap::SetAllParameters()", "GeomSolids0002",
624  FatalException, message);
625  }
626  fCubicVolume= 0.;
627  fSurfaceArea= 0.;
628  fpPolyhedron = 0;
629  fDz=pDz;
630  fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
631  fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
632 
633  fDy1=pDy1;
634  fDx1=pDx1;
635  fDx2=pDx2;
636  fTalpha1=std::tan(pAlp1);
637 
638  fDy2=pDy2;
639  fDx3=pDx3;
640  fDx4=pDx4;
641  fTalpha2=std::tan(pAlp2);
642 
643  MakePlanes();
644 }
645 
647 //
648 // Checking of coplanarity
649 
650 G4bool G4Trap::MakePlanes()
651 {
652  G4bool good = true;
653 
654  G4ThreeVector pt[8] ;
655 
656  pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
657  -fDz*fTthetaSphi-fDy1,-fDz);
658  pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
659  -fDz*fTthetaSphi-fDy1,-fDz);
660  pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
661  -fDz*fTthetaSphi+fDy1,-fDz);
662  pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
663  -fDz*fTthetaSphi+fDy1,-fDz);
664  pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
665  +fDz*fTthetaSphi-fDy2,+fDz);
666  pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
667  +fDz*fTthetaSphi-fDy2,+fDz);
668  pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
669  +fDz*fTthetaSphi+fDy2,+fDz);
670  pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
671  +fDz*fTthetaSphi+fDy2,+fDz);
672 
673  // Bottom side with normal approx. -Y
674  //
675  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]) ;
676  if (!good)
677  {
678  std::ostringstream message;
679  message << "Face at ~-Y not planar for Solid: " << GetName();
680  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
681  FatalException, message);
682  }
683 
684  // Top side with normal approx. +Y
685  //
686  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
687  if (!good)
688  {
689  std::ostringstream message;
690  message << "Face at ~+Y not planar for Solid: " << GetName();
691  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
692  FatalException, message);
693  }
694 
695  // Front side with normal approx. -X
696  //
697  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
698  if (!good)
699  {
700  std::ostringstream message;
701  message << "Face at ~-X not planar for Solid: " << GetName();
702  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
703  FatalException, message);
704  }
705 
706  // Back side iwth normal approx. +X
707  //
708  good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
709  if ( !good )
710  {
711  std::ostringstream message;
712  message << "Face at ~+X not planar for Solid: " << GetName();
713  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
714  FatalException, message);
715  }
716 
717  return good;
718 }
719 
721 //
722 // Calculate the coef's of the plane p1->p2->p3->p4->p1
723 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed from
724 // infront of the plane (i.e. from normal direction).
725 //
726 // Return true if the ThreeVectors are coplanar + set coef;s
727 // false if ThreeVectors are not coplanar
728 
729 G4bool G4Trap::MakePlane( const G4ThreeVector& p1,
730  const G4ThreeVector& p2,
731  const G4ThreeVector& p3,
732  const G4ThreeVector& p4,
733  TrapSidePlane& plane )
734 {
735  G4double a, b, c, sd;
736  G4ThreeVector v12, v13, v14, Vcross;
737 
738  G4bool good;
739 
740  v12 = p2 - p1;
741  v13 = p3 - p1;
742  v14 = p4 - p1;
743  Vcross = v12.cross(v13);
744 
745  if (std::fabs(Vcross.dot(v14)/(Vcross.mag()*v14.mag())) > kCoplanar_Tolerance)
746  {
747  good = false;
748  }
749  else
750  {
751  // a,b,c correspond to the x/y/z components of the
752  // normal vector to the plane
753 
754  // a = (p2.y()-p1.y())*(p1.z()+p2.z())+(p3.y()-p2.y())*(p2.z()+p3.z());
755  // a += (p4.y()-p3.y())*(p3.z()+p4.z())+(p1.y()-p4.y())*(p4.z()+p1.z()); // ?
756  // b = (p2.z()-p1.z())*(p1.x()+p2.x())+(p3.z()-p2.z())*(p2.x()+p3.x());
757  // b += (p4.z()-p3.z())*(p3.x()+p4.x())+(p1.z()-p4.z())*(p4.x()+p1.x()); // ?
758  // c = (p2.x()-p1.x())*(p1.y()+p2.y())+(p3.x()-p2.x())*(p2.y()+p3.y());
759  // c += (p4.x()-p3.x())*(p3.y()+p4.y())+(p1.x()-p4.x())*(p4.y()+p1.y()); // ?
760 
761  // Let create diagonals 4-2 and 3-1 than (4-2)x(3-1) provides
762  // vector perpendicular to the plane directed to outside !!!
763  // and a,b,c, = f(1,2,3,4) external relative to trap normal
764 
765  a = +(p4.y() - p2.y())*(p3.z() - p1.z())
766  - (p3.y() - p1.y())*(p4.z() - p2.z());
767 
768  b = -(p4.x() - p2.x())*(p3.z() - p1.z())
769  + (p3.x() - p1.x())*(p4.z() - p2.z());
770 
771  c = +(p4.x() - p2.x())*(p3.y() - p1.y())
772  - (p3.x() - p1.x())*(p4.y() - p2.y());
773 
774  sd = std::sqrt( a*a + b*b + c*c ); // so now vector plane.(a,b,c) is unit
775 
776  if( sd > 0 )
777  {
778  plane.a = a/sd;
779  plane.b = b/sd;
780  plane.c = c/sd;
781  }
782  else
783  {
784  std::ostringstream message;
785  message << "Invalid parameters: norm.mod() <= 0, for Solid: "
786  << GetName();
787  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
788  FatalException, message) ;
789  }
790  // Calculate D: p1 in in plane so D=-n.p1.Vect()
791 
792  plane.d = -( plane.a*p1.x() + plane.b*p1.y() + plane.c*p1.z() );
793 
794  good = true;
795  }
796  return good;
797 }
798 
800 //
801 // Dispatch to parameterisation for replication mechanism dimension
802 // computation & modification.
803 
804 void G4Trap::ComputeDimensions( G4VPVParameterisation* p,
805  const G4int n,
806  const G4VPhysicalVolume* pRep )
807 {
808  p->ComputeDimensions(*this,n,pRep);
809 }
810 
811 
813 //
814 // Calculate extent under transform and specified limit
815 
816 G4bool G4Trap::CalculateExtent( const EAxis pAxis,
817  const G4VoxelLimits& pVoxelLimit,
818  const G4AffineTransform& pTransform,
819  G4double& pMin, G4double& pMax) const
820 {
821  G4double xMin, xMax, yMin, yMax, zMin, zMax;
822  G4bool flag;
823 
824  if (!pTransform.IsRotated())
825  {
826  // Special case handling for unrotated trapezoids
827  // Compute z/x/y/ mins and maxs respecting limits, with early returns
828  // if outside limits. Then switch() on pAxis
829 
830  G4int i ;
831  G4double xoffset;
832  G4double yoffset;
833  G4double zoffset;
834  G4double temp[8] ; // some points for intersection with zMin/zMax
835  G4ThreeVector pt[8]; // vertices after translation
836 
837  xoffset=pTransform.NetTranslation().x();
838  yoffset=pTransform.NetTranslation().y();
839  zoffset=pTransform.NetTranslation().z();
840 
841  pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
842  yoffset-fDz*fTthetaSphi-fDy1,zoffset-fDz);
843  pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
844  yoffset-fDz*fTthetaSphi-fDy1,zoffset-fDz);
845  pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
846  yoffset-fDz*fTthetaSphi+fDy1,zoffset-fDz);
847  pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
848  yoffset-fDz*fTthetaSphi+fDy1,zoffset-fDz);
849  pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
850  yoffset+fDz*fTthetaSphi-fDy2,zoffset+fDz);
851  pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
852  yoffset+fDz*fTthetaSphi-fDy2,zoffset+fDz);
853  pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
854  yoffset+fDz*fTthetaSphi+fDy2,zoffset+fDz);
855  pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
856  yoffset+fDz*fTthetaSphi+fDy2,zoffset+fDz);
857  zMin=zoffset-fDz;
858  zMax=zoffset+fDz;
859 
860  if ( pVoxelLimit.IsZLimited() )
861  {
862  if ( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance)
863  || (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
864  {
865  return false;
866  }
867  else
868  {
869  if ( zMin < pVoxelLimit.GetMinZExtent() )
870  {
871  zMin = pVoxelLimit.GetMinZExtent() ;
872  }
873  if ( zMax > pVoxelLimit.GetMaxZExtent() )
874  {
875  zMax = pVoxelLimit.GetMaxZExtent() ;
876  }
877  }
878  }
879  temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMin-pt[0].z())
880  /(pt[4].z()-pt[0].z()) ;
881  temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMax-pt[0].z())
882  /(pt[4].z()-pt[0].z()) ;
883  temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMin-pt[2].z())
884  /(pt[6].z()-pt[2].z()) ;
885  temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMax-pt[2].z())
886  /(pt[6].z()-pt[2].z()) ;
887 
888  yMax = yoffset - std::fabs(fDz*fTthetaSphi) - fDy1 - fDy2 ;
889  yMin = -yMax ;
890 
891  for( i = 0 ; i < 4 ; i++ )
892  {
893  if( temp[i] > yMax ) yMax = temp[i] ;
894  if( temp[i] < yMin ) yMin = temp[i] ;
895  }
896  if ( pVoxelLimit.IsYLimited() )
897  {
898  if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance)
899  || (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
900  {
901  return false;
902  }
903  else
904  {
905  if ( yMin < pVoxelLimit.GetMinYExtent() )
906  {
907  yMin = pVoxelLimit.GetMinYExtent() ;
908  }
909  if ( yMax > pVoxelLimit.GetMaxYExtent() )
910  {
911  yMax = pVoxelLimit.GetMaxYExtent() ;
912  }
913  }
914  }
915  temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())
916  *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
917  temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())
918  *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
919  temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())
920  *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
921  temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())
922  *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
923  temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())
924  *(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
925  temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())
926  *(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
927  temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())
928  *(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
929  temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())
930  *(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
931 
932  xMax = xoffset - std::fabs(fDz*fTthetaCphi) - fDx1 - fDx2 -fDx3 - fDx4 ;
933  xMin = -xMax ;
934 
935  for( i = 0 ; i < 8 ; i++ )
936  {
937  if( temp[i] > xMax) xMax = temp[i] ;
938  if( temp[i] < xMin) xMin = temp[i] ;
939  }
940  if (pVoxelLimit.IsXLimited()) // xMax/Min = f(yMax/Min) ?
941  {
942  if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance)
943  || (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
944  {
945  return false;
946  }
947  else
948  {
949  if ( xMin < pVoxelLimit.GetMinXExtent() )
950  {
951  xMin = pVoxelLimit.GetMinXExtent() ;
952  }
953  if ( xMax > pVoxelLimit.GetMaxXExtent() )
954  {
955  xMax = pVoxelLimit.GetMaxXExtent() ;
956  }
957  }
958  }
959  switch (pAxis)
960  {
961  case kXAxis:
962  pMin=xMin;
963  pMax=xMax;
964  break;
965 
966  case kYAxis:
967  pMin=yMin;
968  pMax=yMax;
969  break;
970 
971  case kZAxis:
972  pMin=zMin;
973  pMax=zMax;
974  break;
975 
976  default:
977  break;
978  }
979  pMin -= kCarTolerance;
980  pMax += kCarTolerance;
981 
982  flag = true;
983  }
984  else // General rotated case -
985  {
986  G4bool existsAfterClip = false ;
987  G4ThreeVectorList* vertices;
988  pMin = +kInfinity;
989  pMax = -kInfinity;
990 
991  // Calculate rotated vertex coordinates. Operator 'new' is called
992 
993  vertices = CreateRotatedVertices(pTransform);
994 
995  xMin = +kInfinity; yMin = +kInfinity; zMin = +kInfinity;
996  xMax = -kInfinity; yMax = -kInfinity; zMax = -kInfinity;
997 
998  for( G4int nv = 0 ; nv < 8 ; nv++ )
999  {
1000  if( (*vertices)[nv].x() > xMax ) xMax = (*vertices)[nv].x();
1001  if( (*vertices)[nv].y() > yMax ) yMax = (*vertices)[nv].y();
1002  if( (*vertices)[nv].z() > zMax ) zMax = (*vertices)[nv].z();
1003 
1004  if( (*vertices)[nv].x() < xMin ) xMin = (*vertices)[nv].x();
1005  if( (*vertices)[nv].y() < yMin ) yMin = (*vertices)[nv].y();
1006  if( (*vertices)[nv].z() < zMin ) zMin = (*vertices)[nv].z();
1007  }
1008  if ( pVoxelLimit.IsZLimited() )
1009  {
1010  if ( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance)
1011  || (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
1012  {
1013  delete vertices ; // 'new' in the function called
1014  return false;
1015  }
1016  else
1017  {
1018  if ( zMin < pVoxelLimit.GetMinZExtent() )
1019  {
1020  zMin = pVoxelLimit.GetMinZExtent() ;
1021  }
1022  if ( zMax > pVoxelLimit.GetMaxZExtent() )
1023  {
1024  zMax = pVoxelLimit.GetMaxZExtent() ;
1025  }
1026  }
1027  }
1028  if ( pVoxelLimit.IsYLimited() )
1029  {
1030  if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance)
1031  || (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
1032  {
1033  delete vertices ; // 'new' in the function called
1034  return false;
1035  }
1036  else
1037  {
1038  if ( yMin < pVoxelLimit.GetMinYExtent() )
1039  {
1040  yMin = pVoxelLimit.GetMinYExtent() ;
1041  }
1042  if ( yMax > pVoxelLimit.GetMaxYExtent() )
1043  {
1044  yMax = pVoxelLimit.GetMaxYExtent() ;
1045  }
1046  }
1047  }
1048  if ( pVoxelLimit.IsXLimited() )
1049  {
1050  if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance)
1051  || (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
1052  {
1053  delete vertices ; // 'new' in the function called
1054  return false ;
1055  }
1056  else
1057  {
1058  if ( xMin < pVoxelLimit.GetMinXExtent() )
1059  {
1060  xMin = pVoxelLimit.GetMinXExtent() ;
1061  }
1062  if ( xMax > pVoxelLimit.GetMaxXExtent() )
1063  {
1064  xMax = pVoxelLimit.GetMaxXExtent() ;
1065  }
1066  }
1067  }
1068  switch (pAxis)
1069  {
1070  case kXAxis:
1071  pMin=xMin;
1072  pMax=xMax;
1073  break;
1074 
1075  case kYAxis:
1076  pMin=yMin;
1077  pMax=yMax;
1078  break;
1079 
1080  case kZAxis:
1081  pMin=zMin;
1082  pMax=zMax;
1083  break;
1084 
1085  default:
1086  break;
1087  }
1088  if ( (pMin != kInfinity) || (pMax != -kInfinity) )
1089  {
1090  existsAfterClip=true;
1091 
1092  // Add tolerance to avoid precision troubles
1093  //
1094  pMin -= kCarTolerance ;
1095  pMax += kCarTolerance ;
1096  }
1097  delete vertices ; // 'new' in the function called
1098  flag = existsAfterClip ;
1099  }
1100  return flag;
1101 }
1102 
1103 
1105 //
1106 // Return whether point inside/outside/on surface, using tolerance
1107 
1108 EInside G4Trap::Inside( const G4ThreeVector& p ) const
1109 {
1110  EInside in;
1111  G4double Dist;
1112  G4int i;
1113  if ( std::fabs(p.z()) <= fDz-kCarTolerance*0.5)
1114  {
1115  in = kInside;
1116 
1117  for ( i = 0;i < 4;i++ )
1118  {
1119  Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
1120  +fPlanes[i].c*p.z() + fPlanes[i].d;
1121 
1122  if (Dist > kCarTolerance*0.5) return in = kOutside;
1123  else if (Dist > -kCarTolerance*0.5) in = kSurface;
1124 
1125  }
1126  }
1127  else if (std::fabs(p.z()) <= fDz+kCarTolerance*0.5)
1128  {
1129  in = kSurface;
1130 
1131  for ( i = 0; i < 4; i++ )
1132  {
1133  Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
1134  +fPlanes[i].c*p.z() + fPlanes[i].d;
1135 
1136  if (Dist > kCarTolerance*0.5) return in = kOutside;
1137  }
1138  }
1139  else in = kOutside;
1140 
1141  return in;
1142 }
1143 
1145 //
1146 // Calculate side nearest to p, and return normal
1147 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
1148 
1149 G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const
1150 {
1151  G4int i, noSurfaces = 0;
1152  G4double dist, distz, distx, disty, distmx, distmy, safe = kInfinity;
1153  G4double delta = 0.5*kCarTolerance;
1154  G4ThreeVector norm, sumnorm(0.,0.,0.);
1155 
1156  for (i = 0; i < 4; i++)
1157  {
1158  dist = std::fabs(fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
1159  + fPlanes[i].c*p.z() + fPlanes[i].d);
1160  if ( dist < safe )
1161  {
1162  safe = dist;
1163  }
1164  }
1165  distz = std::fabs( std::fabs( p.z() ) - fDz );
1166 
1167  distmy = std::fabs( fPlanes[0].a*p.x() + fPlanes[0].b*p.y()
1168  + fPlanes[0].c*p.z() + fPlanes[0].d );
1169 
1170  disty = std::fabs( fPlanes[1].a*p.x() + fPlanes[1].b*p.y()
1171  + fPlanes[1].c*p.z() + fPlanes[1].d );
1172 
1173  distmx = std::fabs( fPlanes[2].a*p.x() + fPlanes[2].b*p.y()
1174  + fPlanes[2].c*p.z() + fPlanes[2].d );
1175 
1176  distx = std::fabs( fPlanes[3].a*p.x() + fPlanes[3].b*p.y()
1177  + fPlanes[3].c*p.z() + fPlanes[3].d );
1178 
1179  G4ThreeVector nX = G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1180  G4ThreeVector nmX = G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1181  G4ThreeVector nY = G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1182  G4ThreeVector nmY = G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1183  G4ThreeVector nZ = G4ThreeVector(0.,0.,1.0);
1184 
1185  if (distx <= delta)
1186  {
1187  noSurfaces ++;
1188  sumnorm += nX;
1189  }
1190  if (distmx <= delta)
1191  {
1192  noSurfaces ++;
1193  sumnorm += nmX;
1194  }
1195  if (disty <= delta)
1196  {
1197  noSurfaces ++;
1198  sumnorm += nY;
1199  }
1200  if (distmy <= delta)
1201  {
1202  noSurfaces ++;
1203  sumnorm += nmY;
1204  }
1205  if (distz <= delta)
1206  {
1207  noSurfaces ++;
1208  if ( p.z() >= 0.) sumnorm += nZ;
1209  else sumnorm -= nZ;
1210  }
1211  if ( noSurfaces == 0 )
1212  {
1213 #ifdef G4CSGDEBUG
1214  G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
1215  JustWarning, "Point p is not on surface !?" );
1216 #endif
1217  norm = ApproxSurfaceNormal(p);
1218  }
1219  else if ( noSurfaces == 1 ) norm = sumnorm;
1220  else norm = sumnorm.unit();
1221  return norm;
1222 }
1223 
1225 //
1226 // Algorithm for SurfaceNormal() following the original specification
1227 // for points not on the surface
1228 
1229 G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
1230 {
1231  G4double safe=kInfinity,Dist,safez;
1232  G4int i,imin=0;
1233  for (i=0;i<4;i++)
1234  {
1235  Dist=std::fabs(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1236  +fPlanes[i].c*p.z()+fPlanes[i].d);
1237  if (Dist<safe)
1238  {
1239  safe=Dist;
1240  imin=i;
1241  }
1242  }
1243  safez=std::fabs(std::fabs(p.z())-fDz);
1244  if (safe<safez)
1245  {
1246  return G4ThreeVector(fPlanes[imin].a,fPlanes[imin].b,fPlanes[imin].c);
1247  }
1248  else
1249  {
1250  if (p.z()>0)
1251  {
1252  return G4ThreeVector(0,0,1);
1253  }
1254  else
1255  {
1256  return G4ThreeVector(0,0,-1);
1257  }
1258  }
1259 }
1260 
1262 //
1263 // Calculate distance to shape from outside - return kInfinity if no intersection
1264 //
1265 // ALGORITHM:
1266 // For each component, calculate pair of minimum and maximum intersection
1267 // values for which the particle is in the extent of the shape
1268 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
1269 
1270 G4double G4Trap::DistanceToIn( const G4ThreeVector& p,
1271  const G4ThreeVector& v ) const
1272 {
1273 
1274  G4double snxt; // snxt = default return value
1275  G4double max,smax,smin;
1276  G4double pdist,Comp,vdist;
1277  G4int i;
1278  //
1279  // Z Intersection range
1280  //
1281  if ( v.z() > 0 )
1282  {
1283  max = fDz - p.z() ;
1284  if (max > 0.5*kCarTolerance)
1285  {
1286  smax = max/v.z();
1287  smin = (-fDz-p.z())/v.z();
1288  }
1289  else
1290  {
1291  return snxt=kInfinity;
1292  }
1293  }
1294  else if (v.z() < 0 )
1295  {
1296  max = - fDz - p.z() ;
1297  if (max < -0.5*kCarTolerance )
1298  {
1299  smax=max/v.z();
1300  smin=(fDz-p.z())/v.z();
1301  }
1302  else
1303  {
1304  return snxt=kInfinity;
1305  }
1306  }
1307  else
1308  {
1309  if (std::fabs(p.z())<fDz - 0.5*kCarTolerance) // Inside was <=fDz
1310  {
1311  smin=0;
1312  smax=kInfinity;
1313  }
1314  else
1315  {
1316  return snxt=kInfinity;
1317  }
1318  }
1319 
1320  for (i=0;i<4;i++)
1321  {
1322  pdist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1323  +fPlanes[i].c*p.z()+fPlanes[i].d;
1324  Comp=fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
1325  if ( pdist >= -0.5*kCarTolerance ) // was >0
1326  {
1327  //
1328  // Outside the plane -> this is an extent entry distance
1329  //
1330  if (Comp >= 0) // was >0
1331  {
1332  return snxt=kInfinity ;
1333  }
1334  else
1335  {
1336  vdist=-pdist/Comp;
1337  if (vdist>smin)
1338  {
1339  if (vdist<smax)
1340  {
1341  smin = vdist;
1342  }
1343  else
1344  {
1345  return snxt=kInfinity;
1346  }
1347  }
1348  }
1349  }
1350  else
1351  {
1352  //
1353  // Inside the plane -> couble be an extent exit distance (smax)
1354  //
1355  if (Comp>0) // Will leave extent
1356  {
1357  vdist=-pdist/Comp;
1358  if (vdist<smax)
1359  {
1360  if (vdist>smin)
1361  {
1362  smax=vdist;
1363  }
1364  else
1365  {
1366  return snxt=kInfinity;
1367  }
1368  }
1369  }
1370  }
1371  }
1372  //
1373  // Checks in non z plane intersections ensure smin<smax
1374  //
1375  if (smin >=0 )
1376  {
1377  snxt = smin ;
1378  }
1379  else
1380  {
1381  snxt = 0 ;
1382  }
1383  return snxt;
1384 }
1385 
1387 //
1388 // Calculate exact shortest distance to any boundary from outside
1389 // This is the best fast estimation of the shortest distance to trap
1390 // - Returns 0 is ThreeVector inside
1391 
1392 G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const
1393 {
1394  G4double safe=0.0,Dist;
1395  G4int i;
1396  safe=std::fabs(p.z())-fDz;
1397  for (i=0;i<4;i++)
1398  {
1399  Dist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1400  +fPlanes[i].c*p.z()+fPlanes[i].d;
1401  if (Dist > safe) safe=Dist;
1402  }
1403  if (safe<0) safe=0;
1404  return safe;
1405 }
1406 
1408 //
1409 // Calculate distance to surface of shape from inside
1410 // Calculate distance to x/y/z planes - smallest is exiting distance
1411 
1412 G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
1413  const G4bool calcNorm,
1414  G4bool *validNorm, G4ThreeVector *n) const
1415 {
1416  Eside side = kUndef;
1417  G4double snxt; // snxt = return value
1418  G4double pdist,Comp,vdist,max;
1419  //
1420  // Z Intersections
1421  //
1422  if (v.z()>0)
1423  {
1424  max=fDz-p.z();
1425  if (max>kCarTolerance/2)
1426  {
1427  snxt=max/v.z();
1428  side=kPZ;
1429  }
1430  else
1431  {
1432  if (calcNorm)
1433  {
1434  *validNorm=true;
1435  *n=G4ThreeVector(0,0,1);
1436  }
1437  return snxt=0;
1438  }
1439  }
1440  else if (v.z()<0)
1441  {
1442  max=-fDz-p.z();
1443  if (max<-kCarTolerance/2)
1444  {
1445  snxt=max/v.z();
1446  side=kMZ;
1447  }
1448  else
1449  {
1450  if (calcNorm)
1451  {
1452  *validNorm=true;
1453  *n=G4ThreeVector(0,0,-1);
1454  }
1455  return snxt=0;
1456  }
1457  }
1458  else
1459  {
1460  snxt=kInfinity;
1461  }
1462 
1463  //
1464  // Intersections with planes[0] (expanded because of setting enum)
1465  //
1466  pdist=fPlanes[0].a*p.x()+fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
1467  Comp=fPlanes[0].a*v.x()+fPlanes[0].b*v.y()+fPlanes[0].c*v.z();
1468  if (pdist>0)
1469  {
1470  // Outside the plane
1471  if (Comp>0)
1472  {
1473  // Leaving immediately
1474  if (calcNorm)
1475  {
1476  *validNorm=true;
1477  *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1478  }
1479  return snxt=0;
1480  }
1481  }
1482  else if (pdist<-kCarTolerance/2)
1483  {
1484  // Inside the plane
1485  if (Comp>0)
1486  {
1487  // Will leave extent
1488  vdist=-pdist/Comp;
1489  if (vdist<snxt)
1490  {
1491  snxt=vdist;
1492  side=ks0;
1493  }
1494  }
1495  }
1496  else
1497  {
1498  // On surface
1499  if (Comp>0)
1500  {
1501  if (calcNorm)
1502  {
1503  *validNorm=true;
1504  *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1505  }
1506  return snxt=0;
1507  }
1508  }
1509 
1510  //
1511  // Intersections with planes[1] (expanded because of setting enum)
1512  //
1513  pdist=fPlanes[1].a*p.x()+fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
1514  Comp=fPlanes[1].a*v.x()+fPlanes[1].b*v.y()+fPlanes[1].c*v.z();
1515  if (pdist>0)
1516  {
1517  // Outside the plane
1518  if (Comp>0)
1519  {
1520  // Leaving immediately
1521  if (calcNorm)
1522  {
1523  *validNorm=true;
1524  *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1525  }
1526  return snxt=0;
1527  }
1528  }
1529  else if (pdist<-kCarTolerance/2)
1530  {
1531  // Inside the plane
1532  if (Comp>0)
1533  {
1534  // Will leave extent
1535  vdist=-pdist/Comp;
1536  if (vdist<snxt)
1537  {
1538  snxt=vdist;
1539  side=ks1;
1540  }
1541  }
1542  }
1543  else
1544  {
1545  // On surface
1546  if (Comp>0)
1547  {
1548  if (calcNorm)
1549  {
1550  *validNorm=true;
1551  *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1552  }
1553  return snxt=0;
1554  }
1555  }
1556 
1557  //
1558  // Intersections with planes[2] (expanded because of setting enum)
1559  //
1560  pdist=fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d;
1561  Comp=fPlanes[2].a*v.x()+fPlanes[2].b*v.y()+fPlanes[2].c*v.z();
1562  if (pdist>0)
1563  {
1564  // Outside the plane
1565  if (Comp>0)
1566  {
1567  // Leaving immediately
1568  if (calcNorm)
1569  {
1570  *validNorm=true;
1571  *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1572  }
1573  return snxt=0;
1574  }
1575  }
1576  else if (pdist<-kCarTolerance/2)
1577  {
1578  // Inside the plane
1579  if (Comp>0)
1580  {
1581  // Will leave extent
1582  vdist=-pdist/Comp;
1583  if (vdist<snxt)
1584  {
1585  snxt=vdist;
1586  side=ks2;
1587  }
1588  }
1589  }
1590  else
1591  {
1592  // On surface
1593  if (Comp>0)
1594  {
1595  if (calcNorm)
1596  {
1597  *validNorm=true;
1598  *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1599  }
1600  return snxt=0;
1601  }
1602  }
1603 
1604  //
1605  // Intersections with planes[3] (expanded because of setting enum)
1606  //
1607  pdist=fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d;
1608  Comp=fPlanes[3].a*v.x()+fPlanes[3].b*v.y()+fPlanes[3].c*v.z();
1609  if (pdist>0)
1610  {
1611  // Outside the plane
1612  if (Comp>0)
1613  {
1614  // Leaving immediately
1615  if (calcNorm)
1616  {
1617  *validNorm=true;
1618  *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1619  }
1620  return snxt=0;
1621  }
1622  }
1623  else if (pdist<-kCarTolerance/2)
1624  {
1625  // Inside the plane
1626  if (Comp>0)
1627  {
1628  // Will leave extent
1629  vdist=-pdist/Comp;
1630  if (vdist<snxt)
1631  {
1632  snxt=vdist;
1633  side=ks3;
1634  }
1635  }
1636  }
1637  else
1638  {
1639  // On surface
1640  if (Comp>0)
1641  {
1642  if (calcNorm)
1643  {
1644  *validNorm=true;
1645  *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1646  }
1647  return snxt=0;
1648  }
1649  }
1650 
1651  // set normal
1652  if (calcNorm)
1653  {
1654  *validNorm=true;
1655  switch(side)
1656  {
1657  case ks0:
1658  *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1659  break;
1660  case ks1:
1661  *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1662  break;
1663  case ks2:
1664  *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1665  break;
1666  case ks3:
1667  *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1668  break;
1669  case kMZ:
1670  *n=G4ThreeVector(0,0,-1);
1671  break;
1672  case kPZ:
1673  *n=G4ThreeVector(0,0,1);
1674  break;
1675  default:
1676  G4cout << G4endl;
1677  DumpInfo();
1678  std::ostringstream message;
1679  G4int oldprc = message.precision(16);
1680  message << "Undefined side for valid surface normal to solid."
1681  << G4endl
1682  << "Position:" << G4endl << G4endl
1683  << "p.x() = " << p.x()/mm << " mm" << G4endl
1684  << "p.y() = " << p.y()/mm << " mm" << G4endl
1685  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1686  << "Direction:" << G4endl << G4endl
1687  << "v.x() = " << v.x() << G4endl
1688  << "v.y() = " << v.y() << G4endl
1689  << "v.z() = " << v.z() << G4endl << G4endl
1690  << "Proposed distance :" << G4endl << G4endl
1691  << "snxt = " << snxt/mm << " mm" << G4endl;
1692  message.precision(oldprc);
1693  G4Exception("G4Trap::DistanceToOut(p,v,..)","GeomSolids1002",
1694  JustWarning, message);
1695  break;
1696  }
1697  }
1698  return snxt;
1699 }
1700 
1702 //
1703 // Calculate exact shortest distance to any boundary from inside
1704 // - Returns 0 is ThreeVector outside
1705 
1706 G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const
1707 {
1708  G4double safe=0.0,Dist;
1709  G4int i;
1710 
1711 #ifdef G4CSGDEBUG
1712  if( Inside(p) == kOutside )
1713  {
1714  G4int oldprc = G4cout.precision(16) ;
1715  G4cout << G4endl ;
1716  DumpInfo();
1717  G4cout << "Position:" << G4endl << G4endl ;
1718  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1719  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1720  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1721  G4cout.precision(oldprc) ;
1722  G4Exception("G4Trap::DistanceToOut(p)",
1723  "GeomSolids1002", JustWarning, "Point p is outside !?" );
1724  }
1725 #endif
1726 
1727  safe=fDz-std::fabs(p.z());
1728  if (safe<0) safe=0;
1729  else
1730  {
1731  for (i=0;i<4;i++)
1732  {
1733  Dist=-(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1734  +fPlanes[i].c*p.z()+fPlanes[i].d);
1735  if (Dist<safe) safe=Dist;
1736  }
1737  if (safe<0) safe=0;
1738  }
1739  return safe;
1740 }
1741 
1743 //
1744 // Create a List containing the transformed vertices
1745 // Ordering [0-3] -fDz cross section
1746 // [4-7] +fDz cross section such that [0] is below [4],
1747 // [1] below [5] etc.
1748 // Note:
1749 // Caller has deletion resposibility
1750 
1751 G4ThreeVectorList*
1752 G4Trap::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
1753 {
1754  G4ThreeVectorList *vertices;
1755  vertices=new G4ThreeVectorList();
1756  if (vertices)
1757  {
1758  vertices->reserve(8);
1759  G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
1760  -fDz*fTthetaSphi-fDy1,-fDz);
1761  G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
1762  -fDz*fTthetaSphi-fDy1,-fDz);
1763  G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
1764  -fDz*fTthetaSphi+fDy1,-fDz);
1765  G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
1766  -fDz*fTthetaSphi+fDy1,-fDz);
1767  G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
1768  +fDz*fTthetaSphi-fDy2,+fDz);
1769  G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
1770  +fDz*fTthetaSphi-fDy2,+fDz);
1771  G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
1772  +fDz*fTthetaSphi+fDy2,+fDz);
1773  G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
1774  +fDz*fTthetaSphi+fDy2,+fDz);
1775 
1776  vertices->push_back(pTransform.TransformPoint(vertex0));
1777  vertices->push_back(pTransform.TransformPoint(vertex1));
1778  vertices->push_back(pTransform.TransformPoint(vertex2));
1779  vertices->push_back(pTransform.TransformPoint(vertex3));
1780  vertices->push_back(pTransform.TransformPoint(vertex4));
1781  vertices->push_back(pTransform.TransformPoint(vertex5));
1782  vertices->push_back(pTransform.TransformPoint(vertex6));
1783  vertices->push_back(pTransform.TransformPoint(vertex7));
1784  }
1785  else
1786  {
1787  DumpInfo();
1788  G4Exception("G4Trap::CreateRotatedVertices()",
1789  "GeomSolids0003", FatalException,
1790  "Error in allocation of vertices. Out of memory !");
1791  }
1792  return vertices;
1793 }
1794 
1796 //
1797 // GetEntityType
1798 
1799 G4GeometryType G4Trap::GetEntityType() const
1800 {
1801  return G4String("G4Trap");
1802 }
1803 
1805 //
1806 // Make a clone of the object
1807 //
1808 G4VSolid* G4Trap::Clone() const
1809 {
1810  return new G4Trap(*this);
1811 }
1812 
1814 //
1815 // Stream object contents to an output stream
1816 
1817 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1818 {
1819  G4int oldprc = os.precision(16);
1820  os << "-----------------------------------------------------------\n"
1821  << " *** Dump for solid - " << GetName() << " ***\n"
1822  << " ===================================================\n"
1823  << " Solid type: G4Trap\n"
1824  << " Parameters: \n"
1825  << " half length Z: " << fDz/mm << " mm \n"
1826  << " half length Y of face -fDz: " << fDy1/mm << " mm \n"
1827  << " half length X of side -fDy1, face -fDz: " << fDx1/mm << " mm \n"
1828  << " half length X of side +fDy1, face -fDz: " << fDx2/mm << " mm \n"
1829  << " half length Y of face +fDz: " << fDy2/mm << " mm \n"
1830  << " half length X of side -fDy2, face +fDz: " << fDx3/mm << " mm \n"
1831  << " half length X of side +fDy2, face +fDz: " << fDx4/mm << " mm \n"
1832  << " std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree << " degrees \n"
1833  << " std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree << " degrees \n"
1834  << " std::tan(alpha), -fDz: " << fTalpha1/degree << " degrees \n"
1835  << " std::tan(alpha), +fDz: " << fTalpha2/degree << " degrees \n"
1836  << " trap side plane equations:\n"
1837  << " " << fPlanes[0].a << " X + " << fPlanes[0].b << " Y + "
1838  << fPlanes[0].c << " Z + " << fPlanes[0].d << " = 0\n"
1839  << " " << fPlanes[1].a << " X + " << fPlanes[1].b << " Y + "
1840  << fPlanes[1].c << " Z + " << fPlanes[1].d << " = 0\n"
1841  << " " << fPlanes[2].a << " X + " << fPlanes[2].b << " Y + "
1842  << fPlanes[2].c << " Z + " << fPlanes[2].d << " = 0\n"
1843  << " " << fPlanes[3].a << " X + " << fPlanes[3].b << " Y + "
1844  << fPlanes[3].c << " Z + " << fPlanes[3].d << " = 0\n"
1845  << "-----------------------------------------------------------\n";
1846  os.precision(oldprc);
1847 
1848  return os;
1849 }
1850 
1852 //
1853 // GetPointOnPlane
1854 //
1855 // Auxiliary method for Get Point on Surface
1856 
1857 G4ThreeVector G4Trap::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1,
1858  G4ThreeVector p2, G4ThreeVector p3,
1859  G4double& area) const
1860 {
1861  G4double lambda1, lambda2, chose, aOne, aTwo;
1862  G4ThreeVector t, u, v, w, Area, normal;
1863 
1864  t = p1 - p0;
1865  u = p2 - p1;
1866  v = p3 - p2;
1867  w = p0 - p3;
1868 
1869  Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
1870  w.z()*v.x() - w.x()*v.z(),
1871  w.x()*v.y() - w.y()*v.x());
1872 
1873  aOne = 0.5*Area.mag();
1874 
1875  Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
1876  t.z()*u.x() - t.x()*u.z(),
1877  t.x()*u.y() - t.y()*u.x());
1878 
1879  aTwo = 0.5*Area.mag();
1880 
1881  area = aOne + aTwo;
1882 
1883  chose = RandFlat::shoot(0.,aOne+aTwo);
1884 
1885  if( (chose>=0.) && (chose < aOne) )
1886  {
1887  lambda1 = RandFlat::shoot(0.,1.);
1888  lambda2 = RandFlat::shoot(0.,lambda1);
1889  return (p2+lambda1*v+lambda2*w);
1890  }
1891 
1892  // else
1893 
1894  lambda1 = RandFlat::shoot(0.,1.);
1895  lambda2 = RandFlat::shoot(0.,lambda1);
1896 
1897  return (p0+lambda1*t+lambda2*u);
1898 }
1899 
1901 //
1902 // GetPointOnSurface
1903 
1904 G4ThreeVector G4Trap::GetPointOnSurface() const
1905 {
1906  G4double aOne, aTwo, aThree, aFour, aFive, aSix, chose;
1907  G4ThreeVector One, Two, Three, Four, Five, Six, test;
1908  G4ThreeVector pt[8];
1909 
1910  pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
1911  -fDz*fTthetaSphi-fDy1,-fDz);
1912  pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
1913  -fDz*fTthetaSphi-fDy1,-fDz);
1914  pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
1915  -fDz*fTthetaSphi+fDy1,-fDz);
1916  pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
1917  -fDz*fTthetaSphi+fDy1,-fDz);
1918  pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
1919  +fDz*fTthetaSphi-fDy2,+fDz);
1920  pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
1921  +fDz*fTthetaSphi-fDy2,+fDz);
1922  pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
1923  +fDz*fTthetaSphi+fDy2,+fDz);
1924  pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
1925  +fDz*fTthetaSphi+fDy2,+fDz);
1926 
1927  // make sure we provide the points in a clockwise fashion
1928 
1929  One = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
1930  Two = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
1931  Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
1932  Four = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour);
1933  Five = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
1934  Six = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
1935 
1936  chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
1937  if( (chose>=0.) && (chose<aOne) )
1938  { return One; }
1939  else if( (chose>=aOne) && (chose<aOne+aTwo) )
1940  { return Two; }
1941  else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) )
1942  { return Three; }
1943  else if( (chose>=aOne+aTwo+aThree) && (chose<aOne+aTwo+aThree+aFour) )
1944  { return Four; }
1945  else if( (chose>=aOne+aTwo+aThree+aFour)
1946  && (chose<aOne+aTwo+aThree+aFour+aFive) )
1947  { return Five; }
1948  return Six;
1949 }
1950 
1952 //
1953 // Methods for visualisation
1954 
1955 void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1956 {
1957  scene.AddSolid (*this);
1958 }
1959 
1960 G4Polyhedron* G4Trap::CreatePolyhedron () const
1961 {
1962  G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1963  G4double alpha1 = std::atan(fTalpha1);
1964  G4double alpha2 = std::atan(fTalpha2);
1965  G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi+fTthetaSphi*fTthetaSphi));
1966 
1967  return new G4PolyhedronTrap(fDz, theta, phi,
1968  fDy1, fDx1, fDx2, alpha1,
1969  fDy2, fDx3, fDx4, alpha2);
1970 }
1971 
1972 G4NURBS* G4Trap::CreateNURBS () const
1973 {
1974  // return new G4NURBSbox (fDx, fDy, fDz);
1975  return 0 ;
1976 }
Generated on Sat May 25 2013 14:33:14 for Geant4 by   doxygen 1.8.4