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