11 #include <simulation/kernel/ExtCylSurfaceTarget.h>
12 #include <G4GeometryTolerance.hh>
13 #include <geomdefs.hh>
14 #include <G4Plane3D.hh>
15 #include <framework/logging/Logger.h>
19 using namespace Belle2::Simulation;
21 ExtCylSurfaceTarget::ExtCylSurfaceTarget(
const G4double& radius,
23 const G4double& zmax) :
27 m_tolerance(1000.0 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance())
29 theType = G4ErrorTarget_CylindricalSurface;
30 B2DEBUG(1,
"Simulation::ExtCylSurfaceTarget created with radius "
31 <<
m_radius <<
" zmin " << zmin <<
" zmax " << zmax);
39 const G4ThreeVector& dir)
const
41 if (fabs(dir.mag() - 1.0) > 1.0E-10) {
42 B2FATAL(
"Simulation::ExtCylSurfaceTarget::GetDistanceFromPoint() direction is not a unit vector: " << dir);
50 G4double dirz = dir.z();
51 if (dirz < -1.0E-10) {
52 dist = fmin(dist, (
m_zmin - point.z()) / dirz);
53 }
else if (dirz > 1.0E-10) {
54 dist = fmin(dist, (
m_zmax - point.z()) / dirz);
57 B2DEBUG(300,
"Simulation::ExtCylSurfaceTarget::GetDistanceFromPoint(): Global point "
58 << point <<
" dir " << dir <<
" distance " << dist);
66 G4double dist =
m_radius - point.perp();
67 dist = fmin(dist, point.z() -
m_zmin);
68 dist = fmin(dist,
m_zmax - point.z());
70 B2DEBUG(300,
"Simulation::ExtCylSurfaceTarget::GetDistanceFromPoint(): Global point "
71 << point <<
" minimum distance " << dist);
77 const G4ThreeVector& localDir)
const
80 G4double eqa = localDir.x() * localDir.x() + localDir.y() * localDir.y();
81 G4double eqb = 2.0 * (localPoint.x() * localDir.x() + localPoint.y() * localDir.y());
82 G4double eqc = localPoint.x() * localPoint.x() + localPoint.y() * localPoint.y()
84 G4double eqaInside = (localPoint.perp() >
m_radius) ? -eqa : eqa;
86 G4ThreeVector intersection = localPoint;
87 if (eqaInside > 0.0) {
88 intersection += localDir * ((-eqb + sqrt(eqb * eqb - 4.0 * eqa * eqc)) / (2.0 * eqa));
89 }
else if (eqaInside < 0.0) {
90 intersection += localDir * ((-eqb - sqrt(eqb * eqb - 4.0 * eqa * eqc)) / (2.0 * eqa));
93 intersection -= localDir * (eqc / eqb);
95 B2WARNING(
"Simulation::ExtCylSurfaceTarget::IntersectLocal(): localPoint "
96 << localPoint <<
" localDir " << localDir <<
" does not intersect with cylinder");
97 intersection = localDir * kInfinity;
101 B2DEBUG(300,
"Simulation::ExtCylSurfaceTarget::IntersectLocal(): localPoint "
102 << localPoint <<
" localDir " << localDir <<
" radial intersection " << intersection.perp());
111 B2ERROR(
"Simulation::ExtCylSurfaceTarget::GetTangentPlane(): point "
112 << point <<
" is not at surface; radial distance is " << point.perp() -
m_radius);
115 G4Normal3D normal(point);
117 return G4Plane3D(normal, point);