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