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