Belle II Software  release-05-01-25
SensorInfo.cc
1 /*
2  * SensorInfo.cc
3  *
4  * Created on: Dec 19, 2013
5  * Author: kvasnicka
6  */
7 
8 #include <pxd/geometry/SensorInfo.h>
9 #include <framework/gearbox/Unit.h>
10 #include <framework/gearbox/Const.h>
11 #include <framework/geometry/BFieldManager.h>
12 #include <vxd/geometry/GeoCache.h>
13 
14 using namespace std;
15 using namespace Belle2;
16 using namespace Belle2::PXD;
17 
18 
19 double SensorInfo::getElectronMobility(double E) const
20 {
21  // Electron parameters - maximum velocity, critical intensity, beta factor
22  static double vmElec = 1.53 * pow(m_temperature, -0.87) * 1.E9 * Unit::cm
23  / Unit::s;
24  static double EcElec = 1.01 * pow(m_temperature, +1.55) * Unit::V
25  / Unit::cm;
26  static double betaElec = 2.57 * pow(m_temperature, +0.66) * 1.E-2;
27 
28  return (vmElec / EcElec * 1.
29  / (pow(1. + pow((fabs(E) / EcElec), betaElec), (1. / betaElec))));
30 }
31 
32 const TVector3 SensorInfo::getEField(const TVector3& point) const
33 {
34  // FIXME: Get rid of the gateDepth
35  double depletionVoltage = 0.5 * Unit::e * m_bulkDoping
36  / Const::permSi * m_thickness * m_thickness;
37  double gateZ = 0.5 * m_thickness - m_gateDepth;
38  double Ez = 2 * depletionVoltage * (point.Z() - gateZ) / m_thickness
39  / m_thickness;
40  TVector3 E(0, 0, Ez);
41  return E;
42 }
43 
44 const TVector3 SensorInfo::getBField(const TVector3& point) const
45 {
46  TVector3 pointGlobal = pointToGlobal(point, true);
47  TVector3 bGlobal = BFieldManager::getField(pointGlobal);
48  TVector3 bLocal = vectorToLocal(bGlobal, true);
49  return bLocal;
50 }
51 
52 const TVector3 SensorInfo::getDriftVelocity(const TVector3& E,
53  const TVector3& B) const
54 {
55  // Set mobility parameters
56  double mobility = -getElectronMobility(E.Mag());
57  double mobilityH = m_hallFactor * mobility;
58  // Calculate products
59  TVector3 EcrossB = E.Cross(B);
60  TVector3 BEdotB = E.Dot(B) * B;
61  TVector3 v = mobility * E + mobility * mobilityH * EcrossB
62  + mobility * mobilityH * mobilityH * BEdotB;
63  v *= 1.0 / (1.0 + mobilityH * mobilityH * B.Mag2());
64  return v;
65 }
66 
67 int SensorInfo::getPixelKind(const VxdID sensorID, double v) const
68 {
69  const SensorInfo& Info = dynamic_cast<const SensorInfo&>(VXD::GeoCache::get(sensorID));
70  int i_pixelKind = Info.getVPitchID(v);
71  if (Info.getID().getLayerNumber() == 2) i_pixelKind += 4;
72  if (Info.getID().getSensorNumber() == 2) i_pixelKind += 2;
73  return i_pixelKind;
74 }
75 
76 int SensorInfo::getPixelKindNew(const VxdID& sensorID, int vID) const
77 {
78  const SensorInfo& Info = dynamic_cast<const SensorInfo&>(VXD::GeoCache::get(sensorID));
79  double v = Info.getVCellPosition(vID);
80  double vPitch = Info.getVPitch(v);
81  int i_pixelKind = 0;
82 
83  if (std::fabs(vPitch - 0.0055) < 0.0001)
84  i_pixelKind = 0;
85  else if (std::fabs(vPitch - 0.0060) < 0.0001)
86  i_pixelKind = 1;
87  else if (std::fabs(vPitch - 0.0070) < 0.0001)
88  i_pixelKind = 2;
89  else if (std::fabs(vPitch - 0.0085) < 0.0001)
90  i_pixelKind = 3;
91  else
92  B2FATAL("Unexpected pixel vPitch.");
93  return i_pixelKind;
94 }
95 
96 
97 const TVector3 SensorInfo::getLorentzShift(double u, double v) const
98 {
99  // Constants for the 5-point Gauss quadrature formula
100  const double alphaGL = 1.0 / 3.0 * sqrt(5.0 + 2.0 * sqrt(10.0 / 7.0));
101  const double betaGL = 1.0 / 3.0 * sqrt(5.0 - 2.0 * sqrt(10.0 / 7.0));
102  const double walpha = (322 - 13.0 * sqrt(70)) / 900;
103  const double wbeta = (322 + 13.0 * sqrt(70)) / 900;
104  const double distanceToPlane = 0.5 * m_thickness - m_gateDepth;
105  const double midpoint = 0.5 * distanceToPlane;
106  const double h = 0.5 * distanceToPlane;
107  const double weightGL[5] = {
108  h * walpha, h * wbeta, h * 128.0 / 225.0, h * wbeta, h* walpha
109  };
110  const double zKnots[5] = {
111  midpoint - alphaGL * h, midpoint - betaGL * h, midpoint, midpoint + betaGL * h, midpoint + alphaGL* h
112  };
113  // Integrate v/v_z from z = 0 to z = distanceToPlane
114  TVector3 position(u, v, 0);
115  TVector3 currentBField = getBField(position);
116  for (int iz = 0; iz < 5; ++iz) {
117  // This is OK as long as E only depends on z
118  TVector3 currentEField = getEField(TVector3(0, 0, zKnots[iz]));
119  TVector3 velocity = getDriftVelocity(currentEField, currentBField);
120  position += weightGL[iz] / velocity.Z() * velocity;
121  } // for knots
122  position.SetZ(0);
123  position.SetX(position.X() - u);
124  position.SetY(position.Y() - v);
125  return position;
126 }
127 
128 void SensorInfo::cook()
129 {
130  m_iup = m_uCells / m_width;
131  m_up = m_width / m_uCells;
132 
133  m_vsplit = m_length * (m_splitLength - 0.5);
134 
135  m_vp = m_length * m_splitLength / m_vCells;
136  m_ivp = 1 / m_vp;
137 
138  m_vp2 = m_length * (1 - m_splitLength) / m_vCells2;
139  m_ivp2 = 1 / m_vp2;
140 
141  m_hxIG = 0.5 * m_up - m_clearBorderSmallPitch;
142  m_mIGL = 0.5 * (m_vp - m_sourceBorderLargePitch + m_drainBorderLargePitch);
143  m_sIGL = 0.5 * (m_vp - m_sourceBorderLargePitch - m_drainBorderLargePitch);
144  m_mIGS = 0.5 * (m_vp2 - m_sourceBorderSmallPitch + m_drainBorderSmallPitch);
145  m_sIGS = 0.5 * (m_vp2 - m_sourceBorderSmallPitch - m_drainBorderSmallPitch);
146 }
147 
148 int SensorInfo::getTrappedID(double x, double y) const
149 {
150  double huCells = 0.5 * m_uCells;
151  double ix = floor(x * m_iup + huCells);
152  int jx = ix;
153  double x0 = (ix + 0.5 - huCells) * m_up;
154 
155  if (fabs(x - x0) < m_hxIG) {
156  if ((unsigned)jx >= (unsigned)m_uCells) return -1;
157  double ys = y - m_vsplit;
158  if (ys >= 0) {
159  double iy = floor(ys * m_ivp2);
160  int jy = iy;
161  iy = jy / 2;
162  double y0 = iy * m_vp2 * 2 + m_vp2;
163  double yl = fabs(ys - y0);
164  if (fabs(yl - m_mIGS) < m_sIGS) {
165  if ((unsigned)jy >= (unsigned)m_vCells2) return -1;
166  return jx + m_uCells * (jy + m_vCells);
167  }
168  return -1;
169  } else {
170  ys = y + 0.5 * m_length;
171  double iy = floor(ys * m_ivp);
172  int jy = iy;
173  iy = jy / 2;
174  double y0 = iy * m_vp * 2 + m_vp;
175  double yl = fabs(ys - y0);
176  if (fabs(yl - m_mIGL) < m_sIGL) {
177  if ((unsigned)jy >= (unsigned)m_vCells) return -1;
178  return jx + m_uCells * jy;
179  }
180  return -1;
181  }
182  }
183  return -1;
184 }
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::VXD::SensorInfoBase::getID
VxdID getID() const
Return the ID of the Sensor.
Definition: SensorInfoBase.h:85
Belle2::VXD::SensorInfoBase::getVPitchID
int getVPitchID(double v=0) const
Return the pitch ID of the sensor.
Definition: SensorInfoBase.h:167
Belle2::VXD::SensorInfoBase::getVPitch
double getVPitch(double v=0) const
Return the pitch of the sensor.
Definition: SensorInfoBase.h:150
Belle2::PXD::SensorInfo
Specific implementation of SensorInfo for PXD Sensors which provides additional pixel specific inform...
Definition: SensorInfo.h:34
Belle2::PXD
Namespace to encapsulate code needed for simulation and reconstrucion of the PXD.
Definition: PXDCalibrationUtilities.h:28
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::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:110
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:106