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