Geant4  10.00.p03
VUSolid.cc
Go to the documentation of this file.
1 
2 #include "VUSolid.hh"
3 
5 // "Universal" Solid Interface
6 // Authors: J. Apostolakis, G. Cosmo, M. Gayer, A. Gheata, A. Munnich, T. Nikitina (CERN)
7 //
8 // Created: 25 May 2011
9 //
11 
12 double VUSolid::fgTolerance = 1.0E-9; // cartesian tolerance; to be changed (for U was 1e-8, but we keep Geant4)
13 double VUSolid::frTolerance = 1.0E-9; // radial tolerance; to be changed
14 
15 double VUSolid::faTolerance = 1.0E-9; // angular tolerance; to be changed
16 
17 //______________________________________________________________________________
19 {
20 
21 }
22 
23 //______________________________________________________________________________
24 VUSolid::VUSolid(const std::string &name) :
25 fName(name)
26 {
27  // Named constructor
28  SetName(name);
29 }
30 
31 //______________________________________________________________________________
33 {
34 
35 }
36 
37 /*
38 int UIntersectingCone::LineHitsCone2( const UVector3 &p,
39  const UVector3 &v,
40  double &s1, double &s2 )
41 {
42  double x0 = p.x, y0 = p.y, z0 = p.z;
43  double tx = v.x, ty = v.y, tz = v.z;
44 
45  // Special case which might not be so rare: B = 0 (precisely)
46  //
47  if (B==0)
48  {
49  if (std::fabs(tz) < 1/UUtils::kInfinity) { return 0; }
50 
51  s1 = (A-z0)/tz;
52  return 1;
53  }
54 
55  double B2 = B*B;
56 
57  double a = tz*tz - B2*(tx*tx + ty*ty);
58  double b = 2*( (z0-A)*tz - B2*(x0*tx + y0*ty) );
59  double c = UUtils::sqr(z0-A) - B2*( x0*x0 + y0*y0 );
60 
61  double radical = b*b - 4*a*c;
62 
63  if (radical < -1E-6*std::fabs(b)) { return 0; } // No solution
64 
65  if (radical < 1E-6*std::fabs(b))
66  {
67  //
68  // The radical is roughly zero: check for special, very rare, cases
69  //
70  if (std::fabs(a) > 1/UUtils::kInfinity)
71  {
72  if ( std::fabs(x0*ty - y0*tx) < std::fabs(1E-6/B) )
73  {
74  s1 = -0.5*b/a;
75  return 1;
76  }
77  return 0;
78  }
79  }
80  else
81  {
82  radical = std::sqrt(radical);
83  }
84 
85  if (a < -1/UUtils::kInfinity)
86  {
87  double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
88  sa = q/a;
89  sb = c/q;
90  if (sa < sb) { s1 = sa; s2 = sb; } else { s1 = sb; s2 = sa; }
91  if ((z0 + (s1)*tz - A)/B < 0) { return 0; }
92  return 2;
93  }
94  else if (a > 1/UUtils::kInfinity)
95  {
96  double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
97  sa = q/a;
98  sb = c/q;
99  s1 = (tz*B > 0)^(sa > sb) ? sb : sa;
100  return 1;
101  }
102  else if (std::fabs(b) < 1/UUtils::kInfinity)
103  {
104  return 0;
105  }
106  else
107  {
108  s1 = -c/b;
109  if ((z0 + (s1)*tz - A)/B < 0) { return 0; }
110  return 1;
111  }
112 }
113 int UIntersectingCone::LineHitsCone2( const UVector3 &p,
114  const UVector3 &v,
115  double &s1, double &s2 )
116 {
117  double x0 = p.x, y0 = p.y, z0 = p.z;
118  double tx = v.x, ty = v.y, tz = v.z;
119 
120  // Special case which might not be so rare: B = 0 (precisely)
121  //
122  if (B==0)
123  {
124  if (std::fabs(tz) < 1/UUtils::kInfinity) { return 0; }
125 
126  s1 = (A-z0)/tz;
127  return 1;
128  }
129 
130  double B2 = B*B;
131 
132  double a = tz*tz - B2*(tx*tx + ty*ty);
133  double b = 2*( (z0-A)*tz - B2*(x0*tx + y0*ty) );
134  double c = UUtils::sqr(z0-A) - B2*( x0*x0 + y0*y0 );
135 
136  double radical = b*b - 4*a*c;
137 
138  if (radical < -1E-6*std::fabs(b)) { return 0; } // No solution
139 
140  if (radical < 1E-6*std::fabs(b))
141  {
142  //
143  // The radical is roughly zero: check for special, very rare, cases
144  //
145  if (std::fabs(a) > 1/UUtils::kInfinity)
146  {
147  if ( std::fabs(x0*ty - y0*tx) < std::fabs(1E-6/B) )
148  {
149  s1 = -0.5*b/a;
150  return 1;
151  }
152  return 0;
153  }
154  }
155  else
156  {
157  radical = std::sqrt(radical);
158  }
159 
160  if (a < -1/UUtils::kInfinity)
161  {
162  double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
163  sa = q/a;
164  sb = c/q;
165  if (sa < sb) { s1 = sa; s2 = sb; } else { s1 = sb; s2 = sa; }
166  if ((z0 + (s1)*tz - A)/B < 0) { return 0; }
167  return 2;
168  }
169  else if (a > 1/UUtils::kInfinity)
170  {
171  double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
172  sa = q/a;
173  sb = c/q;
174  s1 = (tz*B > 0)^(sa > sb) ? sb : sa;
175  return 1;
176  }
177  else if (std::fabs(b) < 1/UUtils::kInfinity)
178  {
179  return 0;
180  }
181  else
182  {
183  s1 = -c/b;
184  if ((z0 + (s1)*tz - A)/B < 0) { return 0; }
185  return 1;
186  }
187 }
188 */
189 
191 //
192 // Returns an estimation of the solid volume in internal units.
193 // The number of statistics and error accuracy is fixed.
194 // This method may be overloaded by derived classes to compute the
195 // exact geometrical quantity for solids where this is possible.
196 // or anyway to cache the computed value.
197 // This implementation does NOT cache the computed value.
198 
200 {
201  int cubVolStatistics = 1000000;
202  double cubVolEpsilon = 0.001;
203  return EstimateCubicVolume(cubVolStatistics, cubVolEpsilon);
204 }
205 
207 //
208 // Calculate cubic volume based on Inside() method.
209 // Accuracy is limited by the second argument or the statistics
210 // expressed by the first argument.
211 // Implementation is courtesy of Vasiliki Despoina Mitsou,
212 // University of Athens.
213 
214 double VUSolid::EstimateCubicVolume(int nStat, double epsilon) const
215 {
216  int iInside=0;
217  double px,py,pz,volume;
218  UVector3 min,max;
219  UVector3 p;
221 
222  // values needed for CalculateExtent signature
223 
224  // min max extents of pSolid along X,Y,Z
225 
226  this->Extent(min,max);
227 
228  // limits
229 
230  if(nStat < 100) nStat = 100;
231  if(epsilon > 0.01) epsilon = 0.01;
232 
233  for(int i = 0; i < nStat; i++ )
234  {
235  px = min.x+(max.x-min.x)*UUtils::Random();
236  py = min.y+(max.y-min.y)*UUtils::Random();
237  pz = min.z+(max.z-min.z)*UUtils::Random();
238  p = UVector3(px,py,pz);
239  in = this->Inside(p);
240  if(in != eOutside) iInside++;
241  }
242  volume = (max.x-min.x)*(max.y-min.y)*(max.z-min.z)*iInside/nStat;
243  return volume;
244 }
245 
247 //
248 // Returns an estimation of the solid surface area in internal units.
249 // The number of statistics and error accuracy is fixed.
250 // This method may be overloaded by derived classes to compute the
251 // exact geometrical quantity for solids where this is possible.
252 // or anyway to cache the computed value.
253 // This implementation does NOT cache the computed value.
254 
256 {
257  int stat = 1000000;
258  double ell = -1.;
259  return EstimateSurfaceArea(stat,ell);
260 }
261 
263 //
264 // Estimate surface area based on Inside(), DistanceToIn(), and
265 // DistanceToOut() methods. Accuracy is limited by the statistics
266 // defined by the first argument. Implemented by Mikhail Kosov.
267 
268 double VUSolid::EstimateSurfaceArea(int nStat, double ell) const
269 {
270  int inside=0;
271  double px,py,pz,surf;
272  UVector3 min,max;
273  UVector3 p;
275 
276  // values needed for CalculateExtent signature
277 
278  // min max extents of pSolid along X,Y,Z
279 
280  this->Extent(min,max);
281 
282  // limits
283 
284  if(nStat < 100) { nStat = 100; }
285 
286  double dX=max.x-min.x;
287  double dY=max.y-min.y;
288  double dZ=max.z-min.z;
289  if(ell<=0.) // Automatic definition of skin thickness
290  {
291  double minval=dX;
292  if(dY<dX) { minval=dY; }
293  if(dZ<minval) { minval=dZ; }
294  ell=.01*minval;
295  }
296 
297  double dd=2*ell;
298  min.x-=ell; min.y-=ell; min.z-=ell; dX+=dd; dY+=dd; dZ+=dd;
299 
300  for(int i = 0; i < nStat; i++ )
301  {
302  px = min.x+dX*UUtils::Random();
303  py = min.y+dY*UUtils::Random();
304  pz = min.z+dZ*UUtils::Random();
305  p = UVector3(px,py,pz);
306  in = this->Inside(p);
307  if(in != eOutside)
308  {
309  if (SafetyFromInside(p)<ell) { inside++; }
310  }
311  else if(SafetyFromOutside(p)<ell) { inside++; }
312  }
313  // @@ The conformal correction can be upgraded
314  surf = dX*dY*dZ*inside/dd/nStat;
315  return surf;
316 }
317 
318 void VUSolid::ExtentAxis(EAxisType aAxis, double &aMin, double &aMax) const
319 // Returns the minimum and maximum extent along the specified Cartesian axis
320 {
321  // Returns extent of the solid along a given cartesian axis
322  if (aAxis >= 0 && aAxis <= 2)
323  {
324  UVector3 min, max;
325  Extent(min,max);
326  aMin = min[aAxis]; aMax = max[aAxis];
327  }
328 #ifdef USPECSDEBUG
329  else
330  cout << "Extent: unknown axis" << aAxis << std::endl;
331 #endif
332 }
333 
335 {
337 }
339 {
341 }
343 {
345 }
346 
347 
VUSolid()
Definition: VUSolid.cc:18
void SetRadTolerance(double eps)
Definition: VUSolid.cc:338
static double frTolerance
Definition: VUSolid.hh:31
virtual double SurfaceArea()=0
Definition: VUSolid.cc:255
G4String fName
Definition: G4AttUtils.hh:55
void SetCarTolerance(double eps)
Definition: VUSolid.cc:334
virtual ~VUSolid()
Definition: VUSolid.cc:32
static double fgTolerance
Definition: VUSolid.hh:30
G4String name
Definition: TRTMaterials.hh:40
void SetAngTolerance(double eps)
Definition: VUSolid.cc:342
double EstimateCubicVolume(int nStat, double epsilon) const
Definition: VUSolid.cc:214
static const G4double eps
virtual EnumInside Inside(const UVector3 &aPoint) const =0
double EstimateSurfaceArea(int nStat, double ell) const
Definition: VUSolid.cc:268
double x
Definition: UVector3.hh:136
void SetName(const std::string &aName)
Definition: VUSolid.hh:104
static double faTolerance
Definition: VUSolid.hh:32
EAxisType
Definition: VUSolid.hh:27
EnumInside
Definition: VUSolid.hh:23
virtual double SafetyFromInside(const UVector3 &aPoint, bool aAccurate=false) const =0
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual void ExtentAxis(EAxisType aAxis, double &aMin, double &aMax) const
Definition: VUSolid.cc:318
virtual double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const =0
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
virtual void Extent(UVector3 &aMin, UVector3 &aMax) const =0
double z
Definition: UVector3.hh:138
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
double y
Definition: UVector3.hh:137
virtual double Capacity()=0
Definition: VUSolid.cc:199