3 #include <CLHEP/Matrix/Matrix.h>
4 #include <CLHEP/Matrix/Vector.h>
11 #define TRACKING_RADIUS_MAX 60.0
15 :
MyTrajectory(level), nparams(5), params(nparams), BConst(B_in), delta(nparams, 0.0001), bfield(NULL),
22 :
MyTrajectory(level), nparams(5), params(nparams), delta(nparams, 0.0001), bfield(bfield_in),
32 double xp0, z0, theta0, phi0, ptinv;
43 double phi,
double ptinv) {
44 if (
debug_level >= 3) cout <<
"starting B-field swim, xp0 = " << xp0
45 <<
" z0 = " << z0 <<
" theta = " << theta
46 <<
" phi = " << phi <<
" ptinv = " << ptinv << endl;
54 const double c = 29.9792458;
55 double m=0.139, dt=0.01;
57 HepVector r_0(3), r_1(3), r_2(3), k(3), v_0(3);
58 HepMatrix A(3,3,0),
I(3,3,1), M(3,3), onePlusA(3,3);
62 double p = 1.0/ptinv/
sin(theta);
63 double E =
sqrt(m*m + p*p);
65 beta =
sqrt(1.0 - 1.0/(gamma*gamma));
67 if (ptinv > 0.0) {e = 1.0;}
else {e = -1.0;}
69 k_1 = (e*1.60217653e-19)*(dt*1.0e-9)/(gamma*(m*1.78266181e-27));
77 double sinTheta =
sin(theta);
78 double cosTheta = cos(theta);
79 double sinPhi =
sin(phi);
80 double cosPhi = cos(phi);
81 v_0(1) = v*sinTheta*cosPhi;
82 v_0(2) = v*sinTheta*sinPhi;
85 HepLorentzVector *thisVector;
89 thisVector =
new HepLorentzVector(r_0(1), r_0(2), r_0(3), 0.0);
90 traj.push_back(thisVector);
92 thisVector =
new HepLorentzVector(r_1(1), r_1(2), r_1(3), cdt);
93 traj.push_back(thisVector);
94 double ctime = 2.0*cdt;
96 for (
int istep = 0; istep < 2000; istep++) {
98 if (
debug_level > 2) cout <<
"MyTrajectoryBfield::swim B vector:" << B;
108 M = onePlusA.inverse(ierr)*(
I - A);
109 r_2 = r_1 + M*(r_1 - r_0);
111 if (istep%300 == 0) cout << istep <<
' ' << r_2(1) <<
' '
112 << r_2(2) <<
' ' << r_2(3) << endl;
114 thisVector =
new HepLorentzVector(r_2(1), r_2(2), r_2(3), ctime);
115 traj.push_back(thisVector);
116 r_0 = r_1; r_1 = r_2;
MyTrajectoryBfield(const DMagneticFieldMap *bfield, int level=1)
HepVector getField(HepVector &r)
void swim(const HepVector ¶m)
unsigned int getNumberOfParams()
vector< HepLorentzVector * > traj
const DMagneticFieldMap * bfield
virtual void GetField(const DVector3 &pos, DVector3 &Bout) const =0
#define TRACKING_RADIUS_MAX