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)"
0058            " is_brl=%d, is_pix=%d, is_stereo=%d, has_charge=%d, q_bin=%.2f\n",
0059            m_layer_id, m_rin, m_rout, m_zmin, m_zmax,
0060            is_barrel(), m_is_pixel, m_is_stereo, m_has_charge, 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_sizeof_moduleshape = sizeof(ModuleShape);
0137       int f_n_layers = -1;
0138
0139       GeomFileHeader() = default;
0140
0141       constexpr static int s_magic = 0xB00F;
0142       constexpr static int s_version = 3;
0143     };
0144
0145     template <typename T>
0146     int write_std_vec(FILE* fp, const std::vector<T>& vec, int el_size = 0) {
0147       int n = vec.size();
0148       fwrite(&n, sizeof(int), 1, fp);
0149       if (el_size == 0) {
0150         fwrite(vec.data(), sizeof(T), n, fp);
0151       } else {
0152         for (int i = 0; i < n; ++i)
0153           fwrite(&vec[i], el_size, 1, fp);
0154       }
0155       return n;
0156     }
0157
0158     template <typename T>
0159     int read_std_vec(FILE* fp, std::vector<T>& vec, int el_size = 0) {
0160       int n;
0161       fread(&n, sizeof(int), 1, fp);
0162       vec.resize(n);
0163       if (el_size == 0) {
0164         fread(vec.data(), sizeof(T), n, fp);
0165       } else {
0166         for (int i = 0; i < n; ++i)
0167           fread(&vec[i], el_size, 1, fp);
0168       }
0169       return n;
0170     }
0171
0172     void assert_sizeof_match(int size_on_file, int size_of_class, const char* class_name) {
0173       if (size_on_file != size_of_class) {
0174         fprintf(stderr,
0175                 "sizeof(%s) on file (%d) different from current value (%d).\n",
0176                 class_name,
0177                 size_on_file,
0178                 size_of_class);
0179         throw std::runtime_error("class sizeof mismatch for " + std::string(class_name));
0180       }
0181     }
0182   }  // namespace
0183
0184   void TrackerInfo::write_bin_file(const std::string& fname) const {
0185     FILE* fp = fopen(fname.c_str(), "w");
0186     if (!fp) {
0187       fprintf(stderr,
0188               "TrackerInfo::write_bin_file error opening file '%s', errno=%d: '%s'",
0189               fname.c_str(),
0190               errno,
0191               strerror(errno));
0192       throw std::runtime_error("Filed opening file in TrackerInfo::write_bin_file");
0193     }
0194     GeomFileHeader fh;
0195     fh.f_n_layers = n_layers();
0196     fwrite(&fh, sizeof(GeomFileHeader), 1, fp);
0197
0198     write_std_vec(fp, m_layers, (int)offsetof(LayerInfo, m_final_member_for_streaming));
0199     write_std_vec(fp, m_barrel);
0200     write_std_vec(fp, m_ecap_pos);
0201     write_std_vec(fp, m_ecap_neg);
0202
0203     for (int l = 0; l < fh.f_n_layers; ++l) {
0204       write_std_vec(fp, m_layers[l].m_modules);
0205       write_std_vec(fp, m_layers[l].m_shapes);
0206     }
0207
0208     fwrite(&m_mat_range_z, 4 * sizeof(float), 1, fp);
0209     fwrite(&m_mat_vec, 2 * sizeof(int), 1, fp);
0210     write_std_vec(fp, m_mat_vec.vector());
0211
0212     fclose(fp);
0213   }
0214
0215   void TrackerInfo::read_bin_file(const std::string& fname) {
0216     FILE* fp = fopen(fname.c_str(), "r");
0217     if (!fp) {
0218       fprintf(stderr,
0219               "TrackerInfo::read_bin_file error opening file '%s', errno=%d: '%s'\n",
0220               fname.c_str(),
0221               errno,
0222               strerror(errno));
0223       throw std::runtime_error("Failed opening file in TrackerInfo::read_bin_file");
0224     }
0225     GeomFileHeader fh;
0226     fread(&fh, sizeof(GeomFileHeader), 1, fp);
0227
0228     if (fh.f_magic != GeomFileHeader::s_magic) {
0229       fprintf(stderr, "Incompatible input file (wrong magick).\n");
0230       throw std::runtime_error("Filed opening file in TrackerInfo::read_bin_file");
0231     }
0232     if (fh.f_format_version != GeomFileHeader::s_version) {
0233       fprintf(stderr,
0234               "Unsupported file version %d. Supported version is %d.\n",
0235               fh.f_format_version,
0236               GeomFileHeader::s_version);
0237       throw std::runtime_error("Unsupported file version in TrackerInfo::read_bin_file");
0238     }
0239     assert_sizeof_match(fh.f_sizeof_trackerinfo, sizeof(TrackerInfo), "TrackerInfo");
0240     assert_sizeof_match(fh.f_sizeof_layerinfo, sizeof(LayerInfo), "LayerInfo");
0241     assert_sizeof_match(fh.f_sizeof_moduleinfo, sizeof(ModuleInfo), "ModuleInfo");
0242     assert_sizeof_match(fh.f_sizeof_moduleshape, sizeof(ModuleShape), "ModuleShape");
0243
0244     printf("Opened TrackerInfoGeom file '%s', format version %d, n_layers %d\n",
0245            fname.c_str(),
0246            fh.f_format_version,
0247            fh.f_n_layers);
0248
0249     read_std_vec(fp, m_layers, (int)offsetof(LayerInfo, m_final_member_for_streaming));
0250     read_std_vec(fp, m_barrel);
0251     read_std_vec(fp, m_ecap_pos);
0252     read_std_vec(fp, m_ecap_neg);
0253
0254     for (int l = 0; l < fh.f_n_layers; ++l) {
0255       LayerInfo& li = m_layers[l];
0256       int nm = read_std_vec(fp, li.m_modules);
0257
0258       li.m_detid2sid.clear();
0259       for (int m = 0; m < nm; ++m) {
0260         li.m_detid2sid.insert({li.m_modules[m].detid, m});
0261       }
0262
0263       read_std_vec(fp, li.m_shapes);
0264     }
0265
0266     fread(&m_mat_range_z, 4 * sizeof(float), 1, fp);
0267     fread(&m_mat_vec, 2 * sizeof(int), 1, fp);
0268     read_std_vec(fp, m_mat_vec.vector());
0269
0270     fclose(fp);
0271   }
0272
0273   void TrackerInfo::print_tracker(int level, int precision) const {
0274     // clang-format off
0275     if (level > 0) {
0276       for (int i = 0; i < n_layers(); ++i) {
0277         const LayerInfo& li = layer(i);
0278         li.print_layer();
0279         if (level > 1) {
0280           printf("  Detailed module list N=%d, N_shapes=%d\n", li.n_modules(), li.n_shapes());
0281           for (int j = 0; j < li.n_shapes(); ++j) {
0282             const ModuleShape &ms = li.module_shape(j);
0283             const int w = precision + 1;
0284             printf("    Shape id=%u: dx1=%.*f, dx2=%.*f, dy=%.*f, dz=%.*f\n",
0285                    j, w, ms.dx1, w, ms.dx2, w, ms.dy, w, ms.dz);
0286           }
0287           for (int j = 0; j < li.n_modules(); ++j) {
0288             const ModuleInfo& mi = li.module_info(j);
0289             auto* p = mi.pos.Array();
0290             auto* z = mi.zdir.Array();
0291             auto* x = mi.xdir.Array();
0292             const int w = precision;
0293             printf("    Module id=%u: detid=0x%x shapeid=%u pos=%.*f,%.*f,%.*f, "
0294                    "norm=%.*f,%.*f,%.*f, phi=%.*f,%.*f,%.*f\n",
0295                    j, mi.detid, mi.shapeid, w, p[0], w, p[1], w, p[2],
0296                    w, z[0], w, z[1], w, z[2], w, x[0], w, x[1], w, x[2]);
0297           }
0298           printf("\n");
0299         }
0300       }
0301     }
0302     // clang-format on
0303   }
0304 }  // end namespace mkfit