26 if (
traj.size() != 0) {
27 for (vector<HepLorentzVector*>::iterator iVector =
traj.begin();
28 iVector !=
traj.end();
42 double stepLength = 1.0;
44 HepLorentzVector step;
45 step.setX(stepLength*
sin(theta)*cos(phi));
46 step.setY(stepLength*
sin(theta)*
sin(phi));
47 step.setZ(stepLength*cos(theta));
48 step.setT(stepLength);
49 HepLorentzVector* thisVector;
50 thisVector =
new HepLorentzVector(startingPoint(1), startingPoint(2),
51 startingPoint(3), t0);
52 traj.push_back(thisVector);
53 HepLorentzVector lastVector(*thisVector);
54 for (
int i = 0; i < 600; i++) {
55 thisVector =
new HepLorentzVector();
56 *thisVector = lastVector + step;
57 traj.push_back(thisVector);
58 lastVector = *thisVector;
64 HepVector startingPoint(3);
65 startingPoint(1) = startingVector(1);
66 startingPoint(2) = startingVector(2);
67 startingPoint(3) = 0.0;
68 double theta = startingVector(3);
69 double phi = startingVector(4);
70 swim(startingPoint, theta, phi);
75 int size = params.size();
76 HepVector hepParams(size);
77 for (
int i = 0; i <
size; i++) {
78 hepParams(i + 1) = params[i];
87 HepLorentzVector *point;
88 for (vector<const DMCTrackHit*>::iterator imchit = mctrackhits.begin();
89 imchit != mctrackhits.end();
93 double r = mchit->
r, phi = mchit->
phi, z = mchit->
z;
94 point =
new HepLorentzVector(r*cos(phi), r*
sin(phi), z);
95 traj.push_back(point);
102 cout <<
"### MyTrajectory print out begin ###" << endl;
103 for (vector<HepLorentzVector*>::iterator iVector =
traj.begin();
104 iVector !=
traj.end();
106 cout << (**iVector).x() <<
" " << (**iVector).y() <<
" " << (**iVector).z() <<
" " << (**iVector).t()
109 cout <<
"### MyTrajectory print out end ###" << endl;
117 Hep3Vector
delta, point3(point(1), point(2), point(3)), trajPoint;
118 trajPoint =
traj[trajIndex]->getV();
119 delta = point3 - trajPoint;
120 if (
debug_level >= 4) cout <<
"point3 = " << point3
121 <<
"traj point = " << trajPoint
122 <<
"delta = " << delta
123 <<
"delta.mag = " << delta.mag() << endl;
132 int iBefore = 0, iAfter =
traj.size() - 1, iTry;
133 double zBefore =
traj[iBefore]->z();
134 double zAfter =
traj[iAfter]->z();
136 if (z < zBefore || z > zAfter) {
139 while (iAfter - iBefore > 1) {
141 + (int)((
double)(iAfter - iBefore)*(z - zBefore)/(zAfter - zBefore)
143 if (
debug_level > 3) cout << iBefore <<
' ' << iTry
144 <<
' ' << iAfter << endl;
145 if (iBefore == iTry) iTry++;
146 if (iAfter == iTry) iTry--;
147 if (
debug_level > 3) cout << iBefore <<
' ' << iTry
148 <<
' ' << iAfter << endl;
149 zTry =
traj[iTry]->z();
150 if (
debug_level > 3) cout <<
"zTry = " << zTry << endl;
153 zAfter =
traj[iAfter]->z();
154 }
else if (z > zTry) {
156 zBefore =
traj[iBefore]->z();
159 zBefore =
traj[iBefore]->z();
161 zAfter =
traj[iAfter]->z();
163 if (
debug_level > 3) cout << z <<
' ' << zBefore <<
' ' << zTry
164 <<
' ' << zAfter << endl;
166 double frac, otherfrac;
167 frac = (z - zBefore)/(zAfter - zBefore);
168 otherfrac = 1.0 - frac;
169 if (
debug_level > 3) cout << frac <<
' ' << otherfrac << endl;
171 <<
' '<<
traj[iAfter]->x() << endl;
173 <<
' '<<
traj[iAfter]->y() << endl;
174 x = frac*
traj[iAfter]->x() + otherfrac*
traj[iBefore]->x();
175 y = frac*
traj[iAfter]->y() + otherfrac*
traj[iBefore]->y();
176 t = frac*
traj[iAfter]->t() + otherfrac*
traj[iBefore]->t();
177 if (
debug_level > 3) cout <<
"MyTrajectory::getXYT, x, y, t = " << x <<
" " << y
183 double &xMinFrac,
double &yMin)
const {
185 a = 0.5*(yPlus - 2.0*yZero + yMinus);
186 b = 0.5*(yPlus - yMinus);
188 yMin = -b*b/(4.0*a) + c;
189 if (yMin < 0.0) yMin = 0.0;
190 xMinFrac = -b/(2.0*a);
195 if (
traj.size() != 0) {
203 HepVector trajPoint(3);
204 for (
unsigned int i = 0; i <
traj.size(); i++) {
205 trajPoint = *(
traj[i]);
206 *trajFile << tag <<
" " << i + 1 <<
" " << trajPoint(1) <<
" "
207 << trajPoint(2) <<
" " << trajPoint(3) << endl;
int getXYT(double z, double &x, double &y, double &ct) const
double doca(HepVector point)
void swimMC(vector< const DMCTrackHit * > &mctrackhits)
vector< HepLorentzVector * > traj
DetectorSystem_t system
particle type
double dist(HepVector &point, int trajIndex) const
MyTrajectory(int level=1)
float z
coordinates of hit in cm and rad
void para_min(double yMinus, double yZero, double yPlus, double &xMinFrac, double &yMin) const
virtual unsigned int getNumberOfParams()
void dump_ascii(ostream *trajFile, int tag)
virtual void swim(HepVector startingPoint, double theta, double phi)
int primary
primary track=1 not primary track=0
vector< HepLorentzVector * > * getTrajectory()