SND@LHC Software
Loading...
Searching...
No Matches
MeasuredStateOnPlane.cc
Go to the documentation of this file.
1/* Copyright 2008-2010, Technische Universitaet Muenchen,
2 Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
3
4 This file is part of GENFIT.
5
6 GENFIT is free software: you can redistribute it and/or modify
7 it under the terms of the GNU Lesser General Public License as published
8 by the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 GENFIT is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU Lesser General Public License for more details.
15
16 You should have received a copy of the GNU Lesser General Public License
17 along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18*/
19
21#include "AbsTrackRep.h"
22#include "Exception.h"
23#include "Tools.h"
24
25#include <cassert>
26#include <iostream>
27
28namespace genfit {
29
30void MeasuredStateOnPlane::Print(Option_t*) const {
31 std::cout << "genfit::MeasuredStateOnPlane ";
32 std::cout << "my address " << this << " my plane's address " << this->sharedPlane_.get() << "; use count: " << sharedPlane_.use_count() << std::endl;
33 std::cout << " state vector: "; state_.Print();
34 std::cout << " covariance matrix: "; cov_.Print();
35 if (sharedPlane_ != NULL) {
36 std::cout << " defined in plane "; sharedPlane_->Print();
37 TVector3 pos, mom;
38 TMatrixDSym cov(6,6);
39 getRep()->getPosMomCov(*this, pos, mom, cov);
40 std::cout << " 3D position: "; pos.Print();
41 std::cout << " 3D momentum: "; mom.Print();
42 //std::cout << " 6D covariance: "; cov.Print();
43 }
44}
45
46void MeasuredStateOnPlane::blowUpCov(double blowUpFac, bool resetOffDiagonals) {
47
48 if (resetOffDiagonals) {
49 unsigned int dim = cov_.GetNcols();
50 for (unsigned int i=0; i<dim; ++i) {
51 for (unsigned int j=0; j<dim; ++j) {
52 if (i != j)
53 cov_(i,j) = 0; // reset off-diagonals
54 else
55 cov_(i,j) *= blowUpFac; // blow up diagonals
56 }
57 }
58 }
59 else
60 cov_ *= blowUpFac;
61
62}
63
64
66 // check if both states are defined in the same plane
67 if (forwardState.getPlane() != backwardState.getPlane()) {
68 Exception e("KalmanFitterInfo::calcAverageState: forwardState and backwardState are not defined in the same plane.", __LINE__,__FILE__);
69 throw e;
70 }
71
72 // This code is called a lot, so some effort has gone into reducing
73 // the number of temporary objects being constructed.
74
75#if 0
76 // For ease of understanding, here's a version of the code that is a
77 // few percent slower, depending on the exact use-case, but which
78 // makes the math more explicit:
79 TMatrixDSym fCovInv, bCovInv, smoothed_cov;
80 tools::invertMatrix(forwardState.getCov(), fCovInv);
81 tools::invertMatrix(backwardState.getCov(), bCovInv);
82
83 tools::invertMatrix(fCovInv + bCovInv, smoothed_cov); // one temporary TMatrixDSym
84
85 MeasuredStateOnPlane retVal(forwardState);
86 retVal.setState(smoothed_cov*(fCovInv*forwardState.getState() + bCovInv*backwardState.getState())); // four temporary TVectorD's
87 retVal.setCov(smoothed_cov);
88 return retVal;
89#endif
90
91 static TMatrixDSym fCovInv, bCovInv; // Static to avoid re-constructing for every call
92 tools::invertMatrix(forwardState.getCov(), fCovInv);
93 tools::invertMatrix(backwardState.getCov(), bCovInv);
94
95 // Using a StateOnPlane here turned out at least as fast as
96 // constructing a MeasuredStateOnPlane here, and then resetting its
97 // covariance matrix below, even though it means another copy in the
98 // return statement.
99 StateOnPlane sop(forwardState); // copies auxInfo, plane, rep in the process
100 // Using 'static' + subsequent assignment is measurably slower.
101 // Surprisingly, using only TVectorD
102 // sop(forwardState.getState()) with according
103 // changes below measured slower.
104
105 sop.getState() *= fCovInv;
106 fCovInv += bCovInv;
107 tools::invertMatrix(fCovInv); // This is now the covariance of the average.
108 sop.getState() += bCovInv*backwardState.getState(); // one temporary TVectorD
109 sop.getState() *= fCovInv;
110
111 return MeasuredStateOnPlane(sop, fCovInv);
112}
113
114
115} /* End of namespace genfit */
virtual void getPosMomCov(const MeasuredStateOnPlane &state, TVector3 &pos, TVector3 &mom, TMatrixDSym &cov) const =0
Translates MeasuredStateOnPlane into 3D position, momentum and 6x6 covariance.
Exception class for error handling in GENFIT (provides storage for diagnostic information)
Definition Exception.h:48
StateOnPlane with additional covariance matrix.
void setCov(const TMatrixDSym &cov)
const TMatrixDSym & getCov() const
void blowUpCov(double blowUpFac, bool resetOffDiagonals=true)
virtual void Print(Option_t *option="") const
A state with arbitrary dimension defined in a DetPlane.
const TVectorD & getState() const
SharedPlanePtr sharedPlane_
void setState(const TVectorD &state)
const AbsTrackRep * getRep() const
const SharedPlanePtr & getPlane() const
void invertMatrix(const TMatrixDSym &mat, TMatrixDSym &inv, double *determinant=NULL)
Invert a matrix, throwing an Exception when inversion fails. Optional calculation of determinant.
Definition Tools.cc:38
Matrix inversion tools.
Definition AbsBField.h:29
MeasuredStateOnPlane calcAverageState(const MeasuredStateOnPlane &forwardState, const MeasuredStateOnPlane &backwardState)
Calculate weighted average between two MeasuredStateOnPlanes.