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