Geant4  10.00.p02
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 81636 2014-06-04 09:06:08Z 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 
63 using namespace CLHEP;
64 
66 //
67 // Accuracy of coplanarity
68 
69 const G4double kCoplanar_Tolerance = 1E-4 ;
70 
72 //
73 // Private enum: Not for external use
74 
75 enum Eside {kUndef,ks0,ks1,ks2,ks3,kPZ,kMZ};
76 
78 //
79 // Constructor - check and set half-widths as well as angles:
80 // final check of coplanarity
81 
82 G4Trap::G4Trap( const G4String& pName,
83  G4double pDz,
84  G4double pTheta, G4double pPhi,
85  G4double pDy1, G4double pDx1, G4double pDx2,
86  G4double pAlp1,
87  G4double pDy2, G4double pDx3, G4double pDx4,
88  G4double pAlp2)
89  : G4CSGSolid(pName)
90 {
91  if ( pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 ||
92  pDx2 <= 0 || pDy2 <= 0 || pDx3 <= 0 || pDx4 <= 0 )
93  {
94  std::ostringstream message;
95  message << "Invalid length parameters for Solid: " << GetName() << G4endl
96  << " X - "
97  << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
98  << " Y - " << pDy1 << ", " << pDy2 << G4endl
99  << " Z - " << pDz;
100  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
101  FatalException, message);
102  }
103 
104  fDz=pDz;
105  fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
106  fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
107 
108  fDy1=pDy1;
109  fDx1=pDx1;
110  fDx2=pDx2;
111  fTalpha1=std::tan(pAlp1);
112 
113  fDy2=pDy2;
114  fDx3=pDx3;
115  fDx4=pDx4;
116  fTalpha2=std::tan(pAlp2);
117 
118  MakePlanes();
119 }
120 
122 //
123 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
124 // which are its vertices. Checking of planarity with preparation of
125 // fPlanes[] and than calculation of other members
126 
127 G4Trap::G4Trap( const G4String& pName,
128  const G4ThreeVector pt[8] )
129  : G4CSGSolid(pName)
130 {
131  G4bool good;
132 
133  // Start with check of centering - the center of gravity trap line
134  // should cross the origin of frame
135 
136  if (!( pt[0].z() < 0
137  && pt[0].z() == pt[1].z() && pt[0].z() == pt[2].z()
138  && pt[0].z() == pt[3].z()
139  && pt[4].z() > 0
140  && pt[4].z() == pt[5].z() && pt[4].z() == pt[6].z()
141  && pt[4].z() == pt[7].z()
142  && std::fabs( pt[0].z() + pt[4].z() ) < kCarTolerance
143  && pt[0].y() == pt[1].y() && pt[2].y() == pt[3].y()
144  && pt[4].y() == pt[5].y() && pt[6].y() == pt[7].y()
145  && std::fabs( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) < kCarTolerance
146  && std::fabs( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() +
147  pt[2].x() + pt[3].x() + pt[6].x() + pt[7].x() ) < kCarTolerance ) )
148  {
149  std::ostringstream message;
150  message << "Invalid vertice coordinates for Solid: " << GetName();
151  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
152  FatalException, message);
153  }
154 
155  // Bottom side with normal approx. -Y
156 
157  good = MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
158 
159  if (!good)
160  {
161  DumpInfo();
162  G4Exception("G4Trap::G4Trap()", "GeomSolids0002", FatalException,
163  "Face at ~-Y not planar.");
164  }
165 
166  // Top side with normal approx. +Y
167 
168  good = MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
169 
170  if (!good)
171  {
172  std::ostringstream message;
173  message << "Face at ~+Y not planar for Solid: " << GetName();
174  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
175  FatalException, message);
176  }
177 
178  // Front side with normal approx. -X
179 
180  good = MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
181 
182  if (!good)
183  {
184  std::ostringstream message;
185  message << "Face at ~-X not planar for Solid: " << GetName();
186  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
187  FatalException, message);
188  }
189 
190  // Back side iwth normal approx. +X
191 
192  good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
193  if (!good)
194  {
195  std::ostringstream message;
196  message << "Face at ~+X not planar for Solid: " << GetName();
197  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
198  FatalException, message);
199  }
200  fDz = (pt[7]).z() ;
201 
202  fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5;
203  fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5;
204  fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5;
205  fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1;
206 
207  fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5;
208  fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5;
209  fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5;
210  fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2;
211 
212  fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz;
213  fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;
214 }
215 
217 //
218 // Constructor for Right Angular Wedge from STEP
219 
220 G4Trap::G4Trap( const G4String& pName,
221  G4double pZ,
222  G4double pY,
223  G4double pX, G4double pLTX )
224  : G4CSGSolid(pName)
225 {
226  G4bool good;
227 
228  if ( pZ<=0 || pY<=0 || pX<=0 || pLTX<=0 || pLTX>pX )
229  {
230  std::ostringstream message;
231  message << "Invalid length parameters for Solid: " << GetName();
232  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
233  FatalException, message);
234  }
235 
236  fDz = 0.5*pZ ;
237  fTthetaCphi = 0 ;
238  fTthetaSphi = 0 ;
239 
240  fDy1 = 0.5*pY;
241  fDx1 = 0.5*pX ;
242  fDx2 = 0.5*pLTX;
243  fTalpha1 = 0.5*(pLTX - pX)/pY;
244 
245  fDy2 = fDy1 ;
246  fDx3 = fDx1;
247  fDx4 = fDx2 ;
248  fTalpha2 = fTalpha1 ;
249 
250  G4ThreeVector pt[8] ;
251 
252  pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
253  -fDz*fTthetaSphi-fDy1,-fDz);
254  pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
255  -fDz*fTthetaSphi-fDy1,-fDz);
256  pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
257  -fDz*fTthetaSphi+fDy1,-fDz);
258  pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
259  -fDz*fTthetaSphi+fDy1,-fDz);
260  pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
261  +fDz*fTthetaSphi-fDy2,+fDz);
262  pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
263  +fDz*fTthetaSphi-fDy2,+fDz);
264  pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
265  +fDz*fTthetaSphi+fDy2,+fDz);
266  pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
267  +fDz*fTthetaSphi+fDy2,+fDz);
268 
269  // Bottom side with normal approx. -Y
270  //
271  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
272  if (!good)
273  {
274  std::ostringstream message;
275  message << "Face at ~-Y not planar for Solid: " << GetName();
276  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
277  FatalException, message);
278  }
279 
280  // Top side with normal approx. +Y
281  //
282  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
283  if (!good)
284  {
285  std::ostringstream message;
286  message << "Face at ~+Y not planar for Solid: " << GetName();
287  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
288  FatalException, message);
289  }
290 
291  // Front side with normal approx. -X
292  //
293  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
294  if (!good)
295  {
296  std::ostringstream message;
297  message << "Face at ~-X not planar for Solid: " << GetName();
298  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
299  FatalException, message);
300  }
301 
302  // Back side iwth normal approx. +X
303  //
304  good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
305  if (!good)
306  {
307  std::ostringstream message;
308  message << "Face at ~+X not planar for Solid: " << GetName();
309  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
310  FatalException, message);
311  }
312 }
313 
315 //
316 // Constructor for G4Trd
317 
318 G4Trap::G4Trap( const G4String& pName,
319  G4double pDx1, G4double pDx2,
320  G4double pDy1, G4double pDy2,
321  G4double pDz )
322  : G4CSGSolid(pName)
323 {
324  G4bool good;
325 
326  if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 )
327  {
328  std::ostringstream message;
329  message << "Invalid length parameters for Solid: " << GetName();
330  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
331  FatalException, message);
332  }
333 
334  fDz = pDz;
335  fTthetaCphi = 0 ;
336  fTthetaSphi = 0 ;
337 
338  fDy1 = pDy1 ;
339  fDx1 = pDx1 ;
340  fDx2 = pDx1 ;
341  fTalpha1 = 0 ;
342 
343  fDy2 = pDy2 ;
344  fDx3 = pDx2 ;
345  fDx4 = pDx2 ;
346  fTalpha2 = 0 ;
347 
348  G4ThreeVector pt[8] ;
349 
350  pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
351  -fDz*fTthetaSphi-fDy1,-fDz);
352  pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
353  -fDz*fTthetaSphi-fDy1,-fDz);
354  pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
355  -fDz*fTthetaSphi+fDy1,-fDz);
356  pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
357  -fDz*fTthetaSphi+fDy1,-fDz);
358  pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
359  +fDz*fTthetaSphi-fDy2,+fDz);
360  pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
361  +fDz*fTthetaSphi-fDy2,+fDz);
362  pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
363  +fDz*fTthetaSphi+fDy2,+fDz);
364  pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
365  +fDz*fTthetaSphi+fDy2,+fDz);
366 
367  // Bottom side with normal approx. -Y
368  //
369  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
370  if (!good)
371  {
372  std::ostringstream message;
373  message << "Face at ~-Y not planar for Solid: " << GetName();
374  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
375  FatalException, message);
376  }
377 
378  // Top side with normal approx. +Y
379  //
380  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
381  if (!good)
382  {
383  std::ostringstream message;
384  message << "Face at ~+Y not planar for Solid: " << GetName();
385  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
386  FatalException, message);
387  }
388 
389  // Front side with normal approx. -X
390  //
391  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
392  if (!good)
393  {
394  std::ostringstream message;
395  message << "Face at ~-X not planar for Solid: " << GetName();
396  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
397  FatalException, message);
398  }
399 
400  // Back side iwth normal approx. +X
401  //
402  good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
403  if (!good)
404  {
405  std::ostringstream message;
406  message << "Face at ~+X not planar for Solid: " << GetName();
407  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
408  FatalException, message);
409  }
410 }
411 
413 //
414 // Constructor for G4Para
415 
416 G4Trap::G4Trap( const G4String& pName,
417  G4double pDx, G4double pDy,
418  G4double pDz,
419  G4double pAlpha,
420  G4double pTheta, G4double pPhi )
421  : G4CSGSolid(pName)
422 {
423  G4bool good;
424 
425  if ( pDz<=0 || pDy<=0 || pDx<=0 )
426  {
427  std::ostringstream message;
428  message << "Invalid length parameters for Solid: " << GetName();
429  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
430  FatalException, message);
431  }
432 
433  fDz = pDz ;
434  fTthetaCphi = std::tan(pTheta)*std::cos(pPhi) ;
435  fTthetaSphi = std::tan(pTheta)*std::sin(pPhi) ;
436 
437  fDy1 = pDy ;
438  fDx1 = pDx ;
439  fDx2 = pDx ;
440  fTalpha1 = std::tan(pAlpha) ;
441 
442  fDy2 = pDy ;
443  fDx3 = pDx ;
444  fDx4 = pDx ;
445  fTalpha2 = fTalpha1 ;
446 
447  G4ThreeVector pt[8] ;
448 
449  pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
450  -fDz*fTthetaSphi-fDy1,-fDz);
451  pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
452  -fDz*fTthetaSphi-fDy1,-fDz);
453  pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
454  -fDz*fTthetaSphi+fDy1,-fDz);
455  pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
456  -fDz*fTthetaSphi+fDy1,-fDz);
457  pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
458  +fDz*fTthetaSphi-fDy2,+fDz);
459  pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
460  +fDz*fTthetaSphi-fDy2,+fDz);
461  pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
462  +fDz*fTthetaSphi+fDy2,+fDz);
463  pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
464  +fDz*fTthetaSphi+fDy2,+fDz);
465 
466  // Bottom side with normal approx. -Y
467  //
468  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
469  if (!good)
470  {
471  std::ostringstream message;
472  message << "Face at ~-Y not planar for Solid: " << GetName();
473  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
474  FatalException, message);
475  }
476 
477  // Top side with normal approx. +Y
478  //
479  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
480  if (!good)
481  {
482  std::ostringstream message;
483  message << "Face at ~+Y not planar for Solid: " << GetName();
484  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
485  FatalException, message);
486  }
487 
488  // Front side with normal approx. -X
489  //
490  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
491  if (!good)
492  {
493  std::ostringstream message;
494  message << "Face at ~-X not planar for Solid: " << GetName();
495  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
496  FatalException, message);
497  }
498 
499  // Back side iwth normal approx. +X
500  //
501  good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
502  if (!good)
503  {
504  std::ostringstream message;
505  message << "Face at ~+X not planar for Solid: " << GetName();
506  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
507  FatalException, message);
508  }
509 }
510 
512 //
513 // Nominal constructor for G4Trap whose parameters are to be set by
514 // a G4VParamaterisation later. Check and set half-widths as well as
515 // angles: final check of coplanarity
516 
517 G4Trap::G4Trap( const G4String& pName )
518  : G4CSGSolid (pName), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
519  fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
520  fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
521 {
522  MakePlanes();
523 }
524 
526 //
527 // Fake default constructor - sets only member data and allocates memory
528 // for usage restricted to object persistency.
529 //
530 G4Trap::G4Trap( __void__& a )
531  : G4CSGSolid(a), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
532  fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
533  fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
534 {
535  MakePlanes();
536 }
537 
539 //
540 // Destructor
541 
542 G4Trap::~G4Trap()
543 {
544 }
545 
547 //
548 // Copy constructor
549 
550 G4Trap::G4Trap(const G4Trap& rhs)
551  : G4CSGSolid(rhs), fDz(rhs.fDz),
552  fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
553  fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
554  fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
555 {
556  for (size_t i=0; i<4; ++i)
557  {
558  fPlanes[i].a = rhs.fPlanes[i].a;
559  fPlanes[i].b = rhs.fPlanes[i].b;
560  fPlanes[i].c = rhs.fPlanes[i].c;
561  fPlanes[i].d = rhs.fPlanes[i].d;
562  }
563  fpPolyhedron = GetPolyhedron();
564 }
565 
567 //
568 // Assignment operator
569 
570 G4Trap& G4Trap::operator = (const G4Trap& rhs)
571 {
572  // Check assignment to self
573  //
574  if (this == &rhs) { return *this; }
575 
576  // Copy base class data
577  //
578  G4CSGSolid::operator=(rhs);
579 
580  // Copy data
581  //
582  fDz = rhs.fDz;
583  fTthetaCphi = rhs.fTthetaCphi; fTthetaSphi = rhs.fTthetaSphi;
584  fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
585  fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
586  for (size_t i=0; i<4; ++i)
587  {
588  fPlanes[i].a = rhs.fPlanes[i].a;
589  fPlanes[i].b = rhs.fPlanes[i].b;
590  fPlanes[i].c = rhs.fPlanes[i].c;
591  fPlanes[i].d = rhs.fPlanes[i].d;
592  }
593  fpPolyhedron = GetPolyhedron();
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  delete fpPolyhedron; 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 }
G4Trap::fDz
G4double fDz
Definition: G4Trap.hh:278
G4Trap::fDx1
G4double fDx1
Definition: G4Trap.hh:279
TrapSidePlane::b
G4double b
Definition: G4Trap.hh:103
G4VSolid::GetName
G4String GetName() const
G4Trap::fTthetaCphi
G4double fTthetaCphi
Definition: G4Trap.hh:278
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:67
G4Polyhedron.hh
G4Trap::GetPointOnPlane
G4ThreeVector GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3, G4double &area) const
Definition: G4Trap.cc:1857
G4Trap::GetEntityType
G4GeometryType GetEntityType() const
Definition: G4Trap.cc:1799
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
G4VoxelLimits::GetMinYExtent
G4double GetMinYExtent() const
G4VGraphicsScene.hh
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4Trap::G4Trap
G4Trap(const G4String &pName, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pAlp1, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlp2)
Definition: G4Trap.cc:82
G4VoxelLimits.hh
ks1
Definition: G4Trap.cc:75
G4Trap.hh
kMZ
Definition: G4Trap.cc:75
z
G4double z
Definition: TRTMaterials.hh:39
G4VoxelLimits::IsYLimited
G4bool IsYLimited() const
G4Trap::operator=
G4Trap & operator=(const G4Trap &rhs)
Definition: G4Trap.cc:570
G4AffineTransform::IsRotated
G4bool IsRotated() const
G4CSGSolid::fCubicVolume
G4double fCubicVolume
Definition: G4CSGSolid.hh:78
G4AffineTransform::NetTranslation
G4ThreeVector NetTranslation() const
G4VoxelLimits::IsXLimited
G4bool IsXLimited() const
G4Trap::fDx2
G4double fDx2
Definition: G4Trap.hh:279
a
G4double a
Definition: TRTMaterials.hh:39
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
TrapSidePlane::a
G4double a
Definition: G4Trap.hh:103
G4Trap::fTalpha2
G4double fTalpha2
Definition: G4Trap.hh:280
ks0
Definition: G4Trap.cc:75
smax
const G4int smax
Definition: G4UrQMD1_3Interface.hh:71
G4Trap::fDx4
G4double fDx4
Definition: G4Trap.hh:280
G4VSolid::DumpInfo
void DumpInfo() const
G4Trap::fDy2
G4double fDy2
Definition: G4Trap.hh:280
CLHEP::normal
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4VoxelLimits::GetMaxXExtent
G4double GetMaxXExtent() const
G4Trap::DistanceToOut
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Trap.cc:1412
G4VoxelLimits::GetMinZExtent
G4double GetMinZExtent() const
G4CSGSolid::fSurfaceArea
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79
G4cout
G4GLOB_DLL std::ostream G4cout
G4Trap::~G4Trap
virtual ~G4Trap()
Definition: G4Trap.cc:542
G4CSGSolid::fpPolyhedron
G4Polyhedron * fpPolyhedron
Definition: G4CSGSolid.hh:80
G4Trap::Clone
G4VSolid * Clone() const
Definition: G4Trap.cc:1808
G4bool
bool G4bool
Definition: G4Types.hh:79
TrapSidePlane::d
G4double d
Definition: G4Trap.hh:103
ks2
Definition: G4Trap.cc:75
G4CSGSolid::GetPolyhedron
virtual G4Polyhedron * GetPolyhedron() const
Definition: G4CSGSolid.cc:124
G4Trap::SetAllParameters
void SetAllParameters(G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pAlp1, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlp2)
Definition: G4Trap.cc:603
G4Trap::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
Definition: G4Trap.cc:816
G4VPVParameterisation.hh
G4Trap::fTalpha1
G4double fTalpha1
Definition: G4Trap.hh:279
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
Randomize.hh
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4Trap::CreateRotatedVertices
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Trap.cc:1752
G4Trap::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Trap.cc:804
G4Trap::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Trap.cc:1149
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
kPZ
Definition: G4Trap.cc:75
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:95
G4Trap::fTthetaSphi
G4double fTthetaSphi
Definition: G4Trap.hh:278
G4VoxelLimits::GetMinXExtent
G4double GetMinXExtent() const
G4AffineTransform::TransformPoint
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4Trap::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Trap.cc:1817
G4AffineTransform.hh
G4VoxelLimits::GetMaxZExtent
G4double GetMaxZExtent() const
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
EInside
EInside
Definition: geomdefs.hh:58
G4Trap::fDx3
G4double fDx3
Definition: G4Trap.hh:280
EAxis
EAxis
Definition: geomdefs.hh:54
FatalException
Definition: G4ExceptionSeverity.hh:60
G4Trap::ApproxSurfaceNormal
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Trap.cc:1229
degree
static const double degree
Definition: G4SIunits.hh:125
G4Trap::fDy1
G4double fDy1
Definition: G4Trap.hh:279
ks3
Definition: G4Trap.cc:75
TrapSidePlane::c
G4double c
Definition: G4Trap.hh:103
G4Trap::fPlanes
TrapSidePlane fPlanes[4]
Definition: G4Trap.hh:281
G4Trap::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const
Definition: G4Trap.cc:1904
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VoxelLimits::GetMaxYExtent
G4double GetMaxYExtent() const
Eside
Eside
Definition: G4Trap.cc:75
kCoplanar_Tolerance
const G4double kCoplanar_Tolerance
Definition: G4Trap.cc:69
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:305
kUndef
Definition: G4Trap.cc:75
G4Trap::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4Trap.cc:1960
G4double
double G4double
Definition: G4Types.hh:76
G4Trap::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Trap.cc:1270
G4Trap::MakePlanes
G4bool MakePlanes()
Definition: G4Trap.cc:650
G4Trap::MakePlane
G4bool MakePlane(const G4ThreeVector &p1, const G4ThreeVector &p2, const G4ThreeVector &p3, const G4ThreeVector &p4, TrapSidePlane &plane)
Definition: G4Trap.cc:729
G4Trap::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Trap.cc:1108
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:82
G4VoxelLimits::IsZLimited
G4bool IsZLimited() const
mm
static const double mm
Definition: G4SIunits.hh:102
G4Trap::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Trap.cc:1955