source: trunk/source/geometry/solids/CSG/src/G4Trap.cc@ 1044

Last change on this file since 1044 was 850, checked in by garnier, 17 years ago

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