Geant4  10.01.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::CheckPhiAngles(double sPhi, double dPhi)
110 {
111  CheckDPhiAngle(dPhi);
112  if ((fDPhi < 2 * UUtils::kPi) && (sPhi))
113  {
114  CheckSPhiAngle(sPhi);
115  }
116  InitializeTrigonometry();
117 }
118 
119 inline
120 void UCons::SetInnerRadiusMinusZ(double Rmin1)
121 {
122  fRmin1 = Rmin1 ;
123  Initialize();
124 }
125 
126 inline
127 void UCons::SetOuterRadiusMinusZ(double Rmax1)
128 {
129  fRmax1 = Rmax1 ;
130  Initialize();
131 }
132 
133 inline
134 void UCons::SetInnerRadiusPlusZ(double Rmin2)
135 {
136  fRmin2 = Rmin2 ;
137  Initialize();
138 }
139 
140 inline
141 void UCons::SetOuterRadiusPlusZ(double Rmax2)
142 {
143  fRmax2 = Rmax2 ;
144  Initialize();
145 }
146 
147 inline
148 void UCons::SetZHalfLength(double newDz)
149 {
150  fDz = newDz ;
151  Initialize();
152 }
153 
154 inline
155 void UCons::SetStartPhiAngle(double newSPhi, bool compute)
156 {
157  // Flag 'compute' can be used to explicitely avoid recomputation of
158  // trigonometry in case SetDeltaPhiAngle() is invoked afterwards
159 
160  CheckSPhiAngle(newSPhi);
161  fPhiFullCone = false;
162  if (compute)
163  {
164  InitializeTrigonometry();
165  }
166  Initialize();
167 }
168 
169 void UCons::SetDeltaPhiAngle(double newDPhi)
170 {
171  CheckPhiAngles(fSPhi, newDPhi);
172  Initialize();
173 }
174 
175 // Old access methods ...
176 
177 inline
178 double UCons::GetRmin1() const
179 {
180  return GetInnerRadiusMinusZ();
181 }
182 
183 inline
184 double UCons::GetRmax1() const
185 {
186  return GetOuterRadiusMinusZ();
187 }
188 
189 inline
190 double UCons::GetRmin2() const
191 {
192  return GetInnerRadiusPlusZ();
193 }
194 
195 inline
196 double UCons::GetRmax2() const
197 {
198  return GetOuterRadiusPlusZ();
199 }
200 
201 inline
202 double UCons::GetDz() const
203 {
204  return GetZHalfLength();
205 }
206 
207 inline
208 double UCons::GetSPhi() const
209 {
210  return GetStartPhiAngle();
211 }
212 
213 inline
214 double UCons::GetDPhi() const
215 {
216  return GetDeltaPhiAngle();
217 }
218 
219 inline
220 double UCons::Capacity()
221 {
222  if (fCubicVolume != 0.)
223  {
224  ;
225  }
226  else
227  {
228  double Rmean, rMean, deltaR, deltar;
229 
230  Rmean = 0.5 * (fRmax1 + fRmax2);
231  deltaR = fRmax1 - fRmax2;
232 
233  rMean = 0.5 * (fRmin1 + fRmin2);
234  deltar = fRmin1 - fRmin2;
235  fCubicVolume = fDPhi * fDz * (Rmean * Rmean - rMean * rMean
236  + (deltaR * deltaR - deltar * deltar) / 12);
237  }
238  return fCubicVolume;
239 }
240 
241 inline
242 double UCons::SurfaceArea()
243 {
244  if (fSurfaceArea != 0.)
245  {
246  ;
247  }
248  else
249  {
250  double mmin, mmax, dmin, dmax;
251 
252  mmin = (fRmin1 + fRmin2) * 0.5;
253  mmax = (fRmax1 + fRmax2) * 0.5;
254  dmin = (fRmin2 - fRmin1);
255  dmax = (fRmax2 - fRmax1);
256 
257  fSurfaceArea = fDPhi * (mmin * std::sqrt(dmin * dmin + 4 * fDz * fDz)
258  + mmax * std::sqrt(dmax * dmax + 4 * fDz * fDz)
259  + 0.5 * (fRmax1 * fRmax1 - fRmin1 * fRmin1
260  + fRmax2 * fRmax2 - fRmin2 * fRmin2));
261  if (!fPhiFullCone)
262  {
263  fSurfaceArea = fSurfaceArea + 4 * fDz * (mmax - mmin);
264  }
265  }
266  return fSurfaceArea;
267 }
268 
269 inline
270 double UCons::SafetyToPhi(const UVector3& p,
271  const double rho, bool& outside) const
272 {
273  double cosPsi, safePhi = 0.0;
274  outside = false;
275 
276  cosPsi = (p.x() * cosCPhi + p.y() * sinCPhi) / rho;
277 
278  if (cosPsi < std::cos(fDPhi * 0.5)) // Point lies outside phi range
279  {
280  outside = true;
281  if ((p.y() * cosCPhi - p.x() * sinCPhi) <= 0.0)
282  {
283  safePhi = std::fabs(p.x() * std::sin(fSPhi) - p.y() * std::cos(fSPhi));
284  }
285  else
286  {
287  safePhi = std::fabs(p.x() * sinEPhi - p.y() * cosEPhi);
288  }
289  }
290  return safePhi;
291 }
292 
293 inline
294 double UCons::SafetyFromInsideR(const UVector3& p,
295  const double rho, bool) const
296 {
297  double safe = 0.0, safeR1, safeR2, safePhi;
298  double pRMin;
299  double pRMax;
300 
301  if (fRmin1 || fRmin2)
302  {
303  pRMin = tanRMin * p.z() + (fRmin1 + fRmin2) * 0.5;
304  safeR1 = (rho - pRMin) / secRMin;
305  }
306  else
307  {
308  safeR1 = UUtils::kInfinity;
309  }
310 
311  pRMax = tanRMax * p.z() + (fRmax1 + fRmax2) * 0.5;
312  safeR2 = (pRMax - rho) / secRMax;
313 
314  if (safeR1 < safeR2)
315  {
316  safe = safeR1;
317  }
318  else
319  {
320  safe = safeR2;
321  }
322 
323  // Check if phi divided, Calc distances closest phi plane
324  //
325  if (!fPhiFullCone)
326  {
327  // Above/below central phi of UCons?
328 
329  if ((p.y() * cosCPhi - p.x() * sinCPhi) <= 0)
330  {
331  safePhi = -(p.x() * sinSPhi - p.y() * cosSPhi);
332  }
333  else
334  {
335  safePhi = (p.x() * sinEPhi - p.y() * cosEPhi);
336  }
337  if (safePhi < safe)
338  {
339  safe = safePhi;
340  }
341  }
342 
343  if (safe < 0)
344  {
345  safe = 0;
346  }
347 
348  return safe;
349 }
350 
351 inline
352 double UCons::SafetyFromOutsideR(const UVector3& p,
353  const double rho, bool) const
354 {
355  double safe = 0.0, safeR1, safeR2;
356  double safePhi;
357  double pRMin, pRMax;
358  bool outside;
359  if (fRmin1 || fRmin2)
360  {
361  pRMin = tanRMin * p.z() + (fRmin1 + fRmin2) * 0.5;
362  safeR1 = (rho-pRMin ) / secRMin;
363 
364  pRMax = tanRMax * p.z() + (fRmax1 + fRmax2) * 0.5;
365  safeR2 = (rho - pRMax) / secRMax;
366 
367  if (safeR1 > safeR2)
368  {
369  safe = safeR1;
370  }
371  else
372  {
373  safe = safeR2;
374  }
375  }
376  else
377  {
378  pRMax = tanRMax * p.z() + (fRmax1 + fRmax2) * 0.5;
379  safe = (rho - pRMax) / secRMax;
380  }
381  if (!fPhiFullCone)
382  {
383  safePhi=SafetyToPhi(p,rho,outside);
384  if ((outside) && (safePhi > safe))
385  {
386  safe = safePhi;
387  }
388  }
389 
390  if (safe < 0.0)
391  {
392  safe = 0.0;
393  }
394  return safe; // not accurate safety
395 }
396 
397 inline
398 VUSolid::EnumInside UCons::Inside(const UVector3& p) const
399 {
400  double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2;
401  VUSolid::EnumInside in;
402  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
403  static const double halfRadTolerance = kRadTolerance * 0.5;
404  static const double halfAngTolerance = kAngTolerance * 0.5;
405 
406  if (std::fabs(p.z()) > fDz + halfCarTolerance)
407  {
408  return in = eOutside;
409  }
410  else if (std::fabs(p.z()) >= fDz - halfCarTolerance)
411  {
412  in = eSurface;
413  }
414  else
415  {
416  in = eInside;
417  }
418  r2 = p.x() * p.x() + p.y() * p.y();
419  rl = 0.5 * (fRmin2 * (p.z() + fDz) + fRmin1 * (fDz - p.z())) / fDz;
420  rh = 0.5 * (fRmax2 * (p.z() + fDz) + fRmax1 * (fDz - p.z())) / fDz;
421 
422  tolRMin = rl - halfRadTolerance;
423  if (tolRMin < 0)
424  {
425  tolRMin = 0;
426  }
427  tolRMax = rh + halfRadTolerance;
428 
429  if ((r2 < tolRMin * tolRMin) || (r2 > tolRMax * tolRMax))
430  {
431  return in = eOutside;
432  }
433  if (rl)
434  {
435  tolRMin = rl + halfRadTolerance;
436  }
437  else
438  {
439  tolRMin = 0.0;
440  }
441 
442  tolRMax = rh - halfRadTolerance;
443  if (in == eInside) // else it's eSurface already
444  {
445  if ((r2 < tolRMin * tolRMin) || (r2 >= tolRMax * tolRMax))
446  {
447  in = eSurface;
448  }
449  }
450  if (!fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)))
451  {
452  pPhi = std::atan2(p.y(), p.x());
453  if (pPhi < fSPhi - halfAngTolerance)
454  {
455  pPhi += 2 * UUtils::kPi;
456  }
457  else if (pPhi > fSPhi + fDPhi + halfAngTolerance)
458  {
459  pPhi -= 2 * UUtils::kPi;
460  }
461 
462  if ((pPhi < fSPhi - halfAngTolerance) ||
463  (pPhi > fSPhi + fDPhi + halfAngTolerance))
464  {
465  return in = eOutside;
466  }
467  else if (in == eInside) // else it's eSurface anyway already
468  {
469  if ((pPhi < fSPhi + halfAngTolerance) ||
470  (pPhi > fSPhi + fDPhi - halfAngTolerance))
471  {
472  in = eSurface;
473  }
474  }
475  }
476  else if (!fPhiFullCone)
477  {
478  in = eSurface;
479  }
480 
481  return in;
482 }