00001
00002
00003
00004
00005
00006
00007 #include "BDSDecStepper.hh"
00008 #include "G4ThreeVector.hh"
00009 #include "G4TransportationManager.hh"
00010
00011 extern G4double BDSLocalRadiusOfCurvature;
00012
00013 BDSDecStepper::BDSDecStepper(G4Mag_EqRhs *EqRhs)
00014 : G4MagIntegratorStepper(EqRhs,6),
00015
00016 itsBQuadPrime(0.0), itsDist(0.0)
00017 {
00018 fPtrMagEqOfMot = EqRhs;
00019 }
00020
00021
00022 void BDSDecStepper::AdvanceHelix( const G4double yIn[],
00023 G4ThreeVector,
00024 G4double h,
00025 G4double yDec[])
00026 {
00027 const G4double *pIn = yIn+3;
00028 G4ThreeVector v0= G4ThreeVector( pIn[0], pIn[1], pIn[2]);
00029 G4ThreeVector InitMomDir=v0.unit();
00030
00031 G4ThreeVector GlobalPosition= G4ThreeVector( yIn[0], yIn[1], yIn[2]);
00032 G4double InitMag=v0.mag();
00033 G4double kappa= -fPtrMagEqOfMot->FCof()*itsBQuadPrime/InitMag;
00034
00035
00036
00037 if(fabs(kappa)<1.e-20)
00038 {
00039 G4ThreeVector positionMove = (h/InitMag) * v0;
00040
00041 yDec[0] = yIn[0] + positionMove.x();
00042 yDec[1] = yIn[1] + positionMove.y();
00043 yDec[2] = yIn[2] + positionMove.z();
00044
00045 yDec[3] = v0.x();
00046 yDec[4] = v0.y();
00047 yDec[5] = v0.z();
00048
00049 itsDist=0;
00050 }
00051 else
00052 {
00053 G4Navigator* DecNavigator=
00054 G4TransportationManager::GetTransportationManager()->
00055 GetNavigatorForTracking();
00056
00057 G4AffineTransform LocalAffine=DecNavigator->GetLocalToGlobalTransform();
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067 G4AffineTransform GlobalAffine=DecNavigator->GetGlobalToLocalTransform();
00068 G4ThreeVector LocalR=GlobalAffine.TransformPoint(GlobalPosition);
00069 G4ThreeVector LocalRp=GlobalAffine.TransformAxis(InitMomDir);
00070
00071
00072
00073 G4double x0=LocalR.x();
00074 G4double y0=LocalR.y();
00075
00076 G4double x02My02=(x0*x0-y0*y0);
00077
00078 G4double xp=LocalRp.x();
00079 G4double yp=LocalRp.y();
00080 G4double zp=LocalRp.z();
00081
00082 G4double Bx=4.0*x0*y0*(x02My02);
00083 G4double By=pow(x0,4)-6.0*x0*x0*y0*y0+pow(y0,4);
00084
00085
00086 G4ThreeVector LocalRpp;
00087
00088 LocalRpp.setX( zp*By);
00089 LocalRpp.setY(-zp*Bx);
00090 LocalRpp.setZ( xp*By - yp*Bx);
00091
00092 LocalRpp*=kappa/24;
00093
00094
00095 G4double R_1 = LocalRpp.mag();
00096 if(R_1>0.)
00097 {
00098
00099 BDSLocalRadiusOfCurvature=1/R_1;
00100
00101
00102 G4double h2=h*h;
00103 itsDist= h2*R_1/8;
00104
00105 G4double dx=LocalRp.x()*h + LocalRpp.x()*h2/2;
00106 G4double dy=LocalRp.y()*h + LocalRpp.y()*h2/2;
00107
00108 G4double dz=sqrt(h2*(1.-h2*R_1*R_1/12)-dx*dx-dy*dy);
00109
00110 G4double ScaleFac=(dx*dx+dy*dy+dz*dz)/h2;
00111 if(ScaleFac>1.0000001)
00112 {
00113 ScaleFac=sqrt(ScaleFac);
00114 dx/=ScaleFac;
00115 dy/=ScaleFac;
00116 dz/=ScaleFac;
00117 }
00118
00119 LocalR.setX(LocalR.x()+dx);
00120 LocalR.setY(LocalR.y()+dy);
00121 LocalR.setZ(LocalR.z()+dz);
00122
00123 LocalRp = LocalRp+ h*LocalRpp;
00124 }
00125 else
00126 {LocalR += h*LocalRp;}
00127
00128 GlobalPosition=LocalAffine.TransformPoint(LocalR);
00129 G4ThreeVector GlobalTangent=LocalAffine.TransformAxis(LocalRp)*InitMag;
00130
00131 yDec[0] = GlobalPosition.x();
00132 yDec[1] = GlobalPosition.y();
00133 yDec[2] = GlobalPosition.z();
00134
00135 yDec[3] = GlobalTangent.x();
00136 yDec[4] = GlobalTangent.y();
00137 yDec[5] = GlobalTangent.z();
00138 }
00139 }
00140
00141 void BDSDecStepper::Stepper( const G4double yInput[],
00142 const G4double[],
00143 const G4double hstep,
00144 G4double yOut[],
00145 G4double yErr[] )
00146 {
00147 const G4int nvar = 6 ;
00148
00149 G4int i;
00150 for(i=0;i<nvar;i++) yErr[i]=0;
00151 AdvanceHelix(yInput,(G4ThreeVector)0,hstep,yOut);
00152 }
00153
00154 G4double BDSDecStepper::DistChord() const
00155 {
00156 return itsDist;
00157 }
00158
00159 BDSDecStepper::~BDSDecStepper()
00160 {}