Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // G4ErrorCylSurfaceTarget class implementatio << 26 // 27 // << 27 // $Id: G4ErrorCylSurfaceTarget.cc,v 1.4 2007/06/20 12:50:48 arce Exp $ 28 // Created: P.Arce, September 2004 << 28 // GEANT4 tag $Name: geant4-09-03-patch-02 $ >> 29 // >> 30 // >> 31 // -------------------------------------------------------------------- >> 32 // GEANT 4 class implementation file 29 // ------------------------------------------- 33 // -------------------------------------------------------------------- 30 34 31 #include "G4ErrorCylSurfaceTarget.hh" 35 #include "G4ErrorCylSurfaceTarget.hh" 32 #include "G4GeometryTolerance.hh" 36 #include "G4GeometryTolerance.hh" 33 37 34 #ifdef G4VERBOSE 38 #ifdef G4VERBOSE 35 #include "G4ErrorPropagatorData.hh" // for ver << 39 #include "G4ErrorPropagatorData.hh" //for verbosity checking 36 #endif 40 #endif 37 41 38 #include "geomdefs.hh" 42 #include "geomdefs.hh" 39 #include "G4Normal3D.hh" << 40 #include "G4Plane3D.hh" 43 #include "G4Plane3D.hh" 41 44 42 //-------------------------------------------- 45 //--------------------------------------------------------------------- 43 46 44 G4ErrorCylSurfaceTarget:: 47 G4ErrorCylSurfaceTarget:: 45 G4ErrorCylSurfaceTarget( const G4double& radiu 48 G4ErrorCylSurfaceTarget( const G4double& radius, 46 const G4ThreeVector& 49 const G4ThreeVector& trans, 47 const G4RotationMatri 50 const G4RotationMatrix& rotm ) 48 : fradius(radius) 51 : fradius(radius) 49 { 52 { 50 theType = G4ErrorTarget_CylindricalSurface; 53 theType = G4ErrorTarget_CylindricalSurface; 51 54 52 ftransform = G4AffineTransform( rotm.inverse 55 ftransform = G4AffineTransform( rotm.inverse(), -trans ); 53 #ifdef G4VERBOSE 56 #ifdef G4VERBOSE 54 if(G4ErrorPropagatorData::verbose() >= 2 ) << 57 if(G4ErrorPropagatorData::verbose() >= 2 ) { 55 { << 56 Dump( " $$$ creating G4ErrorCylSurfaceTarg 58 Dump( " $$$ creating G4ErrorCylSurfaceTarget "); 57 } 59 } 58 #endif 60 #endif 59 } 61 } 60 62 61 //-------------------------------------------- 63 //--------------------------------------------------------------------- 62 64 63 G4ErrorCylSurfaceTarget:: 65 G4ErrorCylSurfaceTarget:: 64 G4ErrorCylSurfaceTarget( const G4double& radiu 66 G4ErrorCylSurfaceTarget( const G4double& radius, 65 const G4AffineTransfo 67 const G4AffineTransform& trans ) 66 : fradius(radius), ftransform(trans.Inverse( 68 : fradius(radius), ftransform(trans.Inverse()) 67 { 69 { 68 theType = G4ErrorTarget_CylindricalSurface; 70 theType = G4ErrorTarget_CylindricalSurface; 69 71 70 #ifdef G4VERBOSE 72 #ifdef G4VERBOSE 71 if(G4ErrorPropagatorData::verbose() >= 2 ) 73 if(G4ErrorPropagatorData::verbose() >= 2 ) 72 { 74 { 73 Dump( " $$$ creating G4ErrorCylSurfaceTarg 75 Dump( " $$$ creating G4ErrorCylSurfaceTarget "); 74 } 76 } 75 #endif 77 #endif 76 } 78 } 77 79 78 //-------------------------------------------- 80 //--------------------------------------------------------------------- 79 81 >> 82 G4ErrorCylSurfaceTarget::~G4ErrorCylSurfaceTarget() >> 83 { >> 84 } >> 85 >> 86 //--------------------------------------------------------------------- >> 87 80 G4double G4ErrorCylSurfaceTarget:: 88 G4double G4ErrorCylSurfaceTarget:: 81 GetDistanceFromPoint( const G4ThreeVector& poi 89 GetDistanceFromPoint( const G4ThreeVector& point, 82 const G4ThreeVector& dir 90 const G4ThreeVector& dir ) const 83 { 91 { 84 if( dir.mag() == 0. ) 92 if( dir.mag() == 0. ) 85 { 93 { 86 G4Exception("G4ErrorCylSurfaceTarget::GetD 94 G4Exception("G4ErrorCylSurfaceTarget::GetDistanceFromPoint()", 87 "GeomMgt0003", FatalException, << 95 "InvalidSetup", FatalException, "Direction is zero !"); 88 } 96 } 89 97 90 //----- Get intersection point 98 //----- Get intersection point 91 G4ThreeVector localPoint = ftransform.Transf 99 G4ThreeVector localPoint = ftransform.TransformPoint( point ); 92 G4ThreeVector localDir = ftransform.Transfor 100 G4ThreeVector localDir = ftransform.TransformAxis( dir ); 93 G4ThreeVector inters = IntersectLocal(localP 101 G4ThreeVector inters = IntersectLocal(localPoint, localDir); 94 102 95 G4double dist = (localPoint-inters).mag(); 103 G4double dist = (localPoint-inters).mag(); 96 104 97 #ifdef G4VERBOSE 105 #ifdef G4VERBOSE 98 if(G4ErrorPropagatorData::verbose() >= 3 ) 106 if(G4ErrorPropagatorData::verbose() >= 3 ) 99 { 107 { 100 G4cout << " G4ErrorCylSurfaceTarget::GetDi 108 G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint():" << G4endl 101 << " Global point " << point << " d 109 << " Global point " << point << " dir " << dir << G4endl 102 << " Intersection " << inters << G4 110 << " Intersection " << inters << G4endl 103 << " Distance " << dist << G4endl; 111 << " Distance " << dist << G4endl; 104 Dump( " CylSurface: " ); 112 Dump( " CylSurface: " ); 105 } 113 } 106 #endif 114 #endif 107 115 108 return dist; 116 return dist; 109 } 117 } 110 118 111 119 112 //-------------------------------------------- 120 //--------------------------------------------------------------------- 113 121 114 G4double G4ErrorCylSurfaceTarget:: 122 G4double G4ErrorCylSurfaceTarget:: 115 GetDistanceFromPoint( const G4ThreeVector& poi 123 GetDistanceFromPoint( const G4ThreeVector& point ) const 116 { 124 { 117 G4ThreeVector localPoint = ftransform.Transf 125 G4ThreeVector localPoint = ftransform.TransformPoint( point ); 118 126 119 #ifdef G4VERBOSE 127 #ifdef G4VERBOSE 120 if(G4ErrorPropagatorData::verbose() >= 3 ) 128 if(G4ErrorPropagatorData::verbose() >= 3 ) 121 { 129 { 122 G4cout << " G4ErrorCylSurfaceTarget::GetDi 130 G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint:" << G4endl 123 << " Global point " << point << G4e 131 << " Global point " << point << G4endl 124 << " Distance " << fradius - localP 132 << " Distance " << fradius - localPoint.perp() << G4endl; 125 Dump( " CylSurface: " ); 133 Dump( " CylSurface: " ); 126 } 134 } 127 #endif 135 #endif 128 136 129 return fradius - localPoint.perp(); 137 return fradius - localPoint.perp(); 130 } 138 } 131 139 132 //-------------------------------------------- 140 //--------------------------------------------------------------------- 133 141 134 G4ThreeVector G4ErrorCylSurfaceTarget:: 142 G4ThreeVector G4ErrorCylSurfaceTarget:: 135 IntersectLocal( const G4ThreeVector& localPoin 143 IntersectLocal( const G4ThreeVector& localPoint, 136 const G4ThreeVector& localDir 144 const G4ThreeVector& localDir ) const 137 { 145 { 138 G4double eqa = localDir.x()*localDir.x()+loc 146 G4double eqa = localDir.x()*localDir.x()+localDir.y()*localDir.y(); 139 G4double eqb = 2*(localPoint.x()*localDir.x( 147 G4double eqb = 2*(localPoint.x()*localDir.x()+localPoint.y()*localDir.y()); 140 G4double eqc = -fradius*fradius+localPoint.x 148 G4double eqc = -fradius*fradius+localPoint.x()*localPoint.x() 141 +localPoint.y()*localPoint.y( 149 +localPoint.y()*localPoint.y(); 142 G4int inside = (localPoint.perp() > fradius) 150 G4int inside = (localPoint.perp() > fradius) ? -1 : 1; 143 G4double lambda; 151 G4double lambda; 144 152 145 if( eqa*inside > 0. ) 153 if( eqa*inside > 0. ) 146 { 154 { 147 lambda = (-eqb + std::sqrt(eqb*eqb-4*eqa*e 155 lambda = (-eqb + std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa); 148 } 156 } 149 else if( eqa*inside < 0. ) 157 else if( eqa*inside < 0. ) 150 { 158 { 151 lambda = (-eqb - std::sqrt(eqb*eqb-4*eqa*e 159 lambda = (-eqb - std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa); 152 } 160 } 153 else 161 else 154 { 162 { 155 if( eqb != 0. ) 163 if( eqb != 0. ) 156 { 164 { 157 lambda = -eqc/eqb; 165 lambda = -eqc/eqb; 158 } 166 } 159 else 167 else 160 { 168 { 161 std::ostringstream message; << 169 G4cerr << "WARNING - G4ErrorCylSurfaceTarget::IntersectLocal()" << G4endl 162 message << "Intersection not possible !" << 170 << " Point: " << localPoint << ", direction: " 163 << " Point: " << localP << 171 << localDir << G4endl; 164 << localDir; << 165 Dump( " CylSurface: " ); 172 Dump( " CylSurface: " ); 166 G4Exception("G4ErrorCylSurfaceTarget::In 173 G4Exception("G4ErrorCylSurfaceTarget::IntersectLocal()", 167 "GeomMgt1002", JustWarning, << 174 "IllegalCondition", >> 175 JustWarning, "Intersection not possible !"); 168 lambda = kInfinity; 176 lambda = kInfinity; 169 } 177 } 170 } 178 } 171 179 172 G4ThreeVector inters = localPoint + lambda*l 180 G4ThreeVector inters = localPoint + lambda*localDir/localDir.mag(); 173 181 174 #ifdef G4VERBOSE 182 #ifdef G4VERBOSE 175 if(G4ErrorPropagatorData::verbose() >= 4 ) { 183 if(G4ErrorPropagatorData::verbose() >= 4 ) { 176 G4cout << " G4ErrorCylSurfaceTarget::Inter << 184 G4cout << " G4ErrorCylSurfaceTarget::IntersectLocal " << inters << " " << inters.perp() << " localPoint " << localPoint << " localDir " << localDir << G4endl; 177 << inters.perp() << " localPoint " << 178 << localDir << G4endl; << 179 } 185 } 180 #endif 186 #endif 181 187 182 return inters; 188 return inters; 183 } 189 } 184 190 185 //-------------------------------------------- 191 //--------------------------------------------------------------------- 186 192 187 G4Plane3D G4ErrorCylSurfaceTarget:: 193 G4Plane3D G4ErrorCylSurfaceTarget:: 188 GetTangentPlane( const G4ThreeVector& point ) 194 GetTangentPlane( const G4ThreeVector& point ) const 189 { 195 { 190 G4ThreeVector localPoint = ftransform.Transf 196 G4ThreeVector localPoint = ftransform.TransformPoint( point ); 191 197 192 // check that point is at cylinder surface 198 // check that point is at cylinder surface 193 // 199 // 194 if( std::fabs( localPoint.perp() - fradius ) 200 if( std::fabs( localPoint.perp() - fradius ) 195 > 1000.*G4GeometryTolerance::GetInstance 201 > 1000.*G4GeometryTolerance::GetInstance()->GetSurfaceTolerance() ) 196 { 202 { 197 std::ostringstream message; << 203 G4cerr << "WARNING - G4ErrorCylSurfaceTarget::GetTangentPlane()" << G4endl 198 message << "Local point not at surface !" << 204 << " Point: " << point << ", local: " << localPoint 199 << " Point: " << point << << 205 << G4endl 200 << G4endl << 206 << " is not at surface, but far away by: " 201 << " is not at surface, b << 207 << localPoint.perp() - fradius << " !" << G4endl; 202 << localPoint.perp() - fradius << << 203 G4Exception("G4ErrorCylSurfaceTarget::GetT 208 G4Exception("G4ErrorCylSurfaceTarget::GetTangentPlane()", 204 "GeomMgt1002", JustWarning, me << 209 "IllegalCondition", JustWarning, >> 210 "Local point not at surface !"); 205 } 211 } 206 212 207 G4Normal3D normal = localPoint - ftransform. 213 G4Normal3D normal = localPoint - ftransform.NetTranslation(); 208 214 209 return G4Plane3D( normal, point ); 215 return G4Plane3D( normal, point ); 210 } 216 } 211 217 212 218 213 //-------------------------------------------- 219 //--------------------------------------------------------------------- 214 220 215 void G4ErrorCylSurfaceTarget::Dump( const G4St 221 void G4ErrorCylSurfaceTarget::Dump( const G4String& msg ) const 216 { 222 { 217 G4cout << msg << " radius " << fradius 223 G4cout << msg << " radius " << fradius 218 << " centre " << ftransform.Ne 224 << " centre " << ftransform.NetTranslation() 219 << " rotation " << ftransform. 225 << " rotation " << ftransform.NetRotation() << G4endl; 220 } 226 } 221 227