11 #include <alignment/reconstruction/AlignableCDCRecoHit.h>
13 #include <alignment/GlobalLabel.h>
14 #include <alignment/GlobalDerivatives.h>
15 #include <cdc/dataobjects/WireID.h>
16 #include <cdc/dbobjects/CDCTimeZeros.h>
17 #include <cdc/dbobjects/CDCTimeWalks.h>
18 #include <cdc/geometry/CDCGeometryPar.h>
23 using namespace alignment;
25 bool AlignableCDCRecoHit::s_enableTrackT0LocalDerivative =
true;
26 bool AlignableCDCRecoHit::s_enableWireSaggingGlobalDerivative =
false;
27 bool AlignableCDCRecoHit::s_enableWireByWireAlignmentGlobalDerivatives =
false;
29 std::pair<std::vector<int>, TMatrixD> AlignableCDCRecoHit::globalDerivatives(
const genfit::StateOnPlane* sop)
31 GlobalDerivatives globals;
32 unsigned short LR = (int(m_leftRight) > 0.) ? 1 : 0;
34 const TVector3& mom = sop->getMom();
35 const TVector3& wirePositon = sop->getPlane()->getO();
36 const unsigned short layer = getWireID().getICLayer();
39 const double alpha = cdcgeo.
getAlpha(wirePositon, mom);
40 const double theta = cdcgeo.
getTheta(mom);
41 const TVectorD& stateOnPlane = sop->getState();
42 const double driftLengthEstimate = std::abs(stateOnPlane(3));
43 const double driftTime = cdcgeo.
getDriftTime(driftLengthEstimate, layer, LR, alpha, theta);
44 const double driftVelocity = cdcgeo.
getDriftV(driftTime, layer, LR, alpha, theta);
49 if (driftTime > 20 && driftTime < 200 && fabs(driftVelocity) < 1.0e-2) {
51 GlobalLabel::construct<CDCTimeZeros>(getWireID(), 0),
52 driftVelocity *
double(
int(m_leftRight))
58 GlobalLabel::construct<CDCTimeWalks>(CDCGeometryPar::Instance().getBoardID(getWireID()), 0),
59 driftVelocity * 1. / sqrt(m_adcCount) *
double(
int(m_leftRight))
64 auto tdir = sop->getDir();
65 auto ndir = sop->getPlane()->getNormal();
66 auto udir = sop->getPlane()->getU();
67 auto vdir = sop->getPlane()->getV();
68 auto pos = sop->getPos();
70 auto tn = tdir[0] * ndir[0] + tdir[1] * ndir[1] + tdir[2] * ndir [2];
72 auto drdm = TMatrixD(3, 3);
74 for (
int i = 0; i < 3; ++i) {
75 for (
int j = 0; j < 3; ++j) {
76 drdm(i, j) -= tdir[i] * ndir[j] / tn;
80 auto dmdg = TMatrixD(3, 6);
82 dmdg[0][0] = 1.; dmdg[0][4] = -pos[2]; dmdg[0][5] = pos[1];
83 dmdg[1][1] = 1.; dmdg[1][3] = pos[2]; dmdg[1][5] = -pos[0];
84 dmdg[2][2] = 1.; dmdg[2][3] = -pos[1]; dmdg[2][4] = pos[0];
86 auto drldrg = TMatrixD(3, 3);
87 for (
int i = 0; i < 3; ++i) {
88 drldrg(0, i) = udir[i];
89 drldrg(1, i) = vdir[i];
90 drldrg(2, i) = ndir[i];
93 auto drldg = drldrg * (drdm * dmdg);
96 const double zWireM = s_cdcGeometryTranslator->getWireBackwardPosition(getWireID(), CDCGeometryPar::c_Aligned)[2];
97 const double zWireP = s_cdcGeometryTranslator->getWireForwardPosition(getWireID(), CDCGeometryPar::c_Aligned)[2];
99 const double zRel = std::max(0., std::min(1., (pos[2] - zWireM) / (zWireP - zWireM)));
103 auto layerID =
WireID(getWireID().getICLayer(), 511);
107 GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerX),
113 GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerY),
119 GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerPhi),
126 GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerDx),
132 GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerDy),
138 GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerDPhi),
144 if (s_enableWireByWireAlignmentGlobalDerivatives) {
147 double zRelM = fabs(1. - zRel);
149 double zRelP = fabs(zRel - 0.);
153 GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireBwdX),
158 GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireFwdX),
164 GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireBwdY),
169 GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireFwdY),
178 if (s_enableWireSaggingGlobalDerivative) {
180 GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireTension),
181 drldg(0, 1) * 4.0 * zRel * (1.0 - zRel)
191 if (!s_enableTrackT0LocalDerivative)
194 unsigned short LR = (int(m_leftRight) > 0.) ? 1 : 0;
196 const TVector3& mom = sop->getMom();
197 const TVector3& wirePositon = sop->getPlane()->getO();
198 const unsigned short layer = getWireID().getICLayer();
201 const double alpha = cdcgeo.
getAlpha(wirePositon, mom);
202 const double theta = cdcgeo.
getTheta(mom);
203 const TVectorD& stateOnPlane = sop->getState();
204 const double driftLengthEstimate = std::abs(stateOnPlane(3));
205 const double driftTime = cdcgeo.
getDriftTime(driftLengthEstimate, layer, LR, alpha, theta);
206 const double driftVelocity = cdcgeo.
getDriftV(driftTime, layer, LR, alpha, theta);
208 TMatrixD locals(2, 1);
209 if (driftTime > 20 && driftTime < 200 && fabs(driftVelocity) < 1.0e-2) {
210 locals(0, 0) = - double(
int(m_leftRight)) * driftVelocity;