| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4TwistTrapAlphaSide implementation
27 //
28 // Author: 18/03/2005 - O.Link (Oliver.Link@ce
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,
43 G4double pDz,
44 G4double pTheta,
45 G4double pPhi,
46 G4double pDy1,
47 G4double pDx1,
48 G4double pDx2,
49 G4double pDy2,
50 G4double pDx3,
51 G4double pDx4,
52 G4double pAlph,
53 G4double AngleSide
54 )
55 : G4VTwistSurface(name)
56 {
57 fAxis[0] = kYAxis; // in local coordinate
58 fAxis[1] = kZAxis;
59 fAxisMin[0] = -kInfinity ; // Y Axis bounda
60 fAxisMax[0] = kInfinity ; // depends on
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 de
93
94 fdeltaX = 2 * fDz * std::tan(fTheta) * std::
95 // dx in surface equation
96 fdeltaY = 2 * fDz * std::tan(fTheta) * std::
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( __
113 : G4VTwistSurface(a)
114 {
115 }
116
117
118 //============================================
119 //* destructor -------------------------------
120
121 G4TwistTrapAlphaSide::~G4TwistTrapAlphaSide()
122
123
124 //============================================
125 //* GetNormal --------------------------------
126
127 G4ThreeVector
128 G4TwistTrapAlphaSide::GetNormal(const G4ThreeV
129 G4bool i
130 {
131 // GetNormal returns a normal vector at a s
132 // to surface) point at tmpxx.
133 // If isGlobal=true, it returns the normal
134 //
135
136 G4ThreeVector xx;
137 if (isGlobal)
138 {
139 xx = ComputeLocalPoint(tmpxx);
140 if ((xx - fCurrentNormal.p).mag() < 0.5
141 {
142 return ComputeGlobalDirection(fCurren
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
158
159 G4ThreeVector normal = NormAng(phi,u) ; /
160
161 #ifdef G4TWISTDEBUG
162 G4cout << "normal vector = " << normal <<
163 G4cout << "phi = " << phi << " , u = " << u
164 #endif
165
166 if (isGlobal)
167 {
168 fCurrentNormal.normal = ComputeGlobalDir
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
183 const
184
185
186
187
188
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 solut
197
198 fCurStatWithV.ResetfDone(validate, &gp, &gv)
199
200 if (fCurStatWithV.IsDone())
201 {
202 for (G4int i=0; i<fCurStatWithV.GetNXX();
203 {
204 gxx[i] = fCurStatWithV.GetXX(i);
205 distance[i] = fCurStatWithV.GetDistance(
206 areacode[i] = fCurStatWithV.GetAreacode(
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, kInfini
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 << G4e
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 intersecti
243
244 G4double L = 2*fDz ;
245
246 G4double phixz = fPhiTwist * ( p.x() * v.z()
247 G4double phiyz = fPhiTwist * ( p.y() * v.z()
248
249
250 // special case vz = 0
251
252 if ( v.z() == 0. )
253 {
254 if ( std::fabs(p.z()) <= L ) // intersec
255 {
256 phi = p.z() * fPhiTwist / L ; // phi is
257 u = (fDy1*(4*(-(fdeltaY*phi*v.x()) + fPh
258 + fdeltaX*phi*v.y() - fPhi
259 + ((fDx3plus1 + fDx4plus2)*fP
260 + 2*(fDx3minus1 + fDx4minus2
261 *(v.y()*std::cos(phi) - v
262 /(fPhiTwist*(4*fDy1* v.x() - (fa1md1
263 *std::cos(phi) + fPhiTwist
264 + 4*fDy1*(fTAlph*v.x()
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 t
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]
280 areacode[
281 0, valida
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)*fPh
294 c[6] = -57600*
295 fDy1*(4*fDy1*(phiyz + 2*fDz*v.x() + fTAl
296 - 2*fDy1*(2*fdeltaX + fDx3minus1 + f
297 - 2*fdeltaY*fTAlph)*v.z()
298 + fa1md1*(phixz - 2*fDz*v.y() + fdel
299 c[5] = 4800*
300 fDy1*(fa1md1*(-5*phiyz - 24*fDz*v.x() +
301 fDy1*(20*phixz - 4*(5*fTAlph*phiyz
302 + 24*fDz*v.y()) + (48*fdelt
303 *fPhiTwist
304 c[4] = 4800*
305 fDy1*(fa1md1*(phixz - 10*fDz*v.y() + 5*f
306 + 2*fDy1*(2*phiyz + 20*fDz*v.x()
307 + (-10*fdeltaX + fDx3minus1
308 + 2*fTAlph*(phixz - 10*fDz*
309 c[3] = -96*
310 fDy1*(-(fa1md1*(phiyz + 100*fDz*v.x() -
311 + fDy1*(4*phixz - 400*fDz*v.y()
312 + (200*fdeltaY - (fDx3plus1
313 - 4*fTAlph*(phiyz + 100*fDz*
314 c[2] = 32*
315 fDy1*(4*fDy1*(7*fTAlph*phixz + 7*phiyz -
316 + 6*fDy1*(2*fdeltaX+fDx3minus1+fDx4m
317 + fa1md1*(7*phixz + 6*fDz*v.y() - 3*
318 c[1] = -8*
319 fDy1*(fa1md1*(-9*phiyz - 56*fDz*v.x() +
320 + 4*fDy1*(9*phixz - 9*fTAlph*phiyz -
321 - 56*fDz*v.y() + 28*(fdeltaY
322 c[0] = 72*
323 fDy1*(fa1md1*(2*fDz*v.y() - fdeltaY*v.z(
324 + fDy1*(-8*fDz*v.x() + 8*fDz*fTAlph*
325 + 4*fdeltaX*v.z() - 4*fdeltaY*
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
342 {
343 if ( si[i]==0.0 ) // only real solution
344 {
345 #ifdef G4TWISTDEBUG
346 G4cout << "Solution " << i << " : " <<
347 #endif
348 phi = std::fmod(srd[i] , pihalf) ;
349 u = (fDy1*(4*(phiyz + 2*fDz*phi*v.y(
350 - ((fDx3plus1 + fDx4plus2
351 + 2*(fDx3minus1 + fDx4m
352 /(fPhiTwist*v.z()*(4*fDy1*std::co
353 + (fa1md1 + 4*fDy
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
361
362 #ifdef G4TWISTDEBUG
363 G4cout << "solution " << i << " = " <<
364 #endif
365 } // end if real solution
366 } // end loop i
367 } // end general case
368
369 nxx = (G4int)xbuf.size() ; // save the numb
370
371 G4ThreeVector xxonsurface ; // point on
372 G4ThreeVector surfacenormal ; // normal ve
373 G4double deltaX; // distance between inters
374 G4double theta; // angle between track and
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[
383 << ", uR = " << xbuf[k].u << G4endl
384 #endif
385
386 phi = k.phi ; // get the stored values fo
387 u = k.u ;
388
389 IsConverged = false ; // no convergence
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, xxo
397 deltaX = ( tmpxx - xxonsurface ).mag() ;
398 theta = std::fabs(std::acos(v*surfacenor
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 << ", distanc
411 << ", " << deltaX << G4endl ;
412 G4cout << "X = " << tmpxx << G4endl ;
413 #endif
414
415 GetPhiUAtX(tmpxx, phi, u) ;
416 // the new point xx is accepted and ph
417
418 #ifdef G4TWISTDEBUG
419 G4cout << "approximated phi = " << phi <
420 #endif
421
422 if ( deltaX <= factor*ctol ) { IsConverg
423
424 } // end iterative loop (i)
425
426 if ( std::fabs(tmpdist)<ctol ) { tmpdist =
427
428 #ifdef G4TWISTDEBUG
429 G4cout << "refined solution " << phi << "
430 G4cout << "distance = " << tmpdist << G4en
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 = tru
452 }
453 }
454 else // kDontValidate
455 {
456 G4Exception("G4TwistTrapAlphaSide::Dis
457 "GeomSolids0001", FatalExc
458 "Feature NOT implemented !
459 }
460 }
461 else
462 {
463 tmpdist = kInfinity; // no convergen
464 tmpisvalid = false ; // solution is
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 (vari
475
476 std::sort(xbuf.begin() , xbuf.end(), Distanc
477
478 #ifdef G4TWISTDEBUG
479 G4cout << G4endl << "list xbuf after sorting
480 G4cout << G4endl << G4endl ;
481 #endif
482
483 // erase identical intersection (within kCar
484 //
485 xbuf.erase( std::unique(xbuf.begin(), xbuf.e
486 xbuf.end() );
487
488
489 // add guesses
490 //
491 auto nxxtmp = (G4int)xbuf.size() ;
492
493 if ( nxxtmp<2 || IsParallel ) // positive
494 {
495
496 #ifdef G4TWISTDEBUG
497 G4cout << "add guess at +z/2 .. " << G4end
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
510
511 #ifdef G4TWISTDEBUG
512 G4cout << "add guess at -z/2 .. " << G4end
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
525
526 for ( std::size_t k = nxxtmp ; k<xbuf.size
527 {
528
529 #ifdef G4TWISTDEBUG
530 G4cout << "Solution " << k << " : "
531 << "reconstructed phiR = " << xbu
532 << ", uR = " << xbuf[k].u << G4en
533 #endif
534
535 phi = xbuf[k].phi ; // get the stored v
536 u = xbuf[k].u ;
537
538 IsConverged = false ; // no convergenc
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, x
545 deltaX = ( tmpxx - xxonsurface ).mag()
546 theta = std::fabs(std::acos(v*surfacen
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 << ", dista
558 << ", " << deltaX << G4endl
559 << "X = " << tmpxx << G4endl ;
560 #endif
561
562 GetPhiUAtX(tmpxx, phi, u) ;
563 // the new point xx is accepted and
564
565 #ifdef G4TWISTDEBUG
566 G4cout << "approximated phi = " << phi
567 #endif
568
569 if ( deltaX <= factor*ctol ) { IsConve
570
571 } // end iterative loop (i)
572
573 if ( std::fabs(tmpdist)<ctol ) { tmpdis
574
575 #ifdef G4TWISTDEBUG
576 G4cout << "refined solution " << phi <<
577 G4cout << "distance = " << tmpdist << G4
578 G4cout << "local X = " << tmpxx << G4end
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 =
591 }
592 }
593 else if (validate == kValidateWithoutT
594 {
595 tmpareacode = GetAreaCode(tmpxx, fal
596 if (IsInside(tmpareacode))
597 {
598 if (tmpdist >= 0) { tmpisvalid =
599 }
600 }
601 else // kDontValidate
602 {
603 G4Exception("G4TwistedBoxSide::Dista
604 "GeomSolids0001", FatalE
605 "Feature NOT implemented
606 }
607 }
608 else
609 {
610 tmpdist = kInfinity; // no converg
611 tmpisvalid = false ; // solution i
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(), Distanc
626
627 // erase identical intersection (within kCar
628 xbuf.erase( std::unique(xbuf.begin(), xbuf.e
629 xbuf.end() );
630
631 #ifdef G4TWISTDEBUG
632 G4cout << G4endl << "list xbuf after sorting
633 G4cout << G4endl << G4endl ;
634 #endif
635
636 nxx = (G4int)xbuf.size() ; // determine nu
637
638 for ( G4int i = 0 ; i<(G4int)xbuf.size() ; +
639 {
640 distance[i] = xbuf[i].distance;
641 gxx[i] = ComputeGlobalPoint(xbuf[i].x
642 areacode[i] = xbuf[i].areacode ;
643 isvalid[i] = xbuf[i].isvalid ;
644
645 fCurStatWithV.SetCurrentStatus(i, gxx[i],
646 isvalid[i
647 #ifdef G4TWISTDEBUG
648 G4cout << "element Nr. " << i
649 << ", local Intersection = " << xbu
650 << ", distance = " << xbuf[i].dista
651 << ", u = " << xbuf[i].u
652 << ", phi = " << xbuf[i].phi
653 << ", isvalid = " << xbuf[i].isvali
654 << G4endl ;
655 #endif
656
657 } // end for( i ) loop
658
659 #ifdef G4TWISTDEBUG
660 G4cout << "G4TwistTrapAlphaSide finished " <
661 G4cout << nxx << " possible physical solutio
662 for ( G4int k= 0 ; k< nxx ; k++ )
663 {
664 G4cout << "global intersection Point found
665 G4cout << "distance = " << distance[k] <<
666 G4cout << "isvalid = " << isvalid[k] << G4
667 }
668 #endif
669
670 return nxx ;
671 }
672
673
674 //============================================
675 //* DistanceToSurface ------------------------
676
677 G4int
678 G4TwistTrapAlphaSide::DistanceToSurface(const
679
680
681
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, kInf
704 }
705 }
706
707 G4ThreeVector p = ComputeLocalPoint(gp);
708 G4ThreeVector xx; // intersection point
709 G4ThreeVector xxonsurface ; // interpolated
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,xxonsurfa
725 deltaX = ( xx - xxonsurface ).mag() ;
726
727 #ifdef G4TWISTDEBUG
728 G4cout << "i = " << i << ", distance = "
729 << ", " << deltaX << G4endl
730 << "X = " << xx << G4endl ;
731 #endif
732
733 // the new point xx is accepted and phi/p
734 //
735 GetPhiUAtX(xx, phiR, uR) ;
736
737 if ( deltaX <= ctol ) { break ; }
738 }
739
740 // check validity of solution ( valid phi,p
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] =
752
753 // end of validity
754
755 #ifdef G4TWISTDEBUG
756 G4cout << "refined solution " << phiR << "
757 G4cout << "distance = " << distance[0] << G
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: " << g
766 G4cout << "distance = " << distance[0] << G
767 #endif
768
769 fCurStat.SetCurrentStatus(0, gxx[0], distan
770 isvalid, 1, kDontV
771 return 1;
772 }
773
774
775 //============================================
776 //* GetAreaCode ------------------------------
777
778 G4int
779 G4TwistTrapAlphaSide::GetAreaCode(const G4Thre
780 {
781 // We must use the function in local coordi
782 // See the description of DistanceToSurface
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 << G
795 G4cout << "GetAreaCode: yprime = " << yprim
796 G4cout << "Intervall is " << fYAxisMin << "
797 #endif
798
799 G4int areacode = sInside;
800
801 if (fAxis[0] == kYAxis && fAxis[1] == kZAxi
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 | sA
814 if (yprime <= fYAxisMin - ctol) {
815
816 }
817 else if (yprime > fYAxisMax - ctol)
818 {
819 areacode |= (sAxis0 & (sAxisY | sA
820 if (yprime >= fYAxisMax + ctol) {
821 }
822
823 // test boundary of z-axis
824
825 if (xx.z() < fAxisMin[zaxis] + ctol)
826 {
827 areacode |= (sAxis1 & (sAxisZ | sA
828
829 if ((areacode & sBoundary) != 0)
830 { areacode |= sCorner; }
831
832 else
833 { areacode |= sBoundary; }
834 if (xx.z() <= fAxisMin[zaxis] - ct
835 }
836 else if (xx.z() > fAxisMax[zaxis] - c
837 {
838 areacode |= (sAxis1 & (sAxisZ | sA
839
840 if ((areacode & sBoundary) != 0)
841 { areacode |= sCorner; }
842 else
843 { areacode |= sBoundary; }
844 if (xx.z() >= fAxisMax[zaxis] + ct
845 }
846
847 // if isoutside = true, clear inside
848 // if not on boundary, add axis infor
849
850 if (isoutside)
851 {
852 G4int tmpareacode = areacode & (~s
853 areacode = tmpareacode;
854 }
855 else if ((areacode & sBoundary) != sB
856 {
857 areacode |= (sAxis0 & sAxisY) | (s
858 }
859
860 }
861 else
862 {
863 // boundary of y-axis
864
865 if (yprime < fYAxisMin )
866 {
867 areacode |= (sAxis0 & (sAxisY | sA
868 }
869 else if (yprime > fYAxisMax)
870 {
871 areacode |= (sAxis0 & (sAxisY | sA
872 }
873
874 // boundary of z-axis
875
876 if (xx.z() < fAxisMin[zaxis])
877 {
878 areacode |= (sAxis1 & (sAxisZ | sA
879 if ((areacode & sBoundary) != 0)
880 { areacode |= sCorner; }
881 else
882 { areacode |= sBoundary; }
883 }
884 else if (xx.z() > fAxisMax[zaxis])
885 {
886 areacode |= (sAxis1 & (sAxisZ | sA
887 if ((areacode & sBoundary) != 0)
888 { areacode |= sCorner; }
889 else
890 { areacode |= sBoundary; }
891 }
892
893 if ((areacode & sBoundary) != sBounda
894 {
895 areacode |= (sAxis0 & sAxisY) | (s
896 }
897 }
898 return areacode;
899 }
900 else
901 {
902 G4Exception("G4TwistTrapAlphaSide::GetAr
903 "GeomSolids0001", FatalExcep
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
925 - fDy1*std::sin(fPhiTwist/2.);
926 y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/
927 + (-fDx1 + fDy1*fTAlph)*std::sin(fPhiTwi
928 z = -fDz ;
929
930 // G4cout << "SetCorners: " << x << ",
931
932 SetCorner(sC0Min1Min, x, y, z);
933
934 // corner of Axis0max and Axis1min
935 //
936 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std
937 + fDy1*std::sin(fPhiTwist/2.);
938 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/
939 - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwis
940 z = -fDz ;
941
942 // G4cout << "SetCorners: " << x << ",
943
944 SetCorner(sC0Max1Min, x, y, z);
945
946 // corner of Axis0max and Axis1max
947 //
948 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std:
949 - fDy2*std::sin(fPhiTwist/2.);
950 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2
951 + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwis
952 z = fDz ;
953
954 // G4cout << "SetCorners: " << x << ",
955
956 SetCorner(sC0Max1Max, x, y, z);
957
958 // corner of Axis0min and Axis1max
959 x = fdeltaX/2. + (fDx3 - fDy2*fTAlph)*std:
960 + fDy2*std::sin(fPhiTwist/2.) ;
961 y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2
962 + (fDx3 - fDy2*fTAlph)*std::sin(fPhiTwis
963 z = fDz ;
964
965 // G4cout << "SetCorners: " << x << ",
966
967 SetCorner(sC0Min1Max, x, y, z);
968
969 }
970 else
971 {
972 G4Exception("G4TwistTrapAlphaSide::SetCorn
973 "GeomSolids0001", FatalExcepti
974 "Method NOT implemented !");
975 }
976 }
977
978 //============================================
979 //* SetBoundaries() --------------------------
980
981 void G4TwistTrapAlphaSide::SetBoundaries()
982 {
983 // Set direction-unit vector of boundary-li
984 //
985
986 G4ThreeVector direction;
987
988 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis
989 {
990 // sAxis0 & sAxisMin
991 direction = GetCorner(sC0Min1Max) - GetCor
992 direction = direction.unit();
993 SetBoundary(sAxis0 & (sAxisY | sAxisMin),
994 GetCorner(sC0Min1Min), sAxisZ)
995
996 // sAxis0 & sAxisMax
997 direction = GetCorner(sC0Max1Max) - GetCor
998 direction = direction.unit();
999 SetBoundary(sAxis0 & (sAxisY | sAxisMax),
1000 GetCorner(sC0Max1Min), sAxisZ
1001
1002 // sAxis1 & sAxisMin
1003 direction = GetCorner(sC0Max1Min) - GetCo
1004 direction = direction.unit();
1005 SetBoundary(sAxis1 & (sAxisZ | sAxisMin),
1006 GetCorner(sC0Min1Min), sAxisY
1007
1008 // sAxis1 & sAxisMax
1009 direction = GetCorner(sC0Max1Max) - GetCo
1010 direction = direction.unit();
1011 SetBoundary(sAxis1 & (sAxisZ | sAxisMax),
1012 GetCorner(sC0Min1Max), sAxisY
1013
1014 }
1015 else
1016 {
1017 G4Exception("G4TwistTrapAlphaSide::SetCor
1018 "GeomSolids0001", FatalExcept
1019 "Feature NOT implemented !");
1020 }
1021 }
1022
1023 //===========================================
1024 //* GetPhiUAtX ------------------------------
1025
1026 void
1027 G4TwistTrapAlphaSide::GetPhiUAtX( const G4Thr
1028 {
1029 // find closest point XX on surface for a g
1030 // X0 is a point on the surface, d is the
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 -
1037 - 4*(fDx3plus1 + fDx4plus2)*fDy1
1038 - 2*(2*fDx1*fDx1 - 2*fDx2*fDx2 +
1039 + 4*(fDx3minus1 + fDx4minus2)
1040 - 4*(fa1md1*(fdeltaX*phi - fPhiT
1041 + 4*fDy1*(fdeltaY*phi + fdeltaX*
1042 - fPhiTwist*(fTAlph*p.x(
1043 - 4*(fa1md1*fdeltaY*phi - 4*fdel
1044 + 4*fdeltaY*fDy1*fTAlph*phi +
1045 - fPhiTwist*(fa1md1 + 4*fDy1*
1046 /(fDy1* fPhiTwist*((std::fabs(((fa1md1
1047 /fDy1 - 4*std::sin(
1048 *(std::fabs(((fa1md1
1049 /fDy1 - 4*std::sin(
1050 + (std::fabs(4*std::cos(phi)
1051 + ((fa1md1 + 4*fDy1*fTAlph)*s
1052 * (std::fabs(4*std::cos(phi)
1053 + ((fa1md1 + 4*fDy1*fTAlph)*s
1054 }
1055
1056 //===========================================
1057 //* ProjectPoint ----------------------------
1058
1059 G4ThreeVector
1060 G4TwistTrapAlphaSide::ProjectPoint(const G4Th
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
1079
1080 G4ThreeVector xx = SurfacePoint(phi,u) ;
1081 // transform back to cartesian coordinate
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, G4i
1098 G4int faces[
1099 {
1100
1101 G4double phi ;
1102 G4double b ;
1103
1104 G4double z, u ; // the two parameters f
1105 G4ThreeVector p ; // a point on the surfac
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) ; // surf
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 w
1129 {
1130 nface = GetFace(i,j,k,n,iside) ;
1131 faces[nface][0] = GetEdgeVisibility(i
1132 * (GetNode(i ,j ,k,
1133 faces[nface][1] = GetEdgeVisibility(i
1134 * (GetNode(i ,j+1,k,
1135 faces[nface][2] = GetEdgeVisibility(i
1136 * (GetNode(i+1,j+1,k,
1137 faces[nface][3] = GetEdgeVisibility(i
1138 * (GetNode(i+1,j ,k,
1139 }
1140 }
1141 }
1142 }
1143