File indexing completed on 2025-01-18 03:42:17
0001 #ifndef RecoTracker_MkFitCore_interface_Hit_h
0002 #define RecoTracker_MkFitCore_interface_Hit_h
0003
0004 #include "RecoTracker/MkFitCore/interface/Config.h"
0005 #include "RecoTracker/MkFitCore/interface/MatrixSTypes.h"
0006
0007 #include <cmath>
0008 #include <vdt/atan2.h>
0009 #include <vdt/cos.h>
0010 #include <vdt/log.h>
0011 #include <vdt/sin.h>
0012 #include <vdt/sincos.h>
0013 #include <vdt/tan.h>
0014 #include <vector>
0015 #include <string_view>
0016
0017 namespace mkfit {
0018
0019 template <typename T>
0020 inline T sqr(T x) {
0021 return x * x;
0022 }
0023 template <typename T>
0024 inline T cube(T x) {
0025 return x * x * x;
0026 }
0027
0028 inline float squashPhiGeneral(float phi) {
0029 return phi - floor(0.5 * Const::InvPI * (phi + Const::PI)) * Const::TwoPI;
0030 }
0031
0032 inline float squashPhiMinimal(float phi) {
0033 return phi >= Const::PI ? phi - Const::TwoPI : (phi < -Const::PI ? phi + Const::TwoPI : phi);
0034 }
0035
0036 inline float getRad2(float x, float y) { return x * x + y * y; }
0037
0038 inline float getInvRad2(float x, float y) { return 1.0f / (x * x + y * y); }
0039
0040 inline float getPhi(float x, float y) { return vdt::fast_atan2f(y, x); }
0041
0042 inline float getTheta(float r, float z) { return vdt::fast_atan2f(r, z); }
0043
0044 inline float getEta(float r, float z) { return -1.0f * vdt::fast_logf(vdt::fast_tanf(getTheta(r, z) / 2.0f)); }
0045
0046 inline float getEta(float theta) { return -1.0f * vdt::fast_logf(vdt::fast_tanf(theta / 2.0f)); }
0047
0048 inline float getEta(float x, float y, float z) {
0049 const float theta = vdt::fast_atan2f(std::sqrt(x * x + y * y), z);
0050 return -1.0f * vdt::fast_logf(vdt::fast_tanf(theta / 2.0f));
0051 }
0052
0053 inline float getHypot(float x, float y) { return std::sqrt(x * x + y * y); }
0054
0055 inline float getRadErr2(float x, float y, float exx, float eyy, float exy) {
0056 return (x * x * exx + y * y * eyy + 2.0f * x * y * exy) / getRad2(x, y);
0057 }
0058
0059 inline float getInvRadErr2(float x, float y, float exx, float eyy, float exy) {
0060 return (x * x * exx + y * y * eyy + 2.0f * x * y * exy) / cube(getRad2(x, y));
0061 }
0062
0063 inline float getPhiErr2(float x, float y, float exx, float eyy, float exy) {
0064 const float rad2 = getRad2(x, y);
0065 return (y * y * exx + x * x * eyy - 2.0f * x * y * exy) / (rad2 * rad2);
0066 }
0067
0068 inline float getThetaErr2(
0069 float x, float y, float z, float exx, float eyy, float ezz, float exy, float exz, float eyz) {
0070 const float rad2 = getRad2(x, y);
0071 const float rad = std::sqrt(rad2);
0072 const float hypot2 = rad2 + z * z;
0073 const float dthetadx = x * z / (rad * hypot2);
0074 const float dthetady = y * z / (rad * hypot2);
0075 const float dthetadz = -rad / hypot2;
0076 return dthetadx * dthetadx * exx + dthetady * dthetady * eyy + dthetadz * dthetadz * ezz +
0077 2.0f * dthetadx * dthetady * exy + 2.0f * dthetadx * dthetadz * exz + 2.0f * dthetady * dthetadz * eyz;
0078 }
0079
0080 inline float getEtaErr2(float x, float y, float z, float exx, float eyy, float ezz, float exy, float exz, float eyz) {
0081 const float rad2 = getRad2(x, y);
0082 const float detadx = -x / (rad2 * std::sqrt(1 + rad2 / (z * z)));
0083 const float detady = -y / (rad2 * std::sqrt(1 + rad2 / (z * z)));
0084 const float detadz = 1.0f / (z * std::sqrt(1 + rad2 / (z * z)));
0085 return detadx * detadx * exx + detady * detady * eyy + detadz * detadz * ezz + 2.0f * detadx * detady * exy +
0086 2.0f * detadx * detadz * exz + 2.0f * detady * detadz * eyz;
0087 }
0088
0089 inline float getPxPxErr2(float ipt, float phi, float vipt, float vphi) {
0090 const float iipt2 = 1.0f / (ipt * ipt);
0091 float cosP;
0092 float sinP;
0093 vdt::fast_sincosf(phi, sinP, cosP);
0094 return iipt2 * (iipt2 * cosP * cosP * vipt + sinP * sinP * vphi);
0095 }
0096
0097 inline float getPyPyErr2(float ipt, float phi, float vipt, float vphi) {
0098 const float iipt2 = 1.0f / (ipt * ipt);
0099 float cosP;
0100 float sinP;
0101 vdt::fast_sincosf(phi, sinP, cosP);
0102 return iipt2 * (iipt2 * sinP * sinP * vipt + cosP * cosP * vphi);
0103 }
0104
0105 inline float getPzPzErr2(float ipt, float theta, float vipt, float vtheta) {
0106 const float iipt2 = 1.0f / (ipt * ipt);
0107 const float cotT = 1.0f / vdt::fast_tanf(theta);
0108 const float cscT = 1.0f / vdt::fast_sinf(theta);
0109 return iipt2 * (iipt2 * cotT * cotT * vipt + cscT * cscT * cscT * cscT * vtheta);
0110 }
0111
0112 struct MCHitInfo {
0113 MCHitInfo() {}
0114 MCHitInfo(int track, int layer, int ithlayerhit, int mcHitID)
0115 : mcTrackID_(track), layer_(layer), ithLayerHit_(ithlayerhit), mcHitID_(mcHitID) {}
0116
0117 int mcTrackID_;
0118 int layer_;
0119 int ithLayerHit_;
0120 int mcHitID_;
0121
0122 int mcTrackID() const { return mcTrackID_; }
0123 int layer() const { return layer_; }
0124 int mcHitID() const { return mcHitID_; }
0125 static void reset();
0126 };
0127 typedef std::vector<MCHitInfo> MCHitInfoVec;
0128
0129 struct MeasurementState {
0130 public:
0131 MeasurementState() {}
0132 MeasurementState(const SVector3& p, const SMatrixSym33& e) : pos_(p), err_(e) {}
0133
0134 const SVector3& parameters() const { return pos_; }
0135 const SMatrixSym33& errors() const { return err_; }
0136
0137 SVector3 pos_;
0138 SMatrixSym33 err_;
0139 };
0140
0141 class Hit {
0142 public:
0143 Hit() : mcHitID_(-1) {}
0144
0145 Hit(const SVector3& position, const SMatrixSym33& error, int mcHitID = -1)
0146 : state_(position, error), mcHitID_(mcHitID) {}
0147
0148 const SVector3& position() const { return state_.parameters(); }
0149 const SVector3& parameters() const { return state_.parameters(); }
0150 const SMatrixSym33& error() const { return state_.err_; }
0151
0152 const float* posArray() const { return state_.pos_.Array(); }
0153 const float* errArray() const { return state_.err_.Array(); }
0154
0155
0156 SVector3& parameters_nc() { return state_.pos_; }
0157 SMatrixSym33& error_nc() { return state_.err_; }
0158
0159 float r() const {
0160 return sqrtf(state_.parameters().At(0) * state_.parameters().At(0) +
0161 state_.parameters().At(1) * state_.parameters().At(1));
0162 }
0163 float x() const { return state_.parameters().At(0); }
0164 float y() const { return state_.parameters().At(1); }
0165 float z() const { return state_.parameters().At(2); }
0166 float exx() const { return state_.errors().At(0, 0); }
0167 float eyy() const { return state_.errors().At(1, 1); }
0168 float ezz() const { return state_.errors().At(2, 2); }
0169 float phi() const { return getPhi(state_.parameters().At(0), state_.parameters().At(1)); }
0170 float eta() const {
0171 return getEta(state_.parameters().At(0), state_.parameters().At(1), state_.parameters().At(2));
0172 }
0173 float ephi() const { return getPhiErr2(x(), y(), exx(), eyy(), state_.errors().At(0, 1)); }
0174 float eeta() const {
0175 return getEtaErr2(x(),
0176 y(),
0177 z(),
0178 exx(),
0179 eyy(),
0180 ezz(),
0181 state_.errors().At(0, 1),
0182 state_.errors().At(0, 2),
0183 state_.errors().At(1, 2));
0184 }
0185
0186 const MeasurementState& measurementState() const { return state_; }
0187
0188 int mcHitID() const { return mcHitID_; }
0189 int layer(const MCHitInfoVec& globalMCHitInfo) const { return globalMCHitInfo[mcHitID_].layer(); }
0190 int mcTrackID(const MCHitInfoVec& globalMCHitInfo) const { return globalMCHitInfo[mcHitID_].mcTrackID(); }
0191
0192
0193
0194 static constexpr int kHitMissIdx = -1;
0195 static constexpr int kHitStopIdx = -2;
0196 static constexpr int kHitEdgeIdx = -3;
0197 static constexpr int kHitMaxClusterIdx = -5;
0198 static constexpr int kHitInGapIdx = -7;
0199 static constexpr int kHitCCCFilterIdx = -9;
0200
0201 static constexpr int kMinChargePerCM = 1620;
0202 static constexpr int kChargePerCMBits = 8;
0203 static constexpr int kDetIdInLayerBits = 14;
0204 static constexpr int kClusterSizeBits = 5;
0205 static constexpr int kMaxClusterSize = (1 << kClusterSizeBits) - 1;
0206
0207 struct PackedData {
0208 unsigned int detid_in_layer : kDetIdInLayerBits;
0209 unsigned int charge_pcm : kChargePerCMBits;
0210 unsigned int span_rows : kClusterSizeBits;
0211 unsigned int span_cols : kClusterSizeBits;
0212
0213 PackedData() : detid_in_layer(0), charge_pcm(0), span_rows(0), span_cols(0) {}
0214
0215 void set_charge_pcm(int cpcm) {
0216 if (cpcm < kMinChargePerCM)
0217 charge_pcm = 0;
0218 else
0219 charge_pcm = std::min((1 << kChargePerCMBits) - 1, ((cpcm - kMinChargePerCM) >> 3) + 1);
0220 }
0221 unsigned int get_charge_pcm() const {
0222 if (charge_pcm == 0)
0223 return 0;
0224 else
0225 return ((charge_pcm - 1) << 3) + kMinChargePerCM;
0226 }
0227 };
0228
0229 unsigned int detIDinLayer() const { return pdata_.detid_in_layer; }
0230 unsigned int chargePerCM() const { return pdata_.get_charge_pcm(); }
0231 unsigned int spanRows() const { return pdata_.span_rows + 1; }
0232 unsigned int spanCols() const { return pdata_.span_cols + 1; }
0233
0234 static unsigned int minChargePerCM() { return kMinChargePerCM; }
0235 static unsigned int maxChargePerCM() { return kMinChargePerCM + (((1 << kChargePerCMBits) - 2) << 3); }
0236 static unsigned int maxSpan() { return kMaxClusterSize; }
0237
0238 void setupAsPixel(unsigned int id, int rows, int cols) {
0239 pdata_.detid_in_layer = id;
0240 pdata_.charge_pcm = (1 << kChargePerCMBits) - 1;
0241 pdata_.span_rows = std::min(kMaxClusterSize, rows - 1);
0242 pdata_.span_cols = std::min(kMaxClusterSize, cols - 1);
0243 }
0244
0245 void setupAsStrip(unsigned int id, int cpcm, int rows) {
0246 pdata_.detid_in_layer = id;
0247 pdata_.set_charge_pcm(cpcm);
0248 pdata_.span_rows = std::min(kMaxClusterSize, rows - 1);
0249 }
0250
0251 private:
0252 MeasurementState state_;
0253 int mcHitID_;
0254 PackedData pdata_;
0255 };
0256
0257 typedef std::vector<Hit> HitVec;
0258
0259 struct HitOnTrack {
0260 int index : 24;
0261 int layer : 8;
0262
0263 HitOnTrack() : index(-1), layer(-1) {}
0264 HitOnTrack(int i, int l) : index(i), layer(l) {}
0265
0266 bool operator<(const HitOnTrack o) const {
0267 if (layer != o.layer)
0268 return layer < o.layer;
0269 return index < o.index;
0270 }
0271 };
0272
0273 typedef std::vector<HitOnTrack> HoTVec;
0274
0275 void print(std::string_view label, const MeasurementState& s);
0276
0277 }
0278 #endif