Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4VTwistedFaceted.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 // $Id: G4VTwistedFaceted.cc 99781 2016-10-05 10:18:54Z gcosmo $
27 //
28 //
29 // --------------------------------------------------------------------
30 // GEANT 4 class source file
31 //
32 //
33 // G4VTwistedFaceted.cc
34 //
35 // Author:
36 //
37 // 04-Nov-2004 - O.Link (Oliver.Link@cern.ch)
38 //
39 // --------------------------------------------------------------------
40 
41 #include "G4VTwistedFaceted.hh"
42 
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4VoxelLimits.hh"
46 #include "G4AffineTransform.hh"
47 #include "G4BoundingEnvelope.hh"
48 #include "G4SolidExtentList.hh"
49 #include "G4ClippablePolygon.hh"
50 #include "G4VPVParameterisation.hh"
51 #include "G4GeometryTolerance.hh"
52 #include "meshdefs.hh"
53 
54 #include "G4VGraphicsScene.hh"
55 #include "G4Polyhedron.hh"
56 #include "G4VisExtent.hh"
57 
58 #include "Randomize.hh"
59 
60 #include "G4AutoLock.hh"
61 
62 namespace
63 {
64  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
65 }
66 
67 //=====================================================================
68 //* constructors ------------------------------------------------------
69 
71 G4VTwistedFaceted( const G4String &pname, // Name of instance
72  G4double PhiTwist, // twist angle
73  G4double pDz, // half z length
74  G4double pTheta, // direction between end planes
75  G4double pPhi, // defined by polar and azim. angles
76  G4double pDy1, // half y length at -pDz
77  G4double pDx1, // half x length at -pDz,-pDy
78  G4double pDx2, // half x length at -pDz,+pDy
79  G4double pDy2, // half y length at +pDz
80  G4double pDx3, // half x length at +pDz,-pDy
81  G4double pDx4, // half x length at +pDz,+pDy
82  G4double pAlph // tilt angle
83  )
84  : G4VSolid(pname), fRebuildPolyhedron(false), fpPolyhedron(0),
85  fLowerEndcap(0), fUpperEndcap(0), fSide0(0),
86  fSide90(0), fSide180(0), fSide270(0),
87  fSurfaceArea(0.)
88 {
89 
90  G4double pDytmp ;
91  G4double fDxUp ;
92  G4double fDxDown ;
93 
94  fDx1 = pDx1 ;
95  fDx2 = pDx2 ;
96  fDx3 = pDx3 ;
97  fDx4 = pDx4 ;
98  fDy1 = pDy1 ;
99  fDy2 = pDy2 ;
100  fDz = pDz ;
101 
102  G4double kAngTolerance
104 
105  // maximum values
106  //
107  fDxDown = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ;
108  fDxUp = ( fDx3 > fDx4 ? fDx3 : fDx4 ) ;
109  fDx = ( fDxUp > fDxDown ? fDxUp : fDxDown ) ;
110  fDy = ( fDy1 > fDy2 ? fDy1 : fDy2 ) ;
111 
112  // planarity check
113  //
114  if ( fDx1 != fDx2 && fDx3 != fDx4 )
115  {
116  pDytmp = fDy1 * ( fDx3 - fDx4 ) / ( fDx1 - fDx2 ) ;
117  if ( std::fabs(pDytmp - fDy2) > kCarTolerance )
118  {
119  std::ostringstream message;
120  message << "Not planar surface in untwisted Trapezoid: "
121  << GetName() << G4endl
122  << "fDy2 is " << fDy2 << " but should be "
123  << pDytmp << ".";
124  G4Exception("G4VTwistedFaceted::G4VTwistedFaceted()", "GeomSolids0002",
125  FatalErrorInArgument, message);
126  }
127  }
128 
129 #ifdef G4TWISTDEBUG
130  if ( fDx1 == fDx2 && fDx3 == fDx4 )
131  {
132  G4cout << "Trapezoid is a box" << G4endl ;
133  }
134 
135 #endif
136 
137  if ( ( fDx1 == fDx2 && fDx3 != fDx4 ) || ( fDx1 != fDx2 && fDx3 == fDx4 ) )
138  {
139  std::ostringstream message;
140  message << "Not planar surface in untwisted Trapezoid: "
141  << GetName() << G4endl
142  << "One endcap is rectangular, the other is a trapezoid." << G4endl
143  << "For planarity reasons they have to be rectangles or trapezoids "
144  << "on both sides.";
145  G4Exception("G4VTwistedFaceted::G4VTwistedFaceted()", "GeomSolids0002",
146  FatalErrorInArgument, message);
147  }
148 
149  // twist angle
150  //
151  fPhiTwist = PhiTwist ;
152 
153  // tilt angle
154  //
155  fAlph = pAlph ;
156  fTAlph = std::tan(fAlph) ;
157 
158  fTheta = pTheta ;
159  fPhi = pPhi ;
160 
161  // dx in surface equation
162  //
163  fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
164 
165  // dy in surface equation
166  //
167  fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
168 
169  if ( ! ( ( fDx1 > 2*kCarTolerance)
170  && ( fDx2 > 2*kCarTolerance)
171  && ( fDx3 > 2*kCarTolerance)
172  && ( fDx4 > 2*kCarTolerance)
173  && ( fDy1 > 2*kCarTolerance)
174  && ( fDy2 > 2*kCarTolerance)
175  && ( fDz > 2*kCarTolerance)
176  && ( std::fabs(fPhiTwist) > 2*kAngTolerance )
177  && ( std::fabs(fPhiTwist) < pi/2 )
178  && ( std::fabs(fAlph) < pi/2 )
179  && ( fTheta < pi/2 && fTheta >= 0 ) )
180  )
181  {
182  std::ostringstream message;
183  message << "Invalid dimensions. Too small, or twist angle too big: "
184  << GetName() << G4endl
185  << "fDx 1-4 = " << fDx1/cm << ", " << fDx2/cm << ", "
186  << fDx3/cm << ", " << fDx4/cm << " cm" << G4endl
187  << "fDy 1-2 = " << fDy1/cm << ", " << fDy2/cm << ", "
188  << " cm" << G4endl
189  << "fDz = " << fDz/cm << " cm" << G4endl
190  << " twistangle " << fPhiTwist/deg << " deg" << G4endl
191  << " phi,theta = " << fPhi/deg << ", " << fTheta/deg << " deg";
192  G4Exception("G4TwistedTrap::G4VTwistedFaceted()",
193  "GeomSolids0002", FatalErrorInArgument, message);
194  }
195  CreateSurfaces();
196  fCubicVolume = 2 * fDz * ( ( fDx1 + fDx2 ) * fDy1 + ( fDx3 + fDx4 ) * fDy2 );
197 }
198 
199 
200 //=====================================================================
201 //* Fake default constructor ------------------------------------------
202 
204  : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0),
205  fTheta(0.), fPhi(0.), fDy1(0.),
206  fDx1(0.), fDx2(0.), fDy2(0.), fDx3(0.), fDx4(0.),
207  fDz(0.), fDx(0.), fDy(0.), fAlph(0.),
208  fTAlph(0.), fdeltaX(0.), fdeltaY(0.), fPhiTwist(0.),
209  fLowerEndcap(0), fUpperEndcap(0), fSide0(0), fSide90(0), fSide180(0),
210  fSide270(0), fCubicVolume(0.), fSurfaceArea(0.)
211 {
212 }
213 
214 
215 //=====================================================================
216 //* destructor --------------------------------------------------------
217 
219 {
220  if (fLowerEndcap) { delete fLowerEndcap ; }
221  if (fUpperEndcap) { delete fUpperEndcap ; }
222 
223  if (fSide0) { delete fSide0 ; }
224  if (fSide90) { delete fSide90 ; }
225  if (fSide180) { delete fSide180 ; }
226  if (fSide270) { delete fSide270 ; }
227  if (fpPolyhedron) { delete fpPolyhedron; fpPolyhedron = 0; }
228 }
229 
230 
231 //=====================================================================
232 //* Copy constructor --------------------------------------------------
233 
235  : G4VSolid(rhs), fRebuildPolyhedron(false), fpPolyhedron(0),
236  fTheta(rhs.fTheta), fPhi(rhs.fPhi),
237  fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fDy2(rhs.fDy2),
238  fDx3(rhs.fDx3), fDx4(rhs.fDx4), fDz(rhs.fDz), fDx(rhs.fDx), fDy(rhs.fDy),
239  fAlph(rhs.fAlph), fTAlph(rhs.fTAlph), fdeltaX(rhs.fdeltaX),
240  fdeltaY(rhs.fdeltaY), fPhiTwist(rhs.fPhiTwist), fLowerEndcap(0),
241  fUpperEndcap(0), fSide0(0), fSide90(0), fSide180(0), fSide270(0),
242  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
243  fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
244  fLastDistanceToIn(rhs.fLastDistanceToIn),
245  fLastDistanceToOut(rhs.fLastDistanceToOut),
246  fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
247  fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
248 {
249  CreateSurfaces();
250 }
251 
252 
253 //=====================================================================
254 //* Assignment operator -----------------------------------------------
255 
257 {
258  // Check assignment to self
259  //
260  if (this == &rhs) { return *this; }
261 
262  // Copy base class data
263  //
264  G4VSolid::operator=(rhs);
265 
266  // Copy data
267  //
268  fTheta = rhs.fTheta; fPhi = rhs.fPhi;
269  fDy1= rhs.fDy1; fDx1= rhs.fDx1; fDx2= rhs.fDx2; fDy2= rhs.fDy2;
270  fDx3= rhs.fDx3; fDx4= rhs.fDx4; fDz= rhs.fDz; fDx= rhs.fDx; fDy= rhs.fDy;
271  fAlph= rhs.fAlph; fTAlph= rhs.fTAlph; fdeltaX= rhs.fdeltaX;
272  fdeltaY= rhs.fdeltaY; fPhiTwist= rhs.fPhiTwist; fLowerEndcap= 0;
273  fUpperEndcap= 0; fSide0= 0; fSide90= 0; fSide180= 0; fSide270= 0;
274  fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
275  fRebuildPolyhedron = false;
276  delete fpPolyhedron; fpPolyhedron= 0;
277  fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
278  fLastDistanceToIn= rhs.fLastDistanceToIn;
279  fLastDistanceToOut= rhs.fLastDistanceToOut;
280  fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
281  fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
282 
283  CreateSurfaces();
284 
285  return *this;
286 }
287 
288 
289 //=====================================================================
290 //* ComputeDimensions -------------------------------------------------
291 
293  const G4int ,
294  const G4VPhysicalVolume* )
295 {
296  G4Exception("G4VTwistedFaceted::ComputeDimensions()",
297  "GeomSolids0001", FatalException,
298  "G4VTwistedFaceted does not support Parameterisation.");
299 }
300 
301 
302 //=====================================================================
303 //* Extent ------------------------------------------------------------
304 
306  G4ThreeVector &pMax) const
307 {
308  G4double maxRad = std::sqrt(fDx*fDx + fDy*fDy);
309  pMin.set(-maxRad,-maxRad,-fDz);
310  pMax.set( maxRad, maxRad, fDz);
311 }
312 
313 
314 //=====================================================================
315 //* CalculateExtent ---------------------------------------------------
316 
317 G4bool
319  const G4VoxelLimits &pVoxelLimit,
320  const G4AffineTransform &pTransform,
321  G4double &pMin,
322  G4double &pMax ) const
323 {
324  G4ThreeVector bmin, bmax;
325 
326  // Get bounding box
327  Extent(bmin,bmax);
328 
329  // Find extent
330  G4BoundingEnvelope bbox(bmin,bmax);
331  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
332 }
333 
334 
335 //=====================================================================
336 //* Inside ------------------------------------------------------------
337 
339 {
340 
341  G4ThreeVector *tmpp;
342  EInside *tmpin;
343  if (fLastInside.p == p) {
344  return fLastInside.inside;
345  } else {
346  tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
347  tmpin = const_cast<EInside*>(&(fLastInside.inside));
348  tmpp->set(p.x(), p.y(), p.z());
349  }
350 
351  *tmpin = kOutside ;
352 
353  G4double phi = p.z()/(2*fDz) * fPhiTwist ; // rotate the point to z=0
354  G4double cphi = std::cos(-phi) ;
355  G4double sphi = std::sin(-phi) ;
356 
357  G4double px = p.x() + fdeltaX * ( -phi/fPhiTwist) ; // shift
358  G4double py = p.y() + fdeltaY * ( -phi/fPhiTwist) ;
359  G4double pz = p.z() ;
360 
361  G4double posx = px * cphi - py * sphi ; // rotation
362  G4double posy = px * sphi + py * cphi ;
363  G4double posz = pz ;
364 
365  G4double xMin = Xcoef(posy,phi,fTAlph) - 2*Xcoef(posy,phi,0.) ;
366  G4double xMax = Xcoef(posy,phi,fTAlph) ;
367 
368  G4double yMax = GetValueB(phi)/2. ; // b(phi)/2 is limit
369  G4double yMin = -yMax ;
370 
371 #ifdef G4TWISTDEBUG
372 
373  G4cout << "inside called: p = " << p << G4endl ;
374  G4cout << "fDx1 = " << fDx1 << G4endl ;
375  G4cout << "fDx2 = " << fDx2 << G4endl ;
376  G4cout << "fDx3 = " << fDx3 << G4endl ;
377  G4cout << "fDx4 = " << fDx4 << G4endl ;
378 
379  G4cout << "fDy1 = " << fDy1 << G4endl ;
380  G4cout << "fDy2 = " << fDy2 << G4endl ;
381 
382  G4cout << "fDz = " << fDz << G4endl ;
383 
384  G4cout << "Tilt angle alpha = " << fAlph << G4endl ;
385  G4cout << "phi,theta = " << fPhi << " , " << fTheta << G4endl ;
386 
387  G4cout << "Twist angle = " << fPhiTwist << G4endl ;
388 
389  G4cout << "posx = " << posx << G4endl ;
390  G4cout << "posy = " << posy << G4endl ;
391  G4cout << "xMin = " << xMin << G4endl ;
392  G4cout << "xMax = " << xMax << G4endl ;
393  G4cout << "yMin = " << yMin << G4endl ;
394  G4cout << "yMax = " << yMax << G4endl ;
395 
396 #endif
397 
398 
399  if ( posx <= xMax - kCarTolerance*0.5
400  && posx >= xMin + kCarTolerance*0.5 )
401  {
402  if ( posy <= yMax - kCarTolerance*0.5
403  && posy >= yMin + kCarTolerance*0.5 )
404  {
405  if (std::fabs(posz) <= fDz - kCarTolerance*0.5 ) *tmpin = kInside ;
406  else if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) *tmpin = kSurface ;
407  }
408  else if ( posy <= yMax + kCarTolerance*0.5
409  && posy >= yMin - kCarTolerance*0.5 )
410  {
411  if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) *tmpin = kSurface ;
412  }
413  }
414  else if ( posx <= xMax + kCarTolerance*0.5
415  && posx >= xMin - kCarTolerance*0.5 )
416  {
417  if ( posy <= yMax + kCarTolerance*0.5
418  && posy >= yMin - kCarTolerance*0.5 )
419  {
420  if (std::fabs(posz) <= fDz + kCarTolerance*0.5) *tmpin = kSurface ;
421  }
422  }
423 
424 #ifdef G4TWISTDEBUG
425  G4cout << "inside = " << fLastInside.inside << G4endl ;
426 #endif
427 
428  return fLastInside.inside;
429 
430 }
431 
432 
433 //=====================================================================
434 //* SurfaceNormal -----------------------------------------------------
435 
437 {
438  //
439  // return the normal unit vector to the Hyperbolical Surface at a point
440  // p on (or nearly on) the surface
441  //
442  // Which of the three or four surfaces are we closest to?
443  //
444 
445  if (fLastNormal.p == p)
446  {
447  return fLastNormal.vec;
448  }
449 
450  G4ThreeVector *tmpp = const_cast<G4ThreeVector*>(&(fLastNormal.p));
451  G4ThreeVector *tmpnormal = const_cast<G4ThreeVector*>(&(fLastNormal.vec));
452  G4VTwistSurface **tmpsurface = const_cast<G4VTwistSurface**>(fLastNormal.surface);
453  tmpp->set(p.x(), p.y(), p.z());
454 
455  G4double distance = kInfinity;
456 
457  G4VTwistSurface *surfaces[6];
458 
459  surfaces[0] = fSide0 ;
460  surfaces[1] = fSide90 ;
461  surfaces[2] = fSide180 ;
462  surfaces[3] = fSide270 ;
463  surfaces[4] = fLowerEndcap;
464  surfaces[5] = fUpperEndcap;
465 
466  G4ThreeVector xx;
467  G4ThreeVector bestxx;
468  G4int i;
469  G4int besti = -1;
470  for (i=0; i< 6; i++)
471  {
472  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
473  if (tmpdistance < distance)
474  {
475  distance = tmpdistance;
476  bestxx = xx;
477  besti = i;
478  }
479  }
480 
481  tmpsurface[0] = surfaces[besti];
482  *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
483 
484  return fLastNormal.vec;
485 }
486 
487 
488 //=====================================================================
489 //* DistanceToIn (p, v) -----------------------------------------------
490 
492  const G4ThreeVector& v ) const
493 {
494 
495  // DistanceToIn (p, v):
496  // Calculate distance to surface of shape from `outside'
497  // along with the v, allowing for tolerance.
498  // The function returns kInfinity if no intersection or
499  // just grazing within tolerance.
500 
501  //
502  // checking last value
503  //
504 
505  G4ThreeVector *tmpp;
506  G4ThreeVector *tmpv;
507  G4double *tmpdist;
508  if (fLastDistanceToInWithV.p == p && fLastDistanceToInWithV.vec == v)
509  {
510  return fLastDistanceToIn.value;
511  }
512  else
513  {
514  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
515  tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
516  tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
517  tmpp->set(p.x(), p.y(), p.z());
518  tmpv->set(v.x(), v.y(), v.z());
519  }
520 
521  //
522  // Calculate DistanceToIn(p,v)
523  //
524 
525  EInside currentside = Inside(p);
526 
527  if (currentside == kInside)
528  {
529  }
530  else if (currentside == kSurface)
531  {
532  // particle is just on a boundary.
533  // if the particle is entering to the volume, return 0
534  //
536  if (normal*v < 0)
537  {
538  *tmpdist = 0;
539  return fLastDistanceToInWithV.value;
540  }
541  }
542 
543  // now, we can take smallest positive distance.
544 
545  // Initialize
546  //
547  G4double distance = kInfinity;
548 
549  // Find intersections and choose nearest one
550  //
551  G4VTwistSurface *surfaces[6];
552 
553  surfaces[0] = fSide0;
554  surfaces[1] = fSide90 ;
555  surfaces[2] = fSide180 ;
556  surfaces[3] = fSide270 ;
557  surfaces[4] = fLowerEndcap;
558  surfaces[5] = fUpperEndcap;
559 
560  G4ThreeVector xx;
561  G4ThreeVector bestxx;
562  G4int i;
563  for (i=0; i < 6 ; i++)
564  {
565 #ifdef G4TWISTDEBUG
566  G4cout << G4endl << "surface " << i << ": " << G4endl << G4endl ;
567 #endif
568  G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
569 #ifdef G4TWISTDEBUG
570  G4cout << "Solid DistanceToIn : distance = " << tmpdistance << G4endl ;
571  G4cout << "intersection point = " << xx << G4endl ;
572 #endif
573  if (tmpdistance < distance)
574  {
575  distance = tmpdistance;
576  bestxx = xx;
577  }
578  }
579 
580 #ifdef G4TWISTDEBUG
581  G4cout << "best distance = " << distance << G4endl ;
582 #endif
583 
584  *tmpdist = distance;
585  // timer.Stop();
586  return fLastDistanceToInWithV.value;
587 }
588 
589 
590 //=====================================================================
591 //* DistanceToIn (p) --------------------------------------------------
592 
594 {
595  // DistanceToIn(p):
596  // Calculate distance to surface of shape from `outside',
597  // allowing for tolerance
598  //
599 
600  //
601  // checking last value
602  //
603 
604  G4ThreeVector *tmpp;
605  G4double *tmpdist;
606  if (fLastDistanceToIn.p == p)
607  {
608  return fLastDistanceToIn.value;
609  }
610  else
611  {
612  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
613  tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
614  tmpp->set(p.x(), p.y(), p.z());
615  }
616 
617  //
618  // Calculate DistanceToIn(p)
619  //
620 
621  EInside currentside = Inside(p);
622 
623  switch (currentside)
624  {
625  case (kInside) :
626  {
627  }
628 
629  case (kSurface) :
630  {
631  *tmpdist = 0.;
632  return fLastDistanceToIn.value;
633  }
634 
635  case (kOutside) :
636  {
637  // Initialize
638  //
639  G4double distance = kInfinity;
640 
641  // Find intersections and choose nearest one
642  //
643  G4VTwistSurface *surfaces[6];
644 
645  surfaces[0] = fSide0;
646  surfaces[1] = fSide90 ;
647  surfaces[2] = fSide180 ;
648  surfaces[3] = fSide270 ;
649  surfaces[4] = fLowerEndcap;
650  surfaces[5] = fUpperEndcap;
651 
652  G4int i;
653  G4ThreeVector xx;
654  G4ThreeVector bestxx;
655  for (i=0; i< 6; i++)
656  {
657  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
658  if (tmpdistance < distance)
659  {
660  distance = tmpdistance;
661  bestxx = xx;
662  }
663  }
664  *tmpdist = distance;
665  return fLastDistanceToIn.value;
666  }
667 
668  default :
669  {
670  G4Exception("G4VTwistedFaceted::DistanceToIn(p)", "GeomSolids0003",
671  FatalException, "Unknown point location!");
672  }
673  } // switch end
674 
675  return 0;
676 }
677 
678 
679 //=====================================================================
680 //* DistanceToOut (p, v) ----------------------------------------------
681 
682 G4double
684  const G4ThreeVector& v,
685  const G4bool calcNorm,
686  G4bool *validNorm,
687  G4ThreeVector *norm ) const
688 {
689  // DistanceToOut (p, v):
690  // Calculate distance to surface of shape from `inside'
691  // along with the v, allowing for tolerance.
692  // The function returns kInfinity if no intersection or
693  // just grazing within tolerance.
694 
695  //
696  // checking last value
697  //
698 
699  G4ThreeVector *tmpp;
700  G4ThreeVector *tmpv;
701  G4double *tmpdist;
702  if (fLastDistanceToOutWithV.p == p && fLastDistanceToOutWithV.vec == v )
703  {
704  return fLastDistanceToOutWithV.value;
705  }
706  else
707  {
708  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
709  tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
710  tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
711  tmpp->set(p.x(), p.y(), p.z());
712  tmpv->set(v.x(), v.y(), v.z());
713  }
714 
715  //
716  // Calculate DistanceToOut(p,v)
717  //
718 
719  EInside currentside = Inside(p);
720 
721  if (currentside == kOutside)
722  {
723  }
724  else if (currentside == kSurface)
725  {
726  // particle is just on a boundary.
727  // if the particle is exiting from the volume, return 0
728  //
730  G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
731  if (normal*v > 0)
732  {
733  if (calcNorm)
734  {
735  *norm = (blockedsurface->GetNormal(p, true));
736  *validNorm = blockedsurface->IsValidNorm();
737  }
738  *tmpdist = 0.;
739  // timer.Stop();
740  return fLastDistanceToOutWithV.value;
741  }
742  }
743 
744  // now, we can take smallest positive distance.
745 
746  // Initialize
747  G4double distance = kInfinity;
748 
749  // find intersections and choose nearest one.
750  G4VTwistSurface *surfaces[6];
751 
752  surfaces[0] = fSide0;
753  surfaces[1] = fSide90 ;
754  surfaces[2] = fSide180 ;
755  surfaces[3] = fSide270 ;
756  surfaces[4] = fLowerEndcap;
757  surfaces[5] = fUpperEndcap;
758 
759  G4int i;
760  G4int besti = -1;
761  G4ThreeVector xx;
762  G4ThreeVector bestxx;
763  for (i=0; i< 6 ; i++) {
764  G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
765  if (tmpdistance < distance)
766  {
767  distance = tmpdistance;
768  bestxx = xx;
769  besti = i;
770  }
771  }
772 
773  if (calcNorm)
774  {
775  if (besti != -1)
776  {
777  *norm = (surfaces[besti]->GetNormal(p, true));
778  *validNorm = surfaces[besti]->IsValidNorm();
779  }
780  }
781 
782  *tmpdist = distance;
783  // timer.Stop();
784  return fLastDistanceToOutWithV.value;
785 }
786 
787 
788 //=====================================================================
789 //* DistanceToOut (p) -------------------------------------------------
790 
792 {
793  // DistanceToOut(p):
794  // Calculate distance to surface of shape from `inside',
795  // allowing for tolerance
796 
797  //
798  // checking last value
799  //
800 
801  G4ThreeVector *tmpp;
802  G4double *tmpdist;
803 
804  if (fLastDistanceToOut.p == p)
805  {
806  return fLastDistanceToOut.value;
807  }
808  else
809  {
810  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
811  tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
812  tmpp->set(p.x(), p.y(), p.z());
813  }
814 
815  //
816  // Calculate DistanceToOut(p)
817  //
818 
819  EInside currentside = Inside(p);
820  G4double retval = kInfinity;
821 
822  switch (currentside)
823  {
824  case (kOutside) :
825  {
826 #ifdef G4SPECSDEBUG
827  G4int oldprc = G4cout.precision(16) ;
828  G4cout << G4endl ;
829  DumpInfo();
830  G4cout << "Position:" << G4endl << G4endl ;
831  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
832  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
833  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
834  G4cout.precision(oldprc) ;
835  G4Exception("G4VTwistedFaceted::DistanceToOut(p)", "GeomSolids1002",
836  JustWarning, "Point p is outside !?" );
837 #endif
838  break;
839  }
840  case (kSurface) :
841  {
842  *tmpdist = 0.;
843  retval = fLastDistanceToOut.value;
844  break;
845  }
846 
847  case (kInside) :
848  {
849  // Initialize
850  //
851  G4double distance = kInfinity;
852 
853  // find intersections and choose nearest one
854  //
855  G4VTwistSurface *surfaces[6];
856 
857  surfaces[0] = fSide0;
858  surfaces[1] = fSide90 ;
859  surfaces[2] = fSide180 ;
860  surfaces[3] = fSide270 ;
861  surfaces[4] = fLowerEndcap;
862  surfaces[5] = fUpperEndcap;
863 
864  G4int i;
865  G4ThreeVector xx;
866  G4ThreeVector bestxx;
867  for (i=0; i< 6; i++)
868  {
869  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
870  if (tmpdistance < distance)
871  {
872  distance = tmpdistance;
873  bestxx = xx;
874  }
875  }
876  *tmpdist = distance;
877 
878  retval = fLastDistanceToOut.value;
879  break;
880  }
881 
882  default :
883  {
884  G4Exception("G4VTwistedFaceted::DistanceToOut(p)", "GeomSolids0003",
885  FatalException, "Unknown point location!");
886  break;
887  }
888  } // switch end
889 
890  return retval;
891 }
892 
893 
894 //=====================================================================
895 //* StreamInfo --------------------------------------------------------
896 
897 std::ostream& G4VTwistedFaceted::StreamInfo(std::ostream& os) const
898 {
899  //
900  // Stream object contents to an output stream
901  //
902  G4int oldprc = os.precision(16);
903  os << "-----------------------------------------------------------\n"
904  << " *** Dump for solid - " << GetName() << " ***\n"
905  << " ===================================================\n"
906  << " Solid type: G4VTwistedFaceted\n"
907  << " Parameters: \n"
908  << " polar angle theta = " << fTheta/degree << " deg" << G4endl
909  << " azimuthal angle phi = " << fPhi/degree << " deg" << G4endl
910  << " tilt angle alpha = " << fAlph/degree << " deg" << G4endl
911  << " TWIST angle = " << fPhiTwist/degree << " deg" << G4endl
912  << " Half length along y (lower endcap) = " << fDy1/cm << " cm"
913  << G4endl
914  << " Half length along x (lower endcap, bottom) = " << fDx1/cm << " cm"
915  << G4endl
916  << " Half length along x (lower endcap, top) = " << fDx2/cm << " cm"
917  << G4endl
918  << " Half length along y (upper endcap) = " << fDy2/cm << " cm"
919  << G4endl
920  << " Half length along x (upper endcap, bottom) = " << fDx3/cm << " cm"
921  << G4endl
922  << " Half length along x (upper endcap, top) = " << fDx4/cm << " cm"
923  << G4endl
924  << "-----------------------------------------------------------\n";
925  os.precision(oldprc);
926 
927  return os;
928 }
929 
930 
931 //=====================================================================
932 //* DiscribeYourselfTo ------------------------------------------------
933 
935 {
936  scene.AddSolid (*this);
937 }
938 
939 
940 //=====================================================================
941 //* GetExtent ---------------------------------------------------------
942 
944 {
945  G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
946 
947  return G4VisExtent(-maxRad, maxRad ,
948  -maxRad, maxRad ,
949  -fDz, fDz );
950 }
951 
952 
953 //=====================================================================
954 //* CreateSurfaces ----------------------------------------------------
955 
956 void G4VTwistedFaceted::CreateSurfaces()
957 {
958 
959  // create 6 surfaces of TwistedTub.
960 
961  if ( fDx1 == fDx2 && fDx3 == fDx4 ) // special case : Box
962  {
963  fSide0 = new G4TwistBoxSide("0deg", fPhiTwist, fDz, fTheta, fPhi,
964  fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 0.*deg);
965  fSide180 = new G4TwistBoxSide("180deg", fPhiTwist, fDz, fTheta, fPhi+pi,
966  fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 180.*deg);
967  }
968  else // default general case
969  {
970  fSide0 = new G4TwistTrapAlphaSide("0deg" ,fPhiTwist, fDz, fTheta,
971  fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
972  fSide180 = new G4TwistTrapAlphaSide("180deg", fPhiTwist, fDz, fTheta,
973  fPhi+pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
974  }
975 
976  // create parallel sides
977  //
978  fSide90 = new G4TwistTrapParallelSide("90deg", fPhiTwist, fDz, fTheta,
979  fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
980  fSide270 = new G4TwistTrapParallelSide("270deg", fPhiTwist, fDz, fTheta,
981  fPhi+pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
982 
983  // create endcaps
984  //
985  fUpperEndcap = new G4TwistTrapFlatSide("UpperCap",fPhiTwist, fDx3, fDx4, fDy2,
986  fDz, fAlph, fPhi, fTheta, 1 );
987  fLowerEndcap = new G4TwistTrapFlatSide("LowerCap",fPhiTwist, fDx1, fDx2, fDy1,
988  fDz, fAlph, fPhi, fTheta, -1 );
989 
990  // Set neighbour surfaces
991 
992  fSide0->SetNeighbours( fSide270 , fLowerEndcap , fSide90 , fUpperEndcap );
993  fSide90->SetNeighbours( fSide0 , fLowerEndcap , fSide180 , fUpperEndcap );
994  fSide180->SetNeighbours(fSide90 , fLowerEndcap , fSide270 , fUpperEndcap );
995  fSide270->SetNeighbours(fSide180 , fLowerEndcap , fSide0 , fUpperEndcap );
996  fUpperEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
997  fLowerEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
998 
999 }
1000 
1001 
1002 //=====================================================================
1003 //* GetEntityType -----------------------------------------------------
1004 
1006 {
1007  return G4String("G4VTwistedFaceted");
1008 }
1009 
1010 
1011 //=====================================================================
1012 //* GetPolyhedron -----------------------------------------------------
1013 
1015 {
1016  if (!fpPolyhedron ||
1020  {
1021  G4AutoLock l(&polyhedronMutex);
1022  delete fpPolyhedron;
1024  fRebuildPolyhedron = false;
1025  l.unlock();
1026  }
1027 
1028  return fpPolyhedron;
1029 }
1030 
1031 
1032 //=====================================================================
1033 //* GetPointInSolid ---------------------------------------------------
1034 
1036 {
1037 
1038 
1039  // this routine is only used for a test
1040  // can be deleted ...
1041 
1042  if ( z == fDz ) z -= 0.1*fDz ;
1043  if ( z == -fDz ) z += 0.1*fDz ;
1044 
1045  G4double phi = z/(2*fDz)*fPhiTwist ;
1046 
1047  return G4ThreeVector(fdeltaX * phi/fPhiTwist, fdeltaY * phi/fPhiTwist, z ) ;
1048 }
1049 
1050 
1051 //=====================================================================
1052 //* GetPointOnSurface -------------------------------------------------
1053 
1055 {
1056 
1057  G4double phi = G4RandFlat::shoot(-fPhiTwist/2.,fPhiTwist/2.);
1058  G4double u , umin, umax ; // variable for twisted surfaces
1059  G4double y ; // variable for flat surface (top and bottom)
1060 
1061  // Compute the areas. Attention: Only correct for trapezoids
1062  // where the twisting is done along the z-axis. In the general case
1063  // the computed surface area is more difficult. However this simplification
1064  // does not affect the tracking through the solid.
1065 
1066  G4double a1 = fSide0->GetSurfaceArea();
1067  G4double a2 = fSide90->GetSurfaceArea();
1068  G4double a3 = fSide180->GetSurfaceArea() ;
1069  G4double a4 = fSide270->GetSurfaceArea() ;
1070  G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1071  G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1072 
1073 #ifdef G4TWISTDEBUG
1074  G4cout << "Surface 0 deg = " << a1 << G4endl ;
1075  G4cout << "Surface 90 deg = " << a2 << G4endl ;
1076  G4cout << "Surface 180 deg = " << a3 << G4endl ;
1077  G4cout << "Surface 270 deg = " << a4 << G4endl ;
1078  G4cout << "Surface Lower = " << a5 << G4endl ;
1079  G4cout << "Surface Upper = " << a6 << G4endl ;
1080 #endif
1081 
1082  G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1083 
1084  if(chose < a1)
1085  {
1086 
1087  umin = fSide0->GetBoundaryMin(phi) ;
1088  umax = fSide0->GetBoundaryMax(phi) ;
1089  u = G4RandFlat::shoot(umin,umax) ;
1090 
1091  return fSide0->SurfacePoint(phi, u, true) ; // point on 0deg surface
1092  }
1093 
1094  else if( (chose >= a1) && (chose < a1 + a2 ) )
1095  {
1096 
1097  umin = fSide90->GetBoundaryMin(phi) ;
1098  umax = fSide90->GetBoundaryMax(phi) ;
1099 
1100  u = G4RandFlat::shoot(umin,umax) ;
1101 
1102  return fSide90->SurfacePoint(phi, u, true); // point on 90deg surface
1103  }
1104 
1105  else if( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1106  {
1107 
1108  umin = fSide180->GetBoundaryMin(phi) ;
1109  umax = fSide180->GetBoundaryMax(phi) ;
1110  u = G4RandFlat::shoot(umin,umax) ;
1111 
1112  return fSide180->SurfacePoint(phi, u, true); // point on 180 deg surface
1113  }
1114 
1115  else if( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1116  {
1117 
1118  umin = fSide270->GetBoundaryMin(phi) ;
1119  umax = fSide270->GetBoundaryMax(phi) ;
1120  u = G4RandFlat::shoot(umin,umax) ;
1121 
1122  return fSide270->SurfacePoint(phi, u, true); // point on 270 deg surface
1123  }
1124 
1125  else if( (chose >= a1 + a2 + a3 + a4 ) && (chose < a1 + a2 + a3 + a4 + a5 ) )
1126  {
1127 
1128  y = G4RandFlat::shoot(-fDy1,fDy1) ;
1129  umin = fLowerEndcap->GetBoundaryMin(y) ;
1130  umax = fLowerEndcap->GetBoundaryMax(y) ;
1131  u = G4RandFlat::shoot(umin,umax) ;
1132 
1133  return fLowerEndcap->SurfacePoint(u,y,true); // point on lower endcap
1134  }
1135  else {
1136 
1137  y = G4RandFlat::shoot(-fDy2,fDy2) ;
1138  umin = fUpperEndcap->GetBoundaryMin(y) ;
1139  umax = fUpperEndcap->GetBoundaryMax(y) ;
1140  u = G4RandFlat::shoot(umin,umax) ;
1141 
1142  return fUpperEndcap->SurfacePoint(u,y,true) ; // point on upper endcap
1143 
1144  }
1145 }
1146 
1147 
1148 //=====================================================================
1149 //* CreatePolyhedron --------------------------------------------------
1150 
1152 {
1153  // number of meshes
1154  const G4int k =
1155  G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2;
1156  const G4int n = k;
1157 
1158  const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
1159  const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
1160 
1161  G4Polyhedron *ph=new G4Polyhedron;
1162  typedef G4double G4double3[3];
1163  typedef G4int G4int4[4];
1164  G4double3* xyz = new G4double3[nnodes]; // number of nodes
1165  G4int4* faces = new G4int4[nfaces] ; // number of faces
1166 
1167  fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
1168  fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
1169  fSide270->GetFacets(k,n,xyz,faces,2) ;
1170  fSide0->GetFacets(k,n,xyz,faces,3) ;
1171  fSide90->GetFacets(k,n,xyz,faces,4) ;
1172  fSide180->GetFacets(k,n,xyz,faces,5) ;
1173 
1174  ph->createPolyhedron(nnodes,nfaces,xyz,faces);
1175 
1176  return ph;
1177 }
void set(double x, double y, double z)
G4String GetName() const
ThreeVector shoot(const G4int Ap, const G4int Af)
virtual void ComputeDimensions(G4VPVParameterisation *, const G4int, const G4VPhysicalVolume *)
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
G4VTwistedFaceted(const G4String &pname, G4double PhiTwist, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlph)
CLHEP::Hep3Vector G4ThreeVector
double x() const
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal)=0
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
virtual G4GeometryType GetEntityType() const
const char * p
Definition: xmltok.h:285
G4ThreeVector GetPointOnSurface() const
G4ThreeVector GetPointInSolid(G4double z) const
virtual G4Polyhedron * CreatePolyhedron() const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
virtual void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4VTwistedFaceted & operator=(const G4VTwistedFaceted &rhs)
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
G4bool IsValidNorm() const
double z() const
void DumpInfo() const
virtual G4double DistanceToIn(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double Xcoef(G4double u, G4double phi, G4double ftg) const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
virtual G4double GetSurfaceArea()=0
G4GLOB_DLL std::ostream G4cout
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
void SetNeighbours(G4VTwistSurface *axis0min, G4VTwistSurface *axis1min, G4VTwistSurface *axis0max, G4VTwistSurface *axis1max)
virtual G4double DistanceTo(const G4ThreeVector &gp, G4ThreeVector &gxx)
static constexpr double degree
Definition: G4SIunits.hh:144
virtual G4double DistanceToOut(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)=0
bool G4bool
Definition: G4Types.hh:79
static constexpr double cm
Definition: G4SIunits.hh:119
G4double GetValueB(G4double phi) const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
virtual G4double GetBoundaryMax(G4double)=0
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)=0
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcnorm=false, G4bool *validnorm=0, G4ThreeVector *n=0) const
virtual G4double GetBoundaryMin(G4double)=0
G4int G4Mutex
Definition: G4Threading.hh:173
virtual G4Polyhedron * GetPolyhedron() const
static G4int GetNumberOfRotationSteps()
EInside
Definition: geomdefs.hh:58
virtual G4VisExtent GetExtent() const
virtual EInside Inside(const G4ThreeVector &p) const
EAxis
Definition: geomdefs.hh:54
double y() const
virtual void DescribeYourselfTo(G4VGraphicsScene &scene) const
#define G4endl
Definition: G4ios.hh:61
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:111
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
static constexpr double deg
Definition: G4SIunits.hh:152
virtual std::ostream & StreamInfo(std::ostream &os) const
G4Polyhedron * fpPolyhedron
G4double GetAngularTolerance() const
static G4GeometryTolerance * GetInstance()