source: trunk/source/geometry/solids/specific/src/G4TwistTrapAlphaSide.cc @ 1058

Last change on this file since 1058 was 1058, checked in by garnier, 15 years ago

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