Belle II Software  release-05-01-25
PXDRecoHit.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  * Benjamin Schwenker *
8  * *
9  * This software is provided "as is" without any warranty. *
10  **************************************************************************/
11 
12 #include <framework/logging/Logger.h>
13 #include <pxd/reconstruction/PXDRecoHit.h>
14 #include <pxd/reconstruction/PXDClusterPositionEstimator.h>
15 #include <pxd/reconstruction/PXDGainCalibrator.h>
16 #include <pxd/dataobjects/PXDTrueHit.h>
17 #include <pxd/dataobjects/PXDCluster.h>
18 #include <pxd/geometry/SensorInfo.h>
19 #include <vxd/geometry/SensorPlane.h>
20 #include <vxd/geometry/GeoCache.h>
21 
22 #include <genfit/DetPlane.h>
23 #include <TVector3.h>
24 #include <TRandom.h>
25 
26 using namespace std;
27 using namespace Belle2;
28 
29 PXDRecoHit::PXDRecoHit():
30  genfit::PlanarMeasurement(HIT_DIMENSIONS), m_sensorID(0), m_trueHit(0), m_cluster(0),
31  m_energyDep(0)//, m_energyDepError(0)
32 {}
33 
34 PXDRecoHit::PXDRecoHit(const PXDTrueHit* hit, const genfit::TrackCandHit*, float sigmaU, float sigmaV):
35  genfit::PlanarMeasurement(HIT_DIMENSIONS), m_sensorID(0), m_trueHit(hit), m_cluster(0),
36  m_energyDep(0)//, m_energyDepError(0)
37 {
38  if (!gRandom) B2FATAL("gRandom not initialized, please set up gRandom first");
39 
40  // Set the sensor UID
41  m_sensorID = hit->getSensorID();
42 
43  //If no error is given, estimate the error by dividing the pixel size by sqrt(12)
44  if (sigmaU < 0 || sigmaV < 0) {
45  const PXD::SensorInfo& geometry = dynamic_cast<const PXD::SensorInfo&>(VXD::GeoCache::get(m_sensorID));
46  sigmaU = geometry.getUPitch(hit->getV()) / sqrt(12);
47  sigmaV = geometry.getVPitch(hit->getV()) / sqrt(12);
48  }
49 
50  // Set positions
51  rawHitCoords_(0) = gRandom->Gaus(hit->getU(), sigmaU);
52  rawHitCoords_(1) = gRandom->Gaus(hit->getV(), sigmaV);
53  // Set the error covariance matrix
54  rawHitCov_(0, 0) = sigmaU * sigmaU;
55  rawHitCov_(0, 1) = 0;
56  rawHitCov_(1, 0) = 0;
57  rawHitCov_(1, 1) = sigmaV * sigmaV;
58  // Set physical parameters
59  m_energyDep = hit->getEnergyDep();
60  // Setup geometry information
61  setDetectorPlane();
62 }
63 
64 PXDRecoHit::PXDRecoHit(const PXDCluster* hit, float sigmaU, float sigmaV, float covUV):
65  genfit::PlanarMeasurement(HIT_DIMENSIONS), m_sensorID(0), m_trueHit(0), m_cluster(hit),
66  m_energyDep(0)//, m_energyDepError(0)
67 {
68  // Set the sensor UID
69  m_sensorID = hit->getSensorID();
70  // Set positions
71  rawHitCoords_(0) = hit->getU();
72  rawHitCoords_(1) = hit->getV();
73  // Set the error covariance matrix
74  rawHitCov_(0, 0) = sigmaU * sigmaU;
75  rawHitCov_(0, 1) = covUV;
76  rawHitCov_(1, 0) = covUV;
77  rawHitCov_(1, 1) = sigmaV * sigmaV;
78  // Set physical parameters
79  const PXD::SensorInfo& SensorInfo = dynamic_cast<const PXD::SensorInfo&>(VXD::GeoCache::get(m_sensorID));
80  auto ADUToEnergy = PXD::PXDGainCalibrator::getInstance().getADUToEnergy(m_sensorID, SensorInfo.getUCellID(hit->getU()),
81  SensorInfo.getVCellID(hit->getV()));
82  m_energyDep = hit->getCharge() * ADUToEnergy;
83  //m_energyDepError = 0;
84  // Setup geometry information
85  setDetectorPlane();
86 }
87 
88 
89 PXDRecoHit::PXDRecoHit(const PXDCluster* hit, const genfit::TrackCandHit*):
90  genfit::PlanarMeasurement(HIT_DIMENSIONS), m_sensorID(0), m_trueHit(0), m_cluster(hit),
91  m_energyDep(0)//, m_energyDepError(0)
92 {
93  // Set the sensor UID
94  m_sensorID = hit->getSensorID();
95  // Set positions
96  rawHitCoords_(0) = hit->getU();
97  rawHitCoords_(1) = hit->getV();
98  // Set the error covariance matrix
99  rawHitCov_(0, 0) = hit->getUSigma() * hit->getUSigma();
100  rawHitCov_(0, 1) = hit->getRho() * hit->getUSigma() * hit->getVSigma();
101  rawHitCov_(1, 0) = hit->getRho() * hit->getUSigma() * hit->getVSigma();
102  rawHitCov_(1, 1) = hit->getVSigma() * hit->getVSigma();
103  // Set physical parameters
104  const PXD::SensorInfo& SensorInfo = dynamic_cast<const PXD::SensorInfo&>(VXD::GeoCache::get(m_sensorID));
105  auto ADUToEnergy = PXD::PXDGainCalibrator::getInstance().getADUToEnergy(m_sensorID, SensorInfo.getUCellID(hit->getU()),
106  SensorInfo.getVCellID(hit->getV()));
107  m_energyDep = hit->getCharge() * ADUToEnergy;
108  //m_energyDepError = 0;
109  // Setup geometry information
110  setDetectorPlane();
111 }
112 
113 genfit::AbsMeasurement* PXDRecoHit::clone() const
114 {
115  return new PXDRecoHit(*this);
116 }
117 
118 
119 void PXDRecoHit::setDetectorPlane()
120 {
121  // Construct a finite detector plane and set it.
122  const PXD::SensorInfo& geometry = dynamic_cast<const PXD::SensorInfo&>(VXD::GeoCache::get(m_sensorID));
123 
124  // Construct vectors o, u, v
125  TVector3 uLocal(1, 0, 0);
126  TVector3 vLocal(0, 1, 0);
127  TVector3 origin = geometry.pointToGlobal(TVector3(0, 0, 0), true);
128  TVector3 uGlobal = geometry.vectorToGlobal(uLocal, true);
129  TVector3 vGlobal = geometry.vectorToGlobal(vLocal, true);
130 
131  //Construct the detector plane
132  VXD::SensorPlane* finitePlane = new VXD::SensorPlane(m_sensorID, 20.0, 20.0);
133  genfit::SharedPlanePtr detPlane(new genfit::DetPlane(origin, uGlobal, vGlobal, finitePlane));
134  setPlane(detPlane, m_sensorID);
135 }
136 
137 
138 float PXDRecoHit::getShapeLikelyhood(const genfit::StateOnPlane& state) const
139 {
140  // We need an associated cluster
141  if (this->getCluster()) {
142  // Likelyhood depends on the fitted incidence angles into the sensor
143  const TVectorD& state5 = state.getState();
144  return PXD::PXDClusterPositionEstimator::getInstance().getShapeLikelyhood(*this->getCluster(), state5[1], state5[2]);
145  }
146  // If we reach here, we can do no better than return zero
147  return 0;
148 }
149 
150 TVectorD PXDRecoHit::applyPlanarDeformation(TVectorD hitCoords, std::vector<double> planarParameters,
151  const genfit::StateOnPlane& state) const
152 {
153  // Legendre parametrization of deformation
154  auto L1 = [](double x) {return x;};
155  auto L2 = [](double x) {return (3 * x * x - 1) / 2;};
156  auto L3 = [](double x) {return (5 * x * x * x - 3 * x) / 2;};
157  auto L4 = [](double x) {return (35 * x * x * x * x - 30 * x * x + 3) / 8;};
158 
159  const PXD::SensorInfo& geometry = dynamic_cast<const PXD::SensorInfo&>(VXD::GeoCache::get(m_sensorID));
160 
161  double u = hitCoords[0];
162  double v = hitCoords[1];
163  double width = geometry.getWidth(v); // Width of sensor (U side)
164  double length = geometry.getLength(); // Length of sensor (V side)
165  u = u * 2 / width; // Legendre parametrization required U in (-1, 1)
166  v = v * 2 / length; // Legendre parametrization required V in (-1, 1)
167 
168  /* Planar deformation using Legendre parametrization
169  w(u, v) = L_{31} * L2(u) + L_{32} * L1(u) * L1(v) + L_{33} * L2(v) +
170  L_{41} * L3(u) + L_{42} * L2(u) * L1(v) + L_{43} * L1(u) * L2(v) + L_{44} * L3(v) +
171  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); */
172  double dw =
173  planarParameters[0] * L2(u) + planarParameters[1] * L1(u) * L1(v) + planarParameters[2] * L2(v) +
174  planarParameters[3] * L3(u) + planarParameters[4] * L2(u) * L1(v) + planarParameters[5] * L1(u) * L2(v) + planarParameters[6] * L3(
175  v) +
176  planarParameters[7] * L4(u) + planarParameters[8] * L3(u) * L1(v) + planarParameters[9] * L2(u) * L2(v) + planarParameters[10] * L1(
177  u) * L3(v) + planarParameters[11] * L4(v);
178 
179  double du_dw = state.getState()[1]; // slope in U direction
180  double dv_dw = state.getState()[2]; // slope in V direction
181 
182  u = u * width / 2; // from Legendre to Local parametrization
183  v = v * length / 2; // from Legendre to Local parametrization
184 
185  TVectorD pos(2);
186 
187  pos[0] = u + dw * du_dw;
188  pos[1] = v + dw * dv_dw;
189 
190  return pos;
191 }
192 
193 std::vector<genfit::MeasurementOnPlane*> PXDRecoHit::constructMeasurementsOnPlane(const genfit::StateOnPlane& state) const
194 {
195  // Track-based update only takes place when the RecoHit has an associated cluster
196  if (this->getCluster()) {
197  // Check if we can correct position coordinates based on track info
198  const TVectorD& state5 = state.getState();
199  auto offset = PXD::PXDClusterPositionEstimator::getInstance().getClusterOffset(*this->getCluster(), state5[1], state5[2]);
200 
201  if (offset != nullptr) {
202  // Found a valid offset, lets apply it
203  const Belle2::VxdID& sensorID = (*this->getCluster()).getSensorID();
204  const Belle2::PXD::SensorInfo& Info = dynamic_cast<const Belle2::PXD::SensorInfo&>(VXD::GeoCache::get(sensorID));
205  double posU = Info.getUCellPosition((*this->getCluster()).getUStart());
206  double posV = Info.getVCellPosition((*this->getCluster()).getVStart());
207 
208  TVectorD hitCoords(2);
209  hitCoords(0) = posU + offset->getU();
210  hitCoords(1) = posV + offset->getV();
211  TMatrixDSym hitCov(2);
212  hitCov(0, 0) = offset->getUSigma2();
213  hitCov(0, 1) = offset->getUVCovariance();
214  hitCov(1, 0) = offset->getUVCovariance();
215  hitCov(1, 1) = offset->getVSigma2();
216 
217  // Apply planar deformation
218  TVectorD pos = applyPlanarDeformation(hitCoords, VXD::GeoCache::get(m_sensorID).getSurfaceParameters(), state);
219 
220  return std::vector<genfit::MeasurementOnPlane*>(1, new genfit::MeasurementOnPlane(
221  pos, hitCov, state.getPlane(), state.getRep(), this->constructHMatrix(state.getRep())
222  ));
223  }
224  }
225 
226  // Apply planar deformation
227  TVectorD pos = applyPlanarDeformation(rawHitCoords_, VXD::GeoCache::get(m_sensorID).getSurfaceParameters(), state);
228 
229  // If we reach here, we can do no better than what we have
230  return std::vector<genfit::MeasurementOnPlane*>(1, new genfit::MeasurementOnPlane(
231  pos, rawHitCov_, state.getPlane(), state.getRep(), this->constructHMatrix(state.getRep())
232  ));
233 }
234 
235 
236 
genfit::SharedPlanePtr
std::shared_ptr< genfit::DetPlane > SharedPlanePtr
Shared Pointer to a DetPlane.
Definition: SharedPlanePtr.h:40
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::VXD::SensorInfoBase::getUCellID
int getUCellID(double u, double v=0, bool clamp=false) const
Return the corresponding pixel/strip ID of a given u coordinate.
Definition: SensorInfoBase.h:203
Belle2::PXDTrueHit
Class PXDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition: PXDTrueHit.h:35
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
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
genfit::TrackCandHit
Hit object for use in TrackCand.
Definition: TrackCandHit.h:34
Belle2::PXDRecoHit::m_sensorID
unsigned short m_sensorID
Unique sensor identifier.
Definition: PXDRecoHit.h:145
Belle2::PXDRecoHit::clone
genfit::AbsMeasurement * clone() const override
Creating a deep copy of this hit.
Definition: PXDRecoHit.cc:113
Belle2::PXD::PXDClusterPositionEstimator::getInstance
static PXDClusterPositionEstimator & getInstance()
Main (and only) way to access the PXDClusterPositionEstimator.
Definition: PXDClusterPositionEstimator.cc:51
genfit::AbsMeasurement
Contains the measurement and covariance in raw detector coordinates.
Definition: AbsMeasurement.h:42
Belle2::PXDRecoHit::PXDRecoHit
PXDRecoHit()
Default constructor for ROOT IO.
Definition: PXDRecoHit.cc:29
Belle2::PXD::SensorInfo
Specific implementation of SensorInfo for PXD Sensors which provides additional pixel specific inform...
Definition: SensorInfo.h:34
Belle2::PXD::PXDGainCalibrator::getADUToEnergy
float getADUToEnergy(VxdID id, unsigned int uid, unsigned int vid) const
Get conversion factor from ADU to energy.
Definition: PXDGainCalibrator.cc:45
Belle2::VXD::SensorInfoBase::getVCellID
int getVCellID(double v, bool clamp=false) const
Return the corresponding pixel/strip ID of a given v coordinate.
Definition: SensorInfoBase.h:213
Belle2::VXD::SensorInfoBase::getVCellPosition
double getVCellPosition(int vID) const
Return the position of a specific strip/pixel in v direction.
Definition: SensorInfoBase.h:191
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::PXD::PXDClusterPositionEstimator::getShapeLikelyhood
float getShapeLikelyhood(const PXDCluster &cluster, double tu, double tv) const
Return cluster shape likelyhood.
Definition: PXDClusterPositionEstimator.cc:73
Belle2::PXD::PXDClusterPositionEstimator::getClusterOffset
const PXDClusterOffsetPar * getClusterOffset(const PXDCluster &cluster, double tu, double tv) const
Return pointer to cluster offsets, can be nullptr.
Definition: PXDClusterPositionEstimator.cc:58
Belle2::PXDRecoHit::setDetectorPlane
void setDetectorPlane()
Set up Detector plane information.
Definition: PXDRecoHit.cc:119
Belle2::PXDRecoHit::getSensorID
VxdID getSensorID() const
Get the compact ID.
Definition: PXDRecoHit.h:111
Belle2::PXDCluster
The PXD Cluster class This class stores all information about reconstructed PXD clusters The position...
Definition: PXDCluster.h:41
Belle2::PXDRecoHit::m_energyDep
float m_energyDep
transient member (not written out during streaming)
Definition: PXDRecoHit.h:150
Belle2::PXDRecoHit::getCluster
const PXDCluster * getCluster() const
Get pointer to the Cluster used when creating this RecoHit, can be NULL if created from something els...
Definition: PXDRecoHit.h:116
Belle2::PXDRecoHit::constructMeasurementsOnPlane
virtual std::vector< genfit::MeasurementOnPlane * > constructMeasurementsOnPlane(const genfit::StateOnPlane &state) const override
Methods that actually interface to Genfit.
Definition: PXDRecoHit.cc:193
Belle2::PXD::PXDGainCalibrator::getInstance
static PXDGainCalibrator & getInstance()
Main (and only) way to access the PXDGainCalibrator.
Definition: PXDGainCalibrator.cc:29
genfit::DetPlane
Detector plane.
Definition: DetPlane.h:59
Belle2::PXDRecoHit::getShapeLikelyhood
float getShapeLikelyhood(const genfit::StateOnPlane &state) const
Get deposited energy error.
Definition: PXDRecoHit.cc:138
Belle2::VXD::SensorInfoBase::getUCellPosition
double getUCellPosition(int uID, int vID=-1) const
Return the position of a specific strip/pixel in u direction.
Definition: SensorInfoBase.h:179
Belle2::PXDRecoHit::applyPlanarDeformation
TVectorD applyPlanarDeformation(TVectorD hitCoords, std::vector< double > planarParameters, const genfit::StateOnPlane &state) const
Apply planar deformation of sensors.
Definition: PXDRecoHit.cc:150
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