00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #include "G4ErrorCylSurfaceTarget.hh"
00035 #include "G4GeometryTolerance.hh"
00036
00037 #ifdef G4VERBOSE
00038 #include "G4ErrorPropagatorData.hh"
00039 #endif
00040
00041 #include "geomdefs.hh"
00042 #include "G4Plane3D.hh"
00043
00044
00045
00046 G4ErrorCylSurfaceTarget::
00047 G4ErrorCylSurfaceTarget( const G4double& radius,
00048 const G4ThreeVector& trans,
00049 const G4RotationMatrix& rotm )
00050 : fradius(radius)
00051 {
00052 theType = G4ErrorTarget_CylindricalSurface;
00053
00054 ftransform = G4AffineTransform( rotm.inverse(), -trans );
00055 #ifdef G4VERBOSE
00056 if(G4ErrorPropagatorData::verbose() >= 2 ) {
00057 Dump( " $$$ creating G4ErrorCylSurfaceTarget ");
00058 }
00059 #endif
00060 }
00061
00062
00063
00064 G4ErrorCylSurfaceTarget::
00065 G4ErrorCylSurfaceTarget( const G4double& radius,
00066 const G4AffineTransform& trans )
00067 : fradius(radius), ftransform(trans.Inverse())
00068 {
00069 theType = G4ErrorTarget_CylindricalSurface;
00070
00071 #ifdef G4VERBOSE
00072 if(G4ErrorPropagatorData::verbose() >= 2 )
00073 {
00074 Dump( " $$$ creating G4ErrorCylSurfaceTarget ");
00075 }
00076 #endif
00077 }
00078
00079
00080
00081 G4ErrorCylSurfaceTarget::~G4ErrorCylSurfaceTarget()
00082 {
00083 }
00084
00085
00086
00087 G4double G4ErrorCylSurfaceTarget::
00088 GetDistanceFromPoint( const G4ThreeVector& point,
00089 const G4ThreeVector& dir ) const
00090 {
00091 if( dir.mag() == 0. )
00092 {
00093 G4Exception("G4ErrorCylSurfaceTarget::GetDistanceFromPoint()",
00094 "GeomMgt0003", FatalException, "Direction is zero !");
00095 }
00096
00097
00098 G4ThreeVector localPoint = ftransform.TransformPoint( point );
00099 G4ThreeVector localDir = ftransform.TransformAxis( dir );
00100 G4ThreeVector inters = IntersectLocal(localPoint, localDir);
00101
00102 G4double dist = (localPoint-inters).mag();
00103
00104 #ifdef G4VERBOSE
00105 if(G4ErrorPropagatorData::verbose() >= 3 )
00106 {
00107 G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint():" << G4endl
00108 << " Global point " << point << " dir " << dir << G4endl
00109 << " Intersection " << inters << G4endl
00110 << " Distance " << dist << G4endl;
00111 Dump( " CylSurface: " );
00112 }
00113 #endif
00114
00115 return dist;
00116 }
00117
00118
00119
00120
00121 G4double G4ErrorCylSurfaceTarget::
00122 GetDistanceFromPoint( const G4ThreeVector& point ) const
00123 {
00124 G4ThreeVector localPoint = ftransform.TransformPoint( point );
00125
00126 #ifdef G4VERBOSE
00127 if(G4ErrorPropagatorData::verbose() >= 3 )
00128 {
00129 G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint:" << G4endl
00130 << " Global point " << point << G4endl
00131 << " Distance " << fradius - localPoint.perp() << G4endl;
00132 Dump( " CylSurface: " );
00133 }
00134 #endif
00135
00136 return fradius - localPoint.perp();
00137 }
00138
00139
00140
00141 G4ThreeVector G4ErrorCylSurfaceTarget::
00142 IntersectLocal( const G4ThreeVector& localPoint,
00143 const G4ThreeVector& localDir ) const
00144 {
00145 G4double eqa = localDir.x()*localDir.x()+localDir.y()*localDir.y();
00146 G4double eqb = 2*(localPoint.x()*localDir.x()+localPoint.y()*localDir.y());
00147 G4double eqc = -fradius*fradius+localPoint.x()*localPoint.x()
00148 +localPoint.y()*localPoint.y();
00149 G4int inside = (localPoint.perp() > fradius) ? -1 : 1;
00150 G4double lambda;
00151
00152 if( eqa*inside > 0. )
00153 {
00154 lambda = (-eqb + std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa);
00155 }
00156 else if( eqa*inside < 0. )
00157 {
00158 lambda = (-eqb - std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa);
00159 }
00160 else
00161 {
00162 if( eqb != 0. )
00163 {
00164 lambda = -eqc/eqb;
00165 }
00166 else
00167 {
00168 std::ostringstream message;
00169 message << "Intersection not possible !" << G4endl
00170 << " Point: " << localPoint << ", direction: "
00171 << localDir;
00172 Dump( " CylSurface: " );
00173 G4Exception("G4ErrorCylSurfaceTarget::IntersectLocal()",
00174 "GeomMgt1002", JustWarning, message);
00175 lambda = kInfinity;
00176 }
00177 }
00178
00179 G4ThreeVector inters = localPoint + lambda*localDir/localDir.mag();
00180
00181 #ifdef G4VERBOSE
00182 if(G4ErrorPropagatorData::verbose() >= 4 ) {
00183 G4cout << " G4ErrorCylSurfaceTarget::IntersectLocal " << inters << " "
00184 << inters.perp() << " localPoint " << localPoint << " localDir "
00185 << localDir << G4endl;
00186 }
00187 #endif
00188
00189 return inters;
00190 }
00191
00192
00193
00194 G4Plane3D G4ErrorCylSurfaceTarget::
00195 GetTangentPlane( const G4ThreeVector& point ) const
00196 {
00197 G4ThreeVector localPoint = ftransform.TransformPoint( point );
00198
00199
00200
00201 if( std::fabs( localPoint.perp() - fradius )
00202 > 1000.*G4GeometryTolerance::GetInstance()->GetSurfaceTolerance() )
00203 {
00204 std::ostringstream message;
00205 message << "Local point not at surface !" << G4endl
00206 << " Point: " << point << ", local: " << localPoint
00207 << G4endl
00208 << " is not at surface, but far away by: "
00209 << localPoint.perp() - fradius << " !";
00210 G4Exception("G4ErrorCylSurfaceTarget::GetTangentPlane()",
00211 "GeomMgt1002", JustWarning, message);
00212 }
00213
00214 G4Normal3D normal = localPoint - ftransform.NetTranslation();
00215
00216 return G4Plane3D( normal, point );
00217 }
00218
00219
00220
00221
00222 void G4ErrorCylSurfaceTarget::Dump( const G4String& msg ) const
00223 {
00224 G4cout << msg << " radius " << fradius
00225 << " centre " << ftransform.NetTranslation()
00226 << " rotation " << ftransform.NetRotation() << G4endl;
00227 }