measurements/src/WirePointMeasurement.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 "WirePointMeasurement.h"
0021
0022 #include <Exception.h>
0023 #include <RKTrackRep.h>
0024 #include <HMatrixUV.h>
0025
0026 #include <cassert>
0027 #include <algorithm>
0028
0029
0030 namespace genfit {
0031
0032 WirePointMeasurement::WirePointMeasurement(int nDim)
0033   : WireMeasurement(nDim)
0034 {
0035   assert(nDim >= 8);
0036 }
0037
0038 WirePointMeasurement::WirePointMeasurement(const TVectorD& rawHitCoords, const TMatrixDSym& rawHitCov, int detId, int hitId, TrackPoint* trackPoint)
0039   : WireMeasurement(rawHitCoords, rawHitCov, detId, hitId, trackPoint)
0040 {
0041   assert(rawHitCoords_.GetNrows() >= 8);
0042 }
0043
0044
0045 SharedPlanePtr WirePointMeasurement::constructPlane(const StateOnPlane& state) const {
0046
0047   // copy state. Neglect covariance.
0048   StateOnPlane st(state);
0049
0050   TVector3 wire1(rawHitCoords_(0), rawHitCoords_(1), rawHitCoords_(2));
0051   TVector3 wire2(rawHitCoords_(3), rawHitCoords_(4), rawHitCoords_(5));
0052
0053   // unit vector along the wire (V)
0054   TVector3 wireDirection = wire2 - wire1;
0055   wireDirection.SetMag(1.);
0056
0057   // point of closest approach
0058   const AbsTrackRep* rep = state.getRep();
0059   rep->extrapolateToLine(st, wire1, wireDirection);
0060   //const TVector3& poca = rep->getPos(st);
0061   TVector3 dirInPoca = rep->getMom(st);
0062   dirInPoca.SetMag(1.);
0063   //const TVector3& pocaOnWire = wire1 + wireDirection.Dot(poca - wire1)*wireDirection;
0064
0065   // check if direction is parallel to wire
0066   if (fabs(wireDirection.Angle(dirInPoca)) < 0.01){
0067     Exception exc("WireMeasurement::detPlane(): Cannot construct detector plane, direction is parallel to wire", __LINE__,__FILE__);
0068     throw exc;
0069   }
0070
0071   // construct orthogonal vector
0072   TVector3 U = dirInPoca.Cross(wireDirection);
0073   // U.SetMag(1.); automatically assured
0074
0075   return SharedPlanePtr(new DetPlane(wire1, U, wireDirection));
0076 }
0077
0078
0079 std::vector<MeasurementOnPlane*> WirePointMeasurement::constructMeasurementsOnPlane(const StateOnPlane& state) const
0080 {
0081   MeasurementOnPlane* mopR = new MeasurementOnPlane(TVectorD(2),
0082        TMatrixDSym(2),
0083        state.getPlane(), state.getRep(), constructHMatrix(state.getRep()));
0084
0085   mopR->getState()(0) = rawHitCoords_(6);
0086   mopR->getState()(1) = rawHitCoords_(7);
0087
0088   mopR->getCov()(0,0) = rawHitCov_(6,6);
0089   mopR->getCov()(1,0) = rawHitCov_(7,6);
0090   mopR->getCov()(0,1) = rawHitCov_(6,7);
0091   mopR->getCov()(1,1) = rawHitCov_(7,7);
0092
0093
0094   MeasurementOnPlane* mopL = new MeasurementOnPlane(*mopR);
0095   mopL->getState()(0) *= -1;
0096
0097   // set left/right weights
0098   if (leftRight_ < 0) {
0099     mopL->setWeight(1);
0100     mopR->setWeight(0);
0101   }
0102   else if (leftRight_ > 0) {
0103     mopL->setWeight(0);
0104     mopR->setWeight(1);
0105   }
0106   else {
0107     double val = 0.5 * pow(std::max(0., 1 - rawHitCoords_(6)/maxDistance_), 2.);
0108     mopL->setWeight(val);
0109     mopR->setWeight(val);
0110   }
0111
0112   std::vector<MeasurementOnPlane*> retVal;
0113   retVal.push_back(mopL);
0114   retVal.push_back(mopR);
0115   return retVal;
0116 }
0117
0118 const AbsHMatrix* WirePointMeasurement::constructHMatrix(const AbsTrackRep* rep) const {
0119   if (dynamic_cast<const RKTrackRep*>(rep) == nullptr) {
0120     Exception exc("WirePointMeasurement default implementation can only handle state vectors of type RKTrackRep!", __LINE__,__FILE__);
0121     throw exc;
0122   }
0123
0124   return new HMatrixUV();
0125 }
0126
0127 } /* End of namespace genfit */