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