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