CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
GlobalCoordsObtainer Class Reference

#include <GlobalCoordsObtainer.h>

Public Member Functions

void generate_luts ()
 
std::vector< double > get_global_coordinates (uint32_t, int, int, int)
 
 GlobalCoordsObtainer (const edm::ParameterSet &pset)
 
 ~GlobalCoordsObtainer ()
 

Private Member Functions

std::map< int, lut_valuecalc_atan_lut (int, int, double, double, double, int, int, int, int, int)
 
int from_two_comp (int val, int size)
 
int to_two_comp (int val, int size)
 

Private Attributes

bool cmssw_for_global_
 
std::vector< global_constantglobal_constants
 
edm::FileInPath global_coords_filename_
 
std::map< uint32_t, lut_groupluts
 

Detailed Description

Definition at line 54 of file GlobalCoordsObtainer.h.

Constructor & Destructor Documentation

GlobalCoordsObtainer::GlobalCoordsObtainer ( const edm::ParameterSet pset)

Definition at line 9 of file GlobalCoordsObtainer.cc.

References Exception, edm::ParameterSet::getParameter(), geometryCSVtoXML::line, global_constant_per_sl::perp, perp(), DetId::rawId(), AlCaHLTBitMon_QueryRunRegistry::string, and global_constant_per_sl::x_phi0.

9  {
10  global_coords_filename_ = pset.getParameter<edm::FileInPath>("global_coords_filename");
11  std::ifstream ifin3(global_coords_filename_.fullPath());
12 
13  if (ifin3.fail()) {
14  throw cms::Exception("Missing Input File")
15  << "GlobalCoordsObtainer::GlobalCoordsObtainer() - Cannot find " << global_coords_filename_.fullPath();
16  }
17 
18  int wh, st, se, sl;
19  double perp, x_phi0;
21 
22  global_constant_per_sl sl1_constants;
23  global_constant_per_sl sl3_constants;
24 
25  while (ifin3.good()) {
26  ifin3 >> wh >> st >> se >> sl >> perp >> x_phi0;
27 
28  if (sl == 1) {
29  sl1_constants.perp = perp;
30  sl1_constants.x_phi0 = x_phi0;
31  } else {
32  sl3_constants.perp = perp;
33  sl3_constants.x_phi0 = x_phi0;
34 
35  DTChamberId ChId(wh, st, se);
36  global_constants.push_back({ChId.rawId(), sl1_constants, sl3_constants});
37  }
38  }
39 }
edm::FileInPath global_coords_filename_
std::vector< global_constant > global_constants
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
T perp() const
Magnitude of transverse component.
std::string fullPath() const
Definition: FileInPath.cc:161
GlobalCoordsObtainer::~GlobalCoordsObtainer ( )

Definition at line 41 of file GlobalCoordsObtainer.cc.

41 {}

Member Function Documentation

std::map< int, lut_value > GlobalCoordsObtainer::calc_atan_lut ( int  msb_num,
int  lsb_num,
double  in_res,
double  abscissa_0,
double  out_res,
int  a_extra_bits,
int  b_extra_bits,
int  a_size,
int  b_size,
int  sgn 
)
private

Definition at line 96 of file GlobalCoordsObtainer.cc.

References a, b, cms::cuda::bs, funct::pow(), and mathSSE::sqrt().

