MkFitCore/src/TrackerInfo.cc

0001 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0002 #include "RecoTracker/MkFitCore/interface/PropagationConfig.h"
0003
0004 #include <cassert>
0005 #include <cstring>
0006
0007 namespace mkfit {
0008
0009   //==============================================================================
0010   // PropagationConfig
0011   //==============================================================================
0012
0013   void PropagationConfig::apply_tracker_info(const TrackerInfo* ti) {
0014     finding_inter_layer_pflags.tracker_info = ti;
0015     finding_intra_layer_pflags.tracker_info = ti;
0016     backward_fit_pflags.tracker_info = ti;
0017     forward_fit_pflags.tracker_info = ti;
0018     seed_fit_pflags.tracker_info = ti;
0019     pca_prop_pflags.tracker_info = ti;
0020   }
0021
0022   //==============================================================================
0023   // LayerInfo
0024   //==============================================================================
0025
0026   void LayerInfo::set_limits(float r1, float r2, float z1, float z2) {
0027     m_rin = r1;
0028     m_rout = r2;
0029     m_zmin = z1;
0030     m_zmax = z2;
0031   }
0032
0033   void LayerInfo::extend_limits(float r, float z) {
0034     if (z > m_zmax)
0035       m_zmax = z;
0036     if (z < m_zmin)
0037       m_zmin = z;
0038     if (r > m_rout)
0039       m_rout = r;
0040     if (r < m_rin)
0041       m_rin = r;
0042   }
0043
0044   void LayerInfo::set_r_in_out(float r1, float r2) {
0045     m_rin = r1;
0046     m_rout = r2;
0047   }
0048
0049   void LayerInfo::set_r_hole_range(float rh1, float rh2) {
0050     m_has_r_range_hole = true;
0051     m_hole_r_min = rh1;
0052     m_hole_r_max = rh2;
0053   }
0054
0055   void LayerInfo::print_layer() const {
0056     // clang-format off
0057     printf("Layer %2d  r(%7.4f, %7.4f) z(% 9.4f, % 9.4f) is_brl=%d, is_pix=%d, is_stereo=%d, q_bin=%.2f\n",
0058            m_layer_id,
0059            m_rin, m_rout, m_zmin, m_zmax,
0060            is_barrel(), m_is_pixel, m_is_stereo, m_q_bin);
0061     if (m_has_r_range_hole)
0062       printf("          has_r_range_hole: %.2f -> %.2f, dr: %f\n", m_hole_r_min, m_hole_r_max, m_hole_r_max - m_hole_r_min);
0063     // clang-format on
0064   }
0065
0066   //==============================================================================
0067   // TrackerInfo
0068   //==============================================================================
0069
0070   void TrackerInfo::reserve_layers(int n_brl, int n_ec_pos, int n_ec_neg) {
0071     m_layers.reserve(n_brl + n_ec_pos + n_ec_neg);
0072     m_barrel.reserve(n_brl);
0073     m_ecap_pos.reserve(n_ec_pos);
0074     m_ecap_neg.reserve(n_ec_neg);
0075   }
0076
0077   void TrackerInfo::create_layers(int n_brl, int n_ec_pos, int n_ec_neg) {
0078     reserve_layers(n_brl, n_ec_pos, n_ec_neg);
0079     for (int i = 0; i < n_brl; ++i)
0080       new_barrel_layer();
0081     for (int i = 0; i < n_ec_pos; ++i)
0082       new_ecap_pos_layer();
0083     for (int i = 0; i < n_ec_neg; ++i)
0084       new_ecap_neg_layer();
0085   }
0086
0087   int TrackerInfo::new_layer(LayerInfo::LayerType_e type) {
0088     int l = (int)m_layers.size();
0089     m_layers.emplace_back(LayerInfo(l, type));
0090     return l;
0091   }
0092
0093   LayerInfo& TrackerInfo::new_barrel_layer() {
0094     m_barrel.push_back(new_layer(LayerInfo::Barrel));
0095     return m_layers.back();
0096   }
0097
0098   LayerInfo& TrackerInfo::new_ecap_pos_layer() {
0099     m_ecap_pos.push_back(new_layer(LayerInfo::EndCapPos));
0100     return m_layers.back();
0101   }
0102
0103   LayerInfo& TrackerInfo::new_ecap_neg_layer() {
0104     m_ecap_neg.push_back(new_layer(LayerInfo::EndCapNeg));
0105     return m_layers.back();
0106   }
0107
0108   int TrackerInfo::n_total_modules() const {
0109     int nm = 0;
0110     for (auto& l : m_layers)
0111       nm += l.n_modules();
0112     return nm;
0113   }
0114
0115   //==============================================================================
0116   // Material
0117
0118   void TrackerInfo::create_material(int nBinZ, float rngZ, int nBinR, float rngR) {
0119     m_mat_range_z = rngZ;
0120     m_mat_range_r = rngR;
0121     m_mat_fac_z = nBinZ / m_mat_range_z;
0122     m_mat_fac_r = nBinR / m_mat_range_r;
0123
0124     m_mat_vec.rerect(nBinZ, nBinR);
0125   }
0126
0127   //==============================================================================
0128
0129   namespace {
0130     struct GeomFileHeader {
0131       int f_magic = s_magic;
0132       int f_format_version = s_version;
0133       int f_sizeof_trackerinfo = sizeof(TrackerInfo);
0134       int f_sizeof_layerinfo = sizeof(LayerInfo);
0135       int f_sizeof_moduleinfo = sizeof(ModuleInfo);
0136       int f_n_layers = -1;
0137
0138       GeomFileHeader() = default;
0139
0140       constexpr static int s_magic = 0xB00F;
0141       constexpr static int s_version = 2;
0142     };
0143
0144     template <typename T>
0145     int write_std_vec(FILE* fp, const std::vector<T>& vec, int el_size = 0) {
0146       int n = vec.size();
0147       fwrite(&n, sizeof(int), 1, fp);
0148       if (el_size == 0) {
0149         fwrite(vec.data(), sizeof(T), n, fp);
0150       } else {
0151         for (int i = 0; i < n; ++i)
0152           fwrite(&vec[i], el_size, 1, fp);
0153       }
0154       return n;
0155     }
0156
0157     template <typename T>
0158     int read_std_vec(FILE* fp, std::vector<T>& vec, int el_size = 0) {
0159       int n;
0160       fread(&n, sizeof(int), 1, fp);
0161       vec.resize(n);
0162       if (el_size == 0) {
0163         fread(vec.data(), sizeof(T), n, fp);
0164       } else {
0165         for (int i = 0; i < n; ++i)
0166           fread(&vec[i], el_size, 1, fp);
0167       }
0168       return n;
0169     }
0170   }  // namespace
0171
0172   void TrackerInfo::write_bin_file(const std::string& fname) const {
0173     FILE* fp = fopen(fname.c_str(), "w");
0174     if (!fp) {
0175       fprintf(stderr,
0176               "TrackerInfo::write_bin_file error opening file '%s', errno=%d: '%s'",
0177               fname.c_str(),
0178               errno,
0179               strerror(errno));
0180       throw std::runtime_error("Filed opening file in TrackerInfo::write_bin_file");
0181     }
0182     GeomFileHeader fh;
0183     fh.f_n_layers = n_layers();
0184     fwrite(&fh, sizeof(GeomFileHeader), 1, fp);
0185
0186     write_std_vec(fp, m_layers, (int)(offsetof(LayerInfo, m_is_pixel)) + 1);
0187     write_std_vec(fp, m_barrel);
0188     write_std_vec(fp, m_ecap_pos);
0189     write_std_vec(fp, m_ecap_neg);
0190
0191     for (int l = 0; l < fh.f_n_layers; ++l) {
0192       write_std_vec(fp, m_layers[l].m_modules);
0193     }
0194
0195     fwrite(&m_mat_range_z, 4 * sizeof(float), 1, fp);
0196     fwrite(&m_mat_vec, 2 * sizeof(int), 1, fp);
0197     write_std_vec(fp, m_mat_vec.vector());
0198
0199     fclose(fp);
0200   }
0201
0202   void TrackerInfo::read_bin_file(const std::string& fname) {
0203     FILE* fp = fopen(fname.c_str(), "r");
0204     if (!fp) {
0205       fprintf(stderr,
0206               "TrackerInfo::read_bin_file error opening file '%s', errno=%d: '%s'\n",
0207               fname.c_str(),
0208               errno,
0209               strerror(errno));
0210       throw std::runtime_error("Failed opening file in TrackerInfo::read_bin_file");
0211     }
0212     GeomFileHeader fh;
0213     fread(&fh, sizeof(GeomFileHeader), 1, fp);
0214
0215     if (fh.f_magic != GeomFileHeader::s_magic) {
0216       fprintf(stderr, "Incompatible input file (wrong magick).\n");
0217       throw std::runtime_error("Filed opening file in TrackerInfo::read_bin_file");
0218     }
0219     if (fh.f_format_version != GeomFileHeader::s_version) {
0220       fprintf(stderr,
0221               "Unsupported file version %d. Supported version is %d.\n",
0222               fh.f_format_version,
0223               GeomFileHeader::s_version);
0224       throw std::runtime_error("Unsupported file version in TrackerInfo::read_bin_file");
0225     }
0226     if (fh.f_sizeof_trackerinfo != sizeof(TrackerInfo)) {
0227       fprintf(stderr,
0228               "sizeof(TrackerInfo) on file (%d) different from current value (%d).\n",
0229               fh.f_sizeof_trackerinfo,
0230               (int)sizeof(TrackerInfo));
0231       throw std::runtime_error("sizeof(TrackerInfo) mismatch in TrackerInfo::read_bin_file");
0232     }
0233     if (fh.f_sizeof_layerinfo != sizeof(LayerInfo)) {
0234       fprintf(stderr,
0235               "sizeof(LayerInfo) on file (%d) different from current value (%d).\n",
0236               fh.f_sizeof_layerinfo,
0237               (int)sizeof(LayerInfo));
0238       throw std::runtime_error("sizeof(LayerInfo) mismatch in TrackerInfo::read_bin_file");
0239     }
0240     if (fh.f_sizeof_moduleinfo != sizeof(ModuleInfo)) {
0241       fprintf(stderr,
0242               "sizeof(ModuleInfo) on file (%d) different from current value (%d).\n",
0243               fh.f_sizeof_moduleinfo,
0244               (int)sizeof(ModuleInfo));
0245       throw std::runtime_error("sizeof(ModuleInfo) mismatch in TrackerInfo::read_bin_file");
0246     }
0247
0248     printf("Opened TrackerInfoGeom file '%s', format version %d, n_layers %d\n",
0249            fname.c_str(),
0250            fh.f_format_version,
0251            fh.f_n_layers);
0252
0253     read_std_vec(fp, m_layers, (int)(offsetof(LayerInfo, m_is_pixel)) + 1);
0254     read_std_vec(fp, m_barrel);
0255     read_std_vec(fp, m_ecap_pos);
0256     read_std_vec(fp, m_ecap_neg);
0257
0258     for (int l = 0; l < fh.f_n_layers; ++l) {
0259       LayerInfo& li = m_layers[l];
0260       int nm = read_std_vec(fp, li.m_modules);
0261
0262       li.m_detid2sid.clear();
0263       for (int m = 0; m < nm; ++m) {
0264         li.m_detid2sid.insert({li.m_modules[m].detid, m});
0265       }
0266     }
0267
0268     fread(&m_mat_range_z, 4 * sizeof(float), 1, fp);
0269     fread(&m_mat_vec, 2 * sizeof(int), 1, fp);
0270     read_std_vec(fp, m_mat_vec.vector());
0271
0272     fclose(fp);
0273   }
0274
0275   void TrackerInfo::print_tracker(int level) const {
0276     if (level > 0) {
0277       for (int i = 0; i < n_layers(); ++i) {
0278         const LayerInfo& li = layer(i);
0279         li.print_layer();
0280         if (level > 1) {
0281           printf("  Detailed module list N=%d\n", li.n_modules());
0282           for (int j = 0; j < li.n_modules(); ++j) {
0283             const ModuleInfo& mi = li.module_info(j);
0284             auto* p = mi.pos.Array();
0285             auto* z = mi.zdir.Array();
0286             auto* x = mi.xdir.Array();
0287             // clang-format off
0288             printf("Layer %d, mid=%u; detid=0x%x pos=%.3f,%.3f,%.3f, "
0289                   "norm=%.3f,%.3f,%.3f, phi=%.3f,%.3f,%.3f\n",
0290                   i, j, mi.detid, p[0], p[1], p[2],
0291                   z[0], z[1], z[2], x[0], x[1], x[2]);
0292             // clang-format on
0293           }
0294           printf("\n");
0295         }
0296       }
0297     }
0298   }
0299 }  // end namespace mkfit