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 <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) {  // ipt = 1/pT, v = variance
0090     const float iipt2 = 1.0f / (ipt * ipt);                                 //iipt = 1/(1/pT) = pT
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) {  // ipt = 1/pT, v = variance
0098     const float iipt2 = 1.0f / (ipt * ipt);                                 //iipt = 1/(1/pT) = pT
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) {  // ipt = 1/pT, v = variance
0106     const float iipt2 = 1.0f / (ipt * ipt);                                     //iipt = 1/(1/pT) = pT
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     // Non-const versions needed for CopyOut of Matriplex.
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     // Indices for "fake" hit addition
0193     // Only if fake_hit_idx==-1, then count as missing hit for candidate score
0194     static constexpr int kHitMissIdx = -1;        //  hit is missed
0195     static constexpr int kHitStopIdx = -2;        //  track is stopped
0196     static constexpr int kHitEdgeIdx = -3;        //  track not in sensitive region of detector
0197     static constexpr int kHitMaxClusterIdx = -5;  //  hit cluster size > maxClusterSize
0198     static constexpr int kHitInGapIdx = -7;       //  track passing through inactive module
0199     static constexpr int kHitCCCFilterIdx = -9;   //  hit filtered via CCC (unused)
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;  // MIMI see set/get funcs; applicable for phase-0/1
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 }  // end namespace mkfit
0278 #endif