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