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