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