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