Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // Implementation for G4Cons class 26 // Implementation for G4Cons class 27 // 27 // 28 // ~1994 P.Kent: Created, as main part of t 28 // ~1994 P.Kent: Created, as main part of the geometry prototype 29 // 13.09.96 V.Grichine: Review and final modif 29 // 13.09.96 V.Grichine: Review and final modifications 30 // 03.05.05 V.Grichine: SurfaceNormal(p) accor 30 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal 31 // 12.10.09 T.Nikitina: Added to DistanceToIn( 31 // 12.10.09 T.Nikitina: Added to DistanceToIn(p,v) check on the direction in 32 // case of point on surfa 32 // case of point on surface 33 // 03.10.16 E.Tcherniaev: use G4BoundingEnvelo 33 // 03.10.16 E.Tcherniaev: use G4BoundingEnvelope for CalculateExtent(), 34 // removed CreateRotatedV 34 // removed CreateRotatedVertices() 35 // ------------------------------------------- 35 // -------------------------------------------------------------------- 36 36 37 #include "G4Cons.hh" 37 #include "G4Cons.hh" 38 38 39 #if !defined(G4GEOM_USE_UCONS) 39 #if !defined(G4GEOM_USE_UCONS) 40 40 41 #include "G4GeomTools.hh" 41 #include "G4GeomTools.hh" 42 #include "G4VoxelLimits.hh" 42 #include "G4VoxelLimits.hh" 43 #include "G4AffineTransform.hh" 43 #include "G4AffineTransform.hh" 44 #include "G4BoundingEnvelope.hh" 44 #include "G4BoundingEnvelope.hh" 45 #include "G4GeometryTolerance.hh" 45 #include "G4GeometryTolerance.hh" 46 46 47 #include "G4VPVParameterisation.hh" 47 #include "G4VPVParameterisation.hh" 48 48 49 #include "meshdefs.hh" 49 #include "meshdefs.hh" 50 50 51 #include "Randomize.hh" 51 #include "Randomize.hh" 52 52 53 #include "G4VGraphicsScene.hh" 53 #include "G4VGraphicsScene.hh" 54 54 55 using namespace CLHEP; 55 using namespace CLHEP; 56 56 57 ////////////////////////////////////////////// 57 //////////////////////////////////////////////////////////////////////// 58 // 58 // 59 // Private enums: Not for external use 59 // Private enums: Not for external use 60 60 61 namespace 61 namespace 62 { 62 { 63 // used by DistanceToOut() 63 // used by DistanceToOut() 64 // 64 // 65 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kP 65 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ}; 66 66 67 // used by ApproxSurfaceNormal() 67 // used by ApproxSurfaceNormal() 68 // 68 // 69 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ} 69 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ}; 70 } 70 } 71 71 72 ////////////////////////////////////////////// 72 ////////////////////////////////////////////////////////////////////////// 73 // 73 // 74 // constructor - check parameters, convert ang 74 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI 75 // - note if pDPhi>2PI then reset 75 // - note if pDPhi>2PI then reset to 2PI 76 76 77 G4Cons::G4Cons( const G4String& pName, 77 G4Cons::G4Cons( const G4String& pName, 78 G4double pRmin1, G4doub 78 G4double pRmin1, G4double pRmax1, 79 G4double pRmin2, G4doub 79 G4double pRmin2, G4double pRmax2, 80 G4double pDz, 80 G4double pDz, 81 G4double pSPhi, G4double 81 G4double pSPhi, G4double pDPhi) 82 : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2( 82 : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2), 83 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), 83 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.) 84 { 84 { 85 kRadTolerance = G4GeometryTolerance::GetInst 85 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 86 kAngTolerance = G4GeometryTolerance::GetInst 86 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 87 87 88 halfCarTolerance=kCarTolerance*0.5; 88 halfCarTolerance=kCarTolerance*0.5; 89 halfRadTolerance=kRadTolerance*0.5; 89 halfRadTolerance=kRadTolerance*0.5; 90 halfAngTolerance=kAngTolerance*0.5; 90 halfAngTolerance=kAngTolerance*0.5; 91 91 92 // Check z-len 92 // Check z-len 93 // 93 // 94 if ( pDz < 0 ) 94 if ( pDz < 0 ) 95 { 95 { 96 std::ostringstream message; 96 std::ostringstream message; 97 message << "Invalid Z half-length for Soli 97 message << "Invalid Z half-length for Solid: " << GetName() << G4endl 98 << " hZ = " << pDz; 98 << " hZ = " << pDz; 99 G4Exception("G4Cons::G4Cons()", "GeomSolid 99 G4Exception("G4Cons::G4Cons()", "GeomSolids0002", 100 FatalException, message); 100 FatalException, message); 101 } 101 } 102 102 103 // Check radii 103 // Check radii 104 // 104 // 105 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || 105 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0)) 106 { 106 { 107 std::ostringstream message; 107 std::ostringstream message; 108 message << "Invalid values of radii for So 108 message << "Invalid values of radii for Solid: " << GetName() << G4endl 109 << " pRmin1 = " << pRmin1 < 109 << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2 110 << ", pRmax1 = " << pRmax1 << ", p 110 << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2; 111 G4Exception("G4Cons::G4Cons()", "GeomSolid 111 G4Exception("G4Cons::G4Cons()", "GeomSolids0002", 112 FatalException, message) ; 112 FatalException, message) ; 113 } 113 } 114 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fR 114 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; } 115 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fR 115 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; } 116 116 117 // Check angles 117 // Check angles 118 // 118 // 119 CheckPhiAngles(pSPhi, pDPhi); 119 CheckPhiAngles(pSPhi, pDPhi); 120 } 120 } 121 121 122 ////////////////////////////////////////////// 122 /////////////////////////////////////////////////////////////////////// 123 // 123 // 124 // Fake default constructor - sets only member 124 // Fake default constructor - sets only member data and allocates memory 125 // for usage restri 125 // for usage restricted to object persistency. 126 // 126 // 127 G4Cons::G4Cons( __void__& a ) 127 G4Cons::G4Cons( __void__& a ) 128 : G4CSGSolid(a) << 128 : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.), >> 129 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.), >> 130 fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhi(0.), >> 131 cosHDPhiOT(0.), cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), >> 132 sinEPhi(0.), cosEPhi(0.), >> 133 halfCarTolerance(0.), halfRadTolerance(0.), halfAngTolerance(0.) 129 { 134 { 130 } 135 } 131 136 132 ////////////////////////////////////////////// 137 /////////////////////////////////////////////////////////////////////// 133 // 138 // 134 // Destructor 139 // Destructor 135 140 136 G4Cons::~G4Cons() = default; << 141 G4Cons::~G4Cons() >> 142 { >> 143 } 137 144 138 ////////////////////////////////////////////// 145 ////////////////////////////////////////////////////////////////////////// 139 // 146 // 140 // Copy constructor 147 // Copy constructor 141 148 142 G4Cons::G4Cons(const G4Cons&) = default; << 149 G4Cons::G4Cons(const G4Cons& rhs) >> 150 : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance), >> 151 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2), >> 152 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi), >> 153 fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi), >> 154 cosHDPhi(rhs.cosHDPhi), cosHDPhiOT(rhs.cosHDPhiOT), >> 155 cosHDPhiIT(rhs.cosHDPhiIT), sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), >> 156 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone), >> 157 halfCarTolerance(rhs.halfCarTolerance), >> 158 halfRadTolerance(rhs.halfRadTolerance), >> 159 halfAngTolerance(rhs.halfAngTolerance) >> 160 { >> 161 } 143 162 144 ////////////////////////////////////////////// 163 ////////////////////////////////////////////////////////////////////////// 145 // 164 // 146 // Assignment operator 165 // Assignment operator 147 166 148 G4Cons& G4Cons::operator = (const G4Cons& rhs) 167 G4Cons& G4Cons::operator = (const G4Cons& rhs) 149 { 168 { 150 // Check assignment to self 169 // Check assignment to self 151 // 170 // 152 if (this == &rhs) { return *this; } 171 if (this == &rhs) { return *this; } 153 172 154 // Copy base class data 173 // Copy base class data 155 // 174 // 156 G4CSGSolid::operator=(rhs); 175 G4CSGSolid::operator=(rhs); 157 176 158 // Copy data 177 // Copy data 159 // 178 // 160 kRadTolerance = rhs.kRadTolerance; 179 kRadTolerance = rhs.kRadTolerance; 161 kAngTolerance = rhs.kAngTolerance; 180 kAngTolerance = rhs.kAngTolerance; 162 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2; 181 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2; 163 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2; 182 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2; 164 fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = r 183 fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; 165 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPh 184 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi; cosHDPhi = rhs.cosHDPhi; 166 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r 185 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT; 167 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh 186 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi; 168 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh 187 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi; 169 fPhiFullCone = rhs.fPhiFullCone; 188 fPhiFullCone = rhs.fPhiFullCone; 170 halfCarTolerance = rhs.halfCarTolerance; 189 halfCarTolerance = rhs.halfCarTolerance; 171 halfRadTolerance = rhs.halfRadTolerance; 190 halfRadTolerance = rhs.halfRadTolerance; 172 halfAngTolerance = rhs.halfAngTolerance; 191 halfAngTolerance = rhs.halfAngTolerance; 173 192 174 return *this; 193 return *this; 175 } 194 } 176 195 177 ////////////////////////////////////////////// 196 ///////////////////////////////////////////////////////////////////// 178 // 197 // 179 // Return whether point inside/outside/on surf 198 // Return whether point inside/outside/on surface 180 199 181 EInside G4Cons::Inside(const G4ThreeVector& p) 200 EInside G4Cons::Inside(const G4ThreeVector& p) const 182 { 201 { 183 G4double r2, rl, rh, pPhi, tolRMin, tolRMax; 202 G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ; 184 EInside in; 203 EInside in; 185 204 186 if (std::fabs(p.z()) > fDz + halfCarToleranc 205 if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; } 187 else if(std::fabs(p.z()) >= fDz - halfCarTol 206 else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; } 188 else 207 else { in = kInside; } 189 208 190 r2 = p.x()*p.x() + p.y()*p.y() ; 209 r2 = p.x()*p.x() + p.y()*p.y() ; 191 rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz 210 rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ; 192 rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z 211 rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz; 193 212 194 // rh2 = rh*rh; 213 // rh2 = rh*rh; 195 214 196 tolRMin = rl - halfRadTolerance; 215 tolRMin = rl - halfRadTolerance; 197 if ( tolRMin < 0 ) { tolRMin = 0; } 216 if ( tolRMin < 0 ) { tolRMin = 0; } 198 tolRMax = rh + halfRadTolerance; 217 tolRMax = rh + halfRadTolerance; 199 218 200 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tol 219 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; } 201 220 202 if (rl != 0.0) { tolRMin = rl + halfRadToler << 221 if (rl) { tolRMin = rl + halfRadTolerance; } 203 else { tolRMin = 0.0; } 222 else { tolRMin = 0.0; } 204 tolRMax = rh - halfRadTolerance; 223 tolRMax = rh - halfRadTolerance; 205 224 206 if (in == kInside) // else it's kSurface alr 225 if (in == kInside) // else it's kSurface already 207 { 226 { 208 if ( (r2 < tolRMin*tolRMin) || (r2 >= tol 227 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; } 209 } 228 } 210 if ( !fPhiFullCone && ((p.x() != 0.0) || (p. 229 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) ) 211 { 230 { 212 pPhi = std::atan2(p.y(),p.x()) ; 231 pPhi = std::atan2(p.y(),p.x()) ; 213 232 214 if ( pPhi < fSPhi - halfAngTolerance ) 233 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; } 215 else if ( pPhi > fSPhi + fDPhi + halfAngTo 234 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; } 216 235 217 if ( (pPhi < fSPhi - halfAngTolerance) || 236 if ( (pPhi < fSPhi - halfAngTolerance) || 218 (pPhi > fSPhi + fDPhi + halfAngTolera 237 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; } 219 238 220 else if (in == kInside) // else it's kSur 239 else if (in == kInside) // else it's kSurface anyway already 221 { 240 { 222 if ( (pPhi < fSPhi + halfAngTolerance) 241 if ( (pPhi < fSPhi + halfAngTolerance) || 223 (pPhi > fSPhi + fDPhi - halfAngTol 242 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; } 224 } 243 } 225 } 244 } 226 else if ( !fPhiFullCone ) { in = kSurface; 245 else if ( !fPhiFullCone ) { in = kSurface; } 227 246 228 return in ; 247 return in ; 229 } 248 } 230 249 231 ////////////////////////////////////////////// 250 ///////////////////////////////////////////////////////////////////////// 232 // 251 // 233 // Dispatch to parameterisation for replicatio 252 // Dispatch to parameterisation for replication mechanism dimension 234 // computation & modification. 253 // computation & modification. 235 254 236 void G4Cons::ComputeDimensions( G4VPVPara 255 void G4Cons::ComputeDimensions( G4VPVParameterisation* p, 237 const G4int 256 const G4int n, 238 const G4VPhysic 257 const G4VPhysicalVolume* pRep ) 239 { 258 { 240 p->ComputeDimensions(*this,n,pRep) ; 259 p->ComputeDimensions(*this,n,pRep) ; 241 } 260 } 242 261 243 ////////////////////////////////////////////// 262 /////////////////////////////////////////////////////////////////////// 244 // 263 // 245 // Get bounding box 264 // Get bounding box 246 265 247 void G4Cons::BoundingLimits(G4ThreeVector& pMi 266 void G4Cons::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const 248 { 267 { 249 G4double rmin = std::min(GetInnerRadiusMinus 268 G4double rmin = std::min(GetInnerRadiusMinusZ(),GetInnerRadiusPlusZ()); 250 G4double rmax = std::max(GetOuterRadiusMinus 269 G4double rmax = std::max(GetOuterRadiusMinusZ(),GetOuterRadiusPlusZ()); 251 G4double dz = GetZHalfLength(); 270 G4double dz = GetZHalfLength(); 252 271 253 // Find bounding box 272 // Find bounding box 254 // 273 // 255 if (GetDeltaPhiAngle() < twopi) 274 if (GetDeltaPhiAngle() < twopi) 256 { 275 { 257 G4TwoVector vmin,vmax; 276 G4TwoVector vmin,vmax; 258 G4GeomTools::DiskExtent(rmin,rmax, 277 G4GeomTools::DiskExtent(rmin,rmax, 259 GetSinStartPhi(),G 278 GetSinStartPhi(),GetCosStartPhi(), 260 GetSinEndPhi(),Get 279 GetSinEndPhi(),GetCosEndPhi(), 261 vmin,vmax); 280 vmin,vmax); 262 pMin.set(vmin.x(),vmin.y(),-dz); 281 pMin.set(vmin.x(),vmin.y(),-dz); 263 pMax.set(vmax.x(),vmax.y(), dz); 282 pMax.set(vmax.x(),vmax.y(), dz); 264 } 283 } 265 else 284 else 266 { 285 { 267 pMin.set(-rmax,-rmax,-dz); 286 pMin.set(-rmax,-rmax,-dz); 268 pMax.set( rmax, rmax, dz); 287 pMax.set( rmax, rmax, dz); 269 } 288 } 270 289 271 // Check correctness of the bounding box 290 // Check correctness of the bounding box 272 // 291 // 273 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 292 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 274 { 293 { 275 std::ostringstream message; 294 std::ostringstream message; 276 message << "Bad bounding box (min >= max) 295 message << "Bad bounding box (min >= max) for solid: " 277 << GetName() << " !" 296 << GetName() << " !" 278 << "\npMin = " << pMin 297 << "\npMin = " << pMin 279 << "\npMax = " << pMax; 298 << "\npMax = " << pMax; 280 G4Exception("G4Cons::BoundingLimits()", "G 299 G4Exception("G4Cons::BoundingLimits()", "GeomMgt0001", 281 JustWarning, message); 300 JustWarning, message); 282 DumpInfo(); 301 DumpInfo(); 283 } 302 } 284 } 303 } 285 304 286 ////////////////////////////////////////////// 305 /////////////////////////////////////////////////////////////////////// 287 // 306 // 288 // Calculate extent under transform and specif 307 // Calculate extent under transform and specified limit 289 308 290 G4bool G4Cons::CalculateExtent( const EAxis 309 G4bool G4Cons::CalculateExtent( const EAxis pAxis, 291 const G4VoxelL 310 const G4VoxelLimits& pVoxelLimit, 292 const G4Affine 311 const G4AffineTransform& pTransform, 293 G4double 312 G4double& pMin, 294 G4double 313 G4double& pMax ) const 295 { 314 { 296 G4ThreeVector bmin, bmax; 315 G4ThreeVector bmin, bmax; 297 G4bool exist; 316 G4bool exist; 298 317 299 // Get bounding box 318 // Get bounding box 300 BoundingLimits(bmin,bmax); 319 BoundingLimits(bmin,bmax); 301 320 302 // Check bounding box 321 // Check bounding box 303 G4BoundingEnvelope bbox(bmin,bmax); 322 G4BoundingEnvelope bbox(bmin,bmax); 304 #ifdef G4BBOX_EXTENT 323 #ifdef G4BBOX_EXTENT 305 if (true) return bbox.CalculateExtent(pAxis, 324 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 306 #endif 325 #endif 307 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 326 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax)) 308 { 327 { 309 return exist = pMin < pMax; << 328 return exist = (pMin < pMax) ? true : false; 310 } 329 } 311 330 312 // Get parameters of the solid 331 // Get parameters of the solid 313 G4double rmin1 = GetInnerRadiusMinusZ(); 332 G4double rmin1 = GetInnerRadiusMinusZ(); 314 G4double rmax1 = GetOuterRadiusMinusZ(); 333 G4double rmax1 = GetOuterRadiusMinusZ(); 315 G4double rmin2 = GetInnerRadiusPlusZ(); 334 G4double rmin2 = GetInnerRadiusPlusZ(); 316 G4double rmax2 = GetOuterRadiusPlusZ(); 335 G4double rmax2 = GetOuterRadiusPlusZ(); 317 G4double dz = GetZHalfLength(); 336 G4double dz = GetZHalfLength(); 318 G4double dphi = GetDeltaPhiAngle(); 337 G4double dphi = GetDeltaPhiAngle(); 319 338 320 // Find bounding envelope and calculate exte 339 // Find bounding envelope and calculate extent 321 // 340 // 322 const G4int NSTEPS = 24; // numbe 341 const G4int NSTEPS = 24; // number of steps for whole circle 323 G4double astep = twopi/NSTEPS; // max a 342 G4double astep = twopi/NSTEPS; // max angle for one step 324 G4int ksteps = (dphi <= astep) ? 1 : (G4i 343 G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1; 325 G4double ang = dphi/ksteps; 344 G4double ang = dphi/ksteps; 326 345 327 G4double sinHalf = std::sin(0.5*ang); 346 G4double sinHalf = std::sin(0.5*ang); 328 G4double cosHalf = std::cos(0.5*ang); 347 G4double cosHalf = std::cos(0.5*ang); 329 G4double sinStep = 2.*sinHalf*cosHalf; 348 G4double sinStep = 2.*sinHalf*cosHalf; 330 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 349 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 331 G4double rext1 = rmax1/cosHalf; 350 G4double rext1 = rmax1/cosHalf; 332 G4double rext2 = rmax2/cosHalf; 351 G4double rext2 = rmax2/cosHalf; 333 352 334 // bounding envelope for full cone without h 353 // bounding envelope for full cone without hole consists of two polygons, 335 // in other cases it is a sequence of quadri 354 // in other cases it is a sequence of quadrilaterals 336 if (rmin1 == 0 && rmin2 == 0 && dphi == twop 355 if (rmin1 == 0 && rmin2 == 0 && dphi == twopi) 337 { 356 { 338 G4double sinCur = sinHalf; 357 G4double sinCur = sinHalf; 339 G4double cosCur = cosHalf; 358 G4double cosCur = cosHalf; 340 359 341 G4ThreeVectorList baseA(NSTEPS),baseB(NSTE 360 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS); 342 for (G4int k=0; k<NSTEPS; ++k) 361 for (G4int k=0; k<NSTEPS; ++k) 343 { 362 { 344 baseA[k].set(rext1*cosCur,rext1*sinCur,- 363 baseA[k].set(rext1*cosCur,rext1*sinCur,-dz); 345 baseB[k].set(rext2*cosCur,rext2*sinCur, 364 baseB[k].set(rext2*cosCur,rext2*sinCur, dz); 346 365 347 G4double sinTmp = sinCur; 366 G4double sinTmp = sinCur; 348 sinCur = sinCur*cosStep + cosCur*sinStep 367 sinCur = sinCur*cosStep + cosCur*sinStep; 349 cosCur = cosCur*cosStep - sinTmp*sinStep 368 cosCur = cosCur*cosStep - sinTmp*sinStep; 350 } 369 } 351 std::vector<const G4ThreeVectorList *> pol 370 std::vector<const G4ThreeVectorList *> polygons(2); 352 polygons[0] = &baseA; 371 polygons[0] = &baseA; 353 polygons[1] = &baseB; 372 polygons[1] = &baseB; 354 G4BoundingEnvelope benv(bmin,bmax,polygons 373 G4BoundingEnvelope benv(bmin,bmax,polygons); 355 exist = benv.CalculateExtent(pAxis,pVoxelL 374 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 356 } 375 } 357 else 376 else 358 { 377 { 359 G4double sinStart = GetSinStartPhi(); 378 G4double sinStart = GetSinStartPhi(); 360 G4double cosStart = GetCosStartPhi(); 379 G4double cosStart = GetCosStartPhi(); 361 G4double sinEnd = GetSinEndPhi(); 380 G4double sinEnd = GetSinEndPhi(); 362 G4double cosEnd = GetCosEndPhi(); 381 G4double cosEnd = GetCosEndPhi(); 363 G4double sinCur = sinStart*cosHalf + cos 382 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf; 364 G4double cosCur = cosStart*cosHalf - sin 383 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf; 365 384 366 // set quadrilaterals 385 // set quadrilaterals 367 G4ThreeVectorList pols[NSTEPS+2]; 386 G4ThreeVectorList pols[NSTEPS+2]; 368 for (G4int k=0; k<ksteps+2; ++k) pols[k].r 387 for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4); 369 pols[0][0].set(rmin2*cosStart,rmin2*sinSta 388 pols[0][0].set(rmin2*cosStart,rmin2*sinStart, dz); 370 pols[0][1].set(rmin1*cosStart,rmin1*sinSta 389 pols[0][1].set(rmin1*cosStart,rmin1*sinStart,-dz); 371 pols[0][2].set(rmax1*cosStart,rmax1*sinSta 390 pols[0][2].set(rmax1*cosStart,rmax1*sinStart,-dz); 372 pols[0][3].set(rmax2*cosStart,rmax2*sinSta 391 pols[0][3].set(rmax2*cosStart,rmax2*sinStart, dz); 373 for (G4int k=1; k<ksteps+1; ++k) 392 for (G4int k=1; k<ksteps+1; ++k) 374 { 393 { 375 pols[k][0].set(rmin2*cosCur,rmin2*sinCur 394 pols[k][0].set(rmin2*cosCur,rmin2*sinCur, dz); 376 pols[k][1].set(rmin1*cosCur,rmin1*sinCur 395 pols[k][1].set(rmin1*cosCur,rmin1*sinCur,-dz); 377 pols[k][2].set(rext1*cosCur,rext1*sinCur 396 pols[k][2].set(rext1*cosCur,rext1*sinCur,-dz); 378 pols[k][3].set(rext2*cosCur,rext2*sinCur 397 pols[k][3].set(rext2*cosCur,rext2*sinCur, dz); 379 398 380 G4double sinTmp = sinCur; 399 G4double sinTmp = sinCur; 381 sinCur = sinCur*cosStep + cosCur*sinStep 400 sinCur = sinCur*cosStep + cosCur*sinStep; 382 cosCur = cosCur*cosStep - sinTmp*sinStep 401 cosCur = cosCur*cosStep - sinTmp*sinStep; 383 } 402 } 384 pols[ksteps+1][0].set(rmin2*cosEnd,rmin2*s 403 pols[ksteps+1][0].set(rmin2*cosEnd,rmin2*sinEnd, dz); 385 pols[ksteps+1][1].set(rmin1*cosEnd,rmin1*s 404 pols[ksteps+1][1].set(rmin1*cosEnd,rmin1*sinEnd,-dz); 386 pols[ksteps+1][2].set(rmax1*cosEnd,rmax1*s 405 pols[ksteps+1][2].set(rmax1*cosEnd,rmax1*sinEnd,-dz); 387 pols[ksteps+1][3].set(rmax2*cosEnd,rmax2*s 406 pols[ksteps+1][3].set(rmax2*cosEnd,rmax2*sinEnd, dz); 388 407 389 // set envelope and calculate extent 408 // set envelope and calculate extent 390 std::vector<const G4ThreeVectorList *> pol 409 std::vector<const G4ThreeVectorList *> polygons; 391 polygons.resize(ksteps+2); 410 polygons.resize(ksteps+2); 392 for (G4int k=0; k<ksteps+2; ++k) polygons[ 411 for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k]; 393 G4BoundingEnvelope benv(bmin,bmax,polygons 412 G4BoundingEnvelope benv(bmin,bmax,polygons); 394 exist = benv.CalculateExtent(pAxis,pVoxelL 413 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 395 } 414 } 396 return exist; 415 return exist; 397 } 416 } 398 417 399 ////////////////////////////////////////////// 418 //////////////////////////////////////////////////////////////////////// 400 // 419 // 401 // Return unit normal of surface closest to p 420 // Return unit normal of surface closest to p 402 // - note if point on z axis, ignore phi divid 421 // - note if point on z axis, ignore phi divided sides 403 // - unsafe if point close to z axis a rmin=0 422 // - unsafe if point close to z axis a rmin=0 - no explicit checks 404 423 405 G4ThreeVector G4Cons::SurfaceNormal( const G4T 424 G4ThreeVector G4Cons::SurfaceNormal( const G4ThreeVector& p) const 406 { 425 { 407 G4int noSurfaces = 0; 426 G4int noSurfaces = 0; 408 G4double rho, pPhi; 427 G4double rho, pPhi; 409 G4double distZ, distRMin, distRMax; 428 G4double distZ, distRMin, distRMax; 410 G4double distSPhi = kInfinity, distEPhi = kI 429 G4double distSPhi = kInfinity, distEPhi = kInfinity; 411 G4double tanRMin, secRMin, pRMin, widRMin; 430 G4double tanRMin, secRMin, pRMin, widRMin; 412 G4double tanRMax, secRMax, pRMax, widRMax; 431 G4double tanRMax, secRMax, pRMax, widRMax; 413 432 414 G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = 433 G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.); 415 G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe; 434 G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe; 416 435 417 distZ = std::fabs(std::fabs(p.z()) - fDz); 436 distZ = std::fabs(std::fabs(p.z()) - fDz); 418 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) 437 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()); 419 438 420 tanRMin = (fRmin2 - fRmin1)*0.5/fDz; 439 tanRMin = (fRmin2 - fRmin1)*0.5/fDz; 421 secRMin = std::sqrt(1 + tanRMin*tanRMin); 440 secRMin = std::sqrt(1 + tanRMin*tanRMin); 422 pRMin = rho - p.z()*tanRMin; 441 pRMin = rho - p.z()*tanRMin; 423 widRMin = fRmin2 - fDz*tanRMin; 442 widRMin = fRmin2 - fDz*tanRMin; 424 distRMin = std::fabs(pRMin - widRMin)/secRMi 443 distRMin = std::fabs(pRMin - widRMin)/secRMin; 425 444 426 tanRMax = (fRmax2 - fRmax1)*0.5/fDz; 445 tanRMax = (fRmax2 - fRmax1)*0.5/fDz; 427 secRMax = std::sqrt(1+tanRMax*tanRMax); 446 secRMax = std::sqrt(1+tanRMax*tanRMax); 428 pRMax = rho - p.z()*tanRMax; 447 pRMax = rho - p.z()*tanRMax; 429 widRMax = fRmax2 - fDz*tanRMax; 448 widRMax = fRmax2 - fDz*tanRMax; 430 distRMax = std::fabs(pRMax - widRMax)/secRMa 449 distRMax = std::fabs(pRMax - widRMax)/secRMax; 431 450 432 if (!fPhiFullCone) // Protected against (0 451 if (!fPhiFullCone) // Protected against (0,0,z) 433 { 452 { 434 if ( rho != 0.0 ) << 453 if ( rho ) 435 { 454 { 436 pPhi = std::atan2(p.y(),p.x()); 455 pPhi = std::atan2(p.y(),p.x()); 437 456 438 if (pPhi < fSPhi-halfCarTolerance) 457 if (pPhi < fSPhi-halfCarTolerance) { pPhi += twopi; } 439 else if (pPhi > fSPhi+fDPhi+halfCarToler 458 else if (pPhi > fSPhi+fDPhi+halfCarTolerance) { pPhi -= twopi; } 440 459 441 distSPhi = std::fabs( pPhi - fSPhi ); 460 distSPhi = std::fabs( pPhi - fSPhi ); 442 distEPhi = std::fabs( pPhi - fSPhi - fDP 461 distEPhi = std::fabs( pPhi - fSPhi - fDPhi ); 443 } 462 } 444 else if( ((fRmin1) == 0.0) || ((fRmin2) == << 463 else if( !(fRmin1) || !(fRmin2) ) 445 { 464 { 446 distSPhi = 0.; 465 distSPhi = 0.; 447 distEPhi = 0.; 466 distEPhi = 0.; 448 } 467 } 449 nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0 468 nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0 ); 450 nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0 469 nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0 ); 451 } 470 } 452 if ( rho > halfCarTolerance ) 471 if ( rho > halfCarTolerance ) 453 { 472 { 454 nR = G4ThreeVector(p.x()/rho/secRMax, p.y( 473 nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax); 455 if ((fRmin1 != 0.0) || (fRmin2 != 0.0)) << 474 if (fRmin1 || fRmin2) 456 { 475 { 457 nr = G4ThreeVector(-p.x()/rho/secRMin,-p 476 nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin); 458 } 477 } 459 } 478 } 460 479 461 if( distRMax <= halfCarTolerance ) 480 if( distRMax <= halfCarTolerance ) 462 { 481 { 463 ++noSurfaces; 482 ++noSurfaces; 464 sumnorm += nR; 483 sumnorm += nR; 465 } 484 } 466 if( ((fRmin1 != 0.0) || (fRmin2 != 0.0)) && << 485 if( (fRmin1 || fRmin2) && (distRMin <= halfCarTolerance) ) 467 { 486 { 468 ++noSurfaces; 487 ++noSurfaces; 469 sumnorm += nr; 488 sumnorm += nr; 470 } 489 } 471 if( !fPhiFullCone ) 490 if( !fPhiFullCone ) 472 { 491 { 473 if (distSPhi <= halfAngTolerance) 492 if (distSPhi <= halfAngTolerance) 474 { 493 { 475 ++noSurfaces; 494 ++noSurfaces; 476 sumnorm += nPs; 495 sumnorm += nPs; 477 } 496 } 478 if (distEPhi <= halfAngTolerance) 497 if (distEPhi <= halfAngTolerance) 479 { 498 { 480 ++noSurfaces; 499 ++noSurfaces; 481 sumnorm += nPe; 500 sumnorm += nPe; 482 } 501 } 483 } 502 } 484 if (distZ <= halfCarTolerance) 503 if (distZ <= halfCarTolerance) 485 { 504 { 486 ++noSurfaces; 505 ++noSurfaces; 487 if ( p.z() >= 0.) { sumnorm += nZ; } 506 if ( p.z() >= 0.) { sumnorm += nZ; } 488 else { sumnorm -= nZ; } 507 else { sumnorm -= nZ; } 489 } 508 } 490 if ( noSurfaces == 0 ) 509 if ( noSurfaces == 0 ) 491 { 510 { 492 #ifdef G4CSGDEBUG 511 #ifdef G4CSGDEBUG 493 G4Exception("G4Cons::SurfaceNormal(p)", "G 512 G4Exception("G4Cons::SurfaceNormal(p)", "GeomSolids1002", 494 JustWarning, "Point p is not o 513 JustWarning, "Point p is not on surface !?" ); 495 #endif 514 #endif 496 norm = ApproxSurfaceNormal(p); 515 norm = ApproxSurfaceNormal(p); 497 } 516 } 498 else if ( noSurfaces == 1 ) { norm = sumnor 517 else if ( noSurfaces == 1 ) { norm = sumnorm; } 499 else { norm = sumnor 518 else { norm = sumnorm.unit(); } 500 519 501 return norm ; 520 return norm ; 502 } 521 } 503 522 504 ////////////////////////////////////////////// 523 //////////////////////////////////////////////////////////////////////////// 505 // 524 // 506 // Algorithm for SurfaceNormal() following the 525 // Algorithm for SurfaceNormal() following the original specification 507 // for points not on the surface 526 // for points not on the surface 508 527 509 G4ThreeVector G4Cons::ApproxSurfaceNormal( con 528 G4ThreeVector G4Cons::ApproxSurfaceNormal( const G4ThreeVector& p ) const 510 { 529 { 511 ENorm side ; 530 ENorm side ; 512 G4ThreeVector norm ; 531 G4ThreeVector norm ; 513 G4double rho, phi ; 532 G4double rho, phi ; 514 G4double distZ, distRMin, distRMax, distSPhi 533 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ; 515 G4double tanRMin, secRMin, pRMin, widRMin ; 534 G4double tanRMin, secRMin, pRMin, widRMin ; 516 G4double tanRMax, secRMax, pRMax, widRMax ; 535 G4double tanRMax, secRMax, pRMax, widRMax ; 517 536 518 distZ = std::fabs(std::fabs(p.z()) - fDz) ; 537 distZ = std::fabs(std::fabs(p.z()) - fDz) ; 519 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) 538 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 520 539 521 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 540 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 522 secRMin = std::sqrt(1 + tanRMin*tanRMin) ; 541 secRMin = std::sqrt(1 + tanRMin*tanRMin) ; 523 pRMin = rho - p.z()*tanRMin ; 542 pRMin = rho - p.z()*tanRMin ; 524 widRMin = fRmin2 - fDz*tanRMin ; 543 widRMin = fRmin2 - fDz*tanRMin ; 525 distRMin = std::fabs(pRMin - widRMin)/secRMi 544 distRMin = std::fabs(pRMin - widRMin)/secRMin ; 526 545 527 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 546 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 528 secRMax = std::sqrt(1+tanRMax*tanRMax) ; 547 secRMax = std::sqrt(1+tanRMax*tanRMax) ; 529 pRMax = rho - p.z()*tanRMax ; 548 pRMax = rho - p.z()*tanRMax ; 530 widRMax = fRmax2 - fDz*tanRMax ; 549 widRMax = fRmax2 - fDz*tanRMax ; 531 distRMax = std::fabs(pRMax - widRMax)/secRMa 550 distRMax = std::fabs(pRMax - widRMax)/secRMax ; 532 551 533 if (distRMin < distRMax) // First minimum 552 if (distRMin < distRMax) // First minimum 534 { 553 { 535 if (distZ < distRMin) 554 if (distZ < distRMin) 536 { 555 { 537 distMin = distZ ; 556 distMin = distZ ; 538 side = kNZ ; 557 side = kNZ ; 539 } 558 } 540 else 559 else 541 { 560 { 542 distMin = distRMin ; 561 distMin = distRMin ; 543 side = kNRMin ; 562 side = kNRMin ; 544 } 563 } 545 } 564 } 546 else 565 else 547 { 566 { 548 if (distZ < distRMax) 567 if (distZ < distRMax) 549 { 568 { 550 distMin = distZ ; 569 distMin = distZ ; 551 side = kNZ ; 570 side = kNZ ; 552 } 571 } 553 else 572 else 554 { 573 { 555 distMin = distRMax ; 574 distMin = distRMax ; 556 side = kNRMax ; 575 side = kNRMax ; 557 } 576 } 558 } 577 } 559 if ( !fPhiFullCone && (rho != 0.0) ) // Pro << 578 if ( !fPhiFullCone && rho ) // Protected against (0,0,z) 560 { 579 { 561 phi = std::atan2(p.y(),p.x()) ; 580 phi = std::atan2(p.y(),p.x()) ; 562 581 563 if (phi < 0) { phi += twopi; } 582 if (phi < 0) { phi += twopi; } 564 583 565 if (fSPhi < 0) { distSPhi = std::fabs(phi 584 if (fSPhi < 0) { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; } 566 else { distSPhi = std::fabs(phi 585 else { distSPhi = std::fabs(phi - fSPhi)*rho; } 567 586 568 distEPhi = std::fabs(phi - fSPhi - fDPhi)* 587 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ; 569 588 570 // Find new minimum 589 // Find new minimum 571 590 572 if (distSPhi < distEPhi) 591 if (distSPhi < distEPhi) 573 { 592 { 574 if (distSPhi < distMin) { side = kNSPhi 593 if (distSPhi < distMin) { side = kNSPhi; } 575 } 594 } 576 else 595 else 577 { 596 { 578 if (distEPhi < distMin) { side = kNEPhi 597 if (distEPhi < distMin) { side = kNEPhi; } 579 } 598 } 580 } 599 } 581 switch (side) 600 switch (side) 582 { 601 { 583 case kNRMin: // Inner radius 602 case kNRMin: // Inner radius 584 { 603 { 585 rho *= secRMin ; 604 rho *= secRMin ; 586 norm = G4ThreeVector(-p.x()/rho, -p.y()/ 605 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ; 587 break ; 606 break ; 588 } 607 } 589 case kNRMax: // Outer radius 608 case kNRMax: // Outer radius 590 { 609 { 591 rho *= secRMax ; 610 rho *= secRMax ; 592 norm = G4ThreeVector(p.x()/rho, p.y()/rh 611 norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ; 593 break ; 612 break ; 594 } 613 } 595 case kNZ: // +/- dz 614 case kNZ: // +/- dz 596 { 615 { 597 if (p.z() > 0) { norm = G4ThreeVector(0 616 if (p.z() > 0) { norm = G4ThreeVector(0,0,1); } 598 else { norm = G4ThreeVector(0 617 else { norm = G4ThreeVector(0,0,-1); } 599 break ; 618 break ; 600 } 619 } 601 case kNSPhi: 620 case kNSPhi: 602 { 621 { 603 norm = G4ThreeVector(sinSPhi, -cosSPhi, 622 norm = G4ThreeVector(sinSPhi, -cosSPhi, 0) ; 604 break ; 623 break ; 605 } 624 } 606 case kNEPhi: 625 case kNEPhi: 607 { 626 { 608 norm = G4ThreeVector(-sinEPhi, cosEPhi, 627 norm = G4ThreeVector(-sinEPhi, cosEPhi, 0) ; 609 break ; 628 break ; 610 } 629 } 611 default: // Should never reach th 630 default: // Should never reach this case... 612 { 631 { 613 DumpInfo(); 632 DumpInfo(); 614 G4Exception("G4Cons::ApproxSurfaceNormal 633 G4Exception("G4Cons::ApproxSurfaceNormal()", 615 "GeomSolids1002", JustWarnin 634 "GeomSolids1002", JustWarning, 616 "Undefined side for valid su 635 "Undefined side for valid surface normal to solid."); 617 break ; 636 break ; 618 } 637 } 619 } 638 } 620 return norm ; 639 return norm ; 621 } 640 } 622 641 623 ////////////////////////////////////////////// 642 //////////////////////////////////////////////////////////////////////// 624 // 643 // 625 // Calculate distance to shape from outside, a 644 // Calculate distance to shape from outside, along normalised vector 626 // - return kInfinity if no intersection, or i 645 // - return kInfinity if no intersection, or intersection distance <= tolerance 627 // 646 // 628 // - Compute the intersection with the z plane 647 // - Compute the intersection with the z planes 629 // - if at valid r, phi, return 648 // - if at valid r, phi, return 630 // 649 // 631 // -> If point is outside cone, compute inters 650 // -> If point is outside cone, compute intersection with rmax1*0.5 632 // - if at valid phi,z return 651 // - if at valid phi,z return 633 // - if inside outer cone, handle case 652 // - if inside outer cone, handle case when on tolerant outer cone 634 // boundary and heading inwards(->0 t 653 // boundary and heading inwards(->0 to in) 635 // 654 // 636 // -> Compute intersection with inner cone, ta 655 // -> Compute intersection with inner cone, taking largest +ve root 637 // - if valid (in z,phi), save intersct 656 // - if valid (in z,phi), save intersction 638 // 657 // 639 // -> If phi segmented, compute intersectio 658 // -> If phi segmented, compute intersections with phi half planes 640 // - return smallest of valid phi inter 659 // - return smallest of valid phi intersections and 641 // inner radius intersection 660 // inner radius intersection 642 // 661 // 643 // NOTE: 662 // NOTE: 644 // - `if valid' implies tolerant checking of i 663 // - `if valid' implies tolerant checking of intersection points 645 // - z, phi intersection from Tubs 664 // - z, phi intersection from Tubs 646 665 647 G4double G4Cons::DistanceToIn( const G4ThreeVe 666 G4double G4Cons::DistanceToIn( const G4ThreeVector& p, 648 const G4ThreeVe 667 const G4ThreeVector& v ) const 649 { 668 { 650 G4double snxt = kInfinity ; // snxt = d 669 G4double snxt = kInfinity ; // snxt = default return value 651 const G4double dRmax = 50*(fRmax1+fRmax2);// 670 const G4double dRmax = 50*(fRmax1+fRmax2);// 100*(Rmax1+Rmax2)/2. 652 671 653 G4double tanRMax,secRMax,rMaxAv,rMaxOAv ; / 672 G4double tanRMax,secRMax,rMaxAv,rMaxOAv ; // Data for cones 654 G4double tanRMin,secRMin,rMinAv,rMinOAv ; 673 G4double tanRMin,secRMin,rMinAv,rMinOAv ; 655 G4double rout,rin ; 674 G4double rout,rin ; 656 675 657 G4double tolORMin,tolORMin2,tolIRMin,tolIRMi 676 G4double tolORMin,tolORMin2,tolIRMin,tolIRMin2 ; // `generous' radii squared 658 G4double tolORMax2,tolIRMax,tolIRMax2 ; 677 G4double tolORMax2,tolIRMax,tolIRMax2 ; 659 G4double tolODz,tolIDz ; 678 G4double tolODz,tolIDz ; 660 679 661 G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2, 680 G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ; // Intersection point vars 662 681 663 G4double t1,t2,t3,b,c,d ; // Quadratic so 682 G4double t1,t2,t3,b,c,d ; // Quadratic solver variables 664 G4double nt1,nt2,nt3 ; 683 G4double nt1,nt2,nt3 ; 665 G4double Comp ; 684 G4double Comp ; 666 685 667 G4ThreeVector Normal; 686 G4ThreeVector Normal; 668 687 669 // Cone Precalcs 688 // Cone Precalcs 670 689 671 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 690 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 672 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ; 691 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ; 673 rMinAv = (fRmin1 + fRmin2)*0.5 ; 692 rMinAv = (fRmin1 + fRmin2)*0.5 ; 674 693 675 if (rMinAv > halfRadTolerance) 694 if (rMinAv > halfRadTolerance) 676 { 695 { 677 rMinOAv = rMinAv - halfRadTolerance ; 696 rMinOAv = rMinAv - halfRadTolerance ; 678 } 697 } 679 else 698 else 680 { 699 { 681 rMinOAv = 0.0 ; 700 rMinOAv = 0.0 ; 682 } 701 } 683 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 702 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 684 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 703 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 685 rMaxAv = (fRmax1 + fRmax2)*0.5 ; 704 rMaxAv = (fRmax1 + fRmax2)*0.5 ; 686 rMaxOAv = rMaxAv + halfRadTolerance ; 705 rMaxOAv = rMaxAv + halfRadTolerance ; 687 706 688 // Intersection with z-surfaces 707 // Intersection with z-surfaces 689 708 690 tolIDz = fDz - halfCarTolerance ; 709 tolIDz = fDz - halfCarTolerance ; 691 tolODz = fDz + halfCarTolerance ; 710 tolODz = fDz + halfCarTolerance ; 692 711 693 if (std::fabs(p.z()) >= tolIDz) 712 if (std::fabs(p.z()) >= tolIDz) 694 { 713 { 695 if ( p.z()*v.z() < 0 ) // at +Z going i 714 if ( p.z()*v.z() < 0 ) // at +Z going in -Z or visa versa 696 { 715 { 697 sd = (std::fabs(p.z()) - fDz)/std::fabs( 716 sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance 698 717 699 if( sd < 0.0 ) { sd = 0.0; } 718 if( sd < 0.0 ) { sd = 0.0; } // negative dist -> zero 700 719 701 xi = p.x() + sd*v.x() ; // Intersecti 720 xi = p.x() + sd*v.x() ; // Intersection coords 702 yi = p.y() + sd*v.y() ; 721 yi = p.y() + sd*v.y() ; 703 rhoi2 = xi*xi + yi*yi ; 722 rhoi2 = xi*xi + yi*yi ; 704 723 705 // Check validity of intersection 724 // Check validity of intersection 706 // Calculate (outer) tolerant radi^2 at 725 // Calculate (outer) tolerant radi^2 at intersecion 707 726 708 if (v.z() > 0) 727 if (v.z() > 0) 709 { 728 { 710 tolORMin = fRmin1 - halfRadTolerance* 729 tolORMin = fRmin1 - halfRadTolerance*secRMin ; 711 tolIRMin = fRmin1 + halfRadTolerance* 730 tolIRMin = fRmin1 + halfRadTolerance*secRMin ; 712 tolIRMax = fRmax1 - halfRadTolerance* 731 tolIRMax = fRmax1 - halfRadTolerance*secRMin ; 713 // tolORMax2 = (fRmax1 + halfRadTolera 732 // tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)* 714 // (fRmax1 + halfRadTolera 733 // (fRmax1 + halfRadTolerance*secRMax) ; 715 } 734 } 716 else 735 else 717 { 736 { 718 tolORMin = fRmin2 - halfRadTolerance* 737 tolORMin = fRmin2 - halfRadTolerance*secRMin ; 719 tolIRMin = fRmin2 + halfRadTolerance* 738 tolIRMin = fRmin2 + halfRadTolerance*secRMin ; 720 tolIRMax = fRmax2 - halfRadTolerance* 739 tolIRMax = fRmax2 - halfRadTolerance*secRMin ; 721 // tolORMax2 = (fRmax2 + halfRadTolera 740 // tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)* 722 // (fRmax2 + halfRadTolera 741 // (fRmax2 + halfRadTolerance*secRMax) ; 723 } 742 } 724 if ( tolORMin > 0 ) 743 if ( tolORMin > 0 ) 725 { 744 { 726 // tolORMin2 = tolORMin*tolORMin ; 745 // tolORMin2 = tolORMin*tolORMin ; 727 tolIRMin2 = tolIRMin*tolIRMin ; 746 tolIRMin2 = tolIRMin*tolIRMin ; 728 } 747 } 729 else 748 else 730 { 749 { 731 // tolORMin2 = 0.0 ; 750 // tolORMin2 = 0.0 ; 732 tolIRMin2 = 0.0 ; 751 tolIRMin2 = 0.0 ; 733 } 752 } 734 if ( tolIRMax > 0 ) { tolIRMax2 = tolIR 753 if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; } 735 else { tolIRMax2 = 0.0; 754 else { tolIRMax2 = 0.0; } 736 755 737 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= t 756 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) ) 738 { 757 { 739 if ( !fPhiFullCone && (rhoi2 != 0.0) ) << 758 if ( !fPhiFullCone && rhoi2 ) 740 { 759 { 741 // Psi = angle made with central (av 760 // Psi = angle made with central (average) phi of shape 742 761 743 cosPsi = (xi*cosCPhi + yi*sinCPhi)/s 762 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 744 763 745 if (cosPsi >= cosHDPhiIT) { return 764 if (cosPsi >= cosHDPhiIT) { return sd; } 746 } 765 } 747 else 766 else 748 { 767 { 749 return sd; 768 return sd; 750 } 769 } 751 } 770 } 752 } 771 } 753 else // On/outside extent, and heading aw 772 else // On/outside extent, and heading away -> cannot intersect 754 { 773 { 755 return snxt ; 774 return snxt ; 756 } 775 } 757 } 776 } 758 777 759 // ----> Can not intersect z surfaces 778 // ----> Can not intersect z surfaces 760 779 761 780 762 // Intersection with outer cone (possible retu 781 // Intersection with outer cone (possible return) and 763 // inner cone (must also che 782 // inner cone (must also check phi) 764 // 783 // 765 // Intersection point (xi,yi,zi) on line x=p.x 784 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. 766 // 785 // 767 // Intersects with x^2+y^2=(a*z+b)^2 786 // Intersects with x^2+y^2=(a*z+b)^2 768 // 787 // 769 // where a=tanRMax or tanRMin 788 // where a=tanRMax or tanRMin 770 // b=rMaxAv or rMinAv 789 // b=rMaxAv or rMinAv 771 // 790 // 772 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a 791 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ; 773 // t1 t2 792 // t1 t2 t3 774 // 793 // 775 // \--------u-------/ \-----------v---- 794 // \--------u-------/ \-----------v----------/ \---------w--------/ 776 // 795 // 777 796 778 t1 = 1.0 - v.z()*v.z() ; 797 t1 = 1.0 - v.z()*v.z() ; 779 t2 = p.x()*v.x() + p.y()*v.y() ; 798 t2 = p.x()*v.x() + p.y()*v.y() ; 780 t3 = p.x()*p.x() + p.y()*p.y() ; 799 t3 = p.x()*p.x() + p.y()*p.y() ; 781 rin = tanRMin*p.z() + rMinAv ; 800 rin = tanRMin*p.z() + rMinAv ; 782 rout = tanRMax*p.z() + rMaxAv ; 801 rout = tanRMax*p.z() + rMaxAv ; 783 802 784 // Outer Cone Intersection 803 // Outer Cone Intersection 785 // Must be outside/on outer cone for valid i 804 // Must be outside/on outer cone for valid intersection 786 805 787 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ; 806 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ; 788 nt2 = t2 - tanRMax*v.z()*rout ; 807 nt2 = t2 - tanRMax*v.z()*rout ; 789 nt3 = t3 - rout*rout ; 808 nt3 = t3 - rout*rout ; 790 809 791 if (std::fabs(nt1) > kRadTolerance) // Equa 810 if (std::fabs(nt1) > kRadTolerance) // Equation quadratic => 2 roots 792 { 811 { 793 b = nt2/nt1; 812 b = nt2/nt1; 794 c = nt3/nt1; 813 c = nt3/nt1; 795 d = b*b-c ; 814 d = b*b-c ; 796 if ( (nt3 > rout*rout*kRadTolerance*kRadTo 815 if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax) 797 || (rout < 0) ) 816 || (rout < 0) ) 798 { 817 { 799 // If outside real cone (should be rho-r 818 // If outside real cone (should be rho-rout>kRadTolerance*0.5 800 // NOT rho^2 etc) saves a std::sqrt() at 819 // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy 801 820 802 if (d >= 0) 821 if (d >= 0) 803 { 822 { 804 823 805 if ((rout < 0) && (nt3 <= 0)) 824 if ((rout < 0) && (nt3 <= 0)) 806 { 825 { 807 // Inside `shadow cone' with -ve rad 826 // Inside `shadow cone' with -ve radius 808 // -> 2nd root could be on real cone 827 // -> 2nd root could be on real cone 809 828 810 if (b>0) { sd = c/(-b-std::sqrt(d)); 829 if (b>0) { sd = c/(-b-std::sqrt(d)); } 811 else { sd = -b + std::sqrt(d); 830 else { sd = -b + std::sqrt(d); } 812 } 831 } 813 else 832 else 814 { 833 { 815 if ((b <= 0) && (c >= 0)) // both >= 834 if ((b <= 0) && (c >= 0)) // both >=0, try smaller root 816 { 835 { 817 sd=c/(-b+std::sqrt(d)); 836 sd=c/(-b+std::sqrt(d)); 818 } 837 } 819 else 838 else 820 { 839 { 821 if ( c <= 0 ) // second >=0 840 if ( c <= 0 ) // second >=0 822 { 841 { 823 sd = -b + std::sqrt(d) ; 842 sd = -b + std::sqrt(d) ; 824 if((sd<0.0) && (sd>-halfRadToler << 843 if((sd<0) & (sd>-halfRadTolerance)) sd=0; 825 } 844 } 826 else // both negative, travel awa 845 else // both negative, travel away 827 { 846 { 828 return kInfinity ; 847 return kInfinity ; 829 } 848 } 830 } 849 } 831 } 850 } 832 if ( sd >= 0 ) // If 'forwards'. Chec 851 if ( sd >= 0 ) // If 'forwards'. Check z intersection 833 { 852 { 834 if ( sd>dRmax ) // Avoid rounding er 853 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on 835 { // 64 bits systems. 854 { // 64 bits systems. Split long distances and recompute 836 G4double fTerm = sd-std::fmod(sd,d 855 G4double fTerm = sd-std::fmod(sd,dRmax); 837 sd = fTerm + DistanceToIn(p+fTerm* 856 sd = fTerm + DistanceToIn(p+fTerm*v,v); 838 } 857 } 839 zi = p.z() + sd*v.z() ; 858 zi = p.z() + sd*v.z() ; 840 859 841 if (std::fabs(zi) <= tolODz) 860 if (std::fabs(zi) <= tolODz) 842 { 861 { 843 // Z ok. Check phi intersection if 862 // Z ok. Check phi intersection if reqd 844 863 845 if ( fPhiFullCone ) { return sd; 864 if ( fPhiFullCone ) { return sd; } 846 else 865 else 847 { 866 { 848 xi = p.x() + sd*v.x() ; 867 xi = p.x() + sd*v.x() ; 849 yi = p.y() + sd*v.y() ; 868 yi = p.y() + sd*v.y() ; 850 ri = rMaxAv + zi*tanRMax ; 869 ri = rMaxAv + zi*tanRMax ; 851 cosPsi = (xi*cosCPhi + yi*sinCPh 870 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 852 871 853 if ( cosPsi >= cosHDPhiIT ) { r 872 if ( cosPsi >= cosHDPhiIT ) { return sd; } 854 } 873 } 855 } 874 } 856 } // end if (sd>0) 875 } // end if (sd>0) 857 } 876 } 858 } 877 } 859 else 878 else 860 { 879 { 861 // Inside outer cone 880 // Inside outer cone 862 // check not inside, and heading through 881 // check not inside, and heading through G4Cons (-> 0 to in) 863 882 864 if ( ( t3 > (rin + halfRadTolerance*sec 883 if ( ( t3 > (rin + halfRadTolerance*secRMin)* 865 (rin + halfRadTolerance*sec 884 (rin + halfRadTolerance*secRMin) ) 866 && (nt2 < 0) && (d >= 0) && (std::fabs 885 && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) ) 867 { 886 { 868 // Inside cones, delta r -ve, inside z 887 // Inside cones, delta r -ve, inside z extent 869 // Point is on the Surface => check Di 888 // Point is on the Surface => check Direction using Normal.dot(v) 870 889 871 xi = p.x() ; 890 xi = p.x() ; 872 yi = p.y() ; 891 yi = p.y() ; 873 risec = std::sqrt(xi*xi + yi*yi)*secR 892 risec = std::sqrt(xi*xi + yi*yi)*secRMax ; 874 Normal = G4ThreeVector(xi/risec,yi/ris 893 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ; 875 if ( !fPhiFullCone ) 894 if ( !fPhiFullCone ) 876 { 895 { 877 cosPsi = (p.x()*cosCPhi + p.y()*sinC 896 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ; 878 if ( cosPsi >= cosHDPhiIT ) 897 if ( cosPsi >= cosHDPhiIT ) 879 { 898 { 880 if ( Normal.dot(v) <= 0 ) { retur 899 if ( Normal.dot(v) <= 0 ) { return 0.0; } 881 } 900 } 882 } 901 } 883 else 902 else 884 { 903 { 885 if ( Normal.dot(v) <= 0 ) { return 904 if ( Normal.dot(v) <= 0 ) { return 0.0; } 886 } 905 } 887 } 906 } 888 } 907 } 889 } 908 } 890 else // Single root case 909 else // Single root case 891 { 910 { 892 if ( std::fabs(nt2) > kRadTolerance ) 911 if ( std::fabs(nt2) > kRadTolerance ) 893 { 912 { 894 sd = -0.5*nt3/nt2 ; 913 sd = -0.5*nt3/nt2 ; 895 914 896 if ( sd < 0 ) { return kInfinity; } / 915 if ( sd < 0 ) { return kInfinity; } // travel away 897 else // sd >= 0, If 'forwards'. Check 916 else // sd >= 0, If 'forwards'. Check z intersection 898 { 917 { 899 zi = p.z() + sd*v.z() ; 918 zi = p.z() + sd*v.z() ; 900 919 901 if ((std::fabs(zi) <= tolODz) && (nt2 920 if ((std::fabs(zi) <= tolODz) && (nt2 < 0)) 902 { 921 { 903 // Z ok. Check phi intersection if r 922 // Z ok. Check phi intersection if reqd 904 923 905 if ( fPhiFullCone ) { return sd; } 924 if ( fPhiFullCone ) { return sd; } 906 else 925 else 907 { 926 { 908 xi = p.x() + sd*v.x() ; 927 xi = p.x() + sd*v.x() ; 909 yi = p.y() + sd*v.y() ; 928 yi = p.y() + sd*v.y() ; 910 ri = rMaxAv + zi*tanRMax ; 929 ri = rMaxAv + zi*tanRMax ; 911 cosPsi = (xi*cosCPhi + yi*sinCPhi) 930 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 912 931 913 if (cosPsi >= cosHDPhiIT) { retur 932 if (cosPsi >= cosHDPhiIT) { return sd; } 914 } 933 } 915 } 934 } 916 } 935 } 917 } 936 } 918 else // travel || cone surface from it 937 else // travel || cone surface from its origin 919 { 938 { 920 sd = kInfinity ; 939 sd = kInfinity ; 921 } 940 } 922 } 941 } 923 942 924 // Inner Cone Intersection 943 // Inner Cone Intersection 925 // o Space is divided into 3 areas: 944 // o Space is divided into 3 areas: 926 // 1) Radius greater than real inner cone 945 // 1) Radius greater than real inner cone & imaginary cone & outside 927 // tolerance 946 // tolerance 928 // 2) Radius less than inner or imaginary 947 // 2) Radius less than inner or imaginary cone & outside tolarance 929 // 3) Within tolerance of real or imaginar 948 // 3) Within tolerance of real or imaginary cones 930 // - Extra checks needed for 3's inters 949 // - Extra checks needed for 3's intersections 931 // => lots of duplicated code 950 // => lots of duplicated code 932 951 933 if (rMinAv != 0.0) << 952 if (rMinAv) 934 { 953 { 935 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) 954 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ; 936 nt2 = t2 - tanRMin*v.z()*rin ; 955 nt2 = t2 - tanRMin*v.z()*rin ; 937 nt3 = t3 - rin*rin ; 956 nt3 = t3 - rin*rin ; 938 957 939 if ( nt1 != 0.0 ) << 958 if ( nt1 ) 940 { 959 { 941 if ( nt3 > rin*kRadTolerance*secRMin ) 960 if ( nt3 > rin*kRadTolerance*secRMin ) 942 { 961 { 943 // At radius greater than real & imagi 962 // At radius greater than real & imaginary cones 944 // -> 2nd root, with zi check 963 // -> 2nd root, with zi check 945 964 946 b = nt2/nt1 ; 965 b = nt2/nt1 ; 947 c = nt3/nt1 ; 966 c = nt3/nt1 ; 948 d = b*b-c ; 967 d = b*b-c ; 949 if (d >= 0) // > 0 968 if (d >= 0) // > 0 950 { 969 { 951 if(b>0){sd = c/( -b-std::sqrt(d));} 970 if(b>0){sd = c/( -b-std::sqrt(d));} 952 else {sd = -b + std::sqrt(d) ;} 971 else {sd = -b + std::sqrt(d) ;} 953 972 954 if ( sd >= 0 ) // > 0 973 if ( sd >= 0 ) // > 0 955 { 974 { 956 if ( sd>dRmax ) // Avoid rounding 975 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on 957 { // 64 bits systems 976 { // 64 bits systems. Split long distance and recompute 958 G4double fTerm = sd-std::fmod(sd 977 G4double fTerm = sd-std::fmod(sd,dRmax); 959 sd = fTerm + DistanceToIn(p+fTer 978 sd = fTerm + DistanceToIn(p+fTerm*v,v); 960 } 979 } 961 zi = p.z() + sd*v.z() ; 980 zi = p.z() + sd*v.z() ; 962 981 963 if ( std::fabs(zi) <= tolODz ) 982 if ( std::fabs(zi) <= tolODz ) 964 { 983 { 965 if ( !fPhiFullCone ) 984 if ( !fPhiFullCone ) 966 { 985 { 967 xi = p.x() + sd*v.x() ; 986 xi = p.x() + sd*v.x() ; 968 yi = p.y() + sd*v.y() ; 987 yi = p.y() + sd*v.y() ; 969 ri = rMinAv + zi*tanRMin ; 988 ri = rMinAv + zi*tanRMin ; 970 cosPsi = (xi*cosCPhi + yi*sinC 989 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 971 990 972 if (cosPsi >= cosHDPhiIT) 991 if (cosPsi >= cosHDPhiIT) 973 { 992 { 974 if ( sd > halfRadTolerance ) 993 if ( sd > halfRadTolerance ) { snxt=sd; } 975 else 994 else 976 { 995 { 977 // Calculate a normal vect 996 // Calculate a normal vector in order to check Direction 978 997 979 risec = std::sqrt(xi*xi + 998 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 980 Normal = G4ThreeVector(-xi 999 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin); 981 if ( Normal.dot(v) <= 0 ) 1000 if ( Normal.dot(v) <= 0 ) { snxt = sd; } 982 } 1001 } 983 } 1002 } 984 } 1003 } 985 else 1004 else 986 { 1005 { 987 if ( sd > halfRadTolerance ) 1006 if ( sd > halfRadTolerance ) { return sd; } 988 else 1007 else 989 { 1008 { 990 // Calculate a normal vector 1009 // Calculate a normal vector in order to check Direction 991 1010 992 xi = p.x() + sd*v.x() ; 1011 xi = p.x() + sd*v.x() ; 993 yi = p.y() + sd*v.y() ; 1012 yi = p.y() + sd*v.y() ; 994 risec = std::sqrt(xi*xi + y 1013 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 995 Normal = G4ThreeVector(-xi/r 1014 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ; 996 if ( Normal.dot(v) <= 0 ) { 1015 if ( Normal.dot(v) <= 0 ) { return sd; } 997 } 1016 } 998 } 1017 } 999 } 1018 } 1000 } 1019 } 1001 } 1020 } 1002 } 1021 } 1003 else if ( nt3 < -rin*kRadTolerance*sec 1022 else if ( nt3 < -rin*kRadTolerance*secRMin ) 1004 { 1023 { 1005 // Within radius of inner cone (real 1024 // Within radius of inner cone (real or imaginary) 1006 // -> Try 2nd root, with checking int 1025 // -> Try 2nd root, with checking intersection is with real cone 1007 // -> If check fails, try 1st root, a 1026 // -> If check fails, try 1st root, also checking intersection is 1008 // on real cone 1027 // on real cone 1009 1028 1010 b = nt2/nt1 ; 1029 b = nt2/nt1 ; 1011 c = nt3/nt1 ; 1030 c = nt3/nt1 ; 1012 d = b*b - c ; 1031 d = b*b - c ; 1013 1032 1014 if ( d >= 0 ) // > 0 1033 if ( d >= 0 ) // > 0 1015 { 1034 { 1016 if (b>0) { sd = c/(-b-std::sqrt(d)) 1035 if (b>0) { sd = c/(-b-std::sqrt(d)); } 1017 else { sd = -b + std::sqrt(d); 1036 else { sd = -b + std::sqrt(d); } 1018 zi = p.z() + sd*v.z() ; 1037 zi = p.z() + sd*v.z() ; 1019 ri = rMinAv + zi*tanRMin ; 1038 ri = rMinAv + zi*tanRMin ; 1020 1039 1021 if ( ri > 0 ) 1040 if ( ri > 0 ) 1022 { 1041 { 1023 if ( (sd >= 0) && (std::fabs(zi) 1042 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd > 0 1024 { 1043 { 1025 if ( sd>dRmax ) // Avoid roundi 1044 if ( sd>dRmax ) // Avoid rounding errors due to precision issues 1026 { // seen on 64 b 1045 { // seen on 64 bits systems. Split and recompute 1027 G4double fTerm = sd-std::fmod 1046 G4double fTerm = sd-std::fmod(sd,dRmax); 1028 sd = fTerm + DistanceToIn(p+f 1047 sd = fTerm + DistanceToIn(p+fTerm*v,v); 1029 } 1048 } 1030 if ( !fPhiFullCone ) 1049 if ( !fPhiFullCone ) 1031 { 1050 { 1032 xi = p.x() + sd*v.x() ; 1051 xi = p.x() + sd*v.x() ; 1033 yi = p.y() + sd*v.y() ; 1052 yi = p.y() + sd*v.y() ; 1034 cosPsi = (xi*cosCPhi + yi*sin 1053 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 1035 1054 1036 if (cosPsi >= cosHDPhiOT) 1055 if (cosPsi >= cosHDPhiOT) 1037 { 1056 { 1038 if ( sd > halfRadTolerance 1057 if ( sd > halfRadTolerance ) { snxt=sd; } 1039 else 1058 else 1040 { 1059 { 1041 // Calculate a normal vec 1060 // Calculate a normal vector in order to check Direction 1042 1061 1043 risec = std::sqrt(xi*xi 1062 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 1044 Normal = G4ThreeVector(-x 1063 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin); 1045 if ( Normal.dot(v) <= 0 ) 1064 if ( Normal.dot(v) <= 0 ) { snxt = sd; } 1046 } 1065 } 1047 } 1066 } 1048 } 1067 } 1049 else 1068 else 1050 { 1069 { 1051 if( sd > halfRadTolerance ) 1070 if( sd > halfRadTolerance ) { return sd; } 1052 else 1071 else 1053 { 1072 { 1054 // Calculate a normal vecto 1073 // Calculate a normal vector in order to check Direction 1055 1074 1056 xi = p.x() + sd*v.x() ; 1075 xi = p.x() + sd*v.x() ; 1057 yi = p.y() + sd*v.y() ; 1076 yi = p.y() + sd*v.y() ; 1058 risec = std::sqrt(xi*xi + 1077 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 1059 Normal = G4ThreeVector(-xi/ 1078 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ; 1060 if ( Normal.dot(v) <= 0 ) 1079 if ( Normal.dot(v) <= 0 ) { return sd; } 1061 } 1080 } 1062 } 1081 } 1063 } 1082 } 1064 } 1083 } 1065 else 1084 else 1066 { 1085 { 1067 if (b>0) { sd = -b - std::sqrt(d) 1086 if (b>0) { sd = -b - std::sqrt(d); } 1068 else { sd = c/(-b+std::sqrt(d 1087 else { sd = c/(-b+std::sqrt(d)); } 1069 zi = p.z() + sd*v.z() ; 1088 zi = p.z() + sd*v.z() ; 1070 ri = rMinAv + zi*tanRMin ; 1089 ri = rMinAv + zi*tanRMin ; 1071 1090 1072 if ( (sd >= 0) && (ri > 0) && (st 1091 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0 1073 { 1092 { 1074 if ( sd>dRmax ) // Avoid roundi 1093 if ( sd>dRmax ) // Avoid rounding errors due to precision issues 1075 { // seen on 64 b 1094 { // seen on 64 bits systems. Split and recompute 1076 G4double fTerm = sd-std::fmod 1095 G4double fTerm = sd-std::fmod(sd,dRmax); 1077 sd = fTerm + DistanceToIn(p+f 1096 sd = fTerm + DistanceToIn(p+fTerm*v,v); 1078 } 1097 } 1079 if ( !fPhiFullCone ) 1098 if ( !fPhiFullCone ) 1080 { 1099 { 1081 xi = p.x() + sd*v.x() ; 1100 xi = p.x() + sd*v.x() ; 1082 yi = p.y() + sd*v.y() ; 1101 yi = p.y() + sd*v.y() ; 1083 cosPsi = (xi*cosCPhi + yi*sin 1102 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 1084 1103 1085 if (cosPsi >= cosHDPhiIT) 1104 if (cosPsi >= cosHDPhiIT) 1086 { 1105 { 1087 if ( sd > halfRadTolerance 1106 if ( sd > halfRadTolerance ) { snxt=sd; } 1088 else 1107 else 1089 { 1108 { 1090 // Calculate a normal vec 1109 // Calculate a normal vector in order to check Direction 1091 1110 1092 risec = std::sqrt(xi*xi 1111 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 1093 Normal = G4ThreeVector(-x 1112 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin); 1094 if ( Normal.dot(v) <= 0 ) 1113 if ( Normal.dot(v) <= 0 ) { snxt = sd; } 1095 } 1114 } 1096 } 1115 } 1097 } 1116 } 1098 else 1117 else 1099 { 1118 { 1100 if ( sd > halfRadTolerance ) 1119 if ( sd > halfRadTolerance ) { return sd; } 1101 else 1120 else 1102 { 1121 { 1103 // Calculate a normal vecto 1122 // Calculate a normal vector in order to check Direction 1104 1123 1105 xi = p.x() + sd*v.x() ; 1124 xi = p.x() + sd*v.x() ; 1106 yi = p.y() + sd*v.y() ; 1125 yi = p.y() + sd*v.y() ; 1107 risec = std::sqrt(xi*xi + 1126 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 1108 Normal = G4ThreeVector(-xi/ 1127 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ; 1109 if ( Normal.dot(v) <= 0 ) 1128 if ( Normal.dot(v) <= 0 ) { return sd; } 1110 } 1129 } 1111 } 1130 } 1112 } 1131 } 1113 } 1132 } 1114 } 1133 } 1115 } 1134 } 1116 else 1135 else 1117 { 1136 { 1118 // Within kRadTol*0.5 of inner cone ( 1137 // Within kRadTol*0.5 of inner cone (real OR imaginary) 1119 // ----> Check not travelling through 1138 // ----> Check not travelling through (=>0 to in) 1120 // ----> if not: 1139 // ----> if not: 1121 // -2nd root with validity check 1140 // -2nd root with validity check 1122 1141 1123 if ( std::fabs(p.z()) <= tolODz ) 1142 if ( std::fabs(p.z()) <= tolODz ) 1124 { 1143 { 1125 if ( nt2 > 0 ) 1144 if ( nt2 > 0 ) 1126 { 1145 { 1127 // Inside inner real cone, headin 1146 // Inside inner real cone, heading outwards, inside z range 1128 1147 1129 if ( !fPhiFullCone ) 1148 if ( !fPhiFullCone ) 1130 { 1149 { 1131 cosPsi = (p.x()*cosCPhi + p.y() 1150 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ; 1132 1151 1133 if (cosPsi >= cosHDPhiIT) { re 1152 if (cosPsi >= cosHDPhiIT) { return 0.0; } 1134 } 1153 } 1135 else { return 0.0; } 1154 else { return 0.0; } 1136 } 1155 } 1137 else 1156 else 1138 { 1157 { 1139 // Within z extent, but not trave 1158 // Within z extent, but not travelling through 1140 // -> 2nd root or kInfinity if 1s 1159 // -> 2nd root or kInfinity if 1st root on imaginary cone 1141 1160 1142 b = nt2/nt1 ; 1161 b = nt2/nt1 ; 1143 c = nt3/nt1 ; 1162 c = nt3/nt1 ; 1144 d = b*b - c ; 1163 d = b*b - c ; 1145 1164 1146 if ( d >= 0 ) // > 0 1165 if ( d >= 0 ) // > 0 1147 { 1166 { 1148 if (b>0) { sd = -b - std::sqrt( 1167 if (b>0) { sd = -b - std::sqrt(d); } 1149 else { sd = c/(-b+std::sqrt 1168 else { sd = c/(-b+std::sqrt(d)); } 1150 zi = p.z() + sd*v.z() ; 1169 zi = p.z() + sd*v.z() ; 1151 ri = rMinAv + zi*tanRMin ; 1170 ri = rMinAv + zi*tanRMin ; 1152 1171 1153 if ( ri > 0 ) // 2nd root 1172 if ( ri > 0 ) // 2nd root 1154 { 1173 { 1155 if (b>0) { sd = c/(-b-std::sq 1174 if (b>0) { sd = c/(-b-std::sqrt(d)); } 1156 else { sd = -b + std::sqr 1175 else { sd = -b + std::sqrt(d); } 1157 1176 1158 zi = p.z() + sd*v.z() ; 1177 zi = p.z() + sd*v.z() ; 1159 1178 1160 if ( (sd >= 0) && (std::fabs( 1179 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0 1161 { 1180 { 1162 if ( sd>dRmax ) // Avoid ro 1181 if ( sd>dRmax ) // Avoid rounding errors due to precision issue 1163 { // seen on 1182 { // seen on 64 bits systems. Split and recompute 1164 G4double fTerm = sd-std:: 1183 G4double fTerm = sd-std::fmod(sd,dRmax); 1165 sd = fTerm + DistanceToIn 1184 sd = fTerm + DistanceToIn(p+fTerm*v,v); 1166 } 1185 } 1167 if ( !fPhiFullCone ) 1186 if ( !fPhiFullCone ) 1168 { 1187 { 1169 xi = p.x() + sd*v.x() 1188 xi = p.x() + sd*v.x() ; 1170 yi = p.y() + sd*v.y() 1189 yi = p.y() + sd*v.y() ; 1171 ri = rMinAv + zi*tanR 1190 ri = rMinAv + zi*tanRMin ; 1172 cosPsi = (xi*cosCPhi + yi 1191 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 1173 1192 1174 if ( cosPsi >= cosHDPhiIT 1193 if ( cosPsi >= cosHDPhiIT ) { snxt = sd; } 1175 } 1194 } 1176 else { return sd; } 1195 else { return sd; } 1177 } 1196 } 1178 } 1197 } 1179 else { return kInfinity; } 1198 else { return kInfinity; } 1180 } 1199 } 1181 } 1200 } 1182 } 1201 } 1183 else // 2nd root 1202 else // 2nd root 1184 { 1203 { 1185 b = nt2/nt1 ; 1204 b = nt2/nt1 ; 1186 c = nt3/nt1 ; 1205 c = nt3/nt1 ; 1187 d = b*b - c ; 1206 d = b*b - c ; 1188 1207 1189 if ( d > 0 ) 1208 if ( d > 0 ) 1190 { 1209 { 1191 if (b>0) { sd = c/(-b-std::sqrt(d 1210 if (b>0) { sd = c/(-b-std::sqrt(d)); } 1192 else { sd = -b + std::sqrt(d) 1211 else { sd = -b + std::sqrt(d) ; } 1193 zi = p.z() + sd*v.z() ; 1212 zi = p.z() + sd*v.z() ; 1194 1213 1195 if ( (sd >= 0) && (std::fabs(zi) 1214 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0 1196 { 1215 { 1197 if ( sd>dRmax ) // Avoid roundi 1216 if ( sd>dRmax ) // Avoid rounding errors due to precision issues 1198 { // seen on 64 b 1217 { // seen on 64 bits systems. Split and recompute 1199 G4double fTerm = sd-std::fmod 1218 G4double fTerm = sd-std::fmod(sd,dRmax); 1200 sd = fTerm + DistanceToIn(p+f 1219 sd = fTerm + DistanceToIn(p+fTerm*v,v); 1201 } 1220 } 1202 if ( !fPhiFullCone ) 1221 if ( !fPhiFullCone ) 1203 { 1222 { 1204 xi = p.x() + sd*v.x(); 1223 xi = p.x() + sd*v.x(); 1205 yi = p.y() + sd*v.y(); 1224 yi = p.y() + sd*v.y(); 1206 ri = rMinAv + zi*tanRMin 1225 ri = rMinAv + zi*tanRMin ; 1207 cosPsi = (xi*cosCPhi + yi*sin 1226 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri; 1208 1227 1209 if (cosPsi >= cosHDPhiIT) { 1228 if (cosPsi >= cosHDPhiIT) { snxt = sd; } 1210 } 1229 } 1211 else { return sd; } 1230 else { return sd; } 1212 } 1231 } 1213 } 1232 } 1214 } 1233 } 1215 } 1234 } 1216 } 1235 } 1217 } 1236 } 1218 1237 1219 // Phi segment intersection 1238 // Phi segment intersection 1220 // 1239 // 1221 // o Tolerant of points inside phi planes b 1240 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5 1222 // 1241 // 1223 // o NOTE: Large duplication of code betwee 1242 // o NOTE: Large duplication of code between sphi & ephi checks 1224 // -> only diffs: sphi -> ephi, Com 1243 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane 1225 // intersection check <=0 -> >=0 1244 // intersection check <=0 -> >=0 1226 // -> Should use some form of loop 1245 // -> Should use some form of loop Construct 1227 1246 1228 if ( !fPhiFullCone ) 1247 if ( !fPhiFullCone ) 1229 { 1248 { 1230 // First phi surface (starting phi) 1249 // First phi surface (starting phi) 1231 1250 1232 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 1251 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 1233 1252 1234 if ( Comp < 0 ) // Component in outwar 1253 if ( Comp < 0 ) // Component in outwards normal dirn 1235 { 1254 { 1236 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) 1255 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ; 1237 1256 1238 if (Dist < halfCarTolerance) 1257 if (Dist < halfCarTolerance) 1239 { 1258 { 1240 sd = Dist/Comp ; 1259 sd = Dist/Comp ; 1241 1260 1242 if ( sd < snxt ) 1261 if ( sd < snxt ) 1243 { 1262 { 1244 if ( sd < 0 ) { sd = 0.0; } 1263 if ( sd < 0 ) { sd = 0.0; } 1245 1264 1246 zi = p.z() + sd*v.z() ; 1265 zi = p.z() + sd*v.z() ; 1247 1266 1248 if ( std::fabs(zi) <= tolODz ) 1267 if ( std::fabs(zi) <= tolODz ) 1249 { 1268 { 1250 xi = p.x() + sd*v.x() ; 1269 xi = p.x() + sd*v.x() ; 1251 yi = p.y() + sd*v.y() ; 1270 yi = p.y() + sd*v.y() ; 1252 rhoi2 = xi*xi + yi*yi ; 1271 rhoi2 = xi*xi + yi*yi ; 1253 tolORMin2 = (rMinOAv + zi*tanRMin 1272 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ; 1254 tolORMax2 = (rMaxOAv + zi*tanRMax 1273 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ; 1255 1274 1256 if ( (rhoi2 >= tolORMin2) && (rho 1275 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) ) 1257 { 1276 { 1258 // z and r intersections good - 1277 // z and r intersections good - check intersecting with 1259 // correct half-plane 1278 // correct half-plane 1260 1279 1261 if ((yi*cosCPhi - xi*sinCPhi) < 1280 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; } 1262 } 1281 } 1263 } 1282 } 1264 } 1283 } 1265 } 1284 } 1266 } 1285 } 1267 1286 1268 // Second phi surface (Ending phi) 1287 // Second phi surface (Ending phi) 1269 1288 1270 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi 1289 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ; 1271 1290 1272 if ( Comp < 0 ) // Component in outward 1291 if ( Comp < 0 ) // Component in outwards normal dirn 1273 { 1292 { 1274 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) 1293 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ; 1275 if (Dist < halfCarTolerance) 1294 if (Dist < halfCarTolerance) 1276 { 1295 { 1277 sd = Dist/Comp ; 1296 sd = Dist/Comp ; 1278 1297 1279 if ( sd < snxt ) 1298 if ( sd < snxt ) 1280 { 1299 { 1281 if ( sd < 0 ) { sd = 0.0; } 1300 if ( sd < 0 ) { sd = 0.0; } 1282 1301 1283 zi = p.z() + sd*v.z() ; 1302 zi = p.z() + sd*v.z() ; 1284 1303 1285 if (std::fabs(zi) <= tolODz) 1304 if (std::fabs(zi) <= tolODz) 1286 { 1305 { 1287 xi = p.x() + sd*v.x() ; 1306 xi = p.x() + sd*v.x() ; 1288 yi = p.y() + sd*v.y() ; 1307 yi = p.y() + sd*v.y() ; 1289 rhoi2 = xi*xi + yi*yi ; 1308 rhoi2 = xi*xi + yi*yi ; 1290 tolORMin2 = (rMinOAv + zi*tanRMin 1309 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ; 1291 tolORMax2 = (rMaxOAv + zi*tanRMax 1310 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ; 1292 1311 1293 if ( (rhoi2 >= tolORMin2) && (rho 1312 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) ) 1294 { 1313 { 1295 // z and r intersections good - 1314 // z and r intersections good - check intersecting with 1296 // correct half-plane 1315 // correct half-plane 1297 1316 1298 if ( (yi*cosCPhi - xi*sinCPhi) 1317 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; } 1299 } 1318 } 1300 } 1319 } 1301 } 1320 } 1302 } 1321 } 1303 } 1322 } 1304 } 1323 } 1305 if (snxt < halfCarTolerance) { snxt = 0.; 1324 if (snxt < halfCarTolerance) { snxt = 0.; } 1306 1325 1307 return snxt ; 1326 return snxt ; 1308 } 1327 } 1309 1328 1310 ///////////////////////////////////////////// 1329 ////////////////////////////////////////////////////////////////////////////// 1311 // 1330 // 1312 // Calculate distance (<= actual) to closest 1331 // Calculate distance (<= actual) to closest surface of shape from outside 1313 // - Calculate distance to z, radial planes 1332 // - Calculate distance to z, radial planes 1314 // - Only to phi planes if outside phi extent 1333 // - Only to phi planes if outside phi extent 1315 // - Return 0 if point inside 1334 // - Return 0 if point inside 1316 1335 1317 G4double G4Cons::DistanceToIn(const G4ThreeVe 1336 G4double G4Cons::DistanceToIn(const G4ThreeVector& p) const 1318 { 1337 { 1319 G4double safe=0.0, rho, safeR1, safeR2, saf 1338 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ; 1320 G4double tanRMin, secRMin, pRMin ; 1339 G4double tanRMin, secRMin, pRMin ; 1321 G4double tanRMax, secRMax, pRMax ; 1340 G4double tanRMax, secRMax, pRMax ; 1322 1341 1323 rho = std::sqrt(p.x()*p.x() + p.y()*p.y() 1342 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 1324 safeZ = std::fabs(p.z()) - fDz ; 1343 safeZ = std::fabs(p.z()) - fDz ; 1325 1344 1326 if ( (fRmin1 != 0.0) || (fRmin2 != 0.0) ) << 1345 if ( fRmin1 || fRmin2 ) 1327 { 1346 { 1328 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 1347 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 1329 secRMin = std::sqrt(1.0 + tanRMin*tanRMin 1348 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ; 1330 pRMin = tanRMin*p.z() + (fRmin1 + fRmin 1349 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ; 1331 safeR1 = (pRMin - rho)/secRMin ; 1350 safeR1 = (pRMin - rho)/secRMin ; 1332 1351 1333 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 1352 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 1334 secRMax = std::sqrt(1.0 + tanRMax*tanRMax 1353 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 1335 pRMax = tanRMax*p.z() + (fRmax1 + fRmax 1354 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ; 1336 safeR2 = (rho - pRMax)/secRMax ; 1355 safeR2 = (rho - pRMax)/secRMax ; 1337 1356 1338 if ( safeR1 > safeR2) { safe = safeR1; } 1357 if ( safeR1 > safeR2) { safe = safeR1; } 1339 else { safe = safeR2; } 1358 else { safe = safeR2; } 1340 } 1359 } 1341 else 1360 else 1342 { 1361 { 1343 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 1362 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 1344 secRMax = std::sqrt(1.0 + tanRMax*tanRMax 1363 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 1345 pRMax = tanRMax*p.z() + (fRmax1 + fRmax 1364 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ; 1346 safe = (rho - pRMax)/secRMax ; 1365 safe = (rho - pRMax)/secRMax ; 1347 } 1366 } 1348 if ( safeZ > safe ) { safe = safeZ; } 1367 if ( safeZ > safe ) { safe = safeZ; } 1349 1368 1350 if ( !fPhiFullCone && (rho != 0.0) ) << 1369 if ( !fPhiFullCone && rho ) 1351 { 1370 { 1352 // Psi=angle from central phi to point 1371 // Psi=angle from central phi to point 1353 1372 1354 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/ 1373 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ; 1355 1374 1356 if ( cosPsi < cosHDPhi ) // Point lies ou 1375 if ( cosPsi < cosHDPhi ) // Point lies outside phi range 1357 { 1376 { 1358 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 1377 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 ) 1359 { 1378 { 1360 safePhi = std::fabs(p.x()*sinSPhi-p.y 1379 safePhi = std::fabs(p.x()*sinSPhi-p.y()*cosSPhi); 1361 } 1380 } 1362 else 1381 else 1363 { 1382 { 1364 safePhi = std::fabs(p.x()*sinEPhi-p.y 1383 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi); 1365 } 1384 } 1366 if ( safePhi > safe ) { safe = safePhi 1385 if ( safePhi > safe ) { safe = safePhi; } 1367 } 1386 } 1368 } 1387 } 1369 if ( safe < 0.0 ) { safe = 0.0; } 1388 if ( safe < 0.0 ) { safe = 0.0; } 1370 1389 1371 return safe ; 1390 return safe ; 1372 } 1391 } 1373 1392 1374 ///////////////////////////////////////////// 1393 /////////////////////////////////////////////////////////////// 1375 // 1394 // 1376 // Calculate distance to surface of shape fro 1395 // Calculate distance to surface of shape from 'inside', allowing for tolerance 1377 // - Only Calc rmax intersection if no valid 1396 // - Only Calc rmax intersection if no valid rmin intersection 1378 1397 1379 G4double G4Cons::DistanceToOut( const G4Three 1398 G4double G4Cons::DistanceToOut( const G4ThreeVector& p, 1380 const G4Three 1399 const G4ThreeVector& v, 1381 const G4bool 1400 const G4bool calcNorm, 1382 G4bool* 1401 G4bool* validNorm, 1383 G4Three 1402 G4ThreeVector* n) const 1384 { 1403 { 1385 ESide side = kNull, sider = kNull, sidephi 1404 ESide side = kNull, sider = kNull, sidephi = kNull; 1386 1405 1387 G4double snxt,srd,sphi,pdist ; 1406 G4double snxt,srd,sphi,pdist ; 1388 1407 1389 G4double tanRMax, secRMax, rMaxAv ; // Dat 1408 G4double tanRMax, secRMax, rMaxAv ; // Data for outer cone 1390 G4double tanRMin, secRMin, rMinAv ; // Dat 1409 G4double tanRMin, secRMin, rMinAv ; // Data for inner cone 1391 1410 1392 G4double t1, t2, t3, rout, rin, nt1, nt2, n 1411 G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ; 1393 G4double b, c, d, sr2, sr3 ; 1412 G4double b, c, d, sr2, sr3 ; 1394 1413 1395 // Vars for intersection within tolerance 1414 // Vars for intersection within tolerance 1396 1415 1397 ESide sidetol = kNull ; 1416 ESide sidetol = kNull ; 1398 G4double slentol = kInfinity ; 1417 G4double slentol = kInfinity ; 1399 1418 1400 // Vars for phi intersection: 1419 // Vars for phi intersection: 1401 1420 1402 G4double pDistS, compS, pDistE, compE, sphi 1421 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ; 1403 G4double zi, ri, deltaRoi2 ; 1422 G4double zi, ri, deltaRoi2 ; 1404 1423 1405 // Z plane intersection 1424 // Z plane intersection 1406 1425 1407 if ( v.z() > 0.0 ) 1426 if ( v.z() > 0.0 ) 1408 { 1427 { 1409 pdist = fDz - p.z() ; 1428 pdist = fDz - p.z() ; 1410 1429 1411 if (pdist > halfCarTolerance) 1430 if (pdist > halfCarTolerance) 1412 { 1431 { 1413 snxt = pdist/v.z() ; 1432 snxt = pdist/v.z() ; 1414 side = kPZ ; 1433 side = kPZ ; 1415 } 1434 } 1416 else 1435 else 1417 { 1436 { 1418 if (calcNorm) 1437 if (calcNorm) 1419 { 1438 { 1420 *n = G4ThreeVector(0,0,1) ; 1439 *n = G4ThreeVector(0,0,1) ; 1421 *validNorm = true ; 1440 *validNorm = true ; 1422 } 1441 } 1423 return snxt = 0.0; 1442 return snxt = 0.0; 1424 } 1443 } 1425 } 1444 } 1426 else if ( v.z() < 0.0 ) 1445 else if ( v.z() < 0.0 ) 1427 { 1446 { 1428 pdist = fDz + p.z() ; 1447 pdist = fDz + p.z() ; 1429 1448 1430 if ( pdist > halfCarTolerance) 1449 if ( pdist > halfCarTolerance) 1431 { 1450 { 1432 snxt = -pdist/v.z() ; 1451 snxt = -pdist/v.z() ; 1433 side = kMZ ; 1452 side = kMZ ; 1434 } 1453 } 1435 else 1454 else 1436 { 1455 { 1437 if ( calcNorm ) 1456 if ( calcNorm ) 1438 { 1457 { 1439 *n = G4ThreeVector(0,0,-1) ; 1458 *n = G4ThreeVector(0,0,-1) ; 1440 *validNorm = true ; 1459 *validNorm = true ; 1441 } 1460 } 1442 return snxt = 0.0 ; 1461 return snxt = 0.0 ; 1443 } 1462 } 1444 } 1463 } 1445 else // Travel perpendicular to z axis 1464 else // Travel perpendicular to z axis 1446 { 1465 { 1447 snxt = kInfinity ; 1466 snxt = kInfinity ; 1448 side = kNull ; 1467 side = kNull ; 1449 } 1468 } 1450 1469 1451 // Radial Intersections 1470 // Radial Intersections 1452 // 1471 // 1453 // Intersection with outer cone (possible r 1472 // Intersection with outer cone (possible return) and 1454 // inner cone (must also 1473 // inner cone (must also check phi) 1455 // 1474 // 1456 // Intersection point (xi,yi,zi) on line x= 1475 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. 1457 // 1476 // 1458 // Intersects with x^2+y^2=(a*z+b)^2 1477 // Intersects with x^2+y^2=(a*z+b)^2 1459 // 1478 // 1460 // where a=tanRMax or tanRMin 1479 // where a=tanRMax or tanRMin 1461 // b=rMaxAv or rMinAv 1480 // b=rMaxAv or rMinAv 1462 // 1481 // 1463 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*v 1482 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ; 1464 // t1 t2 1483 // t1 t2 t3 1465 // 1484 // 1466 // \--------u-------/ \-----------v- 1485 // \--------u-------/ \-----------v----------/ \---------w--------/ 1467 1486 1468 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 1487 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 1469 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) 1488 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 1470 rMaxAv = (fRmax1 + fRmax2)*0.5 ; 1489 rMaxAv = (fRmax1 + fRmax2)*0.5 ; 1471 1490 1472 1491 1473 t1 = 1.0 - v.z()*v.z() ; // since v 1492 t1 = 1.0 - v.z()*v.z() ; // since v normalised 1474 t2 = p.x()*v.x() + p.y()*v.y() ; 1493 t2 = p.x()*v.x() + p.y()*v.y() ; 1475 t3 = p.x()*p.x() + p.y()*p.y() ; 1494 t3 = p.x()*p.x() + p.y()*p.y() ; 1476 rout = tanRMax*p.z() + rMaxAv ; 1495 rout = tanRMax*p.z() + rMaxAv ; 1477 1496 1478 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) 1497 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ; 1479 nt2 = t2 - tanRMax*v.z()*rout ; 1498 nt2 = t2 - tanRMax*v.z()*rout ; 1480 nt3 = t3 - rout*rout ; 1499 nt3 = t3 - rout*rout ; 1481 1500 1482 if (v.z() > 0.0) 1501 if (v.z() > 0.0) 1483 { 1502 { 1484 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3 1503 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3 1485 - fRmax2*(fRmax2 + kRadTolera 1504 - fRmax2*(fRmax2 + kRadTolerance*secRMax); 1486 } 1505 } 1487 else if (v.z() < 0.0) 1506 else if (v.z() < 0.0) 1488 { 1507 { 1489 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3 1508 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3 1490 - fRmax1*(fRmax1 + kRadTolera 1509 - fRmax1*(fRmax1 + kRadTolerance*secRMax); 1491 } 1510 } 1492 else 1511 else 1493 { 1512 { 1494 deltaRoi2 = 1.0; 1513 deltaRoi2 = 1.0; 1495 } 1514 } 1496 1515 1497 if ( (nt1 != 0.0) && (deltaRoi2 > 0.0) ) << 1516 if ( nt1 && (deltaRoi2 > 0.0) ) 1498 { 1517 { 1499 // Equation quadratic => 2 roots : second 1518 // Equation quadratic => 2 roots : second root must be leaving 1500 1519 1501 b = nt2/nt1 ; 1520 b = nt2/nt1 ; 1502 c = nt3/nt1 ; 1521 c = nt3/nt1 ; 1503 d = b*b - c ; 1522 d = b*b - c ; 1504 1523 1505 if ( d >= 0 ) 1524 if ( d >= 0 ) 1506 { 1525 { 1507 // Check if on outer cone & heading out 1526 // Check if on outer cone & heading outwards 1508 // NOTE: Should use rho-rout>-kRadToler 1527 // NOTE: Should use rho-rout>-kRadTolerance*0.5 1509 1528 1510 if (nt3 > -halfRadTolerance && nt2 >= 0 1529 if (nt3 > -halfRadTolerance && nt2 >= 0 ) 1511 { 1530 { 1512 if (calcNorm) 1531 if (calcNorm) 1513 { 1532 { 1514 risec = std::sqrt(t3)*secRMax 1533 risec = std::sqrt(t3)*secRMax ; 1515 *validNorm = true ; 1534 *validNorm = true ; 1516 *n = G4ThreeVector(p.x()/ri 1535 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax); 1517 } 1536 } 1518 return snxt=0 ; 1537 return snxt=0 ; 1519 } 1538 } 1520 else 1539 else 1521 { 1540 { 1522 sider = kRMax ; 1541 sider = kRMax ; 1523 if (b>0) { srd = -b - std::sqrt(d); 1542 if (b>0) { srd = -b - std::sqrt(d); } 1524 else { srd = c/(-b+std::sqrt(d)) 1543 else { srd = c/(-b+std::sqrt(d)) ; } 1525 1544 1526 zi = p.z() + srd*v.z() ; 1545 zi = p.z() + srd*v.z() ; 1527 ri = tanRMax*zi + rMaxAv ; 1546 ri = tanRMax*zi + rMaxAv ; 1528 1547 1529 if ((ri >= 0) && (-halfRadTolerance < 1548 if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance)) 1530 { 1549 { 1531 // An intersection within the toler 1550 // An intersection within the tolerance 1532 // we will Store it in case it is 1551 // we will Store it in case it is good - 1533 // 1552 // 1534 slentol = srd ; 1553 slentol = srd ; 1535 sidetol = kRMax ; 1554 sidetol = kRMax ; 1536 } 1555 } 1537 if ( (ri < 0) || (srd < halfRadTolera 1556 if ( (ri < 0) || (srd < halfRadTolerance) ) 1538 { 1557 { 1539 // Safety: if both roots -ve ensure 1558 // Safety: if both roots -ve ensure that srd cannot `win' 1540 // distance to out 1559 // distance to out 1541 1560 1542 if (b>0) { sr2 = c/(-b-std::sqrt(d) 1561 if (b>0) { sr2 = c/(-b-std::sqrt(d)); } 1543 else { sr2 = -b + std::sqrt(d); 1562 else { sr2 = -b + std::sqrt(d); } 1544 zi = p.z() + sr2*v.z() ; 1563 zi = p.z() + sr2*v.z() ; 1545 ri = tanRMax*zi + rMaxAv ; 1564 ri = tanRMax*zi + rMaxAv ; 1546 1565 1547 if ((ri >= 0) && (sr2 > halfRadTole 1566 if ((ri >= 0) && (sr2 > halfRadTolerance)) 1548 { 1567 { 1549 srd = sr2; 1568 srd = sr2; 1550 } 1569 } 1551 else 1570 else 1552 { 1571 { 1553 srd = kInfinity ; 1572 srd = kInfinity ; 1554 1573 1555 if( (-halfRadTolerance <= sr2) && 1574 if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) ) 1556 { 1575 { 1557 // An intersection within the t 1576 // An intersection within the tolerance. 1558 // Storing it in case it is goo 1577 // Storing it in case it is good. 1559 1578 1560 slentol = sr2 ; 1579 slentol = sr2 ; 1561 sidetol = kRMax ; 1580 sidetol = kRMax ; 1562 } 1581 } 1563 } 1582 } 1564 } 1583 } 1565 } 1584 } 1566 } 1585 } 1567 else 1586 else 1568 { 1587 { 1569 // No intersection with outer cone & no 1588 // No intersection with outer cone & not parallel 1570 // -> already outside, no intersection 1589 // -> already outside, no intersection 1571 1590 1572 if ( calcNorm ) 1591 if ( calcNorm ) 1573 { 1592 { 1574 risec = std::sqrt(t3)*secRMax; 1593 risec = std::sqrt(t3)*secRMax; 1575 *validNorm = true; 1594 *validNorm = true; 1576 *n = G4ThreeVector(p.x()/rise 1595 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax); 1577 } 1596 } 1578 return snxt = 0.0 ; 1597 return snxt = 0.0 ; 1579 } 1598 } 1580 } 1599 } 1581 else if ( (nt2 != 0.0) && (deltaRoi2 > 0.0) << 1600 else if ( nt2 && (deltaRoi2 > 0.0) ) 1582 { 1601 { 1583 // Linear case (only one intersection) => 1602 // Linear case (only one intersection) => point outside outer cone 1584 1603 1585 if ( calcNorm ) 1604 if ( calcNorm ) 1586 { 1605 { 1587 risec = std::sqrt(t3)*secRMax; 1606 risec = std::sqrt(t3)*secRMax; 1588 *validNorm = true; 1607 *validNorm = true; 1589 *n = G4ThreeVector(p.x()/risec, 1608 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax); 1590 } 1609 } 1591 return snxt = 0.0 ; 1610 return snxt = 0.0 ; 1592 } 1611 } 1593 else 1612 else 1594 { 1613 { 1595 // No intersection -> parallel to outer c 1614 // No intersection -> parallel to outer cone 1596 // => Z or inner cone intersection 1615 // => Z or inner cone intersection 1597 1616 1598 srd = kInfinity ; 1617 srd = kInfinity ; 1599 } 1618 } 1600 1619 1601 // Check possible intersection within toler 1620 // Check possible intersection within tolerance 1602 1621 1603 if ( slentol <= halfCarTolerance ) 1622 if ( slentol <= halfCarTolerance ) 1604 { 1623 { 1605 // An intersection within the tolerance w 1624 // An intersection within the tolerance was found. 1606 // We must accept it only if the momentum 1625 // We must accept it only if the momentum points outwards. 1607 // 1626 // 1608 // G4ThreeVector ptTol ; // The point of 1627 // G4ThreeVector ptTol ; // The point of the intersection 1609 // ptTol= p + slentol*v ; 1628 // ptTol= p + slentol*v ; 1610 // ri=tanRMax*zi+rMaxAv ; 1629 // ri=tanRMax*zi+rMaxAv ; 1611 // 1630 // 1612 // Calculate a normal vector, as below 1631 // Calculate a normal vector, as below 1613 1632 1614 xi = p.x() + slentol*v.x(); 1633 xi = p.x() + slentol*v.x(); 1615 yi = p.y() + slentol*v.y(); 1634 yi = p.y() + slentol*v.y(); 1616 risec = std::sqrt(xi*xi + yi*yi)*secRMax; 1635 risec = std::sqrt(xi*xi + yi*yi)*secRMax; 1617 G4ThreeVector Normal = G4ThreeVector(xi/r 1636 G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax); 1618 1637 1619 if ( Normal.dot(v) > 0 ) // We will le 1638 if ( Normal.dot(v) > 0 ) // We will leave the Cone immediatelly 1620 { 1639 { 1621 if ( calcNorm ) 1640 if ( calcNorm ) 1622 { 1641 { 1623 *n = Normal.unit() ; 1642 *n = Normal.unit() ; 1624 *validNorm = true ; 1643 *validNorm = true ; 1625 } 1644 } 1626 return snxt = 0.0 ; 1645 return snxt = 0.0 ; 1627 } 1646 } 1628 else // On the surface, but not heading o 1647 else // On the surface, but not heading out so we ignore this intersection 1629 { // 1648 { // (as it is within tolerance). 1630 slentol = kInfinity ; 1649 slentol = kInfinity ; 1631 } 1650 } 1632 } 1651 } 1633 1652 1634 // Inner Cone intersection 1653 // Inner Cone intersection 1635 1654 1636 if ( (fRmin1 != 0.0) || (fRmin2 != 0.0) ) << 1655 if ( fRmin1 || fRmin2 ) 1637 { 1656 { 1638 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 1657 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 1639 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v 1658 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ; 1640 1659 1641 if ( nt1 != 0.0 ) << 1660 if ( nt1 ) 1642 { 1661 { 1643 secRMin = std::sqrt(1.0 + tanRMin*tanRM 1662 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ; 1644 rMinAv = (fRmin1 + fRmin2)*0.5 ; 1663 rMinAv = (fRmin1 + fRmin2)*0.5 ; 1645 rin = tanRMin*p.z() + rMinAv ; 1664 rin = tanRMin*p.z() + rMinAv ; 1646 nt2 = t2 - tanRMin*v.z()*rin ; 1665 nt2 = t2 - tanRMin*v.z()*rin ; 1647 nt3 = t3 - rin*rin ; 1666 nt3 = t3 - rin*rin ; 1648 1667 1649 // Equation quadratic => 2 roots : firs 1668 // Equation quadratic => 2 roots : first root must be leaving 1650 1669 1651 b = nt2/nt1 ; 1670 b = nt2/nt1 ; 1652 c = nt3/nt1 ; 1671 c = nt3/nt1 ; 1653 d = b*b - c ; 1672 d = b*b - c ; 1654 1673 1655 if ( d >= 0.0 ) 1674 if ( d >= 0.0 ) 1656 { 1675 { 1657 // NOTE: should be rho-rin<kRadTolera 1676 // NOTE: should be rho-rin<kRadTolerance*0.5, 1658 // but using squared versions f 1677 // but using squared versions for efficiency 1659 1678 1660 if (nt3 < kRadTolerance*(rin + kRadTo 1679 if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25)) 1661 { 1680 { 1662 if ( nt2 < 0.0 ) 1681 if ( nt2 < 0.0 ) 1663 { 1682 { 1664 if (calcNorm) { *validNorm = fal 1683 if (calcNorm) { *validNorm = false; } 1665 return snxt = 0.0; 1684 return snxt = 0.0; 1666 } 1685 } 1667 } 1686 } 1668 else 1687 else 1669 { 1688 { 1670 if (b>0) { sr2 = -b - std::sqrt(d); 1689 if (b>0) { sr2 = -b - std::sqrt(d); } 1671 else { sr2 = c/(-b+std::sqrt(d) 1690 else { sr2 = c/(-b+std::sqrt(d)); } 1672 zi = p.z() + sr2*v.z() ; 1691 zi = p.z() + sr2*v.z() ; 1673 ri = tanRMin*zi + rMinAv ; 1692 ri = tanRMin*zi + rMinAv ; 1674 1693 1675 if( (ri>=0.0)&&(-halfRadTolerance<= 1694 if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) ) 1676 { 1695 { 1677 // An intersection within the tol 1696 // An intersection within the tolerance 1678 // storing it in case it is good. 1697 // storing it in case it is good. 1679 1698 1680 slentol = sr2 ; 1699 slentol = sr2 ; 1681 sidetol = kRMax ; 1700 sidetol = kRMax ; 1682 } 1701 } 1683 if( (ri<0) || (sr2 < halfRadToleran 1702 if( (ri<0) || (sr2 < halfRadTolerance) ) 1684 { 1703 { 1685 if (b>0) { sr3 = c/(-b-std::sqrt( 1704 if (b>0) { sr3 = c/(-b-std::sqrt(d)); } 1686 else { sr3 = -b + std::sqrt(d 1705 else { sr3 = -b + std::sqrt(d) ; } 1687 1706 1688 // Safety: if both roots -ve ensu 1707 // Safety: if both roots -ve ensure that srd cannot `win' 1689 // distancetoout 1708 // distancetoout 1690 1709 1691 if ( sr3 > halfRadTolerance ) 1710 if ( sr3 > halfRadTolerance ) 1692 { 1711 { 1693 if( sr3 < srd ) 1712 if( sr3 < srd ) 1694 { 1713 { 1695 zi = p.z() + sr3*v.z() ; 1714 zi = p.z() + sr3*v.z() ; 1696 ri = tanRMin*zi + rMinAv ; 1715 ri = tanRMin*zi + rMinAv ; 1697 1716 1698 if ( ri >= 0.0 ) 1717 if ( ri >= 0.0 ) 1699 { 1718 { 1700 srd=sr3 ; 1719 srd=sr3 ; 1701 sider=kRMin ; 1720 sider=kRMin ; 1702 } 1721 } 1703 } 1722 } 1704 } 1723 } 1705 else if ( sr3 > -halfRadTolerance 1724 else if ( sr3 > -halfRadTolerance ) 1706 { 1725 { 1707 // Intersection in tolerance. S 1726 // Intersection in tolerance. Store to check if it's good 1708 1727 1709 slentol = sr3 ; 1728 slentol = sr3 ; 1710 sidetol = kRMin ; 1729 sidetol = kRMin ; 1711 } 1730 } 1712 } 1731 } 1713 else if ( (sr2 < srd) && (sr2 > hal 1732 else if ( (sr2 < srd) && (sr2 > halfCarTolerance) ) 1714 { 1733 { 1715 srd = sr2 ; 1734 srd = sr2 ; 1716 sider = kRMin ; 1735 sider = kRMin ; 1717 } 1736 } 1718 else if (sr2 > -halfCarTolerance) 1737 else if (sr2 > -halfCarTolerance) 1719 { 1738 { 1720 // Intersection in tolerance. Sto 1739 // Intersection in tolerance. Store to check if it's good 1721 1740 1722 slentol = sr2 ; 1741 slentol = sr2 ; 1723 sidetol = kRMin ; 1742 sidetol = kRMin ; 1724 } 1743 } 1725 if( slentol <= halfCarTolerance ) 1744 if( slentol <= halfCarTolerance ) 1726 { 1745 { 1727 // An intersection within the tol 1746 // An intersection within the tolerance was found. 1728 // We must accept it only if the 1747 // We must accept it only if the momentum points outwards. 1729 1748 1730 G4ThreeVector Normal ; 1749 G4ThreeVector Normal ; 1731 1750 1732 // Calculate a normal vector, as 1751 // Calculate a normal vector, as below 1733 1752 1734 xi = p.x() + slentol*v.x() ; 1753 xi = p.x() + slentol*v.x() ; 1735 yi = p.y() + slentol*v.y() ; 1754 yi = p.y() + slentol*v.y() ; 1736 if( sidetol==kRMax ) 1755 if( sidetol==kRMax ) 1737 { 1756 { 1738 risec = std::sqrt(xi*xi + yi*y 1757 risec = std::sqrt(xi*xi + yi*yi)*secRMax ; 1739 Normal = G4ThreeVector(xi/risec 1758 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ; 1740 } 1759 } 1741 else 1760 else 1742 { 1761 { 1743 risec = std::sqrt(xi*xi + yi*y 1762 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 1744 Normal = G4ThreeVector(-xi/rise 1763 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ; 1745 } 1764 } 1746 if( Normal.dot(v) > 0 ) 1765 if( Normal.dot(v) > 0 ) 1747 { 1766 { 1748 // We will leave the cone immed 1767 // We will leave the cone immediately 1749 1768 1750 if( calcNorm ) 1769 if( calcNorm ) 1751 { 1770 { 1752 *n = Normal.unit() ; 1771 *n = Normal.unit() ; 1753 *validNorm = true ; 1772 *validNorm = true ; 1754 } 1773 } 1755 return snxt = 0.0 ; 1774 return snxt = 0.0 ; 1756 } 1775 } 1757 else 1776 else 1758 { 1777 { 1759 // On the surface, but not head 1778 // On the surface, but not heading out so we ignore this 1760 // intersection (as it is withi 1779 // intersection (as it is within tolerance). 1761 1780 1762 slentol = kInfinity ; 1781 slentol = kInfinity ; 1763 } 1782 } 1764 } 1783 } 1765 } 1784 } 1766 } 1785 } 1767 } 1786 } 1768 } 1787 } 1769 1788 1770 // Linear case => point outside inner cone 1789 // Linear case => point outside inner cone ---> outer cone intersect 1771 // 1790 // 1772 // Phi Intersection 1791 // Phi Intersection 1773 1792 1774 if ( !fPhiFullCone ) 1793 if ( !fPhiFullCone ) 1775 { 1794 { 1776 // add angle calculation with correction 1795 // add angle calculation with correction 1777 // of the difference in domain of atan2 a 1796 // of the difference in domain of atan2 and Sphi 1778 1797 1779 vphi = std::atan2(v.y(),v.x()) ; 1798 vphi = std::atan2(v.y(),v.x()) ; 1780 1799 1781 if ( vphi < fSPhi - halfAngTolerance ) 1800 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; } 1782 else if ( vphi > fSPhi + fDPhi + halfAngT 1801 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; } 1783 1802 1784 if ( (p.x() != 0.0) || (p.y() != 0.0) ) << 1803 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later) 1785 { 1804 { 1786 // pDist -ve when inside 1805 // pDist -ve when inside 1787 1806 1788 pDistS = p.x()*sinSPhi - p.y()*cosSPhi 1807 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; 1789 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi 1808 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ; 1790 1809 1791 // Comp -ve when in direction of outwar 1810 // Comp -ve when in direction of outwards normal 1792 1811 1793 compS = -sinSPhi*v.x() + cosSPhi*v.y() 1812 compS = -sinSPhi*v.x() + cosSPhi*v.y() ; 1794 compE = sinEPhi*v.x() - cosEPhi*v.y() ; 1813 compE = sinEPhi*v.x() - cosEPhi*v.y() ; 1795 1814 1796 sidephi = kNull ; 1815 sidephi = kNull ; 1797 1816 1798 if( ( (fDPhi <= pi) && ( (pDistS <= hal 1817 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance) 1799 && (pDistE <= hal 1818 && (pDistE <= halfCarTolerance) ) ) 1800 || ( (fDPhi > pi) && ((pDistS <= h << 1819 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance) 1801 || (pDistE <= << 1820 && (pDistE > halfCarTolerance) ) ) ) 1802 { 1821 { 1803 // Inside both phi *full* planes 1822 // Inside both phi *full* planes 1804 if ( compS < 0 ) 1823 if ( compS < 0 ) 1805 { 1824 { 1806 sphi = pDistS/compS ; 1825 sphi = pDistS/compS ; 1807 if (sphi >= -halfCarTolerance) 1826 if (sphi >= -halfCarTolerance) 1808 { 1827 { 1809 xi = p.x() + sphi*v.x() ; 1828 xi = p.x() + sphi*v.x() ; 1810 yi = p.y() + sphi*v.y() ; 1829 yi = p.y() + sphi*v.y() ; 1811 1830 1812 // Check intersecting with correc 1831 // Check intersecting with correct half-plane 1813 // (if not -> no intersect) 1832 // (if not -> no intersect) 1814 // 1833 // 1815 if ( (std::fabs(xi)<=kCarToleranc 1834 if ( (std::fabs(xi)<=kCarTolerance) 1816 && (std::fabs(yi)<=kCarToleranc 1835 && (std::fabs(yi)<=kCarTolerance) ) 1817 { 1836 { 1818 sidephi= kSPhi; 1837 sidephi= kSPhi; 1819 if ( ( fSPhi-halfAngTolerance < 1838 if ( ( fSPhi-halfAngTolerance <= vphi ) 1820 && ( fSPhi+fDPhi+halfAngToler 1839 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) ) 1821 { 1840 { 1822 sphi = kInfinity; 1841 sphi = kInfinity; 1823 } 1842 } 1824 } 1843 } 1825 else 1844 else 1826 if ( (yi*cosCPhi-xi*sinCPhi)>=0 ) 1845 if ( (yi*cosCPhi-xi*sinCPhi)>=0 ) 1827 { 1846 { 1828 sphi = kInfinity ; 1847 sphi = kInfinity ; 1829 } 1848 } 1830 else 1849 else 1831 { 1850 { 1832 sidephi = kSPhi ; 1851 sidephi = kSPhi ; 1833 if ( pDistS > -halfCarTolerance 1852 if ( pDistS > -halfCarTolerance ) 1834 { 1853 { 1835 sphi = 0.0 ; // Leave by sphi 1854 sphi = 0.0 ; // Leave by sphi immediately 1836 } 1855 } 1837 } 1856 } 1838 } 1857 } 1839 else 1858 else 1840 { 1859 { 1841 sphi = kInfinity ; 1860 sphi = kInfinity ; 1842 } 1861 } 1843 } 1862 } 1844 else 1863 else 1845 { 1864 { 1846 sphi = kInfinity ; 1865 sphi = kInfinity ; 1847 } 1866 } 1848 1867 1849 if ( compE < 0 ) 1868 if ( compE < 0 ) 1850 { 1869 { 1851 sphi2 = pDistE/compE ; 1870 sphi2 = pDistE/compE ; 1852 1871 1853 // Only check further if < starting 1872 // Only check further if < starting phi intersection 1854 // 1873 // 1855 if ( (sphi2 > -halfCarTolerance) && 1874 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) ) 1856 { 1875 { 1857 xi = p.x() + sphi2*v.x() ; 1876 xi = p.x() + sphi2*v.x() ; 1858 yi = p.y() + sphi2*v.y() ; 1877 yi = p.y() + sphi2*v.y() ; 1859 1878 1860 // Check intersecting with correc 1879 // Check intersecting with correct half-plane 1861 1880 1862 if ( (std::fabs(xi)<=kCarToleranc 1881 if ( (std::fabs(xi)<=kCarTolerance) 1863 && (std::fabs(yi)<=kCarToleranc 1882 && (std::fabs(yi)<=kCarTolerance) ) 1864 { 1883 { 1865 // Leaving via ending phi 1884 // Leaving via ending phi 1866 1885 1867 if( (fSPhi-halfAngTolerance > v << 1886 if(!( (fSPhi-halfAngTolerance <= vphi) 1868 || (fSPhi+fDPhi+halfAngToler << 1887 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) ) 1869 { 1888 { 1870 sidephi = kEPhi ; 1889 sidephi = kEPhi ; 1871 if ( pDistE <= -halfCarTolera 1890 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; } 1872 else 1891 else { sphi = 0.0; } 1873 } 1892 } 1874 } 1893 } 1875 else // Check intersecting with c 1894 else // Check intersecting with correct half-plane 1876 if ( yi*cosCPhi-xi*sinCPhi >= 0 ) 1895 if ( yi*cosCPhi-xi*sinCPhi >= 0 ) 1877 { 1896 { 1878 // Leaving via ending phi 1897 // Leaving via ending phi 1879 1898 1880 sidephi = kEPhi ; 1899 sidephi = kEPhi ; 1881 if ( pDistE <= -halfCarToleranc 1900 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; } 1882 else 1901 else { sphi = 0.0; } 1883 } 1902 } 1884 } 1903 } 1885 } 1904 } 1886 } 1905 } 1887 else 1906 else 1888 { 1907 { 1889 sphi = kInfinity ; 1908 sphi = kInfinity ; 1890 } 1909 } 1891 } 1910 } 1892 else 1911 else 1893 { 1912 { 1894 // On z axis + travel not || to z axis 1913 // On z axis + travel not || to z axis -> if phi of vector direction 1895 // within phi of shape, Step limited by 1914 // within phi of shape, Step limited by rmax, else Step =0 1896 1915 1897 if ( (fSPhi-halfAngTolerance <= vphi) 1916 if ( (fSPhi-halfAngTolerance <= vphi) 1898 && (vphi <= fSPhi+fDPhi+halfAngTolera 1917 && (vphi <= fSPhi+fDPhi+halfAngTolerance) ) 1899 { 1918 { 1900 sphi = kInfinity ; 1919 sphi = kInfinity ; 1901 } 1920 } 1902 else 1921 else 1903 { 1922 { 1904 sidephi = kSPhi ; // arbitrary 1923 sidephi = kSPhi ; // arbitrary 1905 sphi = 0.0 ; 1924 sphi = 0.0 ; 1906 } 1925 } 1907 } 1926 } 1908 if ( sphi < snxt ) // Order intersecttio 1927 if ( sphi < snxt ) // Order intersecttions 1909 { 1928 { 1910 snxt = sphi ; 1929 snxt = sphi ; 1911 side = sidephi ; 1930 side = sidephi ; 1912 } 1931 } 1913 } 1932 } 1914 if ( srd < snxt ) // Order intersections 1933 if ( srd < snxt ) // Order intersections 1915 { 1934 { 1916 snxt = srd ; 1935 snxt = srd ; 1917 side = sider ; 1936 side = sider ; 1918 } 1937 } 1919 if (calcNorm) 1938 if (calcNorm) 1920 { 1939 { 1921 switch(side) 1940 switch(side) 1922 { // Note: returned v 1941 { // Note: returned vector not normalised 1923 case kRMax: // (divide by frmax 1942 case kRMax: // (divide by frmax for unit vector) 1924 xi = p.x() + snxt*v.x() ; 1943 xi = p.x() + snxt*v.x() ; 1925 yi = p.y() + snxt*v.y() ; 1944 yi = p.y() + snxt*v.y() ; 1926 risec = std::sqrt(xi*xi + yi*yi) 1945 risec = std::sqrt(xi*xi + yi*yi)*secRMax ; 1927 *n = G4ThreeVector(xi/risec,y 1946 *n = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ; 1928 *validNorm = true ; 1947 *validNorm = true ; 1929 break ; 1948 break ; 1930 case kRMin: 1949 case kRMin: 1931 *validNorm = false ; // Rmin is inco 1950 *validNorm = false ; // Rmin is inconvex 1932 break ; 1951 break ; 1933 case kSPhi: 1952 case kSPhi: 1934 if ( fDPhi <= pi ) 1953 if ( fDPhi <= pi ) 1935 { 1954 { 1936 *n = G4ThreeVector(sinSPhi, 1955 *n = G4ThreeVector(sinSPhi, -cosSPhi, 0); 1937 *validNorm = true ; 1956 *validNorm = true ; 1938 } 1957 } 1939 else 1958 else 1940 { 1959 { 1941 *validNorm = false ; 1960 *validNorm = false ; 1942 } 1961 } 1943 break ; 1962 break ; 1944 case kEPhi: 1963 case kEPhi: 1945 if ( fDPhi <= pi ) 1964 if ( fDPhi <= pi ) 1946 { 1965 { 1947 *n = G4ThreeVector(-sinEPhi, cosEPh 1966 *n = G4ThreeVector(-sinEPhi, cosEPhi, 0); 1948 *validNorm = true ; 1967 *validNorm = true ; 1949 } 1968 } 1950 else 1969 else 1951 { 1970 { 1952 *validNorm = false ; 1971 *validNorm = false ; 1953 } 1972 } 1954 break ; 1973 break ; 1955 case kPZ: 1974 case kPZ: 1956 *n = G4ThreeVector(0,0,1) ; 1975 *n = G4ThreeVector(0,0,1) ; 1957 *validNorm = true ; 1976 *validNorm = true ; 1958 break ; 1977 break ; 1959 case kMZ: 1978 case kMZ: 1960 *n = G4ThreeVector(0,0,-1) ; 1979 *n = G4ThreeVector(0,0,-1) ; 1961 *validNorm = true ; 1980 *validNorm = true ; 1962 break ; 1981 break ; 1963 default: 1982 default: 1964 G4cout << G4endl ; 1983 G4cout << G4endl ; 1965 DumpInfo(); 1984 DumpInfo(); 1966 std::ostringstream message; 1985 std::ostringstream message; 1967 G4long oldprc = message.precision(16) 1986 G4long oldprc = message.precision(16) ; 1968 message << "Undefined side for valid 1987 message << "Undefined side for valid surface normal to solid." 1969 << G4endl 1988 << G4endl 1970 << "Position:" << G4endl << 1989 << "Position:" << G4endl << G4endl 1971 << "p.x() = " << p.x()/mm < 1990 << "p.x() = " << p.x()/mm << " mm" << G4endl 1972 << "p.y() = " << p.y()/mm < 1991 << "p.y() = " << p.y()/mm << " mm" << G4endl 1973 << "p.z() = " << p.z()/mm < 1992 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl 1974 << "pho at z = " << std::sq 1993 << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm 1975 << " mm" << G4endl << G4endl 1994 << " mm" << G4endl << G4endl ; 1976 if( p.x() != 0. || p.y() != 0.) 1995 if( p.x() != 0. || p.y() != 0.) 1977 { 1996 { 1978 message << "point phi = " << std 1997 message << "point phi = " << std::atan2(p.y(),p.x())/degree 1979 << " degrees" << G4endl << 1998 << " degrees" << G4endl << G4endl ; 1980 } 1999 } 1981 message << "Direction:" << G4endl << 2000 message << "Direction:" << G4endl << G4endl 1982 << "v.x() = " << v.x() << G 2001 << "v.x() = " << v.x() << G4endl 1983 << "v.y() = " << v.y() << G 2002 << "v.y() = " << v.y() << G4endl 1984 << "v.z() = " << v.z() << G 2003 << "v.z() = " << v.z() << G4endl<< G4endl 1985 << "Proposed distance :" << G 2004 << "Proposed distance :" << G4endl<< G4endl 1986 << "snxt = " << snxt/mm << 2005 << "snxt = " << snxt/mm << " mm" << G4endl ; 1987 message.precision(oldprc) ; 2006 message.precision(oldprc) ; 1988 G4Exception("G4Cons::DistanceToOut(p, 2007 G4Exception("G4Cons::DistanceToOut(p,v,..)","GeomSolids1002", 1989 JustWarning, message) ; 2008 JustWarning, message) ; 1990 break ; 2009 break ; 1991 } 2010 } 1992 } 2011 } 1993 if (snxt < halfCarTolerance) { snxt = 0.; 2012 if (snxt < halfCarTolerance) { snxt = 0.; } 1994 2013 1995 return snxt ; 2014 return snxt ; 1996 } 2015 } 1997 2016 1998 ///////////////////////////////////////////// 2017 ////////////////////////////////////////////////////////////////// 1999 // 2018 // 2000 // Calculate distance (<=actual) to closest s 2019 // Calculate distance (<=actual) to closest surface of shape from inside 2001 2020 2002 G4double G4Cons::DistanceToOut(const G4ThreeV 2021 G4double G4Cons::DistanceToOut(const G4ThreeVector& p) const 2003 { 2022 { 2004 G4double safe=0.0, rho, safeR1, safeR2, saf 2023 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi; 2005 G4double tanRMin, secRMin, pRMin; 2024 G4double tanRMin, secRMin, pRMin; 2006 G4double tanRMax, secRMax, pRMax; 2025 G4double tanRMax, secRMax, pRMax; 2007 2026 2008 #ifdef G4CSGDEBUG 2027 #ifdef G4CSGDEBUG 2009 if( Inside(p) == kOutside ) 2028 if( Inside(p) == kOutside ) 2010 { 2029 { 2011 G4int oldprc=G4cout.precision(16) ; 2030 G4int oldprc=G4cout.precision(16) ; 2012 G4cout << G4endl ; 2031 G4cout << G4endl ; 2013 DumpInfo(); 2032 DumpInfo(); 2014 G4cout << "Position:" << G4endl << G4end 2033 G4cout << "Position:" << G4endl << G4endl ; 2015 G4cout << "p.x() = " << p.x()/mm << " m 2034 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; 2016 G4cout << "p.y() = " << p.y()/mm << " m 2035 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; 2017 G4cout << "p.z() = " << p.z()/mm << " m 2036 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; 2018 G4cout << "pho at z = " << std::sqrt( p 2037 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm 2019 << " mm" << G4endl << G4endl ; 2038 << " mm" << G4endl << G4endl ; 2020 if( (p.x() != 0.) || (p.x() != 0.) ) 2039 if( (p.x() != 0.) || (p.x() != 0.) ) 2021 { 2040 { 2022 G4cout << "point phi = " << std::atan 2041 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree 2023 << " degrees" << G4endl << G4end 2042 << " degrees" << G4endl << G4endl ; 2024 } 2043 } 2025 G4cout.precision(oldprc) ; 2044 G4cout.precision(oldprc) ; 2026 G4Exception("G4Cons::DistanceToOut(p)", " 2045 G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002", 2027 JustWarning, "Point p is outs 2046 JustWarning, "Point p is outside !?" ); 2028 } 2047 } 2029 #endif 2048 #endif 2030 2049 2031 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) 2050 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 2032 safeZ = fDz - std::fabs(p.z()) ; 2051 safeZ = fDz - std::fabs(p.z()) ; 2033 2052 2034 if ((fRmin1 != 0.0) || (fRmin2 != 0.0)) << 2053 if (fRmin1 || fRmin2) 2035 { 2054 { 2036 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 2055 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 2037 secRMin = std::sqrt(1.0 + tanRMin*tanRMin 2056 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ; 2038 pRMin = tanRMin*p.z() + (fRmin1 + fRmin 2057 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ; 2039 safeR1 = (rho - pRMin)/secRMin ; 2058 safeR1 = (rho - pRMin)/secRMin ; 2040 } 2059 } 2041 else 2060 else 2042 { 2061 { 2043 safeR1 = kInfinity ; 2062 safeR1 = kInfinity ; 2044 } 2063 } 2045 2064 2046 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 2065 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 2047 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) 2066 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 2048 pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0 2067 pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ; 2049 safeR2 = (pRMax - rho)/secRMax ; 2068 safeR2 = (pRMax - rho)/secRMax ; 2050 2069 2051 if (safeR1 < safeR2) { safe = safeR1; } 2070 if (safeR1 < safeR2) { safe = safeR1; } 2052 else { safe = safeR2; } 2071 else { safe = safeR2; } 2053 if (safeZ < safe) { safe = safeZ ; } 2072 if (safeZ < safe) { safe = safeZ ; } 2054 2073 2055 // Check if phi divided, Calc distances clo 2074 // Check if phi divided, Calc distances closest phi plane 2056 2075 2057 if (!fPhiFullCone) 2076 if (!fPhiFullCone) 2058 { 2077 { 2059 // Above/below central phi of G4Cons? 2078 // Above/below central phi of G4Cons? 2060 2079 2061 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 2080 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 ) 2062 { 2081 { 2063 safePhi = -(p.x()*sinSPhi - p.y()*cosSP 2082 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ; 2064 } 2083 } 2065 else 2084 else 2066 { 2085 { 2067 safePhi = (p.x()*sinEPhi - p.y()*cosEPh 2086 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ; 2068 } 2087 } 2069 if (safePhi < safe) { safe = safePhi; } 2088 if (safePhi < safe) { safe = safePhi; } 2070 } 2089 } 2071 if ( safe < 0 ) { safe = 0; } 2090 if ( safe < 0 ) { safe = 0; } 2072 2091 2073 return safe ; 2092 return safe ; 2074 } 2093 } 2075 2094 2076 ///////////////////////////////////////////// 2095 ////////////////////////////////////////////////////////////////////////// 2077 // 2096 // 2078 // GetEntityType 2097 // GetEntityType 2079 2098 2080 G4GeometryType G4Cons::GetEntityType() const 2099 G4GeometryType G4Cons::GetEntityType() const 2081 { 2100 { 2082 return {"G4Cons"}; << 2101 return G4String("G4Cons"); 2083 } 2102 } 2084 2103 2085 ///////////////////////////////////////////// 2104 ////////////////////////////////////////////////////////////////////////// 2086 // 2105 // 2087 // Make a clone of the object 2106 // Make a clone of the object 2088 // 2107 // 2089 G4VSolid* G4Cons::Clone() const 2108 G4VSolid* G4Cons::Clone() const 2090 { 2109 { 2091 return new G4Cons(*this); 2110 return new G4Cons(*this); 2092 } 2111 } 2093 2112 2094 ///////////////////////////////////////////// 2113 ////////////////////////////////////////////////////////////////////////// 2095 // 2114 // 2096 // Stream object contents to an output stream 2115 // Stream object contents to an output stream 2097 2116 2098 std::ostream& G4Cons::StreamInfo(std::ostream 2117 std::ostream& G4Cons::StreamInfo(std::ostream& os) const 2099 { 2118 { 2100 G4long oldprc = os.precision(16); 2119 G4long oldprc = os.precision(16); 2101 os << "------------------------------------ 2120 os << "-----------------------------------------------------------\n" 2102 << " *** Dump for solid - " << GetNam 2121 << " *** Dump for solid - " << GetName() << " ***\n" 2103 << " ================================ 2122 << " ===================================================\n" 2104 << " Solid type: G4Cons\n" 2123 << " Solid type: G4Cons\n" 2105 << " Parameters: \n" 2124 << " Parameters: \n" 2106 << " inside -fDz radius: " << fRmin1 2125 << " inside -fDz radius: " << fRmin1/mm << " mm \n" 2107 << " outside -fDz radius: " << fRmax1 2126 << " outside -fDz radius: " << fRmax1/mm << " mm \n" 2108 << " inside +fDz radius: " << fRmin2 2127 << " inside +fDz radius: " << fRmin2/mm << " mm \n" 2109 << " outside +fDz radius: " << fRmax2 2128 << " outside +fDz radius: " << fRmax2/mm << " mm \n" 2110 << " half length in Z : " << fDz/mm 2129 << " half length in Z : " << fDz/mm << " mm \n" 2111 << " starting angle of segment: " << f 2130 << " starting angle of segment: " << fSPhi/degree << " degrees \n" 2112 << " delta angle of segment : " << f 2131 << " delta angle of segment : " << fDPhi/degree << " degrees \n" 2113 << "------------------------------------ 2132 << "-----------------------------------------------------------\n"; 2114 os.precision(oldprc); 2133 os.precision(oldprc); 2115 2134 2116 return os; 2135 return os; 2117 } 2136 } 2118 2137 2119 2138 2120 2139 2121 ///////////////////////////////////////////// 2140 ///////////////////////////////////////////////////////////////////////// 2122 // 2141 // 2123 // GetPointOnSurface 2142 // GetPointOnSurface 2124 2143 2125 G4ThreeVector G4Cons::GetPointOnSurface() con 2144 G4ThreeVector G4Cons::GetPointOnSurface() const 2126 { 2145 { 2127 // declare working variables 2146 // declare working variables 2128 // 2147 // 2129 G4double rone = (fRmax1-fRmax2)/(2.*fDz); 2148 G4double rone = (fRmax1-fRmax2)/(2.*fDz); 2130 G4double rtwo = (fRmin1-fRmin2)/(2.*fDz); 2149 G4double rtwo = (fRmin1-fRmin2)/(2.*fDz); 2131 G4double qone = (fRmax1 == fRmax2) ? 0. : f 2150 G4double qone = (fRmax1 == fRmax2) ? 0. : fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2); 2132 G4double qtwo = (fRmin1 == fRmin2) ? 0. : f 2151 G4double qtwo = (fRmin1 == fRmin2) ? 0. : fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2); 2133 2152 2134 G4double slin = std::hypot(fRmin1-fRmin2, 2153 G4double slin = std::hypot(fRmin1-fRmin2, 2.*fDz); 2135 G4double slout = std::hypot(fRmax1-fRmax2, 2154 G4double slout = std::hypot(fRmax1-fRmax2, 2.*fDz); 2136 G4double Aone = 0.5*fDPhi*(fRmax2 + fRmax 2155 G4double Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*slout; // outer surface 2137 G4double Atwo = 0.5*fDPhi*(fRmin2 + fRmin 2156 G4double Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*slin; // inner surface 2138 G4double Athree = 0.5*fDPhi*(fRmax1*fRmax1- 2157 G4double Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1); // base at -Dz 2139 G4double Afour = 0.5*fDPhi*(fRmax2*fRmax2- 2158 G4double Afour = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2); // base at +Dz 2140 G4double Afive = fDz*(fRmax1-fRmin1+fRmax2 2159 G4double Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2); // phi section 2141 2160 2142 G4double phi = G4RandFlat::shoot(fSPhi,f 2161 G4double phi = G4RandFlat::shoot(fSPhi,fSPhi+fDPhi); 2143 G4double cosu = std::cos(phi); 2162 G4double cosu = std::cos(phi); 2144 G4double sinu = std::sin(phi); 2163 G4double sinu = std::sin(phi); 2145 G4double rRand1 = GetRadiusInRing(fRmin1, f 2164 G4double rRand1 = GetRadiusInRing(fRmin1, fRmax1); 2146 G4double rRand2 = GetRadiusInRing(fRmin2, f 2165 G4double rRand2 = GetRadiusInRing(fRmin2, fRmax2); 2147 2166 2148 if ( (fSPhi == 0.) && fPhiFullCone ) { Afi 2167 if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; } 2149 G4double chose = G4RandFlat::shoot(0.,Aone 2168 G4double chose = G4RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive); 2150 2169 2151 if( (chose >= 0.) && (chose < Aone) ) // ou 2170 if( (chose >= 0.) && (chose < Aone) ) // outer surface 2152 { 2171 { 2153 if(fRmax1 != fRmax2) 2172 if(fRmax1 != fRmax2) 2154 { 2173 { 2155 G4double zRand = G4RandFlat::shoot(-1.* 2174 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz); 2156 return { rone*cosu*(qone-zRand), rone*s << 2175 return G4ThreeVector (rone*cosu*(qone-zRand), >> 2176 rone*sinu*(qone-zRand), zRand); 2157 } 2177 } 2158 else 2178 else 2159 { 2179 { 2160 return { fRmax1*cosu, fRmax2*sinu, G4Ra << 2180 return G4ThreeVector(fRmax1*cosu, fRmax2*sinu, >> 2181 G4RandFlat::shoot(-1.*fDz,fDz)); 2161 } 2182 } 2162 } 2183 } 2163 else if( (chose >= Aone) && (chose < Aone + 2184 else if( (chose >= Aone) && (chose < Aone + Atwo) ) // inner surface 2164 { 2185 { 2165 if(fRmin1 != fRmin2) 2186 if(fRmin1 != fRmin2) 2166 { 2187 { 2167 G4double zRand = G4RandFlat::shoot(-1.* 2188 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz); 2168 return { rtwo*cosu*(qtwo-zRand), rtwo*s << 2189 return G4ThreeVector (rtwo*cosu*(qtwo-zRand), >> 2190 rtwo*sinu*(qtwo-zRand), zRand); 2169 } 2191 } 2170 else 2192 else 2171 { 2193 { 2172 return { fRmin1*cosu, fRmin2*sinu, G4Ra << 2194 return G4ThreeVector(fRmin1*cosu, fRmin2*sinu, >> 2195 G4RandFlat::shoot(-1.*fDz,fDz)); 2173 } 2196 } 2174 } 2197 } 2175 else if( (chose >= Aone + Atwo) && (chose < 2198 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) ) // base at -Dz 2176 { 2199 { 2177 return {rRand1*cosu, rRand1*sinu, -1*fDz} << 2200 return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz); 2178 } 2201 } 2179 else if( (chose >= Aone + Atwo + Athree) 2202 else if( (chose >= Aone + Atwo + Athree) 2180 && (chose < Aone + Atwo + Athree + Af 2203 && (chose < Aone + Atwo + Athree + Afour) ) // base at +Dz 2181 { 2204 { 2182 return { rRand2*cosu, rRand2*sinu, fDz }; << 2205 return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz); 2183 } 2206 } 2184 else if( (chose >= Aone + Atwo + Athree + A 2207 else if( (chose >= Aone + Atwo + Athree + Afour) // SPhi section 2185 && (chose < Aone + Atwo + Athree + Af 2208 && (chose < Aone + Atwo + Athree + Afour + Afive) ) 2186 { 2209 { 2187 G4double zRand = G4RandFlat::shoot(-1.*f 2210 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz); 2188 rRand1 = G4RandFlat::shoot(fRmin2-((zRand 2211 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2), 2189 fRmax2-((zRand 2212 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2)); 2190 return { rRand1*cosSPhi, rRand1*sinSPhi, << 2213 return G4ThreeVector (rRand1*cosSPhi, >> 2214 rRand1*sinSPhi, zRand); 2191 } 2215 } 2192 else // SPhi+DPhi section 2216 else // SPhi+DPhi section 2193 { 2217 { 2194 G4double zRand = G4RandFlat::shoot(-1.*f 2218 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz); 2195 rRand1 = G4RandFlat::shoot(fRmin2-((zRand 2219 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2), 2196 fRmax2-((zRand 2220 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2)); 2197 return { rRand1*cosEPhi, rRand1*sinEPhi, << 2221 return G4ThreeVector (rRand1*cosEPhi, >> 2222 rRand1*sinEPhi, zRand); 2198 } 2223 } 2199 } 2224 } 2200 2225 2201 ///////////////////////////////////////////// 2226 ////////////////////////////////////////////////////////////////////////// 2202 // 2227 // 2203 // Methods for visualisation 2228 // Methods for visualisation 2204 2229 2205 void G4Cons::DescribeYourselfTo (G4VGraphicsS 2230 void G4Cons::DescribeYourselfTo (G4VGraphicsScene& scene) const 2206 { 2231 { 2207 scene.AddSolid (*this); 2232 scene.AddSolid (*this); 2208 } 2233 } 2209 2234 2210 G4Polyhedron* G4Cons::CreatePolyhedron () con 2235 G4Polyhedron* G4Cons::CreatePolyhedron () const 2211 { 2236 { 2212 return new G4PolyhedronCons(fRmin1,fRmax1,f 2237 return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi); 2213 } 2238 } 2214 2239 2215 #endif 2240 #endif 2216 2241