11 #include <cdc/dataobjects/CDCRecoHit.h>
15 #include <genfit/WireTrackCandHit.h>
16 #include <genfit/AbsTrackRep.h>
17 #include <genfit/RKTrackRep.h>
18 #include <genfit/TrackPoint.h>
19 #include <genfit/AbsFitterInfo.h>
20 #include <genfit/Exception.h>
21 #include <genfit/HMatrixU.h>
27 std::unique_ptr<ADCCountTranslatorBase> CDCRecoHit::s_adcCountTranslator = 0;
28 std::unique_ptr<CDCGeometryTranslatorBase> CDCRecoHit::s_cdcGeometryTranslator = 0;
29 std::unique_ptr<TDCCountTranslatorBase> CDCRecoHit::s_tdcCountTranslator = 0;
30 bool CDCRecoHit::s_useTrackTime =
false;
32 bool CDCRecoHit::s_cosmics =
false;
39 bool useTrackTime,
bool cosmics)
41 s_adcCountTranslator.reset(adcCountTranslator);
42 s_cdcGeometryTranslator.reset(cdcGeometryTranslator);
43 s_tdcCountTranslator.reset(tdcCountTranslator);
44 s_useTrackTime = useTrackTime;
49 CDCRecoHit::CDCRecoHit()
50 :
genfit::AbsMeasurement(1),
51 m_tdcCount(0), m_adcCount(0), m_wireID(
WireID()), m_cdcHit(NULL), m_leftRight(0)
56 :
genfit::AbsMeasurement(1), m_cdcHit(cdcHit), m_leftRight(0)
59 B2FATAL(
"Can't produce CDCRecoHits without setting of the translators.");
69 if (aTrackCandHitPtr) {
70 signed char lrInfo = aTrackCandHitPtr->getLeftRightResolution();
71 B2DEBUG(250,
"l/r: " <<
int(lrInfo));
97 TVector3 noSagWireDirection = noSagWire2 - noSagWire1;
98 noSagWireDirection.SetMag(1.);
103 const TVector3& noSagPoca = rep->
getPos(st);
105 double zPOCA = (noSagWire1.Z()
106 + noSagWireDirection.Dot(noSagPoca - noSagWire1) * noSagWireDirection.Z());
113 TVector3 wireDirection = wire2 - wire1;
114 wireDirection.SetMag(1.);
118 const TVector3& poca = rep->
getPos(st);
119 TVector3 dirInPoca = rep->
getMom(st);
120 dirInPoca.SetMag(1.);
121 const TVector3& pocaOnWire = wire1 + wireDirection.Dot(poca - wire1) * wireDirection;
128 if (fabs(wireDirection.Angle(dirInPoca)) < 0.01) {
129 genfit::Exception exc(
"CDCRecoHit::constructPlane(): Cannot construct detector plane, direction is parallel to wire", __LINE__,
135 const TVector3& U = wireDirection.Cross(dirInPoca);
144 double z = state.getPos().Z();
145 const TVector3& p = state.getMom();
154 const double wx = state.getPlane()->getO().x();
155 const double wy = state.getPlane()->getO().y();
156 const double px = p.x();
157 const double py = p.y();
158 const double cross = wx * py - wy * px;
159 const double dot = wx * px + wy * py;
160 double alpha = atan2(cross, dot);
161 double theta = atan2(p.Perp(), p.z());
179 if (atan2(py, px) > 0.) {
202 static TVectorD m(1);
203 static TMatrixDSym cov(1);
236 double weight = 0.5 * pow(std::max(0., 1 - (mR + mL) / 2 / rMax), 2);
237 mopL->setWeight(weight);
238 mopR->setWeight(weight);
243 if (mL < 0. || mR < 0.) {
244 B2DEBUG(150,
"Ignoring hit with negative drift time.");
256 B2FATAL(
"CDCRecoHit can only handle state vectors of type RKTrackRep!");
265 double z = state.getPos().Z();
266 const TVector3& p = state.getMom();
275 const double wx = state.getPlane()->getO().x();
276 const double wy = state.getPlane()->getO().y();
277 const double px = p.x();
278 const double py = p.y();
279 const double cross = wx * py - wy * px;
280 const double dot = wx * px + wy * py;
281 double alpha = atan2(cross, dot);
282 double theta = atan2(p.Perp(), p.z());
301 auto fL = [&](
const double & t) ->
double {
305 auto fR = [&](
const double & t) ->
double {
318 double rightShort[2], rightFull[2];
319 double leftShort[2], leftFull[2];
320 const double defaultStepT = 1e-3 * trackTime;
323 double temp = trackTime + defaultStepT / 2;
337 stepT = 2 * (temp - trackTime);
339 rightShort[0] = fL(trackTime + .5 * stepT);
340 rightShort[1] = fR(trackTime + .5 * stepT);
343 leftShort[0] = fL(trackTime - .5 * stepT);
344 leftShort[1] = fR(trackTime - .5 * stepT);
347 rightFull[0] = fL(trackTime + stepT);
348 rightFull[1] = fR(trackTime + stepT);
351 leftFull[0] = fL(trackTime - stepT);
352 leftFull[1] = fR(trackTime - stepT);
358 double derivShort[2];
359 for (
size_t j = 0; j < 2; ++j) {
360 derivFull[j] = (rightFull[j] - leftFull[j]) / (2.*stepT);
361 derivShort[j] = (rightShort[j] - leftShort[j]) / stepT;
364 return { +(4.*derivShort[0] - derivFull[0]) / 3.,
365 -(4.*derivShort[1] - derivFull[1]) / 3.};
373 bool usePlaneFromFit)
377 B2ERROR(
"No genfit::TrackPoint for CDCRecoHit.");
382 B2DEBUG(200,
"No genfit::AbsFitterInfo for this CDCRecoHit.");
394 ? mop.getPlane()->getO()
398 pointingVector = fittedPoca - pocaOnWire;
400 trackDir = mop.getMom();