Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
WireMeasurement.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file WireMeasurement.cc
1 /* Copyright 2008-2010, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
23 #include "WireMeasurement.h"
24 #include "IO.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(nDim), maxDistance_(2), leftRight_(0)
41 {
42  assert(nDim >= 7);
43 }
44 
45 WireMeasurement::WireMeasurement(const TVectorD& rawHitCoords, const TMatrixDSym& rawHitCov, int detId, int hitId, TrackPoint* trackPoint)
46  : AbsMeasurement(rawHitCoords, rawHitCov, detId, hitId, trackPoint), maxDistance_(2), leftRight_(0)
47 {
48  assert(rawHitCoords_.GetNrows() >= 7);
49 }
50 
52 
53  // copy state. Neglect covariance.
54  StateOnPlane st(state);
55 
56  TVector3 wire1(rawHitCoords_(0), rawHitCoords_(1), rawHitCoords_(2));
57  TVector3 wire2(rawHitCoords_(3), rawHitCoords_(4), rawHitCoords_(5));
58 
59  // unit vector along the wire (V)
60  TVector3 wireDirection = wire2 - wire1;
61  wireDirection.SetMag(1.);
62 
63  // point of closest approach
64  const AbsTrackRep* rep = state.getRep();
65  rep->extrapolateToLine(st, wire1, wireDirection);
66  const TVector3& poca = rep->getPos(st);
67  TVector3 dirInPoca = rep->getMom(st);
68  dirInPoca.SetMag(1.);
69  const TVector3& pocaOnWire = wire1 + wireDirection.Dot(poca - wire1)*wireDirection;
70 
71  // check if direction is parallel to wire
72  if (fabs(wireDirection.Angle(dirInPoca)) < 0.01){
73  Exception exc("WireMeasurement::detPlane(): Cannot construct detector plane, direction is parallel to wire", __LINE__,__FILE__);
74  throw exc;
75  }
76 
77  // construct orthogonal vector
78  TVector3 U = dirInPoca.Cross(wireDirection);
79  // U.SetMag(1.); automatically assured
80 
81  return SharedPlanePtr(new DetPlane(pocaOnWire, U, wireDirection));
82 }
83 
84 
85 std::vector<MeasurementOnPlane*> WireMeasurement::constructMeasurementsOnPlane(const StateOnPlane& state) const
86 {
87  double mR = rawHitCoords_(6);
88  double mL = -mR;
89  double V = rawHitCov_(6,6);
90 
91  MeasurementOnPlane* mopL = new MeasurementOnPlane(TVectorD(1, &mL),
92  TMatrixDSym(1, &V),
93  state.getPlane(), state.getRep(), constructHMatrix(state.getRep()));
94 
95  MeasurementOnPlane* mopR = new MeasurementOnPlane(TVectorD(1, &mR),
96  TMatrixDSym(1, &V),
97  state.getPlane(), state.getRep(), constructHMatrix(state.getRep()));
98 
99  // set left/right weights
100  if (leftRight_ < 0) {
101  mopL->setWeight(1);
102  mopR->setWeight(0);
103  }
104  else if (leftRight_ > 0) {
105  mopL->setWeight(0);
106  mopR->setWeight(1);
107  }
108  else {
109  double val = 0.5 * pow(std::max(0., 1 - mR/maxDistance_), 2.);
110  mopL->setWeight(val);
111  mopR->setWeight(val);
112  }
113 
114  std::vector<MeasurementOnPlane*> retVal;
115  retVal.push_back(mopL);
116  retVal.push_back(mopR);
117  return retVal;
118 }
119 
121  if (dynamic_cast<const RKTrackRep*>(rep) == nullptr) {
122  Exception exc("WireMeasurement default implementation can only handle state vectors of type RKTrackRep!", __LINE__,__FILE__);
123  throw exc;
124  }
125 
126  return new HMatrixU();
127 }
128 
130  if (lr==0) leftRight_ = 0;
131  else if (lr<0) leftRight_ = -1;
132  else leftRight_ = 1;
133 }
134 
135 
136 } /* End of namespace genfit */