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 // G4CutTubs implementation << 27 // 26 // 28 // 01.06.11 T.Nikitina - Derived from G4Tubs << 27 // >> 28 // >> 29 // class G4CutTubs >> 30 // >> 31 // History: >> 32 // 29 // 30.10.16 E.Tcherniaev - reimplemented Calcu 33 // 30.10.16 E.Tcherniaev - reimplemented CalculateExtent(), 30 // removed CreateRotat << 34 // removed CreateRotatedVetices() 31 // ------------------------------------------- << 35 // 05.04.12 M.Kelsey - GetPointOnSurface() throw flat in sqrt(r) >> 36 // 01.06.11 T.Nikitina - Derived from G4Tubs >> 37 // >> 38 ///////////////////////////////////////////////////////////////////////// 32 39 33 #include "G4CutTubs.hh" 40 #include "G4CutTubs.hh" 34 41 35 #if !defined(G4GEOM_USE_UCTUBS) 42 #if !defined(G4GEOM_USE_UCTUBS) 36 43 37 #include "G4GeomTools.hh" 44 #include "G4GeomTools.hh" 38 #include "G4VoxelLimits.hh" 45 #include "G4VoxelLimits.hh" 39 #include "G4AffineTransform.hh" 46 #include "G4AffineTransform.hh" 40 #include "G4GeometryTolerance.hh" 47 #include "G4GeometryTolerance.hh" 41 #include "G4BoundingEnvelope.hh" 48 #include "G4BoundingEnvelope.hh" 42 49 43 #include "G4VPVParameterisation.hh" 50 #include "G4VPVParameterisation.hh" 44 #include "G4QuickRand.hh" << 45 51 46 #include "G4VGraphicsScene.hh" << 52 #include "Randomize.hh" 47 #include "G4Polyhedron.hh" << 48 53 49 #include "G4AutoLock.hh" << 54 #include "meshdefs.hh" 50 55 51 namespace << 56 #include "G4VGraphicsScene.hh" 52 { << 53 G4Mutex zminmaxMutex = G4MUTEX_INITIALIZER; << 54 } << 55 57 56 using namespace CLHEP; 58 using namespace CLHEP; 57 59 58 ////////////////////////////////////////////// 60 ///////////////////////////////////////////////////////////////////////// 59 // 61 // 60 // Constructor - check parameters, convert ang 62 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI 61 // - note if pdphi>2PI then reset 63 // - note if pdphi>2PI then reset to 2PI 62 64 63 G4CutTubs::G4CutTubs( const G4String &pName, 65 G4CutTubs::G4CutTubs( const G4String &pName, 64 G4double pRMin, G4double 66 G4double pRMin, G4double pRMax, 65 G4double pDz, 67 G4double pDz, 66 G4double pSPhi, G4double 68 G4double pSPhi, G4double pDPhi, 67 G4ThreeVector pLowNorm,G 69 G4ThreeVector pLowNorm,G4ThreeVector pHighNorm ) 68 : G4CSGSolid(pName), fRMin(pRMin), fRMax(pRM << 70 : G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0) 69 fSPhi(0.), fDPhi(0.), fZMin(0.), fZMax(0.) << 70 { 71 { 71 kRadTolerance = G4GeometryTolerance::GetInst 72 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 72 kAngTolerance = G4GeometryTolerance::GetInst 73 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 73 74 74 halfCarTolerance = kCarTolerance*0.5; 75 halfCarTolerance = kCarTolerance*0.5; 75 halfRadTolerance = kRadTolerance*0.5; 76 halfRadTolerance = kRadTolerance*0.5; 76 halfAngTolerance = kAngTolerance*0.5; 77 halfAngTolerance = kAngTolerance*0.5; 77 78 78 if (pDz<=0) // Check z-len 79 if (pDz<=0) // Check z-len 79 { 80 { 80 std::ostringstream message; 81 std::ostringstream message; 81 message << "Negative Z half-length (" << p 82 message << "Negative Z half-length (" << pDz << ") in solid: " << GetName(); 82 G4Exception("G4CutTubs::G4CutTubs()", "Geo 83 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002", FatalException, message); 83 } 84 } 84 if ( (pRMin >= pRMax) || (pRMin < 0) ) // Ch 85 if ( (pRMin >= pRMax) || (pRMin < 0) ) // Check radii 85 { 86 { 86 std::ostringstream message; 87 std::ostringstream message; 87 message << "Invalid values for radii in so 88 message << "Invalid values for radii in solid: " << GetName() 88 << G4endl 89 << G4endl 89 << " pRMin = " << pRMin << 90 << " pRMin = " << pRMin << ", pRMax = " << pRMax; 90 G4Exception("G4CutTubs::G4CutTubs()", "Geo 91 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002", FatalException, message); 91 } 92 } 92 93 93 // Check angles << 94 // Check angles 94 // << 95 // 95 CheckPhiAngles(pSPhi, pDPhi); 96 CheckPhiAngles(pSPhi, pDPhi); 96 97 97 // Check on Cutted Planes Normals 98 // Check on Cutted Planes Normals 98 // If there is NO CUT, propose to use G4Tubs 99 // If there is NO CUT, propose to use G4Tubs instead 99 // 100 // 100 if ( ( pLowNorm.x() == 0.0) && ( pLowNorm.y( << 101 if ( ( !pLowNorm.x()) && ( !pLowNorm.y()) 101 && ( pHighNorm.x() == 0.0) && (pHighNorm.y << 102 && ( !pHighNorm.x()) && (!pHighNorm.y()) ) 102 { 103 { 103 std::ostringstream message; 104 std::ostringstream message; 104 message << "Inexisting Low/High Normal to 105 message << "Inexisting Low/High Normal to Z plane or Parallel to Z." 105 << G4endl 106 << G4endl 106 << "Normals to Z plane are " << pL << 107 << "Normals to Z plane are (" << pLowNorm <<" and " 107 << pHighNorm << " in solid: " << G << 108 << pHighNorm << ") in solid: " << GetName(); 108 G4Exception("G4CutTubs::G4CutTubs()", "Geo 109 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids1001", 109 JustWarning, message, "Should 110 JustWarning, message, "Should use G4Tubs!"); 110 } 111 } 111 112 112 // If Normal is (0,0,0),means parallel to R, 113 // If Normal is (0,0,0),means parallel to R, give it value of (0,0,+/-1) 113 // << 114 // 114 if (pLowNorm.mag2() == 0.) { pLowNorm.setZ( 115 if (pLowNorm.mag2() == 0.) { pLowNorm.setZ(-1.); } 115 if (pHighNorm.mag2()== 0.) { pHighNorm.setZ 116 if (pHighNorm.mag2()== 0.) { pHighNorm.setZ(1.); } 116 117 117 // Given Normals to Cut Planes have to be an << 118 // Given Normals to Cut Planes have to be an unit vectors. 118 // Normalize if it is needed. 119 // Normalize if it is needed. 119 // 120 // 120 if (pLowNorm.mag2() != 1.) { pLowNorm = pL 121 if (pLowNorm.mag2() != 1.) { pLowNorm = pLowNorm.unit(); } 121 if (pHighNorm.mag2()!= 1.) { pHighNorm = pH 122 if (pHighNorm.mag2()!= 1.) { pHighNorm = pHighNorm.unit(); } 122 123 123 // Normals to cutted planes have to point ou 124 // Normals to cutted planes have to point outside Solid 124 // 125 // 125 if( (pLowNorm.mag2() != 0.) && (pHighNorm.ma 126 if( (pLowNorm.mag2() != 0.) && (pHighNorm.mag2()!= 0. ) ) 126 { 127 { 127 if( ( pLowNorm.z()>= 0. ) || ( pHighNorm.z 128 if( ( pLowNorm.z()>= 0. ) || ( pHighNorm.z() <= 0.)) 128 { 129 { 129 std::ostringstream message; 130 std::ostringstream message; 130 message << "Invalid Low or High Normal t 131 message << "Invalid Low or High Normal to Z plane; " 131 "has to point outside Solid." 132 "has to point outside Solid." << G4endl 132 << "Invalid Norm to Z plane (" < 133 << "Invalid Norm to Z plane (" << pLowNorm << " or " 133 << pHighNorm << ") in solid: " < 134 << pHighNorm << ") in solid: " << GetName(); 134 G4Exception("G4CutTubs::G4CutTubs()", "G 135 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002", 135 FatalException, message); 136 FatalException, message); 136 } 137 } 137 } 138 } 138 fLowNorm = pLowNorm; 139 fLowNorm = pLowNorm; 139 fHighNorm = pHighNorm; 140 fHighNorm = pHighNorm; 140 141 141 // Check intersection of cut planes, they MU << 142 // Check Intersection of cut planes. They MUST NOT Intersect 142 // each other inside the lateral surface << 143 // 143 // 144 if(IsCrossingCutPlanes()) << 144 // This check has been disabled as too strict. 145 { << 145 // See problem report #1887 146 std::ostringstream message; << 146 // 147 message << "Invalid normals to Z plane in << 147 // if(IsCrossingCutPlanes()) 148 << "Cut planes are crossing inside << 148 // { 149 << " Solid type: G4CutTubs\n" << 149 // std::ostringstream message; 150 << " Parameters: \n" << 150 // message << "Invalid Low or High Normal to Z plane; " 151 << " inner radius : " << fRMin/ << 151 // << "Crossing Cutted Planes." << G4endl 152 << " outer radius : " << fRMax/ << 152 // << "Invalid Norm to Z plane (" << pLowNorm << " and " 153 << " half length Z: " << fDz/mm << 153 // << pHighNorm << ") in solid: " << GetName(); 154 << " starting phi : " << fSPhi/ << 154 // G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002", 155 << " delta phi : " << fDPhi/ << 155 // FatalException, message); 156 << " low Norm : " << fLowNo << 156 // } 157 << " high Norm : " << fHighN << 158 G4Exception("G4CutTubs::G4CutTubs()", "Geo << 159 FatalException, message); << 160 } << 161 } 157 } 162 158 163 ////////////////////////////////////////////// 159 /////////////////////////////////////////////////////////////////////// 164 // 160 // 165 // Fake default constructor - sets only member 161 // Fake default constructor - sets only member data and allocates memory 166 // for usage restri 162 // for usage restricted to object persistency. 167 // 163 // 168 G4CutTubs::G4CutTubs( __void__& a ) 164 G4CutTubs::G4CutTubs( __void__& a ) 169 : G4CSGSolid(a) << 165 : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.), >> 166 fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.), >> 167 sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.), >> 168 sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.), >> 169 fPhiFullCutTube(false), >> 170 halfCarTolerance(0.), halfRadTolerance(0.), halfAngTolerance(0.), >> 171 fLowNorm(G4ThreeVector()), fHighNorm(G4ThreeVector()) 170 { 172 { 171 } 173 } 172 174 173 ////////////////////////////////////////////// 175 ////////////////////////////////////////////////////////////////////////// 174 // 176 // 175 // Destructor 177 // Destructor 176 178 177 G4CutTubs::~G4CutTubs() = default; << 179 G4CutTubs::~G4CutTubs() >> 180 { >> 181 } 178 182 179 ////////////////////////////////////////////// 183 ////////////////////////////////////////////////////////////////////////// 180 // 184 // 181 // Copy constructor 185 // Copy constructor 182 186 183 G4CutTubs::G4CutTubs(const G4CutTubs&) = defau << 187 G4CutTubs::G4CutTubs(const G4CutTubs& rhs) >> 188 : G4CSGSolid(rhs), >> 189 kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance), >> 190 fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz), >> 191 fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi), >> 192 sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi), >> 193 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT), >> 194 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), >> 195 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), >> 196 fPhiFullCutTube(rhs.fPhiFullCutTube), >> 197 halfCarTolerance(rhs.halfCarTolerance), >> 198 halfRadTolerance(rhs.halfRadTolerance), >> 199 halfAngTolerance(rhs.halfAngTolerance), >> 200 fLowNorm(rhs.fLowNorm), fHighNorm(rhs.fHighNorm) >> 201 { >> 202 } 184 203 185 ////////////////////////////////////////////// 204 ////////////////////////////////////////////////////////////////////////// 186 // 205 // 187 // Assignment operator 206 // Assignment operator 188 207 189 G4CutTubs& G4CutTubs::operator = (const G4CutT << 208 G4CutTubs& G4CutTubs::operator = (const G4CutTubs& rhs) 190 { 209 { 191 // Check assignment to self 210 // Check assignment to self 192 // 211 // 193 if (this == &rhs) { return *this; } 212 if (this == &rhs) { return *this; } 194 213 195 // Copy base class data 214 // Copy base class data 196 // 215 // 197 G4CSGSolid::operator=(rhs); 216 G4CSGSolid::operator=(rhs); 198 217 199 // Copy data 218 // Copy data 200 // 219 // 201 kRadTolerance = rhs.kRadTolerance; kAngTole 220 kRadTolerance = rhs.kRadTolerance; kAngTolerance = rhs.kAngTolerance; 202 fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = 221 fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = rhs.fDz; 203 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; 222 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; 204 fZMin = rhs.fZMin; fZMax = rhs.fZMax; << 205 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPh 223 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi; 206 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r 224 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT; 207 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh 225 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi; 208 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh 226 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi; 209 fPhiFullCutTube = rhs.fPhiFullCutTube; 227 fPhiFullCutTube = rhs.fPhiFullCutTube; 210 halfCarTolerance = rhs.halfCarTolerance; 228 halfCarTolerance = rhs.halfCarTolerance; 211 halfRadTolerance = rhs.halfRadTolerance; 229 halfRadTolerance = rhs.halfRadTolerance; 212 halfAngTolerance = rhs.halfAngTolerance; 230 halfAngTolerance = rhs.halfAngTolerance; 213 fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fH 231 fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fHighNorm; 214 232 215 return *this; 233 return *this; 216 } 234 } 217 235 218 ////////////////////////////////////////////// 236 ////////////////////////////////////////////////////////////////////////// 219 // 237 // 220 // Get volume << 221 << 222 G4double G4CutTubs::GetCubicVolume() << 223 { << 224 constexpr G4int nphi = 200, nrho = 100; << 225 << 226 if (fCubicVolume == 0.) << 227 { << 228 // get parameters << 229 G4double rmin = GetInnerRadius(); << 230 G4double rmax = GetOuterRadius(); << 231 G4double dz = GetZHalfLength(); << 232 G4double sphi = GetStartPhiAngle(); << 233 G4double dphi = GetDeltaPhiAngle(); << 234 << 235 // calculate volume << 236 G4double volume = dz*dphi*(rmax*rmax - rmi << 237 if (dphi < twopi) // make recalculation << 238 { << 239 // set values for calculation of h - dis << 240 // opposite points on bases << 241 G4ThreeVector nbot = GetLowNorm(); << 242 G4ThreeVector ntop = GetHighNorm(); << 243 G4double nx = nbot.x()/nbot.z() - ntop.x << 244 G4double ny = nbot.y()/nbot.z() - ntop.y << 245 << 246 // compute volume by integration << 247 G4double delrho = (rmax - rmin)/nrho; << 248 G4double delphi = dphi/nphi; << 249 volume = 0.; << 250 for (G4int irho=0; irho<nrho; ++irho) << 251 { << 252 G4double r1 = rmin + delrho*irho; << 253 G4double r2 = rmin + delrho*(irho + 1 << 254 G4double rho = 0.5*(r1 + r2); << 255 G4double sector = 0.5*delphi*(r2*r2 - << 256 for (G4int iphi=0; iphi<nphi; ++iphi) << 257 { << 258 G4double phi = sphi + delphi*(iphi + << 259 G4double h = nx*rho*std::cos(phi) + << 260 volume += sector*h; << 261 } << 262 } << 263 } << 264 fCubicVolume = volume; << 265 } << 266 return fCubicVolume; << 267 } << 268 << 269 ////////////////////////////////////////////// << 270 // << 271 // Get surface area << 272 << 273 G4double G4CutTubs::GetSurfaceArea() << 274 { << 275 constexpr G4int nphi = 400; << 276 << 277 if (fSurfaceArea == 0.) << 278 { << 279 // get parameters << 280 G4double rmin = GetInnerRadius(); << 281 G4double rmax = GetOuterRadius(); << 282 G4double dz = GetZHalfLength(); << 283 G4double sphi = GetStartPhiAngle(); << 284 G4double dphi = GetDeltaPhiAngle(); << 285 G4ThreeVector nbot = GetLowNorm(); << 286 G4ThreeVector ntop = GetHighNorm(); << 287 << 288 // calculate lateral surface area << 289 G4double sinner = 2.*dz*dphi*rmin; << 290 G4double souter = 2.*dz*dphi*rmax; << 291 if (dphi < twopi) // make recalculation << 292 { << 293 // set values for calculation of h - dis << 294 // opposite points on bases << 295 G4double nx = nbot.x()/nbot.z() - ntop.x << 296 G4double ny = nbot.y()/nbot.z() - ntop.y << 297 << 298 // compute lateral surface area by integ << 299 G4double delphi = dphi/nphi; << 300 sinner = 0.; << 301 souter = 0.; << 302 for (G4int iphi=0; iphi<nphi; ++iphi) << 303 { << 304 G4double phi = sphi + delphi*(iphi + 0 << 305 G4double cosphi = std::cos(phi); << 306 G4double sinphi = std::sin(phi); << 307 sinner += rmin*(nx*cosphi + ny*sinphi) << 308 souter += rmax*(nx*cosphi + ny*sinphi) << 309 } << 310 sinner *= delphi*rmin; << 311 souter *= delphi*rmax; << 312 } << 313 // set surface area << 314 G4double scut = (dphi == twopi) ? 0. : 2. << 315 G4double szero = 0.5*dphi*(rmax*rmax - rmi << 316 G4double slow = szero/std::abs(nbot.z()); << 317 G4double shigh = szero/std::abs(ntop.z()); << 318 fSurfaceArea = sinner + souter + 2.*scut + << 319 } << 320 return fSurfaceArea; << 321 } << 322 << 323 ////////////////////////////////////////////// << 324 // << 325 // Get bounding box 238 // Get bounding box 326 239 327 void G4CutTubs::BoundingLimits(G4ThreeVector& 240 void G4CutTubs::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const 328 { 241 { 329 G4double rmin = GetInnerRadius(); 242 G4double rmin = GetInnerRadius(); 330 G4double rmax = GetOuterRadius(); 243 G4double rmax = GetOuterRadius(); 331 G4double dz = GetZHalfLength(); 244 G4double dz = GetZHalfLength(); 332 G4double dphi = GetDeltaPhiAngle(); 245 G4double dphi = GetDeltaPhiAngle(); 333 246 334 G4double sinSphi = GetSinStartPhi(); << 247 G4double sinSphi = GetSinStartPhi(); 335 G4double cosSphi = GetCosStartPhi(); << 248 G4double cosSphi = GetCosStartPhi(); 336 G4double sinEphi = GetSinEndPhi(); << 249 G4double sinEphi = GetSinEndPhi(); 337 G4double cosEphi = GetCosEndPhi(); << 250 G4double cosEphi = GetCosEndPhi(); 338 251 339 G4ThreeVector norm; 252 G4ThreeVector norm; 340 G4double mag, topx, topy, dists, diste; << 253 G4double mag, topx, topy, dists, diste; 341 G4bool iftop; 254 G4bool iftop; 342 255 343 // Find Zmin 256 // Find Zmin 344 // 257 // 345 G4double zmin; 258 G4double zmin; 346 norm = GetLowNorm(); 259 norm = GetLowNorm(); 347 mag = std::sqrt(norm.x()*norm.x() + norm.y( 260 mag = std::sqrt(norm.x()*norm.x() + norm.y()*norm.y()); 348 topx = (mag == 0) ? 0 : -rmax*norm.x()/mag; << 261 topx = (mag == 0) ? 0 : -rmax*norm.x()/mag; 349 topy = (mag == 0) ? 0 : -rmax*norm.y()/mag; << 262 topy = (mag == 0) ? 0 : -rmax*norm.y()/mag; 350 dists = sinSphi*topx - cosSphi*topy; 263 dists = sinSphi*topx - cosSphi*topy; 351 diste = -sinEphi*topx + cosEphi*topy; 264 diste = -sinEphi*topx + cosEphi*topy; 352 if (dphi > pi) 265 if (dphi > pi) 353 { 266 { 354 iftop = true; 267 iftop = true; 355 if (dists > 0 && diste > 0)iftop = false; 268 if (dists > 0 && diste > 0)iftop = false; 356 } 269 } 357 else 270 else 358 { 271 { 359 iftop = false; 272 iftop = false; 360 if (dists <= 0 && diste <= 0) iftop = true 273 if (dists <= 0 && diste <= 0) iftop = true; 361 } 274 } 362 if (iftop) 275 if (iftop) 363 { 276 { 364 zmin = -(norm.x()*topx + norm.y()*topy)/no 277 zmin = -(norm.x()*topx + norm.y()*topy)/norm.z() - dz; 365 } 278 } 366 else 279 else 367 { 280 { 368 G4double z1 = -rmin*(norm.x()*cosSphi + no << 281 G4double z1 = -rmin*(norm.x()*cosSphi + norm.y()*sinSphi)/norm.z() - dz; 369 G4double z2 = -rmin*(norm.x()*cosEphi + no << 282 G4double z2 = -rmin*(norm.x()*cosEphi + norm.y()*sinEphi)/norm.z() - dz; 370 G4double z3 = -rmax*(norm.x()*cosSphi + no << 283 G4double z3 = -rmax*(norm.x()*cosSphi + norm.y()*sinSphi)/norm.z() - dz; 371 G4double z4 = -rmax*(norm.x()*cosEphi + no << 284 G4double z4 = -rmax*(norm.x()*cosEphi + norm.y()*sinEphi)/norm.z() - dz; 372 zmin = std::min(std::min(std::min(z1,z2),z 285 zmin = std::min(std::min(std::min(z1,z2),z3),z4); 373 } 286 } 374 287 375 // Find Zmax 288 // Find Zmax 376 // 289 // 377 G4double zmax; 290 G4double zmax; 378 norm = GetHighNorm(); 291 norm = GetHighNorm(); 379 mag = std::sqrt(norm.x()*norm.x() + norm.y( 292 mag = std::sqrt(norm.x()*norm.x() + norm.y()*norm.y()); 380 topx = (mag == 0) ? 0 : -rmax*norm.x()/mag; << 293 topx = (mag == 0) ? 0 : -rmax*norm.x()/mag; 381 topy = (mag == 0) ? 0 : -rmax*norm.y()/mag; << 294 topy = (mag == 0) ? 0 : -rmax*norm.y()/mag; 382 dists = sinSphi*topx - cosSphi*topy; 295 dists = sinSphi*topx - cosSphi*topy; 383 diste = -sinEphi*topx + cosEphi*topy; 296 diste = -sinEphi*topx + cosEphi*topy; 384 if (dphi > pi) 297 if (dphi > pi) 385 { 298 { 386 iftop = true; 299 iftop = true; 387 if (dists > 0 && diste > 0) iftop = false; 300 if (dists > 0 && diste > 0) iftop = false; 388 } 301 } 389 else 302 else 390 { 303 { 391 iftop = false; 304 iftop = false; 392 if (dists <= 0 && diste <= 0) iftop = true 305 if (dists <= 0 && diste <= 0) iftop = true; 393 } 306 } 394 if (iftop) 307 if (iftop) 395 { 308 { 396 zmax = -(norm.x()*topx + norm.y()*topy)/no 309 zmax = -(norm.x()*topx + norm.y()*topy)/norm.z() + dz; 397 } 310 } 398 else 311 else 399 { 312 { 400 G4double z1 = -rmin*(norm.x()*cosSphi + no << 313 G4double z1 = -rmin*(norm.x()*cosSphi + norm.y()*sinSphi)/norm.z() + dz; 401 G4double z2 = -rmin*(norm.x()*cosEphi + no << 314 G4double z2 = -rmin*(norm.x()*cosEphi + norm.y()*sinEphi)/norm.z() + dz; 402 G4double z3 = -rmax*(norm.x()*cosSphi + no << 315 G4double z3 = -rmax*(norm.x()*cosSphi + norm.y()*sinSphi)/norm.z() + dz; 403 G4double z4 = -rmax*(norm.x()*cosEphi + no << 316 G4double z4 = -rmax*(norm.x()*cosEphi + norm.y()*sinEphi)/norm.z() + dz; 404 zmax = std::max(std::max(std::max(z1,z2),z 317 zmax = std::max(std::max(std::max(z1,z2),z3),z4); 405 } 318 } 406 319 407 // Find bounding box 320 // Find bounding box 408 // 321 // 409 if (dphi < twopi) 322 if (dphi < twopi) 410 { 323 { 411 G4TwoVector vmin,vmax; 324 G4TwoVector vmin,vmax; 412 G4GeomTools::DiskExtent(rmin,rmax, 325 G4GeomTools::DiskExtent(rmin,rmax, 413 GetSinStartPhi(),G 326 GetSinStartPhi(),GetCosStartPhi(), 414 GetSinEndPhi(),Get 327 GetSinEndPhi(),GetCosEndPhi(), 415 vmin,vmax); 328 vmin,vmax); 416 pMin.set(vmin.x(),vmin.y(), zmin); 329 pMin.set(vmin.x(),vmin.y(), zmin); 417 pMax.set(vmax.x(),vmax.y(), zmax); 330 pMax.set(vmax.x(),vmax.y(), zmax); 418 } 331 } 419 else 332 else 420 { 333 { 421 pMin.set(-rmax,-rmax, zmin); 334 pMin.set(-rmax,-rmax, zmin); 422 pMax.set( rmax, rmax, zmax); 335 pMax.set( rmax, rmax, zmax); 423 } 336 } 424 337 425 // Check correctness of the bounding box 338 // Check correctness of the bounding box 426 // 339 // 427 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 340 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 428 { 341 { 429 std::ostringstream message; 342 std::ostringstream message; 430 message << "Bad bounding box (min >= max) 343 message << "Bad bounding box (min >= max) for solid: " 431 << GetName() << " !" 344 << GetName() << " !" 432 << "\npMin = " << pMin 345 << "\npMin = " << pMin 433 << "\npMax = " << pMax; 346 << "\npMax = " << pMax; 434 G4Exception("G4CutTubs::BoundingLimits()", 347 G4Exception("G4CutTubs::BoundingLimits()", "GeomMgt0001", 435 JustWarning, message); 348 JustWarning, message); 436 DumpInfo(); 349 DumpInfo(); 437 } 350 } 438 } 351 } 439 352 440 ////////////////////////////////////////////// 353 ////////////////////////////////////////////////////////////////////////// 441 // 354 // 442 // Calculate extent under transform and specif 355 // Calculate extent under transform and specified limit 443 356 444 G4bool G4CutTubs::CalculateExtent( const EAxis 357 G4bool G4CutTubs::CalculateExtent( const EAxis pAxis, 445 const G4Vox 358 const G4VoxelLimits& pVoxelLimit, 446 const G4Aff 359 const G4AffineTransform& pTransform, 447 G4dou << 360 G4double& pMin, 448 G4dou 361 G4double& pMax ) const 449 { 362 { 450 G4ThreeVector bmin, bmax; 363 G4ThreeVector bmin, bmax; 451 G4bool exist; 364 G4bool exist; 452 365 453 // Get bounding box 366 // Get bounding box 454 BoundingLimits(bmin,bmax); 367 BoundingLimits(bmin,bmax); 455 368 456 // Check bounding box 369 // Check bounding box 457 G4BoundingEnvelope bbox(bmin,bmax); 370 G4BoundingEnvelope bbox(bmin,bmax); 458 #ifdef G4BBOX_EXTENT 371 #ifdef G4BBOX_EXTENT 459 return bbox.CalculateExtent(pAxis,pVoxelLimi << 372 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 460 #endif 373 #endif 461 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 374 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax)) 462 { 375 { 463 return exist = pMin < pMax; << 376 return exist = (pMin < pMax) ? true : false; 464 } 377 } 465 378 466 // Get parameters of the solid 379 // Get parameters of the solid 467 G4double rmin = GetInnerRadius(); 380 G4double rmin = GetInnerRadius(); 468 G4double rmax = GetOuterRadius(); 381 G4double rmax = GetOuterRadius(); 469 G4double dphi = GetDeltaPhiAngle(); 382 G4double dphi = GetDeltaPhiAngle(); 470 G4double zmin = bmin.z(); 383 G4double zmin = bmin.z(); 471 G4double zmax = bmax.z(); 384 G4double zmax = bmax.z(); 472 385 473 // Find bounding envelope and calculate exte 386 // Find bounding envelope and calculate extent 474 // 387 // 475 const G4int NSTEPS = 24; // numbe 388 const G4int NSTEPS = 24; // number of steps for whole circle 476 G4double astep = twopi/NSTEPS; // max a 389 G4double astep = twopi/NSTEPS; // max angle for one step 477 G4int ksteps = (dphi <= astep) ? 1 : (G4i 390 G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1; 478 G4double ang = dphi/ksteps; 391 G4double ang = dphi/ksteps; 479 392 480 G4double sinHalf = std::sin(0.5*ang); 393 G4double sinHalf = std::sin(0.5*ang); 481 G4double cosHalf = std::cos(0.5*ang); 394 G4double cosHalf = std::cos(0.5*ang); 482 G4double sinStep = 2.*sinHalf*cosHalf; 395 G4double sinStep = 2.*sinHalf*cosHalf; 483 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 396 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 484 G4double rext = rmax/cosHalf; 397 G4double rext = rmax/cosHalf; 485 398 486 // bounding envelope for full cylinder consi 399 // bounding envelope for full cylinder consists of two polygons, 487 // in other cases it is a sequence of quadri 400 // in other cases it is a sequence of quadrilaterals 488 if (rmin == 0 && dphi == twopi) 401 if (rmin == 0 && dphi == twopi) 489 { 402 { 490 G4double sinCur = sinHalf; 403 G4double sinCur = sinHalf; 491 G4double cosCur = cosHalf; 404 G4double cosCur = cosHalf; 492 405 493 G4ThreeVectorList baseA(NSTEPS),baseB(NSTE 406 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS); 494 for (G4int k=0; k<NSTEPS; ++k) 407 for (G4int k=0; k<NSTEPS; ++k) 495 { 408 { 496 baseA[k].set(rext*cosCur,rext*sinCur,zmi 409 baseA[k].set(rext*cosCur,rext*sinCur,zmin); 497 baseB[k].set(rext*cosCur,rext*sinCur,zma 410 baseB[k].set(rext*cosCur,rext*sinCur,zmax); 498 411 499 G4double sinTmp = sinCur; 412 G4double sinTmp = sinCur; 500 sinCur = sinCur*cosStep + cosCur*sinStep 413 sinCur = sinCur*cosStep + cosCur*sinStep; 501 cosCur = cosCur*cosStep - sinTmp*sinStep 414 cosCur = cosCur*cosStep - sinTmp*sinStep; 502 } 415 } 503 std::vector<const G4ThreeVectorList *> pol 416 std::vector<const G4ThreeVectorList *> polygons(2); 504 polygons[0] = &baseA; 417 polygons[0] = &baseA; 505 polygons[1] = &baseB; 418 polygons[1] = &baseB; 506 G4BoundingEnvelope benv(bmin,bmax,polygons 419 G4BoundingEnvelope benv(bmin,bmax,polygons); 507 exist = benv.CalculateExtent(pAxis,pVoxelL 420 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 508 } 421 } 509 else 422 else 510 { 423 { 511 G4double sinStart = GetSinStartPhi(); 424 G4double sinStart = GetSinStartPhi(); 512 G4double cosStart = GetCosStartPhi(); 425 G4double cosStart = GetCosStartPhi(); 513 G4double sinEnd = GetSinEndPhi(); 426 G4double sinEnd = GetSinEndPhi(); 514 G4double cosEnd = GetCosEndPhi(); 427 G4double cosEnd = GetCosEndPhi(); 515 G4double sinCur = sinStart*cosHalf + cos 428 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf; 516 G4double cosCur = cosStart*cosHalf - sin 429 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf; 517 430 518 // set quadrilaterals 431 // set quadrilaterals 519 G4ThreeVectorList pols[NSTEPS+2]; 432 G4ThreeVectorList pols[NSTEPS+2]; 520 for (G4int k=0; k<ksteps+2; ++k) pols[k].r 433 for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4); 521 pols[0][0].set(rmin*cosStart,rmin*sinStart 434 pols[0][0].set(rmin*cosStart,rmin*sinStart,zmax); 522 pols[0][1].set(rmin*cosStart,rmin*sinStart 435 pols[0][1].set(rmin*cosStart,rmin*sinStart,zmin); 523 pols[0][2].set(rmax*cosStart,rmax*sinStart 436 pols[0][2].set(rmax*cosStart,rmax*sinStart,zmin); 524 pols[0][3].set(rmax*cosStart,rmax*sinStart 437 pols[0][3].set(rmax*cosStart,rmax*sinStart,zmax); 525 for (G4int k=1; k<ksteps+1; ++k) 438 for (G4int k=1; k<ksteps+1; ++k) 526 { 439 { 527 pols[k][0].set(rmin*cosCur,rmin*sinCur,z 440 pols[k][0].set(rmin*cosCur,rmin*sinCur,zmax); 528 pols[k][1].set(rmin*cosCur,rmin*sinCur,z 441 pols[k][1].set(rmin*cosCur,rmin*sinCur,zmin); 529 pols[k][2].set(rext*cosCur,rext*sinCur,z 442 pols[k][2].set(rext*cosCur,rext*sinCur,zmin); 530 pols[k][3].set(rext*cosCur,rext*sinCur,z 443 pols[k][3].set(rext*cosCur,rext*sinCur,zmax); 531 444 532 G4double sinTmp = sinCur; 445 G4double sinTmp = sinCur; 533 sinCur = sinCur*cosStep + cosCur*sinStep 446 sinCur = sinCur*cosStep + cosCur*sinStep; 534 cosCur = cosCur*cosStep - sinTmp*sinStep 447 cosCur = cosCur*cosStep - sinTmp*sinStep; 535 } 448 } 536 pols[ksteps+1][0].set(rmin*cosEnd,rmin*sin 449 pols[ksteps+1][0].set(rmin*cosEnd,rmin*sinEnd,zmax); 537 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sin 450 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sinEnd,zmin); 538 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sin 451 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sinEnd,zmin); 539 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sin 452 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sinEnd,zmax); 540 453 541 // set envelope and calculate extent 454 // set envelope and calculate extent 542 std::vector<const G4ThreeVectorList *> pol 455 std::vector<const G4ThreeVectorList *> polygons; 543 polygons.resize(ksteps+2); 456 polygons.resize(ksteps+2); 544 for (G4int k=0; k<ksteps+2; ++k) { polygon << 457 for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k]; 545 G4BoundingEnvelope benv(bmin,bmax,polygons 458 G4BoundingEnvelope benv(bmin,bmax,polygons); 546 exist = benv.CalculateExtent(pAxis,pVoxelL 459 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 547 } 460 } 548 return exist; 461 return exist; 549 } 462 } 550 463 551 ////////////////////////////////////////////// 464 ////////////////////////////////////////////////////////////////////////// 552 // 465 // 553 // Return whether point inside/outside/on surf 466 // Return whether point inside/outside/on surface 554 467 555 EInside G4CutTubs::Inside( const G4ThreeVector 468 EInside G4CutTubs::Inside( const G4ThreeVector& p ) const 556 { 469 { 557 G4ThreeVector vZ = G4ThreeVector(0,0,fDz); 470 G4ThreeVector vZ = G4ThreeVector(0,0,fDz); 558 EInside in = kInside; 471 EInside in = kInside; 559 472 560 // Check the lower cut plane 473 // Check the lower cut plane 561 // 474 // 562 G4double zinLow =(p+vZ).dot(fLowNorm); 475 G4double zinLow =(p+vZ).dot(fLowNorm); 563 if (zinLow > halfCarTolerance) { return kOu 476 if (zinLow > halfCarTolerance) { return kOutside; } 564 477 565 // Check the higher cut plane 478 // Check the higher cut plane 566 // 479 // 567 G4double zinHigh = (p-vZ).dot(fHighNorm); 480 G4double zinHigh = (p-vZ).dot(fHighNorm); 568 if (zinHigh > halfCarTolerance) { return kO 481 if (zinHigh > halfCarTolerance) { return kOutside; } 569 482 570 // Check radius 483 // Check radius 571 // 484 // 572 G4double r2 = p.x()*p.x() + p.y()*p.y() ; 485 G4double r2 = p.x()*p.x() + p.y()*p.y() ; 573 486 574 G4double tolRMin = fRMin - halfRadTolerance; 487 G4double tolRMin = fRMin - halfRadTolerance; 575 G4double tolRMax = fRMax + halfRadTolerance; 488 G4double tolRMax = fRMax + halfRadTolerance; 576 if ( tolRMin < 0 ) { tolRMin = 0; } 489 if ( tolRMin < 0 ) { tolRMin = 0; } 577 490 578 if (r2 < tolRMin*tolRMin || r2 > tolRMax*tol 491 if (r2 < tolRMin*tolRMin || r2 > tolRMax*tolRMax) { return kOutside; } 579 492 580 // Check Phi cut 493 // Check Phi cut 581 // 494 // 582 if(!fPhiFullCutTube) 495 if(!fPhiFullCutTube) 583 { 496 { 584 if ((tolRMin == 0) && (std::fabs(p.x()) <= 497 if ((tolRMin == 0) && (std::fabs(p.x()) <= halfCarTolerance) 585 && (std::fabs(p.y()) <= 498 && (std::fabs(p.y()) <= halfCarTolerance)) 586 { 499 { 587 return kSurface; 500 return kSurface; 588 } 501 } 589 502 590 G4double phi0 = std::atan2(p.y(),p.x()); 503 G4double phi0 = std::atan2(p.y(),p.x()); 591 G4double phi1 = phi0 - twopi; 504 G4double phi1 = phi0 - twopi; 592 G4double phi2 = phi0 + twopi; 505 G4double phi2 = phi0 + twopi; 593 506 594 in = kOutside; 507 in = kOutside; 595 G4double sphi = fSPhi - halfAngTolerance; 508 G4double sphi = fSPhi - halfAngTolerance; 596 G4double ephi = sphi + fDPhi + kAngToleran 509 G4double ephi = sphi + fDPhi + kAngTolerance; 597 if ((phi0 >= sphi && phi0 <= ephi) || 510 if ((phi0 >= sphi && phi0 <= ephi) || 598 (phi1 >= sphi && phi1 <= ephi) || 511 (phi1 >= sphi && phi1 <= ephi) || 599 (phi2 >= sphi && phi2 <= ephi)) in = 512 (phi2 >= sphi && phi2 <= ephi)) in = kSurface; 600 if (in == kOutside) { return kOutside; } 513 if (in == kOutside) { return kOutside; } 601 514 602 sphi += kAngTolerance; 515 sphi += kAngTolerance; 603 ephi -= kAngTolerance; 516 ephi -= kAngTolerance; 604 if ((phi0 >= sphi && phi0 <= ephi) || 517 if ((phi0 >= sphi && phi0 <= ephi) || 605 (phi1 >= sphi && phi1 <= ephi) || 518 (phi1 >= sphi && phi1 <= ephi) || 606 (phi2 >= sphi && phi2 <= ephi)) in = 519 (phi2 >= sphi && phi2 <= ephi)) in = kInside; 607 if (in == kSurface) { return kSurface; } 520 if (in == kSurface) { return kSurface; } 608 } 521 } 609 522 610 // Check on the Surface for Z 523 // Check on the Surface for Z 611 // 524 // 612 if ((zinLow >= -halfCarTolerance) || (zinHig 525 if ((zinLow >= -halfCarTolerance) || (zinHigh >= -halfCarTolerance)) 613 { 526 { 614 return kSurface; 527 return kSurface; 615 } 528 } 616 529 617 // Check on the Surface for R 530 // Check on the Surface for R 618 // 531 // 619 if (fRMin != 0.0) { tolRMin = fRMin + halfRa << 532 if (fRMin) { tolRMin = fRMin + halfRadTolerance; } 620 else { tolRMin = 0; } 533 else { tolRMin = 0; } 621 tolRMax = fRMax - halfRadTolerance; 534 tolRMax = fRMax - halfRadTolerance; 622 if (((r2 <= tolRMin*tolRMin) || (r2 >= tolRM 535 if (((r2 <= tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax)) && 623 (r2 >= halfRadTolerance*halfRadToleranc 536 (r2 >= halfRadTolerance*halfRadTolerance)) 624 { 537 { 625 return kSurface; 538 return kSurface; 626 } 539 } 627 540 628 return in; 541 return in; 629 } 542 } 630 543 631 ////////////////////////////////////////////// 544 /////////////////////////////////////////////////////////////////////////// 632 // 545 // 633 // Return unit normal of surface closest to p 546 // Return unit normal of surface closest to p 634 // - note if point on z axis, ignore phi divid 547 // - note if point on z axis, ignore phi divided sides 635 // - unsafe if point close to z axis a rmin=0 548 // - unsafe if point close to z axis a rmin=0 - no explicit checks 636 549 637 G4ThreeVector G4CutTubs::SurfaceNormal( const 550 G4ThreeVector G4CutTubs::SurfaceNormal( const G4ThreeVector& p ) const 638 { 551 { 639 G4int noSurfaces = 0; 552 G4int noSurfaces = 0; 640 G4double rho, pPhi; 553 G4double rho, pPhi; 641 G4double distZLow,distZHigh, distRMin, distR 554 G4double distZLow,distZHigh, distRMin, distRMax; 642 G4double distSPhi = kInfinity, distEPhi = kI 555 G4double distSPhi = kInfinity, distEPhi = kInfinity; 643 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 556 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 644 557 645 G4ThreeVector norm, sumnorm(0.,0.,0.); 558 G4ThreeVector norm, sumnorm(0.,0.,0.); 646 G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0); 559 G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0); 647 G4ThreeVector nR, nPs, nPe; 560 G4ThreeVector nR, nPs, nPe; 648 561 649 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()); 562 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()); 650 563 651 distRMin = std::fabs(rho - fRMin); 564 distRMin = std::fabs(rho - fRMin); 652 distRMax = std::fabs(rho - fRMax); 565 distRMax = std::fabs(rho - fRMax); 653 566 654 // dist to Low Cut 567 // dist to Low Cut 655 // 568 // 656 distZLow =std::fabs((p+vZ).dot(fLowNorm)); 569 distZLow =std::fabs((p+vZ).dot(fLowNorm)); 657 << 570 658 // dist to High Cut 571 // dist to High Cut 659 // 572 // 660 distZHigh = std::fabs((p-vZ).dot(fHighNorm)) 573 distZHigh = std::fabs((p-vZ).dot(fHighNorm)); 661 574 662 if (!fPhiFullCutTube) // Protected agains << 575 if (!fPhiFullCutTube) // Protected against (0,0,z) 663 { 576 { 664 if ( rho > halfCarTolerance ) 577 if ( rho > halfCarTolerance ) 665 { 578 { 666 pPhi = std::atan2(p.y(),p.x()); 579 pPhi = std::atan2(p.y(),p.x()); 667 << 580 668 if(pPhi < fSPhi- halfCarTolerance) 581 if(pPhi < fSPhi- halfCarTolerance) { pPhi += twopi; } 669 else if(pPhi > fSPhi+fDPhi+ halfCarToler 582 else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; } 670 583 671 distSPhi = std::fabs(pPhi - fSPhi); << 584 distSPhi = std::fabs(pPhi - fSPhi); 672 distEPhi = std::fabs(pPhi - fSPhi - fDPh << 585 distEPhi = std::fabs(pPhi - fSPhi - fDPhi); 673 } 586 } 674 else if( fRMin == 0.0 ) << 587 else if( !fRMin ) 675 { 588 { 676 distSPhi = 0.; << 589 distSPhi = 0.; 677 distEPhi = 0.; << 590 distEPhi = 0.; 678 } 591 } 679 nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0 << 592 nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0); 680 nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0 << 593 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0); 681 } 594 } 682 if ( rho > halfCarTolerance ) { nR = G4Three 595 if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); } 683 596 684 if( distRMax <= halfCarTolerance ) << 597 if( distRMax <= halfCarTolerance ) 685 { 598 { 686 ++noSurfaces; << 599 noSurfaces ++; 687 sumnorm += nR; 600 sumnorm += nR; 688 } 601 } 689 if( (fRMin != 0.0) && (distRMin <= halfCarTo << 602 if( fRMin && (distRMin <= halfCarTolerance) ) 690 { 603 { 691 ++noSurfaces; << 604 noSurfaces ++; 692 sumnorm -= nR; 605 sumnorm -= nR; 693 } 606 } 694 if( fDPhi < twopi ) << 607 if( fDPhi < twopi ) 695 { 608 { 696 if (distSPhi <= halfAngTolerance) << 609 if (distSPhi <= halfAngTolerance) 697 { 610 { 698 ++noSurfaces; << 611 noSurfaces ++; 699 sumnorm += nPs; 612 sumnorm += nPs; 700 } 613 } 701 if (distEPhi <= halfAngTolerance) << 614 if (distEPhi <= halfAngTolerance) 702 { 615 { 703 ++noSurfaces; << 616 noSurfaces ++; 704 sumnorm += nPe; 617 sumnorm += nPe; 705 } 618 } 706 } 619 } 707 if (distZLow <= halfCarTolerance) << 620 if (distZLow <= halfCarTolerance) 708 { 621 { 709 ++noSurfaces; << 622 noSurfaces ++; 710 sumnorm += fLowNorm; 623 sumnorm += fLowNorm; 711 } 624 } 712 if (distZHigh <= halfCarTolerance) << 625 if (distZHigh <= halfCarTolerance) 713 { 626 { 714 ++noSurfaces; << 627 noSurfaces ++; 715 sumnorm += fHighNorm; 628 sumnorm += fHighNorm; 716 } 629 } 717 if ( noSurfaces == 0 ) 630 if ( noSurfaces == 0 ) 718 { 631 { 719 #ifdef G4CSGDEBUG 632 #ifdef G4CSGDEBUG 720 G4Exception("G4CutTubs::SurfaceNormal(p)", 633 G4Exception("G4CutTubs::SurfaceNormal(p)", "GeomSolids1002", 721 JustWarning, "Point p is not o 634 JustWarning, "Point p is not on surface !?" ); 722 G4int oldprc = G4cout.precision(20); 635 G4int oldprc = G4cout.precision(20); 723 G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<< 636 G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); " 724 << G4endl << G4endl; 637 << G4endl << G4endl; 725 G4cout.precision(oldprc) ; 638 G4cout.precision(oldprc) ; 726 #endif << 639 #endif 727 norm = ApproxSurfaceNormal(p); 640 norm = ApproxSurfaceNormal(p); 728 } 641 } 729 else if ( noSurfaces == 1 ) { norm = sumnor 642 else if ( noSurfaces == 1 ) { norm = sumnorm; } 730 else { norm = sumnor 643 else { norm = sumnorm.unit(); } 731 644 732 return norm; 645 return norm; 733 } 646 } 734 647 735 ////////////////////////////////////////////// 648 ///////////////////////////////////////////////////////////////////////////// 736 // 649 // 737 // Algorithm for SurfaceNormal() following the 650 // Algorithm for SurfaceNormal() following the original specification 738 // for points not on the surface 651 // for points not on the surface 739 652 740 G4ThreeVector G4CutTubs::ApproxSurfaceNormal( 653 G4ThreeVector G4CutTubs::ApproxSurfaceNormal( const G4ThreeVector& p ) const 741 { 654 { 742 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ} 655 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ}; 743 656 744 ENorm side ; 657 ENorm side ; 745 G4ThreeVector norm ; 658 G4ThreeVector norm ; 746 G4double rho, phi ; 659 G4double rho, phi ; 747 G4double distZLow,distZHigh,distZ; 660 G4double distZLow,distZHigh,distZ; 748 G4double distRMin, distRMax, distSPhi, distE 661 G4double distRMin, distRMax, distSPhi, distEPhi, distMin ; 749 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 662 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 750 663 751 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 664 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 752 665 753 distRMin = std::fabs(rho - fRMin) ; 666 distRMin = std::fabs(rho - fRMin) ; 754 distRMax = std::fabs(rho - fRMax) ; 667 distRMax = std::fabs(rho - fRMax) ; 755 668 756 //dist to Low Cut 669 //dist to Low Cut 757 // 670 // 758 distZLow =std::fabs((p+vZ).dot(fLowNorm)); 671 distZLow =std::fabs((p+vZ).dot(fLowNorm)); 759 672 760 //dist to High Cut 673 //dist to High Cut 761 // 674 // 762 distZHigh = std::fabs((p-vZ).dot(fHighNorm)) 675 distZHigh = std::fabs((p-vZ).dot(fHighNorm)); 763 distZ=std::min(distZLow,distZHigh); 676 distZ=std::min(distZLow,distZHigh); 764 677 765 if (distRMin < distRMax) // First minimum 678 if (distRMin < distRMax) // First minimum 766 { 679 { 767 if ( distZ < distRMin ) 680 if ( distZ < distRMin ) 768 { 681 { 769 distMin = distZ ; 682 distMin = distZ ; 770 side = kNZ ; 683 side = kNZ ; 771 } 684 } 772 else 685 else 773 { 686 { 774 distMin = distRMin ; 687 distMin = distRMin ; 775 side = kNRMin ; 688 side = kNRMin ; 776 } 689 } 777 } 690 } 778 else 691 else 779 { 692 { 780 if ( distZ < distRMax ) 693 if ( distZ < distRMax ) 781 { 694 { 782 distMin = distZ ; 695 distMin = distZ ; 783 side = kNZ ; 696 side = kNZ ; 784 } 697 } 785 else 698 else 786 { 699 { 787 distMin = distRMax ; 700 distMin = distRMax ; 788 side = kNRMax ; 701 side = kNRMax ; 789 } 702 } 790 } << 703 } 791 if (!fPhiFullCutTube && (rho != 0.0) ) // << 704 if (!fPhiFullCutTube && rho ) // Protected against (0,0,z) 792 { 705 { 793 phi = std::atan2(p.y(),p.x()) ; 706 phi = std::atan2(p.y(),p.x()) ; 794 707 795 if ( phi < 0 ) { phi += twopi; } 708 if ( phi < 0 ) { phi += twopi; } 796 709 797 if ( fSPhi < 0 ) 710 if ( fSPhi < 0 ) 798 { 711 { 799 distSPhi = std::fabs(phi - (fSPhi + twop 712 distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ; 800 } 713 } 801 else 714 else 802 { 715 { 803 distSPhi = std::fabs(phi - fSPhi)*rho ; 716 distSPhi = std::fabs(phi - fSPhi)*rho ; 804 } 717 } 805 distEPhi = std::fabs(phi - fSPhi - fDPhi)* 718 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ; 806 << 719 807 if (distSPhi < distEPhi) // Find new minim 720 if (distSPhi < distEPhi) // Find new minimum 808 { 721 { 809 if ( distSPhi < distMin ) 722 if ( distSPhi < distMin ) 810 { 723 { 811 side = kNSPhi ; 724 side = kNSPhi ; 812 } 725 } 813 } 726 } 814 else 727 else 815 { 728 { 816 if ( distEPhi < distMin ) 729 if ( distEPhi < distMin ) 817 { 730 { 818 side = kNEPhi ; 731 side = kNEPhi ; 819 } 732 } 820 } 733 } 821 } << 734 } 822 switch ( side ) 735 switch ( side ) 823 { 736 { 824 case kNRMin : // Inner radius 737 case kNRMin : // Inner radius 825 { << 738 { 826 norm = G4ThreeVector(-p.x()/rho, -p.y()/ 739 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ; 827 break ; 740 break ; 828 } 741 } 829 case kNRMax : // Outer radius 742 case kNRMax : // Outer radius 830 { << 743 { 831 norm = G4ThreeVector(p.x()/rho, p.y()/rh 744 norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ; 832 break ; 745 break ; 833 } 746 } 834 case kNZ : // + or - dz 747 case kNZ : // + or - dz 835 { << 748 { 836 if ( distZHigh > distZLow ) { norm = fH 749 if ( distZHigh > distZLow ) { norm = fHighNorm ; } 837 else { norm = fL 750 else { norm = fLowNorm; } 838 break ; 751 break ; 839 } 752 } 840 case kNSPhi: 753 case kNSPhi: 841 { 754 { 842 norm = G4ThreeVector(sinSPhi, -cosSPhi, << 755 norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ; 843 break ; 756 break ; 844 } 757 } 845 case kNEPhi: 758 case kNEPhi: 846 { 759 { 847 norm = G4ThreeVector(-sinEPhi, cosEPhi, << 760 norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ; 848 break; 761 break; 849 } 762 } 850 default: // Should never reach this c 763 default: // Should never reach this case ... 851 { 764 { 852 DumpInfo(); 765 DumpInfo(); 853 G4Exception("G4CutTubs::ApproxSurfaceNor 766 G4Exception("G4CutTubs::ApproxSurfaceNormal()", 854 "GeomSolids1002", JustWarnin 767 "GeomSolids1002", JustWarning, 855 "Undefined side for valid su 768 "Undefined side for valid surface normal to solid."); 856 break ; 769 break ; 857 } << 770 } 858 } << 771 } 859 return norm; 772 return norm; 860 } 773 } 861 774 862 ////////////////////////////////////////////// 775 //////////////////////////////////////////////////////////////////// 863 // 776 // 864 // 777 // 865 // Calculate distance to shape from outside, a 778 // Calculate distance to shape from outside, along normalised vector 866 // - return kInfinity if no intersection, or i 779 // - return kInfinity if no intersection, or intersection distance <= tolerance 867 // 780 // 868 // - Compute the intersection with the z plane << 781 // - Compute the intersection with the z planes 869 // - if at valid r, phi, return 782 // - if at valid r, phi, return 870 // 783 // 871 // -> If point is outer outer radius, compute 784 // -> If point is outer outer radius, compute intersection with rmax 872 // - if at valid phi,z return 785 // - if at valid phi,z return 873 // 786 // 874 // -> Compute intersection with inner radius, 787 // -> Compute intersection with inner radius, taking largest +ve root 875 // - if valid (in z,phi), save intersct 788 // - if valid (in z,phi), save intersction 876 // 789 // 877 // -> If phi segmented, compute intersectio 790 // -> If phi segmented, compute intersections with phi half planes 878 // - return smallest of valid phi inter 791 // - return smallest of valid phi intersections and 879 // inner radius intersection 792 // inner radius intersection 880 // 793 // 881 // NOTE: 794 // NOTE: 882 // - 'if valid' implies tolerant checking of i 795 // - 'if valid' implies tolerant checking of intersection points 883 796 884 G4double G4CutTubs::DistanceToIn( const G4Thre 797 G4double G4CutTubs::DistanceToIn( const G4ThreeVector& p, 885 const G4Thre 798 const G4ThreeVector& v ) const 886 { 799 { 887 G4double snxt = kInfinity ; // snxt = d 800 G4double snxt = kInfinity ; // snxt = default return value 888 G4double tolORMin2, tolIRMax2 ; // 'generou 801 G4double tolORMin2, tolIRMax2 ; // 'generous' radii squared 889 G4double tolORMax2, tolIRMin2; 802 G4double tolORMax2, tolIRMin2; 890 const G4double dRmax = 100.*fRMax; 803 const G4double dRmax = 100.*fRMax; 891 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 804 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 892 << 805 893 // Intersection point variables 806 // Intersection point variables 894 // 807 // 895 G4double Dist, sd=0, xi, yi, zi, rho2, inum, 808 G4double Dist, sd=0, xi, yi, zi, rho2, inum, iden, cosPsi, Comp,calf ; 896 G4double t1, t2, t3, b, c, d ; // Quadra << 809 G4double t1, t2, t3, b, c, d ; // Quadratic solver variables 897 G4double distZLow,distZHigh; 810 G4double distZLow,distZHigh; 898 // Calculate tolerant rmin and rmax 811 // Calculate tolerant rmin and rmax 899 812 900 if (fRMin > kRadTolerance) 813 if (fRMin > kRadTolerance) 901 { 814 { 902 tolORMin2 = (fRMin - halfRadTolerance)*(fR 815 tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ; 903 tolIRMin2 = (fRMin + halfRadTolerance)*(fR 816 tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ; 904 } 817 } 905 else 818 else 906 { 819 { 907 tolORMin2 = 0.0 ; 820 tolORMin2 = 0.0 ; 908 tolIRMin2 = 0.0 ; 821 tolIRMin2 = 0.0 ; 909 } 822 } 910 tolORMax2 = (fRMax + halfRadTolerance)*(fRMa 823 tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ; 911 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMa 824 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ; 912 825 913 // Intersection with ZCut surfaces 826 // Intersection with ZCut surfaces 914 827 915 // dist to Low Cut 828 // dist to Low Cut 916 // 829 // 917 distZLow =(p+vZ).dot(fLowNorm); 830 distZLow =(p+vZ).dot(fLowNorm); 918 831 919 // dist to High Cut 832 // dist to High Cut 920 // 833 // 921 distZHigh = (p-vZ).dot(fHighNorm); 834 distZHigh = (p-vZ).dot(fHighNorm); 922 835 923 if ( distZLow >= -halfCarTolerance ) 836 if ( distZLow >= -halfCarTolerance ) 924 { 837 { 925 calf = v.dot(fLowNorm); 838 calf = v.dot(fLowNorm); 926 if (calf<0) 839 if (calf<0) 927 { 840 { 928 sd = -distZLow/calf; 841 sd = -distZLow/calf; 929 if(sd < 0.0) { sd = 0.0; } 842 if(sd < 0.0) { sd = 0.0; } 930 843 931 xi = p.x() + sd*v.x() ; 844 xi = p.x() + sd*v.x() ; // Intersection coords 932 yi = p.y() + sd*v.y() ; 845 yi = p.y() + sd*v.y() ; 933 rho2 = xi*xi + yi*yi ; 846 rho2 = xi*xi + yi*yi ; 934 847 935 // Check validity of intersection 848 // Check validity of intersection 936 849 937 if ((tolIRMin2 <= rho2) && (rho2 <= tolI 850 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2)) 938 { 851 { 939 if (!fPhiFullCutTube && (rho2 != 0.0)) << 852 if (!fPhiFullCutTube && rho2) 940 { 853 { 941 // Psi = angle made with central (av 854 // Psi = angle made with central (average) phi of shape 942 // 855 // 943 inum = xi*cosCPhi + yi*sinCPhi ; 856 inum = xi*cosCPhi + yi*sinCPhi ; 944 iden = std::sqrt(rho2) ; 857 iden = std::sqrt(rho2) ; 945 cosPsi = inum/iden ; 858 cosPsi = inum/iden ; 946 if (cosPsi >= cosHDPhiIT) { return 859 if (cosPsi >= cosHDPhiIT) { return sd ; } 947 } 860 } 948 else 861 else 949 { 862 { 950 return sd ; 863 return sd ; 951 } 864 } 952 } 865 } 953 } 866 } 954 else 867 else 955 { 868 { 956 if ( sd<halfCarTolerance ) 869 if ( sd<halfCarTolerance ) 957 { 870 { 958 if(calf>=0) { sd=kInfinity; } 871 if(calf>=0) { sd=kInfinity; } 959 return sd ; // On/outside extent, and 872 return sd ; // On/outside extent, and heading away 960 } // -> cannot intersect 873 } // -> cannot intersect 961 } 874 } 962 } 875 } 963 876 964 if(distZHigh >= -halfCarTolerance ) 877 if(distZHigh >= -halfCarTolerance ) 965 { 878 { 966 calf = v.dot(fHighNorm); 879 calf = v.dot(fHighNorm); 967 if (calf<0) 880 if (calf<0) 968 { 881 { 969 sd = -distZHigh/calf; 882 sd = -distZHigh/calf; 970 883 971 if(sd < 0.0) { sd = 0.0; } 884 if(sd < 0.0) { sd = 0.0; } 972 885 973 xi = p.x() + sd*v.x() ; 886 xi = p.x() + sd*v.x() ; // Intersection coords 974 yi = p.y() + sd*v.y() ; 887 yi = p.y() + sd*v.y() ; 975 rho2 = xi*xi + yi*yi ; 888 rho2 = xi*xi + yi*yi ; 976 889 977 // Check validity of intersection 890 // Check validity of intersection 978 891 979 if ((tolIRMin2 <= rho2) && (rho2 <= tolI 892 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2)) 980 { 893 { 981 if (!fPhiFullCutTube && (rho2 != 0.0)) << 894 if (!fPhiFullCutTube && rho2) 982 { 895 { 983 // Psi = angle made with central (av 896 // Psi = angle made with central (average) phi of shape 984 // 897 // 985 inum = xi*cosCPhi + yi*sinCPhi ; 898 inum = xi*cosCPhi + yi*sinCPhi ; 986 iden = std::sqrt(rho2) ; 899 iden = std::sqrt(rho2) ; 987 cosPsi = inum/iden ; 900 cosPsi = inum/iden ; 988 if (cosPsi >= cosHDPhiIT) { return 901 if (cosPsi >= cosHDPhiIT) { return sd ; } 989 } 902 } 990 else 903 else 991 { 904 { 992 return sd ; 905 return sd ; 993 } 906 } 994 } 907 } 995 } 908 } 996 else 909 else 997 { 910 { 998 if ( sd<halfCarTolerance ) 911 if ( sd<halfCarTolerance ) 999 { << 912 { 1000 if(calf>=0) { sd=kInfinity; } 913 if(calf>=0) { sd=kInfinity; } 1001 return sd ; // On/outside extent, an 914 return sd ; // On/outside extent, and heading away 1002 } // -> cannot intersect 915 } // -> cannot intersect 1003 } 916 } 1004 } 917 } 1005 918 1006 // -> Can not intersect z surfaces 919 // -> Can not intersect z surfaces 1007 // 920 // 1008 // Intersection with rmax (possible return) 921 // Intersection with rmax (possible return) and rmin (must also check phi) 1009 // 922 // 1010 // Intersection point (xi,yi,zi) on line x= 923 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. 1011 // 924 // 1012 // Intersects with x^2+y^2=R^2 925 // Intersects with x^2+y^2=R^2 1013 // 926 // 1014 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v 927 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0 1015 // t1 t2 928 // t1 t2 t3 1016 929 1017 t1 = 1.0 - v.z()*v.z() ; 930 t1 = 1.0 - v.z()*v.z() ; 1018 t2 = p.x()*v.x() + p.y()*v.y() ; 931 t2 = p.x()*v.x() + p.y()*v.y() ; 1019 t3 = p.x()*p.x() + p.y()*p.y() ; 932 t3 = p.x()*p.x() + p.y()*p.y() ; 1020 if ( t1 > 0 ) // Check not || to z a 933 if ( t1 > 0 ) // Check not || to z axis 1021 { 934 { 1022 b = t2/t1 ; 935 b = t2/t1 ; 1023 c = t3 - fRMax*fRMax ; 936 c = t3 - fRMax*fRMax ; 1024 << 937 1025 if ((t3 >= tolORMax2) && (t2<0)) // Thi 938 if ((t3 >= tolORMax2) && (t2<0)) // This also handles the tangent case 1026 { 939 { 1027 // Try outer cylinder intersection, c=( 940 // Try outer cylinder intersection, c=(t3-fRMax*fRMax)/t1; 1028 941 1029 c /= t1 ; 942 c /= t1 ; 1030 d = b*b - c ; 943 d = b*b - c ; 1031 944 1032 if (d >= 0) // If real root 945 if (d >= 0) // If real root 1033 { 946 { 1034 sd = c/(-b+std::sqrt(d)); 947 sd = c/(-b+std::sqrt(d)); 1035 if (sd >= 0) // If 'forwards' 948 if (sd >= 0) // If 'forwards' 1036 { 949 { 1037 if ( sd>dRmax ) // Avoid rounding e 950 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on 1038 { // 64 bits systems. 951 { // 64 bits systems. Split long distances and recompute 1039 G4double fTerm = sd-std::fmod(sd, 952 G4double fTerm = sd-std::fmod(sd,dRmax); 1040 sd = fTerm + DistanceToIn(p+fTerm 953 sd = fTerm + DistanceToIn(p+fTerm*v,v); 1041 } << 954 } 1042 // Check z intersection 955 // Check z intersection 1043 // 956 // 1044 zi = p.z() + sd*v.z() ; 957 zi = p.z() + sd*v.z() ; 1045 xi = p.x() + sd*v.x() ; 958 xi = p.x() + sd*v.x() ; 1046 yi = p.y() + sd*v.y() ; 959 yi = p.y() + sd*v.y() ; 1047 if ((-xi*fLowNorm.x()-yi*fLowNorm.y 960 if ((-xi*fLowNorm.x()-yi*fLowNorm.y() 1048 -(zi+fDz)*fLowNorm.z())>-halfC 961 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance) 1049 { 962 { 1050 if ((-xi*fHighNorm.x()-yi*fHighNo 963 if ((-xi*fHighNorm.x()-yi*fHighNorm.y() 1051 +(fDz-zi)*fHighNorm.z())>-ha 964 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance) 1052 { 965 { 1053 // Z ok. Check phi intersection 966 // Z ok. Check phi intersection if reqd 1054 // 967 // 1055 if (fPhiFullCutTube) 968 if (fPhiFullCutTube) 1056 { 969 { 1057 return sd ; 970 return sd ; 1058 } 971 } 1059 else 972 else 1060 { 973 { 1061 xi = p.x() + sd*v.x() ; 974 xi = p.x() + sd*v.x() ; 1062 yi = p.y() + sd*v.y() ; 975 yi = p.y() + sd*v.y() ; 1063 cosPsi = (xi*cosCPhi + yi*sin 976 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ; 1064 if (cosPsi >= cosHDPhiIT) { 977 if (cosPsi >= cosHDPhiIT) { return sd ; } 1065 } 978 } 1066 } // end if std::fabs(zi) 979 } // end if std::fabs(zi) 1067 } 980 } 1068 } // end if (sd>=0) 981 } // end if (sd>=0) 1069 } // end if (d>=0) 982 } // end if (d>=0) 1070 } // end if (r>=fRMax) 983 } // end if (r>=fRMax) 1071 else << 984 else 1072 { 985 { 1073 // Inside outer radius : 986 // Inside outer radius : 1074 // check not inside, and heading throug 987 // check not inside, and heading through tubs (-> 0 to in) 1075 if ((t3 > tolIRMin2) && (t2 < 0) 988 if ((t3 > tolIRMin2) && (t2 < 0) 1076 && (std::fabs(p.z()) <= std::fabs(GetC 989 && (std::fabs(p.z()) <= std::fabs(GetCutZ(p))-halfCarTolerance )) 1077 { 990 { 1078 // Inside both radii, delta r -ve, in 991 // Inside both radii, delta r -ve, inside z extent 1079 992 1080 if (!fPhiFullCutTube) 993 if (!fPhiFullCutTube) 1081 { 994 { 1082 inum = p.x()*cosCPhi + p.y()*sinC 995 inum = p.x()*cosCPhi + p.y()*sinCPhi ; 1083 iden = std::sqrt(t3) ; 996 iden = std::sqrt(t3) ; 1084 cosPsi = inum/iden ; 997 cosPsi = inum/iden ; 1085 if (cosPsi >= cosHDPhiIT) 998 if (cosPsi >= cosHDPhiIT) 1086 { 999 { 1087 // In the old version, the small 1000 // In the old version, the small negative tangent for the point 1088 // on surface was not taken in ac 1001 // on surface was not taken in account, and returning 0.0 ... 1089 // New version: check the tangent << 1002 // New version: check the tangent for the point on surface and 1090 // if no intersection, return kIn 1003 // if no intersection, return kInfinity, if intersection instead 1091 // return sd. 1004 // return sd. 1092 // 1005 // 1093 c = t3-fRMax*fRMax; << 1006 c = t3-fRMax*fRMax; 1094 if ( c<=0.0 ) 1007 if ( c<=0.0 ) 1095 { 1008 { 1096 return 0.0; 1009 return 0.0; 1097 } 1010 } 1098 else 1011 else 1099 { 1012 { 1100 c = c/t1 ; 1013 c = c/t1 ; 1101 d = b*b-c; 1014 d = b*b-c; 1102 if ( d>=0.0 ) 1015 if ( d>=0.0 ) 1103 { 1016 { 1104 snxt = c/(-b+std::sqrt(d)); / 1017 snxt = c/(-b+std::sqrt(d)); // using safe solution 1105 / << 1018 // for quadratic equation 1106 if ( snxt < halfCarTolerance 1019 if ( snxt < halfCarTolerance ) { snxt=0; } 1107 return snxt ; 1020 return snxt ; 1108 } << 1021 } 1109 else 1022 else 1110 { 1023 { 1111 return kInfinity; 1024 return kInfinity; 1112 } 1025 } 1113 } 1026 } 1114 } << 1027 } 1115 } 1028 } 1116 else 1029 else 1117 { << 1030 { 1118 // In the old version, the small ne 1031 // In the old version, the small negative tangent for the point 1119 // on surface was not taken in acco 1032 // on surface was not taken in account, and returning 0.0 ... 1120 // New version: check the tangent f << 1033 // New version: check the tangent for the point on surface and 1121 // if no intersection, return kInfi 1034 // if no intersection, return kInfinity, if intersection instead 1122 // return sd. 1035 // return sd. 1123 // 1036 // 1124 c = t3 - fRMax*fRMax; << 1037 c = t3 - fRMax*fRMax; 1125 if ( c<=0.0 ) 1038 if ( c<=0.0 ) 1126 { 1039 { 1127 return 0.0; 1040 return 0.0; 1128 } 1041 } 1129 else 1042 else 1130 { 1043 { 1131 c = c/t1 ; 1044 c = c/t1 ; 1132 d = b*b-c; 1045 d = b*b-c; 1133 if ( d>=0.0 ) 1046 if ( d>=0.0 ) 1134 { 1047 { 1135 snxt= c/(-b+std::sqrt(d)); // u 1048 snxt= c/(-b+std::sqrt(d)); // using safe solution 1136 // f << 1049 // for quadratic equation 1137 if ( snxt < halfCarTolerance ) 1050 if ( snxt < halfCarTolerance ) { snxt=0; } 1138 return snxt ; 1051 return snxt ; 1139 } << 1052 } 1140 else 1053 else 1141 { 1054 { 1142 return kInfinity; 1055 return kInfinity; 1143 } 1056 } 1144 } 1057 } 1145 } // end if (!fPhiFullCutTube) 1058 } // end if (!fPhiFullCutTube) 1146 } // end if (t3>tolIRMin2) 1059 } // end if (t3>tolIRMin2) 1147 } // end if (Inside Outer Radius) << 1060 } // end if (Inside Outer Radius) 1148 << 1061 1149 if ( fRMin != 0.0 ) // Try inner cylin << 1062 if ( fRMin ) // Try inner cylinder intersection 1150 { 1063 { 1151 c = (t3 - fRMin*fRMin)/t1 ; 1064 c = (t3 - fRMin*fRMin)/t1 ; 1152 d = b*b - c ; 1065 d = b*b - c ; 1153 if ( d >= 0.0 ) // If real root 1066 if ( d >= 0.0 ) // If real root 1154 { 1067 { 1155 // Always want 2nd root - we are outs 1068 // Always want 2nd root - we are outside and know rmax Hit was bad 1156 // - If on surface of rmin also need 1069 // - If on surface of rmin also need farthest root 1157 << 1070 1158 sd =( b > 0. )? c/(-b - std::sqrt(d)) 1071 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d)); 1159 if (sd >= -10*halfCarTolerance) // c 1072 if (sd >= -10*halfCarTolerance) // check forwards 1160 { 1073 { 1161 // Check z intersection 1074 // Check z intersection 1162 // 1075 // 1163 if (sd < 0.0) { sd = 0.0; } 1076 if (sd < 0.0) { sd = 0.0; } 1164 if (sd>dRmax) // Avoid rounding err 1077 if (sd>dRmax) // Avoid rounding errors due to precision issues seen 1165 { // 64 bits systems. S 1078 { // 64 bits systems. Split long distances and recompute 1166 G4double fTerm = sd-std::fmod(sd, 1079 G4double fTerm = sd-std::fmod(sd,dRmax); 1167 sd = fTerm + DistanceToIn(p+fTerm 1080 sd = fTerm + DistanceToIn(p+fTerm*v,v); 1168 } << 1081 } 1169 zi = p.z() + sd*v.z() ; 1082 zi = p.z() + sd*v.z() ; 1170 xi = p.x() + sd*v.x() ; 1083 xi = p.x() + sd*v.x() ; 1171 yi = p.y() + sd*v.y() ; 1084 yi = p.y() + sd*v.y() ; 1172 if ((-xi*fLowNorm.x()-yi*fLowNorm.y 1085 if ((-xi*fLowNorm.x()-yi*fLowNorm.y() 1173 -(zi+fDz)*fLowNorm.z())>-halfC 1086 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance) 1174 { 1087 { 1175 if ((-xi*fHighNorm.x()-yi*fHighNo 1088 if ((-xi*fHighNorm.x()-yi*fHighNorm.y() 1176 +(fDz-zi)*fHighNorm.z())>-ha 1089 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance) 1177 { 1090 { 1178 // Z ok. Check phi 1091 // Z ok. Check phi 1179 // 1092 // 1180 if ( fPhiFullCutTube ) 1093 if ( fPhiFullCutTube ) 1181 { 1094 { 1182 return sd ; << 1095 return sd ; 1183 } 1096 } 1184 else 1097 else 1185 { 1098 { 1186 cosPsi = (xi*cosCPhi + yi*sin 1099 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ; 1187 if (cosPsi >= cosHDPhiIT) 1100 if (cosPsi >= cosHDPhiIT) 1188 { 1101 { 1189 // Good inner radius isect 1102 // Good inner radius isect 1190 // - but earlier phi isect 1103 // - but earlier phi isect still possible 1191 // 1104 // 1192 snxt = sd ; 1105 snxt = sd ; 1193 } 1106 } 1194 } 1107 } 1195 } // end if std::fabs(zi) 1108 } // end if std::fabs(zi) 1196 } 1109 } 1197 } // end if (sd>=0) 1110 } // end if (sd>=0) 1198 } // end if (d>=0) 1111 } // end if (d>=0) 1199 } // end if (fRMin) 1112 } // end if (fRMin) 1200 } 1113 } 1201 1114 1202 // Phi segment intersection 1115 // Phi segment intersection 1203 // 1116 // 1204 // o Tolerant of points inside phi planes b 1117 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5 1205 // 1118 // 1206 // o NOTE: Large duplication of code betwee 1119 // o NOTE: Large duplication of code between sphi & ephi checks 1207 // -> only diffs: sphi -> ephi, Com 1120 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane 1208 // intersection check <=0 -> >=0 1121 // intersection check <=0 -> >=0 1209 // -> use some form of loop Constru 1122 // -> use some form of loop Construct ? 1210 // 1123 // 1211 if ( !fPhiFullCutTube ) 1124 if ( !fPhiFullCutTube ) 1212 { 1125 { 1213 // First phi surface (Starting phi) 1126 // First phi surface (Starting phi) 1214 // 1127 // 1215 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 1128 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 1216 << 1129 1217 if ( Comp < 0 ) // Component in outwards 1130 if ( Comp < 0 ) // Component in outwards normal dirn 1218 { 1131 { 1219 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) 1132 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ; 1220 1133 1221 if ( Dist < halfCarTolerance ) 1134 if ( Dist < halfCarTolerance ) 1222 { 1135 { 1223 sd = Dist/Comp ; 1136 sd = Dist/Comp ; 1224 1137 1225 if (sd < snxt) 1138 if (sd < snxt) 1226 { 1139 { 1227 if ( sd < 0 ) { sd = 0.0; } 1140 if ( sd < 0 ) { sd = 0.0; } 1228 zi = p.z() + sd*v.z() ; 1141 zi = p.z() + sd*v.z() ; 1229 xi = p.x() + sd*v.x() ; 1142 xi = p.x() + sd*v.x() ; 1230 yi = p.y() + sd*v.y() ; 1143 yi = p.y() + sd*v.y() ; 1231 if ((-xi*fLowNorm.x()-yi*fLowNorm.y 1144 if ((-xi*fLowNorm.x()-yi*fLowNorm.y() 1232 -(zi+fDz)*fLowNorm.z())>-halfC 1145 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance) 1233 { 1146 { 1234 if ((-xi*fHighNorm.x()-yi*fHighNo 1147 if ((-xi*fHighNorm.x()-yi*fHighNorm.y() 1235 +(fDz-zi)*fHighNorm.z())>-ha << 1148 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance) 1236 { 1149 { 1237 rho2 = xi*xi + yi*yi ; 1150 rho2 = xi*xi + yi*yi ; 1238 if ( ( (rho2 >= tolIRMin2) && ( 1151 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) ) 1239 || ( (rho2 > tolORMin2) && ( 1152 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2) 1240 && ( v.y()*cosSPhi - v.x()* 1153 && ( v.y()*cosSPhi - v.x()*sinSPhi > 0 ) 1241 && ( v.x()*cosSPhi + v.y()* 1154 && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 ) ) 1242 || ( (rho2 > tolIRMax2) && (r 1155 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2) 1243 && (v.y()*cosSPhi - v.x()*s 1156 && (v.y()*cosSPhi - v.x()*sinSPhi > 0) 1244 && (v.x()*cosSPhi + v.y()*s 1157 && (v.x()*cosSPhi + v.y()*sinSPhi < 0) ) ) 1245 { 1158 { 1246 // z and r intersections good 1159 // z and r intersections good 1247 // - check intersecting with 1160 // - check intersecting with correct half-plane 1248 // 1161 // 1249 if ((yi*cosCPhi-xi*sinCPhi) < 1162 if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = sd; } 1250 } 1163 } 1251 } //two Z conditions 1164 } //two Z conditions 1252 } 1165 } 1253 } 1166 } 1254 } << 1167 } 1255 } 1168 } 1256 << 1169 1257 // Second phi surface (Ending phi) 1170 // Second phi surface (Ending phi) 1258 // 1171 // 1259 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ; 1172 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ; 1260 << 1173 1261 if (Comp < 0 ) // Component in outwards 1174 if (Comp < 0 ) // Component in outwards normal dirn 1262 { 1175 { 1263 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) 1176 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ; 1264 1177 1265 if ( Dist < halfCarTolerance ) 1178 if ( Dist < halfCarTolerance ) 1266 { 1179 { 1267 sd = Dist/Comp ; 1180 sd = Dist/Comp ; 1268 1181 1269 if (sd < snxt) 1182 if (sd < snxt) 1270 { 1183 { 1271 if ( sd < 0 ) { sd = 0; } 1184 if ( sd < 0 ) { sd = 0; } 1272 zi = p.z() + sd*v.z() ; 1185 zi = p.z() + sd*v.z() ; 1273 xi = p.x() + sd*v.x() ; 1186 xi = p.x() + sd*v.x() ; 1274 yi = p.y() + sd*v.y() ; 1187 yi = p.y() + sd*v.y() ; 1275 if ((-xi*fLowNorm.x()-yi*fLowNorm.y 1188 if ((-xi*fLowNorm.x()-yi*fLowNorm.y() 1276 -(zi+fDz)*fLowNorm.z())>-halfC 1189 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance) 1277 { 1190 { 1278 if ((-xi*fHighNorm.x()-yi*fHighNo 1191 if ((-xi*fHighNorm.x()-yi*fHighNorm.y() 1279 +(fDz-zi)*fHighNorm.z())>-ha 1192 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance) 1280 { 1193 { 1281 xi = p.x() + sd*v.x() ; 1194 xi = p.x() + sd*v.x() ; 1282 yi = p.y() + sd*v.y() ; 1195 yi = p.y() + sd*v.y() ; 1283 rho2 = xi*xi + yi*yi ; 1196 rho2 = xi*xi + yi*yi ; 1284 if ( ( (rho2 >= tolIRMin2) && ( 1197 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) ) 1285 || ( (rho2 > tolORMin2) && 1198 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2) 1286 && (v.x()*sinEPhi - v.y() 1199 && (v.x()*sinEPhi - v.y()*cosEPhi > 0) 1287 && (v.x()*cosEPhi + v.y() 1200 && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) ) 1288 || ( (rho2 > tolIRMax2) && 1201 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2) 1289 && (v.x()*sinEPhi - v.y() 1202 && (v.x()*sinEPhi - v.y()*cosEPhi > 0) 1290 && (v.x()*cosEPhi + v.y() 1203 && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) ) 1291 { 1204 { 1292 // z and r intersections good 1205 // z and r intersections good 1293 // - check intersecting with 1206 // - check intersecting with correct half-plane 1294 // 1207 // 1295 if ( (yi*cosCPhi-xi*sinCPhi) 1208 if ( (yi*cosCPhi-xi*sinCPhi) >= -halfCarTolerance ) 1296 { 1209 { 1297 snxt = sd; 1210 snxt = sd; 1298 } 1211 } 1299 } //?? >=-halfCarTolerance 1212 } //?? >=-halfCarTolerance 1300 } 1213 } 1301 } // two Z conditions 1214 } // two Z conditions 1302 } 1215 } 1303 } 1216 } 1304 } // Comp < 0 1217 } // Comp < 0 1305 } // !fPhiFullTube << 1218 } // !fPhiFullTube 1306 if ( snxt<halfCarTolerance ) { snxt=0; } 1219 if ( snxt<halfCarTolerance ) { snxt=0; } 1307 1220 1308 return snxt ; 1221 return snxt ; 1309 } 1222 } 1310 << 1223 1311 ///////////////////////////////////////////// 1224 ////////////////////////////////////////////////////////////////// 1312 // 1225 // 1313 // Calculate distance to shape from outside, 1226 // Calculate distance to shape from outside, along normalised vector 1314 // - return kInfinity if no intersection, or 1227 // - return kInfinity if no intersection, or intersection distance <= tolerance 1315 // 1228 // 1316 // - Compute the intersection with the z plan << 1229 // - Compute the intersection with the z planes 1317 // - if at valid r, phi, return 1230 // - if at valid r, phi, return 1318 // 1231 // 1319 // -> If point is outer outer radius, compute 1232 // -> If point is outer outer radius, compute intersection with rmax 1320 // - if at valid phi,z return 1233 // - if at valid phi,z return 1321 // 1234 // 1322 // -> Compute intersection with inner radius, 1235 // -> Compute intersection with inner radius, taking largest +ve root 1323 // - if valid (in z,phi), save intersc 1236 // - if valid (in z,phi), save intersction 1324 // 1237 // 1325 // -> If phi segmented, compute intersecti 1238 // -> If phi segmented, compute intersections with phi half planes 1326 // - return smallest of valid phi inte 1239 // - return smallest of valid phi intersections and 1327 // inner radius intersection 1240 // inner radius intersection 1328 // 1241 // 1329 // NOTE: 1242 // NOTE: 1330 // - Precalculations for phi trigonometry are 1243 // - Precalculations for phi trigonometry are Done `just in time' 1331 // - `if valid' implies tolerant checking of 1244 // - `if valid' implies tolerant checking of intersection points 1332 // Calculate distance (<= actual) to closes 1245 // Calculate distance (<= actual) to closest surface of shape from outside 1333 // - Calculate distance to z, radial planes 1246 // - Calculate distance to z, radial planes 1334 // - Only to phi planes if outside phi extent 1247 // - Only to phi planes if outside phi extent 1335 // - Return 0 if point inside 1248 // - Return 0 if point inside 1336 1249 1337 G4double G4CutTubs::DistanceToIn( const G4Thr 1250 G4double G4CutTubs::DistanceToIn( const G4ThreeVector& p ) const 1338 { 1251 { 1339 G4double safRMin,safRMax,safZLow,safZHigh,s 1252 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi; 1340 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 1253 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 1341 1254 1342 // Distance to R 1255 // Distance to R 1343 // 1256 // 1344 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) 1257 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 1345 1258 1346 safRMin = fRMin- rho ; 1259 safRMin = fRMin- rho ; 1347 safRMax = rho - fRMax ; 1260 safRMax = rho - fRMax ; 1348 1261 1349 // Distances to ZCut(Low/High) 1262 // Distances to ZCut(Low/High) 1350 1263 1351 // Dist to Low Cut 1264 // Dist to Low Cut 1352 // 1265 // 1353 safZLow = (p+vZ).dot(fLowNorm); 1266 safZLow = (p+vZ).dot(fLowNorm); 1354 1267 1355 // Dist to High Cut 1268 // Dist to High Cut 1356 // 1269 // 1357 safZHigh = (p-vZ).dot(fHighNorm); 1270 safZHigh = (p-vZ).dot(fHighNorm); 1358 1271 1359 safe = std::max(safZLow,safZHigh); 1272 safe = std::max(safZLow,safZHigh); 1360 1273 1361 if ( safRMin > safe ) { safe = safRMin; } 1274 if ( safRMin > safe ) { safe = safRMin; } 1362 if ( safRMax> safe ) { safe = safRMax; } 1275 if ( safRMax> safe ) { safe = safRMax; } 1363 1276 1364 // Distance to Phi 1277 // Distance to Phi 1365 // 1278 // 1366 if ( (!fPhiFullCutTube) && ((rho) != 0.0) ) << 1279 if ( (!fPhiFullCutTube) && (rho) ) 1367 { 1280 { 1368 // Psi=angle from central phi to point 1281 // Psi=angle from central phi to point 1369 // 1282 // 1370 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi) 1283 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ; 1371 << 1284 1372 if ( cosPsi < cosHDPhi ) << 1285 if ( cosPsi < std::cos(fDPhi*0.5) ) 1373 { 1286 { 1374 // Point lies outside phi range 1287 // Point lies outside phi range 1375 << 1288 1376 if ( (p.y()*cosCPhi - p.x()*sinCPhi) < 1289 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 ) 1377 { 1290 { 1378 safePhi = std::fabs(p.x()*sinSPhi - 1291 safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ; 1379 } 1292 } 1380 else 1293 else 1381 { 1294 { 1382 safePhi = std::fabs(p.x()*sinEPhi - 1295 safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ; 1383 } 1296 } 1384 if ( safePhi > safe ) { safe = safePh 1297 if ( safePhi > safe ) { safe = safePhi; } 1385 } 1298 } 1386 } 1299 } 1387 if ( safe < 0 ) { safe = 0; } 1300 if ( safe < 0 ) { safe = 0; } 1388 1301 1389 return safe ; 1302 return safe ; 1390 } 1303 } 1391 1304 1392 ///////////////////////////////////////////// 1305 ////////////////////////////////////////////////////////////////////////////// 1393 // 1306 // 1394 // Calculate distance to surface of shape fro 1307 // Calculate distance to surface of shape from `inside', allowing for tolerance 1395 // - Only Calc rmax intersection if no valid 1308 // - Only Calc rmax intersection if no valid rmin intersection 1396 1309 1397 G4double G4CutTubs::DistanceToOut( const G4Th 1310 G4double G4CutTubs::DistanceToOut( const G4ThreeVector& p, 1398 const G4Th 1311 const G4ThreeVector& v, 1399 const G4bo 1312 const G4bool calcNorm, 1400 G4bo << 1313 G4bool *validNorm, 1401 G4Th << 1314 G4ThreeVector *n ) const 1402 { 1315 { 1403 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,k 1316 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ}; 1404 1317 1405 ESide side=kNull , sider=kNull, sidephi=kNu 1318 ESide side=kNull , sider=kNull, sidephi=kNull ; 1406 G4double snxt=kInfinity, srd=kInfinity,sz=k 1319 G4double snxt=kInfinity, srd=kInfinity,sz=kInfinity, sphi=kInfinity ; 1407 G4double deltaR, t1, t2, t3, b, c, d2, roMi 1320 G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ; 1408 G4double distZLow,distZHigh,calfH,calfL; 1321 G4double distZLow,distZHigh,calfH,calfL; 1409 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 1322 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 1410 << 1323 1411 // Vars for phi intersection: 1324 // Vars for phi intersection: 1412 // 1325 // 1413 G4double pDistS, compS, pDistE, compE, sphi 1326 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ; 1414 << 1327 1415 // Z plane intersection 1328 // Z plane intersection 1416 // Distances to ZCut(Low/High) 1329 // Distances to ZCut(Low/High) 1417 1330 1418 // dist to Low Cut 1331 // dist to Low Cut 1419 // 1332 // 1420 distZLow =(p+vZ).dot(fLowNorm); 1333 distZLow =(p+vZ).dot(fLowNorm); 1421 1334 1422 // dist to High Cut 1335 // dist to High Cut 1423 // 1336 // 1424 distZHigh = (p-vZ).dot(fHighNorm); 1337 distZHigh = (p-vZ).dot(fHighNorm); 1425 1338 1426 calfH = v.dot(fHighNorm); 1339 calfH = v.dot(fHighNorm); 1427 calfL = v.dot(fLowNorm); 1340 calfL = v.dot(fLowNorm); 1428 1341 1429 if (calfH > 0 ) 1342 if (calfH > 0 ) 1430 { 1343 { 1431 if ( distZHigh < halfCarTolerance ) 1344 if ( distZHigh < halfCarTolerance ) 1432 { 1345 { 1433 snxt = -distZHigh/calfH ; 1346 snxt = -distZHigh/calfH ; 1434 side = kPZ ; 1347 side = kPZ ; 1435 } 1348 } 1436 else 1349 else 1437 { 1350 { 1438 if (calcNorm) 1351 if (calcNorm) 1439 { 1352 { 1440 *n = G4ThreeVector(0,0,1) ; 1353 *n = G4ThreeVector(0,0,1) ; 1441 *validNorm = true ; 1354 *validNorm = true ; 1442 } 1355 } 1443 return snxt = 0 ; 1356 return snxt = 0 ; 1444 } 1357 } 1445 } 1358 } 1446 if ( calfL>0) 1359 if ( calfL>0) 1447 { 1360 { 1448 << 1361 1449 if ( distZLow < halfCarTolerance ) 1362 if ( distZLow < halfCarTolerance ) 1450 { 1363 { 1451 sz = -distZLow/calfL ; 1364 sz = -distZLow/calfL ; 1452 if(sz<snxt){ 1365 if(sz<snxt){ 1453 snxt=sz; 1366 snxt=sz; 1454 side = kMZ ; 1367 side = kMZ ; 1455 } 1368 } 1456 << 1369 1457 } 1370 } 1458 else 1371 else 1459 { 1372 { 1460 if (calcNorm) 1373 if (calcNorm) 1461 { 1374 { 1462 *n = G4ThreeVector(0,0,-1) ; 1375 *n = G4ThreeVector(0,0,-1) ; 1463 *validNorm = true ; 1376 *validNorm = true ; 1464 } 1377 } 1465 return snxt = 0.0 ; 1378 return snxt = 0.0 ; 1466 } 1379 } 1467 } 1380 } 1468 if((calfH<=0)&&(calfL<=0)) 1381 if((calfH<=0)&&(calfL<=0)) 1469 { 1382 { 1470 snxt = kInfinity ; // Travel perpendic 1383 snxt = kInfinity ; // Travel perpendicular to z axis 1471 side = kNull; 1384 side = kNull; 1472 } 1385 } 1473 // Radial Intersections 1386 // Radial Intersections 1474 // 1387 // 1475 // Find intersection with cylinders at rmax 1388 // Find intersection with cylinders at rmax/rmin 1476 // Intersection point (xi,yi,zi) on line x= 1389 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. 1477 // 1390 // 1478 // Intersects with x^2+y^2=R^2 1391 // Intersects with x^2+y^2=R^2 1479 // 1392 // 1480 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v 1393 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0 1481 // 1394 // 1482 // t1 t2 1395 // t1 t2 t3 1483 1396 1484 t1 = 1.0 - v.z()*v.z() ; // since v 1397 t1 = 1.0 - v.z()*v.z() ; // since v normalised 1485 t2 = p.x()*v.x() + p.y()*v.y() ; 1398 t2 = p.x()*v.x() + p.y()*v.y() ; 1486 t3 = p.x()*p.x() + p.y()*p.y() ; 1399 t3 = p.x()*p.x() + p.y()*p.y() ; 1487 1400 1488 if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fR 1401 if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fRMax*fRMax; } 1489 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t 1402 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; } // radius^2 on +-fDz 1490 1403 1491 if ( t1 > 0 ) // Check not parallel 1404 if ( t1 > 0 ) // Check not parallel 1492 { 1405 { 1493 // Calculate srd, r exit distance 1406 // Calculate srd, r exit distance 1494 << 1407 1495 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax 1408 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) ) 1496 { 1409 { 1497 // Delta r not negative => leaving via 1410 // Delta r not negative => leaving via rmax 1498 1411 1499 deltaR = t3 - fRMax*fRMax ; 1412 deltaR = t3 - fRMax*fRMax ; 1500 1413 1501 // NOTE: Should use rho-fRMax<-kRadTole 1414 // NOTE: Should use rho-fRMax<-kRadTolerance*0.5 1502 // - avoid sqrt for efficiency 1415 // - avoid sqrt for efficiency 1503 1416 1504 if ( deltaR < -kRadTolerance*fRMax ) 1417 if ( deltaR < -kRadTolerance*fRMax ) 1505 { 1418 { 1506 b = t2/t1 ; 1419 b = t2/t1 ; 1507 c = deltaR/t1 ; 1420 c = deltaR/t1 ; 1508 d2 = b*b-c; 1421 d2 = b*b-c; 1509 if( d2 >= 0 ) { srd = c/( -b - std::s 1422 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); } 1510 else { srd = 0.; } 1423 else { srd = 0.; } 1511 sider = kRMax ; 1424 sider = kRMax ; 1512 } 1425 } 1513 else 1426 else 1514 { 1427 { 1515 // On tolerant boundary & heading out 1428 // On tolerant boundary & heading outwards (or perpendicular to) 1516 // outer radial surface -> leaving im 1429 // outer radial surface -> leaving immediately 1517 1430 1518 if ( calcNorm ) << 1431 if ( calcNorm ) 1519 { 1432 { 1520 *n = G4ThreeVector(p.x()/fR 1433 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ; 1521 *validNorm = true ; 1434 *validNorm = true ; 1522 } 1435 } 1523 return snxt = 0 ; // Leaving by rmax 1436 return snxt = 0 ; // Leaving by rmax immediately 1524 } 1437 } 1525 } << 1438 } 1526 else if ( t2 < 0. ) // i.e. t2 < 0; Poss 1439 else if ( t2 < 0. ) // i.e. t2 < 0; Possible rmin intersection 1527 { 1440 { 1528 roMin2 = t3 - t2*t2/t1 ; // min ro2 of << 1441 roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement 1529 1442 1530 if ( (fRMin != 0.0) && (roMin2 < fRMin* << 1443 if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) ) 1531 { 1444 { 1532 deltaR = t3 - fRMin*fRMin ; 1445 deltaR = t3 - fRMin*fRMin ; 1533 b = t2/t1 ; 1446 b = t2/t1 ; 1534 c = deltaR/t1 ; 1447 c = deltaR/t1 ; 1535 d2 = b*b - c ; 1448 d2 = b*b - c ; 1536 1449 1537 if ( d2 >= 0 ) // Leaving via rmin 1450 if ( d2 >= 0 ) // Leaving via rmin 1538 { 1451 { 1539 // NOTE: SHould use rho-rmin>kRadTo 1452 // NOTE: SHould use rho-rmin>kRadTolerance*0.5 1540 // - avoid sqrt for efficiency 1453 // - avoid sqrt for efficiency 1541 1454 1542 if (deltaR > kRadTolerance*fRMin) 1455 if (deltaR > kRadTolerance*fRMin) 1543 { 1456 { 1544 srd = c/(-b+std::sqrt(d2)); << 1457 srd = c/(-b+std::sqrt(d2)); 1545 sider = kRMin ; 1458 sider = kRMin ; 1546 } 1459 } 1547 else 1460 else 1548 { 1461 { 1549 if ( calcNorm ) { *validNorm = fa 1462 if ( calcNorm ) { *validNorm = false; } // Concave side 1550 return snxt = 0.0; 1463 return snxt = 0.0; 1551 } 1464 } 1552 } 1465 } 1553 else // No rmin intersect -> must 1466 else // No rmin intersect -> must be rmax intersect 1554 { 1467 { 1555 deltaR = t3 - fRMax*fRMax ; 1468 deltaR = t3 - fRMax*fRMax ; 1556 c = deltaR/t1 ; 1469 c = deltaR/t1 ; 1557 d2 = b*b-c; 1470 d2 = b*b-c; 1558 if( d2 >=0. ) 1471 if( d2 >=0. ) 1559 { 1472 { 1560 srd = -b + std::sqrt(d2) ; 1473 srd = -b + std::sqrt(d2) ; 1561 sider = kRMax ; 1474 sider = kRMax ; 1562 } 1475 } 1563 else // Case: On the border+t2<kRad 1476 else // Case: On the border+t2<kRadTolerance 1564 // (v is perpendicular t 1477 // (v is perpendicular to the surface) 1565 { 1478 { 1566 if (calcNorm) 1479 if (calcNorm) 1567 { 1480 { 1568 *n = G4ThreeVector(p.x()/fRMax, 1481 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ; 1569 *validNorm = true ; 1482 *validNorm = true ; 1570 } 1483 } 1571 return snxt = 0.0; 1484 return snxt = 0.0; 1572 } 1485 } 1573 } 1486 } 1574 } 1487 } 1575 else if ( roi2 > fRMax*(fRMax + kRadTol 1488 else if ( roi2 > fRMax*(fRMax + kRadTolerance) ) 1576 // No rmin intersect -> must be rm 1489 // No rmin intersect -> must be rmax intersect 1577 { 1490 { 1578 deltaR = t3 - fRMax*fRMax ; 1491 deltaR = t3 - fRMax*fRMax ; 1579 b = t2/t1 ; 1492 b = t2/t1 ; 1580 c = deltaR/t1; 1493 c = deltaR/t1; 1581 d2 = b*b-c; 1494 d2 = b*b-c; 1582 if( d2 >= 0 ) 1495 if( d2 >= 0 ) 1583 { 1496 { 1584 srd = -b + std::sqrt(d2) ; 1497 srd = -b + std::sqrt(d2) ; 1585 sider = kRMax ; 1498 sider = kRMax ; 1586 } 1499 } 1587 else // Case: On the border+t2<kRadTo 1500 else // Case: On the border+t2<kRadTolerance 1588 // (v is perpendicular to 1501 // (v is perpendicular to the surface) 1589 { 1502 { 1590 if (calcNorm) 1503 if (calcNorm) 1591 { 1504 { 1592 *n = G4ThreeVector(p.x()/fRMax,p. 1505 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ; 1593 *validNorm = true ; 1506 *validNorm = true ; 1594 } 1507 } 1595 return snxt = 0.0; 1508 return snxt = 0.0; 1596 } 1509 } 1597 } 1510 } 1598 } 1511 } 1599 // Phi Intersection 1512 // Phi Intersection 1600 1513 1601 if ( !fPhiFullCutTube ) 1514 if ( !fPhiFullCutTube ) 1602 { 1515 { 1603 // add angle calculation with correctio << 1516 // add angle calculation with correction 1604 // of the difference in domain of atan2 1517 // of the difference in domain of atan2 and Sphi 1605 // 1518 // 1606 vphi = std::atan2(v.y(),v.x()) ; 1519 vphi = std::atan2(v.y(),v.x()) ; 1607 << 1520 1608 if ( vphi < fSPhi - halfAngTolerance ) 1521 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; } 1609 else if ( vphi > fSPhi + fDPhi + halfAn 1522 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; } 1610 1523 1611 1524 1612 if ( (p.x() != 0.0) || (p.y() != 0.0) ) << 1525 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later) 1613 { 1526 { 1614 // pDist -ve when inside 1527 // pDist -ve when inside 1615 1528 1616 pDistS = p.x()*sinSPhi - p.y()*cosSPh 1529 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; 1617 pDistE = -p.x()*sinEPhi + p.y()*cosEP 1530 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ; 1618 1531 1619 // Comp -ve when in direction of outw 1532 // Comp -ve when in direction of outwards normal 1620 1533 1621 compS = -sinSPhi*v.x() + cosSPhi*v. 1534 compS = -sinSPhi*v.x() + cosSPhi*v.y() ; 1622 compE = sinEPhi*v.x() - cosEPhi*v. 1535 compE = sinEPhi*v.x() - cosEPhi*v.y() ; 1623 << 1536 1624 sidephi = kNull; 1537 sidephi = kNull; 1625 << 1538 1626 if( ( (fDPhi <= pi) && ( (pDistS <= h 1539 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance) 1627 && (pDistE <= h 1540 && (pDistE <= halfCarTolerance) ) ) 1628 || ( (fDPhi > pi) && ((pDistS <= h << 1541 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance) 1629 || (pDistE <= << 1542 && (pDistE > halfCarTolerance) ) ) ) 1630 { 1543 { 1631 // Inside both phi *full* planes 1544 // Inside both phi *full* planes 1632 << 1545 1633 if ( compS < 0 ) 1546 if ( compS < 0 ) 1634 { 1547 { 1635 sphi = pDistS/compS ; 1548 sphi = pDistS/compS ; 1636 << 1549 1637 if (sphi >= -halfCarTolerance) 1550 if (sphi >= -halfCarTolerance) 1638 { 1551 { 1639 xi = p.x() + sphi*v.x() ; 1552 xi = p.x() + sphi*v.x() ; 1640 yi = p.y() + sphi*v.y() ; 1553 yi = p.y() + sphi*v.y() ; 1641 << 1554 1642 // Check intersecting with corr 1555 // Check intersecting with correct half-plane 1643 // (if not -> no intersect) 1556 // (if not -> no intersect) 1644 // 1557 // 1645 if( (std::fabs(xi)<=kCarToleran 1558 if( (std::fabs(xi)<=kCarTolerance) 1646 && (std::fabs(yi)<=kCarToleran 1559 && (std::fabs(yi)<=kCarTolerance) ) 1647 { 1560 { 1648 sidephi = kSPhi; 1561 sidephi = kSPhi; 1649 if (((fSPhi-halfAngTolerance) 1562 if (((fSPhi-halfAngTolerance)<=vphi) 1650 &&((fSPhi+fDPhi+halfAngTol 1563 &&((fSPhi+fDPhi+halfAngTolerance)>=vphi)) 1651 { 1564 { 1652 sphi = kInfinity; 1565 sphi = kInfinity; 1653 } 1566 } 1654 } 1567 } 1655 else if ( yi*cosCPhi-xi*sinCPhi 1568 else if ( yi*cosCPhi-xi*sinCPhi >=0 ) 1656 { 1569 { 1657 sphi = kInfinity ; 1570 sphi = kInfinity ; 1658 } 1571 } 1659 else 1572 else 1660 { 1573 { 1661 sidephi = kSPhi ; 1574 sidephi = kSPhi ; 1662 if ( pDistS > -halfCarToleran 1575 if ( pDistS > -halfCarTolerance ) 1663 { 1576 { 1664 sphi = 0.0 ; // Leave by sp 1577 sphi = 0.0 ; // Leave by sphi immediately 1665 } << 1578 } 1666 } << 1579 } 1667 } 1580 } 1668 else 1581 else 1669 { 1582 { 1670 sphi = kInfinity ; 1583 sphi = kInfinity ; 1671 } 1584 } 1672 } 1585 } 1673 else 1586 else 1674 { 1587 { 1675 sphi = kInfinity ; 1588 sphi = kInfinity ; 1676 } 1589 } 1677 1590 1678 if ( compE < 0 ) 1591 if ( compE < 0 ) 1679 { 1592 { 1680 sphi2 = pDistE/compE ; 1593 sphi2 = pDistE/compE ; 1681 << 1594 1682 // Only check further if < starti 1595 // Only check further if < starting phi intersection 1683 // 1596 // 1684 if ( (sphi2 > -halfCarTolerance) 1597 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) ) 1685 { 1598 { 1686 xi = p.x() + sphi2*v.x() ; 1599 xi = p.x() + sphi2*v.x() ; 1687 yi = p.y() + sphi2*v.y() ; 1600 yi = p.y() + sphi2*v.y() ; 1688 << 1601 1689 if ((std::fabs(xi)<=kCarToleran 1602 if ((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance)) 1690 { 1603 { 1691 // Leaving via ending phi 1604 // Leaving via ending phi 1692 // 1605 // 1693 if( (fSPhi-halfAngTolerance > << 1606 if( !((fSPhi-halfAngTolerance <= vphi) 1694 ||(fSPhi+fDPhi+halfAngTo << 1607 &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) ) 1695 { 1608 { 1696 sidephi = kEPhi ; 1609 sidephi = kEPhi ; 1697 if ( pDistE <= -halfCarTole 1610 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; } 1698 else 1611 else { sphi = 0.0 ; } 1699 } 1612 } 1700 } << 1613 } 1701 else // Check intersecting w << 1614 else // Check intersecting with correct half-plane 1702 1615 1703 if ( (yi*cosCPhi-xi*sinCPhi) >= 1616 if ( (yi*cosCPhi-xi*sinCPhi) >= 0) 1704 { 1617 { 1705 // Leaving via ending phi 1618 // Leaving via ending phi 1706 // 1619 // 1707 sidephi = kEPhi ; 1620 sidephi = kEPhi ; 1708 if ( pDistE <= -halfCarTolera 1621 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; } 1709 else 1622 else { sphi = 0.0 ; } 1710 } 1623 } 1711 } 1624 } 1712 } 1625 } 1713 } 1626 } 1714 else 1627 else 1715 { 1628 { 1716 sphi = kInfinity ; 1629 sphi = kInfinity ; 1717 } 1630 } 1718 } 1631 } 1719 else 1632 else 1720 { 1633 { 1721 // On z axis + travel not || to z axi 1634 // On z axis + travel not || to z axis -> if phi of vector direction 1722 // within phi of shape, Step limited 1635 // within phi of shape, Step limited by rmax, else Step =0 1723 << 1636 1724 if ( (fSPhi - halfAngTolerance <= vph 1637 if ( (fSPhi - halfAngTolerance <= vphi) 1725 && (vphi <= fSPhi + fDPhi + halfAn 1638 && (vphi <= fSPhi + fDPhi + halfAngTolerance ) ) 1726 { 1639 { 1727 sphi = kInfinity ; 1640 sphi = kInfinity ; 1728 } 1641 } 1729 else 1642 else 1730 { 1643 { 1731 sidephi = kSPhi ; // arbitrary << 1644 sidephi = kSPhi ; // arbitrary 1732 sphi = 0.0 ; 1645 sphi = 0.0 ; 1733 } 1646 } 1734 } 1647 } 1735 if (sphi < snxt) // Order intersecttio 1648 if (sphi < snxt) // Order intersecttions 1736 { 1649 { 1737 snxt = sphi ; 1650 snxt = sphi ; 1738 side = sidephi ; 1651 side = sidephi ; 1739 } 1652 } 1740 } 1653 } 1741 if (srd < snxt) // Order intersections 1654 if (srd < snxt) // Order intersections 1742 { 1655 { 1743 snxt = srd ; 1656 snxt = srd ; 1744 side = sider ; 1657 side = sider ; 1745 } 1658 } 1746 } 1659 } 1747 if (calcNorm) 1660 if (calcNorm) 1748 { 1661 { 1749 switch(side) 1662 switch(side) 1750 { 1663 { 1751 case kRMax: 1664 case kRMax: 1752 // Note: returned vector not normalis 1665 // Note: returned vector not normalised 1753 // (divide by fRMax for unit vector) 1666 // (divide by fRMax for unit vector) 1754 // 1667 // 1755 xi = p.x() + snxt*v.x() ; 1668 xi = p.x() + snxt*v.x() ; 1756 yi = p.y() + snxt*v.y() ; 1669 yi = p.y() + snxt*v.y() ; 1757 *n = G4ThreeVector(xi/fRMax,yi/fRMax, 1670 *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ; 1758 *validNorm = true ; 1671 *validNorm = true ; 1759 break ; 1672 break ; 1760 1673 1761 case kRMin: 1674 case kRMin: 1762 *validNorm = false ; // Rmin is inco 1675 *validNorm = false ; // Rmin is inconvex 1763 break ; 1676 break ; 1764 1677 1765 case kSPhi: 1678 case kSPhi: 1766 if ( fDPhi <= pi ) 1679 if ( fDPhi <= pi ) 1767 { 1680 { 1768 *n = G4ThreeVector(sinSPhi, 1681 *n = G4ThreeVector(sinSPhi,-cosSPhi,0) ; 1769 *validNorm = true ; 1682 *validNorm = true ; 1770 } 1683 } 1771 else 1684 else 1772 { 1685 { 1773 *validNorm = false ; 1686 *validNorm = false ; 1774 } 1687 } 1775 break ; 1688 break ; 1776 1689 1777 case kEPhi: 1690 case kEPhi: 1778 if (fDPhi <= pi) 1691 if (fDPhi <= pi) 1779 { 1692 { 1780 *n = G4ThreeVector(-sinEPhi,cosEPhi 1693 *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ; 1781 *validNorm = true ; 1694 *validNorm = true ; 1782 } 1695 } 1783 else 1696 else 1784 { 1697 { 1785 *validNorm = false ; 1698 *validNorm = false ; 1786 } 1699 } 1787 break ; 1700 break ; 1788 1701 1789 case kPZ: 1702 case kPZ: 1790 *n = fHighNorm ; 1703 *n = fHighNorm ; 1791 *validNorm = true ; 1704 *validNorm = true ; 1792 break ; 1705 break ; 1793 1706 1794 case kMZ: 1707 case kMZ: 1795 *n = fLowNorm ; 1708 *n = fLowNorm ; 1796 *validNorm = true ; 1709 *validNorm = true ; 1797 break ; 1710 break ; 1798 1711 1799 default: 1712 default: 1800 G4cout << G4endl ; 1713 G4cout << G4endl ; 1801 DumpInfo(); 1714 DumpInfo(); 1802 std::ostringstream message; 1715 std::ostringstream message; 1803 G4long oldprc = message.precision(16) << 1716 G4int oldprc = message.precision(16); 1804 message << "Undefined side for valid 1717 message << "Undefined side for valid surface normal to solid." 1805 << G4endl 1718 << G4endl 1806 << "Position:" << G4endl << 1719 << "Position:" << G4endl << G4endl 1807 << "p.x() = " << p.x()/mm < 1720 << "p.x() = " << p.x()/mm << " mm" << G4endl 1808 << "p.y() = " << p.y()/mm < 1721 << "p.y() = " << p.y()/mm << " mm" << G4endl 1809 << "p.z() = " << p.z()/mm < 1722 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl 1810 << "Direction:" << G4endl << 1723 << "Direction:" << G4endl << G4endl 1811 << "v.x() = " << v.x() << G 1724 << "v.x() = " << v.x() << G4endl 1812 << "v.y() = " << v.y() << G 1725 << "v.y() = " << v.y() << G4endl 1813 << "v.z() = " << v.z() << G 1726 << "v.z() = " << v.z() << G4endl << G4endl 1814 << "Proposed distance :" << G 1727 << "Proposed distance :" << G4endl << G4endl 1815 << "snxt = " << snxt/mm << 1728 << "snxt = " << snxt/mm << " mm" << G4endl ; 1816 message.precision(oldprc) ; 1729 message.precision(oldprc) ; 1817 G4Exception("G4CutTubs::DistanceToOut 1730 G4Exception("G4CutTubs::DistanceToOut(p,v,..)", "GeomSolids1002", 1818 JustWarning, message); 1731 JustWarning, message); 1819 break ; 1732 break ; 1820 } 1733 } 1821 } 1734 } 1822 if ( snxt<halfCarTolerance ) { snxt=0 ; } 1735 if ( snxt<halfCarTolerance ) { snxt=0 ; } 1823 return snxt ; 1736 return snxt ; 1824 } 1737 } 1825 1738 1826 ///////////////////////////////////////////// 1739 ////////////////////////////////////////////////////////////////////////// 1827 // 1740 // 1828 // Calculate distance (<=actual) to closest s 1741 // Calculate distance (<=actual) to closest surface of shape from inside 1829 1742 1830 G4double G4CutTubs::DistanceToOut( const G4Th 1743 G4double G4CutTubs::DistanceToOut( const G4ThreeVector& p ) const 1831 { 1744 { 1832 G4double safRMin,safRMax,safZLow,safZHigh,s 1745 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho; 1833 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 1746 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 1834 1747 1835 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) 1748 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; // Distance to R 1836 1749 1837 safRMin = rho - fRMin ; 1750 safRMin = rho - fRMin ; 1838 safRMax = fRMax - rho ; 1751 safRMax = fRMax - rho ; 1839 1752 1840 // Distances to ZCut(Low/High) 1753 // Distances to ZCut(Low/High) 1841 1754 1842 // Dist to Low Cut 1755 // Dist to Low Cut 1843 // 1756 // 1844 safZLow = std::fabs((p+vZ).dot(fLowNorm)); 1757 safZLow = std::fabs((p+vZ).dot(fLowNorm)); 1845 1758 1846 // Dist to High Cut 1759 // Dist to High Cut 1847 // 1760 // 1848 safZHigh = std::fabs((p-vZ).dot(fHighNorm)) 1761 safZHigh = std::fabs((p-vZ).dot(fHighNorm)); 1849 safe = std::min(safZLow,safZHigh); 1762 safe = std::min(safZLow,safZHigh); 1850 1763 1851 if ( safRMin < safe ) { safe = safRMin; } 1764 if ( safRMin < safe ) { safe = safRMin; } 1852 if ( safRMax< safe ) { safe = safRMax; } 1765 if ( safRMax< safe ) { safe = safRMax; } 1853 1766 1854 // Check if phi divided, Calc distances clo 1767 // Check if phi divided, Calc distances closest phi plane 1855 // 1768 // 1856 if ( !fPhiFullCutTube ) 1769 if ( !fPhiFullCutTube ) 1857 { 1770 { 1858 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 ) 1771 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 ) 1859 { 1772 { 1860 safePhi = -(p.x()*sinSPhi - p.y()*cosSP 1773 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ; 1861 } 1774 } 1862 else 1775 else 1863 { 1776 { 1864 safePhi = (p.x()*sinEPhi - p.y()*cosEPh 1777 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ; 1865 } 1778 } 1866 if (safePhi < safe) { safe = safePhi ; } 1779 if (safePhi < safe) { safe = safePhi ; } 1867 } 1780 } 1868 if ( safe < 0 ) { safe = 0; } 1781 if ( safe < 0 ) { safe = 0; } 1869 1782 1870 return safe ; 1783 return safe ; 1871 } 1784 } 1872 1785 1873 ///////////////////////////////////////////// 1786 ////////////////////////////////////////////////////////////////////////// 1874 // 1787 // 1875 // Stream object contents to an output stream 1788 // Stream object contents to an output stream 1876 1789 1877 G4GeometryType G4CutTubs::GetEntityType() con 1790 G4GeometryType G4CutTubs::GetEntityType() const 1878 { 1791 { 1879 return {"G4CutTubs"}; << 1792 return G4String("G4CutTubs"); 1880 } 1793 } 1881 1794 1882 ///////////////////////////////////////////// 1795 ////////////////////////////////////////////////////////////////////////// 1883 // 1796 // 1884 // Make a clone of the object 1797 // Make a clone of the object 1885 // 1798 // 1886 G4VSolid* G4CutTubs::Clone() const 1799 G4VSolid* G4CutTubs::Clone() const 1887 { 1800 { 1888 return new G4CutTubs(*this); 1801 return new G4CutTubs(*this); 1889 } 1802 } 1890 1803 1891 ///////////////////////////////////////////// 1804 ////////////////////////////////////////////////////////////////////////// 1892 // 1805 // 1893 // Stream object contents to an output stream 1806 // Stream object contents to an output stream 1894 1807 1895 std::ostream& G4CutTubs::StreamInfo( std::ost 1808 std::ostream& G4CutTubs::StreamInfo( std::ostream& os ) const 1896 { 1809 { 1897 G4long oldprc = os.precision(16); << 1810 G4int oldprc = os.precision(16); 1898 os << "------------------------------------ 1811 os << "-----------------------------------------------------------\n" 1899 << " *** Dump for solid - " << GetNam 1812 << " *** Dump for solid - " << GetName() << " ***\n" 1900 << " ================================ 1813 << " ===================================================\n" 1901 << " Solid type: G4CutTubs\n" 1814 << " Solid type: G4CutTubs\n" 1902 << " Parameters: \n" 1815 << " Parameters: \n" 1903 << " inner radius : " << fRMin/mm << 1816 << " inner radius : " << fRMin/mm << " mm \n" 1904 << " outer radius : " << fRMax/mm << 1817 << " outer radius : " << fRMax/mm << " mm \n" 1905 << " half length Z: " << fDz/mm << " 1818 << " half length Z: " << fDz/mm << " mm \n" 1906 << " starting phi : " << fSPhi/degree 1819 << " starting phi : " << fSPhi/degree << " degrees \n" 1907 << " delta phi : " << fDPhi/degree 1820 << " delta phi : " << fDPhi/degree << " degrees \n" 1908 << " low Norm : " << fLowNorm << 1821 << " low Norm : " << fLowNorm << " \n" 1909 << " high Norm : " <<fHighNorm 1822 << " high Norm : " <<fHighNorm << " \n" 1910 << "------------------------------------ 1823 << "-----------------------------------------------------------\n"; 1911 os.precision(oldprc); 1824 os.precision(oldprc); 1912 1825 1913 return os; 1826 return os; 1914 } 1827 } 1915 1828 1916 ///////////////////////////////////////////// 1829 ///////////////////////////////////////////////////////////////////////// 1917 // 1830 // 1918 // GetPointOnSurface 1831 // GetPointOnSurface 1919 1832 1920 G4ThreeVector G4CutTubs::GetPointOnSurface() 1833 G4ThreeVector G4CutTubs::GetPointOnSurface() const 1921 { 1834 { 1922 // Set min and max z << 1835 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose, 1923 if (fZMin == 0. && fZMax == 0.) << 1836 aOne, aTwo, aThr, aFou; >> 1837 G4double rRand; >> 1838 >> 1839 aOne = 2.*fDz*fDPhi*fRMax; >> 1840 aTwo = 2.*fDz*fDPhi*fRMin; >> 1841 aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin); >> 1842 aFou = 2.*fDz*(fRMax-fRMin); >> 1843 >> 1844 phi = G4RandFlat::shoot(fSPhi, fSPhi+fDPhi); >> 1845 cosphi = std::cos(phi); >> 1846 sinphi = std::sin(phi); >> 1847 >> 1848 rRand = GetRadiusInRing(fRMin,fRMax); >> 1849 >> 1850 if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; } >> 1851 >> 1852 chose = G4RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou); >> 1853 >> 1854 if( (chose >=0) && (chose < aOne) ) >> 1855 { >> 1856 xRand = fRMax*cosphi; >> 1857 yRand = fRMax*sinphi; >> 1858 zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)), >> 1859 GetCutZ(G4ThreeVector(xRand,yRand,fDz))); >> 1860 return G4ThreeVector (xRand, yRand, zRand); >> 1861 } >> 1862 else if( (chose >= aOne) && (chose < aOne + aTwo) ) >> 1863 { >> 1864 xRand = fRMin*cosphi; >> 1865 yRand = fRMin*sinphi; >> 1866 zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)), >> 1867 GetCutZ(G4ThreeVector(xRand,yRand,fDz))); >> 1868 return G4ThreeVector (xRand, yRand, zRand); >> 1869 } >> 1870 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) ) >> 1871 { >> 1872 xRand = rRand*cosphi; >> 1873 yRand = rRand*sinphi; >> 1874 zRand = GetCutZ(G4ThreeVector(xRand,yRand,fDz)); >> 1875 return G4ThreeVector (xRand, yRand, zRand); >> 1876 } >> 1877 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) ) >> 1878 { >> 1879 xRand = rRand*cosphi; >> 1880 yRand = rRand*sinphi; >> 1881 zRand = GetCutZ(G4ThreeVector(xRand,yRand,-fDz)); >> 1882 return G4ThreeVector (xRand, yRand, zRand); >> 1883 } >> 1884 else if( (chose >= aOne + aTwo + 2.*aThr) >> 1885 && (chose < aOne + aTwo + 2.*aThr + aFou) ) >> 1886 { >> 1887 xRand = rRand*std::cos(fSPhi); >> 1888 yRand = rRand*std::sin(fSPhi); >> 1889 zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)), >> 1890 GetCutZ(G4ThreeVector(xRand,yRand,fDz))); >> 1891 return G4ThreeVector (xRand, yRand, zRand); >> 1892 } >> 1893 else 1924 { 1894 { 1925 G4AutoLock l(&zminmaxMutex); << 1895 xRand = rRand*std::cos(fSPhi+fDPhi); 1926 G4ThreeVector bmin, bmax; << 1896 yRand = rRand*std::sin(fSPhi+fDPhi); 1927 BoundingLimits(bmin,bmax); << 1897 zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)), 1928 fZMin = bmin.z(); << 1898 GetCutZ(G4ThreeVector(xRand,yRand,fDz))); 1929 fZMax = bmax.z(); << 1899 return G4ThreeVector (xRand, yRand, zRand); 1930 l.unlock(); << 1900 } 1931 } << 1932 << 1933 // Set parameters << 1934 G4double hmax = fZMax - fZMin; << 1935 G4double sphi = fSPhi; << 1936 G4double dphi = fDPhi; << 1937 G4double rmin = fRMin; << 1938 G4double rmax = fRMax; << 1939 G4double rrmax = rmax*rmax; << 1940 G4double rrmin = rmin*rmin; << 1941 << 1942 G4ThreeVector nbot = GetLowNorm(); << 1943 G4ThreeVector ntop = GetHighNorm(); << 1944 << 1945 // Set array of surface areas << 1946 G4double sbase = 0.5*dphi*(rrmax - rrmin); << 1947 G4double sbot = sbase/std::abs(nbot.z()); << 1948 G4double stop = sbase/std::abs(ntop.z()); << 1949 G4double scut = (dphi == twopi) ? 0. : hmax << 1950 G4double ssurf[6] = { scut, scut, sbot, sto << 1951 ssurf[1] += ssurf[0]; << 1952 ssurf[2] += ssurf[1]; << 1953 ssurf[3] += ssurf[2]; << 1954 ssurf[4] += ssurf[3]; << 1955 ssurf[5] += ssurf[4]; << 1956 << 1957 constexpr G4int ntry = 100000; << 1958 for (G4int i=0; i<ntry; ++i) << 1959 { << 1960 // Select surface << 1961 G4double select = ssurf[5]*G4QuickRand(); << 1962 G4int k = 5; << 1963 k -= (G4int)(select <= ssurf[4]); << 1964 k -= (G4int)(select <= ssurf[3]); << 1965 k -= (G4int)(select <= ssurf[2]); << 1966 k -= (G4int)(select <= ssurf[1]); << 1967 k -= (G4int)(select <= ssurf[0]); << 1968 << 1969 // Generate point on selected surface (re << 1970 G4ThreeVector p(0,0,0); << 1971 switch(k) << 1972 { << 1973 case 0: // cut at start phi << 1974 { << 1975 G4double r = rmin + (rmax - rmin)*G4Q << 1976 p.set(r*cosSPhi, r*sinSPhi, fZMin + h << 1977 break; << 1978 } << 1979 case 1: // cut at end phi << 1980 { << 1981 G4double r = rmin + (rmax - rmin)*G4Q << 1982 p.set(r*cosEPhi, r*sinEPhi, fZMin + h << 1983 break; << 1984 } << 1985 case 2: // base at low z << 1986 { << 1987 G4double r = std::sqrt(rrmin + (rrmax << 1988 G4double phi = sphi + dphi*G4QuickRan << 1989 G4double x = r*std::cos(phi); << 1990 G4double y = r*std::sin(phi); << 1991 G4double z = -fDz - (x*nbot.x() + y*n << 1992 return {x, y, z}; << 1993 } << 1994 case 3: // base at high z << 1995 { << 1996 G4double r = std::sqrt(rrmin + (rrmax << 1997 G4double phi = sphi + dphi*G4QuickRan << 1998 G4double x = r*std::cos(phi); << 1999 G4double y = r*std::sin(phi); << 2000 G4double z = fDz - (x*ntop.x() + y*nt << 2001 return {x, y, z}; << 2002 } << 2003 case 4: // external lateral surface << 2004 { << 2005 G4double phi = sphi + dphi*G4QuickRan << 2006 G4double z = fZMin + hmax*G4QuickRand << 2007 G4double x = rmax*std::cos(phi); << 2008 G4double y = rmax*std::sin(phi); << 2009 p.set(x, y, z); << 2010 break; << 2011 } << 2012 case 5: // internal lateral surface << 2013 { << 2014 G4double phi = sphi + dphi*G4QuickRan << 2015 G4double z = fZMin + hmax*G4QuickRand << 2016 G4double x = rmin*std::cos(phi); << 2017 G4double y = rmin*std::sin(phi); << 2018 p.set(x, y, z); << 2019 break; << 2020 } << 2021 } << 2022 if ((ntop.dot(p) - fDz*ntop.z()) > 0.) co << 2023 if ((nbot.dot(p) + fDz*nbot.z()) > 0.) co << 2024 return p; << 2025 } << 2026 // Just in case, if all attempts to generat << 2027 // Normally should never happen << 2028 G4double x = rmax*std::cos(sphi + 0.5*dphi) << 2029 G4double y = rmax*std::sin(sphi + 0.5*dphi) << 2030 G4double z = fDz - (x*ntop.x() + y*ntop.y() << 2031 return {x, y, z}; << 2032 } 1901 } 2033 1902 2034 ///////////////////////////////////////////// 1903 /////////////////////////////////////////////////////////////////////////// 2035 // 1904 // 2036 // Methods for visualisation 1905 // Methods for visualisation 2037 1906 2038 void G4CutTubs::DescribeYourselfTo ( G4VGraph << 1907 void G4CutTubs::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 2039 { 1908 { 2040 scene.AddSolid (*this) ; 1909 scene.AddSolid (*this) ; 2041 } 1910 } 2042 1911 2043 G4Polyhedron* G4CutTubs::CreatePolyhedron () << 1912 G4Polyhedron* G4CutTubs::CreatePolyhedron () const 2044 { 1913 { 2045 typedef G4double G4double3[3]; 1914 typedef G4double G4double3[3]; 2046 typedef G4int G4int4[4]; 1915 typedef G4int G4int4[4]; 2047 1916 2048 auto ph = new G4Polyhedron; << 1917 G4Polyhedron *ph = new G4Polyhedron; 2049 G4Polyhedron *ph1 = new G4PolyhedronTubs (f 1918 G4Polyhedron *ph1 = new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi); 2050 G4int nn=ph1->GetNoVertices(); 1919 G4int nn=ph1->GetNoVertices(); 2051 G4int nf=ph1->GetNoFacets(); 1920 G4int nf=ph1->GetNoFacets(); 2052 auto xyz = new G4double3[nn]; // number of << 1921 G4double3* xyz = new G4double3[nn]; // number of nodes 2053 auto faces = new G4int4[nf] ; // number of << 1922 G4int4* faces = new G4int4[nf] ; // number of faces 2054 1923 2055 for(G4int i=0; i<nn; ++i) << 1924 for(G4int i=0;i<nn;i++) 2056 { 1925 { 2057 xyz[i][0]=ph1->GetVertex(i+1).x(); 1926 xyz[i][0]=ph1->GetVertex(i+1).x(); 2058 xyz[i][1]=ph1->GetVertex(i+1).y(); 1927 xyz[i][1]=ph1->GetVertex(i+1).y(); 2059 G4double tmpZ=ph1->GetVertex(i+1).z(); 1928 G4double tmpZ=ph1->GetVertex(i+1).z(); 2060 if(tmpZ>=fDz-kCarTolerance) 1929 if(tmpZ>=fDz-kCarTolerance) 2061 { 1930 { 2062 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][ 1931 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],fDz)); 2063 } 1932 } 2064 else if(tmpZ<=-fDz+kCarTolerance) 1933 else if(tmpZ<=-fDz+kCarTolerance) 2065 { 1934 { 2066 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][ 1935 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],-fDz)); 2067 } 1936 } 2068 else 1937 else 2069 { 1938 { 2070 xyz[i][2]=tmpZ; 1939 xyz[i][2]=tmpZ; 2071 } 1940 } 2072 } 1941 } 2073 G4int iNodes[4]; 1942 G4int iNodes[4]; 2074 G4int* iEdge = nullptr; << 1943 G4int *iEdge=0; 2075 G4int n; 1944 G4int n; 2076 for(G4int i=0; i<nf ; ++i) << 1945 for(G4int i=0;i<nf;i++) 2077 { 1946 { 2078 ph1->GetFacet(i+1,n,iNodes,iEdge); 1947 ph1->GetFacet(i+1,n,iNodes,iEdge); 2079 for(G4int k=0; k<n; ++k) << 1948 for(G4int k=0;k<n;k++) 2080 { 1949 { 2081 faces[i][k]=iNodes[k]; 1950 faces[i][k]=iNodes[k]; 2082 } 1951 } 2083 for(G4int k=n; k<4; ++k) << 1952 for(G4int k=n;k<4;k++) 2084 { 1953 { 2085 faces[i][k]=0; 1954 faces[i][k]=0; 2086 } 1955 } 2087 } 1956 } 2088 ph->createPolyhedron(nn,nf,xyz,faces); 1957 ph->createPolyhedron(nn,nf,xyz,faces); 2089 1958 2090 delete [] xyz; 1959 delete [] xyz; 2091 delete [] faces; 1960 delete [] faces; 2092 delete ph1; 1961 delete ph1; 2093 1962 2094 return ph; 1963 return ph; 2095 } 1964 } 2096 1965 2097 // Auxilary Methods for Solid 1966 // Auxilary Methods for Solid 2098 << 1967 2099 ///////////////////////////////////////////// << 1968 /////////////////////////////////////////////////////////////////////////// 2100 // << 1969 // Return true if Cutted planes are crossing 2101 // Check set of points on the outer lateral s << 1970 // Check Intersection Points on OX and OY axes 2102 // if the cut planes are crossing inside the << 2103 // << 2104 1971 2105 G4bool G4CutTubs::IsCrossingCutPlanes() const 1972 G4bool G4CutTubs::IsCrossingCutPlanes() const 2106 { 1973 { 2107 constexpr G4int npoints = 30; << 1974 G4double zXLow1,zXLow2,zYLow1,zYLow2; >> 1975 G4double zXHigh1,zXHigh2,zYHigh1,zYHigh2; >> 1976 >> 1977 zXLow1 = GetCutZ(G4ThreeVector(-fRMax, 0,-fDz)); >> 1978 zXLow2 = GetCutZ(G4ThreeVector( fRMax, 0,-fDz)); >> 1979 zYLow1 = GetCutZ(G4ThreeVector( 0,-fRMax,-fDz)); >> 1980 zYLow2 = GetCutZ(G4ThreeVector( 0, fRMax,-fDz)); >> 1981 zXHigh1 = GetCutZ(G4ThreeVector(-fRMax, 0, fDz)); >> 1982 zXHigh2 = GetCutZ(G4ThreeVector( fRMax, 0, fDz)); >> 1983 zYHigh1 = GetCutZ(G4ThreeVector( 0,-fRMax, fDz)); >> 1984 zYHigh2 = GetCutZ(G4ThreeVector( 0, fRMax, fDz)); >> 1985 if ( (zXLow1>zXHigh1) ||(zXLow2>zXHigh2) >> 1986 || (zYLow1>zYHigh1) ||(zYLow2>zYHigh2)) { return true; } 2108 1987 2109 // set values for calculation of h - distan << 2110 // opposite points on bases << 2111 G4ThreeVector nbot = GetLowNorm(); << 2112 G4ThreeVector ntop = GetHighNorm(); << 2113 if (std::abs(nbot.z()) < kCarTolerance) ret << 2114 if (std::abs(ntop.z()) < kCarTolerance) ret << 2115 G4double nx = nbot.x()/nbot.z() - ntop.x()/ << 2116 G4double ny = nbot.y()/nbot.z() - ntop.y()/ << 2117 << 2118 // check points << 2119 G4double cosphi = GetCosStartPhi(); << 2120 G4double sinphi = GetSinStartPhi(); << 2121 G4double delphi = GetDeltaPhiAngle()/npoint << 2122 G4double cosdel = std::cos(delphi); << 2123 G4double sindel = std::sin(delphi); << 2124 G4double hzero = 2.*GetZHalfLength()/GetOut << 2125 for (G4int i=0; i<npoints+1; ++i) << 2126 { << 2127 G4double h = nx*cosphi + ny*sinphi + hzer << 2128 if (h < 0.) return true; << 2129 G4double sintmp = sinphi; << 2130 sinphi = sintmp*cosdel + cosphi*sindel; << 2131 cosphi = cosphi*cosdel - sintmp*sindel; << 2132 } << 2133 return false; 1988 return false; 2134 } 1989 } 2135 1990 2136 ///////////////////////////////////////////// 1991 /////////////////////////////////////////////////////////////////////////// 2137 // 1992 // 2138 // Return real Z coordinate of point on Cutte 1993 // Return real Z coordinate of point on Cutted +/- fDZ plane 2139 1994 2140 G4double G4CutTubs::GetCutZ(const G4ThreeVect 1995 G4double G4CutTubs::GetCutZ(const G4ThreeVector& p) const 2141 { 1996 { 2142 G4double newz = p.z(); // p.z() should be 1997 G4double newz = p.z(); // p.z() should be either +fDz or -fDz 2143 if (p.z()<0) 1998 if (p.z()<0) 2144 { 1999 { 2145 if(fLowNorm.z()!=0.) 2000 if(fLowNorm.z()!=0.) 2146 { 2001 { 2147 newz = -fDz-(p.x()*fLowNorm.x()+p.y()* 2002 newz = -fDz-(p.x()*fLowNorm.x()+p.y()*fLowNorm.y())/fLowNorm.z(); 2148 } 2003 } 2149 } 2004 } 2150 else 2005 else 2151 { 2006 { 2152 if(fHighNorm.z()!=0.) 2007 if(fHighNorm.z()!=0.) 2153 { 2008 { 2154 newz = fDz-(p.x()*fHighNorm.x()+p.y()* 2009 newz = fDz-(p.x()*fHighNorm.x()+p.y()*fHighNorm.y())/fHighNorm.z(); 2155 } 2010 } 2156 } 2011 } 2157 return newz; 2012 return newz; >> 2013 } >> 2014 >> 2015 /////////////////////////////////////////////////////////////////////////// >> 2016 // >> 2017 // Calculate Min and Max Z for CutZ >> 2018 >> 2019 void G4CutTubs::GetMaxMinZ(G4double& zmin,G4double& zmax)const >> 2020 >> 2021 { >> 2022 G4double phiLow = std::atan2(fLowNorm.y(),fLowNorm.x()); >> 2023 G4double phiHigh= std::atan2(fHighNorm.y(),fHighNorm.x()); >> 2024 >> 2025 G4double xc=0, yc=0,z1; >> 2026 G4double z[8]; >> 2027 G4bool in_range_low = false; >> 2028 G4bool in_range_hi = false; >> 2029 >> 2030 G4int i; >> 2031 for (i=0; i<2; i++) >> 2032 { >> 2033 if (phiLow<0) { phiLow+=twopi; } >> 2034 G4double ddp = phiLow-fSPhi; >> 2035 if (ddp<0) { ddp += twopi; } >> 2036 if (ddp <= fDPhi) >> 2037 { >> 2038 xc = fRMin*std::cos(phiLow); >> 2039 yc = fRMin*std::sin(phiLow); >> 2040 z1 = GetCutZ(G4ThreeVector(xc, yc, -fDz)); >> 2041 xc = fRMax*std::cos(phiLow); >> 2042 yc = fRMax*std::sin(phiLow); >> 2043 z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, -fDz))); >> 2044 if (in_range_low) { zmin = std::min(zmin, z1); } >> 2045 else { zmin = z1; } >> 2046 in_range_low = true; >> 2047 } >> 2048 phiLow += pi; >> 2049 if (phiLow>twopi) { phiLow-=twopi; } >> 2050 } >> 2051 for (i=0; i<2; i++) >> 2052 { >> 2053 if (phiHigh<0) { phiHigh+=twopi; } >> 2054 G4double ddp = phiHigh-fSPhi; >> 2055 if (ddp<0) { ddp += twopi; } >> 2056 if (ddp <= fDPhi) >> 2057 { >> 2058 xc = fRMin*std::cos(phiHigh); >> 2059 yc = fRMin*std::sin(phiHigh); >> 2060 z1 = GetCutZ(G4ThreeVector(xc, yc, fDz)); >> 2061 xc = fRMax*std::cos(phiHigh); >> 2062 yc = fRMax*std::sin(phiHigh); >> 2063 z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, fDz))); >> 2064 if (in_range_hi) { zmax = std::min(zmax, z1); } >> 2065 else { zmax = z1; } >> 2066 in_range_hi = true; >> 2067 } >> 2068 phiHigh += pi; >> 2069 if (phiHigh>twopi) { phiHigh-=twopi; } >> 2070 } >> 2071 >> 2072 xc = fRMin*std::cos(fSPhi); >> 2073 yc = fRMin*std::sin(fSPhi); >> 2074 z[0] = GetCutZ(G4ThreeVector(xc, yc, -fDz)); >> 2075 z[4] = GetCutZ(G4ThreeVector(xc, yc, fDz)); >> 2076 >> 2077 xc = fRMin*std::cos(fSPhi+fDPhi); >> 2078 yc = fRMin*std::sin(fSPhi+fDPhi); >> 2079 z[1] = GetCutZ(G4ThreeVector(xc, yc, -fDz)); >> 2080 z[5] = GetCutZ(G4ThreeVector(xc, yc, fDz)); >> 2081 >> 2082 xc = fRMax*std::cos(fSPhi); >> 2083 yc = fRMax*std::sin(fSPhi); >> 2084 z[2] = GetCutZ(G4ThreeVector(xc, yc, -fDz)); >> 2085 z[6] = GetCutZ(G4ThreeVector(xc, yc, fDz)); >> 2086 >> 2087 xc = fRMax*std::cos(fSPhi+fDPhi); >> 2088 yc = fRMax*std::sin(fSPhi+fDPhi); >> 2089 z[3] = GetCutZ(G4ThreeVector(xc, yc, -fDz)); >> 2090 z[7] = GetCutZ(G4ThreeVector(xc, yc, fDz)); >> 2091 >> 2092 // Find min/max >> 2093 >> 2094 z1=z[0]; >> 2095 for (i = 1; i < 4; i++) >> 2096 { >> 2097 if(z[i] < z[i-1])z1=z[i]; >> 2098 } >> 2099 >> 2100 if (in_range_low) >> 2101 { >> 2102 zmin = std::min(zmin, z1); >> 2103 } >> 2104 else >> 2105 { >> 2106 zmin = z1; >> 2107 } >> 2108 z1=z[4]; >> 2109 for (i = 1; i < 4; i++) >> 2110 { >> 2111 if(z[4+i] > z[4+i-1]) { z1=z[4+i]; } >> 2112 } >> 2113 >> 2114 if (in_range_hi) { zmax = std::max(zmax, z1); } >> 2115 else { zmax = z1; } 2158 } 2116 } 2159 #endif 2117 #endif 2160 2118