measurements/src/WireMeasurementNew.cc

0001 /* Copyright 2008-2010, Technische Universitaet Muenchen,
0002              2014, Ludwig-Maximilians-Universität München
0003    Authors: Tobias Schlüter
0004
0005    This file is part of GENFIT.
0006
0007    GENFIT is free software: you can redistribute it and/or modify
0008    it under the terms of the GNU Lesser General Public License as published
0009    by the Free Software Foundation, either version 3 of the License, or
0010    (at your option) any later version.
0011
0012    GENFIT is distributed in the hope that it will be useful,
0013    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0015    GNU Lesser General Public License for more details.
0016
0017    You should have received a copy of the GNU Lesser General Public License
0018    along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
0019 */
0020 /** @addtogroup genfit
0021  * @{
0022  */
0023
0024 #include "WireMeasurementNew.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 WireMeasurementNew::WireMeasurementNew()
0040   : AbsMeasurement(1), maxDistance_(2), leftRight_(0)
0041 {
0042   memset(wireEndPoint1_, 0, 3*sizeof(double));
0043   memset(wireEndPoint2_, 0, 3*sizeof(double));
0044 }
0045
0046 WireMeasurementNew::WireMeasurementNew(double driftDistance, double driftDistanceError, const TVector3& endPoint1, const TVector3& endPoint2, int detId, int hitId, TrackPoint* trackPoint)
0047   : AbsMeasurement(1), maxDistance_(2), leftRight_(0)
0048 {
0049   TVectorD coords(1);
0050   coords[0] = driftDistance;
0051   this->setRawHitCoords(coords);
0052
0053   TMatrixDSym cov(1);
0054   cov(0,0) = driftDistanceError*driftDistanceError;
0055   this->setRawHitCov(cov);
0056
0057   this->setWireEndPoints(endPoint1, endPoint2);
0058
0059   this->setDetId(detId);
0060   this->setHitId(hitId);
0061   this->setTrackPoint(trackPoint);
0062 }
0063
0064 SharedPlanePtr WireMeasurementNew::constructPlane(const StateOnPlane& state) const {
0065
0066   // copy state. Neglect covariance.
0067   StateOnPlane st(state);
0068
0069   TVector3 wire1(wireEndPoint1_);
0070   TVector3 wire2(wireEndPoint2_);
0071
0072   // unit vector along the wire (V)
0073   TVector3 wireDirection = wire2 - wire1;
0074   wireDirection.SetMag(1.);
0075
0076   // point of closest approach
0077   const AbsTrackRep* rep = state.getRep();
0078   rep->extrapolateToLine(st, wire1, wireDirection);
0079   const TVector3& poca = rep->getPos(st);
0080   TVector3 dirInPoca = rep->getMom(st);
0081   dirInPoca.SetMag(1.);
0082   const TVector3& pocaOnWire = wire1 + wireDirection.Dot(poca - wire1)*wireDirection;
0083
0084   // check if direction is parallel to wire
0085   if (fabs(wireDirection.Angle(dirInPoca)) < 0.01){
0086     Exception exc("WireMeasurementNew::detPlane(): Cannot construct detector plane, direction is parallel to wire", __LINE__,__FILE__);
0087     throw exc;
0088   }
0089
0090   // construct orthogonal vector
0091   TVector3 U = wireDirection.Cross(dirInPoca);
0092   // U.SetMag(1.); automatically assured
0093
0094   return SharedPlanePtr(new DetPlane(pocaOnWire, U, wireDirection));
0095 }
0096
0097
0098 std::vector<MeasurementOnPlane*> WireMeasurementNew::constructMeasurementsOnPlane(const StateOnPlane& state) const
0099 {
0100   double mR = getRawHitCoords()(0);
0101   double mL = -mR;
0102
0103   MeasurementOnPlane* mopL = new MeasurementOnPlane(TVectorD(1, &mL),
0104        getRawHitCov(),
0105        state.getPlane(), state.getRep(), constructHMatrix(state.getRep()));
0106
0107   MeasurementOnPlane* mopR = new MeasurementOnPlane(TVectorD(1, &mR),
0108        getRawHitCov(),
0109        state.getPlane(), state.getRep(), constructHMatrix(state.getRep()));
0110
0111   // set left/right weights
0112   if (leftRight_ < 0) {
0113     mopL->setWeight(1);
0114     mopR->setWeight(0);
0115   }
0116   else if (leftRight_ > 0) {
0117     mopL->setWeight(0);
0118     mopR->setWeight(1);
0119   }
0120   else {
0121     double val = 0.5 * pow(std::max(0., 1 - mR/maxDistance_), 2.);
0122     mopL->setWeight(val);
0123     mopR->setWeight(val);
0124   }
0125
0126   std::vector<MeasurementOnPlane*> retVal;
0127   retVal.push_back(mopL);
0128   retVal.push_back(mopR);
0129   return retVal;
0130 }
0131
0132 const AbsHMatrix* WireMeasurementNew::constructHMatrix(const AbsTrackRep* rep) const {
0133   if (dynamic_cast<const RKTrackRep*>(rep) == nullptr) {
0134     Exception exc("WireMeasurementNew default implementation can only handle state vectors of type RKTrackRep!", __LINE__,__FILE__);
0135     throw exc;
0136   }
0137
0138   return new HMatrixU();
0139 }
0140
0141 void WireMeasurementNew::setWireEndPoints(const TVector3& endPoint1, const TVector3& endPoint2)
0142 {
0143   wireEndPoint1_[0] = endPoint1.X();
0144   wireEndPoint1_[1] = endPoint1.Y();
0145   wireEndPoint1_[2] = endPoint1.Z();
0146
0147   wireEndPoint2_[0] = endPoint2.X();
0148   wireEndPoint2_[1] = endPoint2.Y();
0149   wireEndPoint2_[2] = endPoint2.Z();
0150 }
0151
0152 void WireMeasurementNew::setLeftRightResolution(int lr){
0153   if (lr==0) leftRight_ = 0;
0154   else if (lr<0) leftRight_ = -1;
0155   else leftRight_ = 1;
0156 }
0157
0158
0159 } /* End of namespace genfit */
0160 /** @} */