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