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