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