| 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 // G4GenericTrap implementation << 26 //
>> 27 // $Id: G4GenericTrap.cc,v 1.21 2010-11-26 13:30:26 gcosmo Exp $
>> 28 // GEANT4 tag $Name: not supported by cvs2svn $
>> 29 //
>> 30 //
>> 31 // --------------------------------------------------------------------
>> 32 // GEANT 4 class source file
>> 33 //
>> 34 // G4GenericTrap.cc
27 // 35 //
28 // Authors: 36 // Authors:
29 // Tatiana Nikitina, CERN; Ivana Hrivnacova, 37 // Tatiana Nikitina, CERN; Ivana Hrivnacova, IPN Orsay
30 // Adapted from Root Arb8 implementation by << 38 // Adapted from Root Arb8 implementation by Andrei Gheata, CERN
31 // 39 //
32 // 27.05.2024 - Evgueni Tcherniaev, complete << 40 // History :
>> 41 // 04 August 2011 T.Nikitina Add SetReferences() and InvertFacets()
>> 42 // to CreatePolyhedron() for Visualisation of Boolean
33 // ------------------------------------------- 43 // --------------------------------------------------------------------
34 44
35 #include "G4GenericTrap.hh" <<
36 <<
37 #if !defined(G4GEOM_USE_UGENERICTRAP) <<
38 <<
39 #include <iomanip> 45 #include <iomanip>
40 46
41 #include "G4PhysicalConstants.hh" << 47 #include "G4GenericTrap.hh"
42 #include "G4SystemOfUnits.hh" << 48 #include "G4TessellatedSolid.hh"
43 #include "G4VoxelLimits.hh" << 49 #include "G4TriangularFacet.hh"
44 #include "G4AffineTransform.hh" << 50 #include "G4QuadrangularFacet.hh"
45 #include "G4BoundingEnvelope.hh" << 51 #include "Randomize.hh"
46 #include "G4QuickRand.hh" <<
47 <<
48 #include "G4GeomTools.hh" <<
49 52
50 #include "G4VGraphicsScene.hh" 53 #include "G4VGraphicsScene.hh"
51 #include "G4Polyhedron.hh" 54 #include "G4Polyhedron.hh"
>> 55 #include "G4PolyhedronArbitrary.hh"
>> 56 #include "G4NURBS.hh"
>> 57 #include "G4NURBSbox.hh"
52 #include "G4VisExtent.hh" 58 #include "G4VisExtent.hh"
53 59
54 #include "G4AutoLock.hh" << 60 const G4int G4GenericTrap::fgkNofVertices = 8;
>> 61 const G4double G4GenericTrap::fgkTolerance = 1E-3;
55 62
56 namespace << 63 // --------------------------------------------------------------------
57 { <<
58 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZ <<
59 G4Mutex lateralareaMutex = G4MUTEX_INITIALIZ <<
60 } <<
61 64
62 ////////////////////////////////////////////// << 65 G4GenericTrap::G4GenericTrap( const G4String& name, G4double hz,
63 // << 66 const std::vector<G4TwoVector>& vertices )
64 // Constructor << 67 : G4VSolid(name),
65 // << 68 fpPolyhedron(0),
66 G4GenericTrap::G4GenericTrap(const G4String& n << 69 fDz(hz),
67 const std::vector << 70 fVertices(),
68 : G4VSolid(name) << 71 fIsTwisted(false),
69 { << 72 fTessellatedSolid(0),
70 halfTolerance = 0.5*kCarTolerance; << 73 fMinBBoxVector(G4ThreeVector(0,0,0)),
71 CheckParameters(halfZ, vertices); << 74 fMaxBBoxVector(G4ThreeVector(0,0,0)),
72 ComputeLateralSurfaces(); << 75 fVisSubdivisions(0),
73 ComputeBoundingBox(); << 76 fSurfaceArea(0.),
74 ComputeScratchLength(); << 77 fCubicVolume(0.)
>> 78
>> 79 {
>> 80 // General constructor
>> 81 const G4double min_length=5*1.e-6;
>> 82 G4double length = 0.;
>> 83 G4int k=0;
>> 84 G4String errorDescription = "InvalidSetup in \" ";
>> 85 errorDescription += name;
>> 86 errorDescription += "\"";
>> 87
>> 88 // Check vertices size
>> 89
>> 90 if ( G4int(vertices.size()) != fgkNofVertices )
>> 91 {
>> 92 G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002",
>> 93 FatalErrorInArgument, "Number of vertices != 8");
>> 94 }
>> 95
>> 96 // Check dZ
>> 97 //
>> 98 if (hz < kCarTolerance)
>> 99 {
>> 100 G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002",
>> 101 FatalErrorInArgument, "dZ is too small or negative");
>> 102 }
>> 103
>> 104 // Check Ordering and Copy vertices
>> 105 //
>> 106 if(CheckOrder(vertices))
>> 107 {
>> 108 for (G4int i=0; i<fgkNofVertices; ++i) {fVertices.push_back(vertices[i]);}
>> 109 }
>> 110 else
>> 111 {
>> 112 for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[3-i]);}
>> 113 for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[7-i]);}
>> 114 }
>> 115
>> 116 // Check length of segments and Adjust
>> 117 //
>> 118 for (G4int j=0; j < 2; j++)
>> 119 {
>> 120 for (G4int i=1; i<4; ++i)
>> 121 {
>> 122 k = j*4+i;
>> 123 length = (fVertices[k]-fVertices[k-1]).mag();
>> 124 if ( ( length < min_length) && ( length > kCarTolerance ) )
>> 125 {
>> 126 std::ostringstream message;
>> 127 message << "Length segment is too small." << G4endl
>> 128 << "Distance between " << fVertices[k-1] << " and "
>> 129 << fVertices[k] << " is only " << length << " mm !";
>> 130 G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids1001",
>> 131 JustWarning, message, "Vertices will be collapsed.");
>> 132 fVertices[k]=fVertices[k-1];
>> 133 }
>> 134 }
>> 135 }
>> 136
>> 137 // Compute Twist
>> 138 //
>> 139 for( G4int i=0; i<4; i++) { fTwist[i]=0.; }
>> 140 fIsTwisted = ComputeIsTwisted();
>> 141
>> 142 // Compute Bounding Box
>> 143 //
>> 144 ComputeBBox();
>> 145
>> 146 // If not twisted - create tessellated solid
>> 147 // (an alternative implementation for testing)
>> 148 //
>> 149 #ifdef G4TESS_TEST
>> 150 if ( !fIsTwisted ) { fTessellatedSolid = CreateTessellatedSolid(); }
>> 151 #endif
75 } 152 }
76 153
77 ////////////////////////////////////////////// << 154 // --------------------------------------------------------------------
78 // << 155
79 // Fake default constructor - sets only member << 156 G4GenericTrap::G4GenericTrap( __void__& a )
80 // for usage restri << 157 : G4VSolid(a),
81 G4GenericTrap::G4GenericTrap(__void__& a) << 158 fpPolyhedron(0),
82 : G4VSolid(a) << 159 fDz(0.),
>> 160 fVertices(),
>> 161 fIsTwisted(false),
>> 162 fTessellatedSolid(0),
>> 163 fMinBBoxVector(G4ThreeVector(0,0,0)),
>> 164 fMaxBBoxVector(G4ThreeVector(0,0,0)),
>> 165 fVisSubdivisions(0),
>> 166 fSurfaceArea(0.),
>> 167 fCubicVolume(0.)
83 { 168 {
>> 169 // Fake default constructor - sets only member data and allocates memory
>> 170 // for usage restricted to object persistency.
>> 171
>> 172 for (size_t i=0; i<4; ++i) { fTwist[i]=0.; }
84 } 173 }
85 174
86 ////////////////////////////////////////////// << 175 // --------------------------------------------------------------------
87 // <<
88 // Destructor <<
89 176
90 G4GenericTrap::~G4GenericTrap() 177 G4GenericTrap::~G4GenericTrap()
91 { 178 {
>> 179 // Destructor
>> 180 delete fTessellatedSolid;
92 } 181 }
93 182
94 ////////////////////////////////////////////// << 183 // --------------------------------------------------------------------
95 // << 184
96 // Copy constructor <<
97 // <<
98 G4GenericTrap::G4GenericTrap(const G4GenericTr 185 G4GenericTrap::G4GenericTrap(const G4GenericTrap& rhs)
99 : G4VSolid(rhs), 186 : G4VSolid(rhs),
100 halfTolerance(rhs.halfTolerance), fScratch << 187 fpPolyhedron(0), fDz(rhs.fDz), fVertices(rhs.fVertices),
101 fDz(rhs.fDz), fVertices(rhs.fVertices), fI << 188 fIsTwisted(rhs.fIsTwisted), fTessellatedSolid(0),
102 fMinBBox(rhs.fMinBBox), fMaxBBox(rhs.fMaxB << 189 fMinBBoxVector(rhs.fMinBBoxVector), fMaxBBoxVector(rhs.fMaxBBoxVector),
103 fVisSubdivisions(rhs.fVisSubdivisions), 190 fVisSubdivisions(rhs.fVisSubdivisions),
104 fSurfaceArea(rhs.fSurfaceArea), fCubicVolu << 191 fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume)
105 { 192 {
106 for (auto i = 0; i < 5; ++i) { fTwist[i] = r << 193 for (size_t i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; }
107 for (auto i = 0; i < 4; ++i) { fDelta[i] = r << 194 #ifdef G4TESS_TEST
108 for (auto i = 0; i < 8; ++i) { fPlane[i] = r << 195 if (rhs.fTessellatedSolid && !fIsTwisted )
109 for (auto i = 0; i < 4; ++i) { fSurf[i] = rh << 196 { fTessellatedSolid = CreateTessellatedSolid(); }
110 for (auto i = 0; i < 4; ++i) { fArea[i] = rh << 197 #endif
111 } 198 }
112 199
113 ////////////////////////////////////////////// << 200 // --------------------------------------------------------------------
114 // << 201
115 // Assignment << 202 G4GenericTrap& G4GenericTrap::operator = (const G4GenericTrap& rhs)
116 // <<
117 G4GenericTrap& G4GenericTrap::operator=(const <<
118 { 203 {
119 // Check assignment to self << 204 // Check assignment to self
120 if (this == &rhs) { return *this; } << 205 //
>> 206 if (this == &rhs) { return *this; }
>> 207
>> 208 // Copy base class data
>> 209 //
>> 210 G4VSolid::operator=(rhs);
>> 211
>> 212 // Copy data
>> 213 //
>> 214 fpPolyhedron = 0; fDz = rhs.fDz; fVertices = rhs.fVertices;
>> 215 fIsTwisted = rhs.fIsTwisted; fTessellatedSolid = 0;
>> 216 fMinBBoxVector = rhs.fMinBBoxVector; fMaxBBoxVector = rhs.fMaxBBoxVector;
>> 217 fVisSubdivisions = rhs.fVisSubdivisions;
>> 218 fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
>> 219
>> 220 for (size_t i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; }
>> 221 #ifdef G4TESS_TEST
>> 222 if (rhs.fTessellatedSolid && !fIsTwisted )
>> 223 { delete fTessellatedSolid; fTessellatedSolid = CreateTessellatedSolid(); }
>> 224 #endif
121 225
122 // Copy base class data << 226 return *this;
123 G4VSolid::operator=(rhs); << 227 }
124 228
125 // Copy data << 229 // --------------------------------------------------------------------
126 halfTolerance = rhs.halfTolerance; fScratch <<
127 fDz = rhs.fDz; fVertices = rhs.fVertices; fI <<
128 fMinBBox = rhs.fMinBBox; fMaxBBox = rhs.fMax <<
129 fVisSubdivisions = rhs.fVisSubdivisions; <<
130 fSurfaceArea = rhs.fSurfaceArea; fCubicVolum <<
131 230
132 for (auto i = 0; i < 8; ++i) { fVertices[i] << 231 EInside
133 for (auto i = 0; i < 5; ++i) { fTwist[i] = r << 232 G4GenericTrap::InsidePolygone(const G4ThreeVector& p,
134 for (auto i = 0; i < 4; ++i) { fDelta[i] = r << 233 const std::vector<G4TwoVector>& poly) const
135 for (auto i = 0; i < 8; ++i) { fPlane[i] = r << 234 {
136 for (auto i = 0; i < 4; ++i) { fSurf[i] = rh << 235 static const G4double halfCarTolerance=kCarTolerance*0.5;
137 for (auto i = 0; i < 4; ++i) { fArea[i] = rh << 236 EInside in = kInside;
>> 237 G4double cross, len2;
>> 238 G4int count=0;
138 239
139 fRebuildPolyhedron = false; << 240 for (G4int i = 0; i < 4; i++)
140 delete fpPolyhedron; fpPolyhedron = nullptr; << 241 {
>> 242 G4int j = (i+1) % 4;
141 243
142 return *this; << 244 cross = (p.x()-poly[i].x())*(poly[j].y()-poly[i].y())-
143 } << 245 (p.y()-poly[i].y())*(poly[j].x()-poly[i].x());
144 246
145 ////////////////////////////////////////////// << 247 len2=(poly[i]-poly[j]).mag2();
146 // << 248 if (len2 > kCarTolerance)
147 // Returns position of the point (inside/outsi << 249 {
148 // << 250 if(cross*cross<=len2*halfCarTolerance*halfCarTolerance) // Surface check
149 EInside G4GenericTrap::Inside(const G4ThreeVec << 251 {
150 { << 252 G4double test;
151 G4double px = p.x(), py = p.y(), pz = p.z(); << 253 G4int k,l;
>> 254 if ( poly[i].y() > poly[j].y() ) { k=i; l=j; }
>> 255 else { k=j; l=i; }
152 256
153 // intersect edges by z = pz plane << 257 if ( poly[k].x() != poly[l].x() )
154 G4TwoVector pp[4]; << 258 {
155 G4double z = (pz + fDz); << 259 test = (p.x()-poly[l].x())/(poly[k].x()-poly[l].x())
156 for (auto i = 0; i < 4; ++i) { pp[i] = fVert << 260 * (poly[k].y()-poly[l].y())+poly[l].y();
>> 261 }
>> 262 else
>> 263 {
>> 264 test = p.y();
>> 265 }
157 266
158 // estimate distance to the solid << 267 // Check if point is Inside Segment
159 G4double dist = std::abs(pz) - fDz; << 268 //
160 for (auto i = 0; i < 4; ++i) << 269 if( (test>=(poly[l].y()-halfCarTolerance))
161 { << 270 && (test<=(poly[k].y()+halfCarTolerance)) )
162 if (fTwist[i] == 0.) << 271 {
163 { << 272 return kSurface;
164 G4double dd = fSurf[i].D*px + fSurf[i].E << 273 }
165 dist = std::max(dd, dist); << 274 else
>> 275 {
>> 276 return kOutside;
>> 277 }
>> 278 }
>> 279 else if (cross<0.) { return kOutside; }
166 } 280 }
167 else 281 else
168 { 282 {
169 auto j = (i + 1)%4; << 283 count++;
170 G4TwoVector a = pp[i]; <<
171 G4TwoVector b = pp[j]; <<
172 G4double dx = b.x() - a.x(); <<
173 G4double dy = b.y() - a.y(); <<
174 G4double leng = std::sqrt(dx*dx + dy*dy) <<
175 G4double dd = (dx*(py - a.y()) - dy*(px <<
176 dist = std::max(dd, dist); <<
177 } 284 }
178 } 285 }
179 return (dist > halfTolerance) ? kOutside : <<
180 ((dist > -halfTolerance) ? kSurface : kIns <<
181 } <<
182 286
183 ////////////////////////////////////////////// << 287 // All collapsed vertices, Tet like
184 // << 288 //
185 // Return unit normal to surface at p << 289 if(count==4)
186 // << 290 {
187 G4ThreeVector G4GenericTrap::SurfaceNormal( co << 291 if ( (std::fabs(p.x()-poly[0].x())+std::fabs(p.y()-poly[0].y())) > halfCarTolerance )
188 { <<
189 G4double halfToleranceSquared = halfToleranc <<
190 G4double px = p.x(), py = p.y(), pz = p.z(); <<
191 G4double nx = 0, ny = 0, nz = 0; <<
192 <<
193 // intersect edges by z = pz plane <<
194 G4TwoVector pp[4]; <<
195 G4double tz = (pz + fDz); <<
196 for (auto i = 0; i < 4; ++i) { pp[i] = fVert <<
197 <<
198 // check bottom and top faces <<
199 G4double dz = std::abs(pz) - fDz; <<
200 nz = std::copysign(G4double(std::abs(dz) <= <<
201 <<
202 // check lateral faces <<
203 for (auto i = 0; i < 4; ++i) <<
204 { <<
205 if (fTwist[i] == 0.) <<
206 { <<
207 G4double dd = fSurf[i].D*px + fSurf[i].E <<
208 if (std::abs(dd) <= halfTolerance) <<
209 { <<
210 nx += fSurf[i].D; <<
211 ny += fSurf[i].E; <<
212 nz += fSurf[i].F; <<
213 } <<
214 } <<
215 else <<
216 { 292 {
217 auto j = (i + 1)%4; << 293 in=kOutside;
218 G4TwoVector a = pp[i]; <<
219 G4TwoVector b = pp[j]; <<
220 G4double dx = b.x() - a.x(); <<
221 G4double dy = b.y() - a.y(); <<
222 G4double ll = dx*dx + dy*dy; <<
223 G4double dd = dx*(py - a.y()) - dy*(px - <<
224 if (dd*dd <= halfToleranceSquared*ll) <<
225 { <<
226 G4double x = fSurf[i].A*pz + fSurf[i]. <<
227 G4double y = fSurf[i].B*pz + fSurf[i]. <<
228 G4double z = fSurf[i].A*px + fSurf[i]. <<
229 G4double inv = 1./std::sqrt(x*x + y*y <<
230 nx += x*inv; <<
231 ny += y*inv; <<
232 nz += z*inv; <<
233 } <<
234 } 294 }
235 } 295 }
236 << 296 return in;
237 // return normal <<
238 G4double mag2 = nx*nx + ny*ny + nz*nz; <<
239 if (mag2 == 0.) return ApproxSurfaceNormal(p <<
240 G4double mag = std::sqrt(mag2); <<
241 G4double inv = (mag == 1.) ? 1. : 1./mag; <<
242 return { nx*inv, ny*inv, nz*inv }; <<
243 } 297 }
244 298
245 ////////////////////////////////////////////// << 299 // --------------------------------------------------------------------
246 // << 300
247 // Find surface nearest to the point and retur << 301 EInside G4GenericTrap::Inside(const G4ThreeVector& p) const
248 // Normally this method should not be called <<
249 // <<
250 G4ThreeVector G4GenericTrap::ApproxSurfaceNorm <<
251 { 302 {
252 #ifdef G4SPECSDEBUG << 303 // Test if point is inside this shape
253 std::ostringstream message; <<
254 message.precision(16); <<
255 message << "Point p is not on surface of sol <<
256 << "Position: " << p << " is " <<
257 << ((Inside(p) == kInside) ? "inside <<
258 StreamInfo(message); <<
259 G4Exception("G4GenericTrap::ApproxSurfaceNor <<
260 JustWarning, message ); <<
261 #endif <<
262 G4double px = p.x(), py = p.y(), pz = p.z(); <<
263 G4double dist = -DBL_MAX; <<
264 auto iside = 0; <<
265 <<
266 // intersect edges by z = pz plane <<
267 G4TwoVector pp[4]; <<
268 G4double tz = (pz + fDz); <<
269 for (auto i = 0; i < 4; ++i) { pp[i] = fVert <<
270 304
271 // check lateral faces << 305 #ifdef G4TESS_TEST
272 for (auto i = 0; i < 4; ++i) << 306 if ( fTessellatedSolid )
>> 307 {
>> 308 return fTessellatedSolid->Inside(p);
>> 309 }
>> 310 #endif
>> 311
>> 312 static const G4double halfCarTolerance=kCarTolerance*0.5;
>> 313 EInside innew=kOutside;
>> 314 std::vector<G4TwoVector> xy;
>> 315
>> 316 if (std::fabs(p.z()) <= fDz+halfCarTolerance) // First check Z range
273 { 317 {
274 if (fTwist[i] == 0.) << 318 // Compute intersection between Z plane containing point and the shape
>> 319 //
>> 320 G4double cf = 0.5*(fDz-p.z())/fDz;
>> 321 for (G4int i=0; i<4; i++)
275 { 322 {
276 G4double d = fSurf[i].D*px + fSurf[i].E* << 323 xy.push_back(fVertices[i+4]+cf*( fVertices[i]-fVertices[i+4]));
277 if (d > dist) { dist = d; iside = i; } <<
278 } 324 }
279 else <<
280 { <<
281 auto j = (i + 1)%4; <<
282 G4TwoVector a = pp[i]; <<
283 G4TwoVector b = pp[j]; <<
284 G4double dx = b.x() - a.x(); <<
285 G4double dy = b.y() - a.y(); <<
286 G4double leng = std::sqrt(dx*dx + dy*dy) <<
287 G4double d = (dx*(py - a.y()) - dy*(px - <<
288 if (d > dist) { dist = d; iside = i; } <<
289 } <<
290 } <<
291 // check bottom and top faces <<
292 G4double distz = std::abs(pz) - fDz; <<
293 if (distz >= dist) return { 0., 0., std::cop <<
294 <<
295 G4double x = fSurf[iside].A*pz + fSurf[iside <<
296 G4double y = fSurf[iside].B*pz + fSurf[iside <<
297 G4double z = fSurf[iside].A*px + fSurf[iside <<
298 G4double inv = 1./std::sqrt(x*x + y*y + z*z) <<
299 return { x*inv, y*inv, z*inv }; <<
300 } <<
301 325
302 ////////////////////////////////////////////// << 326 innew=InsidePolygone(p,xy);
303 // << 327
304 // Compute distance to the surface from outsid << 328 if( (innew==kInside) || (innew==kSurface) )
305 // return kInfinity if no intersection << 329 {
306 // << 330 if(std::fabs(p.z()) > fDz-halfCarTolerance) { innew=kSurface; }
307 G4double G4GenericTrap::DistanceToIn(const G4T << 331 }
308 const G4T << 332 }
>> 333 return innew;
>> 334 }
>> 335
>> 336 // --------------------------------------------------------------------
>> 337
>> 338 G4ThreeVector G4GenericTrap::SurfaceNormal( const G4ThreeVector& p ) const
309 { 339 {
310 G4double px = p.x(), py = p.y(), pz = p.z(); << 340 // Calculate side nearest to p, and return normal
311 G4double vx = v.x(), vy = v.y(), vz = v.z(); << 341 // If two sides are equidistant, sum of the Normal is returned
312 342
313 // Find intersections with the bounding poly << 343 #ifdef G4TESS_TEST
314 // << 344 if ( fTessellatedSolid )
315 if (std::abs(pz) - fDz >= -halfTolerance && << 345 {
316 G4double invz = (vz == 0) ? kInfinity : -1./ << 346 return fTessellatedSolid->SurfaceNormal(p);
317 G4double dz = std::copysign(fDz,invz); << 347 }
318 G4double xin = (pz - dz)*invz; << 348 #endif
319 G4double xout = (pz + dz)*invz; << 349
320 << 350 static const G4double halfCarTolerance=kCarTolerance*0.5;
321 // Check plane faces << 351
322 for (auto k = 0; k < 4; ++k) << 352 G4ThreeVector lnorm, sumnorm(0.,0.,0.), apprnorm(0.,0.,1.),
323 { << 353 p0, p1, p2, r1, r2, r3, r4;
324 if (fTwist[k] != 0) continue; // skip twis << 354 G4int noSurfaces = 0;
325 const G4GenericTrapPlane& surf = fPlane[2* << 355 G4double distxy,distz;
326 G4double cosa = surf.A*vx + surf.B*vy + su << 356 G4bool zPlusSide=false;
327 G4double dist = surf.A*px + surf.B*py + su << 357
328 if (dist >= -halfTolerance) << 358 distz = fDz-std::fabs(p.z());
329 { << 359 if (distz < halfCarTolerance)
330 if (cosa >= 0.) { return kInfinity; } // << 360 {
331 G4double tmp = -dist/cosa; << 361 if(p.z()>0)
332 xin = std::max(tmp, xin); << 362 {
>> 363 zPlusSide=true;
>> 364 sumnorm=G4ThreeVector(0,0,1);
333 } 365 }
334 else 366 else
335 { 367 {
336 if (cosa > 0) { xout = std::min(-dist/co << 368 sumnorm=G4ThreeVector(0,0,-1);
337 if (cosa < 0) { xin = std::max(-dist/cos <<
338 } 369 }
339 } << 370 noSurfaces ++;
>> 371 }
340 372
341 // Check planes around twisted faces, each t << 373 // Check lateral planes
342 G4double tin = xin; << 374 //
343 G4double tout = xout; << 375 std:: vector<G4TwoVector> vertices;
344 for (auto i = 0; i < 4; ++i) << 376 G4double cf = 0.5*(fDz-p.z())/fDz;
>> 377 for (G4int i=0; i<4; i++)
345 { 378 {
346 if (fTwist[i] == 0) continue; // skip plan << 379 vertices.push_back(fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]));
>> 380 }
347 381
348 // check intersection with 1st bounding pl << 382 // Compute distance for lateral planes
349 const G4GenericTrapPlane& surf1 = fPlane[2 << 383 //
350 G4double cosa = surf1.A*vx + surf1.B*vy + << 384 for (G4int s=0; s<4; s++)
351 G4double dist = surf1.A*px + surf1.B*py + << 385 {
352 if (dist >= -halfTolerance) << 386 p0=G4ThreeVector(vertices[s].x(),vertices[s].y(),p.z());
>> 387 if(zPlusSide)
353 { 388 {
354 if (cosa >= 0.) { return kInfinity; } // << 389 p1=G4ThreeVector(fVertices[s].x(),fVertices[s].y(),-fDz);
355 G4double tmp = -dist/cosa; <<
356 tin = std::max(tmp, tin); <<
357 } 390 }
358 else 391 else
359 { 392 {
360 if (cosa > 0) { tout = std::min(-dist/co << 393 p1=G4ThreeVector(fVertices[s+4].x(),fVertices[s+4].y(),fDz);
361 if (cosa < 0) { tin = std::max(-dist/cos <<
362 } 394 }
>> 395 p2=G4ThreeVector(vertices[(s+1)%4].x(),vertices[(s+1)%4].y(),p.z());
363 396
364 // check intersection with 2nd bounding pl << 397 // Collapsed vertices
365 const G4GenericTrapPlane& surf2 = fPlane[2 << 398 //
366 cosa = surf2.A*vx + surf2.B*vy + surf2.C*v << 399 if ( (p2-p0).mag2() < kCarTolerance )
367 dist = surf2.A*px + surf2.B*py + surf2.C*p << 400 {
368 if (dist >= -halfTolerance) << 401 if ( std::fabs(p.z()+fDz) > kCarTolerance )
369 { << 402 {
370 if (cosa >= 0.) { return kInfinity; } // << 403 p2=G4ThreeVector(fVertices[(s+1)%4].x(),fVertices[(s+1)%4].y(),-fDz);
371 G4double tmp = -dist/cosa; << 404 }
372 tin = std::max(tmp, tin); << 405 else
>> 406 {
>> 407 p2=G4ThreeVector(fVertices[(s+1)%4+4].x(),fVertices[(s+1)%4+4].y(),fDz);
>> 408 }
373 } 409 }
374 else << 410 lnorm = (p1-p0).cross(p2-p0);
>> 411 lnorm = lnorm.unit();
>> 412 if(zPlusSide) { lnorm=-lnorm; }
>> 413
>> 414 // Adjust Normal for Twisted Surface
>> 415 //
>> 416 if ( (fIsTwisted) && (GetTwistAngle(s)!=0) )
>> 417 {
>> 418 G4double normP=(p2-p0).mag();
>> 419 if(normP)
>> 420 {
>> 421 G4double proj=(p-p0).dot(p2-p0)/normP;
>> 422 if(proj<0) { proj=0; }
>> 423 if(proj>normP) { proj=normP; }
>> 424 G4int j=(s+1)%4;
>> 425 r1=G4ThreeVector(fVertices[s+4].x(),fVertices[s+4].y(),fDz);
>> 426 r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
>> 427 r3=G4ThreeVector(fVertices[s].x(),fVertices[s].y(),-fDz);
>> 428 r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
>> 429 r1=r1+proj*(r2-r1)/normP;
>> 430 r3=r3+proj*(r4-r3)/normP;
>> 431 r2=r1-r3;
>> 432 r4=r2.cross(p2-p0); r4=r4.unit();
>> 433 lnorm=r4;
>> 434 }
>> 435 } // End if fIsTwisted
>> 436
>> 437 distxy=std::fabs((p0-p).dot(lnorm));
>> 438 if ( distxy<halfCarTolerance )
375 { 439 {
376 if (cosa > 0) { tout = std::min(-dist/co << 440 noSurfaces ++;
377 if (cosa < 0) { tin = std::max(-dist/cos << 441
>> 442 // Negative sign for Normal is taken for Outside Normal
>> 443 //
>> 444 sumnorm=sumnorm+lnorm;
378 } 445 }
>> 446
>> 447 // For ApproxSurfaceNormal
>> 448 //
>> 449 if (distxy<distz)
>> 450 {
>> 451 distz=distxy;
>> 452 apprnorm=lnorm;
>> 453 }
>> 454 } // End for loop
>> 455
>> 456 // Calculate final Normal, add Normal in the Corners and Touching Sides
>> 457 //
>> 458 if ( noSurfaces == 0 )
>> 459 {
>> 460 G4Exception("G4GenericTrap::SurfaceNormal(p)", "GeomSolids1002",
>> 461 JustWarning, "Point p is not on surface !?" );
>> 462 sumnorm=apprnorm;
>> 463 // Add Approximative Surface Normal Calculation?
379 } 464 }
380 if (tout - tin <= halfTolerance) { return kI << 465 else if ( noSurfaces == 1 ) { sumnorm = sumnorm; }
>> 466 else { sumnorm = sumnorm.unit(); }
381 467
382 // if point is outside of the bounding box << 468 return sumnorm ;
383 // then move it to the surface of the boundi << 469 }
384 G4double t0 = 0., x0 = px, y0 = py, z0 = pz; <<
385 if (x0 < fMinBBox.x() - halfTolerance || <<
386 y0 < fMinBBox.y() - halfTolerance || <<
387 z0 < fMinBBox.z() - halfTolerance || <<
388 x0 > fMaxBBox.x() + halfTolerance || <<
389 y0 > fMaxBBox.y() + halfTolerance || <<
390 z0 > fMaxBBox.z() + halfTolerance) <<
391 { <<
392 t0 = tin; <<
393 x0 += vx*t0; <<
394 y0 += vy*t0; <<
395 z0 += vz*t0; <<
396 } <<
397 470
398 // Check intersections with twisted faces << 471 // --------------------------------------------------------------------
>> 472
>> 473 G4ThreeVector G4GenericTrap::NormalToPlane( const G4ThreeVector& p,
>> 474 const G4int ipl ) const
>> 475 {
>> 476 // Return normal to given lateral plane ipl
>> 477
>> 478 #ifdef G4TESS_TEST
>> 479 if ( fTessellatedSolid )
>> 480 {
>> 481 return fTessellatedSolid->SurfaceNormal(p);
>> 482 }
>> 483 #endif
>> 484
>> 485 static const G4double halfCarTolerance=kCarTolerance*0.5;
>> 486 G4ThreeVector lnorm, norm(0.,0.,0.), p0,p1,p2;
>> 487
>> 488 G4double distz = fDz-p.z();
>> 489 G4int i=ipl; // current plane index
>> 490
>> 491 G4TwoVector u,v;
>> 492 G4ThreeVector r1,r2,r3,r4;
>> 493 G4double cf = 0.5*(fDz-p.z())/fDz;
>> 494 G4int j=(i+1)%4;
>> 495
>> 496 u=fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]);
>> 497 v=fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]);
>> 498
>> 499 // Compute cross product
399 // 500 //
400 G4double ttin[2] = { DBL_MAX, DBL_MAX }; << 501 p0=G4ThreeVector(u.x(),u.y(),p.z());
401 G4double ttout[2] = { tout, tout }; << 502
>> 503 if (std::fabs(distz)<halfCarTolerance)
>> 504 {
>> 505 p1=G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz);distz=-1;}
>> 506 else
>> 507 {
>> 508 p1=G4ThreeVector(fVertices[i+4].x(),fVertices[i+4].y(),fDz);
>> 509 }
>> 510 p2=G4ThreeVector(v.x(),v.y(),p.z());
402 511
403 if (tin - xin < halfTolerance) ttin[0] = xin << 512 // Collapsed vertices
404 if (xout - tout < halfTolerance) ttout[0] = << 513 //
405 G4double tminimal = tin - halfTolerance; << 514 if ( (p2-p0).mag2() < kCarTolerance )
406 G4double tmaximal = tout + halfTolerance; << 515 {
407 << 516 if ( std::fabs(p.z()+fDz) > halfCarTolerance )
408 constexpr G4double EPSILON = 1000.*DBL_EPSIL << 517 {
409 for (auto i = 0; i < 4; ++i) << 518 p2=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
410 { <<
411 if (fTwist[i] == 0) continue; // skip plan <<
412 <<
413 // twisted face, solve quadratic equation <<
414 G4double ABC = fSurf[i].A*vx + fSurf[i].B <<
415 G4double ABCF = fSurf[i].A*x0 + fSurf[i].B <<
416 G4double A = ABC*vz; <<
417 G4double B = 0.5*(fSurf[i].D*vx + fSurf[i] <<
418 G4double C = fSurf[i].D*x0 + fSurf[i].E*y0 <<
419 if (std::abs(A) <= EPSILON) <<
420 { <<
421 // case 1: track is parallel to the surf <<
422 if (std::abs(B) <= EPSILON) <<
423 { <<
424 // check position of the track relativ <<
425 // for the reason of precision it's be <<
426 auto j = (i + 1)%4; <<
427 G4double z = (z0 + fDz); <<
428 G4TwoVector a = fVertices[i] + fDelta[ <<
429 G4TwoVector b = fVertices[j] + fDelta[ <<
430 G4double dx = b.x() - a.x(); <<
431 G4double dy = b.y() - a.y(); <<
432 G4double leng = std::sqrt(dx*dx + dy*d <<
433 G4double dist = dx*(y0 - a.y()) - dy*( <<
434 if (dist >= -halfTolerance*leng) { ret <<
435 continue; <<
436 } <<
437 <<
438 // case 2: single root <<
439 G4double tmp = t0 - 0.5*C/B; <<
440 // compute normal at the intersection po <<
441 G4double x = px + vx*tmp; <<
442 G4double y = py + vy*tmp; <<
443 G4double z = pz + vz*tmp; <<
444 const G4GenericTrapSurface& surf = fSurf <<
445 G4double nx = surf.A*z + surf.D; <<
446 G4double ny = surf.B*z + surf.E; <<
447 G4double nz = surf.A*x + surf.B*y + 2.*s <<
448 <<
449 if (nx*vx + ny*vy + nz*vz >= 0.) // poin <<
450 { <<
451 auto k = (i == 3) ? 0 : i + 1; <<
452 G4double tz = (pz + fDz); <<
453 G4TwoVector a = fVertices[i] + fDelta[ <<
454 G4TwoVector b = fVertices[k] + fDelta[ <<
455 G4double dx = b.x() - a.x(); <<
456 G4double dy = b.y() - a.y(); <<
457 G4double leng = std::sqrt(dx*dx + dy*d <<
458 G4double dist = dx*(py - a.y()) - dy*( <<
459 if (dist >= -halfTolerance*leng) { ret <<
460 <<
461 if (tmp < tminimal || tmp > tmaximal) <<
462 if (std::abs(tmp - ttout[0]) < halfTol <<
463 if (tmp < ttout[0]) <<
464 { <<
465 ttout[1] = ttout[0]; <<
466 ttout[0] = tmp; <<
467 } <<
468 else { ttout[1] = std::min(tmp, ttout[ <<
469 } <<
470 else // point is flying to inside <<
471 { <<
472 if (tmp < tminimal || tmp > tmaximal) <<
473 if (std::abs(tmp - ttin[0]) < halfTole <<
474 if (tmp < ttin[0]) <<
475 { <<
476 ttin[1] = ttin[0]; <<
477 ttin[0] = tmp; <<
478 } <<
479 else { ttin[1] = std::min(tmp, ttin[1] <<
480 } <<
481 continue; <<
482 } 519 }
483 << 520 else
484 // case 3: scratch or no intersection <<
485 G4double D = B*B - A*C; <<
486 if (D < 0.25*fScratch*fScratch*A*A) <<
487 { 521 {
488 if (A > 0) return kInfinity; << 522 p2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
489 continue; <<
490 } 523 }
>> 524 }
>> 525 lnorm=-(p1-p0).cross(p2-p0);
>> 526 if (distz>-halfCarTolerance) { lnorm=-lnorm.unit(); }
>> 527 else { lnorm=lnorm.unit(); }
>> 528
>> 529 // Adjust Normal for Twisted Surface
>> 530 //
>> 531 if( (fIsTwisted) && (GetTwistAngle(ipl)!=0) )
>> 532 {
>> 533 G4double normP=(p2-p0).mag();
>> 534 if(normP)
>> 535 {
>> 536 G4double proj=(p-p0).dot(p2-p0)/normP;
>> 537 if (proj<0) { proj=0; }
>> 538 if (proj>normP) { proj=normP; }
>> 539
>> 540 r1=G4ThreeVector(fVertices[i+4].x(),fVertices[i+4].y(),fDz);
>> 541 r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
>> 542 r3=G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz);
>> 543 r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
>> 544 r1=r1+proj*(r2-r1)/normP;
>> 545 r3=r3+proj*(r4-r3)/normP;
>> 546 r2=r1-r3;
>> 547 r4=r2.cross(p2-p0);r4=r4.unit();
>> 548 lnorm=r4;
>> 549 }
>> 550 } // End if fIsTwisted
491 551
492 // case 4: two intersection points << 552 return lnorm;
493 G4double tmp = -B - std::copysign(std::sqr << 553 }
494 G4double t1 = tmp/A + t0; <<
495 G4double t2 = C/tmp + t0; <<
496 G4double tsurfin = std::min(t1, t2); <<
497 G4double tsurfout = std::max(t1, t2); <<
498 if (A < 0) std::swap(tsurfin, tsurfout); <<
499 554
500 if (tsurfin >= tminimal && tsurfin <= tmax << 555 // --------------------------------------------------------------------
>> 556
>> 557 G4double G4GenericTrap::DistToPlane(const G4ThreeVector& p,
>> 558 const G4ThreeVector& v,
>> 559 const G4int ipl) const
>> 560 {
>> 561 // Computes distance to plane ipl :
>> 562 // ipl=0 : points 0,4,1,5
>> 563 // ipl=1 : points 1,5,2,6
>> 564 // ipl=2 : points 2,6,3,7
>> 565 // ipl=3 : points 3,7,0,4
>> 566
>> 567 static const G4double halfCarTolerance=0.5*kCarTolerance;
>> 568 G4double xa,xb,xc,xd,ya,yb,yc,yd;
>> 569
>> 570 G4int j = (ipl+1)%4;
>> 571
>> 572 xa=fVertices[ipl].x();
>> 573 ya=fVertices[ipl].y();
>> 574 xb=fVertices[ipl+4].x();
>> 575 yb=fVertices[ipl+4].y();
>> 576 xc=fVertices[j].x();
>> 577 yc=fVertices[j].y();
>> 578 xd=fVertices[4+j].x();
>> 579 yd=fVertices[4+j].y();
>> 580
>> 581 G4double dz2 =0.5/fDz;
>> 582 G4double tx1 =dz2*(xb-xa);
>> 583 G4double ty1 =dz2*(yb-ya);
>> 584 G4double tx2 =dz2*(xd-xc);
>> 585 G4double ty2 =dz2*(yd-yc);
>> 586 G4double dzp =fDz+p.z();
>> 587 G4double xs1 =xa+tx1*dzp;
>> 588 G4double ys1 =ya+ty1*dzp;
>> 589 G4double xs2 =xc+tx2*dzp;
>> 590 G4double ys2 =yc+ty2*dzp;
>> 591 G4double dxs =xs2-xs1;
>> 592 G4double dys =ys2-ys1;
>> 593 G4double dtx =tx2-tx1;
>> 594 G4double dty =ty2-ty1;
>> 595
>> 596 G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
>> 597 G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
>> 598 + tx1*ys2-tx2*ys1)*v.z();
>> 599 G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
>> 600 G4double s=kInfinity;
>> 601 G4double x1,x2,y1,y2,xp,yp,zi;
>> 602
>> 603 if (std::fabs(a)<kCarTolerance)
>> 604 {
>> 605 if (std::fabs(b)<kCarTolerance) { return kInfinity; }
>> 606 s=-c/b;
>> 607
>> 608 // Check if Point is on the Surface
>> 609
>> 610 if (s>-halfCarTolerance)
>> 611 {
>> 612 if (s<halfCarTolerance)
>> 613 {
>> 614 if (NormalToPlane(p,ipl).dot(v)<=0)
>> 615 { if(Inside(p) != kOutside) { return 0.; } }
>> 616 else
>> 617 { return kInfinity; }
>> 618 }
>> 619
>> 620 // Check the Intersection
>> 621 //
>> 622 zi=p.z()+s*v.z();
>> 623 if (std::fabs(zi)<fDz)
>> 624 {
>> 625 x1=xs1+tx1*v.z()*s;
>> 626 x2=xs2+tx2*v.z()*s;
>> 627 xp=p.x()+s*v.x();
>> 628 y1=ys1+ty1*v.z()*s;
>> 629 y2=ys2+ty2*v.z()*s;
>> 630 yp=p.y()+s*v.y();
>> 631 zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
>> 632 if (zi<=halfCarTolerance) { return s; }
>> 633 }
>> 634 }
>> 635 return kInfinity;
>> 636 }
>> 637 G4double d=b*b-4*a*c;
>> 638 if (d>=0)
>> 639 {
>> 640 if (a>0) { s=0.5*(-b-std::sqrt(d))/a; }
>> 641 else { s=0.5*(-b+std::sqrt(d))/a; }
>> 642
>> 643 // Check if Point is on the Surface
>> 644 //
>> 645 if (s>-halfCarTolerance)
501 { 646 {
502 if (std::abs(tsurfin - ttin[0]) >= halfT << 647 if(s<halfCarTolerance)
503 { 648 {
504 if (tsurfin < ttin[0]) << 649 if (NormalToPlane(p,ipl).dot(v)<=0)
>> 650 {
>> 651 if(Inside(p)!= kOutside) { return 0.; }
>> 652 }
>> 653 else // Check second root; return kInfinity
505 { 654 {
506 ttin[1] = ttin[0]; << 655 if (a>0) { s=0.5*(-b+std::sqrt(d))/a; }
507 ttin[0] = tsurfin; << 656 else { s=0.5*(-b-std::sqrt(d))/a; }
>> 657 if (s<=halfCarTolerance) { return kInfinity; }
508 } 658 }
509 else { ttin[1] = std::min(tsurfin, tti << 659 }
>> 660 // Check the Intersection
>> 661 //
>> 662 zi=p.z()+s*v.z();
>> 663 if (std::fabs(zi)<fDz)
>> 664 {
>> 665 x1=xs1+tx1*v.z()*s;
>> 666 x2=xs2+tx2*v.z()*s;
>> 667 xp=p.x()+s*v.x();
>> 668 y1=ys1+ty1*v.z()*s;
>> 669 y2=ys2+ty2*v.z()*s;
>> 670 yp=p.y()+s*v.y();
>> 671 zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
>> 672 if (zi<=halfCarTolerance) { return s; }
510 } 673 }
511 } 674 }
512 if (tsurfout >= tminimal && tsurfout <= tm << 675 if (a>0) { s=0.5*(-b+std::sqrt(d))/a; }
>> 676 else { s=0.5*(-b-std::sqrt(d))/a; }
>> 677
>> 678 // Check if Point is on the Surface
>> 679 //
>> 680 if (s>-halfCarTolerance)
513 { 681 {
514 if (std::abs(tsurfout - ttout[0]) >= hal << 682 if(s<halfCarTolerance)
515 { 683 {
516 if (tsurfout < ttout[0]) << 684 if (NormalToPlane(p,ipl).dot(v)<=0)
517 { 685 {
518 ttout[1] = ttout[0]; << 686 if(Inside(p) != kOutside) { return 0.; }
519 ttout[0] = tsurfout; <<
520 } 687 }
521 else { ttout[1] = std::min(tsurfout, t << 688 else // Check second root; return kInfinity.
>> 689 {
>> 690 if (a>0) { s=0.5*(-b-std::sqrt(d))/a; }
>> 691 else { s=0.5*(-b+std::sqrt(d))/a; }
>> 692 if (s<=halfCarTolerance) { return kInfinity; }
>> 693 }
>> 694 }
>> 695 // Check the Intersection
>> 696 //
>> 697 zi=p.z()+s*v.z();
>> 698 if (std::fabs(zi)<fDz)
>> 699 {
>> 700 x1=xs1+tx1*v.z()*s;
>> 701 x2=xs2+tx2*v.z()*s;
>> 702 xp=p.x()+s*v.x();
>> 703 y1=ys1+ty1*v.z()*s;
>> 704 y2=ys2+ty2*v.z()*s;
>> 705 yp=p.y()+s*v.y();
>> 706 zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
>> 707 if (zi<=halfCarTolerance) { return s; }
522 } 708 }
523 } 709 }
524 } 710 }
>> 711 return kInfinity;
>> 712 }
525 713
526 // Compute distance to In << 714 // --------------------------------------------------------------------
>> 715
>> 716 G4double G4GenericTrap::DistanceToIn(const G4ThreeVector& p,
>> 717 const G4ThreeVector& v) const
>> 718 {
>> 719 #ifdef G4TESS_TEST
>> 720 if ( fTessellatedSolid )
>> 721 {
>> 722 return fTessellatedSolid->DistanceToIn(p, v);
>> 723 }
>> 724 #endif
>> 725
>> 726 static const G4double halfCarTolerance=kCarTolerance*0.5;
>> 727
>> 728 G4double dist[5];
>> 729 G4ThreeVector n;
>> 730
>> 731 // Check lateral faces
527 // 732 //
528 if (ttin[0] == DBL_MAX) { return kInfinity; << 733 G4int i;
>> 734 for (i=0; i<4; i++)
>> 735 {
>> 736 dist[i]=DistToPlane(p, v, i);
>> 737 }
529 738
530 // single entry point << 739 // Check Z planes
531 if (ttin[1] == DBL_MAX) << 740 //
>> 741 dist[4]=kInfinity;
>> 742 if (std::fabs(p.z())>fDz-halfCarTolerance)
532 { 743 {
533 G4double distin = ttin[0]; << 744 if (v.z())
534 G4double distout = (distin >= ttout[0] - h << 745 {
535 G4double dist = (distout <= halfTolerance << 746 G4ThreeVector pt;
536 return (dist < halfTolerance) ? 0. : dist; << 747 if (p.z()>0)
>> 748 {
>> 749 dist[4] = (fDz-p.z())/v.z();
>> 750 }
>> 751 else
>> 752 {
>> 753 dist[4] = (-fDz-p.z())/v.z();
>> 754 }
>> 755 if (dist[4]<-halfCarTolerance)
>> 756 {
>> 757 dist[4]=kInfinity;
>> 758 }
>> 759 else
>> 760 {
>> 761 if(dist[4]<halfCarTolerance)
>> 762 {
>> 763 if(p.z()>0) { n=G4ThreeVector(0,0,1); }
>> 764 else { n=G4ThreeVector(0,0,-1); }
>> 765 if (n.dot(v)<0) { dist[4]=0.; }
>> 766 else { dist[4]=kInfinity; }
>> 767 }
>> 768 pt=p+dist[4]*v;
>> 769 if (Inside(pt)==kOutside) { dist[4]=kInfinity; }
>> 770 }
>> 771 }
>> 772 }
>> 773 G4double distmin = dist[0];
>> 774 for (i=1;i<5;i++)
>> 775 {
>> 776 if (dist[i] < distmin) { distmin = dist[i]; }
537 } 777 }
538 778
539 // two entry points << 779 if (distmin<halfCarTolerance) { distmin=0.; }
540 if (ttin[1] < ttout[0]) << 780
>> 781 return distmin;
>> 782 }
>> 783
>> 784 // --------------------------------------------------------------------
>> 785
>> 786 G4double G4GenericTrap::DistanceToIn(const G4ThreeVector& p) const
>> 787 {
>> 788 // Computes the closest distance from given point to this shape
>> 789
>> 790 #ifdef G4TESS_TEST
>> 791 if ( fTessellatedSolid )
541 { 792 {
542 G4double distin = ttin[1], distout = ttout << 793 return fTessellatedSolid->DistanceToIn(p);
543 G4double dist = (distout <= halfTolerance << 794 }
544 return (dist < halfTolerance) ? 0. : dist; << 795 #endif
>> 796
>> 797 G4double safz = std::fabs(p.z())-fDz;
>> 798 if(safz<0) { safz=0; }
>> 799
>> 800 G4int iseg;
>> 801 G4double safe = safz;
>> 802 G4double safxy = safz;
>> 803
>> 804 for (iseg=0; iseg<4; iseg++)
>> 805 {
>> 806 safxy = SafetyToFace(p,iseg);
>> 807 if (safxy>safe) { safe=safxy; }
545 } 808 }
546 809
547 // check 1st pair of in-out points << 810 return safe;
548 G4double distin1 = ttin[0], distout1 = ttout << 811 }
549 G4double dist1 = (distout1 <= halfTolerance << 812
550 if (dist1 != kInfinity) { return (dist1 < ha << 813 // --------------------------------------------------------------------
>> 814
>> 815 G4double
>> 816 G4GenericTrap::SafetyToFace(const G4ThreeVector& p, const G4int iseg) const
>> 817 {
>> 818 // Estimate distance to lateral plane defined by segment iseg in range [0,3]
>> 819 // Might be negative: plane seen only from inside
>> 820
>> 821 G4ThreeVector p1,norm;
>> 822 G4double safe;
551 823
552 // check 2nd pair of in-out points << 824 p1=G4ThreeVector(fVertices[iseg].x(),fVertices[iseg].y(),-fDz);
553 G4double distin2 = ttin[1], distout2 = ttout << 825
554 G4double dist2 = (distout2 <= halfTolerance << 826 norm=NormalToPlane(p,iseg);
555 return (dist2 < halfTolerance) ? 0. : dist2; << 827 safe = (p-p1).dot(norm); // Can be negative
>> 828
>> 829 return safe;
556 } 830 }
557 831
558 ////////////////////////////////////////////// << 832 // --------------------------------------------------------------------
559 // << 833
560 // Estimate safety distance to the surface fro << 834 G4double
561 // << 835 G4GenericTrap::DistToTriangle(const G4ThreeVector& p,
562 G4double G4GenericTrap::DistanceToIn(const G4T << 836 const G4ThreeVector& v, const G4int ipl) const
563 { 837 {
564 G4double px = p.x(), py = p.y(), pz = p.z(); << 838 static const G4double halfCarTolerance=kCarTolerance*0.5;
565 839
566 // intersect edges by z = pz plane << 840 G4double xa=fVertices[ipl].x();
567 G4TwoVector pp[4]; << 841 G4double ya=fVertices[ipl].y();
568 G4double z = (pz + fDz); << 842 G4double xb=fVertices[ipl+4].x();
569 for (auto i = 0; i < 4; ++i) { pp[i] = fVert << 843 G4double yb=fVertices[ipl+4].y();
>> 844 G4int j=(ipl+1)%4;
>> 845 G4double xc=fVertices[j].x();
>> 846 G4double yc=fVertices[j].y();
>> 847 G4double zab=2*fDz;
>> 848 G4double zac=0;
>> 849
>> 850 if ( (std::fabs(xa-xc)+std::fabs(ya-yc)) < halfCarTolerance )
>> 851 {
>> 852 xc=fVertices[j+4].x();
>> 853 yc=fVertices[j+4].y();
>> 854 zac=2*fDz;
>> 855 zab=2*fDz;
570 856
571 // estimate distance to the solid << 857 //Line case
572 G4double dist = std::abs(pz) - fDz; << 858 //
573 for (auto i = 0; i < 4; ++i) << 859 if ( (std::fabs(xb-xc)+std::fabs(yb-yc)) < halfCarTolerance )
>> 860 {
>> 861 return kInfinity;
>> 862 }
>> 863 }
>> 864 G4double a=(yb-ya)*zac-(yc-ya)*zab;
>> 865 G4double b=(xc-xa)*zab-(xb-xa)*zac;
>> 866 G4double c=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
>> 867 G4double d=-xa*a-ya*b+fDz*c;
>> 868 G4double t=a*v.x()+b*v.y()+c*v.z();
>> 869
>> 870 if (t!=0)
>> 871 {
>> 872 t=-(a*p.x()+b*p.y()+c*p.z()+d)/t;
>> 873 }
>> 874 if ( (t<halfCarTolerance) && (t>-halfCarTolerance) )
574 { 875 {
575 if (fTwist[i] == 0.) << 876 if (NormalToPlane(p,ipl).dot(v)<kCarTolerance)
576 { 877 {
577 G4double dd = fSurf[i].D*px + fSurf[i].E << 878 t=kInfinity;
578 dist = std::max(dd, dist); <<
579 } 879 }
580 else 880 else
581 { 881 {
582 // comptute distance to the wedge (two p << 882 t=0;
583 auto j = (i + 1)%4; <<
584 G4double cx = pp[j].x() - pp[i].x(); <<
585 G4double cy = pp[j].y() - pp[i].y(); <<
586 G4double d = (pp[i].x() - px)*cy + (py - <<
587 G4ThreeVector na(-cy, cx, fDelta[i].x()* <<
588 G4ThreeVector nb(-cy, cx, fDelta[j].x()* <<
589 G4double amag2 = na.mag2(); <<
590 G4double bmag2 = nb.mag2(); <<
591 G4double distab = (amag2 > bmag2) ? d/st <<
592 dist = std::max(distab, dist); <<
593 } 883 }
594 } 884 }
595 return (dist > 0.) ? dist : 0.; // return sa << 885 if (Inside(p+v*t) != kSurface) { t=kInfinity; }
596 } << 886
>> 887 return t;
>> 888 }
>> 889
>> 890 // --------------------------------------------------------------------
597 891
598 ////////////////////////////////////////////// <<
599 // <<
600 // Calculate distance to surface from inside <<
601 // <<
602 G4double G4GenericTrap::DistanceToOut(const G4 892 G4double G4GenericTrap::DistanceToOut(const G4ThreeVector& p,
603 const G4 893 const G4ThreeVector& v,
604 const G4 894 const G4bool calcNorm,
605 G4 895 G4bool* validNorm,
606 G4 896 G4ThreeVector* n) const
607 { 897 {
608 G4double px = p.x(), py = p.y(), pz = p.z(); << 898 #ifdef G4TESS_TEST
609 G4double vx = v.x(), vy = v.y(), vz = v.z(); << 899 if ( fTessellatedSolid )
>> 900 {
>> 901 return fTessellatedSolid->DistanceToOut(p, v, calcNorm, validNorm, n);
>> 902 }
>> 903 #endif
610 904
611 // Check intersections with plane faces << 905 static const G4double halfCarTolerance=kCarTolerance*0.5;
612 // << 906
613 if ((std::abs(pz) - fDz) >= -halfTolerance & << 907 G4double distmin;
>> 908 G4bool lateral_cross = false;
>> 909 ESide side = kUndefined;
>> 910
>> 911 if (calcNorm) { *validNorm=true; } // All normals are valid
>> 912
>> 913 if (v.z() < 0)
614 { 914 {
615 if (calcNorm) << 915 distmin=(-fDz-p.z())/v.z();
616 { << 916 if (calcNorm) { side=kMZ; *n=G4ThreeVector(0,0,-1); }
617 *validNorm = true; <<
618 n->set(0., 0., std::copysign(1., pz)); <<
619 } <<
620 return 0.; <<
621 } 917 }
622 G4double tout = (vz == 0) ? DBL_MAX : (std:: << 918 else
623 G4int iface = (vz < 0) ? -4 : -2; // little <<
624 <<
625 for (auto i = 0; i < 4; ++i) <<
626 { 919 {
627 if (fTwist[i] != 0) continue; // skip twis << 920 if (v.z() > 0)
628 const G4GenericTrapPlane& surf = fPlane[2* << 921 {
629 G4double cosa = surf.A*vx + surf.B*vy + su << 922 distmin = (fDz-p.z())/v.z();
630 if (cosa > 0) << 923 if (calcNorm) { side=kPZ; *n=G4ThreeVector(0,0,1); }
631 { << 924 }
632 G4double dist = surf.A*px + surf.B*py + << 925 else { distmin = kInfinity; }
633 if (dist >= -halfTolerance) << 926 }
>> 927
>> 928 G4double dz2 =0.5/fDz;
>> 929 G4double xa,xb,xc,xd;
>> 930 G4double ya,yb,yc,yd;
>> 931
>> 932 for (G4int ipl=0; ipl<4; ipl++)
>> 933 {
>> 934 G4int j = (ipl+1)%4;
>> 935 xa=fVertices[ipl].x();
>> 936 ya=fVertices[ipl].y();
>> 937 xb=fVertices[ipl+4].x();
>> 938 yb=fVertices[ipl+4].y();
>> 939 xc=fVertices[j].x();
>> 940 yc=fVertices[j].y();
>> 941 xd=fVertices[4+j].x();
>> 942 yd=fVertices[4+j].y();
>> 943
>> 944 if ( ((std::fabs(xb-xd)+std::fabs(yb-yd))<halfCarTolerance)
>> 945 || ((std::fabs(xa-xc)+std::fabs(ya-yc))<halfCarTolerance) )
>> 946 {
>> 947 G4double s=DistToTriangle(p,v,ipl) ;
>> 948 if ( (s>=0) && (s<distmin) )
>> 949 {
>> 950 distmin=s;
>> 951 lateral_cross=true;
>> 952 side=ESide(ipl+1);
>> 953 }
>> 954 continue;
>> 955 }
>> 956 G4double tx1 =dz2*(xb-xa);
>> 957 G4double ty1 =dz2*(yb-ya);
>> 958 G4double tx2 =dz2*(xd-xc);
>> 959 G4double ty2 =dz2*(yd-yc);
>> 960 G4double dzp =fDz+p.z();
>> 961 G4double xs1 =xa+tx1*dzp;
>> 962 G4double ys1 =ya+ty1*dzp;
>> 963 G4double xs2 =xc+tx2*dzp;
>> 964 G4double ys2 =yc+ty2*dzp;
>> 965 G4double dxs =xs2-xs1;
>> 966 G4double dys =ys2-ys1;
>> 967 G4double dtx =tx2-tx1;
>> 968 G4double dty =ty2-ty1;
>> 969 G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
>> 970 G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
>> 971 + tx1*ys2-tx2*ys1)*v.z();
>> 972 G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
>> 973 G4double s=kInfinity;
>> 974
>> 975 if (std::fabs(a) < kCarTolerance)
>> 976 {
>> 977 if (std::fabs(b) < kCarTolerance) { continue; }
>> 978 s=-c/b;
>> 979
>> 980 // Check for Point on the Surface
>> 981 //
>> 982 if ((s > -halfCarTolerance) && (s < distmin))
634 { 983 {
635 if (calcNorm) << 984 if (s < halfCarTolerance)
636 { 985 {
637 *validNorm = true; << 986 if (NormalToPlane(p,ipl).dot(v)<0.) { continue; }
638 n->set(surf.A, surf.B, surf.C); <<
639 } 987 }
640 return 0.; << 988 distmin =s;
641 } << 989 lateral_cross=true;
642 G4double tmp = -dist/cosa; << 990 side=ESide(ipl+1);
643 if (tout > tmp) { tout = tmp; iface = i; << 991 }
>> 992 continue;
644 } 993 }
645 } << 994 G4double d=b*b-4*a*c;
>> 995 if (d >= 0.)
>> 996 {
>> 997 if (a > 0) { s=0.5*(-b-std::sqrt(d))/a; }
>> 998 else { s=0.5*(-b+std::sqrt(d))/a; }
646 999
647 // Check intersections with twisted faces << 1000 // Check for Point on the Surface
648 // << 1001 //
649 constexpr G4double EPSILON = 1000.*DBL_EPSIL << 1002 if (s > -halfCarTolerance )
650 for (auto i = 0; i < 4; ++i) << 1003 {
651 { << 1004 if (s < distmin)
652 if (fTwist[i] == 0) continue; // skip plan <<
653 <<
654 // twisted face, solve quadratic equation <<
655 const G4GenericTrapSurface& surf = fSurf[i <<
656 G4double ABC = surf.A*vx + surf.B*vy + su <<
657 G4double ABCF = surf.A*px + surf.B*py + su <<
658 G4double A = ABC*vz; <<
659 G4double B = 0.5*(surf.D*vx + surf.E*vy + <<
660 G4double C = surf.D*px + surf.E*py + ABCF* <<
661 <<
662 if (std::abs(A) <= EPSILON) <<
663 { <<
664 // case 1: track is parallel to the surf <<
665 if (std::abs(B) <= EPSILON) { continue; <<
666 <<
667 // case 2: single root <<
668 G4double tmp = -0.5*C/B; <<
669 // compute normal at intersection point <<
670 G4double x = px + vx*tmp; <<
671 G4double y = py + vy*tmp; <<
672 G4double z = pz + vz*tmp; <<
673 G4double nx = surf.A*z + surf.D; <<
674 G4double ny = surf.B*z + surf.E; <<
675 G4double nz = surf.A*x + surf.B*y + 2.*s <<
676 <<
677 if (nx*vx + ny*vy + nz*vz > 0.) // point <<
678 { <<
679 auto k = (i + 1)%4; <<
680 G4double tz = (pz + fDz); <<
681 G4TwoVector a = fVertices[i] + fDelta[ <<
682 G4TwoVector b = fVertices[k] + fDelta[ <<
683 G4double dx = b.x() - a.x(); <<
684 G4double dy = b.y() - a.y(); <<
685 G4double leng = std::sqrt(dx*dx + dy*d <<
686 G4double dist = dx*(py - a.y()) - dy*( <<
687 if (dist >= -halfTolerance*leng) // po <<
688 { 1005 {
689 if (calcNorm) << 1006 if(s < halfCarTolerance)
690 { 1007 {
691 *validNorm = false; << 1008 if (NormalToPlane(p,ipl).dot(v)<0.) // Check second root
692 G4double normx = surf.A*pz + surf. << 1009 {
693 G4double normy = surf.B*pz + surf. << 1010 if (a > 0) { s=0.5*(-b+std::sqrt(d))/a; }
694 G4double normz = surf.A*px + surf. << 1011 else { s=0.5*(-b-std::sqrt(d))/a; }
695 G4double inv = 1./std::sqrt(normx* << 1012 if (( s > halfCarTolerance) && (s < distmin))
696 n->set(normx*inv, normy*inv, normz << 1013 {
>> 1014 distmin=s;
>> 1015 lateral_cross = true;
>> 1016 side=ESide(ipl+1);
>> 1017 }
>> 1018 continue;
>> 1019 }
697 } 1020 }
698 return 0.; << 1021 distmin = s;
>> 1022 lateral_cross = true;
>> 1023 side=ESide(ipl+1);
699 } 1024 }
700 if (tout > tmp) { tout = tmp; iface = <<
701 } 1025 }
702 continue; << 1026 else
703 } <<
704 <<
705 // case 3: scratch or no intersection <<
706 G4double D = B*B - A*C; <<
707 if (D < 0.25*fScratch*fScratch*A*A) <<
708 { <<
709 // check position of the point <<
710 auto j = (i + 1)%4; <<
711 G4double tz = pz + fDz; <<
712 G4TwoVector a = fVertices[i] + fDelta[i] <<
713 G4TwoVector b = fVertices[j] + fDelta[j] <<
714 G4double dx = b.x() - a.x(); <<
715 G4double dy = b.y() - a.y(); <<
716 G4double leng = std::sqrt(dx*dx + dy*dy) <<
717 G4double dist = dx*(py - a.y()) - dy*(px <<
718 if (dist <= -halfTolerance*leng) { cont <<
719 if (A > 0 || dist > halfTolerance*leng) <<
720 { 1027 {
721 if (calcNorm) << 1028 if (a > 0) { s=0.5*(-b+std::sqrt(d))/a; }
>> 1029 else { s=0.5*(-b-std::sqrt(d))/a; }
>> 1030
>> 1031 // Check for Point on the Surface
>> 1032 //
>> 1033 if ((s > -halfCarTolerance) && (s < distmin))
722 { 1034 {
723 *validNorm = false; << 1035 if (s < halfCarTolerance)
724 G4double nx = surf.A*pz + surf.D; << 1036 {
725 G4double ny = surf.B*pz + surf.E; << 1037 if (NormalToPlane(p,ipl).dot(v)<0.) // Check second root
726 G4double nz = surf.A*px + surf.B*py << 1038 {
727 G4double inv = 1./std::sqrt(nx*nx + << 1039 if (a > 0) { s=0.5*(-b-std::sqrt(d))/a; }
728 n->set(nx*inv, ny*inv, nz*inv); << 1040 else { s=0.5*(-b+std::sqrt(d))/a; }
729 } << 1041 if ( ( s > halfCarTolerance) && (s < distmin) )
730 return 0.; << 1042 {
>> 1043 distmin=s;
>> 1044 lateral_cross = true;
>> 1045 side=ESide(ipl+1);
>> 1046 }
>> 1047 continue;
>> 1048 }
>> 1049 }
>> 1050 distmin =s;
>> 1051 lateral_cross = true;
>> 1052 side=ESide(ipl+1);
>> 1053 }
731 } 1054 }
732 continue; <<
733 } 1055 }
>> 1056 }
>> 1057 if (!lateral_cross) // Make sure that track crosses the top or bottom
>> 1058 {
>> 1059 if (distmin >= kInfinity) { distmin=kCarTolerance; }
>> 1060 G4ThreeVector pt=p+distmin*v;
>> 1061
>> 1062 // Check if propagated point is in the polygon
>> 1063 //
>> 1064 G4int i=0;
>> 1065 if (v.z()>0.) { i=4; }
>> 1066 std::vector<G4TwoVector> xy;
>> 1067 for ( G4int j=0; j<4; j++) { xy.push_back(fVertices[i+j]); }
>> 1068
>> 1069 // Check Inside
>> 1070 //
>> 1071 if (InsidePolygone(pt,xy)==kOutside)
>> 1072 {
>> 1073 if(calcNorm)
>> 1074 {
>> 1075 if (v.z()>0) {side= kPZ; *n = G4ThreeVector(0,0,1);}
>> 1076 else { side=kMZ; *n = G4ThreeVector(0,0,-1);}
>> 1077 }
>> 1078 return 0.;
>> 1079 }
>> 1080 else
>> 1081 {
>> 1082 if(v.z()>0) {side=kPZ;}
>> 1083 else {side=kMZ;}
>> 1084 }
>> 1085 }
>> 1086
>> 1087 if (calcNorm)
>> 1088 {
>> 1089 G4ThreeVector pt=p+v*distmin;
>> 1090 switch (side)
>> 1091 {
>> 1092 case kXY0:
>> 1093 *n=NormalToPlane(pt,0);
>> 1094 break;
>> 1095 case kXY1:
>> 1096 *n=NormalToPlane(pt,1);
>> 1097 break;
>> 1098 case kXY2:
>> 1099 *n=NormalToPlane(pt,2);
>> 1100 break;
>> 1101 case kXY3:
>> 1102 *n=NormalToPlane(pt,3);
>> 1103 break;
>> 1104 case kPZ:
>> 1105 *n=G4ThreeVector(0,0,1);
>> 1106 break;
>> 1107 case kMZ:
>> 1108 *n=G4ThreeVector(0,0,-1);
>> 1109 break;
>> 1110 default:
>> 1111 DumpInfo();
>> 1112 std::ostringstream message;
>> 1113 G4int oldprc = message.precision(16);
>> 1114 message << "Undefined side for valid surface normal to solid." << G4endl
>> 1115 << "Position:" << G4endl
>> 1116 << " p.x() = " << p.x()/mm << " mm" << G4endl
>> 1117 << " p.y() = " << p.y()/mm << " mm" << G4endl
>> 1118 << " p.z() = " << p.z()/mm << " mm" << G4endl
>> 1119 << "Direction:" << G4endl
>> 1120 << " v.x() = " << v.x() << G4endl
>> 1121 << " v.y() = " << v.y() << G4endl
>> 1122 << " v.z() = " << v.z() << G4endl
>> 1123 << "Proposed distance :" << G4endl
>> 1124 << " distmin = " << distmin/mm << " mm";
>> 1125 message.precision(oldprc);
>> 1126 G4Exception("G4GenericTrap::DistanceToOut(p,v,..)",
>> 1127 "GeomSolids1002", JustWarning, message);
>> 1128 break;
>> 1129 }
>> 1130 }
>> 1131
>> 1132 if (distmin<halfCarTolerance) { distmin=0.; }
>> 1133
>> 1134 return distmin;
>> 1135 }
>> 1136
>> 1137 // --------------------------------------------------------------------
>> 1138
>> 1139 G4double G4GenericTrap::DistanceToOut(const G4ThreeVector& p) const
>> 1140 {
>> 1141
>> 1142 #ifdef G4TESS_TEST
>> 1143 if ( fTessellatedSolid )
>> 1144 {
>> 1145 return fTessellatedSolid->DistanceToOut(p);
>> 1146 }
>> 1147 #endif
734 1148
735 // case 4: two intersection points << 1149 G4double safz = fDz-std::fabs(p.z());
736 G4double tmp = -B - std::copysign(std::sqr << 1150 if (safz<0) { safz = 0; }
737 G4double t1 = tmp / A; << 1151
738 G4double t2 = C / tmp; << 1152 G4double safe = safz;
739 G4double tmin = std::min(t1, t2); << 1153 G4double safxy = safz;
740 G4double tmax = std::max(t1, t2); << 1154
741 << 1155 for (G4int iseg=0; iseg<4; iseg++)
742 if (A < 0) // concave profile: tmin(out) - << 1156 {
743 { << 1157 safxy = std::fabs(SafetyToFace(p,iseg));
744 if (std::abs(tmax) < std::abs(tmin)) con << 1158 if (safxy < safe) { safe = safxy; }
745 if (tmin <= halfTolerance) // point is o << 1159 }
746 { << 1160
747 G4double t = 0.5*(tmin + tmax); << 1161 return safe;
748 G4double x = px + vx*t; << 1162 }
749 G4double y = py + vy*t; << 1163
750 G4double z = pz + vz*t; << 1164 // --------------------------------------------------------------------
751 << 1165
752 auto j = (i + 1)%4; << 1166 G4bool G4GenericTrap::CalculateExtent(const EAxis pAxis,
753 G4double tz = z + fDz; << 1167 const G4VoxelLimits& pVoxelLimit,
754 G4TwoVector a = fVertices[i] + fDelta[ << 1168 const G4AffineTransform& pTransform,
755 G4TwoVector b = fVertices[j] + fDelta[ << 1169 G4double& pMin, G4double& pMax) const
756 G4double dx = b.x() - a.x(); << 1170 {
757 G4double dy = b.y() - a.y(); << 1171 #ifdef G4TESS_TEST
758 G4double leng = std::sqrt(dx*dx + dy*d << 1172 if ( fTessellatedSolid )
759 G4double dist = dx*(y - a.y()) - dy*(x << 1173 {
760 if (dist <= -halfTolerance*leng) cont << 1174 return fTessellatedSolid->CalculateExtent(pAxis, pVoxelLimit,
761 if (calcNorm) << 1175 pTransform, pMin, pMax);
>> 1176 }
>> 1177 #endif
>> 1178
>> 1179 // Computes bounding vectors for a shape
>> 1180 //
>> 1181 G4double Dx,Dy;
>> 1182 G4ThreeVector minVec = GetMinimumBBox();
>> 1183 G4ThreeVector maxVec = GetMaximumBBox();
>> 1184 Dx = 0.5*(maxVec.x()- minVec.x());
>> 1185 Dy = 0.5*(maxVec.y()- minVec.y());
>> 1186
>> 1187 if (!pTransform.IsRotated())
>> 1188 {
>> 1189 // Special case handling for unrotated shapes
>> 1190 // Compute x/y/z mins and maxs respecting limits, with early returns
>> 1191 // if outside limits. Then switch() on pAxis
>> 1192 //
>> 1193 G4double xoffset,xMin,xMax;
>> 1194 G4double yoffset,yMin,yMax;
>> 1195 G4double zoffset,zMin,zMax;
>> 1196
>> 1197 xoffset=pTransform.NetTranslation().x();
>> 1198 xMin=xoffset-Dx;
>> 1199 xMax=xoffset+Dx;
>> 1200 if (pVoxelLimit.IsXLimited())
>> 1201 {
>> 1202 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
>> 1203 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
>> 1204 {
>> 1205 return false;
>> 1206 }
>> 1207 else
>> 1208 {
>> 1209 if (xMin<pVoxelLimit.GetMinXExtent())
762 { 1210 {
763 *validNorm = false; << 1211 xMin=pVoxelLimit.GetMinXExtent();
764 G4double nx = surf.A*pz + surf.D; << 1212 }
765 G4double ny = surf.B*pz + surf.E; << 1213 if (xMax>pVoxelLimit.GetMaxXExtent())
766 G4double nz = surf.A*px + surf.B*py << 1214 {
767 G4double inv = 1./std::sqrt(nx*nx + << 1215 xMax=pVoxelLimit.GetMaxXExtent();
768 n->set(nx*inv, ny*inv, nz*inv); <<
769 } 1216 }
770 return 0.; <<
771 } 1217 }
772 if (tout > tmin) { tout = tmin; iface = <<
773 } 1218 }
774 else // convex profile: tmin(in) -> tmax(o << 1219
>> 1220 yoffset=pTransform.NetTranslation().y();
>> 1221 yMin=yoffset-Dy;
>> 1222 yMax=yoffset+Dy;
>> 1223 if (pVoxelLimit.IsYLimited())
775 { 1224 {
776 if (tmax < halfTolerance) // point is on << 1225 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
>> 1226 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
777 { 1227 {
778 if (calcNorm) << 1228 return false;
>> 1229 }
>> 1230 else
>> 1231 {
>> 1232 if (yMin<pVoxelLimit.GetMinYExtent())
>> 1233 {
>> 1234 yMin=pVoxelLimit.GetMinYExtent();
>> 1235 }
>> 1236 if (yMax>pVoxelLimit.GetMaxYExtent())
779 { 1237 {
780 *validNorm = false; << 1238 yMax=pVoxelLimit.GetMaxYExtent();
781 G4double nx = surf.A*pz + surf.D; <<
782 G4double ny = surf.B*pz + surf.E; <<
783 G4double nz = surf.A*px + surf.B*py <<
784 G4double inv = 1./std::sqrt(nx*nx + <<
785 n->set(nx*inv, ny*inv, nz*inv); <<
786 } 1239 }
787 return 0.; <<
788 } 1240 }
789 if (tout > tmax) { tout = tmax; iface = <<
790 } 1241 }
791 } <<
792 1242
793 // Compute normal, if required, and return d << 1243 zoffset=pTransform.NetTranslation().z();
794 // << 1244 zMin=zoffset-fDz;
795 if (tout < halfTolerance) tout = 0.; << 1245 zMax=zoffset+fDz;
796 if (calcNorm) << 1246 if (pVoxelLimit.IsZLimited())
797 { <<
798 if (iface < 0) <<
799 { <<
800 *validNorm = true; <<
801 n->set(0, 0, iface + 3); // little trick <<
802 } <<
803 else <<
804 { 1247 {
805 const G4GenericTrapSurface& surf = fSurf << 1248 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
806 if (fTwist[iface] == 0) << 1249 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
807 { 1250 {
808 *validNorm = true; << 1251 return false;
809 n->set(surf.D, surf.E, surf.F); <<
810 } 1252 }
811 else 1253 else
812 { 1254 {
813 *validNorm = false; << 1255 if (zMin<pVoxelLimit.GetMinZExtent())
814 G4double x = px + vx*tout; << 1256 {
815 G4double y = py + vy*tout; << 1257 zMin=pVoxelLimit.GetMinZExtent();
816 G4double z = pz + vz*tout; << 1258 }
817 G4double nx = surf.A*z + surf.D; << 1259 if (zMax>pVoxelLimit.GetMaxZExtent())
818 G4double ny = surf.B*z + surf.E; << 1260 {
819 G4double nz = surf.A*x + surf.B*y + 2. << 1261 zMax=pVoxelLimit.GetMaxZExtent();
820 G4double inv = 1./std::sqrt(nx*nx + ny << 1262 }
821 n->set(nx*inv, ny*inv, nz*inv); <<
822 } 1263 }
823 } 1264 }
>> 1265
>> 1266 switch (pAxis)
>> 1267 {
>> 1268 case kXAxis:
>> 1269 pMin = xMin;
>> 1270 pMax = xMax;
>> 1271 break;
>> 1272 case kYAxis:
>> 1273 pMin = yMin;
>> 1274 pMax = yMax;
>> 1275 break;
>> 1276 case kZAxis:
>> 1277 pMin = zMin;
>> 1278 pMax = zMax;
>> 1279 break;
>> 1280 default:
>> 1281 break;
>> 1282 }
>> 1283 pMin-=kCarTolerance;
>> 1284 pMax+=kCarTolerance;
>> 1285
>> 1286 return true;
824 } 1287 }
825 return tout; << 1288 else
826 } << 1289 {
>> 1290 // General rotated case - create and clip mesh to boundaries
827 1291
828 ////////////////////////////////////////////// << 1292 G4bool existsAfterClip=false;
829 // << 1293 G4ThreeVectorList *vertices;
830 // Estimate safety distance to the surface fro <<
831 // <<
832 G4double G4GenericTrap::DistanceToOut(const G4 <<
833 { <<
834 G4double px = p.x(), py = p.y(), pz = p.z(); <<
835 1294
836 // intersect edges by z = pz plane << 1295 pMin=+kInfinity;
837 G4TwoVector pp[4]; << 1296 pMax=-kInfinity;
838 G4double z = (pz + fDz); <<
839 for (auto i = 0; i < 4; ++i) { pp[i] = fVert <<
840 1297
841 // estimate distance to the solid << 1298 // Calculate rotated vertex coordinates
842 G4double dist = std::abs(pz) - fDz; << 1299 //
843 for (auto i = 0; i < 4; ++i) << 1300 vertices=CreateRotatedVertices(pTransform);
844 { << 1301 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
845 if (fTwist[i] == 0.) << 1302 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
846 { << 1303 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
847 G4double dd = fSurf[i].D*px + fSurf[i].E << 1304
848 dist = std::max(dd, dist); << 1305 if ( (pMin!=kInfinity) || (pMax!=-kInfinity) )
>> 1306 {
>> 1307 existsAfterClip=true;
>> 1308
>> 1309 // Add 2*tolerance to avoid precision troubles
>> 1310 //
>> 1311 pMin-=kCarTolerance;
>> 1312 pMax+=kCarTolerance;
849 } 1313 }
850 else 1314 else
851 { 1315 {
852 // comptute distance to the wedge (two p << 1316 // Check for case where completely enveloping clipping volume.
853 auto j = (i + 1)%4; << 1317 // If point inside then we are confident that the solid completely
854 G4double cx = pp[j].x() - pp[i].x(); << 1318 // envelopes the clipping volume. Hence set min/max extents according
855 G4double cy = pp[j].y() - pp[i].y(); << 1319 // to clipping volume extents along the specified axis.
856 G4double d = (pp[i].x() - px)*cy + (py - << 1320 //
857 G4ThreeVector na(-cy, cx, fDelta[i].x()* << 1321 G4ThreeVector clipCentre(
858 G4ThreeVector nb(-cy, cx, fDelta[j].x()* << 1322 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
859 G4double amag2 = na.mag2(); << 1323 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
860 G4double bmag2 = nb.mag2(); << 1324 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
861 G4double distab = (amag2 > bmag2) ? d/st << 1325
862 dist = std::max(distab, dist); << 1326 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
>> 1327 {
>> 1328 existsAfterClip=true;
>> 1329 pMin=pVoxelLimit.GetMinExtent(pAxis);
>> 1330 pMax=pVoxelLimit.GetMaxExtent(pAxis);
>> 1331 }
863 } 1332 }
>> 1333 delete vertices;
>> 1334 return existsAfterClip;
864 } 1335 }
865 return (dist < 0.) ? -dist : 0.; // return s <<
866 } 1336 }
867 1337
868 ////////////////////////////////////////////// << 1338 // --------------------------------------------------------------------
869 // <<
870 // Return bounding box <<
871 // <<
872 void G4GenericTrap::BoundingLimits(G4ThreeVect <<
873 G4ThreeVect <<
874 { <<
875 pMin = fMinBBox; <<
876 pMax = fMaxBBox; <<
877 } <<
878 <<
879 ////////////////////////////////////////////// <<
880 // <<
881 // Calculate extent under transform and specif <<
882 // <<
883 G4bool <<
884 G4GenericTrap::CalculateExtent(const EAxis pAx <<
885 const G4VoxelLi <<
886 const G4AffineT <<
887 G4double& <<
888 { <<
889 G4ThreeVector bmin, bmax; <<
890 G4bool exist; <<
891 <<
892 // Check bounding box (bbox) <<
893 // <<
894 BoundingLimits(bmin,bmax); <<
895 G4BoundingEnvelope bbox(bmin,bmax); <<
896 #ifdef G4BBOX_EXTENT <<
897 return bbox.CalculateExtent(pAxis,pVoxelLimi <<
898 #endif <<
899 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox <<
900 { <<
901 return exist = pMin < pMax; <<
902 } <<
903 1339
904 // Set bounding envelope (benv) and calculat << 1340 G4ThreeVectorList*
905 // << 1341 G4GenericTrap::CreateRotatedVertices(const G4AffineTransform& pTransform) const
906 // To build the bounding envelope with plane << 1342 {
907 // the trapezoid is subdivided in two triang << 1343 // Create a List containing the transformed vertices
908 // duplication of vertices in the bases in a << 1344 // Ordering [0-3] -fDz cross section
909 // a convex polyhedron (some faces of the en << 1345 // [4-7] +fDz cross section such that [0] is below [4],
910 // << 1346 // [1] below [5] etc.
911 G4double dz = GetZHalfLength(); << 1347 // Note: caller has deletion responsibility
912 G4ThreeVectorList baseA(8), baseB(8); << 1348
913 for (G4int i = 0; i < 4; ++i) << 1349 G4ThreeVector Min = GetMinimumBBox();
914 { << 1350 G4ThreeVector Max = GetMaximumBBox();
915 G4TwoVector va = GetVertex(i); << 1351
916 G4TwoVector vb = GetVertex(i+4); << 1352 G4ThreeVectorList *vertices;
917 baseA[2*i].set(va.x(), va.y(),-dz); << 1353 vertices=new G4ThreeVectorList();
918 baseB[2*i].set(vb.x(), vb.y(), dz); << 1354
>> 1355 if (vertices)
>> 1356 {
>> 1357 vertices->reserve(8);
>> 1358 G4ThreeVector vertex0(Min.x(),Min.y(),Min.z());
>> 1359 G4ThreeVector vertex1(Max.x(),Min.y(),Min.z());
>> 1360 G4ThreeVector vertex2(Max.x(),Max.y(),Min.z());
>> 1361 G4ThreeVector vertex3(Min.x(),Max.y(),Min.z());
>> 1362 G4ThreeVector vertex4(Min.x(),Min.y(),Max.z());
>> 1363 G4ThreeVector vertex5(Max.x(),Min.y(),Max.z());
>> 1364 G4ThreeVector vertex6(Max.x(),Max.y(),Max.z());
>> 1365 G4ThreeVector vertex7(Min.x(),Max.y(),Max.z());
>> 1366
>> 1367 vertices->push_back(pTransform.TransformPoint(vertex0));
>> 1368 vertices->push_back(pTransform.TransformPoint(vertex1));
>> 1369 vertices->push_back(pTransform.TransformPoint(vertex2));
>> 1370 vertices->push_back(pTransform.TransformPoint(vertex3));
>> 1371 vertices->push_back(pTransform.TransformPoint(vertex4));
>> 1372 vertices->push_back(pTransform.TransformPoint(vertex5));
>> 1373 vertices->push_back(pTransform.TransformPoint(vertex6));
>> 1374 vertices->push_back(pTransform.TransformPoint(vertex7));
919 } 1375 }
920 for (G4int i = 0; i < 4; ++i) << 1376 else
921 { 1377 {
922 G4int k1 = 2*i, k2 = (2*i + 2)%8; << 1378 G4Exception("G4GenericTrap::CreateRotatedVertices()", "FatalError",
923 G4double ax = (baseA[k2].x() - baseA[k1].x << 1379 FatalException, "Out of memory - Cannot allocate vertices!");
924 G4double ay = (baseA[k2].y() - baseA[k1].y <<
925 G4double bx = (baseB[k2].x() - baseB[k1].x <<
926 G4double by = (baseB[k2].y() - baseB[k1].y <<
927 G4double znorm = ax*by - ay*bx; <<
928 baseA[k1+1] = (znorm < 0.0) ? baseA[k2] : <<
929 baseB[k1+1] = (znorm < 0.0) ? baseB[k1] : <<
930 } 1380 }
931 << 1381 return vertices;
932 std::vector<const G4ThreeVectorList *> polyg <<
933 polygons[0] = &baseA; <<
934 polygons[1] = &baseB; <<
935 <<
936 G4BoundingEnvelope benv(bmin, bmax, polygons <<
937 exist = benv.CalculateExtent(pAxis,pVoxelLim <<
938 return exist; <<
939 } 1382 }
>> 1383
>> 1384 // --------------------------------------------------------------------
940 1385
941 ////////////////////////////////////////////// <<
942 // <<
943 // Return type of the solid <<
944 // <<
945 G4GeometryType G4GenericTrap::GetEntityType() 1386 G4GeometryType G4GenericTrap::GetEntityType() const
946 { 1387 {
947 return { "G4GenericTrap" }; << 1388 return G4String("G4GenericTrap");
948 } <<
949 <<
950 ////////////////////////////////////////////// <<
951 // <<
952 // Return true if not twisted, false otherwise <<
953 // <<
954 G4bool G4GenericTrap::IsFaceted() const <<
955 { <<
956 return (!fIsTwisted); <<
957 } 1389 }
>> 1390
>> 1391 // --------------------------------------------------------------------
958 1392
959 ////////////////////////////////////////////// <<
960 // <<
961 // Make a clone of the solid <<
962 // <<
963 G4VSolid* G4GenericTrap::Clone() const 1393 G4VSolid* G4GenericTrap::Clone() const
964 { 1394 {
965 return new G4GenericTrap(*this); 1395 return new G4GenericTrap(*this);
966 } 1396 }
967 1397
968 ////////////////////////////////////////////// << 1398 // --------------------------------------------------------------------
969 // << 1399
970 // Write parameters of the solid to the specif <<
971 // <<
972 std::ostream& G4GenericTrap::StreamInfo(std::o 1400 std::ostream& G4GenericTrap::StreamInfo(std::ostream& os) const
973 { 1401 {
974 G4long oldprc = os.precision(16); << 1402 G4int oldprc = os.precision(16);
975 os << "------------------------------------- 1403 os << "-----------------------------------------------------------\n"
976 << " *** Dump for solid - " << GetName << 1404 << " *** Dump for solid - " << GetName() << " *** \n"
977 << " ================================= << 1405 << " =================================================== \n"
978 << "Solid geometry type: " << GetEntityTy << 1406 << " Solid geometry type: " << GetEntityType() << G4endl
979 << " half length Z: " << fDz/mm << "\n" << 1407 << " half length Z: " << fDz/mm << " mm \n"
980 << " list of vertices:\n"; 1408 << " list of vertices:\n";
981 for (G4int i = 0; i < 8; ++i) << 1409
>> 1410 for ( G4int i=0; i<fgkNofVertices; ++i )
982 { 1411 {
983 os << " #" << i << " " << fVertices[i] << 1412 os << std::setw(5) << "#" << i
>> 1413 << " vx = " << fVertices[i].x()/mm << " mm"
>> 1414 << " vy = " << fVertices[i].y()/mm << " mm" << G4endl;
984 } 1415 }
985 os << "------------------------------------- <<
986 os.precision(oldprc); 1416 os.precision(oldprc);
>> 1417
987 return os; 1418 return os;
988 } << 1419 }
>> 1420
>> 1421 // --------------------------------------------------------------------
989 1422
990 ////////////////////////////////////////////// <<
991 // <<
992 // Pick up a random point on the surface <<
993 // <<
994 G4ThreeVector G4GenericTrap::GetPointOnSurface 1423 G4ThreeVector G4GenericTrap::GetPointOnSurface() const
995 { 1424 {
996 if (fArea[0] + fArea[1] + fArea[2] + fArea[3 << 1425
997 { << 1426 #ifdef G4TESS_TEST
998 G4AutoLock l(&lateralareaMutex); << 1427 if ( fTessellatedSolid )
999 fArea[0] = GetLateralFaceArea(0); << 1428 {
1000 fArea[1] = GetLateralFaceArea(1); << 1429 return fTessellatedSolid->GetPointOnSurface();
1001 fArea[2] = GetLateralFaceArea(2); << 1430 }
1002 fArea[3] = GetLateralFaceArea(3); << 1431 #endif
1003 l.unlock(); <<
1004 } <<
1005 <<
1006 constexpr G4int iface[6][4] = <<
1007 { {0,1,2,3}, {0,4,5,1}, {1,5,6,2}, {2,6,7 <<
1008 <<
1009 G4bool isTwisted[6] = {false}; <<
1010 for (auto i = 0; i < 4; ++i) { isTwisted[i <<
1011 <<
1012 G4double ssurf[6]; <<
1013 G4TwoVector A = fVertices[3] - fVertices[1] <<
1014 G4TwoVector B = fVertices[2] - fVertices[0] <<
1015 G4TwoVector C = fVertices[7] - fVertices[5] <<
1016 G4TwoVector D = fVertices[6] - fVertices[4] <<
1017 ssurf[0] = (A.x()*B.y() - A.y()*B.x())*0.5; <<
1018 ssurf[1] = ssurf[0] + fArea[0]; <<
1019 ssurf[2] = ssurf[1] + fArea[1]; <<
1020 ssurf[3] = ssurf[2] + fArea[2]; <<
1021 ssurf[4] = ssurf[3] + fArea[3]; <<
1022 ssurf[5] = ssurf[4] + (C.x()*D.y() - C.y()* <<
1023 <<
1024 G4double select = ssurf[5]*G4QuickRand(); <<
1025 G4int k; <<
1026 for (k = 0; k < 5; ++k) { if (select <= ssu <<
1027 <<
1028 G4int i0 = iface[k][0]; <<
1029 G4int i1 = iface[k][1]; <<
1030 G4int i2 = iface[k][2]; <<
1031 G4int i3 = iface[k][3]; <<
1032 G4ThreeVector pp[4]; <<
1033 pp[0].set(fVertices[i0].x(), fVertices[i0]. <<
1034 pp[1].set(fVertices[i1].x(), fVertices[i1]. <<
1035 pp[2].set(fVertices[i2].x(), fVertices[i2]. <<
1036 pp[3].set(fVertices[i3].x(), fVertices[i3]. <<
1037 1432
1038 G4ThreeVector point; 1433 G4ThreeVector point;
1039 if (isTwisted[k]) << 1434 G4TwoVector u,v,w;
1040 { // twisted face, rejection sampling << 1435 G4double rand,area,chose,cf,lambda0,lambda1,alfa,beta,zp;
1041 G4double maxlength = std::max((pp[2] - pp << 1436 G4int ipl,j;
1042 point = (pp[0] + pp[1] + pp[2] + pp[3])*0 << 1437
1043 for (auto i = 0; i < 10000; ++i) << 1438 std::vector<G4ThreeVector> vertices;
1044 { << 1439 for (G4int i=0; i<4;i++)
1045 G4double u = G4QuickRand(); << 1440 {
1046 G4ThreeVector p0 = pp[0] + (pp[1] - pp[ << 1441 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
1047 G4ThreeVector p1 = pp[3] + (pp[2] - pp[ <<
1048 G4double v = G4QuickRand()*(maxlength/( <<
1049 if (v >= 1.) continue; <<
1050 point = p0 + (p1 - p0)*v; <<
1051 break; <<
1052 } <<
1053 } 1442 }
1054 else << 1443 for (G4int i=4; i<8;i++)
1055 { // plane face << 1444 {
1056 G4double u = G4QuickRand(); << 1445 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
1057 G4double v = G4QuickRand(); << 1446 }
1058 if (u + v > 1.) { u = 1. - u; v = 1. - v; << 1447
1059 G4double ss = (((pp[2] - pp[0]).cross(pp[ << 1448 // Surface Area of Planes(only estimation for twisted)
1060 G4int j = (select > ssurf[k] - ss) ? 3 : << 1449 //
1061 point = pp[0] + (pp[2] - pp[0])*u + (pp[j << 1450 G4double Surface0=GetFaceSurfaceArea(vertices[0],vertices[1],
>> 1451 vertices[2],vertices[3]);//-fDz plane
>> 1452 G4double Surface1=GetFaceSurfaceArea(vertices[0],vertices[1],
>> 1453 vertices[5],vertices[4]);// Lat plane
>> 1454 G4double Surface2=GetFaceSurfaceArea(vertices[3],vertices[0],
>> 1455 vertices[4],vertices[7]);// Lat plane
>> 1456 G4double Surface3=GetFaceSurfaceArea(vertices[2],vertices[3],
>> 1457 vertices[7],vertices[6]);// Lat plane
>> 1458 G4double Surface4=GetFaceSurfaceArea(vertices[2],vertices[1],
>> 1459 vertices[5],vertices[6]);// Lat plane
>> 1460 G4double Surface5=GetFaceSurfaceArea(vertices[4],vertices[5],
>> 1461 vertices[6],vertices[7]);// fDz plane
>> 1462 rand = G4UniformRand();
>> 1463 area = Surface0+Surface1+Surface2+Surface3+Surface4+Surface5;
>> 1464 chose = rand*area;
>> 1465
>> 1466 if ( ( chose < Surface0)
>> 1467 || ( chose > (Surface0+Surface1+Surface2+Surface3+Surface4)) )
>> 1468 { // fDz or -fDz Plane
>> 1469 ipl = G4int(G4UniformRand()*4);
>> 1470 j = (ipl+1)%4;
>> 1471 if(chose < Surface0)
>> 1472 {
>> 1473 zp = -fDz;
>> 1474 u = fVertices[ipl]; v = fVertices[j];
>> 1475 w = fVertices[(ipl+3)%4];
>> 1476 }
>> 1477 else
>> 1478 {
>> 1479 zp = fDz;
>> 1480 u = fVertices[ipl+4]; v = fVertices[j+4];
>> 1481 w = fVertices[(ipl+3)%4+4];
>> 1482 }
>> 1483 alfa = G4UniformRand();
>> 1484 beta = G4UniformRand();
>> 1485 lambda1=alfa*beta;
>> 1486 lambda0=alfa-lambda1;
>> 1487 v = v-u;
>> 1488 w = w-u;
>> 1489 v = u+lambda0*v+lambda1*w;
>> 1490 }
>> 1491 else // Lateral Plane Twisted or Not
>> 1492 {
>> 1493 if (chose < Surface0+Surface1) { ipl=0; }
>> 1494 else if (chose < Surface0+Surface1+Surface2) { ipl=1; }
>> 1495 else if (chose < Surface0+Surface1+Surface2+Surface3) { ipl=2; }
>> 1496 else { ipl=3; }
>> 1497 j = (ipl+1)%4;
>> 1498 zp = -fDz+G4UniformRand()*2*fDz;
>> 1499 cf = 0.5*(fDz-zp)/fDz;
>> 1500 u = fVertices[ipl+4]+cf*( fVertices[ipl]-fVertices[ipl+4]);
>> 1501 v = fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]);
>> 1502 v = u+(v-u)*G4UniformRand();
1062 } 1503 }
>> 1504 point=G4ThreeVector(v.x(),v.y(),zp);
>> 1505
1063 return point; 1506 return point;
1064 } 1507 }
1065 1508
1066 ///////////////////////////////////////////// << 1509 // --------------------------------------------------------------------
1067 // << 1510
1068 // Return volume of the solid <<
1069 // <<
1070 G4double G4GenericTrap::GetCubicVolume() 1511 G4double G4GenericTrap::GetCubicVolume()
1071 { 1512 {
1072 if (fCubicVolume == 0.0) << 1513 if(fCubicVolume != 0.) {;}
1073 { << 1514 else { fCubicVolume = G4VSolid::GetCubicVolume(); }
1074 // diagonals <<
1075 G4TwoVector A = fVertices[3] - fVertices[ <<
1076 G4TwoVector B = fVertices[2] - fVertices[ <<
1077 G4TwoVector C = fVertices[7] - fVertices[ <<
1078 G4TwoVector D = fVertices[6] - fVertices[ <<
1079 <<
1080 // kross products <<
1081 G4double AB = A.x()*B.y() - A.y()*B.x(); <<
1082 G4double CD = C.x()*D.y() - C.y()*D.x(); <<
1083 G4double AD = A.x()*D.y() - A.y()*D.x(); <<
1084 G4double CB = C.x()*B.y() - C.y()*B.x(); <<
1085 <<
1086 fCubicVolume = ((AB + CD)/3. + (AD + CB)/ <<
1087 } <<
1088 return fCubicVolume; 1515 return fCubicVolume;
1089 } 1516 }
1090 1517
1091 ///////////////////////////////////////////// << 1518 // --------------------------------------------------------------------
1092 // <<
1093 // Compute lateral face area <<
1094 // <<
1095 G4double G4GenericTrap::GetLateralFaceArea(G4 <<
1096 { <<
1097 G4int i1 = iface, i2 = (i1 + 1)%4, i3 = i1 <<
1098 1519
1099 // plane face << 1520 G4double G4GenericTrap::GetSurfaceArea()
1100 if (fTwist[iface] == 0.0) << 1521 {
>> 1522 if(fSurfaceArea != 0.) {;}
>> 1523 else
1101 { 1524 {
1102 G4ThreeVector p1(fVertices[i1].x(), fVert << 1525 std::vector<G4ThreeVector> vertices;
1103 G4ThreeVector p2(fVertices[i2].x(), fVert << 1526 for (G4int i=0; i<4;i++)
1104 G4ThreeVector p3(fVertices[i3].x(), fVert << 1527 {
1105 G4ThreeVector p4(fVertices[i4].x(), fVert << 1528 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
1106 return ((p4 - p1).cross(p3 - p2)).mag()*0 << 1529 }
>> 1530 for (G4int i=4; i<8;i++)
>> 1531 {
>> 1532 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
>> 1533 }
>> 1534
>> 1535 // Surface Area of Planes(only estimation for twisted)
>> 1536 //
>> 1537 G4double fSurface0=GetFaceSurfaceArea(vertices[0],vertices[1],
>> 1538 vertices[2],vertices[3]);//-fDz plane
>> 1539 G4double fSurface1=GetFaceSurfaceArea(vertices[0],vertices[1],
>> 1540 vertices[5],vertices[4]);// Lat plane
>> 1541 G4double fSurface2=GetFaceSurfaceArea(vertices[3],vertices[0],
>> 1542 vertices[4],vertices[7]);// Lat plane
>> 1543 G4double fSurface3=GetFaceSurfaceArea(vertices[2],vertices[3],
>> 1544 vertices[7],vertices[6]);// Lat plane
>> 1545 G4double fSurface4=GetFaceSurfaceArea(vertices[2],vertices[1],
>> 1546 vertices[5],vertices[6]);// Lat plane
>> 1547 G4double fSurface5=GetFaceSurfaceArea(vertices[4],vertices[5],
>> 1548 vertices[6],vertices[7]);// fDz plane
>> 1549
>> 1550 // Total Surface Area
>> 1551 //
>> 1552 if(!fIsTwisted)
>> 1553 {
>> 1554 fSurfaceArea = fSurface0+fSurface1+fSurface2
>> 1555 + fSurface3+fSurface4+fSurface5;
>> 1556 }
>> 1557 else
>> 1558 {
>> 1559 fSurfaceArea = G4VSolid::GetSurfaceArea();
>> 1560 }
1107 } 1561 }
>> 1562 return fSurfaceArea;
>> 1563 }
1108 1564
1109 // twisted face << 1565 // --------------------------------------------------------------------
1110 constexpr G4int NSTEP = 250; << 1566
1111 constexpr G4double dt = 1./NSTEP; << 1567 G4double G4GenericTrap::GetFaceSurfaceArea(const G4ThreeVector& p0,
>> 1568 const G4ThreeVector& p1,
>> 1569 const G4ThreeVector& p2,
>> 1570 const G4ThreeVector& p3) const
>> 1571 {
>> 1572 // Auxiliary method for Get Surface Area of Face
>> 1573
>> 1574 G4double aOne, aTwo;
>> 1575 G4ThreeVector t, u, v, w, Area, normal;
>> 1576
>> 1577 t = p2 - p1;
>> 1578 u = p0 - p1;
>> 1579 v = p2 - p3;
>> 1580 w = p0 - p3;
>> 1581
>> 1582 Area = w.cross(v);
>> 1583 aOne = 0.5*Area.mag();
>> 1584
>> 1585 Area = t.cross(u);
>> 1586 aTwo = 0.5*Area.mag();
>> 1587
>> 1588 return aOne + aTwo;
>> 1589 }
>> 1590
>> 1591 // --------------------------------------------------------------------
1112 1592
1113 G4double x21 = fVertices[i2].x() - fVertice << 1593 G4bool G4GenericTrap::ComputeIsTwisted()
1114 G4double y21 = fVertices[i2].y() - fVertice << 1594 {
1115 G4double x31 = fVertices[i3].x() - fVertice << 1595 // Computes tangents of twist angles (angles between projections on XY plane
1116 G4double y31 = fVertices[i3].y() - fVertice << 1596 // of corresponding -dz +dz edges).
1117 G4double x42 = fVertices[i4].x() - fVertice <<
1118 G4double y42 = fVertices[i4].y() - fVertice <<
1119 G4double x43 = fVertices[i4].x() - fVertice <<
1120 G4double y43 = fVertices[i4].y() - fVertice <<
1121 1597
1122 G4double A = x21*y43 - y21*x43; << 1598 G4bool twisted = false;
1123 G4double B0 = x21*y31 - y21*x31; << 1599 G4double dx1, dy1, dx2, dy2;
1124 G4double B1 = x42*y31 - y42*x31; << 1600 G4int nv = fgkNofVertices/2;
1125 G4double HH = 4*fDz*fDz; <<
1126 G4double invAA = 1./(A*A); <<
1127 G4double sqrtAA = 2.*std::abs(A); <<
1128 G4double invSqrtAA = 1./sqrtAA; <<
1129 1601
1130 G4double area = 0.; << 1602 for ( G4int i=0; i<4; i++ )
1131 for (G4int i = 0; i < NSTEP; ++i) <<
1132 { 1603 {
1133 G4double t = (i + 0.5)*dt; << 1604 dx1 = fVertices[(i+1)%nv].x()-fVertices[i].x();
1134 G4double I = y21 + (y43 - y21)*t; << 1605 dy1 = fVertices[(i+1)%nv].y()-fVertices[i].y();
1135 G4double J = x21 + (x43 - x21)*t; << 1606 if ( (dx1 == 0) && (dy1 == 0) ) { continue; }
1136 G4double IIJJ = HH*(I*I + J*J); << 1607
1137 G4double B = B1*t + B0; << 1608 dx2 = fVertices[nv+(i+1)%nv].x()-fVertices[nv+i].x();
>> 1609 dy2 = fVertices[nv+(i+1)%nv].y()-fVertices[nv+i].y();
1138 1610
1139 G4double aa = A*A; << 1611 if ( dx2 == 0 && dy2 == 0 ) { continue; }
1140 G4double bb = 2.*A*B; << 1612 G4double twist_angle = std::fabs(dy1*dx2 - dx1*dy2);
1141 G4double cc = IIJJ + B*B; << 1613 if ( twist_angle < fgkTolerance ) { continue; }
>> 1614 twisted = true;
>> 1615 SetTwistAngle(i,twist_angle);
1142 1616
1143 G4double R1 = std::sqrt(aa + bb + cc); << 1617 // Check on big angles, potentially navigation problem
1144 G4double R0 = std::sqrt(cc); <<
1145 G4double log1 = std::log(std::abs(sqrtAA* <<
1146 G4double log0 = std::log(std::abs(sqrtAA* <<
1147 1618
1148 area += 0.5*R1 + 0.25*bb*invAA*(R1 - R0) << 1619 twist_angle = std::acos( (dx1*dx2 + dy1*dy2)
>> 1620 / (std::sqrt(dx1*dx1+dy1*dy1)
>> 1621 * std::sqrt(dx2*dx2+dy2*dy2)) );
>> 1622
>> 1623 if ( std::fabs(twist_angle) > 0.5*pi+kCarTolerance )
>> 1624 {
>> 1625 std::ostringstream message;
>> 1626 message << "Twisted Angle is bigger than 90 degrees - " << GetName()
>> 1627 << G4endl
>> 1628 << " Potential problem of malformed Solid !" << G4endl
>> 1629 << " TwistANGLE = " << twist_angle
>> 1630 << "*rad for lateral plane N= " << i;
>> 1631 G4Exception("G4GenericTrap::ComputeIsTwisted()", "GeomSolids1002",
>> 1632 JustWarning, message);
>> 1633 }
1149 } 1634 }
1150 return area*dt; << 1635
>> 1636 return twisted;
1151 } 1637 }
1152 1638
1153 ///////////////////////////////////////////// << 1639 // --------------------------------------------------------------------
1154 // << 1640
1155 // Return surface area of the solid << 1641 G4bool G4GenericTrap::CheckOrder(const std::vector<G4TwoVector>& vertices) const
1156 // <<
1157 G4double G4GenericTrap::GetSurfaceArea() <<
1158 { 1642 {
1159 if (fSurfaceArea == 0.0) << 1643 // Test if the vertices are in a clockwise order, if not reorder them.
>> 1644 // Also test if they're well defined without crossing opposite segments
>> 1645
>> 1646 G4bool clockwise_order=true;
>> 1647 G4double sum1 = 0.;
>> 1648 G4double sum2 = 0.;
>> 1649 G4int j;
>> 1650
>> 1651 for (G4int i=0; i<4; i++)
1160 { 1652 {
1161 G4TwoVector A = fVertices[3] - fVertices[ << 1653 j = (i+1)%4;
1162 G4TwoVector B = fVertices[2] - fVertices[ << 1654 sum1 += vertices[i].x()*vertices[j].y() - vertices[j].x()*vertices[i].y();
1163 G4TwoVector C = fVertices[7] - fVertices[ << 1655 sum2 += vertices[i+4].x()*vertices[j+4].y()
1164 G4TwoVector D = fVertices[6] - fVertices[ << 1656 - vertices[j+4].x()*vertices[i+4].y();
1165 G4double S_bot = (A.x()*B.y() - A.y()*B.x <<
1166 G4double S_top = (C.x()*D.y() - C.y()*D.x <<
1167 fArea[0] = GetLateralFaceArea(0); <<
1168 fArea[1] = GetLateralFaceArea(1); <<
1169 fArea[2] = GetLateralFaceArea(2); <<
1170 fArea[3] = GetLateralFaceArea(3); <<
1171 fSurfaceArea = S_bot + S_top + fArea[0] + <<
1172 } 1657 }
1173 return fSurfaceArea; << 1658 if (sum1*sum2 < -fgkTolerance)
>> 1659 {
>> 1660 std::ostringstream message;
>> 1661 message << "Lower/upper faces defined with opposite clockwise - "
>> 1662 << GetName();
>> 1663 G4Exception("G4GenericTrap::CheckOrder()", "GeomSolids0002",
>> 1664 FatalException, message);
>> 1665 }
>> 1666
>> 1667 if ((sum1 > 0.)||(sum2 > 0.))
>> 1668 {
>> 1669 std::ostringstream message;
>> 1670 message << "Vertices must be defined in clockwise XY planes - "
>> 1671 << GetName();
>> 1672 G4Exception("G4GenericTrap::CheckOrder()", "GeomSolids1001",
>> 1673 JustWarning,message, "Re-ordering...");
>> 1674 clockwise_order = false;
>> 1675 }
>> 1676
>> 1677 // Check for illegal crossings
>> 1678 //
>> 1679 G4bool illegal_cross = false;
>> 1680 illegal_cross = IsSegCrossingZ(vertices[0],vertices[4],
>> 1681 vertices[1],vertices[5]);
>> 1682
>> 1683 if (!illegal_cross)
>> 1684 {
>> 1685 illegal_cross = IsSegCrossingZ(vertices[2],vertices[6],
>> 1686 vertices[3],vertices[7]);
>> 1687 }
>> 1688 // +/- dZ planes
>> 1689 if (!illegal_cross)
>> 1690 {
>> 1691 illegal_cross = IsSegCrossing(vertices[0],vertices[1],
>> 1692 vertices[2],vertices[3]);
>> 1693 }
>> 1694 if (!illegal_cross)
>> 1695 {
>> 1696 illegal_cross = IsSegCrossing(vertices[0],vertices[3],
>> 1697 vertices[1],vertices[2]);
>> 1698 }
>> 1699 if (!illegal_cross)
>> 1700 {
>> 1701 illegal_cross = IsSegCrossing(vertices[4],vertices[5],
>> 1702 vertices[6],vertices[7]);
>> 1703 }
>> 1704 if (!illegal_cross)
>> 1705 {
>> 1706 illegal_cross = IsSegCrossing(vertices[4],vertices[7],
>> 1707 vertices[5],vertices[6]);
>> 1708 }
>> 1709
>> 1710 if (illegal_cross)
>> 1711 {
>> 1712 std::ostringstream message;
>> 1713 message << "Malformed polygone with opposite sides - " << GetName();
>> 1714 G4Exception("G4GenericTrap::CheckOrderAndSetup()",
>> 1715 "GeomSolids0002", FatalException, message);
>> 1716 }
>> 1717 return clockwise_order;
1174 } 1718 }
1175 1719
1176 ///////////////////////////////////////////// << 1720 // --------------------------------------------------------------------
1177 // << 1721
1178 // Check parameters of the solid << 1722 void G4GenericTrap::ReorderVertices(std::vector<G4ThreeVector>& vertices) const
1179 // <<
1180 void <<
1181 G4GenericTrap::CheckParameters(G4double halfZ <<
1182 const std::vec <<
1183 { 1723 {
1184 // Check number of vertices << 1724 // Reorder the vector of vertices
1185 // << 1725
1186 if (vertices.size() != 8) << 1726 std::vector<G4ThreeVector> oldVertices(vertices);
>> 1727
>> 1728 for ( G4int i=0; i < G4int(oldVertices.size()); ++i )
1187 { 1729 {
1188 std::ostringstream message; << 1730 vertices[i] = oldVertices[oldVertices.size()-1-i];
1189 message << "Number of vertices is " << ve << 1731 }
1190 << " instead of 8 for Solid: " << << 1732 }
1191 G4Exception("G4GenericTrap::CheckParamete << 1733
1192 FatalException, message); << 1734 // --------------------------------------------------------------------
1193 } <<
1194 1735
1195 // Check dZ << 1736 G4bool
1196 // << 1737 G4GenericTrap::IsSegCrossing(const G4TwoVector& a, const G4TwoVector& b,
1197 if ((fDz = halfZ) < 2.*kCarTolerance) << 1738 const G4TwoVector& c, const G4TwoVector& d) const
>> 1739 {
>> 1740 // Check if segments [A,B] and [C,D] are crossing
>> 1741
>> 1742 G4bool stand1 = false;
>> 1743 G4bool stand2 = false;
>> 1744 G4double dx1,dx2,xm=0.,ym=0.,a1=0.,a2=0.,b1=0.,b2=0.;
>> 1745 dx1=(b-a).x();
>> 1746 dx2=(d-c).x();
>> 1747
>> 1748 if( std::fabs(dx1) < fgkTolerance ) { stand1 = true; }
>> 1749 if( std::fabs(dx2) < fgkTolerance ) { stand2 = true; }
>> 1750 if (!stand1)
>> 1751 {
>> 1752 a1 = (b.x()*a.y()-a.x()*b.y())/dx1;
>> 1753 b1 = (b-a).y()/dx1;
>> 1754 }
>> 1755 if (!stand2)
>> 1756 {
>> 1757 a2 = (d.x()*c.y()-c.x()*d.y())/dx2;
>> 1758 b2 = (d-c).y()/dx2;
>> 1759 }
>> 1760 if (stand1 && stand2)
>> 1761 {
>> 1762 // Segments parallel and vertical
>> 1763 //
>> 1764 if (std::fabs(a.x()-c.x())<fgkTolerance)
>> 1765 {
>> 1766 // Check if segments are overlapping
>> 1767 //
>> 1768 if ( ((c.y()-a.y())*(c.y()-b.y())<-fgkTolerance)
>> 1769 || ((d.y()-a.y())*(d.y()-b.y())<-fgkTolerance)
>> 1770 || ((a.y()-c.y())*(a.y()-d.y())<-fgkTolerance)
>> 1771 || ((b.y()-c.y())*(b.y()-d.y())<-fgkTolerance) ) { return true; }
>> 1772
>> 1773 return false;
>> 1774 }
>> 1775 // Different x values
>> 1776 //
>> 1777 return false;
>> 1778 }
>> 1779
>> 1780 if (stand1) // First segment vertical
1198 { 1781 {
1199 std::ostringstream message; << 1782 xm = a.x();
1200 message << "Z-dimension is too small or n << 1783 ym = a2+b2*xm;
1201 << " for Solid: " << GetName() << <<
1202 G4Exception("G4GenericTrap::CheckParamete <<
1203 FatalException, message); <<
1204 } 1784 }
>> 1785 else
>> 1786 {
>> 1787 if (stand2) // Second segment vertical
>> 1788 {
>> 1789 xm = c.x();
>> 1790 ym = a1+b1*xm;
>> 1791 }
>> 1792 else // Normal crossing
>> 1793 {
>> 1794 if (std::fabs(b1-b2) < fgkTolerance)
>> 1795 {
>> 1796 // Parallel segments, are they aligned
>> 1797 //
>> 1798 if (std::fabs(c.y()-(a1+b1*c.x())) > fgkTolerance) { return false; }
>> 1799
>> 1800 // Aligned segments, are they overlapping
>> 1801 //
>> 1802 if ( ((c.x()-a.x())*(c.x()-b.x())<-fgkTolerance)
>> 1803 || ((d.x()-a.x())*(d.x()-b.x())<-fgkTolerance)
>> 1804 || ((a.x()-c.x())*(a.x()-d.x())<-fgkTolerance)
>> 1805 || ((b.x()-c.x())*(b.x()-d.x())<-fgkTolerance) ) { return true; }
1205 1806
1206 // Check order of vertices << 1807 return false;
>> 1808 }
>> 1809 xm = (a1-a2)/(b2-b1);
>> 1810 ym = (a1*b2-a2*b1)/(b2-b1);
>> 1811 }
>> 1812 }
>> 1813
>> 1814 // Check if crossing point is both between A,B and C,D
>> 1815 //
>> 1816 G4double check = (xm-a.x())*(xm-b.x())+(ym-a.y())*(ym-b.y());
>> 1817 if (check > -fgkTolerance) { return false; }
>> 1818 check = (xm-c.x())*(xm-d.x())+(ym-c.y())*(ym-d.y());
>> 1819 if (check > -fgkTolerance) { return false; }
>> 1820
>> 1821 return true;
>> 1822 }
>> 1823
>> 1824 // --------------------------------------------------------------------
>> 1825
>> 1826 G4bool
>> 1827 G4GenericTrap::IsSegCrossingZ(const G4TwoVector& a, const G4TwoVector& b,
>> 1828 const G4TwoVector& c, const G4TwoVector& d) const
>> 1829 {
>> 1830 // Check if segments [A,B] and [C,D] are crossing when
>> 1831 // A and C are on -dZ and B and D are on +dZ
>> 1832
>> 1833 // Calculate the Intersection point between two lines in 3D
1207 // 1834 //
1208 for (auto i = 0; i < 8; ++i) { fVertices.at << 1835 G4ThreeVector temp1,temp2;
>> 1836 G4ThreeVector v1,v2,p1,p2,p3,p4,dv;
>> 1837 G4double s,det;
>> 1838 p1=G4ThreeVector(a.x(),a.y(),-fDz);
>> 1839 p2=G4ThreeVector(c.x(),c.y(),-fDz);
>> 1840 p3=G4ThreeVector(b.x(),b.y(),fDz);
>> 1841 p4=G4ThreeVector(d.x(),d.y(),fDz);
>> 1842 v1=p3-p1;
>> 1843 v2=p4-p2;
>> 1844 dv=p2-p1;
1209 1845
1210 // Bottom polygon area << 1846 // In case of Collapsed Vertices No crossing
1211 G4TwoVector diag1 = fVertices[2] - fVertice << 1847 //
1212 G4TwoVector diag2 = fVertices[3] - fVertice << 1848 if( (std::fabs(dv.x()) < kCarTolerance )&&
1213 G4double ldiagbot = std::max(diag1.mag(), d << 1849 (std::fabs(dv.y()) < kCarTolerance ) ) { return false; }
1214 G4double zbot = diag1.x()*diag2.y() - diag1 << 1850
1215 if (std::abs(zbot) < ldiagbot*kCarTolerance << 1851 if( (std::fabs((p4-p3).x()) < kCarTolerance )&&
1216 << 1852 (std::fabs((p4-p3).y()) < kCarTolerance ) ) { return false; }
1217 // Top polygon area << 1853
1218 G4TwoVector diag3 = fVertices[6] - fVertice << 1854 // First estimate if Intersection is possible( if det is 0)
1219 G4TwoVector diag4 = fVertices[7] - fVertice << 1855 //
1220 G4double ldiagtop = std::max(diag3.mag(), d << 1856 det = dv.x()*v1.y()*v2.z()+dv.y()*v1.z()*v2.x()
1221 G4double ztop = diag3.x()*diag4.y() - diag3 << 1857 - dv.x()*v1.z()*v2.y()-dv.y()*v1.x()*v2.z();
1222 if (std::abs(ztop) < ldiagtop*kCarTolerance << 1858
>> 1859 if (std::fabs(det)<kCarTolerance) //Intersection
>> 1860 {
>> 1861 temp1 = v1.cross(v2);
>> 1862 temp2 = (p2-p1).cross(v2);
>> 1863 if (temp1.dot(temp2) < 0) { return false; } // intersection negative
>> 1864 s = temp1.mag();
>> 1865
>> 1866 if ( s < kCarTolerance ) { return false; } // parallel lines
>> 1867 s = ((dv).cross(v2)).mag()/s;
>> 1868
>> 1869 if(s < 1.-kCarTolerance) { return true; }
>> 1870 }
>> 1871 return false;
>> 1872 }
>> 1873
>> 1874 // --------------------------------------------------------------------
>> 1875
>> 1876 G4VFacet*
>> 1877 G4GenericTrap::MakeDownFacet(const std::vector<G4ThreeVector>& fromVertices,
>> 1878 G4int ind1, G4int ind2, G4int ind3) const
>> 1879 {
>> 1880 // Create a triangular facet from the polygon points given by indices
>> 1881 // forming the down side ( the normal goes in -z)
>> 1882 // Do not create facet if 2 vertices are the same
>> 1883
>> 1884 if ( (fromVertices[ind1] == fromVertices[ind2]) ||
>> 1885 (fromVertices[ind2] == fromVertices[ind3]) ||
>> 1886 (fromVertices[ind1] == fromVertices[ind3]) ) { return 0; }
>> 1887
>> 1888 std::vector<G4ThreeVector> vertices;
>> 1889 vertices.push_back(fromVertices[ind1]);
>> 1890 vertices.push_back(fromVertices[ind2]);
>> 1891 vertices.push_back(fromVertices[ind3]);
>> 1892
>> 1893 // first vertex most left
>> 1894 //
>> 1895 G4ThreeVector cross=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
1223 1896
1224 if (zbot*ztop < 0.) << 1897 if ( cross.z() > 0.0 )
1225 { 1898 {
>> 1899 // Should not happen, as vertices should have been reordered at this stage
>> 1900
1226 std::ostringstream message; 1901 std::ostringstream message;
1227 message << "Vertices of the bottom and to << 1902 message << "Vertices in wrong order - " << GetName();
1228 StreamInfo(message); << 1903 G4Exception("G4GenericTrap::MakeDownFacet", "GeomSolids0002",
1229 G4Exception("G4GenericTrap::CheckParamete <<
1230 FatalException, message); 1904 FatalException, message);
1231 } 1905 }
1232 if ((zbot > 0.) || (ztop > 0.)) << 1906
1233 { << 1907 return new G4TriangularFacet(vertices[0], vertices[1], vertices[2], ABSOLUTE);
1234 std::swap(fVertices[1], fVertices[3]); << 1908 }
1235 std::swap(fVertices[5], fVertices[7]); <<
1236 std::ostringstream message; <<
1237 message << "Vertices re-ordered in Solid: <<
1238 << "Vertices of the bottom and to <<
1239 G4Exception("G4GenericTrap::CheckParamete <<
1240 JustWarning, message); <<
1241 } <<
1242 1909
1243 // Check for degeneracy << 1910 // --------------------------------------------------------------------
>> 1911
>> 1912 G4VFacet*
>> 1913 G4GenericTrap::MakeUpFacet(const std::vector<G4ThreeVector>& fromVertices,
>> 1914 G4int ind1, G4int ind2, G4int ind3) const
>> 1915 {
>> 1916 // Create a triangular facet from the polygon points given by indices
>> 1917 // forming the upper side ( z>0 )
>> 1918
>> 1919 // Do not create facet if 2 vertices are the same
1244 // 1920 //
1245 G4int ndegenerated = 0; << 1921 if ( (fromVertices[ind1] == fromVertices[ind2]) ||
1246 for (auto i = 0; i < 4; ++i) << 1922 (fromVertices[ind2] == fromVertices[ind3]) ||
1247 { << 1923 (fromVertices[ind1] == fromVertices[ind3]) ) { return 0; }
1248 auto j = (i + 1)%4; << 1924
1249 G4double lbot = (fVertices[j] - fVertices << 1925 std::vector<G4ThreeVector> vertices;
1250 G4double ltop = (fVertices[j + 4] - fVert << 1926 vertices.push_back(fromVertices[ind1]);
1251 ndegenerated += (std::max(lbot, ltop) < k << 1927 vertices.push_back(fromVertices[ind2]);
1252 } << 1928 vertices.push_back(fromVertices[ind3]);
1253 if (ndegenerated > 1 || << 1929
1254 GetCubicVolume() < fDz*std::max(ldiagbo << 1930 // First vertex most left
>> 1931 //
>> 1932 G4ThreeVector cross=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
>> 1933
>> 1934 if ( cross.z() < 0.0 )
1255 { 1935 {
>> 1936 // Should not happen, as vertices should have been reordered at this stage
>> 1937
1256 std::ostringstream message; 1938 std::ostringstream message;
1257 message << "Degenerated solid !\n"; << 1939 message << "Vertices in wrong order - " << GetName();
1258 StreamInfo(message); << 1940 G4Exception("G4GenericTrap::MakeUpFacet", "GeomSolids0002",
1259 G4Exception("G4GenericTrap::CheckParamete <<
1260 FatalException, message); 1941 FatalException, message);
1261 } 1942 }
>> 1943
>> 1944 return new G4TriangularFacet(vertices[0], vertices[1], vertices[2], ABSOLUTE);
>> 1945 }
1262 1946
1263 // Check that the polygons are convex << 1947 // --------------------------------------------------------------------
1264 // << 1948
1265 G4bool isConvex = true; << 1949 G4VFacet*
1266 for (auto i = 0; i < 4; ++i) << 1950 G4GenericTrap::MakeSideFacet(const G4ThreeVector& downVertex0,
>> 1951 const G4ThreeVector& downVertex1,
>> 1952 const G4ThreeVector& upVertex1,
>> 1953 const G4ThreeVector& upVertex0) const
>> 1954 {
>> 1955 // Creates a triangular facet from the polygon points given by indices
>> 1956 // forming the upper side ( z>0 )
>> 1957
>> 1958 if ( (downVertex0 == downVertex1) && (upVertex0 == upVertex1) )
1267 { 1959 {
1268 auto j = (i + 1)%4; << 1960 return 0;
1269 auto k = (j + 1)%4; <<
1270 G4TwoVector edge1 = fVertices[j] - fVerti <<
1271 G4TwoVector edge2 = fVertices[k] - fVerti <<
1272 isConvex = ((edge1.x()*edge2.y() - edge1. <<
1273 if (!isConvex) break; <<
1274 G4TwoVector edge3 = fVertices[j + 4] - fV <<
1275 G4TwoVector edge4 = fVertices[k + 4] - fV <<
1276 isConvex = ((edge3.x()*edge4.y() - edge3. <<
1277 if (!isConvex) break; <<
1278 } 1961 }
1279 if (!isConvex) << 1962
>> 1963 if ( downVertex0 == downVertex1 )
1280 { 1964 {
1281 std::ostringstream message; << 1965 return new G4TriangularFacet(downVertex0, upVertex1, upVertex0, ABSOLUTE);
1282 message << "The bottom and top faces must <<
1283 StreamInfo(message); <<
1284 G4Exception("G4GenericTrap::CheckParamete <<
1285 FatalException, message); <<
1286 } 1966 }
1287 } <<
1288 1967
1289 ///////////////////////////////////////////// << 1968 if ( upVertex0 == upVertex1 )
1290 // <<
1291 // Compute surface equations and twist angles <<
1292 // <<
1293 void G4GenericTrap::ComputeLateralSurfaces() <<
1294 { <<
1295 for (auto i = 0; i < 4; ++i) <<
1296 { 1969 {
1297 auto j = (i + 1)%4; << 1970 return new G4TriangularFacet(downVertex0, downVertex1, upVertex0, ABSOLUTE);
1298 G4ThreeVector p1(fVertices[j].x(), fVerti <<
1299 G4ThreeVector p2(fVertices[i].x(), fVerti <<
1300 G4ThreeVector p3(fVertices[j + 4].x(), fV <<
1301 G4ThreeVector p4(fVertices[i + 4].x(), fV <<
1302 G4ThreeVector ebot = p2 - p1; <<
1303 G4ThreeVector etop = p4 - p3; <<
1304 G4double lbot = ebot.mag(); <<
1305 G4double ltop = etop.mag(); <<
1306 G4double zcross = ebot.x()*etop.y() - ebo <<
1307 G4double eps = kCarTolerance*std::max(lbo <<
1308 if (std::min(lbot, ltop) < kCarTolerance <<
1309 { // plane surface: Dx + Ey + Fz + G = 0 <<
1310 G4ThreeVector normal; <<
1311 if (std::max(lbot, ltop) < kCarToleranc <<
1312 { <<
1313 auto k = (j + 1)%4; <<
1314 auto l = (k + 1)%4; <<
1315 G4TwoVector vl = fVertices[l] + fVert <<
1316 G4TwoVector vi = fVertices[i] + fVert <<
1317 G4TwoVector vj = fVertices[j] + fVert <<
1318 G4TwoVector vk = fVertices[k] + fVert <<
1319 G4TwoVector vij = (vi - vl).unit() + <<
1320 G4ThreeVector epar = (p4 + p3 - p2 - <<
1321 G4ThreeVector eort = epar.cross(G4Thr <<
1322 normal = (eort.cross(epar)).unit(); <<
1323 } <<
1324 else <<
1325 { <<
1326 normal = ((p4 - p1).cross(p3 - p2)).u <<
1327 } <<
1328 fSurf[i].D = fPlane[2*i].A = fPlane[2*i <<
1329 fSurf[i].E = fPlane[2*i].B = fPlane[2*i <<
1330 fSurf[i].F = fPlane[2*i].C = fPlane[2*i <<
1331 fSurf[i].G = fPlane[2*i].D = fPlane[2*i <<
1332 } <<
1333 else <<
1334 { // hyperbolic paraboloid: Axz + Byz + C <<
1335 fIsTwisted = true; <<
1336 G4double angle = std::acos(ebot.dot(eto <<
1337 if (angle > CLHEP::halfpi) <<
1338 { <<
1339 std::ostringstream message; <<
1340 message << "Twist on " << angle/CLHEP <<
1341 << " degrees, should not be m <<
1342 StreamInfo(message); <<
1343 G4Exception("G4GenericTrap::ComputeLa <<
1344 FatalException, message); <<
1345 } <<
1346 fTwist[i] = std::copysign(angle, zcross <<
1347 // set equation of twisted surface (hyp <<
1348 fSurf[i].A = 2.*fDz*(p4.y() - p3.y() - <<
1349 fSurf[i].B =-2.*fDz*(p4.x() - p3.x() - <<
1350 fSurf[i].C = ((p4.x() - p2.x())*(p3.y() <<
1351 fSurf[i].D = 2.*fDz*fDz*(p4.y() - p3.y( <<
1352 fSurf[i].E =-2.*fDz*fDz*(p4.x() - p3.x( <<
1353 fSurf[i].F = 2.*fDz*(p4.x()*p3.y() - p3 <<
1354 fSurf[i].G = fDz*fDz*((p4.x() + p2.x()) <<
1355 G4double magnitude = G4ThreeVector(fSu <<
1356 if (magnitude < kCarTolerance) continue <<
1357 fSurf[i].A /= magnitude; <<
1358 fSurf[i].B /= magnitude; <<
1359 fSurf[i].C /= magnitude; <<
1360 fSurf[i].D /= magnitude; <<
1361 fSurf[i].E /= magnitude; <<
1362 fSurf[i].F /= magnitude; <<
1363 fSurf[i].G /= magnitude; <<
1364 // set planes of bounding polyhedron <<
1365 G4ThreeVector normal1, normal2; <<
1366 G4ThreeVector c1, c2; <<
1367 if (fTwist[i] < 0.) <<
1368 { <<
1369 normal1 = ((p2 - p1).cross(p4 - p1)). <<
1370 normal2 = ((p3 - p4).cross(p1 - p4)). <<
1371 c1 = p1; <<
1372 c2 = p4; <<
1373 } <<
1374 else <<
1375 { <<
1376 normal1 = ((p3 - p2).cross(p1 - p2)). <<
1377 normal2 = ((p2 - p3).cross(p4 - p3)). <<
1378 c1 = p2; <<
1379 c2 = p3; <<
1380 } <<
1381 fPlane[2*i].A = normal1.x(); <<
1382 fPlane[2*i].B = normal1.y(); <<
1383 fPlane[2*i].C = normal1.z(); <<
1384 fPlane[2*i].D = -normal1.dot(c1); <<
1385 fPlane[2*i + 1].A = normal2.x(); <<
1386 fPlane[2*i + 1].B = normal2.y(); <<
1387 fPlane[2*i + 1].C = normal2.z(); <<
1388 fPlane[2*i + 1].D = -normal2.dot(c2); <<
1389 } <<
1390 fDelta[i] = (fVertices[i + 4] - fVertices <<
1391 } 1971 }
1392 } <<
1393 1972
1394 ///////////////////////////////////////////// << 1973 return new G4QuadrangularFacet(downVertex0, downVertex1,
1395 // << 1974 upVertex1, upVertex0, ABSOLUTE);
1396 // Set bounding box << 1975 }
1397 // << 1976
1398 void G4GenericTrap::ComputeBoundingBox() << 1977 // --------------------------------------------------------------------
>> 1978
>> 1979 G4TessellatedSolid* G4GenericTrap::CreateTessellatedSolid() const
1399 { 1980 {
1400 G4double minX, maxX, minY, maxY; << 1981 // 3D vertices
1401 minX = maxX = fVertices[0].x(); << 1982 //
1402 minY = maxY = fVertices[0].y(); << 1983 G4int nv = fgkNofVertices/2;
1403 for (auto i = 1; i < 8; ++i) << 1984 std::vector<G4ThreeVector> downVertices;
1404 { << 1985 for ( G4int i=0; i<nv; i++ )
1405 minX = std::min(minX, fVertices[i].x()); << 1986 {
1406 maxX = std::max(maxX, fVertices[i].x()); << 1987 downVertices.push_back(G4ThreeVector(fVertices[i].x(),
1407 minY = std::min(minY, fVertices[i].y()); << 1988 fVertices[i].y(), -fDz));
1408 maxY = std::max(maxY, fVertices[i].y()); << 1989 }
>> 1990
>> 1991 std::vector<G4ThreeVector> upVertices;
>> 1992 for ( G4int i=nv; i<2*nv; i++ )
>> 1993 {
>> 1994 upVertices.push_back(G4ThreeVector(fVertices[i].x(),
>> 1995 fVertices[i].y(), fDz));
1409 } 1996 }
1410 fMinBBox = G4ThreeVector(minX, minY,-fDz); << 1997
1411 fMaxBBox = G4ThreeVector(maxX, maxY, fDz); << 1998 // Reorder vertices if they are not ordered anti-clock wise
>> 1999 //
>> 2000 G4ThreeVector cross
>> 2001 = (downVertices[1]-downVertices[0]).cross(downVertices[2]-downVertices[1]);
>> 2002 G4ThreeVector cross1
>> 2003 = (upVertices[1]-upVertices[0]).cross(upVertices[2]-upVertices[1]);
>> 2004 if ( (cross.z() > 0.0) || (cross1.z() > 0.0) )
>> 2005 {
>> 2006 ReorderVertices(downVertices);
>> 2007 ReorderVertices(upVertices);
>> 2008 }
>> 2009
>> 2010 G4TessellatedSolid* tessellatedSolid = new G4TessellatedSolid(GetName());
>> 2011
>> 2012 G4VFacet* facet = 0;
>> 2013 facet = MakeDownFacet(downVertices, 0, 1, 2);
>> 2014 if (facet) { tessellatedSolid->AddFacet( facet ); }
>> 2015 facet = MakeDownFacet(downVertices, 0, 2, 3);
>> 2016 if (facet) { tessellatedSolid->AddFacet( facet ); }
>> 2017 facet = MakeUpFacet(upVertices, 0, 2, 1);
>> 2018 if (facet) { tessellatedSolid->AddFacet( facet ); }
>> 2019 facet = MakeUpFacet(upVertices, 0, 3, 2);
>> 2020 if (facet) { tessellatedSolid->AddFacet( facet ); }
1412 2021
1413 // Check correctness of the bounding box << 2022 // The quadrangular sides
1414 // 2023 //
1415 if (minX >= maxX || minY >= maxY || -fDz >= << 2024 for ( G4int i = 0; i < nv; ++i )
1416 { 2025 {
1417 std::ostringstream message; << 2026 G4int j = (i+1) % nv;
1418 message << "Bad bounding box (min >= max) << 2027 facet = MakeSideFacet(downVertices[j], downVertices[i],
1419 << GetName() << " !" << 2028 upVertices[i], upVertices[j]);
1420 << "\npMin = " << fMinBBox << 2029
1421 << "\npMax = " << fMaxBBox; << 2030 if ( facet ) { tessellatedSolid->AddFacet( facet ); }
1422 G4Exception("G4GenericTrap::ComputeBoundi <<
1423 JustWarning, message); <<
1424 DumpInfo(); <<
1425 } 2031 }
1426 } <<
1427 2032
1428 ///////////////////////////////////////////// << 2033 tessellatedSolid->SetSolidClosed(true);
1429 // << 2034
1430 // Set max length of a scratch << 2035 return tessellatedSolid;
1431 // << 2036 }
1432 void G4GenericTrap::ComputeScratchLength() << 2037
>> 2038 // --------------------------------------------------------------------
>> 2039
>> 2040 void G4GenericTrap::ComputeBBox()
1433 { 2041 {
1434 G4double scratch = kInfinity; << 2042 // Computes bounding box for a shape.
1435 for (auto i = 0; i < 4; ++i) << 2043
>> 2044 G4double minX, maxX, minY, maxY;
>> 2045 minX = maxX = fVertices[0].x();
>> 2046 minY = maxY = fVertices[0].y();
>> 2047
>> 2048 for (G4int i=1; i< fgkNofVertices; i++)
1436 { 2049 {
1437 if (fTwist[i] == 0.) continue; // skip pl << 2050 if (minX>fVertices[i].x()) { minX=fVertices[i].x(); }
1438 auto k = (i + 1)%4; << 2051 if (maxX<fVertices[i].x()) { maxX=fVertices[i].x(); }
1439 G4ThreeVector p1(fVertices[i].x(), fVerti << 2052 if (minY>fVertices[i].y()) { minY=fVertices[i].y(); }
1440 G4ThreeVector p2(fVertices[k].x(), fVerti << 2053 if (maxY<fVertices[i].y()) { maxY=fVertices[i].y(); }
1441 G4ThreeVector p3(fVertices[i + 4].x(), fV <<
1442 G4ThreeVector p4(fVertices[k + 4].x(), fV <<
1443 G4ThreeVector p0 = (p1 + p2 + p3 + p4)*0. <<
1444 G4ThreeVector norm = SurfaceNormal(p0); <<
1445 G4ThreeVector pp[2]; // points inside and <<
1446 pp[0] = p0 - norm * halfTolerance; <<
1447 pp[1] = p0 + norm * halfTolerance; <<
1448 G4ThreeVector vv[2]; // unit vectors alon <<
1449 vv[0] = (p4 - p1).unit(); <<
1450 vv[1] = (p3 - p2).unit(); <<
1451 // find intersection points and compute t <<
1452 for (auto ip = 0; ip < 2; ++ip) <<
1453 { <<
1454 G4double px = pp[ip].x(); <<
1455 G4double py = pp[ip].y(); <<
1456 G4double pz = pp[ip].z(); <<
1457 for (auto iv = 0; iv < 2; ++iv) <<
1458 { <<
1459 G4double vx = vv[iv].x(); <<
1460 G4double vy = vv[iv].y(); <<
1461 G4double vz = vv[iv].z(); <<
1462 // solve quadratic equation <<
1463 G4double ABC = fSurf[i].A*vx + fSurf <<
1464 G4double ABCF = fSurf[i].A*px + fSurf <<
1465 G4double A = ABC*vz; <<
1466 G4double B = 0.5*(fSurf[i].D*vx + fSu <<
1467 G4double C = fSurf[i].D*px + fSurf[i] <<
1468 G4double D = B*B - A*C; <<
1469 if (D < 0) continue; <<
1470 G4double leng = 2.*std::sqrt(D)/std:: <<
1471 scratch = std::min(leng, scratch); <<
1472 } <<
1473 } <<
1474 } 2054 }
1475 fScratch = std::max(kCarTolerance, scratch) << 2055 fMinBBoxVector = G4ThreeVector(minX,minY,-fDz);
>> 2056 fMaxBBoxVector = G4ThreeVector(maxX,maxY, fDz);
1476 } 2057 }
1477 2058
1478 ///////////////////////////////////////////// << 2059 // --------------------------------------------------------------------
1479 // << 2060
1480 // GetPolyhedron <<
1481 // <<
1482 G4Polyhedron* G4GenericTrap::GetPolyhedron () 2061 G4Polyhedron* G4GenericTrap::GetPolyhedron () const
1483 { 2062 {
1484 if ( (fpPolyhedron == nullptr) << 2063
1485 || fRebuildPolyhedron << 2064 #ifdef G4TESS_TEST
>> 2065 if ( fTessellatedSolid )
>> 2066 {
>> 2067 return fTessellatedSolid->GetPolyhedron();
>> 2068 }
>> 2069 #endif
>> 2070
>> 2071 if ( (!fpPolyhedron)
1486 || (fpPolyhedron->GetNumberOfRotationStep 2072 || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1487 fpPolyhedron->GetNumberOfRotationStep 2073 fpPolyhedron->GetNumberOfRotationSteps()) )
1488 { 2074 {
1489 G4AutoLock l(&polyhedronMutex); <<
1490 delete fpPolyhedron; 2075 delete fpPolyhedron;
1491 fpPolyhedron = CreatePolyhedron(); 2076 fpPolyhedron = CreatePolyhedron();
1492 fRebuildPolyhedron = false; <<
1493 l.unlock(); <<
1494 } 2077 }
1495 return fpPolyhedron; 2078 return fpPolyhedron;
1496 } << 2079 }
>> 2080
>> 2081 // --------------------------------------------------------------------
1497 2082
1498 ///////////////////////////////////////////// <<
1499 // <<
1500 // Method for visualisation <<
1501 // <<
1502 void G4GenericTrap::DescribeYourselfTo(G4VGra 2083 void G4GenericTrap::DescribeYourselfTo(G4VGraphicsScene& scene) const
1503 { 2084 {
>> 2085
>> 2086 #ifdef G4TESS_TEST
>> 2087 if ( fTessellatedSolid )
>> 2088 {
>> 2089 return fTessellatedSolid->DescribeYourselfTo(scene);
>> 2090 }
>> 2091 #endif
>> 2092
1504 scene.AddSolid(*this); 2093 scene.AddSolid(*this);
1505 } 2094 }
1506 2095
1507 ///////////////////////////////////////////// << 2096 // --------------------------------------------------------------------
1508 // << 2097
1509 // Return VisExtent <<
1510 // <<
1511 G4VisExtent G4GenericTrap::GetExtent() const 2098 G4VisExtent G4GenericTrap::GetExtent() const
1512 { 2099 {
1513 return { fMinBBox.x(), fMaxBBox.x(), << 2100 // Computes bounding vectors for the shape
1514 fMinBBox.y(), fMaxBBox.y(), << 2101
1515 fMinBBox.z(), fMaxBBox.z() }; << 2102 #ifdef G4TESS_TEST
1516 } << 2103 if ( fTessellatedSolid )
>> 2104 {
>> 2105 return fTessellatedSolid->GetExtent();
>> 2106 }
>> 2107 #endif
>> 2108
>> 2109 G4double Dx,Dy;
>> 2110 G4ThreeVector minVec = GetMinimumBBox();
>> 2111 G4ThreeVector maxVec = GetMaximumBBox();
>> 2112 Dx = 0.5*(maxVec.x()- minVec.x());
>> 2113 Dy = 0.5*(maxVec.y()- minVec.y());
>> 2114
>> 2115 return G4VisExtent (-Dx, Dx, -Dy, Dy, -fDz, fDz);
>> 2116 }
1517 2117
1518 // ------------------------------------------ 2118 // --------------------------------------------------------------------
1519 2119
1520 G4Polyhedron* G4GenericTrap::CreatePolyhedron 2120 G4Polyhedron* G4GenericTrap::CreatePolyhedron() const
1521 { 2121 {
>> 2122
>> 2123 #ifdef G4TESS_TEST
>> 2124 if ( fTessellatedSolid )
>> 2125 {
>> 2126 return fTessellatedSolid->CreatePolyhedron();
>> 2127 }
>> 2128 #endif
>> 2129
1522 // Approximation of Twisted Side 2130 // Approximation of Twisted Side
1523 // Construct extra Points, if Twisted Side 2131 // Construct extra Points, if Twisted Side
1524 // 2132 //
1525 G4Polyhedron* polyhedron; << 2133 G4PolyhedronArbitrary* polyhedron;
1526 G4int nVertices, nFacets; << 2134 size_t nVertices, nFacets;
1527 2135
1528 G4int subdivisions = 0; << 2136 G4int subdivisions=0;
1529 if (fIsTwisted) << 2137 G4int i;
>> 2138 if(fIsTwisted)
1530 { 2139 {
1531 if (GetVisSubdivisions() != 0) << 2140 if ( GetVisSubdivisions()!= 0 )
1532 { 2141 {
1533 subdivisions = GetVisSubdivisions(); << 2142 subdivisions=GetVisSubdivisions();
1534 } 2143 }
1535 else 2144 else
1536 { 2145 {
1537 // Estimation of Number of Subdivisions 2146 // Estimation of Number of Subdivisions for smooth visualisation
1538 // 2147 //
1539 G4double maxTwist = 0.; << 2148 G4double maxTwist=0.;
1540 for(G4int i = 0; i < 4; ++i) << 2149 for(i=0; i<4; i++)
1541 { 2150 {
1542 if (GetTwistAngle(i) > maxTwist) { ma << 2151 if(GetTwistAngle(i)>maxTwist) { maxTwist=GetTwistAngle(i); }
1543 } 2152 }
1544 2153
1545 // Computes bounding vectors for the sh 2154 // Computes bounding vectors for the shape
1546 // 2155 //
1547 G4double Dx, Dy; << 2156 G4double Dx,Dy;
1548 Dx = 0.5*(fMaxBBox.x() - fMinBBox.y()); << 2157 G4ThreeVector minVec = GetMinimumBBox();
1549 Dy = 0.5*(fMaxBBox.y() - fMinBBox.y()); << 2158 G4ThreeVector maxVec = GetMaximumBBox();
1550 if (Dy > Dx) { Dx = Dy; } << 2159 Dx = 0.5*(maxVec.x()- minVec.y());
>> 2160 Dy = 0.5*(maxVec.y()- minVec.y());
>> 2161 if (Dy > Dx) { Dx=Dy; }
>> 2162
>> 2163 subdivisions=8*G4int(maxTwist/(Dx*Dx*Dx)*fDz);
>> 2164 if (subdivisions<4) { subdivisions=4; }
>> 2165 if (subdivisions>30) { subdivisions=30; }
>> 2166 }
>> 2167 }
>> 2168 G4int sub4=4*subdivisions;
>> 2169 nVertices = 8+subdivisions*4;
>> 2170 nFacets = 6+subdivisions*4;
>> 2171 G4double cf=1./(subdivisions+1);
>> 2172 polyhedron = new G4PolyhedronArbitrary (nVertices, nFacets);
1551 2173
1552 subdivisions = 8*G4int(maxTwist/(Dx*Dx* << 2174 // Add Vertex
1553 if (subdivisions < 4) { subdivisions = <<
1554 if (subdivisions > 30) { subdivisions = <<
1555 } <<
1556 } <<
1557 G4int sub4 = 4*subdivisions; <<
1558 nVertices = 8 + subdivisions*4; <<
1559 nFacets = 6 + subdivisions*4; <<
1560 G4double cf = 1./(subdivisions + 1); <<
1561 polyhedron = new G4Polyhedron(nVertices, nF <<
1562 <<
1563 // Set vertices <<
1564 // 2175 //
1565 G4int icur = 0; << 2176 for (i=0;i<4;i++)
1566 for (G4int i = 0; i < 4; ++i) <<
1567 { 2177 {
1568 G4ThreeVector v(fVertices[i].x(),fVertice << 2178 polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
1569 polyhedron->SetVertex(++icur, v); << 2179 fVertices[i].y(),-fDz));
1570 } 2180 }
1571 for (G4int i = 0; i < subdivisions; ++i) << 2181 for( i=0;i<subdivisions;i++)
1572 { 2182 {
1573 for (G4int j = 0; j < 4; ++j) << 2183 for(G4int j=0;j<4;j++)
1574 { 2184 {
1575 G4TwoVector u = fVertices[j]+cf*(i+1)*( << 2185 G4TwoVector u=fVertices[j]+cf*(i+1)*( fVertices[j+4]-fVertices[j]);
1576 G4ThreeVector v(u.x(),u.y(),-fDz+cf*2*f << 2186 polyhedron->AddVertex(G4ThreeVector(u.x(),u.y(),-fDz+cf*2*fDz*(i+1)));
1577 polyhedron->SetVertex(++icur, v); << 2187 }
1578 } <<
1579 } 2188 }
1580 for (G4int i = 4; i < 8; ++i) << 2189 for (i=4;i<8;i++)
1581 { 2190 {
1582 G4ThreeVector v(fVertices[i].x(),fVertice << 2191 polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
1583 polyhedron->SetVertex(++icur, v); << 2192 fVertices[i].y(),fDz));
1584 } 2193 }
1585 2194
1586 // Set facets << 2195 // Add Facets
1587 // 2196 //
1588 icur = 0; << 2197 polyhedron->AddFacet(1,4,3,2); //Z-plane
1589 polyhedron->SetFacet(++icur, 1, 4, 3, 2); / << 2198 for (i=0;i<subdivisions+1;i++)
1590 for (G4int i = 0; i < subdivisions + 1; ++i <<
1591 { 2199 {
1592 G4int is = i*4; << 2200 G4int is=i*4;
1593 polyhedron->SetFacet(++icur, 5+is, 8+is, << 2201 polyhedron->AddFacet(5+is,8+is,4+is,1+is);
1594 polyhedron->SetFacet(++icur, 8+is, 7+is, << 2202 polyhedron->AddFacet(8+is,7+is,3+is,4+is);
1595 polyhedron->SetFacet(++icur, 7+is, 6+is, << 2203 polyhedron->AddFacet(7+is,6+is,2+is,3+is);
1596 polyhedron->SetFacet(++icur, 6+is, 5+is, << 2204 polyhedron->AddFacet(6+is,5+is,1+is,2+is);
1597 } 2205 }
1598 polyhedron->SetFacet(++icur, 5+sub4, 6+sub4 << 2206 polyhedron->AddFacet(5+sub4,6+sub4,7+sub4,8+sub4); //Z-plane
1599 2207
1600 polyhedron->SetReferences(); 2208 polyhedron->SetReferences();
1601 polyhedron->InvertFacets(); 2209 polyhedron->InvertFacets();
1602 2210
1603 return polyhedron; << 2211 return (G4Polyhedron*) polyhedron;
1604 } 2212 }
1605 2213
1606 ///////////////////////////////////////////// << 2214 // --------------------------------------------------------------------
1607 // <<
1608 // Print out a warning if A has an unexpected <<
1609 // <<
1610 void G4GenericTrap::WarningSignA(const G4Stri <<
1611 const G4Thre <<
1612 { <<
1613 std::ostringstream message; <<
1614 message.precision(16); <<
1615 message << icase << " in " << GetName() << <<
1616 << " p" << p << "\n" <<
1617 << " v" << v << "\n" <<
1618 << " A = " << A << " (is " <<
1619 << ((A < 0) ? "negative, instead of <<
1620 StreamInfo(message); <<
1621 const G4String function = "G4GenericTrap::D <<
1622 G4Exception(function, "GeomSolids1002", Jus <<
1623 } <<
1624 2215
1625 ///////////////////////////////////////////// << 2216 G4NURBS* G4GenericTrap::CreateNURBS() const
1626 // << 2217 {
1627 // Print out a warning if B has an unexpected << 2218 #ifdef G4TESS_TEST
1628 // << 2219 if ( fTessellatedSolid )
1629 void G4GenericTrap::WarningSignB(const G4Stri << 2220 {
1630 G4double f, << 2221 return fTessellatedSolid->CreateNURBS();
1631 const G4Thre << 2222 }
1632 { << 2223 #endif
1633 std::ostringstream message; <<
1634 message.precision(16); <<
1635 message << icase << " in " << GetName() << <<
1636 << " p" << p << "\n" <<
1637 << " v" << v << "\n" <<
1638 << " f = " << f << " B = " << B < <<
1639 << ((B < 0) ? "negative, instead of <<
1640 StreamInfo(message); <<
1641 const G4String function = "G4GenericTrap::D <<
1642 G4Exception(function, "GeomSolids1002", Jus <<
1643 } <<
1644 2224
1645 ///////////////////////////////////////////// << 2225 // Computes bounding vectors for the shape
1646 // << 2226 //
1647 // Print out a warning in DistanceToIn(p,v) << 2227 G4double Dx,Dy;
1648 // << 2228 G4ThreeVector minVec = GetMinimumBBox();
1649 void G4GenericTrap::WarningDistanceToIn(G4int << 2229 G4ThreeVector maxVec = GetMaximumBBox();
1650 G4dou << 2230 Dx = 0.5*(maxVec.x()- minVec.y());
1651 const << 2231 Dy = 0.5*(maxVec.y()- minVec.y());
1652 { <<
1653 G4String check = ""; <<
1654 if (ttin[1] != DBL_MAX) <<
1655 { <<
1656 G4double tcheck = 0.5*(ttout[0] + ttin[1] <<
1657 if (Inside(p + v*tcheck) != kOutside) che <<
1658 } <<
1659 <<
1660 auto position = Inside(p); <<
1661 std::ostringstream message; <<
1662 message.precision(16); <<
1663 message << k << "_Unexpected sequence of in <<
1664 << " position = " << ((position = <<
1665 << " p" << p << "\n" <<
1666 << " v" << v << "\n" <<
1667 << " range : [" << tmin << ", <<
1668 << " ttin[2] : " <<
1669 << ((ttin[0] == DBL_MAX) ? kInfinit <<
1670 << ((ttin[1] == DBL_MAX) ? kInfinit <<
1671 << " ttout[2] : " <<
1672 << ((ttout[0] == DBL_MAX) ? kInfini <<
1673 << ((ttout[1] == DBL_MAX) ? kInfini <<
1674 StreamInfo(message); <<
1675 G4Exception("G4GenericTrap::DistanceToIn(p, <<
1676 JustWarning, message ); <<
1677 } <<
1678 2232
1679 ///////////////////////////////////////////// << 2233 return new G4NURBSbox (Dx, Dy, fDz);
1680 // << 2234 }
1681 // Print out a warning in DistanceToOut(p,v) <<
1682 // <<
1683 void G4GenericTrap::WarningDistanceToOut(cons <<
1684 cons <<
1685 G4do <<
1686 { <<
1687 auto position = Inside(p); <<
1688 std::ostringstream message; <<
1689 message.precision(16); <<
1690 message << "Unexpected final tout = " << to <<
1691 << " position = " << ((position = <<
1692 << " p" << p << "\n" <<
1693 << " v" << v << "\n"; <<
1694 StreamInfo(message); <<
1695 G4Exception("G4GenericTrap::DistanceToOut(p <<
1696 JustWarning, message ); <<
1697 } <<
1698 2235
1699 #endif << 2236 // --------------------------------------------------------------------
1700 2237