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 // G4GenericPolycone implementation << 27 // 26 // 28 // Authors: T.Nikitina, G.Cosmo - CERN << 27 // $Id: G4GenericPolycone.cc 72131 2013-07-11 07:12:24Z gcosmo $ >> 28 // >> 29 // >> 30 // -------------------------------------------------------------------- >> 31 // GEANT 4 class source file >> 32 // >> 33 // >> 34 // G4GenericPolycone.cc >> 35 // >> 36 // Implementation of a CSG polycone >> 37 // 29 // ------------------------------------------- 38 // -------------------------------------------------------------------- 30 39 31 #include "G4GenericPolycone.hh" 40 #include "G4GenericPolycone.hh" 32 41 33 #if !defined(G4GEOM_USE_UGENERICPOLYCONE) 42 #if !defined(G4GEOM_USE_UGENERICPOLYCONE) 34 43 35 #include "G4PolyconeSide.hh" 44 #include "G4PolyconeSide.hh" 36 #include "G4PolyPhiFace.hh" 45 #include "G4PolyPhiFace.hh" 37 46 38 #include "G4GeomTools.hh" << 47 #include "Randomize.hh" 39 #include "G4VoxelLimits.hh" << 40 #include "G4AffineTransform.hh" << 41 #include "G4BoundingEnvelope.hh" << 42 << 43 #include "G4QuickRand.hh" << 44 48 45 #include "G4Polyhedron.hh" 49 #include "G4Polyhedron.hh" 46 #include "G4EnclosingCylinder.hh" 50 #include "G4EnclosingCylinder.hh" 47 #include "G4ReduciblePolygon.hh" 51 #include "G4ReduciblePolygon.hh" 48 #include "G4VPVParameterisation.hh" 52 #include "G4VPVParameterisation.hh" 49 53 50 namespace << 51 { << 52 G4Mutex surface_elementsMutex = G4MUTEX_INIT << 53 } << 54 << 55 using namespace CLHEP; 54 using namespace CLHEP; 56 55 >> 56 // 57 // Constructor (generic parameters) 57 // Constructor (generic parameters) 58 // 58 // 59 G4GenericPolycone::G4GenericPolycone( const G4 << 59 G4GenericPolycone::G4GenericPolycone( const G4String& name, 60 G4double phiStar 60 G4double phiStart, 61 G4double phiTota 61 G4double phiTotal, 62 G4int numRZ, 62 G4int numRZ, 63 const G4double r[], 63 const G4double r[], 64 const G4double z[] ) 64 const G4double z[] ) 65 : G4VCSGfaceted( name ) 65 : G4VCSGfaceted( name ) 66 { << 66 { 67 << 67 68 auto rz = new G4ReduciblePolygon( r, z, numR << 68 G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ ); 69 << 69 70 Create( phiStart, phiTotal, rz ); 70 Create( phiStart, phiTotal, rz ); 71 << 71 72 // Set original_parameters struct for consis 72 // Set original_parameters struct for consistency 73 // 73 // 74 //SetOriginalParameters(rz); 74 //SetOriginalParameters(rz); 75 75 76 delete rz; 76 delete rz; 77 } 77 } 78 78 >> 79 // 79 // Create 80 // Create 80 // 81 // 81 // Generic create routine, called by each cons 82 // Generic create routine, called by each constructor after 82 // conversion of arguments 83 // conversion of arguments 83 // 84 // 84 void G4GenericPolycone::Create( G4double phiSt 85 void G4GenericPolycone::Create( G4double phiStart, 85 G4double phiTo << 86 G4double phiTotal, 86 G4ReduciblePol << 87 G4ReduciblePolygon *rz ) 87 { 88 { 88 // 89 // 89 // Perform checks of rz values 90 // Perform checks of rz values 90 // 91 // 91 if (rz->Amin() < 0.0) 92 if (rz->Amin() < 0.0) 92 { 93 { 93 std::ostringstream message; 94 std::ostringstream message; 94 message << "Illegal input parameters - " < 95 message << "Illegal input parameters - " << GetName() << G4endl 95 << " All R values must be > 96 << " All R values must be >= 0 !"; 96 G4Exception("G4GenericPolycone::Create()", 97 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002", 97 FatalErrorInArgument, message) 98 FatalErrorInArgument, message); 98 } 99 } 99 << 100 100 G4double rzArea = rz->Area(); 101 G4double rzArea = rz->Area(); 101 if (rzArea < -kCarTolerance) 102 if (rzArea < -kCarTolerance) 102 { << 103 rz->ReverseOrder(); 103 rz->ReverseOrder(); 104 } << 104 105 else if (rzArea < kCarTolerance) << 105 else if (rzArea < -kCarTolerance) 106 { 106 { 107 std::ostringstream message; 107 std::ostringstream message; 108 message << "Illegal input parameters - " < 108 message << "Illegal input parameters - " << GetName() << G4endl 109 << " R/Z cross section is z 109 << " R/Z cross section is zero or near zero: " << rzArea; 110 G4Exception("G4GenericPolycone::Create()", 110 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002", 111 FatalErrorInArgument, message) 111 FatalErrorInArgument, message); 112 } 112 } 113 << 113 114 if ( (!rz->RemoveDuplicateVertices( kCarTole 114 if ( (!rz->RemoveDuplicateVertices( kCarTolerance )) 115 || (!rz->RemoveRedundantVertices( kCarTole << 115 || (!rz->RemoveRedundantVertices( kCarTolerance )) ) 116 { 116 { 117 std::ostringstream message; 117 std::ostringstream message; 118 message << "Illegal input parameters - " < 118 message << "Illegal input parameters - " << GetName() << G4endl 119 << " Too few unique R/Z val 119 << " Too few unique R/Z values !"; 120 G4Exception("G4GenericPolycone::Create()", 120 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002", 121 FatalErrorInArgument, message) 121 FatalErrorInArgument, message); 122 } 122 } 123 123 124 if (rz->CrossesItself(1/kInfinity)) << 124 if (rz->CrossesItself(1/kInfinity)) 125 { 125 { 126 std::ostringstream message; 126 std::ostringstream message; 127 message << "Illegal input parameters - " < 127 message << "Illegal input parameters - " << GetName() << G4endl 128 << " R/Z segments cross !"; 128 << " R/Z segments cross !"; 129 G4Exception("G4GenericPolycone::Create()", 129 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002", 130 FatalErrorInArgument, message) 130 FatalErrorInArgument, message); 131 } 131 } 132 132 133 numCorner = rz->NumVertices(); 133 numCorner = rz->NumVertices(); 134 134 135 // 135 // 136 // Phi opening? Account for some possible ro 136 // Phi opening? Account for some possible roundoff, and interpret 137 // nonsense value as representing no phi ope 137 // nonsense value as representing no phi opening 138 // 138 // 139 if (phiTotal <= 0 || phiTotal > twopi-1E-10) 139 if (phiTotal <= 0 || phiTotal > twopi-1E-10) 140 { 140 { 141 phiIsOpen = false; 141 phiIsOpen = false; 142 startPhi = 0; 142 startPhi = 0; 143 endPhi = twopi; 143 endPhi = twopi; 144 } 144 } 145 else 145 else 146 { 146 { 147 phiIsOpen = true; 147 phiIsOpen = true; 148 << 148 149 // 149 // 150 // Convert phi into our convention 150 // Convert phi into our convention 151 // 151 // 152 startPhi = phiStart; 152 startPhi = phiStart; 153 while( startPhi < 0 ) // Loop checking, << 153 while( startPhi < 0 ) startPhi += twopi; 154 startPhi += twopi; << 154 155 << 156 endPhi = phiStart+phiTotal; 155 endPhi = phiStart+phiTotal; 157 while( endPhi < startPhi ) // Loop chec << 156 while( endPhi < startPhi ) endPhi += twopi; 158 endPhi += twopi; << 159 } 157 } 160 << 158 161 // 159 // 162 // Allocate corner array. << 160 // Allocate corner array. 163 // 161 // 164 corners = new G4PolyconeSideRZ[numCorner]; 162 corners = new G4PolyconeSideRZ[numCorner]; 165 163 166 // 164 // 167 // Copy corners 165 // Copy corners 168 // 166 // 169 G4ReduciblePolygonIterator iterRZ(rz); 167 G4ReduciblePolygonIterator iterRZ(rz); 170 << 168 171 G4PolyconeSideRZ* next = corners; << 169 G4PolyconeSideRZ *next = corners; 172 iterRZ.Begin(); 170 iterRZ.Begin(); 173 do // Loop checking, 13.08.2015, G.Cosmo << 171 do 174 { 172 { 175 next->r = iterRZ.GetA(); 173 next->r = iterRZ.GetA(); 176 next->z = iterRZ.GetB(); 174 next->z = iterRZ.GetB(); 177 } while( ++next, iterRZ.Next() ); 175 } while( ++next, iterRZ.Next() ); 178 << 176 179 // 177 // 180 // Allocate face pointer array 178 // Allocate face pointer array 181 // 179 // 182 numFace = phiIsOpen ? numCorner+2 : numCorne 180 numFace = phiIsOpen ? numCorner+2 : numCorner; 183 faces = new G4VCSGface*[numFace]; 181 faces = new G4VCSGface*[numFace]; 184 << 182 185 // 183 // 186 // Construct conical faces 184 // Construct conical faces 187 // 185 // 188 // But! Don't construct a face if both point 186 // But! Don't construct a face if both points are at zero radius! 189 // 187 // 190 G4PolyconeSideRZ *corner = corners, 188 G4PolyconeSideRZ *corner = corners, 191 *prev = corners + numCorner 189 *prev = corners + numCorner-1, 192 *nextNext; 190 *nextNext; 193 G4VCSGface** face = faces; << 191 G4VCSGface **face = faces; 194 do // Loop checking, 13.08.2015, G.Cosmo << 192 do 195 { 193 { 196 next = corner+1; 194 next = corner+1; 197 if (next >= corners+numCorner) next = corn 195 if (next >= corners+numCorner) next = corners; 198 nextNext = next+1; 196 nextNext = next+1; 199 if (nextNext >= corners+numCorner) nextNex 197 if (nextNext >= corners+numCorner) nextNext = corners; 200 << 198 201 if (corner->r < 1/kInfinity && next->r < 1 199 if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue; 202 << 200 203 // 201 // 204 // We must decide here if we can dare decl 202 // We must decide here if we can dare declare one of our faces 205 // as having a "valid" normal (i.e. allBeh 203 // as having a "valid" normal (i.e. allBehind = true). This 206 // is never possible if the face faces "in 204 // is never possible if the face faces "inward" in r. 207 // 205 // 208 G4bool allBehind; 206 G4bool allBehind; 209 if (corner->z > next->z) 207 if (corner->z > next->z) 210 { 208 { 211 allBehind = false; 209 allBehind = false; 212 } 210 } 213 else 211 else 214 { 212 { 215 // 213 // 216 // Otherwise, it is only true if the lin 214 // Otherwise, it is only true if the line passing 217 // through the two points of the segment 215 // through the two points of the segment do not 218 // split the r/z cross section 216 // split the r/z cross section 219 // 217 // 220 allBehind = !rz->BisectedBy( corner->r, 218 allBehind = !rz->BisectedBy( corner->r, corner->z, 221 next->r, next->z, kCarToleran 219 next->r, next->z, kCarTolerance ); 222 } 220 } 223 << 221 224 *face++ = new G4PolyconeSide( prev, corner 222 *face++ = new G4PolyconeSide( prev, corner, next, nextNext, 225 startPhi, endPhi-startPhi, phi 223 startPhi, endPhi-startPhi, phiIsOpen, allBehind ); 226 } while( prev=corner, corner=next, corner > 224 } while( prev=corner, corner=next, corner > corners ); 227 << 225 228 if (phiIsOpen) 226 if (phiIsOpen) 229 { 227 { 230 // 228 // 231 // Construct phi open edges 229 // Construct phi open edges 232 // 230 // 233 *face++ = new G4PolyPhiFace( rz, startPhi, 231 *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi ); 234 *face++ = new G4PolyPhiFace( rz, endPhi, 232 *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi ); 235 } 233 } 236 << 234 237 // 235 // 238 // We might have dropped a face or two: reca 236 // We might have dropped a face or two: recalculate numFace 239 // 237 // 240 numFace = (G4int)(face-faces); << 238 numFace = face-faces; 241 << 239 242 // 240 // 243 // Make enclosingCylinder 241 // Make enclosingCylinder 244 // 242 // 245 enclosingCylinder = 243 enclosingCylinder = 246 new G4EnclosingCylinder( rz, phiIsOpen, ph 244 new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal ); 247 } 245 } 248 246 >> 247 >> 248 // 249 // Fake default constructor - sets only member 249 // Fake default constructor - sets only member data and allocates memory 250 // for usage restri 250 // for usage restricted to object persistency. 251 // 251 // 252 G4GenericPolycone::G4GenericPolycone( __void__ 252 G4GenericPolycone::G4GenericPolycone( __void__& a ) 253 : G4VCSGfaceted(a), startPhi(0.), endPhi(0.) << 253 : G4VCSGfaceted(a), startPhi(0.), endPhi(0.), phiIsOpen(false), >> 254 numCorner(0), corners(0), enclosingCylinder(0) 254 { 255 { 255 } 256 } 256 257 >> 258 >> 259 // 257 // Destructor 260 // Destructor 258 // 261 // 259 G4GenericPolycone::~G4GenericPolycone() 262 G4GenericPolycone::~G4GenericPolycone() 260 { 263 { 261 delete [] corners; 264 delete [] corners; 262 delete enclosingCylinder; 265 delete enclosingCylinder; 263 delete fElements; << 264 delete fpPolyhedron; << 265 corners = nullptr; << 266 enclosingCylinder = nullptr; << 267 fElements = nullptr; << 268 fpPolyhedron = nullptr; << 269 } 266 } 270 267 >> 268 >> 269 // 271 // Copy constructor 270 // Copy constructor 272 // 271 // 273 G4GenericPolycone::G4GenericPolycone( const G4 << 272 G4GenericPolycone::G4GenericPolycone( const G4GenericPolycone &source ) 274 : G4VCSGfaceted( source ) 273 : G4VCSGfaceted( source ) 275 { 274 { 276 CopyStuff( source ); 275 CopyStuff( source ); 277 } 276 } 278 277 >> 278 >> 279 // 279 // Assignment operator 280 // Assignment operator 280 // 281 // 281 G4GenericPolycone& << 282 const G4GenericPolycone &G4GenericPolycone::operator=( const G4GenericPolycone &source ) 282 G4GenericPolycone::operator=( const G4GenericP << 283 { 283 { 284 if (this == &source) return *this; 284 if (this == &source) return *this; 285 << 285 286 G4VCSGfaceted::operator=( source ); 286 G4VCSGfaceted::operator=( source ); 287 << 287 288 delete [] corners; 288 delete [] corners; 289 // if (original_parameters) delete original_ 289 // if (original_parameters) delete original_parameters; 290 << 290 291 delete enclosingCylinder; 291 delete enclosingCylinder; 292 << 292 293 CopyStuff( source ); 293 CopyStuff( source ); 294 << 294 295 return *this; 295 return *this; 296 } 296 } 297 297 >> 298 >> 299 // 298 // CopyStuff 300 // CopyStuff 299 // 301 // 300 void G4GenericPolycone::CopyStuff( const G4Gen << 302 void G4GenericPolycone::CopyStuff( const G4GenericPolycone &source ) 301 { 303 { 302 // 304 // 303 // Simple stuff 305 // Simple stuff 304 // 306 // 305 startPhi = source.startPhi; 307 startPhi = source.startPhi; 306 endPhi = source.endPhi; 308 endPhi = source.endPhi; 307 phiIsOpen = source.phiIsOpen; << 309 phiIsOpen = source.phiIsOpen; 308 numCorner = source.numCorner; << 310 numCorner = source.numCorner; 309 311 310 // 312 // 311 // The corner array 313 // The corner array 312 // 314 // 313 corners = new G4PolyconeSideRZ[numCorner]; 315 corners = new G4PolyconeSideRZ[numCorner]; 314 << 316 315 G4PolyconeSideRZ *corn = corners, 317 G4PolyconeSideRZ *corn = corners, 316 *sourceCorn = source.corners; 318 *sourceCorn = source.corners; 317 do // Loop checking, 13.08.2015, G.Cosmo << 319 do 318 { 320 { 319 *corn = *sourceCorn; 321 *corn = *sourceCorn; 320 } while( ++sourceCorn, ++corn < corners+numC 322 } while( ++sourceCorn, ++corn < corners+numCorner ); 321 << 323 322 // 324 // 323 // Enclosing cylinder 325 // Enclosing cylinder 324 // 326 // 325 enclosingCylinder = new G4EnclosingCylinder( 327 enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder ); 326 << 327 // << 328 // Surface elements << 329 // << 330 delete fElements; << 331 fElements = nullptr; << 332 << 333 // Polyhedron << 334 // << 335 fRebuildPolyhedron = false; << 336 delete fpPolyhedron; << 337 fpPolyhedron = nullptr; << 338 } 328 } 339 329 >> 330 >> 331 // 340 // Reset 332 // Reset 341 // 333 // 342 G4bool G4GenericPolycone::Reset() 334 G4bool G4GenericPolycone::Reset() 343 { 335 { 344 std::ostringstream message; << 336 345 message << "Solid " << GetName() << " built << 337 std::ostringstream message; 346 << G4endl << "Not applicable to the << 338 message << "Solid " << GetName() << " built using generic construct." 347 G4Exception("G4GenericPolycone::Reset()", "G << 339 << G4endl << "Not applicable to the generic construct !"; 348 JustWarning, message, "Parameter << 340 G4Exception("G4GenericPolycone::Reset()", "GeomSolids1001", 349 return true; << 341 JustWarning, message, "Parameters NOT resetted."); >> 342 return 1; >> 343 350 } 344 } 351 345 >> 346 >> 347 // 352 // Inside 348 // Inside 353 // 349 // 354 // This is an override of G4VCSGfaceted::Insid 350 // This is an override of G4VCSGfaceted::Inside, created in order 355 // to speed things up by first checking with G 351 // to speed things up by first checking with G4EnclosingCylinder. 356 // 352 // 357 EInside G4GenericPolycone::Inside( const G4Thr << 353 EInside G4GenericPolycone::Inside( const G4ThreeVector &p ) const 358 { 354 { 359 // 355 // 360 // Quick test 356 // Quick test 361 // 357 // 362 if (enclosingCylinder->MustBeOutside(p)) ret 358 if (enclosingCylinder->MustBeOutside(p)) return kOutside; 363 359 364 // 360 // 365 // Long answer 361 // Long answer 366 // 362 // 367 return G4VCSGfaceted::Inside(p); 363 return G4VCSGfaceted::Inside(p); 368 } 364 } 369 365 >> 366 >> 367 // 370 // DistanceToIn 368 // DistanceToIn 371 // 369 // 372 // This is an override of G4VCSGfaceted::Insid 370 // This is an override of G4VCSGfaceted::Inside, created in order 373 // to speed things up by first checking with G 371 // to speed things up by first checking with G4EnclosingCylinder. 374 // 372 // 375 G4double G4GenericPolycone::DistanceToIn( cons << 373 G4double G4GenericPolycone::DistanceToIn( const G4ThreeVector &p, 376 cons << 374 const G4ThreeVector &v ) const 377 { 375 { 378 // 376 // 379 // Quick test 377 // Quick test 380 // 378 // 381 if (enclosingCylinder->ShouldMiss(p,v)) 379 if (enclosingCylinder->ShouldMiss(p,v)) 382 return kInfinity; 380 return kInfinity; 383 << 381 384 // 382 // 385 // Long answer 383 // Long answer 386 // 384 // 387 return G4VCSGfaceted::DistanceToIn( p, v ); 385 return G4VCSGfaceted::DistanceToIn( p, v ); 388 } 386 } 389 387 >> 388 >> 389 // 390 // DistanceToIn 390 // DistanceToIn 391 // 391 // 392 G4double G4GenericPolycone::DistanceToIn( cons << 392 G4double G4GenericPolycone::DistanceToIn( const G4ThreeVector &p ) const 393 { 393 { 394 return G4VCSGfaceted::DistanceToIn(p); 394 return G4VCSGfaceted::DistanceToIn(p); 395 } 395 } 396 396 397 // BoundingLimits << 397 398 // 398 // 399 // Get bounding box << 399 // ComputeDimensions 400 // 400 // 401 void << 401 /*void G4GenericPolycone::ComputeDimensions( G4VPVParameterisation* p, 402 G4GenericPolycone::BoundingLimits(G4ThreeVecto << 402 const G4int n, 403 G4ThreeVecto << 403 const G4VPhysicalVolume* pRep ) 404 { 404 { 405 G4double rmin = kInfinity, rmax = -kInfinity << 405 p->ComputeDimensions(*this,n,pRep); 406 G4double zmin = kInfinity, zmax = -kInfinity << 407 << 408 for (G4int i=0; i<GetNumRZCorner(); ++i) << 409 { << 410 G4PolyconeSideRZ corner = GetCorner(i); << 411 if (corner.r < rmin) rmin = corner.r; << 412 if (corner.r > rmax) rmax = corner.r; << 413 if (corner.z < zmin) zmin = corner.z; << 414 if (corner.z > zmax) zmax = corner.z; << 415 } << 416 << 417 if (IsOpen()) << 418 { << 419 G4TwoVector vmin,vmax; << 420 G4GeomTools::DiskExtent(rmin,rmax, << 421 GetSinStartPhi(),G << 422 GetSinEndPhi(),Get << 423 vmin,vmax); << 424 pMin.set(vmin.x(),vmin.y(),zmin); << 425 pMax.set(vmax.x(),vmax.y(),zmax); << 426 } << 427 else << 428 { << 429 pMin.set(-rmax,-rmax, zmin); << 430 pMax.set( rmax, rmax, zmax); << 431 } << 432 << 433 // Check correctness of the bounding box << 434 // << 435 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 436 { << 437 std::ostringstream message; << 438 message << "Bad bounding box (min >= max) << 439 << GetName() << " !" << 440 << "\npMin = " << pMin << 441 << "\npMax = " << pMax; << 442 G4Exception("GenericG4Polycone::BoundingLi << 443 JustWarning, message); << 444 DumpInfo(); << 445 } << 446 } 406 } 447 << 407 */ 448 // CalculateExtent << 449 // << 450 // Calculate extent under transform and specif << 451 // 408 // 452 G4bool << 453 G4GenericPolycone::CalculateExtent(const EAxis << 454 const G4Vox << 455 const G4Aff << 456 G4dou << 457 { << 458 G4ThreeVector bmin, bmax; << 459 G4bool exist; << 460 << 461 // Check bounding box (bbox) << 462 // << 463 BoundingLimits(bmin,bmax); << 464 G4BoundingEnvelope bbox(bmin,bmax); << 465 #ifdef G4BBOX_EXTENT << 466 return bbox.CalculateExtent(pAxis,pVoxelLimi << 467 #endif << 468 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 469 { << 470 return exist = pMin < pMax; << 471 } << 472 << 473 // To find the extent, RZ contour of the pol << 474 // in triangles. The extent is calculated as << 475 // all sub-polycones formed by rotation of t << 476 // << 477 G4TwoVectorList contourRZ; << 478 G4TwoVectorList triangles; << 479 G4double eminlim = pVoxelLimit.GetMinExtent( << 480 G4double emaxlim = pVoxelLimit.GetMaxExtent( << 481 << 482 // get RZ contour, ensure anticlockwise orde << 483 for (G4int i=0; i<GetNumRZCorner(); ++i) << 484 { << 485 G4PolyconeSideRZ corner = GetCorner(i); << 486 contourRZ.emplace_back(corner.r,corner.z); << 487 } << 488 G4double area = G4GeomTools::PolygonArea(con << 489 if (area < 0.) std::reverse(contourRZ.begin( << 490 << 491 // triangulate RZ countour << 492 if (!G4GeomTools::TriangulatePolygon(contour << 493 { << 494 std::ostringstream message; << 495 message << "Triangulation of RZ contour ha << 496 << GetName() << " !" << 497 << "\nExtent has been calculated u << 498 G4Exception("G4GenericPolycone::CalculateE << 499 "GeomMgt1002", JustWarning, me << 500 return bbox.CalculateExtent(pAxis,pVoxelLi << 501 } << 502 << 503 // set trigonometric values << 504 const G4int NSTEPS = 24; // numbe << 505 G4double astep = twopi/NSTEPS; // max a << 506 << 507 G4double sphi = GetStartPhi(); << 508 G4double ephi = GetEndPhi(); << 509 G4double dphi = IsOpen() ? ephi-sphi : two << 510 G4int ksteps = (dphi <= astep) ? 1 : (G4i << 511 G4double ang = dphi/ksteps; << 512 << 513 G4double sinHalf = std::sin(0.5*ang); << 514 G4double cosHalf = std::cos(0.5*ang); << 515 G4double sinStep = 2.*sinHalf*cosHalf; << 516 G4double cosStep = 1. - 2.*sinHalf*sinHalf; << 517 << 518 G4double sinStart = GetSinStartPhi(); << 519 G4double cosStart = GetCosStartPhi(); << 520 G4double sinEnd = GetSinEndPhi(); << 521 G4double cosEnd = GetCosEndPhi(); << 522 << 523 // define vectors and arrays << 524 std::vector<const G4ThreeVectorList *> polyg << 525 polygons.resize(ksteps+2); << 526 G4ThreeVectorList pols[NSTEPS+2]; << 527 for (G4int k=0; k<ksteps+2; ++k) pols[k].res << 528 for (G4int k=0; k<ksteps+2; ++k) polygons[k] << 529 G4double r0[6],z0[6]; // contour with origin << 530 G4double r1[6]; // shifted radii of ex << 531 << 532 // main loop along triangles << 533 pMin = kInfinity; << 534 pMax =-kInfinity; << 535 G4int ntria = (G4int)triangles.size()/3; << 536 for (G4int i=0; i<ntria; ++i) << 537 { << 538 G4int i3 = i*3; << 539 for (G4int k=0; k<3; ++k) << 540 { << 541 G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3; << 542 G4int k2 = k*2; << 543 // set contour with original edges of tr << 544 r0[k2+0] = triangles[e0].x(); z0[k2+0] = << 545 r0[k2+1] = triangles[e1].x(); z0[k2+1] = << 546 // set shifted radii << 547 r1[k2+0] = r0[k2+0]; << 548 r1[k2+1] = r0[k2+1]; << 549 if (z0[k2+1] - z0[k2+0] <= 0) continue; << 550 r1[k2+0] /= cosHalf; << 551 r1[k2+1] /= cosHalf; << 552 } << 553 << 554 // rotate countour, set sequence of 6-side << 555 G4double sinCur = sinStart*cosHalf + cosSt << 556 G4double cosCur = cosStart*cosHalf - sinSt << 557 for (G4int j=0; j<6; ++j) << 558 { << 559 pols[0][j].set(r0[j]*cosStart,r0[j]*sinS << 560 } << 561 for (G4int k=1; k<ksteps+1; ++k) << 562 { << 563 for (G4int j=0; j<6; ++j) << 564 { << 565 pols[k][j].set(r1[j]*cosCur,r1[j]*sinC << 566 } << 567 G4double sinTmp = sinCur; << 568 sinCur = sinCur*cosStep + cosCur*sinStep << 569 cosCur = cosCur*cosStep - sinTmp*sinStep << 570 } << 571 for (G4int j=0; j<6; ++j) << 572 { << 573 pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j] << 574 } << 575 << 576 // set sub-envelope and adjust extent << 577 G4double emin,emax; << 578 G4BoundingEnvelope benv(polygons); << 579 if (!benv.CalculateExtent(pAxis,pVoxelLimi << 580 if (emin < pMin) pMin = emin; << 581 if (emax > pMax) pMax = emax; << 582 if (eminlim > pMin && emaxlim < pMax) retu << 583 } << 584 return (pMin < pMax); << 585 } << 586 << 587 // GetEntityType 409 // GetEntityType 588 // 410 // 589 G4GeometryType G4GenericPolycone::GetEntityTy 411 G4GeometryType G4GenericPolycone::GetEntityType() const 590 { 412 { 591 return {"G4GenericPolycone"}; << 413 return G4String("G4GenericPolycone"); 592 } 414 } 593 415 594 // Clone << 416 595 // 417 // 596 // Make a clone of the object 418 // Make a clone of the object 597 // 419 // 598 G4VSolid* G4GenericPolycone::Clone() const 420 G4VSolid* G4GenericPolycone::Clone() const 599 { 421 { 600 return new G4GenericPolycone(*this); 422 return new G4GenericPolycone(*this); 601 } 423 } 602 424 603 // StreamInfo << 604 // 425 // 605 // Stream object contents to an output stream 426 // Stream object contents to an output stream 606 // 427 // 607 std::ostream& G4GenericPolycone::StreamInfo( s 428 std::ostream& G4GenericPolycone::StreamInfo( std::ostream& os ) const 608 { 429 { 609 G4long oldprc = os.precision(16); << 430 G4int oldprc = os.precision(16); 610 os << "------------------------------------- 431 os << "-----------------------------------------------------------\n" 611 << " *** Dump for solid - " << GetName 432 << " *** Dump for solid - " << GetName() << " ***\n" 612 << " ================================= 433 << " ===================================================\n" 613 << " Solid type: G4GenericPolycone\n" 434 << " Solid type: G4GenericPolycone\n" 614 << " Parameters: \n" 435 << " Parameters: \n" 615 << " starting phi angle : " << startPh 436 << " starting phi angle : " << startPhi/degree << " degrees \n" 616 << " ending phi angle : " << endPhi/ 437 << " ending phi angle : " << endPhi/degree << " degrees \n"; 617 G4int i=0; 438 G4int i=0; 618 << 439 619 os << " number of RZ points: " << numCorn 440 os << " number of RZ points: " << numCorner << "\n" 620 << " RZ values (corners): \n 441 << " RZ values (corners): \n"; 621 for (i=0; i<numCorner; i++) 442 for (i=0; i<numCorner; i++) 622 { 443 { 623 os << " " 444 os << " " 624 << corners[i].r << ", " << corners[i 445 << corners[i].r << ", " << corners[i].z << "\n"; 625 } 446 } 626 os << "------------------------------------- 447 os << "-----------------------------------------------------------\n"; 627 os.precision(oldprc); 448 os.precision(oldprc); 628 449 629 return os; 450 return os; 630 } 451 } 631 452 632 ////////////////////////////////////////////// << 633 // << 634 // Return volume << 635 453 636 G4double G4GenericPolycone::GetCubicVolume() << 454 >> 455 // >> 456 // GetPointOnSurface >> 457 // >> 458 G4ThreeVector G4GenericPolycone::GetPointOnSurface() const 637 { 459 { 638 if (fCubicVolume == 0.) << 460 return GetPointOnSurfaceGeneric(); 639 { << 461 640 G4double total = 0.; << 641 G4int nrz = GetNumRZCorner(); << 642 G4PolyconeSideRZ a = GetCorner(nrz - 1); << 643 for (G4int i=0; i<nrz; ++i) << 644 { << 645 G4PolyconeSideRZ b = GetCorner(i); << 646 total += (b.r*b.r + b.r*a.r + a.r*a.r)*( << 647 a = b; << 648 } << 649 fCubicVolume = std::abs(total)*(GetEndPhi( << 650 } << 651 return fCubicVolume; << 652 } 462 } 653 463 654 ////////////////////////////////////////////// << 655 // 464 // 656 // Return surface area << 465 // CreatePolyhedron 657 << 466 // 658 G4double G4GenericPolycone::GetSurfaceArea() << 467 G4Polyhedron* G4GenericPolycone::CreatePolyhedron() const 659 { << 468 { 660 if (fSurfaceArea == 0.) << 469 // The following code prepares for: 661 { << 470 // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces, 662 // phi cut area << 471 // const double xyz[][3], 663 G4int nrz = GetNumRZCorner(); << 472 // const int faces_vec[][4]) 664 G4double scut = 0.; << 473 // Here is an extract from the header file HepPolyhedron.h: 665 if (IsOpen()) << 474 /** >> 475 * Creates user defined polyhedron. >> 476 * This function allows to the user to define arbitrary polyhedron. >> 477 * The faces of the polyhedron should be either triangles or planar >> 478 * quadrilateral. Nodes of a face are defined by indexes pointing to >> 479 * the elements in the xyz array. Numeration of the elements in the >> 480 * array starts from 1 (like in fortran). The indexes can be positive >> 481 * or negative. Negative sign means that the corresponding edge is >> 482 * invisible. The normal of the face should be directed to exterior >> 483 * of the polyhedron. >> 484 * >> 485 * @param Nnodes number of nodes >> 486 * @param Nfaces number of faces >> 487 * @param xyz nodes >> 488 * @param faces_vec faces (quadrilaterals or triangles) >> 489 * @return status of the operation - is non-zero in case of problem >> 490 */ >> 491 const G4int numSide = >> 492 G4int(G4Polyhedron::GetNumberOfRotationSteps() >> 493 * (endPhi - startPhi) / twopi) + 1; >> 494 G4int nNodes; >> 495 G4int nFaces; >> 496 typedef G4double double3[3]; >> 497 double3* xyz; >> 498 typedef G4int int4[4]; >> 499 int4* faces_vec; >> 500 if (phiIsOpen) 666 { 501 { 667 G4PolyconeSideRZ a = GetCorner(nrz - 1); << 502 // Triangulate open ends. Simple ear-chopping algorithm... 668 for (G4int i=0; i<nrz; ++i) << 503 // I'm not sure how robust this algorithm is (J.Allison). >> 504 // >> 505 std::vector<G4bool> chopped(numCorner, false); >> 506 std::vector<G4int*> triQuads; >> 507 G4int remaining = numCorner; >> 508 G4int iStarter = 0; >> 509 while (remaining >= 3) 669 { 510 { 670 G4PolyconeSideRZ b = GetCorner(i); << 511 // Find unchopped corners... 671 scut += a.r*b.z - a.z*b.r; << 512 // 672 a = b; << 513 G4int A = -1, B = -1, C = -1; >> 514 G4int iStepper = iStarter; >> 515 do >> 516 { >> 517 if (A < 0) { A = iStepper; } >> 518 else if (B < 0) { B = iStepper; } >> 519 else if (C < 0) { C = iStepper; } >> 520 do >> 521 { >> 522 if (++iStepper >= numCorner) { iStepper = 0; } >> 523 } >> 524 while (chopped[iStepper]); >> 525 } >> 526 while (C < 0 && iStepper != iStarter); >> 527 >> 528 // Check triangle at B is pointing outward (an "ear"). >> 529 // Sign of z cross product determines... >> 530 // >> 531 G4double BAr = corners[A].r - corners[B].r; >> 532 G4double BAz = corners[A].z - corners[B].z; >> 533 G4double BCr = corners[C].r - corners[B].r; >> 534 G4double BCz = corners[C].z - corners[B].z; >> 535 if (BAr * BCz - BAz * BCr < kCarTolerance) >> 536 { >> 537 G4int* tq = new G4int[3]; >> 538 tq[0] = A + 1; >> 539 tq[1] = B + 1; >> 540 tq[2] = C + 1; >> 541 triQuads.push_back(tq); >> 542 chopped[B] = true; >> 543 --remaining; >> 544 } >> 545 else >> 546 { >> 547 do >> 548 { >> 549 if (++iStarter >= numCorner) { iStarter = 0; } >> 550 } >> 551 while (chopped[iStarter]); >> 552 } >> 553 } >> 554 // Transfer to faces... >> 555 // >> 556 nNodes = (numSide + 1) * numCorner; >> 557 nFaces = numSide * numCorner + 2 * triQuads.size(); >> 558 faces_vec = new int4[nFaces]; >> 559 G4int iface = 0; >> 560 G4int addition = numCorner * numSide; >> 561 G4int d = numCorner - 1; >> 562 for (G4int iEnd = 0; iEnd < 2; ++iEnd) >> 563 { >> 564 for (size_t i = 0; i < triQuads.size(); ++i) >> 565 { >> 566 // Negative for soft/auxiliary/normally invisible edges... >> 567 // >> 568 G4int a, b, c; >> 569 if (iEnd == 0) >> 570 { >> 571 a = triQuads[i][0]; >> 572 b = triQuads[i][1]; >> 573 c = triQuads[i][2]; >> 574 } >> 575 else >> 576 { >> 577 a = triQuads[i][0] + addition; >> 578 b = triQuads[i][2] + addition; >> 579 c = triQuads[i][1] + addition; >> 580 } >> 581 G4int ab = std::abs(b - a); >> 582 G4int bc = std::abs(c - b); >> 583 G4int ca = std::abs(a - c); >> 584 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a; >> 585 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b; >> 586 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c; >> 587 faces_vec[iface][3] = 0; >> 588 ++iface; >> 589 } 673 } 590 } 674 scut = std::abs(scut); << 675 } << 676 // lateral surface area << 677 G4double slat = 0; << 678 G4PolyconeSideRZ a = GetCorner(nrz - 1); << 679 for (G4int i=0; i<nrz; ++i) << 680 { << 681 G4PolyconeSideRZ b = GetCorner(i); << 682 G4double h = std::sqrt((b.r - a.r)*(b.r << 683 slat += (b.r + a.r)*h; << 684 a = b; << 685 } << 686 slat *= (GetEndPhi() - GetStartPhi())/2.; << 687 fSurfaceArea = scut + slat; << 688 } << 689 return fSurfaceArea; << 690 } << 691 591 692 ////////////////////////////////////////////// << 592 // Continue with sides... 693 // << 694 // Set vector of surface elements, auxiliary m << 695 // random points on surface << 696 593 697 void G4GenericPolycone::SetSurfaceElements() c << 594 xyz = new double3[nNodes]; 698 { << 595 const G4double dPhi = (endPhi - startPhi) / numSide; 699 fElements = new std::vector<G4GenericPolycon << 596 G4double phi = startPhi; 700 G4double sarea = 0.; << 597 G4int ixyz = 0; 701 G4int nrz = GetNumRZCorner(); << 598 for (G4int iSide = 0; iSide < numSide; ++iSide) >> 599 { >> 600 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner) >> 601 { >> 602 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi); >> 603 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi); >> 604 xyz[ixyz][2] = corners[iCorner].z; >> 605 if (iSide == 0) // startPhi >> 606 { >> 607 if (iCorner < numCorner - 1) >> 608 { >> 609 faces_vec[iface][0] = ixyz + 1; >> 610 faces_vec[iface][1] = -(ixyz + numCorner + 1); >> 611 faces_vec[iface][2] = ixyz + numCorner + 2; >> 612 faces_vec[iface][3] = ixyz + 2; >> 613 } >> 614 else >> 615 { >> 616 faces_vec[iface][0] = ixyz + 1; >> 617 faces_vec[iface][1] = -(ixyz + numCorner + 1); >> 618 faces_vec[iface][2] = ixyz + 2; >> 619 faces_vec[iface][3] = ixyz - numCorner + 2; >> 620 } >> 621 } >> 622 else if (iSide == numSide - 1) // endPhi >> 623 { >> 624 if (iCorner < numCorner - 1) >> 625 { >> 626 faces_vec[iface][0] = ixyz + 1; >> 627 faces_vec[iface][1] = ixyz + numCorner + 1; >> 628 faces_vec[iface][2] = ixyz + numCorner + 2; >> 629 faces_vec[iface][3] = -(ixyz + 2); >> 630 } >> 631 else >> 632 { >> 633 faces_vec[iface][0] = ixyz + 1; >> 634 faces_vec[iface][1] = ixyz + numCorner + 1; >> 635 faces_vec[iface][2] = ixyz + 2; >> 636 faces_vec[iface][3] = -(ixyz - numCorner + 2); >> 637 } >> 638 } >> 639 else >> 640 { >> 641 if (iCorner < numCorner - 1) >> 642 { >> 643 faces_vec[iface][0] = ixyz + 1; >> 644 faces_vec[iface][1] = -(ixyz + numCorner + 1); >> 645 faces_vec[iface][2] = ixyz + numCorner + 2; >> 646 faces_vec[iface][3] = -(ixyz + 2); >> 647 } >> 648 else >> 649 { >> 650 faces_vec[iface][0] = ixyz + 1; >> 651 faces_vec[iface][1] = -(ixyz + numCorner + 1); >> 652 faces_vec[iface][2] = ixyz + 2; >> 653 faces_vec[iface][3] = -(ixyz - numCorner + 2); >> 654 } >> 655 } >> 656 ++iface; >> 657 ++ixyz; >> 658 } >> 659 phi += dPhi; >> 660 } 702 661 703 // set lateral surface elements << 662 // Last corners... 704 G4double dphi = GetEndPhi() - GetStartPhi(); << 705 G4int ia = nrz - 1; << 706 for (G4int ib=0; ib<nrz; ++ib) << 707 { << 708 G4PolyconeSideRZ a = GetCorner(ia); << 709 G4PolyconeSideRZ b = GetCorner(ib); << 710 G4GenericPolycone::surface_element selem; << 711 selem.i0 = ia; << 712 selem.i1 = ib; << 713 selem.i2 = -1; << 714 ia = ib; << 715 if (a.r == 0. && b.r == 0.) continue; << 716 G4double h = std::sqrt((b.r - a.r)*(b.r - << 717 sarea += 0.5*dphi*(b.r + a.r)*h; << 718 selem.area = sarea; << 719 fElements->push_back(selem); << 720 } << 721 663 722 // set elements for phi cuts << 664 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner) 723 if (IsOpen()) << 665 { 724 { << 666 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi); 725 G4TwoVectorList contourRZ; << 667 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi); 726 std::vector<G4int> triangles; << 668 xyz[ixyz][2] = corners[iCorner].z; 727 for (G4int i=0; i<nrz; ++i) << 669 ++ixyz; 728 { << 670 } 729 G4PolyconeSideRZ corner = GetCorner(i); << 730 contourRZ.emplace_back(corner.r, corner. << 731 } 671 } 732 G4GeomTools::TriangulatePolygon(contourRZ, << 672 else // !phiIsOpen - i.e., a complete 360 degrees. 733 auto ntria = (G4int)triangles.size(); << 734 for (G4int i=0; i<ntria; i+=3) << 735 { 673 { 736 G4GenericPolycone::surface_element selem << 674 nNodes = numSide * numCorner; 737 selem.i0 = triangles[i]; << 675 nFaces = numSide * numCorner;; 738 selem.i1 = triangles[i+1]; << 676 xyz = new double3[nNodes]; 739 selem.i2 = triangles[i+2]; << 677 faces_vec = new int4[nFaces]; 740 G4PolyconeSideRZ a = GetCorner(selem.i0) << 678 const G4double dPhi = (endPhi - startPhi) / numSide; 741 G4PolyconeSideRZ b = GetCorner(selem.i1) << 679 G4double phi = startPhi; 742 G4PolyconeSideRZ c = GetCorner(selem.i2) << 680 G4int ixyz = 0, iface = 0; 743 G4double stria = << 681 for (G4int iSide = 0; iSide < numSide; ++iSide) 744 std::abs(G4GeomTools::TriangleArea(a.r << 682 { 745 sarea += stria; << 683 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner) 746 selem.area = sarea; << 684 { 747 fElements->push_back(selem); // start ph << 685 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi); 748 sarea += stria; << 686 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi); 749 selem.area = sarea; << 687 xyz[ixyz][2] = corners[iCorner].z; 750 selem.i0 += nrz; << 688 751 fElements->push_back(selem); // end phi << 689 if (iSide < numSide - 1) 752 } << 690 { 753 } << 691 if (iCorner < numCorner - 1) 754 } << 692 { 755 << 693 faces_vec[iface][0] = ixyz + 1; 756 ////////////////////////////////////////////// << 694 faces_vec[iface][1] = -(ixyz + numCorner + 1); 757 // << 695 faces_vec[iface][2] = ixyz + numCorner + 2; 758 // Generate random point on surface << 696 faces_vec[iface][3] = -(ixyz + 2); 759 << 697 } 760 G4ThreeVector G4GenericPolycone::GetPointOnSur << 698 else 761 { << 699 { 762 // Set surface elements << 700 faces_vec[iface][0] = ixyz + 1; 763 if (fElements == nullptr) << 701 faces_vec[iface][1] = -(ixyz + numCorner + 1); 764 { << 702 faces_vec[iface][2] = ixyz + 2; 765 G4AutoLock l(&surface_elementsMutex); << 703 faces_vec[iface][3] = -(ixyz - numCorner + 2); 766 SetSurfaceElements(); << 704 } 767 l.unlock(); << 705 } 768 } << 706 else // Last side joins ends... 769 << 707 { 770 // Select surface element << 708 if (iCorner < numCorner - 1) 771 G4GenericPolycone::surface_element selem; << 709 { 772 selem = fElements->back(); << 710 faces_vec[iface][0] = ixyz + 1; 773 G4double select = selem.area*G4QuickRand(); << 711 faces_vec[iface][1] = -(ixyz + numCorner - nFaces + 1); 774 auto it = std::lower_bound(fElements->begin( << 712 faces_vec[iface][2] = ixyz + numCorner - nFaces + 2; 775 [](const G4Generi << 713 faces_vec[iface][3] = -(ixyz + 2); 776 -> G4bool { retur << 714 } 777 << 715 else 778 // Generate random point << 716 { 779 G4double r = 0, z = 0, phi = 0; << 717 faces_vec[iface][0] = ixyz + 1; 780 G4double u = G4QuickRand(); << 718 faces_vec[iface][1] = -(ixyz - nFaces + numCorner + 1); 781 G4double v = G4QuickRand(); << 719 faces_vec[iface][2] = ixyz - nFaces + 2; 782 G4int i0 = (*it).i0; << 720 faces_vec[iface][3] = -(ixyz - numCorner + 2); 783 G4int i1 = (*it).i1; << 721 } 784 G4int i2 = (*it).i2; << 722 } 785 if (i2 < 0) // lateral surface << 723 ++ixyz; 786 { << 724 ++iface; 787 G4PolyconeSideRZ p0 = GetCorner(i0); << 725 } 788 G4PolyconeSideRZ p1 = GetCorner(i1); << 726 phi += dPhi; 789 if (p1.r < p0.r) << 727 } 790 { << 791 p0 = GetCorner(i1); << 792 p1 = GetCorner(i0); << 793 } 728 } 794 if (p1.r - p0.r < kCarTolerance) // cylind << 729 G4Polyhedron* polyhedron = new G4Polyhedron; 795 { << 730 G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec); 796 r = (p1.r - p0.r)*u + p0.r; << 731 delete [] faces_vec; 797 z = (p1.z - p0.z)*u + p0.z; << 732 delete [] xyz; >> 733 if (problem) >> 734 { >> 735 std::ostringstream message; >> 736 message << "Problem creating G4Polyhedron for: " << GetName(); >> 737 G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002", >> 738 JustWarning, message); >> 739 delete polyhedron; >> 740 return 0; 798 } 741 } 799 else // conical surface << 742 else 800 { 743 { 801 r = std::sqrt(p1.r*p1.r*u + p0.r*p0.r*(1 << 744 return polyhedron; 802 z = p0.z + (p1.z - p0.z)*(r - p0.r)/(p1. << 803 } 745 } 804 phi = (GetEndPhi() - GetStartPhi())*v + Ge << 805 } << 806 else // phi cut << 807 { << 808 G4int nrz = GetNumRZCorner(); << 809 phi = (i0 < nrz) ? GetStartPhi() : GetEndP << 810 if (i0 >= nrz) { i0 -= nrz; } << 811 G4PolyconeSideRZ p0 = GetCorner(i0); << 812 G4PolyconeSideRZ p1 = GetCorner(i1); << 813 G4PolyconeSideRZ p2 = GetCorner(i2); << 814 if (u + v > 1.) { u = 1. - u; v = 1. - v; << 815 r = (p1.r - p0.r)*u + (p2.r - p0.r)*v + p << 816 z = (p1.z - p0.z)*u + (p2.z - p0.z)*v + p << 817 } << 818 return {r*std::cos(phi), r*std::sin(phi), z} << 819 } << 820 << 821 ////////////////////////////////////////////// << 822 // << 823 // CreatePolyhedron << 824 << 825 G4Polyhedron* G4GenericPolycone::CreatePolyhed << 826 { << 827 std::vector<G4TwoVector> rz(numCorner); << 828 for (G4int i = 0; i < numCorner; ++i) << 829 rz[i].set(corners[i].r, corners[i].z); << 830 return new G4PolyhedronPcon(startPhi, endPhi << 831 } 746 } 832 747 833 #endif 748 #endif 834 749