Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
WireMeasurementNew.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file WireMeasurementNew.cc
1 /* Copyright 2008-2010, Technische Universitaet Muenchen,
2  2014, Ludwig-Maximilians-Universität München
3  Authors: Tobias Schlüter
4 
5  This file is part of GENFIT.
6 
7  GENFIT is free software: you can redistribute it and/or modify
8  it under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  GENFIT is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU Lesser General Public License for more details.
16 
17  You should have received a copy of the GNU Lesser General Public License
18  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
19 */
24 #include "WireMeasurementNew.h"
25 
26 #include <cmath>
27 #include <algorithm>
28 
29 #include <Exception.h>
30 #include <RKTrackRep.h>
31 #include <HMatrixU.h>
32 
33 #include <cassert>
34 
35 
36 namespace genfit {
37 
38 
40  : AbsMeasurement(1), maxDistance_(2), leftRight_(0)
41 {
42  memset(wireEndPoint1_, 0, 3*sizeof(double));
43  memset(wireEndPoint2_, 0, 3*sizeof(double));
44 }
45 
46 WireMeasurementNew::WireMeasurementNew(double driftDistance, double driftDistanceError, const TVector3& endPoint1, const TVector3& endPoint2, int detId, int hitId, TrackPoint* trackPoint)
47  : AbsMeasurement(1), maxDistance_(2), leftRight_(0)
48 {
49  TVectorD coords(1);
50  coords[0] = driftDistance;
51  this->setRawHitCoords(coords);
52 
53  TMatrixDSym cov(1);
54  cov(0,0) = driftDistanceError*driftDistanceError;
55  this->setRawHitCov(cov);
56 
57  this->setWireEndPoints(endPoint1, endPoint2);
58 
59  this->setDetId(detId);
60  this->setHitId(hitId);
61  this->setTrackPoint(trackPoint);
62 }
63 
65 
66  // copy state. Neglect covariance.
67  StateOnPlane st(state);
68 
69  TVector3 wire1(wireEndPoint1_);
70  TVector3 wire2(wireEndPoint2_);
71 
72  // unit vector along the wire (V)
73  TVector3 wireDirection = wire2 - wire1;
74  wireDirection.SetMag(1.);
75 
76  // point of closest approach
77  const AbsTrackRep* rep = state.getRep();
78  rep->extrapolateToLine(st, wire1, wireDirection);
79  const TVector3& poca = rep->getPos(st);
80  TVector3 dirInPoca = rep->getMom(st);
81  dirInPoca.SetMag(1.);
82  const TVector3& pocaOnWire = wire1 + wireDirection.Dot(poca - wire1)*wireDirection;
83 
84  // check if direction is parallel to wire
85  if (fabs(wireDirection.Angle(dirInPoca)) < 0.01){
86  Exception exc("WireMeasurementNew::detPlane(): Cannot construct detector plane, direction is parallel to wire", __LINE__,__FILE__);
87  throw exc;
88  }
89 
90  // construct orthogonal vector
91  TVector3 U = wireDirection.Cross(dirInPoca);
92  // U.SetMag(1.); automatically assured
93 
94  return SharedPlanePtr(new DetPlane(pocaOnWire, U, wireDirection));
95 }
96 
97 
98 std::vector<MeasurementOnPlane*> WireMeasurementNew::constructMeasurementsOnPlane(const StateOnPlane& state) const
99 {
100  double mR = getRawHitCoords()(0);
101  double mL = -mR;
102 
103  MeasurementOnPlane* mopL = new MeasurementOnPlane(TVectorD(1, &mL),
104  getRawHitCov(),
105  state.getPlane(), state.getRep(), constructHMatrix(state.getRep()));
106 
107  MeasurementOnPlane* mopR = new MeasurementOnPlane(TVectorD(1, &mR),
108  getRawHitCov(),
109  state.getPlane(), state.getRep(), constructHMatrix(state.getRep()));
110 
111  // set left/right weights
112  if (leftRight_ < 0) {
113  mopL->setWeight(1);
114  mopR->setWeight(0);
115  }
116  else if (leftRight_ > 0) {
117  mopL->setWeight(0);
118  mopR->setWeight(1);
119  }
120  else {
121  double val = 0.5 * pow(std::max(0., 1 - mR/maxDistance_), 2.);
122  mopL->setWeight(val);
123  mopR->setWeight(val);
124  }
125 
126  std::vector<MeasurementOnPlane*> retVal;
127  retVal.push_back(mopL);
128  retVal.push_back(mopR);
129  return retVal;
130 }
131 
133  if (dynamic_cast<const RKTrackRep*>(rep) == nullptr) {
134  Exception exc("WireMeasurementNew default implementation can only handle state vectors of type RKTrackRep!", __LINE__,__FILE__);
135  throw exc;
136  }
137 
138  return new HMatrixU();
139 }
140 
141 void WireMeasurementNew::setWireEndPoints(const TVector3& endPoint1, const TVector3& endPoint2)
142 {
143  wireEndPoint1_[0] = endPoint1.X();
144  wireEndPoint1_[1] = endPoint1.Y();
145  wireEndPoint1_[2] = endPoint1.Z();
146 
147  wireEndPoint2_[0] = endPoint2.X();
148  wireEndPoint2_[1] = endPoint2.Y();
149  wireEndPoint2_[2] = endPoint2.Z();
150 }
151 
153  if (lr==0) leftRight_ = 0;
154  else if (lr<0) leftRight_ = -1;
155  else leftRight_ = 1;
156 }
157 
158 
159 } /* End of namespace genfit */