Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TwistTrapAlphaSide.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4TwistTrapAlphaSide.cc 67011 2013-01-29 16:17:41Z 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 
51 G4TwistTrapAlphaSide::
52 G4TwistTrapAlphaSide(const G4String &name,
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 
124 G4TwistTrapAlphaSide::G4TwistTrapAlphaSide( __void__& a )
125  : G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
126  fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
127  fAngleSide(0.), 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 
137 G4TwistTrapAlphaSide::~G4TwistTrapAlphaSide()
138 {
139 }
140 
141 
142 //=====================================================================
143 //* GetNormal ---------------------------------------------------------
144 
145 G4ThreeVector
146 G4TwistTrapAlphaSide::GetNormal(const G4ThreeVector &tmpxx,
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 
154  G4ThreeVector xx;
155  if (isGlobal)
156  {
157  xx = ComputeLocalPoint(tmpxx);
158  if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
159  {
160  return ComputeGlobalDirection(fCurrentNormal.normal);
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  {
186  fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
187  }
188  else
189  {
190  fCurrentNormal.normal = normal.unit();
191  }
192 
193  return fCurrentNormal.normal;
194 }
195 
196 //=====================================================================
197 //* DistanceToSurface -------------------------------------------------
198 
199 G4int
200 G4TwistTrapAlphaSide::DistanceToSurface(const G4ThreeVector &gp,
201  const G4ThreeVector &gv,
202  G4ThreeVector gxx[],
203  G4double distance[],
204  G4int areacode[],
205  G4bool isvalid[],
206  EValidate validate)
207 {
208  static const G4double ctol = 0.5 * kCarTolerance;
209  static const G4double pihalf = pi/2 ;
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 (register int 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 (register int 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 
240  G4ThreeVector p = ComputeLocalPoint(gp);
241  G4ThreeVector v = ComputeLocalDirection(gv);
242 
243 #ifdef G4TWISTDEBUG
244  G4cout << "Local point p = " << p << G4endl ;
245  G4cout << "Local direction v = " << v << G4endl ;
246 #endif
247 
248  G4double phi,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())
277  + ((fDx3plus1 + fDx4plus2)*fPhiTwist
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;
294  gxx[0].set(kInfinity,kInfinity,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())
334  + 6*fDy1*(2*fdeltaX+fDx3minus1+fDx4minus2-2*fdeltaY*fTAlph)*v.z()
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 (register int 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())
368  - ((fDx3plus1 + fDx4plus2)*fPhiTwist
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 ( register 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 ( register int 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 ( register 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 ( register int 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 ( register 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
696 G4TwistTrapAlphaSide::DistanceToSurface(const G4ThreeVector &gp,
697  G4ThreeVector gxx[],
698  G4double distance[],
699  G4int areacode[])
700 {
701  static const G4double ctol = 0.5 * kCarTolerance;
702 
703  fCurStat.ResetfDone(kDontValidate, &gp);
704 
705  if (fCurStat.IsDone())
706  {
707  for (register int 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 (register int i=0; i<G4VSURFACENXX; i++)
718  {
719  distance[i] = kInfinity;
720  areacode[i] = sOutside;
721  gxx[i].set(kInfinity, kInfinity, kInfinity);
722  }
723  }
724 
725  G4ThreeVector p = ComputeLocalPoint(gp);
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 ( register int 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
797 G4TwistTrapAlphaSide::GetAreaCode(const G4ThreeVector &xx, G4bool withTol)
798 {
799  // We must use the function in local coordinate system.
800  // See the description of DistanceToSurface(p,v).
801 
802  static 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 
930 void G4TwistTrapAlphaSide::SetCorners()
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 
999 void G4TwistTrapAlphaSide::SetBoundaries()
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,
1012  GetCorner(sC0Min1Min), sAxisZ) ;
1013 
1014  // sAxis0 & sAxisMax
1015  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
1016  direction = direction.unit();
1017  SetBoundary(sAxis0 & (sAxisY | sAxisMax), direction,
1018  GetCorner(sC0Max1Min), sAxisZ);
1019 
1020  // sAxis1 & sAxisMin
1021  direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
1022  direction = direction.unit();
1023  SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
1024  GetCorner(sC0Min1Min), sAxisY);
1025 
1026  // sAxis1 & sAxisMax
1027  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
1028  direction = direction.unit();
1029  SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
1030  GetCorner(sC0Min1Max), sAxisY);
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
1045 G4TwistTrapAlphaSide::GetPhiUAtX( G4ThreeVector p, G4double &phi, G4double &u )
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 
1077 G4ThreeVector
1078 G4TwistTrapAlphaSide::ProjectPoint(const G4ThreeVector &p, G4bool isglobal)
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
1115 G4TwistTrapAlphaSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
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 ( register int 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 ( register int 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 }
Generated on Sat May 25 2013 14:33:16 for Geant4 by   doxygen 1.8.4