Belle II Software  release-05-01-25
SVDRecoHit.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Peter Kvasnicka, Martin Ritter, Moritz Nadler *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <framework/logging/Logger.h>
12 #include <svd/reconstruction/SVDRecoHit.h>
13 #include <svd/geometry/SensorInfo.h>
14 #include <vxd/geometry/SensorPlane.h>
15 #include <vxd/geometry/GeoCache.h>
16 
17 #include <genfit/DetPlane.h>
18 #include <TVector3.h>
19 #include <TRandom.h>
20 #include <cmath>
21 
22 using namespace std;
23 using namespace Belle2;
24 
25 SVDRecoHit::SVDRecoHit():
26  genfit::PlanarMeasurement(HIT_DIMENSIONS), m_sensorID(0), m_trueHit(0),
27  m_cluster(0), m_isU(0), m_energyDep(0), m_rotationPhi(0)
28 {
29  setStripV(!m_isU);
30 }
31 
32 SVDRecoHit::SVDRecoHit(const SVDTrueHit* hit, bool uDirection, float sigma):
33  genfit::PlanarMeasurement(HIT_DIMENSIONS), m_sensorID(0), m_trueHit(hit),
34  m_cluster(0), m_isU(uDirection), m_energyDep(0), m_rotationPhi(0)
35 {
36  setStripV(!m_isU);
37 
38  // Smear the coordinate when constructing from a TrueHit.
39  if (!gRandom) B2FATAL("gRandom not initialized, please set up gRandom first");
40 
41  // Set the sensor UID
42  m_sensorID = hit->getSensorID();
43 
44  //If no error is given, estimate the error to pitch/sqrt(12)
45  if (sigma < 0) {
46  const SVD::SensorInfo& geometry = dynamic_cast<const SVD::SensorInfo&>(VXD::GeoCache::get(m_sensorID));
47  sigma = (m_isU) ? geometry.getUPitch(hit->getV()) / sqrt(12) : geometry.getVPitch() / sqrt(12);
48  }
49 
50  // Set positions
51  rawHitCoords_(0) = (m_isU) ? gRandom->Gaus(hit->getU(), sigma) : gRandom->Gaus(hit->getV(), sigma);
52  // Set the error covariance matrix
53  rawHitCov_(0, 0) = sigma * sigma;
54  // Set physical parameters
55  m_energyDep = hit->getEnergyDep();
56  // Setup geometry information
58 }
59 
61  genfit::PlanarMeasurement(HIT_DIMENSIONS), m_sensorID(0), m_trueHit(0),
62  m_cluster(hit), m_energyDep(0), m_rotationPhi(0)
63 {
64  // Set the sensor UID
65  m_sensorID = hit->getSensorID();
66  m_isU = hit->isUCluster();
67 
68  setStripV(!m_isU);
69 
70  // Determine if we have a wedge sensor.
71  const SVD::SensorInfo& geometry = dynamic_cast<const SVD::SensorInfo&>(VXD::GeoCache::get(m_sensorID));
72 
73  bool isWedgeU = m_isU && (geometry.getBackwardWidth() > geometry.getForwardWidth());
74 
75  // Set positions
76  rawHitCoords_(0) = hit->getPosition();
77  if (isWedgeU) {
78  // For u coordinate in a wedge sensor, the position line is not u = const.
79  // We have to rotate the coordinate system to achieve this.
80  m_rotationPhi = atan2((geometry.getBackwardWidth() - geometry.getForwardWidth()) / geometry.getWidth(0) * hit->getPosition(),
81  geometry.getLength());
82  }
83  // Set the error covariance matrix (this does not scale with position)
84  rawHitCov_(0, 0) = hit->getPositionSigma() * hit->getPositionSigma();
85  // Set physical parameters
86  m_energyDep = hit->getCharge();
87  // Setup geometry information
89 }
90 
92 {
93  // Construct a finite detector plane and set it.
94  const SVD::SensorInfo& geometry = dynamic_cast<const SVD::SensorInfo&>(VXD::GeoCache::get(m_sensorID));
95  bool isWedgeU = m_isU && (geometry.getBackwardWidth() > geometry.getForwardWidth());
96 
97  // Construct vectors o, u, v
98  TVector3 uLocal(1, 0, 0);
99  TVector3 vLocal(0, 1, 0);
100  TVector3 origin = geometry.pointToGlobal(TVector3(0, 0, 0), true);
101  TVector3 uGlobal = geometry.vectorToGlobal(uLocal, true);
102  TVector3 vGlobal = geometry.vectorToGlobal(vLocal, true);
103 
104  //Construct the detector plane
105  VXD::SensorPlane* finitePlane = new VXD::SensorPlane(m_sensorID, 20.0, 20.0);
106  if (isWedgeU) finitePlane->setRotation(m_rotationPhi);
107  genfit::SharedPlanePtr detPlane(new genfit::DetPlane(origin, uGlobal, vGlobal, finitePlane));
108  setPlane(detPlane, m_sensorID);
109 }
110 
112 {
113  return new SVDRecoHit(*this);
114 }
115 
116 TVectorD SVDRecoHit::applyPlanarDeformation(TVectorD rawHit, std::vector<double> planarParameters,
117  const genfit::StateOnPlane& state) const
118 {
119  // Legendre parametrization of deformation
120  auto L1 = [](double x) {return x;};
121  auto L2 = [](double x) {return (3 * x * x - 1) / 2;};
122  auto L3 = [](double x) {return (5 * x * x * x - 3 * x) / 2;};
123  auto L4 = [](double x) {return (35 * x * x * x * x - 30 * x * x + 3) / 8;};
124 
125  const SVD::SensorInfo& geometry = dynamic_cast<const SVD::SensorInfo&>(VXD::GeoCache::get(m_sensorID));
126 
127  double u = 0;
128  double v = 0;
129  double length = 0;
130  double width = 0;
131 
132  if (m_isU) {
133  u = rawHit[0]; // U coordinate of hit
134  v = state.getState()(4); // V coordinate of hit
135  width = geometry.getWidth(v); // Width of sensor (U side) is function of V (slanted)
136  length = geometry.getLength(); // Length of sensor (V side)
137 
138  u = u * 2 / width; // Legendre parametrization required U in (-1, 1)
139  v = v * 2 / length; // Legendre parametrization required V in (-1, 1)
140 
141  } else {
142  v = rawHit[0]; // V coordinate of hit
143  u = state.getState()(3); // U coordinate of hit
144  length = geometry.getLength(); // Length of sensor (V side) is fuction of V (slanted)
145  width = geometry.getWidth(v); // Width of sensor (U side)
146 
147  v = v * 2 / length; // Legendre parametrization required V in (-1, 1)
148  u = u * 2 / width; // Legendre parametrization required U in (-1, 1)
149 
150  }
151 
152  /* Planar deformation using Legendre parametrization
153  w(u, v) = L_{31} * L2(u) + L_{32} * L1(u) * L1(v) + L_{33} * L2(v) +
154  L_{41} * L3(u) + L_{42} * L2(u) * L1(v) + L_{43} * L1(u) * L2(v) + L_{44} * L3(v) +
155  L_{51} * L4(u) + L_{52} * L3(u) * L1(v) + L_{53} * L2(u) * L2(v) + L_{54} * L1(u) * L3(v) + L_{55} * L4(v); */
156  double dw =
157  planarParameters[0] * L2(u) + planarParameters[1] * L1(u) * L1(v) + planarParameters[2] * L2(v) +
158  planarParameters[3] * L3(u) + planarParameters[4] * L2(u) * L1(v) + planarParameters[5] * L1(u) * L2(v) + planarParameters[6] * L3(
159  v) +
160  planarParameters[7] * L4(u) + planarParameters[8] * L3(u) * L1(v) + planarParameters[9] * L2(u) * L2(v) + planarParameters[10] * L1(
161  u) * L3(v) + planarParameters[11] * L4(v);
162 
163  double du_dw = state.getState()[1]; // slope in U direction
164  double dv_dw = state.getState()[2]; // slope in V direction
165 
166  u = u * width / 2; // from Legendre to Local parametrization
167  v = v * length / 2; // from Legendre to Local parametrization
168 
169  TVectorD pos(1);
170 
171  if (m_isU) {
172  pos[0] = u + dw * du_dw;
173  } else {
174  pos[0] = v + dw * dv_dw;
175  }
176 
177  return pos;
178 }
179 
180 std::vector<genfit::MeasurementOnPlane*> SVDRecoHit::constructMeasurementsOnPlane(const genfit::StateOnPlane& state) const
181 {
182  if (!m_isU || m_rotationPhi == 0.0) {
183 
184  // Apply planar deformation to rectangular sensor or V coordinate of slanted sensor
185  TVectorD pos = applyPlanarDeformation(rawHitCoords_, VXD::GeoCache::get(m_sensorID).getSurfaceParameters(), state);
186 
187  return std::vector<genfit::MeasurementOnPlane*>(1, new genfit::MeasurementOnPlane(pos, rawHitCov_, state.getPlane(),
188  state.getRep(), this->constructHMatrix(state.getRep())));
189  }
190 
191  // Wedged sensor: the measured coordinate in U depends on V and the
192  // rotation angle. Namely, it needs to be scaled.
193  double u = rawHitCoords_(0);
194  double v = state.getState()(4);
195  double uPrime = u - v * tan(m_rotationPhi);
196  double scale = uPrime / u;
197 
198  TVectorD coords(1);
199  coords(0) = uPrime;
200 
201  // Apply planar deformation to U coordinate of slanted sensor
202  TVectorD pos = applyPlanarDeformation(coords, VXD::GeoCache::get(m_sensorID).getSurfaceParameters(), state);
203 
204  TMatrixDSym cov(scale * scale * rawHitCov_);
205 
206  return std::vector<genfit::MeasurementOnPlane*>(1, new genfit::MeasurementOnPlane(pos, cov, state.getPlane(), state.getRep(),
207  this->constructHMatrix(state.getRep())));
208 }
genfit::SharedPlanePtr
std::shared_ptr< genfit::DetPlane > SharedPlanePtr
Shared Pointer to a DetPlane.
Definition: SharedPlanePtr.h:40
Belle2::SVDRecoHit::applyPlanarDeformation
TVectorD applyPlanarDeformation(TVectorD rawHit, std::vector< double > planarParameters, const genfit::StateOnPlane &state) const
Apply planar deformation of sensors.
Definition: SVDRecoHit.cc:116
Belle2::SVDRecoHit::constructMeasurementsOnPlane
virtual std::vector< genfit::MeasurementOnPlane * > constructMeasurementsOnPlane(const genfit::StateOnPlane &state) const override
Methods that actually interface to Genfit.
Definition: SVDRecoHit.cc:180
Belle2::VXD::GeoCache::get
static const SensorInfoBase & get(Belle2::VxdID id)
Return a reference to the SensorInfo of a given SensorID.
Definition: GeoCache.h:141
Belle2::SVDRecoHit::m_energyDep
float m_energyDep
deposited energy.
Definition: SVDRecoHit.h:131
genfit::StateOnPlane
A state with arbitrary dimension defined in a DetPlane.
Definition: StateOnPlane.h:47
genfit
Defines for I/O streams used for error and debug printing.
Definition: AlignablePXDRecoHit.h:19
Belle2::SVDTrueHit
Class SVDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition: SVDTrueHit.h:35
genfit::TrackCandHit
Hit object for use in TrackCand.
Definition: TrackCandHit.h:34
Belle2::SVD::SensorInfo
Specific implementation of SensorInfo for SVD Sensors which provides additional sensor specific infor...
Definition: SensorInfo.h:35
Belle2::SVDRecoHit::m_sensorID
unsigned short m_sensorID
Unique sensor identifier.
Definition: SVDRecoHit.h:125
Belle2::VXD::SensorPlane::setRotation
void setRotation(double phi)
Set plane rotation angle.
Definition: SensorPlane.h:57
Belle2::SVDRecoHit::m_rotationPhi
float m_rotationPhi
angle of the plane rotation, for u in wedge sensors.
Definition: SVDRecoHit.h:133
genfit::AbsMeasurement
Contains the measurement and covariance in raw detector coordinates.
Definition: AbsMeasurement.h:42
genfit::PlanarMeasurement::setStripV
void setStripV(bool v=true)
Use if the coordinate for 1D hits measured in V direction.
Definition: PlanarMeasurement.h:70
Belle2::SVDRecoHit::m_isU
bool m_isU
transient member (not written out during streaming)
Definition: SVDRecoHit.h:130
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::SVDCluster
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:38
genfit::DetPlane
Detector plane.
Definition: DetPlane.h:59
Belle2::SVDRecoHit::SVDRecoHit
SVDRecoHit()
Default constructor for ROOT IO.
Definition: SVDRecoHit.cc:25
Belle2::SVDRecoHit::setDetectorPlane
void setDetectorPlane()
Set up Detector plane information.
Definition: SVDRecoHit.cc:91
genfit::MeasurementOnPlane
Measured coordinates on a plane.
Definition: MeasurementOnPlane.h:46
Belle2::VXD::SensorPlane
A Finite plane of one VXD Sensor.
Definition: SensorPlane.h:44
Belle2::SVDRecoHit::clone
genfit::AbsMeasurement * clone() const override
Creating a deep copy of this hit.
Definition: SVDRecoHit.cc:111