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