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