DTTriggerPhase2/src/MuonPathSLFitter.cc

0001 #include "L1Trigger/DTTriggerPhase2/interface/MuonPathSLFitter.h"
0002 #include <cmath>
0003 #include <memory>
0004
0005 using namespace edm;
0006 using namespace std;
0007 using namespace cmsdt;
0008 // ============================================================================
0009 // Constructors and destructor
0010 // ============================================================================
0011 MuonPathSLFitter::MuonPathSLFitter(const ParameterSet &pset,
0012                                    edm::ConsumesCollector &iC,
0013                                    std::shared_ptr<GlobalCoordsObtainer> &globalcoordsobtainer)
0014     : MuonPathFitter(pset, iC, globalcoordsobtainer) {
0015   if (debug_)
0016     LogDebug("MuonPathSLFitter") << "MuonPathSLFitter: constructor";
0017
0018   //shift theta
0019   int rawId;
0020   double shift;
0021   shift_theta_filename_ = pset.getParameter<edm::FileInPath>("shift_theta_filename");
0022   std::ifstream ifin4(shift_theta_filename_.fullPath());
0023   if (ifin4.fail()) {
0024     throw cms::Exception("Missing Input File")
0025         << "MuonPathSLFitter::MuonPathSLFitter() -  Cannot find " << shift_theta_filename_.fullPath();
0026   }
0027
0028   while (ifin4.good()) {
0029     ifin4 >> rawId >> shift;
0030     shiftthetainfo_[rawId] = shift;
0031   }
0032
0033   // LUTs
0034   sl1_filename_ = pset.getParameter<edm::FileInPath>("lut_sl1");
0035   sl2_filename_ = pset.getParameter<edm::FileInPath>("lut_sl2");
0036   sl3_filename_ = pset.getParameter<edm::FileInPath>("lut_sl3");
0037
0038   fillLuts();
0039
0040   setChi2Th(pset.getParameter<double>("chi2Th"));
0041   setTanPhiTh(pset.getParameter<double>("tanPhiTh"));
0042 }
0043
0044 MuonPathSLFitter::~MuonPathSLFitter() {
0045   if (debug_)
0046     LogDebug("MuonPathSLFitter") << "MuonPathSLFitter: destructor";
0047 }
0048
0049 // ============================================================================
0050 // Main methods (initialise, run, finish)
0051 // ============================================================================
0052 void MuonPathSLFitter::initialise(const edm::EventSetup &iEventSetup) {
0053   if (debug_)
0054     LogDebug("MuonPathSLFitter") << "MuonPathSLFitter::initialiase";
0055
0056   auto geom = iEventSetup.getHandle(dtGeomH);
0057   dtGeo_ = &(*geom);
0058 }
0059
0060 void MuonPathSLFitter::run(edm::Event &iEvent,
0061                            const edm::EventSetup &iEventSetup,
0062                            MuonPathPtrs &muonpaths,
0063                            std::vector<lat_vector> &lateralities,
0064                            std::vector<metaPrimitive> &metaPrimitives) {
0065   if (debug_)
0066     LogDebug("MuonPathSLFitter") << "MuonPathSLFitter: run";
0067
0068   // fit per SL (need to allow for multiple outputs for a single mpath)
0069   // for (auto &muonpath : muonpaths) {
0070   if (!muonpaths.empty()) {
0071     auto muonpath = muonpaths[0];
0072     int rawId = muonpath->primitive(0)->cameraId();
0073     if (muonpath->primitive(0)->cameraId() == -1) {
0074       rawId = muonpath->primitive(1)->cameraId();
0075     }
0076     const DTLayerId dtLId(rawId);
0077     max_drift_tdc = maxdriftinfo_[dtLId.wheel() + 2][dtLId.station() - 1][dtLId.sector() - 1];
0078   }
0079
0080   for (size_t i = 0; i < muonpaths.size(); i++) {
0081     auto muonpath = muonpaths[i];
0082     auto lats = lateralities[i];
0083     analyze(muonpath, lats, metaPrimitives);
0084   }
0085 }
0086
0087 void MuonPathSLFitter::finish() {
0088   if (debug_)
0089     LogDebug("MuonPathSLFitter") << "MuonPathSLFitter: finish";
0090 };
0091
0092 //------------------------------------------------------------------
0093 //--- Metodos privados
0094 //------------------------------------------------------------------
0095
0096 void MuonPathSLFitter::analyze(MuonPathPtr &inMPath,
0097                                lat_vector lat_combs,
0098                                std::vector<cmsdt::metaPrimitive> &metaPrimitives) {
0099   auto sl = inMPath->primitive(0)->superLayerId();  // 0, 1, 2
0100
0101   int selected_lay = 1;
0102   if (inMPath->primitive(0)->cameraId() > 0)
0103     selected_lay = 0;
0104
0105   int dumLayId = inMPath->primitive(selected_lay)->cameraId();
0106   auto dtDumlayerId = DTLayerId(dumLayId);
0107   DTSuperLayerId MuonPathSLId(dtDumlayerId.wheel(), dtDumlayerId.station(), dtDumlayerId.sector(), sl + 1);
0108
0109   DTChamberId ChId(MuonPathSLId.wheel(), MuonPathSLId.station(), MuonPathSLId.sector());
0110
0111   DTSuperLayerId MuonPathSL1Id(dtDumlayerId.wheel(), dtDumlayerId.station(), dtDumlayerId.sector(), 1);
0112   DTSuperLayerId MuonPathSL2Id(dtDumlayerId.wheel(), dtDumlayerId.station(), dtDumlayerId.sector(), 2);
0113   DTSuperLayerId MuonPathSL3Id(dtDumlayerId.wheel(), dtDumlayerId.station(), dtDumlayerId.sector(), 3);
0114   DTWireId wireIdSL1(MuonPathSL1Id, 2, 1);
0115   DTWireId wireIdSL2(MuonPathSL2Id, 2, 1);
0116   DTWireId wireIdSL3(MuonPathSL3Id, 2, 1);
0117   auto sl_shift_cm = shiftinfo_[wireIdSL1.rawId()] - shiftinfo_[wireIdSL3.rawId()];
0118
0119   fit_common_in_t fit_common_in;
0120
0121   // 8-element vectors, for the 8 layers. As here we are fitting one SL only, we leave the other SL values as dummy ones
0122   fit_common_in.hits = {};
0123   fit_common_in.hits_valid = {};
0124
0125   int quality = 3;
0126   if (inMPath->missingLayer() != -1)
0127     quality = 1;
0128
0129   int minISL = 1;
0130   int maxISL = 3;
0131   if (sl < 1) {
0132     minISL = 0;
0133     maxISL = 2;
0134   }
0135
0136   for (int isl = minISL; isl < maxISL; isl++) {
0137     for (int i = 0; i < NUM_LAYERS; i++) {
0138       if (isl == sl && inMPath->missingLayer() != i) {
0139         // Include both valid and non-valid hits. Non-valid values can be whatever, leaving all as -1 to make debugging easier.
0140         auto ti = inMPath->primitive(i)->tdcTimeStamp();
0141         if (ti != -1)
0142           ti = (int)round(((float)TIME_TO_TDC_COUNTS / (float)LHC_CLK_FREQ) * ti);
0143         auto wi = inMPath->primitive(i)->channelId();
0144         auto ly = inMPath->primitive(i)->layerId();
0145         // int layId = inMPath->primitive(i)->cameraId();
0146         // auto dtlayerId = DTLayerId(layId);
0147         // auto wireId = DTWireId(dtlayerId, wi + 1); // wire start from 1, mixer groups them starting from 0
0148         // int rawId = wireId.rawId();
0149         // wp in tdc counts (still in floating point)
0150         int wp_semicells = (wi - SL1_CELLS_OFFSET) * 2 + 1;
0151         if (ly % 2 == 1)
0152           wp_semicells -= 1;
0153         if (isl == 2)
0154           wp_semicells -= (int)round((sl_shift_cm * 10) / CELL_SEMILENGTH);
0155
0156         float wp_tdc = wp_semicells * max_drift_tdc;
0157         int wp = (int)((long int)(round(wp_tdc * std::pow(2, WIREPOS_WIDTH))) / (int)std::pow(2, WIREPOS_WIDTH));
0158         fit_common_in.hits.push_back({ti, wi, ly, wp});
0159         // fill valids as well
0160         if (inMPath->missingLayer() == i)
0161           fit_common_in.hits_valid.push_back(0);
0162         else
0163           fit_common_in.hits_valid.push_back(1);
0164       } else {
0165         fit_common_in.hits.push_back({-1, -1, -1, -1});
0166         fit_common_in.hits_valid.push_back(0);
0167       }
0168     }
0169   }
0170
0171   int smallest_time = 999999, tmp_coarse_wirepos_1 = -1, tmp_coarse_wirepos_3 = -1;
0172   // coarse_bctr is the 12 MSB of the smallest tdc
0173   for (int isl = 0; isl < 3; isl++) {
0174     if (isl != sl)
0175       continue;
0176     int myisl = isl < 2 ? 0 : 1;
0177     for (size_t i = 0; i < NUM_LAYERS; i++) {
0178       if (fit_common_in.hits_valid[NUM_LAYERS * myisl + i] == 0)
0179         continue;
0180       else if (fit_common_in.hits[NUM_LAYERS * myisl + i].ti < smallest_time)
0181         smallest_time = fit_common_in.hits[NUM_LAYERS * myisl + i].ti;
0182     }
0183     if (fit_common_in.hits_valid[NUM_LAYERS * myisl + 0] == 1)
0184       tmp_coarse_wirepos_1 = fit_common_in.hits[NUM_LAYERS * myisl + 0].wp;
0185     else
0186       tmp_coarse_wirepos_1 = fit_common_in.hits[NUM_LAYERS * myisl + 1].wp;
0187     if (fit_common_in.hits_valid[NUM_LAYERS * myisl + 3] == 1)
0188       tmp_coarse_wirepos_3 = fit_common_in.hits[NUM_LAYERS * myisl + 3].wp;
0189     else
0190       tmp_coarse_wirepos_3 = fit_common_in.hits[NUM_LAYERS * myisl + 2].wp;
0191
0192     tmp_coarse_wirepos_1 = tmp_coarse_wirepos_1 >> WIREPOS_NORM_LSB_IGNORED;
0193     tmp_coarse_wirepos_3 = tmp_coarse_wirepos_3 >> WIREPOS_NORM_LSB_IGNORED;
0194   }
0195   fit_common_in.coarse_bctr = smallest_time >> (WIDTH_FULL_TIME - WIDTH_COARSED_TIME);
0196   fit_common_in.coarse_wirepos = (tmp_coarse_wirepos_1 + tmp_coarse_wirepos_3) >> 1;
0197
0198   for (auto &lat_comb : lat_combs) {
0199     if (lat_comb[0] == 0 && lat_comb[1] == 0 && lat_comb[2] == 0 && lat_comb[3] == 0)
0200       continue;
0201     fit_common_in.lateralities.clear();
0202
0203     auto rom_addr = get_rom_addr(inMPath, lat_comb);
0204     coeffs_t coeffs;
0205     if (sl == 0) {
0206       coeffs =
0207           RomDataConvert(lut_sl1[rom_addr], COEFF_WIDTH_SL_T0, COEFF_WIDTH_SL_POSITION, COEFF_WIDTH_SL_SLOPE, 0, 3);
0208     } else if (sl == 1) {
0209       coeffs =
0210           RomDataConvert(lut_sl2[rom_addr], COEFF_WIDTH_SL_T0, COEFF_WIDTH_SL2_POSITION, COEFF_WIDTH_SL_SLOPE, 0, 3);
0211     } else {
0212       coeffs =
0213           RomDataConvert(lut_sl3[rom_addr], COEFF_WIDTH_SL_T0, COEFF_WIDTH_SL_POSITION, COEFF_WIDTH_SL_SLOPE, 4, 7);
0214     }
0215     // Filling lateralities
0216     int minISL = 1;
0217     int maxISL = 3;
0218     if (sl < 1) {
0219       minISL = 0;
0220       maxISL = 2;
0221     }
0222
0223     for (int isl = minISL; isl < maxISL; isl++) {
0224       for (size_t i = 0; i < NUM_LAYERS; i++) {
0225         if (isl == sl) {
0226           fit_common_in.lateralities.push_back(lat_comb[i]);
0227         } else
0228           fit_common_in.lateralities.push_back(-1);
0229       }
0230     }
0231     fit_common_in.coeffs = coeffs;
0232
0233     auto fit_common_out = fit(fit_common_in,
0234                               XI_SL_WIDTH,
0235                               COEFF_WIDTH_SL_T0,
0236                               sl == 1 ? COEFF_WIDTH_SL2_POSITION : COEFF_WIDTH_SL_POSITION,
0237                               COEFF_WIDTH_SL_SLOPE,
0238                               PRECISSION_SL_T0,
0239                               PRECISSION_SL_POSITION,
0240                               PRECISSION_SL_SLOPE,
0241                               PROD_RESIZE_SL_T0,
0242                               sl == 1 ? PROD_RESIZE_SL2_POSITION : PROD_RESIZE_SL_POSITION,
0243                               PROD_RESIZE_SL_SLOPE,
0244                               max_drift_tdc,
0245                               sl + 1);
0246
0247     if (fit_common_out.valid_fit == 1) {
0248       float t0_f = ((float)fit_common_out.t0) * (float)LHC_CLK_FREQ / (float)TIME_TO_TDC_COUNTS;
0249
0250       float slope_f = -fit_common_out.slope * ((float)CELL_SEMILENGTH / max_drift_tdc) * (1) / (CELL_SEMIHEIGHT * 16.);
0251       if (sl != 1 && std::abs(slope_f) > tanPhiTh_)
0252         continue;
0253
0254       DTWireId wireId(MuonPathSLId, 2, 1);
0255       float pos_ch_f = (float)(fit_common_out.position) * ((float)CELL_SEMILENGTH / (float)max_drift_tdc) / 10;
0256       pos_ch_f += (SL1_CELLS_OFFSET * CELL_LENGTH) / 10.;
0257       if (sl != 1)
0258         pos_ch_f += shiftinfo_[wireIdSL1.rawId()];
0259       else if (sl == 1)
0260         pos_ch_f += shiftthetainfo_[wireIdSL2.rawId()];
0261
0262       float pos_sl_f = pos_ch_f - (sl - 1) * slope_f * VERT_PHI1_PHI3 / 2;
0263       float chi2_f = fit_common_out.chi2 * std::pow(((float)CELL_SEMILENGTH / (float)max_drift_tdc), 2) / 100;
0264
0265       // obtention of global coordinates using luts
0266       int pos = (int)(10 * (pos_sl_f - shiftinfo_[wireId.rawId()]) * INCREASED_RES_POS_POW);
0267       int slope = (int)(-slope_f * INCREASED_RES_SLOPE_POW);
0268
0269       auto global_coords = globalcoordsobtainer_->get_global_coordinates(ChId.rawId(), sl + 1, pos, slope);
0270       float phi = global_coords[0];
0271       float phiB = global_coords[1];
0272
0273       // obtention of global coordinates using cmssw geometry
0274       double z = 0;
0275       if (ChId.station() == 3 or ChId.station() == 4) {
0276         z = Z_SHIFT_MB4;
0277       }
0278       GlobalPoint jm_x_cmssw_global = dtGeo_->chamber(ChId)->toGlobal(LocalPoint(pos_sl_f, 0., z));
0279       int thisec = ChId.sector();
0280       if (thisec == 13)
0281         thisec = 4;
0282       if (thisec == 14)
0283         thisec = 10;
0284       float phi_cmssw = jm_x_cmssw_global.phi() - PHI_CONV * (thisec - 1);
0285       float psi = atan(slope_f);
0286       float phiB_cmssw = hasPosRF(ChId.wheel(), ChId.sector()) ? psi - phi_cmssw : -psi - phi_cmssw;
0287       if (sl == 0)
0288         metaPrimitives.emplace_back(metaPrimitive({MuonPathSLId.rawId(),
0289                                                    t0_f,
0290                                                    (double)(fit_common_out.position),
0291                                                    (double)fit_common_out.slope,
0292                                                    phi,
0293                                                    phiB,
0294                                                    phi_cmssw,
0295                                                    phiB_cmssw,
0296                                                    chi2_f,
0297                                                    quality,
0298                                                    inMPath->primitive(0)->channelId(),
0299                                                    inMPath->primitive(0)->tdcTimeStamp(),
0300                                                    lat_comb[0],
0301                                                    inMPath->primitive(1)->channelId(),
0302                                                    inMPath->primitive(1)->tdcTimeStamp(),
0303                                                    lat_comb[1],
0304                                                    inMPath->primitive(2)->channelId(),
0305                                                    inMPath->primitive(2)->tdcTimeStamp(),
0306                                                    lat_comb[2],
0307                                                    inMPath->primitive(3)->channelId(),
0308                                                    inMPath->primitive(3)->tdcTimeStamp(),
0309                                                    lat_comb[3],
0310                                                    -1,
0311                                                    -1,
0312                                                    -1,
0313                                                    -1,
0314                                                    -1,
0315                                                    -1,
0316                                                    -1,
0317                                                    -1,
0318                                                    -1,
0319                                                    -1,
0320                                                    -1,
0321                                                    -1,
0322                                                    -1}));
0323       else if (sl == 2)
0324         metaPrimitives.emplace_back(metaPrimitive({MuonPathSLId.rawId(),
0325                                                    t0_f,
0326                                                    (double)(fit_common_out.position),
0327                                                    (double)fit_common_out.slope,
0328                                                    phi,
0329                                                    phiB,
0330                                                    phi_cmssw,
0331                                                    phiB_cmssw,
0332                                                    chi2_f,
0333                                                    quality,
0334                                                    -1,
0335                                                    -1,
0336                                                    -1,
0337                                                    -1,
0338                                                    -1,
0339                                                    -1,
0340                                                    -1,
0341                                                    -1,
0342                                                    -1,
0343                                                    -1,
0344                                                    -1,
0345                                                    -1,
0346                                                    inMPath->primitive(0)->channelId(),
0347                                                    inMPath->primitive(0)->tdcTimeStamp(),
0348                                                    lat_comb[0],
0349                                                    inMPath->primitive(1)->channelId(),
0350                                                    inMPath->primitive(1)->tdcTimeStamp(),
0351                                                    lat_comb[1],
0352                                                    inMPath->primitive(2)->channelId(),
0353                                                    inMPath->primitive(2)->tdcTimeStamp(),
0354                                                    lat_comb[2],
0355                                                    inMPath->primitive(3)->channelId(),
0356                                                    inMPath->primitive(3)->tdcTimeStamp(),
0357                                                    lat_comb[3],
0358                                                    -1}));
0359       else if (sl == 1) {
0360         // fw-like calculation
0361         DTLayerId SL2_layer2Id(MuonPathSLId, 2);
0362         double z_shift = shiftthetainfo_[SL2_layer2Id.rawId()];
0363         double pos_cm =
0364             pos / 10 / INCREASED_RES_POS_POW;  // fixed to have z_shift and the position in the same units (MC)
0365         double jm_y = hasPosRF(MuonPathSLId.wheel(), MuonPathSLId.sector()) ? z_shift - pos_cm : z_shift + pos_cm;
0366
0367         phi = jm_y;
0368
0369         // Fixed sign of k (MC)
0370         double k_fromfw = hasPosRF(MuonPathSLId.wheel(), MuonPathSLId.sector()) ? slope_f : -slope_f;
0371         phiB = k_fromfw;
0372
0373         // cmssw-like calculation
0374         LocalPoint wire1_in_layer(dtGeo_->layer(SL2_layer2Id)->specificTopology().wirePosition(1), 0, -0.65);
0375         GlobalPoint wire1_in_global = dtGeo_->layer(SL2_layer2Id)->toGlobal(wire1_in_layer);
0376         LocalPoint wire1_in_sl = dtGeo_->superLayer(MuonPathSLId)->toLocal(wire1_in_global);
0377         double x_shift = wire1_in_sl.x();
0378         jm_y = (dtGeo_->superLayer(MuonPathSLId)
0379                     ->toGlobal(LocalPoint(double(pos) / (10 * pow(2, INCREASED_RES_POS)) + x_shift, 0., 0)))
0380                    .z();
0381         phi_cmssw = jm_y;
0382
0383         double x_cmssw = (dtGeo_->superLayer(MuonPathSLId)
0384                               ->toGlobal(LocalPoint(double(pos) / (10 * pow(2, INCREASED_RES_POS)) + x_shift, 0., 0)))
0385                              .x();
0386         double y_cmssw = (dtGeo_->superLayer(MuonPathSLId)
0387                               ->toGlobal(LocalPoint(double(pos) / (10 * pow(2, INCREASED_RES_POS)) + x_shift, 0., 0)))
0388                              .y();
0389         double r_cmssw = sqrt(x_cmssw * x_cmssw + y_cmssw * y_cmssw);
0390         double k_cmssw = jm_y / r_cmssw;
0391         phiB_cmssw = k_cmssw;
0392
0393         metaPrimitives.emplace_back(metaPrimitive({MuonPathSLId.rawId(),
0394                                                    t0_f,
0395                                                    (double)(fit_common_out.position),
0396                                                    (double)fit_common_out.slope,
0397                                                    phi,
0398                                                    phiB,
0399                                                    phi_cmssw,
0400                                                    phiB_cmssw,
0401                                                    chi2_f,
0402                                                    quality,
0403                                                    inMPath->primitive(0)->channelId(),
0404                                                    inMPath->primitive(0)->tdcTimeStamp(),
0405                                                    lat_comb[0],
0406                                                    inMPath->primitive(1)->channelId(),
0407                                                    inMPath->primitive(1)->tdcTimeStamp(),
0408                                                    lat_comb[1],
0409                                                    inMPath->primitive(2)->channelId(),
0410                                                    inMPath->primitive(2)->tdcTimeStamp(),
0411                                                    lat_comb[2],
0412                                                    inMPath->primitive(3)->channelId(),
0413                                                    inMPath->primitive(3)->tdcTimeStamp(),
0414                                                    lat_comb[3],
0415                                                    -1,
0416                                                    -1,
0417                                                    -1,
0418                                                    -1,
0419                                                    -1,
0420                                                    -1,
0421                                                    -1,
0422                                                    -1,
0423                                                    -1,
0424                                                    -1,
0425                                                    -1,
0426                                                    -1,
0427                                                    -1}));
0428       }
0429     }  // (fit_common_out.valid_fit == 1)
0430   }    // loop in lat_combs
0431   return;
0432 }
0433
0434 void MuonPathSLFitter::fillLuts() {
0435   std::ifstream ifinsl1(sl1_filename_.fullPath());
0436   std::string line;
0437   while (ifinsl1.good()) {
0438     ifinsl1 >> line;
0439
0440     std::vector<int> myNumbers;
0441     for (size_t i = 0; i < line.size(); i++) {
0442       // This converts the char into an int and pushes it into vec
0443       myNumbers.push_back(line[i] - '0');  // The digits will be in the same order as before
0444     }
0445     std::reverse(myNumbers.begin(), myNumbers.end());
0446     lut_sl1.push_back(myNumbers);
0447   }
0448
0449   std::ifstream ifinsl2(sl2_filename_.fullPath());
0450   while (ifinsl2.good()) {
0451     ifinsl2 >> line;
0452
0453     std::vector<int> myNumbers;
0454     for (size_t i = 0; i < line.size(); i++) {
0455       // This converts the char into an int and pushes it into vec
0456       myNumbers.push_back(line[i] - '0');  // The digits will be in the same order as before
0457     }
0458     std::reverse(myNumbers.begin(), myNumbers.end());
0459     lut_sl2.push_back(myNumbers);
0460   }
0461
0462   std::ifstream ifinsl3(sl3_filename_.fullPath());
0463   while (ifinsl3.good()) {
0464     ifinsl3 >> line;
0465
0466     std::vector<int> myNumbers;
0467     for (size_t i = 0; i < line.size(); i++) {
0468       // This converts the char into an int and pushes it into vec
0469       myNumbers.push_back(line[i] - '0');  // The digits will be in the same order as before
0470     }
0471     std::reverse(myNumbers.begin(), myNumbers.end());
0472     lut_sl3.push_back(myNumbers);
0473   }
0474
0475   return;
0476 }
0477
0478 int MuonPathSLFitter::get_rom_addr(MuonPathPtr &inMPath, latcomb lats) {
0479   std::vector<int> rom_addr;
0480   auto missing_layer = inMPath->missingLayer();
0481   if (missing_layer == -1) {
0482     rom_addr.push_back(1);
0483     rom_addr.push_back(0);
0484   } else {
0485     if (missing_layer == 0) {
0486       rom_addr.push_back(0);
0487       rom_addr.push_back(0);
0488     } else if (missing_layer == 1) {
0489       rom_addr.push_back(0);
0490       rom_addr.push_back(1);
0491     } else if (missing_layer == 2) {
0492       rom_addr.push_back(1);
0493       rom_addr.push_back(0);
0494     } else {  // missing_layer == 3
0495       rom_addr.push_back(1);
0496       rom_addr.push_back(1);
0497     }
0498   }
0499   for (size_t ilat = 0; ilat < lats.size(); ilat++) {
0500     if ((int)ilat == missing_layer)  // only applies to 3-hit, as in 4-hit missL=-1
0501       continue;
0502     auto lat = lats[ilat];
0503     if (lat == -1)
0504       lat = 0;
0505     rom_addr.push_back(lat);
0506   }
0507   std::reverse(rom_addr.begin(), rom_addr.end());
0508   return vhdl_unsigned_to_int(rom_addr);
0509 }