ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ErrorCylSurfaceTarget.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ErrorCylSurfaceTarget.cc
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 // G4ErrorCylSurfaceTarget class implementation
27 //
28 // Created: P.Arce, September 2004
29 // --------------------------------------------------------------------
30 
32 #include "G4GeometryTolerance.hh"
33 
34 #ifdef G4VERBOSE
35 #include "G4ErrorPropagatorData.hh" // for verbosity checking
36 #endif
37 
38 #include "geomdefs.hh"
39 #include "G4Plane3D.hh"
40 
41 //---------------------------------------------------------------------
42 
45  const G4ThreeVector& trans,
46  const G4RotationMatrix& rotm )
47  : fradius(radius)
48 {
50 
51  ftransform = G4AffineTransform( rotm.inverse(), -trans );
52 #ifdef G4VERBOSE
54  {
55  Dump( " $$$ creating G4ErrorCylSurfaceTarget ");
56  }
57 #endif
58 }
59 
60 //---------------------------------------------------------------------
61 
64  const G4AffineTransform& trans )
65  : fradius(radius), ftransform(trans.Inverse())
66 {
68 
69 #ifdef G4VERBOSE
71  {
72  Dump( " $$$ creating G4ErrorCylSurfaceTarget ");
73  }
74 #endif
75 }
76 
77 //---------------------------------------------------------------------
78 
80 {
81 }
82 
83 //---------------------------------------------------------------------
84 
87  const G4ThreeVector& dir ) const
88 {
89  if( dir.mag() == 0. )
90  {
91  G4Exception("G4ErrorCylSurfaceTarget::GetDistanceFromPoint()",
92  "GeomMgt0003", FatalException, "Direction is zero !");
93  }
94 
95  //----- Get intersection point
96  G4ThreeVector localPoint = ftransform.TransformPoint( point );
97  G4ThreeVector localDir = ftransform.TransformAxis( dir );
98  G4ThreeVector inters = IntersectLocal(localPoint, localDir);
99 
100  G4double dist = (localPoint-inters).mag();
101 
102 #ifdef G4VERBOSE
104  {
105  G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint():" << G4endl
106  << " Global point " << point << " dir " << dir << G4endl
107  << " Intersection " << inters << G4endl
108  << " Distance " << dist << G4endl;
109  Dump( " CylSurface: " );
110  }
111 #endif
112 
113  return dist;
114 }
115 
116 
117 //---------------------------------------------------------------------
118 
121 {
122  G4ThreeVector localPoint = ftransform.TransformPoint( point );
123 
124 #ifdef G4VERBOSE
126  {
127  G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint:" << G4endl
128  << " Global point " << point << G4endl
129  << " Distance " << fradius - localPoint.perp() << G4endl;
130  Dump( " CylSurface: " );
131  }
132 #endif
133 
134  return fradius - localPoint.perp();
135 }
136 
137 //---------------------------------------------------------------------
138 
140 IntersectLocal( const G4ThreeVector& localPoint,
141  const G4ThreeVector& localDir ) const
142 {
143  G4double eqa = localDir.x()*localDir.x()+localDir.y()*localDir.y();
144  G4double eqb = 2*(localPoint.x()*localDir.x()+localPoint.y()*localDir.y());
145  G4double eqc = -fradius*fradius+localPoint.x()*localPoint.x()
146  +localPoint.y()*localPoint.y();
147  G4int inside = (localPoint.perp() > fradius) ? -1 : 1;
149 
150  if( eqa*inside > 0. )
151  {
152  lambda = (-eqb + std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa);
153  }
154  else if( eqa*inside < 0. )
155  {
156  lambda = (-eqb - std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa);
157  }
158  else
159  {
160  if( eqb != 0. )
161  {
162  lambda = -eqc/eqb;
163  }
164  else
165  {
166  std::ostringstream message;
167  message << "Intersection not possible !" << G4endl
168  << " Point: " << localPoint << ", direction: "
169  << localDir;
170  Dump( " CylSurface: " );
171  G4Exception("G4ErrorCylSurfaceTarget::IntersectLocal()",
172  "GeomMgt1002", JustWarning, message);
173  lambda = kInfinity;
174  }
175  }
176 
177  G4ThreeVector inters = localPoint + lambda*localDir/localDir.mag();
178 
179 #ifdef G4VERBOSE
180  if(G4ErrorPropagatorData::verbose() >= 4 ) {
181  G4cout << " G4ErrorCylSurfaceTarget::IntersectLocal " << inters << " "
182  << inters.perp() << " localPoint " << localPoint << " localDir "
183  << localDir << G4endl;
184  }
185 #endif
186 
187  return inters;
188 }
189 
190 //---------------------------------------------------------------------
191 
194 {
195  G4ThreeVector localPoint = ftransform.TransformPoint( point );
196 
197  // check that point is at cylinder surface
198  //
199  if( std::fabs( localPoint.perp() - fradius )
201  {
202  std::ostringstream message;
203  message << "Local point not at surface !" << G4endl
204  << " Point: " << point << ", local: " << localPoint
205  << G4endl
206  << " is not at surface, but far away by: "
207  << localPoint.perp() - fradius << " !";
208  G4Exception("G4ErrorCylSurfaceTarget::GetTangentPlane()",
209  "GeomMgt1002", JustWarning, message);
210  }
211 
212  G4Normal3D normal = localPoint - ftransform.NetTranslation();
213 
214  return G4Plane3D( normal, point );
215 }
216 
217 
218 //---------------------------------------------------------------------
219 
220 void G4ErrorCylSurfaceTarget::Dump( const G4String& msg ) const
221 {
222  G4cout << msg << " radius " << fradius
223  << " centre " << ftransform.NetTranslation()
224  << " rotation " << ftransform.NetRotation() << G4endl;
225 }