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