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