Belle II Software  release-05-01-25
AlignableCDCRecoHit.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2012, 2015 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Tadeas Bilka *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <alignment/reconstruction/AlignableCDCRecoHit.h>
12 
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>
19 
20 using namespace std;
21 using namespace Belle2;
22 using namespace CDC;
23 using namespace alignment;
24 
25 bool AlignableCDCRecoHit::s_enableTrackT0LocalDerivative = true;
26 bool AlignableCDCRecoHit::s_enableWireSaggingGlobalDerivative = false;
27 bool AlignableCDCRecoHit::s_enableWireByWireAlignmentGlobalDerivatives = false;
28 
29 std::pair<std::vector<int>, TMatrixD> AlignableCDCRecoHit::globalDerivatives(const genfit::StateOnPlane* sop)
30 {
31  GlobalDerivatives globals;
32  unsigned short LR = (int(m_leftRight) > 0.) ? 1 : 0;
33 
34  const TVector3& mom = sop->getMom();
35  const TVector3& wirePositon = sop->getPlane()->getO();
36  const unsigned short layer = getWireID().getICLayer();
37 
38  CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
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);
45 
46  // CDC Calibration -------------------------------------------------
47 
48  // T0 calibration (per wire) TODO check sign!!!
49  if (driftTime > 20 && driftTime < 200 && fabs(driftVelocity) < 1.0e-2) {
50  globals.add(
51  GlobalLabel::construct<CDCTimeZeros>(getWireID(), 0),
52  driftVelocity * double(int(m_leftRight))
53  );
54  }
55 
56  // Time-walk calibration (per board) TODO checksign!!!
57  globals.add(
58  GlobalLabel::construct<CDCTimeWalks>(CDCGeometryPar::Instance().getBoardID(getWireID()), 0),
59  driftVelocity * 1. / sqrt(m_adcCount) * double(int(m_leftRight))
60  );
61 
62  // CDC Alignment ---------------------------------------------------
63 
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();
69 
70  auto tn = tdir[0] * ndir[0] + tdir[1] * ndir[1] + tdir[2] * ndir [2];
71  // d residual / d measurement
72  auto drdm = TMatrixD(3, 3);
73  drdm.UnitMatrix();
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;
77  }
78  }
79  // d measurement / d global rigid body param
80  auto dmdg = TMatrixD(3, 6);
81  dmdg.Zero();
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];
85  // d local residual / d global residual
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];
91  }
92  // d local residual / d global rigid body param
93  auto drldg = drldrg * (drdm * dmdg);
94 
95  // wire ends
96  const double zWireM = s_cdcGeometryTranslator->getWireBackwardPosition(getWireID(), CDCGeometryPar::c_Aligned)[2];
97  const double zWireP = s_cdcGeometryTranslator->getWireForwardPosition(getWireID(), CDCGeometryPar::c_Aligned)[2];
98  // relative Z position [0..1]
99  const double zRel = std::max(0., std::min(1., (pos[2] - zWireM) / (zWireP - zWireM)));
100 
101  // Layer alignment
102  // wire 511 = no wire (0 is reserved for existing wires) - this has to be compatible with code in CDCGeometryPar::setWirPosAlignParams
103  auto layerID = WireID(getWireID().getICLayer(), 511);
104 
105  // Alignment of layer X (bwd)
106  globals.add(
107  GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerX),
108  drldg(0, 0)
109  );
110 
111  // Alignment of layer Y (bwd)
112  globals.add(
113  GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerY),
114  drldg(0, 1)
115  );
116 
117  // Alignment of layer rotation (gamma) (bwd)
118  globals.add(
119  GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerPhi),
120  drldg(0, 5)
121  );
122 
123  // Difference between wire ends (end plates)
124  // Alignment of layer dX, dX = foward - backward endplate
125  globals.add(
126  GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerDx),
127  drldg(0, 0) * zRel
128  );
129 
130  // Alignment of layer dY, dY = foward - backward endplate
131  globals.add(
132  GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerDy),
133  drldg(0, 1) * zRel
134  );
135 
136  // Alignment of layer rotation difference d(gamma or phi), dPhi = foward - backward endplate
137  globals.add(
138  GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerDPhi),
139  drldg(0, 5) * zRel
140  );
141 
142  //WARNING: experimental (disabled by default)
143  // Wire-by-wire alignment
144  if (s_enableWireByWireAlignmentGlobalDerivatives) {
145  // How much shift (in X or Y) on BWD wire-end will change the residual at estimated track crossing
146  // with the wire (at relative z-position on wire = zRel)
147  double zRelM = fabs(1. - zRel);
148  // Same as above but for FWD wire-end (residual at zRel = zRel * delta(X or Y at FWD wire-end)
149  double zRelP = fabs(zRel - 0.);
150 
151  // Alignment of wires X in global coords at BWD wire-end
152  globals.add(
153  GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireBwdX),
154  drldg(0, 0) * zRelM
155  );
156  // Alignment of wires X in global coords at FWD wire-end
157  globals.add(
158  GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireFwdX),
159  drldg(0, 0) * zRelP
160  );
161 
162  // Alignment of wires Y in global coords at BWD wire-end
163  globals.add(
164  GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireBwdY),
165  drldg(0, 1) * zRelM
166  );
167  // Alignment of wires Y in global coords at FWD wire-end
168  globals.add(
169  GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireFwdY),
170  drldg(0, 1) * zRelP
171  );
172  }
173 
174  //WARNING: experimental (disabled by default)
175  //TODO: missing some factors (we do not align directly the wire-sag coefficient, but
176  // wire tension ... coef = pi * ro * r * r / 8 / tension
177  //TODO: need to get these numbers from CDCGeometryPar!
178  if (s_enableWireSaggingGlobalDerivative) {
179  globals.add(
180  GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireTension),
181  drldg(0, 1) * 4.0 * zRel * (1.0 - zRel)
182  );
183  }
184 
185  return globals;
186 }
187 
188 
189 TMatrixD AlignableCDCRecoHit::localDerivatives(const genfit::StateOnPlane* sop)
190 {
191  if (!s_enableTrackT0LocalDerivative)
192  return TMatrixD();
193 
194  unsigned short LR = (int(m_leftRight) > 0.) ? 1 : 0;
195 
196  const TVector3& mom = sop->getMom();
197  const TVector3& wirePositon = sop->getPlane()->getO();
198  const unsigned short layer = getWireID().getICLayer();
199 
200  CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
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);
207 
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;
211  locals(1, 0) = 0.; // insesitive coordinate along wire
212  }
213 
214  return locals;
215 }
Belle2::WireID
Class to identify a wire inside the CDC.
Definition: WireID.h:44
genfit::StateOnPlane
A state with arbitrary dimension defined in a DetPlane.
Definition: StateOnPlane.h:47
Belle2::CDC::CDCGeometryPar
The Class for CDC Geometry Parameters.
Definition: CDCGeometryPar.h:75
Belle2::CDC::CDCGeometryPar::getTheta
double getTheta(const TVector3 &momentum) const
Returns track incident angle (theta in rad.).
Definition: CDCGeometryPar.cc:2674
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDC::CDCGeometryPar::getDriftTime
double getDriftTime(double dist, unsigned short layer, unsigned short lr, double alpha, double theta) const
Return the drift time to the sense wire.
Definition: CDCGeometryPar.cc:2457
Belle2::CDC::CDCGeometryPar::getDriftV
double getDriftV(double dt, unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI) const
Get the realistic drift velocity.
Definition: CDCGeometryPar.cc:2006
Belle2::CDC::CDCGeometryPar::getAlpha
double getAlpha(const TVector3 &posOnWire, const TVector3 &momentum) const
Returns track incident angle in rphi plane (alpha in rad.).
Definition: CDCGeometryPar.cc:2661