source: trunk/source/geometry/solids/specific/src/G4TwistedTubs.cc @ 1294

Last change on this file since 1294 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
File size: 39.3 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: G4TwistedTubs.cc,v 1.24 2007/05/18 07:39:56 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-03 $
29//
30//
31// --------------------------------------------------------------------
32// GEANT 4 class source file
33//
34//
35// G4TwistTubsSide.cc
36//
37// Author:
38//   01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
39//
40// History:
41//   13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
42//                 from original version in Jupiter-2.5.02 application.
43// --------------------------------------------------------------------
44
45#include "G4TwistedTubs.hh"
46
47#include "G4GeometryTolerance.hh"
48#include "G4VoxelLimits.hh"
49#include "G4AffineTransform.hh"
50#include "G4SolidExtentList.hh"
51#include "G4ClippablePolygon.hh"
52#include "G4VPVParameterisation.hh"
53#include "meshdefs.hh"
54
55#include "G4VGraphicsScene.hh"
56#include "G4Polyhedron.hh"
57#include "G4VisExtent.hh"
58#include "G4NURBS.hh"
59#include "G4NURBStube.hh"
60#include "G4NURBScylinder.hh"
61#include "G4NURBStubesector.hh"
62
63#include "Randomize.hh"
64
65//=====================================================================
66//* constructors ------------------------------------------------------
67
68G4TwistedTubs::G4TwistedTubs(const G4String &pname,
69                                   G4double  twistedangle,
70                                   G4double  endinnerrad,
71                                   G4double  endouterrad,
72                                   G4double  halfzlen,
73                                   G4double  dphi)
74   : G4VSolid(pname), fDPhi(dphi), 
75     fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
76     fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
77     fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
78{
79   if (endinnerrad < DBL_MIN)
80   {
81      G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
82                  FatalException, "Invalid end-inner-radius!");
83   }
84           
85   G4double sinhalftwist = std::sin(0.5 * twistedangle);
86
87   G4double endinnerradX = endinnerrad * sinhalftwist;
88   G4double innerrad     = std::sqrt( endinnerrad * endinnerrad
89                                 - endinnerradX * endinnerradX );
90
91   G4double endouterradX = endouterrad * sinhalftwist;
92   G4double outerrad     = std::sqrt( endouterrad * endouterrad
93                                 - endouterradX * endouterradX );
94   
95   // temporary treatment!!
96   SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
97   CreateSurfaces();
98}
99
100G4TwistedTubs::G4TwistedTubs(const G4String &pname,     
101                                   G4double  twistedangle,   
102                                   G4double  endinnerrad, 
103                                   G4double  endouterrad, 
104                                   G4double  halfzlen,
105                                   G4int     nseg,
106                                   G4double  totphi)
107   : G4VSolid(pname),
108     fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
109     fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
110     fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
111{
112
113   if (!nseg)
114   {
115     G4cerr << "ERROR - G4TwistedTubs::G4TwistedTubs()" << G4endl
116            << "        Invalid nseg. nseg = " << nseg << G4endl;
117   }
118   if (totphi == DBL_MIN || endinnerrad < DBL_MIN)
119   {
120      G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
121                  FatalException, "Invalid total-phi or end-inner-radius!");
122   }   
123         
124   G4double sinhalftwist = std::sin(0.5 * twistedangle);
125
126   G4double endinnerradX = endinnerrad * sinhalftwist;
127   G4double innerrad     = std::sqrt( endinnerrad * endinnerrad
128                                 - endinnerradX * endinnerradX );
129
130   G4double endouterradX = endouterrad * sinhalftwist;
131   G4double outerrad     = std::sqrt( endouterrad * endouterrad
132                                 - endouterradX * endouterradX );
133   
134   // temporary treatment!!
135   fDPhi = totphi / nseg;
136   SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
137   CreateSurfaces();
138}
139
140G4TwistedTubs::G4TwistedTubs(const G4String &pname,
141                                   G4double  twistedangle,
142                                   G4double  innerrad,
143                                   G4double  outerrad,
144                                   G4double  negativeEndz,
145                                   G4double  positiveEndz,
146                                   G4double  dphi)
147   : G4VSolid(pname), fDPhi(dphi),
148     fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
149     fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
150     fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
151{
152   if (innerrad < DBL_MIN)
153   {
154      G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
155                  FatalException, "Invalid end-inner-radius!");
156   }
157                 
158   SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
159   CreateSurfaces();
160}
161
162G4TwistedTubs::G4TwistedTubs(const G4String &pname,     
163                                   G4double  twistedangle,   
164                                   G4double  innerrad, 
165                                   G4double  outerrad, 
166                                   G4double  negativeEndz,
167                                   G4double  positiveEndz,
168                                   G4int     nseg,
169                                   G4double  totphi)
170   : G4VSolid(pname),
171     fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
172     fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
173     fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
174{
175   if (!nseg)
176   {
177     G4cerr << "ERROR - G4TwistedTubs::G4TwistedTubs()" << G4endl
178            << "        Invalid nseg. nseg = " << nseg << G4endl;
179   }
180   if (totphi == DBL_MIN || innerrad < DBL_MIN)
181   {
182      G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
183                  FatalException, "Invalid total-phi or end-inner-radius!");
184   }
185           
186   fDPhi = totphi / nseg;
187   SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
188   CreateSurfaces();
189}
190
191//=====================================================================
192//* Fake default constructor ------------------------------------------
193
194G4TwistedTubs::G4TwistedTubs( __void__& a )
195  : G4VSolid(a), fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
196    fFormerTwisted(0), fInnerHype(0), fOuterHype(0), fCubicVolume(0.),
197    fSurfaceArea(0.), fpPolyhedron(0)
198{
199}
200
201//=====================================================================
202//* destructor --------------------------------------------------------
203
204G4TwistedTubs::~G4TwistedTubs()
205{
206   if (fLowerEndcap)   { delete fLowerEndcap;   }
207   if (fUpperEndcap)   { delete fUpperEndcap;   }
208   if (fLatterTwisted) { delete fLatterTwisted; }
209   if (fFormerTwisted) { delete fFormerTwisted; }
210   if (fInnerHype)     { delete fInnerHype;     }
211   if (fOuterHype)     { delete fOuterHype;     }
212   if (fpPolyhedron)   { delete fpPolyhedron;   }
213}
214
215//=====================================================================
216//* ComputeDimensions -------------------------------------------------
217
218void G4TwistedTubs::ComputeDimensions(G4VPVParameterisation* /* p */ ,
219                                      const G4int            /* n  */ ,
220                                      const G4VPhysicalVolume* /* pRep */ )
221{
222  G4Exception("G4TwistedTubs::ComputeDimensions()",
223              "NotSupported", FatalException,
224              "G4TwistedTubs does not support Parameterisation.");
225}
226
227
228//=====================================================================
229//* CalculateExtent ---------------------------------------------------
230
231G4bool G4TwistedTubs::CalculateExtent( const EAxis              axis,
232                                       const G4VoxelLimits     &voxelLimit,
233                                       const G4AffineTransform &transform,
234                                             G4double          &min,
235                                             G4double          &max ) const
236{
237
238  G4SolidExtentList  extentList( axis, voxelLimit );
239  G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
240                             fEndOuterRadius[0] : fEndOuterRadius[1]);
241  G4double maxEndInnerRad = (fEndInnerRadius[0] > fEndInnerRadius[1] ?
242                             fEndInnerRadius[0] : fEndInnerRadius[1]);
243  G4double maxphi         = (std::fabs(fEndPhi[0]) > std::fabs(fEndPhi[1]) ?
244                             std::fabs(fEndPhi[0]) : std::fabs(fEndPhi[1]));
245   
246  //
247  // Choose phi size of our segment(s) based on constants as
248  // defined in meshdefs.hh
249  //
250  // G4int numPhi = kMaxMeshSections;
251  G4double sigPhi = 2*maxphi + fDPhi;
252  G4double rFudge = 1.0/std::cos(0.5*sigPhi);
253  G4double fudgeEndOuterRad = rFudge * maxEndOuterRad;
254 
255  //
256  // We work around in phi building polygons along the way.
257  // As a reasonable compromise between accuracy and
258  // complexity (=cpu time), the following facets are chosen:
259  //
260  //   1. If fOuterRadius/maxEndOuterRad > 0.95, approximate
261  //      the outer surface as a cylinder, and use one
262  //      rectangular polygon (0-1) to build its mesh.
263  //
264  //      Otherwise, use two trapazoidal polygons that
265  //      meet at z = 0 (0-4-1)
266  //
267  //   2. If there is no inner surface, then use one
268  //      polygon for each entire endcap.  (0) and (1)
269  //
270  //      Otherwise, use a trapazoidal polygon for each
271  //      phi segment of each endcap.    (0-2) and (1-3)
272  //
273  //   3. For the inner surface, if fInnerRadius/maxEndInnerRad > 0.95,
274  //      approximate the inner surface as a cylinder of
275  //      radius fInnerRadius and use one rectangular polygon
276  //      to build each phi segment of its mesh.   (2-3)
277  //
278  //      Otherwise, use one rectangular polygon centered
279  //      at z = 0 (5-6) and two connecting trapazoidal polygons
280  //      for each phi segment (2-5) and (3-6).
281  //
282
283  G4bool splitOuter = (fOuterRadius/maxEndOuterRad < 0.95);
284  G4bool splitInner = (fInnerRadius/maxEndInnerRad < 0.95);
285 
286  //
287  // Vertex assignments (v and w arrays)
288  // [0] and [1] are mandatory
289  // the rest are optional
290  //
291  //     +                     -
292  //      [0]------[4]------[1]      <--- outer radius
293  //       |                 |       
294  //       |                 |       
295  //      [2]---[5]---[6]---[3]      <--- inner radius
296  //
297
298  G4ClippablePolygon endPoly1, endPoly2;
299 
300  G4double phimax   = maxphi + 0.5*fDPhi;
301  G4double phimin   = - phimax;
302
303  G4ThreeVector v0, v1, v2, v3, v4, v5, v6;   // -ve phi verticies for polygon
304  G4ThreeVector w0, w1, w2, w3, w4, w5, w6;   // +ve phi verticies for polygon
305
306  //
307  // decide verticies of -ve phi boundary
308  //
309 
310  G4double cosPhi = std::cos(phimin);
311  G4double sinPhi = std::sin(phimin);
312
313  // Outer hyperbolic surface 
314
315  v0 = transform.TransformPoint( 
316       G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 
317                     + fZHalfLength));
318  v1 = transform.TransformPoint( 
319       G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 
320                     - fZHalfLength));
321  if (splitOuter)
322  {
323     v4 = transform.TransformPoint(
324          G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 0));
325  }
326 
327  // Inner hyperbolic surface 
328
329  G4double zInnerSplit = 0.;
330  if (splitInner)
331  {
332     v2 = transform.TransformPoint( 
333          G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, 
334                        + fZHalfLength));
335     v3 = transform.TransformPoint( 
336          G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, 
337                        - fZHalfLength));
338     
339     // Find intersection of tangential line of inner
340     // surface at z = fZHalfLength and line r=fInnerRadius.
341     G4double dr = fZHalfLength * fTanInnerStereo2;
342     G4double dz = maxEndInnerRad;
343     zInnerSplit = fZHalfLength + (fInnerRadius - maxEndInnerRad) * dz / dr;
344
345     // Build associated vertices
346     v5 = transform.TransformPoint( 
347          G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, 
348                        + zInnerSplit));
349     v6 = transform.TransformPoint( 
350          G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, 
351                        - zInnerSplit));
352  }
353  else
354  {
355     v2 = transform.TransformPoint(
356          G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, 
357                        + fZHalfLength));
358     v3 = transform.TransformPoint(
359          G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, 
360                        - fZHalfLength));
361  }
362
363  //
364  // decide vertices of +ve phi boundary
365  //
366
367  cosPhi = std::cos(phimax);
368  sinPhi = std::sin(phimax);
369 
370  // Outer hyperbolic surface 
371 
372  w0 = transform.TransformPoint(
373       G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
374                     + fZHalfLength));
375  w1 = transform.TransformPoint(
376       G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
377                     - fZHalfLength));
378  if (splitOuter)
379  {
380     G4double r = rFudge*fOuterRadius;
381     
382     w4 = transform.TransformPoint(G4ThreeVector( r*cosPhi, r*sinPhi, 0 ));
383     
384     AddPolyToExtent( v0, v4, w4, w0, voxelLimit, axis, extentList );
385     AddPolyToExtent( v4, v1, w1, w4, voxelLimit, axis, extentList );
386  }
387  else
388  {
389     AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList );
390  }
391 
392  // Inner hyperbolic surface
393 
394  if (splitInner)
395  {
396     w2 = transform.TransformPoint(
397          G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, 
398                        + fZHalfLength));
399     w3 = transform.TransformPoint(
400          G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, 
401                        - fZHalfLength));
402         
403     w5 = transform.TransformPoint(
404          G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
405                        + zInnerSplit));
406     w6 = transform.TransformPoint(
407          G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
408                        - zInnerSplit));
409                                   
410     AddPolyToExtent( v3, v6, w6, w3, voxelLimit, axis, extentList );
411     AddPolyToExtent( v6, v5, w5, w6, voxelLimit, axis, extentList );
412     AddPolyToExtent( v5, v2, w2, w5, voxelLimit, axis, extentList );
413     
414  }
415  else
416  {
417     w2 = transform.TransformPoint(
418          G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
419                        + fZHalfLength));
420     w3 = transform.TransformPoint(
421          G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
422                        - fZHalfLength));
423
424     AddPolyToExtent( v3, v2, w2, w3, voxelLimit, axis, extentList );
425  }
426
427  //
428  // Endplate segments
429  //
430  AddPolyToExtent( v1, v3, w3, w1, voxelLimit, axis, extentList );
431  AddPolyToExtent( v2, v0, w0, w2, voxelLimit, axis, extentList );
432 
433  //
434  // Return min/max value
435  //
436  return extentList.GetExtent( min, max );
437}
438
439
440//=====================================================================
441//* AddPolyToExtent ---------------------------------------------------
442
443void G4TwistedTubs::AddPolyToExtent( const G4ThreeVector &v0,
444                                     const G4ThreeVector &v1,
445                                     const G4ThreeVector &w1,
446                                     const G4ThreeVector &w0,
447                                     const G4VoxelLimits &voxelLimit,
448                                     const EAxis          axis,
449                                     G4SolidExtentList   &extentList ) 
450{
451    // Utility function for CalculateExtent
452    //
453    G4ClippablePolygon phiPoly;
454
455    phiPoly.AddVertexInOrder( v0 );
456    phiPoly.AddVertexInOrder( v1 );
457    phiPoly.AddVertexInOrder( w1 );
458    phiPoly.AddVertexInOrder( w0 );
459
460    if (phiPoly.PartialClip( voxelLimit, axis ))
461    {
462        phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
463        extentList.AddSurface( phiPoly );
464    }
465}
466
467
468//=====================================================================
469//* Inside ------------------------------------------------------------
470
471EInside G4TwistedTubs::Inside(const G4ThreeVector& p) const
472{
473
474   static const G4double halftol
475     = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
476   // static G4int timerid = -1;
477   // G4Timer timer(timerid, "G4TwistedTubs", "Inside");
478   // timer.Start();
479   
480
481
482   G4ThreeVector *tmpp;
483   EInside       *tmpinside;
484   if (fLastInside.p == p)
485   {
486     return fLastInside.inside;
487   }
488   else
489   {
490      tmpp      = const_cast<G4ThreeVector*>(&(fLastInside.p));
491      tmpinside = const_cast<EInside*>(&(fLastInside.inside));
492      tmpp->set(p.x(), p.y(), p.z());
493   }
494   
495   EInside  outerhypearea = ((G4TwistTubsHypeSide *)fOuterHype)->Inside(p);
496   G4double innerhyperho  = ((G4TwistTubsHypeSide *)fInnerHype)->GetRhoAtPZ(p);
497   G4double distanceToOut = p.getRho() - innerhyperho; // +ve: inside
498
499   if ((outerhypearea == kOutside) || (distanceToOut < -halftol))
500   {
501      *tmpinside = kOutside;
502   }
503   else if (outerhypearea == kSurface)
504   {
505      *tmpinside = kSurface;
506   }
507   else
508   {
509      if (distanceToOut <= halftol)
510      {
511         *tmpinside = kSurface;
512      }
513      else
514      {
515         *tmpinside = kInside;
516      }
517   }
518
519   return fLastInside.inside;
520}
521
522//=====================================================================
523//* SurfaceNormal -----------------------------------------------------
524
525G4ThreeVector G4TwistedTubs::SurfaceNormal(const G4ThreeVector& p) const
526{
527   //
528   // return the normal unit vector to the Hyperbolical Surface at a point
529   // p on (or nearly on) the surface
530   //
531   // Which of the three or four surfaces are we closest to?
532   //
533
534   if (fLastNormal.p == p)
535   {
536      return fLastNormal.vec;
537   }   
538   G4ThreeVector *tmpp          =
539     const_cast<G4ThreeVector*>(&(fLastNormal.p));
540   G4ThreeVector *tmpnormal     =
541     const_cast<G4ThreeVector*>(&(fLastNormal.vec));
542   G4VTwistSurface **tmpsurface =
543     const_cast<G4VTwistSurface**>(fLastNormal.surface);
544   tmpp->set(p.x(), p.y(), p.z());
545
546   G4double      distance = kInfinity;
547
548   G4VTwistSurface *surfaces[6];
549   surfaces[0] = fLatterTwisted;
550   surfaces[1] = fFormerTwisted;
551   surfaces[2] = fInnerHype;
552   surfaces[3] = fOuterHype;
553   surfaces[4] = fLowerEndcap;
554   surfaces[5] = fUpperEndcap;
555
556   G4ThreeVector xx;
557   G4ThreeVector bestxx;
558   G4int i;
559   G4int besti = -1;
560   for (i=0; i< 6; i++)
561   {
562      G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
563      if (tmpdistance < distance)
564      {
565         distance = tmpdistance;
566         bestxx = xx;
567         besti = i; 
568      }
569   }
570
571   tmpsurface[0] = surfaces[besti];
572   *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
573   
574   return fLastNormal.vec;
575}
576
577//=====================================================================
578//* DistanceToIn (p, v) -----------------------------------------------
579
580G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p,
581                                      const G4ThreeVector& v ) const
582{
583
584   // DistanceToIn (p, v):
585   // Calculate distance to surface of shape from `outside'
586   // along with the v, allowing for tolerance.
587   // The function returns kInfinity if no intersection or
588   // just grazing within tolerance.
589
590   //
591   // checking last value
592   //
593   
594   G4ThreeVector *tmpp;
595   G4ThreeVector *tmpv;
596   G4double      *tmpdist;
597   if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v))
598   {
599     return fLastDistanceToIn.value;
600   }
601   else
602   {
603      tmpp    = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
604      tmpv    = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
605      tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
606      tmpp->set(p.x(), p.y(), p.z());
607      tmpv->set(v.x(), v.y(), v.z());
608   }
609
610   //
611   // Calculate DistanceToIn(p,v)
612   //
613   
614   EInside currentside = Inside(p);
615
616   if (currentside == kInside)
617   {
618   }
619   else
620   {
621     if (currentside == kSurface)
622     {
623       // particle is just on a boundary.
624       // If the particle is entering to the volume, return 0.
625       //
626       G4ThreeVector normal = SurfaceNormal(p);
627       if (normal*v < 0)
628       {
629         *tmpdist = 0;
630         return fLastDistanceToInWithV.value;
631       } 
632     }
633   }
634     
635   // now, we can take smallest positive distance.
636   
637   // Initialize
638   //
639   G4double      distance = kInfinity;   
640
641   // find intersections and choose nearest one.
642   //
643   G4VTwistSurface *surfaces[6];
644   surfaces[0] = fLowerEndcap;
645   surfaces[1] = fUpperEndcap;
646   surfaces[2] = fLatterTwisted;
647   surfaces[3] = fFormerTwisted;
648   surfaces[4] = fInnerHype;
649   surfaces[5] = fOuterHype;
650   
651   G4ThreeVector xx;
652   G4ThreeVector bestxx;
653   G4int i;
654   G4int besti = -1;
655   for (i=0; i< 6; i++)
656   {
657      G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
658      if (tmpdistance < distance)
659      {
660         distance = tmpdistance;
661         bestxx = xx;
662         besti = i;
663      }
664   }
665   *tmpdist = distance;
666
667   return fLastDistanceToInWithV.value;
668}
669 
670//=====================================================================
671//* DistanceToIn (p) --------------------------------------------------
672
673G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p) const
674{
675   // DistanceToIn(p):
676   // Calculate distance to surface of shape from `outside',
677   // allowing for tolerance
678   
679   //
680   // checking last value
681   //
682   
683   G4ThreeVector *tmpp;
684   G4double      *tmpdist;
685   if (fLastDistanceToIn.p == p)
686   {
687     return fLastDistanceToIn.value;
688   }
689   else
690   {
691      tmpp    = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
692      tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
693      tmpp->set(p.x(), p.y(), p.z());
694   }
695
696   //
697   // Calculate DistanceToIn(p)
698   //
699   
700   EInside currentside = Inside(p);
701
702   switch (currentside)
703   {
704      case (kInside) :
705      {}
706      case (kSurface) :
707      {
708         *tmpdist = 0.;
709         return fLastDistanceToIn.value;
710      }
711      case (kOutside) :
712      {
713         // Initialize
714         G4double      distance = kInfinity;   
715
716         // find intersections and choose nearest one.
717         G4VTwistSurface *surfaces[6];
718         surfaces[0] = fLowerEndcap;
719         surfaces[1] = fUpperEndcap;
720         surfaces[2] = fLatterTwisted;
721         surfaces[3] = fFormerTwisted;
722         surfaces[4] = fInnerHype;
723         surfaces[5] = fOuterHype;
724
725         G4int i;
726         G4int besti = -1;
727         G4ThreeVector xx;
728         G4ThreeVector bestxx;
729         for (i=0; i< 6; i++)
730         {
731            G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
732            if (tmpdistance < distance)
733            {
734               distance = tmpdistance;
735               bestxx = xx;
736               besti = i;
737            }
738         }
739         *tmpdist = distance;
740         return fLastDistanceToIn.value;
741      }
742      default :
743      {
744         G4Exception("G4TwistedTubs::DistanceToIn(p)", "InvalidCondition",
745                     FatalException, "Unknown point location!");
746      }
747   } // switch end
748
749   return kInfinity;
750}
751
752//=====================================================================
753//* DistanceToOut (p, v) ----------------------------------------------
754
755G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p,
756                                       const G4ThreeVector& v,
757                                       const G4bool calcNorm,
758                                       G4bool *validNorm,
759                                       G4ThreeVector *norm ) const
760{
761   // DistanceToOut (p, v):
762   // Calculate distance to surface of shape from `inside'
763   // along with the v, allowing for tolerance.
764   // The function returns kInfinity if no intersection or
765   // just grazing within tolerance.
766
767   //
768   // checking last value
769   //
770   
771   G4ThreeVector *tmpp;
772   G4ThreeVector *tmpv;
773   G4double      *tmpdist;
774   if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) )
775   {
776     return fLastDistanceToOutWithV.value;
777   }
778   else
779   {
780      tmpp    = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
781      tmpv    = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
782      tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
783      tmpp->set(p.x(), p.y(), p.z());
784      tmpv->set(v.x(), v.y(), v.z());
785   }
786
787   //
788   // Calculate DistanceToOut(p,v)
789   //
790   
791   EInside currentside = Inside(p);
792
793   if (currentside == kOutside)
794   {
795   }
796   else
797   {
798     if (currentside == kSurface)
799     {
800       // particle is just on a boundary.
801       // If the particle is exiting from the volume, return 0.
802       //
803       G4ThreeVector normal = SurfaceNormal(p);
804       G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
805       if (normal*v > 0)
806       {
807         if (calcNorm)
808         {
809           *norm = (blockedsurface->GetNormal(p, true));
810           *validNorm = blockedsurface->IsValidNorm();
811         }
812         *tmpdist = 0.;
813         return fLastDistanceToOutWithV.value;
814       }
815     }
816   }
817     
818   // now, we can take smallest positive distance.
819   
820   // Initialize
821   //
822   G4double      distance = kInfinity;
823       
824   // find intersections and choose nearest one.
825   //
826   G4VTwistSurface *surfaces[6];
827   surfaces[0] = fLatterTwisted;
828   surfaces[1] = fFormerTwisted;
829   surfaces[2] = fInnerHype;
830   surfaces[3] = fOuterHype;
831   surfaces[4] = fLowerEndcap;
832   surfaces[5] = fUpperEndcap;
833   
834   G4int i;
835   G4int besti = -1;
836   G4ThreeVector xx;
837   G4ThreeVector bestxx;
838   for (i=0; i< 6; i++)
839   {
840      G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
841      if (tmpdistance < distance)
842      {
843         distance = tmpdistance;
844         bestxx = xx; 
845         besti = i;
846      }
847   }
848
849   if (calcNorm)
850   {
851      if (besti != -1)
852      {
853         *norm = (surfaces[besti]->GetNormal(p, true));
854         *validNorm = surfaces[besti]->IsValidNorm();
855      }
856   }
857
858   *tmpdist = distance;
859
860   return fLastDistanceToOutWithV.value;
861}
862
863
864//=====================================================================
865//* DistanceToOut (p) ----------------------------------------------
866
867G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p ) const
868{
869   // DistanceToOut(p):
870   // Calculate distance to surface of shape from `inside',
871   // allowing for tolerance
872
873   //
874   // checking last value
875   //
876   
877   G4ThreeVector *tmpp;
878   G4double      *tmpdist;
879   if (fLastDistanceToOut.p == p)
880   {
881      return fLastDistanceToOut.value;
882   }
883   else
884   {
885      tmpp    = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
886      tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
887      tmpp->set(p.x(), p.y(), p.z());
888   }
889   
890   //
891   // Calculate DistanceToOut(p)
892   //
893   
894   EInside currentside = Inside(p);
895
896   switch (currentside)
897   {
898      case (kOutside) :
899      {
900      }
901      case (kSurface) :
902      {
903        *tmpdist = 0.;
904         return fLastDistanceToOut.value;
905      }
906      case (kInside) :
907      {
908         // Initialize
909         G4double      distance = kInfinity;
910   
911         // find intersections and choose nearest one.
912         G4VTwistSurface *surfaces[6];
913         surfaces[0] = fLatterTwisted;
914         surfaces[1] = fFormerTwisted;
915         surfaces[2] = fInnerHype;
916         surfaces[3] = fOuterHype;
917         surfaces[4] = fLowerEndcap;
918         surfaces[5] = fUpperEndcap;
919
920         G4int i;
921         G4int besti = -1;
922         G4ThreeVector xx;
923         G4ThreeVector bestxx;
924         for (i=0; i< 6; i++)
925         {
926            G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
927            if (tmpdistance < distance)
928            {
929               distance = tmpdistance;
930               bestxx = xx;
931               besti = i;
932            }
933         }
934         *tmpdist = distance;
935   
936         return fLastDistanceToOut.value;
937      }
938      default :
939      {
940         G4Exception("G4TwistedTubs::DistanceToOut(p)", "InvalidCondition",
941                     FatalException, "Unknown point location!");
942      }
943   } // switch end
944
945   return 0;
946}
947
948//=====================================================================
949//* StreamInfo --------------------------------------------------------
950
951std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
952{
953  //
954  // Stream object contents to an output stream
955  //
956  os << "-----------------------------------------------------------\n"
957     << "    *** Dump for solid - " << GetName() << " ***\n"
958     << "    ===================================================\n"
959     << " Solid type: G4TwistedTubs\n"
960     << " Parameters: \n"
961     << "    -ve end Z              : " << fEndZ[0]/mm << " mm \n"
962     << "    +ve end Z              : " << fEndZ[1]/mm << " mm \n"
963     << "    inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
964     << "    inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
965     << "    outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
966     << "    outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
967     << "    inner radius (z=0)     : " << fInnerRadius/mm << " mm \n"
968     << "    outer radius (z=0)     : " << fOuterRadius/mm << " mm \n"
969     << "    twisted angle          : " << fPhiTwist/degree << " degrees \n"
970     << "    inner stereo angle     : " << fInnerStereo/degree << " degrees \n"
971     << "    outer stereo angle     : " << fOuterStereo/degree << " degrees \n"
972     << "    phi-width of a piece   : " << fDPhi/degree << " degrees \n"
973     << "-----------------------------------------------------------\n";
974
975  return os;
976}
977
978
979//=====================================================================
980//* DiscribeYourselfTo ------------------------------------------------
981
982void G4TwistedTubs::DescribeYourselfTo (G4VGraphicsScene& scene) const 
983{
984  scene.AddSolid (*this);
985}
986
987//=====================================================================
988//* GetExtent ---------------------------------------------------------
989
990G4VisExtent G4TwistedTubs::GetExtent() const 
991{
992  // Define the sides of the box into which the G4Tubs instance would fit.
993
994  G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
995  return G4VisExtent( -maxEndOuterRad, maxEndOuterRad, 
996                      -maxEndOuterRad, maxEndOuterRad, 
997                      -fZHalfLength, fZHalfLength );
998}
999
1000//=====================================================================
1001//* CreatePolyhedron --------------------------------------------------
1002
1003G4Polyhedron* G4TwistedTubs::CreatePolyhedron () const 
1004{
1005  // number of meshes
1006  //
1007  G4double dA = std::max(fDPhi,fPhiTwist);
1008  const G4int m =
1009    G4int(G4Polyhedron::GetNumberOfRotationSteps() * dA / twopi) + 2;
1010  const G4int n =
1011    G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2;
1012
1013  const G4int nnodes = 4*(m-1)*(n-2) + 2*m*m ;
1014  const G4int nfaces = 4*(m-1)*(n-1) + 2*(m-1)*(m-1) ;
1015
1016  G4Polyhedron *ph=new G4Polyhedron;
1017  typedef G4double G4double3[3];
1018  typedef G4int G4int4[4];
1019  G4double3* xyz = new G4double3[nnodes];  // number of nodes
1020  G4int4*  faces = new G4int4[nfaces] ;    // number of faces
1021  fLowerEndcap->GetFacets(m,m,xyz,faces,0) ;
1022  fUpperEndcap->GetFacets(m,m,xyz,faces,1) ;
1023  fInnerHype->GetFacets(m,n,xyz,faces,2) ;
1024  fFormerTwisted->GetFacets(m,n,xyz,faces,3) ;
1025  fOuterHype->GetFacets(m,n,xyz,faces,4) ;
1026  fLatterTwisted->GetFacets(m,n,xyz,faces,5) ;
1027
1028  ph->createPolyhedron(nnodes,nfaces,xyz,faces);
1029
1030  delete[] xyz;
1031  delete[] faces;
1032
1033  return ph;
1034}
1035
1036//=====================================================================
1037//* CreateNUBS --------------------------------------------------------
1038
1039G4NURBS* G4TwistedTubs::CreateNURBS () const 
1040{
1041   G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
1042   G4double maxEndInnerRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
1043   return new G4NURBStube(maxEndInnerRad, maxEndOuterRad, fZHalfLength); 
1044   // Tube for now!!!
1045}
1046
1047//=====================================================================
1048//* GetPolyhedron -----------------------------------------------------
1049
1050G4Polyhedron* G4TwistedTubs::GetPolyhedron () const
1051{
1052  if ((!fpPolyhedron) ||
1053      (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1054       fpPolyhedron->GetNumberOfRotationSteps()))
1055  {
1056    delete fpPolyhedron;
1057    fpPolyhedron = CreatePolyhedron();
1058  }
1059  return fpPolyhedron;
1060}
1061
1062//=====================================================================
1063//* CreateSurfaces ----------------------------------------------------
1064
1065void G4TwistedTubs::CreateSurfaces() 
1066{
1067   // create 6 surfaces of TwistedTub
1068
1069   G4ThreeVector x0(0, 0, fEndZ[0]);
1070   G4ThreeVector n (0, 0, -1);
1071
1072   fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
1073                                    fEndInnerRadius, fEndOuterRadius,
1074                                    fDPhi, fEndPhi, fEndZ, -1) ;
1075
1076   fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap", 
1077                                    fEndInnerRadius, fEndOuterRadius,
1078                                    fDPhi, fEndPhi, fEndZ, 1) ;
1079
1080   G4RotationMatrix    rotHalfDPhi;
1081   rotHalfDPhi.rotateZ(0.5*fDPhi);
1082
1083   fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
1084                                         fEndInnerRadius, fEndOuterRadius,
1085                                         fDPhi, fEndPhi, fEndZ, 
1086                                         fInnerRadius, fOuterRadius, fKappa,
1087                                         1 ) ; 
1088   fFormerTwisted = new G4TwistTubsSide("FormerTwisted", 
1089                                         fEndInnerRadius, fEndOuterRadius,
1090                                         fDPhi, fEndPhi, fEndZ, 
1091                                         fInnerRadius, fOuterRadius, fKappa,
1092                                         -1 ) ; 
1093
1094   fInnerHype = new G4TwistTubsHypeSide("InnerHype",
1095                                        fEndInnerRadius, fEndOuterRadius,
1096                                        fDPhi, fEndPhi, fEndZ, 
1097                                        fInnerRadius, fOuterRadius,fKappa,
1098                                        fTanInnerStereo, fTanOuterStereo, -1) ;
1099   fOuterHype = new G4TwistTubsHypeSide("OuterHype", 
1100                                        fEndInnerRadius, fEndOuterRadius,
1101                                        fDPhi, fEndPhi, fEndZ, 
1102                                        fInnerRadius, fOuterRadius,fKappa,
1103                                        fTanInnerStereo, fTanOuterStereo, 1) ;
1104
1105
1106   // set neighbour surfaces
1107   //
1108   fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
1109                               fOuterHype, fFormerTwisted);
1110   fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
1111                               fOuterHype, fFormerTwisted);
1112   fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
1113                                 fOuterHype, fUpperEndcap);
1114   fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
1115                                 fOuterHype, fUpperEndcap);
1116   fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
1117                             fFormerTwisted, fUpperEndcap);
1118   fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
1119                             fFormerTwisted, fUpperEndcap);
1120}
1121
1122
1123//=====================================================================
1124//* GetEntityType -----------------------------------------------------
1125
1126G4GeometryType G4TwistedTubs::GetEntityType() const
1127{
1128  return G4String("G4TwistedTubs");
1129}
1130
1131//=====================================================================
1132//* GetCubicVolume ----------------------------------------------------
1133
1134G4double G4TwistedTubs::GetCubicVolume()
1135{
1136  if(fCubicVolume != 0.) {;}
1137  else   { fCubicVolume = fDPhi*fZHalfLength*(fOuterRadius*fOuterRadius
1138                                             -fInnerRadius*fInnerRadius); }
1139  return fCubicVolume;
1140}
1141
1142//=====================================================================
1143//* GetSurfaceArea ----------------------------------------------------
1144
1145G4double G4TwistedTubs::GetSurfaceArea()
1146{
1147  if(fSurfaceArea != 0.) {;}
1148  else   { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
1149  return fSurfaceArea;
1150}
1151
1152//=====================================================================
1153//* GetPointOnSurface -------------------------------------------------
1154
1155G4ThreeVector G4TwistedTubs::GetPointOnSurface() const
1156{
1157
1158  G4double  z = CLHEP::RandFlat::shoot(fEndZ[0],fEndZ[1]);
1159  G4double phi , phimin, phimax ;
1160  G4double x   , xmin,   xmax ;
1161  G4double r   , rmin,   rmax ;
1162
1163  G4double a1 = fOuterHype->GetSurfaceArea() ;
1164  G4double a2 = fInnerHype->GetSurfaceArea() ;
1165  G4double a3 = fLatterTwisted->GetSurfaceArea() ;
1166  G4double a4 = fFormerTwisted->GetSurfaceArea() ;
1167  G4double a5 = fLowerEndcap->GetSurfaceArea()  ;
1168  G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1169
1170  G4double chose = CLHEP::RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1171
1172  if(chose < a1)
1173  {
1174
1175    phimin = fOuterHype->GetBoundaryMin(z) ;
1176    phimax = fOuterHype->GetBoundaryMax(z) ;
1177    phi = CLHEP::RandFlat::shoot(phimin,phimax) ;
1178
1179    return fOuterHype->SurfacePoint(phi,z,true) ;
1180
1181  }
1182  else if ( (chose >= a1) && (chose < a1 + a2 ) )
1183  {
1184
1185    phimin = fInnerHype->GetBoundaryMin(z) ;
1186    phimax = fInnerHype->GetBoundaryMax(z) ;
1187    phi = CLHEP::RandFlat::shoot(phimin,phimax) ;
1188
1189    return fInnerHype->SurfacePoint(phi,z,true) ;
1190
1191  }
1192  else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) ) 
1193  {
1194
1195    xmin = fLatterTwisted->GetBoundaryMin(z) ; 
1196    xmax = fLatterTwisted->GetBoundaryMax(z) ;
1197    x = CLHEP::RandFlat::shoot(xmin,xmax) ;
1198   
1199    return fLatterTwisted->SurfacePoint(x,z,true) ;
1200
1201  }
1202  else if ( (chose >= a1 + a2 + a3  ) && (chose < a1 + a2 + a3 + a4  ) )
1203  {
1204
1205    xmin = fFormerTwisted->GetBoundaryMin(z) ; 
1206    xmax = fFormerTwisted->GetBoundaryMax(z) ;
1207    x = CLHEP::RandFlat::shoot(xmin,xmax) ;
1208
1209    return fFormerTwisted->SurfacePoint(x,z,true) ;
1210 
1211  }
1212  else if( (chose >= a1 + a2 + a3 + a4  )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1213  {
1214
1215    rmin = GetEndInnerRadius(0) ;
1216    rmax = GetEndOuterRadius(0) ;
1217    r = CLHEP::RandFlat::shoot(rmin,rmax) ;
1218
1219    phimin = fLowerEndcap->GetBoundaryMin(r) ; 
1220    phimax = fLowerEndcap->GetBoundaryMax(r) ;
1221    phi    = CLHEP::RandFlat::shoot(phimin,phimax) ;
1222
1223    return fLowerEndcap->SurfacePoint(phi,r,true) ;
1224
1225  }
1226  else
1227  {
1228    rmin = GetEndInnerRadius(1) ;
1229    rmax = GetEndOuterRadius(1) ;
1230    r = CLHEP::RandFlat::shoot(rmin,rmax) ;
1231
1232    phimin = fUpperEndcap->GetBoundaryMin(r) ; 
1233    phimax = fUpperEndcap->GetBoundaryMax(r) ;
1234    phi    = CLHEP::RandFlat::shoot(phimin,phimax) ;
1235
1236    return fUpperEndcap->SurfacePoint(phi,r,true) ;
1237  }
1238}
Note: See TracBrowser for help on using the repository browser.