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