| 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 // G4TwistTrapAlphaSide implementation 26 // G4TwistTrapAlphaSide implementation
27 // 27 //
28 // Author: 18/03/2005 - O.Link (Oliver.Link@ce 28 // Author: 18/03/2005 - O.Link (Oliver.Link@cern.ch)
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include <cmath> 31 #include <cmath>
32 32
33 #include "G4TwistTrapAlphaSide.hh" 33 #include "G4TwistTrapAlphaSide.hh"
34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh"
35 #include "G4JTPolynomialSolver.hh" 35 #include "G4JTPolynomialSolver.hh"
36 36
37 //============================================ 37 //=====================================================================
38 //* constructors ----------------------------- 38 //* constructors ------------------------------------------------------
39 39
40 G4TwistTrapAlphaSide:: 40 G4TwistTrapAlphaSide::
41 G4TwistTrapAlphaSide(const G4String& name, 41 G4TwistTrapAlphaSide(const G4String& name,
42 G4double PhiTwist, 42 G4double PhiTwist, // twist angle
43 G4double pDz, 43 G4double pDz, // half z lenght
44 G4double pTheta, 44 G4double pTheta, // direction between end planes
45 G4double pPhi, 45 G4double pPhi, // by polar and azimutal angles
46 G4double pDy1, 46 G4double pDy1, // half y length at -pDz
47 G4double pDx1, 47 G4double pDx1, // half x length at -pDz,-pDy
48 G4double pDx2, 48 G4double pDx2, // half x length at -pDz,+pDy
49 G4double pDy2, 49 G4double pDy2, // half y length at +pDz
50 G4double pDx3, 50 G4double pDx3, // half x length at +pDz,-pDy
51 G4double pDx4, 51 G4double pDx4, // half x length at +pDz,+pDy
52 G4double pAlph, 52 G4double pAlph, // tilt angle at +pDz
53 G4double AngleSide 53 G4double AngleSide // parity
54 ) 54 )
55 : G4VTwistSurface(name) 55 : G4VTwistSurface(name)
56 { 56 {
57 fAxis[0] = kYAxis; // in local coordinate 57 fAxis[0] = kYAxis; // in local coordinate system
58 fAxis[1] = kZAxis; 58 fAxis[1] = kZAxis;
59 fAxisMin[0] = -kInfinity ; // Y Axis bounda 59 fAxisMin[0] = -kInfinity ; // Y Axis boundary
60 fAxisMax[0] = kInfinity ; // depends on 60 fAxisMax[0] = kInfinity ; // depends on z !!
61 fAxisMin[1] = -pDz ; // Z Axis boundary 61 fAxisMin[1] = -pDz ; // Z Axis boundary
62 fAxisMax[1] = pDz ; 62 fAxisMax[1] = pDz ;
63 63
64 fDx1 = pDx1 ; 64 fDx1 = pDx1 ;
65 fDx2 = pDx2 ; 65 fDx2 = pDx2 ;
66 fDx3 = pDx3 ; 66 fDx3 = pDx3 ;
67 fDx4 = pDx4 ; 67 fDx4 = pDx4 ;
68 68
69 fDy1 = pDy1 ; 69 fDy1 = pDy1 ;
70 fDy2 = pDy2 ; 70 fDy2 = pDy2 ;
71 71
72 fDz = pDz ; 72 fDz = pDz ;
73 73
74 fAlph = pAlph ; 74 fAlph = pAlph ;
75 fTAlph = std::tan(fAlph) ; 75 fTAlph = std::tan(fAlph) ;
76 76
77 fTheta = pTheta ; 77 fTheta = pTheta ;
78 fPhi = pPhi ; 78 fPhi = pPhi ;
79 79
80 // precalculate frequently used parameters 80 // precalculate frequently used parameters
81 fDx4plus2 = fDx4 + fDx2 ; 81 fDx4plus2 = fDx4 + fDx2 ;
82 fDx4minus2 = fDx4 - fDx2 ; 82 fDx4minus2 = fDx4 - fDx2 ;
83 fDx3plus1 = fDx3 + fDx1 ; 83 fDx3plus1 = fDx3 + fDx1 ;
84 fDx3minus1 = fDx3 - fDx1 ; 84 fDx3minus1 = fDx3 - fDx1 ;
85 fDy2plus1 = fDy2 + fDy1 ; 85 fDy2plus1 = fDy2 + fDy1 ;
86 fDy2minus1 = fDy2 - fDy1 ; 86 fDy2minus1 = fDy2 - fDy1 ;
87 87
88 fa1md1 = 2*fDx2 - 2*fDx1 ; 88 fa1md1 = 2*fDx2 - 2*fDx1 ;
89 fa2md2 = 2*fDx4 - 2*fDx3 ; 89 fa2md2 = 2*fDx4 - 2*fDx3 ;
90 90
91 fPhiTwist = PhiTwist ; // dphi << 91 fPhiTwist = PhiTwist ; // dphi
92 fAngleSide = AngleSide ; // 0,90,180,270 de 92 fAngleSide = AngleSide ; // 0,90,180,270 deg
93 93
94 fdeltaX = 2 * fDz * std::tan(fTheta) * std:: 94 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi);
95 // dx in surface equation 95 // dx in surface equation
96 fdeltaY = 2 * fDz * std::tan(fTheta) * std:: 96 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi);
97 // dy in surface equation 97 // dy in surface equation
98 98
99 fRot.rotateZ( AngleSide ) ; 99 fRot.rotateZ( AngleSide ) ;
100 100
101 fTrans.set(0, 0, 0); // No Translation 101 fTrans.set(0, 0, 0); // No Translation
102 fIsValidNorm = false; 102 fIsValidNorm = false;
103 103
104 SetCorners() ; 104 SetCorners() ;
105 SetBoundaries() ; 105 SetBoundaries() ;
106 } 106 }
107 107
108 108
109 //============================================ 109 //=====================================================================
110 //* Fake default constructor ----------------- 110 //* Fake default constructor ------------------------------------------
111 111
112 G4TwistTrapAlphaSide::G4TwistTrapAlphaSide( __ 112 G4TwistTrapAlphaSide::G4TwistTrapAlphaSide( __void__& a )
113 : G4VTwistSurface(a) << 113 : G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
>> 114 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
>> 115 fAngleSide(0.), fDx4plus2(0.), fDx4minus2(0.), fDx3plus1(0.), fDx3minus1(0.),
>> 116 fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.), fa2md2(0.), fdeltaX(0.),
>> 117 fdeltaY(0.)
114 { 118 {
115 } 119 }
116 120
117 121
118 //============================================ 122 //=====================================================================
119 //* destructor ------------------------------- 123 //* destructor --------------------------------------------------------
120 124
121 G4TwistTrapAlphaSide::~G4TwistTrapAlphaSide() << 125 G4TwistTrapAlphaSide::~G4TwistTrapAlphaSide()
>> 126 {
>> 127 }
122 128
123 129
124 //============================================ 130 //=====================================================================
125 //* GetNormal -------------------------------- 131 //* GetNormal ---------------------------------------------------------
126 132
127 G4ThreeVector 133 G4ThreeVector
128 G4TwistTrapAlphaSide::GetNormal(const G4ThreeV 134 G4TwistTrapAlphaSide::GetNormal(const G4ThreeVector& tmpxx,
129 G4bool i 135 G4bool isGlobal)
130 { 136 {
131 // GetNormal returns a normal vector at a s 137 // GetNormal returns a normal vector at a surface (or very close
132 // to surface) point at tmpxx. 138 // to surface) point at tmpxx.
133 // If isGlobal=true, it returns the normal 139 // If isGlobal=true, it returns the normal in global coordinate.
134 // 140 //
135 141
136 G4ThreeVector xx; 142 G4ThreeVector xx;
137 if (isGlobal) 143 if (isGlobal)
138 { 144 {
139 xx = ComputeLocalPoint(tmpxx); 145 xx = ComputeLocalPoint(tmpxx);
140 if ((xx - fCurrentNormal.p).mag() < 0.5 146 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
141 { 147 {
142 return ComputeGlobalDirection(fCurren 148 return ComputeGlobalDirection(fCurrentNormal.normal);
143 } 149 }
144 } 150 }
145 else 151 else
146 { 152 {
147 xx = tmpxx; 153 xx = tmpxx;
148 if (xx == fCurrentNormal.p) 154 if (xx == fCurrentNormal.p)
149 { 155 {
150 return fCurrentNormal.normal; 156 return fCurrentNormal.normal;
151 } 157 }
152 } 158 }
153 159
154 G4double phi ; 160 G4double phi ;
155 G4double u ; 161 G4double u ;
156 162
157 GetPhiUAtX(xx,phi,u) ; // phi,u for point 163 GetPhiUAtX(xx,phi,u) ; // phi,u for point xx close to surface
158 164
159 G4ThreeVector normal = NormAng(phi,u) ; / 165 G4ThreeVector normal = NormAng(phi,u) ; // the normal vector at phi,u
160 166
161 #ifdef G4TWISTDEBUG 167 #ifdef G4TWISTDEBUG
162 G4cout << "normal vector = " << normal << 168 G4cout << "normal vector = " << normal << G4endl ;
163 G4cout << "phi = " << phi << " , u = " << u 169 G4cout << "phi = " << phi << " , u = " << u << G4endl ;
164 #endif 170 #endif
165 171
166 if (isGlobal) 172 if (isGlobal)
167 { 173 {
168 fCurrentNormal.normal = ComputeGlobalDir 174 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
169 } 175 }
170 else 176 else
171 { 177 {
172 fCurrentNormal.normal = normal.unit(); 178 fCurrentNormal.normal = normal.unit();
173 } 179 }
174 180
175 return fCurrentNormal.normal; 181 return fCurrentNormal.normal;
176 } 182 }
177 183
178 //============================================ 184 //=====================================================================
179 //* DistanceToSurface ------------------------ 185 //* DistanceToSurface -------------------------------------------------
180 186
181 G4int 187 G4int
182 G4TwistTrapAlphaSide::DistanceToSurface(const 188 G4TwistTrapAlphaSide::DistanceToSurface(const G4ThreeVector& gp,
183 const 189 const G4ThreeVector& gv,
184 190 G4ThreeVector gxx[],
185 191 G4double distance[],
186 192 G4int areacode[],
187 193 G4bool isvalid[],
188 194 EValidate validate)
189 { 195 {
190 static const G4double pihalf = pi/2 ; 196 static const G4double pihalf = pi/2 ;
191 const G4double ctol = 0.5 * kCarTolerance; 197 const G4double ctol = 0.5 * kCarTolerance;
192 198
193 G4bool IsParallel = false ; 199 G4bool IsParallel = false ;
194 G4bool IsConverged = false ; 200 G4bool IsConverged = false ;
195 201
196 G4int nxx = 0 ; // number of physical solut 202 G4int nxx = 0 ; // number of physical solutions
197 203
198 fCurStatWithV.ResetfDone(validate, &gp, &gv) 204 fCurStatWithV.ResetfDone(validate, &gp, &gv);
199 205
200 if (fCurStatWithV.IsDone()) 206 if (fCurStatWithV.IsDone())
201 { 207 {
202 for (G4int i=0; i<fCurStatWithV.GetNXX(); 208 for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
203 { 209 {
204 gxx[i] = fCurStatWithV.GetXX(i); 210 gxx[i] = fCurStatWithV.GetXX(i);
205 distance[i] = fCurStatWithV.GetDistance( 211 distance[i] = fCurStatWithV.GetDistance(i);
206 areacode[i] = fCurStatWithV.GetAreacode( 212 areacode[i] = fCurStatWithV.GetAreacode(i);
207 isvalid[i] = fCurStatWithV.IsValid(i); 213 isvalid[i] = fCurStatWithV.IsValid(i);
208 } 214 }
209 return fCurStatWithV.GetNXX(); 215 return fCurStatWithV.GetNXX();
210 } 216 }
211 else // initialise 217 else // initialise
212 { 218 {
213 for (G4int j=0; j<G4VSURFACENXX ; ++j) 219 for (G4int j=0; j<G4VSURFACENXX ; ++j)
214 { 220 {
215 distance[j] = kInfinity; 221 distance[j] = kInfinity;
216 areacode[j] = sOutside; 222 areacode[j] = sOutside;
217 isvalid[j] = false; 223 isvalid[j] = false;
218 gxx[j].set(kInfinity, kInfinity, kInfini 224 gxx[j].set(kInfinity, kInfinity, kInfinity);
219 } 225 }
220 } 226 }
221 227
222 G4ThreeVector p = ComputeLocalPoint(gp); 228 G4ThreeVector p = ComputeLocalPoint(gp);
223 G4ThreeVector v = ComputeLocalDirection(gv); 229 G4ThreeVector v = ComputeLocalDirection(gv);
224 230
225 #ifdef G4TWISTDEBUG 231 #ifdef G4TWISTDEBUG
226 G4cout << "Local point p = " << p << G4endl 232 G4cout << "Local point p = " << p << G4endl ;
227 G4cout << "Local direction v = " << v << G4e 233 G4cout << "Local direction v = " << v << G4endl ;
228 #endif 234 #endif
229 235
230 G4double phi,u ; // parameters 236 G4double phi,u ; // parameters
231 237
232 // temporary variables 238 // temporary variables
233 239
234 G4double tmpdist = kInfinity ; 240 G4double tmpdist = kInfinity ;
235 G4ThreeVector tmpxx; 241 G4ThreeVector tmpxx;
236 G4int tmpareacode = sOutside ; 242 G4int tmpareacode = sOutside ;
237 G4bool tmpisvalid = false ; 243 G4bool tmpisvalid = false ;
238 244
239 std::vector<Intersection> xbuf ; 245 std::vector<Intersection> xbuf ;
240 Intersection xbuftmp ; 246 Intersection xbuftmp ;
241 247
242 // prepare some variables for the intersecti 248 // prepare some variables for the intersection finder
243 249
244 G4double L = 2*fDz ; 250 G4double L = 2*fDz ;
245 251
246 G4double phixz = fPhiTwist * ( p.x() * v.z() 252 G4double phixz = fPhiTwist * ( p.x() * v.z() - p.z() * v.x() ) ;
247 G4double phiyz = fPhiTwist * ( p.y() * v.z() 253 G4double phiyz = fPhiTwist * ( p.y() * v.z() - p.z() * v.y() ) ;
248 254
249 255
250 // special case vz = 0 256 // special case vz = 0
251 257
252 if ( v.z() == 0. ) 258 if ( v.z() == 0. )
253 { << 259 {
254 if ( std::fabs(p.z()) <= L ) // intersec 260 if ( std::fabs(p.z()) <= L ) // intersection possible in z
255 { 261 {
256 phi = p.z() * fPhiTwist / L ; // phi is 262 phi = p.z() * fPhiTwist / L ; // phi is determined by the z-position
257 u = (fDy1*(4*(-(fdeltaY*phi*v.x()) + fPh 263 u = (fDy1*(4*(-(fdeltaY*phi*v.x()) + fPhiTwist*p.y()*v.x()
258 + fdeltaX*phi*v.y() - fPhi 264 + fdeltaX*phi*v.y() - fPhiTwist*p.x()*v.y())
259 + ((fDx3plus1 + fDx4plus2)*fP 265 + ((fDx3plus1 + fDx4plus2)*fPhiTwist
260 + 2*(fDx3minus1 + fDx4minus2 266 + 2*(fDx3minus1 + fDx4minus2)*phi)
261 *(v.y()*std::cos(phi) - v 267 *(v.y()*std::cos(phi) - v.x()*std::sin(phi))))
262 /(fPhiTwist*(4*fDy1* v.x() - (fa1md1 268 /(fPhiTwist*(4*fDy1* v.x() - (fa1md1 + 4*fDy1*fTAlph)*v.y())
263 *std::cos(phi) + fPhiTwist 269 *std::cos(phi) + fPhiTwist*(fa1md1*v.x()
264 + 4*fDy1*(fTAlph*v.x() 270 + 4*fDy1*(fTAlph*v.x() + v.y()))*std::sin(phi));
265 xbuftmp.phi = phi ; 271 xbuftmp.phi = phi ;
266 xbuftmp.u = u ; 272 xbuftmp.u = u ;
267 xbuftmp.areacode = sOutside ; 273 xbuftmp.areacode = sOutside ;
268 xbuftmp.distance = kInfinity ; 274 xbuftmp.distance = kInfinity ;
269 xbuftmp.isvalid = false ; 275 xbuftmp.isvalid = false ;
270 276
271 xbuf.push_back(xbuftmp) ; // store it t 277 xbuf.push_back(xbuftmp) ; // store it to xbuf
272 } 278 }
273 else // no intersection possible 279 else // no intersection possible
274 { 280 {
275 distance[0] = kInfinity; 281 distance[0] = kInfinity;
276 gxx[0].set(kInfinity,kInfinity,kInfinity 282 gxx[0].set(kInfinity,kInfinity,kInfinity);
277 isvalid[0] = false ; 283 isvalid[0] = false ;
278 areacode[0] = sOutside ; 284 areacode[0] = sOutside ;
279 fCurStatWithV.SetCurrentStatus(0, gxx[0] 285 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
280 areacode[ 286 areacode[0], isvalid[0],
281 0, valida 287 0, validate, &gp, &gv);
282 return 0; 288 return 0;
283 } // end std::fabs(p.z() <= L 289 } // end std::fabs(p.z() <= L
284 } // end v.z() == 0 290 } // end v.z() == 0
285 else // general solution for non-zero vz 291 else // general solution for non-zero vz
286 { 292 {
287 293
288 G4double c[8],srd[7],si[7] ; 294 G4double c[8],srd[7],si[7] ;
289 295
290 c[7] = 57600* 296 c[7] = 57600*
291 fDy1*(fa1md1*phiyz + 297 fDy1*(fa1md1*phiyz +
292 fDy1*(-4*phixz + 4*fTAlph*phiyz 298 fDy1*(-4*phixz + 4*fTAlph*phiyz
293 + (fDx3plus1 + fDx4plus2)*fPh 299 + (fDx3plus1 + fDx4plus2)*fPhiTwist*v.z())) ;
294 c[6] = -57600* 300 c[6] = -57600*
295 fDy1*(4*fDy1*(phiyz + 2*fDz*v.x() + fTAl 301 fDy1*(4*fDy1*(phiyz + 2*fDz*v.x() + fTAlph*(phixz - 2*fDz*v.y()))
296 - 2*fDy1*(2*fdeltaX + fDx3minus1 + f 302 - 2*fDy1*(2*fdeltaX + fDx3minus1 + fDx4minus2
297 - 2*fdeltaY*fTAlph)*v.z() 303 - 2*fdeltaY*fTAlph)*v.z()
298 + fa1md1*(phixz - 2*fDz*v.y() + fdel 304 + fa1md1*(phixz - 2*fDz*v.y() + fdeltaY*v.z()));
299 c[5] = 4800* 305 c[5] = 4800*
300 fDy1*(fa1md1*(-5*phiyz - 24*fDz*v.x() + 306 fDy1*(fa1md1*(-5*phiyz - 24*fDz*v.x() + 12*fdeltaX*v.z()) +
301 fDy1*(20*phixz - 4*(5*fTAlph*phiyz 307 fDy1*(20*phixz - 4*(5*fTAlph*phiyz + 24*fDz*fTAlph*v.x()
302 + 24*fDz*v.y()) + (48*fdelt 308 + 24*fDz*v.y()) + (48*fdeltaY + (fDx3plus1 + fDx4plus2)
303 *fPhiTwist 309 *fPhiTwist + 48*fdeltaX*fTAlph)*v.z()));
304 c[4] = 4800* 310 c[4] = 4800*
305 fDy1*(fa1md1*(phixz - 10*fDz*v.y() + 5*f 311 fDy1*(fa1md1*(phixz - 10*fDz*v.y() + 5*fdeltaY*v.z())
306 + 2*fDy1*(2*phiyz + 20*fDz*v.x() 312 + 2*fDy1*(2*phiyz + 20*fDz*v.x()
307 + (-10*fdeltaX + fDx3minus1 313 + (-10*fdeltaX + fDx3minus1 + fDx4minus2)*v.z()
308 + 2*fTAlph*(phixz - 10*fDz* 314 + 2*fTAlph*(phixz - 10*fDz*v.y() + 5*fdeltaY*v.z())));
309 c[3] = -96* 315 c[3] = -96*
310 fDy1*(-(fa1md1*(phiyz + 100*fDz*v.x() - 316 fDy1*(-(fa1md1*(phiyz + 100*fDz*v.x() - 50*fdeltaX*v.z()))
311 + fDy1*(4*phixz - 400*fDz*v.y() 317 + fDy1*(4*phixz - 400*fDz*v.y()
312 + (200*fdeltaY - (fDx3plus1 318 + (200*fdeltaY - (fDx3plus1 + fDx4plus2)*fPhiTwist)*v.z()
313 - 4*fTAlph*(phiyz + 100*fDz* 319 - 4*fTAlph*(phiyz + 100*fDz*v.x() - 50*fdeltaX*v.z())));
314 c[2] = 32* 320 c[2] = 32*
315 fDy1*(4*fDy1*(7*fTAlph*phixz + 7*phiyz - 321 fDy1*(4*fDy1*(7*fTAlph*phixz + 7*phiyz - 6*fDz*v.x() + 6*fDz*fTAlph*v.y())
316 + 6*fDy1*(2*fdeltaX+fDx3minus1+fDx4m 322 + 6*fDy1*(2*fdeltaX+fDx3minus1+fDx4minus2-2*fdeltaY*fTAlph)*v.z()
317 + fa1md1*(7*phixz + 6*fDz*v.y() - 3* 323 + fa1md1*(7*phixz + 6*fDz*v.y() - 3*fdeltaY*v.z()));
318 c[1] = -8* 324 c[1] = -8*
319 fDy1*(fa1md1*(-9*phiyz - 56*fDz*v.x() + 325 fDy1*(fa1md1*(-9*phiyz - 56*fDz*v.x() + 28*fdeltaX*v.z())
320 + 4*fDy1*(9*phixz - 9*fTAlph*phiyz - 326 + 4*fDy1*(9*phixz - 9*fTAlph*phiyz - 56*fDz*fTAlph*v.x()
321 - 56*fDz*v.y() + 28*(fdeltaY 327 - 56*fDz*v.y() + 28*(fdeltaY + fdeltaX*fTAlph)*v.z()));
322 c[0] = 72* 328 c[0] = 72*
323 fDy1*(fa1md1*(2*fDz*v.y() - fdeltaY*v.z( 329 fDy1*(fa1md1*(2*fDz*v.y() - fdeltaY*v.z())
324 + fDy1*(-8*fDz*v.x() + 8*fDz*fTAlph* 330 + fDy1*(-8*fDz*v.x() + 8*fDz*fTAlph*v.y()
325 + 4*fdeltaX*v.z() - 4*fdeltaY* 331 + 4*fdeltaX*v.z() - 4*fdeltaY*fTAlph*v.z()));
326 332
327 #ifdef G4TWISTDEBUG 333 #ifdef G4TWISTDEBUG
328 G4cout << "coef = " << c[0] << " " 334 G4cout << "coef = " << c[0] << " "
329 << c[1] << " " 335 << c[1] << " "
330 << c[2] << " " 336 << c[2] << " "
331 << c[3] << " " 337 << c[3] << " "
332 << c[4] << " " 338 << c[4] << " "
333 << c[5] << " " 339 << c[5] << " "
334 << c[6] << " " 340 << c[6] << " "
335 << c[7] << G4endl ; 341 << c[7] << G4endl ;
336 #endif 342 #endif
337 343
338 G4JTPolynomialSolver trapEq ; 344 G4JTPolynomialSolver trapEq ;
339 G4int num = trapEq.FindRoots(c,7,srd,si); 345 G4int num = trapEq.FindRoots(c,7,srd,si);
340 346
341 for (G4int i = 0 ; i<num ; i++ ) // loop 347 for (G4int i = 0 ; i<num ; i++ ) // loop over all math solutions
342 { 348 {
343 if ( si[i]==0.0 ) // only real solution 349 if ( si[i]==0.0 ) // only real solutions
344 { 350 {
345 #ifdef G4TWISTDEBUG 351 #ifdef G4TWISTDEBUG
346 G4cout << "Solution " << i << " : " << 352 G4cout << "Solution " << i << " : " << srd[i] << G4endl ;
347 #endif 353 #endif
348 phi = std::fmod(srd[i] , pihalf) ; 354 phi = std::fmod(srd[i] , pihalf) ;
349 u = (fDy1*(4*(phiyz + 2*fDz*phi*v.y( 355 u = (fDy1*(4*(phiyz + 2*fDz*phi*v.y() - fdeltaY*phi*v.z())
350 - ((fDx3plus1 + fDx4plus2 356 - ((fDx3plus1 + fDx4plus2)*fPhiTwist
351 + 2*(fDx3minus1 + fDx4m 357 + 2*(fDx3minus1 + fDx4minus2)*phi)*v.z()*std::sin(phi)))
352 /(fPhiTwist*v.z()*(4*fDy1*std::co 358 /(fPhiTwist*v.z()*(4*fDy1*std::cos(phi)
353 + (fa1md1 + 4*fDy 359 + (fa1md1 + 4*fDy1*fTAlph)*std::sin(phi)));
354 xbuftmp.phi = phi ; 360 xbuftmp.phi = phi ;
355 xbuftmp.u = u ; 361 xbuftmp.u = u ;
356 xbuftmp.areacode = sOutside ; 362 xbuftmp.areacode = sOutside ;
357 xbuftmp.distance = kInfinity ; 363 xbuftmp.distance = kInfinity ;
358 xbuftmp.isvalid = false ; 364 xbuftmp.isvalid = false ;
359 365
360 xbuf.push_back(xbuftmp) ; // store it 366 xbuf.push_back(xbuftmp) ; // store it to xbuf
361 367
362 #ifdef G4TWISTDEBUG 368 #ifdef G4TWISTDEBUG
363 G4cout << "solution " << i << " = " << 369 G4cout << "solution " << i << " = " << phi << " , " << u << G4endl ;
364 #endif 370 #endif
365 } // end if real solution 371 } // end if real solution
366 } // end loop i 372 } // end loop i
367 } // end general case 373 } // end general case
368 374
369 nxx = (G4int)xbuf.size() ; // save the numb << 375 nxx = xbuf.size() ; // save the number of solutions
370 376
371 G4ThreeVector xxonsurface ; // point on 377 G4ThreeVector xxonsurface ; // point on surface
372 G4ThreeVector surfacenormal ; // normal ve 378 G4ThreeVector surfacenormal ; // normal vector
373 G4double deltaX; // distance between inters 379 G4double deltaX; // distance between intersection point and point on surface
374 G4double theta; // angle between track and 380 G4double theta; // angle between track and surfacenormal
375 G4double factor; // a scaling factor 381 G4double factor; // a scaling factor
376 G4int maxint=30; // number of iterations 382 G4int maxint=30; // number of iterations
377 383
378 for (auto & k : xbuf) << 384 for ( size_t k = 0 ; k<xbuf.size() ; ++k )
379 { 385 {
380 #ifdef G4TWISTDEBUG 386 #ifdef G4TWISTDEBUG
381 G4cout << "Solution " << k << " : " 387 G4cout << "Solution " << k << " : "
382 << "reconstructed phiR = " << xbuf[ 388 << "reconstructed phiR = " << xbuf[k].phi
383 << ", uR = " << xbuf[k].u << G4endl 389 << ", uR = " << xbuf[k].u << G4endl ;
384 #endif 390 #endif
385 391
386 phi = k.phi ; // get the stored values fo << 392 phi = xbuf[k].phi ; // get the stored values for phi and u
387 u = k.u ; << 393 u = xbuf[k].u ;
388 394
389 IsConverged = false ; // no convergence 395 IsConverged = false ; // no convergence at the beginning
390 396
391 for ( G4int i = 1 ; i<maxint ; ++i ) 397 for ( G4int i = 1 ; i<maxint ; ++i )
392 { 398 {
393 xxonsurface = SurfacePoint(phi,u) ; 399 xxonsurface = SurfacePoint(phi,u) ;
394 surfacenormal = NormAng(phi,u) ; 400 surfacenormal = NormAng(phi,u) ;
395 401
396 tmpdist = DistanceToPlaneWithV(p, v, xxo 402 tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx);
397 deltaX = ( tmpxx - xxonsurface ).mag() ; 403 deltaX = ( tmpxx - xxonsurface ).mag() ;
398 theta = std::fabs(std::acos(v*surfacenor 404 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
399 if ( theta < 0.001 ) 405 if ( theta < 0.001 )
400 { 406 {
401 factor = 50 ; 407 factor = 50 ;
402 IsParallel = true ; 408 IsParallel = true ;
403 } 409 }
404 else 410 else
405 { 411 {
406 factor = 1 ; 412 factor = 1 ;
407 } 413 }
408 414
409 #ifdef G4TWISTDEBUG 415 #ifdef G4TWISTDEBUG
410 G4cout << "Step i = " << i << ", distanc 416 G4cout << "Step i = " << i << ", distance = " << tmpdist
411 << ", " << deltaX << G4endl ; 417 << ", " << deltaX << G4endl ;
412 G4cout << "X = " << tmpxx << G4endl ; 418 G4cout << "X = " << tmpxx << G4endl ;
413 #endif 419 #endif
414 420
415 GetPhiUAtX(tmpxx, phi, u) ; 421 GetPhiUAtX(tmpxx, phi, u) ;
416 // the new point xx is accepted and ph 422 // the new point xx is accepted and phi/u replaced
417 423
418 #ifdef G4TWISTDEBUG 424 #ifdef G4TWISTDEBUG
419 G4cout << "approximated phi = " << phi < 425 G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ;
420 #endif 426 #endif
421 427
422 if ( deltaX <= factor*ctol ) { IsConverg 428 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
423 429
424 } // end iterative loop (i) 430 } // end iterative loop (i)
425 431
426 if ( std::fabs(tmpdist)<ctol ) { tmpdist = 432 if ( std::fabs(tmpdist)<ctol ) { tmpdist = 0 ; }
427 433
428 #ifdef G4TWISTDEBUG 434 #ifdef G4TWISTDEBUG
429 G4cout << "refined solution " << phi << " 435 G4cout << "refined solution " << phi << " , " << u << G4endl ;
430 G4cout << "distance = " << tmpdist << G4en 436 G4cout << "distance = " << tmpdist << G4endl ;
431 G4cout << "local X = " << tmpxx << G4endl 437 G4cout << "local X = " << tmpxx << G4endl ;
432 #endif 438 #endif
433 439
434 tmpisvalid = false ; // init 440 tmpisvalid = false ; // init
435 441
436 if ( IsConverged ) 442 if ( IsConverged )
437 { 443 {
438 if (validate == kValidateWithTol) 444 if (validate == kValidateWithTol)
439 { 445 {
440 tmpareacode = GetAreaCode(tmpxx); 446 tmpareacode = GetAreaCode(tmpxx);
441 if (!IsOutside(tmpareacode)) 447 if (!IsOutside(tmpareacode))
442 { 448 {
443 if (tmpdist >= 0) tmpisvalid = true; 449 if (tmpdist >= 0) tmpisvalid = true;
444 } 450 }
445 } 451 }
446 else if (validate == kValidateWithoutTol 452 else if (validate == kValidateWithoutTol)
447 { 453 {
448 tmpareacode = GetAreaCode(tmpxx, false 454 tmpareacode = GetAreaCode(tmpxx, false);
449 if (IsInside(tmpareacode)) 455 if (IsInside(tmpareacode))
450 { 456 {
451 if (tmpdist >= 0) { tmpisvalid = tru 457 if (tmpdist >= 0) { tmpisvalid = true; }
452 } 458 }
453 } 459 }
454 else // kDontValidate 460 else // kDontValidate
455 { 461 {
456 G4Exception("G4TwistTrapAlphaSide::Dis 462 G4Exception("G4TwistTrapAlphaSide::DistanceToSurface()",
457 "GeomSolids0001", FatalExc 463 "GeomSolids0001", FatalException,
458 "Feature NOT implemented ! 464 "Feature NOT implemented !");
459 } 465 }
460 } 466 }
461 else 467 else
462 { 468 {
463 tmpdist = kInfinity; // no convergen 469 tmpdist = kInfinity; // no convergence after 10 steps
464 tmpisvalid = false ; // solution is 470 tmpisvalid = false ; // solution is not vaild
465 } 471 }
466 472
467 // store the found values 473 // store the found values
468 // 474 //
469 k.xx = tmpxx ; << 475 xbuf[k].xx = tmpxx ;
470 k.distance = tmpdist ; << 476 xbuf[k].distance = tmpdist ;
471 k.areacode = tmpareacode ; << 477 xbuf[k].areacode = tmpareacode ;
472 k.isvalid = tmpisvalid ; << 478 xbuf[k].isvalid = tmpisvalid ;
473 479
474 } // end loop over physical solutions (vari 480 } // end loop over physical solutions (variable k)
475 481
476 std::sort(xbuf.begin() , xbuf.end(), Distanc 482 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting
477 483
478 #ifdef G4TWISTDEBUG 484 #ifdef G4TWISTDEBUG
479 G4cout << G4endl << "list xbuf after sorting 485 G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
480 G4cout << G4endl << G4endl ; 486 G4cout << G4endl << G4endl ;
481 #endif 487 #endif
482 488
483 // erase identical intersection (within kCar 489 // erase identical intersection (within kCarTolerance)
484 // 490 //
485 xbuf.erase( std::unique(xbuf.begin(), xbuf.e 491 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ),
486 xbuf.end() ); 492 xbuf.end() );
487 493
488 494
489 // add guesses 495 // add guesses
490 // 496 //
491 auto nxxtmp = (G4int)xbuf.size() ; << 497 G4int nxxtmp = xbuf.size() ;
492 498
493 if ( nxxtmp<2 || IsParallel ) // positive 499 if ( nxxtmp<2 || IsParallel ) // positive end
494 { 500 {
495 501
496 #ifdef G4TWISTDEBUG 502 #ifdef G4TWISTDEBUG
497 G4cout << "add guess at +z/2 .. " << G4end 503 G4cout << "add guess at +z/2 .. " << G4endl ;
498 #endif 504 #endif
499 505
500 phi = fPhiTwist/2 ; 506 phi = fPhiTwist/2 ;
501 u = 0 ; 507 u = 0 ;
502 508
503 xbuftmp.phi = phi ; 509 xbuftmp.phi = phi ;
504 xbuftmp.u = u ; 510 xbuftmp.u = u ;
505 xbuftmp.areacode = sOutside ; 511 xbuftmp.areacode = sOutside ;
506 xbuftmp.distance = kInfinity ; 512 xbuftmp.distance = kInfinity ;
507 xbuftmp.isvalid = false ; 513 xbuftmp.isvalid = false ;
508 514
509 xbuf.push_back(xbuftmp) ; // store it to 515 xbuf.push_back(xbuftmp) ; // store it to xbuf
510 516
511 #ifdef G4TWISTDEBUG 517 #ifdef G4TWISTDEBUG
512 G4cout << "add guess at -z/2 .. " << G4end 518 G4cout << "add guess at -z/2 .. " << G4endl ;
513 #endif 519 #endif
514 520
515 phi = -fPhiTwist/2 ; 521 phi = -fPhiTwist/2 ;
516 u = 0 ; 522 u = 0 ;
517 523
518 xbuftmp.phi = phi ; 524 xbuftmp.phi = phi ;
519 xbuftmp.u = u ; 525 xbuftmp.u = u ;
520 xbuftmp.areacode = sOutside ; 526 xbuftmp.areacode = sOutside ;
521 xbuftmp.distance = kInfinity ; 527 xbuftmp.distance = kInfinity ;
522 xbuftmp.isvalid = false ; 528 xbuftmp.isvalid = false ;
523 529
524 xbuf.push_back(xbuftmp) ; // store it to 530 xbuf.push_back(xbuftmp) ; // store it to xbuf
525 531
526 for ( std::size_t k = nxxtmp ; k<xbuf.size << 532 for ( size_t k = nxxtmp ; k<xbuf.size() ; ++k )
527 { 533 {
528 534
529 #ifdef G4TWISTDEBUG 535 #ifdef G4TWISTDEBUG
530 G4cout << "Solution " << k << " : " 536 G4cout << "Solution " << k << " : "
531 << "reconstructed phiR = " << xbu 537 << "reconstructed phiR = " << xbuf[k].phi
532 << ", uR = " << xbuf[k].u << G4en 538 << ", uR = " << xbuf[k].u << G4endl ;
533 #endif 539 #endif
534 540
535 phi = xbuf[k].phi ; // get the stored v 541 phi = xbuf[k].phi ; // get the stored values for phi and u
536 u = xbuf[k].u ; 542 u = xbuf[k].u ;
537 543
538 IsConverged = false ; // no convergenc 544 IsConverged = false ; // no convergence at the beginning
539 545
540 for ( G4int i = 1 ; i<maxint ; ++i ) 546 for ( G4int i = 1 ; i<maxint ; ++i )
541 { 547 {
542 xxonsurface = SurfacePoint(phi,u) ; 548 xxonsurface = SurfacePoint(phi,u) ;
543 surfacenormal = NormAng(phi,u) ; 549 surfacenormal = NormAng(phi,u) ;
544 tmpdist = DistanceToPlaneWithV(p, v, x 550 tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx);
545 deltaX = ( tmpxx - xxonsurface ).mag() 551 deltaX = ( tmpxx - xxonsurface ).mag();
546 theta = std::fabs(std::acos(v*surfacen 552 theta = std::fabs(std::acos(v*surfacenormal) - pihalf);
547 if ( theta < 0.001 ) 553 if ( theta < 0.001 )
548 { 554 {
549 factor = 50 ; 555 factor = 50 ;
550 } 556 }
551 else 557 else
552 { 558 {
553 factor = 1 ; 559 factor = 1 ;
554 } 560 }
555 561
556 #ifdef G4TWISTDEBUG 562 #ifdef G4TWISTDEBUG
557 G4cout << "Step i = " << i << ", dista 563 G4cout << "Step i = " << i << ", distance = " << tmpdist
558 << ", " << deltaX << G4endl 564 << ", " << deltaX << G4endl
559 << "X = " << tmpxx << G4endl ; 565 << "X = " << tmpxx << G4endl ;
560 #endif 566 #endif
561 567
562 GetPhiUAtX(tmpxx, phi, u) ; 568 GetPhiUAtX(tmpxx, phi, u) ;
563 // the new point xx is accepted and 569 // the new point xx is accepted and phi/u replaced
564 570
565 #ifdef G4TWISTDEBUG 571 #ifdef G4TWISTDEBUG
566 G4cout << "approximated phi = " << phi 572 G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ;
567 #endif 573 #endif
568 574
569 if ( deltaX <= factor*ctol ) { IsConve 575 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
570 576
571 } // end iterative loop (i) 577 } // end iterative loop (i)
572 578
573 if ( std::fabs(tmpdist)<ctol ) { tmpdis 579 if ( std::fabs(tmpdist)<ctol ) { tmpdist = 0; }
574 580
575 #ifdef G4TWISTDEBUG 581 #ifdef G4TWISTDEBUG
576 G4cout << "refined solution " << phi << 582 G4cout << "refined solution " << phi << " , " << u << G4endl ;
577 G4cout << "distance = " << tmpdist << G4 583 G4cout << "distance = " << tmpdist << G4endl ;
578 G4cout << "local X = " << tmpxx << G4end 584 G4cout << "local X = " << tmpxx << G4endl ;
579 #endif 585 #endif
580 586
581 tmpisvalid = false ; // init 587 tmpisvalid = false ; // init
582 588
583 if ( IsConverged ) 589 if ( IsConverged )
584 { 590 {
585 if (validate == kValidateWithTol) 591 if (validate == kValidateWithTol)
586 { 592 {
587 tmpareacode = GetAreaCode(tmpxx); 593 tmpareacode = GetAreaCode(tmpxx);
588 if (!IsOutside(tmpareacode)) 594 if (!IsOutside(tmpareacode))
589 { 595 {
590 if (tmpdist >= 0) { tmpisvalid = 596 if (tmpdist >= 0) { tmpisvalid = true; }
591 } 597 }
592 } 598 }
593 else if (validate == kValidateWithoutT 599 else if (validate == kValidateWithoutTol)
594 { 600 {
595 tmpareacode = GetAreaCode(tmpxx, fal 601 tmpareacode = GetAreaCode(tmpxx, false);
596 if (IsInside(tmpareacode)) 602 if (IsInside(tmpareacode))
597 { 603 {
598 if (tmpdist >= 0) { tmpisvalid = 604 if (tmpdist >= 0) { tmpisvalid = true; }
599 } 605 }
600 } 606 }
601 else // kDontValidate 607 else // kDontValidate
602 { 608 {
603 G4Exception("G4TwistedBoxSide::Dista 609 G4Exception("G4TwistedBoxSide::DistanceToSurface()",
604 "GeomSolids0001", FatalE 610 "GeomSolids0001", FatalException,
605 "Feature NOT implemented 611 "Feature NOT implemented !");
606 } 612 }
607 } 613 }
608 else 614 else
609 { 615 {
610 tmpdist = kInfinity; // no converg 616 tmpdist = kInfinity; // no convergence after 10 steps
611 tmpisvalid = false ; // solution i 617 tmpisvalid = false ; // solution is not vaild
612 } 618 }
613 619
614 // store the found values 620 // store the found values
615 // 621 //
616 xbuf[k].xx = tmpxx ; 622 xbuf[k].xx = tmpxx ;
617 xbuf[k].distance = tmpdist ; 623 xbuf[k].distance = tmpdist ;
618 xbuf[k].areacode = tmpareacode ; 624 xbuf[k].areacode = tmpareacode ;
619 xbuf[k].isvalid = tmpisvalid ; 625 xbuf[k].isvalid = tmpisvalid ;
620 626
621 } // end loop over physical solutions 627 } // end loop over physical solutions
622 } // end less than 2 solutions 628 } // end less than 2 solutions
623 629
624 // sort again 630 // sort again
625 std::sort(xbuf.begin() , xbuf.end(), Distanc 631 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting
626 632
627 // erase identical intersection (within kCar 633 // erase identical intersection (within kCarTolerance)
628 xbuf.erase( std::unique(xbuf.begin(), xbuf.e 634 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) ,
629 xbuf.end() ); 635 xbuf.end() );
630 636
631 #ifdef G4TWISTDEBUG 637 #ifdef G4TWISTDEBUG
632 G4cout << G4endl << "list xbuf after sorting 638 G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
633 G4cout << G4endl << G4endl ; 639 G4cout << G4endl << G4endl ;
634 #endif 640 #endif
635 641
636 nxx = (G4int)xbuf.size() ; // determine nu << 642 nxx = xbuf.size() ; // determine number of solutions again.
637 643
638 for ( G4int i = 0 ; i<(G4int)xbuf.size() ; + << 644 for ( size_t i = 0 ; i<xbuf.size() ; ++i )
639 { 645 {
640 distance[i] = xbuf[i].distance; 646 distance[i] = xbuf[i].distance;
641 gxx[i] = ComputeGlobalPoint(xbuf[i].x 647 gxx[i] = ComputeGlobalPoint(xbuf[i].xx);
642 areacode[i] = xbuf[i].areacode ; 648 areacode[i] = xbuf[i].areacode ;
643 isvalid[i] = xbuf[i].isvalid ; 649 isvalid[i] = xbuf[i].isvalid ;
644 650
645 fCurStatWithV.SetCurrentStatus(i, gxx[i], 651 fCurStatWithV.SetCurrentStatus(i, gxx[i], distance[i], areacode[i],
646 isvalid[i 652 isvalid[i], nxx, validate, &gp, &gv);
647 #ifdef G4TWISTDEBUG 653 #ifdef G4TWISTDEBUG
648 G4cout << "element Nr. " << i 654 G4cout << "element Nr. " << i
649 << ", local Intersection = " << xbu 655 << ", local Intersection = " << xbuf[i].xx
650 << ", distance = " << xbuf[i].dista 656 << ", distance = " << xbuf[i].distance
651 << ", u = " << xbuf[i].u 657 << ", u = " << xbuf[i].u
652 << ", phi = " << xbuf[i].phi 658 << ", phi = " << xbuf[i].phi
653 << ", isvalid = " << xbuf[i].isvali 659 << ", isvalid = " << xbuf[i].isvalid
654 << G4endl ; 660 << G4endl ;
655 #endif 661 #endif
656 662
657 } // end for( i ) loop 663 } // end for( i ) loop
658 664
659 #ifdef G4TWISTDEBUG 665 #ifdef G4TWISTDEBUG
660 G4cout << "G4TwistTrapAlphaSide finished " < 666 G4cout << "G4TwistTrapAlphaSide finished " << G4endl ;
661 G4cout << nxx << " possible physical solutio 667 G4cout << nxx << " possible physical solutions found" << G4endl ;
662 for ( G4int k= 0 ; k< nxx ; k++ ) 668 for ( G4int k= 0 ; k< nxx ; k++ )
663 { 669 {
664 G4cout << "global intersection Point found 670 G4cout << "global intersection Point found: " << gxx[k] << G4endl ;
665 G4cout << "distance = " << distance[k] << 671 G4cout << "distance = " << distance[k] << G4endl ;
666 G4cout << "isvalid = " << isvalid[k] << G4 672 G4cout << "isvalid = " << isvalid[k] << G4endl ;
667 } 673 }
668 #endif 674 #endif
669 675
670 return nxx ; 676 return nxx ;
671 } 677 }
672 678
673 679
674 //============================================ 680 //=====================================================================
675 //* DistanceToSurface ------------------------ 681 //* DistanceToSurface -------------------------------------------------
676 682
677 G4int 683 G4int
678 G4TwistTrapAlphaSide::DistanceToSurface(const 684 G4TwistTrapAlphaSide::DistanceToSurface(const G4ThreeVector& gp,
679 685 G4ThreeVector gxx[],
680 686 G4double distance[],
681 687 G4int areacode[])
682 { 688 {
683 const G4double ctol = 0.5 * kCarTolerance; 689 const G4double ctol = 0.5 * kCarTolerance;
684 690
685 fCurStat.ResetfDone(kDontValidate, &gp); 691 fCurStat.ResetfDone(kDontValidate, &gp);
686 692
687 if (fCurStat.IsDone()) 693 if (fCurStat.IsDone())
688 { 694 {
689 for (G4int i=0; i<fCurStat.GetNXX(); ++i 695 for (G4int i=0; i<fCurStat.GetNXX(); ++i)
690 { 696 {
691 gxx[i] = fCurStat.GetXX(i); 697 gxx[i] = fCurStat.GetXX(i);
692 distance[i] = fCurStat.GetDistance(i) 698 distance[i] = fCurStat.GetDistance(i);
693 areacode[i] = fCurStat.GetAreacode(i) 699 areacode[i] = fCurStat.GetAreacode(i);
694 } 700 }
695 return fCurStat.GetNXX(); 701 return fCurStat.GetNXX();
696 } 702 }
697 else // initialize 703 else // initialize
698 { 704 {
699 for (G4int i=0; i<G4VSURFACENXX; ++i) 705 for (G4int i=0; i<G4VSURFACENXX; ++i)
700 { 706 {
701 distance[i] = kInfinity; 707 distance[i] = kInfinity;
702 areacode[i] = sOutside; 708 areacode[i] = sOutside;
703 gxx[i].set(kInfinity, kInfinity, kInf 709 gxx[i].set(kInfinity, kInfinity, kInfinity);
704 } 710 }
705 } 711 }
706 712
707 G4ThreeVector p = ComputeLocalPoint(gp); 713 G4ThreeVector p = ComputeLocalPoint(gp);
708 G4ThreeVector xx; // intersection point 714 G4ThreeVector xx; // intersection point
709 G4ThreeVector xxonsurface ; // interpolated 715 G4ThreeVector xxonsurface ; // interpolated intersection point
710 716
711 // the surfacenormal at that surface point 717 // the surfacenormal at that surface point
712 // 718 //
713 G4double phiR = 0 ; 719 G4double phiR = 0 ;
714 G4double uR = 0 ; 720 G4double uR = 0 ;
715 721
716 G4ThreeVector surfacenormal ; 722 G4ThreeVector surfacenormal ;
717 G4double deltaX, uMax ; 723 G4double deltaX, uMax ;
718 G4double halfphi = 0.5*fPhiTwist ; 724 G4double halfphi = 0.5*fPhiTwist ;
719 725
720 for ( G4int i = 1 ; i<20 ; ++i ) 726 for ( G4int i = 1 ; i<20 ; ++i )
721 { 727 {
722 xxonsurface = SurfacePoint(phiR,uR) ; 728 xxonsurface = SurfacePoint(phiR,uR) ;
723 surfacenormal = NormAng(phiR,uR) ; 729 surfacenormal = NormAng(phiR,uR) ;
724 distance[0] = DistanceToPlane(p,xxonsurfa 730 distance[0] = DistanceToPlane(p,xxonsurface,surfacenormal,xx); // new XX
725 deltaX = ( xx - xxonsurface ).mag() ; 731 deltaX = ( xx - xxonsurface ).mag() ;
726 732
727 #ifdef G4TWISTDEBUG 733 #ifdef G4TWISTDEBUG
728 G4cout << "i = " << i << ", distance = " 734 G4cout << "i = " << i << ", distance = " << distance[0]
729 << ", " << deltaX << G4endl 735 << ", " << deltaX << G4endl
730 << "X = " << xx << G4endl ; 736 << "X = " << xx << G4endl ;
731 #endif 737 #endif
732 738
733 // the new point xx is accepted and phi/p 739 // the new point xx is accepted and phi/psi replaced
734 // 740 //
735 GetPhiUAtX(xx, phiR, uR) ; 741 GetPhiUAtX(xx, phiR, uR) ;
736 742
737 if ( deltaX <= ctol ) { break ; } 743 if ( deltaX <= ctol ) { break ; }
738 } 744 }
739 745
740 // check validity of solution ( valid phi,p 746 // check validity of solution ( valid phi,psi )
741 747
742 uMax = GetBoundaryMax(phiR) ; 748 uMax = GetBoundaryMax(phiR) ;
743 749
744 if ( phiR > halfphi ) { phiR = halfphi ; 750 if ( phiR > halfphi ) { phiR = halfphi ; }
745 if ( phiR < -halfphi ) { phiR = -halfphi ; 751 if ( phiR < -halfphi ) { phiR = -halfphi ; }
746 if ( uR > uMax ) { uR = uMax ; } 752 if ( uR > uMax ) { uR = uMax ; }
747 if ( uR < -uMax ) { uR = -uMax ; } 753 if ( uR < -uMax ) { uR = -uMax ; }
748 754
749 xxonsurface = SurfacePoint(phiR,uR) ; 755 xxonsurface = SurfacePoint(phiR,uR) ;
750 distance[0] = ( p - xx ).mag() ; 756 distance[0] = ( p - xx ).mag() ;
751 if ( distance[0] <= ctol ) { distance[0] = 757 if ( distance[0] <= ctol ) { distance[0] = 0 ; }
752 758
753 // end of validity 759 // end of validity
754 760
755 #ifdef G4TWISTDEBUG 761 #ifdef G4TWISTDEBUG
756 G4cout << "refined solution " << phiR << " 762 G4cout << "refined solution " << phiR << " , " << uR << " , " << G4endl ;
757 G4cout << "distance = " << distance[0] << G 763 G4cout << "distance = " << distance[0] << G4endl ;
758 G4cout << "X = " << xx << G4endl ; 764 G4cout << "X = " << xx << G4endl ;
759 #endif 765 #endif
760 766
761 G4bool isvalid = true; 767 G4bool isvalid = true;
762 gxx[0] = ComputeGlobalPoint(xx); 768 gxx[0] = ComputeGlobalPoint(xx);
763 769
764 #ifdef G4TWISTDEBUG 770 #ifdef G4TWISTDEBUG
765 G4cout << "intersection Point found: " << g 771 G4cout << "intersection Point found: " << gxx[0] << G4endl ;
766 G4cout << "distance = " << distance[0] << G 772 G4cout << "distance = " << distance[0] << G4endl ;
767 #endif 773 #endif
768 774
769 fCurStat.SetCurrentStatus(0, gxx[0], distan 775 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
770 isvalid, 1, kDontV 776 isvalid, 1, kDontValidate, &gp);
771 return 1; 777 return 1;
772 } 778 }
773 779
774 780
775 //============================================ 781 //=====================================================================
776 //* GetAreaCode ------------------------------ 782 //* GetAreaCode -------------------------------------------------------
777 783
778 G4int 784 G4int
779 G4TwistTrapAlphaSide::GetAreaCode(const G4Thre 785 G4TwistTrapAlphaSide::GetAreaCode(const G4ThreeVector& xx, G4bool withTol)
780 { 786 {
781 // We must use the function in local coordi 787 // We must use the function in local coordinate system.
782 // See the description of DistanceToSurface 788 // See the description of DistanceToSurface(p,v).
783 789
784 const G4double ctol = 0.5 * kCarTolerance; 790 const G4double ctol = 0.5 * kCarTolerance;
785 791
786 G4double phi ; 792 G4double phi ;
787 G4double yprime ; 793 G4double yprime ;
788 GetPhiUAtX(xx, phi,yprime ) ; 794 GetPhiUAtX(xx, phi,yprime ) ;
789 795
790 G4double fYAxisMax = GetBoundaryMax(phi) ; 796 G4double fYAxisMax = GetBoundaryMax(phi) ;
791 G4double fYAxisMin = GetBoundaryMin(phi) ; 797 G4double fYAxisMin = GetBoundaryMin(phi) ;
792 798
793 #ifdef G4TWISTDEBUG 799 #ifdef G4TWISTDEBUG
794 G4cout << "GetAreaCode: phi = " << phi << G 800 G4cout << "GetAreaCode: phi = " << phi << G4endl ;
795 G4cout << "GetAreaCode: yprime = " << yprim 801 G4cout << "GetAreaCode: yprime = " << yprime << G4endl ;
796 G4cout << "Intervall is " << fYAxisMin << " 802 G4cout << "Intervall is " << fYAxisMin << " to " << fYAxisMax << G4endl ;
797 #endif 803 #endif
798 804
799 G4int areacode = sInside; 805 G4int areacode = sInside;
800 806
801 if (fAxis[0] == kYAxis && fAxis[1] == kZAxi 807 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
802 { 808 {
803 G4int zaxis = 1; 809 G4int zaxis = 1;
804 810
805 if (withTol) 811 if (withTol)
806 { 812 {
807 G4bool isoutside = false; 813 G4bool isoutside = false;
808 814
809 // test boundary of yaxis 815 // test boundary of yaxis
810 816
811 if (yprime < fYAxisMin + ctol) 817 if (yprime < fYAxisMin + ctol)
812 { 818 {
813 areacode |= (sAxis0 & (sAxisY | sA 819 areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
814 if (yprime <= fYAxisMin - ctol) { 820 if (yprime <= fYAxisMin - ctol) { isoutside = true; }
815 821
816 } 822 }
817 else if (yprime > fYAxisMax - ctol) 823 else if (yprime > fYAxisMax - ctol)
818 { 824 {
819 areacode |= (sAxis0 & (sAxisY | sA 825 areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
820 if (yprime >= fYAxisMax + ctol) { 826 if (yprime >= fYAxisMax + ctol) { isoutside = true; }
821 } 827 }
822 828
823 // test boundary of z-axis 829 // test boundary of z-axis
824 830
825 if (xx.z() < fAxisMin[zaxis] + ctol) 831 if (xx.z() < fAxisMin[zaxis] + ctol)
826 { 832 {
827 areacode |= (sAxis1 & (sAxisZ | sA 833 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
828 834
829 if ((areacode & sBoundary) != 0) << 835 if (areacode & sBoundary) // xx is on the corner
830 { areacode |= sCorner; } 836 { areacode |= sCorner; }
831 837
832 else 838 else
833 { areacode |= sBoundary; } 839 { areacode |= sBoundary; }
834 if (xx.z() <= fAxisMin[zaxis] - ct 840 if (xx.z() <= fAxisMin[zaxis] - ctol) { isoutside = true; }
835 } 841 }
836 else if (xx.z() > fAxisMax[zaxis] - c 842 else if (xx.z() > fAxisMax[zaxis] - ctol)
837 { 843 {
838 areacode |= (sAxis1 & (sAxisZ | sA 844 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
839 845
840 if ((areacode & sBoundary) != 0) << 846 if (areacode & sBoundary) // xx is on the corner
841 { areacode |= sCorner; } 847 { areacode |= sCorner; }
842 else 848 else
843 { areacode |= sBoundary; } 849 { areacode |= sBoundary; }
844 if (xx.z() >= fAxisMax[zaxis] + ct 850 if (xx.z() >= fAxisMax[zaxis] + ctol) { isoutside = true; }
845 } 851 }
846 852
847 // if isoutside = true, clear inside 853 // if isoutside = true, clear inside bit.
848 // if not on boundary, add axis infor 854 // if not on boundary, add axis information.
849 855
850 if (isoutside) 856 if (isoutside)
851 { 857 {
852 G4int tmpareacode = areacode & (~s 858 G4int tmpareacode = areacode & (~sInside);
853 areacode = tmpareacode; 859 areacode = tmpareacode;
854 } 860 }
855 else if ((areacode & sBoundary) != sB 861 else if ((areacode & sBoundary) != sBoundary)
856 { 862 {
857 areacode |= (sAxis0 & sAxisY) | (s 863 areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
858 } 864 }
859 865
860 } 866 }
861 else 867 else
862 { 868 {
863 // boundary of y-axis 869 // boundary of y-axis
864 870
865 if (yprime < fYAxisMin ) 871 if (yprime < fYAxisMin )
866 { 872 {
867 areacode |= (sAxis0 & (sAxisY | sA 873 areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
868 } 874 }
869 else if (yprime > fYAxisMax) 875 else if (yprime > fYAxisMax)
870 { 876 {
871 areacode |= (sAxis0 & (sAxisY | sA 877 areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
872 } 878 }
873 879
874 // boundary of z-axis 880 // boundary of z-axis
875 881
876 if (xx.z() < fAxisMin[zaxis]) 882 if (xx.z() < fAxisMin[zaxis])
877 { 883 {
878 areacode |= (sAxis1 & (sAxisZ | sA 884 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
879 if ((areacode & sBoundary) != 0) << 885 if (areacode & sBoundary) // xx is on the corner
880 { areacode |= sCorner; } 886 { areacode |= sCorner; }
881 else 887 else
882 { areacode |= sBoundary; } 888 { areacode |= sBoundary; }
883 } 889 }
884 else if (xx.z() > fAxisMax[zaxis]) 890 else if (xx.z() > fAxisMax[zaxis])
885 { 891 {
886 areacode |= (sAxis1 & (sAxisZ | sA 892 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
887 if ((areacode & sBoundary) != 0) << 893 if (areacode & sBoundary) // xx is on the corner
888 { areacode |= sCorner; } 894 { areacode |= sCorner; }
889 else 895 else
890 { areacode |= sBoundary; } 896 { areacode |= sBoundary; }
891 } 897 }
892 898
893 if ((areacode & sBoundary) != sBounda 899 if ((areacode & sBoundary) != sBoundary)
894 { 900 {
895 areacode |= (sAxis0 & sAxisY) | (s 901 areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
896 } 902 }
897 } 903 }
898 return areacode; 904 return areacode;
899 } 905 }
900 else 906 else
901 { 907 {
902 G4Exception("G4TwistTrapAlphaSide::GetAr 908 G4Exception("G4TwistTrapAlphaSide::GetAreaCode()",
903 "GeomSolids0001", FatalExcep 909 "GeomSolids0001", FatalException,
904 "Feature NOT implemented !") 910 "Feature NOT implemented !");
905 } 911 }
906 return areacode; 912 return areacode;
907 } 913 }
908 914
909 //============================================ 915 //=====================================================================
910 //* SetCorners() ----------------------------- 916 //* SetCorners() ------------------------------------------------------
911 917
912 void G4TwistTrapAlphaSide::SetCorners() 918 void G4TwistTrapAlphaSide::SetCorners()
913 { 919 {
914 920
915 // Set Corner points in local coodinate. 921 // Set Corner points in local coodinate.
916 922
917 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis 923 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
918 { 924 {
919 925
920 G4double x, y, z; 926 G4double x, y, z;
921 927
922 // corner of Axis0min and Axis1min 928 // corner of Axis0min and Axis1min
923 // 929 //
924 x = -fdeltaX/2. + (fDx1 - fDy1*fTAlph)*std 930 x = -fdeltaX/2. + (fDx1 - fDy1*fTAlph)*std::cos(fPhiTwist/2.)
925 - fDy1*std::sin(fPhiTwist/2.); 931 - fDy1*std::sin(fPhiTwist/2.);
926 y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/ 932 y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/2.)
927 + (-fDx1 + fDy1*fTAlph)*std::sin(fPhiTwi 933 + (-fDx1 + fDy1*fTAlph)*std::sin(fPhiTwist/2.);
928 z = -fDz ; 934 z = -fDz ;
929 935
930 // G4cout << "SetCorners: " << x << ", 936 // G4cout << "SetCorners: " << x << ", " << y << ", " << z << G4endl ;
931 937
932 SetCorner(sC0Min1Min, x, y, z); 938 SetCorner(sC0Min1Min, x, y, z);
933 939
934 // corner of Axis0max and Axis1min 940 // corner of Axis0max and Axis1min
935 // 941 //
936 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std 942 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.)
937 + fDy1*std::sin(fPhiTwist/2.); 943 + fDy1*std::sin(fPhiTwist/2.);
938 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/ 944 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.)
939 - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwis 945 - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.);
940 z = -fDz ; 946 z = -fDz ;
941 947
942 // G4cout << "SetCorners: " << x << ", 948 // G4cout << "SetCorners: " << x << ", " << y << ", " << z << G4endl ;
943 949
944 SetCorner(sC0Max1Min, x, y, z); 950 SetCorner(sC0Max1Min, x, y, z);
945 951
946 // corner of Axis0max and Axis1max 952 // corner of Axis0max and Axis1max
947 // 953 //
948 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std: 954 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.)
949 - fDy2*std::sin(fPhiTwist/2.); 955 - fDy2*std::sin(fPhiTwist/2.);
950 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2 956 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.)
951 + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwis 957 + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.);
952 z = fDz ; 958 z = fDz ;
953 959
954 // G4cout << "SetCorners: " << x << ", 960 // G4cout << "SetCorners: " << x << ", " << y << ", " << z << G4endl ;
955 961
956 SetCorner(sC0Max1Max, x, y, z); 962 SetCorner(sC0Max1Max, x, y, z);
957 963
958 // corner of Axis0min and Axis1max 964 // corner of Axis0min and Axis1max
959 x = fdeltaX/2. + (fDx3 - fDy2*fTAlph)*std: 965 x = fdeltaX/2. + (fDx3 - fDy2*fTAlph)*std::cos(fPhiTwist/2.)
960 + fDy2*std::sin(fPhiTwist/2.) ; 966 + fDy2*std::sin(fPhiTwist/2.) ;
961 y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2 967 y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2.)
962 + (fDx3 - fDy2*fTAlph)*std::sin(fPhiTwis 968 + (fDx3 - fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
963 z = fDz ; 969 z = fDz ;
964 970
965 // G4cout << "SetCorners: " << x << ", 971 // G4cout << "SetCorners: " << x << ", " << y << ", " << z << G4endl ;
966 972
967 SetCorner(sC0Min1Max, x, y, z); 973 SetCorner(sC0Min1Max, x, y, z);
968 974
969 } 975 }
970 else 976 else
971 { 977 {
972 G4Exception("G4TwistTrapAlphaSide::SetCorn 978 G4Exception("G4TwistTrapAlphaSide::SetCorners()",
973 "GeomSolids0001", FatalExcepti 979 "GeomSolids0001", FatalException,
974 "Method NOT implemented !"); 980 "Method NOT implemented !");
975 } 981 }
976 } 982 }
977 983
978 //============================================ 984 //=====================================================================
979 //* SetBoundaries() -------------------------- 985 //* SetBoundaries() ---------------------------------------------------
980 986
981 void G4TwistTrapAlphaSide::SetBoundaries() 987 void G4TwistTrapAlphaSide::SetBoundaries()
982 { 988 {
983 // Set direction-unit vector of boundary-li 989 // Set direction-unit vector of boundary-lines in local coodinate.
984 // 990 //
985 991
986 G4ThreeVector direction; 992 G4ThreeVector direction;
987 993
988 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis 994 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
989 { 995 {
990 // sAxis0 & sAxisMin 996 // sAxis0 & sAxisMin
991 direction = GetCorner(sC0Min1Max) - GetCor 997 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
992 direction = direction.unit(); 998 direction = direction.unit();
993 SetBoundary(sAxis0 & (sAxisY | sAxisMin), 999 SetBoundary(sAxis0 & (sAxisY | sAxisMin), direction,
994 GetCorner(sC0Min1Min), sAxisZ) 1000 GetCorner(sC0Min1Min), sAxisZ) ;
995 1001
996 // sAxis0 & sAxisMax 1002 // sAxis0 & sAxisMax
997 direction = GetCorner(sC0Max1Max) - GetCor 1003 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
998 direction = direction.unit(); 1004 direction = direction.unit();
999 SetBoundary(sAxis0 & (sAxisY | sAxisMax), 1005 SetBoundary(sAxis0 & (sAxisY | sAxisMax), direction,
1000 GetCorner(sC0Max1Min), sAxisZ 1006 GetCorner(sC0Max1Min), sAxisZ);
1001 1007
1002 // sAxis1 & sAxisMin 1008 // sAxis1 & sAxisMin
1003 direction = GetCorner(sC0Max1Min) - GetCo 1009 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
1004 direction = direction.unit(); 1010 direction = direction.unit();
1005 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), 1011 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
1006 GetCorner(sC0Min1Min), sAxisY 1012 GetCorner(sC0Min1Min), sAxisY);
1007 1013
1008 // sAxis1 & sAxisMax 1014 // sAxis1 & sAxisMax
1009 direction = GetCorner(sC0Max1Max) - GetCo 1015 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
1010 direction = direction.unit(); 1016 direction = direction.unit();
1011 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), 1017 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
1012 GetCorner(sC0Min1Max), sAxisY 1018 GetCorner(sC0Min1Max), sAxisY);
1013 1019
1014 } 1020 }
1015 else 1021 else
1016 { 1022 {
1017 G4Exception("G4TwistTrapAlphaSide::SetCor 1023 G4Exception("G4TwistTrapAlphaSide::SetCorners()",
1018 "GeomSolids0001", FatalExcept 1024 "GeomSolids0001", FatalException,
1019 "Feature NOT implemented !"); 1025 "Feature NOT implemented !");
1020 } 1026 }
1021 } 1027 }
1022 1028
1023 //=========================================== 1029 //=====================================================================
1024 //* GetPhiUAtX ------------------------------ 1030 //* GetPhiUAtX --------------------------------------------------------
1025 1031
1026 void 1032 void
1027 G4TwistTrapAlphaSide::GetPhiUAtX( const G4Thr << 1033 G4TwistTrapAlphaSide::GetPhiUAtX( G4ThreeVector p, G4double& phi, G4double& u )
1028 { 1034 {
1029 // find closest point XX on surface for a g 1035 // find closest point XX on surface for a given point p
1030 // X0 is a point on the surface, d is the 1036 // X0 is a point on the surface, d is the direction
1031 // ( both for a fixed z = pz) 1037 // ( both for a fixed z = pz)
1032 1038
1033 // phi is given by the z coordinate of p 1039 // phi is given by the z coordinate of p
1034 1040
1035 phi = p.z()/(2*fDz)*fPhiTwist ; 1041 phi = p.z()/(2*fDz)*fPhiTwist ;
1036 u = (fPhiTwist*(2*fDx1*fDx1 - 2*fDx2*fDx2 - 1042 u = (fPhiTwist*(2*fDx1*fDx1 - 2*fDx2*fDx2 - fa1md1*(fDx3 + fDx4)
1037 - 4*(fDx3plus1 + fDx4plus2)*fDy1 1043 - 4*(fDx3plus1 + fDx4plus2)*fDy1*fTAlph)
1038 - 2*(2*fDx1*fDx1 - 2*fDx2*fDx2 + 1044 - 2*(2*fDx1*fDx1 - 2*fDx2*fDx2 + fa1md1*(fDx3 + fDx4)
1039 + 4*(fDx3minus1 + fDx4minus2) 1045 + 4*(fDx3minus1 + fDx4minus2)*fDy1*fTAlph)*phi
1040 - 4*(fa1md1*(fdeltaX*phi - fPhiT 1046 - 4*(fa1md1*(fdeltaX*phi - fPhiTwist*p.x())
1041 + 4*fDy1*(fdeltaY*phi + fdeltaX* 1047 + 4*fDy1*(fdeltaY*phi + fdeltaX*fTAlph*phi
1042 - fPhiTwist*(fTAlph*p.x( 1048 - fPhiTwist*(fTAlph*p.x() + p.y())))*std::cos(phi)
1043 - 4*(fa1md1*fdeltaY*phi - 4*fdel 1049 - 4*(fa1md1*fdeltaY*phi - 4*fdeltaX*fDy1*phi
1044 + 4*fdeltaY*fDy1*fTAlph*phi + 1050 + 4*fdeltaY*fDy1*fTAlph*phi + 4*fDy1*fPhiTwist*p.x()
1045 - fPhiTwist*(fa1md1 + 4*fDy1* 1051 - fPhiTwist*(fa1md1 + 4*fDy1*fTAlph)*p.y())*std::sin(phi))
1046 /(fDy1* fPhiTwist*((std::fabs(((fa1md1 1052 /(fDy1* fPhiTwist*((std::fabs(((fa1md1 + 4*fDy1*fTAlph)*std::cos(phi))
1047 /fDy1 - 4*std::sin( 1053 /fDy1 - 4*std::sin(phi)))
1048 *(std::fabs(((fa1md1 1054 *(std::fabs(((fa1md1 + 4*fDy1*fTAlph)*std::cos(phi))
1049 /fDy1 - 4*std::sin( 1055 /fDy1 - 4*std::sin(phi)))
1050 + (std::fabs(4*std::cos(phi) 1056 + (std::fabs(4*std::cos(phi)
1051 + ((fa1md1 + 4*fDy1*fTAlph)*s 1057 + ((fa1md1 + 4*fDy1*fTAlph)*std::sin(phi))/fDy1))
1052 * (std::fabs(4*std::cos(phi) 1058 * (std::fabs(4*std::cos(phi)
1053 + ((fa1md1 + 4*fDy1*fTAlph)*s 1059 + ((fa1md1 + 4*fDy1*fTAlph)*std::sin(phi))/fDy1)))) ;
1054 } 1060 }
1055 1061
1056 //=========================================== 1062 //=====================================================================
1057 //* ProjectPoint ---------------------------- 1063 //* ProjectPoint ------------------------------------------------------
1058 1064
1059 G4ThreeVector 1065 G4ThreeVector
1060 G4TwistTrapAlphaSide::ProjectPoint(const G4Th 1066 G4TwistTrapAlphaSide::ProjectPoint(const G4ThreeVector& p, G4bool isglobal)
1061 { 1067 {
1062 // Get Rho at p.z() on Hyperbolic Surface. 1068 // Get Rho at p.z() on Hyperbolic Surface.
1063 1069
1064 G4ThreeVector tmpp; 1070 G4ThreeVector tmpp;
1065 if (isglobal) 1071 if (isglobal)
1066 { 1072 {
1067 tmpp = fRot.inverse()*p - fTrans; 1073 tmpp = fRot.inverse()*p - fTrans;
1068 } 1074 }
1069 else 1075 else
1070 { 1076 {
1071 tmpp = p; 1077 tmpp = p;
1072 } 1078 }
1073 1079
1074 G4double phi ; 1080 G4double phi ;
1075 G4double u ; 1081 G4double u ;
1076 1082
1077 GetPhiUAtX( tmpp, phi, u ) ; 1083 GetPhiUAtX( tmpp, phi, u ) ;
1078 // calculate (phi, u) for a point p close 1084 // calculate (phi, u) for a point p close the surface
1079 1085
1080 G4ThreeVector xx = SurfacePoint(phi,u) ; 1086 G4ThreeVector xx = SurfacePoint(phi,u) ;
1081 // transform back to cartesian coordinate 1087 // transform back to cartesian coordinates
1082 1088
1083 if (isglobal) 1089 if (isglobal)
1084 { 1090 {
1085 return (fRot * xx + fTrans); 1091 return (fRot * xx + fTrans);
1086 } 1092 }
1087 else 1093 else
1088 { 1094 {
1089 return xx; 1095 return xx;
1090 } 1096 }
1091 } 1097 }
1092 1098
1093 //=========================================== 1099 //=====================================================================
1094 //* GetFacets ------------------------------- 1100 //* GetFacets ---------------------------------------------------------
1095 1101
1096 void 1102 void
1097 G4TwistTrapAlphaSide::GetFacets( G4int k, G4i 1103 G4TwistTrapAlphaSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
1098 G4int faces[ 1104 G4int faces[][4], G4int iside )
1099 { 1105 {
1100 1106
1101 G4double phi ; 1107 G4double phi ;
1102 G4double b ; 1108 G4double b ;
1103 1109
1104 G4double z, u ; // the two parameters f 1110 G4double z, u ; // the two parameters for the surface equation
1105 G4ThreeVector p ; // a point on the surfac 1111 G4ThreeVector p ; // a point on the surface, given by (z,u)
1106 1112
1107 G4int nnode ; 1113 G4int nnode ;
1108 G4int nface ; 1114 G4int nface ;
1109 1115
1110 // calculate the (n-1)*(k-1) vertices 1116 // calculate the (n-1)*(k-1) vertices
1111 1117
1112 for ( G4int i = 0 ; i<n ; ++i ) 1118 for ( G4int i = 0 ; i<n ; ++i )
1113 { 1119 {
1114 z = -fDz+i*(2.*fDz)/(n-1) ; 1120 z = -fDz+i*(2.*fDz)/(n-1) ;
1115 phi = z*fPhiTwist/(2*fDz) ; 1121 phi = z*fPhiTwist/(2*fDz) ;
1116 b = GetValueB(phi) ; 1122 b = GetValueB(phi) ;
1117 1123
1118 for ( G4int j = 0 ; j<k ; ++j ) 1124 for ( G4int j = 0 ; j<k ; ++j )
1119 { 1125 {
1120 nnode = GetNode(i,j,k,n,iside) ; 1126 nnode = GetNode(i,j,k,n,iside) ;
1121 u = -b/2 +j*b/(k-1) ; 1127 u = -b/2 +j*b/(k-1) ;
1122 p = SurfacePoint(phi,u,true) ; // surf 1128 p = SurfacePoint(phi,u,true) ; // surface point in global coordinates
1123 1129
1124 xyz[nnode][0] = p.x() ; 1130 xyz[nnode][0] = p.x() ;
1125 xyz[nnode][1] = p.y() ; 1131 xyz[nnode][1] = p.y() ;
1126 xyz[nnode][2] = p.z() ; 1132 xyz[nnode][2] = p.z() ;
1127 1133
1128 if ( i<n-1 && j<k-1 ) // conterclock w 1134 if ( i<n-1 && j<k-1 ) // conterclock wise filling
1129 { 1135 {
1130 nface = GetFace(i,j,k,n,iside) ; 1136 nface = GetFace(i,j,k,n,iside) ;
1131 faces[nface][0] = GetEdgeVisibility(i 1137 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1)
1132 * (GetNode(i ,j ,k, 1138 * (GetNode(i ,j ,k,n,iside)+1) ; // f77 numbering
1133 faces[nface][1] = GetEdgeVisibility(i 1139 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1)
1134 * (GetNode(i ,j+1,k, 1140 * (GetNode(i ,j+1,k,n,iside)+1) ;
1135 faces[nface][2] = GetEdgeVisibility(i 1141 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1)
1136 * (GetNode(i+1,j+1,k, 1142 * (GetNode(i+1,j+1,k,n,iside)+1) ;
1137 faces[nface][3] = GetEdgeVisibility(i 1143 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1)
1138 * (GetNode(i+1,j ,k, 1144 * (GetNode(i+1,j ,k,n,iside)+1) ;
1139 } 1145 }
1140 } 1146 }
1141 } 1147 }
1142 } 1148 }
1143 1149