Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Para.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4Para.cc 101121 2016-11-07 09:18:01Z gcosmo $
28 //
29 // class G4Para
30 //
31 // Implementation for G4Para class
32 //
33 // History:
34 //
35 // 23.09.16 E.Tcherniaev: added Extent(pmin,pmax),
36 // use G4BoundingEnvelope for CalculateExtent(),
37 // removed CreateRotatedVertices()
38 // 23.10.05 V.Grichine: bug fixed in DistanceToOut(p,v,...) for the v.x()<0 case
39 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal
40 // 30.11.04 V.Grichine: modifications in SurfaceNormal for edges/vertices and
41 // in constructor with vertices
42 // 14.02.02 V.Grichine: bug fixed in Inside according to proposal of D.Wright
43 // 18.11.99 V.Grichine: kUndef was added to ESide
44 // 31.10.96 V.Grichine: Modifications according G4Box/Tubs before to commit
45 // 21.03.95 P.Kent: Modified for `tolerant' geom
46 //
48 
49 #include "G4Para.hh"
50 
51 #include "G4VoxelLimits.hh"
52 #include "G4AffineTransform.hh"
53 #include "G4BoundingEnvelope.hh"
54 #include "Randomize.hh"
55 
56 #include "G4VPVParameterisation.hh"
57 
58 #include "G4VGraphicsScene.hh"
59 
60 using namespace CLHEP;
61 
62 // Private enum: Not for external use
63 
65 
66 // used internally for normal routine
67 
68 enum ENSide {kNZ,kNX,kNY};
69 
71 //
72 // Constructor - check and set half-widths
73 
75  G4double pAlpha, G4double pTheta, G4double pPhi )
76 {
77  if ( pDx > 0 && pDy > 0 && pDz > 0 )
78  {
79  fDx = pDx;
80  fDy = pDy;
81  fDz = pDz;
82  fTalpha = std::tan(pAlpha);
83  fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
84  fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
85  }
86  else
87  {
88  std::ostringstream message;
89  message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
90  << " pDx, pDy, pDz = "
91  << pDx << ", " << pDy << ", " << pDz;
92  G4Exception("G4Para::SetAllParameters()", "GeomSolids0002",
93  FatalException, message);
94  }
95  fCubicVolume = 0.;
96  fSurfaceArea = 0.;
97 }
98 
100 //
101 
103  G4double pDx, G4double pDy, G4double pDz,
104  G4double pAlpha, G4double pTheta, G4double pPhi)
105  : G4CSGSolid(pName)
106 {
107  if ((pDx<=0) || (pDy<=0) || (pDz<=0))
108  {
109  std::ostringstream message;
110  message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
111  << " pDx, pDy, pDz = "
112  << pDx << ", " << pDy << ", " << pDz;
113  G4Exception("G4Para::G4Para()", "GeomSolids0002",
114  FatalException, message);
115  }
116  SetAllParameters( pDx, pDy, pDz, pAlpha, pTheta, pPhi);
117 }
118 
120 //
121 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
122 // which are its vertices. Checking of planarity with preparation of
123 // fPlanes[] and than calculation of other members
124 
125 G4Para::G4Para( const G4String& pName,
126  const G4ThreeVector pt[8] )
127  : G4CSGSolid(pName)
128 {
129  if (!( pt[0].z()<0 && pt[0].z()==pt[1].z() && pt[0].z()==pt[2].z() &&
130  pt[0].z()==pt[3].z() && pt[4].z()>0 && pt[4].z()==pt[5].z() &&
131  pt[4].z()==pt[6].z() && pt[4].z()==pt[7].z() &&
132  (pt[0].z()+pt[4].z())==0 &&
133  pt[0].y()==pt[1].y() && pt[2].y()==pt[3].y() &&
134  pt[4].y()==pt[5].y() && pt[6].y()==pt[7].y() &&
135  ( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) == 0 &&
136  ( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() ) == 0) )
137  {
138  std::ostringstream message;
139  message << "Invalid vertice coordinates for Solid: " << GetName();
140  G4Exception("G4Para::G4Para()", "GeomSolids0002",
141  FatalException, message);
142  }
143  fDx = ((pt[3]).x()-(pt[2]).x())*0.5;
144  fDy = ((pt[2]).y()-(pt[1]).y())*0.5;
145  fDz = (pt[7]).z();
146  fTalpha = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy ;
147  fTthetaCphi = ((pt[4]).x()+fDy*fTalpha+fDx)/fDz ;
148  fTthetaSphi = ((pt[4]).y()+fDy)/fDz ;
149  fCubicVolume = 0.;
150  fSurfaceArea = 0.;
151 }
152 
154 //
155 // Fake default constructor - sets only member data and allocates memory
156 // for usage restricted to object persistency.
157 //
158 G4Para::G4Para( __void__& a )
159  : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.),
160  fTalpha(0.), fTthetaCphi(0.), fTthetaSphi(0.)
161 {
162 }
163 
165 //
166 
168 {
169 }
170 
172 //
173 // Copy constructor
174 
176  : G4CSGSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),
177  fTalpha(rhs.fTalpha), fTthetaCphi(rhs.fTthetaCphi),
178  fTthetaSphi(rhs.fTthetaSphi)
179 {
180 }
181 
183 //
184 // Assignment operator
185 
187 {
188  // Check assignment to self
189  //
190  if (this == &rhs) { return *this; }
191 
192  // Copy base class data
193  //
195 
196  // Copy data
197  //
198  fDx = rhs.fDx; fDy = rhs.fDy; fDz = rhs.fDz;
199  fTalpha = rhs.fTalpha; fTthetaCphi = rhs.fTthetaCphi;
200  fTthetaSphi = rhs.fTthetaSphi;
201 
202  return *this;
203 }
204 
206 //
207 // Dispatch to parameterisation for replication mechanism dimension
208 // computation & modification.
209 
211  const G4int n,
212  const G4VPhysicalVolume* pRep )
213 {
214  p->ComputeDimensions(*this,n,pRep);
215 }
216 
218 //
219 // Get bounding box
220 
222 {
223  G4double dz = GetZHalfLength();
224  G4double dx = GetXHalfLength();
225  G4double dy = GetYHalfLength();
226 
227  G4double x0 = dz*fTthetaCphi;
228  G4double x1 = dy*GetTanAlpha();
229  G4double xmin =
230  std::min(
231  std::min(
232  std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0+x1-dx);
233  G4double xmax =
234  std::max(
235  std::max(
236  std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0+x1+dx);
237 
238  G4double y0 = dz*fTthetaSphi;
239  G4double ymin = std::min(-y0-dy,y0-dy);
240  G4double ymax = std::max(-y0+dy,y0+dy);
241 
242  pMin.set(xmin,ymin,-dz);
243  pMax.set(xmax,ymax, dz);
244 
245  // Check correctness of the bounding box
246  //
247  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
248  {
249  std::ostringstream message;
250  message << "Bad bounding box (min >= max) for solid: "
251  << GetName() << " !"
252  << "\npMin = " << pMin
253  << "\npMax = " << pMax;
254  G4Exception("G4Para::Extent()", "GeomMgt0001", JustWarning, message);
255  DumpInfo();
256  }
257 }
258 
260 //
261 // Calculate extent under transform and specified limit
262 
264  const G4VoxelLimits& pVoxelLimit,
265  const G4AffineTransform& pTransform,
266  G4double& pMin, G4double& pMax ) const
267 {
268  G4ThreeVector bmin, bmax;
269  G4bool exist;
270 
271  // Check bounding box (bbox)
272  //
273  Extent(bmin,bmax);
274  G4BoundingEnvelope bbox(bmin,bmax);
275 #ifdef G4BBOX_EXTENT
276  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
277 #endif
278  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
279  {
280  return exist = (pMin < pMax) ? true : false;
281  }
282 
283  // Set bounding envelope (benv) and calculate extent
284  //
285  G4double dz = GetZHalfLength();
286  G4double dx = GetXHalfLength();
287  G4double dy = GetYHalfLength();
288 
289  G4double x0 = dz*fTthetaCphi;
290  G4double x1 = dy*GetTanAlpha();
291  G4double y0 = dz*fTthetaSphi;
292 
293  G4ThreeVectorList baseA(4), baseB(4);
294  baseA[0].set(-x0-x1-dx,-y0-dy,-dz);
295  baseA[1].set(-x0-x1+dx,-y0-dy,-dz);
296  baseA[2].set(-x0+x1+dx,-y0+dy,-dz);
297  baseA[3].set(-x0+x1-dx,-y0+dy,-dz);
298 
299  baseB[0].set(+x0-x1-dx, y0-dy, dz);
300  baseB[1].set(+x0-x1+dx, y0-dy, dz);
301  baseB[2].set(+x0+x1+dx, y0+dy, dz);
302  baseB[3].set(+x0+x1-dx, y0+dy, dz);
303 
304  std::vector<const G4ThreeVectorList *> polygons(2);
305  polygons[0] = &baseA;
306  polygons[1] = &baseB;
307 
308  G4BoundingEnvelope benv(bmin,bmax,polygons);
309  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
310  return exist;
311 }
312 
314 //
315 // Check in p is inside/on surface/outside solid
316 
318 {
319  G4double xt, yt, yt1;
320  EInside in = kOutside;
321 
322  yt1 = p.y() - fTthetaSphi*p.z();
323  yt = std::fabs(yt1) ;
324 
325  // xt = std::fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt );
326 
327  xt = std::fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt1 );
328 
329  if ( std::fabs( p.z() ) <= fDz - kCarTolerance*0.5)
330  {
331  if (yt <= fDy - kCarTolerance*0.5)
332  {
333  if ( xt <= fDx - kCarTolerance*0.5 ) in = kInside;
334  else if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
335  }
336  else if ( yt <= fDy + kCarTolerance*0.5)
337  {
338  if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
339  }
340  }
341  else if ( std::fabs(p.z()) <= fDz + kCarTolerance*0.5 )
342  {
343  if ( yt <= fDy + kCarTolerance*0.5)
344  {
345  if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
346  }
347  }
348  return in;
349 }
350 
352 //
353 // Calculate side nearest to p, and return normal
354 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
355 
357 {
358  G4ThreeVector norm, sumnorm(0.,0.,0.);
359  G4int noSurfaces = 0;
360  G4double distx,disty,distz;
361  G4double newpx,newpy,xshift;
362  G4double calpha,salpha; // Sin/Cos(alpha) - needed to recalc G4Parameter
363  G4double tntheta,cosntheta; // tan and cos of normal's theta component
364  G4double ycomp;
365  G4double delta = 0.5*kCarTolerance;
366 
367  newpx = p.x()-fTthetaCphi*p.z();
368  newpy = p.y()-fTthetaSphi*p.z();
369 
370  calpha = 1/std::sqrt(1+fTalpha*fTalpha);
371  salpha = -calpha*fTalpha; // NOTE: using MINUS std::sin(alpha)
372 
373  // xshift = newpx*calpha+newpy*salpha;
374  xshift = newpx - newpy*fTalpha;
375 
376  // distx = std::fabs(std::fabs(xshift)-fDx*calpha);
377  distx = std::fabs(std::fabs(xshift)-fDx);
378  disty = std::fabs(std::fabs(newpy)-fDy);
379  distz = std::fabs(std::fabs(p.z())-fDz);
380 
381  tntheta = fTthetaCphi*calpha + fTthetaSphi*salpha;
382  cosntheta = 1/std::sqrt(1+tntheta*tntheta);
383  ycomp = 1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
384 
385  G4ThreeVector nX = G4ThreeVector( calpha*cosntheta,
386  salpha*cosntheta,
387  -tntheta*cosntheta);
388  G4ThreeVector nY = G4ThreeVector( 0, ycomp,-fTthetaSphi*ycomp);
389  G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0);
390 
391  if (distx <= delta)
392  {
393  noSurfaces ++;
394  if ( xshift >= 0.) {sumnorm += nX;}
395  else {sumnorm -= nX;}
396  }
397  if (disty <= delta)
398  {
399  noSurfaces ++;
400  if ( newpy >= 0.) {sumnorm += nY;}
401  else {sumnorm -= nY;}
402  }
403  if (distz <= delta)
404  {
405  noSurfaces ++;
406  if ( p.z() >= 0.) {sumnorm += nZ;}
407  else {sumnorm -= nZ;}
408  }
409  if ( noSurfaces == 0 )
410  {
411 #ifdef G4CSGDEBUG
412  G4Exception("G4Para::SurfaceNormal(p)", "GeomSolids1002",
413  JustWarning, "Point p is not on surface !?" );
414 #endif
415  norm = ApproxSurfaceNormal(p);
416  }
417  else if ( noSurfaces == 1 ) {norm = sumnorm;}
418  else {norm = sumnorm.unit();}
419 
420  return norm;
421 }
422 
423 
425 //
426 // Algorithm for SurfaceNormal() following the original specification
427 // for points not on the surface
428 
429 G4ThreeVector G4Para::ApproxSurfaceNormal( const G4ThreeVector& p ) const
430 {
431  ENSide side;
432  G4ThreeVector norm;
433  G4double distx,disty,distz;
434  G4double newpx,newpy,xshift;
435  G4double calpha,salpha; // Sin/Cos(alpha) - needed to recalc G4Parameter
436  G4double tntheta,cosntheta; // tan and cos of normal's theta component
437  G4double ycomp;
438 
439  newpx=p.x()-fTthetaCphi*p.z();
440  newpy=p.y()-fTthetaSphi*p.z();
441 
442  calpha=1/std::sqrt(1+fTalpha*fTalpha);
443  salpha=-calpha*fTalpha; // NOTE: actually use MINUS std::sin(alpha)
444 
445  xshift=newpx*calpha+newpy*salpha;
446 
447  distx=std::fabs(std::fabs(xshift)-fDx*calpha);
448  disty=std::fabs(std::fabs(newpy)-fDy);
449  distz=std::fabs(std::fabs(p.z())-fDz);
450 
451  if (distx<disty)
452  {
453  if (distx<distz) {side=kNX;}
454  else {side=kNZ;}
455  }
456  else
457  {
458  if (disty<distz) {side=kNY;}
459  else {side=kNZ;}
460  }
461 
462  switch (side)
463  {
464  case kNX:
465  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
466  if (xshift<0)
467  {
468  cosntheta=-1/std::sqrt(1+tntheta*tntheta);
469  }
470  else
471  {
472  cosntheta=1/std::sqrt(1+tntheta*tntheta);
473  }
474  norm=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
475  break;
476  case kNY:
477  if (newpy<0)
478  {
479  ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
480  }
481  else
482  {
483  ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
484  }
485  norm=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
486  break;
487  case kNZ: // Closest to Z
488  if (p.z()>=0)
489  {
490  norm=G4ThreeVector(0,0,1);
491  }
492  else
493  {
494  norm=G4ThreeVector(0,0,-1);
495  }
496  break;
497  }
498  return norm;
499 }
500 
502 //
503 // Calculate distance to shape from outside
504 // - return kInfinity if no intersection
505 //
506 // ALGORITHM:
507 // For each component, calculate pair of minimum and maximum intersection
508 // values for which the particle is in the extent of the shape
509 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
510 // - Z plane intersectin uses tolerance
511 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
512 // (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
513 // cases)
514 // - Note: XZ and YZ planes each divide space into four regions,
515 // characterised by ss1 ss2
516 
518  const G4ThreeVector& v ) const
519 {
520  G4double snxt; // snxt = default return value
521  G4double smin,smax;
522  G4double tmin,tmax;
523  G4double yt,vy,xt,vx;
524  G4double max;
525  //
526  // Z Intersection range
527  //
528  if (v.z()>0)
529  {
530  max=fDz-p.z();
531  if (max>kCarTolerance*0.5)
532  {
533  smax=max/v.z();
534  smin=(-fDz-p.z())/v.z();
535  }
536  else
537  {
538  return snxt=kInfinity;
539  }
540  }
541  else if (v.z()<0)
542  {
543  max=-fDz-p.z();
544  if (max<-kCarTolerance*0.5)
545  {
546  smax=max/v.z();
547  smin=(fDz-p.z())/v.z();
548  }
549  else
550  {
551  return snxt=kInfinity;
552  }
553  }
554  else
555  {
556  if (std::fabs(p.z())<=fDz) // Inside
557  {
558  smin=0;
559  smax=kInfinity;
560  }
561  else
562  {
563  return snxt=kInfinity;
564  }
565  }
566 
567  //
568  // Y G4Parallel planes intersection
569  //
570 
571  yt=p.y()-fTthetaSphi*p.z();
572  vy=v.y()-fTthetaSphi*v.z();
573 
574  if (vy>0)
575  {
576  max=fDy-yt;
577  if (max>kCarTolerance*0.5)
578  {
579  tmax=max/vy;
580  tmin=(-fDy-yt)/vy;
581  }
582  else
583  {
584  return snxt=kInfinity;
585  }
586  }
587  else if (vy<0)
588  {
589  max=-fDy-yt;
590  if (max<-kCarTolerance*0.5)
591  {
592  tmax=max/vy;
593  tmin=(fDy-yt)/vy;
594  }
595  else
596  {
597  return snxt=kInfinity;
598  }
599  }
600  else
601  {
602  if (std::fabs(yt)<=fDy)
603  {
604  tmin=0;
605  tmax=kInfinity;
606  }
607  else
608  {
609  return snxt=kInfinity;
610  }
611  }
612 
613  // Re-Calc valid intersection range
614  //
615  if (tmin>smin) smin=tmin;
616  if (tmax<smax) smax=tmax;
617  if (smax<=smin)
618  {
619  return snxt=kInfinity;
620  }
621  else
622  {
623  //
624  // X G4Parallel planes intersection
625  //
626  xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
627  vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
628  if (vx>0)
629  {
630  max=fDx-xt;
631  if (max>kCarTolerance*0.5)
632  {
633  tmax=max/vx;
634  tmin=(-fDx-xt)/vx;
635  }
636  else
637  {
638  return snxt=kInfinity;
639  }
640  }
641  else if (vx<0)
642  {
643  max=-fDx-xt;
644  if (max<-kCarTolerance*0.5)
645  {
646  tmax=max/vx;
647  tmin=(fDx-xt)/vx;
648  }
649  else
650  {
651  return snxt=kInfinity;
652  }
653  }
654  else
655  {
656  if (std::fabs(xt)<=fDx)
657  {
658  tmin=0;
659  tmax=kInfinity;
660  }
661  else
662  {
663  return snxt=kInfinity;
664  }
665  }
666  if (tmin>smin) smin=tmin;
667  if (tmax<smax) smax=tmax;
668  }
669 
670  if (smax>0&&smin<smax)
671  {
672  if (smin>0)
673  {
674  snxt=smin;
675  }
676  else
677  {
678  snxt=0;
679  }
680  }
681  else
682  {
683  snxt=kInfinity;
684  }
685  return snxt;
686 }
687 
689 //
690 // Calculate exact shortest distance to any boundary from outside
691 // - Returns 0 is point inside
692 
694 {
695  G4double safe=0.0;
696  G4double distz1,distz2,disty1,disty2,distx1,distx2;
697  G4double trany,cosy,tranx,cosx;
698 
699  // Z planes
700  //
701  distz1=p.z()-fDz;
702  distz2=-fDz-p.z();
703  if (distz1>distz2)
704  {
705  safe=distz1;
706  }
707  else
708  {
709  safe=distz2;
710  }
711 
712  trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
713 
714  // Transformed x into `box' system
715  //
716  cosy=1.0/std::sqrt(1.0+fTthetaSphi*fTthetaSphi);
717  disty1=(trany-fDy)*cosy;
718  disty2=(-fDy-trany)*cosy;
719 
720  if (disty1>safe) safe=disty1;
721  if (disty2>safe) safe=disty2;
722 
723  tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
724  cosx=1.0/std::sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
725  distx1=(tranx-fDx)*cosx;
726  distx2=(-fDx-tranx)*cosx;
727 
728  if (distx1>safe) safe=distx1;
729  if (distx2>safe) safe=distx2;
730 
731  if (safe<0) safe=0;
732  return safe;
733 }
734 
736 //
737 // Calculate distance to surface of shape from inside
738 // Calculate distance to x/y/z planes - smallest is exiting distance
739 
741  const G4bool calcNorm,
742  G4bool *validNorm, G4ThreeVector *n) const
743 {
744  ESide side = kUndef;
745  G4double snxt; // snxt = return value
746  G4double max,tmax;
747  G4double yt,vy,xt,vx;
748 
749  G4double ycomp,calpha,salpha,tntheta,cosntheta;
750 
751  //
752  // Z Intersections
753  //
754 
755  if (v.z()>0)
756  {
757  max=fDz-p.z();
758  if (max>kCarTolerance*0.5)
759  {
760  snxt=max/v.z();
761  side=kPZ;
762  }
763  else
764  {
765  if (calcNorm)
766  {
767  *validNorm=true;
768  *n=G4ThreeVector(0,0,1);
769  }
770  return snxt=0;
771  }
772  }
773  else if (v.z()<0)
774  {
775  max=-fDz-p.z();
776  if (max<-kCarTolerance*0.5)
777  {
778  snxt=max/v.z();
779  side=kMZ;
780  }
781  else
782  {
783  if (calcNorm)
784  {
785  *validNorm=true;
786  *n=G4ThreeVector(0,0,-1);
787  }
788  return snxt=0;
789  }
790  }
791  else
792  {
793  snxt=kInfinity;
794  }
795 
796  //
797  // Y plane intersection
798  //
799 
800  yt=p.y()-fTthetaSphi*p.z();
801  vy=v.y()-fTthetaSphi*v.z();
802 
803  if (vy>0)
804  {
805  max=fDy-yt;
806  if (max>kCarTolerance*0.5)
807  {
808  tmax=max/vy;
809  if (tmax<snxt)
810  {
811  snxt=tmax;
812  side=kPY;
813  }
814  }
815  else
816  {
817  if (calcNorm)
818  {
819  *validNorm=true; // Leaving via plus Y
820  ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
821  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
822  }
823  return snxt=0;
824  }
825  }
826  else if (vy<0)
827  {
828  max=-fDy-yt;
829  if (max<-kCarTolerance*0.5)
830  {
831  tmax=max/vy;
832  if (tmax<snxt)
833  {
834  snxt=tmax;
835  side=kMY;
836  }
837  }
838  else
839  {
840  if (calcNorm)
841  {
842  *validNorm=true; // Leaving via minus Y
843  ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
844  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
845  }
846  return snxt=0;
847  }
848  }
849 
850  //
851  // X plane intersection
852  //
853 
854  xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
855  vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
856  if (vx>0)
857  {
858  max=fDx-xt;
859  if (max>kCarTolerance*0.5)
860  {
861  tmax=max/vx;
862  if (tmax<snxt)
863  {
864  snxt=tmax;
865  side=kPX;
866  }
867  }
868  else
869  {
870  if (calcNorm)
871  {
872  *validNorm=true; // Leaving via plus X
873  calpha=1/std::sqrt(1+fTalpha*fTalpha);
874  salpha=-calpha*fTalpha; // NOTE: actually use MINUS std::sin(alpha)
875  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
876  cosntheta=1/std::sqrt(1+tntheta*tntheta);
877  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
878  }
879  return snxt=0;
880  }
881  }
882  else if (vx<0)
883  {
884  max=-fDx-xt;
885  if (max<-kCarTolerance*0.5)
886  {
887  tmax=max/vx;
888  if (tmax<snxt)
889  {
890  snxt=tmax;
891  side=kMX;
892  }
893  }
894  else
895  {
896  if (calcNorm)
897  {
898  *validNorm=true; // Leaving via minus X
899  calpha=1/std::sqrt(1+fTalpha*fTalpha);
900  salpha=-calpha*fTalpha; // NOTE: actually use MINUS std::sin(alpha)
901  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
902  cosntheta=-1/std::sqrt(1+tntheta*tntheta);
903  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
904  }
905  return snxt=0;
906  }
907  }
908 
909  if (calcNorm)
910  {
911  *validNorm=true;
912  switch (side)
913  {
914  case kMZ:
915  *n=G4ThreeVector(0,0,-1);
916  break;
917  case kPZ:
918  *n=G4ThreeVector(0,0,1);
919  break;
920  case kMY:
921  ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
922  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
923  break;
924  case kPY:
925  ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
926  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
927  break;
928  case kMX:
929  calpha=1/std::sqrt(1+fTalpha*fTalpha);
930  salpha=-calpha*fTalpha; // NOTE: actually use MINUS std::sin(alpha)
931  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
932  cosntheta=-1/std::sqrt(1+tntheta*tntheta);
933  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
934  break;
935  case kPX:
936  calpha=1/std::sqrt(1+fTalpha*fTalpha);
937  salpha=-calpha*fTalpha; // NOTE: actually use MINUS std::sin(alpha)
938  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
939  cosntheta=1/std::sqrt(1+tntheta*tntheta);
940  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
941  break;
942  default:
943  DumpInfo();
944  G4Exception("G4Para::DistanceToOut(p,v,..)",
945  "GeomSolids1002",JustWarning,
946  "Undefined side for valid surface normal to solid.");
947  break;
948  }
949  }
950  return snxt;
951 }
952 
954 //
955 // Calculate exact shortest distance to any boundary from inside
956 // - Returns 0 is point outside
957 
959 {
960  G4double safe=0.0;
961  G4double distz1,distz2,disty1,disty2,distx1,distx2;
962  G4double trany,cosy,tranx,cosx;
963 
964 #ifdef G4CSGDEBUG
965  if( Inside(p) == kOutside )
966  {
967  G4int oldprc = G4cout.precision(16) ;
968  G4cout << G4endl ;
969  DumpInfo();
970  G4cout << "Position:" << G4endl << G4endl ;
971  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
972  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
973  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
974  G4cout.precision(oldprc) ;
975  G4Exception("G4Para::DistanceToOut(p)", "GeomSolids1002",
976  JustWarning, "Point p is outside !?" );
977  }
978 #endif
979 
980  // Z planes
981  //
982  distz1=fDz-p.z();
983  distz2=fDz+p.z();
984  if (distz1<distz2)
985  {
986  safe=distz1;
987  }
988  else
989  {
990  safe=distz2;
991  }
992 
993  trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
994 
995  // Transformed x into `box' system
996  //
997  cosy=1.0/std::sqrt(1.0+fTthetaSphi*fTthetaSphi);
998  disty1=(fDy-trany)*cosy;
999  disty2=(fDy+trany)*cosy;
1000 
1001  if (disty1<safe) safe=disty1;
1002  if (disty2<safe) safe=disty2;
1003 
1004  tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
1005  cosx=1.0/std::sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
1006  distx1=(fDx-tranx)*cosx;
1007  distx2=(fDx+tranx)*cosx;
1008 
1009  if (distx1<safe) safe=distx1;
1010  if (distx2<safe) safe=distx2;
1011 
1012  if (safe<0) safe=0;
1013  return safe;
1014 }
1015 
1017 //
1018 // GetEntityType
1019 
1021 {
1022  return G4String("G4Para");
1023 }
1024 
1026 //
1027 // Make a clone of the object
1028 //
1030 {
1031  return new G4Para(*this);
1032 }
1033 
1035 //
1036 // Stream object contents to an output stream
1037 
1038 std::ostream& G4Para::StreamInfo( std::ostream& os ) const
1039 {
1040  G4int oldprc = os.precision(16);
1041  os << "-----------------------------------------------------------\n"
1042  << " *** Dump for solid - " << GetName() << " ***\n"
1043  << " ===================================================\n"
1044  << " Solid type: G4Para\n"
1045  << " Parameters: \n"
1046  << " half length X: " << fDx/mm << " mm \n"
1047  << " half length Y: " << fDy/mm << " mm \n"
1048  << " half length Z: " << fDz/mm << " mm \n"
1049  << " std::tan(alpha) : " << fTalpha/degree << " degrees \n"
1050  << " std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree
1051  << " degrees \n"
1052  << " std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree
1053  << " degrees \n"
1054  << "-----------------------------------------------------------\n";
1055  os.precision(oldprc);
1056 
1057  return os;
1058 }
1059 
1061 //
1062 // GetPointOnPlane
1063 // Auxiliary method for Get Point on Surface
1064 //
1065 
1066 G4ThreeVector G4Para::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1,
1067  G4ThreeVector p2, G4ThreeVector p3,
1068  G4double& area) const
1069 {
1070  G4double lambda1, lambda2, chose, aOne, aTwo;
1071  G4ThreeVector t, u, v, w, Area, normal;
1072 
1073  t = p1 - p0;
1074  u = p2 - p1;
1075  v = p3 - p2;
1076  w = p0 - p3;
1077 
1078  Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
1079  w.z()*v.x() - w.x()*v.z(),
1080  w.x()*v.y() - w.y()*v.x());
1081 
1082  aOne = 0.5*Area.mag();
1083 
1084  Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
1085  t.z()*u.x() - t.x()*u.z(),
1086  t.x()*u.y() - t.y()*u.x());
1087 
1088  aTwo = 0.5*Area.mag();
1089 
1090  area = aOne + aTwo;
1091 
1092  chose = G4RandFlat::shoot(0.,aOne+aTwo);
1093 
1094  if( (chose>=0.) && (chose < aOne) )
1095  {
1096  lambda1 = G4RandFlat::shoot(0.,1.);
1097  lambda2 = G4RandFlat::shoot(0.,lambda1);
1098  return (p2+lambda1*v+lambda2*w);
1099  }
1100 
1101  // else
1102 
1103  lambda1 = G4RandFlat::shoot(0.,1.);
1104  lambda2 = G4RandFlat::shoot(0.,lambda1);
1105  return (p0+lambda1*t+lambda2*u);
1106 }
1107 
1109 //
1110 // GetPointOnSurface
1111 //
1112 // Return a point (G4ThreeVector) randomly and uniformly
1113 // selected on the solid surface
1114 
1116 {
1117  G4ThreeVector One, Two, Three, Four, Five, Six;
1118  G4ThreeVector pt[8] ;
1119  G4double chose, aOne, aTwo, aThree, aFour, aFive, aSix;
1120 
1121  pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha-fDx,
1122  -fDz*fTthetaSphi-fDy, -fDz);
1123  pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha+fDx,
1124  -fDz*fTthetaSphi-fDy, -fDz);
1125  pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha-fDx,
1126  -fDz*fTthetaSphi+fDy, -fDz);
1127  pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha+fDx,
1128  -fDz*fTthetaSphi+fDy, -fDz);
1129  pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha-fDx,
1130  +fDz*fTthetaSphi-fDy, +fDz);
1131  pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha+fDx,
1132  +fDz*fTthetaSphi-fDy, +fDz);
1133  pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha-fDx,
1134  +fDz*fTthetaSphi+fDy, +fDz);
1135  pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha+fDx,
1136  +fDz*fTthetaSphi+fDy, +fDz);
1137 
1138  // make sure we provide the points in a clockwise fashion
1139 
1140  One = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
1141  Two = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
1142  Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
1143  Four = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour);
1144  Five = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
1145  Six = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
1146 
1147  chose = G4RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
1148 
1149  if( (chose>=0.) && (chose<aOne) )
1150  { return One; }
1151  else if(chose>=aOne && chose<aOne+aTwo)
1152  { return Two; }
1153  else if(chose>=aOne+aTwo && chose<aOne+aTwo+aThree)
1154  { return Three; }
1155  else if(chose>=aOne+aTwo+aThree && chose<aOne+aTwo+aThree+aFour)
1156  { return Four; }
1157  else if(chose>=aOne+aTwo+aThree+aFour && chose<aOne+aTwo+aThree+aFour+aFive)
1158  { return Five; }
1159  return Six;
1160 }
1161 
1163 //
1164 // Methods for visualisation
1165 
1167 {
1168  scene.AddSolid (*this);
1169 }
1170 
1172 {
1173  G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1174  G4double alpha = std::atan(fTalpha);
1175  G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
1176  +fTthetaSphi*fTthetaSphi));
1177 
1178  return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
1179 }
void set(double x, double y, double z)
G4String GetName() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Para.cc:740
Definition: G4Para.cc:64
ThreeVector shoot(const G4int Ap, const G4int Af)
G4ThreeVector GetPointOnSurface() const
Definition: G4Para.cc:1115
Definition: G4Cons.cc:76
Definition: G4Para.hh:77
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Para.cc:210
EInside Inside(const G4ThreeVector &p) const
Definition: G4Para.cc:317
static constexpr double mm
Definition: G4SIunits.hh:115
void SetAllParameters(G4double pDx, G4double pDy, G4double pDz, G4double pAlpha, G4double pTheta, G4double pPhi)
Definition: G4Para.cc:74
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
Definition: G4Para.cc:263
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
double x() const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const char * p
Definition: xmltok.h:285
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Para.cc:356
Definition: G4Para.cc:64
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Para.cc:517
G4Para(const G4String &pName, G4double pDx, G4double pDy, G4double pDz, G4double pAlpha, G4double pTheta, G4double pPhi)
Definition: G4Para.cc:102
G4double fCubicVolume
Definition: G4CSGSolid.hh:78
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Para.cc:1038
tuple x
Definition: test.py:50
G4Para & operator=(const G4Para &rhs)
Definition: G4Para.cc:186
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
const G4int smax
double z() const
G4Polyhedron * CreatePolyhedron() const
Definition: G4Para.cc:1171
void DumpInfo() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
Definition: G4Para.cc:68
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79
G4GLOB_DLL std::ostream G4cout
Definition: G4Para.cc:68
static constexpr double degree
Definition: G4SIunits.hh:144
Definition: G4Para.cc:64
bool G4bool
Definition: G4Types.hh:79
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
Definition: G4Para.cc:64
const G4int n
G4GeometryType GetEntityType() const
Definition: G4Para.cc:1020
std::vector< G4ThreeVector > G4ThreeVectorList
tuple v
Definition: test.py:18
G4double GetXHalfLength() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
ENSide
Definition: G4Para.cc:68
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Para.cc:1166
G4double GetTanAlpha() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
Hep3Vector unit() const
double y() const
Definition: G4Cons.cc:76
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetZHalfLength() const
tuple z
Definition: test.py:28
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4VSolid * Clone() const
Definition: G4Para.cc:1029
ESide
Definition: G4Cons.cc:76
virtual ~G4Para()
Definition: G4Para.cc:167
double G4double
Definition: G4Types.hh:76
static const G4double alpha
Definition: G4Para.cc:64
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:91
double mag() const
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4Para.cc:221
G4double GetYHalfLength() const
Definition: G4Cons.cc:80