DTTriggerPhase2/src/MuonPathCorFitter.cc

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 // Constructors and destructor
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   // LUTs
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 // Main methods (initialise, run, finish)
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       }  // looping over the 0 -> N-1 BX groups
0120     }    // looping over the 1 -> N BX groups
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   // moving position from SL RF to chamber RF
0145   float pos_ch_sl1_f = mp_sl1.x;
0146   float pos_ch_sl3_f = mp_sl3.x;
0147
0148   // translating into tdc counts
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 //--- Metodos privados
0174 //------------------------------------------------------------------
0175
0176 void MuonPathCorFitter::analyze(mp_group mp, std::vector<cmsdt::metaPrimitive>& metaPrimitives) {
0177   //FIXME
0178   DTSuperLayerId MuonPathSLId(mp[0].rawId);  // SL1
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   // 8-element vectors, for the 8 layers. As here we are fitting one SL only, we leave the other SL values as dummy ones
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         // Include both valid and non-valid hits. Non-valid values can be whatever, leaving all as -1 to make debugging easier.
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)  // SL3
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         // wp in tdc counts (still in floating point)
0243         fit_common_in.hits.push_back({ti, wi, ly, wp});
0244         // fill valids as well
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   // coarse_bctr is the 12 MSB of the smallest tdc
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   // Filling lateralities
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     // obtention of global coordinates using luts
0336     int pos = (int)(10 * (pos_ch_f - shiftinfo_[wireId.rawId()]) * INCREASED_RES_POS_POW);
0337     int slope = (int)(-slope_f * INCREASED_RES_SLOPE_POW);
0338     auto global_coords = globalcoordsobtainer_->get_global_coordinates(ChId.rawId(), 0, pos, slope);
0339     float phi = global_coords[0];
0340     float phiB = global_coords[1];
0341
0342     // obtention of global coordinates using cmssw geometry
0343     double z = 0;
0344     if (ChId.station() == 3 or ChId.station() == 4) {
0345       z += Z_SHIFT_MB4;
0346     }
0347     GlobalPoint jm_x_cmssw_global = dtGeo_->chamber(ChId)->toGlobal(LocalPoint(pos_ch_f, 0., z));
0348     int thisec = ChId.sector();
0349     if (thisec == 13)
0350       thisec = 4;
0351     if (thisec == 14)
0352       thisec = 10;
0353     float phi_cmssw = jm_x_cmssw_global.phi() - PHI_CONV * (thisec - 1);
0354     float psi = atan(slope_f);
0355     float phiB_cmssw = hasPosRF(ChId.wheel(), ChId.sector()) ? psi - phi_cmssw : -psi - phi_cmssw;
0356     metaPrimitives.emplace_back(metaPrimitive({MuonPathSLId.rawId(),
0357                                                t0_f,
0358                                                (double)fit_common_out.position,
0359                                                (double)fit_common_out.slope,
0360                                                phi,
0361                                                phiB,
0362                                                phi_cmssw,
0363                                                phiB_cmssw,
0364                                                chi2_f,
0365                                                quality,
0366                                                mp[0].wi1,
0367                                                mp[0].tdc1,
0368                                                mp[0].lat1,
0369                                                mp[0].wi2,
0370                                                mp[0].tdc2,
0371                                                mp[0].lat2,
0372                                                mp[0].wi3,
0373                                                mp[0].tdc3,
0374                                                mp[0].lat3,
0375                                                mp[0].wi4,
0376                                                mp[0].tdc4,
0377                                                mp[0].lat4,
0378                                                mp[1].wi5,
0379                                                mp[1].tdc5,
0380                                                mp[1].lat5,
0381                                                mp[1].wi6,
0382                                                mp[1].tdc6,
0383                                                mp[1].lat6,
0384                                                mp[1].wi7,
0385                                                mp[1].tdc7,
0386                                                mp[1].lat7,
0387                                                mp[1].wi8,
0388                                                mp[1].tdc8,
0389                                                mp[1].lat8,
0390                                                -1}));
0391   }
0392   return;
0393 }
0394
0395 void MuonPathCorFitter::fillLuts() {
0396   std::ifstream ifin2sl(both_sl_filename_.fullPath());
0397   std::string line;
0398   while (ifin2sl.good()) {
0399     ifin2sl >> line;
0400
0401     std::vector<int> myNumbers;
0402     for (size_t i = 0; i < line.size(); i++) {
0403       // This converts the char into an int and pushes it into vec
0404       myNumbers.push_back(line[i] - '0');  // The digits will be in the same order as before
0405     }
0406     std::reverse(myNumbers.begin(), myNumbers.end());
0407     lut_2sl.push_back(myNumbers);
0408   }
0409
0410   return;
0411 }
0412
0413 int MuonPathCorFitter::get_rom_addr(mp_group mps, std::vector<int> missing_hits) {
0414   std::vector<int> lats = {
0415       mps[0].lat1, mps[0].lat2, mps[0].lat3, mps[0].lat4, mps[1].lat5, mps[1].lat6, mps[1].lat7, mps[1].lat8};
0416
0417   std::vector<int> rom_addr;
0418   if (missing_hits.size() == 1)
0419     rom_addr.push_back(1);
0420   else
0421     rom_addr.push_back(0);
0422
0423   if (missing_hits.size() == 1) {  // 7 layers fit
0424     if (missing_hits[0] < 4)
0425       rom_addr.push_back(0);  // First SL has 4 hits (1) or 3 (0)
0426     else
0427       rom_addr.push_back(1);
0428     if (missing_hits[0] % 4 == 0) {
0429       rom_addr.push_back(0);
0430       rom_addr.push_back(0);
0431     } else if (missing_hits[0] % 4 == 1) {
0432       rom_addr.push_back(0);
0433       rom_addr.push_back(1);
0434     } else if (missing_hits[0] % 4 == 2) {
0435       rom_addr.push_back(1);
0436       rom_addr.push_back(0);
0437     } else {  // missing_hits[0] == 3
0438       rom_addr.push_back(1);
0439       rom_addr.push_back(1);
0440     }
0441     for (size_t ilat = 0; ilat < lats.size(); ilat++) {
0442       if ((int)ilat == missing_hits[0])  // only applies to 3-hit, as in 4-hit missL=-1
0443         continue;
0444       auto lat = lats[ilat];
0445       if (lat == -1)
0446         lat = 0;
0447       rom_addr.push_back(lat);
0448     }
0449
0450   } else if (missing_hits.empty()) {  // 8 layers fit
0451     for (size_t ilat = 0; ilat < lats.size(); ilat++) {
0452       auto lat = lats[ilat];
0453       if (lat == -1)
0454         lat = 0;
0455       rom_addr.push_back(lat);
0456     }
0457     auto lat = lats[NUM_LAYERS + 3];
0458     if (lat == -1)
0459       lat = 0;
0460     rom_addr.push_back(lat);
0461     rom_addr.push_back(lat);
0462
0463   } else {  // 6 layers fit
0464     for (int i = missing_hits.size() - 1; i >= 0; i--) {
0465       if (missing_hits[i] % 4 == 0) {
0466         rom_addr.push_back(0);
0467         rom_addr.push_back(0);
0468       } else if (missing_hits[i] % 4 == 1) {
0469         rom_addr.push_back(0);
0470         rom_addr.push_back(1);
0471       } else if (missing_hits[i] % 4 == 2) {
0472         rom_addr.push_back(1);
0473         rom_addr.push_back(0);
0474       } else {  // missing_hits[i] % 4 == 3
0475         rom_addr.push_back(1);
0476         rom_addr.push_back(1);
0477       }
0478     }
0479     for (size_t ilat = 0; ilat < lats.size(); ilat++) {
0480       if ((int)ilat == missing_hits[0] || (int)ilat == missing_hits[1])  // only applies to 3-hit, as in 4-hit missL=-1
0481         continue;
0482       auto lat = lats[ilat];
0483       if (lat == -1)
0484         lat = 0;
0485       rom_addr.push_back(lat);
0486     }
0487   }
0488   std::reverse(rom_addr.begin(), rom_addr.end());
0489   return vhdl_unsigned_to_int(rom_addr);
0490 }