#include <TrkHelixUtils.h>
Public Member Functions | |
TrkHelixUtils () | |
Static Public Member Functions | |
static TrkExchangePar | helixFromMom (const HepPoint3D &vertex, const Hep3Vector &p, double sign, const BField &) |
static TrkExchangePar | helixFromMomErr (const BesPointErr &vertex, const BesVectorErr &p, const HepMatrix &cxp, double sign, const BField &) |
static NeutParams | lineFromMomErr (const BesPointErr &vertex, const BesVectorErr &p, const HepMatrix &cxp, double sign, const BField &) |
static HepMatrix | jacobianExtrapolate (const TrkExchangePar &, double fltNew) |
static HepSymMatrix | extrapolateCov (TrkExchangePar &, double fltNew) |
static double | fltToRad (const TrkExchangePar &hel, double rad) |
Private Member Functions | |
TrkHelixUtils & | operator= (const TrkHelixUtils &) |
TrkHelixUtils (const TrkHelixUtils &) |
Definition at line 37 of file TrkHelixUtils.h.
TrkHelixUtils::TrkHelixUtils | ( | ) | [inline] |
TrkHelixUtils::TrkHelixUtils | ( | const TrkHelixUtils & | ) | [private] |
HepSymMatrix TrkHelixUtils::extrapolateCov | ( | TrkExchangePar & | , | |
double | fltNew | |||
) | [static] |
Definition at line 98 of file TrkHelixUtils.cxx.
References TrkExchangePar::covariance(), and jacobianExtrapolate().
00099 { 00100 //---------------------------------------------------------------------- 00101 00102 return par.covariance().similarity(jacobianExtrapolate(par, fltNew)); 00103 }
double TrkHelixUtils::fltToRad | ( | const TrkExchangePar & | hel, | |
double | rad | |||
) | [static] |
Definition at line 333 of file TrkHelixUtils.cxx.
References TrkExchangePar::d0(), TrkExchangePar::omega(), Constants::pi, and TrkExchangePar::tanDip().
Referenced by MdcSegInfoSterO::calcStereo().
00333 { 00334 //------------------------------------------------------------------------ 00335 double d0 = hel.d0(); 00336 double omega = hel.omega(); 00337 double tanDip = hel.tanDip(); 00338 //std::cout<< "omega "<<omega << std::endl; 00339 // GS To be able to use this with straight lines: 00340 // if( fabs(omega) < 0.0001 ) omega=copysign(0.0001,omega); // assume the pt= 1000 GeV 00341 if( fabs(omega) < 0.0001 ) omega = (omega<0.0) ? -0.0001 : 0.0001 ; 00342 00343 double stuff = ( rad*rad - d0*d0 ) / (1 + omega * d0); 00344 if (stuff <= 0.0) return 0.; 00345 //std::cout<< "in TrkHelixUtils::fltToRad "<< stuff << std::endl; 00346 if (omega == 0.) return sqrt(stuff); 00347 double sinAngle = 0.5 * omega * sqrt(stuff); 00348 double dist2d = 0; 00349 if (sinAngle < -1.0 || sinAngle > 1.0) { 00350 dist2d = Constants::pi / fabs(omega); 00351 } else { 00352 dist2d = 2. * asin(sinAngle) / omega; 00353 } 00354 //std::cout<< "in TrkHelixUtils::fltToRad dist2d "<< dist2d << std::endl; 00355 return dist2d * sqrt(1. + tanDip*tanDip); 00356 }
TrkExchangePar TrkHelixUtils::helixFromMom | ( | const HepPoint3D & | vertex, | |
const Hep3Vector & | p, | |||
double | sign, | |||
const BField & | ||||
) | [static] |
Definition at line 106 of file TrkHelixUtils.cxx.
References BField::bFieldNominal(), BField::cmTeslaToGeVc, cos(), and sin().
00107 { 00108 //---------------------------------------------------------------------- 00109 //As documented in 00110 //R.~Luchsinger and C.~Grab, Comp.~Phys.~Comm.~{\bf 76}, 263-280 (1993). 00111 //Equation 14 on page 270. 00112 //Note: We use the opposite sign for d0. 00113 //We use tandip and not the dip angle itself. 00114 00115 register double phip,gamma,rho,pt; 00116 static HepVector pars(5); 00117 00118 double px = pmom.x(); 00119 double py = pmom.y(); 00120 pt=sqrt(px*px+py*py); 00121 00122 if (pt < 1.e-7) pt = 1.e-7; // hack to avoid pt=0 tracks 00123 if (fabs(px) < 1.e-7) px = (px<0.0) ? -1.e-7 : 1.e-7; // hack to avoid pt=0 tracks 00124 00125 double Bval = fieldmap.bFieldNominal(); 00126 00127 pars(3)=-BField::cmTeslaToGeVc*Bval*sign/pt; //omega 00128 pars(5)=pmom.z()/pt; //tandip 00129 00130 rho=1./pars(3); 00131 phip=atan2(py,px); 00132 double cosphip=cos(phip); double sinphip=sin(phip); 00133 00134 gamma=atan((pos.x()*cosphip+pos.y()*sinphip)/ 00135 (pos.x()*sinphip-pos.y()*cosphip-rho)); // NOTE: do NOT use atan2 here! 00136 00137 00138 pars(2)=BesAngle(phip+gamma).rad(); //phi0 00139 pars(1)=(rho+pos.y()*cosphip-pos.x()*sinphip)/cos(gamma)-rho; //d0 00140 pars(4)=pos.z()+rho*gamma*pars(5); //z0 00141 00142 return TrkExchangePar(pars); 00143 }
TrkExchangePar TrkHelixUtils::helixFromMomErr | ( | const BesPointErr & | vertex, | |
const BesVectorErr & | p, | |||
const HepMatrix & | cxp, | |||
double | sign, | |||
const BField & | ||||
) | [static] |
Definition at line 146 of file TrkHelixUtils.cxx.
References BField::bFieldNominal(), BField::cmTeslaToGeVc, cos(), BesVectorErr::covMatrix(), genRecEmupikp::i, ganga-rec::j, DifArray::jacobian(), and boss::pos.
Referenced by TrkCompTrk::TrkCompTrk().
00148 { 00149 //---------------------------------------------------------------------- 00150 00151 DifNumber xDF(pos.x(),1,6), yDF(pos.y(),2,6), zDF(pos.z(),3,6); 00152 DifNumber pxDF(pmom.x(), 4, 6); 00153 DifNumber pyDF(pmom.y(), 5, 6); 00154 DifNumber pzDF(pmom.z(), 6, 6); 00155 00156 static DifNumber phip, cosphip, sinphip, gamma; 00157 static DifArray pars(5,6); 00158 00159 DifNumber invpt = pxDF; 00160 invpt *= pxDF; 00161 invpt += pyDF*pyDF; 00162 invpt = sqrt(invpt); 00163 00164 if (invpt < 1.e-7) invpt = 1.e-7; // hack to avoid pt=0 tracks 00165 if (fabs(pxDF) < 1.e-7) pxDF = pxDF<0?-1e-7:1e-7; // hack to avoid pt=0 tracks 00166 invpt = 1./invpt; 00167 00168 //omega 00169 double Bval = fieldmap.bFieldNominal(); 00170 // if (fabs(Bval) < 1.e-7) { // hack to avoid B=0 tracks (very far Z) 00171 // Bval = 1.e-7; 00172 // } 00173 00174 // pars(3) = -BField::cmTeslaToGeVc*Bval*sign/pt; 00175 pars(3) = invpt; 00176 pars(3) *= sign; 00177 pars(3) *= Bval; 00178 pars(3) *= -BField::cmTeslaToGeVc; 00179 00180 // pars(5) = pzDF / pt; //tandip 00181 pars(5) = pzDF; 00182 pars(5) *= invpt; 00183 00184 DifNumber rho = 1./pars[2]; 00185 phip=atan2(pyDF,pxDF); 00186 phip.cosAndSin(cosphip,sinphip); 00187 00188 // pars(1) = (rho + yDF*cosphip - xDF*sinphip)/cos(gamma) - rho;//d0 00189 pars(1) = rho; 00190 pars(1) += yDF*cosphip; 00191 pars(1) -= xDF*sinphip; // continued below ... 00192 00193 gamma=atan((xDF*cosphip+yDF*sinphip)/ -pars(1)); // NOTE: do NOT use atan2 here 00194 00195 pars(1) /= cos(gamma); 00196 pars(1) -= rho; 00197 00198 // pars(2) = phip+gamma; //phi0 00199 pars(2) = phip; 00200 pars(2) += gamma; 00201 // pars(2).mod(0., Constants::twoPi); 00202 00203 00204 // pars(4) = zDF + rho*gamma*pars[4]; //z0 00205 pars(4) = pars[4]; // weird 00206 pars(4) *= gamma; 00207 pars(4) *= rho; 00208 pars(4) += zDF; 00209 00210 00211 // Get error matrix for position and momentum 00212 00213 static HepSymMatrix posandmomErr(6); 00214 static HepVector parsVec(5); 00215 00216 int i; 00217 for (i = 1; i <= 3; ++i) { 00218 int j; 00219 for (j = 1; j <= i; ++j) { 00220 // with "fast" access, row should be >= col 00221 posandmomErr.fast(i,j) = pos.covMatrix().fast(i,j); 00222 posandmomErr.fast(i+3,j+3) = pmom.covMatrix().fast(i,j); 00223 } 00224 for (j = 1; j <= 3; ++j) { 00225 posandmomErr.fast(j+3,i) = cxp(i,j); 00226 } 00227 } 00228 for (i = 1; i <= 5; ++i) { 00229 // make the array of DifNums into a HepVector 00230 // (needed for TrkExchangePar init) 00231 parsVec(i) = pars(i).number(); 00232 } 00233 // Now calculate error on the helix pars--the real calculation 00234 00235 return TrkExchangePar(parsVec,posandmomErr.similarity(pars.jacobian()) ); 00236 }
HepMatrix TrkHelixUtils::jacobianExtrapolate | ( | const TrkExchangePar & | , | |
double | fltNew | |||
) | [static] |
Definition at line 33 of file TrkHelixUtils.cxx.
References cos(), TrkExchangePar::omega(), sin(), and TrkExchangePar::tanDip().
Referenced by extrapolateCov().
00034 { 00035 //------------------------------------------------------------------------ 00036 00037 //---------------------------------------------------------------------- 00038 // Compute and return the jacobian that takes the covariance matrix 00039 // from fltOld to fltNew 00040 // 00041 // "fltLen" -- the signed 3-d pathlength a particle travels 00042 // along the orbit starting from the point on the orbit 00043 // close to the origin. 00044 // 00045 // Each element in this matrix is a partial derivative of an orbit 00046 // parameter at some reference point to an orbit parameter at 00047 // the fit point. What is kept fixed in taking these partial 00048 // derivatives is that the orbit parameters are those at the point 00049 // on the orbit that is closest in the x-y plane to the reference point. 00050 // This transform matrix has the property that transform(-arclength) 00051 // is the inverse of transform(arclength). 00052 // (repaired by Gerry Lynch, I think -- sfs) 00053 //---------------------------------------------------------------------- 00054 00055 HepMatrix transform(5, 5, 0); 00056 00057 double tanDip = par.tanDip(); 00058 double omega = par.omega(); 00059 // Convert to 2-d arclength 00060 double darc = fltNew / sqrt(1. + tanDip * tanDip); 00061 double dphi = omega * darc; 00062 double cosDphi = cos(dphi); 00063 double sinDphi = sin(dphi); 00064 00065 // partials of d0 00066 00067 transform[0][0] = cosDphi; 00068 transform[0][1] = sinDphi / omega; 00069 transform[0][2] = (1.0-cosDphi) / (omega * omega); 00070 00071 // partials of phi0 00072 00073 transform[1][0] = -sinDphi * omega; 00074 transform[1][1] = cosDphi; 00075 transform[1][2] = sinDphi / omega; 00076 00077 // partials of omega 00078 00079 transform[2][2] = 1.0; 00080 00081 // partials of z0 00082 00083 transform[3][0] = -tanDip * sinDphi; 00084 transform[3][1] = -tanDip * (1.0-cosDphi) / omega; 00085 transform[3][2] = -tanDip * (dphi-sinDphi) / (omega*omega); 00086 transform[3][3] = 1.; 00087 transform[3][4] = darc; 00088 00089 // partials of tanDip 00090 00091 transform[4][4] = 1.; 00092 00093 return transform; 00094 }
NeutParams TrkHelixUtils::lineFromMomErr | ( | const BesPointErr & | vertex, | |
const BesVectorErr & | p, | |||
const HepMatrix & | cxp, | |||
double | sign, | |||
const BField & | ||||
) | [static] |
Definition at line 238 of file TrkHelixUtils.cxx.
References BesVectorErr::covMatrix(), genRecEmupikp::i, ganga-rec::j, and boss::pos.
Referenced by TrkCompTrk::TrkCompTrk().
00240 { 00241 //---------------------------------------------------------------------- 00242 00243 DifNumber xDF(pos.x(),1,6), yDF(pos.y(),2,6), zDF(pos.z(),3,6); 00244 DifNumber pxDF(pmom.x(), 4, 6); 00245 DifNumber pyDF(pmom.y(), 5, 6); 00246 DifNumber pzDF(pmom.z(), 6, 6); 00247 00248 static DifArray pars(6,6); 00249 00250 DifNumber pt = pxDF; 00251 pt *= pxDF; 00252 pt += pyDF*pyDF; 00253 00254 if (pt < 1.e-14) pt = 1.e-14; // hack to avoid pt=0 tracks 00255 if (fabs(pxDF) < 1.e-7) pxDF = pxDF<0?-1e-7:1e-7; // hack to avoid pt=0 tracks 00256 00257 // pars(3) = sqrt(pt*pt+pzDF*pzDF); // Magnitude of p 00258 pars(3) = pzDF; 00259 pars(3) *= pzDF; 00260 pars(3) += pt; 00261 pars(3) = sqrt(pars(3)); 00262 00263 pt = sqrt(pt); 00264 DifNumber invpt = 1./pt; 00265 00266 //Next lines modified by Eugenio Paoloni 19-XI-98 00267 00268 // DifNumber pvxDF=pxDF/pt; // versor along pt x 00269 DifNumber pvxDF=pxDF; 00270 pvxDF *= invpt; 00271 00272 // DifNumber pvyDF=pyDF/pt; // and y component 00273 DifNumber pvyDF=pyDF; 00274 pvyDF *= invpt; 00275 00276 00277 // pars(5) = pzDF / pt; //tandip 00278 pars(5) = pzDF; 00279 pars(5) *= invpt; 00280 00281 pars(2) = atan2(pyDF,pxDF); //phi0 00282 00283 // pars(1) = yDF*pvxDF - xDF*pvyDF;//d0 00284 pars(1) = yDF; 00285 pars(1) *= pvxDF; 00286 pars(1) -= xDF*pvyDF; 00287 00288 // pars(4) = zDF -pars(5)*(xDF*pvxDF+yDF*pvyDF) ; //z0 00289 pars(4) = xDF; 00290 pars(4) *= pvxDF; 00291 pars(4) += yDF*pvyDF; 00292 00293 // insert this in the middle for speed 00294 // pars(6) = (xDF*pvxDF+yDF*pvyDF)*pars(3)/pt;//s0 00295 pars(6) = pars(4); 00296 pars(6) *= pars(3); 00297 pars(6) *= invpt; 00298 00299 00300 pars(4) *= -pars(5); 00301 pars(4) += zDF; 00302 00303 00304 00305 00306 // Get error matrix for position and momentum 00307 00308 static HepSymMatrix posandmomErr(6); 00309 static HepVector parsVec(6); 00310 00311 int i; 00312 for (i = 1; i <= 3; ++i) { 00313 int j; 00314 for (j = 1; j <= i; ++j) { 00315 // with "fast" access, row should be >= col 00316 posandmomErr.fast(i,j) = pos.covMatrix().fast(i,j); 00317 posandmomErr.fast(i+3,j+3) = pmom.covMatrix().fast(i,j); 00318 } 00319 for (j = 1; j <= 3; ++j) { 00320 posandmomErr.fast(j+3,i) = cxp(i,j); 00321 } 00322 } 00323 for (i = 1; i <= 6; ++i) { 00324 // make the array of DifNums into a HepVector 00325 // (needed for TrkExchangePar init) 00326 parsVec(i) = pars(i).number(); 00327 } 00328 // Now calculate error on the helix pars--the real calculation 00329 return NeutParams(parsVec,posandmomErr.similarity(pars.jacobian()) ); 00330 }
TrkHelixUtils& TrkHelixUtils::operator= | ( | const TrkHelixUtils & | ) | [private] |