105  {
106  // Calculates the coefficients needed to calculate the arc-tan function in fw
107  // by doing a piece-wise linear approximation.
108  // In fw, input (x) and output (y) are integers, these conversions are needed
109  // t = x*in_res - abscissa_0
110  // phi = arctan(t)
111  // y = phi/out_res
112  // => y = arctan(x*in_res - abcissa_0)/out_res
113  // The linear function is approximated as
114  // y = a*x_lsb + b
115  // Where a, b = func(x_msb) are the coefficients stored in the lut
116 
117  // a is stored as unsigned, b as signed, with their respective sizes a_size, b_size,
118  // previously shifted left by a_extra_bits and b_extra_bits, respectively
119 
120  long int a_min = -std::pow(2, (a_size - 1));
121  long int a_max = std::pow(2, (a_size - 1)) - 1;
122  long int b_min = -std::pow(2, (b_size - 1));
123  long int b_max = std::pow(2, (b_size - 1)) - 1;
124 
125  std::map<int, lut_value> lut;
126 
127  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++) {
128  int x1 = ((x_msb) << lsb_num);
129  int x2 = ((x_msb + 1) << lsb_num) - 1;
130 
131  double t1 = x1 * in_res - abscissa_0;
132  double t2 = x2 * in_res - abscissa_0;
133 
134  double phi1 = sgn * atan(t1);
135  double phi2 = sgn * atan(t2);
136 
137  double y1 = phi1 / out_res;
138  double y2 = phi2 / out_res;
139 
140  // we want to find a, b so that the error in the extremes is the same as the error in the center
141  // so the error in the extremes will be the same, so the "a" is determined by those points
142  double a = (y2 - y1) / (x2 - x1);
143 
144  // by derivating the error function and equaling to 0, you get this is the point
145  // towards the interval center with the highest error
146  // err_f = y - (a*x+b) = sgn*arctan(x*in_res - abcissa_0)/out_res - (a*x+b)
147  // d(err_f)/dx = sgn*(1/(1+(x*in_res - abcissa_0)^2))*in_res/out_res - a
148  // d(err_f)/dx = 0 => x_max_err = (sqrt(in_res/out_res/a-1) + abscissa_0)/in_res
149  // There is sign ambiguity in the sqrt operation. The sqrt has units of t (adimensional).
150  // It is resolved by setting the sqrt to have the same sign as (t1+t2)/2
151 
152  double t_max_err = sqrt(sgn * in_res / out_res / a - 1);
153  if ((t1 + t2) < 0) {
154  t_max_err *= -1;
155  }
156 
157  double x_max_err = (t_max_err + abscissa_0) / in_res;
158  double phi_max_err = sgn * atan(t_max_err);
159  double y_max_err = phi_max_err / out_res;
160 
161  // once you have the point of max error, the "b" parameter is chosen as the average between
162  // those two numbers, which makes the error at the center be equal in absolute value
163  // to the error in the extremes
164  // units: rad
165 
166  double b = (y1 + y_max_err - a * (x_max_err - x1)) / 2;
167 
168  // increase b in 1/2 of y_lsb, so that fw truncate operation on the of the result
169  // is equivalent to a round function instead of a floor function
170  b += 0.5;
171 
172  // shift left and round
173  long int a_int = (long int)(round(a * (pow(2, a_extra_bits))));
174  long int b_int = (long int)(round(b * (pow(2, b_extra_bits))));
175 
176  // tuck a, b constants into the bit size of the output (un)signed integer
177  std::vector<long int> as = {a_min, a_int, a_max};
178  std::vector<long int> bs = {b_min, b_int, b_max};
179 
180  std::sort(as.begin(), as.end());
181  std::sort(bs.begin(), bs.end());
182 
183  a_int = as.at(1);
184  b_int = bs.at(1);
185 
186  // // convert a, b to two's complement
187  // auto a_signed = a_int % (long int) (pow(2, a_size));
188  // auto b_signed = b_int % (long int) (pow(2, b_size));
189 
190  // convert x_msb to two's complement signed
191  int index = to_two_comp(x_msb, msb_num);
192  lut[index] = {a_int, b_int};
193  }
194  return lut;
195 }
float sgn(float val)
Definition: FWPFMaths.cc:9
int to_two_comp(int val, int size)
T sqrt(T t)
Definition: SSEVec.h:19
double b
Definition: hdecay.h:118
double a
Definition: hdecay.h:119
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
int GlobalCoordsObtainer::from_two_comp ( int  val,
int  size 
)
inlineprivate

Definition at line 71 of file GlobalCoordsObtainer.h.

References hgcalPerformanceValidation::val.

71 { return val - ((2 * val) & (1 << size)); }
tuple size
Write out results.
void GlobalCoordsObtainer::generate_luts ( )

Definition at line 43 of file GlobalCoordsObtainer.cc.

References global_constant::chid, global_constant_per_sl::perp, funct::pow(), DTChamberId::sector(), FWPFMaths::sgn(), global_constant::sl1, global_constant::sl3, DTChamberId::wheel(), and global_constant_per_sl::x_phi0.

