2 #define M_CHARGED_PI 0.13957018
5 #define TRACKING_RADIUS_MAX 60.0
6 #define TIME_STEP_SIZE 0.03
20 void grkuta_(
float* charge,
float* step,
float* vect,
float* vout);
28 double xp0, z0, theta0, phi0, ptinv;
49 double sinTheta =
sin(theta);
50 double cosTheta = cos(theta);
51 double sinPhi =
sin(phi);
52 double cosPhi = cos(phi);
53 double xStart = xp0*sinPhi;
54 double yStart = -xp0*cosPhi;
57 HepLorentzVector* thisVectorPtr;
58 float charge, step, vect[7], vout[7];
59 if (ptinv > 0.0) {charge = 1.0;}
else {charge = -1.0;}
60 vect[0] = xStart; vect[1] = yStart; vect[2] = zStart;
61 vect[3] = sinTheta*cosPhi; vect[4] = sinTheta*sinPhi; vect[5] = cosTheta;
62 double ptot = abs(1.0/(ptinv*sinTheta));
66 double beta =
sqrt(1.0 - 1.0/(gamma*gamma));
68 double ctStep =
c*tStep;
69 step = (float)(beta*ctStep);
71 globalGrkutaPtr =
this;
72 thisVectorPtr =
new HepLorentzVector;
73 thisVectorPtr->setX(vect[0]);
74 thisVectorPtr->setY(vect[1]);
75 thisVectorPtr->setZ(vect[2]);
76 thisVectorPtr->setT(ct);
77 traj.push_back(thisVectorPtr);
78 for (
int i = 0; i < 2000; i++) {
79 grkuta_(&charge, &step, vect, vout);
81 thisVectorPtr =
new HepLorentzVector;
82 thisVectorPtr->setX(vout[0]);
83 thisVectorPtr->setY(vout[1]);
84 thisVectorPtr->setZ(vout[2]);
85 thisVectorPtr->setT(ct);
86 if (
debug_level > 3) cout << setprecision(14) <<
"MyTrajectoryGrcuta::swim: i = " << i <<
" LorentzVector = " << *thisVectorPtr << endl;
87 traj.push_back(thisVectorPtr);
89 for (
int j = 0; j < 7; j++) {vect[j] = vout[j];}
95 const double kGperT = 10.0;
100 HepVector B = globalGrkutaPtr->
getField(r);
void gufld_(float *x, float *f)
HepVector getField(HepVector &r)
MyTrajectoryGrkuta * globalGrkutaPtr
#define TRACKING_RADIUS_MAX
vector< HepLorentzVector * > traj
MyTrajectoryGrkuta(const DMagneticFieldMap *bfield, int level=1)
void swim(const HepVector ¶m)
int grkuta_(double *CHARGE, double *STEP, double *VECT, double *VOUT, const DMagneticFieldMap *bfield)