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