43  {
44  for (auto& global_constant : global_constants) {
45  int sgn = 1;
47  // typical hasPosRF function
48  if (ChId.wheel() > 0 || (ChId.wheel() == 0 && ChId.sector() % 4 > 1)) {
49  sgn = -1;
50  }
51 
52  auto phi1 = calc_atan_lut(12,
53  6,
54  (1. / 16) / (global_constant.sl1.perp * 10),
56  1. / std::pow(2, 17),
57  10,
58  3,
59  12,
60  20,
61  sgn);
62 
63  auto phi3 = calc_atan_lut(12,
64  6,
65  (1. / 16) / (global_constant.sl3.perp * 10),
67  1. / std::pow(2, 17),
68  10,
69  3,
70  12,
71  20,
72  sgn);
73 
74  double max_x_phi0 = global_constant.sl1.x_phi0;
75  if (global_constant.sl3.x_phi0 > max_x_phi0) {
76  max_x_phi0 = global_constant.sl3.x_phi0;
77  }
78 
79  auto phic = calc_atan_lut(12,
80  6,
81  (1. / 16) / ((global_constant.sl1.perp + global_constant.sl3.perp) / .2),
82  max_x_phi0 / ((global_constant.sl1.perp + global_constant.sl3.perp) / 2),
83  1. / std::pow(2, 17),
84  10,
85  3,
86  12,
87  20,
88  sgn);
89 
90  auto phib = calc_atan_lut(9, 6, 1. / 4096, 0., 4. / std::pow(2, 13), 10, 3, 10, 16, sgn);
91 
92  luts[global_constant.chid] = {phic, phi1, phi3, phib};
93  }
94 }
float sgn(float val)
Definition: FWPFMaths.cc:9
global_constant_per_sl sl1
std::vector< global_constant > global_constants
std::map< uint32_t, lut_group > luts
global_constant_per_sl sl3
std::map< int, lut_value > calc_atan_lut(int, int, double, double, double, int, int, int, int, int)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
std::vector< double > GlobalCoordsObtainer::get_global_coordinates ( uint32_t  chid,
int  sl,
int  x,
int  tanpsi 
)

Definition at line 197 of file GlobalCoordsObtainer.cc.

References cmsdt::PHI_B_SHL_BITS, cmsdt::PHI_LUT_ADDR_WIDTH, cmsdt::PHI_MULT_SHR_BITS, cmsdt::PHI_PHIB_RES_DIFF_BITS, cmsdt::PHI_SIZE, cmsdt::PHIB_B_SHL_BITS, cmsdt::PHIB_LUT_ADDR_WIDTH, cmsdt::PHIB_MULT_SHR_BITS, cmsdt::PHIB_SIZE, funct::pow(), cmsdt::TANPSI_SIZE, trackerHitRTTI::vector, and cmsdt::X_SIZE.

