DTTriggerPhase2/src/MuonPathFitter.cc

0001 #include "L1Trigger/DTTriggerPhase2/interface/MuonPathFitter.h"
0002 #include <cmath>
0003 #include <memory>
0004
0005 using namespace edm;
0006 using namespace std;
0007 using namespace cmsdt;
0008
0009 // ============================================================================
0010 // Constructors and destructor
0011 // ============================================================================
0012 MuonPathFitter::MuonPathFitter(const ParameterSet &pset,
0013                                edm::ConsumesCollector &iC,
0014                                std::shared_ptr<GlobalCoordsObtainer> &globalcoordsobtainer)
0015     : MuonPathAnalyzer(pset, iC), debug_(pset.getUntrackedParameter<bool>("debug")) {
0016   if (debug_)
0017     LogDebug("MuonPathFitter") << "MuonPathAnalyzer: constructor";
0018
0019   //shift phi
0020   int rawId;
0021   shift_filename_ = pset.getParameter<edm::FileInPath>("shift_filename");
0022   std::ifstream ifin3(shift_filename_.fullPath());
0023   double shift;
0024   if (ifin3.fail()) {
0025     throw cms::Exception("Missing Input File")
0026         << "MuonPathFitter::MuonPathFitter() -  Cannot find " << shift_filename_.fullPath();
0027   }
0028   while (ifin3.good()) {
0029     ifin3 >> rawId >> shift;
0030     shiftinfo_[rawId] = shift;
0031   }
0032
0033   int wh, st, se, maxdrift;
0034   maxdrift_filename_ = pset.getParameter<edm::FileInPath>("maxdrift_filename");
0035   std::ifstream ifind(maxdrift_filename_.fullPath());
0036   if (ifind.fail()) {
0037     throw cms::Exception("Missing Input File")
0038         << "MPSLFilter::MPSLFilter() -  Cannot find " << maxdrift_filename_.fullPath();
0039   }
0040   while (ifind.good()) {
0041     ifind >> wh >> st >> se >> maxdrift;
0042     maxdriftinfo_[wh][st][se] = maxdrift;
0043   }
0044
0045   dtGeomH = iC.esConsumes<DTGeometry, MuonGeometryRecord, edm::Transition::BeginRun>();
0046   globalcoordsobtainer_ = globalcoordsobtainer;
0047 }
0048
0049 MuonPathFitter::~MuonPathFitter() {
0050   if (debug_)
0051     LogDebug("MuonPathFitter") << "MuonPathAnalyzer: destructor";
0052 }
0053
0054 // ============================================================================
0055 // Main methods (initialise, run, finish)
0056 // ============================================================================
0057
0058 //------------------------------------------------------------------
0059 //--- Metodos privados
0060 //------------------------------------------------------------------
0061
0062 fit_common_out_t MuonPathFitter::fit(fit_common_in_t fit_common_in,
0063                                      int XI_WIDTH,
0064                                      int COEFF_WIDTH_T0,
0065                                      int COEFF_WIDTH_POSITION,
0066                                      int COEFF_WIDTH_SLOPE,
0067                                      int PRECISSION_T0,
0068                                      int PRECISSION_POSITION,
0069                                      int PRECISSION_SLOPE,
0070                                      int PROD_RESIZE_T0,
0071                                      int PROD_RESIZE_POSITION,
0072                                      int PROD_RESIZE_SLOPE,
0073                                      int MAX_DRIFT_TDC,
0074                                      int sl) {
0075   const int PARTIALS_PRECISSION = 4;
0076   const int PARTIALS_SHR_T0 = PRECISSION_T0 - PARTIALS_PRECISSION;
0077   const int PARTIALS_SHR_POSITION = PRECISSION_POSITION - PARTIALS_PRECISSION;
0078   const int PARTIALS_SHR_SLOPE = PRECISSION_SLOPE - PARTIALS_PRECISSION;
0079   const int PARTIALS_WIDTH_T0 = PROD_RESIZE_T0 - PARTIALS_SHR_T0;
0080   const int PARTIALS_WIDTH_POSITION = PROD_RESIZE_POSITION - PARTIALS_SHR_POSITION;
0081   const int PARTIALS_WIDTH_SLOPE = PROD_RESIZE_SLOPE - PARTIALS_SHR_SLOPE;
0082
0083   const int WIDTH_TO_PREC = 11 + PARTIALS_PRECISSION;
0084   const int WIDTH_SLOPE_PREC = 14 + PARTIALS_PRECISSION;
0085   const int WIDTH_POSITION_PREC = WIDTH_SLOPE_PREC + 1;
0086
0087   const int SEMICHAMBER_H_PRECISSION = 13 + PARTIALS_PRECISSION;
0088   const float SEMICHAMBER_H_REAL = ((235. / 2.) / (16. * 6.5)) * std::pow(2, SEMICHAMBER_H_PRECISSION);
0089   const int SEMICHAMBER_H = (int)SEMICHAMBER_H_REAL;  // signed(SEMICHAMBER_H_WIDTH-1 downto 0)
0090
0091   const int SEMICHAMBER_RES_SHR = SEMICHAMBER_H_PRECISSION;
0092
0093   const int LYRANDAHALF_RES_SHR = 4;
0094
0095   const int CHI2_CALC_RES_BITS = 7;
0096
0097   /*******************************
0098             clock cycle 1
0099   *******************************/
0100   std::vector<int> normalized_times;
0101   std::vector<int> normalized_wirepos;
0102
0103   for (int i = 0; i < 2 * NUM_LAYERS; i++) {
0104     // normalized times
0105     // this should be resized to an unsigned of 10 bits (max drift time ~508 TDC counts, using 9+1 to include tolerance)
0106     // leaving it as an integer for now
0107     // we are obtaining the difference as the difference in BX + the LS bits from the hit time
0108
0109     if (fit_common_in.hits_valid[i] == 1) {
0110       int dif_bx = (fit_common_in.hits[i].ti >> (WIDTH_FULL_TIME - WIDTH_COARSED_TIME)) - fit_common_in.coarse_bctr;
0111
0112       int tmp_norm_time = (dif_bx << (WIDTH_FULL_TIME - WIDTH_COARSED_TIME)) +
0113                           (fit_common_in.hits[i].ti % (int)std::pow(2, WIDTH_FULL_TIME - WIDTH_COARSED_TIME));
0114       // resize test
0115       // this has implications in the FW (reducing number of bits).
0116       // we keep here the int as it is, but we do the same check done in the fw
0117       std::vector<int> tmp_dif_bx_vector;
0118       vhdl_int_to_unsigned(dif_bx, tmp_dif_bx_vector);
0119       vhdl_resize_unsigned(tmp_dif_bx_vector, 12);
0120       if (!vhdl_resize_unsigned_ok(tmp_dif_bx_vector, WIDTH_DIFBX))
0121         return fit_common_out_t();
0122
0123       normalized_times.push_back(tmp_norm_time);
0124       int tmp_wirepos = fit_common_in.hits[i].wp - (fit_common_in.coarse_wirepos << WIREPOS_NORM_LSB_IGNORED);
0125       // resize test
0126       std::vector<int> tmp_wirepos_vector;
0127       vhdl_int_to_signed(tmp_wirepos, tmp_wirepos_vector);
0128       vhdl_resize_signed(tmp_wirepos_vector, WIREPOS_WIDTH);
0129       if (!vhdl_resize_signed_ok(tmp_wirepos_vector, XI_WIDTH))
0130         return fit_common_out_t();
0131
0132       normalized_wirepos.push_back(tmp_wirepos);
0133     } else {  // dummy hit
0134       normalized_times.push_back(-1);
0135       normalized_wirepos.push_back(-1);
0136     }
0137   }
0138
0139   /*******************************
0140             clock cycle 2
0141   *******************************/
0142
0143   std::vector<int> xi_arr;
0144   // min and max times are computed throught several clk cycles in the fw,
0145   // here we compute it at once
0146   int min_hit_time = 999999, max_hit_time = 0;
0147   for (int i = 0; i < 2 * NUM_LAYERS; i++) {
0148     if (fit_common_in.hits_valid[i] == 1) {
0149       // calculate xi array
0150       auto tmp_xi_incr = normalized_wirepos[i];
0151       tmp_xi_incr += (-1 + 2 * fit_common_in.lateralities[i]) * normalized_times[i];
0152
0153       // resize test
0154       std::vector<int> tmp_xi_incr_vector;
0155       vhdl_int_to_signed(tmp_xi_incr, tmp_xi_incr_vector);
0156       vhdl_resize_signed(tmp_xi_incr_vector, XI_WIDTH + 1);
0157       if (!vhdl_resize_signed_ok(tmp_xi_incr_vector, XI_WIDTH))
0158         return fit_common_out_t();
0159       xi_arr.push_back(tmp_xi_incr);
0160
0161       // calculate min and max times
0162       if (normalized_times[i] < min_hit_time) {
0163         min_hit_time = normalized_times[i];
0164       }
0165       if (normalized_times[i] > max_hit_time) {
0166         max_hit_time = normalized_times[i];
0167       }
0168     } else {
0169       xi_arr.push_back(-1);
0170     }
0171   }
0172
0173   /*******************************
0174             clock cycle 3
0175   *******************************/
0176
0177   std::vector<int> products_t0;
0178   std::vector<int> products_position;
0179   std::vector<int> products_slope;
0180   for (int i = 0; i < 2 * NUM_LAYERS; i++) {
0181     if (fit_common_in.hits_valid[i] == 0) {
0182       products_t0.push_back(-1);
0183       products_position.push_back(-1);
0184       products_slope.push_back(-1);
0185     } else {
0186       products_t0.push_back(xi_arr[i] * vhdl_signed_to_int(fit_common_in.coeffs.t0[i]));
0187       products_position.push_back(xi_arr[i] * vhdl_signed_to_int(fit_common_in.coeffs.position[i]));
0188       products_slope.push_back(xi_arr[i] * vhdl_signed_to_int(fit_common_in.coeffs.slope[i]));
0189     }
0190   }
0191
0192   /*******************************
0193             clock cycle 4
0194   *******************************/
0195   // Do the 8 element sums
0196   int t0_prec = 0, position_prec = 0, slope_prec = 0;
0197   for (int i = 0; i < 2 * NUM_LAYERS; i++) {
0198     if (fit_common_in.hits_valid[i] == 0) {
0199       continue;
0200     } else {
0201       t0_prec += products_t0[i] >> PARTIALS_SHR_T0;
0202       position_prec += products_position[i] >> PARTIALS_SHR_POSITION;
0203       slope_prec += products_slope[i] >> PARTIALS_SHR_SLOPE;
0204     }
0205   }
0206
0207   /*******************************
0208             clock cycle 5
0209   *******************************/
0210   // Do resize tests for the computed sums with full precision
0211   std::vector<int> t0_prec_vector, position_prec_vector, slope_prec_vector;
0212   vhdl_int_to_signed(t0_prec, t0_prec_vector);
0213
0214   vhdl_resize_signed(t0_prec_vector, PARTIALS_WIDTH_T0);
0215   if (!vhdl_resize_signed_ok(t0_prec_vector, WIDTH_TO_PREC))
0216     return fit_common_out_t();
0217
0218   vhdl_int_to_signed(position_prec, position_prec_vector);
0219   vhdl_resize_signed(position_prec_vector, PARTIALS_WIDTH_POSITION);
0220   if (!vhdl_resize_signed_ok(position_prec_vector, WIDTH_POSITION_PREC))
0221     return fit_common_out_t();
0222
0223   vhdl_int_to_signed(slope_prec, slope_prec_vector);
0224   vhdl_resize_signed(slope_prec_vector, PARTIALS_WIDTH_SLOPE);
0225   if (!vhdl_resize_signed_ok(slope_prec_vector, WIDTH_SLOPE_PREC))
0226     return fit_common_out_t();
0227
0228   /*******************************
0229             clock cycle 6
0230   *******************************/
0231   // Round the fitting parameters to the final resolution;
0232   // in vhdl something more sofisticated is done, here we do a float division, round
0233   // and cast again to integer
0234
0235   int norm_t0 = ((t0_prec >> (PARTIALS_PRECISSION - 1)) + 1) >> 1;
0236   int norm_position = ((position_prec >> (PARTIALS_PRECISSION - 1)) + 1) >> 1;
0237   int norm_slope = ((slope_prec >> (PARTIALS_PRECISSION - 1)) + 1) >> 1;
0238
0239   // Calculate the (-xi) + pos (+/-) t0, which only is lacking the slope term to become the residuals
0240   std::vector<int> res_partials_arr;
0241   for (int i = 0; i < 2 * NUM_LAYERS; i++) {
0242     if (fit_common_in.hits_valid[i] == 0) {
0243       res_partials_arr.push_back(-1);
0244     } else {
0245       int tmp_position_prec = position_prec - (xi_arr[i] << PARTIALS_PRECISSION);
0246       // rounding
0247       tmp_position_prec += std::pow(2, PARTIALS_PRECISSION - 1);
0248
0249       tmp_position_prec += (-1 + 2 * fit_common_in.lateralities[i]) * t0_prec;
0250       res_partials_arr.push_back(tmp_position_prec);
0251     }
0252   }
0253
0254   // calculate the { slope x semichamber, slope x 1.5 layers, slope x 0.5 layers }
0255   // these 3 values are later combined with different signs to get the slope part
0256   // of the residual for each of the layers.
0257   int slope_x_halfchamb = (((long int)slope_prec * (long int)SEMICHAMBER_H)) >> SEMICHAMBER_RES_SHR;
0258   if (sl == 2)
0259     slope_x_halfchamb = 0;
0260   int slope_x_3semicells = (slope_prec * 3) >> LYRANDAHALF_RES_SHR;
0261   int slope_x_1semicell = (slope_prec * 1) >> LYRANDAHALF_RES_SHR;
0262
0263   /*******************************
0264             clock cycle 7
0265   *******************************/
0266   // Complete the residuals calculation by constructing the slope term (1/2)
0267   for (int i = 0; i < 2 * NUM_LAYERS; i++) {
0268     if (fit_common_in.hits_valid[i] == 1) {
0269       if (i % 4 == 0)
0270         res_partials_arr[i] -= slope_x_3semicells;
0271       else if (i % 4 == 1)
0272         res_partials_arr[i] -= slope_x_1semicell;
0273       else if (i % 4 == 2)
0274         res_partials_arr[i] += slope_x_1semicell;
0275       else
0276         res_partials_arr[i] += slope_x_3semicells;
0277     }
0278   }
0279
0280   /*******************************
0281             clock cycle 8
0282   *******************************/
0283   // Complete the residuals calculation by constructing the slope term (2/2)
0284   std::vector<int> residuals, position_prec_arr;
0285   for (int i = 0; i < 2 * NUM_LAYERS; i++) {
0286     if (fit_common_in.hits_valid[i] == 0) {
0287       residuals.push_back(-1);
0288       position_prec_arr.push_back(-1);
0289     } else {
0290       int tmp_position_prec = res_partials_arr[i];
0291       tmp_position_prec += (-1 + 2 * (int)(i >= NUM_LAYERS)) * slope_x_halfchamb;
0292       position_prec_arr.push_back(tmp_position_prec);
0293       residuals.push_back(abs(tmp_position_prec >> PARTIALS_PRECISSION));
0294     }
0295   }
0296
0297   // minimum and maximum fit t0
0298   int min_t0 = max_hit_time - MAX_DRIFT_TDC - T0_CUT_TOLERANCE;
0299   int max_t0 = min_hit_time + T0_CUT_TOLERANCE;
0300
0301   /*******************************
0302             clock cycle 9
0303   *******************************/
0304   // Prepare addition of coarse_offset to T0 (T0 de-normalization)
0305   int t0_fine = norm_t0 & (int)(std::pow(2, 5) - 1);
0306   int t0_bx_sign = ((int)(norm_t0 < 0)) * 1;
0307   int t0_bx_abs = abs(norm_t0 >> 5);
0308
0309   // De-normalize Position and slope
0310   int position = (fit_common_in.coarse_wirepos << WIREPOS_NORM_LSB_IGNORED) + norm_position;
0311   int slope = norm_slope;
0312
0313   // Apply T0 cuts
0314   if (norm_t0 < min_t0)
0315     return fit_common_out_t();
0316   if (norm_t0 > max_t0)
0317     return fit_common_out_t();
0318
0319   // square the residuals
0320   std::vector<int> squared_residuals;
0321   for (int i = 0; i < 2 * NUM_LAYERS; i++) {
0322     if (fit_common_in.hits_valid[i] == 0) {
0323       squared_residuals.push_back(-1);
0324     } else {
0325       squared_residuals.push_back(residuals[i] * residuals[i]);
0326     }
0327   }
0328
0329   // check for residuals overflow
0330   for (int i = 0; i < 2 * NUM_LAYERS; i++) {
0331     if (fit_common_in.hits_valid[i] == 1) {
0332       std::vector<int> tmp_vector;
0333       int tmp_position_prec = (position_prec_arr[i] >> PARTIALS_PRECISSION);
0334       vhdl_int_to_signed(tmp_position_prec, tmp_vector);
0335       vhdl_resize_signed(tmp_vector, WIDTH_POSITION_PREC);
0336       if (!vhdl_resize_signed_ok(tmp_vector, CHI2_CALC_RES_BITS + 1))
0337         return fit_common_out_t();
0338       // Commented for now, maybe later we need to do something here
0339       // if ((tmp_position_prec / (int) std::pow(2, CHI2_CALC_RES_BITS)) > 0)
0340       // return fit_common_out_t();
0341     }
0342   }
0343
0344   /*******************************
0345         clock cycle 10, 11, 12
0346   *******************************/
0347   int t0 = t0_fine;
0348   t0 += (fit_common_in.coarse_bctr - (-1 + 2 * t0_bx_sign) * t0_bx_abs) * (int)std::pow(2, 5);
0349
0350   int chi2 = 0;
0351   for (int i = 0; i < 2 * NUM_LAYERS; i++) {
0352     if (fit_common_in.hits_valid[i] == 1) {
0353       chi2 += squared_residuals[i];
0354     }
0355   }
0356
0357   // Impose the thresholds
0358
0359   if (chi2 / 16 >= (int)round(chi2Th_ * (std::pow((float)MAX_DRIFT_TDC / ((float)CELL_SEMILENGTH / 10.), 2)) / 16.))
0360     return fit_common_out_t();
0361
0362   fit_common_out_t fit_common_out;
0363   fit_common_out.position = position;
0364   fit_common_out.slope = slope;
0365   fit_common_out.t0 = t0;
0366   fit_common_out.chi2 = chi2;
0367   fit_common_out.valid_fit = 1;
0368
0369   return fit_common_out;
0370 }
0371
0372 coeffs_t MuonPathFitter::RomDataConvert(std::vector<int> slv,
0373                                         short COEFF_WIDTH_T0,
0374                                         short COEFF_WIDTH_POSITION,
0375                                         short COEFF_WIDTH_SLOPE,
0376                                         short LOLY,
0377                                         short HILY) {
0378   coeffs_t res;
0379   int ctr = 0;
0380   for (int i = LOLY; i <= HILY; i++) {
0381     res.t0[i] = vhdl_slice(slv, COEFF_WIDTH_T0 + ctr - 1, ctr);
0382     vhdl_resize_unsigned(res.t0[i], GENERIC_COEFF_WIDTH);
0383     res.t0[i] = vhdl_slice(res.t0[i], COEFF_WIDTH_T0 - 1, 0);
0384     ctr += COEFF_WIDTH_T0;
0385   }
0386   for (int i = LOLY; i <= HILY; i++) {
0387     res.position[i] = vhdl_slice(slv, COEFF_WIDTH_POSITION + ctr - 1, ctr);
0388     vhdl_resize_unsigned(res.position[i], GENERIC_COEFF_WIDTH);
0389     res.position[i] = vhdl_slice(res.position[i], COEFF_WIDTH_POSITION - 1, 0);
0390     ctr += COEFF_WIDTH_POSITION;
0391   }
0392   for (int i = LOLY; i <= HILY; i++) {
0393     res.slope[i] = vhdl_slice(slv, COEFF_WIDTH_SLOPE + ctr - 1, ctr);
0394     vhdl_resize_unsigned(res.slope[i], GENERIC_COEFF_WIDTH);
0395     res.slope[i] = vhdl_slice(res.slope[i], COEFF_WIDTH_SLOPE - 1, 0);
0396     ctr += COEFF_WIDTH_SLOPE;
0397   }
0398   return res;
0399 }