MarlinTPC  1.2.0
TrackRecoTools.h
1 #ifndef TrackRecoTools_h
2 #define TrackRecoTools_h 1
3 
4 // C++
5 #include <string>
6 #include <cmath>
7 #include <algorithm>
8 #include <CLHEP/Vector/ThreeVector.h>
9 
10 // Marlin
11 #include "marlin/Processor.h"
12 
13 // LCIO
14 #include "lcio.h"
15 #include <EVENT/TrackerHit.h>
16 #include <EVENT/Track.h>
17 
18 using namespace marlintpc;
19 
20 namespace tracktools{
21 
22 static bool compareHit(const std::pair<TrackerHit*, double> & pair1,
23  const std::pair<TrackerHit*, double> & pair2)
24  {
25  return pair1.second < pair2.second;
26  }
27 
28 static double CalculateS(EVENT::Track* currentTrack, EVENT::TrackerHit * currentHit)
29 {
30 
31  double pca_x = ( -1.0 * currentTrack->getD0() ) * std::sin(currentTrack->getPhi());
32  double pca_y = ( currentTrack->getD0() ) * std::cos(currentTrack->getPhi());
33 
34  double s = 0.;
35  const double* HitPosition=currentHit->getPosition();
36 
37  if(fabs(currentTrack->getOmega()) < 0.000000001)
38  {
39  s = std::sqrt(
40  static_cast<long double>(( HitPosition[0] - pca_x ) * ( HitPosition[0] - pca_x )
41  + ( HitPosition[1] - pca_y ) * ( HitPosition[1] - pca_y )));
42 
43  //make sure the sign of s is right
44  if(( ( ( -1.25 * M_PI ) <= currentTrack->getPhi() ) && ( currentTrack->getPhi() <= ( -0.75 * M_PI ) ) )
45  || ( ( ( -0.25 * M_PI ) <= currentTrack->getPhi() ) && ( currentTrack->getPhi() <= ( 0.25 * M_PI ) ) )
46  || ( ( ( 0.75 * M_PI ) <= currentTrack->getPhi() ) && currentTrack->getPhi() <= ( 1.25 * M_PI ) ))
47  {
48  if(pca_x < HitPosition[0])
49  {
50  s = ( -1.0 ) * s;
51  }
52  }
53  else
54  {
55  if(pca_y < HitPosition[1])
56  {
57  s = ( -1.0 ) * s;
58  }
59  }
60  s = (-1.)*s; //Needed to make tracks look correctly in CEDViewer
61  }
62  else
63  {
64  double r = 1. / currentTrack->getOmega();
65  //center of circle
66  double xm, ym;
67  xm = (r - currentTrack->getD0()) * sin(currentTrack->getPhi());
68  ym = -(r - currentTrack->getD0()) * cos(currentTrack->getPhi());
69  s = r * ( atan2(pca_y - ym , pca_x - xm)
70  -
71  atan2(HitPosition[1]- ym , HitPosition[0] - xm)
72  );
73 
74  double s_1, s_2;
75  s_1 = s;
76  if(s < 0)
77  {
78  s_2 = s + fabs(2.*M_PI/currentTrack->getOmega());
79  }
80  else
81  {
82  s_2 = s - fabs(2.*M_PI/currentTrack->getOmega());
83  }
84 
85  if(fabs(s_1) <= fabs(s_2)) s = s_1;
86  else s = s_2;
87 
88  }
89  return s;
90 }
91 
96 //static double calculateMaxTrackLength(EVENT::Track* track)
97 //{
98 // typedef std::vector< std::pair<EVENT::TrackerHit*, double> > HitDouble;
99 // HitDouble HitPairVec;
100 // // int hitNr=0;
101 //
102 // // first calculate some track parameters from track:
103 // double phi = track->getPhi(); // just to improve readability
104 //
105 // // unit vector from reference point to point of closest approach (pca) (xy plane only)
106 // CLHEP::Hep3Vector n_pca( -sin( phi ), cos ( phi ) );
107 //
108 // // reference point in xy projection
109 // CLHEP::Hep3Vector reference_point( track->getReferencePoint ()[0], track->getReferencePoint ()[1],
110 // track->getReferencePoint ()[2]);
111 //
112 // // point of closest approach (3D point which is closest in the xy plane)
113 // CLHEP::Hep3Vector pca = reference_point
114 // + track->getD0() * n_pca // the xy part
115 // + track->getZ0() * CLHEP::Hep3Vector(0,0,1); // the z part
116 //
117 // streamlog_out(DEBUG4) << "PCA is " << pca[0] << " , " << pca[1] << " , " << pca[2] << std::endl;
118 //
119 // // psi is the angle of the pca wrt the helix center
120 // // (Note: in the paper from Lasiuk et al. from 1998 this is Phi_0,
121 // // while Psi is the angle in track direktion.
122 // // As Phi is the angle in track direction in LCIO,
123 // // Phi and psi are exchanged wrt. the Lasiuk paper.)
124 // double psi = track->getPhi() + ( track->getOmega()<0 ? -1 : 1 ) * M_PI_2;
125 //
126 // // wrap psi to -2pi .. 2pi to prevent the while loop running for ages
127 // // in case there are nonsense values for the track's phi
128 // if ( std::fabs(psi) > 20 )
129 // {
130 // streamlog_out(WARNING3) << "Psi is " << psi << ", check phi of the input tracks!"
131 // << std::endl;
132 // }
133 //
134 // psi = std::fmod(psi, 2*M_PI);
135 //
136 // TrackerHitVec hitVec = track->getTrackerHits();
137 // float sumZ=0;
138 // for ( TrackerHitVec::iterator hitIter = hitVec.begin(); hitIter < hitVec.end(); hitIter++ )
139 // {
140 // TrackerHit* trackerHit = *hitIter;
141 // // A Hep3Vector for the test hit
142 // CLHEP::Hep3Vector hit(trackerHit->getPosition()[0],trackerHit->getPosition()[1],
143 // trackerHit->getPosition()[2] );
144 //
145 // sumZ+=hit[2];
146 // double s; // the track length in the xy-plane (projection)
147 // if ( track->getOmega() != 0. )
148 // {
149 // double co = std::fabs(track->getOmega())*(hit[0] - pca[0]) + cos(psi);
150 // double si = std::fabs(track->getOmega())*(hit[1] - pca[1]) + sin(psi);
151 // streamlog_out(DEBUG2)<< "TrackFitterBase: co "<<co <<" : si " << si
152 // << " : psi " << psi << std::endl;
153 //
154 // double phi = ( co == 0 ? ( si < 0 ? - M_PI_2 : M_PI_2 ) : ( co > 0 ? - atan(si/co) :
155 // M_PI - atan(si/co) ) );
156 //
157 // // wrap phi to be in +- pi wrt. -psi (minus because the pca is at -psi wrt helix centre)
158 // while ( -psi - phi > M_PI ) phi+=2*M_PI ;
159 // while ( -psi - phi < -M_PI ) phi-=2*M_PI ;
160 //
161 // s = (-psi - phi) / track->getOmega();
162 // // s = ( atan(si/co) + ( co < 0 ? M_PI : 0 ) - psi ) / track->getOmega();
163 //
164 // // //This is the method from Lasiuk et al. It does not work because
165 // // // the atan only works for +- pi/2, but we need a full circle
166 // // s = atan( ( (hit[1] - pca[1] ) * sin(phi) + (hit[0] -pca[0]) * cos(phi) ) /
167 // // ( 1/std::fabs( track->getOmega() ) +
168 // // (hit[0] - pca[0]) * sin(phi) -
169 // // (hit[1] - pca[1]) * cos(phi)
170 // // )
171 // // )
172 // // / track->getOmega() ;
173 // streamlog_out(DEBUG3)<< "TrackFitterBase: Hit "<<hit[0] <<" : " <<hit[1] <<" : "<<hit[2]
174 // <<" ; s is " << s << std::endl;
175 // }
176 // else
177 // {
178 // s = (hit[1] - pca[1] ) * sin(phi) + (hit[0] -pca[0]) * cos(phi); // tracklength in xy projection
179 //
180 // streamlog_out(DEBUG3)<< "TrackFitterBase: Hit "<<hit[0] <<" : " <<hit[1] <<" : "<<hit[2]
181 // <<" ; s is " << s << std::endl;
182 // }
183 //
184 // s = s * std::sqrt( 1 + track->getTanLambda()*track->getTanLambda()); // here comes the third dimension
185 // std::pair<TrackerHit*, double> HitPair( trackerHit , s);
186 // HitPairVec.push_back( HitPair );
187 // }// for hitIter
188 //
189 // sort( HitPairVec.begin(), HitPairVec.end(), compareHit );
190 // sort( HitPairVec.begin(), HitPairVec.end(), compareHit );
191 //
192 // //effective length of track is distance between outer hits times n/(n-1):
193 // //distance between outer points on track
194 // double l = fabs( (*HitPairVec.begin()).second - (*(HitPairVec.end()-1)).second );
195 // //effective length
196 // double l_eff = l*hitVec.size()/(hitVec.size()-1);
197 // return l_eff;
198 //}
199 }
200 
201 #endif
Generated on Fri Dec 2 2016 12:42:18 for MarlinTPC by   doxygen 1.8.6