197  {
198  // Depending on the type of primitive (SL1, SL3 or correlated), choose the
199  // appropriate input data (x, tanpsi) from the input primitive data structure
200  // and the corresponding phi-lut from the 3 available options
201  auto phi_lut = &luts[chid].phic;
202  if (sl == 1) {
203  phi_lut = &luts[chid].phi1;
204  } else if (sl == 3) {
205  phi_lut = &luts[chid].phi3;
206  }
207 
208  auto phib_lut = &luts[chid].phib;
209 
210  // x and slope are given in two's complement in fw
211  x = to_two_comp(x, X_SIZE);
212  tanpsi = to_two_comp(tanpsi, TANPSI_SIZE);
213 
214  // Slice x and tanpsi
215  // Both x and tanpsi are represented in vhdl as signed, this means their values
216  // are stored as two's complement.
217 
218  // The MSB part is going to be used to index the luts and obtain a and b parameters
219  // Converting the upper part of the signed to an integer (with sign).
220 
221  int x_msb = x >> (X_SIZE - PHI_LUT_ADDR_WIDTH);
222  x_msb = from_two_comp(x_msb, PHI_LUT_ADDR_WIDTH);
223 
224  int tanpsi_msb = tanpsi >> (TANPSI_SIZE - PHIB_LUT_ADDR_WIDTH);
225  tanpsi_msb = from_two_comp(tanpsi_msb, PHIB_LUT_ADDR_WIDTH);
226 
227  x_msb = x >> (X_SIZE - PHI_LUT_ADDR_WIDTH);
228  x_msb = from_two_comp(x_msb, PHI_LUT_ADDR_WIDTH);
229 
230  tanpsi_msb = tanpsi >> (TANPSI_SIZE - PHIB_LUT_ADDR_WIDTH);
231  tanpsi_msb = from_two_comp(tanpsi_msb, PHIB_LUT_ADDR_WIDTH);
232 
233  // The LSB part can be sliced right away because it must yield a positive integer
234  int x_lsb = x & (int)(std::pow(2, (X_SIZE - PHI_LUT_ADDR_WIDTH)) - 1);
235  int tanpsi_lsb = tanpsi & (int)(std::pow(2, (TANPSI_SIZE - PHIB_LUT_ADDR_WIDTH)) - 1);
236 
237  // Index the luts wiht the MSB parts already calculated
238  auto phi_lut_q = phi_lut->at(to_two_comp(x_msb, PHI_LUT_ADDR_WIDTH));
239  auto phib_lut_q = phib_lut->at(to_two_comp(tanpsi_msb, PHIB_LUT_ADDR_WIDTH));
240 
241  // Separate this data into the coefficients a and b
242  auto phi_lut_a = phi_lut_q.a;
243  auto phi_lut_b = phi_lut_q.b;
244  auto phib_lut_a = phib_lut_q.a;
245  auto phib_lut_b = phib_lut_q.b;
246 
247  // Do the linear piece-wise operations
248  // At this point all variables that can be negative have already been converted
249  // so will yield negative values when necessary
250  int phi_uncut = (phi_lut_b << PHI_B_SHL_BITS) + x_lsb * phi_lut_a;
251  int psi_uncut = (phib_lut_b << PHIB_B_SHL_BITS) + tanpsi_lsb * phib_lut_a;
252 
253  // Trim phi to its final size
254  int phi = (phi_uncut >> PHI_MULT_SHR_BITS);
255 
256  // Calculate phi_bending from the uncut version of phi and psi, and the trim it to size
257  int phib_uncut = psi_uncut - (phi_uncut >> (PHI_PHIB_RES_DIFF_BITS + PHI_MULT_SHR_BITS - PHIB_MULT_SHR_BITS));
258  int phib = (phib_uncut >> PHIB_MULT_SHR_BITS);
259 
260  double phi_f = (double)phi / pow(2, PHI_SIZE);
261  double phib_f = (double)phib / pow(2, PHIB_SIZE);
262 
263  return std::vector({phi_f, phib_f});
264 }
constexpr int PHI_LUT_ADDR_WIDTH
Definition: constants.h:293
constexpr int PHIB_SIZE
Definition: constants.h:291
constexpr int PHI_PHIB_RES_DIFF_BITS
Definition: constants.h:305
constexpr int PHIB_B_SHL_BITS
Definition: constants.h:300
int to_two_comp(int val, int size)
constexpr int PHI_B_SHL_BITS
Definition: constants.h:294
constexpr int PHI_SIZE
Definition: constants.h:290
std::map< uint32_t, lut_group > luts
constexpr int PHIB_LUT_ADDR_WIDTH
Definition: constants.h:299
constexpr int PHI_MULT_SHR_BITS
Definition: constants.h:295
constexpr int TANPSI_SIZE
Definition: constants.h:289
int from_two_comp(int val, int size)
constexpr int PHIB_MULT_SHR_BITS
Definition: constants.h:301
constexpr int X_SIZE
Definition: constants.h:288
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
int GlobalCoordsObtainer::to_two_comp ( int  val,
int  size 
)
inlineprivate

Definition at line 65 of file GlobalCoordsObtainer.h.

References funct::pow(), and hgcalPerformanceValidation::val.

65  {
66  if (val >= 0)
67  return val;
68  return std::pow(2, size) + val;
69  }
tuple size
Write out results.
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

Member Data Documentation

bool GlobalCoordsObtainer::cmssw_for_global_
private

Definition at line 74 of file GlobalCoordsObtainer.h.

std::vector<global_constant> GlobalCoordsObtainer::global_constants
private

Definition at line 76 of file GlobalCoordsObtainer.h.

edm::FileInPath GlobalCoordsObtainer::global_coords_filename_
private

Definition at line 75 of file GlobalCoordsObtainer.h.

std::map<uint32_t, lut_group> GlobalCoordsObtainer::luts
private

Definition at line 77 of file GlobalCoordsObtainer.h.