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