MkFitCore/interface/TrackerInfo.h

0001 #ifndef RecoTracker_MkFitCore_interface_TrackerInfo_h
0002 #define RecoTracker_MkFitCore_interface_TrackerInfo_h
0003
0004 #include "RecoTracker/MkFitCore/interface/MatrixSTypes.h"
0005 #include "RecoTracker/MkFitCore/interface/PropagationConfig.h"
0006 #include "RecoTracker/MkFitCore/interface/Config.h"
0007 #include <string>
0008 #include <unordered_map>
0009 #include <vector>
0010 #include <unordered_map>
0011
0012 namespace mkfit {
0013
0014   //==============================================================================
0015   // WSR -- WithinSensitiveRegion state
0016
0017   enum WithinSensitiveRegion_e { WSR_Undef = -1, WSR_Inside = 0, WSR_Edge, WSR_Outside, WSR_Failed };
0018
0019   struct WSR_Result {
0020     // Could also store XHitSize count equivalent here : 16;
0021     WithinSensitiveRegion_e m_wsr : 8;
0022     bool m_in_gap : 8;
0023
0024     WSR_Result() : m_wsr(WSR_Undef), m_in_gap(false) {}
0025
0026     WSR_Result(WithinSensitiveRegion_e wsr, bool in_gap) : m_wsr(wsr), m_in_gap(in_gap) {}
0027   };
0028
0029   //==============================================================================
0030
0031   struct ModuleShape {
0032     float dx1;  // half-extent along x (bottom edge for trap)
0033     float dx2;  // 0 for rect; half-extent along x, top edge for trap
0034     float dy;   // half extent along y / less precise direction
0035     float dz;   // half thickness in z
0036
0037     // round to 1 mum precision
0038     float rmu(float x) { return std::round(1e4f * x) * 1e-4f; }
0039     void round_assign(float x1, float x2, float y, float z) {
0040       dx1 = rmu(x1);
0041       dx2 = rmu(x2);
0042       dy = rmu(y);
0043       dz = rmu(z);
0044     }
0045
0046     bool is_rect() const { return dx2 == 0.f; }
0047     bool is_trap() const { return dx2 != 0.f; }
0048
0049     bool operator==(const ModuleShape& s) const { return dx1 == s.dx1 && dx2 == s.dx2 && dy == s.dy && dz == s.dz; }
0050   };
0051
0052   struct ModuleInfo {
0053     SVector3 pos;
0054     SVector3 zdir;  // normal to module plane
0055     SVector3 xdir;  // the precise / "phi" direction
0056     unsigned int detid;
0057     unsigned short shapeid;
0058
0059     ModuleInfo() = default;
0060     ModuleInfo(SVector3 p, SVector3 zd, SVector3 xd, unsigned int did, unsigned short sid)
0061         : pos(p), zdir(zd), xdir(xd), detid(did), shapeid(sid) {}
0062
0063     SVector3 calc_ydir() const {
0064       return {zdir[1] * xdir[2] - zdir[2] * xdir[1],
0065               zdir[2] * xdir[0] - zdir[0] * xdir[2],
0066               zdir[0] * zdir[1] - zdir[1] * xdir[0]};
0067     }
0068   };
0069
0070   //==============================================================================
0071
0072   class LayerInfo {
0073     friend class TrackerInfo;
0074
0075   public:
0076     enum LayerType_e { Undef = -1, Barrel = 0, EndCapPos = 1, EndCapNeg = 2 };
0077
0078     LayerInfo() = default;
0079     LayerInfo(int lid, LayerType_e type) : m_layer_id(lid), m_layer_type(type) {}
0080
0081     void set_layer_type(LayerType_e t) { m_layer_type = t; }
0082     void set_limits(float r1, float r2, float z1, float z2);
0083     void extend_limits(float r, float z);
0084     void set_r_in_out(float r1, float r2);
0085     void set_propagate_to(float pto) { m_propagate_to = pto; }
0086     void set_r_hole_range(float rh1, float rh2);
0087     void set_q_bin(float qb) { m_q_bin = qb; }
0088     void set_subdet(int sd) { m_subdet = sd; }
0089     void set_is_pixel(bool p) { m_is_pixel = p; }
0090     void set_is_stereo(bool s) { m_is_stereo = s; }
0091     void set_has_charge(bool c) { m_has_charge = c; }
0092
0093     int layer_id() const { return m_layer_id; }
0094     LayerType_e layer_type() const { return m_layer_type; }
0095     float rin() const { return m_rin; }
0096     float rout() const { return m_rout; }
0097     float r_mean() const { return 0.5f * (m_rin + m_rout); }
0098     float zmin() const { return m_zmin; }
0099     float zmax() const { return m_zmax; }
0100     float z_mean() const { return 0.5f * (m_zmin + m_zmax); }
0101     float propagate_to() const { return m_propagate_to; }
0102     float q_bin() const { return m_q_bin; }
0103
0104     int subdet() const { return m_subdet; }
0105     bool is_barrel() const { return m_layer_type == Barrel; }
0106     bool is_pixel() const { return m_is_pixel; }
0107     bool is_stereo() const { return m_is_stereo; }
0108     bool has_charge() const { return m_has_charge; }
0109
0110     bool is_within_z_limits(float z) const { return z > m_zmin && z < m_zmax; }
0111     bool is_within_r_limits(float r) const { return r > m_rin && r < m_rout; }
0112     bool is_within_q_limits(float q) const { return is_barrel() ? is_within_z_limits(q) : is_within_r_limits(q); }
0113
0114     bool is_in_r_hole(float r) const { return m_has_r_range_hole ? is_in_r_hole_no_check(r) : false; }
0115
0116     WSR_Result is_within_z_sensitive_region(float z, float dz) const {
0117       if (z > m_zmax + dz || z < m_zmin - dz)
0118         return WSR_Result(WSR_Outside, false);
0119       if (z < m_zmax - dz && z > m_zmin + dz)
0120         return WSR_Result(WSR_Inside, false);
0121       return WSR_Result(WSR_Edge, false);
0122     }
0123
0124     WSR_Result is_within_r_sensitive_region(float r, float dr) const {
0125       if (r > m_rout + dr || r < m_rin - dr)
0126         return WSR_Result(WSR_Outside, false);
0127       if (r < m_rout - dr && r > m_rin + dr) {
0128         if (m_has_r_range_hole) {
0129           if (r < m_hole_r_max - dr && r > m_hole_r_min + dr)
0130             return WSR_Result(WSR_Outside, true);
0131           if (r < m_hole_r_max + dr && r > m_hole_r_min - dr)
0132             return WSR_Result(WSR_Edge, true);
0133         }
0134         return WSR_Result(WSR_Inside, false);
0135       }
0136       return WSR_Result(WSR_Edge, false);
0137     }
0138
0139     void print_layer() const;
0140
0141     // module & detid interface
0142     void reserve_modules(int nm) { m_modules.reserve(nm); }
0143     unsigned int register_module(ModuleInfo&& mi) {
0144       unsigned int pos = m_modules.size();
0145       m_modules.emplace_back(mi);
0146       m_detid2sid[mi.detid] = pos;
0147       return pos;
0148     }
0149     unsigned int shrink_modules() {
0150       m_modules.shrink_to_fit();
0151       return m_modules.size() - 1;
0152     }
0153
0154     void resize_shapes(int ns) { m_shapes.resize(ns); }
0155     void register_shape(const ModuleShape& ms, unsigned short sid) { m_shapes[sid] = ms; }
0156
0157     unsigned int short_id(unsigned int detid) const { return m_detid2sid.at(detid); }
0158     int n_modules() const { return m_modules.size(); }
0159     int n_shapes() const { return m_shapes.size(); }
0160     const ModuleInfo& module_info(unsigned int sid) const { return m_modules[sid]; }
0161     const ModuleShape& module_shape(unsigned short msid) const { return m_shapes[msid]; }
0162
0163   private:
0164     bool is_in_r_hole_no_check(float r) const { return r > m_hole_r_min && r < m_hole_r_max; }
0165
0166     int m_layer_id = -1;
0167     LayerType_e m_layer_type = Undef;
0168     int m_subdet = -1;  // sub-detector id, not used in core mkFit
0169
0170     float m_rin = 0, m_rout = 0, m_zmin = 0, m_zmax = 0;
0171     float m_propagate_to = 0;
0172
0173     float m_q_bin = 0;                         // > 0 - bin width, < 0 - number of bins
0174     float m_hole_r_min = 0, m_hole_r_max = 0;  // This could be turned into std::function when needed.
0175     bool m_has_r_range_hole = false;
0176     bool m_is_stereo = false;
0177     bool m_is_pixel = false;
0178     bool m_has_charge = true;
0179     // NOTE: offset of the last element is used in write/read_bin file.
0180     bool m_final_member_for_streaming = false;
0181
0182     std::unordered_map<unsigned int, unsigned int> m_detid2sid;
0183     std::vector<ModuleInfo> m_modules;
0184     std::vector<ModuleShape> m_shapes;
0185   };
0186
0187   //==============================================================================
0188
0189   template <typename T>
0190   class rectvec {
0191   public:
0192     rectvec(int n1 = 0, int n2 = 0) : m_n1(n1), m_n2(n2), m_vec(n1 * n2) {}
0193
0194     void rerect(int n1, int n2) {
0195       m_n1 = n1;
0196       m_n2 = n2;
0197       m_vec.resize(n1 * n2);
0198     }
0199
0200     const T& operator()(int i1, int i2) const { return m_vec[i1 * m_n2 + i2]; }
0201     T& operator()(int i1, int i2) { return m_vec[i1 * m_n2 + i2]; }
0202
0203     const T* operator[](int i1) const { return &m_vec[i1 * m_n2]; }
0204     T* operator[](int i1) { return &m_vec[i1 * m_n2]; }
0205
0206     const std::vector<T>& vector() const { return m_vec; }
0207     std::vector<T>& vector() { return m_vec; }
0208
0209     int n1() const { return m_n1; }
0210     int n2() const { return m_n2; }
0211     bool check_idcs(int i1, int i2) const { return i1 >= 0 && i1 < m_n1 && i2 >= 0 && i2 < m_n2; }
0212
0213   private:
0214     int m_n1, m_n2;
0215     std::vector<T> m_vec;
0216   };
0217
0218   class TrackerInfo {
0219   public:
0220     enum EtaRegion {
0221       Reg_Begin = 0,
0222       Reg_Endcap_Neg = 0,
0223       Reg_Transition_Neg,
0224       Reg_Barrel,
0225       Reg_Transition_Pos,
0226       Reg_Endcap_Pos,
0227       Reg_End,
0228       Reg_Count = Reg_End
0229     };
0230     struct Material {
0231       float bbxi{0}, radl{0};
0232     };
0233
0234     void reserve_layers(int n_brl, int n_ec_pos, int n_ec_neg);
0235     void create_layers(int n_brl, int n_ec_pos, int n_ec_neg);
0236     LayerInfo& new_barrel_layer();
0237     LayerInfo& new_ecap_pos_layer();
0238     LayerInfo& new_ecap_neg_layer();
0239
0240     int n_layers() const { return m_layers.size(); }
0241     const LayerInfo& layer(int l) const { return m_layers[l]; }
0242     LayerInfo& layer_nc(int l) { return m_layers[l]; }
0243
0244     int n_total_modules() const;
0245
0246     const LayerInfo& operator[](int l) const { return m_layers[l]; }
0247
0248     const LayerInfo& outer_barrel_layer() const { return m_layers[m_barrel.back()]; }
0249
0250     const std::vector<int>& barrel_layers() const { return m_barrel; }
0251     const std::vector<int>& endcap_pos_layers() const { return m_ecap_pos; }
0252     const std::vector<int>& endcap_neg_layers() const { return m_ecap_neg; }
0253
0254     const PropagationConfig& prop_config() const { return m_prop_config; }
0255     PropagationConfig& prop_config_nc() { return m_prop_config; }
0256
0257     void write_bin_file(const std::string& fname) const;
0258     void read_bin_file(const std::string& fname);
0259     void print_tracker(int level, int precision = 3) const;
0260
0261     void create_material(int nBinZ, float rngZ, int nBinR, float rngR);
0262     int mat_nbins_z() const { return m_mat_vec.n1(); }
0263     int mat_nbins_r() const { return m_mat_vec.n2(); }
0264     float mat_range_z() const { return m_mat_range_z; }
0265     float mat_range_r() const { return m_mat_range_r; }
0266     int mat_bin_z(float z) const { return z * m_mat_fac_z; }
0267     int mat_bin_r(float r) const { return r * m_mat_fac_r; }
0268     bool check_bins(int bz, int br) const { return m_mat_vec.check_idcs(bz, br); }
0269
0270     float material_bbxi(int binZ, int binR) const { return m_mat_vec(binZ, binR).bbxi; }
0271     float material_radl(int binZ, int binR) const { return m_mat_vec(binZ, binR).radl; }
0272     float& material_bbxi(int binZ, int binR) { return m_mat_vec(binZ, binR).bbxi; }
0273     float& material_radl(int binZ, int binR) { return m_mat_vec(binZ, binR).radl; }
0274
0275     Material material_checked(float z, float r) const {
0276       const int zbin = mat_bin_z(z), rbin = mat_bin_r(r);
0277       return check_bins(zbin, rbin) ? m_mat_vec(zbin, rbin) : Material();
0278     }
0279
0280   private:
0281     int new_layer(LayerInfo::LayerType_e type);
0282
0283     std::vector<LayerInfo> m_layers;
0284
0285     std::vector<int> m_barrel;
0286     std::vector<int> m_ecap_pos;
0287     std::vector<int> m_ecap_neg;
0288
0289     float m_mat_range_z, m_mat_range_r;
0290     float m_mat_fac_z, m_mat_fac_r;
0291     rectvec<Material> m_mat_vec;
0292
0293     PropagationConfig m_prop_config;
0294   };
0295
0296 }  // end namespace mkfit
0297 #endif