Geant4_10
G4ErrorCylSurfaceTarget.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4ErrorCylSurfaceTarget.cc 66356 2012-12-18 09:02:32Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class implementation file
32 // --------------------------------------------------------------------
33 
35 #include "G4GeometryTolerance.hh"
36 
37 #ifdef G4VERBOSE
38 #include "G4ErrorPropagatorData.hh" //for verbosity checking
39 #endif
40 
41 #include "geomdefs.hh"
42 #include "G4Plane3D.hh"
43 
44 //---------------------------------------------------------------------
45 
48  const G4ThreeVector& trans,
49  const G4RotationMatrix& rotm )
50  : fradius(radius)
51 {
53 
54  ftransform = G4AffineTransform( rotm.inverse(), -trans );
55 #ifdef G4VERBOSE
56  if(G4ErrorPropagatorData::verbose() >= 2 ) {
57  Dump( " $$$ creating G4ErrorCylSurfaceTarget ");
58  }
59 #endif
60 }
61 
62 //---------------------------------------------------------------------
63 
66  const G4AffineTransform& trans )
67  : fradius(radius), ftransform(trans.Inverse())
68 {
70 
71 #ifdef G4VERBOSE
73  {
74  Dump( " $$$ creating G4ErrorCylSurfaceTarget ");
75  }
76 #endif
77 }
78 
79 //---------------------------------------------------------------------
80 
82 {
83 }
84 
85 //---------------------------------------------------------------------
86 
89  const G4ThreeVector& dir ) const
90 {
91  if( dir.mag() == 0. )
92  {
93  G4Exception("G4ErrorCylSurfaceTarget::GetDistanceFromPoint()",
94  "GeomMgt0003", FatalException, "Direction is zero !");
95  }
96 
97  //----- Get intersection point
98  G4ThreeVector localPoint = ftransform.TransformPoint( point );
99  G4ThreeVector localDir = ftransform.TransformAxis( dir );
100  G4ThreeVector inters = IntersectLocal(localPoint, localDir);
101 
102  G4double dist = (localPoint-inters).mag();
103 
104 #ifdef G4VERBOSE
106  {
107  G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint():" << G4endl
108  << " Global point " << point << " dir " << dir << G4endl
109  << " Intersection " << inters << G4endl
110  << " Distance " << dist << G4endl;
111  Dump( " CylSurface: " );
112  }
113 #endif
114 
115  return dist;
116 }
117 
118 
119 //---------------------------------------------------------------------
120 
122 GetDistanceFromPoint( const G4ThreeVector& point ) const
123 {
124  G4ThreeVector localPoint = ftransform.TransformPoint( point );
125 
126 #ifdef G4VERBOSE
128  {
129  G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint:" << G4endl
130  << " Global point " << point << G4endl
131  << " Distance " << fradius - localPoint.perp() << G4endl;
132  Dump( " CylSurface: " );
133  }
134 #endif
135 
136  return fradius - localPoint.perp();
137 }
138 
139 //---------------------------------------------------------------------
140 
142 IntersectLocal( const G4ThreeVector& localPoint,
143  const G4ThreeVector& localDir ) const
144 {
145  G4double eqa = localDir.x()*localDir.x()+localDir.y()*localDir.y();
146  G4double eqb = 2*(localPoint.x()*localDir.x()+localPoint.y()*localDir.y());
147  G4double eqc = -fradius*fradius+localPoint.x()*localPoint.x()
148  +localPoint.y()*localPoint.y();
149  G4int inside = (localPoint.perp() > fradius) ? -1 : 1;
151 
152  if( eqa*inside > 0. )
153  {
154  lambda = (-eqb + std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa);
155  }
156  else if( eqa*inside < 0. )
157  {
158  lambda = (-eqb - std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa);
159  }
160  else
161  {
162  if( eqb != 0. )
163  {
164  lambda = -eqc/eqb;
165  }
166  else
167  {
168  std::ostringstream message;
169  message << "Intersection not possible !" << G4endl
170  << " Point: " << localPoint << ", direction: "
171  << localDir;
172  Dump( " CylSurface: " );
173  G4Exception("G4ErrorCylSurfaceTarget::IntersectLocal()",
174  "GeomMgt1002", JustWarning, message);
175  lambda = kInfinity;
176  }
177  }
178 
179  G4ThreeVector inters = localPoint + lambda*localDir/localDir.mag();
180 
181 #ifdef G4VERBOSE
182  if(G4ErrorPropagatorData::verbose() >= 4 ) {
183  G4cout << " G4ErrorCylSurfaceTarget::IntersectLocal " << inters << " "
184  << inters.perp() << " localPoint " << localPoint << " localDir "
185  << localDir << G4endl;
186  }
187 #endif
188 
189  return inters;
190 }
191 
192 //---------------------------------------------------------------------
193 
195 GetTangentPlane( const G4ThreeVector& point ) const
196 {
197  G4ThreeVector localPoint = ftransform.TransformPoint( point );
198 
199  // check that point is at cylinder surface
200  //
201  if( std::fabs( localPoint.perp() - fradius )
203  {
204  std::ostringstream message;
205  message << "Local point not at surface !" << G4endl
206  << " Point: " << point << ", local: " << localPoint
207  << G4endl
208  << " is not at surface, but far away by: "
209  << localPoint.perp() - fradius << " !";
210  G4Exception("G4ErrorCylSurfaceTarget::GetTangentPlane()",
211  "GeomMgt1002", JustWarning, message);
212  }
213 
214  G4Normal3D normal = localPoint - ftransform.NetTranslation();
215 
216  return G4Plane3D( normal, point );
217 }
218 
219 
220 //---------------------------------------------------------------------
221 
222 void G4ErrorCylSurfaceTarget::Dump( const G4String& msg ) const
223 {
224  G4cout << msg << " radius " << fradius
225  << " centre " << ftransform.NetTranslation()
226  << " rotation " << ftransform.NetRotation() << G4endl;
227 }
double x() const
G4double GetSurfaceTolerance() const
G4ThreeVector NetTranslation() const
int G4int
Definition: G4Types.hh:78
HepRotation inverse() const
G4ErrorCylSurfaceTarget(const G4double &radius, const G4ThreeVector &trans=G4ThreeVector(), const G4RotationMatrix &rotm=G4RotationMatrix())
G4GLOB_DLL std::ostream G4cout
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4RotationMatrix NetRotation() const
G4ErrorTargetType theType
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
virtual G4Plane3D GetTangentPlane(const G4ThreeVector &point) const
virtual G4ThreeVector IntersectLocal(const G4ThreeVector &point, const G4ThreeVector &direc) const
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
double y() const
HepGeom::Plane3D< G4double > G4Plane3D
Definition: G4Plane3D.hh:37
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
virtual void Dump(const G4String &msg) const
double perp() const
double mag() const
virtual G4double GetDistanceFromPoint(const G4ThreeVector &point, const G4ThreeVector &direc) const
static G4GeometryTolerance * GetInstance()
TDirectory * dir
Definition: macro.C:5