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