Geant4  10.00.p03
UCons.icc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UCons.icc
12 //
13 // Implementation of inline methods of UCons
14 //
15 // 19.10.12 Marek Gayer
16 // Created from original implementation in Geant4
17 // --------------------------------------------------------------------
18 
19 inline
20 double UCons::GetInnerRadiusMinusZ() const
21 {
22  return fRmin1 ;
23 }
24 
25 inline
26 double UCons::GetOuterRadiusMinusZ() const
27 {
28  return fRmax1 ;
29 }
30 
31 inline
32 double UCons::GetInnerRadiusPlusZ() const
33 {
34  return fRmin2 ;
35 }
36 
37 inline
38 double UCons::GetOuterRadiusPlusZ() const
39 {
40  return fRmax2 ;
41 }
42 
43 inline
44 double UCons::GetZHalfLength() const
45 {
46  return fDz ;
47 }
48 
49 inline
50 double UCons::GetStartPhiAngle() const
51 {
52  return fSPhi ;
53 }
54 
55 inline
56 double UCons::GetDeltaPhiAngle() const
57 {
58  return fDPhi;
59 }
60 
61 inline
62 void UCons::Initialize()
63 {
64  fCubicVolume = 0.;
65  fSurfaceArea = 0.;
66 
67  tanRMin = (fRmin2 - fRmin1) * 0.5 / fDz;
68  secRMin = std::sqrt(1.0 + tanRMin * tanRMin);
69 
70  tanRMax = (fRmax2 - fRmax1) * 0.5 / fDz;
71  secRMax = std::sqrt(1.0 + tanRMax * tanRMax);
72 }
73 
74 inline
75 void UCons::InitializeTrigonometry()
76 {
77  double hDPhi = 0.5 * fDPhi; // half delta phi
78  double cPhi = fSPhi + hDPhi;
79  double ePhi = fSPhi + fDPhi;
80 
81  sinCPhi = std::sin(cPhi);
82  cosCPhi = std::cos(cPhi);
83  cosHDPhiIT = std::cos(hDPhi - 0.5 * kAngTolerance); // inner/outer tol half dphi
84  cosHDPhiOT = std::cos(hDPhi + 0.5 * kAngTolerance);
85  sinSPhi = std::sin(fSPhi);
86  cosSPhi = std::cos(fSPhi);
87  sinEPhi = std::sin(ePhi);
88  cosEPhi = std::cos(ePhi);
89 }
90 
91 inline void UCons::CheckSPhiAngle(double sPhi)
92 {
93  // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
94 
95  if (sPhi < 0)
96  {
97  fSPhi = 2 * UUtils::kPi - std::fmod(std::fabs(sPhi), 2 * UUtils::kPi);
98  }
99  else
100  {
101  fSPhi = std::fmod(sPhi, 2 * UUtils::kPi) ;
102  }
103  if (fSPhi + fDPhi > 2 * UUtils::kPi)
104  {
105  fSPhi -= 2 * UUtils::kPi ;
106  }
107 }
108 
109 inline void UCons::CheckDPhiAngle(double dPhi)
110 {
111  fPhiFullCone = true;
112  if (dPhi >= 2 * UUtils::kPi - kAngTolerance * 0.5)
113  {
114  fDPhi = 2 * UUtils::kPi;
115  fSPhi = 0;
116  }
117  else
118  {
119  fPhiFullCone = false;
120  if (dPhi > 0)
121  {
122  fDPhi = dPhi;
123  }
124  else
125  {
126  std::ostringstream message;
127  message << "Invalid dphi." << std::endl
128  << "Negative or zero delta-Phi (" << dPhi << ") in solid: "
129  << GetName();
130  UUtils::Exception("UCons::CheckDPhiAngle()", "GeomSolids0002",
131  FatalErrorInArguments, 1, message.str().c_str());
132  }
133  }
134 }
135 
136 inline void UCons::CheckPhiAngles(double sPhi, double dPhi)
137 {
138  CheckDPhiAngle(dPhi);
139  if ((fDPhi < 2 * UUtils::kPi) && (sPhi))
140  {
141  CheckSPhiAngle(sPhi);
142  }
143  InitializeTrigonometry();
144 }
145 
146 inline
147 void UCons::SetInnerRadiusMinusZ(double Rmin1)
148 {
149  fRmin1 = Rmin1 ;
150  Initialize();
151 }
152 
153 inline
154 void UCons::SetOuterRadiusMinusZ(double Rmax1)
155 {
156  fRmax1 = Rmax1 ;
157  Initialize();
158 }
159 
160 inline
161 void UCons::SetInnerRadiusPlusZ(double Rmin2)
162 {
163  fRmin2 = Rmin2 ;
164  Initialize();
165 }
166 
167 inline
168 void UCons::SetOuterRadiusPlusZ(double Rmax2)
169 {
170  fRmax2 = Rmax2 ;
171  Initialize();
172 }
173 
174 inline
175 void UCons::SetZHalfLength(double newDz)
176 {
177  fDz = newDz ;
178  Initialize();
179 }
180 
181 inline
182 void UCons::SetStartPhiAngle(double newSPhi, bool compute)
183 {
184  // Flag 'compute' can be used to explicitely avoid recomputation of
185  // trigonometry in case SetDeltaPhiAngle() is invoked afterwards
186 
187  CheckSPhiAngle(newSPhi);
188  fPhiFullCone = false;
189  if (compute)
190  {
191  InitializeTrigonometry();
192  }
193  Initialize();
194 }
195 
196 void UCons::SetDeltaPhiAngle(double newDPhi)
197 {
198  CheckPhiAngles(fSPhi, newDPhi);
199  Initialize();
200 }
201 
202 // Old access methods ...
203 
204 inline
205 double UCons::GetRmin1() const
206 {
207  return GetInnerRadiusMinusZ();
208 }
209 
210 inline
211 double UCons::GetRmax1() const
212 {
213  return GetOuterRadiusMinusZ();
214 }
215 
216 inline
217 double UCons::GetRmin2() const
218 {
219  return GetInnerRadiusPlusZ();
220 }
221 
222 inline
223 double UCons::GetRmax2() const
224 {
225  return GetOuterRadiusPlusZ();
226 }
227 
228 inline
229 double UCons::GetDz() const
230 {
231  return GetZHalfLength();
232 }
233 
234 inline
235 double UCons::GetSPhi() const
236 {
237  return GetStartPhiAngle();
238 }
239 
240 inline
241 double UCons::GetDPhi() const
242 {
243  return GetDeltaPhiAngle();
244 }
245 
246 inline
247 double UCons::Capacity()
248 {
249  if (fCubicVolume != 0.)
250  {
251  ;
252  }
253  else
254  {
255  double Rmean, rMean, deltaR, deltar;
256 
257  Rmean = 0.5 * (fRmax1 + fRmax2);
258  deltaR = fRmax1 - fRmax2;
259 
260  rMean = 0.5 * (fRmin1 + fRmin2);
261  deltar = fRmin1 - fRmin2;
262  fCubicVolume = fDPhi * fDz * (Rmean * Rmean - rMean * rMean
263  + (deltaR * deltaR - deltar * deltar) / 12);
264  }
265  return fCubicVolume;
266 }
267 
268 inline
269 double UCons::SurfaceArea()
270 {
271  if (fSurfaceArea != 0.)
272  {
273  ;
274  }
275  else
276  {
277  double mmin, mmax, dmin, dmax;
278 
279  mmin = (fRmin1 + fRmin2) * 0.5;
280  mmax = (fRmax1 + fRmax2) * 0.5;
281  dmin = (fRmin2 - fRmin1);
282  dmax = (fRmax2 - fRmax1);
283 
284  fSurfaceArea = fDPhi * (mmin * std::sqrt(dmin * dmin + 4 * fDz * fDz)
285  + mmax * std::sqrt(dmax * dmax + 4 * fDz * fDz)
286  + 0.5 * (fRmax1 * fRmax1 - fRmin1 * fRmin1
287  + fRmax2 * fRmax2 - fRmin2 * fRmin2));
288  if (!fPhiFullCone)
289  {
290  fSurfaceArea = fSurfaceArea + 4 * fDz * (mmax - mmin);
291  }
292  }
293  return fSurfaceArea;
294 }
295 
296 inline
297 double UCons::SafetyToPhi(const UVector3& p,
298  const double rho, bool& outside) const
299 {
300  double cosPsi, safePhi = 0.0;
301  outside = false;
302 
303  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / rho;
304 
305  if (cosPsi < std::cos(fDPhi * 0.5)) // Point lies outside phi range
306  {
307  outside = true;
308  if ((p.y * cosCPhi - p.x * sinCPhi) <= 0.0)
309  {
310  safePhi = std::fabs(p.x * std::sin(fSPhi) - p.y * std::cos(fSPhi));
311  }
312  else
313  {
314  safePhi = std::fabs(p.x * sinEPhi - p.y * cosEPhi);
315  }
316  }
317  return safePhi;
318 }
319 
320 inline
321 double UCons::SafetyFromInsideR(const UVector3& p,
322  const double rho, bool) const
323 {
324  double safe = 0.0, safeR1, safeR2, safePhi;
325  double pRMin;
326  double pRMax;
327 
328  if (fRmin1 || fRmin2)
329  {
330  pRMin = tanRMin * p.z + (fRmin1 + fRmin2) * 0.5;
331  safeR1 = (rho - pRMin) / secRMin;
332  }
333  else
334  {
335  safeR1 = UUtils::kInfinity;
336  }
337 
338  pRMax = tanRMax * p.z + (fRmax1 + fRmax2) * 0.5;
339  safeR2 = (pRMax - rho) / secRMax;
340 
341  if (safeR1 < safeR2)
342  {
343  safe = safeR1;
344  }
345  else
346  {
347  safe = safeR2;
348  }
349 
350  // Check if phi divided, Calc distances closest phi plane
351  //
352  if (!fPhiFullCone)
353  {
354  // Above/below central phi of UCons?
355 
356  if ((p.y * cosCPhi - p.x * sinCPhi) <= 0)
357  {
358  safePhi = -(p.x * sinSPhi - p.y * cosSPhi);
359  }
360  else
361  {
362  safePhi = (p.x * sinEPhi - p.y * cosEPhi);
363  }
364  if (safePhi < safe)
365  {
366  safe = safePhi;
367  }
368  }
369 
370  if (safe < 0)
371  {
372  safe = 0;
373  }
374 
375  return safe;
376 }
377 
378 inline
379 double UCons::SafetyFromOutsideR(const UVector3& p,
380  const double rho, bool) const
381 {
382  double safe = 0.0, safeR1, safeR2;
383  double safePhi;
384  double pRMin, pRMax;
385  bool outside;
386  if (fRmin1 || fRmin2)
387  {
388  pRMin = tanRMin * p.z + (fRmin1 + fRmin2) * 0.5;
389  safeR1 = (rho-pRMin ) / secRMin;
390 
391  pRMax = tanRMax * p.z + (fRmax1 + fRmax2) * 0.5;
392  safeR2 = (rho - pRMax) / secRMax;
393 
394  if (safeR1 > safeR2)
395  {
396  safe = safeR1;
397  }
398  else
399  {
400  safe = safeR2;
401  }
402  }
403  else
404  {
405  pRMax = tanRMax * p.z + (fRmax1 + fRmax2) * 0.5;
406  safe = (rho - pRMax) / secRMax;
407  }
408  if (!fPhiFullCone)
409  {
410  safePhi=SafetyToPhi(p,rho,outside);
411  if ((outside) && (safePhi > safe))
412  {
413  safe = safePhi;
414  }
415  }
416 
417  if (safe < 0.0)
418  {
419  safe = 0.0;
420  }
421  return safe; // not accurate safety
422 }
423 
424 inline
425 VUSolid::EnumInside UCons::Inside(const UVector3& p) const
426 {
427  double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2;
428  VUSolid::EnumInside in;
429  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
430  static const double halfRadTolerance = kRadTolerance * 0.5;
431  static const double halfAngTolerance = kAngTolerance * 0.5;
432 
433  if (std::fabs(p.z) > fDz + halfCarTolerance)
434  {
435  return in = eOutside;
436  }
437  else if (std::fabs(p.z) >= fDz - halfCarTolerance)
438  {
439  in = eSurface;
440  }
441  else
442  {
443  in = eInside;
444  }
445  r2 = p.x * p.x + p.y * p.y;
446  rl = 0.5 * (fRmin2 * (p.z + fDz) + fRmin1 * (fDz - p.z)) / fDz;
447  rh = 0.5 * (fRmax2 * (p.z + fDz) + fRmax1 * (fDz - p.z)) / fDz;
448 
449  tolRMin = rl - halfRadTolerance;
450  if (tolRMin < 0)
451  {
452  tolRMin = 0;
453  }
454  tolRMax = rh + halfRadTolerance;
455 
456  if ((r2 < tolRMin * tolRMin) || (r2 > tolRMax * tolRMax))
457  {
458  return in = eOutside;
459  }
460  if (rl)
461  {
462  tolRMin = rl + halfRadTolerance;
463  }
464  else
465  {
466  tolRMin = 0.0;
467  }
468 
469  tolRMax = rh - halfRadTolerance;
470  if (in == eInside) // else it's eSurface already
471  {
472  if ((r2 < tolRMin * tolRMin) || (r2 >= tolRMax * tolRMax))
473  {
474  in = eSurface;
475  }
476  }
477  if (!fPhiFullCone && ((p.x != 0.0) || (p.y != 0.0)))
478  {
479  pPhi = std::atan2(p.y, p.x);
480  if (pPhi < fSPhi - halfAngTolerance)
481  {
482  pPhi += 2 * UUtils::kPi;
483  }
484  else if (pPhi > fSPhi + fDPhi + halfAngTolerance)
485  {
486  pPhi -= 2 * UUtils::kPi;
487  }
488 
489  if ((pPhi < fSPhi - halfAngTolerance) ||
490  (pPhi > fSPhi + fDPhi + halfAngTolerance))
491  {
492  return in = eOutside;
493  }
494  else if (in == eInside) // else it's eSurface anyway already
495  {
496  if ((pPhi < fSPhi + halfAngTolerance) ||
497  (pPhi > fSPhi + fDPhi - halfAngTolerance))
498  {
499  in = eSurface;
500  }
501  }
502  }
503  else if (!fPhiFullCone)
504  {
505  in = eSurface;
506  }
507 
508  return in;
509 }