Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4Trap.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/CSG/src/G4Trap.cc (Version 11.3.0) and /geometry/solids/CSG/src/G4Trap.cc (Version ReleaseNotes)


** Warning: Cannot open xref database.

  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 // Implementation for G4Trap class                
 27 //                                                
 28 // 21.03.95 P.Kent: Modified for `tolerant' ge    
 29 // 09.09.96 V.Grichine: Final modifications be    
 30 // 08.12.97 J.Allison: Added "nominal" constru    
 31 // 28.04.05 V.Grichine: new SurfaceNormal acco    
 32 // 18.04.17 E.Tcherniaev: complete revision, s    
 33 // -------------------------------------------    
 34                                                   
 35 #include "G4Trap.hh"                              
 36                                                   
 37 #if !defined(G4GEOM_USE_UTRAP)                    
 38                                                   
 39 #include "globals.hh"                             
 40 #include "G4GeomTools.hh"                         
 41                                                   
 42 #include "G4VoxelLimits.hh"                       
 43 #include "G4AffineTransform.hh"                   
 44 #include "G4BoundingEnvelope.hh"                  
 45                                                   
 46 #include "G4VPVParameterisation.hh"               
 47                                                   
 48 #include "G4QuickRand.hh"                         
 49                                                   
 50 #include "G4VGraphicsScene.hh"                    
 51 #include "G4Polyhedron.hh"                        
 52                                                   
 53 using namespace CLHEP;                            
 54                                                   
 55 //////////////////////////////////////////////    
 56 //                                                
 57 // Constructor - check and set half-widths as     
 58 // final check of coplanarity                     
 59                                                   
 60 G4Trap::G4Trap( const G4String& pName,            
 61                       G4double pDz,               
 62                       G4double pTheta, G4doubl    
 63                       G4double pDy1, G4double     
 64                       G4double pAlp1,             
 65                       G4double pDy2, G4double     
 66                       G4double pAlp2 )            
 67   : G4CSGSolid(pName), halfCarTolerance(0.5*kC    
 68 {                                                 
 69   fDz = pDz;                                      
 70   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi    
 71   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi    
 72                                                   
 73   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp    
 74   fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp    
 75                                                   
 76   CheckParameters();                              
 77   MakePlanes();                                   
 78 }                                                 
 79                                                   
 80 //////////////////////////////////////////////    
 81 //                                                
 82 // Constructor - Design of trapezoid based on     
 83 // which are its vertices. Checking of planari    
 84 // fPlanes[] and than calculation of other mem    
 85                                                   
 86 G4Trap::G4Trap( const G4String& pName,            
 87                 const G4ThreeVector pt[8] )       
 88   : G4CSGSolid(pName), halfCarTolerance(0.5*kC    
 89 {                                                 
 90   // Start with check of centering - the cente    
 91   // should cross the origin of frame             
 92   //                                              
 93   if (  pt[0].z() >= 0                            
 94         || pt[0].z() != pt[1].z()                 
 95         || pt[0].z() != pt[2].z()                 
 96         || pt[0].z() != pt[3].z()                 
 97                                                   
 98         || pt[4].z() <= 0                         
 99         || pt[4].z() != pt[5].z()                 
100         || pt[4].z() != pt[6].z()                 
101         || pt[4].z() != pt[7].z()                 
102                                                   
103         || std::fabs( pt[0].z() + pt[4].z() )     
104                                                   
105         || pt[0].y() != pt[1].y()                 
106         || pt[2].y() != pt[3].y()                 
107         || pt[4].y() != pt[5].y()                 
108         || pt[6].y() != pt[7].y()                 
109                                                   
110         || std::fabs(pt[0].y()+pt[2].y()+pt[4]    
111         || std::fabs(pt[0].x()+pt[1].x()+pt[4]    
112                      pt[2].x()+pt[3].x()+pt[6]    
113   {                                               
114     std::ostringstream message;                   
115     message << "Invalid vertice coordinates fo    
116     G4Exception("G4Trap::G4Trap()", "GeomSolid    
117                 FatalException, message);         
118   }                                               
119                                                   
120   // Set parameters                               
121   //                                              
122   fDz = (pt[7]).z();                              
123                                                   
124   fDy1     = ((pt[2]).y()-(pt[1]).y())*0.5;       
125   fDx1     = ((pt[1]).x()-(pt[0]).x())*0.5;       
126   fDx2     = ((pt[3]).x()-(pt[2]).x())*0.5;       
127   fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).    
128                                                   
129   fDy2     = ((pt[6]).y()-(pt[5]).y())*0.5;       
130   fDx3     = ((pt[5]).x()-(pt[4]).x())*0.5;       
131   fDx4     = ((pt[7]).x()-(pt[6]).x())*0.5;       
132   fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).    
133                                                   
134   fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx    
135   fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;           
136                                                   
137   CheckParameters();                              
138   MakePlanes(pt);                                 
139 }                                                 
140                                                   
141 //////////////////////////////////////////////    
142 //                                                
143 // Constructor for Right Angular Wedge from ST    
144                                                   
145 G4Trap::G4Trap( const G4String& pName,            
146                       G4double pZ,                
147                       G4double pY,                
148                       G4double pX, G4double pL    
149   : G4CSGSolid(pName), halfCarTolerance(0.5*kC    
150 {                                                 
151   fDz  = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi     
152   fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLT    
153   fDy2 = fDy1;   fDx3 = fDx1;   fDx4 = fDx2;      
154                                                   
155   CheckParameters();                              
156   MakePlanes();                                   
157 }                                                 
158                                                   
159 //////////////////////////////////////////////    
160 //                                                
161 // Constructor for G4Trd                          
162                                                   
163 G4Trap::G4Trap( const G4String& pName,            
164                       G4double pDx1,  G4double    
165                       G4double pDy1,  G4double    
166                       G4double pDz )              
167   : G4CSGSolid(pName), halfCarTolerance(0.5*kC    
168 {                                                 
169   fDz  = pDz;  fTthetaCphi = 0; fTthetaSphi =     
170   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalp    
171   fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalp    
172                                                   
173   CheckParameters();                              
174   MakePlanes();                                   
175 }                                                 
176                                                   
177 //////////////////////////////////////////////    
178 //                                                
179 // Constructor for G4Para                         
180                                                   
181 G4Trap::G4Trap( const G4String& pName,            
182                       G4double pDx, G4double p    
183                       G4double pDz,               
184                       G4double pAlpha,            
185                       G4double pTheta, G4doubl    
186   : G4CSGSolid(pName), halfCarTolerance(0.5*kC    
187 {                                                 
188   fDz = pDz;                                      
189   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi    
190   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi    
191                                                   
192   fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1    
193   fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2    
194                                                   
195   CheckParameters();                              
196   MakePlanes();                                   
197 }                                                 
198                                                   
199 //////////////////////////////////////////////    
200 //                                                
201 // Nominal constructor for G4Trap whose parame    
202 // a G4VParamaterisation later.  Check and set    
203 // angles: final check of coplanarity             
204                                                   
205 G4Trap::G4Trap( const G4String& pName )           
206   : G4CSGSolid (pName), halfCarTolerance(0.5*k    
207     fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),    
208     fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.)    
209     fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)    
210 {                                                 
211   MakePlanes();                                   
212 }                                                 
213                                                   
214 //////////////////////////////////////////////    
215 //                                                
216 // Fake default constructor - sets only member    
217 //                            for usage restri    
218 //                                                
219 G4Trap::G4Trap( __void__& a )                     
220   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo    
221     fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),    
222     fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.)    
223     fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)    
224 {                                                 
225   MakePlanes();                                   
226 }                                                 
227                                                   
228 //////////////////////////////////////////////    
229 //                                                
230 // Destructor                                     
231                                                   
232 G4Trap::~G4Trap() = default;                      
233                                                   
234 //////////////////////////////////////////////    
235 //                                                
236 // Copy constructor                               
237                                                   
238 G4Trap::G4Trap(const G4Trap& rhs)                 
239   : G4CSGSolid(rhs), halfCarTolerance(rhs.half    
240     fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi)    
241     fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.f    
242     fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.f    
243 {                                                 
244   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs    
245   for (G4int i=0; i<6; ++i) { fAreas[i] = rhs.    
246   fTrapType = rhs.fTrapType;                      
247 }                                                 
248                                                   
249 //////////////////////////////////////////////    
250 //                                                
251 // Assignment operator                            
252                                                   
253 G4Trap& G4Trap::operator = (const G4Trap& rhs)    
254 {                                                 
255   // Check assignment to self                     
256   //                                              
257   if (this == &rhs)  { return *this; }            
258                                                   
259   // Copy base class data                         
260   //                                              
261   G4CSGSolid::operator=(rhs);                     
262                                                   
263   // Copy data                                    
264   //                                              
265   halfCarTolerance = rhs.halfCarTolerance;        
266   fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi    
267   fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs    
268   fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs    
269   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs    
270   for (G4int i=0; i<6; ++i) { fAreas[i] = rhs.    
271   fTrapType = rhs.fTrapType;                      
272   return *this;                                   
273 }                                                 
274                                                   
275 //////////////////////////////////////////////    
276 //                                                
277 // Set all parameters, as for constructor - ch    
278 // as well as angles: final check of coplanari    
279                                                   
280 void G4Trap::SetAllParameters ( G4double pDz,     
281                                 G4double pThet    
282                                 G4double pPhi,    
283                                 G4double pDy1,    
284                                 G4double pDx1,    
285                                 G4double pDx2,    
286                                 G4double pAlp1    
287                                 G4double pDy2,    
288                                 G4double pDx3,    
289                                 G4double pDx4,    
290                                 G4double pAlp2    
291 {                                                 
292   // Reset data of the base class                 
293   fCubicVolume = 0;                               
294   fSurfaceArea = 0;                               
295   fRebuildPolyhedron = true;                      
296                                                   
297   // Set parameters                               
298   fDz = pDz;                                      
299   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi    
300   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi    
301                                                   
302   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp    
303   fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp    
304                                                   
305   CheckParameters();                              
306   MakePlanes();                                   
307 }                                                 
308                                                   
309 //////////////////////////////////////////////    
310 //                                                
311 // Check length parameters                        
312                                                   
313 void G4Trap::CheckParameters()                    
314 {                                                 
315   if (fDz<=0 ||                                   
316       fDy1<=0 || fDx1<=0 || fDx2<=0 ||            
317       fDy2<=0 || fDx3<=0 || fDx4<=0)              
318   {                                               
319     std::ostringstream message;                   
320     message << "Invalid Length Parameters for     
321             << "\n  X - " <<fDx1<<", "<<fDx2<<    
322             << "\n  Y - " <<fDy1<<", "<<fDy2      
323             << "\n  Z - " <<fDz;                  
324     G4Exception("G4Trap::CheckParameters()", "    
325                 FatalException, message);         
326   }                                               
327 }                                                 
328                                                   
329 //////////////////////////////////////////////    
330 //                                                
331 // Compute vertices and set side planes           
332                                                   
333 void G4Trap::MakePlanes()                         
334 {                                                 
335   G4double DzTthetaCphi = fDz*fTthetaCphi;        
336   G4double DzTthetaSphi = fDz*fTthetaSphi;        
337   G4double Dy1Talpha1   = fDy1*fTalpha1;          
338   G4double Dy2Talpha2   = fDy2*fTalpha2;          
339                                                   
340   G4ThreeVector pt[8] =                           
341   {                                               
342     G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx    
343     G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx    
344     G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx    
345     G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx    
346     G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx    
347     G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx    
348     G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx    
349     G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx    
350   };                                              
351                                                   
352   MakePlanes(pt);                                 
353 }                                                 
354                                                   
355 //////////////////////////////////////////////    
356 //                                                
357 // Set side planes, check planarity               
358                                                   
359 void G4Trap::MakePlanes(const G4ThreeVector pt    
360 {                                                 
361   constexpr G4int iface[4][4] = { {0,4,5,1}, {    
362   const static G4String side[4] = { "~-Y", "~+    
363                                                   
364   for (G4int i=0; i<4; ++i)                       
365   {                                               
366     if (MakePlane(pt[iface[i][0]],                
367                   pt[iface[i][1]],                
368                   pt[iface[i][2]],                
369                   pt[iface[i][3]],                
370                   fPlanes[i])) continue;          
371                                                   
372     // Non planar side face                       
373     G4ThreeVector normal(fPlanes[i].a,fPlanes[    
374     G4double dmax = 0;                            
375     for (G4int k=0; k<4; ++k)                     
376     {                                             
377       G4double dist = normal.dot(pt[iface[i][k    
378       if (std::abs(dist) > std::abs(dmax)) dma    
379     }                                             
380     std::ostringstream message;                   
381     message << "Side face " << side[i] << " is    
382             << GetName() << "\nDiscrepancy: "     
383     StreamInfo(message);                          
384     G4Exception("G4Trap::MakePlanes()", "GeomS    
385                 FatalException, message);         
386   }                                               
387                                                   
388   // Re-compute parameters                        
389   SetCachedValues();                              
390 }                                                 
391                                                   
392 //////////////////////////////////////////////    
393 //                                                
394 // Calculate the coef's of the plane p1->p2->p    
395 // where the ThreeVectors 1-4 are in anti-cloc    
396 // from infront of the plane (i.e. from normal    
397 //                                                
398 // Return true if the points are coplanar, fal    
399                                                   
400 G4bool G4Trap::MakePlane( const G4ThreeVector&    
401                           const G4ThreeVector&    
402                           const G4ThreeVector&    
403                           const G4ThreeVector&    
404                                 TrapSidePlane&    
405 {                                                 
406   G4ThreeVector normal = ((p4 - p2).cross(p3 -    
407   if (std::abs(normal.x()) < DBL_EPSILON) norm    
408   if (std::abs(normal.y()) < DBL_EPSILON) norm    
409   if (std::abs(normal.z()) < DBL_EPSILON) norm    
410   normal = normal.unit();                         
411                                                   
412   G4ThreeVector centre = (p1 + p2 + p3 + p4)*0    
413   plane.a =  normal.x();                          
414   plane.b =  normal.y();                          
415   plane.c =  normal.z();                          
416   plane.d = -normal.dot(centre);                  
417                                                   
418   // compute distances and check planarity        
419   G4double d1 = std::abs(normal.dot(p1) + plan    
420   G4double d2 = std::abs(normal.dot(p2) + plan    
421   G4double d3 = std::abs(normal.dot(p3) + plan    
422   G4double d4 = std::abs(normal.dot(p4) + plan    
423   G4double dmax = std::max(std::max(std::max(d    
424                                                   
425   return dmax <= 1000 * kCarTolerance;            
426 }                                                 
427                                                   
428 //////////////////////////////////////////////    
429 //                                                
430 // Recompute parameters using planes              
431                                                   
432 void G4Trap::SetCachedValues()                    
433 {                                                 
434   // Set indeces                                  
435   constexpr  G4int iface[6][4] =                  
436       { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,    
437                                                   
438   // Get vertices                                 
439   G4ThreeVector pt[8];                            
440   GetVertices(pt);                                
441                                                   
442   // Set face areas                               
443   for (G4int i=0; i<6; ++i)                       
444   {                                               
445     fAreas[i] = G4GeomTools::QuadAreaNormal(pt    
446                                             pt    
447                                             pt    
448                                             pt    
449   }                                               
450   for (G4int i=1; i<6; ++i) { fAreas[i] += fAr    
451                                                   
452   // Define type of trapezoid                     
453   fTrapType = 0;                                  
454   if (fPlanes[0].b  == -1 && fPlanes[1].b == 1    
455       std::abs(fPlanes[0].a) < DBL_EPSILON &&     
456       std::abs(fPlanes[0].c) < DBL_EPSILON &&     
457       std::abs(fPlanes[1].a) < DBL_EPSILON &&     
458       std::abs(fPlanes[1].c) < DBL_EPSILON)       
459   {                                               
460     fTrapType = 1; // YZ section is a rectangl    
461     if (std::abs(fPlanes[2].a + fPlanes[3].a)     
462         std::abs(fPlanes[2].c - fPlanes[3].c)     
463         fPlanes[2].b == 0 &&                      
464         fPlanes[3].b == 0)                        
465     {                                             
466       fTrapType = 2; // ... and XZ section is     
467       fPlanes[2].a = -fPlanes[3].a;               
468       fPlanes[2].c =  fPlanes[3].c;               
469     }                                             
470     if (std::abs(fPlanes[2].a + fPlanes[3].a)     
471         std::abs(fPlanes[2].b - fPlanes[3].b)     
472         fPlanes[2].c == 0 &&                      
473         fPlanes[3].c == 0)                        
474     {                                             
475       fTrapType = 3; // ... and XY section is     
476       fPlanes[2].a = -fPlanes[3].a;               
477       fPlanes[2].b =  fPlanes[3].b;               
478     }                                             
479   }                                               
480 }                                                 
481                                                   
482 //////////////////////////////////////////////    
483 //                                                
484 // Get volume                                     
485                                                   
486 G4double G4Trap::GetCubicVolume()                 
487 {                                                 
488   if (fCubicVolume == 0)                          
489   {                                               
490     G4ThreeVector pt[8];                          
491     GetVertices(pt);                              
492                                                   
493     G4double dz  = pt[4].z() - pt[0].z();         
494     G4double dy1 = pt[2].y() - pt[0].y();         
495     G4double dx1 = pt[1].x() - pt[0].x();         
496     G4double dx2 = pt[3].x() - pt[2].x();         
497     G4double dy2 = pt[6].y() - pt[4].y();         
498     G4double dx3 = pt[5].x() - pt[4].x();         
499     G4double dx4 = pt[7].x() - pt[6].x();         
500                                                   
501     fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(d    
502                     (dx4 + dx3 - dx2 - dx1)*(d    
503   }                                               
504   return fCubicVolume;                            
505 }                                                 
506                                                   
507 //////////////////////////////////////////////    
508 //                                                
509 // Get surface area                               
510                                                   
511 G4double G4Trap::GetSurfaceArea()                 
512 {                                                 
513   if (fSurfaceArea == 0)                          
514   {                                               
515     G4ThreeVector pt[8];                          
516     G4int iface [6][4] =                          
517       { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,    
518                                                   
519     GetVertices(pt);                              
520     for (const auto & i : iface)                  
521     {                                             
522       fSurfaceArea += G4GeomTools::QuadAreaNor    
523                                                   
524                                                   
525                                                   
526     }                                             
527   }                                               
528   return fSurfaceArea;                            
529 }                                                 
530                                                   
531 //////////////////////////////////////////////    
532 //                                                
533 // Dispatch to parameterisation for replicatio    
534 // computation & modification.                    
535                                                   
536 void G4Trap::ComputeDimensions(       G4VPVPar    
537                                 const G4int n,    
538                                 const G4VPhysi    
539 {                                                 
540   p->ComputeDimensions(*this,n,pRep);             
541 }                                                 
542                                                   
543 //////////////////////////////////////////////    
544 //                                                
545 // Get bounding box                               
546                                                   
547 void G4Trap::BoundingLimits(G4ThreeVector& pMi    
548 {                                                 
549   G4ThreeVector pt[8];                            
550   GetVertices(pt);                                
551                                                   
552   G4double xmin = kInfinity, xmax = -kInfinity    
553   G4double ymin = kInfinity, ymax = -kInfinity    
554   for (const auto & i : pt)                       
555   {                                               
556     G4double x = i.x();                           
557     if (x < xmin) xmin = x;                       
558     if (x > xmax) xmax = x;                       
559     G4double y = i.y();                           
560     if (y < ymin) ymin = y;                       
561     if (y > ymax) ymax = y;                       
562   }                                               
563                                                   
564   G4double dz   = GetZHalfLength();               
565   pMin.set(xmin,ymin,-dz);                        
566   pMax.set(xmax,ymax, dz);                        
567                                                   
568   // Check correctness of the bounding box        
569   //                                              
570   if (pMin.x() >= pMax.x() || pMin.y() >= pMax    
571   {                                               
572     std::ostringstream message;                   
573     message << "Bad bounding box (min >= max)     
574             << GetName() << " !"                  
575             << "\npMin = " << pMin                
576             << "\npMax = " << pMax;               
577     G4Exception("G4Trap::BoundingLimits()", "G    
578                 JustWarning, message);            
579     DumpInfo();                                   
580   }                                               
581 }                                                 
582                                                   
583 //////////////////////////////////////////////    
584 //                                                
585 // Calculate extent under transform and specif    
586                                                   
587 G4bool G4Trap::CalculateExtent( const EAxis pA    
588                                 const G4VoxelL    
589                                 const G4Affine    
590                                       G4double    
591 {                                                 
592   G4ThreeVector bmin, bmax;                       
593   G4bool exist;                                   
594                                                   
595   // Check bounding box (bbox)                    
596   //                                              
597   BoundingLimits(bmin,bmax);                      
598   G4BoundingEnvelope bbox(bmin,bmax);             
599 #ifdef G4BBOX_EXTENT                              
600   return bbox.CalculateExtent(pAxis,pVoxelLimi    
601 #endif                                            
602   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox    
603   {                                               
604     return exist = pMin < pMax;                   
605   }                                               
606                                                   
607   // Set bounding envelope (benv) and calculat    
608   //                                              
609   G4ThreeVector pt[8];                            
610   GetVertices(pt);                                
611                                                   
612   G4ThreeVectorList baseA(4), baseB(4);           
613   baseA[0] = pt[0];                               
614   baseA[1] = pt[1];                               
615   baseA[2] = pt[3];                               
616   baseA[3] = pt[2];                               
617                                                   
618   baseB[0] = pt[4];                               
619   baseB[1] = pt[5];                               
620   baseB[2] = pt[7];                               
621   baseB[3] = pt[6];                               
622                                                   
623   std::vector<const G4ThreeVectorList *> polyg    
624   polygons[0] = &baseA;                           
625   polygons[1] = &baseB;                           
626                                                   
627   G4BoundingEnvelope benv(bmin,bmax,polygons);    
628   exist = benv.CalculateExtent(pAxis,pVoxelLim    
629   return exist;                                   
630 }                                                 
631                                                   
632 //////////////////////////////////////////////    
633 //                                                
634 // Return whether point is inside/outside/on_s    
635                                                   
636 EInside G4Trap::Inside( const G4ThreeVector& p    
637 {                                                 
638   switch (fTrapType)                              
639   {                                               
640     case 0: // General case                       
641     {                                             
642       G4double dz = std::abs(p.z())-fDz;          
643       G4double dy1 = fPlanes[0].b*p.y()+fPlane    
644       G4double dy2 = fPlanes[1].b*p.y()+fPlane    
645       G4double dy = std::max(dz,std::max(dy1,d    
646                                                   
647       G4double dx1 = fPlanes[2].a*p.x()+fPlane    
648                    + fPlanes[2].c*p.z()+fPlane    
649       G4double dx2 = fPlanes[3].a*p.x()+fPlane    
650                    + fPlanes[3].c*p.z()+fPlane    
651       G4double dist = std::max(dy,std::max(dx1    
652                                                   
653       return (dist > halfCarTolerance) ? kOuts    
654         ((dist > -halfCarTolerance) ? kSurface    
655     }                                             
656     case 1: // YZ section is a rectangle          
657     {                                             
658       G4double dz = std::abs(p.z())-fDz;          
659       G4double dy = std::max(dz,std::abs(p.y()    
660       G4double dx1 = fPlanes[2].a*p.x()+fPlane    
661                    + fPlanes[2].c*p.z()+fPlane    
662       G4double dx2 = fPlanes[3].a*p.x()+fPlane    
663                    + fPlanes[3].c*p.z()+fPlane    
664       G4double dist = std::max(dy,std::max(dx1    
665                                                   
666       return (dist > halfCarTolerance) ? kOuts    
667         ((dist > -halfCarTolerance) ? kSurface    
668     }                                             
669     case 2: // YZ section is a rectangle and      
670     {       // XZ section is an isosceles trap    
671       G4double dz = std::abs(p.z())-fDz;          
672       G4double dy = std::max(dz,std::abs(p.y()    
673       G4double dx = fPlanes[3].a*std::abs(p.x(    
674                   + fPlanes[3].c*p.z()+fPlanes    
675       G4double dist = std::max(dy,dx);            
676                                                   
677       return (dist > halfCarTolerance) ? kOuts    
678         ((dist > -halfCarTolerance) ? kSurface    
679     }                                             
680     case 3: // YZ section is a rectangle and      
681     {       // XY section is an isosceles trap    
682       G4double dz = std::abs(p.z())-fDz;          
683       G4double dy = std::max(dz,std::abs(p.y()    
684       G4double dx = fPlanes[3].a*std::abs(p.x(    
685                   + fPlanes[3].b*p.y()+fPlanes    
686       G4double dist = std::max(dy,dx);            
687                                                   
688       return (dist > halfCarTolerance) ? kOuts    
689         ((dist > -halfCarTolerance) ? kSurface    
690     }                                             
691   }                                               
692   return kOutside;                                
693 }                                                 
694                                                   
695 //////////////////////////////////////////////    
696 //                                                
697 // Determine side, and return corresponding no    
698                                                   
699 G4ThreeVector G4Trap::SurfaceNormal( const G4T    
700 {                                                 
701   G4double nx = 0, ny = 0, nz = 0;                
702   G4double dz = std::abs(p.z()) - fDz;            
703   nz = std::copysign(G4double(std::abs(dz) <=     
704                                                   
705   switch (fTrapType)                              
706   {                                               
707     case 0: // General case                       
708     {                                             
709       for (G4int i=0; i<2; ++i)                   
710       {                                           
711         G4double dy = fPlanes[i].b*p.y() + fPl    
712         if (std::abs(dy) > halfCarTolerance) c    
713         ny  = fPlanes[i].b;                       
714         nz += fPlanes[i].c;                       
715         break;                                    
716       }                                           
717       for (G4int i=2; i<4; ++i)                   
718       {                                           
719         G4double dx = fPlanes[i].a*p.x() +        
720                       fPlanes[i].b*p.y() + fPl    
721         if (std::abs(dx) > halfCarTolerance) c    
722         nx  = fPlanes[i].a;                       
723         ny += fPlanes[i].b;                       
724         nz += fPlanes[i].c;                       
725         break;                                    
726       }                                           
727       break;                                      
728     }                                             
729     case 1: // YZ section - rectangle             
730     {                                             
731       G4double dy = std::abs(p.y()) + fPlanes[    
732       ny = std::copysign(G4double(std::abs(dy)    
733       for (G4int i=2; i<4; ++i)                   
734       {                                           
735         G4double dx = fPlanes[i].a*p.x() +        
736                       fPlanes[i].b*p.y() + fPl    
737         if (std::abs(dx) > halfCarTolerance) c    
738         nx  = fPlanes[i].a;                       
739         ny += fPlanes[i].b;                       
740         nz += fPlanes[i].c;                       
741         break;                                    
742       }                                           
743       break;                                      
744     }                                             
745     case 2: // YZ section - rectangle, XZ sect    
746     {                                             
747       G4double dy = std::abs(p.y()) + fPlanes[    
748       ny = std::copysign(G4double(std::abs(dy)    
749       G4double dx = fPlanes[3].a*std::abs(p.x(    
750                     fPlanes[3].c*p.z() + fPlan    
751       G4double k = std::abs(dx) <= halfCarTole    
752       nx  = std::copysign(k, p.x())*fPlanes[3]    
753       nz += k*fPlanes[3].c;                       
754       break;                                      
755     }                                             
756     case 3: // YZ section - rectangle, XY sect    
757     {                                             
758       G4double dy = std::abs(p.y()) + fPlanes[    
759       ny = std::copysign(G4double(std::abs(dy)    
760       G4double dx = fPlanes[3].a*std::abs(p.x(    
761                     fPlanes[3].b*p.y() + fPlan    
762       G4double k = std::abs(dx) <= halfCarTole    
763       nx  = std::copysign(k, p.x())*fPlanes[3]    
764       ny += k*fPlanes[3].b;                       
765       break;                                      
766     }                                             
767   }                                               
768                                                   
769   // Return normal                                
770   //                                              
771   G4double mag2 = nx*nx + ny*ny + nz*nz;          
772   if (mag2 == 1)      return { nx,ny,nz };        
773   else if (mag2 != 0) return G4ThreeVector(nx,    
774   else                                            
775   {                                               
776     // Point is not on the surface                
777     //                                            
778 #ifdef G4CSGDEBUG                                 
779     std::ostringstream message;                   
780     G4long oldprc = message.precision(16);        
781     message << "Point p is not on surface (!?)    
782             << GetName() << G4endl;               
783     message << "Position:\n";                     
784     message << "   p.x() = " << p.x()/mm << "     
785     message << "   p.y() = " << p.y()/mm << "     
786     message << "   p.z() = " << p.z()/mm << "     
787     G4cout.precision(oldprc) ;                    
788     G4Exception("G4Trap::SurfaceNormal(p)", "G    
789                 JustWarning, message );           
790     DumpInfo();                                   
791 #endif                                            
792     return ApproxSurfaceNormal(p);                
793   }                                               
794 }                                                 
795                                                   
796 //////////////////////////////////////////////    
797 //                                                
798 // Algorithm for SurfaceNormal() following the    
799 // for points not on the surface                  
800                                                   
801 G4ThreeVector G4Trap::ApproxSurfaceNormal( con    
802 {                                                 
803   G4double dist = -DBL_MAX;                       
804   G4int iside = 0;                                
805   for (G4int i=0; i<4; ++i)                       
806   {                                               
807     G4double d = fPlanes[i].a*p.x() +             
808                  fPlanes[i].b*p.y() +             
809                  fPlanes[i].c*p.z() + fPlanes[    
810     if (d > dist) { dist = d; iside = i; }        
811   }                                               
812                                                   
813   G4double distz = std::abs(p.z()) - fDz;         
814   if (dist > distz)                               
815     return { fPlanes[iside].a, fPlanes[iside].    
816   else                                            
817     return { 0, 0, (G4double)((p.z() < 0) ? -1    
818 }                                                 
819                                                   
820 //////////////////////////////////////////////    
821 //                                                
822 // Calculate distance to shape from outside       
823 //  - return kInfinity if no intersection         
824                                                   
825 G4double G4Trap::DistanceToIn(const G4ThreeVec    
826                               const G4ThreeVec    
827 {                                                 
828   // Z intersections                              
829   //                                              
830   if ((std::abs(p.z()) - fDz) >= -halfCarToler    
831     return kInfinity;                             
832   G4double invz = (-v.z() == 0) ? DBL_MAX : -1    
833   G4double dz = (invz < 0) ? fDz : -fDz;          
834   G4double tzmin = (p.z() + dz)*invz;             
835   G4double tzmax = (p.z() - dz)*invz;             
836                                                   
837   // Y intersections                              
838   //                                              
839   G4double tymin = 0, tymax = DBL_MAX;            
840   G4int i = 0;                                    
841   for ( ; i<2; ++i)                               
842   {                                               
843     G4double cosa = fPlanes[i].b*v.y() + fPlan    
844     G4double dist = fPlanes[i].b*p.y() + fPlan    
845     if (dist >= -halfCarTolerance)                
846     {                                             
847       if (cosa >= 0) return kInfinity;            
848       G4double tmp  = -dist/cosa;                 
849       if (tymin < tmp) tymin = tmp;               
850     }                                             
851     else if (cosa > 0)                            
852     {                                             
853       G4double tmp  = -dist/cosa;                 
854       if (tymax > tmp) tymax = tmp;               
855     }                                             
856   }                                               
857                                                   
858   // Z intersections                              
859   //                                              
860   G4double txmin = 0, txmax = DBL_MAX;            
861   for ( ; i<4; ++i)                               
862   {                                               
863     G4double cosa = fPlanes[i].a*v.x()+fPlanes    
864     G4double dist = fPlanes[i].a*p.x()+fPlanes    
865                     fPlanes[i].d;                 
866     if (dist >= -halfCarTolerance)                
867     {                                             
868       if (cosa >= 0) return kInfinity;            
869       G4double tmp  = -dist/cosa;                 
870       if (txmin < tmp) txmin = tmp;               
871     }                                             
872     else if (cosa > 0)                            
873     {                                             
874       G4double tmp  = -dist/cosa;                 
875       if (txmax > tmp) txmax = tmp;               
876     }                                             
877   }                                               
878                                                   
879   // Find distance                                
880   //                                              
881   G4double tmin = std::max(std::max(txmin,tymi    
882   G4double tmax = std::min(std::min(txmax,tyma    
883                                                   
884   if (tmax <= tmin + halfCarTolerance) return     
885   return (tmin < halfCarTolerance ) ? 0. : tmi    
886 }                                                 
887                                                   
888 //////////////////////////////////////////////    
889 //                                                
890 // Calculate exact shortest distance to any bo    
891 // This is the best fast estimation of the sho    
892 // - return 0 if point is inside                  
893                                                   
894 G4double G4Trap::DistanceToIn( const G4ThreeVe    
895 {                                                 
896   switch (fTrapType)                              
897   {                                               
898     case 0: // General case                       
899     {                                             
900       G4double dz = std::abs(p.z())-fDz;          
901       G4double dy1 = fPlanes[0].b*p.y()+fPlane    
902       G4double dy2 = fPlanes[1].b*p.y()+fPlane    
903       G4double dy = std::max(dz,std::max(dy1,d    
904                                                   
905       G4double dx1 = fPlanes[2].a*p.x()+fPlane    
906                    + fPlanes[2].c*p.z()+fPlane    
907       G4double dx2 = fPlanes[3].a*p.x()+fPlane    
908                    + fPlanes[3].c*p.z()+fPlane    
909       G4double dist = std::max(dy,std::max(dx1    
910       return (dist > 0) ? dist : 0.;              
911     }                                             
912     case 1: // YZ section is a rectangle          
913     {                                             
914       G4double dz = std::abs(p.z())-fDz;          
915       G4double dy = std::max(dz,std::abs(p.y()    
916       G4double dx1 = fPlanes[2].a*p.x()+fPlane    
917                    + fPlanes[2].c*p.z()+fPlane    
918       G4double dx2 = fPlanes[3].a*p.x()+fPlane    
919                    + fPlanes[3].c*p.z()+fPlane    
920       G4double dist = std::max(dy,std::max(dx1    
921       return (dist > 0) ? dist : 0.;              
922     }                                             
923     case 2: // YZ section is a rectangle and      
924     {       // XZ section is an isosceles trap    
925       G4double dz = std::abs(p.z())-fDz;          
926       G4double dy = std::max(dz,std::abs(p.y()    
927       G4double dx = fPlanes[3].a*std::abs(p.x(    
928                   + fPlanes[3].c*p.z()+fPlanes    
929       G4double dist = std::max(dy,dx);            
930       return (dist > 0) ? dist : 0.;              
931     }                                             
932     case 3: // YZ section is a rectangle and      
933     {       // XY section is an isosceles trap    
934       G4double dz = std::abs(p.z())-fDz;          
935       G4double dy = std::max(dz,std::abs(p.y()    
936       G4double dx = fPlanes[3].a*std::abs(p.x(    
937                   + fPlanes[3].b*p.y()+fPlanes    
938       G4double dist = std::max(dy,dx);            
939       return (dist > 0) ? dist : 0.;              
940     }                                             
941   }                                               
942   return 0.;                                      
943 }                                                 
944                                                   
945 //////////////////////////////////////////////    
946 //                                                
947 // Calculate distance to surface of shape from    
948 // find normal at exit point, if required         
949 // - when leaving the surface, return 0           
950                                                   
951 G4double G4Trap::DistanceToOut(const G4ThreeVe    
952                                const G4bool ca    
953                                      G4bool* v    
954 {                                                 
955   // Z intersections                              
956   //                                              
957   if ((std::abs(p.z()) - fDz) >= -halfCarToler    
958   {                                               
959     if (calcNorm)                                 
960     {                                             
961       *validNorm = true;                          
962       n->set(0, 0, (p.z() < 0) ? -1 : 1);         
963     }                                             
964     return 0;                                     
965   }                                               
966   G4double vz = v.z();                            
967   G4double tmax = (vz == 0) ? DBL_MAX : (std::    
968   G4int iside = (vz < 0) ? -4 : -2; // little     
969                                                   
970   // Y intersections                              
971   //                                              
972   G4int i = 0;                                    
973   for ( ; i<2; ++i)                               
974   {                                               
975     G4double cosa = fPlanes[i].b*v.y() + fPlan    
976     if (cosa > 0)                                 
977     {                                             
978       G4double dist = fPlanes[i].b*p.y() + fPl    
979       if (dist >= -halfCarTolerance)              
980       {                                           
981         if (calcNorm)                             
982         {                                         
983           *validNorm = true;                      
984           n->set(0, fPlanes[i].b, fPlanes[i].c    
985         }                                         
986         return 0;                                 
987       }                                           
988       G4double tmp = -dist/cosa;                  
989       if (tmax > tmp) { tmax = tmp; iside = i;    
990     }                                             
991   }                                               
992                                                   
993   // X intersections                              
994   //                                              
995   for ( ; i<4; ++i)                               
996   {                                               
997     G4double cosa = fPlanes[i].a*v.x()+fPlanes    
998     if (cosa > 0)                                 
999     {                                             
1000       G4double dist = fPlanes[i].a*p.x() +       
1001                       fPlanes[i].b*p.y() + fP    
1002       if (dist >= -halfCarTolerance)             
1003       {                                          
1004         if (calcNorm)                            
1005         {                                        
1006            *validNorm = true;                    
1007            n->set(fPlanes[i].a, fPlanes[i].b,    
1008         }                                        
1009         return 0;                                
1010       }                                          
1011       G4double tmp = -dist/cosa;                 
1012       if (tmax > tmp) { tmax = tmp; iside = i    
1013     }                                            
1014   }                                              
1015                                                  
1016   // Set normal, if required, and return dist    
1017   //                                             
1018   if (calcNorm)                                  
1019   {                                              
1020     *validNorm = true;                           
1021     if (iside < 0)                               
1022       n->set(0, 0, iside + 3); // (-4+3)=-1,     
1023     else                                         
1024       n->set(fPlanes[iside].a, fPlanes[iside]    
1025   }                                              
1026   return tmax;                                   
1027 }                                                
1028                                                  
1029 /////////////////////////////////////////////    
1030 //                                               
1031 // Calculate exact shortest distance to any b    
1032 // - Returns 0 is ThreeVector outside            
1033                                                  
1034 G4double G4Trap::DistanceToOut( const G4Three    
1035 {                                                
1036 #ifdef G4CSGDEBUG                                
1037   if( Inside(p) == kOutside )                    
1038   {                                              
1039     std::ostringstream message;                  
1040     G4long oldprc = message.precision(16);       
1041     message << "Point p is outside (!?) of so    
1042     message << "Position:\n";                    
1043     message << "   p.x() = " << p.x()/mm << "    
1044     message << "   p.y() = " << p.y()/mm << "    
1045     message << "   p.z() = " << p.z()/mm << "    
1046     G4cout.precision(oldprc);                    
1047     G4Exception("G4Trap::DistanceToOut(p)", "    
1048                 JustWarning, message );          
1049     DumpInfo();                                  
1050   }                                              
1051 #endif                                           
1052   switch (fTrapType)                             
1053   {                                              
1054     case 0: // General case                      
1055     {                                            
1056       G4double dz = std::abs(p.z())-fDz;         
1057       G4double dy1 = fPlanes[0].b*p.y()+fPlan    
1058       G4double dy2 = fPlanes[1].b*p.y()+fPlan    
1059       G4double dy = std::max(dz,std::max(dy1,    
1060                                                  
1061       G4double dx1 = fPlanes[2].a*p.x()+fPlan    
1062                    + fPlanes[2].c*p.z()+fPlan    
1063       G4double dx2 = fPlanes[3].a*p.x()+fPlan    
1064                    + fPlanes[3].c*p.z()+fPlan    
1065       G4double dist = std::max(dy,std::max(dx    
1066       return (dist < 0) ? -dist : 0.;            
1067     }                                            
1068     case 1: // YZ section is a rectangle         
1069     {                                            
1070       G4double dz = std::abs(p.z())-fDz;         
1071       G4double dy = std::max(dz,std::abs(p.y(    
1072       G4double dx1 = fPlanes[2].a*p.x()+fPlan    
1073                    + fPlanes[2].c*p.z()+fPlan    
1074       G4double dx2 = fPlanes[3].a*p.x()+fPlan    
1075                    + fPlanes[3].c*p.z()+fPlan    
1076       G4double dist = std::max(dy,std::max(dx    
1077       return (dist < 0) ? -dist : 0.;            
1078     }                                            
1079     case 2: // YZ section is a rectangle and     
1080     {       // XZ section is an isosceles tra    
1081       G4double dz = std::abs(p.z())-fDz;         
1082       G4double dy = std::max(dz,std::abs(p.y(    
1083       G4double dx = fPlanes[3].a*std::abs(p.x    
1084                   + fPlanes[3].c*p.z()+fPlane    
1085       G4double dist = std::max(dy,dx);           
1086       return (dist < 0) ? -dist : 0.;            
1087     }                                            
1088     case 3: // YZ section is a rectangle and     
1089     {       // XY section is an isosceles tra    
1090       G4double dz = std::abs(p.z())-fDz;         
1091       G4double dy = std::max(dz,std::abs(p.y(    
1092       G4double dx = fPlanes[3].a*std::abs(p.x    
1093                   + fPlanes[3].b*p.y()+fPlane    
1094       G4double dist = std::max(dy,dx);           
1095       return (dist < 0) ? -dist : 0.;            
1096     }                                            
1097   }                                              
1098   return 0.;                                     
1099 }                                                
1100                                                  
1101 /////////////////////////////////////////////    
1102 //                                               
1103 // GetEntityType                                 
1104                                                  
1105 G4GeometryType G4Trap::GetEntityType() const     
1106 {                                                
1107   return {"G4Trap"};                             
1108 }                                                
1109                                                  
1110 /////////////////////////////////////////////    
1111 //                                               
1112 // IsFaceted                                     
1113                                                  
1114 G4bool G4Trap::IsFaceted() const                 
1115 {                                                
1116   return true;                                   
1117 }                                                
1118                                                  
1119 /////////////////////////////////////////////    
1120 //                                               
1121 // Make a clone of the object                    
1122 //                                               
1123 G4VSolid* G4Trap::Clone() const                  
1124 {                                                
1125   return new G4Trap(*this);                      
1126 }                                                
1127                                                  
1128 /////////////////////////////////////////////    
1129 //                                               
1130 // Stream object contents to an output stream    
1131                                                  
1132 std::ostream& G4Trap::StreamInfo( std::ostrea    
1133 {                                                
1134   G4double phi    = GetPhi();                    
1135   G4double theta  = GetTheta();                  
1136   G4double alpha1 = GetAlpha1();                 
1137   G4double alpha2 = GetAlpha2();                 
1138                                                  
1139   G4long oldprc = os.precision(16);              
1140   os << "------------------------------------    
1141      << "    *** Dump for solid: " << GetName    
1142      << "    ================================    
1143      << " Solid type: G4Trap\n"                  
1144      << " Parameters:\n"                         
1145      << "    half length Z: " << fDz/mm << "     
1146      << "    half length Y, face -Dz: " << fD    
1147      << "    half length X, face -Dz, side -D    
1148      << "    half length X, face -Dz, side +D    
1149      << "    half length Y, face +Dz: " << fD    
1150      << "    half length X, face +Dz, side -D    
1151      << "    half length X, face +Dz, side +D    
1152      << "    theta: " << theta/degree << " de    
1153      << "    phi:   " << phi/degree << " degr    
1154      << "    alpha, face -Dz: " << alpha1/deg    
1155      << "    alpha, face +Dz: " << alpha2/deg    
1156      << "------------------------------------    
1157   os.precision(oldprc);                          
1158                                                  
1159   return os;                                     
1160 }                                                
1161                                                  
1162 /////////////////////////////////////////////    
1163 //                                               
1164 // Compute vertices from planes                  
1165                                                  
1166 void G4Trap::GetVertices(G4ThreeVector pt[8])    
1167 {                                                
1168   for (G4int i=0; i<8; ++i)                      
1169   {                                              
1170     G4int iy = (i==0 || i==1 || i==4 || i==5)    
1171     G4int ix = (i==0 || i==2 || i==4 || i==6)    
1172     G4double z = (i < 4) ? -fDz : fDz;           
1173     G4double y = -(fPlanes[iy].c*z + fPlanes[    
1174     G4double x = -(fPlanes[ix].b*y + fPlanes[    
1175                    + fPlanes[ix].d)/fPlanes[i    
1176     pt[i].set(x,y,z);                            
1177   }                                              
1178 }                                                
1179                                                  
1180 /////////////////////////////////////////////    
1181 //                                               
1182 // Generate random point on the surface          
1183                                                  
1184 G4ThreeVector G4Trap::GetPointOnSurface() con    
1185 {                                                
1186   // Set indeces                                 
1187   constexpr G4int iface [6][4] =                 
1188     { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6    
1189                                                  
1190   // Set vertices                                
1191   G4ThreeVector pt[8];                           
1192   GetVertices(pt);                               
1193                                                  
1194   // Select face                                 
1195   //                                             
1196   G4double select = fAreas[5]*G4QuickRand();     
1197   G4int k = 5;                                   
1198   k -= (G4int)(select <= fAreas[4]);             
1199   k -= (G4int)(select <= fAreas[3]);             
1200   k -= (G4int)(select <= fAreas[2]);             
1201   k -= (G4int)(select <= fAreas[1]);             
1202   k -= (G4int)(select <= fAreas[0]);             
1203                                                  
1204   // Select sub-triangle                         
1205   //                                             
1206   G4int i0 = iface[k][0];                        
1207   G4int i1 = iface[k][1];                        
1208   G4int i2 = iface[k][2];                        
1209   G4int i3 = iface[k][3];                        
1210   G4double s2 = G4GeomTools::TriangleAreaNorm    
1211   if (select > fAreas[k] - s2) i0 = i2;          
1212                                                  
1213   // Generate point                              
1214   //                                             
1215   G4double u = G4QuickRand();                    
1216   G4double v = G4QuickRand();                    
1217   if (u + v > 1.) { u = 1. - u; v = 1. - v; }    
1218   return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3    
1219 }                                                
1220                                                  
1221 /////////////////////////////////////////////    
1222 //                                               
1223 // Methods for visualisation                     
1224                                                  
1225 void G4Trap::DescribeYourselfTo ( G4VGraphics    
1226 {                                                
1227   scene.AddSolid (*this);                        
1228 }                                                
1229                                                  
1230 G4Polyhedron* G4Trap::CreatePolyhedron () con    
1231 {                                                
1232   G4double phi = std::atan2(fTthetaSphi, fTth    
1233   G4double alpha1 = std::atan(fTalpha1);         
1234   G4double alpha2 = std::atan(fTalpha2);         
1235   G4double theta = std::atan(std::sqrt(fTthet    
1236                                       +fTthet    
1237                                                  
1238   return new G4PolyhedronTrap(fDz, theta, phi    
1239                               fDy1, fDx1, fDx    
1240                               fDy2, fDx3, fDx    
1241 }                                                
1242                                                  
1243 #endif                                           
1244