1 #ifndef _MYTRAJECTORY_H_
2 #define _MYTRAJECTORY_H_
6 #include <CLHEP/Vector/LorentzVector.h>
10 #define DIST_BIG 10000.0
11 #define MAX_ITERATIONS 100000
20 virtual void swim(HepVector startingPoint,
double theta,
double phi);
21 virtual void swim(
const HepVector &startingVector);
22 void swim(
const vector<double> &startingStdVector);
23 void swimMC(vector<const DMCTrackHit*> &mctrackhits);
25 vector<HepLorentzVector*>* getTrajectory();
26 template<
class C>
double doca(C& spaceObject, HepLorentzVector &poca)
const;
28 virtual unsigned int getNumberOfParams();
29 double dist(HepVector& point,
int trajIndex)
const;
30 double dist(
DLine& line,
int trajIndex)
const;
31 void para_min(
double yMinus,
double yZero,
double yPlus,
double &xMinFrac,
36 int getXYT(
double z,
double &
x,
double &
y,
double &ct)
const;
37 void dump_ascii(ostream *trajFile,
int tag);
40 vector<HepLorentzVector*>
traj;
49 unsigned int ilo, ihi, imid;
51 ihi = traj.size() - 1;
52 double distlo, disthi, distmid, distdown, distup;
53 distlo = dist(spaceObject, ilo);
54 disthi = dist(spaceObject, ihi);
56 cout <<
"MyTrajectory::doca: end point of trajectory too far away " << distlo <<
" " << disthi << endl;
60 if (distlo < disthi) {
65 distmid = dist(spaceObject, imid);
66 if (debug_level >= 4) cout << setprecision(14) <<
"MyTrajectory::doca: initialize " << ilo <<
" " << imid <<
" " << ihi <<
" " << distlo <<
" " << distmid <<
" " << disthi << endl;
68 cout <<
"MyTrajectory::doca: distmid is not a number: " << ilo <<
" " << imid <<
" " << ihi <<
" " << distlo <<
" " << distmid <<
" " << disthi << endl;
72 if (distmid >= distlo || distmid >= disthi) {
73 cout <<
"MyTrajectory::doca: bad initialization of doca search: " << ilo <<
" " << imid <<
" " << ihi <<
" " << distlo <<
" " << distmid <<
" " << disthi << endl;
77 double x1, x2, y1, y2, xnew, distnew = 0;
80 while (imid - ilo > 1 || ihi - imid > 1) {
81 x1 = (double)(imid - ilo);
82 x2 = -(double)(ihi - imid);
83 y1 = distmid - distlo;
84 y2 = distmid - disthi;
91 inew = (
unsigned int)(xnew + 0.5);
92 if (inew < ilo || inew > ihi) {
93 cout <<
"MyTrajectory::doca: jumped out of the miminum bracket: " << ilo <<
" " << inew <<
" " << ihi << endl;
97 distnew = dist(spaceObject, inew);
98 if (debug_level >= 4) cout <<
"MyTrajectory::doca: xnew = " << xnew <<
" inew = " << inew <<
" distnew = " << distnew << endl;
101 distup = dist(spaceObject, imid + 1);
102 distdown = dist(spaceObject, imid - 1);
103 if (distup > distmid && distdown > distmid) {
108 }
else if (distup > distmid) {
109 if (ihi - imid > 1) {
111 disthi = dist(spaceObject, ihi);
114 distmid = dist(spaceObject, imid);
117 if (imid - ilo > 1) {
119 distlo = dist(spaceObject, ilo);
122 distmid = dist(spaceObject, imid);
127 if (distnew < distmid) {
148 if (debug_level >= 4) cout <<
"MyTrajectory::doca, after housekeeping: " << ilo <<
" " << imid <<
" " << ihi <<
" " << distlo <<
" " << distmid <<
" " << disthi << endl;
151 cout <<
"MyTrajectory::doca: too many iterations " << ilo <<
" " << imid <<
" " << ihi <<
" " << distlo <<
" " << distmid <<
" " << disthi << endl;
156 double dlon, dmidn, dhin, dinterp;
157 dlon = dist(spaceObject, ilo);
158 dmidn = dist(spaceObject, imid);
159 dhin = dist(spaceObject, ihi);
160 double xMinFrac, yMin;
161 para_min(dlon*dlon, dmidn*dmidn, dhin*dhin, xMinFrac, yMin);
162 dinterp =
sqrt(yMin);
163 if (xMinFrac > 0.0) {
164 poca = xMinFrac*(*traj[ihi]) + (1.0 - xMinFrac)*(*traj[imid]);
166 poca = -xMinFrac*(*traj[ilo]) + (1.0 + xMinFrac)*(*traj[imid]);
168 if (debug_level >= 4) cout <<
"doca final calc: ilo,imid,ihi = " << ilo <<
',' << imid <<
',' << ihi <<
" dists = " << dlon <<
"," << dmidn <<
"," << dhin <<
" dinterp = " << dinterp << endl;
173 #endif // _MYTRAJECTORY_H_
double doca(C &spaceObject, HepLorentzVector &poca) const
vector< HepLorentzVector * > traj
virtual vector< double > getDelta()