Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:54:21

0001 #include "L1Trigger/DTTriggerPhase2/interface/GlobalCoordsObtainer.h"
0002 
0003 using namespace edm;
0004 using namespace cmsdt;
0005 
0006 // ============================================================================
0007 // Constructors and destructor
0008 // ============================================================================
0009 GlobalCoordsObtainer::GlobalCoordsObtainer(const ParameterSet& pset) {
0010   global_coords_filename_ = pset.getParameter<edm::FileInPath>("global_coords_filename");
0011   std::ifstream ifin3(global_coords_filename_.fullPath());
0012 
0013   if (ifin3.fail()) {
0014     throw cms::Exception("Missing Input File")
0015         << "GlobalCoordsObtainer::GlobalCoordsObtainer() -  Cannot find " << global_coords_filename_.fullPath();
0016   }
0017 
0018   int wh, st, se, sl;
0019   double perp, x_phi0;
0020   std::string line;
0021 
0022   global_constant_per_sl sl1_constants;
0023   global_constant_per_sl sl3_constants;
0024 
0025   while (ifin3.good()) {
0026     ifin3 >> wh >> st >> se >> sl >> perp >> x_phi0;
0027 
0028     if (sl == 1) {
0029       sl1_constants.perp = perp;
0030       sl1_constants.x_phi0 = x_phi0;
0031     } else {
0032       sl3_constants.perp = perp;
0033       sl3_constants.x_phi0 = x_phi0;
0034 
0035       DTChamberId ChId(wh, st, se);
0036       global_constants.push_back({ChId.rawId(), sl1_constants, sl3_constants});
0037     }
0038   }
0039 }
0040 
0041 GlobalCoordsObtainer::~GlobalCoordsObtainer() {}
0042 
0043 void GlobalCoordsObtainer::generate_luts() {
0044   for (auto& global_constant : global_constants) {
0045     int sgn = 1;
0046     DTChamberId ChId(global_constant.chid);
0047     // typical hasPosRF function
0048     if (ChId.wheel() > 0 || (ChId.wheel() == 0 && ChId.sector() % 4 > 1)) {
0049       sgn = -1;
0050     }
0051 
0052     auto phi1 = calc_atan_lut(12,
0053                               6,
0054                               (1. / 16) / (global_constant.sl1.perp * 10),
0055                               global_constant.sl1.x_phi0 / global_constant.sl1.perp,
0056                               1. / std::pow(2, 17),
0057                               10,
0058                               3,
0059                               12,
0060                               20,
0061                               sgn);
0062 
0063     auto phi3 = calc_atan_lut(12,
0064                               6,
0065                               (1. / 16) / (global_constant.sl3.perp * 10),
0066                               global_constant.sl3.x_phi0 / global_constant.sl3.perp,
0067                               1. / std::pow(2, 17),
0068                               10,
0069                               3,
0070                               12,
0071                               20,
0072                               sgn);
0073 
0074     double max_x_phi0 = global_constant.sl1.x_phi0;
0075     if (global_constant.sl3.x_phi0 > max_x_phi0) {
0076       max_x_phi0 = global_constant.sl3.x_phi0;
0077     }
0078 
0079     auto phic = calc_atan_lut(12,
0080                               6,
0081                               (1. / 16) / ((global_constant.sl1.perp + global_constant.sl3.perp) / .2),
0082                               max_x_phi0 / ((global_constant.sl1.perp + global_constant.sl3.perp) / 2),
0083                               1. / std::pow(2, 17),
0084                               10,
0085                               3,
0086                               12,
0087                               20,
0088                               sgn);
0089 
0090     auto phib = calc_atan_lut(9, 6, 1. / 4096, 0., 4. / std::pow(2, 13), 10, 3, 10, 16, sgn);
0091 
0092     luts[global_constant.chid] = {phic, phi1, phi3, phib};
0093   }
0094 }
0095 
0096 std::map<int, lut_value> GlobalCoordsObtainer::calc_atan_lut(int msb_num,
0097                                                              int lsb_num,
0098                                                              double in_res,
0099                                                              double abscissa_0,
0100                                                              double out_res,
0101                                                              int a_extra_bits,
0102                                                              int b_extra_bits,
0103                                                              int a_size,
0104                                                              int b_size,
0105                                                              int sgn) {
0106   // Calculates the coefficients needed to calculate the arc-tan function in fw
0107   // by doing a piece-wise linear approximation.
0108   // In fw, input (x) and output (y) are integers, these conversions are needed
0109   //   t = x*in_res - abscissa_0
0110   //   phi = arctan(t)
0111   //   y = phi/out_res
0112   //   => y = arctan(x*in_res - abcissa_0)/out_res
0113   // The linear function is approximated as
0114   //   y = a*x_lsb + b
0115   // Where a, b = func(x_msb) are the coefficients stored in the lut
0116 
0117   // a is stored as unsigned, b as signed, with their respective sizes a_size, b_size,
0118   // previously shifted left by a_extra_bits and b_extra_bits, respectively
0119 
0120   long int a_min = -std::pow(2, (a_size - 1));
0121   long int a_max = std::pow(2, (a_size - 1)) - 1;
0122   long int b_min = -std::pow(2, (b_size - 1));
0123   long int b_max = std::pow(2, (b_size - 1)) - 1;
0124 
0125   std::map<int, lut_value> lut;
0126 
0127   for (long int x_msb = -(long int)std::pow(2, msb_num - 1); x_msb < (long int)std::pow(2, msb_num - 1); x_msb++) {
0128     int x1 = ((x_msb) << lsb_num);
0129     int x2 = ((x_msb + 1) << lsb_num) - 1;
0130 
0131     double t1 = x1 * in_res - abscissa_0;
0132     double t2 = x2 * in_res - abscissa_0;
0133 
0134     double phi1 = sgn * atan(t1);
0135     double phi2 = sgn * atan(t2);
0136 
0137     double y1 = phi1 / out_res;
0138     double y2 = phi2 / out_res;
0139 
0140     // we want to find a, b so that the error in the extremes is the same as the error in the center
0141     // so the error in the extremes will be the same, so the "a" is determined by those points
0142     double a = (y2 - y1) / (x2 - x1);
0143 
0144     // by derivating the error function and equaling to 0, you get this is the point
0145     // towards the interval center with the highest error
0146     // err_f = y - (a*x+b) = sgn*arctan(x*in_res - abcissa_0)/out_res - (a*x+b)
0147     // d(err_f)/dx = sgn*(1/(1+(x*in_res - abcissa_0)^2))*in_res/out_res - a
0148     // d(err_f)/dx = 0 => x_max_err = (sqrt(in_res/out_res/a-1) + abscissa_0)/in_res
0149     // There is sign ambiguity in the sqrt operation. The sqrt has units of t (adimensional).
0150     // It is resolved by setting the sqrt to have the same sign as (t1+t2)/2
0151 
0152     double t_max_err = sqrt(sgn * in_res / out_res / a - 1);
0153     if ((t1 + t2) < 0) {
0154       t_max_err *= -1;
0155     }
0156 
0157     double x_max_err = (t_max_err + abscissa_0) / in_res;
0158     double phi_max_err = sgn * atan(t_max_err);
0159     double y_max_err = phi_max_err / out_res;
0160 
0161     // once you have the point of max error, the "b" parameter is chosen as the average between
0162     // those two numbers, which makes the error at the center be equal in absolute value
0163     // to the error in the extremes
0164     // units: rad
0165 
0166     double b = (y1 + y_max_err - a * (x_max_err - x1)) / 2;
0167 
0168     // increase b in 1/2 of y_lsb, so that fw truncate operation on the of the result
0169     // is equivalent to a round function instead of a floor function
0170     b += 0.5;
0171 
0172     // shift left and round
0173     long int a_int = (long int)(round(a * (pow(2, a_extra_bits))));
0174     long int b_int = (long int)(round(b * (pow(2, b_extra_bits))));
0175 
0176     // tuck a, b constants into the bit size of the output (un)signed integer
0177     std::vector<long int> as = {a_min, a_int, a_max};
0178     std::vector<long int> bs = {b_min, b_int, b_max};
0179 
0180     std::sort(as.begin(), as.end());
0181     std::sort(bs.begin(), bs.end());
0182 
0183     a_int = as.at(1);
0184     b_int = bs.at(1);
0185 
0186     // // convert a, b to two's complement
0187     // auto a_signed = a_int % (long int) (pow(2, a_size));
0188     // auto b_signed = b_int % (long int) (pow(2, b_size));
0189 
0190     // convert x_msb to two's complement signed
0191     int index = to_two_comp(x_msb, msb_num);
0192     lut[index] = {a_int, b_int};
0193   }
0194   return lut;
0195 }
0196 
0197 std::vector<double> GlobalCoordsObtainer::get_global_coordinates(uint32_t chid, int sl, int x, int tanpsi) {
0198   // Depending on the type of primitive (SL1, SL3 or correlated), choose the
0199   // appropriate input data (x, tanpsi) from the input primitive data structure
0200   // and the corresponding phi-lut from the 3 available options
0201   auto phi_lut = &luts[chid].phic;
0202   if (sl == 1) {
0203     phi_lut = &luts[chid].phi1;
0204   } else if (sl == 3) {
0205     phi_lut = &luts[chid].phi3;
0206   }
0207 
0208   auto phib_lut = &luts[chid].phib;
0209 
0210   // x and slope are given in two's complement in fw
0211   x = to_two_comp(x, X_SIZE);
0212   tanpsi = to_two_comp(tanpsi, TANPSI_SIZE);
0213 
0214   // Slice x and tanpsi
0215   // Both x and tanpsi are represented in vhdl as signed, this means their values
0216   // are stored as two's complement.
0217 
0218   // The MSB part is going to be used to index the luts and obtain a and b parameters
0219   // Converting the upper part of the signed to an integer (with sign).
0220 
0221   int x_msb = x >> (X_SIZE - PHI_LUT_ADDR_WIDTH);
0222   x_msb = from_two_comp(x_msb, PHI_LUT_ADDR_WIDTH);
0223 
0224   int tanpsi_msb = tanpsi >> (TANPSI_SIZE - PHIB_LUT_ADDR_WIDTH);
0225   tanpsi_msb = from_two_comp(tanpsi_msb, PHIB_LUT_ADDR_WIDTH);
0226 
0227   x_msb = x >> (X_SIZE - PHI_LUT_ADDR_WIDTH);
0228   x_msb = from_two_comp(x_msb, PHI_LUT_ADDR_WIDTH);
0229 
0230   tanpsi_msb = tanpsi >> (TANPSI_SIZE - PHIB_LUT_ADDR_WIDTH);
0231   tanpsi_msb = from_two_comp(tanpsi_msb, PHIB_LUT_ADDR_WIDTH);
0232 
0233   // The LSB part can be sliced right away because it must yield a positive integer
0234   int x_lsb = x & (int)(std::pow(2, (X_SIZE - PHI_LUT_ADDR_WIDTH)) - 1);
0235   int tanpsi_lsb = tanpsi & (int)(std::pow(2, (TANPSI_SIZE - PHIB_LUT_ADDR_WIDTH)) - 1);
0236 
0237   // Index the luts wiht the MSB parts already calculated
0238   auto phi_lut_q = phi_lut->at(to_two_comp(x_msb, PHI_LUT_ADDR_WIDTH));
0239   auto phib_lut_q = phib_lut->at(to_two_comp(tanpsi_msb, PHIB_LUT_ADDR_WIDTH));
0240 
0241   // Separate this data into the coefficients a and b
0242   auto phi_lut_a = phi_lut_q.a;
0243   auto phi_lut_b = phi_lut_q.b;
0244   auto phib_lut_a = phib_lut_q.a;
0245   auto phib_lut_b = phib_lut_q.b;
0246 
0247   // Do the linear piece-wise operations
0248   // At this point all variables that can be negative have already been converted
0249   // so will yield negative values when necessary
0250   int phi_uncut = (phi_lut_b << PHI_B_SHL_BITS) + x_lsb * phi_lut_a;
0251   int psi_uncut = (phib_lut_b << PHIB_B_SHL_BITS) + tanpsi_lsb * phib_lut_a;
0252 
0253   // Trim phi to its final size
0254   int phi = (phi_uncut >> PHI_MULT_SHR_BITS);
0255 
0256   // Calculate phi_bending from the uncut version of phi and psi, and the trim it to size
0257   int phib_uncut = psi_uncut - (phi_uncut >> (PHI_PHIB_RES_DIFF_BITS + PHI_MULT_SHR_BITS - PHIB_MULT_SHR_BITS));
0258   int phib = (phib_uncut >> PHIB_MULT_SHR_BITS);
0259 
0260   double phi_f = (double)phi / pow(2, PHI_SIZE);
0261   double phib_f = (double)phib / pow(2, PHIB_SIZE);
0262 
0263   return std::vector({phi_f, phib_f});
0264 }