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