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

Last change on this file was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
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: geant4-09-04-beta-01 $
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.