source: trunk/source/geometry/solids/specific/src/G4TwistTrapParallelSide.cc@ 1316

Last change on this file since 1316 was 1228, checked in by garnier, 16 years ago

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