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