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