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 26 27 // ------------------------------------------- 27 // -------------------------------------------------------------------- 28 // Implementation for FastAerosolSolid class 28 // Implementation for FastAerosolSolid class 29 // Author: A.Knaian (ara@nklabs.com), N.MacFad 29 // Author: A.Knaian (ara@nklabs.com), N.MacFadden (natemacfadden@gmail.com) 30 // ------------------------------------------- 30 // -------------------------------------------------------------------- 31 31 32 #include "FastAerosolSolid.hh" 32 #include "FastAerosolSolid.hh" 33 33 34 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh" 35 35 36 // calculate extent 36 // calculate extent 37 #include "G4BoundingEnvelope.hh" 37 #include "G4BoundingEnvelope.hh" 38 #include "G4AffineTransform.hh" 38 #include "G4AffineTransform.hh" 39 #include "G4VoxelLimits.hh" 39 #include "G4VoxelLimits.hh" 40 40 41 // visualization 41 // visualization 42 #include "G4VGraphicsScene.hh" 42 #include "G4VGraphicsScene.hh" 43 #include "G4VisExtent.hh" 43 #include "G4VisExtent.hh" 44 44 45 // polyhedron 45 // polyhedron 46 #include "G4AutoLock.hh" 46 #include "G4AutoLock.hh" 47 #include "G4Polyhedron.hh" 47 #include "G4Polyhedron.hh" 48 #include "HepPolyhedronProcessor.h" 48 #include "HepPolyhedronProcessor.h" 49 49 50 namespace 50 namespace 51 { 51 { 52 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER << 52 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER; 53 } 53 } 54 54 55 55 56 ////////////////////////////////////////////// 56 /////////////////////////////////////////////////////////////////////////////// 57 // 57 // 58 // Constructor 58 // Constructor 59 // 59 // 60 FastAerosolSolid::FastAerosolSolid(const G4Str << 60 FastAerosolSolid::FastAerosolSolid(const G4String& pName, >> 61 FastAerosol* pCloud, >> 62 G4VSolid* pDroplet, 61 std::function<G4RotationM 63 std::function<G4RotationMatrix (G4ThreeVector)> pRotation) 62 : G4VSolid(pName), fCloud(pCloud), fDroplet( << 64 : G4VSolid(pName), fCloud(pCloud), fDroplet(pDroplet), fRotation(pRotation), fRebuildPolyhedron(false), fpPolyhedron(0) 63 { 65 { 64 // Get cloud size from fCloud 66 // Get cloud size from fCloud 65 G4ThreeVector cloudPMin, cloudPMax; 67 G4ThreeVector cloudPMin, cloudPMax; 66 fCloud->GetBoundingLimits(cloudPMin, cloudPM 68 fCloud->GetBoundingLimits(cloudPMin, cloudPMax); 67 69 68 fVisDx = cloudPMax.x(); 70 fVisDx = cloudPMax.x(); 69 fVisDy = cloudPMax.y(); 71 fVisDy = cloudPMax.y(); 70 fVisDz = cloudPMax.z(); 72 fVisDz = cloudPMax.z(); 71 73 72 // Check and set droplet radius 74 // Check and set droplet radius 73 G4double pR = fCloud->GetRadius(); 75 G4double pR = fCloud->GetRadius(); 74 // would be nice to add a check to make sure 76 // would be nice to add a check to make sure pDroplet fits in sphere of radius pR 75 fR = pR; 77 fR = pR; 76 78 77 fBulk = fCloud->GetBulk(); 79 fBulk = fCloud->GetBulk(); 78 80 79 farFromCloudDist = fCloud->GetPreSphereR()*f 81 farFromCloudDist = fCloud->GetPreSphereR()*fR; 80 } 82 } 81 83 82 84 83 ////////////////////////////////////////////// 85 /////////////////////////////////////////////////////////////////////////////// 84 // 86 // 85 // Alternative constructor (constant rotation 87 // Alternative constructor (constant rotation function) 86 // 88 // 87 FastAerosolSolid::FastAerosolSolid(const G4Str 89 FastAerosolSolid::FastAerosolSolid(const G4String& pName, 88 FastAerosol* pCloud, 90 FastAerosol* pCloud, 89 G4VSolid* pDroplet): 91 G4VSolid* pDroplet): 90 FastAerosolSolid(pName, pCloud, pDroplet, 92 FastAerosolSolid(pName, pCloud, pDroplet, 91 [](G4ThreeVector) {return G4Rotatio 93 [](G4ThreeVector) {return G4RotationMatrix();}) 92 {} 94 {} 93 95 94 ////////////////////////////////////////////// 96 /////////////////////////////////////////////////////////////////////////////// 95 // 97 // 96 // Fake default constructor - sets only member 98 // Fake default constructor - sets only member data and allocates memory 97 // for usage restri 99 // for usage restricted to object persistency. 98 // 100 // 99 FastAerosolSolid::FastAerosolSolid( __void__& 101 FastAerosolSolid::FastAerosolSolid( __void__& a ) 100 : G4VSolid(a), fCloud(nullptr), fDroplet(nul 102 : G4VSolid(a), fCloud(nullptr), fDroplet(nullptr), 101 fBulk(nullptr), fR(0.), 103 fBulk(nullptr), fR(0.), 102 fVisDx(0.), fVisDy(0.), fVisDz(0.), 104 fVisDx(0.), fVisDy(0.), fVisDz(0.), 103 fCubicVolume(0.), fSurfaceArea(0.), 105 fCubicVolume(0.), fSurfaceArea(0.), 104 farFromCloudDist(0.), 106 farFromCloudDist(0.), 105 fRotation([](G4ThreeVector) {return G4Rotati 107 fRotation([](G4ThreeVector) {return G4RotationMatrix();}), 106 fRebuildPolyhedron(false), fpPolyhedron(null << 108 fRebuildPolyhedron(false), fpPolyhedron(0) 107 { 109 { 108 } 110 } 109 111 110 ////////////////////////////////////////////// 112 /////////////////////////////////////////////////////////////////////////////// 111 // 113 // >> 114 // Destructor >> 115 // >> 116 FastAerosolSolid::~FastAerosolSolid() { >> 117 } >> 118 >> 119 /////////////////////////////////////////////////////////////////////////////// >> 120 // 112 // Copy constructor 121 // Copy constructor 113 // 122 // 114 FastAerosolSolid::FastAerosolSolid(const FastA 123 FastAerosolSolid::FastAerosolSolid(const FastAerosolSolid &rhs) 115 : G4VSolid(rhs), fCloud(rhs.fCloud), fDrople 124 : G4VSolid(rhs), fCloud(rhs.fCloud), fDroplet(rhs.fDroplet), 116 fBulk(rhs.fBulk), fR(rhs.fR), 125 fBulk(rhs.fBulk), fR(rhs.fR), 117 fVisDx(rhs.fVisDx), fVisDy(rhs.fVisDy), fVis 126 fVisDx(rhs.fVisDx), fVisDy(rhs.fVisDy), fVisDz(rhs.fVisDz), 118 fCubicVolume(rhs.fCubicVolume), fSurfaceArea 127 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea), 119 farFromCloudDist(rhs.farFromCloudDist), 128 farFromCloudDist(rhs.farFromCloudDist), 120 fRotation(rhs.fRotation), 129 fRotation(rhs.fRotation), 121 fRebuildPolyhedron(rhs.fRebuildPolyhedron), 130 fRebuildPolyhedron(rhs.fRebuildPolyhedron), fpPolyhedron(rhs.fpPolyhedron) 122 { 131 { 123 } 132 } 124 133 125 ////////////////////////////////////////////// 134 ////////////////////////////////////////////////////////////////////////// 126 // 135 // 127 // Assignment operator 136 // Assignment operator 128 // 137 // 129 FastAerosolSolid &FastAerosolSolid::operator = 138 FastAerosolSolid &FastAerosolSolid::operator = (const FastAerosolSolid &rhs) 130 { 139 { 131 // Check assignment to self 140 // Check assignment to self 132 // 141 // 133 if (this == &rhs) 142 if (this == &rhs) 134 { 143 { 135 return *this; 144 return *this; 136 } 145 } 137 146 138 // Copy base class data 147 // Copy base class data 139 // 148 // 140 G4VSolid::operator=(rhs); 149 G4VSolid::operator=(rhs); 141 150 142 // Copy data 151 // Copy data 143 // 152 // 144 fCloud = rhs.fCloud; 153 fCloud = rhs.fCloud; 145 fDroplet = rhs.fDroplet; 154 fDroplet = rhs.fDroplet; 146 fBulk = rhs.fBulk; 155 fBulk = rhs.fBulk; 147 fR = rhs.fR; 156 fR = rhs.fR; 148 fVisDx = rhs.fVisDx; 157 fVisDx = rhs.fVisDx; 149 fVisDy = rhs.fVisDy; 158 fVisDy = rhs.fVisDy; 150 fVisDz = rhs.fVisDz; 159 fVisDz = rhs.fVisDz; 151 fCubicVolume = rhs.fCubicVolume; 160 fCubicVolume = rhs.fCubicVolume; 152 fSurfaceArea = rhs.fSurfaceArea; 161 fSurfaceArea = rhs.fSurfaceArea; 153 farFromCloudDist = rhs.farFromCloudDist; 162 farFromCloudDist = rhs.farFromCloudDist; 154 fRotation = rhs.fRotation; 163 fRotation = rhs.fRotation; 155 fRebuildPolyhedron = rhs.fRebuildPolyhedron; 164 fRebuildPolyhedron = rhs.fRebuildPolyhedron; 156 fpPolyhedron = rhs.fpPolyhedron; 165 fpPolyhedron = rhs.fpPolyhedron; 157 166 158 167 159 return *this; 168 return *this; 160 } 169 } 161 170 162 ////////////////////////////////////////////// 171 /////////////////////////////////////////////////////////////////////////////// 163 // 172 // 164 // Calculate extent under transform and specif 173 // Calculate extent under transform and specified limit 165 // 174 // 166 G4bool FastAerosolSolid::CalculateExtent(const 175 G4bool FastAerosolSolid::CalculateExtent(const EAxis pAxis, 167 const G4VoxelLimits &pVo 176 const G4VoxelLimits &pVoxelLimit, 168 const G4AffineTransform 177 const G4AffineTransform &pTransform, 169 G4double &pMin, G4double 178 G4double &pMin, G4double &pMax) const 170 { 179 { 171 // Get smallest box to fully contain the clo 180 // Get smallest box to fully contain the cloud of objects, not just the centers 172 // 181 // 173 G4ThreeVector bmin, bmax; 182 G4ThreeVector bmin, bmax; 174 fCloud->GetBoundingLimits(bmin, bmax); 183 fCloud->GetBoundingLimits(bmin, bmax); 175 184 176 // Find extent 185 // Find extent 177 // 186 // 178 G4BoundingEnvelope bbox(bmin, bmax); 187 G4BoundingEnvelope bbox(bmin, bmax); 179 return bbox.CalculateExtent(pAxis, pVoxelLim 188 return bbox.CalculateExtent(pAxis, pVoxelLimit, pTransform, pMin, pMax); 180 } 189 } 181 190 182 ////////////////////////////////////////////// 191 /////////////////////////////////////////////////////////////////////////////// 183 // 192 // 184 // Return whether point inside/outside/on surf 193 // Return whether point inside/outside/on surface 185 // 194 // 186 // This function assumes the cloud has at leas 195 // This function assumes the cloud has at least 1 droplet 187 // 196 // 188 EInside FastAerosolSolid::Inside(const G4Three 197 EInside FastAerosolSolid::Inside(const G4ThreeVector &p) const 189 { 198 { 190 G4ThreeVector center; 199 G4ThreeVector center; 191 G4double closestDistance; 200 G4double closestDistance; 192 201 193 fCloud->GetNearestDroplet(p, center, closest 202 fCloud->GetNearestDroplet(p, center, closestDistance, fR, fDroplet, fRotation); 194 203 195 if (closestDistance==0.0) 204 if (closestDistance==0.0) 196 { 205 { 197 G4RotationMatrix irotm = fRotation(center) 206 G4RotationMatrix irotm = fRotation(center).inverse(); 198 207 199 return fDroplet->Inside( irotm*(p - center 208 return fDroplet->Inside( irotm*(p - center) ); 200 } 209 } 201 else 210 else 202 { 211 { 203 return kOutside; 212 return kOutside; 204 } 213 } 205 } 214 } 206 215 207 ////////////////////////////////////////////// 216 ///////////////////////////////////////////////////////////////////// 208 // 217 // 209 // Return unit normal of surface closest to p 218 // Return unit normal of surface closest to p 210 // 219 // 211 // This function assumes the cloud has at leas 220 // This function assumes the cloud has at least 1 droplet 212 // 221 // 213 G4ThreeVector FastAerosolSolid::SurfaceNormal( 222 G4ThreeVector FastAerosolSolid::SurfaceNormal(const G4ThreeVector &p) const 214 { 223 { 215 G4ThreeVector center; 224 G4ThreeVector center; 216 G4double closestDistance; 225 G4double closestDistance; 217 226 218 fCloud->GetNearestDroplet(p, center, closest 227 fCloud->GetNearestDroplet(p, center, closestDistance, DBL_MAX, fDroplet, fRotation); 219 228 220 G4RotationMatrix rotm = fRotation(center); 229 G4RotationMatrix rotm = fRotation(center); 221 230 222 return rotm*( fDroplet->SurfaceNormal( rotm. 231 return rotm*( fDroplet->SurfaceNormal( rotm.inverse()*(p - center) ) ); 223 } 232 } 224 233 225 ////////////////////////////////////////////// 234 /////////////////////////////////////////////////////////////////////////////// 226 // 235 // 227 // Calculate distance to shape from outside, a 236 // Calculate distance to shape from outside, along normalised vector 228 // 237 // 229 // This CANNOT be an underestimate 238 // This CANNOT be an underestimate 230 // 239 // 231 G4double FastAerosolSolid::DistanceToIn(const 240 G4double FastAerosolSolid::DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const 232 { 241 { 233 G4ThreeVector center; 242 G4ThreeVector center; 234 G4double closestDistance; 243 G4double closestDistance; 235 244 236 if (fCloud->GetNearestDroplet(p, v, center, 245 if (fCloud->GetNearestDroplet(p, v, center, closestDistance, fStepLim, fDroplet, fRotation)) // if we found a droplet within fStepLim of query 237 { 246 { 238 return closestDistance; 247 return closestDistance; 239 } 248 } 240 else if (fCloud->DistanceToCloud(p,v)<DBL_MA 249 else if (fCloud->DistanceToCloud(p,v)<DBL_MAX) // if there is cloud in front of us 241 { 250 { 242 return 1.1*fStepLim; 251 return 1.1*fStepLim; 243 } 252 } 244 else // flying away 253 else // flying away from cloud 245 { 254 { 246 return kInfinity; 255 return kInfinity; 247 } 256 } 248 } 257 } 249 258 250 ////////////////////////////////////////////// 259 ////////////////////////////////////////////////////////////////////// 251 // 260 // 252 // Calculate distance (<= actual) to closest s 261 // Calculate distance (<= actual) to closest surface of shape from outside 253 // 262 // 254 // This function assumes the cloud has at leas 263 // This function assumes the cloud has at least 1 droplet 255 // 264 // 256 // This can be an underestimate 265 // This can be an underestimate 257 // 266 // 258 G4double FastAerosolSolid::DistanceToIn(const 267 G4double FastAerosolSolid::DistanceToIn(const G4ThreeVector &p) const 259 { 268 { 260 G4ThreeVector center; 269 G4ThreeVector center; 261 G4double closestDistance; 270 G4double closestDistance; 262 271 263 G4double distanceToCloud = fBulk->DistanceTo 272 G4double distanceToCloud = fBulk->DistanceToIn(p); 264 273 265 if (fBulk->Inside(p)==kOutside && distanceTo 274 if (fBulk->Inside(p)==kOutside && distanceToCloud>=farFromCloudDist) 266 { 275 { 267 return distanceToCloud; 276 return distanceToCloud; 268 } 277 } 269 else if (fCloud->GetNearestDroplet(p, center 278 else if (fCloud->GetNearestDroplet(p, center, closestDistance, fStepLim, fDroplet, fRotation)) // if we found a droplet within fStepLim of query 270 { 279 { 271 return closestDistance; 280 return closestDistance; 272 } 281 } 273 else 282 else 274 { 283 { 275 return 1.1*fStepLim; 284 return 1.1*fStepLim; 276 } 285 } 277 } 286 } 278 287 279 ////////////////////////////////////////////// 288 ////////////////////////////////////////////////////////////////////// 280 // 289 // 281 // Calculate distance (<= actual) to closest s 290 // Calculate distance (<= actual) to closest surface of shape from inside 282 // 291 // 283 // Despite being a vector distance, we find th 292 // Despite being a vector distance, we find the absolutely closest 284 // droplet to our point since we assume that p 293 // droplet to our point since we assume that p is in a droplet and p 285 // could be past the center 294 // could be past the center 286 // 295 // 287 // This CANNOT be an underestimate 296 // This CANNOT be an underestimate 288 // 297 // 289 G4double FastAerosolSolid::DistanceToOut(const 298 G4double FastAerosolSolid::DistanceToOut(const G4ThreeVector &p, 290 const G4ThreeVector &v, 299 const G4ThreeVector &v, 291 const G4bool calcNorm, 300 const G4bool calcNorm, 292 G4bool *validNorm, 301 G4bool *validNorm, 293 G4ThreeVector *n) const 302 G4ThreeVector *n) const 294 { 303 { 295 G4ThreeVector center; 304 G4ThreeVector center; 296 G4double distanceToIn; // should be 0 305 G4double distanceToIn; // should be 0 297 306 298 fCloud->GetNearestDroplet(p, center, distanc 307 fCloud->GetNearestDroplet(p, center, distanceToIn, fR, fDroplet, fRotation); // if we call this function, must be inside and thus must have a droplet within fR 299 308 300 G4RotationMatrix rotm = fRotation(center); 309 G4RotationMatrix rotm = fRotation(center); 301 G4RotationMatrix irotm = rotm.inverse(); 310 G4RotationMatrix irotm = rotm.inverse(); 302 311 303 G4ThreeVector relPos = irotm*(p-center); 312 G4ThreeVector relPos = irotm*(p-center); 304 313 305 if (fDroplet->Inside(relPos) == kOutside) // 314 if (fDroplet->Inside(relPos) == kOutside) // something went wrong... we should be inside 306 { 315 { 307 std::ostringstream message; 316 std::ostringstream message; 308 message << std::setprecision(15) << "The p 317 message << std::setprecision(15) << "The particle at point p = " << p/mm << "mm" 309 << std::setprecision(15) << " called D 318 << std::setprecision(15) << " called DistanceToOut(p,v) and found the closest droplet to be at center = " << center/mm << "mm" 310 << " but p is outside the droplet!"; 319 << " but p is outside the droplet!"; 311 G4Exception("FastAerosolSolid::DistanceToO 320 G4Exception("FastAerosolSolid::DistanceToOut()", "GeomSolids0002", 312 FatalErrorInArgument, message); 321 FatalErrorInArgument, message); 313 } 322 } 314 323 315 G4double dist = fDroplet->DistanceToOut(relP 324 G4double dist = fDroplet->DistanceToOut(relPos, irotm*v, calcNorm, validNorm, n); 316 *n = rotm*(*n); 325 *n = rotm*(*n); 317 *validNorm = false; // even if droplet is co 326 *validNorm = false; // even if droplet is convex, the aerosol isn't 318 327 319 return dist; 328 return dist; 320 } 329 } 321 330 322 ////////////////////////////////////////////// 331 ///////////////////////////////////////////////////////////////////////// 323 // 332 // 324 // Calculate distance (<=actual) to closest su 333 // Calculate distance (<=actual) to closest surface of shape from inside 325 // 334 // 326 // This can be an underestimate 335 // This can be an underestimate 327 // 336 // 328 G4double FastAerosolSolid::DistanceToOut(const 337 G4double FastAerosolSolid::DistanceToOut(const G4ThreeVector &p) const 329 { 338 { 330 G4ThreeVector center; 339 G4ThreeVector center; 331 G4double distanceToIn; // should be 0 340 G4double distanceToIn; // should be 0 332 341 333 fCloud->GetNearestDroplet(p, center, distanc 342 fCloud->GetNearestDroplet(p, center, distanceToIn, fR, fDroplet, fRotation); // if we call this function, must be inside and thus must have a droplet within fR 334 343 335 G4RotationMatrix irotm = fRotation(center).i 344 G4RotationMatrix irotm = fRotation(center).inverse(); 336 G4ThreeVector relPos = irotm*(p-center); 345 G4ThreeVector relPos = irotm*(p-center); 337 346 338 if (fDroplet->Inside(relPos) == kOutside) // 347 if (fDroplet->Inside(relPos) == kOutside) // something went wrong... we should be inside 339 { 348 { 340 std::ostringstream message; 349 std::ostringstream message; 341 message << "The particle at point p = " << 350 message << "The particle at point p = " << p/mm << "mm" 342 << " called DistanceToOut(p) and found 351 << " called DistanceToOut(p) and found the closest droplet to be at center = " << center/mm << "mm" 343 << " but p is outside the droplet!"; 352 << " but p is outside the droplet!"; 344 G4Exception("FastAerosolSolid::DistanceToO 353 G4Exception("FastAerosolSolid::DistanceToOut()", "GeomSolids0002", 345 FatalErrorInArgument, message); 354 FatalErrorInArgument, message); 346 } 355 } 347 356 348 return fDroplet->DistanceToOut(relPos); 357 return fDroplet->DistanceToOut(relPos); 349 } 358 } 350 359 351 ////////////////////////////////////////////// 360 ////////////////////////////////////////////////////////////////////////// 352 // 361 // 353 // G4EntityType 362 // G4EntityType 354 // 363 // 355 G4GeometryType FastAerosolSolid::GetEntityType 364 G4GeometryType FastAerosolSolid::GetEntityType() const 356 { 365 { 357 return G4String("FastAerosolSolid"); 366 return G4String("FastAerosolSolid"); 358 } 367 } 359 368 360 ////////////////////////////////////////////// 369 ////////////////////////////////////////////////////////////////////////// 361 // 370 // 362 // G4EntityType 371 // G4EntityType 363 // 372 // 364 G4VSolid* FastAerosolSolid::Clone() const 373 G4VSolid* FastAerosolSolid::Clone() const 365 { 374 { 366 return new FastAerosolSolid(*this); 375 return new FastAerosolSolid(*this); 367 } 376 } 368 377 369 ////////////////////////////////////////////// 378 ////////////////////////////////////////////////////////////////////////// 370 // 379 // 371 // Stream object contents to an output stream 380 // Stream object contents to an output stream 372 // 381 // 373 std::ostream &FastAerosolSolid::StreamInfo(std 382 std::ostream &FastAerosolSolid::StreamInfo(std::ostream &os) const 374 { 383 { 375 os << "------------------------------------- 384 os << "-----------------------------------------------------------\n" 376 << " *** Dump for solid - " << GetName 385 << " *** Dump for solid - " << GetName() << " ***\n" 377 << " ================================= 386 << " ===================================================\n" 378 << " Solid type: FastAerosolSolid\n" 387 << " Solid type: FastAerosolSolid\n" 379 << " Parameters: \n" 388 << " Parameters: \n" 380 << " numDroplets: " << fCloud->GetNumD 389 << " numDroplets: " << fCloud->GetNumDroplets() << "\n" 381 << " fDroplet type: " << fDroplet->Get 390 << " fDroplet type: " << fDroplet->GetName() << "\n" 382 << " fDroplet parameters: \n"; 391 << " fDroplet parameters: \n"; 383 fDroplet->StreamInfo(os); 392 fDroplet->StreamInfo(os); 384 os << "---------------------------------- 393 os << "-----------------------------------------------------------\n"; 385 return os; 394 return os; 386 } 395 } 387 396 388 ////////////////////////////////////////////// 397 //////////////////////////////////////////////////////////////////////////////// 389 // 398 // 390 // GetPointOnSurface 399 // GetPointOnSurface 391 // 400 // 392 // Currently hardcoded to look at all droplets 401 // Currently hardcoded to look at all droplets, not just the populated ones 393 // 402 // 394 G4ThreeVector FastAerosolSolid::GetPointOnSurf 403 G4ThreeVector FastAerosolSolid::GetPointOnSurface() const 395 { 404 { 396 G4ThreeVector center; 405 G4ThreeVector center; 397 G4double closestDistance; 406 G4double closestDistance; 398 407 399 G4double fDx = fCloud->GetXHalfLength(); 408 G4double fDx = fCloud->GetXHalfLength(); 400 G4double fDy = fCloud->GetYHalfLength(); 409 G4double fDy = fCloud->GetYHalfLength(); 401 G4double fDz = fCloud->GetZHalfLength(); 410 G4double fDz = fCloud->GetZHalfLength(); 402 411 403 G4ThreeVector p(2.0*fDx*G4UniformRand(),2.0* 412 G4ThreeVector p(2.0*fDx*G4UniformRand(),2.0*fDy*G4UniformRand(),2.0*fDz*G4UniformRand()); 404 p -= G4ThreeVector(fDx, fDy, fDz); 413 p -= G4ThreeVector(fDx, fDy, fDz); 405 414 406 fCloud->GetNearestDroplet(p, center, closest 415 fCloud->GetNearestDroplet(p, center, closestDistance, DBL_MAX, fDroplet, fRotation); 407 416 408 return(center + fRotation(center)*fDroplet-> 417 return(center + fRotation(center)*fDroplet->GetPointOnSurface()); 409 } 418 } 410 419 411 420 412 ////////////////////////////////////////////// 421 ///////////////////////////////////////////////////////////////////////////// 413 // 422 // 414 // Methods for visualisation 423 // Methods for visualisation 415 // 424 // 416 void FastAerosolSolid::DescribeYourselfTo (G4V 425 void FastAerosolSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const 417 { 426 { 418 scene.AddSolid(*this); 427 scene.AddSolid(*this); 419 } 428 } 420 429 421 G4VisExtent FastAerosolSolid::GetExtent() cons 430 G4VisExtent FastAerosolSolid::GetExtent() const 422 { 431 { 423 return G4VisExtent (-fVisDx, fVisDx, -fVisDy 432 return G4VisExtent (-fVisDx, fVisDx, -fVisDy, fVisDy, -fVisDz, fVisDz); 424 } 433 } 425 434 426 G4Polyhedron* FastAerosolSolid::CreatePolyhedr 435 G4Polyhedron* FastAerosolSolid::CreatePolyhedron () const 427 { 436 { 428 return fBulk->CreatePolyhedron(); 437 return fBulk->CreatePolyhedron(); 429 } 438 } 430 439 431 440 432 // copied from G4Ellipsoid 441 // copied from G4Ellipsoid 433 G4Polyhedron* FastAerosolSolid::GetPolyhedron 442 G4Polyhedron* FastAerosolSolid::GetPolyhedron () const 434 { 443 { 435 if (!fpPolyhedron || 444 if (!fpPolyhedron || 436 fRebuildPolyhedron || 445 fRebuildPolyhedron || 437 fpPolyhedron->GetNumberOfRotationStepsAtTi 446 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != 438 fpPolyhedron->GetNumberOfRotationSteps()) 447 fpPolyhedron->GetNumberOfRotationSteps()) 439 { 448 { 440 G4AutoLock l(&polyhedronMutex); 449 G4AutoLock l(&polyhedronMutex); 441 delete fpPolyhedron; 450 delete fpPolyhedron; 442 fpPolyhedron = CreatePolyhedron(); 451 fpPolyhedron = CreatePolyhedron(); 443 fRebuildPolyhedron = false; 452 fRebuildPolyhedron = false; 444 l.unlock(); 453 l.unlock(); 445 } 454 } 446 return fpPolyhedron; 455 return fpPolyhedron; 447 } << 456 } 448