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