measurements/src/SpacepointMeasurement.cc

0001 /* Copyright 2008-2010, Technische Universitaet Muenchen,
0002    Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
0003
0004    This file is part of GENFIT.
0005
0006    GENFIT is free software: you can redistribute it and/or modify
0007    it under the terms of the GNU Lesser General Public License as published
0008    by the Free Software Foundation, either version 3 of the License, or
0009    (at your option) any later version.
0010
0011    GENFIT is distributed in the hope that it will be useful,
0012    but WITHOUT ANY WARRANTY; without even the implied warranty of
0013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0014    GNU Lesser General Public License for more details.
0015
0016    You should have received a copy of the GNU Lesser General Public License
0017    along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
0018 */
0019
0020 #include "SpacepointMeasurement.h"
0021
0022 #include "Exception.h"
0023 #include "RKTrackRep.h"
0024 #include "Tools.h"
0025 #include "HMatrixUV.h"
0026 #include "MeasurementOnPlane.h"
0027
0028 #include <cassert>
0029 #include <TBuffer.h>
0030
0031 namespace genfit {
0032
0033 SpacepointMeasurement::SpacepointMeasurement(int nDim)
0034   : AbsMeasurement(nDim), weightedPlaneContruction_(true), G_(3), cutCov_(true)
0035 {
0036   assert(nDim >= 3);
0037
0038   //this->initG();
0039 }
0040
0041 SpacepointMeasurement::SpacepointMeasurement(const TVectorD& rawHitCoords, const TMatrixDSym& rawHitCov, int detId, int hitId, TrackPoint* trackPoint,
0042     bool weightedPlaneContruction, bool cutCov)
0043   : AbsMeasurement(rawHitCoords, rawHitCov, detId, hitId, trackPoint),
0044     weightedPlaneContruction_(weightedPlaneContruction), cutCov_(cutCov)
0045 {
0046   assert(rawHitCoords_.GetNrows() >= 3);
0047
0048   if (weightedPlaneContruction_)
0049     this->initG();
0050 }
0051
0052
0053 SharedPlanePtr SpacepointMeasurement::constructPlane(const StateOnPlane& state) const {
0054
0055   // copy state. Neglect covariance.
0056   StateOnPlane st(state);
0057
0058   const TVector3 point(rawHitCoords_(0), rawHitCoords_(1), rawHitCoords_(2));
0059
0060   if (weightedPlaneContruction_)
0061     st.extrapolateToPoint(point, G_);
0062   else
0063     st.extrapolateToPoint(point);
0064
0065   return st.getPlane();
0066 }
0067
0068
0069 std::vector<MeasurementOnPlane*> SpacepointMeasurement::constructMeasurementsOnPlane(const StateOnPlane& state) const
0070 {
0071   MeasurementOnPlane* mop = new MeasurementOnPlane(TVectorD(2),
0072        TMatrixDSym(3), // will be resized to 2x2 by Similarity later
0073        state.getPlane(), state.getRep(), constructHMatrix(state.getRep()));
0074
0075   TVectorD& m = mop->getState();
0076   TMatrixDSym& V = mop->getCov();
0077
0078   const TVector3& o(state.getPlane()->getO());
0079   const TVector3& u(state.getPlane()->getU());
0080   const TVector3& v(state.getPlane()->getV());
0081
0082   // m
0083   m(0) = (rawHitCoords_(0)-o.X()) * u.X() +
0084          (rawHitCoords_(1)-o.Y()) * u.Y() +
0085          (rawHitCoords_(2)-o.Z()) * u.Z();
0086
0087   m(1) = (rawHitCoords_(0)-o.X()) * v.X() +
0088          (rawHitCoords_(1)-o.Y()) * v.Y() +
0089          (rawHitCoords_(2)-o.Z()) * v.Z();
0090
0091   //
0092   // V
0093   //
0094   V = rawHitCov_;
0095
0096   // jac = dF_i/dx_j = s_unitvec * t_untivec, with s=u,v and t=x,y,z
0097   TMatrixD jac(3,2);
0098   jac(0,0) = u.X();
0099   jac(1,0) = u.Y();
0100   jac(2,0) = u.Z();
0101   jac(0,1) = v.X();
0102   jac(1,1) = v.Y();
0103   jac(2,1) = v.Z();
0104
0105   if (cutCov_) { // cut (default)
0106     tools::invertMatrix(V);
0107     V.SimilarityT(jac);
0108     tools::invertMatrix(V);
0109   }
0110   else { // projection
0111     V.SimilarityT(jac);
0112   }
0113
0114   std::vector<MeasurementOnPlane*> retVal;
0115   retVal.push_back(mop);
0116   return retVal;
0117 }
0118
0119
0120 const AbsHMatrix* SpacepointMeasurement::constructHMatrix(const AbsTrackRep* rep) const {
0121   if (dynamic_cast<const RKTrackRep*>(rep) == nullptr) {
0122     Exception exc("SpacepointMeasurement default implementation can only handle state vectors of type RKTrackRep!", __LINE__,__FILE__);
0123     throw exc;
0124   }
0125
0126   return new HMatrixUV();
0127 }
0128
0129
0130 void SpacepointMeasurement::initG() {
0131   rawHitCov_.GetSub(0, 2, G_);
0132   tools::invertMatrix(G_);
0133 }
0134
0135
0136 // Modified from auto-generated Streamer to account for non-persistent G_
0137 void SpacepointMeasurement::Streamer(TBuffer &R__b)
0138 {
0139    // Stream an object of class genfit::SpacepointMeasurement.
0140
0141    if (R__b.IsReading()) {
0142       R__b.ReadClassBuffer(genfit::SpacepointMeasurement::Class(),this);
0143
0144       if (weightedPlaneContruction_)
0145         this->initG();
0146    } else {
0147       R__b.WriteClassBuffer(genfit::SpacepointMeasurement::Class(),this);
0148    }
0149 }
0150
0151
0152 } /* End of namespace genfit */