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