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