source: trunk/source/geometry/solids/specific/src/G4TwistBoxSide.cc@ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 15 years ago

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