Geant4_10
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 66356 2012-12-18 09:02:32Z 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 }
564 
566 //
567 // Assignment operator
568 
569 G4Trap& G4Trap::operator = (const G4Trap& rhs)
570 {
571  // Check assignment to self
572  //
573  if (this == &rhs) { return *this; }
574 
575  // Copy base class data
576  //
577  G4CSGSolid::operator=(rhs);
578 
579  // Copy data
580  //
581  fDz = rhs.fDz;
582  fTthetaCphi = rhs.fTthetaCphi; fTthetaSphi = rhs.fTthetaSphi;
583  fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
584  fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
585  for (size_t i=0; i<4; ++i)
586  {
587  fPlanes[i].a = rhs.fPlanes[i].a;
588  fPlanes[i].b = rhs.fPlanes[i].b;
589  fPlanes[i].c = rhs.fPlanes[i].c;
590  fPlanes[i].d = rhs.fPlanes[i].d;
591  }
592 
593  return *this;
594 }
595 
597 //
598 // Set all parameters, as for constructor - check and set half-widths
599 // as well as angles: final check of coplanarity
600 
601 void G4Trap::SetAllParameters ( G4double pDz,
602  G4double pTheta,
603  G4double pPhi,
604  G4double pDy1,
605  G4double pDx1,
606  G4double pDx2,
607  G4double pAlp1,
608  G4double pDy2,
609  G4double pDx3,
610  G4double pDx4,
611  G4double pAlp2 )
612 {
613  if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 || pDx3<=0 || pDx4<=0 )
614  {
615  std::ostringstream message;
616  message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
617  << " X - "
618  << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
619  << " Y - " << pDy1 << ", " << pDy2 << G4endl
620  << " Z - " << pDz;
621  G4Exception("G4Trap::SetAllParameters()", "GeomSolids0002",
622  FatalException, message);
623  }
624  fCubicVolume= 0.;
625  fSurfaceArea= 0.;
626  fpPolyhedron = 0;
627  fDz=pDz;
628  fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
629  fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
630 
631  fDy1=pDy1;
632  fDx1=pDx1;
633  fDx2=pDx2;
634  fTalpha1=std::tan(pAlp1);
635 
636  fDy2=pDy2;
637  fDx3=pDx3;
638  fDx4=pDx4;
639  fTalpha2=std::tan(pAlp2);
640 
641  MakePlanes();
642 }
643 
645 //
646 // Checking of coplanarity
647 
648 G4bool G4Trap::MakePlanes()
649 {
650  G4bool good = true;
651 
652  G4ThreeVector pt[8] ;
653 
654  pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
655  -fDz*fTthetaSphi-fDy1,-fDz);
656  pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
657  -fDz*fTthetaSphi-fDy1,-fDz);
658  pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
659  -fDz*fTthetaSphi+fDy1,-fDz);
660  pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
661  -fDz*fTthetaSphi+fDy1,-fDz);
662  pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
663  +fDz*fTthetaSphi-fDy2,+fDz);
664  pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
665  +fDz*fTthetaSphi-fDy2,+fDz);
666  pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
667  +fDz*fTthetaSphi+fDy2,+fDz);
668  pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
669  +fDz*fTthetaSphi+fDy2,+fDz);
670 
671  // Bottom side with normal approx. -Y
672  //
673  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]) ;
674  if (!good)
675  {
676  std::ostringstream message;
677  message << "Face at ~-Y not planar for Solid: " << GetName();
678  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
679  FatalException, message);
680  }
681 
682  // Top side with normal approx. +Y
683  //
684  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
685  if (!good)
686  {
687  std::ostringstream message;
688  message << "Face at ~+Y not planar for Solid: " << GetName();
689  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
690  FatalException, message);
691  }
692 
693  // Front side with normal approx. -X
694  //
695  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
696  if (!good)
697  {
698  std::ostringstream message;
699  message << "Face at ~-X not planar for Solid: " << GetName();
700  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
701  FatalException, message);
702  }
703 
704  // Back side iwth normal approx. +X
705  //
706  good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
707  if ( !good )
708  {
709  std::ostringstream message;
710  message << "Face at ~+X not planar for Solid: " << GetName();
711  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
712  FatalException, message);
713  }
714 
715  return good;
716 }
717 
719 //
720 // Calculate the coef's of the plane p1->p2->p3->p4->p1
721 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed from
722 // infront of the plane (i.e. from normal direction).
723 //
724 // Return true if the ThreeVectors are coplanar + set coef;s
725 // false if ThreeVectors are not coplanar
726 
727 G4bool G4Trap::MakePlane( const G4ThreeVector& p1,
728  const G4ThreeVector& p2,
729  const G4ThreeVector& p3,
730  const G4ThreeVector& p4,
731  TrapSidePlane& plane )
732 {
733  G4double a, b, c, sd;
734  G4ThreeVector v12, v13, v14, Vcross;
735 
736  G4bool good;
737 
738  v12 = p2 - p1;
739  v13 = p3 - p1;
740  v14 = p4 - p1;
741  Vcross = v12.cross(v13);
742 
743  if (std::fabs(Vcross.dot(v14)/(Vcross.mag()*v14.mag())) > kCoplanar_Tolerance)
744  {
745  good = false;
746  }
747  else
748  {
749  // a,b,c correspond to the x/y/z components of the
750  // normal vector to the plane
751 
752  // a = (p2.y()-p1.y())*(p1.z()+p2.z())+(p3.y()-p2.y())*(p2.z()+p3.z());
753  // a += (p4.y()-p3.y())*(p3.z()+p4.z())+(p1.y()-p4.y())*(p4.z()+p1.z()); // ?
754  // b = (p2.z()-p1.z())*(p1.x()+p2.x())+(p3.z()-p2.z())*(p2.x()+p3.x());
755  // b += (p4.z()-p3.z())*(p3.x()+p4.x())+(p1.z()-p4.z())*(p4.x()+p1.x()); // ?
756  // c = (p2.x()-p1.x())*(p1.y()+p2.y())+(p3.x()-p2.x())*(p2.y()+p3.y());
757  // c += (p4.x()-p3.x())*(p3.y()+p4.y())+(p1.x()-p4.x())*(p4.y()+p1.y()); // ?
758 
759  // Let create diagonals 4-2 and 3-1 than (4-2)x(3-1) provides
760  // vector perpendicular to the plane directed to outside !!!
761  // and a,b,c, = f(1,2,3,4) external relative to trap normal
762 
763  a = +(p4.y() - p2.y())*(p3.z() - p1.z())
764  - (p3.y() - p1.y())*(p4.z() - p2.z());
765 
766  b = -(p4.x() - p2.x())*(p3.z() - p1.z())
767  + (p3.x() - p1.x())*(p4.z() - p2.z());
768 
769  c = +(p4.x() - p2.x())*(p3.y() - p1.y())
770  - (p3.x() - p1.x())*(p4.y() - p2.y());
771 
772  sd = std::sqrt( a*a + b*b + c*c ); // so now vector plane.(a,b,c) is unit
773 
774  if( sd > 0 )
775  {
776  plane.a = a/sd;
777  plane.b = b/sd;
778  plane.c = c/sd;
779  }
780  else
781  {
782  std::ostringstream message;
783  message << "Invalid parameters: norm.mod() <= 0, for Solid: "
784  << GetName();
785  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
786  FatalException, message) ;
787  }
788  // Calculate D: p1 in in plane so D=-n.p1.Vect()
789 
790  plane.d = -( plane.a*p1.x() + plane.b*p1.y() + plane.c*p1.z() );
791 
792  good = true;
793  }
794  return good;
795 }
796 
798 //
799 // Dispatch to parameterisation for replication mechanism dimension
800 // computation & modification.
801 
802 void G4Trap::ComputeDimensions( G4VPVParameterisation* p,
803  const G4int n,
804  const G4VPhysicalVolume* pRep )
805 {
806  p->ComputeDimensions(*this,n,pRep);
807 }
808 
809 
811 //
812 // Calculate extent under transform and specified limit
813 
814 G4bool G4Trap::CalculateExtent( const EAxis pAxis,
815  const G4VoxelLimits& pVoxelLimit,
816  const G4AffineTransform& pTransform,
817  G4double& pMin, G4double& pMax) const
818 {
819  G4double xMin, xMax, yMin, yMax, zMin, zMax;
820  G4bool flag;
821 
822  if (!pTransform.IsRotated())
823  {
824  // Special case handling for unrotated trapezoids
825  // Compute z/x/y/ mins and maxs respecting limits, with early returns
826  // if outside limits. Then switch() on pAxis
827 
828  G4int i ;
829  G4double xoffset;
830  G4double yoffset;
831  G4double zoffset;
832  G4double temp[8] ; // some points for intersection with zMin/zMax
833  G4ThreeVector pt[8]; // vertices after translation
834 
835  xoffset=pTransform.NetTranslation().x();
836  yoffset=pTransform.NetTranslation().y();
837  zoffset=pTransform.NetTranslation().z();
838 
839  pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
840  yoffset-fDz*fTthetaSphi-fDy1,zoffset-fDz);
841  pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
842  yoffset-fDz*fTthetaSphi-fDy1,zoffset-fDz);
843  pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
844  yoffset-fDz*fTthetaSphi+fDy1,zoffset-fDz);
845  pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
846  yoffset-fDz*fTthetaSphi+fDy1,zoffset-fDz);
847  pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
848  yoffset+fDz*fTthetaSphi-fDy2,zoffset+fDz);
849  pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
850  yoffset+fDz*fTthetaSphi-fDy2,zoffset+fDz);
851  pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
852  yoffset+fDz*fTthetaSphi+fDy2,zoffset+fDz);
853  pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
854  yoffset+fDz*fTthetaSphi+fDy2,zoffset+fDz);
855  zMin=zoffset-fDz;
856  zMax=zoffset+fDz;
857 
858  if ( pVoxelLimit.IsZLimited() )
859  {
860  if ( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance)
861  || (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
862  {
863  return false;
864  }
865  else
866  {
867  if ( zMin < pVoxelLimit.GetMinZExtent() )
868  {
869  zMin = pVoxelLimit.GetMinZExtent() ;
870  }
871  if ( zMax > pVoxelLimit.GetMaxZExtent() )
872  {
873  zMax = pVoxelLimit.GetMaxZExtent() ;
874  }
875  }
876  }
877  temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMin-pt[0].z())
878  /(pt[4].z()-pt[0].z()) ;
879  temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMax-pt[0].z())
880  /(pt[4].z()-pt[0].z()) ;
881  temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMin-pt[2].z())
882  /(pt[6].z()-pt[2].z()) ;
883  temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMax-pt[2].z())
884  /(pt[6].z()-pt[2].z()) ;
885 
886  yMax = yoffset - std::fabs(fDz*fTthetaSphi) - fDy1 - fDy2 ;
887  yMin = -yMax ;
888 
889  for( i = 0 ; i < 4 ; i++ )
890  {
891  if( temp[i] > yMax ) yMax = temp[i] ;
892  if( temp[i] < yMin ) yMin = temp[i] ;
893  }
894  if ( pVoxelLimit.IsYLimited() )
895  {
896  if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance)
897  || (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
898  {
899  return false;
900  }
901  else
902  {
903  if ( yMin < pVoxelLimit.GetMinYExtent() )
904  {
905  yMin = pVoxelLimit.GetMinYExtent() ;
906  }
907  if ( yMax > pVoxelLimit.GetMaxYExtent() )
908  {
909  yMax = pVoxelLimit.GetMaxYExtent() ;
910  }
911  }
912  }
913  temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())
914  *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
915  temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())
916  *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
917  temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())
918  *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
919  temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())
920  *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
921  temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())
922  *(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
923  temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())
924  *(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
925  temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())
926  *(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
927  temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())
928  *(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
929 
930  xMax = xoffset - std::fabs(fDz*fTthetaCphi) - fDx1 - fDx2 -fDx3 - fDx4 ;
931  xMin = -xMax ;
932 
933  for( i = 0 ; i < 8 ; i++ )
934  {
935  if( temp[i] > xMax) xMax = temp[i] ;
936  if( temp[i] < xMin) xMin = temp[i] ;
937  }
938  if (pVoxelLimit.IsXLimited()) // xMax/Min = f(yMax/Min) ?
939  {
940  if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance)
941  || (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
942  {
943  return false;
944  }
945  else
946  {
947  if ( xMin < pVoxelLimit.GetMinXExtent() )
948  {
949  xMin = pVoxelLimit.GetMinXExtent() ;
950  }
951  if ( xMax > pVoxelLimit.GetMaxXExtent() )
952  {
953  xMax = pVoxelLimit.GetMaxXExtent() ;
954  }
955  }
956  }
957  switch (pAxis)
958  {
959  case kXAxis:
960  pMin=xMin;
961  pMax=xMax;
962  break;
963 
964  case kYAxis:
965  pMin=yMin;
966  pMax=yMax;
967  break;
968 
969  case kZAxis:
970  pMin=zMin;
971  pMax=zMax;
972  break;
973 
974  default:
975  break;
976  }
977  pMin -= kCarTolerance;
978  pMax += kCarTolerance;
979 
980  flag = true;
981  }
982  else // General rotated case -
983  {
984  G4bool existsAfterClip = false ;
985  G4ThreeVectorList* vertices;
986  pMin = +kInfinity;
987  pMax = -kInfinity;
988 
989  // Calculate rotated vertex coordinates. Operator 'new' is called
990 
991  vertices = CreateRotatedVertices(pTransform);
992 
993  xMin = +kInfinity; yMin = +kInfinity; zMin = +kInfinity;
994  xMax = -kInfinity; yMax = -kInfinity; zMax = -kInfinity;
995 
996  for( G4int nv = 0 ; nv < 8 ; nv++ )
997  {
998  if( (*vertices)[nv].x() > xMax ) xMax = (*vertices)[nv].x();
999  if( (*vertices)[nv].y() > yMax ) yMax = (*vertices)[nv].y();
1000  if( (*vertices)[nv].z() > zMax ) zMax = (*vertices)[nv].z();
1001 
1002  if( (*vertices)[nv].x() < xMin ) xMin = (*vertices)[nv].x();
1003  if( (*vertices)[nv].y() < yMin ) yMin = (*vertices)[nv].y();
1004  if( (*vertices)[nv].z() < zMin ) zMin = (*vertices)[nv].z();
1005  }
1006  if ( pVoxelLimit.IsZLimited() )
1007  {
1008  if ( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance)
1009  || (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
1010  {
1011  delete vertices ; // 'new' in the function called
1012  return false;
1013  }
1014  else
1015  {
1016  if ( zMin < pVoxelLimit.GetMinZExtent() )
1017  {
1018  zMin = pVoxelLimit.GetMinZExtent() ;
1019  }
1020  if ( zMax > pVoxelLimit.GetMaxZExtent() )
1021  {
1022  zMax = pVoxelLimit.GetMaxZExtent() ;
1023  }
1024  }
1025  }
1026  if ( pVoxelLimit.IsYLimited() )
1027  {
1028  if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance)
1029  || (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
1030  {
1031  delete vertices ; // 'new' in the function called
1032  return false;
1033  }
1034  else
1035  {
1036  if ( yMin < pVoxelLimit.GetMinYExtent() )
1037  {
1038  yMin = pVoxelLimit.GetMinYExtent() ;
1039  }
1040  if ( yMax > pVoxelLimit.GetMaxYExtent() )
1041  {
1042  yMax = pVoxelLimit.GetMaxYExtent() ;
1043  }
1044  }
1045  }
1046  if ( pVoxelLimit.IsXLimited() )
1047  {
1048  if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance)
1049  || (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
1050  {
1051  delete vertices ; // 'new' in the function called
1052  return false ;
1053  }
1054  else
1055  {
1056  if ( xMin < pVoxelLimit.GetMinXExtent() )
1057  {
1058  xMin = pVoxelLimit.GetMinXExtent() ;
1059  }
1060  if ( xMax > pVoxelLimit.GetMaxXExtent() )
1061  {
1062  xMax = pVoxelLimit.GetMaxXExtent() ;
1063  }
1064  }
1065  }
1066  switch (pAxis)
1067  {
1068  case kXAxis:
1069  pMin=xMin;
1070  pMax=xMax;
1071  break;
1072 
1073  case kYAxis:
1074  pMin=yMin;
1075  pMax=yMax;
1076  break;
1077 
1078  case kZAxis:
1079  pMin=zMin;
1080  pMax=zMax;
1081  break;
1082 
1083  default:
1084  break;
1085  }
1086  if ( (pMin != kInfinity) || (pMax != -kInfinity) )
1087  {
1088  existsAfterClip=true;
1089 
1090  // Add tolerance to avoid precision troubles
1091  //
1092  pMin -= kCarTolerance ;
1093  pMax += kCarTolerance ;
1094  }
1095  delete vertices ; // 'new' in the function called
1096  flag = existsAfterClip ;
1097  }
1098  return flag;
1099 }
1100 
1101 
1103 //
1104 // Return whether point inside/outside/on surface, using tolerance
1105 
1106 EInside G4Trap::Inside( const G4ThreeVector& p ) const
1107 {
1108  EInside in;
1109  G4double Dist;
1110  G4int i;
1111  if ( std::fabs(p.z()) <= fDz-kCarTolerance*0.5)
1112  {
1113  in = kInside;
1114 
1115  for ( i = 0;i < 4;i++ )
1116  {
1117  Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
1118  +fPlanes[i].c*p.z() + fPlanes[i].d;
1119 
1120  if (Dist > kCarTolerance*0.5) return in = kOutside;
1121  else if (Dist > -kCarTolerance*0.5) in = kSurface;
1122 
1123  }
1124  }
1125  else if (std::fabs(p.z()) <= fDz+kCarTolerance*0.5)
1126  {
1127  in = kSurface;
1128 
1129  for ( i = 0; i < 4; i++ )
1130  {
1131  Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
1132  +fPlanes[i].c*p.z() + fPlanes[i].d;
1133 
1134  if (Dist > kCarTolerance*0.5) return in = kOutside;
1135  }
1136  }
1137  else in = kOutside;
1138 
1139  return in;
1140 }
1141 
1143 //
1144 // Calculate side nearest to p, and return normal
1145 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
1146 
1147 G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const
1148 {
1149  G4int i, noSurfaces = 0;
1150  G4double dist, distz, distx, disty, distmx, distmy, safe = kInfinity;
1151  G4double delta = 0.5*kCarTolerance;
1152  G4ThreeVector norm, sumnorm(0.,0.,0.);
1153 
1154  for (i = 0; i < 4; i++)
1155  {
1156  dist = std::fabs(fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
1157  + fPlanes[i].c*p.z() + fPlanes[i].d);
1158  if ( dist < safe )
1159  {
1160  safe = dist;
1161  }
1162  }
1163  distz = std::fabs( std::fabs( p.z() ) - fDz );
1164 
1165  distmy = std::fabs( fPlanes[0].a*p.x() + fPlanes[0].b*p.y()
1166  + fPlanes[0].c*p.z() + fPlanes[0].d );
1167 
1168  disty = std::fabs( fPlanes[1].a*p.x() + fPlanes[1].b*p.y()
1169  + fPlanes[1].c*p.z() + fPlanes[1].d );
1170 
1171  distmx = std::fabs( fPlanes[2].a*p.x() + fPlanes[2].b*p.y()
1172  + fPlanes[2].c*p.z() + fPlanes[2].d );
1173 
1174  distx = std::fabs( fPlanes[3].a*p.x() + fPlanes[3].b*p.y()
1175  + fPlanes[3].c*p.z() + fPlanes[3].d );
1176 
1177  G4ThreeVector nX = G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1178  G4ThreeVector nmX = G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1179  G4ThreeVector nY = G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1180  G4ThreeVector nmY = G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1181  G4ThreeVector nZ = G4ThreeVector(0.,0.,1.0);
1182 
1183  if (distx <= delta)
1184  {
1185  noSurfaces ++;
1186  sumnorm += nX;
1187  }
1188  if (distmx <= delta)
1189  {
1190  noSurfaces ++;
1191  sumnorm += nmX;
1192  }
1193  if (disty <= delta)
1194  {
1195  noSurfaces ++;
1196  sumnorm += nY;
1197  }
1198  if (distmy <= delta)
1199  {
1200  noSurfaces ++;
1201  sumnorm += nmY;
1202  }
1203  if (distz <= delta)
1204  {
1205  noSurfaces ++;
1206  if ( p.z() >= 0.) sumnorm += nZ;
1207  else sumnorm -= nZ;
1208  }
1209  if ( noSurfaces == 0 )
1210  {
1211 #ifdef G4CSGDEBUG
1212  G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
1213  JustWarning, "Point p is not on surface !?" );
1214 #endif
1215  norm = ApproxSurfaceNormal(p);
1216  }
1217  else if ( noSurfaces == 1 ) norm = sumnorm;
1218  else norm = sumnorm.unit();
1219  return norm;
1220 }
1221 
1223 //
1224 // Algorithm for SurfaceNormal() following the original specification
1225 // for points not on the surface
1226 
1227 G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
1228 {
1229  G4double safe=kInfinity,Dist,safez;
1230  G4int i,imin=0;
1231  for (i=0;i<4;i++)
1232  {
1233  Dist=std::fabs(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1234  +fPlanes[i].c*p.z()+fPlanes[i].d);
1235  if (Dist<safe)
1236  {
1237  safe=Dist;
1238  imin=i;
1239  }
1240  }
1241  safez=std::fabs(std::fabs(p.z())-fDz);
1242  if (safe<safez)
1243  {
1244  return G4ThreeVector(fPlanes[imin].a,fPlanes[imin].b,fPlanes[imin].c);
1245  }
1246  else
1247  {
1248  if (p.z()>0)
1249  {
1250  return G4ThreeVector(0,0,1);
1251  }
1252  else
1253  {
1254  return G4ThreeVector(0,0,-1);
1255  }
1256  }
1257 }
1258 
1260 //
1261 // Calculate distance to shape from outside - return kInfinity if no intersection
1262 //
1263 // ALGORITHM:
1264 // For each component, calculate pair of minimum and maximum intersection
1265 // values for which the particle is in the extent of the shape
1266 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
1267 
1268 G4double G4Trap::DistanceToIn( const G4ThreeVector& p,
1269  const G4ThreeVector& v ) const
1270 {
1271 
1272  G4double snxt; // snxt = default return value
1273  G4double max,smax,smin;
1274  G4double pdist,Comp,vdist;
1275  G4int i;
1276  //
1277  // Z Intersection range
1278  //
1279  if ( v.z() > 0 )
1280  {
1281  max = fDz - p.z() ;
1282  if (max > 0.5*kCarTolerance)
1283  {
1284  smax = max/v.z();
1285  smin = (-fDz-p.z())/v.z();
1286  }
1287  else
1288  {
1289  return snxt=kInfinity;
1290  }
1291  }
1292  else if (v.z() < 0 )
1293  {
1294  max = - fDz - p.z() ;
1295  if (max < -0.5*kCarTolerance )
1296  {
1297  smax=max/v.z();
1298  smin=(fDz-p.z())/v.z();
1299  }
1300  else
1301  {
1302  return snxt=kInfinity;
1303  }
1304  }
1305  else
1306  {
1307  if (std::fabs(p.z())<fDz - 0.5*kCarTolerance) // Inside was <=fDz
1308  {
1309  smin=0;
1310  smax=kInfinity;
1311  }
1312  else
1313  {
1314  return snxt=kInfinity;
1315  }
1316  }
1317 
1318  for (i=0;i<4;i++)
1319  {
1320  pdist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1321  +fPlanes[i].c*p.z()+fPlanes[i].d;
1322  Comp=fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
1323  if ( pdist >= -0.5*kCarTolerance ) // was >0
1324  {
1325  //
1326  // Outside the plane -> this is an extent entry distance
1327  //
1328  if (Comp >= 0) // was >0
1329  {
1330  return snxt=kInfinity ;
1331  }
1332  else
1333  {
1334  vdist=-pdist/Comp;
1335  if (vdist>smin)
1336  {
1337  if (vdist<smax)
1338  {
1339  smin = vdist;
1340  }
1341  else
1342  {
1343  return snxt=kInfinity;
1344  }
1345  }
1346  }
1347  }
1348  else
1349  {
1350  //
1351  // Inside the plane -> couble be an extent exit distance (smax)
1352  //
1353  if (Comp>0) // Will leave extent
1354  {
1355  vdist=-pdist/Comp;
1356  if (vdist<smax)
1357  {
1358  if (vdist>smin)
1359  {
1360  smax=vdist;
1361  }
1362  else
1363  {
1364  return snxt=kInfinity;
1365  }
1366  }
1367  }
1368  }
1369  }
1370  //
1371  // Checks in non z plane intersections ensure smin<smax
1372  //
1373  if (smin >=0 )
1374  {
1375  snxt = smin ;
1376  }
1377  else
1378  {
1379  snxt = 0 ;
1380  }
1381  return snxt;
1382 }
1383 
1385 //
1386 // Calculate exact shortest distance to any boundary from outside
1387 // This is the best fast estimation of the shortest distance to trap
1388 // - Returns 0 is ThreeVector inside
1389 
1390 G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const
1391 {
1392  G4double safe=0.0,Dist;
1393  G4int i;
1394  safe=std::fabs(p.z())-fDz;
1395  for (i=0;i<4;i++)
1396  {
1397  Dist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1398  +fPlanes[i].c*p.z()+fPlanes[i].d;
1399  if (Dist > safe) safe=Dist;
1400  }
1401  if (safe<0) safe=0;
1402  return safe;
1403 }
1404 
1406 //
1407 // Calculate distance to surface of shape from inside
1408 // Calculate distance to x/y/z planes - smallest is exiting distance
1409 
1410 G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
1411  const G4bool calcNorm,
1412  G4bool *validNorm, G4ThreeVector *n) const
1413 {
1414  Eside side = kUndef;
1415  G4double snxt; // snxt = return value
1416  G4double pdist,Comp,vdist,max;
1417  //
1418  // Z Intersections
1419  //
1420  if (v.z()>0)
1421  {
1422  max=fDz-p.z();
1423  if (max>kCarTolerance/2)
1424  {
1425  snxt=max/v.z();
1426  side=kPZ;
1427  }
1428  else
1429  {
1430  if (calcNorm)
1431  {
1432  *validNorm=true;
1433  *n=G4ThreeVector(0,0,1);
1434  }
1435  return snxt=0;
1436  }
1437  }
1438  else if (v.z()<0)
1439  {
1440  max=-fDz-p.z();
1441  if (max<-kCarTolerance/2)
1442  {
1443  snxt=max/v.z();
1444  side=kMZ;
1445  }
1446  else
1447  {
1448  if (calcNorm)
1449  {
1450  *validNorm=true;
1451  *n=G4ThreeVector(0,0,-1);
1452  }
1453  return snxt=0;
1454  }
1455  }
1456  else
1457  {
1458  snxt=kInfinity;
1459  }
1460 
1461  //
1462  // Intersections with planes[0] (expanded because of setting enum)
1463  //
1464  pdist=fPlanes[0].a*p.x()+fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
1465  Comp=fPlanes[0].a*v.x()+fPlanes[0].b*v.y()+fPlanes[0].c*v.z();
1466  if (pdist>0)
1467  {
1468  // Outside the plane
1469  if (Comp>0)
1470  {
1471  // Leaving immediately
1472  if (calcNorm)
1473  {
1474  *validNorm=true;
1475  *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1476  }
1477  return snxt=0;
1478  }
1479  }
1480  else if (pdist<-kCarTolerance/2)
1481  {
1482  // Inside the plane
1483  if (Comp>0)
1484  {
1485  // Will leave extent
1486  vdist=-pdist/Comp;
1487  if (vdist<snxt)
1488  {
1489  snxt=vdist;
1490  side=ks0;
1491  }
1492  }
1493  }
1494  else
1495  {
1496  // On surface
1497  if (Comp>0)
1498  {
1499  if (calcNorm)
1500  {
1501  *validNorm=true;
1502  *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1503  }
1504  return snxt=0;
1505  }
1506  }
1507 
1508  //
1509  // Intersections with planes[1] (expanded because of setting enum)
1510  //
1511  pdist=fPlanes[1].a*p.x()+fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
1512  Comp=fPlanes[1].a*v.x()+fPlanes[1].b*v.y()+fPlanes[1].c*v.z();
1513  if (pdist>0)
1514  {
1515  // Outside the plane
1516  if (Comp>0)
1517  {
1518  // Leaving immediately
1519  if (calcNorm)
1520  {
1521  *validNorm=true;
1522  *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1523  }
1524  return snxt=0;
1525  }
1526  }
1527  else if (pdist<-kCarTolerance/2)
1528  {
1529  // Inside the plane
1530  if (Comp>0)
1531  {
1532  // Will leave extent
1533  vdist=-pdist/Comp;
1534  if (vdist<snxt)
1535  {
1536  snxt=vdist;
1537  side=ks1;
1538  }
1539  }
1540  }
1541  else
1542  {
1543  // On surface
1544  if (Comp>0)
1545  {
1546  if (calcNorm)
1547  {
1548  *validNorm=true;
1549  *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1550  }
1551  return snxt=0;
1552  }
1553  }
1554 
1555  //
1556  // Intersections with planes[2] (expanded because of setting enum)
1557  //
1558  pdist=fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d;
1559  Comp=fPlanes[2].a*v.x()+fPlanes[2].b*v.y()+fPlanes[2].c*v.z();
1560  if (pdist>0)
1561  {
1562  // Outside the plane
1563  if (Comp>0)
1564  {
1565  // Leaving immediately
1566  if (calcNorm)
1567  {
1568  *validNorm=true;
1569  *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1570  }
1571  return snxt=0;
1572  }
1573  }
1574  else if (pdist<-kCarTolerance/2)
1575  {
1576  // Inside the plane
1577  if (Comp>0)
1578  {
1579  // Will leave extent
1580  vdist=-pdist/Comp;
1581  if (vdist<snxt)
1582  {
1583  snxt=vdist;
1584  side=ks2;
1585  }
1586  }
1587  }
1588  else
1589  {
1590  // On surface
1591  if (Comp>0)
1592  {
1593  if (calcNorm)
1594  {
1595  *validNorm=true;
1596  *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1597  }
1598  return snxt=0;
1599  }
1600  }
1601 
1602  //
1603  // Intersections with planes[3] (expanded because of setting enum)
1604  //
1605  pdist=fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d;
1606  Comp=fPlanes[3].a*v.x()+fPlanes[3].b*v.y()+fPlanes[3].c*v.z();
1607  if (pdist>0)
1608  {
1609  // Outside the plane
1610  if (Comp>0)
1611  {
1612  // Leaving immediately
1613  if (calcNorm)
1614  {
1615  *validNorm=true;
1616  *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1617  }
1618  return snxt=0;
1619  }
1620  }
1621  else if (pdist<-kCarTolerance/2)
1622  {
1623  // Inside the plane
1624  if (Comp>0)
1625  {
1626  // Will leave extent
1627  vdist=-pdist/Comp;
1628  if (vdist<snxt)
1629  {
1630  snxt=vdist;
1631  side=ks3;
1632  }
1633  }
1634  }
1635  else
1636  {
1637  // On surface
1638  if (Comp>0)
1639  {
1640  if (calcNorm)
1641  {
1642  *validNorm=true;
1643  *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1644  }
1645  return snxt=0;
1646  }
1647  }
1648 
1649  // set normal
1650  if (calcNorm)
1651  {
1652  *validNorm=true;
1653  switch(side)
1654  {
1655  case ks0:
1656  *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1657  break;
1658  case ks1:
1659  *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1660  break;
1661  case ks2:
1662  *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1663  break;
1664  case ks3:
1665  *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1666  break;
1667  case kMZ:
1668  *n=G4ThreeVector(0,0,-1);
1669  break;
1670  case kPZ:
1671  *n=G4ThreeVector(0,0,1);
1672  break;
1673  default:
1674  G4cout << G4endl;
1675  DumpInfo();
1676  std::ostringstream message;
1677  G4int oldprc = message.precision(16);
1678  message << "Undefined side for valid surface normal to solid."
1679  << G4endl
1680  << "Position:" << G4endl << G4endl
1681  << "p.x() = " << p.x()/mm << " mm" << G4endl
1682  << "p.y() = " << p.y()/mm << " mm" << G4endl
1683  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1684  << "Direction:" << G4endl << G4endl
1685  << "v.x() = " << v.x() << G4endl
1686  << "v.y() = " << v.y() << G4endl
1687  << "v.z() = " << v.z() << G4endl << G4endl
1688  << "Proposed distance :" << G4endl << G4endl
1689  << "snxt = " << snxt/mm << " mm" << G4endl;
1690  message.precision(oldprc);
1691  G4Exception("G4Trap::DistanceToOut(p,v,..)","GeomSolids1002",
1692  JustWarning, message);
1693  break;
1694  }
1695  }
1696  return snxt;
1697 }
1698 
1700 //
1701 // Calculate exact shortest distance to any boundary from inside
1702 // - Returns 0 is ThreeVector outside
1703 
1704 G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const
1705 {
1706  G4double safe=0.0,Dist;
1707  G4int i;
1708 
1709 #ifdef G4CSGDEBUG
1710  if( Inside(p) == kOutside )
1711  {
1712  G4int oldprc = G4cout.precision(16) ;
1713  G4cout << G4endl ;
1714  DumpInfo();
1715  G4cout << "Position:" << G4endl << G4endl ;
1716  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1717  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1718  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1719  G4cout.precision(oldprc) ;
1720  G4Exception("G4Trap::DistanceToOut(p)",
1721  "GeomSolids1002", JustWarning, "Point p is outside !?" );
1722  }
1723 #endif
1724 
1725  safe=fDz-std::fabs(p.z());
1726  if (safe<0) safe=0;
1727  else
1728  {
1729  for (i=0;i<4;i++)
1730  {
1731  Dist=-(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1732  +fPlanes[i].c*p.z()+fPlanes[i].d);
1733  if (Dist<safe) safe=Dist;
1734  }
1735  if (safe<0) safe=0;
1736  }
1737  return safe;
1738 }
1739 
1741 //
1742 // Create a List containing the transformed vertices
1743 // Ordering [0-3] -fDz cross section
1744 // [4-7] +fDz cross section such that [0] is below [4],
1745 // [1] below [5] etc.
1746 // Note:
1747 // Caller has deletion resposibility
1748 
1749 G4ThreeVectorList*
1750 G4Trap::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
1751 {
1752  G4ThreeVectorList *vertices;
1753  vertices=new G4ThreeVectorList();
1754  if (vertices)
1755  {
1756  vertices->reserve(8);
1757  G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
1758  -fDz*fTthetaSphi-fDy1,-fDz);
1759  G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
1760  -fDz*fTthetaSphi-fDy1,-fDz);
1761  G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
1762  -fDz*fTthetaSphi+fDy1,-fDz);
1763  G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
1764  -fDz*fTthetaSphi+fDy1,-fDz);
1765  G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
1766  +fDz*fTthetaSphi-fDy2,+fDz);
1767  G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
1768  +fDz*fTthetaSphi-fDy2,+fDz);
1769  G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
1770  +fDz*fTthetaSphi+fDy2,+fDz);
1771  G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
1772  +fDz*fTthetaSphi+fDy2,+fDz);
1773 
1774  vertices->push_back(pTransform.TransformPoint(vertex0));
1775  vertices->push_back(pTransform.TransformPoint(vertex1));
1776  vertices->push_back(pTransform.TransformPoint(vertex2));
1777  vertices->push_back(pTransform.TransformPoint(vertex3));
1778  vertices->push_back(pTransform.TransformPoint(vertex4));
1779  vertices->push_back(pTransform.TransformPoint(vertex5));
1780  vertices->push_back(pTransform.TransformPoint(vertex6));
1781  vertices->push_back(pTransform.TransformPoint(vertex7));
1782  }
1783  else
1784  {
1785  DumpInfo();
1786  G4Exception("G4Trap::CreateRotatedVertices()",
1787  "GeomSolids0003", FatalException,
1788  "Error in allocation of vertices. Out of memory !");
1789  }
1790  return vertices;
1791 }
1792 
1794 //
1795 // GetEntityType
1796 
1797 G4GeometryType G4Trap::GetEntityType() const
1798 {
1799  return G4String("G4Trap");
1800 }
1801 
1803 //
1804 // Make a clone of the object
1805 //
1806 G4VSolid* G4Trap::Clone() const
1807 {
1808  return new G4Trap(*this);
1809 }
1810 
1812 //
1813 // Stream object contents to an output stream
1814 
1815 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1816 {
1817  G4int oldprc = os.precision(16);
1818  os << "-----------------------------------------------------------\n"
1819  << " *** Dump for solid - " << GetName() << " ***\n"
1820  << " ===================================================\n"
1821  << " Solid type: G4Trap\n"
1822  << " Parameters: \n"
1823  << " half length Z: " << fDz/mm << " mm \n"
1824  << " half length Y of face -fDz: " << fDy1/mm << " mm \n"
1825  << " half length X of side -fDy1, face -fDz: " << fDx1/mm << " mm \n"
1826  << " half length X of side +fDy1, face -fDz: " << fDx2/mm << " mm \n"
1827  << " half length Y of face +fDz: " << fDy2/mm << " mm \n"
1828  << " half length X of side -fDy2, face +fDz: " << fDx3/mm << " mm \n"
1829  << " half length X of side +fDy2, face +fDz: " << fDx4/mm << " mm \n"
1830  << " std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree << " degrees \n"
1831  << " std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree << " degrees \n"
1832  << " std::tan(alpha), -fDz: " << fTalpha1/degree << " degrees \n"
1833  << " std::tan(alpha), +fDz: " << fTalpha2/degree << " degrees \n"
1834  << " trap side plane equations:\n"
1835  << " " << fPlanes[0].a << " X + " << fPlanes[0].b << " Y + "
1836  << fPlanes[0].c << " Z + " << fPlanes[0].d << " = 0\n"
1837  << " " << fPlanes[1].a << " X + " << fPlanes[1].b << " Y + "
1838  << fPlanes[1].c << " Z + " << fPlanes[1].d << " = 0\n"
1839  << " " << fPlanes[2].a << " X + " << fPlanes[2].b << " Y + "
1840  << fPlanes[2].c << " Z + " << fPlanes[2].d << " = 0\n"
1841  << " " << fPlanes[3].a << " X + " << fPlanes[3].b << " Y + "
1842  << fPlanes[3].c << " Z + " << fPlanes[3].d << " = 0\n"
1843  << "-----------------------------------------------------------\n";
1844  os.precision(oldprc);
1845 
1846  return os;
1847 }
1848 
1850 //
1851 // GetPointOnPlane
1852 //
1853 // Auxiliary method for Get Point on Surface
1854 
1855 G4ThreeVector G4Trap::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1,
1856  G4ThreeVector p2, G4ThreeVector p3,
1857  G4double& area) const
1858 {
1859  G4double lambda1, lambda2, chose, aOne, aTwo;
1860  G4ThreeVector t, u, v, w, Area, normal;
1861 
1862  t = p1 - p0;
1863  u = p2 - p1;
1864  v = p3 - p2;
1865  w = p0 - p3;
1866 
1867  Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
1868  w.z()*v.x() - w.x()*v.z(),
1869  w.x()*v.y() - w.y()*v.x());
1870 
1871  aOne = 0.5*Area.mag();
1872 
1873  Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
1874  t.z()*u.x() - t.x()*u.z(),
1875  t.x()*u.y() - t.y()*u.x());
1876 
1877  aTwo = 0.5*Area.mag();
1878 
1879  area = aOne + aTwo;
1880 
1881  chose = RandFlat::shoot(0.,aOne+aTwo);
1882 
1883  if( (chose>=0.) && (chose < aOne) )
1884  {
1885  lambda1 = RandFlat::shoot(0.,1.);
1886  lambda2 = RandFlat::shoot(0.,lambda1);
1887  return (p2+lambda1*v+lambda2*w);
1888  }
1889 
1890  // else
1891 
1892  lambda1 = RandFlat::shoot(0.,1.);
1893  lambda2 = RandFlat::shoot(0.,lambda1);
1894 
1895  return (p0+lambda1*t+lambda2*u);
1896 }
1897 
1899 //
1900 // GetPointOnSurface
1901 
1902 G4ThreeVector G4Trap::GetPointOnSurface() const
1903 {
1904  G4double aOne, aTwo, aThree, aFour, aFive, aSix, chose;
1905  G4ThreeVector One, Two, Three, Four, Five, Six, test;
1906  G4ThreeVector pt[8];
1907 
1908  pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
1909  -fDz*fTthetaSphi-fDy1,-fDz);
1910  pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
1911  -fDz*fTthetaSphi-fDy1,-fDz);
1912  pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
1913  -fDz*fTthetaSphi+fDy1,-fDz);
1914  pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
1915  -fDz*fTthetaSphi+fDy1,-fDz);
1916  pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
1917  +fDz*fTthetaSphi-fDy2,+fDz);
1918  pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
1919  +fDz*fTthetaSphi-fDy2,+fDz);
1920  pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
1921  +fDz*fTthetaSphi+fDy2,+fDz);
1922  pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
1923  +fDz*fTthetaSphi+fDy2,+fDz);
1924 
1925  // make sure we provide the points in a clockwise fashion
1926 
1927  One = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
1928  Two = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
1929  Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
1930  Four = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour);
1931  Five = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
1932  Six = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
1933 
1934  chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
1935  if( (chose>=0.) && (chose<aOne) )
1936  { return One; }
1937  else if( (chose>=aOne) && (chose<aOne+aTwo) )
1938  { return Two; }
1939  else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) )
1940  { return Three; }
1941  else if( (chose>=aOne+aTwo+aThree) && (chose<aOne+aTwo+aThree+aFour) )
1942  { return Four; }
1943  else if( (chose>=aOne+aTwo+aThree+aFour)
1944  && (chose<aOne+aTwo+aThree+aFour+aFive) )
1945  { return Five; }
1946  return Six;
1947 }
1948 
1950 //
1951 // Methods for visualisation
1952 
1953 void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1954 {
1955  scene.AddSolid (*this);
1956 }
1957 
1958 G4Polyhedron* G4Trap::CreatePolyhedron () const
1959 {
1960  G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1961  G4double alpha1 = std::atan(fTalpha1);
1962  G4double alpha2 = std::atan(fTalpha2);
1963  G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi+fTthetaSphi*fTthetaSphi));
1964 
1965  return new G4PolyhedronTrap(fDz, theta, phi,
1966  fDy1, fDx1, fDx2, alpha1,
1967  fDy2, fDx3, fDx4, alpha2);
1968 }
TrapSidePlane::b
G4double b
Definition: G4Trap.hh:103
G4VSolid::GetName
G4String GetName() const
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:67
G4Polyhedron.hh
kMZ
Definition: G4Cons.cc:68
G4Trap::GetEntityType
G4GeometryType GetEntityType() const
Definition: G4Trap.cc:1797
test.a
tuple a
Definition: test.py:11
G4VoxelLimits::GetMinYExtent
G4double GetMinYExtent() const
G4VGraphicsScene.hh
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
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
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
G4VoxelLimits.hh
ks1
Definition: G4Trap.cc:75
G4Trap.hh
G4VoxelLimits::IsYLimited
G4bool IsYLimited() const
in
ifstream in
Definition: comparison.C:7
p
const char * p
Definition: xmltok.h:285
norm
Float_t norm
Definition: comparison_ascii.C:71
G4Trap::operator=
G4Trap & operator=(const G4Trap &rhs)
Definition: G4Trap.cc:569
G4AffineTransform::IsRotated
G4bool IsRotated() const
G4CSGSolid::fCubicVolume
G4double fCubicVolume
Definition: G4CSGSolid.hh:78
G4AffineTransform::NetTranslation
G4ThreeVector NetTranslation() const
G4VoxelLimits::IsXLimited
G4bool IsXLimited() const
test.x
tuple x
Definition: test.py:50
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
ks0
Definition: G4Trap.cc:75
smax
const G4int smax
Definition: G4UrQMD1_3Interface.hh:71
CLHEP::Hep3Vector::z
double z() const
G4VSolid::DumpInfo
void DumpInfo() const
y
Double_t y
Definition: plot.C:279
test.b
tuple b
Definition: test.py:12
G4VoxelLimits::GetMaxXExtent
G4double GetMaxXExtent() const
n
Char_t n[5]
Definition: comparison_ascii.C:55
G4Trap::DistanceToOut
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Trap.cc:1410
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
python.hepunit.degree
tuple degree
Definition: hepunit.py:69
G4Trap::Clone
G4VSolid * Clone() const
Definition: G4Trap.cc:1806
G4bool
bool G4bool
Definition: G4Types.hh:79
TrapSidePlane::d
G4double d
Definition: G4Trap.hh:103
ks2
Definition: G4Trap.cc:75
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:601
G4Trap::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
Definition: G4Trap.cc:814
G4VPVParameterisation.hh
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
kUndef
Definition: G4Para.cc:61
Randomize.hh
G4Trap::CreateRotatedVertices
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Trap.cc:1750
test.v
tuple v
Definition: test.py:18
G4Trap::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Trap.cc:802
pt
TMarker * pt
Definition: egs.C:22
G4Trap::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Trap.cc:1147
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:95
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:1815
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
EAxis
EAxis
Definition: geomdefs.hh:54
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::y
double y() const
FatalException
Definition: G4ExceptionSeverity.hh:60
kPZ
Definition: G4Cons.cc:68
test.z
tuple z
Definition: test.py:28
ks3
Definition: G4Trap.cc:75
TrapSidePlane::c
G4double c
Definition: G4Trap.hh:103
G4Trap::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const
Definition: G4Trap.cc:1902
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
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
mcscore.test
def test
Definition: mcscore.py:117
G4Trap::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4Trap.cc:1958
G4double
double G4double
Definition: G4Types.hh:76
G4Trap::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Trap.cc:1268
test.c
tuple c
Definition: test.py:13
G4Trap::MakePlanes
G4bool MakePlanes()
Definition: G4Trap.cc:648
G4Trap::MakePlane
G4bool MakePlane(const G4ThreeVector &p1, const G4ThreeVector &p2, const G4ThreeVector &p3, const G4ThreeVector &p4, TrapSidePlane &plane)
Definition: G4Trap.cc:727
G4Trap::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Trap.cc:1106
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:82
CLHEP::Hep3Vector::mag
double mag() const
G4VoxelLimits::IsZLimited
G4bool IsZLimited() const
G4Trap::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Trap.cc:1953
Generated on Sat Dec 14 2013 14:34:10 for Geant4_10 by   doxygen 1.8.5