Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
44 #include "G4TwistTrapParallelSide.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4JTPolynomialSolver.hh"
47 
48 //=====================================================================
49 //* constructors ------------------------------------------------------
50 
51 G4TwistTrapParallelSide::G4TwistTrapParallelSide(const G4String &name,
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 
122 G4TwistTrapParallelSide::G4TwistTrapParallelSide( __void__& a )
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 
135 G4TwistTrapParallelSide::~G4TwistTrapParallelSide()
136 {
137 }
138 
139 //=====================================================================
140 //* GetNormal ---------------------------------------------------------
141 
142 G4ThreeVector G4TwistTrapParallelSide::GetNormal(const G4ThreeVector &tmpxx,
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 
150  G4ThreeVector xx;
151  if (isGlobal) {
152  xx = ComputeLocalPoint(tmpxx);
153  if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
154  return ComputeGlobalDirection(fCurrentNormal.normal);
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) {
178  fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
179  } else {
180  fCurrentNormal.normal = normal.unit();
181  }
182  return fCurrentNormal.normal;
183 }
184 
185 //=====================================================================
186 //* DistanceToSurface -------------------------------------------------
187 
188 G4int G4TwistTrapParallelSide::DistanceToSurface(const G4ThreeVector &gp,
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 ctol = 0.5 * kCarTolerance;
198  static const G4double pihalf = pi/2 ;
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 
228  G4ThreeVector p = ComputeLocalPoint(gp);
229  G4ThreeVector v = ComputeLocalDirection(gv);
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 
655 G4int G4TwistTrapParallelSide::DistanceToSurface(const G4ThreeVector &gp,
656  G4ThreeVector gxx[],
657  G4double distance[],
658  G4int areacode[])
659 {
660  // to do
661 
662  static const G4double ctol = 0.5 * kCarTolerance;
663 
664  fCurStat.ResetfDone(kDontValidate, &gp);
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 
684  G4ThreeVector p = ComputeLocalPoint(gp);
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  static 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,
927  GetCorner(sC0Min1Min), sAxisZ) ;
928 
929  // sAxis0 & sAxisMax
930  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
931  direction = direction.unit();
932  SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction,
933  GetCorner(sC0Max1Min), sAxisZ);
934 
935  // sAxis1 & sAxisMin
936  direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
937  direction = direction.unit();
938  SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
939  GetCorner(sC0Min1Min), sAxisX);
940 
941  // sAxis1 & sAxisMax
942  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
943  direction = direction.unit();
944  SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
945  GetCorner(sC0Min1Max), sAxisX);
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 }
Generated on Sat May 25 2013 14:33:16 for Geant4 by   doxygen 1.8.4