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