| 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 // Implementation of G4Polyhedra, a CSG polyhe << 26 //
27 // as an inherited class of G4VCSGfaceted. << 27 // $Id: G4Polyhedra.cc 88514 2015-02-25 09:59:10Z gcosmo $
>> 28 //
>> 29 //
>> 30 // --------------------------------------------------------------------
>> 31 // GEANT 4 class source file
>> 32 //
>> 33 //
>> 34 // G4Polyhedra.cc
>> 35 //
>> 36 // Implementation of a CSG polyhedra, as an inherited class of G4VCSGfaceted.
>> 37 //
>> 38 // To be done:
>> 39 // * Cracks: there are probably small cracks in the seams between the
>> 40 // phi face (G4PolyPhiFace) and sides (G4PolyhedraSide) that are not
>> 41 // entirely leakproof. Also, I am not sure all vertices are leak proof.
>> 42 // * Many optimizations are possible, but not implemented.
>> 43 // * Visualization needs to be updated outside of this routine.
28 // 44 //
29 // Utility classes: 45 // Utility classes:
30 // * G4EnclosingCylinder: decided a quick c << 46 // * G4EnclosingCylinder: I decided a quick check of geometry would be a
31 // good idea (for CPU speed). If the quic 47 // good idea (for CPU speed). If the quick check fails, the regular
32 // full-blown G4VCSGfaceted version is in 48 // full-blown G4VCSGfaceted version is invoked.
33 // * G4ReduciblePolygon: Really meant as a 49 // * G4ReduciblePolygon: Really meant as a check of input parameters,
34 // this utility class also "converts" the 50 // this utility class also "converts" the GEANT3-like PGON/PCON
35 // arguments into the newer ones. 51 // arguments into the newer ones.
36 // Both these classes are implemented outside 52 // Both these classes are implemented outside this file because they are
37 // shared with G4Polycone. 53 // shared with G4Polycone.
38 // 54 //
39 // Author: David C. Williams (davidw@scipp.ucs <<
40 // ------------------------------------------- 55 // --------------------------------------------------------------------
41 56
42 #include "G4Polyhedra.hh" 57 #include "G4Polyhedra.hh"
43 58
44 #if !defined(G4GEOM_USE_UPOLYHEDRA) 59 #if !defined(G4GEOM_USE_UPOLYHEDRA)
45 60
46 #include "G4PolyhedraSide.hh" 61 #include "G4PolyhedraSide.hh"
47 #include "G4PolyPhiFace.hh" 62 #include "G4PolyPhiFace.hh"
48 63
49 #include "G4GeomTools.hh" << 64 #include "Randomize.hh"
50 #include "G4VoxelLimits.hh" <<
51 #include "G4AffineTransform.hh" <<
52 #include "G4BoundingEnvelope.hh" <<
53 <<
54 #include "G4QuickRand.hh" <<
55 65
56 #include "G4EnclosingCylinder.hh" 66 #include "G4EnclosingCylinder.hh"
57 #include "G4ReduciblePolygon.hh" 67 #include "G4ReduciblePolygon.hh"
58 #include "G4VPVParameterisation.hh" 68 #include "G4VPVParameterisation.hh"
59 69
60 namespace << 70 #include <sstream>
61 { <<
62 G4Mutex surface_elementsMutex = G4MUTEX_INIT <<
63 } <<
64 71
65 using namespace CLHEP; 72 using namespace CLHEP;
66 73
>> 74 //
67 // Constructor (GEANT3 style parameters) 75 // Constructor (GEANT3 style parameters)
68 // 76 //
69 // GEANT3 PGON radii are specified in the dist 77 // GEANT3 PGON radii are specified in the distance to the norm of each face.
70 // << 78 //
71 G4Polyhedra::G4Polyhedra( const G4String& name << 79 G4Polyhedra::G4Polyhedra( const G4String& name,
72 G4double phiSt 80 G4double phiStart,
73 G4double thePh 81 G4double thePhiTotal,
74 G4int theNumSi << 82 G4int theNumSide,
75 G4int numZPlan 83 G4int numZPlanes,
76 const G4double zPlan 84 const G4double zPlane[],
77 const G4double rInne 85 const G4double rInner[],
78 const G4double rOute 86 const G4double rOuter[] )
79 : G4VCSGfaceted( name ) << 87 : G4VCSGfaceted( name ), genericPgon(false)
80 { 88 {
81 if (theNumSide <= 0) 89 if (theNumSide <= 0)
82 { 90 {
83 std::ostringstream message; 91 std::ostringstream message;
84 message << "Solid must have at least one s 92 message << "Solid must have at least one side - " << GetName() << G4endl
85 << " No sides specified !"; 93 << " No sides specified !";
86 G4Exception("G4Polyhedra::G4Polyhedra()", 94 G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
87 FatalErrorInArgument, message) 95 FatalErrorInArgument, message);
88 } 96 }
89 97
90 // 98 //
91 // Calculate conversion factor from G3 radiu 99 // Calculate conversion factor from G3 radius to G4 radius
92 // 100 //
93 G4double phiTotal = thePhiTotal; 101 G4double phiTotal = thePhiTotal;
94 if ( (phiTotal <=0) || (phiTotal >= twopi*(1 102 if ( (phiTotal <=0) || (phiTotal >= twopi*(1-DBL_EPSILON)) )
95 { phiTotal = twopi; } 103 { phiTotal = twopi; }
96 G4double convertRad = std::cos(0.5*phiTotal/ 104 G4double convertRad = std::cos(0.5*phiTotal/theNumSide);
97 105
98 // 106 //
99 // Some historical stuff 107 // Some historical stuff
100 // 108 //
101 original_parameters = new G4PolyhedraHistori 109 original_parameters = new G4PolyhedraHistorical;
102 << 110
103 original_parameters->numSide = theNumSide; 111 original_parameters->numSide = theNumSide;
104 original_parameters->Start_angle = phiStart; 112 original_parameters->Start_angle = phiStart;
105 original_parameters->Opening_angle = phiTota 113 original_parameters->Opening_angle = phiTotal;
106 original_parameters->Num_z_planes = numZPlan 114 original_parameters->Num_z_planes = numZPlanes;
107 original_parameters->Z_values = new G4double 115 original_parameters->Z_values = new G4double[numZPlanes];
108 original_parameters->Rmin = new G4double[num 116 original_parameters->Rmin = new G4double[numZPlanes];
109 original_parameters->Rmax = new G4double[num 117 original_parameters->Rmax = new G4double[numZPlanes];
110 118
111 for (G4int i=0; i<numZPlanes; ++i) << 119 G4int i;
>> 120 for (i=0; i<numZPlanes; i++)
112 { 121 {
113 if (( i < numZPlanes-1) && ( zPlane[i] == 122 if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
114 { 123 {
115 if( (rInner[i] > rOuter[i+1]) 124 if( (rInner[i] > rOuter[i+1])
116 ||(rInner[i+1] > rOuter[i]) ) 125 ||(rInner[i+1] > rOuter[i]) )
117 { 126 {
118 DumpInfo(); 127 DumpInfo();
119 std::ostringstream message; 128 std::ostringstream message;
120 message << "Cannot create a Polyhedra 129 message << "Cannot create a Polyhedra with no contiguous segments."
121 << G4endl 130 << G4endl
122 << " Segments are not c 131 << " Segments are not contiguous !" << G4endl
123 << " rMin[" << i << "] 132 << " rMin[" << i << "] = " << rInner[i]
124 << " -- rMax[" << i+1 << "] = 133 << " -- rMax[" << i+1 << "] = " << rOuter[i+1] << G4endl
125 << " rMin[" << i+1 << " 134 << " rMin[" << i+1 << "] = " << rInner[i+1]
126 << " -- rMax[" << i << "] = " 135 << " -- rMax[" << i << "] = " << rOuter[i];
127 G4Exception("G4Polyhedra::G4Polyhedra( 136 G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
128 FatalErrorInArgument, mess 137 FatalErrorInArgument, message);
129 } 138 }
130 } 139 }
131 original_parameters->Z_values[i] = zPlane[ 140 original_parameters->Z_values[i] = zPlane[i];
132 original_parameters->Rmin[i] = rInner[i]/c 141 original_parameters->Rmin[i] = rInner[i]/convertRad;
133 original_parameters->Rmax[i] = rOuter[i]/c 142 original_parameters->Rmax[i] = rOuter[i]/convertRad;
134 } 143 }
135 << 144
136 << 145
137 // 146 //
138 // Build RZ polygon using special PCON/PGON 147 // Build RZ polygon using special PCON/PGON GEANT3 constructor
139 // 148 //
140 auto rz = new G4ReduciblePolygon( rInner, rO << 149 G4ReduciblePolygon *rz =
>> 150 new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
141 rz->ScaleA( 1/convertRad ); 151 rz->ScaleA( 1/convertRad );
142 << 152
143 // 153 //
144 // Do the real work 154 // Do the real work
145 // 155 //
146 Create( phiStart, phiTotal, theNumSide, rz ) 156 Create( phiStart, phiTotal, theNumSide, rz );
147 << 157
148 delete rz; 158 delete rz;
149 } 159 }
150 160
>> 161
>> 162 //
151 // Constructor (generic parameters) 163 // Constructor (generic parameters)
152 // 164 //
153 G4Polyhedra::G4Polyhedra( const G4String& name << 165 G4Polyhedra::G4Polyhedra( const G4String& name,
154 G4double phiSt 166 G4double phiStart,
155 G4double phiTo 167 G4double phiTotal,
156 G4int theNu << 168 G4int theNumSide,
157 G4int numRZ 169 G4int numRZ,
158 const G4double r[], 170 const G4double r[],
159 const G4double z[] 171 const G4double z[] )
160 : G4VCSGfaceted( name ), genericPgon(true) 172 : G4VCSGfaceted( name ), genericPgon(true)
161 { << 173 {
162 if (theNumSide <= 0) 174 if (theNumSide <= 0)
163 { 175 {
164 std::ostringstream message; 176 std::ostringstream message;
165 message << "Solid must have at least one s 177 message << "Solid must have at least one side - " << GetName() << G4endl
166 << " No sides specified !"; 178 << " No sides specified !";
167 G4Exception("G4Polyhedra::G4Polyhedra()", 179 G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
168 FatalErrorInArgument, message) 180 FatalErrorInArgument, message);
169 } 181 }
170 182
171 auto rz = new G4ReduciblePolygon( r, z, numR << 183 G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
172 << 184
173 Create( phiStart, phiTotal, theNumSide, rz ) 185 Create( phiStart, phiTotal, theNumSide, rz );
174 << 186
175 // Set original_parameters struct for consis 187 // Set original_parameters struct for consistency
176 // 188 //
177 SetOriginalParameters(rz); 189 SetOriginalParameters(rz);
178 << 190
179 delete rz; 191 delete rz;
180 } 192 }
181 193
>> 194
>> 195 //
182 // Create 196 // Create
183 // 197 //
184 // Generic create routine, called by each cons 198 // Generic create routine, called by each constructor
185 // after conversion of arguments 199 // after conversion of arguments
186 // 200 //
187 void G4Polyhedra::Create( G4double phiStart, 201 void G4Polyhedra::Create( G4double phiStart,
188 G4double phiTotal, 202 G4double phiTotal,
189 G4int theNumSide, << 203 G4int theNumSide,
190 G4ReduciblePolygon* << 204 G4ReduciblePolygon *rz )
191 { 205 {
192 // 206 //
193 // Perform checks of rz values 207 // Perform checks of rz values
194 // 208 //
195 if (rz->Amin() < 0.0) 209 if (rz->Amin() < 0.0)
196 { 210 {
197 std::ostringstream message; 211 std::ostringstream message;
198 message << "Illegal input parameters - " < 212 message << "Illegal input parameters - " << GetName() << G4endl
199 << " All R values must be > 213 << " All R values must be >= 0 !";
200 G4Exception("G4Polyhedra::Create()", "Geom 214 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
201 FatalErrorInArgument, message) 215 FatalErrorInArgument, message);
202 } 216 }
203 217
204 G4double rzArea = rz->Area(); 218 G4double rzArea = rz->Area();
205 if (rzArea < -kCarTolerance) 219 if (rzArea < -kCarTolerance)
206 { <<
207 rz->ReverseOrder(); 220 rz->ReverseOrder();
208 } << 221
209 else if (rzArea < kCarTolerance) << 222 else if (rzArea < -kCarTolerance)
210 { 223 {
211 std::ostringstream message; 224 std::ostringstream message;
212 message << "Illegal input parameters - " < 225 message << "Illegal input parameters - " << GetName() << G4endl
213 << " R/Z cross section is z 226 << " R/Z cross section is zero or near zero: " << rzArea;
214 G4Exception("G4Polyhedra::Create()", "Geom 227 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
215 FatalErrorInArgument, message) 228 FatalErrorInArgument, message);
216 } 229 }
217 << 230
218 if ( (!rz->RemoveDuplicateVertices( kCarTole 231 if ( (!rz->RemoveDuplicateVertices( kCarTolerance ))
219 || (!rz->RemoveRedundantVertices( kCarTole << 232 || (!rz->RemoveRedundantVertices( kCarTolerance )) )
220 { 233 {
221 std::ostringstream message; 234 std::ostringstream message;
222 message << "Illegal input parameters - " < 235 message << "Illegal input parameters - " << GetName() << G4endl
223 << " Too few unique R/Z val 236 << " Too few unique R/Z values !";
224 G4Exception("G4Polyhedra::Create()", "Geom 237 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
225 FatalErrorInArgument, message) 238 FatalErrorInArgument, message);
226 } 239 }
227 240
228 if (rz->CrossesItself( 1/kInfinity )) << 241 if (rz->CrossesItself( 1/kInfinity ))
229 { 242 {
230 std::ostringstream message; 243 std::ostringstream message;
231 message << "Illegal input parameters - " < 244 message << "Illegal input parameters - " << GetName() << G4endl
232 << " R/Z segments cross !"; 245 << " R/Z segments cross !";
233 G4Exception("G4Polyhedra::Create()", "Geom 246 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
234 FatalErrorInArgument, message) 247 FatalErrorInArgument, message);
235 } 248 }
236 249
237 numCorner = rz->NumVertices(); 250 numCorner = rz->NumVertices();
238 251
239 252
240 startPhi = phiStart; 253 startPhi = phiStart;
241 while( startPhi < 0 ) // Loop checking, 1 << 254 while( startPhi < 0 ) startPhi += twopi;
242 startPhi += twopi; <<
243 // 255 //
244 // Phi opening? Account for some possible ro 256 // Phi opening? Account for some possible roundoff, and interpret
245 // nonsense value as representing no phi ope 257 // nonsense value as representing no phi opening
246 // 258 //
247 if ( (phiTotal <= 0) || (phiTotal > twopi*(1 259 if ( (phiTotal <= 0) || (phiTotal > twopi*(1-DBL_EPSILON)) )
248 { 260 {
249 phiIsOpen = false; 261 phiIsOpen = false;
250 endPhi = startPhi + twopi; << 262 endPhi = phiStart+twopi;
251 } 263 }
252 else 264 else
253 { 265 {
254 phiIsOpen = true; 266 phiIsOpen = true;
255 endPhi = startPhi + phiTotal; << 267
>> 268 //
>> 269 // Convert phi into our convention
>> 270 //
>> 271 endPhi = phiStart+phiTotal;
>> 272 while( endPhi < startPhi ) endPhi += twopi;
256 } 273 }
257 << 274
258 // 275 //
259 // Save number sides 276 // Save number sides
260 // 277 //
261 numSide = theNumSide; 278 numSide = theNumSide;
262 << 279
263 // 280 //
264 // Allocate corner array. << 281 // Allocate corner array.
265 // 282 //
266 corners = new G4PolyhedraSideRZ[numCorner]; 283 corners = new G4PolyhedraSideRZ[numCorner];
267 284
268 // 285 //
269 // Copy corners 286 // Copy corners
270 // 287 //
271 G4ReduciblePolygonIterator iterRZ(rz); 288 G4ReduciblePolygonIterator iterRZ(rz);
272 << 289
273 G4PolyhedraSideRZ *next = corners; 290 G4PolyhedraSideRZ *next = corners;
274 iterRZ.Begin(); 291 iterRZ.Begin();
275 do // Loop checking, 13.08.2015, G.Cosmo << 292 do
276 { 293 {
277 next->r = iterRZ.GetA(); 294 next->r = iterRZ.GetA();
278 next->z = iterRZ.GetB(); 295 next->z = iterRZ.GetB();
279 } while( ++next, iterRZ.Next() ); 296 } while( ++next, iterRZ.Next() );
280 << 297
281 // 298 //
282 // Allocate face pointer array 299 // Allocate face pointer array
283 // 300 //
284 numFace = phiIsOpen ? numCorner+2 : numCorne 301 numFace = phiIsOpen ? numCorner+2 : numCorner;
285 faces = new G4VCSGface*[numFace]; 302 faces = new G4VCSGface*[numFace];
286 << 303
287 // 304 //
288 // Construct side faces 305 // Construct side faces
289 // 306 //
290 // To do so properly, we need to keep track 307 // To do so properly, we need to keep track of four successive RZ
291 // corners. 308 // corners.
292 // 309 //
293 // But! Don't construct a face if both point 310 // But! Don't construct a face if both points are at zero radius!
294 // 311 //
295 G4PolyhedraSideRZ* corner = corners, << 312 G4PolyhedraSideRZ *corner = corners,
296 * prev = corners + numCorne << 313 *prev = corners + numCorner-1,
297 * nextNext; << 314 *nextNext;
298 G4VCSGface** face = faces; << 315 G4VCSGface **face = faces;
299 do // Loop checking, 13.08.2015, G.Cosmo << 316 do
300 { 317 {
301 next = corner+1; 318 next = corner+1;
302 if (next >= corners+numCorner) next = corn 319 if (next >= corners+numCorner) next = corners;
303 nextNext = next+1; 320 nextNext = next+1;
304 if (nextNext >= corners+numCorner) nextNex 321 if (nextNext >= corners+numCorner) nextNext = corners;
305 << 322
306 if (corner->r < 1/kInfinity && next->r < 1 323 if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
307 /* 324 /*
308 // We must decide here if we can dare decl 325 // We must decide here if we can dare declare one of our faces
309 // as having a "valid" normal (i.e. allBeh 326 // as having a "valid" normal (i.e. allBehind = true). This
310 // is never possible if the face faces "in 327 // is never possible if the face faces "inward" in r *unless*
311 // we have only one side 328 // we have only one side
312 // 329 //
313 G4bool allBehind; 330 G4bool allBehind;
314 if ((corner->z > next->z) && (numSide > 1) 331 if ((corner->z > next->z) && (numSide > 1))
315 { 332 {
316 allBehind = false; 333 allBehind = false;
317 } 334 }
318 else 335 else
319 { 336 {
320 // 337 //
321 // Otherwise, it is only true if the lin 338 // Otherwise, it is only true if the line passing
322 // through the two points of the segment 339 // through the two points of the segment do not
323 // split the r/z cross section 340 // split the r/z cross section
324 // 341 //
325 allBehind = !rz->BisectedBy( corner->r, 342 allBehind = !rz->BisectedBy( corner->r, corner->z,
326 next->r, ne 343 next->r, next->z, kCarTolerance );
327 } 344 }
328 */ 345 */
329 *face++ = new G4PolyhedraSide( prev, corne 346 *face++ = new G4PolyhedraSide( prev, corner, next, nextNext,
330 numSide, startPhi, endPhi-sta 347 numSide, startPhi, endPhi-startPhi, phiIsOpen );
331 } while( prev=corner, corner=next, corner > 348 } while( prev=corner, corner=next, corner > corners );
332 << 349
333 if (phiIsOpen) 350 if (phiIsOpen)
334 { 351 {
335 // 352 //
336 // Construct phi open edges 353 // Construct phi open edges
337 // 354 //
338 *face++ = new G4PolyPhiFace( rz, startPhi, 355 *face++ = new G4PolyPhiFace( rz, startPhi, phiTotal/numSide, endPhi );
339 *face++ = new G4PolyPhiFace( rz, endPhi, 356 *face++ = new G4PolyPhiFace( rz, endPhi, phiTotal/numSide, startPhi );
340 } 357 }
341 << 358
342 // 359 //
343 // We might have dropped a face or two: reca 360 // We might have dropped a face or two: recalculate numFace
344 // 361 //
345 numFace = (G4int)(face-faces); << 362 numFace = face-faces;
346 << 363
347 // 364 //
348 // Make enclosingCylinder 365 // Make enclosingCylinder
349 // 366 //
350 enclosingCylinder = 367 enclosingCylinder =
351 new G4EnclosingCylinder( rz, phiIsOpen, ph 368 new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
352 } 369 }
353 370
>> 371
>> 372 //
354 // Fake default constructor - sets only member 373 // Fake default constructor - sets only member data and allocates memory
355 // for usage restri 374 // for usage restricted to object persistency.
356 // 375 //
357 G4Polyhedra::G4Polyhedra( __void__& a ) 376 G4Polyhedra::G4Polyhedra( __void__& a )
358 : G4VCSGfaceted(a), startPhi(0.), endPhi(0.) << 377 : G4VCSGfaceted(a), numSide(0), startPhi(0.), endPhi(0.),
>> 378 phiIsOpen(false), genericPgon(false), numCorner(0), corners(0),
>> 379 original_parameters(0), enclosingCylinder(0)
359 { 380 {
360 } 381 }
361 382
>> 383
>> 384 //
362 // Destructor 385 // Destructor
363 // 386 //
364 G4Polyhedra::~G4Polyhedra() 387 G4Polyhedra::~G4Polyhedra()
365 { 388 {
366 delete [] corners; 389 delete [] corners;
367 delete original_parameters; << 390 if (original_parameters) delete original_parameters;
>> 391
368 delete enclosingCylinder; 392 delete enclosingCylinder;
369 delete fElements; <<
370 delete fpPolyhedron; <<
371 corners = nullptr; <<
372 original_parameters = nullptr; <<
373 enclosingCylinder = nullptr; <<
374 fElements = nullptr; <<
375 fpPolyhedron = nullptr; <<
376 } 393 }
377 394
>> 395
>> 396 //
378 // Copy constructor 397 // Copy constructor
379 // 398 //
380 G4Polyhedra::G4Polyhedra( const G4Polyhedra& s << 399 G4Polyhedra::G4Polyhedra( const G4Polyhedra &source )
381 : G4VCSGfaceted( source ) 400 : G4VCSGfaceted( source )
382 { 401 {
383 CopyStuff( source ); 402 CopyStuff( source );
384 } 403 }
385 404
>> 405
>> 406 //
386 // Assignment operator 407 // Assignment operator
387 // 408 //
388 G4Polyhedra &G4Polyhedra::operator=( const G4P << 409 G4Polyhedra &G4Polyhedra::operator=( const G4Polyhedra &source )
389 { 410 {
390 if (this == &source) return *this; 411 if (this == &source) return *this;
391 412
392 G4VCSGfaceted::operator=( source ); 413 G4VCSGfaceted::operator=( source );
393 << 414
394 delete [] corners; 415 delete [] corners;
395 delete original_parameters; << 416 if (original_parameters) delete original_parameters;
>> 417
396 delete enclosingCylinder; 418 delete enclosingCylinder;
397 << 419
398 CopyStuff( source ); 420 CopyStuff( source );
399 << 421
400 return *this; 422 return *this;
401 } 423 }
402 424
>> 425
>> 426 //
403 // CopyStuff 427 // CopyStuff
404 // 428 //
405 void G4Polyhedra::CopyStuff( const G4Polyhedra << 429 void G4Polyhedra::CopyStuff( const G4Polyhedra &source )
406 { 430 {
407 // 431 //
408 // Simple stuff 432 // Simple stuff
409 // 433 //
410 numSide = source.numSide; 434 numSide = source.numSide;
411 startPhi = source.startPhi; 435 startPhi = source.startPhi;
412 endPhi = source.endPhi; 436 endPhi = source.endPhi;
413 phiIsOpen = source.phiIsOpen; 437 phiIsOpen = source.phiIsOpen;
414 numCorner = source.numCorner; 438 numCorner = source.numCorner;
415 genericPgon= source.genericPgon; 439 genericPgon= source.genericPgon;
416 440
417 // 441 //
418 // The corner array 442 // The corner array
419 // 443 //
420 corners = new G4PolyhedraSideRZ[numCorner]; 444 corners = new G4PolyhedraSideRZ[numCorner];
421 << 445
422 G4PolyhedraSideRZ* corn = corners, << 446 G4PolyhedraSideRZ *corn = corners,
423 * sourceCorn = source.corne << 447 *sourceCorn = source.corners;
424 do // Loop checking, 13.08.2015, G.Cosmo << 448 do
425 { 449 {
426 *corn = *sourceCorn; 450 *corn = *sourceCorn;
427 } while( ++sourceCorn, ++corn < corners+numC 451 } while( ++sourceCorn, ++corn < corners+numCorner );
428 << 452
429 // 453 //
430 // Original parameters 454 // Original parameters
431 // 455 //
432 if (source.original_parameters != nullptr) << 456 if (source.original_parameters)
433 { 457 {
434 original_parameters = 458 original_parameters =
435 new G4PolyhedraHistorical( *source.origi 459 new G4PolyhedraHistorical( *source.original_parameters );
436 } 460 }
437 << 461
438 // 462 //
439 // Enclosing cylinder 463 // Enclosing cylinder
440 // 464 //
441 enclosingCylinder = new G4EnclosingCylinder( 465 enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder );
442 466
443 // <<
444 // Surface elements <<
445 // <<
446 delete fElements; <<
447 fElements = nullptr; <<
448 <<
449 // <<
450 // Polyhedron <<
451 // <<
452 fRebuildPolyhedron = false; 467 fRebuildPolyhedron = false;
453 delete fpPolyhedron; << 468 fpPolyhedron = 0;
454 fpPolyhedron = nullptr; <<
455 } 469 }
456 470
>> 471
>> 472 //
457 // Reset 473 // Reset
458 // 474 //
459 // Recalculates and reshapes the solid, given 475 // Recalculates and reshapes the solid, given pre-assigned scaled
460 // original_parameters. 476 // original_parameters.
461 // 477 //
462 G4bool G4Polyhedra::Reset() 478 G4bool G4Polyhedra::Reset()
463 { 479 {
464 if (genericPgon) 480 if (genericPgon)
465 { 481 {
466 std::ostringstream message; 482 std::ostringstream message;
467 message << "Solid " << GetName() << " buil 483 message << "Solid " << GetName() << " built using generic construct."
468 << G4endl << "Not applicable to th 484 << G4endl << "Not applicable to the generic construct !";
469 G4Exception("G4Polyhedra::Reset()", "GeomS 485 G4Exception("G4Polyhedra::Reset()", "GeomSolids1001",
470 JustWarning, message, "Paramet 486 JustWarning, message, "Parameters NOT resetted.");
471 return true; << 487 return 1;
472 } 488 }
473 489
474 // 490 //
475 // Clear old setup 491 // Clear old setup
476 // 492 //
477 G4VCSGfaceted::DeleteStuff(); 493 G4VCSGfaceted::DeleteStuff();
478 delete [] corners; 494 delete [] corners;
479 delete enclosingCylinder; 495 delete enclosingCylinder;
480 delete fElements; <<
481 corners = nullptr; <<
482 fElements = nullptr; <<
483 enclosingCylinder = nullptr; <<
484 496
485 // 497 //
486 // Rebuild polyhedra 498 // Rebuild polyhedra
487 // 499 //
488 auto rz = new G4ReduciblePolygon( original_p << 500 G4ReduciblePolygon *rz =
489 original_p << 501 new G4ReduciblePolygon( original_parameters->Rmin,
490 original_p << 502 original_parameters->Rmax,
491 original_p << 503 original_parameters->Z_values,
>> 504 original_parameters->Num_z_planes );
492 Create( original_parameters->Start_angle, 505 Create( original_parameters->Start_angle,
493 original_parameters->Opening_angle, 506 original_parameters->Opening_angle,
494 original_parameters->numSide, rz ); 507 original_parameters->numSide, rz );
495 delete rz; 508 delete rz;
496 509
497 return false; << 510 return 0;
498 } 511 }
499 512
>> 513
>> 514 //
500 // Inside 515 // Inside
501 // 516 //
502 // This is an override of G4VCSGfaceted::Insid 517 // This is an override of G4VCSGfaceted::Inside, created in order
503 // to speed things up by first checking with G 518 // to speed things up by first checking with G4EnclosingCylinder.
504 // 519 //
505 EInside G4Polyhedra::Inside( const G4ThreeVect << 520 EInside G4Polyhedra::Inside( const G4ThreeVector &p ) const
506 { 521 {
507 // 522 //
508 // Quick test 523 // Quick test
509 // 524 //
510 if (enclosingCylinder->MustBeOutside(p)) ret 525 if (enclosingCylinder->MustBeOutside(p)) return kOutside;
511 526
512 // 527 //
513 // Long answer 528 // Long answer
514 // 529 //
515 return G4VCSGfaceted::Inside(p); 530 return G4VCSGfaceted::Inside(p);
516 } 531 }
517 532
>> 533
>> 534 //
518 // DistanceToIn 535 // DistanceToIn
519 // 536 //
520 // This is an override of G4VCSGfaceted::Insid 537 // This is an override of G4VCSGfaceted::Inside, created in order
521 // to speed things up by first checking with G 538 // to speed things up by first checking with G4EnclosingCylinder.
522 // 539 //
523 G4double G4Polyhedra::DistanceToIn( const G4Th << 540 G4double G4Polyhedra::DistanceToIn( const G4ThreeVector &p,
524 const G4Th << 541 const G4ThreeVector &v ) const
525 { 542 {
526 // 543 //
527 // Quick test 544 // Quick test
528 // 545 //
529 if (enclosingCylinder->ShouldMiss(p,v)) 546 if (enclosingCylinder->ShouldMiss(p,v))
530 return kInfinity; 547 return kInfinity;
531 << 548
532 // 549 //
533 // Long answer 550 // Long answer
534 // 551 //
535 return G4VCSGfaceted::DistanceToIn( p, v ); 552 return G4VCSGfaceted::DistanceToIn( p, v );
536 } 553 }
537 554
>> 555
>> 556 //
538 // DistanceToIn 557 // DistanceToIn
539 // 558 //
540 G4double G4Polyhedra::DistanceToIn( const G4Th << 559 G4double G4Polyhedra::DistanceToIn( const G4ThreeVector &p ) const
541 { 560 {
542 return G4VCSGfaceted::DistanceToIn(p); 561 return G4VCSGfaceted::DistanceToIn(p);
543 } 562 }
544 563
545 // Get bounding box <<
546 // <<
547 void G4Polyhedra::BoundingLimits(G4ThreeVector <<
548 G4ThreeVector <<
549 { <<
550 G4double rmin = kInfinity, rmax = -kInfinity <<
551 G4double zmin = kInfinity, zmax = -kInfinity <<
552 for (G4int i=0; i<GetNumRZCorner(); ++i) <<
553 { <<
554 G4PolyhedraSideRZ corner = GetCorner(i); <<
555 if (corner.r < rmin) rmin = corner.r; <<
556 if (corner.r > rmax) rmax = corner.r; <<
557 if (corner.z < zmin) zmin = corner.z; <<
558 if (corner.z > zmax) zmax = corner.z; <<
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 } <<
612 564
613 // Calculate extent under transform and specif <<
614 // 565 //
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 <<
729 // ComputeDimensions 566 // ComputeDimensions
730 // 567 //
731 void G4Polyhedra::ComputeDimensions( G4V 568 void G4Polyhedra::ComputeDimensions( G4VPVParameterisation* p,
732 const G4i 569 const G4int n,
733 const G4V 570 const G4VPhysicalVolume* pRep )
734 { 571 {
735 p->ComputeDimensions(*this,n,pRep); 572 p->ComputeDimensions(*this,n,pRep);
736 } 573 }
737 574
>> 575
>> 576 //
738 // GetEntityType 577 // GetEntityType
739 // 578 //
740 G4GeometryType G4Polyhedra::GetEntityType() co 579 G4GeometryType G4Polyhedra::GetEntityType() const
741 { 580 {
742 return {"G4Polyhedra"}; << 581 return G4String("G4Polyhedra");
743 } 582 }
744 583
745 // IsFaceted <<
746 // <<
747 G4bool G4Polyhedra::IsFaceted() const <<
748 { <<
749 return true; <<
750 } <<
751 584
>> 585 //
752 // Make a clone of the object 586 // Make a clone of the object
753 // 587 //
754 G4VSolid* G4Polyhedra::Clone() const 588 G4VSolid* G4Polyhedra::Clone() const
755 { 589 {
756 return new G4Polyhedra(*this); 590 return new G4Polyhedra(*this);
757 } 591 }
758 592
>> 593
>> 594 //
759 // Stream object contents to an output stream 595 // Stream object contents to an output stream
760 // 596 //
761 std::ostream& G4Polyhedra::StreamInfo( std::os 597 std::ostream& G4Polyhedra::StreamInfo( std::ostream& os ) const
762 { 598 {
763 G4long oldprc = os.precision(16); << 599 G4int oldprc = os.precision(16);
764 os << "------------------------------------- 600 os << "-----------------------------------------------------------\n"
765 << " *** Dump for solid - " << GetName 601 << " *** Dump for solid - " << GetName() << " ***\n"
766 << " ================================= 602 << " ===================================================\n"
767 << " Solid type: G4Polyhedra\n" 603 << " Solid type: G4Polyhedra\n"
768 << " Parameters: \n" 604 << " Parameters: \n"
769 << " starting phi angle : " << startPh 605 << " starting phi angle : " << startPhi/degree << " degrees \n"
770 << " ending phi angle : " << endPhi/ 606 << " ending phi angle : " << endPhi/degree << " degrees \n"
771 << " number of sides : " << numSide 607 << " number of sides : " << numSide << " \n";
772 G4int i=0; 608 G4int i=0;
773 if (!genericPgon) 609 if (!genericPgon)
774 { 610 {
775 G4int numPlanes = original_parameters->Num 611 G4int numPlanes = original_parameters->Num_z_planes;
776 os << " number of Z planes: " << numPla 612 os << " number of Z planes: " << numPlanes << "\n"
777 << " Z values: \n"; 613 << " Z values: \n";
778 for (i=0; i<numPlanes; ++i) << 614 for (i=0; i<numPlanes; i++)
779 { 615 {
780 os << " Z plane " << i << " 616 os << " Z plane " << i << ": "
781 << original_parameters->Z_values[i] < 617 << original_parameters->Z_values[i] << "\n";
782 } 618 }
783 os << " Tangent distances to 619 os << " Tangent distances to inner surface (Rmin): \n";
784 for (i=0; i<numPlanes; ++i) << 620 for (i=0; i<numPlanes; i++)
785 { 621 {
786 os << " Z plane " << i << " 622 os << " Z plane " << i << ": "
787 << original_parameters->Rmin[i] << "\ 623 << original_parameters->Rmin[i] << "\n";
788 } 624 }
789 os << " Tangent distances to 625 os << " Tangent distances to outer surface (Rmax): \n";
790 for (i=0; i<numPlanes; ++i) << 626 for (i=0; i<numPlanes; i++)
791 { 627 {
792 os << " Z plane " << i << " 628 os << " Z plane " << i << ": "
793 << original_parameters->Rmax[i] << "\ 629 << original_parameters->Rmax[i] << "\n";
794 } 630 }
795 } 631 }
796 os << " number of RZ points: " << numCorn 632 os << " number of RZ points: " << numCorner << "\n"
797 << " RZ values (corners): \n 633 << " RZ values (corners): \n";
798 for (i=0; i<numCorner; ++i) << 634 for (i=0; i<numCorner; i++)
799 { 635 {
800 os << " " 636 os << " "
801 << corners[i].r << ", " << corners[i 637 << corners[i].r << ", " << corners[i].z << "\n";
802 } 638 }
803 os << "------------------------------------- 639 os << "-----------------------------------------------------------\n";
804 os.precision(oldprc); 640 os.precision(oldprc);
805 641
806 return os; 642 return os;
807 } 643 }
808 644
809 ////////////////////////////////////////////// <<
810 // <<
811 // Return volume <<
812 645
813 G4double G4Polyhedra::GetCubicVolume() << 646 //
>> 647 // GetPointOnPlane
>> 648 //
>> 649 // Auxiliary method for get point on surface
>> 650 //
>> 651 G4ThreeVector G4Polyhedra::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1,
>> 652 G4ThreeVector p2, G4ThreeVector p3) const
814 { 653 {
815 if (fCubicVolume == 0.) << 654 G4double lambda1, lambda2, chose,aOne,aTwo;
>> 655 G4ThreeVector t, u, v, w, Area, normal;
>> 656 aOne = 1.;
>> 657 aTwo = 1.;
>> 658
>> 659 t = p1 - p0;
>> 660 u = p2 - p1;
>> 661 v = p3 - p2;
>> 662 w = p0 - p3;
>> 663
>> 664 chose = RandFlat::shoot(0.,aOne+aTwo);
>> 665 if( (chose>=0.) && (chose < aOne) )
816 { 666 {
817 G4double total = 0.; << 667 lambda1 = RandFlat::shoot(0.,1.);
818 G4int nrz = GetNumRZCorner(); << 668 lambda2 = RandFlat::shoot(0.,lambda1);
819 G4PolyhedraSideRZ a = GetCorner(nrz - 1); << 669 return (p2+lambda1*v+lambda2*w);
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 } 670 }
829 return fCubicVolume; << 671
>> 672 lambda1 = RandFlat::shoot(0.,1.);
>> 673 lambda2 = RandFlat::shoot(0.,lambda1);
>> 674 return (p0+lambda1*t+lambda2*u);
830 } 675 }
831 676
832 ////////////////////////////////////////////// << 677
>> 678 //
>> 679 // GetPointOnTriangle
833 // 680 //
834 // Return surface area << 681 // Auxiliary method for get point on surface
>> 682 //
>> 683 G4ThreeVector G4Polyhedra::GetPointOnTriangle(G4ThreeVector p1,
>> 684 G4ThreeVector p2,
>> 685 G4ThreeVector p3) const
>> 686 {
>> 687 G4double lambda1,lambda2;
>> 688 G4ThreeVector v=p3-p1, w=p1-p2;
>> 689
>> 690 lambda1 = RandFlat::shoot(0.,1.);
>> 691 lambda2 = RandFlat::shoot(0.,lambda1);
>> 692
>> 693 return (p2 + lambda1*w + lambda2*v);
>> 694 }
>> 695
835 696
836 G4double G4Polyhedra::GetSurfaceArea() << 697 //
>> 698 // GetPointOnSurface
>> 699 //
>> 700 G4ThreeVector G4Polyhedra::GetPointOnSurface() const
837 { 701 {
838 if (fSurfaceArea == 0.) << 702 if( !genericPgon ) // Polyhedra by faces
839 { 703 {
840 G4double total = 0.; << 704 G4int j, numPlanes = original_parameters->Num_z_planes, Flag=0;
841 G4int nrz = GetNumRZCorner(); << 705 G4double chose, totArea=0., Achose1, Achose2,
842 if (IsOpen()) << 706 rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
>> 707 G4double a, b, l2, rang, totalPhi, ksi,
>> 708 area, aTop=0., aBottom=0., zVal=0.;
>> 709
>> 710 G4ThreeVector p0, p1, p2, p3;
>> 711 std::vector<G4double> aVector1;
>> 712 std::vector<G4double> aVector2;
>> 713 std::vector<G4double> aVector3;
>> 714
>> 715 totalPhi= (phiIsOpen) ? (endPhi-startPhi) : twopi;
>> 716 ksi = totalPhi/numSide;
>> 717 G4double cosksi = std::cos(ksi/2.);
>> 718
>> 719 // Below we generate the areas relevant to our solid
>> 720 //
>> 721 for(j=0; j<numPlanes-1; j++)
843 { 722 {
844 G4PolyhedraSideRZ a = GetCorner(nrz - 1) << 723 a = original_parameters->Rmax[j+1];
845 for (G4int i=0; i<nrz; ++i) << 724 b = original_parameters->Rmax[j];
>> 725 l2 = sqr(original_parameters->Z_values[j]
>> 726 -original_parameters->Z_values[j+1]) + sqr(b-a);
>> 727 area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
>> 728 aVector1.push_back(area);
>> 729 }
>> 730
>> 731 for(j=0; j<numPlanes-1; j++)
>> 732 {
>> 733 a = original_parameters->Rmin[j+1];//*cosksi;
>> 734 b = original_parameters->Rmin[j];//*cosksi;
>> 735 l2 = sqr(original_parameters->Z_values[j]
>> 736 -original_parameters->Z_values[j+1]) + sqr(b-a);
>> 737 area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
>> 738 aVector2.push_back(area);
>> 739 }
>> 740
>> 741 for(j=0; j<numPlanes-1; j++)
>> 742 {
>> 743 if(phiIsOpen == true)
846 { 744 {
847 G4PolyhedraSideRZ b = GetCorner(i); << 745 aVector3.push_back(0.5*(original_parameters->Rmax[j]
848 total += a.r*b.z - a.z*b.r; << 746 -original_parameters->Rmin[j]
849 a = b; << 747 +original_parameters->Rmax[j+1]
>> 748 -original_parameters->Rmin[j+1])
>> 749 *std::fabs(original_parameters->Z_values[j+1]
>> 750 -original_parameters->Z_values[j]));
850 } 751 }
851 total = std::abs(total); << 752 else { aVector3.push_back(0.); }
852 } 753 }
853 G4double alp = (GetEndPhi() - GetStartPhi( << 754
854 G4double cosa = std::cos(alp); << 755 for(j=0; j<numPlanes-1; j++)
855 G4double sina = std::sin(alp); << 756 {
856 G4PolyhedraSideRZ a = GetCorner(nrz - 1); << 757 totArea += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
857 for (G4int i=0; i<nrz; ++i) << 758 }
858 { << 759
859 G4PolyhedraSideRZ b = GetCorner(i); << 760 // Must include top and bottom areas
860 G4ThreeVector p1(a.r, 0, a.z); << 761 //
861 G4ThreeVector p2(a.r*cosa, a.r*sina, a.z << 762 if(original_parameters->Rmax[numPlanes-1] != 0.)
862 G4ThreeVector p3(b.r*cosa, b.r*sina, b.z << 763 {
863 G4ThreeVector p4(b.r, 0, b.z); << 764 a = original_parameters->Rmax[numPlanes-1];
864 total += GetNumSide()*(G4GeomTools::Quad << 765 b = original_parameters->Rmin[numPlanes-1];
865 a = b; << 766 l2 = sqr(a-b);
>> 767 aTop = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
>> 768 }
>> 769
>> 770 if(original_parameters->Rmax[0] != 0.)
>> 771 {
>> 772 a = original_parameters->Rmax[0];
>> 773 b = original_parameters->Rmin[0];
>> 774 l2 = sqr(a-b);
>> 775 aBottom = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
>> 776 }
>> 777
>> 778 Achose1 = 0.;
>> 779 Achose2 = numSide*(aVector1[0]+aVector2[0])+2.*aVector3[0];
>> 780
>> 781 chose = RandFlat::shoot(0.,totArea+aTop+aBottom);
>> 782 if( (chose >= 0.) && (chose < aTop + aBottom) )
>> 783 {
>> 784 chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
>> 785 rang = std::floor((chose-startPhi)/ksi-0.01);
>> 786 if(rang<0) { rang=0; }
>> 787 rang = std::fabs(rang);
>> 788 sinphi1 = std::sin(startPhi+rang*ksi);
>> 789 sinphi2 = std::sin(startPhi+(rang+1)*ksi);
>> 790 cosphi1 = std::cos(startPhi+rang*ksi);
>> 791 cosphi2 = std::cos(startPhi+(rang+1)*ksi);
>> 792 chose = RandFlat::shoot(0., aTop + aBottom);
>> 793 if(chose>=0. && chose<aTop)
>> 794 {
>> 795 rad1 = original_parameters->Rmin[numPlanes-1];
>> 796 rad2 = original_parameters->Rmax[numPlanes-1];
>> 797 zVal = original_parameters->Z_values[numPlanes-1];
>> 798 }
>> 799 else
>> 800 {
>> 801 rad1 = original_parameters->Rmin[0];
>> 802 rad2 = original_parameters->Rmax[0];
>> 803 zVal = original_parameters->Z_values[0];
>> 804 }
>> 805 p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
>> 806 p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
>> 807 p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
>> 808 p3 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
>> 809 return GetPointOnPlane(p0,p1,p2,p3);
>> 810 }
>> 811 else
>> 812 {
>> 813 for (j=0; j<numPlanes-1; j++)
>> 814 {
>> 815 if( ((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes-2) )
>> 816 {
>> 817 Flag = j; break;
>> 818 }
>> 819 Achose1 += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
>> 820 Achose2 = Achose1 + numSide*(aVector1[j+1]+aVector2[j+1])
>> 821 + 2.*aVector3[j+1];
>> 822 }
866 } 823 }
867 fSurfaceArea = total; <<
868 } <<
869 return fSurfaceArea; <<
870 } <<
871 824
872 ////////////////////////////////////////////// << 825 // At this point we have chosen a subsection
873 // << 826 // between to adjacent plane cuts...
874 // Set vector of surface elements, auxiliary m <<
875 // random points on surface <<
876 827
877 void G4Polyhedra::SetSurfaceElements() const << 828 j = Flag;
878 { << 829
879 fElements = new std::vector<G4Polyhedra::sur << 830 totArea = numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
880 G4double total = 0.; << 831 chose = RandFlat::shoot(0.,totArea);
881 G4int nrz = GetNumRZCorner(); << 832
882 << 833 if( (chose>=0.) && (chose<numSide*aVector1[j]) )
883 // set lateral surface elements << 834 {
884 G4double dphi = (GetEndPhi() - GetStartPhi() << 835 chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
885 G4double cosa = std::cos(dphi); << 836 rang = std::floor((chose-startPhi)/ksi-0.01);
886 G4double sina = std::sin(dphi); << 837 if(rang<0) { rang=0; }
887 G4int ia = nrz - 1; << 838 rang = std::fabs(rang);
888 for (G4int ib=0; ib<nrz; ++ib) << 839 rad1 = original_parameters->Rmax[j];
889 { << 840 rad2 = original_parameters->Rmax[j+1];
890 G4PolyhedraSideRZ a = GetCorner(ia); << 841 sinphi1 = std::sin(startPhi+rang*ksi);
891 G4PolyhedraSideRZ b = GetCorner(ib); << 842 sinphi2 = std::sin(startPhi+(rang+1)*ksi);
892 G4Polyhedra::surface_element selem; << 843 cosphi1 = std::cos(startPhi+rang*ksi);
893 selem.i0 = ia; << 844 cosphi2 = std::cos(startPhi+(rang+1)*ksi);
894 selem.i1 = ib; << 845 zVal = original_parameters->Z_values[j];
895 ia = ib; << 846
896 if (a.r == 0. && b.r == 0.) continue; << 847 p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
897 G4ThreeVector p1(a.r, 0, a.z); << 848 p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
898 G4ThreeVector p2(a.r*cosa, a.r*sina, a.z); << 849
899 G4ThreeVector p3(b.r*cosa, b.r*sina, b.z); << 850 zVal = original_parameters->Z_values[j+1];
900 G4ThreeVector p4(b.r, 0, b.z); << 851
901 if (a.r > 0.) << 852 p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
902 { << 853 p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
903 selem.i2 = -1; << 854 return GetPointOnPlane(p0,p1,p2,p3);
904 total += GetNumSide()*(G4GeomTools::Tria << 855 }
905 selem.area = total; << 856 else if ( (chose >= numSide*aVector1[j])
906 fElements->push_back(selem); << 857 && (chose <= numSide*(aVector1[j]+aVector2[j])) )
907 } << 858 {
908 if (b.r > 0.) << 859 chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
909 { << 860 rang = std::floor((chose-startPhi)/ksi-0.01);
910 selem.i2 = -2; << 861 if(rang<0) { rang=0; }
911 total += GetNumSide()*(G4GeomTools::Tria << 862 rang = std::fabs(rang);
912 selem.area = total; << 863 rad1 = original_parameters->Rmin[j];
913 fElements->push_back(selem); << 864 rad2 = original_parameters->Rmin[j+1];
914 } << 865 sinphi1 = std::sin(startPhi+rang*ksi);
915 } << 866 sinphi2 = std::sin(startPhi+(rang+1)*ksi);
916 << 867 cosphi1 = std::cos(startPhi+rang*ksi);
917 // set elements for phi cuts << 868 cosphi2 = std::cos(startPhi+(rang+1)*ksi);
918 if (IsOpen()) << 869 zVal = original_parameters->Z_values[j];
919 { << 870
920 G4TwoVectorList contourRZ; << 871 p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
921 std::vector<G4int> triangles; << 872 p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
922 for (G4int i=0; i<nrz; ++i) << 873
923 { << 874 zVal = original_parameters->Z_values[j+1];
924 G4PolyhedraSideRZ corner = GetCorner(i); << 875
925 contourRZ.emplace_back(corner.r, corner. << 876 p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
926 } << 877 p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
927 G4GeomTools::TriangulatePolygon(contourRZ, << 878 return GetPointOnPlane(p0,p1,p2,p3);
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 } 879 }
948 } <<
949 } <<
950 880
951 ////////////////////////////////////////////// << 881 chose = RandFlat::shoot(0.,2.2);
952 // << 882 if( (chose>=0.) && (chose < 1.) )
953 // Generate random point on surface << 883 {
954 << 884 rang = startPhi;
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 } 885 }
1005 else 886 else
1006 { 887 {
1007 if (iside == nside) --iside; // iside m << 888 rang = endPhi;
1008 G4double phi = iside*dphi + GetStartPhi << 889 }
1009 G4double cosphi = std::cos(phi); << 890
1010 G4double sinphi = std::sin(phi); << 891 cosphi1 = std::cos(rang); rad1 = original_parameters->Rmin[j];
1011 x = p0.x()*cosphi - p0.y()*sinphi; << 892 sinphi1 = std::sin(rang); rad2 = original_parameters->Rmax[j];
1012 y = p0.x()*sinphi + p0.y()*cosphi; << 893
1013 z = p0.z(); << 894 p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
1014 } << 895 original_parameters->Z_values[j]);
1015 } << 896 p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
1016 else // phi cut << 897 original_parameters->Z_values[j]);
1017 { << 898
1018 G4int nrz = GetNumRZCorner(); << 899 rad1 = original_parameters->Rmax[j+1];
1019 G4double phi = (i0 < nrz) ? GetStartPhi() << 900 rad2 = original_parameters->Rmin[j+1];
1020 if (i0 >= nrz) { i0 -= nrz; } << 901
1021 G4PolyhedraSideRZ p0 = GetCorner(i0); << 902 p2 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
1022 G4PolyhedraSideRZ p1 = GetCorner(i1); << 903 original_parameters->Z_values[j+1]);
1023 G4PolyhedraSideRZ p2 = GetCorner(i2); << 904 p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
1024 G4double r = (p1.r - p0.r)*u + (p2.r - p0 << 905 original_parameters->Z_values[j+1]);
1025 x = r*std::cos(phi); << 906 return GetPointOnPlane(p0,p1,p2,p3);
1026 y = r*std::sin(phi); << 907 }
1027 z = (p1.z - p0.z)*u + (p2.z - p0.z)*v + p << 908 else // Generic polyhedra
>> 909 {
>> 910 return GetPointOnSurfaceGeneric();
1028 } 911 }
1029 return {x, y, z}; <<
1030 } 912 }
1031 913
1032 ///////////////////////////////////////////// <<
1033 // 914 //
1034 // CreatePolyhedron 915 // CreatePolyhedron
1035 << 916 //
1036 G4Polyhedron* G4Polyhedra::CreatePolyhedron() 917 G4Polyhedron* G4Polyhedra::CreatePolyhedron() const
1037 { << 918 {
1038 std::vector<G4TwoVector> rz(numCorner); << 919 if (!genericPgon)
1039 for (G4int i = 0; i < numCorner; ++i) << 920 {
1040 rz[i].set(corners[i].r, corners[i].z); << 921 return new G4PolyhedronPgon( original_parameters->Start_angle,
1041 return new G4PolyhedronPgon(startPhi, endPh << 922 original_parameters->Opening_angle,
>> 923 original_parameters->numSide,
>> 924 original_parameters->Num_z_planes,
>> 925 original_parameters->Z_values,
>> 926 original_parameters->Rmin,
>> 927 original_parameters->Rmax);
>> 928 }
>> 929 else
>> 930 {
>> 931 // The following code prepares for:
>> 932 // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
>> 933 // const double xyz[][3],
>> 934 // const int faces_vec[][4])
>> 935 // Here is an extract from the header file HepPolyhedron.h:
>> 936 /**
>> 937 * Creates user defined polyhedron.
>> 938 * This function allows to the user to define arbitrary polyhedron.
>> 939 * The faces of the polyhedron should be either triangles or planar
>> 940 * quadrilateral. Nodes of a face are defined by indexes pointing to
>> 941 * the elements in the xyz array. Numeration of the elements in the
>> 942 * array starts from 1 (like in fortran). The indexes can be positive
>> 943 * or negative. Negative sign means that the corresponding edge is
>> 944 * invisible. The normal of the face should be directed to exterior
>> 945 * of the polyhedron.
>> 946 *
>> 947 * @param Nnodes number of nodes
>> 948 * @param Nfaces number of faces
>> 949 * @param xyz nodes
>> 950 * @param faces_vec faces (quadrilaterals or triangles)
>> 951 * @return status of the operation - is non-zero in case of problem
>> 952 */
>> 953 G4int nNodes;
>> 954 G4int nFaces;
>> 955 typedef G4double double3[3];
>> 956 double3* xyz;
>> 957 typedef G4int int4[4];
>> 958 int4* faces_vec;
>> 959 if (phiIsOpen)
>> 960 {
>> 961 // Triangulate open ends. Simple ear-chopping algorithm...
>> 962 // I'm not sure how robust this algorithm is (J.Allison).
>> 963 //
>> 964 std::vector<G4bool> chopped(numCorner, false);
>> 965 std::vector<G4int*> triQuads;
>> 966 G4int remaining = numCorner;
>> 967 G4int iStarter = 0;
>> 968 while (remaining >= 3)
>> 969 {
>> 970 // Find unchopped corners...
>> 971 //
>> 972 G4int A = -1, B = -1, C = -1;
>> 973 G4int iStepper = iStarter;
>> 974 do
>> 975 {
>> 976 if (A < 0) { A = iStepper; }
>> 977 else if (B < 0) { B = iStepper; }
>> 978 else if (C < 0) { C = iStepper; }
>> 979 do
>> 980 {
>> 981 if (++iStepper >= numCorner) iStepper = 0;
>> 982 }
>> 983 while (chopped[iStepper]);
>> 984 }
>> 985 while (C < 0 && iStepper != iStarter);
>> 986
>> 987 // Check triangle at B is pointing outward (an "ear").
>> 988 // Sign of z cross product determines...
>> 989
>> 990 G4double BAr = corners[A].r - corners[B].r;
>> 991 G4double BAz = corners[A].z - corners[B].z;
>> 992 G4double BCr = corners[C].r - corners[B].r;
>> 993 G4double BCz = corners[C].z - corners[B].z;
>> 994 if (BAr * BCz - BAz * BCr < kCarTolerance)
>> 995 {
>> 996 G4int* tq = new G4int[3];
>> 997 tq[0] = A + 1;
>> 998 tq[1] = B + 1;
>> 999 tq[2] = C + 1;
>> 1000 triQuads.push_back(tq);
>> 1001 chopped[B] = true;
>> 1002 --remaining;
>> 1003 }
>> 1004 else
>> 1005 {
>> 1006 do
>> 1007 {
>> 1008 if (++iStarter >= numCorner) { iStarter = 0; }
>> 1009 }
>> 1010 while (chopped[iStarter]);
>> 1011 }
>> 1012 }
>> 1013
>> 1014 // Transfer to faces...
>> 1015
>> 1016 nNodes = (numSide + 1) * numCorner;
>> 1017 nFaces = numSide * numCorner + 2 * triQuads.size();
>> 1018 faces_vec = new int4[nFaces];
>> 1019 G4int iface = 0;
>> 1020 G4int addition = numCorner * numSide;
>> 1021 G4int d = numCorner - 1;
>> 1022 for (G4int iEnd = 0; iEnd < 2; ++iEnd)
>> 1023 {
>> 1024 for (size_t i = 0; i < triQuads.size(); ++i)
>> 1025 {
>> 1026 // Negative for soft/auxiliary/normally invisible edges...
>> 1027 //
>> 1028 G4int a, b, c;
>> 1029 if (iEnd == 0)
>> 1030 {
>> 1031 a = triQuads[i][0];
>> 1032 b = triQuads[i][1];
>> 1033 c = triQuads[i][2];
>> 1034 }
>> 1035 else
>> 1036 {
>> 1037 a = triQuads[i][0] + addition;
>> 1038 b = triQuads[i][2] + addition;
>> 1039 c = triQuads[i][1] + addition;
>> 1040 }
>> 1041 G4int ab = std::abs(b - a);
>> 1042 G4int bc = std::abs(c - b);
>> 1043 G4int ca = std::abs(a - c);
>> 1044 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
>> 1045 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
>> 1046 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
>> 1047 faces_vec[iface][3] = 0;
>> 1048 ++iface;
>> 1049 }
>> 1050 }
>> 1051
>> 1052 // Continue with sides...
>> 1053
>> 1054 xyz = new double3[nNodes];
>> 1055 const G4double dPhi = (endPhi - startPhi) / numSide;
>> 1056 G4double phi = startPhi;
>> 1057 G4int ixyz = 0;
>> 1058 for (G4int iSide = 0; iSide < numSide; ++iSide)
>> 1059 {
>> 1060 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
>> 1061 {
>> 1062 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
>> 1063 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
>> 1064 xyz[ixyz][2] = corners[iCorner].z;
>> 1065 if (iCorner < numCorner - 1)
>> 1066 {
>> 1067 faces_vec[iface][0] = ixyz + 1;
>> 1068 faces_vec[iface][1] = ixyz + numCorner + 1;
>> 1069 faces_vec[iface][2] = ixyz + numCorner + 2;
>> 1070 faces_vec[iface][3] = ixyz + 2;
>> 1071 }
>> 1072 else
>> 1073 {
>> 1074 faces_vec[iface][0] = ixyz + 1;
>> 1075 faces_vec[iface][1] = ixyz + numCorner + 1;
>> 1076 faces_vec[iface][2] = ixyz + 2;
>> 1077 faces_vec[iface][3] = ixyz - numCorner + 2;
>> 1078 }
>> 1079 ++iface;
>> 1080 ++ixyz;
>> 1081 }
>> 1082 phi += dPhi;
>> 1083 }
>> 1084
>> 1085 // Last corners...
>> 1086
>> 1087 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
>> 1088 {
>> 1089 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
>> 1090 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
>> 1091 xyz[ixyz][2] = corners[iCorner].z;
>> 1092 ++ixyz;
>> 1093 }
>> 1094 }
>> 1095 else // !phiIsOpen - i.e., a complete 360 degrees.
>> 1096 {
>> 1097 nNodes = numSide * numCorner;
>> 1098 nFaces = numSide * numCorner;;
>> 1099 xyz = new double3[nNodes];
>> 1100 faces_vec = new int4[nFaces];
>> 1101 // const G4double dPhi = (endPhi - startPhi) / numSide;
>> 1102 const G4double dPhi = twopi / numSide; // !phiIsOpen endPhi-startPhi = 360 degrees.
>> 1103 G4double phi = startPhi;
>> 1104 G4int ixyz = 0, iface = 0;
>> 1105 for (G4int iSide = 0; iSide < numSide; ++iSide)
>> 1106 {
>> 1107 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
>> 1108 {
>> 1109 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
>> 1110 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
>> 1111 xyz[ixyz][2] = corners[iCorner].z;
>> 1112 if (iSide < numSide - 1)
>> 1113 {
>> 1114 if (iCorner < numCorner - 1)
>> 1115 {
>> 1116 faces_vec[iface][0] = ixyz + 1;
>> 1117 faces_vec[iface][1] = ixyz + numCorner + 1;
>> 1118 faces_vec[iface][2] = ixyz + numCorner + 2;
>> 1119 faces_vec[iface][3] = ixyz + 2;
>> 1120 }
>> 1121 else
>> 1122 {
>> 1123 faces_vec[iface][0] = ixyz + 1;
>> 1124 faces_vec[iface][1] = ixyz + numCorner + 1;
>> 1125 faces_vec[iface][2] = ixyz + 2;
>> 1126 faces_vec[iface][3] = ixyz - numCorner + 2;
>> 1127 }
>> 1128 }
>> 1129 else // Last side joins ends...
>> 1130 {
>> 1131 if (iCorner < numCorner - 1)
>> 1132 {
>> 1133 faces_vec[iface][0] = ixyz + 1;
>> 1134 faces_vec[iface][1] = ixyz + numCorner - nFaces + 1;
>> 1135 faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
>> 1136 faces_vec[iface][3] = ixyz + 2;
>> 1137 }
>> 1138 else
>> 1139 {
>> 1140 faces_vec[iface][0] = ixyz + 1;
>> 1141 faces_vec[iface][1] = ixyz - nFaces + numCorner + 1;
>> 1142 faces_vec[iface][2] = ixyz - nFaces + 2;
>> 1143 faces_vec[iface][3] = ixyz - numCorner + 2;
>> 1144 }
>> 1145 }
>> 1146 ++ixyz;
>> 1147 ++iface;
>> 1148 }
>> 1149 phi += dPhi;
>> 1150 }
>> 1151 }
>> 1152 G4Polyhedron* polyhedron = new G4Polyhedron;
>> 1153 G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
>> 1154 delete [] faces_vec;
>> 1155 delete [] xyz;
>> 1156 if (problem)
>> 1157 {
>> 1158 std::ostringstream message;
>> 1159 message << "Problem creating G4Polyhedron for: " << GetName();
>> 1160 G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
>> 1161 JustWarning, message);
>> 1162 delete polyhedron;
>> 1163 return 0;
>> 1164 }
>> 1165 else
>> 1166 {
>> 1167 return polyhedron;
>> 1168 }
>> 1169 }
1042 } 1170 }
1043 1171
1044 // SetOriginalParameters << 1172
1045 // << 1173 void G4Polyhedra::SetOriginalParameters(G4ReduciblePolygon *rz)
1046 void G4Polyhedra::SetOriginalParameters(G4Red <<
1047 { 1174 {
1048 G4int numPlanes = numCorner; << 1175 G4int numPlanes = (G4int)numCorner;
1049 G4bool isConvertible = true; << 1176 G4bool isConvertible=true;
1050 G4double Zmax=rz->Bmax(); 1177 G4double Zmax=rz->Bmax();
1051 rz->StartWithZMin(); 1178 rz->StartWithZMin();
1052 1179
1053 // Prepare vectors for storage << 1180 // Prepare vectors for storage
1054 // 1181 //
1055 std::vector<G4double> Z; 1182 std::vector<G4double> Z;
1056 std::vector<G4double> Rmin; 1183 std::vector<G4double> Rmin;
1057 std::vector<G4double> Rmax; 1184 std::vector<G4double> Rmax;
1058 1185
1059 G4int countPlanes=1; 1186 G4int countPlanes=1;
1060 G4int icurr=0; 1187 G4int icurr=0;
1061 G4int icurl=0; 1188 G4int icurl=0;
1062 1189
1063 // first plane Z=Z[0] 1190 // first plane Z=Z[0]
1064 // 1191 //
1065 Z.push_back(corners[0].z); 1192 Z.push_back(corners[0].z);
1066 G4double Zprev=Z[0]; 1193 G4double Zprev=Z[0];
1067 if (Zprev == corners[1].z) 1194 if (Zprev == corners[1].z)
1068 { 1195 {
1069 Rmin.push_back(corners[0].r); << 1196 Rmin.push_back(corners[0].r);
1070 Rmax.push_back (corners[1].r);icurr=1; << 1197 Rmax.push_back (corners[1].r);icurr=1;
1071 } 1198 }
1072 else if (Zprev == corners[numPlanes-1].z) 1199 else if (Zprev == corners[numPlanes-1].z)
1073 { 1200 {
1074 Rmin.push_back(corners[numPlanes-1].r); << 1201 Rmin.push_back(corners[numPlanes-1].r);
1075 Rmax.push_back (corners[0].r); 1202 Rmax.push_back (corners[0].r);
1076 icurl=numPlanes-1; << 1203 icurl=numPlanes-1;
1077 } 1204 }
1078 else 1205 else
1079 { 1206 {
1080 Rmin.push_back(corners[0].r); << 1207 Rmin.push_back(corners[0].r);
1081 Rmax.push_back (corners[0].r); 1208 Rmax.push_back (corners[0].r);
1082 } 1209 }
1083 1210
1084 // next planes until last 1211 // next planes until last
1085 // 1212 //
1086 G4int inextr=0, inextl=0; << 1213 G4int inextr=0, inextl=0;
1087 for (G4int i=0; i < numPlanes-2; ++i) << 1214 for (G4int i=0; i < numPlanes-2; i++)
1088 { 1215 {
1089 inextr=1+icurr; 1216 inextr=1+icurr;
1090 inextl=(icurl <= 0)? numPlanes-1 : icurl- 1217 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1091 1218
1092 if((corners[inextr].z >= Zmax) & (corners 1219 if((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax)) { break; }
1093 1220
1094 G4double Zleft = corners[inextl].z; 1221 G4double Zleft = corners[inextl].z;
1095 G4double Zright = corners[inextr].z; 1222 G4double Zright = corners[inextr].z;
1096 if(Zright>Zleft) 1223 if(Zright>Zleft)
1097 { 1224 {
1098 Z.push_back(Zleft); << 1225 Z.push_back(Zleft);
1099 countPlanes++; 1226 countPlanes++;
1100 G4double difZr=corners[inextr].z - corn 1227 G4double difZr=corners[inextr].z - corners[icurr].z;
1101 G4double difZl=corners[inextl].z - corn 1228 G4double difZl=corners[inextl].z - corners[icurl].z;
1102 1229
1103 if(std::fabs(difZl) < kCarTolerance) 1230 if(std::fabs(difZl) < kCarTolerance)
1104 { 1231 {
1105 if(std::fabs(difZr) < kCarTolerance) 1232 if(std::fabs(difZr) < kCarTolerance)
1106 { 1233 {
1107 Rmin.push_back(corners[inextl].r); 1234 Rmin.push_back(corners[inextl].r);
1108 Rmax.push_back(corners[icurr].r); 1235 Rmax.push_back(corners[icurr].r);
1109 } 1236 }
1110 else 1237 else
1111 { 1238 {
1112 Rmin.push_back(corners[inextl].r); 1239 Rmin.push_back(corners[inextl].r);
1113 Rmax.push_back(corners[icurr].r + ( 1240 Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1114 *(corners[ine << 1241 *(corners[inextr].r - corners[icurr].r));
1115 } 1242 }
1116 } 1243 }
1117 else if (difZl >= kCarTolerance) 1244 else if (difZl >= kCarTolerance)
1118 { 1245 {
1119 if(std::fabs(difZr) < kCarTolerance) 1246 if(std::fabs(difZr) < kCarTolerance)
1120 { 1247 {
1121 Rmin.push_back(corners[icurl].r); 1248 Rmin.push_back(corners[icurl].r);
1122 Rmax.push_back(corners[icurr].r); 1249 Rmax.push_back(corners[icurr].r);
1123 } 1250 }
1124 else 1251 else
1125 { 1252 {
1126 Rmin.push_back(corners[icurl].r); 1253 Rmin.push_back(corners[icurl].r);
1127 Rmax.push_back(corners[icurr].r + ( 1254 Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1128 *(corners[ine 1255 *(corners[inextr].r - corners[icurr].r));
1129 } 1256 }
1130 } 1257 }
1131 else 1258 else
1132 { 1259 {
1133 isConvertible=false; break; 1260 isConvertible=false; break;
1134 } 1261 }
1135 icurl=(icurl == 0)? numPlanes-1 : icurl 1262 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1136 } 1263 }
1137 else if(std::fabs(Zright-Zleft)<kCarToler 1264 else if(std::fabs(Zright-Zleft)<kCarTolerance) // Zright=Zleft
1138 { 1265 {
1139 Z.push_back(Zleft); << 1266 Z.push_back(Zleft);
1140 ++countPlanes; << 1267 countPlanes++;
1141 ++icurr; << 1268 icurr++;
1142 1269
1143 icurl=(icurl == 0)? numPlanes-1 : icurl 1270 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1144 1271
1145 Rmin.push_back(corners[inextl].r); << 1272 Rmin.push_back(corners[inextl].r);
1146 Rmax.push_back (corners[inextr].r); 1273 Rmax.push_back (corners[inextr].r);
1147 } 1274 }
1148 else // Zright<Zleft 1275 else // Zright<Zleft
1149 { 1276 {
1150 Z.push_back(Zright); << 1277 Z.push_back(Zright);
1151 ++countPlanes; << 1278 countPlanes++;
1152 1279
1153 G4double difZr=corners[inextr].z - corn 1280 G4double difZr=corners[inextr].z - corners[icurr].z;
1154 G4double difZl=corners[inextl].z - corn 1281 G4double difZl=corners[inextl].z - corners[icurl].z;
1155 if(std::fabs(difZr) < kCarTolerance) 1282 if(std::fabs(difZr) < kCarTolerance)
1156 { 1283 {
1157 if(std::fabs(difZl) < kCarTolerance) 1284 if(std::fabs(difZl) < kCarTolerance)
1158 { 1285 {
1159 Rmax.push_back(corners[inextr].r); 1286 Rmax.push_back(corners[inextr].r);
1160 Rmin.push_back(corners[icurr].r); << 1287 Rmin.push_back(corners[icurr].r);
1161 } << 1288 }
1162 else 1289 else
1163 { 1290 {
1164 Rmin.push_back(corners[icurl].r + ( 1291 Rmin.push_back(corners[icurl].r + (Zright-corners[icurl].z)/difZl
1165 * (corners[in 1292 * (corners[inextl].r - corners[icurl].r));
1166 Rmax.push_back(corners[inextr].r); 1293 Rmax.push_back(corners[inextr].r);
1167 } 1294 }
1168 ++icurr; << 1295 icurr++;
1169 } // plate 1296 } // plate
1170 else if (difZr >= kCarTolerance) 1297 else if (difZr >= kCarTolerance)
1171 { 1298 {
1172 if(std::fabs(difZl) < kCarTolerance) 1299 if(std::fabs(difZl) < kCarTolerance)
1173 { 1300 {
1174 Rmax.push_back(corners[inextr].r); 1301 Rmax.push_back(corners[inextr].r);
1175 Rmin.push_back (corners[icurr].r); << 1302 Rmin.push_back (corners[icurr].r);
1176 } << 1303 }
1177 else 1304 else
1178 { 1305 {
1179 Rmax.push_back(corners[inextr].r); 1306 Rmax.push_back(corners[inextr].r);
1180 Rmin.push_back (corners[icurl].r+(Z 1307 Rmin.push_back (corners[icurl].r+(Zright-corners[icurl].z)/difZl
1181 * (corners[ 1308 * (corners[inextl].r - corners[icurl].r));
1182 } 1309 }
1183 ++icurr; << 1310 icurr++;
1184 } 1311 }
1185 else 1312 else
1186 { 1313 {
1187 isConvertible=false; break; 1314 isConvertible=false; break;
1188 } 1315 }
1189 } 1316 }
1190 } // end for loop 1317 } // end for loop
1191 1318
1192 // last plane Z=Zmax 1319 // last plane Z=Zmax
1193 // 1320 //
1194 Z.push_back(Zmax); 1321 Z.push_back(Zmax);
1195 ++countPlanes; << 1322 countPlanes++;
1196 inextr=1+icurr; 1323 inextr=1+icurr;
1197 inextl=(icurl <= 0)? numPlanes-1 : icurl-1; 1324 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1198 << 1325
1199 Rmax.push_back(corners[inextr].r); << 1326 if(corners[inextr].z==corners[inextl].z)
1200 Rmin.push_back(corners[inextl].r); << 1327 {
>> 1328 Rmax.push_back(corners[inextr].r);
>> 1329 Rmin.push_back(corners[inextl].r);
>> 1330 }
>> 1331 else
>> 1332 {
>> 1333 Rmax.push_back(corners[inextr].r);
>> 1334 Rmin.push_back(corners[inextl].r);
>> 1335 }
1201 1336
1202 // Set original parameters Rmin,Rmax,Z 1337 // Set original parameters Rmin,Rmax,Z
1203 // 1338 //
1204 if(isConvertible) 1339 if(isConvertible)
1205 { 1340 {
1206 original_parameters = new G4PolyhedraHisto 1341 original_parameters = new G4PolyhedraHistorical;
1207 original_parameters->numSide = numSide; 1342 original_parameters->numSide = numSide;
1208 original_parameters->Z_values = new G4doub 1343 original_parameters->Z_values = new G4double[countPlanes];
1209 original_parameters->Rmin = new G4double[c 1344 original_parameters->Rmin = new G4double[countPlanes];
1210 original_parameters->Rmax = new G4double[c 1345 original_parameters->Rmax = new G4double[countPlanes];
1211 << 1346
1212 for(G4int j=0; j < countPlanes; ++j) << 1347 for(G4int j=0; j < countPlanes; j++)
1213 { 1348 {
1214 original_parameters->Z_values[j] = Z[j]; 1349 original_parameters->Z_values[j] = Z[j];
1215 original_parameters->Rmax[j] = Rmax[j]; 1350 original_parameters->Rmax[j] = Rmax[j];
1216 original_parameters->Rmin[j] = Rmin[j]; 1351 original_parameters->Rmin[j] = Rmin[j];
1217 } 1352 }
1218 original_parameters->Start_angle = startPh 1353 original_parameters->Start_angle = startPhi;
1219 original_parameters->Opening_angle = endPh 1354 original_parameters->Opening_angle = endPhi-startPhi;
1220 original_parameters->Num_z_planes = countP 1355 original_parameters->Num_z_planes = countPlanes;
1221 << 1356
1222 } 1357 }
1223 else // Set parameters(r,z) with Rmin==0 a 1358 else // Set parameters(r,z) with Rmin==0 as convention
1224 { 1359 {
1225 #ifdef G4SPECSDEBUG 1360 #ifdef G4SPECSDEBUG
1226 std::ostringstream message; 1361 std::ostringstream message;
1227 message << "Polyhedra " << GetName() << G 1362 message << "Polyhedra " << GetName() << G4endl
1228 << "cannot be converted to Polyhedra wi 1363 << "cannot be converted to Polyhedra with (Rmin,Rmaz,Z) parameters!";
1229 G4Exception("G4Polyhedra::SetOriginalPara 1364 G4Exception("G4Polyhedra::SetOriginalParameters()",
1230 "GeomSolids0002", JustWarning 1365 "GeomSolids0002", JustWarning, message);
1231 #endif 1366 #endif
1232 original_parameters = new G4PolyhedraHist 1367 original_parameters = new G4PolyhedraHistorical;
1233 original_parameters->numSide = numSide; 1368 original_parameters->numSide = numSide;
1234 original_parameters->Z_values = new G4dou 1369 original_parameters->Z_values = new G4double[numPlanes];
1235 original_parameters->Rmin = new G4double[ 1370 original_parameters->Rmin = new G4double[numPlanes];
1236 original_parameters->Rmax = new G4double[ 1371 original_parameters->Rmax = new G4double[numPlanes];
1237 << 1372
1238 for(G4int j=0; j < numPlanes; ++j) << 1373 for(G4int j=0; j < numPlanes; j++)
1239 { 1374 {
1240 original_parameters->Z_values[j] = corn 1375 original_parameters->Z_values[j] = corners[j].z;
1241 original_parameters->Rmax[j] = corners[ 1376 original_parameters->Rmax[j] = corners[j].r;
1242 original_parameters->Rmin[j] = 0.0; 1377 original_parameters->Rmin[j] = 0.0;
1243 } 1378 }
1244 original_parameters->Start_angle = startP 1379 original_parameters->Start_angle = startPhi;
1245 original_parameters->Opening_angle = endP 1380 original_parameters->Opening_angle = endPhi-startPhi;
1246 original_parameters->Num_z_planes = numPl 1381 original_parameters->Num_z_planes = numPlanes;
1247 } 1382 }
>> 1383 //return isConvertible;
1248 } 1384 }
1249 1385
1250 #endif 1386 #endif
1251 1387