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