10 TMatrixT<double> smoothed_state;
11 TMatrixT<double> smoothed_cov;
17 TMatrixT<double> pos_tmp(H * smoothed_state);
18 pos.ResizeTo(pos_tmp);
29 TMatrixT<double> smoothed_state;
30 TMatrixT<double> smoothed_cov;
37 TMatrixT<double> pos_tmp(H * smoothed_state);
38 pos.ResizeTo(pos_tmp);
41 if (ret) *ret =
false;
46 if (pos.GetNrows() != 2 || pos.GetNcols() != 1){
47 GFException exc(
"GFTools::getSmoothedPosXYZ ==> dimension of hit in plane is not 2, cannot calculate (x,y,z) hit position",__LINE__,__FILE__);
52 TVector3 pos3D(plane.
getO());
53 pos3D +=
pos(0,0) * plane.
getU();
54 pos3D +=
pos(1,0) * plane.
getV();
62 TMatrixT<double> smoothed_state;
63 TMatrixT<double> smoothed_cov;
78 TMatrixT<double> pos_tmp(H * smoothed_state);
79 pos.ResizeTo(pos_tmp);
83 if (ret) *ret =
false;
90 if (pos.GetNrows() != 2 || pos.GetNcols() != 1){
92 if (ret) *ret =
false;
101 TVector3 pos3D(plane.
getO());
102 pos3D +=
pos(0,0) * plane.
getU();
103 pos3D +=
pos(1,0) * plane.
getV();
105 if (ret) *ret =
true;
115 TMatrixT<double> smoothed_state, smoothed_cov, auxInfo;
119 if (!
getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, smoothing_plane, auxInfo,
false)) {
120 if (ret) *ret =
false;
124 if(rep->hasAuxInfo()) {
125 rep->setData(smoothed_state, smoothing_plane, &smoothed_cov, &auxInfo);
127 rep->setData(smoothed_state, smoothing_plane, &smoothed_cov);
130 if (ret) *ret =
true;
131 return rep->getMom(smoothing_plane);
137 TMatrixT<double> smoothed_state;
138 TMatrixT<double> smoothed_cov;
139 TMatrixT<double>
pos;
145 TMatrixT<double> pos_tmp(H * smoothed_state);
146 pos.ResizeTo(pos_tmp);
156 TMatrixT<double> smoothed_state;
157 TMatrixT<double> smoothed_cov;
163 TMatrixT<double> cov_tmp(smoothed_cov,TMatrixT<double>::kMultTranspose,H);
164 TMatrixT<double>
cov(H,TMatrixT<double>::kMult,cov_tmp);
172 TMatrixT<double> smoothed_state;
173 TMatrixT<double> smoothed_cov;
179 TMatrixT<double> cov_tmp(smoothed_cov,TMatrixT<double>::kMultTranspose,H);
180 TMatrixT<double>
cov(H,TMatrixT<double>::kMult,cov_tmp);
189 TMatrixT<double> auxInfo;
197 TMatrixT<double> auxInfo;
204 TMatrixT<double> auxInfo;
211 TMatrixT<double> auxInfo;
213 smoothing_plane, auxInfo, extrapolation_allowed);
222 std::cout <<
"Trying to get smoothed hit coordinates from a track without smoothing! Aborting..." << std::endl;
223 TMatrixT<double> failed(1,1);
228 std::cout <<
"Hit number out of bounds while trying to get smoothed hit coordinates! Aborting..." <<std::endl;
229 TMatrixT<double> failed(1,1);
237 TMatrixT<double> fUpSt;
238 TMatrixT<double> fUpCov;
239 TMatrixT<double> fAuxInfo;
240 TMatrixT<double>* fAuxInfoP;
241 TMatrixT<double> bUpSt;
242 TMatrixT<double> bUpCov;
243 TMatrixT<double> bAuxInfo;
244 TMatrixT<double>* bAuxInfoP;
246 TMatrixT<double> fSt;
247 TMatrixT<double> fCov;
248 TMatrixT<double> bSt;
249 TMatrixT<double> bCov;
268 bAuxInfoP = &bAuxInfo;
272 if(bUpSt.GetNrows() == 0) {
printf(
"Ku-0\n");
return false; };
273 rep->
setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
274 rep->
extrapolate(smoothing_plane,smoothed_state,smoothed_cov);
286 fAuxInfoP = &fAuxInfo;
290 if(fUpSt.GetNrows() == 0) {
printf(
"Ku-1\n");
return false; };
291 rep->
setData(fUpSt,fPl,&fUpCov,fAuxInfoP);
292 rep->
extrapolate(smoothing_plane,smoothed_state,smoothed_cov);
303 fAuxInfoP = &fAuxInfo;
304 bAuxInfoP = &bAuxInfo;
306 fAuxInfoP = bAuxInfoP = NULL;
312 if(fUpSt.GetNrows() == 0 || bUpSt.GetNrows() == 0) {
printf(
"Ku-2\n");
return false; };
316 rep->
setData(fUpSt,fPl,&fUpCov,fAuxInfoP);
320 rep->
setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
332 bAuxInfoP = &bAuxInfo;
336 if(bUpSt.GetNrows() == 0)
return false;
337 rep->
setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
338 rep->
extrapolate(smoothing_plane,smoothed_state,smoothed_cov);
354 if(smoothing_plane == bPl) {
363 bAuxInfoP = &bAuxInfo;
367 rep->
setData(bUpSt, bPl, &bUpCov, bAuxInfoP);
373 TMatrixT<double> fCovInvert;
374 TMatrixT<double> bCovInvert;
381 smoothed_state.ResizeTo(smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt));
382 smoothed_state = smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt);
391 std::cout <<
"Trying to get smoothed hit coordinates from a track without smoothing! Aborting..." << std::endl;
392 TMatrixT<double> failed(1,1);
397 std::cout <<
"Hit number out of bounds while trying to get smoothed hit coordinates! Aborting..." <<std::endl;
398 TMatrixT<double> failed(1,1);
403 TMatrixT<double> bUpSt;
404 TMatrixT<double> bUpCov;
405 TMatrixT<double> bAuxInfo;
406 TMatrixT<double>* bAuxInfoP;
434 TMatrixT<double> fSt;
435 TMatrixT<double> fCov;
436 TMatrixT<double> bSt;
437 TMatrixT<double> bCov;
447 bAuxInfoP = &bAuxInfo;
454 if(bUpSt.GetNrows() == 0) {
460 rep->
setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
470 if(smoothing_plane == bPl) {
478 if (!extrapolation_allowed)
return false;
484 bAuxInfoP = &bAuxInfo;
491 if(bUpSt.GetNrows() == 0) {
495 rep->
setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
503 TMatrixT<double> fCovInvert;
504 TMatrixT<double> bCovInvert;
511 smoothed_state.ResizeTo(smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt));
512 smoothed_state = smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt);
529 GFException exc(
"Trying to get tracklength from a track without smoothing!",__LINE__,__FILE__);
534 GFException exc(
"Hit number out of bounds while trying to get tracklength!",__LINE__,__FILE__);
539 if(startHit > endHit) {
540 unsigned int biggerOne = startHit;
546 if (startHit==0 && endHit==0) endHit = trk->
getNumHits()-1;
547 if (startHit == endHit)
return 0.;
549 double totLen(0), fLen(0), bLen(0);
551 for (
unsigned int i=startHit; i!=endHit; ++i){
554 totLen += 0.5*(fLen - bLen);
557 if (inv)
return -totLen;
566 if (!(mat<1.E100) || !(mat>-1.E100)){
567 GFException e(
"cannot invert matrix GFTools::invertMatrix(), entries too big (>1e100)",
573 if (mat.GetNrows() == 2){
574 double det = mat(0,0)*mat(1,1) - mat(1,0)*mat(0,1);
575 if(fabs(det) < 1E-50){
576 GFException e(
"cannot invert matrix GFTools::invertMatrix(), determinant = 0",
582 inv(0,0) = det * mat(1,1);
583 inv(0,1) = -1.*det * mat(0,1);
584 inv(1,0) = -1.*det * mat(1,0);
585 inv(1,1) = det * mat(0,0);
591 TDecompSVD invertAlgo(mat);
593 invertAlgo.SetTol(1E-50);
595 inv = invertAlgo.Invert(status);
597 GFException e(
"cannot invert matrix GFTools::invertMatrix(), status = 0",
606 TMatrixT<double> smoothed_state, smoothed_cov, auxInfo;
613 rep->
setData(smoothed_state, smoothing_plane, &smoothed_cov, &auxInfo);
615 rep->
setData(smoothed_state, smoothing_plane, &smoothed_cov);
622 TMatrixT<double> smoothed_state;
623 TMatrixT<double> smoothed_cov;
626 TMatrixT<double> HT(TMatrixT<double>::kTransposed, H);
627 TMatrixT<double>
pos = H * smoothed_state;
630 TMatrixT<double> res = m -
pos;
631 TMatrixT<double> resT(TMatrixT<double>::kTransposed, res);
632 TMatrixT<double>
R = V - H*smoothed_cov*HT;
633 TMatrixT<double> invR;
635 TMatrixT<double> smoothedChiSqu = resT*invR*res;
636 return smoothedChiSqu(0,0);
642 GFException exc(
"Trying to get closest hit from a track without smoothing!",__LINE__,__FILE__);
647 if (nHits == 1)
return 0;
650 double minDist(1.E99), dist;
651 unsigned int minId(0);
653 if (checkEveryHit || nHits<4){
655 for (
unsigned int i=0; i<nHits; ++i){
657 dist = (pos-hitPos).Mag();
667 if (distFirst <= distLast){
670 for (
unsigned int i=1; i<nHits-1; ++i){
672 dist = (pos-hitPos).Mag();
683 for (
unsigned int i=nHits-2; i>1; --i){
685 dist = (pos-hitPos).Mag();