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 #include "G4DELPHIMagField.hh"
00031 #include "globals.hh"
00032 #include "G4SystemOfUnits.hh"
00033
00034 G4DELPHIMagField::G4DELPHIMagField()
00035 {
00036 }
00037
00039
00040 G4DELPHIMagField::~G4DELPHIMagField()
00041 {
00042 }
00043
00045
00046
00047 void G4DELPHIMagField::GetFieldValue( const G4double yTrack[7],
00048 G4double B[3] ) const
00049 {
00050 G4int i, n = 8 ;
00051 G4double a = 0.001 ;
00052 G4double x = a*yTrack[0], y = a*yTrack[1], z = a*yTrack[2] ;
00053 G4double x2 = x*x, y2 = y*y, z2 = z*z, r2 = x2 + y2 ;
00054 G4double r4 = r2*r2, z4 = z2*z2, r6 = r4*r2, z6 = z4*z2 ;
00055 G4double r8 = r4*r4, z8 = z4*z4, r10 = r8*r2, z10 = z8*z2 ;
00056 G4double rz = z*std::sqrt(r2), r = std::sqrt(r2+a*a) ;
00057 G4double Br ;
00058 G4double P[8], Q[8] ;
00059 static G4double c[8] = {
00060 -9.26e-5, -3.51e-5, 2.94e-6, -1.10e-6,
00061 6.25e-8, -1.77e-8, -6.88e-10, -7.52e-11
00062 } ;
00063 P[0] = 2*rz ;
00064 P[1] = 4*rz*(r2 - 4.0*z2/3.0) ;
00065 P[2] = 6*rz*(r4 - 4*r2*z2 + 1.6*z4) ;
00066 P[3] = 8*rz*(r6 - 8*r4*z2 + 9.6*r2*z4 - 64.0*z6/35.0) ;
00067 P[4] = 10*rz*(r8 - 40.0*r6*z2/3.0 + 32*r4*z4
00068 - 128.0*r2*z6/7.0 + 128.0*z8/63.0);
00069 P[5] = 0 ;
00070 P[6] = 0 ;
00071 P[7] = 0 ;
00072
00073 Q[0] = r2 - 2*z2 ;
00074 Q[1] = r4 - 8*r2*z2 + 8.0*z4/3.0 ;
00075 Q[2] = r6 - 18*r4*z2 + 24*r2*z4 - 3.2*z6 ;
00076 Q[3] = r8 - 32*r6*z2 + 96*r4*z4 - 51.2*r2*z6 +128.0*z8/35.0 ;
00077 Q[4] = r10 - 50*r8*z2 + 800.0*r6*z4/3.0 - 320*r4*z6
00078 + 640.0*r2*z8/7.0 - 256.0*z10/63.0 ;
00079 Q[5] = 0 ;
00080 Q[6] = 0 ;
00081 Q[7] = 0 ;
00082
00083 Br = 0 ;
00084 B[2] = 1.2*tesla ;
00085 for(i=0;i<n;i++)
00086 {
00087 Br += c[i]*P[i] ;
00088 B[2] += c[i]*Q[i] ;
00089 }
00090 B[0] = Br*x/r ;
00091 B[1] = Br*y/r ;
00092 return ;
00093 }
00094
00095