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