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