2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
9 // --------------------------------------------------------------------
13 // Implementation of inline methods of UCons
15 // 19.10.12 Marek Gayer
16 // Created from original implementation in Geant4
17 // --------------------------------------------------------------------
20 double UCons::GetInnerRadiusMinusZ() const
26 double UCons::GetOuterRadiusMinusZ() const
32 double UCons::GetInnerRadiusPlusZ() const
38 double UCons::GetOuterRadiusPlusZ() const
44 double UCons::GetZHalfLength() const
50 double UCons::GetStartPhiAngle() const
56 double UCons::GetDeltaPhiAngle() const
62 void UCons::Initialize()
67 tanRMin = (fRmin2 - fRmin1) * 0.5 / fDz;
68 secRMin = std::sqrt(1.0 + tanRMin * tanRMin);
70 tanRMax = (fRmax2 - fRmax1) * 0.5 / fDz;
71 secRMax = std::sqrt(1.0 + tanRMax * tanRMax);
75 void UCons::InitializeTrigonometry()
77 double hDPhi = 0.5 * fDPhi; // half delta phi
78 double cPhi = fSPhi + hDPhi;
79 double ePhi = fSPhi + fDPhi;
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);
91 inline void UCons::CheckSPhiAngle(double sPhi)
93 // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
97 fSPhi = 2 * UUtils::kPi - std::fmod(std::fabs(sPhi), 2 * UUtils::kPi);
101 fSPhi = std::fmod(sPhi, 2 * UUtils::kPi) ;
103 if (fSPhi + fDPhi > 2 * UUtils::kPi)
105 fSPhi -= 2 * UUtils::kPi ;
109 inline void UCons::CheckDPhiAngle(double dPhi)
112 if (dPhi >= 2 * UUtils::kPi - kAngTolerance * 0.5)
114 fDPhi = 2 * UUtils::kPi;
119 fPhiFullCone = false;
126 std::ostringstream message;
127 message << "Invalid dphi." << std::endl
128 << "Negative or zero delta-Phi (" << dPhi << ") in solid: "
130 UUtils::Exception("UCons::CheckDPhiAngle()", "GeomSolids0002",
131 FatalErrorInArguments, 1, message.str().c_str());
136 inline void UCons::CheckPhiAngles(double sPhi, double dPhi)
138 CheckDPhiAngle(dPhi);
139 if ((fDPhi < 2 * UUtils::kPi) && (sPhi))
141 CheckSPhiAngle(sPhi);
143 InitializeTrigonometry();
147 void UCons::SetInnerRadiusMinusZ(double Rmin1)
154 void UCons::SetOuterRadiusMinusZ(double Rmax1)
161 void UCons::SetInnerRadiusPlusZ(double Rmin2)
168 void UCons::SetOuterRadiusPlusZ(double Rmax2)
175 void UCons::SetZHalfLength(double newDz)
182 void UCons::SetStartPhiAngle(double newSPhi, bool compute)
184 // Flag 'compute' can be used to explicitely avoid recomputation of
185 // trigonometry in case SetDeltaPhiAngle() is invoked afterwards
187 CheckSPhiAngle(newSPhi);
188 fPhiFullCone = false;
191 InitializeTrigonometry();
196 void UCons::SetDeltaPhiAngle(double newDPhi)
198 CheckPhiAngles(fSPhi, newDPhi);
202 // Old access methods ...
205 double UCons::GetRmin1() const
207 return GetInnerRadiusMinusZ();
211 double UCons::GetRmax1() const
213 return GetOuterRadiusMinusZ();
217 double UCons::GetRmin2() const
219 return GetInnerRadiusPlusZ();
223 double UCons::GetRmax2() const
225 return GetOuterRadiusPlusZ();
229 double UCons::GetDz() const
231 return GetZHalfLength();
235 double UCons::GetSPhi() const
237 return GetStartPhiAngle();
241 double UCons::GetDPhi() const
243 return GetDeltaPhiAngle();
247 double UCons::Capacity()
249 if (fCubicVolume != 0.)
255 double Rmean, rMean, deltaR, deltar;
257 Rmean = 0.5 * (fRmax1 + fRmax2);
258 deltaR = fRmax1 - fRmax2;
260 rMean = 0.5 * (fRmin1 + fRmin2);
261 deltar = fRmin1 - fRmin2;
262 fCubicVolume = fDPhi * fDz * (Rmean * Rmean - rMean * rMean
263 + (deltaR * deltaR - deltar * deltar) / 12);
269 double UCons::SurfaceArea()
271 if (fSurfaceArea != 0.)
277 double mmin, mmax, dmin, dmax;
279 mmin = (fRmin1 + fRmin2) * 0.5;
280 mmax = (fRmax1 + fRmax2) * 0.5;
281 dmin = (fRmin2 - fRmin1);
282 dmax = (fRmax2 - fRmax1);
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));
290 fSurfaceArea = fSurfaceArea + 4 * fDz * (mmax - mmin);