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