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