| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // Implementation of G4Polyhedra, a CSG polyhe << 23 //
27 // as an inherited class of G4VCSGfaceted. << 24 // $Id: G4Polyhedra.cc,v 1.2.4.1 2001/06/28 19:09:26 gunter Exp $
>> 25 // GEANT4 tag $Name: $
>> 26 //
>> 27 //
>> 28 // --------------------------------------------------------------------
>> 29 // GEANT 4 class source file
>> 30 //
>> 31 //
>> 32 // G4Polyhedra.cc
>> 33 //
>> 34 // Implementation of a CSG polyhedra, as an inherited class of G4VCSGfaceted.
>> 35 //
>> 36 // To be done:
>> 37 // * Cracks: there are probably small cracks in the seams between the
>> 38 // phi face (G4PolyPhiFace) and sides (G4PolyhedraSide) that are not
>> 39 // entirely leakproof. Also, I am not sure all vertices are leak proof.
>> 40 // * Many optimizations are possible, but not implemented.
>> 41 // * Visualization needs to be updated outside of this routine.
28 // 42 //
29 // Utility classes: 43 // Utility classes:
30 // * G4EnclosingCylinder: decided a quick c << 44 // * G4EnclosingCylinder: I decided a quick check of geometry would be a
31 // good idea (for CPU speed). If the quic 45 // good idea (for CPU speed). If the quick check fails, the regular
32 // full-blown G4VCSGfaceted version is in 46 // full-blown G4VCSGfaceted version is invoked.
33 // * G4ReduciblePolygon: Really meant as a 47 // * G4ReduciblePolygon: Really meant as a check of input parameters,
34 // this utility class also "converts" the 48 // this utility class also "converts" the GEANT3-like PGON/PCON
35 // arguments into the newer ones. 49 // arguments into the newer ones.
36 // Both these classes are implemented outside 50 // Both these classes are implemented outside this file because they are
37 // shared with G4Polycone. 51 // shared with G4Polycone.
38 // 52 //
39 // Author: David C. Williams (davidw@scipp.ucs <<
40 // ------------------------------------------- 53 // --------------------------------------------------------------------
41 54
42 #include "G4Polyhedra.hh" 55 #include "G4Polyhedra.hh"
43 <<
44 #if !defined(G4GEOM_USE_UPOLYHEDRA) <<
45 <<
46 #include "G4PolyhedraSide.hh" 56 #include "G4PolyhedraSide.hh"
47 #include "G4PolyPhiFace.hh" 57 #include "G4PolyPhiFace.hh"
48 58
49 #include "G4GeomTools.hh" << 59 #include "G4Polyhedron.hh"
50 #include "G4VoxelLimits.hh" <<
51 #include "G4AffineTransform.hh" <<
52 #include "G4BoundingEnvelope.hh" <<
53 <<
54 #include "G4QuickRand.hh" <<
55 <<
56 #include "G4EnclosingCylinder.hh" 60 #include "G4EnclosingCylinder.hh"
57 #include "G4ReduciblePolygon.hh" 61 #include "G4ReduciblePolygon.hh"
58 #include "G4VPVParameterisation.hh" <<
59 <<
60 namespace <<
61 { <<
62 G4Mutex surface_elementsMutex = G4MUTEX_INIT <<
63 } <<
64 <<
65 using namespace CLHEP; <<
66 62
>> 63 //
67 // Constructor (GEANT3 style parameters) 64 // Constructor (GEANT3 style parameters)
68 // 65 //
69 // GEANT3 PGON radii are specified in the dist 66 // GEANT3 PGON radii are specified in the distance to the norm of each face.
70 // << 67 //
71 G4Polyhedra::G4Polyhedra( const G4String& name << 68 G4Polyhedra::G4Polyhedra( const G4String& name,
72 G4double phiSt 69 G4double phiStart,
73 G4double thePh 70 G4double thePhiTotal,
74 G4int theNumSi << 71 G4int theNumSide,
75 G4int numZPlan 72 G4int numZPlanes,
76 const G4double zPlan 73 const G4double zPlane[],
77 const G4double rInne 74 const G4double rInner[],
78 const G4double rOute << 75 const G4double rOuter[] ) : G4VCSGfaceted( name )
79 : G4VCSGfaceted( name ) <<
80 { 76 {
81 if (theNumSide <= 0) << 77 if (theNumSide <= 0) G4Exception( "G4Polyhedra:: must have at least one side" );
82 { <<
83 std::ostringstream message; <<
84 message << "Solid must have at least one s <<
85 << " No sides specified !"; <<
86 G4Exception("G4Polyhedra::G4Polyhedra()", <<
87 FatalErrorInArgument, message) <<
88 } <<
89 <<
90 // <<
91 // Calculate conversion factor from G3 radiu <<
92 // <<
93 G4double phiTotal = thePhiTotal; <<
94 if ( (phiTotal <=0) || (phiTotal >= twopi*(1 <<
95 { phiTotal = twopi; } <<
96 G4double convertRad = std::cos(0.5*phiTotal/ <<
97 <<
98 // <<
99 // Some historical stuff <<
100 // <<
101 original_parameters = new G4PolyhedraHistori <<
102 <<
103 original_parameters->numSide = theNumSide; <<
104 original_parameters->Start_angle = phiStart; <<
105 original_parameters->Opening_angle = phiTota <<
106 original_parameters->Num_z_planes = numZPlan <<
107 original_parameters->Z_values = new G4double <<
108 original_parameters->Rmin = new G4double[num <<
109 original_parameters->Rmax = new G4double[num <<
110 <<
111 for (G4int i=0; i<numZPlanes; ++i) <<
112 { <<
113 if (( i < numZPlanes-1) && ( zPlane[i] == <<
114 { <<
115 if( (rInner[i] > rOuter[i+1]) <<
116 ||(rInner[i+1] > rOuter[i]) ) <<
117 { <<
118 DumpInfo(); <<
119 std::ostringstream message; <<
120 message << "Cannot create a Polyhedra <<
121 << G4endl <<
122 << " Segments are not c <<
123 << " rMin[" << i << "] <<
124 << " -- rMax[" << i+1 << "] = <<
125 << " rMin[" << i+1 << " <<
126 << " -- rMax[" << i << "] = " <<
127 G4Exception("G4Polyhedra::G4Polyhedra( <<
128 FatalErrorInArgument, mess <<
129 } <<
130 } <<
131 original_parameters->Z_values[i] = zPlane[ <<
132 original_parameters->Rmin[i] = rInner[i]/c <<
133 original_parameters->Rmax[i] = rOuter[i]/c <<
134 } <<
135 <<
136 <<
137 // <<
138 // Build RZ polygon using special PCON/PGON <<
139 // <<
140 auto rz = new G4ReduciblePolygon( rInner, rO <<
141 rz->ScaleA( 1/convertRad ); <<
142 <<
143 // <<
144 // Do the real work <<
145 // <<
146 Create( phiStart, phiTotal, theNumSide, rz ) <<
147 78
148 delete rz; << 79 //
>> 80 // Calculate conversion factor from G3 radius to G4 radius
>> 81 //
>> 82 G4double phiTotal = thePhiTotal;
>> 83 if (phiTotal <=0 || phiTotal >= 2*M_PI*(1-DBL_EPSILON)) phiTotal = 2*M_PI;
>> 84 G4double convertRad = cos(0.5*phiTotal/theNumSide);
>> 85
>> 86 //
>> 87 // Some historical stuff
>> 88 //
>> 89 original_parameters = new G4PolyhedraHistorical;
>> 90
>> 91 original_parameters->numSide = theNumSide;
>> 92 original_parameters->Start_angle = phiStart;
>> 93 original_parameters->Opening_angle = phiTotal;
>> 94 original_parameters->Num_z_planes = numZPlanes;
>> 95 original_parameters->Z_values = new G4double[numZPlanes];
>> 96 original_parameters->Rmin = new G4double[numZPlanes];
>> 97 original_parameters->Rmax = new G4double[numZPlanes];
>> 98
>> 99 G4int i;
>> 100 for (i=0; i<numZPlanes; i++) {
>> 101 original_parameters->Z_values[i] = zPlane[i];
>> 102 original_parameters->Rmin[i] = rInner[i]/convertRad;
>> 103 original_parameters->Rmax[i] = rOuter[i]/convertRad;
>> 104 }
>> 105
>> 106
>> 107 //
>> 108 // Build RZ polygon using special PCON/PGON GEANT3 constructor
>> 109 //
>> 110 G4ReduciblePolygon *rz = new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
>> 111 rz->ScaleA( 1/convertRad );
>> 112
>> 113 //
>> 114 // Do the real work
>> 115 //
>> 116 Create( phiStart, phiTotal, theNumSide, rz );
>> 117
>> 118 delete rz;
149 } 119 }
150 120
>> 121
>> 122 //
151 // Constructor (generic parameters) 123 // Constructor (generic parameters)
152 // 124 //
153 G4Polyhedra::G4Polyhedra( const G4String& name << 125 G4Polyhedra::G4Polyhedra( const G4String& name,
154 G4double phiSt << 126 G4double phiStart,
155 G4double phiTo << 127 G4double phiTotal,
156 G4int theNu << 128 G4int theNumSide,
157 G4int numRZ << 129 G4int numRZ,
158 const G4double r[], << 130 const G4double r[],
159 const G4double z[] << 131 const G4double z[] ) : G4VCSGfaceted( name )
160 : G4VCSGfaceted( name ), genericPgon(true) << 132 {
161 { << 133 original_parameters = 0;
162 if (theNumSide <= 0) << 134
163 { << 135 G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
164 std::ostringstream message; << 136
165 message << "Solid must have at least one s << 137 Create( phiStart, phiTotal, theNumSide, rz );
166 << " No sides specified !"; << 138
167 G4Exception("G4Polyhedra::G4Polyhedra()", << 139 delete rz;
168 FatalErrorInArgument, message) <<
169 } <<
170 <<
171 auto rz = new G4ReduciblePolygon( r, z, numR <<
172 <<
173 Create( phiStart, phiTotal, theNumSide, rz ) <<
174 <<
175 // Set original_parameters struct for consis <<
176 // <<
177 SetOriginalParameters(rz); <<
178 <<
179 delete rz; <<
180 } 140 }
181 141
>> 142
>> 143 //
182 // Create 144 // Create
183 // 145 //
184 // Generic create routine, called by each cons << 146 // Generic create routine, called by each constructor after conversion of arguments
185 // after conversion of arguments <<
186 // 147 //
187 void G4Polyhedra::Create( G4double phiStart, 148 void G4Polyhedra::Create( G4double phiStart,
188 G4double phiTotal, << 149 G4double phiTotal,
189 G4int theNumSide, << 150 G4int theNumSide,
190 G4ReduciblePolygon* << 151 G4ReduciblePolygon *rz )
191 { << 152 {
192 // << 153 //
193 // Perform checks of rz values << 154 // Perform checks of rz values
194 // << 155 //
195 if (rz->Amin() < 0.0) << 156 if (rz->Amin() < 0.0)
196 { << 157 G4Exception( "G4Polyhedra: Illegal input parameters: All R values must be >= 0" );
197 std::ostringstream message; << 158
198 message << "Illegal input parameters - " < << 159 G4double rzArea = rz->Area();
199 << " All R values must be > << 160 if (rzArea < -kCarTolerance) rz->ReverseOrder();
200 G4Exception("G4Polyhedra::Create()", "Geom << 161
201 FatalErrorInArgument, message) << 162 else if (rzArea < -kCarTolerance)
202 } << 163 G4Exception( "G4Polyhedra: Illegal input parameters: R/Z cross section is zero or near zero" );
203 << 164
204 G4double rzArea = rz->Area(); << 165 if ((!rz->RemoveDuplicateVertices( kCarTolerance )) ||
205 if (rzArea < -kCarTolerance) << 166 (!rz->RemoveRedundantVertices( kCarTolerance )) )
206 { << 167 G4Exception( "G4Polyhedra: Illegal input parameters: Too few unique R/Z values" );
207 rz->ReverseOrder(); << 168
208 } << 169 if (rz->CrossesItself( 1/kInfinity ))
209 else if (rzArea < kCarTolerance) << 170 G4Exception( "G4Polyhedra: Illegal input parameters: R/Z segments cross" );
210 { << 171
211 std::ostringstream message; << 172 numCorner = rz->NumVertices();
212 message << "Illegal input parameters - " < << 173
213 << " R/Z cross section is z << 174
214 G4Exception("G4Polyhedra::Create()", "Geom << 175 startPhi = phiStart;
215 FatalErrorInArgument, message) << 176 while( startPhi < 0 ) startPhi += 2*M_PI;
216 } << 177 //
217 << 178 // Phi opening? Account for some possible roundoff, and interpret
218 if ( (!rz->RemoveDuplicateVertices( kCarTole << 179 // nonsense value as representing no phi opening
219 || (!rz->RemoveRedundantVertices( kCarTole << 180 //
220 { << 181 if (phiTotal <= 0 || phiTotal > 2.0*M_PI*(1-DBL_EPSILON)) {
221 std::ostringstream message; << 182 phiIsOpen = false;
222 message << "Illegal input parameters - " < << 183 endPhi = phiStart+2*M_PI;
223 << " Too few unique R/Z val << 184 }
224 G4Exception("G4Polyhedra::Create()", "Geom << 185 else {
225 FatalErrorInArgument, message) << 186 phiIsOpen = true;
226 } << 187
227 << 188 //
228 if (rz->CrossesItself( 1/kInfinity )) << 189 // Convert phi into our convention
229 { << 190 //
230 std::ostringstream message; << 191 endPhi = phiStart+phiTotal;
231 message << "Illegal input parameters - " < << 192 while( endPhi < startPhi ) endPhi += 2*M_PI;
232 << " R/Z segments cross !"; << 193 }
233 G4Exception("G4Polyhedra::Create()", "Geom << 194
234 FatalErrorInArgument, message) << 195 //
235 } << 196 // Save number sides
236 << 197 //
237 numCorner = rz->NumVertices(); << 198 numSide = theNumSide;
238 << 199
239 << 200 //
240 startPhi = phiStart; << 201 // Allocate corner array.
241 while( startPhi < 0 ) // Loop checking, 1 << 202 //
242 startPhi += twopi; << 203 corners = new G4PolyhedraSideRZ[numCorner];
243 // << 204
244 // Phi opening? Account for some possible ro << 205 //
245 // nonsense value as representing no phi ope << 206 // Copy corners
246 // << 207 //
247 if ( (phiTotal <= 0) || (phiTotal > twopi*(1 << 208 G4ReduciblePolygonIterator iterRZ(rz);
248 { << 209
249 phiIsOpen = false; << 210 G4PolyhedraSideRZ *next = corners;
250 endPhi = startPhi + twopi; << 211 iterRZ.Begin();
251 } << 212 do {
252 else << 213 next->r = iterRZ.GetA();
253 { << 214 next->z = iterRZ.GetB();
254 phiIsOpen = true; << 215 } while( ++next, iterRZ.Next() );
255 endPhi = startPhi + phiTotal; << 216
256 } << 217 //
257 << 218 // Allocate face pointer array
258 // << 219 //
259 // Save number sides << 220 numFace = phiIsOpen ? numCorner+2 : numCorner;
260 // << 221 faces = new G4VCSGface*[numFace];
261 numSide = theNumSide; << 222
262 << 223 //
263 // << 224 // Construct side faces
264 // Allocate corner array. << 225 //
265 // << 226 // To do so properly, we need to keep track of four successive RZ
266 corners = new G4PolyhedraSideRZ[numCorner]; << 227 // corners.
267 << 228 //
268 // << 229 // But! Don't construct a face if both points are at zero radius!
269 // Copy corners << 230 //
270 // << 231 G4PolyhedraSideRZ *corner = corners,
271 G4ReduciblePolygonIterator iterRZ(rz); << 232 *prev = corners + numCorner-1,
272 << 233 *nextNext;
273 G4PolyhedraSideRZ *next = corners; << 234 G4VCSGface **face = faces;
274 iterRZ.Begin(); << 235 do {
275 do // Loop checking, 13.08.2015, G.Cosmo << 236 next = corner+1;
276 { << 237 if (next >= corners+numCorner) next = corners;
277 next->r = iterRZ.GetA(); << 238 nextNext = next+1;
278 next->z = iterRZ.GetB(); << 239 if (nextNext >= corners+numCorner) nextNext = corners;
279 } while( ++next, iterRZ.Next() ); << 240
280 << 241 if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
281 // << 242
282 // Allocate face pointer array << 243 //
283 // << 244 // We must decide here if we can dare declare one of our faces
284 numFace = phiIsOpen ? numCorner+2 : numCorne << 245 // as having a "valid" normal (i.e. allBehind = true). This
285 faces = new G4VCSGface*[numFace]; << 246 // is never possible if the face faces "inward" in r *unless*
286 << 247 // we have only one side
287 // << 248 //
288 // Construct side faces << 249 G4bool allBehind;
289 // << 250 if ((corner->z > next->z) && (numSide > 1)) {
290 // To do so properly, we need to keep track << 251 allBehind = false;
291 // corners. << 252 }
292 // << 253 else {
293 // But! Don't construct a face if both point << 254 //
294 // << 255 // Otherwise, it is only true if the line passing
295 G4PolyhedraSideRZ* corner = corners, << 256 // through the two points of the segment do not
296 * prev = corners + numCorne << 257 // split the r/z cross section
297 * nextNext; << 258 //
298 G4VCSGface** face = faces; << 259 allBehind = !rz->BisectedBy( corner->r, corner->z,
299 do // Loop checking, 13.08.2015, G.Cosmo << 260 next->r, next->z, kCarTolerance );
300 { << 261 }
301 next = corner+1; << 262
302 if (next >= corners+numCorner) next = corn << 263 *face++ = new G4PolyhedraSide( prev, corner, next, nextNext,
303 nextNext = next+1; << 264 numSide, startPhi, endPhi-startPhi, phiIsOpen );
304 if (nextNext >= corners+numCorner) nextNex << 265 } while( prev=corner, corner=next, corner > corners );
305 << 266
306 if (corner->r < 1/kInfinity && next->r < 1 << 267 if (phiIsOpen) {
307 /* << 268 //
308 // We must decide here if we can dare decl << 269 // Construct phi open edges
309 // as having a "valid" normal (i.e. allBeh << 270 //
310 // is never possible if the face faces "in << 271 *face++ = new G4PolyPhiFace( rz, startPhi, phiTotal/numSide, endPhi );
311 // we have only one side << 272 *face++ = new G4PolyPhiFace( rz, endPhi, phiTotal/numSide, startPhi );
312 // << 273 }
313 G4bool allBehind; << 274
314 if ((corner->z > next->z) && (numSide > 1) << 275 //
315 { << 276 // We might have dropped a face or two: recalculate numFace
316 allBehind = false; << 277 //
317 } << 278 numFace = face-faces;
318 else << 279
319 { << 280 //
320 // << 281 // Make enclosingCylinder
321 // Otherwise, it is only true if the lin << 282 //
322 // through the two points of the segment << 283 enclosingCylinder = new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
323 // split the r/z cross section <<
324 // <<
325 allBehind = !rz->BisectedBy( corner->r, <<
326 next->r, ne <<
327 } <<
328 */ <<
329 *face++ = new G4PolyhedraSide( prev, corne <<
330 numSide, startPhi, endPhi-sta <<
331 } while( prev=corner, corner=next, corner > <<
332 <<
333 if (phiIsOpen) <<
334 { <<
335 // <<
336 // Construct phi open edges <<
337 // <<
338 *face++ = new G4PolyPhiFace( rz, startPhi, <<
339 *face++ = new G4PolyPhiFace( rz, endPhi, <<
340 } <<
341 <<
342 // <<
343 // We might have dropped a face or two: reca <<
344 // <<
345 numFace = (G4int)(face-faces); <<
346 <<
347 // <<
348 // Make enclosingCylinder <<
349 // <<
350 enclosingCylinder = <<
351 new G4EnclosingCylinder( rz, phiIsOpen, ph <<
352 } 284 }
353 285
354 // Fake default constructor - sets only member <<
355 // for usage restri <<
356 // <<
357 G4Polyhedra::G4Polyhedra( __void__& a ) <<
358 : G4VCSGfaceted(a), startPhi(0.), endPhi(0.) <<
359 { <<
360 } <<
361 286
>> 287 //
362 // Destructor 288 // Destructor
363 // 289 //
364 G4Polyhedra::~G4Polyhedra() 290 G4Polyhedra::~G4Polyhedra()
365 { 291 {
366 delete [] corners; << 292 delete [] corners;
367 delete original_parameters; << 293 if (original_parameters) delete original_parameters;
368 delete enclosingCylinder; << 294
369 delete fElements; << 295 delete enclosingCylinder;
370 delete fpPolyhedron; <<
371 corners = nullptr; <<
372 original_parameters = nullptr; <<
373 enclosingCylinder = nullptr; <<
374 fElements = nullptr; <<
375 fpPolyhedron = nullptr; <<
376 } 296 }
377 297
>> 298
>> 299 //
378 // Copy constructor 300 // Copy constructor
379 // 301 //
380 G4Polyhedra::G4Polyhedra( const G4Polyhedra& s << 302 G4Polyhedra::G4Polyhedra( const G4Polyhedra &source ) : G4VCSGfaceted( source )
381 : G4VCSGfaceted( source ) <<
382 { 303 {
383 CopyStuff( source ); << 304 CopyStuff( source );
384 } 305 }
385 306
>> 307
>> 308 //
386 // Assignment operator 309 // Assignment operator
387 // 310 //
388 G4Polyhedra &G4Polyhedra::operator=( const G4P << 311 const G4Polyhedra &G4Polyhedra::operator=( const G4Polyhedra &source )
389 { 312 {
390 if (this == &source) return *this; << 313 if (this == &source) return *this;
391 <<
392 G4VCSGfaceted::operator=( source ); <<
393 314
394 delete [] corners; << 315 G4VCSGfaceted::operator=( source );
395 delete original_parameters; << 316
396 delete enclosingCylinder; << 317 delete [] corners;
397 << 318 if (original_parameters) delete original_parameters;
398 CopyStuff( source ); << 319
399 << 320 delete enclosingCylinder;
400 return *this; << 321
>> 322 CopyStuff( source );
>> 323
>> 324 return *this;
401 } 325 }
402 326
403 // CopyStuff <<
404 // <<
405 void G4Polyhedra::CopyStuff( const G4Polyhedra <<
406 { <<
407 // <<
408 // Simple stuff <<
409 // <<
410 numSide = source.numSide; <<
411 startPhi = source.startPhi; <<
412 endPhi = source.endPhi; <<
413 phiIsOpen = source.phiIsOpen; <<
414 numCorner = source.numCorner; <<
415 genericPgon= source.genericPgon; <<
416 <<
417 // <<
418 // The corner array <<
419 // <<
420 corners = new G4PolyhedraSideRZ[numCorner]; <<
421 <<
422 G4PolyhedraSideRZ* corn = corners, <<
423 * sourceCorn = source.corne <<
424 do // Loop checking, 13.08.2015, G.Cosmo <<
425 { <<
426 *corn = *sourceCorn; <<
427 } while( ++sourceCorn, ++corn < corners+numC <<
428 <<
429 // <<
430 // Original parameters <<
431 // <<
432 if (source.original_parameters != nullptr) <<
433 { <<
434 original_parameters = <<
435 new G4PolyhedraHistorical( *source.origi <<
436 } <<
437 <<
438 // <<
439 // Enclosing cylinder <<
440 // <<
441 enclosingCylinder = new G4EnclosingCylinder( <<
442 <<
443 // <<
444 // Surface elements <<
445 // <<
446 delete fElements; <<
447 fElements = nullptr; <<
448 <<
449 // <<
450 // Polyhedron <<
451 // <<
452 fRebuildPolyhedron = false; <<
453 delete fpPolyhedron; <<
454 fpPolyhedron = nullptr; <<
455 } <<
456 327
457 // Reset <<
458 // 328 //
459 // Recalculates and reshapes the solid, given << 329 // CopyStuff
460 // original_parameters. <<
461 // 330 //
462 G4bool G4Polyhedra::Reset() << 331 void G4Polyhedra::CopyStuff( const G4Polyhedra &source )
463 { 332 {
464 if (genericPgon) << 333 //
465 { << 334 // Simple stuff
466 std::ostringstream message; << 335 //
467 message << "Solid " << GetName() << " buil << 336 numSide = source.numSide;
468 << G4endl << "Not applicable to th << 337 startPhi = source.startPhi;
469 G4Exception("G4Polyhedra::Reset()", "GeomS << 338 endPhi = source.endPhi;
470 JustWarning, message, "Paramet << 339 phiIsOpen = source.phiIsOpen;
471 return true; << 340 numCorner = source.numCorner;
472 } << 341
473 << 342 //
474 // << 343 // The corner array
475 // Clear old setup << 344 //
476 // << 345 corners = new G4PolyhedraSideRZ[numCorner];
477 G4VCSGfaceted::DeleteStuff(); << 346
478 delete [] corners; << 347 G4PolyhedraSideRZ *corn = corners,
479 delete enclosingCylinder; << 348 *sourceCorn = source.corners;
480 delete fElements; << 349 do {
481 corners = nullptr; << 350 *corn = *sourceCorn;
482 fElements = nullptr; << 351 } while( ++sourceCorn, ++corn < corners+numCorner );
483 enclosingCylinder = nullptr; << 352
484 << 353 //
485 // << 354 // Original parameters
486 // Rebuild polyhedra << 355 //
487 // << 356 if (source.original_parameters) {
488 auto rz = new G4ReduciblePolygon( original_p << 357 original_parameters = new G4PolyhedraHistorical( *source.original_parameters );
489 original_p << 358 }
490 original_p << 359
491 original_p << 360 //
492 Create( original_parameters->Start_angle, << 361 // Enclosing cylinder
493 original_parameters->Opening_angle, << 362 //
494 original_parameters->numSide, rz ); << 363 enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder );
495 delete rz; <<
496 <<
497 return false; <<
498 } 364 }
499 365
>> 366
>> 367 //
500 // Inside 368 // Inside
501 // 369 //
502 // This is an override of G4VCSGfaceted::Insid << 370 // This is an override of G4VCSGfaceted::Inside, created in order to speed things
503 // to speed things up by first checking with G << 371 // up by first checking with G4EnclosingCylinder.
504 // 372 //
505 EInside G4Polyhedra::Inside( const G4ThreeVect << 373 EInside G4Polyhedra::Inside( const G4ThreeVector &p ) const
506 { 374 {
507 // << 375 //
508 // Quick test << 376 // Quick test
509 // << 377 //
510 if (enclosingCylinder->MustBeOutside(p)) ret << 378 if (enclosingCylinder->MustBeOutside(p)) return kOutside;
511 << 379
512 // << 380 //
513 // Long answer << 381 // Long answer
514 // << 382 //
515 return G4VCSGfaceted::Inside(p); << 383 return G4VCSGfaceted::Inside(p);
516 } 384 }
517 385
518 // DistanceToIn <<
519 // <<
520 // This is an override of G4VCSGfaceted::Insid <<
521 // to speed things up by first checking with G <<
522 // <<
523 G4double G4Polyhedra::DistanceToIn( const G4Th <<
524 const G4Th <<
525 { <<
526 // <<
527 // Quick test <<
528 // <<
529 if (enclosingCylinder->ShouldMiss(p,v)) <<
530 return kInfinity; <<
531 <<
532 // <<
533 // Long answer <<
534 // <<
535 return G4VCSGfaceted::DistanceToIn( p, v ); <<
536 } <<
537 386
>> 387 //
538 // DistanceToIn 388 // DistanceToIn
539 // 389 //
540 G4double G4Polyhedra::DistanceToIn( const G4Th << 390 // This is an override of G4VCSGfaceted::Inside, created in order to speed things
541 { << 391 // up by first checking with G4EnclosingCylinder.
542 return G4VCSGfaceted::DistanceToIn(p); <<
543 } <<
544 <<
545 // Get bounding box <<
546 // 392 //
547 void G4Polyhedra::BoundingLimits(G4ThreeVector << 393 G4double G4Polyhedra::DistanceToIn( const G4ThreeVector &p, const G4ThreeVector &v ) const
548 G4ThreeVector <<
549 { 394 {
550 G4double rmin = kInfinity, rmax = -kInfinity << 395 //
551 G4double zmin = kInfinity, zmax = -kInfinity << 396 // Quick test
552 for (G4int i=0; i<GetNumRZCorner(); ++i) << 397 //
553 { << 398 if (enclosingCylinder->ShouldMiss(p,v)) return kInfinity;
554 G4PolyhedraSideRZ corner = GetCorner(i); << 399
555 if (corner.r < rmin) rmin = corner.r; << 400 //
556 if (corner.r > rmax) rmax = corner.r; << 401 // Long answer
557 if (corner.z < zmin) zmin = corner.z; << 402 //
558 if (corner.z > zmax) zmax = corner.z; << 403 return G4VCSGfaceted::DistanceToIn( p, v );
559 } <<
560 <<
561 G4double sphi = GetStartPhi(); <<
562 G4double ephi = GetEndPhi(); <<
563 G4double dphi = IsOpen() ? ephi-sphi : tw <<
564 G4int ksteps = GetNumSide(); <<
565 G4double astep = dphi/ksteps; <<
566 G4double sinStep = std::sin(astep); <<
567 G4double cosStep = std::cos(astep); <<
568 <<
569 G4double sinCur = GetSinStartPhi(); <<
570 G4double cosCur = GetCosStartPhi(); <<
571 if (!IsOpen()) rmin = 0.; <<
572 G4double xmin = rmin*cosCur, xmax = xmin; <<
573 G4double ymin = rmin*sinCur, ymax = ymin; <<
574 for (G4int k=0; k<ksteps+1; ++k) <<
575 { <<
576 G4double x = rmax*cosCur; <<
577 if (x < xmin) xmin = x; <<
578 if (x > xmax) xmax = x; <<
579 G4double y = rmax*sinCur; <<
580 if (y < ymin) ymin = y; <<
581 if (y > ymax) ymax = y; <<
582 if (rmin > 0) <<
583 { <<
584 G4double xx = rmin*cosCur; <<
585 if (xx < xmin) xmin = xx; <<
586 if (xx > xmax) xmax = xx; <<
587 G4double yy = rmin*sinCur; <<
588 if (yy < ymin) ymin = yy; <<
589 if (yy > ymax) ymax = yy; <<
590 } <<
591 G4double sinTmp = sinCur; <<
592 sinCur = sinCur*cosStep + cosCur*sinStep; <<
593 cosCur = cosCur*cosStep - sinTmp*sinStep; <<
594 } <<
595 pMin.set(xmin,ymin,zmin); <<
596 pMax.set(xmax,ymax,zmax); <<
597 <<
598 // Check correctness of the bounding box <<
599 // <<
600 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
601 { <<
602 std::ostringstream message; <<
603 message << "Bad bounding box (min >= max) <<
604 << GetName() << " !" <<
605 << "\npMin = " << pMin <<
606 << "\npMax = " << pMax; <<
607 G4Exception("G4Polyhedra::BoundingLimits() <<
608 JustWarning, message); <<
609 DumpInfo(); <<
610 } <<
611 } 404 }
612 405
613 // Calculate extent under transform and specif <<
614 // <<
615 G4bool G4Polyhedra::CalculateExtent(const EAxi <<
616 const G4Vo <<
617 const G4Af <<
618 G4double& <<
619 { <<
620 G4ThreeVector bmin, bmax; <<
621 G4bool exist; <<
622 <<
623 // Check bounding box (bbox) <<
624 // <<
625 BoundingLimits(bmin,bmax); <<
626 G4BoundingEnvelope bbox(bmin,bmax); <<
627 #ifdef G4BBOX_EXTENT <<
628 return bbox.CalculateExtent(pAxis,pVoxelLimi <<
629 #endif <<
630 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox <<
631 { <<
632 return exist = pMin < pMax; <<
633 } <<
634 <<
635 // To find the extent, RZ contour of the pol <<
636 // in triangles. The extent is calculated as <<
637 // all sub-polycones formed by rotation of t <<
638 // <<
639 G4TwoVectorList contourRZ; <<
640 G4TwoVectorList triangles; <<
641 std::vector<G4int> iout; <<
642 G4double eminlim = pVoxelLimit.GetMinExtent( <<
643 G4double emaxlim = pVoxelLimit.GetMaxExtent( <<
644 <<
645 // get RZ contour, ensure anticlockwise orde <<
646 for (G4int i=0; i<GetNumRZCorner(); ++i) <<
647 { <<
648 G4PolyhedraSideRZ corner = GetCorner(i); <<
649 contourRZ.emplace_back(corner.r,corner.z); <<
650 } <<
651 G4GeomTools::RemoveRedundantVertices(contour <<
652 G4double area = G4GeomTools::PolygonArea(con <<
653 if (area < 0.) std::reverse(contourRZ.begin( <<
654 <<
655 // triangulate RZ countour <<
656 if (!G4GeomTools::TriangulatePolygon(contour <<
657 { <<
658 std::ostringstream message; <<
659 message << "Triangulation of RZ contour ha <<
660 << GetName() << " !" <<
661 << "\nExtent has been calculated u <<
662 G4Exception("G4Polyhedra::CalculateExtent( <<
663 "GeomMgt1002",JustWarning,mess <<
664 return bbox.CalculateExtent(pAxis,pVoxelLi <<
665 } <<
666 <<
667 // set trigonometric values <<
668 G4double sphi = GetStartPhi(); <<
669 G4double ephi = GetEndPhi(); <<
670 G4double dphi = IsOpen() ? ephi-sphi : t <<
671 G4int ksteps = GetNumSide(); <<
672 G4double astep = dphi/ksteps; <<
673 G4double sinStep = std::sin(astep); <<
674 G4double cosStep = std::cos(astep); <<
675 G4double sinStart = GetSinStartPhi(); <<
676 G4double cosStart = GetCosStartPhi(); <<
677 <<
678 // allocate vector lists <<
679 std::vector<const G4ThreeVectorList *> polyg <<
680 polygons.resize(ksteps+1); <<
681 for (G4int k=0; k<ksteps+1; ++k) <<
682 { <<
683 polygons[k] = new G4ThreeVectorList(3); <<
684 } <<
685 <<
686 // main loop along triangles <<
687 pMin = kInfinity; <<
688 pMax = -kInfinity; <<
689 G4int ntria = (G4int)triangles.size()/3; <<
690 for (G4int i=0; i<ntria; ++i) <<
691 { <<
692 G4double sinCur = sinStart; <<
693 G4double cosCur = cosStart; <<
694 G4int i3 = i*3; <<
695 for (G4int k=0; k<ksteps+1; ++k) // rotate <<
696 { <<
697 auto ptr = const_cast<G4ThreeVectorList* <<
698 auto iter = ptr->begin(); <<
699 iter->set(triangles[i3+0].x()*cosCur, <<
700 triangles[i3+0].x()*sinCur, <<
701 triangles[i3+0].y()); <<
702 iter++; <<
703 iter->set(triangles[i3+1].x()*cosCur, <<
704 triangles[i3+1].x()*sinCur, <<
705 triangles[i3+1].y()); <<
706 iter++; <<
707 iter->set(triangles[i3+2].x()*cosCur, <<
708 triangles[i3+2].x()*sinCur, <<
709 triangles[i3+2].y()); <<
710 <<
711 G4double sinTmp = sinCur; <<
712 sinCur = sinCur*cosStep + cosCur*sinStep <<
713 cosCur = cosCur*cosStep - sinTmp*sinStep <<
714 } <<
715 <<
716 // set sub-envelope and adjust extent <<
717 G4double emin,emax; <<
718 G4BoundingEnvelope benv(polygons); <<
719 if (!benv.CalculateExtent(pAxis,pVoxelLimi <<
720 if (emin < pMin) pMin = emin; <<
721 if (emax > pMax) pMax = emax; <<
722 if (eminlim > pMin && emaxlim < pMax) brea <<
723 } <<
724 // free memory <<
725 for (G4int k=0; k<ksteps+1; ++k) { delete po <<
726 return (pMin < pMax); <<
727 } <<
728 406
>> 407 //
729 // ComputeDimensions 408 // ComputeDimensions
730 // 409 //
731 void G4Polyhedra::ComputeDimensions( G4V << 410 void G4Polyhedra::ComputeDimensions( G4VPVParameterisation* p,
732 const G4i << 411 const G4int n,
733 const G4V << 412 const G4VPhysicalVolume* pRep)
734 { 413 {
735 p->ComputeDimensions(*this,n,pRep); <<
736 } 414 }
737 415
738 // GetEntityType <<
739 // <<
740 G4GeometryType G4Polyhedra::GetEntityType() co <<
741 { <<
742 return {"G4Polyhedra"}; <<
743 } <<
744 416
745 // IsFaceted <<
746 // 417 //
747 G4bool G4Polyhedra::IsFaceted() const << 418 // CreatePolyhedron
748 { <<
749 return true; <<
750 } <<
751 <<
752 // Make a clone of the object <<
753 // 419 //
754 G4VSolid* G4Polyhedra::Clone() const << 420 G4Polyhedron *G4Polyhedra::CreatePolyhedron() const
755 { << 421 {
756 return new G4Polyhedra(*this); << 422 //
757 } << 423 // This has to be fixed in visualization. Fake it for the moment.
>> 424 //
>> 425 if (original_parameters) {
>> 426
>> 427 return new G4PolyhedronPgon( original_parameters->Start_angle,
>> 428 original_parameters->Opening_angle,
>> 429 original_parameters->numSide,
>> 430 original_parameters->Num_z_planes,
>> 431 original_parameters->Z_values,
>> 432 original_parameters->Rmin,
>> 433 original_parameters->Rmax);
>> 434 }
>> 435 else {
>> 436 G4cerr << "G4Polyhedra: visualization of this type of G4Polyhedra is not supported at this time" << G4endl;
>> 437 return 0;
>> 438 }
758 439
759 // Stream object contents to an output stream << 440 }
760 // <<
761 std::ostream& G4Polyhedra::StreamInfo( std::os <<
762 { <<
763 G4long oldprc = os.precision(16); <<
764 os << "------------------------------------- <<
765 << " *** Dump for solid - " << GetName <<
766 << " ================================= <<
767 << " Solid type: G4Polyhedra\n" <<
768 << " Parameters: \n" <<
769 << " starting phi angle : " << startPh <<
770 << " ending phi angle : " << endPhi/ <<
771 << " number of sides : " << numSide <<
772 G4int i=0; <<
773 if (!genericPgon) <<
774 { <<
775 G4int numPlanes = original_parameters->Num <<
776 os << " number of Z planes: " << numPla <<
777 << " Z values: \n"; <<
778 for (i=0; i<numPlanes; ++i) <<
779 { <<
780 os << " Z plane " << i << " <<
781 << original_parameters->Z_values[i] < <<
782 } <<
783 os << " Tangent distances to <<
784 for (i=0; i<numPlanes; ++i) <<
785 { <<
786 os << " Z plane " << i << " <<
787 << original_parameters->Rmin[i] << "\ <<
788 } <<
789 os << " Tangent distances to <<
790 for (i=0; i<numPlanes; ++i) <<
791 { <<
792 os << " Z plane " << i << " <<
793 << original_parameters->Rmax[i] << "\ <<
794 } <<
795 } <<
796 os << " number of RZ points: " << numCorn <<
797 << " RZ values (corners): \n <<
798 for (i=0; i<numCorner; ++i) <<
799 { <<
800 os << " " <<
801 << corners[i].r << ", " << corners[i <<
802 } <<
803 os << "------------------------------------- <<
804 os.precision(oldprc); <<
805 441
806 return os; <<
807 } <<
808 442
809 ////////////////////////////////////////////// <<
810 // 443 //
811 // Return volume << 444 // CreateNURBS
812 <<
813 G4double G4Polyhedra::GetCubicVolume() <<
814 { <<
815 if (fCubicVolume == 0.) <<
816 { <<
817 G4double total = 0.; <<
818 G4int nrz = GetNumRZCorner(); <<
819 G4PolyhedraSideRZ a = GetCorner(nrz - 1); <<
820 for (G4int i=0; i<nrz; ++i) <<
821 { <<
822 G4PolyhedraSideRZ b = GetCorner(i); <<
823 total += (b.r*b.r + b.r*a.r + a.r*a.r)*( <<
824 a = b; <<
825 } <<
826 fCubicVolume = std::abs(total)* <<
827 std::sin((GetEndPhi() - GetStartPhi())/G <<
828 } <<
829 return fCubicVolume; <<
830 } <<
831 <<
832 ////////////////////////////////////////////// <<
833 // 445 //
834 // Return surface area << 446 G4NURBS *G4Polyhedra::CreateNURBS() const
835 <<
836 G4double G4Polyhedra::GetSurfaceArea() <<
837 { 447 {
838 if (fSurfaceArea == 0.) << 448 return 0;
839 { <<
840 G4double total = 0.; <<
841 G4int nrz = GetNumRZCorner(); <<
842 if (IsOpen()) <<
843 { <<
844 G4PolyhedraSideRZ a = GetCorner(nrz - 1) <<
845 for (G4int i=0; i<nrz; ++i) <<
846 { <<
847 G4PolyhedraSideRZ b = GetCorner(i); <<
848 total += a.r*b.z - a.z*b.r; <<
849 a = b; <<
850 } <<
851 total = std::abs(total); <<
852 } <<
853 G4double alp = (GetEndPhi() - GetStartPhi( <<
854 G4double cosa = std::cos(alp); <<
855 G4double sina = std::sin(alp); <<
856 G4PolyhedraSideRZ a = GetCorner(nrz - 1); <<
857 for (G4int i=0; i<nrz; ++i) <<
858 { <<
859 G4PolyhedraSideRZ b = GetCorner(i); <<
860 G4ThreeVector p1(a.r, 0, a.z); <<
861 G4ThreeVector p2(a.r*cosa, a.r*sina, a.z <<
862 G4ThreeVector p3(b.r*cosa, b.r*sina, b.z <<
863 G4ThreeVector p4(b.r, 0, b.z); <<
864 total += GetNumSide()*(G4GeomTools::Quad <<
865 a = b; <<
866 } <<
867 fSurfaceArea = total; <<
868 } <<
869 return fSurfaceArea; <<
870 } 449 }
871 450
872 ////////////////////////////////////////////// <<
873 // <<
874 // Set vector of surface elements, auxiliary m <<
875 // random points on surface <<
876 <<
877 void G4Polyhedra::SetSurfaceElements() const <<
878 { <<
879 fElements = new std::vector<G4Polyhedra::sur <<
880 G4double total = 0.; <<
881 G4int nrz = GetNumRZCorner(); <<
882 <<
883 // set lateral surface elements <<
884 G4double dphi = (GetEndPhi() - GetStartPhi() <<
885 G4double cosa = std::cos(dphi); <<
886 G4double sina = std::sin(dphi); <<
887 G4int ia = nrz - 1; <<
888 for (G4int ib=0; ib<nrz; ++ib) <<
889 { <<
890 G4PolyhedraSideRZ a = GetCorner(ia); <<
891 G4PolyhedraSideRZ b = GetCorner(ib); <<
892 G4Polyhedra::surface_element selem; <<
893 selem.i0 = ia; <<
894 selem.i1 = ib; <<
895 ia = ib; <<
896 if (a.r == 0. && b.r == 0.) continue; <<
897 G4ThreeVector p1(a.r, 0, a.z); <<
898 G4ThreeVector p2(a.r*cosa, a.r*sina, a.z); <<
899 G4ThreeVector p3(b.r*cosa, b.r*sina, b.z); <<
900 G4ThreeVector p4(b.r, 0, b.z); <<
901 if (a.r > 0.) <<
902 { <<
903 selem.i2 = -1; <<
904 total += GetNumSide()*(G4GeomTools::Tria <<
905 selem.area = total; <<
906 fElements->push_back(selem); <<
907 } <<
908 if (b.r > 0.) <<
909 { <<
910 selem.i2 = -2; <<
911 total += GetNumSide()*(G4GeomTools::Tria <<
912 selem.area = total; <<
913 fElements->push_back(selem); <<
914 } <<
915 } <<
916 <<
917 // set elements for phi cuts <<
918 if (IsOpen()) <<
919 { <<
920 G4TwoVectorList contourRZ; <<
921 std::vector<G4int> triangles; <<
922 for (G4int i=0; i<nrz; ++i) <<
923 { <<
924 G4PolyhedraSideRZ corner = GetCorner(i); <<
925 contourRZ.emplace_back(corner.r, corner. <<
926 } <<
927 G4GeomTools::TriangulatePolygon(contourRZ, <<
928 auto ntria = (G4int)triangles.size(); <<
929 for (G4int i=0; i<ntria; i+=3) <<
930 { <<
931 G4Polyhedra::surface_element selem; <<
932 selem.i0 = triangles[i]; <<
933 selem.i1 = triangles[i+1]; <<
934 selem.i2 = triangles[i+2]; <<
935 G4PolyhedraSideRZ a = GetCorner(selem.i0 <<
936 G4PolyhedraSideRZ b = GetCorner(selem.i1 <<
937 G4PolyhedraSideRZ c = GetCorner(selem.i2 <<
938 G4double stria = <<
939 std::abs(G4GeomTools::TriangleArea(a.r <<
940 total += stria; <<
941 selem.area = total; <<
942 fElements->push_back(selem); // start ph <<
943 total += stria; <<
944 selem.area = total; <<
945 selem.i0 += nrz; <<
946 fElements->push_back(selem); // end phi <<
947 } <<
948 } <<
949 } <<
950 451
951 ////////////////////////////////////////////// <<
952 // <<
953 // Generate random point on surface <<
954 452
955 G4ThreeVector G4Polyhedra::GetPointOnSurface() <<
956 { <<
957 // Set surface elements <<
958 if (fElements == nullptr) <<
959 { <<
960 G4AutoLock l(&surface_elementsMutex); <<
961 SetSurfaceElements(); <<
962 l.unlock(); <<
963 } <<
964 <<
965 // Select surface element <<
966 G4Polyhedra::surface_element selem; <<
967 selem = fElements->back(); <<
968 G4double select = selem.area*G4QuickRand(); <<
969 auto it = std::lower_bound(fElements->begin( <<
970 [](const G4Polyhe <<
971 -> G4bool { retur <<
972 <<
973 // Generate random point <<
974 G4double x = 0, y = 0, z = 0; <<
975 G4double u = G4QuickRand(); <<
976 G4double v = G4QuickRand(); <<
977 if (u + v > 1.) { u = 1. - u; v = 1. - v; } <<
978 G4int i0 = (*it).i0; <<
979 G4int i1 = (*it).i1; <<
980 G4int i2 = (*it).i2; <<
981 if (i2 < 0) // lateral surface <<
982 { <<
983 // sample point <<
984 G4int nside = GetNumSide(); <<
985 G4double dphi = (GetEndPhi() - GetStartPhi <<
986 G4double cosa = std::cos(dphi); <<
987 G4double sina = std::sin(dphi); <<
988 G4PolyhedraSideRZ a = GetCorner(i0); <<
989 G4PolyhedraSideRZ b = GetCorner(i1); <<
990 G4ThreeVector p0(a.r, 0, a.z); <<
991 G4ThreeVector p1(b.r, 0, b.z); <<
992 G4ThreeVector p2(b.r*cosa, b.r*sina, b.z); <<
993 if (i2 == -1) p1.set(a.r*cosa, a.r*sina, a <<
994 p0 += (p1 - p0)*u + (p2 - p0)*v; <<
995 // find selected side and rotate point <<
996 G4double scurr = (*it).area; <<
997 G4double sprev = (it == fElements->begin() <<
998 G4int iside = nside*(select - sprev)/(scur <<
999 if (iside == 0 && GetStartPhi() == 0.) <<
1000 { <<
1001 x = p0.x(); <<
1002 y = p0.y(); <<
1003 z = p0.z(); <<
1004 } <<
1005 else <<
1006 { <<
1007 if (iside == nside) --iside; // iside m <<
1008 G4double phi = iside*dphi + GetStartPhi <<
1009 G4double cosphi = std::cos(phi); <<
1010 G4double sinphi = std::sin(phi); <<
1011 x = p0.x()*cosphi - p0.y()*sinphi; <<
1012 y = p0.x()*sinphi + p0.y()*cosphi; <<
1013 z = p0.z(); <<
1014 } <<
1015 } <<
1016 else // phi cut <<
1017 { <<
1018 G4int nrz = GetNumRZCorner(); <<
1019 G4double phi = (i0 < nrz) ? GetStartPhi() <<
1020 if (i0 >= nrz) { i0 -= nrz; } <<
1021 G4PolyhedraSideRZ p0 = GetCorner(i0); <<
1022 G4PolyhedraSideRZ p1 = GetCorner(i1); <<
1023 G4PolyhedraSideRZ p2 = GetCorner(i2); <<
1024 G4double r = (p1.r - p0.r)*u + (p2.r - p0 <<
1025 x = r*std::cos(phi); <<
1026 y = r*std::sin(phi); <<
1027 z = (p1.z - p0.z)*u + (p2.z - p0.z)*v + p <<
1028 } <<
1029 return {x, y, z}; <<
1030 } <<
1031 453
1032 ///////////////////////////////////////////// <<
1033 // 454 //
1034 // CreatePolyhedron << 455 // G4Polyhedra::G4PolyhedraHistorical stuff
1035 << 456 //
1036 G4Polyhedron* G4Polyhedra::CreatePolyhedron() << 457 G4Polyhedra::G4PolyhedraHistorical::~G4PolyhedraHistorical()
1037 { 458 {
1038 std::vector<G4TwoVector> rz(numCorner); << 459 delete [] Z_values;
1039 for (G4int i = 0; i < numCorner; ++i) << 460 delete [] Rmin;
1040 rz[i].set(corners[i].r, corners[i].z); << 461 delete [] Rmax;
1041 return new G4PolyhedronPgon(startPhi, endPh <<
1042 } 462 }
1043 463
1044 // SetOriginalParameters << 464 G4Polyhedra::G4PolyhedraHistorical::G4PolyhedraHistorical( const G4PolyhedraHistorical &source )
1045 // <<
1046 void G4Polyhedra::SetOriginalParameters(G4Red <<
1047 { 465 {
1048 G4int numPlanes = numCorner; << 466 Start_angle = source.Start_angle;
1049 G4bool isConvertible = true; << 467 Opening_angle = source.Opening_angle;
1050 G4double Zmax=rz->Bmax(); << 468 numSide = source.numSide;
1051 rz->StartWithZMin(); << 469 Num_z_planes = source.Num_z_planes;
1052 << 470
1053 // Prepare vectors for storage << 471 Z_values = new G4double[Num_z_planes];
1054 // << 472 Rmin = new G4double[Num_z_planes];
1055 std::vector<G4double> Z; << 473 Rmax = new G4double[Num_z_planes];
1056 std::vector<G4double> Rmin; << 474
1057 std::vector<G4double> Rmax; << 475 G4int i;
1058 << 476 for( i = 0; i < Num_z_planes; i++) {
1059 G4int countPlanes=1; << 477 Z_values[i] = source.Z_values[i];
1060 G4int icurr=0; << 478 Rmin[i] = source.Rmin[i];
1061 G4int icurl=0; << 479 Rmax[i] = source.Rmax[i];
1062 << 480 }
1063 // first plane Z=Z[0] <<
1064 // <<
1065 Z.push_back(corners[0].z); <<
1066 G4double Zprev=Z[0]; <<
1067 if (Zprev == corners[1].z) <<
1068 { <<
1069 Rmin.push_back(corners[0].r); <<
1070 Rmax.push_back (corners[1].r);icurr=1; <<
1071 } <<
1072 else if (Zprev == corners[numPlanes-1].z) <<
1073 { <<
1074 Rmin.push_back(corners[numPlanes-1].r); <<
1075 Rmax.push_back (corners[0].r); <<
1076 icurl=numPlanes-1; <<
1077 } <<
1078 else <<
1079 { <<
1080 Rmin.push_back(corners[0].r); <<
1081 Rmax.push_back (corners[0].r); <<
1082 } <<
1083 <<
1084 // next planes until last <<
1085 // <<
1086 G4int inextr=0, inextl=0; <<
1087 for (G4int i=0; i < numPlanes-2; ++i) <<
1088 { <<
1089 inextr=1+icurr; <<
1090 inextl=(icurl <= 0)? numPlanes-1 : icurl- <<
1091 <<
1092 if((corners[inextr].z >= Zmax) & (corners <<
1093 <<
1094 G4double Zleft = corners[inextl].z; <<
1095 G4double Zright = corners[inextr].z; <<
1096 if(Zright>Zleft) <<
1097 { <<
1098 Z.push_back(Zleft); <<
1099 countPlanes++; <<
1100 G4double difZr=corners[inextr].z - corn <<
1101 G4double difZl=corners[inextl].z - corn <<
1102 <<
1103 if(std::fabs(difZl) < kCarTolerance) <<
1104 { <<
1105 if(std::fabs(difZr) < kCarTolerance) <<
1106 { <<
1107 Rmin.push_back(corners[inextl].r); <<
1108 Rmax.push_back(corners[icurr].r); <<
1109 } <<
1110 else <<
1111 { <<
1112 Rmin.push_back(corners[inextl].r); <<
1113 Rmax.push_back(corners[icurr].r + ( <<
1114 *(corners[ine <<
1115 } <<
1116 } <<
1117 else if (difZl >= kCarTolerance) <<
1118 { <<
1119 if(std::fabs(difZr) < kCarTolerance) <<
1120 { <<
1121 Rmin.push_back(corners[icurl].r); <<
1122 Rmax.push_back(corners[icurr].r); <<
1123 } <<
1124 else <<
1125 { <<
1126 Rmin.push_back(corners[icurl].r); <<
1127 Rmax.push_back(corners[icurr].r + ( <<
1128 *(corners[ine <<
1129 } <<
1130 } <<
1131 else <<
1132 { <<
1133 isConvertible=false; break; <<
1134 } <<
1135 icurl=(icurl == 0)? numPlanes-1 : icurl <<
1136 } <<
1137 else if(std::fabs(Zright-Zleft)<kCarToler <<
1138 { <<
1139 Z.push_back(Zleft); <<
1140 ++countPlanes; <<
1141 ++icurr; <<
1142 <<
1143 icurl=(icurl == 0)? numPlanes-1 : icurl <<
1144 <<
1145 Rmin.push_back(corners[inextl].r); <<
1146 Rmax.push_back (corners[inextr].r); <<
1147 } <<
1148 else // Zright<Zleft <<
1149 { <<
1150 Z.push_back(Zright); <<
1151 ++countPlanes; <<
1152 <<
1153 G4double difZr=corners[inextr].z - corn <<
1154 G4double difZl=corners[inextl].z - corn <<
1155 if(std::fabs(difZr) < kCarTolerance) <<
1156 { <<
1157 if(std::fabs(difZl) < kCarTolerance) <<
1158 { <<
1159 Rmax.push_back(corners[inextr].r); <<
1160 Rmin.push_back(corners[icurr].r); <<
1161 } <<
1162 else <<
1163 { <<
1164 Rmin.push_back(corners[icurl].r + ( <<
1165 * (corners[in <<
1166 Rmax.push_back(corners[inextr].r); <<
1167 } <<
1168 ++icurr; <<
1169 } // plate <<
1170 else if (difZr >= kCarTolerance) <<
1171 { <<
1172 if(std::fabs(difZl) < kCarTolerance) <<
1173 { <<
1174 Rmax.push_back(corners[inextr].r); <<
1175 Rmin.push_back (corners[icurr].r); <<
1176 } <<
1177 else <<
1178 { <<
1179 Rmax.push_back(corners[inextr].r); <<
1180 Rmin.push_back (corners[icurl].r+(Z <<
1181 * (corners[ <<
1182 } <<
1183 ++icurr; <<
1184 } <<
1185 else <<
1186 { <<
1187 isConvertible=false; break; <<
1188 } <<
1189 } <<
1190 } // end for loop <<
1191 <<
1192 // last plane Z=Zmax <<
1193 // <<
1194 Z.push_back(Zmax); <<
1195 ++countPlanes; <<
1196 inextr=1+icurr; <<
1197 inextl=(icurl <= 0)? numPlanes-1 : icurl-1; <<
1198 <<
1199 Rmax.push_back(corners[inextr].r); <<
1200 Rmin.push_back(corners[inextl].r); <<
1201 <<
1202 // Set original parameters Rmin,Rmax,Z <<
1203 // <<
1204 if(isConvertible) <<
1205 { <<
1206 original_parameters = new G4PolyhedraHisto <<
1207 original_parameters->numSide = numSide; <<
1208 original_parameters->Z_values = new G4doub <<
1209 original_parameters->Rmin = new G4double[c <<
1210 original_parameters->Rmax = new G4double[c <<
1211 <<
1212 for(G4int j=0; j < countPlanes; ++j) <<
1213 { <<
1214 original_parameters->Z_values[j] = Z[j]; <<
1215 original_parameters->Rmax[j] = Rmax[j]; <<
1216 original_parameters->Rmin[j] = Rmin[j]; <<
1217 } <<
1218 original_parameters->Start_angle = startPh <<
1219 original_parameters->Opening_angle = endPh <<
1220 original_parameters->Num_z_planes = countP <<
1221 <<
1222 } <<
1223 else // Set parameters(r,z) with Rmin==0 a <<
1224 { <<
1225 #ifdef G4SPECSDEBUG <<
1226 std::ostringstream message; <<
1227 message << "Polyhedra " << GetName() << G <<
1228 << "cannot be converted to Polyhedra wi <<
1229 G4Exception("G4Polyhedra::SetOriginalPara <<
1230 "GeomSolids0002", JustWarning <<
1231 #endif <<
1232 original_parameters = new G4PolyhedraHist <<
1233 original_parameters->numSide = numSide; <<
1234 original_parameters->Z_values = new G4dou <<
1235 original_parameters->Rmin = new G4double[ <<
1236 original_parameters->Rmax = new G4double[ <<
1237 <<
1238 for(G4int j=0; j < numPlanes; ++j) <<
1239 { <<
1240 original_parameters->Z_values[j] = corn <<
1241 original_parameters->Rmax[j] = corners[ <<
1242 original_parameters->Rmin[j] = 0.0; <<
1243 } <<
1244 original_parameters->Start_angle = startP <<
1245 original_parameters->Opening_angle = endP <<
1246 original_parameters->Num_z_planes = numPl <<
1247 } <<
1248 } 481 }
1249 482
1250 #endif <<
1251 483