| 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 // G4TwistTubsHypeSide implementation <<
27 // 26 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu << 27 //
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), << 28 //
30 // from original version in Jupi << 29 // --------------------------------------------------------------------
>> 30 // GEANT 4 class source file
>> 31 //
>> 32 //
>> 33 // G4TwistTubsHypeSide.cc
>> 34 //
>> 35 // Author:
>> 36 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
>> 37 //
>> 38 // History:
>> 39 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
>> 40 // from original version in Jupiter-2.5.02 application.
31 // ------------------------------------------- 41 // --------------------------------------------------------------------
32 42
33 #include "G4TwistTubsHypeSide.hh" 43 #include "G4TwistTubsHypeSide.hh"
34 #include "G4PhysicalConstants.hh" 44 #include "G4PhysicalConstants.hh"
35 #include "G4GeometryTolerance.hh" 45 #include "G4GeometryTolerance.hh"
36 46
37 //============================================ 47 //=====================================================================
38 //* constructors ----------------------------- 48 //* constructors ------------------------------------------------------
39 49
40 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const << 50 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const G4String &name,
41 const << 51 const G4RotationMatrix &rot,
42 const << 52 const G4ThreeVector &tlate,
43 const 53 const G4int handedness,
44 const 54 const G4double kappa,
45 const 55 const G4double tanstereo,
46 const 56 const G4double r0,
47 const 57 const EAxis axis0,
48 const 58 const EAxis axis1,
49 59 G4double axis0min,
50 60 G4double axis1min,
51 61 G4double axis0max,
52 62 G4double axis1max )
53 : G4VTwistSurface(name, rot, tlate, handedne 63 : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
54 axis0min, axis1min, axis0m << 64 axis0min, axis1min, axis0max, axis1max),
55 fKappa(kappa), fTanStereo(tanstereo), 65 fKappa(kappa), fTanStereo(tanstereo),
56 fTan2Stereo(tanstereo*tanstereo), fR0(r0), 66 fTan2Stereo(tanstereo*tanstereo), fR0(r0), fR02(r0*r0), fDPhi(twopi)
57 { 67 {
58 if ( (axis0 == kZAxis) && (axis1 == kPhi) ) 68 if ( (axis0 == kZAxis) && (axis1 == kPhi) )
59 { 69 {
60 G4Exception("G4TwistTubsHypeSide::G4Twis 70 G4Exception("G4TwistTubsHypeSide::G4TwistTubsHypeSide()",
61 "GeomSolids0002", FatalError 71 "GeomSolids0002", FatalErrorInArgument,
62 "Should swap axis0 and axis1 72 "Should swap axis0 and axis1!");
63 } 73 }
>> 74
64 fInside.gp.set(kInfinity, kInfinity, kInfin 75 fInside.gp.set(kInfinity, kInfinity, kInfinity);
65 fInside.inside = kOutside; 76 fInside.inside = kOutside;
66 fIsValidNorm = false; 77 fIsValidNorm = false;
>> 78
67 SetCorners(); 79 SetCorners();
68 SetBoundaries(); 80 SetBoundaries();
>> 81
69 } 82 }
70 83
71 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const << 84 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const G4String &name,
72 << 85 G4double EndInnerRadius[2],
73 << 86 G4double EndOuterRadius[2],
74 << 87 G4double DPhi,
75 << 88 G4double EndPhi[2],
76 << 89 G4double EndZ[2],
77 << 90 G4double InnerRadius,
78 << 91 G4double OuterRadius,
79 << 92 G4double Kappa,
80 << 93 G4double TanInnerStereo,
81 << 94 G4double TanOuterStereo,
82 << 95 G4int handedness)
83 : G4VTwistSurface(name) 96 : G4VTwistSurface(name)
84 { 97 {
85 98
86 fHandedness = handedness; // +z = +ve, -z 99 fHandedness = handedness; // +z = +ve, -z = -ve
87 fAxis[0] = kPhi; 100 fAxis[0] = kPhi;
88 fAxis[1] = kZAxis; 101 fAxis[1] = kZAxis;
89 fAxisMin[0] = kInfinity; // we cann 102 fAxisMin[0] = kInfinity; // we cannot fix boundary min of Phi,
90 fAxisMax[0] = kInfinity; // because 103 fAxisMax[0] = kInfinity; // because it depends on z.
91 fAxisMin[1] = EndZ[0]; 104 fAxisMin[1] = EndZ[0];
92 fAxisMax[1] = EndZ[1]; 105 fAxisMax[1] = EndZ[1];
93 fKappa = Kappa; 106 fKappa = Kappa;
94 fDPhi = DPhi ; 107 fDPhi = DPhi ;
95 108
96 if (handedness < 0) // inner hyperbolic su << 109 if (handedness < 0) { // inner hyperbolic surface
97 { <<
98 fTanStereo = TanInnerStereo; 110 fTanStereo = TanInnerStereo;
99 fR0 = InnerRadius; 111 fR0 = InnerRadius;
100 } << 112 } else { // outer hyperbolic surface
101 else // outer hyperbolic su <<
102 { <<
103 fTanStereo = TanOuterStereo; 113 fTanStereo = TanOuterStereo;
104 fR0 = OuterRadius; 114 fR0 = OuterRadius;
105 } 115 }
106 fTan2Stereo = fTanStereo * fTanStereo; 116 fTan2Stereo = fTanStereo * fTanStereo;
107 fR02 = fR0 * fR0; 117 fR02 = fR0 * fR0;
108 118
109 fTrans.set(0, 0, 0); 119 fTrans.set(0, 0, 0);
110 fIsValidNorm = false; 120 fIsValidNorm = false;
111 121
112 fInside.gp.set(kInfinity, kInfinity, kInfin 122 fInside.gp.set(kInfinity, kInfinity, kInfinity);
113 fInside.inside = kOutside; 123 fInside.inside = kOutside;
114 124
115 SetCorners(EndInnerRadius, EndOuterRadius, 125 SetCorners(EndInnerRadius, EndOuterRadius, DPhi, EndPhi, EndZ) ;
116 126
117 SetBoundaries(); 127 SetBoundaries();
118 } 128 }
119 129
120 //============================================ 130 //=====================================================================
121 //* Fake default constructor ----------------- 131 //* Fake default constructor ------------------------------------------
122 132
123 G4TwistTubsHypeSide::G4TwistTubsHypeSide( __vo 133 G4TwistTubsHypeSide::G4TwistTubsHypeSide( __void__& a )
124 : G4VTwistSurface(a) << 134 : G4VTwistSurface(a), fKappa(0.), fTanStereo(0.), fTan2Stereo(0.),
>> 135 fR0(0.), fR02(0.), fDPhi(0.)
125 { 136 {
126 } 137 }
127 138
128 //============================================ 139 //=====================================================================
129 //* destructor ------------------------------- 140 //* destructor --------------------------------------------------------
130 141
131 G4TwistTubsHypeSide::~G4TwistTubsHypeSide() = << 142 G4TwistTubsHypeSide::~G4TwistTubsHypeSide()
>> 143 {
>> 144 }
132 145
133 //============================================ 146 //=====================================================================
134 //* GetNormal -------------------------------- 147 //* GetNormal ---------------------------------------------------------
135 148
136 G4ThreeVector G4TwistTubsHypeSide::GetNormal(c << 149 G4ThreeVector G4TwistTubsHypeSide::GetNormal(const G4ThreeVector &tmpxx,
137 150 G4bool isGlobal)
138 { 151 {
139 // GetNormal returns a normal vector at a s 152 // GetNormal returns a normal vector at a surface (or very close
140 // to surface) point at tmpxx. 153 // to surface) point at tmpxx.
141 // If isGlobal=true, it returns the normal 154 // If isGlobal=true, it returns the normal in global coordinate.
>> 155 //
142 156
143 G4ThreeVector xx; 157 G4ThreeVector xx;
144 if (isGlobal) << 158 if (isGlobal) {
145 { <<
146 xx = ComputeLocalPoint(tmpxx); 159 xx = ComputeLocalPoint(tmpxx);
147 if ((xx - fCurrentNormal.p).mag() < 0.5 << 160 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
148 { <<
149 return ComputeGlobalDirection(fCurren 161 return ComputeGlobalDirection(fCurrentNormal.normal);
150 } 162 }
151 } << 163 } else {
152 else <<
153 { <<
154 xx = tmpxx; 164 xx = tmpxx;
155 if (xx == fCurrentNormal.p) << 165 if (xx == fCurrentNormal.p) {
156 { <<
157 return fCurrentNormal.normal; 166 return fCurrentNormal.normal;
158 } 167 }
159 } 168 }
160 169
161 fCurrentNormal.p = xx; 170 fCurrentNormal.p = xx;
162 171
163 G4ThreeVector normal( xx.x(), xx.y(), -xx.z 172 G4ThreeVector normal( xx.x(), xx.y(), -xx.z() * fTan2Stereo);
164 normal *= fHandedness; 173 normal *= fHandedness;
165 normal = normal.unit(); 174 normal = normal.unit();
166 175
167 if (isGlobal) << 176 if (isGlobal) {
168 { <<
169 fCurrentNormal.normal = ComputeLocalDire 177 fCurrentNormal.normal = ComputeLocalDirection(normal);
170 } << 178 } else {
171 else <<
172 { <<
173 fCurrentNormal.normal = normal; 179 fCurrentNormal.normal = normal;
174 } 180 }
175 return fCurrentNormal.normal; 181 return fCurrentNormal.normal;
176 } 182 }
177 183
178 //============================================ 184 //=====================================================================
179 //* Inside() --------------------------------- 185 //* Inside() ----------------------------------------------------------
180 186
181 EInside G4TwistTubsHypeSide::Inside(const G4Th << 187 EInside G4TwistTubsHypeSide::Inside(const G4ThreeVector &gp)
182 { 188 {
183 // Inside returns 189 // Inside returns
184 const G4double halftol 190 const G4double halftol
185 = 0.5 * G4GeometryTolerance::GetInstance( 191 = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
186 192
187 if (fInside.gp == gp) << 193 if (fInside.gp == gp) {
188 { <<
189 return fInside.inside; 194 return fInside.inside;
190 } 195 }
191 fInside.gp = gp; 196 fInside.gp = gp;
192 197
193 G4ThreeVector p = ComputeLocalPoint(gp); 198 G4ThreeVector p = ComputeLocalPoint(gp);
194 199
195 200
196 if (p.mag() < DBL_MIN) << 201 if (p.mag() < DBL_MIN) {
197 { <<
198 fInside.inside = kOutside; 202 fInside.inside = kOutside;
199 return fInside.inside; 203 return fInside.inside;
200 } 204 }
201 205
202 G4double rhohype = GetRhoAtPZ(p); 206 G4double rhohype = GetRhoAtPZ(p);
203 G4double distanceToOut = fHandedness * (rho 207 G4double distanceToOut = fHandedness * (rhohype - p.getRho());
204 // +ve : inside 208 // +ve : inside
205 209
206 if (distanceToOut < -halftol) << 210 if (distanceToOut < -halftol) {
207 { << 211
208 fInside.inside = kOutside; 212 fInside.inside = kOutside;
209 } << 213
210 else << 214 } else {
211 { << 215
212 G4int areacode = GetAreaCode(p); 216 G4int areacode = GetAreaCode(p);
213 if (IsOutside(areacode)) << 217 if (IsOutside(areacode)) {
214 { <<
215 fInside.inside = kOutside; 218 fInside.inside = kOutside;
216 } << 219 } else if (IsBoundary(areacode)) {
217 else if (IsBoundary(areacode)) <<
218 { <<
219 fInside.inside = kSurface; 220 fInside.inside = kSurface;
220 } << 221 } else if (IsInside(areacode)) {
221 else if (IsInside(areacode)) << 222 if (distanceToOut <= halftol) {
222 { <<
223 if (distanceToOut <= halftol) <<
224 { <<
225 fInside.inside = kSurface; 223 fInside.inside = kSurface;
226 } << 224 } else {
227 else <<
228 { <<
229 fInside.inside = kInside; 225 fInside.inside = kInside;
230 } 226 }
231 } << 227 } else {
232 else <<
233 { <<
234 G4cout << "WARNING - G4TwistTubsHypeS 228 G4cout << "WARNING - G4TwistTubsHypeSide::Inside()" << G4endl
235 << " Invalid option ! 229 << " Invalid option !" << G4endl
236 << " name, areacode, 230 << " name, areacode, distanceToOut = "
237 << GetName() << ", " << std::h << 231 << GetName() << ", " << std::hex << areacode << std::dec << ", "
238 << std::dec << ", " << distanc << 232 << distanceToOut << G4endl;
239 } 233 }
240 } 234 }
241 << 235
242 return fInside.inside; 236 return fInside.inside;
243 } 237 }
244 238
245 //============================================ 239 //=====================================================================
246 //* DistanceToSurface ------------------------ 240 //* DistanceToSurface -------------------------------------------------
247 241
248 G4int G4TwistTubsHypeSide::DistanceToSurface(c << 242 G4int G4TwistTubsHypeSide::DistanceToSurface(const G4ThreeVector &gp,
249 c << 243 const G4ThreeVector &gv,
250 244 G4ThreeVector gxx[],
251 245 G4double distance[],
252 246 G4int areacode[],
253 247 G4bool isvalid[],
254 248 EValidate validate)
255 { 249 {
>> 250 //
256 // Decide if and where a line intersects wi 251 // Decide if and where a line intersects with a hyperbolic
257 // surface (of infinite extent) 252 // surface (of infinite extent)
258 // 253 //
259 // Arguments: 254 // Arguments:
260 // p - (in) Point on trajectory 255 // p - (in) Point on trajectory
261 // v - (in) Vector along trajecto 256 // v - (in) Vector along trajectory
262 // r2 - (in) Square of radius at z 257 // r2 - (in) Square of radius at z = 0
263 // tan2phi - (in) std::tan(stereo)**2 258 // tan2phi - (in) std::tan(stereo)**2
264 // s - (out) Up to two points of 259 // s - (out) Up to two points of intersection, where the
265 // intersection point i 260 // intersection point is p + s*v, and if there are
266 // two intersections, s 261 // two intersections, s[0] < s[1]. May be negative.
267 // Returns: 262 // Returns:
268 // The number of intersections. If 0, t 263 // The number of intersections. If 0, the trajectory misses.
269 // 264 //
270 // 265 //
271 // Equation of a line: 266 // Equation of a line:
272 // 267 //
273 // x = x0 + s*tx y = y0 + s*ty 268 // x = x0 + s*tx y = y0 + s*ty z = z0 + s*tz
274 // 269 //
275 // Equation of a hyperbolic surface: 270 // Equation of a hyperbolic surface:
276 // 271 //
277 // x**2 + y**2 = r**2 + (z*tanPhi)**2 272 // x**2 + y**2 = r**2 + (z*tanPhi)**2
278 // 273 //
279 // Solution is quadratic: 274 // Solution is quadratic:
280 // 275 //
281 // a*s**2 + b*s + c = 0 276 // a*s**2 + b*s + c = 0
282 // 277 //
283 // where: 278 // where:
284 // 279 //
285 // a = tx**2 + ty**2 - (tz*tanPhi)**2 280 // a = tx**2 + ty**2 - (tz*tanPhi)**2
286 // 281 //
287 // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 282 // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
288 // 283 //
289 // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)* 284 // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
290 // 285 //
291 286
292 fCurStatWithV.ResetfDone(validate, &gp, &gv 287 fCurStatWithV.ResetfDone(validate, &gp, &gv);
293 288
294 if (fCurStatWithV.IsDone()) << 289 if (fCurStatWithV.IsDone()) {
295 { << 290 G4int i;
296 for (G4int i=0; i<fCurStatWithV.GetNXX() << 291 for (i=0; i<fCurStatWithV.GetNXX(); i++) {
297 { <<
298 gxx[i] = fCurStatWithV.GetXX(i); 292 gxx[i] = fCurStatWithV.GetXX(i);
299 distance[i] = fCurStatWithV.GetDistan 293 distance[i] = fCurStatWithV.GetDistance(i);
300 areacode[i] = fCurStatWithV.GetAreaco 294 areacode[i] = fCurStatWithV.GetAreacode(i);
301 isvalid[i] = fCurStatWithV.IsValid(i 295 isvalid[i] = fCurStatWithV.IsValid(i);
302 } 296 }
303 return fCurStatWithV.GetNXX(); 297 return fCurStatWithV.GetNXX();
304 } << 298 } else {
305 else // initialize << 299 // initialize
306 { << 300 G4int i;
307 for (auto i=0; i<2; ++i) << 301 for (i=0; i<2; i++) {
308 { <<
309 distance[i] = kInfinity; 302 distance[i] = kInfinity;
310 areacode[i] = sOutside; 303 areacode[i] = sOutside;
311 isvalid[i] = false; 304 isvalid[i] = false;
312 gxx[i].set(kInfinity, kInfinity, kInf 305 gxx[i].set(kInfinity, kInfinity, kInfinity);
313 } 306 }
314 } 307 }
315 308
316 G4ThreeVector p = ComputeLocalPoint(gp); 309 G4ThreeVector p = ComputeLocalPoint(gp);
317 G4ThreeVector v = ComputeLocalDirection(gv) 310 G4ThreeVector v = ComputeLocalDirection(gv);
318 G4ThreeVector xx[2]; 311 G4ThreeVector xx[2];
319 312
320 // 313 //
321 // special case! p is on origin. 314 // special case! p is on origin.
322 // 315 //
323 316
324 if (p.mag() == 0) << 317 if (p.mag() == 0) {
325 { <<
326 // p is origin. 318 // p is origin.
327 // unique solution of 2-dimension questi 319 // unique solution of 2-dimension question in r-z plane
328 // Equations: 320 // Equations:
329 // r^2 = fR02 + z^2*fTan2Stere0 321 // r^2 = fR02 + z^2*fTan2Stere0
330 // r = beta*z 322 // r = beta*z
331 // where 323 // where
332 // beta = vrho / vz 324 // beta = vrho / vz
333 // Solution (z value of intersection poi 325 // Solution (z value of intersection point):
334 // xxz = +- std::sqrt (fR02 / (beta^2 326 // xxz = +- std::sqrt (fR02 / (beta^2 - fTan2Stereo))
335 // 327 //
336 328
337 G4double vz = v.z(); 329 G4double vz = v.z();
338 G4double absvz = std::fabs(vz); 330 G4double absvz = std::fabs(vz);
339 G4double vrho = v.getRho(); 331 G4double vrho = v.getRho();
340 G4double vslope = vrho/vz; 332 G4double vslope = vrho/vz;
341 G4double vslope2 = vslope * vslope; 333 G4double vslope2 = vslope * vslope;
342 if (vrho == 0 || (vrho/absvz) <= (absvz* << 334 if (vrho == 0 || (vrho/absvz) <= (absvz*std::fabs(fTanStereo)/absvz)) {
343 { <<
344 // vz/vrho is bigger than slope of as 335 // vz/vrho is bigger than slope of asymptonic line
345 distance[0] = kInfinity; 336 distance[0] = kInfinity;
346 fCurStatWithV.SetCurrentStatus(0, gxx 337 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
347 isvali 338 isvalid[0], 0, validate, &gp, &gv);
348 return 0; 339 return 0;
349 } 340 }
350 341
351 if (vz != 0.0) << 342 if (vz) {
352 { <<
353 G4double xxz = std::sqrt(fR02 / (vsl 343 G4double xxz = std::sqrt(fR02 / (vslope2 - fTan2Stereo))
354 * (vz / std::fabs(vz)) 344 * (vz / std::fabs(vz)) ;
355 G4double t = xxz / vz; 345 G4double t = xxz / vz;
356 xx[0].set(t*v.x(), t*v.y(), xxz); 346 xx[0].set(t*v.x(), t*v.y(), xxz);
357 } << 347 } else {
358 else <<
359 { <<
360 // p.z = 0 && v.z =0 348 // p.z = 0 && v.z =0
361 xx[0].set(v.x()*fR0, v.y()*fR0, 0); 349 xx[0].set(v.x()*fR0, v.y()*fR0, 0); // v is a unit vector.
362 } 350 }
363 distance[0] = xx[0].mag(); 351 distance[0] = xx[0].mag();
364 gxx[0] = ComputeGlobalPoint(xx[0]); 352 gxx[0] = ComputeGlobalPoint(xx[0]);
365 353
366 if (validate == kValidateWithTol) << 354 if (validate == kValidateWithTol) {
367 { <<
368 areacode[0] = GetAreaCode(xx[0]); 355 areacode[0] = GetAreaCode(xx[0]);
369 if (!IsOutside(areacode[0])) << 356 if (!IsOutside(areacode[0])) {
370 { <<
371 if (distance[0] >= 0) isvalid[0] = 357 if (distance[0] >= 0) isvalid[0] = true;
372 } 358 }
373 } << 359 } else if (validate == kValidateWithoutTol) {
374 else if (validate == kValidateWithoutTol <<
375 { <<
376 areacode[0] = GetAreaCode(xx[0], fals 360 areacode[0] = GetAreaCode(xx[0], false);
377 if (IsInside(areacode[0])) << 361 if (IsInside(areacode[0])) {
378 { <<
379 if (distance[0] >= 0) isvalid[0] = 362 if (distance[0] >= 0) isvalid[0] = true;
380 } 363 }
381 } << 364 } else { // kDontValidate
382 else // kDontValidate <<
383 { <<
384 areacode[0] = sInside; 365 areacode[0] = sInside;
385 if (distance[0] >= 0) isvalid[0] = 366 if (distance[0] >= 0) isvalid[0] = true;
386 } 367 }
387 368
388 fCurStatWithV.SetCurrentStatus(0, gxx[0] 369 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
389 isvalid[0 << 370 isvalid[0], 1, validate, &gp, &gv);
390 return 1; 371 return 1;
391 } 372 }
392 373
393 // 374 //
394 // special case end. 375 // special case end.
395 // 376 //
396 377
397 G4double a = v.x()*v.x() + v.y()*v.y() - v. 378 G4double a = v.x()*v.x() + v.y()*v.y() - v.z()*v.z()*fTan2Stereo;
398 G4double b = 2.0 << 379 G4double b = 2.0 * ( p.x() * v.x() + p.y() * v.y() - p.z() * v.z() * fTan2Stereo );
399 * ( p.x() * v.x() + p.y() * v.y( <<
400 G4double c = p.x()*p.x() + p.y()*p.y() - fR 380 G4double c = p.x()*p.x() + p.y()*p.y() - fR02 - p.z()*p.z()*fTan2Stereo;
401 G4double D = b*b - 4*a*c; //discri 381 G4double D = b*b - 4*a*c; //discriminant
402 G4int vout = 0; 382 G4int vout = 0;
403 383
404 if (std::fabs(a) < DBL_MIN) << 384 if (std::fabs(a) < DBL_MIN) {
405 { << 385 if (std::fabs(b) > DBL_MIN) { // single solution
406 if (std::fabs(b) > DBL_MIN) / << 386
407 { <<
408 distance[0] = -c/b; 387 distance[0] = -c/b;
409 xx[0] = p + distance[0]*v; 388 xx[0] = p + distance[0]*v;
410 gxx[0] = ComputeGlobalPoint(xx[0]); 389 gxx[0] = ComputeGlobalPoint(xx[0]);
411 390
412 if (validate == kValidateWithTol) << 391 if (validate == kValidateWithTol) {
413 { <<
414 areacode[0] = GetAreaCode(xx[0]); 392 areacode[0] = GetAreaCode(xx[0]);
415 if (!IsOutside(areacode[0])) << 393 if (!IsOutside(areacode[0])) {
416 { <<
417 if (distance[0] >= 0) isvalid[0 394 if (distance[0] >= 0) isvalid[0] = true;
418 } 395 }
419 } << 396 } else if (validate == kValidateWithoutTol) {
420 else if (validate == kValidateWithout <<
421 { <<
422 areacode[0] = GetAreaCode(xx[0], f 397 areacode[0] = GetAreaCode(xx[0], false);
423 if (IsInside(areacode[0])) << 398 if (IsInside(areacode[0])) {
424 { <<
425 if (distance[0] >= 0) isvalid[0 399 if (distance[0] >= 0) isvalid[0] = true;
426 } 400 }
427 } << 401 } else { // kDontValidate
428 else // kDontValidate <<
429 { <<
430 areacode[0] = sInside; 402 areacode[0] = sInside;
431 if (distance[0] >= 0) isvalid[0 403 if (distance[0] >= 0) isvalid[0] = true;
432 } 404 }
433 405
434 fCurStatWithV.SetCurrentStatus(0, gxx 406 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
435 isvali 407 isvalid[0], 1, validate, &gp, &gv);
436 vout = 1; 408 vout = 1;
437 } << 409
438 else << 410 } else {
439 { << 411 // if a=b=0 and c != 0, p is origin and v is parallel to asymptotic line.
440 // if a=b=0 and c != 0, p is origin an << 412 // if a=b=c=0, p is on surface and v is paralell to stereo wire.
441 // if a=b=c=0, p is on surface and v i << 413 // return distance = infinity.
442 // return distance = infinity <<
443 414
444 fCurStatWithV.SetCurrentStatus(0, gxx 415 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
445 isvali 416 isvalid[0], 0, validate, &gp, &gv);
>> 417
446 vout = 0; 418 vout = 0;
447 } 419 }
448 } << 420
449 else if (D > DBL_MIN) // double so << 421 } else if (D > DBL_MIN) { // double solutions
450 { << 422
451 D = std::sqrt(D); 423 D = std::sqrt(D);
452 G4double factor = 0.5/a; 424 G4double factor = 0.5/a;
453 G4double tmpdist[2] = {kInfinity, k 425 G4double tmpdist[2] = {kInfinity, kInfinity};
454 G4ThreeVector tmpxx[2] ; 426 G4ThreeVector tmpxx[2] ;
455 G4int tmpareacode[2] = {sOutside 427 G4int tmpareacode[2] = {sOutside, sOutside};
456 G4bool tmpisvalid[2] = {false, f 428 G4bool tmpisvalid[2] = {false, false};
>> 429 G4int i;
457 430
458 for (auto i=0; i<2; ++i) << 431 for (i=0; i<2; i++) {
459 { <<
460 tmpdist[i] = factor*(-b - D); 432 tmpdist[i] = factor*(-b - D);
461 D = -D; 433 D = -D;
462 tmpxx[i] = p + tmpdist[i]*v; 434 tmpxx[i] = p + tmpdist[i]*v;
463 435
464 if (validate == kValidateWithTol) << 436 if (validate == kValidateWithTol) {
465 { <<
466 tmpareacode[i] = GetAreaCode(tmpxx 437 tmpareacode[i] = GetAreaCode(tmpxx[i]);
467 if (!IsOutside(tmpareacode[i])) << 438 if (!IsOutside(tmpareacode[i])) {
468 { <<
469 if (tmpdist[i] >= 0) tmpisvalid 439 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
470 continue; 440 continue;
471 } 441 }
472 } << 442 } else if (validate == kValidateWithoutTol) {
473 else if (validate == kValidateWithout <<
474 { <<
475 tmpareacode[i] = GetAreaCode(tmpxx 443 tmpareacode[i] = GetAreaCode(tmpxx[i], false);
476 if (IsInside(tmpareacode[i])) << 444 if (IsInside(tmpareacode[i])) {
477 { <<
478 if (tmpdist[i] >= 0) tmpisvalid 445 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
479 continue; 446 continue;
480 } 447 }
481 } << 448 } else { // kDontValidate
482 else // kDontValidate <<
483 { <<
484 tmpareacode[i] = sInside; 449 tmpareacode[i] = sInside;
485 if (tmpdist[i] >= 0) tmpisvalid 450 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
486 continue; 451 continue;
487 } 452 }
488 } 453 }
489 454
490 if (tmpdist[0] <= tmpdist[1]) << 455 if (tmpdist[0] <= tmpdist[1]) {
491 { <<
492 distance[0] = tmpdist[0]; 456 distance[0] = tmpdist[0];
493 distance[1] = tmpdist[1]; 457 distance[1] = tmpdist[1];
494 xx[0] = tmpxx[0]; 458 xx[0] = tmpxx[0];
495 xx[1] = tmpxx[1]; 459 xx[1] = tmpxx[1];
496 gxx[0] = ComputeGlobalPoint(tmp 460 gxx[0] = ComputeGlobalPoint(tmpxx[0]);
497 gxx[1] = ComputeGlobalPoint(tmp 461 gxx[1] = ComputeGlobalPoint(tmpxx[1]);
498 areacode[0] = tmpareacode[0]; 462 areacode[0] = tmpareacode[0];
499 areacode[1] = tmpareacode[1]; 463 areacode[1] = tmpareacode[1];
500 isvalid[0] = tmpisvalid[0]; 464 isvalid[0] = tmpisvalid[0];
501 isvalid[1] = tmpisvalid[1]; 465 isvalid[1] = tmpisvalid[1];
502 } << 466 } else {
503 else <<
504 { <<
505 distance[0] = tmpdist[1]; 467 distance[0] = tmpdist[1];
506 distance[1] = tmpdist[0]; 468 distance[1] = tmpdist[0];
507 xx[0] = tmpxx[1]; 469 xx[0] = tmpxx[1];
508 xx[1] = tmpxx[0]; 470 xx[1] = tmpxx[0];
509 gxx[0] = ComputeGlobalPoint(tmp 471 gxx[0] = ComputeGlobalPoint(tmpxx[1]);
510 gxx[1] = ComputeGlobalPoint(tmp 472 gxx[1] = ComputeGlobalPoint(tmpxx[0]);
511 areacode[0] = tmpareacode[1]; 473 areacode[0] = tmpareacode[1];
512 areacode[1] = tmpareacode[0]; 474 areacode[1] = tmpareacode[0];
513 isvalid[0] = tmpisvalid[1]; 475 isvalid[0] = tmpisvalid[1];
514 isvalid[1] = tmpisvalid[0]; 476 isvalid[1] = tmpisvalid[0];
515 } 477 }
516 478
517 fCurStatWithV.SetCurrentStatus(0, gxx[0] 479 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
518 isvalid[0 480 isvalid[0], 2, validate, &gp, &gv);
519 fCurStatWithV.SetCurrentStatus(1, gxx[1] 481 fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
520 isvalid[1 482 isvalid[1], 2, validate, &gp, &gv);
521 vout = 2; 483 vout = 2;
522 } << 484
523 else << 485 } else {
524 { <<
525 // if D<0, no solution 486 // if D<0, no solution
526 // if D=0, just grazing the surfaces, re 487 // if D=0, just grazing the surfaces, return kInfinity
527 488
528 fCurStatWithV.SetCurrentStatus(0, gxx[0] 489 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
529 isvalid[0 490 isvalid[0], 0, validate, &gp, &gv);
530 vout = 0; 491 vout = 0;
531 } 492 }
532 return vout; 493 return vout;
533 } 494 }
534 495
>> 496
535 //============================================ 497 //=====================================================================
536 //* DistanceToSurface ------------------------ 498 //* DistanceToSurface -------------------------------------------------
537 499
538 G4int G4TwistTubsHypeSide::DistanceToSurface(c << 500 G4int G4TwistTubsHypeSide::DistanceToSurface(const G4ThreeVector &gp,
539 501 G4ThreeVector gxx[],
540 502 G4double distance[],
541 503 G4int areacode[])
542 { 504 {
543 // Find the approximate distance of a poin 505 // Find the approximate distance of a point of a hyperbolic surface.
544 // The distance must be an underestimate. 506 // The distance must be an underestimate.
545 // It will also be nice (although not nece 507 // It will also be nice (although not necessary) that the estimate is
546 // always finite no matter how close the p 508 // always finite no matter how close the point is.
547 // 509 //
548 // We arranged G4Hype::ApproxDistOutside a 510 // We arranged G4Hype::ApproxDistOutside and G4Hype::ApproxDistInside
549 // for this function. See these discriptio 511 // for this function. See these discriptions.
550 512
551 const G4double halftol 513 const G4double halftol
552 = 0.5 * G4GeometryTolerance::GetInstance( 514 = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
553 515
554 fCurStat.ResetfDone(kDontValidate, &gp); 516 fCurStat.ResetfDone(kDontValidate, &gp);
555 517
556 if (fCurStat.IsDone()) << 518 if (fCurStat.IsDone()) {
557 { << 519 for (G4int i=0; i<fCurStat.GetNXX(); i++) {
558 for (G4int i=0; i<fCurStat.GetNXX(); ++i <<
559 { <<
560 gxx[i] = fCurStat.GetXX(i); 520 gxx[i] = fCurStat.GetXX(i);
561 distance[i] = fCurStat.GetDistance(i) 521 distance[i] = fCurStat.GetDistance(i);
562 areacode[i] = fCurStat.GetAreacode(i) 522 areacode[i] = fCurStat.GetAreacode(i);
563 } 523 }
564 return fCurStat.GetNXX(); 524 return fCurStat.GetNXX();
565 } << 525 } else {
566 else // initialize << 526 // initialize
567 { << 527 for (G4int i=0; i<2; i++) {
568 for (auto i=0; i<2; ++i) <<
569 { <<
570 distance[i] = kInfinity; 528 distance[i] = kInfinity;
571 areacode[i] = sOutside; 529 areacode[i] = sOutside;
572 gxx[i].set(kInfinity, kInfinity, kInf 530 gxx[i].set(kInfinity, kInfinity, kInfinity);
573 } 531 }
574 } 532 }
575 533
576 534
577 G4ThreeVector p = ComputeLocalPoint(gp); 535 G4ThreeVector p = ComputeLocalPoint(gp);
578 G4ThreeVector xx; 536 G4ThreeVector xx;
579 537
580 // 538 //
581 // special case! 539 // special case!
582 // If p is on surface, return distance = 0 540 // If p is on surface, return distance = 0 immediatery .
583 // 541 //
584 G4ThreeVector lastgxx[2]; 542 G4ThreeVector lastgxx[2];
585 for (auto i=0; i<2; ++i) << 543 for (G4int i=0; i<2; i++) {
586 { <<
587 lastgxx[i] = fCurStatWithV.GetXX(i); 544 lastgxx[i] = fCurStatWithV.GetXX(i);
588 } 545 }
589 546
590 if ((gp - lastgxx[0]).mag() < halftol || (g << 547 if ((gp - lastgxx[0]).mag() < halftol || (gp - lastgxx[1]).mag() < halftol) {
591 { <<
592 // last winner, or last poststep point i 548 // last winner, or last poststep point is on the surface.
593 xx = p; 549 xx = p;
594 gxx[0] = gp; 550 gxx[0] = gp;
595 distance[0] = 0; 551 distance[0] = 0;
596 552
597 G4bool isvalid = true; 553 G4bool isvalid = true;
598 fCurStat.SetCurrentStatus(0, gxx[0], dis 554 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
599 isvalid, 1, kD 555 isvalid, 1, kDontValidate, &gp);
600 556
601 return 1; 557 return 1;
>> 558
602 } 559 }
603 // 560 //
604 // special case end 561 // special case end
605 // 562 //
606 563
607 G4double prho = p.getRho(); 564 G4double prho = p.getRho();
608 G4double pz = std::fabs(p.z()); 565 G4double pz = std::fabs(p.z()); // use symmetry
609 G4double r1 = std::sqrt(fR02 + pz * 566 G4double r1 = std::sqrt(fR02 + pz * pz * fTan2Stereo);
610 567
611 G4ThreeVector pabsz(p.x(), p.y(), pz); 568 G4ThreeVector pabsz(p.x(), p.y(), pz);
612 569
613 if (prho > r1 + halftol) // p is outside << 570 if (prho > r1 + halftol) { // p is outside of Hyperbolic surface
614 { << 571
615 // First point xx1 572 // First point xx1
616 G4double t = r1 / prho; 573 G4double t = r1 / prho;
617 G4ThreeVector xx1(t * pabsz.x(), t * pab 574 G4ThreeVector xx1(t * pabsz.x(), t * pabsz.y() , pz);
618 575
619 // Second point xx2 576 // Second point xx2
620 G4double z2 = (prho * fTanStereo + pz) / 577 G4double z2 = (prho * fTanStereo + pz) / (1 + fTan2Stereo);
621 G4double r2 = std::sqrt(fR02 + z2 * z2 * 578 G4double r2 = std::sqrt(fR02 + z2 * z2 * fTan2Stereo);
622 t = r2 / prho; 579 t = r2 / prho;
623 G4ThreeVector xx2(t * pabsz.x(), t * pab 580 G4ThreeVector xx2(t * pabsz.x(), t * pabsz.y() , z2);
624 581
625 G4double len = (xx2 - xx1).mag(); 582 G4double len = (xx2 - xx1).mag();
626 if (len < DBL_MIN) << 583 if (len < DBL_MIN) {
627 { <<
628 // xx2 = xx1?? I guess we 584 // xx2 = xx1?? I guess we
629 // must have really bracketed the nor 585 // must have really bracketed the normal
630 distance[0] = (pabsz - xx1).mag(); 586 distance[0] = (pabsz - xx1).mag();
631 xx = xx1; 587 xx = xx1;
632 } << 588 } else {
633 else <<
634 { <<
635 distance[0] = DistanceToLine(pabsz, x 589 distance[0] = DistanceToLine(pabsz, xx1, (xx2 - xx1) , xx);
636 } 590 }
637 591
638 } << 592 } else if (prho < r1 - halftol) { // p is inside of Hyperbolic surface.
639 else if (prho < r1 - halftol) // p is insi << 593
640 { <<
641 // First point xx1 594 // First point xx1
642 G4double t; 595 G4double t;
643 G4ThreeVector xx1; 596 G4ThreeVector xx1;
644 if (prho < DBL_MIN) << 597 if (prho < DBL_MIN) {
645 { <<
646 xx1.set(r1, 0. , pz); 598 xx1.set(r1, 0. , pz);
647 } << 599 } else {
648 else <<
649 { <<
650 t = r1 / prho; 600 t = r1 / prho;
651 xx1.set(t * pabsz.x(), t * pabsz.y() 601 xx1.set(t * pabsz.x(), t * pabsz.y() , pz);
652 } 602 }
653 603
654 // dr, dz is tangential vector of Hyparb 604 // dr, dz is tangential vector of Hyparbolic surface at xx1
655 // dr = r, dz = z*tan2stereo 605 // dr = r, dz = z*tan2stereo
656 G4double dr = pz * fTan2Stereo; 606 G4double dr = pz * fTan2Stereo;
657 G4double dz = r1; 607 G4double dz = r1;
658 G4double tanbeta = dr / dz; 608 G4double tanbeta = dr / dz;
659 G4double pztanbeta = pz * tanbeta; 609 G4double pztanbeta = pz * tanbeta;
660 610
661 // Second point xx2 611 // Second point xx2
662 // xx2 is intersection between x-axis an 612 // xx2 is intersection between x-axis and tangential vector
663 G4double r2 = r1 - pztanbeta; 613 G4double r2 = r1 - pztanbeta;
664 G4ThreeVector xx2; 614 G4ThreeVector xx2;
665 if (prho < DBL_MIN) << 615 if (prho < DBL_MIN) {
666 { <<
667 xx2.set(r2, 0. , 0.); 616 xx2.set(r2, 0. , 0.);
668 } << 617 } else {
669 else <<
670 { <<
671 t = r2 / prho; 618 t = r2 / prho;
672 xx2.set(t * pabsz.x(), t * pabsz.y() 619 xx2.set(t * pabsz.x(), t * pabsz.y() , 0.);
673 } 620 }
674 621
675 G4ThreeVector d = xx2 - xx1; 622 G4ThreeVector d = xx2 - xx1;
676 distance[0] = DistanceToLine(pabsz, xx1, 623 distance[0] = DistanceToLine(pabsz, xx1, d, xx);
677 624
678 } << 625 } else { // p is on Hyperbolic surface.
679 else // p is on Hyperbolic surface. << 626
680 { <<
681 distance[0] = 0; 627 distance[0] = 0;
682 xx.set(p.x(), p.y(), pz); 628 xx.set(p.x(), p.y(), pz);
>> 629
683 } 630 }
684 631
685 if (p.z() < 0) << 632 if (p.z() < 0) {
686 { <<
687 G4ThreeVector tmpxx(xx.x(), xx.y(), -xx. 633 G4ThreeVector tmpxx(xx.x(), xx.y(), -xx.z());
688 xx = tmpxx; 634 xx = tmpxx;
689 } 635 }
690 636
691 gxx[0] = ComputeGlobalPoint(xx); 637 gxx[0] = ComputeGlobalPoint(xx);
692 areacode[0] = sInside; 638 areacode[0] = sInside;
693 G4bool isvalid = true; 639 G4bool isvalid = true;
694 fCurStat.SetCurrentStatus(0, gxx[0], distan 640 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
695 isvalid, 1, kDont 641 isvalid, 1, kDontValidate, &gp);
696 return 1; 642 return 1;
697 } 643 }
698 644
699 //============================================ 645 //=====================================================================
700 //* GetAreaCode ------------------------------ 646 //* GetAreaCode -------------------------------------------------------
701 647
702 G4int G4TwistTubsHypeSide::GetAreaCode(const G << 648 G4int G4TwistTubsHypeSide::GetAreaCode(const G4ThreeVector &xx,
703 G 649 G4bool withTol)
704 { 650 {
705 const G4double ctol = 0.5 * kCarTolerance; 651 const G4double ctol = 0.5 * kCarTolerance;
706 G4int areacode = sInside; 652 G4int areacode = sInside;
707 653
708 if ((fAxis[0] == kPhi && fAxis[1] == kZAxis << 654 if ((fAxis[0] == kPhi && fAxis[1] == kZAxis)) {
709 { << 655 //G4int phiaxis = 0;
710 G4int zaxis = 1; 656 G4int zaxis = 1;
711 657
712 if (withTol) << 658 if (withTol) {
713 { << 659
714 G4bool isoutside = false; 660 G4bool isoutside = false;
715 G4int phiareacode = GetAreaCodeIn 661 G4int phiareacode = GetAreaCodeInPhi(xx);
716 G4bool isoutsideinphi = IsOutside(phi 662 G4bool isoutsideinphi = IsOutside(phiareacode);
717 663
718 // test boundary of phiaxis 664 // test boundary of phiaxis
719 665
720 if ((phiareacode & sAxisMin) == sAxis << 666 if ((phiareacode & sAxisMin) == sAxisMin) {
721 { << 667
722 areacode |= (sAxis0 & (sAxisPhi | 668 areacode |= (sAxis0 & (sAxisPhi | sAxisMin)) | sBoundary;
723 if (isoutsideinphi) isoutside = tr 669 if (isoutsideinphi) isoutside = true;
724 } << 670
725 else if ((phiareacode & sAxisMax) == << 671 } else if ((phiareacode & sAxisMax) == sAxisMax) {
726 { << 672
727 areacode |= (sAxis0 & (sAxisPhi | 673 areacode |= (sAxis0 & (sAxisPhi | sAxisMax)) | sBoundary;
728 if (isoutsideinphi) isoutside = tr 674 if (isoutsideinphi) isoutside = true;
>> 675
729 } 676 }
730 677
731 // test boundary of zaxis 678 // test boundary of zaxis
732 679
733 if (xx.z() < fAxisMin[zaxis] + ctol) << 680 if (xx.z() < fAxisMin[zaxis] + ctol) {
734 { << 681
735 areacode |= (sAxis1 & (sAxisZ | sA 682 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
736 if ((areacode & sBoundary) != 0) << 683 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
737 else areaco 684 else areacode |= sBoundary;
738 685
739 if (xx.z() <= fAxisMin[zaxis] - ct 686 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
740 687
741 } << 688 } else if (xx.z() > fAxisMax[zaxis] - ctol) {
742 else if (xx.z() > fAxisMax[zaxis] - c << 689
743 { <<
744 areacode |= (sAxis1 & (sAxisZ | sA 690 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
745 if ((areacode & sBoundary) != 0) << 691 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
746 else areaco 692 else areacode |= sBoundary;
747 693
748 if (xx.z() >= fAxisMax[zaxis] + ct 694 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
749 } 695 }
750 696
751 // if isoutside = true, clear sInside 697 // if isoutside = true, clear sInside bit.
752 // if not on boundary, add boundary i 698 // if not on boundary, add boundary information.
753 699
754 if (isoutside) << 700 if (isoutside) {
755 { <<
756 G4int tmpareacode = areacode & (~s 701 G4int tmpareacode = areacode & (~sInside);
757 areacode = tmpareacode; 702 areacode = tmpareacode;
758 } << 703 } else if ((areacode & sBoundary) != sBoundary) {
759 else if ((areacode & sBoundary) != sB <<
760 { <<
761 areacode |= (sAxis0 & sAxisPhi) | 704 areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
762 } 705 }
763 706
764 return areacode; 707 return areacode;
765 } << 708
766 else << 709 } else {
767 { << 710
768 G4int phiareacode = GetAreaCodeInPhi( 711 G4int phiareacode = GetAreaCodeInPhi(xx, false);
769 712
770 // test boundary of z-axis 713 // test boundary of z-axis
771 714
772 if (xx.z() < fAxisMin[zaxis]) << 715 if (xx.z() < fAxisMin[zaxis]) {
773 { << 716
774 areacode |= (sAxis1 & (sAxisZ | sA 717 areacode |= (sAxis1 & (sAxisZ | sAxisMin)) | sBoundary;
775 718
776 } << 719 } else if (xx.z() > fAxisMax[zaxis]) {
777 else if (xx.z() > fAxisMax[zaxis]) << 720
778 { <<
779 areacode |= (sAxis1 & (sAxisZ | sA 721 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) | sBoundary;
>> 722
780 } 723 }
781 724
782 // boundary of phi-axis 725 // boundary of phi-axis
783 726
784 if (phiareacode == sAxisMin) << 727 if (phiareacode == sAxisMin) {
785 { << 728
786 areacode |= (sAxis0 & (sAxisPhi | 729 areacode |= (sAxis0 & (sAxisPhi | sAxisMin));
787 if ((areacode & sBoundary) != 0) << 730 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
788 else areaco 731 else areacode |= sBoundary;
789 732
790 } << 733 } else if (phiareacode == sAxisMax) {
791 else if (phiareacode == sAxisMax) << 734
792 { <<
793 areacode |= (sAxis0 & (sAxisPhi | 735 areacode |= (sAxis0 & (sAxisPhi | sAxisMax));
794 if ((areacode & sBoundary) != 0) << 736 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
795 else areaco << 737 else areacode |= sBoundary;
>> 738
796 } 739 }
797 740
798 // if not on boundary, add boundary i 741 // if not on boundary, add boundary information.
799 742
800 if ((areacode & sBoundary) != sBounda << 743 if ((areacode & sBoundary) != sBoundary) {
801 { <<
802 areacode |= (sAxis0 & sAxisPhi) | 744 areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
803 } 745 }
804 return areacode; 746 return areacode;
805 } 747 }
806 } << 748 } else {
807 else <<
808 { <<
809 std::ostringstream message; 749 std::ostringstream message;
810 message << "Feature NOT implemented !" < 750 message << "Feature NOT implemented !" << G4endl
811 << " fAxis[0] = " << fAxi 751 << " fAxis[0] = " << fAxis[0] << G4endl
812 << " fAxis[1] = " << fAxi 752 << " fAxis[1] = " << fAxis[1];
813 G4Exception("G4TwistTubsHypeSide::GetAre 753 G4Exception("G4TwistTubsHypeSide::GetAreaCode()",
814 "GeomSolids0001", FatalExcep 754 "GeomSolids0001", FatalException, message);
815 } 755 }
816 return areacode; 756 return areacode;
817 } 757 }
818 758
819 //============================================ 759 //=====================================================================
820 //* GetAreaCodeInPhi ------------------------- 760 //* GetAreaCodeInPhi --------------------------------------------------
821 761
822 G4int G4TwistTubsHypeSide::GetAreaCodeInPhi(co << 762 G4int G4TwistTubsHypeSide::GetAreaCodeInPhi(const G4ThreeVector &xx,
823 763 G4bool withTol)
824 { 764 {
825 765
826 G4ThreeVector lowerlimit; // lower phi-boun 766 G4ThreeVector lowerlimit; // lower phi-boundary limit at z = xx.z()
827 G4ThreeVector upperlimit; // upper phi-boun 767 G4ThreeVector upperlimit; // upper phi-boundary limit at z = xx.z()
828 lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxis 768 lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxisMin, xx);
829 upperlimit = GetBoundaryAtPZ(sAxis0 & sAxis 769 upperlimit = GetBoundaryAtPZ(sAxis0 & sAxisMax, xx);
830 770
831 G4int areacode = sInside; 771 G4int areacode = sInside;
832 G4bool isoutside = false; 772 G4bool isoutside = false;
833 773
834 if (withTol) << 774 if (withTol) {
835 { << 775
836 if (AmIOnLeftSide(xx, lowerlimit) >= 0) << 776 if (AmIOnLeftSide(xx, lowerlimit) >= 0) { // xx is on lowerlimit
837 { <<
838 areacode |= (sAxisMin | sBoundary); 777 areacode |= (sAxisMin | sBoundary);
839 if (AmIOnLeftSide(xx, lowerlimit) > 0 778 if (AmIOnLeftSide(xx, lowerlimit) > 0) isoutside = true;
840 779
841 } << 780 } else if (AmIOnLeftSide(xx, upperlimit) <= 0) { // xx is on upperlimit
842 else if (AmIOnLeftSide(xx, upperlimit) < <<
843 { <<
844 areacode |= (sAxisMax | sBoundary); 781 areacode |= (sAxisMax | sBoundary);
845 if (AmIOnLeftSide(xx, upperlimit) < 0 782 if (AmIOnLeftSide(xx, upperlimit) < 0) isoutside = true;
846 } 783 }
847 784
848 // if isoutside = true, clear inside bit 785 // if isoutside = true, clear inside bit.
849 786
850 if (isoutside) << 787 if (isoutside) {
851 { <<
852 G4int tmpareacode = areacode & (~sIns 788 G4int tmpareacode = areacode & (~sInside);
853 areacode = tmpareacode; 789 areacode = tmpareacode;
854 } 790 }
855 } << 791
856 else << 792
857 { << 793 } else {
858 if (AmIOnLeftSide(xx, lowerlimit, false) << 794
859 { << 795 if (AmIOnLeftSide(xx, lowerlimit, false) >= 0) {
860 areacode |= (sAxisMin | sBoundary); 796 areacode |= (sAxisMin | sBoundary);
861 } << 797 } else if (AmIOnLeftSide(xx, upperlimit, false) <= 0) {
862 else if (AmIOnLeftSide(xx, upperlimit, f <<
863 { <<
864 areacode |= (sAxisMax | sBoundary); 798 areacode |= (sAxisMax | sBoundary);
865 } 799 }
866 } 800 }
867 801
868 return areacode; 802 return areacode;
>> 803
869 } 804 }
870 805
871 //============================================ 806 //=====================================================================
872 //* SetCorners(EndInnerRadius, EndOuterRadius, 807 //* SetCorners(EndInnerRadius, EndOuterRadius,DPhi,EndPhi,EndZ) -------
873 808
874 void G4TwistTubsHypeSide::SetCorners( G4double << 809 void G4TwistTubsHypeSide::SetCorners(
875 G4double << 810 G4double EndInnerRadius[2],
876 G4double << 811 G4double EndOuterRadius[2],
877 G4double << 812 G4double DPhi,
878 G4double << 813 G4double endPhi[2],
>> 814 G4double endZ[2]
>> 815 )
879 { 816 {
880 // Set Corner points in local coodinate. 817 // Set Corner points in local coodinate.
881 818
882 if (fAxis[0] == kPhi && fAxis[1] == kZAxis) 819 if (fAxis[0] == kPhi && fAxis[1] == kZAxis) {
883 820
>> 821 G4int i;
884 G4double endRad[2]; 822 G4double endRad[2];
885 G4double halfdphi = 0.5*DPhi; 823 G4double halfdphi = 0.5*DPhi;
886 824
887 for (auto i=0; i<2; ++i) // i=0,1 : -ve << 825 for (i=0; i<2; i++) { // i=0,1 : -ve z, +ve z
888 { << 826 endRad[i] = (fHandedness == 1 ? EndOuterRadius[i]
889 endRad[i] = (fHandedness == 1 ? EndOut << 827 : EndInnerRadius[i]);
890 } 828 }
891 829
892 G4int zmin = 0 ; // at -ve z 830 G4int zmin = 0 ; // at -ve z
893 G4int zmax = 1 ; // at +ve z 831 G4int zmax = 1 ; // at +ve z
894 832
895 G4double x, y, z; 833 G4double x, y, z;
896 834
897 // corner of Axis0min and Axis1min 835 // corner of Axis0min and Axis1min
898 x = endRad[zmin]*std::cos(endPhi[zmin] - 836 x = endRad[zmin]*std::cos(endPhi[zmin] - halfdphi);
899 y = endRad[zmin]*std::sin(endPhi[zmin] - 837 y = endRad[zmin]*std::sin(endPhi[zmin] - halfdphi);
900 z = endZ[zmin]; 838 z = endZ[zmin];
901 SetCorner(sC0Min1Min, x, y, z); 839 SetCorner(sC0Min1Min, x, y, z);
902 840
903 // corner of Axis0max and Axis1min 841 // corner of Axis0max and Axis1min
904 x = endRad[zmin]*std::cos(endPhi[zmin] + 842 x = endRad[zmin]*std::cos(endPhi[zmin] + halfdphi);
905 y = endRad[zmin]*std::sin(endPhi[zmin] + 843 y = endRad[zmin]*std::sin(endPhi[zmin] + halfdphi);
906 z = endZ[zmin]; 844 z = endZ[zmin];
907 SetCorner(sC0Max1Min, x, y, z); 845 SetCorner(sC0Max1Min, x, y, z);
908 846
909 // corner of Axis0max and Axis1max 847 // corner of Axis0max and Axis1max
910 x = endRad[zmax]*std::cos(endPhi[zmax] + 848 x = endRad[zmax]*std::cos(endPhi[zmax] + halfdphi);
911 y = endRad[zmax]*std::sin(endPhi[zmax] + 849 y = endRad[zmax]*std::sin(endPhi[zmax] + halfdphi);
912 z = endZ[zmax]; 850 z = endZ[zmax];
913 SetCorner(sC0Max1Max, x, y, z); 851 SetCorner(sC0Max1Max, x, y, z);
914 852
915 // corner of Axis0min and Axis1max 853 // corner of Axis0min and Axis1max
916 x = endRad[zmax]*std::cos(endPhi[zmax] - 854 x = endRad[zmax]*std::cos(endPhi[zmax] - halfdphi);
917 y = endRad[zmax]*std::sin(endPhi[zmax] - 855 y = endRad[zmax]*std::sin(endPhi[zmax] - halfdphi);
918 z = endZ[zmax]; 856 z = endZ[zmax];
919 SetCorner(sC0Min1Max, x, y, z); 857 SetCorner(sC0Min1Max, x, y, z);
920 858
921 } << 859 } else {
922 else <<
923 { <<
924 std::ostringstream message; 860 std::ostringstream message;
925 message << "Feature NOT implemented !" < 861 message << "Feature NOT implemented !" << G4endl
926 << " fAxis[0] = " << fAxi 862 << " fAxis[0] = " << fAxis[0] << G4endl
927 << " fAxis[1] = " << fAxi 863 << " fAxis[1] = " << fAxis[1];
928 G4Exception("G4TwistTubsHypeSide::SetCor 864 G4Exception("G4TwistTubsHypeSide::SetCorners()",
929 "GeomSolids0001", FatalExcep 865 "GeomSolids0001", FatalException, message);
930 } 866 }
931 } 867 }
932 868
>> 869
933 //============================================ 870 //=====================================================================
934 //* SetCorners() ----------------------------- 871 //* SetCorners() ------------------------------------------------------
935 872
936 void G4TwistTubsHypeSide::SetCorners() 873 void G4TwistTubsHypeSide::SetCorners()
937 { 874 {
938 G4Exception("G4TwistTubsHypeSide::SetCorner 875 G4Exception("G4TwistTubsHypeSide::SetCorners()",
939 "GeomSolids0001", FatalExceptio 876 "GeomSolids0001", FatalException,
940 "Method NOT implemented !"); 877 "Method NOT implemented !");
941 } 878 }
942 879
943 //============================================ 880 //=====================================================================
944 //* SetBoundaries() -------------------------- 881 //* SetBoundaries() ---------------------------------------------------
945 882
946 void G4TwistTubsHypeSide::SetBoundaries() 883 void G4TwistTubsHypeSide::SetBoundaries()
947 { 884 {
948 // Set direction-unit vector of phi-boundar 885 // Set direction-unit vector of phi-boundary-lines in local coodinate.
949 // sAxis0 must be kPhi. 886 // sAxis0 must be kPhi.
950 // This fanction set lower phi-boundary and 887 // This fanction set lower phi-boundary and upper phi-boundary.
951 888
952 if (fAxis[0] == kPhi && fAxis[1] == kZAxis) << 889 if (fAxis[0] == kPhi && fAxis[1] == kZAxis) {
953 { << 890
954 G4ThreeVector direction; 891 G4ThreeVector direction;
955 // sAxis0 & sAxisMin 892 // sAxis0 & sAxisMin
956 direction = GetCorner(sC0Min1Max) - GetC 893 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
957 direction = direction.unit(); 894 direction = direction.unit();
958 SetBoundary(sAxis0 & (sAxisPhi | sAxisMi 895 SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
959 GetCorner(sC0Min1Min), sAxi 896 GetCorner(sC0Min1Min), sAxisZ);
960 897
961 // sAxis0 & sAxisMax 898 // sAxis0 & sAxisMax
962 direction = GetCorner(sC0Max1Max) - GetC 899 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
963 direction = direction.unit(); 900 direction = direction.unit();
964 SetBoundary(sAxis0 & (sAxisPhi | sAxisMa 901 SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
965 GetCorner(sC0Max1Min), sAxis 902 GetCorner(sC0Max1Min), sAxisZ);
966 903
967 // sAxis1 & sAxisMin 904 // sAxis1 & sAxisMin
968 direction = GetCorner(sC0Max1Min) - GetC 905 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
969 direction = direction.unit(); 906 direction = direction.unit();
970 SetBoundary(sAxis1 & (sAxisZ | sAxisMin) 907 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
971 GetCorner(sC0Min1Min), sAxi 908 GetCorner(sC0Min1Min), sAxisPhi);
972 909
973 // sAxis1 & sAxisMax 910 // sAxis1 & sAxisMax
974 direction = GetCorner(sC0Max1Max) - GetC 911 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
975 direction = direction.unit(); 912 direction = direction.unit();
976 SetBoundary(sAxis1 & (sAxisZ | sAxisMax) 913 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
977 GetCorner(sC0Min1Max), sAxis 914 GetCorner(sC0Min1Max), sAxisPhi);
978 } << 915 } else {
979 else <<
980 { <<
981 std::ostringstream message; 916 std::ostringstream message;
982 message << "Feature NOT implemented !" < 917 message << "Feature NOT implemented !" << G4endl
983 << " fAxis[0] = " << fAxi 918 << " fAxis[0] = " << fAxis[0] << G4endl
984 << " fAxis[1] = " << fAxi 919 << " fAxis[1] = " << fAxis[1];
985 G4Exception("G4TwistTubsHypeSide::SetBou 920 G4Exception("G4TwistTubsHypeSide::SetBoundaries()",
986 "GeomSolids0001", FatalExcep 921 "GeomSolids0001", FatalException, message);
987 } 922 }
988 } 923 }
989 924
990 //============================================ 925 //=====================================================================
991 //* GetFacets() ------------------------------ 926 //* GetFacets() -------------------------------------------------------
992 927
993 void G4TwistTubsHypeSide::GetFacets( G4int k, 928 void G4TwistTubsHypeSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
994 G4int fac 929 G4int faces[][4], G4int iside )
995 { 930 {
>> 931
996 G4double z ; // the two parameters for t 932 G4double z ; // the two parameters for the surface equation
997 G4double x,xmin,xmax ; 933 G4double x,xmin,xmax ;
998 934
999 G4ThreeVector p ; // a point on the surface 935 G4ThreeVector p ; // a point on the surface, given by (z,u)
1000 936
1001 G4int nnode ; 937 G4int nnode ;
1002 G4int nface ; 938 G4int nface ;
1003 939
1004 // calculate the (n-1)*(k-1) vertices 940 // calculate the (n-1)*(k-1) vertices
1005 941
1006 for ( G4int i = 0 ; i<n ; ++i ) << 942 G4int i,j ;
1007 { << 943
>> 944 for ( i = 0 ; i<n ; i++ ) {
>> 945
1008 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin 946 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
1009 947
1010 for ( G4int j = 0 ; j<k ; ++j ) << 948 for ( j = 0 ; j<k ; j++ )
1011 { 949 {
1012 nnode = GetNode(i,j,k,n,iside) ; 950 nnode = GetNode(i,j,k,n,iside) ;
1013 951
1014 xmin = GetBoundaryMin(z) ; 952 xmin = GetBoundaryMin(z) ;
1015 xmax = GetBoundaryMax(z) ; 953 xmax = GetBoundaryMax(z) ;
1016 954
1017 if (fHandedness < 0) // inner hyperbol << 955 if (fHandedness < 0) { // inner hyperbolic surface
1018 { <<
1019 x = xmin + j*(xmax-xmin)/(k-1) ; 956 x = xmin + j*(xmax-xmin)/(k-1) ;
1020 } << 957 } else { // outer hyperbolic surface
1021 else // outer hyperbol <<
1022 { <<
1023 x = xmax - j*(xmax-xmin)/(k-1) ; 958 x = xmax - j*(xmax-xmin)/(k-1) ;
1024 } 959 }
1025 960
1026 p = SurfacePoint(x,z,true) ; // surfac 961 p = SurfacePoint(x,z,true) ; // surface point in global coord.system
1027 962
1028 xyz[nnode][0] = p.x() ; 963 xyz[nnode][0] = p.x() ;
1029 xyz[nnode][1] = p.y() ; 964 xyz[nnode][1] = p.y() ;
1030 xyz[nnode][2] = p.z() ; 965 xyz[nnode][2] = p.z() ;
1031 966
1032 if ( i<n-1 && j<k-1 ) // clock wise << 967 if ( i<n-1 && j<k-1 ) { // clock wise filling
1033 { << 968
1034 nface = GetFace(i,j,k,n,iside) ; 969 nface = GetFace(i,j,k,n,iside) ;
1035 faces[nface][0] = GetEdgeVisibility(i,j,k,n << 970
1036 * ( GetNode(i ,j ,k << 971 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i ,j ,k,n,iside)+1) ;
1037 faces[nface][1] = GetEdgeVisibility(i,j,k,n << 972 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
1038 * ( GetNode(i+1,j ,k << 973 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1039 faces[nface][2] = GetEdgeVisibility(i,j,k,n << 974 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
1040 * ( GetNode(i+1,j+1,k << 975
1041 faces[nface][3] = GetEdgeVisibility(i,j,k,n <<
1042 * ( GetNode(i ,j+1,k <<
1043 } 976 }
1044 } 977 }
1045 } 978 }
1046 } 979 }
1047 980