MkFitCore/interface/Hit.h

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 <vector>
0009 #include <string_view>
0010
0011 namespace mkfit {
0012
0013   template <typename T>
0014   inline T sqr(T x) {
0015     return x * x;
0016   }
0017   template <typename T>
0018   inline T cube(T x) {
0019     return x * x * x;
0020   }
0021
0022   inline float squashPhiGeneral(float phi) {
0023     return phi - floor(0.5 * Const::InvPI * (phi + Const::PI)) * Const::TwoPI;
0024   }
0025
0026   inline float squashPhiMinimal(float phi) {
0027     return phi >= Const::PI ? phi - Const::TwoPI : (phi < -Const::PI ? phi + Const::TwoPI : phi);
0028   }
0029
0030   inline float getRad2(float x, float y) { return x * x + y * y; }
0031
0032   inline float getInvRad2(float x, float y) { return 1.0f / (x * x + y * y); }
0033
0034   inline float getPhi(float x, float y) { return std::atan2(y, x); }
0035
0036   inline float getTheta(float r, float z) { return std::atan2(r, z); }
0037
0038   inline float getEta(float r, float z) { return -1.0f * std::log(std::tan(getTheta(r, z) / 2.0f)); }
0039
0040   inline float getEta(float theta) { return -1.0f * std::log(std::tan(theta / 2.0f)); }
0041
0042   inline float getEta(float x, float y, float z) {
0043     const float theta = std::atan2(std::sqrt(x * x + y * y), z);
0044     return -1.0f * std::log(std::tan(theta / 2.0f));
0045   }
0046
0047   inline float getHypot(float x, float y) { return std::sqrt(x * x + y * y); }
0048
0049   inline float getRadErr2(float x, float y, float exx, float eyy, float exy) {
0050     return (x * x * exx + y * y * eyy + 2.0f * x * y * exy) / getRad2(x, y);
0051   }
0052
0053   inline float getInvRadErr2(float x, float y, float exx, float eyy, float exy) {
0054     return (x * x * exx + y * y * eyy + 2.0f * x * y * exy) / cube(getRad2(x, y));
0055   }
0056
0057   inline float getPhiErr2(float x, float y, float exx, float eyy, float exy) {
0058     const float rad2 = getRad2(x, y);
0059     return (y * y * exx + x * x * eyy - 2.0f * x * y * exy) / (rad2 * rad2);
0060   }
0061
0062   inline float getThetaErr2(
0063       float x, float y, float z, float exx, float eyy, float ezz, float exy, float exz, float eyz) {
0064     const float rad2 = getRad2(x, y);
0065     const float rad = std::sqrt(rad2);
0066     const float hypot2 = rad2 + z * z;
0067     const float dthetadx = x * z / (rad * hypot2);
0068     const float dthetady = y * z / (rad * hypot2);
0069     const float dthetadz = -rad / hypot2;
0070     return dthetadx * dthetadx * exx + dthetady * dthetady * eyy + dthetadz * dthetadz * ezz +
0071            2.0f * dthetadx * dthetady * exy + 2.0f * dthetadx * dthetadz * exz + 2.0f * dthetady * dthetadz * eyz;
0072   }
0073
0074   inline float getEtaErr2(float x, float y, float z, float exx, float eyy, float ezz, float exy, float exz, float eyz) {
0075     const float rad2 = getRad2(x, y);
0076     const float detadx = -x / (rad2 * std::sqrt(1 + rad2 / (z * z)));
0077     const float detady = -y / (rad2 * std::sqrt(1 + rad2 / (z * z)));
0078     const float detadz = 1.0f / (z * std::sqrt(1 + rad2 / (z * z)));
0079     return detadx * detadx * exx + detady * detady * eyy + detadz * detadz * ezz + 2.0f * detadx * detady * exy +
0080            2.0f * detadx * detadz * exz + 2.0f * detady * detadz * eyz;
0081   }
0082
0083   inline float getPxPxErr2(float ipt, float phi, float vipt, float vphi) {  // ipt = 1/pT, v = variance
0084     const float iipt2 = 1.0f / (ipt * ipt);                                 //iipt = 1/(1/pT) = pT
0085     const float cosP = std::cos(phi);
0086     const float sinP = std::sin(phi);
0087     return iipt2 * (iipt2 * cosP * cosP * vipt + sinP * sinP * vphi);
0088   }
0089
0090   inline float getPyPyErr2(float ipt, float phi, float vipt, float vphi) {  // ipt = 1/pT, v = variance
0091     const float iipt2 = 1.0f / (ipt * ipt);                                 //iipt = 1/(1/pT) = pT
0092     const float cosP = std::cos(phi);
0093     const float sinP = std::sin(phi);
0094     return iipt2 * (iipt2 * sinP * sinP * vipt + cosP * cosP * vphi);
0095   }
0096
0097   inline float getPzPzErr2(float ipt, float theta, float vipt, float vtheta) {  // ipt = 1/pT, v = variance
0098     const float iipt2 = 1.0f / (ipt * ipt);                                     //iipt = 1/(1/pT) = pT
0099     const float cotT = 1.0f / std::tan(theta);
0100     const float cscT = 1.0f / std::sin(theta);
0101     return iipt2 * (iipt2 * cotT * cotT * vipt + cscT * cscT * cscT * cscT * vtheta);
0102   }
0103
0104   struct MCHitInfo {
0105     MCHitInfo() {}
0106     MCHitInfo(int track, int layer, int ithlayerhit, int mcHitID)
0107         : mcTrackID_(track), layer_(layer), ithLayerHit_(ithlayerhit), mcHitID_(mcHitID) {}
0108
0109     int mcTrackID_;
0110     int layer_;
0111     int ithLayerHit_;
0112     int mcHitID_;
0113
0114     int mcTrackID() const { return mcTrackID_; }
0115     int layer() const { return layer_; }
0116     int mcHitID() const { return mcHitID_; }
0117     static void reset();
0118   };
0119   typedef std::vector<MCHitInfo> MCHitInfoVec;
0120
0121   struct MeasurementState {
0122   public:
0123     MeasurementState() {}
0124     MeasurementState(const SVector3& p, const SVector6& e) : pos_(p), err_(e) {}
0125     MeasurementState(const SVector3& p, const SMatrixSym33& e) : pos_(p) {
0126       for (int i = 0; i < 6; ++i)
0127         err_[i] = e.Array()[i];
0128     }
0129     const SVector3& parameters() const { return pos_; }
0130     SMatrixSym33 errors() const {
0131       SMatrixSym33 result;
0132       for (int i = 0; i < 6; ++i)
0133         result.Array()[i] = err_[i];
0134       return result;
0135     }
0136     SVector3 pos_;
0137     SVector6 err_;
0138   };
0139
0140   class Hit {
0141   public:
0142     Hit() : mcHitID_(-1) {}
0143
0144     Hit(const SVector3& position, const SMatrixSym33& error, int mcHitID = -1)
0145         : state_(position, error), mcHitID_(mcHitID) {}
0146
0147     const SVector3& position() const { return state_.parameters(); }
0148     const SVector3& parameters() const { return state_.parameters(); }
0149     const SMatrixSym33 error() const { return state_.errors(); }
0150
0151     const float* posArray() const { return state_.pos_.Array(); }
0152     const float* errArray() const { return state_.err_.Array(); }
0153
0154     // Non-const versions needed for CopyOut of Matriplex.
0155     SVector3& parameters_nc() { return state_.pos_; }
0156     SVector6& error_nc() { return state_.err_; }
0157
0158     float r() const {
0159       return sqrtf(state_.parameters().At(0) * state_.parameters().At(0) +
0160                    state_.parameters().At(1) * state_.parameters().At(1));
0161     }
0162     float x() const { return state_.parameters().At(0); }
0163     float y() const { return state_.parameters().At(1); }
0164     float z() const { return state_.parameters().At(2); }
0165     float exx() const { return state_.errors().At(0, 0); }
0166     float eyy() const { return state_.errors().At(1, 1); }
0167     float ezz() const { return state_.errors().At(2, 2); }
0168     float phi() const { return getPhi(state_.parameters().At(0), state_.parameters().At(1)); }
0169     float eta() const {
0170       return getEta(state_.parameters().At(0), state_.parameters().At(1), state_.parameters().At(2));
0171     }
0172     float ephi() const { return getPhiErr2(x(), y(), exx(), eyy(), state_.errors().At(0, 1)); }
0173     float eeta() const {
0174       return getEtaErr2(x(),
0175                         y(),
0176                         z(),
0177                         exx(),
0178                         eyy(),
0179                         ezz(),
0180                         state_.errors().At(0, 1),
0181                         state_.errors().At(0, 2),
0182                         state_.errors().At(1, 2));
0183     }
0184
0185     const MeasurementState& measurementState() const { return state_; }
0186
0187     int mcHitID() const { return mcHitID_; }
0188     int layer(const MCHitInfoVec& globalMCHitInfo) const { return globalMCHitInfo[mcHitID_].layer(); }
0189     int mcTrackID(const MCHitInfoVec& globalMCHitInfo) const { return globalMCHitInfo[mcHitID_].mcTrackID(); }
0190
0191     // Indices for "fake" hit addition
0192     // Only if fake_hit_idx==-1, then count as missing hit for candidate score
0193     static constexpr int kHitMissIdx = -1;        //  hit is missed
0194     static constexpr int kHitStopIdx = -2;        //  track is stopped
0195     static constexpr int kHitEdgeIdx = -3;        //  track not in sensitive region of detector
0196     static constexpr int kHitMaxClusterIdx = -5;  //  hit cluster size > maxClusterSize
0197     static constexpr int kHitInGapIdx = -7;       //  track passing through inactive module
0198     static constexpr int kHitCCCFilterIdx = -9;   //  hit filtered via CCC (unused)
0199
0200     static constexpr int kMinChargePerCM = 1620;
0201     static constexpr int kChargePerCMBits = 8;
0202     static constexpr int kDetIdInLayerBits = 14;
0203     static constexpr int kClusterSizeBits = 5;
0204     static constexpr int kMaxClusterSize = (1 << kClusterSizeBits) - 1;
0205
0206     struct PackedData {
0207       unsigned int detid_in_layer : kDetIdInLayerBits;
0208       unsigned int charge_pcm : kChargePerCMBits;  // MIMI see set/get funcs; applicable for phase-0/1
0209       unsigned int span_rows : kClusterSizeBits;
0210       unsigned int span_cols : kClusterSizeBits;
0211
0212       PackedData() : detid_in_layer(0), charge_pcm(0), span_rows(0), span_cols(0) {}
0213
0214       void set_charge_pcm(int cpcm) {
0215         if (cpcm < kMinChargePerCM)
0216           charge_pcm = 0;
0217         else
0218           charge_pcm = std::min((1 << kChargePerCMBits) - 1, ((cpcm - kMinChargePerCM) >> 3) + 1);
0219       }
0220       unsigned int get_charge_pcm() const {
0221         if (charge_pcm == 0)
0222           return 0;
0223         else
0224           return ((charge_pcm - 1) << 3) + kMinChargePerCM;
0225       }
0226     };
0227
0228     unsigned int detIDinLayer() const { return pdata_.detid_in_layer; }
0229     unsigned int chargePerCM() const { return pdata_.get_charge_pcm(); }
0230     unsigned int spanRows() const { return pdata_.span_rows + 1; }
0231     unsigned int spanCols() const { return pdata_.span_cols + 1; }
0232
0233     static unsigned int minChargePerCM() { return kMinChargePerCM; }
0234     static unsigned int maxChargePerCM() { return kMinChargePerCM + (((1 << kChargePerCMBits) - 2) << 3); }
0235     static unsigned int maxSpan() { return kMaxClusterSize; }
0236
0237     void setupAsPixel(unsigned int id, int rows, int cols) {
0238       pdata_.detid_in_layer = id;
0239       pdata_.charge_pcm = (1 << kChargePerCMBits) - 1;
0240       pdata_.span_rows = std::min(kMaxClusterSize, rows - 1);
0241       pdata_.span_cols = std::min(kMaxClusterSize, cols - 1);
0242     }
0243
0244     void setupAsStrip(unsigned int id, int cpcm, int rows) {
0245       pdata_.detid_in_layer = id;
0246       pdata_.set_charge_pcm(cpcm);
0247       pdata_.span_rows = std::min(kMaxClusterSize, rows - 1);
0248     }
0249
0250   private:
0251     MeasurementState state_;
0252     int mcHitID_;
0253     PackedData pdata_;
0254   };
0255
0256   typedef std::vector<Hit> HitVec;
0257
0258   struct HitOnTrack {
0259     int index : 24;
0260     int layer : 8;
0261
0262     HitOnTrack() : index(-1), layer(-1) {}
0263     HitOnTrack(int i, int l) : index(i), layer(l) {}
0264
0265     bool operator<(const HitOnTrack o) const {
0266       if (layer != o.layer)
0267         return layer < o.layer;
0268       return index < o.index;
0269     }
0270   };
0271
0272   typedef std::vector<HitOnTrack> HoTVec;
0273
0274   void print(std::string_view label, const MeasurementState& s);
0275
0276   struct DeadRegion {
0277     float phi1, phi2, q1, q2;
0278     DeadRegion(float a1, float a2, float b1, float b2) : phi1(a1), phi2(a2), q1(b1), q2(b2) {}
0279   };
0280   typedef std::vector<DeadRegion> DeadVec;
0281
0282   struct BeamSpot {
0283     float x = 0, y = 0, z = 0;
0284     float sigmaZ = 5;
0285     float beamWidthX = 5e-4, beamWidthY = 5e-4;
0286     float dxdz = 0, dydz = 0;
0287
0288     BeamSpot() = default;
0289     BeamSpot(float ix, float iy, float iz, float is, float ibx, float iby, float idxdz, float idydz)
0290         : x(ix), y(iy), z(iz), sigmaZ(is), beamWidthX(ibx), beamWidthY(iby), dxdz(idxdz), dydz(idydz) {}
0291   };
0292 }  // end namespace mkfit
0293 #endif