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