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