Geant4  10.03
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 101121 2016-11-07 09:18:01Z gcosmo $
28 //
29 // class G4Trap
30 //
31 // Implementation for G4Trap class
32 //
33 // History:
34 //
35 // 23.09.16 E.Tcherniaev: added Extent(pmin,pmax),
36 // use G4BoundingEnvelope for CalculateExtent(),
37 // removed CreateRotatedVertices()
38 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal
39 // 26.04.05 V.Grichine: new SurfaceNormal is default
40 // 19.04.05 V.Grichine: bug fixed in G4Trap("name",G4ThreeVector[8] vp)
41 // 12.12.04 V.Grichine: SurfaceNormal with edges/vertices
42 // 15.11.04 V.Grichine: bug fixed in G4Trap("name",G4ThreeVector[8] vp)
43 // 13.12.99 V.Grichine: bug fixed in DistanceToIn(p,v)
44 // 19.11.99 V.Grichine: kUndef was added to Eside enum
45 // 04.06.99 S.Giani: Fixed CalculateExtent in rotated case.
46 // 08.12.97 J.Allison: Added "nominal" constructor and method SetAllParameters.
47 // 01.11.96 V.Grichine: Costructor for Right Angular Wedge from STEP, G4Trd/Para
48 // 09.09.96 V.Grichine: Final modifications before to commit
49 // 21.03.95 P.Kent: Modified for `tolerant' geometry
50 //
52 
53 #include "G4Trap.hh"
54 
55 #if !defined(G4GEOM_USE_UTRAP)
56 
57 #include "globals.hh"
58 
59 #include "G4VoxelLimits.hh"
60 #include "G4AffineTransform.hh"
61 #include "G4BoundingEnvelope.hh"
62 
63 #include "G4VPVParameterisation.hh"
64 
65 #include "Randomize.hh"
66 
67 #include "G4VGraphicsScene.hh"
68 #include "G4Polyhedron.hh"
69 
70 using namespace CLHEP;
71 
73 //
74 // Accuracy of coplanarity
75 
77 
79 //
80 // Private enum: Not for external use
81 
83 
85 //
86 // Constructor - check and set half-widths as well as angles:
87 // final check of coplanarity
88 
89 G4Trap::G4Trap( const G4String& pName,
90  G4double pDz,
91  G4double pTheta, G4double pPhi,
92  G4double pDy1, G4double pDx1, G4double pDx2,
93  G4double pAlp1,
94  G4double pDy2, G4double pDx3, G4double pDx4,
95  G4double pAlp2)
96  : G4CSGSolid(pName)
97 {
98  if ( pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 ||
99  pDx2 <= 0 || pDy2 <= 0 || pDx3 <= 0 || pDx4 <= 0 )
100  {
101  std::ostringstream message;
102  message << "Invalid length parameters for Solid: " << GetName() << G4endl
103  << " X - "
104  << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
105  << " Y - " << pDy1 << ", " << pDy2 << G4endl
106  << " Z - " << pDz;
107  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
108  FatalException, message);
109  }
110 
111  fDz=pDz;
112  fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
113  fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
114 
115  fDy1=pDy1;
116  fDx1=pDx1;
117  fDx2=pDx2;
118  fTalpha1=std::tan(pAlp1);
119 
120  fDy2=pDy2;
121  fDx3=pDx3;
122  fDx4=pDx4;
123  fTalpha2=std::tan(pAlp2);
124 
125  MakePlanes();
126 }
127 
129 //
130 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
131 // which are its vertices. Checking of planarity with preparation of
132 // fPlanes[] and than calculation of other members
133 
134 G4Trap::G4Trap( const G4String& pName,
135  const G4ThreeVector pt[8] )
136  : G4CSGSolid(pName)
137 {
138  G4bool good;
139 
140  // Start with check of centering - the center of gravity trap line
141  // should cross the origin of frame
142 
143  if (!( pt[0].z() < 0
144  && pt[0].z() == pt[1].z() && pt[0].z() == pt[2].z()
145  && pt[0].z() == pt[3].z()
146  && pt[4].z() > 0
147  && pt[4].z() == pt[5].z() && pt[4].z() == pt[6].z()
148  && pt[4].z() == pt[7].z()
149  && std::fabs( pt[0].z() + pt[4].z() ) < kCarTolerance
150  && pt[0].y() == pt[1].y() && pt[2].y() == pt[3].y()
151  && pt[4].y() == pt[5].y() && pt[6].y() == pt[7].y()
152  && std::fabs( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) < kCarTolerance
153  && std::fabs( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() +
154  pt[2].x() + pt[3].x() + pt[6].x() + pt[7].x() ) < kCarTolerance ) )
155  {
156  std::ostringstream message;
157  message << "Invalid vertice coordinates for Solid: " << GetName();
158  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
159  FatalException, message);
160  }
161 
162  // Bottom side with normal approx. -Y
163 
164  good = MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
165 
166  if (!good)
167  {
168  DumpInfo();
169  G4Exception("G4Trap::G4Trap()", "GeomSolids0002", FatalException,
170  "Face at ~-Y not planar.");
171  }
172 
173  // Top side with normal approx. +Y
174 
175  good = MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
176 
177  if (!good)
178  {
179  std::ostringstream message;
180  message << "Face at ~+Y not planar for Solid: " << GetName();
181  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
182  FatalException, message);
183  }
184 
185  // Front side with normal approx. -X
186 
187  good = MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
188 
189  if (!good)
190  {
191  std::ostringstream message;
192  message << "Face at ~-X not planar for Solid: " << GetName();
193  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
194  FatalException, message);
195  }
196 
197  // Back side iwth normal approx. +X
198 
199  good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
200  if (!good)
201  {
202  std::ostringstream message;
203  message << "Face at ~+X not planar for Solid: " << GetName();
204  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
205  FatalException, message);
206  }
207  fDz = (pt[7]).z() ;
208 
209  fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5;
210  fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5;
211  fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5;
212  fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1;
213 
214  fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5;
215  fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5;
216  fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5;
217  fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2;
218 
219  fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz;
220  fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;
221 }
222 
224 //
225 // Constructor for Right Angular Wedge from STEP
226 
227 G4Trap::G4Trap( const G4String& pName,
228  G4double pZ,
229  G4double pY,
230  G4double pX, G4double pLTX )
231  : G4CSGSolid(pName)
232 {
233  G4bool good;
234 
235  if ( pZ<=0 || pY<=0 || pX<=0 || pLTX<=0 || pLTX>pX )
236  {
237  std::ostringstream message;
238  message << "Invalid length parameters for Solid: " << GetName();
239  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
240  FatalException, message);
241  }
242 
243  fDz = 0.5*pZ ;
244  fTthetaCphi = 0 ;
245  fTthetaSphi = 0 ;
246 
247  fDy1 = 0.5*pY;
248  fDx1 = 0.5*pX ;
249  fDx2 = 0.5*pLTX;
250  fTalpha1 = 0.5*(pLTX - pX)/pY;
251 
252  fDy2 = fDy1 ;
253  fDx3 = fDx1;
254  fDx4 = fDx2 ;
255  fTalpha2 = fTalpha1 ;
256 
257  G4ThreeVector pt[8] ;
258 
260  -fDz*fTthetaSphi-fDy1,-fDz);
262  -fDz*fTthetaSphi-fDy1,-fDz);
264  -fDz*fTthetaSphi+fDy1,-fDz);
266  -fDz*fTthetaSphi+fDy1,-fDz);
268  +fDz*fTthetaSphi-fDy2,+fDz);
270  +fDz*fTthetaSphi-fDy2,+fDz);
272  +fDz*fTthetaSphi+fDy2,+fDz);
274  +fDz*fTthetaSphi+fDy2,+fDz);
275 
276  // Bottom side with normal approx. -Y
277  //
278  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
279  if (!good)
280  {
281  std::ostringstream message;
282  message << "Face at ~-Y not planar for Solid: " << GetName();
283  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
284  FatalException, message);
285  }
286 
287  // Top side with normal approx. +Y
288  //
289  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
290  if (!good)
291  {
292  std::ostringstream message;
293  message << "Face at ~+Y not planar for Solid: " << GetName();
294  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
295  FatalException, message);
296  }
297 
298  // Front side with normal approx. -X
299  //
300  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
301  if (!good)
302  {
303  std::ostringstream message;
304  message << "Face at ~-X not planar for Solid: " << GetName();
305  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
306  FatalException, message);
307  }
308 
309  // Back side iwth normal approx. +X
310  //
311  good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
312  if (!good)
313  {
314  std::ostringstream message;
315  message << "Face at ~+X not planar for Solid: " << GetName();
316  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
317  FatalException, message);
318  }
319 }
320 
322 //
323 // Constructor for G4Trd
324 
325 G4Trap::G4Trap( const G4String& pName,
326  G4double pDx1, G4double pDx2,
327  G4double pDy1, G4double pDy2,
328  G4double pDz )
329  : G4CSGSolid(pName)
330 {
331  G4bool good;
332 
333  if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 )
334  {
335  std::ostringstream message;
336  message << "Invalid length parameters for Solid: " << GetName();
337  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
338  FatalException, message);
339  }
340 
341  fDz = pDz;
342  fTthetaCphi = 0 ;
343  fTthetaSphi = 0 ;
344 
345  fDy1 = pDy1 ;
346  fDx1 = pDx1 ;
347  fDx2 = pDx1 ;
348  fTalpha1 = 0 ;
349 
350  fDy2 = pDy2 ;
351  fDx3 = pDx2 ;
352  fDx4 = pDx2 ;
353  fTalpha2 = 0 ;
354 
355  G4ThreeVector pt[8] ;
356 
358  -fDz*fTthetaSphi-fDy1,-fDz);
360  -fDz*fTthetaSphi-fDy1,-fDz);
362  -fDz*fTthetaSphi+fDy1,-fDz);
364  -fDz*fTthetaSphi+fDy1,-fDz);
366  +fDz*fTthetaSphi-fDy2,+fDz);
368  +fDz*fTthetaSphi-fDy2,+fDz);
370  +fDz*fTthetaSphi+fDy2,+fDz);
372  +fDz*fTthetaSphi+fDy2,+fDz);
373 
374  // Bottom side with normal approx. -Y
375  //
376  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
377  if (!good)
378  {
379  std::ostringstream message;
380  message << "Face at ~-Y not planar for Solid: " << GetName();
381  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
382  FatalException, message);
383  }
384 
385  // Top side with normal approx. +Y
386  //
387  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
388  if (!good)
389  {
390  std::ostringstream message;
391  message << "Face at ~+Y not planar for Solid: " << GetName();
392  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
393  FatalException, message);
394  }
395 
396  // Front side with normal approx. -X
397  //
398  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
399  if (!good)
400  {
401  std::ostringstream message;
402  message << "Face at ~-X not planar for Solid: " << GetName();
403  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
404  FatalException, message);
405  }
406 
407  // Back side iwth normal approx. +X
408  //
409  good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
410  if (!good)
411  {
412  std::ostringstream message;
413  message << "Face at ~+X not planar for Solid: " << GetName();
414  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
415  FatalException, message);
416  }
417 }
418 
420 //
421 // Constructor for G4Para
422 
423 G4Trap::G4Trap( const G4String& pName,
424  G4double pDx, G4double pDy,
425  G4double pDz,
426  G4double pAlpha,
427  G4double pTheta, G4double pPhi )
428  : G4CSGSolid(pName)
429 {
430  G4bool good;
431 
432  if ( pDz<=0 || pDy<=0 || pDx<=0 )
433  {
434  std::ostringstream message;
435  message << "Invalid length parameters for Solid: " << GetName();
436  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
437  FatalException, message);
438  }
439 
440  fDz = pDz ;
441  fTthetaCphi = std::tan(pTheta)*std::cos(pPhi) ;
442  fTthetaSphi = std::tan(pTheta)*std::sin(pPhi) ;
443 
444  fDy1 = pDy ;
445  fDx1 = pDx ;
446  fDx2 = pDx ;
447  fTalpha1 = std::tan(pAlpha) ;
448 
449  fDy2 = pDy ;
450  fDx3 = pDx ;
451  fDx4 = pDx ;
452  fTalpha2 = fTalpha1 ;
453 
454  G4ThreeVector pt[8] ;
455 
457  -fDz*fTthetaSphi-fDy1,-fDz);
459  -fDz*fTthetaSphi-fDy1,-fDz);
461  -fDz*fTthetaSphi+fDy1,-fDz);
463  -fDz*fTthetaSphi+fDy1,-fDz);
465  +fDz*fTthetaSphi-fDy2,+fDz);
467  +fDz*fTthetaSphi-fDy2,+fDz);
469  +fDz*fTthetaSphi+fDy2,+fDz);
471  +fDz*fTthetaSphi+fDy2,+fDz);
472 
473  // Bottom side with normal approx. -Y
474  //
475  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
476  if (!good)
477  {
478  std::ostringstream message;
479  message << "Face at ~-Y not planar for Solid: " << GetName();
480  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
481  FatalException, message);
482  }
483 
484  // Top side with normal approx. +Y
485  //
486  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
487  if (!good)
488  {
489  std::ostringstream message;
490  message << "Face at ~+Y not planar for Solid: " << GetName();
491  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
492  FatalException, message);
493  }
494 
495  // Front side with normal approx. -X
496  //
497  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
498  if (!good)
499  {
500  std::ostringstream message;
501  message << "Face at ~-X not planar for Solid: " << GetName();
502  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
503  FatalException, message);
504  }
505 
506  // Back side iwth normal approx. +X
507  //
508  good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
509  if (!good)
510  {
511  std::ostringstream message;
512  message << "Face at ~+X not planar for Solid: " << GetName();
513  G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
514  FatalException, message);
515  }
516 }
517 
519 //
520 // Nominal constructor for G4Trap whose parameters are to be set by
521 // a G4VParamaterisation later. Check and set half-widths as well as
522 // angles: final check of coplanarity
523 
524 G4Trap::G4Trap( const G4String& pName )
525  : G4CSGSolid (pName), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
526  fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
527  fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
528 {
529  MakePlanes();
530 }
531 
533 //
534 // Fake default constructor - sets only member data and allocates memory
535 // for usage restricted to object persistency.
536 //
537 G4Trap::G4Trap( __void__& a )
538  : G4CSGSolid(a), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
539  fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
540  fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
541 {
542  MakePlanes();
543 }
544 
546 //
547 // Destructor
548 
550 {
551 }
552 
554 //
555 // Copy constructor
556 
558  : G4CSGSolid(rhs), fDz(rhs.fDz),
559  fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
560  fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
561  fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
562 {
563  for (size_t i=0; i<4; ++i)
564  {
565  fPlanes[i].a = rhs.fPlanes[i].a;
566  fPlanes[i].b = rhs.fPlanes[i].b;
567  fPlanes[i].c = rhs.fPlanes[i].c;
568  fPlanes[i].d = rhs.fPlanes[i].d;
569  }
570 }
571 
573 //
574 // Assignment operator
575 
577 {
578  // Check assignment to self
579  //
580  if (this == &rhs) { return *this; }
581 
582  // Copy base class data
583  //
585 
586  // Copy data
587  //
588  fDz = rhs.fDz;
590  fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
591  fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
592  for (size_t i=0; i<4; ++i)
593  {
594  fPlanes[i].a = rhs.fPlanes[i].a;
595  fPlanes[i].b = rhs.fPlanes[i].b;
596  fPlanes[i].c = rhs.fPlanes[i].c;
597  fPlanes[i].d = rhs.fPlanes[i].d;
598  }
599 
600  return *this;
601 }
602 
604 //
605 // Set all parameters, as for constructor - check and set half-widths
606 // as well as angles: final check of coplanarity
607 
609  G4double pTheta,
610  G4double pPhi,
611  G4double pDy1,
612  G4double pDx1,
613  G4double pDx2,
614  G4double pAlp1,
615  G4double pDy2,
616  G4double pDx3,
617  G4double pDx4,
618  G4double pAlp2 )
619 {
620  if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 || pDx3<=0 || pDx4<=0 )
621  {
622  std::ostringstream message;
623  message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
624  << " X - "
625  << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
626  << " Y - " << pDy1 << ", " << pDy2 << G4endl
627  << " Z - " << pDz;
628  G4Exception("G4Trap::SetAllParameters()", "GeomSolids0002",
629  FatalException, message);
630  }
631  fCubicVolume= 0.;
632  fSurfaceArea= 0.;
633  fRebuildPolyhedron = true;
634  fDz=pDz;
635  fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
636  fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
637 
638  fDy1=pDy1;
639  fDx1=pDx1;
640  fDx2=pDx2;
641  fTalpha1=std::tan(pAlp1);
642 
643  fDy2=pDy2;
644  fDx3=pDx3;
645  fDx4=pDx4;
646  fTalpha2=std::tan(pAlp2);
647 
648  MakePlanes();
649 }
650 
652 //
653 // Checking of coplanarity
654 
656 {
657  G4bool good = true;
658 
659  G4ThreeVector pt[8] ;
660 
662  -fDz*fTthetaSphi-fDy1,-fDz);
664  -fDz*fTthetaSphi-fDy1,-fDz);
666  -fDz*fTthetaSphi+fDy1,-fDz);
668  -fDz*fTthetaSphi+fDy1,-fDz);
670  +fDz*fTthetaSphi-fDy2,+fDz);
672  +fDz*fTthetaSphi-fDy2,+fDz);
674  +fDz*fTthetaSphi+fDy2,+fDz);
676  +fDz*fTthetaSphi+fDy2,+fDz);
677 
678  // Bottom side with normal approx. -Y
679  //
680  good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]) ;
681  if (!good)
682  {
683  std::ostringstream message;
684  message << "Face at ~-Y not planar for Solid: " << GetName();
685  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
686  FatalException, message);
687  }
688 
689  // Top side with normal approx. +Y
690  //
691  good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
692  if (!good)
693  {
694  std::ostringstream message;
695  message << "Face at ~+Y not planar for Solid: " << GetName();
696  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
697  FatalException, message);
698  }
699 
700  // Front side with normal approx. -X
701  //
702  good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
703  if (!good)
704  {
705  std::ostringstream message;
706  message << "Face at ~-X not planar for Solid: " << GetName();
707  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
708  FatalException, message);
709  }
710 
711  // Back side iwth normal approx. +X
712  //
713  good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
714  if ( !good )
715  {
716  std::ostringstream message;
717  message << "Face at ~+X not planar for Solid: " << GetName();
718  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
719  FatalException, message);
720  }
721 
722  return good;
723 }
724 
726 //
727 // Calculate the coef's of the plane p1->p2->p3->p4->p1
728 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed from
729 // infront of the plane (i.e. from normal direction).
730 //
731 // Return true if the ThreeVectors are coplanar + set coef;s
732 // false if ThreeVectors are not coplanar
733 
735  const G4ThreeVector& p2,
736  const G4ThreeVector& p3,
737  const G4ThreeVector& p4,
738  TrapSidePlane& plane )
739 {
740  G4double a, b, c, sd;
741  G4ThreeVector v12, v13, v14, Vcross;
742 
743  G4bool good;
744 
745  v12 = p2 - p1;
746  v13 = p3 - p1;
747  v14 = p4 - p1;
748  Vcross = v12.cross(v13);
749 
750  if (std::fabs(Vcross.dot(v14)/(Vcross.mag()*v14.mag())) > kCoplanar_Tolerance)
751  {
752  good = false;
753  }
754  else
755  {
756  // a,b,c correspond to the x/y/z components of the
757  // normal vector to the plane
758 
759  // a = (p2.y()-p1.y())*(p1.z()+p2.z())+(p3.y()-p2.y())*(p2.z()+p3.z());
760  // a += (p4.y()-p3.y())*(p3.z()+p4.z())+(p1.y()-p4.y())*(p4.z()+p1.z()); // ?
761  // b = (p2.z()-p1.z())*(p1.x()+p2.x())+(p3.z()-p2.z())*(p2.x()+p3.x());
762  // b += (p4.z()-p3.z())*(p3.x()+p4.x())+(p1.z()-p4.z())*(p4.x()+p1.x()); // ?
763  // c = (p2.x()-p1.x())*(p1.y()+p2.y())+(p3.x()-p2.x())*(p2.y()+p3.y());
764  // c += (p4.x()-p3.x())*(p3.y()+p4.y())+(p1.x()-p4.x())*(p4.y()+p1.y()); // ?
765 
766  // Let create diagonals 4-2 and 3-1 than (4-2)x(3-1) provides
767  // vector perpendicular to the plane directed to outside !!!
768  // and a,b,c, = f(1,2,3,4) external relative to trap normal
769 
770  a = +(p4.y() - p2.y())*(p3.z() - p1.z())
771  - (p3.y() - p1.y())*(p4.z() - p2.z());
772 
773  b = -(p4.x() - p2.x())*(p3.z() - p1.z())
774  + (p3.x() - p1.x())*(p4.z() - p2.z());
775 
776  c = +(p4.x() - p2.x())*(p3.y() - p1.y())
777  - (p3.x() - p1.x())*(p4.y() - p2.y());
778 
779  sd = std::sqrt( a*a + b*b + c*c ); // so now vector plane.(a,b,c) is unit
780 
781  if( sd > 0 )
782  {
783  plane.a = a/sd;
784  plane.b = b/sd;
785  plane.c = c/sd;
786  }
787  else
788  {
789  std::ostringstream message;
790  message << "Invalid parameters: norm.mod() <= 0, for Solid: "
791  << GetName();
792  G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
793  FatalException, message) ;
794  }
795  // Calculate D: p1 in in plane so D=-n.p1.Vect()
796 
797  plane.d = -( plane.a*p1.x() + plane.b*p1.y() + plane.c*p1.z() );
798 
799  good = true;
800  }
801  return good;
802 }
803 
805 //
806 // Dispatch to parameterisation for replication mechanism dimension
807 // computation & modification.
808 
810  const G4int n,
811  const G4VPhysicalVolume* pRep )
812 {
813  p->ComputeDimensions(*this,n,pRep);
814 }
815 
817 //
818 // Get bounding box
819 
821 {
822  G4double dz = GetZHalfLength();
823  G4double dx1 = GetXHalfLength1();
824  G4double dx2 = GetXHalfLength2();
825  G4double dx3 = GetXHalfLength3();
826  G4double dx4 = GetXHalfLength4();
827  G4double dy1 = GetYHalfLength1();
828  G4double dy2 = GetYHalfLength2();
829 
830  G4double x0 = dz*fTthetaCphi;
831  G4double x1 = dy1*GetTanAlpha1();
832  G4double x2 = dy2*GetTanAlpha2();
833  G4double xmin =
834  std::min(
835  std::min(
836  std::min(-x0-x1-dx1,-x0+x1-dx2),x0-x2-dx3),x0+x2-dx4);
837  G4double xmax =
838  std::max(
839  std::max(
840  std::max(-x0-x1+dx1,-x0+x1+dx2),x0-x2+dx3),x0+x2+dx4);
841 
842  G4double y0 = dz*fTthetaSphi;
843  G4double ymin = std::min(-y0-dy1,y0-dy2);
844  G4double ymax = std::max(-y0+dy1,y0+dy2);
845 
846  pMin.set(xmin,ymin,-dz);
847  pMax.set(xmax,ymax, dz);
848 
849  // Check correctness of the bounding box
850  //
851  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
852  {
853  std::ostringstream message;
854  message << "Bad bounding box (min >= max) for solid: "
855  << GetName() << " !"
856  << "\npMin = " << pMin
857  << "\npMax = " << pMax;
858  G4Exception("G4Trap::Extent()", "GeomMgt0001", JustWarning, message);
859  DumpInfo();
860  }
861 }
862 
864 //
865 // Calculate extent under transform and specified limit
866 
868  const G4VoxelLimits& pVoxelLimit,
869  const G4AffineTransform& pTransform,
870  G4double& pMin, G4double& pMax) const
871 {
872  G4ThreeVector bmin, bmax;
873  G4bool exist;
874 
875  // Check bounding box (bbox)
876  //
877  Extent(bmin,bmax);
878  G4BoundingEnvelope bbox(bmin,bmax);
879 #ifdef G4BBOX_EXTENT
880  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
881 #endif
882  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
883  {
884  return exist = (pMin < pMax) ? true : false;
885  }
886 
887  // Set bounding envelope (benv) and calculate extent
888  //
889  G4double dz = GetZHalfLength();
890  G4double dx1 = GetXHalfLength1();
891  G4double dx2 = GetXHalfLength2();
892  G4double dx3 = GetXHalfLength3();
893  G4double dx4 = GetXHalfLength4();
894  G4double dy1 = GetYHalfLength1();
895  G4double dy2 = GetYHalfLength2();
896 
897  G4double x0 = dz*fTthetaCphi;
898  G4double x1 = dy1*GetTanAlpha1();
899  G4double x2 = dy2*GetTanAlpha2();
900  G4double y0 = dz*fTthetaSphi;
901 
902  G4ThreeVectorList baseA(4), baseB(4);
903  baseA[0].set(-x0-x1-dx1,-y0-dy1,-dz);
904  baseA[1].set(-x0-x1+dx1,-y0-dy1,-dz);
905  baseA[2].set(-x0+x1+dx2,-y0+dy1,-dz);
906  baseA[3].set(-x0+x1-dx2,-y0+dy1,-dz);
907 
908  baseB[0].set( x0-x2-dx3, y0-dy2, dz);
909  baseB[1].set( x0-x2+dx3, y0-dy2, dz);
910  baseB[2].set( x0+x2+dx4, y0+dy2, dz);
911  baseB[3].set( x0+x2-dx4, y0+dy2, dz);
912 
913  std::vector<const G4ThreeVectorList *> polygons(2);
914  polygons[0] = &baseA;
915  polygons[1] = &baseB;
916 
917  G4BoundingEnvelope benv(bmin,bmax,polygons);
918  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
919  return exist;
920 }
921 
923 //
924 // Return whether point inside/outside/on surface, using tolerance
925 
927 {
928  EInside in;
929  G4double Dist;
930  G4int i;
931  if ( std::fabs(p.z()) <= fDz-kCarTolerance*0.5)
932  {
933  in = kInside;
934 
935  for ( i = 0;i < 4;i++ )
936  {
937  Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
938  +fPlanes[i].c*p.z() + fPlanes[i].d;
939 
940  if (Dist > kCarTolerance*0.5) return in = kOutside;
941  else if (Dist > -kCarTolerance*0.5) in = kSurface;
942 
943  }
944  }
945  else if (std::fabs(p.z()) <= fDz+kCarTolerance*0.5)
946  {
947  in = kSurface;
948 
949  for ( i = 0; i < 4; i++ )
950  {
951  Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
952  +fPlanes[i].c*p.z() + fPlanes[i].d;
953 
954  if (Dist > kCarTolerance*0.5) return in = kOutside;
955  }
956  }
957  else in = kOutside;
958 
959  return in;
960 }
961 
963 //
964 // Calculate side nearest to p, and return normal
965 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
966 
968 {
969  G4int i, noSurfaces = 0;
970  G4double dist, distz, distx, disty, distmx, distmy, safe = kInfinity;
971  G4double delta = 0.5*kCarTolerance;
972  G4ThreeVector norm, sumnorm(0.,0.,0.);
973 
974  for (i = 0; i < 4; i++)
975  {
976  dist = std::fabs(fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
977  + fPlanes[i].c*p.z() + fPlanes[i].d);
978  if ( dist < safe )
979  {
980  safe = dist;
981  }
982  }
983  distz = std::fabs( std::fabs( p.z() ) - fDz );
984 
985  distmy = std::fabs( fPlanes[0].a*p.x() + fPlanes[0].b*p.y()
986  + fPlanes[0].c*p.z() + fPlanes[0].d );
987 
988  disty = std::fabs( fPlanes[1].a*p.x() + fPlanes[1].b*p.y()
989  + fPlanes[1].c*p.z() + fPlanes[1].d );
990 
991  distmx = std::fabs( fPlanes[2].a*p.x() + fPlanes[2].b*p.y()
992  + fPlanes[2].c*p.z() + fPlanes[2].d );
993 
994  distx = std::fabs( fPlanes[3].a*p.x() + fPlanes[3].b*p.y()
995  + fPlanes[3].c*p.z() + fPlanes[3].d );
996 
997  G4ThreeVector nX = G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
998  G4ThreeVector nmX = G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
999  G4ThreeVector nY = G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1000  G4ThreeVector nmY = G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1001  G4ThreeVector nZ = G4ThreeVector(0.,0.,1.0);
1002 
1003  if (distx <= delta)
1004  {
1005  noSurfaces ++;
1006  sumnorm += nX;
1007  }
1008  if (distmx <= delta)
1009  {
1010  noSurfaces ++;
1011  sumnorm += nmX;
1012  }
1013  if (disty <= delta)
1014  {
1015  noSurfaces ++;
1016  sumnorm += nY;
1017  }
1018  if (distmy <= delta)
1019  {
1020  noSurfaces ++;
1021  sumnorm += nmY;
1022  }
1023  if (distz <= delta)
1024  {
1025  noSurfaces ++;
1026  if ( p.z() >= 0.) sumnorm += nZ;
1027  else sumnorm -= nZ;
1028  }
1029  if ( noSurfaces == 0 )
1030  {
1031 #ifdef G4CSGDEBUG
1032  G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
1033  JustWarning, "Point p is not on surface !?" );
1034 #endif
1035  norm = ApproxSurfaceNormal(p);
1036  }
1037  else if ( noSurfaces == 1 ) norm = sumnorm;
1038  else norm = sumnorm.unit();
1039  return norm;
1040 }
1041 
1043 //
1044 // Algorithm for SurfaceNormal() following the original specification
1045 // for points not on the surface
1046 
1048 {
1049  G4double safe=kInfinity,Dist,safez;
1050  G4int i,imin=0;
1051  for (i=0;i<4;i++)
1052  {
1053  Dist=std::fabs(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1054  +fPlanes[i].c*p.z()+fPlanes[i].d);
1055  if (Dist<safe)
1056  {
1057  safe=Dist;
1058  imin=i;
1059  }
1060  }
1061  safez=std::fabs(std::fabs(p.z())-fDz);
1062  if (safe<safez)
1063  {
1064  return G4ThreeVector(fPlanes[imin].a,fPlanes[imin].b,fPlanes[imin].c);
1065  }
1066  else
1067  {
1068  if (p.z()>0)
1069  {
1070  return G4ThreeVector(0,0,1);
1071  }
1072  else
1073  {
1074  return G4ThreeVector(0,0,-1);
1075  }
1076  }
1077 }
1078 
1080 //
1081 // Calculate distance to shape from outside - return kInfinity if no intersection
1082 //
1083 // ALGORITHM:
1084 // For each component, calculate pair of minimum and maximum intersection
1085 // values for which the particle is in the extent of the shape
1086 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
1087 
1089  const G4ThreeVector& v ) const
1090 {
1091 
1092  G4double snxt; // snxt = default return value
1093  G4double max,smax,smin;
1094  G4double pdist,Comp,vdist;
1095  G4int i;
1096  //
1097  // Z Intersection range
1098  //
1099  if ( v.z() > 0 )
1100  {
1101  max = fDz - p.z() ;
1102  if (max > 0.5*kCarTolerance)
1103  {
1104  smax = max/v.z();
1105  smin = (-fDz-p.z())/v.z();
1106  }
1107  else
1108  {
1109  return snxt=kInfinity;
1110  }
1111  }
1112  else if (v.z() < 0 )
1113  {
1114  max = - fDz - p.z() ;
1115  if (max < -0.5*kCarTolerance )
1116  {
1117  smax=max/v.z();
1118  smin=(fDz-p.z())/v.z();
1119  }
1120  else
1121  {
1122  return snxt=kInfinity;
1123  }
1124  }
1125  else
1126  {
1127  if (std::fabs(p.z())<fDz - 0.5*kCarTolerance) // Inside was <=fDz
1128  {
1129  smin=0;
1130  smax=kInfinity;
1131  }
1132  else
1133  {
1134  return snxt=kInfinity;
1135  }
1136  }
1137 
1138  for (i=0;i<4;i++)
1139  {
1140  pdist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1141  +fPlanes[i].c*p.z()+fPlanes[i].d;
1142  Comp=fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
1143  if ( pdist >= -0.5*kCarTolerance ) // was >0
1144  {
1145  //
1146  // Outside the plane -> this is an extent entry distance
1147  //
1148  if (Comp >= 0) // was >0
1149  {
1150  return snxt=kInfinity ;
1151  }
1152  else
1153  {
1154  vdist=-pdist/Comp;
1155  if (vdist>smin)
1156  {
1157  if (vdist<smax)
1158  {
1159  smin = vdist;
1160  }
1161  else
1162  {
1163  return snxt=kInfinity;
1164  }
1165  }
1166  }
1167  }
1168  else
1169  {
1170  //
1171  // Inside the plane -> couble be an extent exit distance (smax)
1172  //
1173  if (Comp>0) // Will leave extent
1174  {
1175  vdist=-pdist/Comp;
1176  if (vdist<smax)
1177  {
1178  if (vdist>smin)
1179  {
1180  smax=vdist;
1181  }
1182  else
1183  {
1184  return snxt=kInfinity;
1185  }
1186  }
1187  }
1188  }
1189  }
1190  //
1191  // Checks in non z plane intersections ensure smin<smax
1192  //
1193  if (smin >=0 )
1194  {
1195  snxt = smin ;
1196  }
1197  else
1198  {
1199  snxt = 0 ;
1200  }
1201  return snxt;
1202 }
1203 
1205 //
1206 // Calculate exact shortest distance to any boundary from outside
1207 // This is the best fast estimation of the shortest distance to trap
1208 // - Returns 0 is ThreeVector inside
1209 
1211 {
1212  G4double safe=0.0,Dist;
1213  G4int i;
1214  safe=std::fabs(p.z())-fDz;
1215  for (i=0;i<4;i++)
1216  {
1217  Dist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1218  +fPlanes[i].c*p.z()+fPlanes[i].d;
1219  if (Dist > safe) safe=Dist;
1220  }
1221  if (safe<0) safe=0;
1222  return safe;
1223 }
1224 
1226 //
1227 // Calculate distance to surface of shape from inside
1228 // Calculate distance to x/y/z planes - smallest is exiting distance
1229 
1231  const G4bool calcNorm,
1232  G4bool *validNorm, G4ThreeVector *n) const
1233 {
1234  Eside side = kUndef;
1235  G4double snxt; // snxt = return value
1236  G4double pdist,Comp,vdist,max;
1237  //
1238  // Z Intersections
1239  //
1240  if (v.z()>0)
1241  {
1242  max=fDz-p.z();
1243  if (max>kCarTolerance/2)
1244  {
1245  snxt=max/v.z();
1246  side=kPZ;
1247  }
1248  else
1249  {
1250  if (calcNorm)
1251  {
1252  *validNorm=true;
1253  *n=G4ThreeVector(0,0,1);
1254  }
1255  return snxt=0;
1256  }
1257  }
1258  else if (v.z()<0)
1259  {
1260  max=-fDz-p.z();
1261  if (max<-kCarTolerance/2)
1262  {
1263  snxt=max/v.z();
1264  side=kMZ;
1265  }
1266  else
1267  {
1268  if (calcNorm)
1269  {
1270  *validNorm=true;
1271  *n=G4ThreeVector(0,0,-1);
1272  }
1273  return snxt=0;
1274  }
1275  }
1276  else
1277  {
1278  snxt=kInfinity;
1279  }
1280 
1281  //
1282  // Intersections with planes[0] (expanded because of setting enum)
1283  //
1284  pdist=fPlanes[0].a*p.x()+fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
1285  Comp=fPlanes[0].a*v.x()+fPlanes[0].b*v.y()+fPlanes[0].c*v.z();
1286  if (pdist>0)
1287  {
1288  // Outside the plane
1289  if (Comp>0)
1290  {
1291  // Leaving immediately
1292  if (calcNorm)
1293  {
1294  *validNorm=true;
1295  *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1296  }
1297  return snxt=0;
1298  }
1299  }
1300  else if (pdist<-kCarTolerance/2)
1301  {
1302  // Inside the plane
1303  if (Comp>0)
1304  {
1305  // Will leave extent
1306  vdist=-pdist/Comp;
1307  if (vdist<snxt)
1308  {
1309  snxt=vdist;
1310  side=ks0;
1311  }
1312  }
1313  }
1314  else
1315  {
1316  // On surface
1317  if (Comp>0)
1318  {
1319  if (calcNorm)
1320  {
1321  *validNorm=true;
1322  *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1323  }
1324  return snxt=0;
1325  }
1326  }
1327 
1328  //
1329  // Intersections with planes[1] (expanded because of setting enum)
1330  //
1331  pdist=fPlanes[1].a*p.x()+fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
1332  Comp=fPlanes[1].a*v.x()+fPlanes[1].b*v.y()+fPlanes[1].c*v.z();
1333  if (pdist>0)
1334  {
1335  // Outside the plane
1336  if (Comp>0)
1337  {
1338  // Leaving immediately
1339  if (calcNorm)
1340  {
1341  *validNorm=true;
1342  *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1343  }
1344  return snxt=0;
1345  }
1346  }
1347  else if (pdist<-kCarTolerance/2)
1348  {
1349  // Inside the plane
1350  if (Comp>0)
1351  {
1352  // Will leave extent
1353  vdist=-pdist/Comp;
1354  if (vdist<snxt)
1355  {
1356  snxt=vdist;
1357  side=ks1;
1358  }
1359  }
1360  }
1361  else
1362  {
1363  // On surface
1364  if (Comp>0)
1365  {
1366  if (calcNorm)
1367  {
1368  *validNorm=true;
1369  *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1370  }
1371  return snxt=0;
1372  }
1373  }
1374 
1375  //
1376  // Intersections with planes[2] (expanded because of setting enum)
1377  //
1378  pdist=fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d;
1379  Comp=fPlanes[2].a*v.x()+fPlanes[2].b*v.y()+fPlanes[2].c*v.z();
1380  if (pdist>0)
1381  {
1382  // Outside the plane
1383  if (Comp>0)
1384  {
1385  // Leaving immediately
1386  if (calcNorm)
1387  {
1388  *validNorm=true;
1389  *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1390  }
1391  return snxt=0;
1392  }
1393  }
1394  else if (pdist<-kCarTolerance/2)
1395  {
1396  // Inside the plane
1397  if (Comp>0)
1398  {
1399  // Will leave extent
1400  vdist=-pdist/Comp;
1401  if (vdist<snxt)
1402  {
1403  snxt=vdist;
1404  side=ks2;
1405  }
1406  }
1407  }
1408  else
1409  {
1410  // On surface
1411  if (Comp>0)
1412  {
1413  if (calcNorm)
1414  {
1415  *validNorm=true;
1416  *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1417  }
1418  return snxt=0;
1419  }
1420  }
1421 
1422  //
1423  // Intersections with planes[3] (expanded because of setting enum)
1424  //
1425  pdist=fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d;
1426  Comp=fPlanes[3].a*v.x()+fPlanes[3].b*v.y()+fPlanes[3].c*v.z();
1427  if (pdist>0)
1428  {
1429  // Outside the plane
1430  if (Comp>0)
1431  {
1432  // Leaving immediately
1433  if (calcNorm)
1434  {
1435  *validNorm=true;
1436  *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1437  }
1438  return snxt=0;
1439  }
1440  }
1441  else if (pdist<-kCarTolerance/2)
1442  {
1443  // Inside the plane
1444  if (Comp>0)
1445  {
1446  // Will leave extent
1447  vdist=-pdist/Comp;
1448  if (vdist<snxt)
1449  {
1450  snxt=vdist;
1451  side=ks3;
1452  }
1453  }
1454  }
1455  else
1456  {
1457  // On surface
1458  if (Comp>0)
1459  {
1460  if (calcNorm)
1461  {
1462  *validNorm=true;
1463  *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1464  }
1465  return snxt=0;
1466  }
1467  }
1468 
1469  // set normal
1470  if (calcNorm)
1471  {
1472  *validNorm=true;
1473  switch(side)
1474  {
1475  case ks0:
1476  *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
1477  break;
1478  case ks1:
1479  *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1480  break;
1481  case ks2:
1482  *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
1483  break;
1484  case ks3:
1485  *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
1486  break;
1487  case kMZ:
1488  *n=G4ThreeVector(0,0,-1);
1489  break;
1490  case kPZ:
1491  *n=G4ThreeVector(0,0,1);
1492  break;
1493  default:
1494  G4cout << G4endl;
1495  DumpInfo();
1496  std::ostringstream message;
1497  G4int oldprc = message.precision(16);
1498  message << "Undefined side for valid surface normal to solid."
1499  << G4endl
1500  << "Position:" << G4endl << G4endl
1501  << "p.x() = " << p.x()/mm << " mm" << G4endl
1502  << "p.y() = " << p.y()/mm << " mm" << G4endl
1503  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1504  << "Direction:" << G4endl << G4endl
1505  << "v.x() = " << v.x() << G4endl
1506  << "v.y() = " << v.y() << G4endl
1507  << "v.z() = " << v.z() << G4endl << G4endl
1508  << "Proposed distance :" << G4endl << G4endl
1509  << "snxt = " << snxt/mm << " mm" << G4endl;
1510  message.precision(oldprc);
1511  G4Exception("G4Trap::DistanceToOut(p,v,..)","GeomSolids1002",
1512  JustWarning, message);
1513  break;
1514  }
1515  }
1516  return snxt;
1517 }
1518 
1520 //
1521 // Calculate exact shortest distance to any boundary from inside
1522 // - Returns 0 is ThreeVector outside
1523 
1525 {
1526  G4double safe=0.0,Dist;
1527  G4int i;
1528 
1529 #ifdef G4CSGDEBUG
1530  if( Inside(p) == kOutside )
1531  {
1532  G4int oldprc = G4cout.precision(16) ;
1533  G4cout << G4endl ;
1534  DumpInfo();
1535  G4cout << "Position:" << G4endl << G4endl ;
1536  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1537  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1538  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1539  G4cout.precision(oldprc) ;
1540  G4Exception("G4Trap::DistanceToOut(p)",
1541  "GeomSolids1002", JustWarning, "Point p is outside !?" );
1542  }
1543 #endif
1544 
1545  safe=fDz-std::fabs(p.z());
1546  if (safe<0) safe=0;
1547  else
1548  {
1549  for (i=0;i<4;i++)
1550  {
1551  Dist=-(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1552  +fPlanes[i].c*p.z()+fPlanes[i].d);
1553  if (Dist<safe) safe=Dist;
1554  }
1555  if (safe<0) safe=0;
1556  }
1557  return safe;
1558 }
1559 
1561 //
1562 // GetEntityType
1563 
1565 {
1566  return G4String("G4Trap");
1567 }
1568 
1570 //
1571 // Make a clone of the object
1572 //
1574 {
1575  return new G4Trap(*this);
1576 }
1577 
1579 //
1580 // Stream object contents to an output stream
1581 
1582 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1583 {
1584  G4int oldprc = os.precision(16);
1585  os << "-----------------------------------------------------------\n"
1586  << " *** Dump for solid - " << GetName() << " ***\n"
1587  << " ===================================================\n"
1588  << " Solid type: G4Trap\n"
1589  << " Parameters: \n"
1590  << " half length Z: " << fDz/mm << " mm \n"
1591  << " half length Y of face -fDz: " << fDy1/mm << " mm \n"
1592  << " half length X of side -fDy1, face -fDz: " << fDx1/mm << " mm \n"
1593  << " half length X of side +fDy1, face -fDz: " << fDx2/mm << " mm \n"
1594  << " half length Y of face +fDz: " << fDy2/mm << " mm \n"
1595  << " half length X of side -fDy2, face +fDz: " << fDx3/mm << " mm \n"
1596  << " half length X of side +fDy2, face +fDz: " << fDx4/mm << " mm \n"
1597  << " std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree << " degrees \n"
1598  << " std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree << " degrees \n"
1599  << " std::tan(alpha), -fDz: " << fTalpha1/degree << " degrees \n"
1600  << " std::tan(alpha), +fDz: " << fTalpha2/degree << " degrees \n"
1601  << " trap side plane equations:\n"
1602  << " " << fPlanes[0].a << " X + " << fPlanes[0].b << " Y + "
1603  << fPlanes[0].c << " Z + " << fPlanes[0].d << " = 0\n"
1604  << " " << fPlanes[1].a << " X + " << fPlanes[1].b << " Y + "
1605  << fPlanes[1].c << " Z + " << fPlanes[1].d << " = 0\n"
1606  << " " << fPlanes[2].a << " X + " << fPlanes[2].b << " Y + "
1607  << fPlanes[2].c << " Z + " << fPlanes[2].d << " = 0\n"
1608  << " " << fPlanes[3].a << " X + " << fPlanes[3].b << " Y + "
1609  << fPlanes[3].c << " Z + " << fPlanes[3].d << " = 0\n"
1610  << "-----------------------------------------------------------\n";
1611  os.precision(oldprc);
1612 
1613  return os;
1614 }
1615 
1617 //
1618 // GetPointOnPlane
1619 //
1620 // Auxiliary method for Get Point on Surface
1621 
1624  G4double& area) const
1625 {
1626  G4double lambda1, lambda2, chose, aOne, aTwo;
1627  G4ThreeVector t, u, v, w, Area, normal;
1628 
1629  t = p1 - p0;
1630  u = p2 - p1;
1631  v = p3 - p2;
1632  w = p0 - p3;
1633 
1634  Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
1635  w.z()*v.x() - w.x()*v.z(),
1636  w.x()*v.y() - w.y()*v.x());
1637 
1638  aOne = 0.5*Area.mag();
1639 
1640  Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
1641  t.z()*u.x() - t.x()*u.z(),
1642  t.x()*u.y() - t.y()*u.x());
1643 
1644  aTwo = 0.5*Area.mag();
1645 
1646  area = aOne + aTwo;
1647 
1648  chose = G4RandFlat::shoot(0.,aOne+aTwo);
1649 
1650  if( (chose>=0.) && (chose < aOne) )
1651  {
1652  lambda1 = G4RandFlat::shoot(0.,1.);
1653  lambda2 = G4RandFlat::shoot(0.,lambda1);
1654  return (p2+lambda1*v+lambda2*w);
1655  }
1656 
1657  // else
1658 
1659  lambda1 = G4RandFlat::shoot(0.,1.);
1660  lambda2 = G4RandFlat::shoot(0.,lambda1);
1661 
1662  return (p0+lambda1*t+lambda2*u);
1663 }
1664 
1666 //
1667 // GetPointOnSurface
1668 
1670 {
1671  G4double aOne, aTwo, aThree, aFour, aFive, aSix, chose;
1672  G4ThreeVector One, Two, Three, Four, Five, Six, test;
1673  G4ThreeVector pt[8];
1674 
1676  -fDz*fTthetaSphi-fDy1,-fDz);
1678  -fDz*fTthetaSphi-fDy1,-fDz);
1680  -fDz*fTthetaSphi+fDy1,-fDz);
1682  -fDz*fTthetaSphi+fDy1,-fDz);
1684  +fDz*fTthetaSphi-fDy2,+fDz);
1686  +fDz*fTthetaSphi-fDy2,+fDz);
1688  +fDz*fTthetaSphi+fDy2,+fDz);
1690  +fDz*fTthetaSphi+fDy2,+fDz);
1691 
1692  // make sure we provide the points in a clockwise fashion
1693 
1694  One = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
1695  Two = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
1696  Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
1697  Four = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour);
1698  Five = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
1699  Six = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
1700 
1701  chose = G4RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
1702  if( (chose>=0.) && (chose<aOne) )
1703  { return One; }
1704  else if( (chose>=aOne) && (chose<aOne+aTwo) )
1705  { return Two; }
1706  else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) )
1707  { return Three; }
1708  else if( (chose>=aOne+aTwo+aThree) && (chose<aOne+aTwo+aThree+aFour) )
1709  { return Four; }
1710  else if( (chose>=aOne+aTwo+aThree+aFour)
1711  && (chose<aOne+aTwo+aThree+aFour+aFive) )
1712  { return Five; }
1713  return Six;
1714 }
1715 
1717 //
1718 // Methods for visualisation
1719 
1721 {
1722  scene.AddSolid (*this);
1723 }
1724 
1726 {
1727  G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1728  G4double alpha1 = std::atan(fTalpha1);
1729  G4double alpha2 = std::atan(fTalpha2);
1730  G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi+fTthetaSphi*fTthetaSphi));
1731 
1732  return new G4PolyhedronTrap(fDz, theta, phi,
1733  fDy1, fDx1, fDx2, alpha1,
1734  fDy2, fDx3, fDx4, alpha2);
1735 }
1736 
1737 #endif
G4double fDz
Definition: G4Trap.hh:285
G4double fDx1
Definition: G4Trap.hh:286
G4double b
Definition: G4Trap.hh:103
G4String GetName() const
G4double GetXHalfLength4() const
G4double fTthetaCphi
Definition: G4Trap.hh:285
ThreeVector shoot(const G4int Ap, const G4int Af)
G4ThreeVector GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3, G4double &area) const
Definition: G4Trap.cc:1622
G4GeometryType GetEntityType() const
Definition: G4Trap.cc:1564
static constexpr double mm
Definition: G4SIunits.hh:115
G4double GetYHalfLength2() const
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
G4bool fRebuildPolyhedron
Definition: G4CSGSolid.hh:80
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:89
Definition: G4Trap.cc:82
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
Definition: G4Trap.cc:82
G4double GetZHalfLength() const
G4Trap & operator=(const G4Trap &rhs)
Definition: G4Trap.cc:576
G4double fCubicVolume
Definition: G4CSGSolid.hh:78
G4double fDx2
Definition: G4Trap.hh:286
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
G4double a
Definition: G4Trap.hh:103
G4double fTalpha2
Definition: G4Trap.hh:287
Definition: G4Trap.cc:82
G4double GetXHalfLength2() const
const G4int smax
G4double fDx4
Definition: G4Trap.hh:287
void DumpInfo() const
G4double GetTanAlpha2() const
G4double fDy2
Definition: G4Trap.hh:287
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Trap.cc:1230
G4double GetXHalfLength1() const
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79
G4GLOB_DLL std::ostream G4cout
virtual ~G4Trap()
Definition: G4Trap.cc:549
G4double GetXHalfLength3() const
G4VSolid * Clone() const
Definition: G4Trap.cc:1573
static constexpr double degree
Definition: G4SIunits.hh:144
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4Trap.cc:820
bool G4bool
Definition: G4Types.hh:79
G4double d
Definition: G4Trap.hh:103
Definition: G4Trap.cc:82
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:608
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
Definition: G4Trap.cc:867
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4double fTalpha1
Definition: G4Trap.hh:286
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
const G4int n
std::vector< G4ThreeVector > G4ThreeVectorList
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Trap.cc:809
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Trap.cc:967
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
Definition: G4Trap.cc:82
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4double fTthetaSphi
Definition: G4Trap.hh:285
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Trap.cc:1582
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
G4double fDx3
Definition: G4Trap.hh:287
EAxis
Definition: geomdefs.hh:54
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Trap.cc:1047
G4double fDy1
Definition: G4Trap.hh:286
Definition: G4Trap.cc:82
G4double c
Definition: G4Trap.hh:103
TrapSidePlane fPlanes[4]
Definition: G4Trap.hh:288
G4ThreeVector GetPointOnSurface() const
Definition: G4Trap.cc:1669
#define G4endl
Definition: G4ios.hh:61
Eside
Definition: G4Trap.cc:82
const G4double kCoplanar_Tolerance
Definition: G4Trap.cc:76
G4double kCarTolerance
Definition: G4VSolid.hh:307
Definition: G4Trap.cc:82
G4double GetYHalfLength1() const
G4Polyhedron * CreatePolyhedron() const
Definition: G4Trap.cc:1725
double G4double
Definition: G4Types.hh:76
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Trap.cc:1088
G4bool MakePlanes()
Definition: G4Trap.cc:655
G4bool MakePlane(const G4ThreeVector &p1, const G4ThreeVector &p2, const G4ThreeVector &p3, const G4ThreeVector &p4, TrapSidePlane &plane)
Definition: G4Trap.cc:734
G4double GetTanAlpha1() const
EInside Inside(const G4ThreeVector &p) const
Definition: G4Trap.cc:926
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:91
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Trap.cc:1720