Belle II Software  release-05-01-25
ExtCylSurfaceTarget.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Leo Piilonen *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <simulation/kernel/ExtCylSurfaceTarget.h>
12 #include <G4GeometryTolerance.hh>
13 #include <geomdefs.hh>
14 #include <G4Plane3D.hh>
15 #include <framework/logging/Logger.h>
16 
17 using namespace std;
18 using namespace Belle2;
19 using namespace Belle2::Simulation;
20 
21 ExtCylSurfaceTarget::ExtCylSurfaceTarget(const G4double& radius,
22  const G4double& zmin,
23  const G4double& zmax) :
24  m_radius(radius),
25  m_zmin(zmin),
26  m_zmax(zmax),
27  m_tolerance(1000.0 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance())
28 {
29  theType = G4ErrorTarget_CylindricalSurface;
30  B2DEBUG(1, "Simulation::ExtCylSurfaceTarget created with radius "
31  << m_radius << " zmin " << zmin << " zmax " << zmax);
32 }
33 
35 {
36 }
37 
38 G4double ExtCylSurfaceTarget::GetDistanceFromPoint(const G4ThreeVector& point,
39  const G4ThreeVector& dir) const
40 {
41  if (fabs(dir.mag() - 1.0) > 1.0E-10) {
42  B2FATAL("Simulation::ExtCylSurfaceTarget::GetDistanceFromPoint() direction is not a unit vector: " << dir);
43  }
44 
45  // Get distance to intersection point with the cylinder's curved surface
46  // should be negative if outside! G4double dist = (point - IntersectLocal(point, dir)).mag();
47  G4double dist = (IntersectLocal(point, dir) - point) * dir;
48 
49  // Get intersection point with the plane at either zmin or zmax
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);
55  }
56 
57  B2DEBUG(300, "Simulation::ExtCylSurfaceTarget::GetDistanceFromPoint(): Global point "
58  << point << " dir " << dir << " distance " << dist);
59 
60  return dist;
61 }
62 
63 G4double ExtCylSurfaceTarget::GetDistanceFromPoint(const G4ThreeVector& point) const
64 {
65 
66  G4double dist = m_radius - point.perp();
67  dist = fmin(dist, point.z() - m_zmin);
68  dist = fmin(dist, m_zmax - point.z());
69 
70  B2DEBUG(300, "Simulation::ExtCylSurfaceTarget::GetDistanceFromPoint(): Global point "
71  << point << " minimum distance " << dist);
72 
73  return dist;
74 }
75 
76 G4ThreeVector ExtCylSurfaceTarget::IntersectLocal(const G4ThreeVector& localPoint,
77  const G4ThreeVector& localDir) const
78 {
79  // localDir has already been verified to be a unit vector
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()
83  - m_radius * m_radius;
84  G4double eqaInside = (localPoint.perp() > m_radius) ? -eqa : eqa;
85 
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));
91  } else {
92  if (eqb != 0.0) {
93  intersection -= localDir * (eqc / eqb);
94  } else {
95  B2WARNING("Simulation::ExtCylSurfaceTarget::IntersectLocal(): localPoint "
96  << localPoint << " localDir " << localDir << " does not intersect with cylinder");
97  intersection = localDir * kInfinity;
98  }
99  }
100 
101  B2DEBUG(300, "Simulation::ExtCylSurfaceTarget::IntersectLocal(): localPoint "
102  << localPoint << " localDir " << localDir << " radial intersection " << intersection.perp());
103 
104  return intersection;
105 }
106 
107 G4Plane3D ExtCylSurfaceTarget::GetTangentPlane(const G4ThreeVector& point) const
108 {
109  // check that point is at the cylinder's curved surface
110  if (fabs(point.perp() - m_radius) > m_tolerance) {
111  B2ERROR("Simulation::ExtCylSurfaceTarget::GetTangentPlane(): point "
112  << point << " is not at surface; radial distance is " << point.perp() - m_radius);
113  }
114 
115  G4Normal3D normal(point);
116 
117  return G4Plane3D(normal, point);
118 }
119 
121 void ExtCylSurfaceTarget::Dump(const G4String&) const
122 {
123 }
124 
Belle2::Simulation::ExtCylSurfaceTarget::m_zmax
G4double m_zmax
Cylinder maximum-z coordinate.
Definition: ExtCylSurfaceTarget.h:76
Belle2::Simulation::ExtCylSurfaceTarget::Dump
virtual void Dump(const G4String &msg) const
Dump the cylinder parameters.
Definition: ExtCylSurfaceTarget.cc:121
Belle2::Simulation::ExtCylSurfaceTarget::GetDistanceFromPoint
virtual G4double GetDistanceFromPoint(const G4ThreeVector &point, const G4ThreeVector &direc) const
Get the distance from a point to the cylinder along direc.
Definition: ExtCylSurfaceTarget.cc:38
Belle2::Simulation::ExtCylSurfaceTarget::m_radius
G4double m_radius
Cylinder radius.
Definition: ExtCylSurfaceTarget.h:70
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::Simulation::ExtCylSurfaceTarget::~ExtCylSurfaceTarget
~ExtCylSurfaceTarget()
Destructor.
Definition: ExtCylSurfaceTarget.cc:34
Belle2::Simulation::ExtCylSurfaceTarget::GetTangentPlane
virtual G4Plane3D GetTangentPlane(const G4ThreeVector &point) const
Get the plane tangent to the cylinder at a given point.
Definition: ExtCylSurfaceTarget.cc:107
Belle2::Simulation::ExtCylSurfaceTarget::m_tolerance
G4double m_tolerance
Tolerance for distance between a point and cylinder's curved surface.
Definition: ExtCylSurfaceTarget.h:79
Belle2::Simulation::ExtCylSurfaceTarget::m_zmin
G4double m_zmin
Cylinder minimum-z coordinate.
Definition: ExtCylSurfaceTarget.h:73
Belle2::Simulation::ExtCylSurfaceTarget::IntersectLocal
virtual G4ThreeVector IntersectLocal(const G4ThreeVector &point, const G4ThreeVector &direc) const
Return the intersection of the cylinder with the line defined in local (cylinder) coordinates by poin...
Definition: ExtCylSurfaceTarget.cc:76