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