File indexing completed on 2024-09-07 04:36:47
0001 #include "L1Trigger/DTTriggerPhase2/interface/MuonPathCorFitter.h"
0002 #include <cmath>
0003 #include <memory>
0004
0005 using namespace edm;
0006 using namespace std;
0007 using namespace cmsdt;
0008
0009
0010
0011 MuonPathCorFitter::MuonPathCorFitter(const ParameterSet& pset,
0012 edm::ConsumesCollector& iC,
0013 std::shared_ptr<GlobalCoordsObtainer>& globalcoordsobtainer)
0014 : MuonPathFitter(pset, iC, globalcoordsobtainer), dT0_correlate_TP_(pset.getParameter<double>("dT0_correlate_TP")) {
0015 if (debug_)
0016 LogDebug("MuonPathCorFitter") << "MuonPathCorFitter: constructor";
0017
0018
0019 both_sl_filename_ = pset.getParameter<edm::FileInPath>("lut_2sl");
0020
0021 fillLuts();
0022
0023 setChi2Th(pset.getParameter<double>("chi2corTh"));
0024 setTanPhiTh(pset.getParameter<double>("dTanPsi_correlate_TP"));
0025 }
0026
0027 MuonPathCorFitter::~MuonPathCorFitter() {
0028 if (debug_)
0029 LogDebug("MuonPathCorFitter") << "MuonPathCorFitter: destructor";
0030 }
0031
0032
0033
0034
0035 void MuonPathCorFitter::initialise(const edm::EventSetup& iEventSetup) {
0036 if (debug_)
0037 LogDebug("MuonPathCorFitter") << "MuonPathCorFitter::initialiase";
0038
0039 auto geom = iEventSetup.getHandle(dtGeomH);
0040 dtGeo_ = &(*geom);
0041 }
0042
0043 void MuonPathCorFitter::run(edm::Event& iEvent,
0044 const edm::EventSetup& iEventSetup,
0045 std::vector<cmsdt::metaPrimitive>& inMPaths,
0046 std::vector<cmsdt::metaPrimitive>& outMPaths) {
0047 if (debug_)
0048 LogDebug("MuonPathCorFitter") << "MuonPathCorFitter: run";
0049 if (!inMPaths.empty()) {
0050 int dum_sl_rawid = inMPaths[0].rawId;
0051 DTSuperLayerId dumSlId(dum_sl_rawid);
0052 DTChamberId ChId(dumSlId.wheel(), dumSlId.station(), dumSlId.sector());
0053 max_drift_tdc = maxdriftinfo_[dumSlId.wheel() + 2][dumSlId.station() - 1][dumSlId.sector() - 1];
0054 DTSuperLayerId sl1Id(ChId.rawId(), 1);
0055 DTSuperLayerId sl3Id(ChId.rawId(), 3);
0056
0057 std::map<int, std::vector<metaPrimitive>> SL1metaPrimitivesPerBX;
0058 std::map<int, std::vector<metaPrimitive>> SL3metaPrimitivesPerBX;
0059 for (const auto& metaprimitiveIt : inMPaths) {
0060 int BX = metaprimitiveIt.t0 / 25;
0061 if (metaprimitiveIt.rawId == sl1Id.rawId())
0062 SL1metaPrimitivesPerBX[BX].push_back(metaprimitiveIt);
0063 else if (metaprimitiveIt.rawId == sl3Id.rawId())
0064 SL3metaPrimitivesPerBX[BX].push_back(metaprimitiveIt);
0065 }
0066
0067 std::vector<bx_sl_vector> bxs_to_consider;
0068 bxs_to_consider.reserve(SL1metaPrimitivesPerBX.size());
0069 for (auto& prims_sl1 : SL1metaPrimitivesPerBX)
0070 bxs_to_consider.push_back(bx_sl_vector({prims_sl1.first, prims_sl1.second, 1}));
0071
0072 for (auto& prims_sl3 : SL3metaPrimitivesPerBX)
0073 bxs_to_consider.push_back(bx_sl_vector({prims_sl3.first, prims_sl3.second, 3}));
0074
0075 std::stable_sort(bxs_to_consider.begin(), bxs_to_consider.end(), bxSort);
0076
0077 std::vector<mp_group> mps_q8;
0078 std::vector<mp_group> mps_q7;
0079 std::vector<mp_group> mps_q6;
0080
0081 for (size_t ibx = 1; ibx < bxs_to_consider.size(); ibx++) {
0082 for (size_t ibx2 = 0; ibx2 < ibx; ibx2++) {
0083 if (bxs_to_consider[ibx].sl != bxs_to_consider[ibx2].sl &&
0084 (abs(bxs_to_consider[ibx].bx - bxs_to_consider[ibx2].bx)) <= MAX_BX_FOR_COR) {
0085 int isl1 = 0;
0086 for (auto& prim1 : bxs_to_consider[ibx].mps) {
0087 if (isl1 >= MAX_PRIM_PER_BX_FOR_COR)
0088 break;
0089 int isl2 = 0;
0090 for (auto& prim2 : bxs_to_consider[ibx2].mps) {
0091 if (isl2 >= MAX_PRIM_PER_BX_FOR_COR)
0092 break;
0093 if (bxs_to_consider[ibx].sl == 1) {
0094 if (!canCorrelate(prim1, prim2)) {
0095 continue;
0096 }
0097 if (prim1.quality >= 3 && prim2.quality >= 3)
0098 mps_q8.push_back(mp_group({prim1, prim2}));
0099 else if ((prim1.quality >= 3 && prim2.quality < 3) || (prim1.quality < 3 && prim2.quality >= 3))
0100 mps_q7.push_back(mp_group({prim1, prim2}));
0101 else
0102 mps_q6.push_back(mp_group({prim1, prim2}));
0103 } else {
0104 if (!canCorrelate(prim2, prim1)) {
0105 continue;
0106 }
0107 if (prim2.quality >= 3 && prim1.quality >= 3)
0108 mps_q8.push_back(mp_group({prim2, prim1}));
0109 else if ((prim2.quality >= 3 && prim1.quality < 3) || (prim2.quality < 3 && prim1.quality >= 3))
0110 mps_q7.push_back(mp_group({prim2, prim1}));
0111 else
0112 mps_q6.push_back(mp_group({prim2, prim1}));
0113 }
0114 isl2++;
0115 }
0116 isl1++;
0117 }
0118 }
0119 }
0120 }
0121 int iq = 0;
0122 for (size_t i = 0; i < mps_q8.size(); i++) {
0123 if (iq >= MAX_PRIM_FOR_COR)
0124 break;
0125 analyze(mps_q8[i], outMPaths);
0126 iq += 1;
0127 }
0128 for (size_t i = 0; i < mps_q7.size(); i++) {
0129 if (iq >= MAX_PRIM_FOR_COR)
0130 break;
0131 analyze(mps_q7[i], outMPaths);
0132 iq += 1;
0133 }
0134 for (size_t i = 0; i < mps_q6.size(); i++) {
0135 if (iq >= MAX_PRIM_FOR_COR)
0136 break;
0137 analyze(mps_q6[i], outMPaths);
0138 iq += 1;
0139 }
0140 }
0141 }
0142
0143 bool MuonPathCorFitter::canCorrelate(cmsdt::metaPrimitive mp_sl1, cmsdt::metaPrimitive mp_sl3) {
0144
0145 float pos_ch_sl1_f = mp_sl1.x;
0146 float pos_ch_sl3_f = mp_sl3.x;
0147
0148
0149 int pos_ch_sl1 = int(pos_ch_sl1_f);
0150 int pos_ch_sl3 = int(pos_ch_sl3_f);
0151
0152 int slope_sl1 = (int)mp_sl1.tanPhi;
0153 int slope_sl3 = (int)mp_sl3.tanPhi;
0154
0155 if (abs((slope_sl1 >> WIDTH_POS_SLOPE_CORR) - (slope_sl3 >> WIDTH_POS_SLOPE_CORR)) > 1)
0156 return false;
0157
0158 if (abs((pos_ch_sl1 >> WIDTH_POS_SLOPE_CORR) - (pos_ch_sl3 >> WIDTH_POS_SLOPE_CORR)) > 1)
0159 return false;
0160
0161 if (abs(mp_sl1.t0 - mp_sl3.t0) > dT0_correlate_TP_)
0162 return false;
0163
0164 return true;
0165 }
0166
0167 void MuonPathCorFitter::finish() {
0168 if (debug_)
0169 LogDebug("MuonPathCorFitter") << "MuonPathCorFitter: finish";
0170 };
0171
0172
0173
0174
0175
0176 void MuonPathCorFitter::analyze(mp_group mp, std::vector<cmsdt::metaPrimitive>& metaPrimitives) {
0177
0178 DTSuperLayerId MuonPathSLId(mp[0].rawId);
0179
0180 DTChamberId ChId(MuonPathSLId.wheel(), MuonPathSLId.station(), MuonPathSLId.sector());
0181
0182 DTSuperLayerId MuonPathSL1Id(ChId.wheel(), ChId.station(), ChId.sector(), 1);
0183 DTSuperLayerId MuonPathSL3Id(ChId.wheel(), ChId.station(), ChId.sector(), 3);
0184 DTWireId wireIdSL1(MuonPathSL1Id, 2, 1);
0185 DTWireId wireIdSL3(MuonPathSL3Id, 2, 1);
0186 auto sl_shift_cm = shiftinfo_[wireIdSL1.rawId()] - shiftinfo_[wireIdSL3.rawId()];
0187
0188 fit_common_in_t fit_common_in;
0189
0190
0191 fit_common_in.hits = {};
0192 fit_common_in.hits_valid = {};
0193 short quality = 0;
0194 if (mp[0].quality >= 3 && mp[1].quality >= 3)
0195 quality = 8;
0196 else if ((mp[0].quality >= 3 && mp[1].quality < 3) || (mp[0].quality < 3 && mp[1].quality >= 3))
0197 quality = 7;
0198 else
0199 quality = 6;
0200
0201 std::vector<int> missing_layers;
0202
0203 for (int isl = 0; isl < 2; isl++) {
0204 int wire[4], tdc[4];
0205 if (isl != 1) {
0206 wire[0] = mp[isl].wi1;
0207 tdc[0] = mp[isl].tdc1;
0208 wire[1] = mp[isl].wi2;
0209 tdc[1] = mp[isl].tdc2;
0210 wire[2] = mp[isl].wi3;
0211 tdc[2] = mp[isl].tdc3;
0212 wire[3] = mp[isl].wi4;
0213 tdc[3] = mp[isl].tdc4;
0214 } else {
0215 wire[0] = mp[isl].wi5;
0216 tdc[0] = mp[isl].tdc5;
0217 wire[1] = mp[isl].wi6;
0218 tdc[1] = mp[isl].tdc6;
0219 wire[2] = mp[isl].wi7;
0220 tdc[2] = mp[isl].tdc7;
0221 wire[3] = mp[isl].wi8;
0222 tdc[3] = mp[isl].tdc8;
0223 }
0224
0225 for (int i = 0; i < NUM_LAYERS; i++) {
0226 if (wire[i] != -1) {
0227
0228 auto ti = tdc[i];
0229 if (ti != -1)
0230 ti = (int)round(((float)TIME_TO_TDC_COUNTS / (float)LHC_CLK_FREQ) * ti);
0231 auto wi = wire[i];
0232 auto ly = i;
0233
0234 int wp_semicells = (wi - SL1_CELLS_OFFSET) * 2 + 1;
0235 if (ly % 2 == 1)
0236 wp_semicells -= 1;
0237 if (isl == 1)
0238 wp_semicells -= (int)round((sl_shift_cm * 10) / CELL_SEMILENGTH);
0239 float wp_tdc = wp_semicells * max_drift_tdc;
0240 int wp = (int)((long int)(round(wp_tdc * std::pow(2, WIREPOS_WIDTH))) / (int)std::pow(2, WIREPOS_WIDTH));
0241
0242
0243 fit_common_in.hits.push_back({ti, wi, ly, wp});
0244
0245 fit_common_in.hits_valid.push_back(1);
0246 } else {
0247 missing_layers.push_back(isl * NUM_LAYERS + i);
0248 fit_common_in.hits.push_back({-1, -1, -1, -1});
0249 fit_common_in.hits_valid.push_back(0);
0250 }
0251 }
0252 }
0253
0254 int smallest_time = 999999, tmp_coarse_wirepos_1 = -1, tmp_coarse_wirepos_3 = -1;
0255
0256 for (int isl = 0; isl < 2; isl++) {
0257 for (size_t i = 0; i < NUM_LAYERS; i++) {
0258 if (fit_common_in.hits_valid[NUM_LAYERS * isl + i] == 0)
0259 continue;
0260 else if (fit_common_in.hits[NUM_LAYERS * isl + i].ti < smallest_time)
0261 smallest_time = fit_common_in.hits[NUM_LAYERS * isl + i].ti;
0262 }
0263 }
0264 if (fit_common_in.hits_valid[NUM_LAYERS * 0 + 0] == 1)
0265 tmp_coarse_wirepos_1 = fit_common_in.hits[NUM_LAYERS * 0 + 0].wp;
0266 else
0267 tmp_coarse_wirepos_1 = fit_common_in.hits[NUM_LAYERS * 0 + 1].wp;
0268 if (fit_common_in.hits_valid[NUM_LAYERS * 1 + 3] == 1)
0269 tmp_coarse_wirepos_3 = fit_common_in.hits[NUM_LAYERS * 1 + 3].wp;
0270 else
0271 tmp_coarse_wirepos_3 = fit_common_in.hits[NUM_LAYERS * 1 + 2].wp;
0272
0273 tmp_coarse_wirepos_1 = tmp_coarse_wirepos_1 >> WIREPOS_NORM_LSB_IGNORED;
0274 tmp_coarse_wirepos_3 = tmp_coarse_wirepos_3 >> WIREPOS_NORM_LSB_IGNORED;
0275
0276 fit_common_in.coarse_bctr = smallest_time >> (WIDTH_FULL_TIME - WIDTH_COARSED_TIME);
0277 fit_common_in.coarse_wirepos = (tmp_coarse_wirepos_1 + tmp_coarse_wirepos_3) >> 1;
0278
0279 fit_common_in.lateralities.clear();
0280
0281 auto rom_addr = get_rom_addr(mp, missing_layers);
0282
0283 coeffs_t coeffs =
0284 RomDataConvert(lut_2sl[rom_addr], COEFF_WIDTH_COR_T0, COEFF_WIDTH_COR_POSITION, COEFF_WIDTH_COR_SLOPE, 0, 7);
0285
0286
0287 for (int isl = 0; isl < 2; isl++) {
0288 int lat[4];
0289 if (isl != 1) {
0290 lat[0] = mp[isl].lat1;
0291 lat[1] = mp[isl].lat2;
0292 lat[2] = mp[isl].lat3;
0293 lat[3] = mp[isl].lat4;
0294 } else {
0295 lat[0] = mp[isl].lat5;
0296 lat[1] = mp[isl].lat6;
0297 lat[2] = mp[isl].lat7;
0298 lat[3] = mp[isl].lat8;
0299 }
0300
0301 for (size_t i = 0; i < NUM_LAYERS; i++) {
0302 fit_common_in.lateralities.push_back(lat[i]);
0303 }
0304 }
0305
0306 fit_common_in.coeffs = coeffs;
0307
0308 auto fit_common_out = fit(fit_common_in,
0309 XI_COR_WIDTH,
0310 COEFF_WIDTH_COR_T0,
0311 COEFF_WIDTH_COR_POSITION,
0312 COEFF_WIDTH_COR_SLOPE,
0313 PRECISSION_COR_T0,
0314 PRECISSION_COR_POSITION,
0315 PRECISSION_COR_SLOPE,
0316 PROD_RESIZE_COR_T0,
0317 PROD_RESIZE_COR_POSITION,
0318 PROD_RESIZE_COR_SLOPE,
0319 max_drift_tdc,
0320 0);
0321
0322 if (fit_common_out.valid_fit == 1) {
0323 float t0_f = ((float)fit_common_out.t0) * (float)LHC_CLK_FREQ / (float)TIME_TO_TDC_COUNTS;
0324 float slope_f = -fit_common_out.slope * ((float)CELL_SEMILENGTH / max_drift_tdc) * (1) / (CELL_SEMIHEIGHT * 16.);
0325 if (std::abs(slope_f) > tanPhiTh_)
0326 return;
0327
0328 DTWireId wireId(MuonPathSLId, 2, 1);
0329 float pos_ch_f = (float)(fit_common_out.position) * ((float)CELL_SEMILENGTH / (float)max_drift_tdc) / 10;
0330 pos_ch_f += (SL1_CELLS_OFFSET * CELL_LENGTH) / 10.;
0331 pos_ch_f += shiftinfo_[wireId.rawId()];
0332
0333 float chi2_f = fit_common_out.chi2 * std::pow(((float)CELL_SEMILENGTH / (float)max_drift_tdc), 2) / 100;
0334
0335
0336
0337 LogDebug("MuonPathCorFitter") << "==================== CORRELATED PRIMITIVE =================================";
0338 LogDebug("MuonPathCorFitter") << "WHEEL = " << ChId.wheel();
0339 LogDebug("MuonPathCorFitter") << "SECTOR = " << ChId.sector();
0340 LogDebug("MuonPathCorFitter") << "STATION = " << ChId.station();
0341 LogDebug("MuonPathCorFitter") << "QUALITY = " << quality;
0342 LogDebug("MuonPathCorFitter") << "POSITION = " << (double)fit_common_out.position;
0343 LogDebug("MuonPathCorFitter") << "SLOPE = " << (double)fit_common_out.slope;
0344
0345 auto global_coords =
0346 globalcoordsobtainer_->get_global_coordinates(ChId.rawId(), 0, fit_common_out.position, fit_common_out.slope);
0347 float phi = global_coords[0];
0348 float phiB = global_coords[1];
0349
0350
0351 double z = 0;
0352 if (ChId.station() == 3 or ChId.station() == 4) {
0353 z += Z_SHIFT_MB4;
0354 }
0355 GlobalPoint jm_x_cmssw_global = dtGeo_->chamber(ChId)->toGlobal(LocalPoint(pos_ch_f, 0., z));
0356 int thisec = ChId.sector();
0357 if (thisec == 13)
0358 thisec = 4;
0359 if (thisec == 14)
0360 thisec = 10;
0361 float phi_cmssw = jm_x_cmssw_global.phi() - PHI_CONV * (thisec - 1);
0362 float psi = atan(slope_f);
0363 float phiB_cmssw = hasPosRF(ChId.wheel(), ChId.sector()) ? psi - phi_cmssw : -psi - phi_cmssw;
0364 metaPrimitives.emplace_back(metaPrimitive({MuonPathSLId.rawId(),
0365 t0_f,
0366 (double)fit_common_out.position,
0367 (double)fit_common_out.slope,
0368 phi,
0369 phiB,
0370 phi_cmssw,
0371 phiB_cmssw,
0372 chi2_f,
0373 quality,
0374 mp[0].wi1,
0375 mp[0].tdc1,
0376 mp[0].lat1,
0377 mp[0].wi2,
0378 mp[0].tdc2,
0379 mp[0].lat2,
0380 mp[0].wi3,
0381 mp[0].tdc3,
0382 mp[0].lat3,
0383 mp[0].wi4,
0384 mp[0].tdc4,
0385 mp[0].lat4,
0386 mp[1].wi5,
0387 mp[1].tdc5,
0388 mp[1].lat5,
0389 mp[1].wi6,
0390 mp[1].tdc6,
0391 mp[1].lat6,
0392 mp[1].wi7,
0393 mp[1].tdc7,
0394 mp[1].lat7,
0395 mp[1].wi8,
0396 mp[1].tdc8,
0397 mp[1].lat8,
0398 -1}));
0399 }
0400 return;
0401 }
0402
0403 void MuonPathCorFitter::fillLuts() {
0404 std::ifstream ifin2sl(both_sl_filename_.fullPath());
0405 std::string line;
0406 while (ifin2sl.good()) {
0407 ifin2sl >> line;
0408
0409 std::vector<int> myNumbers;
0410 for (size_t i = 0; i < line.size(); i++) {
0411
0412 myNumbers.push_back(line[i] - '0');
0413 }
0414 std::reverse(myNumbers.begin(), myNumbers.end());
0415 lut_2sl.push_back(myNumbers);
0416 }
0417
0418 return;
0419 }
0420
0421 int MuonPathCorFitter::get_rom_addr(mp_group mps, std::vector<int> missing_hits) {
0422 std::vector<int> lats = {
0423 mps[0].lat1, mps[0].lat2, mps[0].lat3, mps[0].lat4, mps[1].lat5, mps[1].lat6, mps[1].lat7, mps[1].lat8};
0424
0425 std::vector<int> rom_addr;
0426 if (missing_hits.size() == 1)
0427 rom_addr.push_back(1);
0428 else
0429 rom_addr.push_back(0);
0430
0431 if (missing_hits.size() == 1) {
0432 if (missing_hits[0] < 4)
0433 rom_addr.push_back(0);
0434 else
0435 rom_addr.push_back(1);
0436 if (missing_hits[0] % 4 == 0) {
0437 rom_addr.push_back(0);
0438 rom_addr.push_back(0);
0439 } else if (missing_hits[0] % 4 == 1) {
0440 rom_addr.push_back(0);
0441 rom_addr.push_back(1);
0442 } else if (missing_hits[0] % 4 == 2) {
0443 rom_addr.push_back(1);
0444 rom_addr.push_back(0);
0445 } else {
0446 rom_addr.push_back(1);
0447 rom_addr.push_back(1);
0448 }
0449 for (size_t ilat = 0; ilat < lats.size(); ilat++) {
0450 if ((int)ilat == missing_hits[0])
0451 continue;
0452 auto lat = lats[ilat];
0453 if (lat == -1)
0454 lat = 0;
0455 rom_addr.push_back(lat);
0456 }
0457
0458 } else if (missing_hits.empty()) {
0459 for (size_t ilat = 0; ilat < lats.size(); ilat++) {
0460 auto lat = lats[ilat];
0461 if (lat == -1)
0462 lat = 0;
0463 rom_addr.push_back(lat);
0464 }
0465 auto lat = lats[NUM_LAYERS + 3];
0466 if (lat == -1)
0467 lat = 0;
0468 rom_addr.push_back(lat);
0469 rom_addr.push_back(lat);
0470
0471 } else {
0472 for (int i = missing_hits.size() - 1; i >= 0; i--) {
0473 if (missing_hits[i] % 4 == 0) {
0474 rom_addr.push_back(0);
0475 rom_addr.push_back(0);
0476 } else if (missing_hits[i] % 4 == 1) {
0477 rom_addr.push_back(0);
0478 rom_addr.push_back(1);
0479 } else if (missing_hits[i] % 4 == 2) {
0480 rom_addr.push_back(1);
0481 rom_addr.push_back(0);
0482 } else {
0483 rom_addr.push_back(1);
0484 rom_addr.push_back(1);
0485 }
0486 }
0487 for (size_t ilat = 0; ilat < lats.size(); ilat++) {
0488 if ((int)ilat == missing_hits[0] || (int)ilat == missing_hits[1])
0489 continue;
0490 auto lat = lats[ilat];
0491 if (lat == -1)
0492 lat = 0;
0493 rom_addr.push_back(lat);
0494 }
0495 }
0496 std::reverse(rom_addr.begin(), rom_addr.end());
0497 return vhdl_unsigned_to_int(rom_addr);
0498 }