CMS 3D CMS Logo
PtAssignmentEngine2017.cc
Go to the documentation of this file.
1 #include "L1Trigger/L1TMuonEndCap/interface/PtAssignmentEngine2017.h"
2 #include "L1Trigger/L1TMuonEndCap/interface/PtAssignmentEngineAux2017.h"
3 
4 #include <iostream>
5 #include <sstream>
6 
7 const PtAssignmentEngineAux2017& PtAssignmentEngine2017::aux() const {
8  static const PtAssignmentEngineAux2017
9  instance; // KK: arguable design solution, but const qualifier makes it thread-safe anyway
10  return instance;
11 }
12 
13 float PtAssignmentEngine2017::scale_pt(const float pt, const int mode) const {
14  emtf_assert(ptLUTVersion_ >= 6);
15 
16  float pt_xml = -99;
17  float pt_scale = -99;
18 
19  // Scaling to achieve 90% efficency at any given L1 pT threshold
20  // For now, a linear scaling based on SingleMu-quality (modes 11, 13, 14, 15), CSC+RPC tracks
21  // Should maybe scale each mode differently in the future - AWB 31.05.17
22 
23  // TRG = (1.2 + 0.015*TRG) * XML
24  // TRG = 1.2*XML / (1 - 0.015*XML)
25  // TRG / XML = 1.2 / (1 - 0.015*XML)
26 
27  // First "physics" LUTs for 2017, deployed June 7
28  if (ptLUTVersion_ >= 6) {
29  pt_xml = fmin(20., pt); // Maximum scale set by muons with XML pT = 20 GeV (scaled pT ~35 GeV)
30  pt_scale = 1.2 / (1 - 0.015 * pt_xml);
31  }
32 
33  return pt_scale;
34 }
35 
36 float PtAssignmentEngine2017::unscale_pt(const float pt, const int mode) const {
37  emtf_assert(ptLUTVersion_ >= 6);
38 
39  float pt_unscale = -99;
40 
41  // First "physics" LUTs for 2017, deployed June 7
42  if (ptLUTVersion_ >= 6) {
43  pt_unscale = 1 / (1.2 + 0.015 * pt);
44  pt_unscale = fmax(pt_unscale, (1 - 0.015 * 20) / 1.2);
45  }
46 
47  return pt_unscale;
48 }
49 
50 PtAssignmentEngine::address_t PtAssignmentEngine2017::calculate_address(const EMTFTrack& track) const {
51  address_t address = 0;
52 
53  EMTFPtLUT data = track.PtLUT();
54 
55  int mode = track.Mode();
56  int theta = track.Theta_fp();
57  int endcap = track.Endcap();
58  int nHits = (mode / 8) + ((mode % 8) / 4) + ((mode % 4) / 2) + ((mode % 2) / 1);
59  emtf_assert(nHits > 1 && nHits < 5);
60 
61  // 'A' is first station in the track, 'B' the second, etc.
62  int mode_ID = -1;
63  int iA = -1, iB = -1, iC = -1, iD = -1;
64  int iAB, iAC, iAD, iBC, iCD;
65 
66  switch (mode) { // Indices for quantities by station or station pair
67  case 15:
68  mode_ID = 0b1;
69  iA = 0;
70  iB = 1;
71  iC = 2;
72  iD = 3;
73  break;
74  case 14:
75  mode_ID = 0b11;
76  iA = 0;
77  iB = 1;
78  iC = 2;
79  break;
80  case 13:
81  mode_ID = 0b10;
82  iA = 0;
83  iB = 1;
84  iC = 3;
85  break;
86  case 11:
87  mode_ID = 0b01;
88  iA = 0;
89  iB = 2;
90  iC = 3;
91  break;
92  case 7:
93  mode_ID = 0b1;
94  iA = 1;
95  iB = 2;
96  iC = 3;
97  break;
98  case 12:
99  mode_ID = 0b111;
100  iA = 0;
101  iB = 1;
102  break;
103  case 10:
104  mode_ID = 0b110;
105  iA = 0;
106  iB = 2;
107  break;
108  case 9:
109  mode_ID = 0b101;
110  iA = 0;
111  iB = 3;
112  break;
113  case 6:
114  mode_ID = 0b100;
115  iA = 1;
116  iB = 2;
117  break;
118  case 5:
119  mode_ID = 0b011;
120  iA = 1;
121  iB = 3;
122  break;
123  case 3:
124  mode_ID = 0b010;
125  iA = 2;
126  iB = 3;
127  break;
128  default:
129  break;
130  }
131  iAB = (iA >= 0 && iB >= 0) ? iA + iB - (iA == 0) : -1;
132  iAC = (iA >= 0 && iC >= 0) ? iA + iC - (iA == 0) : -1;
133  iAD = (iA >= 0 && iD >= 0) ? 2 : -1;
134  iBC = (iB >= 0 && iC >= 0) ? iB + iC : -1;
135  iCD = (iC >= 0 && iD >= 0) ? 5 : -1;
136 
137  // Fill variable info from pT LUT data
138  int st1_ring2 = data.st1_ring2;
139  if (nHits == 4) {
140  int dPhiAB = data.delta_ph[iAB];
141  int dPhiBC = data.delta_ph[iBC];
142  int dPhiCD = data.delta_ph[iCD];
143  int sPhiAB = data.sign_ph[iAB];
144  int sPhiBC = (data.sign_ph[iBC] == sPhiAB);
145  int sPhiCD = (data.sign_ph[iCD] == sPhiAB);
146  int dTheta = data.delta_th[iAD] * (data.sign_th[iAD] ? 1 : -1);
147  int frA = data.fr[iA];
148  int clctA = data.cpattern[iA];
149  int clctB = data.cpattern[iB];
150  int clctC = data.cpattern[iC];
151  int clctD = data.cpattern[iD];
152 
153  // Convert variables to words for pT LUT address
154  dPhiAB = aux().getNLBdPhiBin(dPhiAB, 7, 512);
155  dPhiBC = aux().getNLBdPhiBin(dPhiBC, 5, 256);
156  dPhiCD = aux().getNLBdPhiBin(dPhiCD, 4, 256);
157  dTheta = aux().getdTheta(dTheta, 2);
158  // Combines track theta, stations with RPC hits, and station 1 bend information
159  int mode15_8b = aux().get8bMode15(theta, st1_ring2, endcap, (sPhiAB == 1 ? 1 : -1), clctA, clctB, clctC, clctD);
160 
161  // Form the pT LUT address
162  address |= (dPhiAB & ((1 << 7) - 1)) << (0);
163  address |= (dPhiBC & ((1 << 5) - 1)) << (0 + 7);
164  address |= (dPhiCD & ((1 << 4) - 1)) << (0 + 7 + 5);
165  address |= (sPhiBC & ((1 << 1) - 1)) << (0 + 7 + 5 + 4);
166  address |= (sPhiCD & ((1 << 1) - 1)) << (0 + 7 + 5 + 4 + 1);
167  address |= (dTheta & ((1 << 2) - 1)) << (0 + 7 + 5 + 4 + 1 + 1);
168  address |= (frA & ((1 << 1) - 1)) << (0 + 7 + 5 + 4 + 1 + 1 + 2);
169  address |= (mode15_8b & ((1 << 8) - 1)) << (0 + 7 + 5 + 4 + 1 + 1 + 2 + 1);
170  address |= (mode_ID & ((1 << 1) - 1)) << (0 + 7 + 5 + 4 + 1 + 1 + 2 + 1 + 8);
171  emtf_assert(address < pow(2, 30) && address >= pow(2, 29));
172  } else if (nHits == 3) {
173  int dPhiAB = data.delta_ph[iAB];
174  int dPhiBC = data.delta_ph[iBC];
175  int sPhiAB = data.sign_ph[iAB];
176  int sPhiBC = (data.sign_ph[iBC] == sPhiAB);
177  int dTheta = data.delta_th[iAC] * (data.sign_th[iAC] ? 1 : -1);
178  int frA = data.fr[iA];
179  int frB = data.fr[iB];
180  int clctA = data.cpattern[iA];
181  int clctB = data.cpattern[iB];
182  int clctC = data.cpattern[iC];
183 
184  // Convert variables to words for pT LUT address
185  dPhiAB = aux().getNLBdPhiBin(dPhiAB, 7, 512);
186  dPhiBC = aux().getNLBdPhiBin(dPhiBC, 5, 256);
187  dTheta = aux().getdTheta(dTheta, 3);
188  // Identifies which stations have RPC hits in 3-station tracks
189  int rpc_2b = aux().get2bRPC(clctA, clctB, clctC); // Have to use un-compressed CLCT words
190  clctA = aux().getCLCT(clctA, endcap, (sPhiAB == 1 ? 1 : -1), 2);
191  theta = aux().getTheta(theta, st1_ring2, 5);
192 
193  // Form the pT LUT address
194  address |= (dPhiAB & ((1 << 7) - 1)) << (0);
195  address |= (dPhiBC & ((1 << 5) - 1)) << (0 + 7);
196  address |= (sPhiBC & ((1 << 1) - 1)) << (0 + 7 + 5);
197  address |= (dTheta & ((1 << 3) - 1)) << (0 + 7 + 5 + 1);
198  address |= (frA & ((1 << 1) - 1)) << (0 + 7 + 5 + 1 + 3);
199  int bit = 0;
200  if (mode != 7) {
201  address |= (frB & ((1 << 1) - 1)) << (0 + 7 + 5 + 1 + 3 + 1);
202  bit = 1;
203  }
204  address |= (clctA & ((1 << 2) - 1)) << (0 + 7 + 5 + 1 + 3 + 1 + bit);
205  address |= (rpc_2b & ((1 << 2) - 1)) << (0 + 7 + 5 + 1 + 3 + 1 + bit + 2);
206  address |= (theta & ((1 << 5) - 1)) << (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2);
207  if (mode != 7) {
208  address |= (mode_ID & ((1 << 2) - 1)) << (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2 + 5);
209  emtf_assert(address < pow(2, 29) && address >= pow(2, 27));
210  } else {
211  address |= (mode_ID & ((1 << 1) - 1)) << (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2 + 5);
212  emtf_assert(address < pow(2, 27) && address >= pow(2, 26));
213  }
214 
215  } else if (nHits == 2) {
216  int dPhiAB = data.delta_ph[iAB];
217  int sPhiAB = data.sign_ph[iAB];
218  int dTheta = data.delta_th[iAB] * (data.sign_th[iAB] ? 1 : -1);
219  int frA = data.fr[iA];
220  int frB = data.fr[iB];
221  int clctA = data.cpattern[iA];
222  int clctB = data.cpattern[iB];
223 
224  // Convert variables to words for pT LUT address
225  dPhiAB = aux().getNLBdPhiBin(dPhiAB, 7, 512);
226  dTheta = aux().getdTheta(dTheta, 3);
227  clctA = aux().getCLCT(clctA, endcap, (sPhiAB == 1 ? 1 : -1), 3);
228  clctB = aux().getCLCT(clctB, endcap, (sPhiAB == 1 ? 1 : -1), 3);
229  theta = aux().getTheta(theta, st1_ring2, 5);
230 
231  // Form the pT LUT address
232  address |= (dPhiAB & ((1 << 7) - 1)) << (0);
233  address |= (dTheta & ((1 << 3) - 1)) << (0 + 7);
234  address |= (frA & ((1 << 1) - 1)) << (0 + 7 + 3);
235  address |= (frB & ((1 << 1) - 1)) << (0 + 7 + 3 + 1);
236  address |= (clctA & ((1 << 3) - 1)) << (0 + 7 + 3 + 1 + 1);
237  address |= (clctB & ((1 << 3) - 1)) << (0 + 7 + 3 + 1 + 1 + 3);
238  address |= (theta & ((1 << 5) - 1)) << (0 + 7 + 3 + 1 + 1 + 3 + 3);
239  address |= (mode_ID & ((1 << 3) - 1)) << (0 + 7 + 3 + 1 + 1 + 3 + 3 + 5);
240  emtf_assert(address < pow(2, 26) && address >= pow(2, 24));
241  }
242  return address;
243 } // End function: PtAssignmentEngine2017::calculate_address()
244 
245 // Calculate XML pT from address
246 float PtAssignmentEngine2017::calculate_pt_xml(const address_t& address) const {
247  // std::cout << "Inside calculate_pt_xml, examining address: ";
248  // for (int j = 0; j < 30; j++) {
249  // std::cout << ((address >> (29 - j)) & 0x1);
250  // if ((j % 5) == 4) std::cout << " ";
251  // }
252  // std::cout << std::endl;
253 
254  float pt_xml = 0.;
255  int nHits = -1, mode = -1;
256 
257  emtf_assert(address < pow(2, 30));
258  if (address >= pow(2, 29)) {
259  nHits = 4;
260  mode = 15;
261  } else if (address >= pow(2, 27)) {
262  nHits = 3;
263  } else if (address >= pow(2, 26)) {
264  nHits = 3;
265  mode = 7;
266  } else if (address >= pow(2, 24)) {
267  nHits = 2;
268  } else
269  return pt_xml;
270 
271  // Variables to unpack from the pT LUT address
272  int mode_ID, theta, dTheta;
273  int dPhiAB, dPhiBC = -1, dPhiCD = -1;
274  int sPhiAB, sPhiBC = -1, sPhiCD = -1;
275  int frA, frB = -1;
276  int clctA, clctB = -1;
277  int rpcA, rpcB, rpcC, rpcD;
278  int endcap = 1;
279  sPhiAB = 1; // Assume positive endcap and dPhiAB for unpacking CLCT bend
280  int mode15_8b = -1; // Combines track theta, stations with RPC hits, and station 1 bend information
281  int rpc_2b = -1; // Identifies which stations have RPC hits in 3-station tracks
282 
283  // Unpack variable words from the pT LUT address
284  if (nHits == 4) {
285  dPhiAB = (address >> (0) & ((1 << 7) - 1));
286  dPhiBC = (address >> (0 + 7) & ((1 << 5) - 1));
287  dPhiCD = (address >> (0 + 7 + 5) & ((1 << 4) - 1));
288  sPhiBC = (address >> (0 + 7 + 5 + 4) & ((1 << 1) - 1));
289  sPhiCD = (address >> (0 + 7 + 5 + 4 + 1) & ((1 << 1) - 1));
290  dTheta = (address >> (0 + 7 + 5 + 4 + 1 + 1) & ((1 << 2) - 1));
291  frA = (address >> (0 + 7 + 5 + 4 + 1 + 1 + 2) & ((1 << 1) - 1));
292  mode15_8b = (address >> (0 + 7 + 5 + 4 + 1 + 1 + 2 + 1) & ((1 << 8) - 1));
293  mode_ID = (address >> (0 + 7 + 5 + 4 + 1 + 1 + 2 + 1 + 8) & ((1 << 1) - 1));
294  emtf_assert(address < pow(2, 30));
295  } else if (nHits == 3) {
296  dPhiAB = (address >> (0) & ((1 << 7) - 1));
297  dPhiBC = (address >> (0 + 7) & ((1 << 5) - 1));
298  sPhiBC = (address >> (0 + 7 + 5) & ((1 << 1) - 1));
299  dTheta = (address >> (0 + 7 + 5 + 1) & ((1 << 3) - 1));
300  frA = (address >> (0 + 7 + 5 + 1 + 3) & ((1 << 1) - 1));
301  int bit = 0;
302  if (mode != 7) {
303  frB = (address >> (0 + 7 + 5 + 1 + 3 + 1) & ((1 << 1) - 1));
304  bit = 1;
305  }
306  clctA = (address >> (0 + 7 + 5 + 1 + 3 + 1 + bit) & ((1 << 2) - 1));
307  rpc_2b = (address >> (0 + 7 + 5 + 1 + 3 + 1 + bit + 2) & ((1 << 2) - 1));
308  theta = (address >> (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2) & ((1 << 5) - 1));
309  if (mode != 7) {
310  mode_ID = (address >> (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2 + 5) & ((1 << 2) - 1));
311  emtf_assert(address < pow(2, 29));
312  } else {
313  mode_ID = (address >> (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2 + 5) & ((1 << 1) - 1));
314  emtf_assert(address < pow(2, 27));
315  }
316  } else if (nHits == 2) {
317  dPhiAB = (address >> (0) & ((1 << 7) - 1));
318  dTheta = (address >> (0 + 7) & ((1 << 3) - 1));
319  frA = (address >> (0 + 7 + 3) & ((1 << 1) - 1));
320  frB = (address >> (0 + 7 + 3 + 1) & ((1 << 1) - 1));
321  clctA = (address >> (0 + 7 + 3 + 1 + 1) & ((1 << 3) - 1));
322  clctB = (address >> (0 + 7 + 3 + 1 + 1 + 3) & ((1 << 3) - 1));
323  theta = (address >> (0 + 7 + 3 + 1 + 1 + 3 + 3) & ((1 << 5) - 1));
324  mode_ID = (address >> (0 + 7 + 3 + 1 + 1 + 3 + 3 + 5) & ((1 << 3) - 1));
325  emtf_assert(address < pow(2, 26));
326  }
327 
328  // Infer track mode (and stations with hits) from mode_ID
329  if (nHits == 3 && mode != 7) {
330  switch (mode_ID) {
331  case 0b11:
332  mode = 14;
333  break;
334  case 0b10:
335  mode = 13;
336  break;
337  case 0b01:
338  mode = 11;
339  break;
340  default:
341  break;
342  }
343  } else if (nHits == 2) {
344  switch (mode_ID) {
345  case 0b111:
346  mode = 12;
347  break;
348  case 0b110:
349  mode = 10;
350  break;
351  case 0b101:
352  mode = 9;
353  break;
354  case 0b100:
355  mode = 6;
356  break;
357  case 0b011:
358  mode = 5;
359  break;
360  case 0b010:
361  mode = 3;
362  break;
363  default:
364  break;
365  }
366  }
367 
368  emtf_assert(mode > 0);
369 
370  // Un-compress words from address
371  // For most variables (e.g. theta, dTheta, CLCT) don't need to unpack, since compressed version was used in training
372  int St1_ring2 = -1;
373  if (nHits == 4) {
374  dPhiAB = aux().getdPhiFromBin(dPhiAB, 7, 512);
375  dPhiBC = aux().getdPhiFromBin(dPhiBC, 5, 256) * (sPhiBC == 1 ? 1 : -1);
376  dPhiCD = aux().getdPhiFromBin(dPhiCD, 4, 256) * (sPhiCD == 1 ? 1 : -1);
377  aux().unpack8bMode15(mode15_8b, theta, St1_ring2, endcap, (sPhiAB == 1 ? 1 : -1), clctA, rpcA, rpcB, rpcC, rpcD);
378 
379  // // Check bit-wise compression / de-compression
380  // emtf_assert( dTheta == aux().getdTheta( aux().unpackdTheta( dTheta, 2), 2) );
381  } else if (nHits == 3) {
382  dPhiAB = aux().getdPhiFromBin(dPhiAB, 7, 512);
383  dPhiBC = aux().getdPhiFromBin(dPhiBC, 5, 256) * (sPhiBC == 1 ? 1 : -1);
384  St1_ring2 = aux().unpackSt1Ring2(theta, 5);
385  aux().unpack2bRPC(rpc_2b, rpcA, rpcB, rpcC);
386 
387  // // Check bit-wise compression / de-compression
388  // emtf_assert( dTheta == aux().getdTheta( aux().unpackdTheta( dTheta, 3), 3) );
389  // emtf_assert( clctA == aux().getCLCT( aux().unpackCLCT( clctA, endcap, (sPhiAB == 1 ? 1 : -1), 2),
390  // endcap, (sPhiAB == 1 ? 1 : -1), 2) );
391  // int theta_unp = theta;
392  // aux().unpackTheta( theta_unp, St1_ring2, 5 );
393  // emtf_assert( theta == aux().getTheta(theta_unp, St1_ring2, 5) );
394  } else if (nHits == 2) {
395  dPhiAB = aux().getdPhiFromBin(dPhiAB, 7, 512);
396  St1_ring2 = aux().unpackSt1Ring2(theta, 5);
397 
398  // // Check bit-wise compression / de-compression
399  // emtf_assert( dTheta == aux().getdTheta( aux().unpackdTheta( dTheta, 3), 3) );
400  // emtf_assert( clctA == aux().getCLCT( aux().unpackCLCT( clctA, endcap, (sPhiAB == 1 ? 1 : -1), 3),
401  // endcap, (sPhiAB == 1 ? 1 : -1), 3) );
402  // emtf_assert( clctB == aux().getCLCT( aux().unpackCLCT( clctB, endcap, (sPhiAB == 1 ? 1 : -1), 3),
403  // endcap, (sPhiAB == 1 ? 1 : -1), 3) );
404  // int theta_unp = theta;
405  // aux().unpackTheta( theta_unp, St1_ring2, 5 );
406  // emtf_assert( theta == aux().getTheta(theta_unp, St1_ring2, 5) );
407  }
408 
409  // Fill vectors of variables for XMLs
410  // KK: sequence of variables here should exaclty match <Variables> block produced by TMVA
411  std::vector<int> predictors;
412 
413  // Variables for input to XMLs
414  int dPhiSum4, dPhiSum4A, dPhiSum3, dPhiSum3A, outStPhi;
415 
416  // Convert words into variables for XMLs
417  if (nHits == 4) {
418  predictors = {
419  theta, St1_ring2, dPhiAB, dPhiBC, dPhiCD, dPhiAB + dPhiBC, dPhiAB + dPhiBC + dPhiCD, dPhiBC + dPhiCD, frA, clctA};
420 
421  aux().calcDeltaPhiSums(dPhiSum4,
422  dPhiSum4A,
423  dPhiSum3,
424  dPhiSum3A,
425  outStPhi,
426  dPhiAB,
427  dPhiAB + dPhiBC,
428  dPhiAB + dPhiBC + dPhiCD,
429  dPhiBC,
430  dPhiBC + dPhiCD,
431  dPhiCD);
432 
433  int tmp[10] = {dPhiSum4, dPhiSum4A, dPhiSum3, dPhiSum3A, outStPhi, dTheta, rpcA, rpcB, rpcC, rpcD};
434  predictors.insert(predictors.end(), tmp, tmp + 10);
435  } else if (nHits == 3) {
436  if (mode == 14)
437  predictors = {theta, St1_ring2, dPhiAB, dPhiBC, dPhiAB + dPhiBC, frA, frB, clctA, dTheta, rpcA, rpcB, rpcC};
438  else if (mode == 13)
439  predictors = {theta, St1_ring2, dPhiAB, dPhiAB + dPhiBC, dPhiBC, frA, frB, clctA, dTheta, rpcA, rpcB, rpcC};
440  else if (mode == 11)
441  predictors = {theta, St1_ring2, dPhiBC, dPhiAB, dPhiAB + dPhiBC, frA, frB, clctA, dTheta, rpcA, rpcB, rpcC};
442  else if (mode == 7)
443  predictors = {theta, dPhiAB, dPhiBC, dPhiAB + dPhiBC, frA, clctA, dTheta, rpcA, rpcB, rpcC};
444  } else if (nHits == 2 && mode >= 8) {
445  predictors = {theta, St1_ring2, dPhiAB, frA, frB, clctA, clctB, dTheta, (clctA == 0), (clctB == 0)};
446  } else if (nHits == 2 && mode < 8) {
447  predictors = {theta, dPhiAB, frA, frB, clctA, clctB, dTheta, (clctA == 0), (clctB == 0)};
448  } else {
449  emtf_assert(false && "Incorrect nHits or mode");
450  }
451 
452  // Retreive pT from XMLs
453  std::vector<double> tree_data(predictors.cbegin(), predictors.cend());
454 
455  auto tree_event = std::make_unique<emtf::Event>();
456  tree_event->predictedValue = 0;
457  tree_event->data = tree_data;
458 
459  // forests_.at(mode).predictEvent(tree_event.get(), 400);
460  emtf::Forest& forest = const_cast<emtf::Forest&>(forests_.at(mode));
461  forest.predictEvent(tree_event.get(), 400);
462 
463  // // Adjust this for different XMLs
464  // float log2_pt = tree_event->predictedValue;
465  // pt_xml = pow(2, fmax(0.0, log2_pt)); // Protect against negative values
466 
467  float inv_pt = tree_event->predictedValue;
468  pt_xml = 1.0 / fmax(0.001, inv_pt); // Protect against negative values
469 
470  return pt_xml;
471 
472 } // End function: float PtAssignmentEngine2017::calculate_pt_xml(const address_t& address)
473 
474 // Calculate XML pT directly from track quantities, without forming an address
475 float PtAssignmentEngine2017::calculate_pt_xml(const EMTFTrack& track) const {
476  float pt_xml = 0.;
477 
478  EMTFPtLUT data = track.PtLUT();
479 
480  auto contain = [](const std::vector<int>& vec, int elem) {
481  return (std::find(vec.begin(), vec.end(), elem) != vec.end());
482  };
483 
484  int endcap = track.Endcap();
485  int mode = track.Mode();
486  int theta = track.Theta_fp();
487  int phi = track.Phi_fp();
488  if (!contain(allowedModes_, mode))
489  return pt_xml;
490 
491  // Which stations have hits
492  int st1 = (mode >= 8);
493  int st2 = ((mode % 8) >= 4);
494  int st3 = ((mode % 4) >= 2);
495  int st4 = ((mode % 2) == 1);
496 
497  // Variables for input to XMLs
498  int dPhi_12, dPhi_13, dPhi_14, dPhi_23, dPhi_24, dPhi_34, dPhiSign;
499  int dPhiSum4, dPhiSum4A, dPhiSum3, dPhiSum3A, outStPhi;
500  int dTh_12, dTh_13, dTh_14, dTh_23, dTh_24, dTh_34;
501  int FR_1, FR_2, FR_3, FR_4;
502  int bend_1, bend_2, bend_3, bend_4;
503  int RPC_1, RPC_2, RPC_3, RPC_4;
504  int St1_ring2 = data.st1_ring2;
505 
506  int ph1 = -99, ph2 = -99, ph3 = -99, ph4 = -99;
507  int th1 = -99, th2 = -99, th3 = -99, th4 = -99;
508  int pat1 = -99, pat2 = -99, pat3 = -99, pat4 = -99;
509 
510  // Compute the original phi and theta coordinates
511  if (st2) {
512  ph2 = phi; // Track phi is from station 2 (if it exists), otherwise 3 or 4
513  th2 = theta; // Likewise for track theta
514  if (st1)
515  ph1 = ph2 - data.delta_ph[0] * (data.sign_ph[0] ? 1 : -1);
516  if (st1)
517  th1 = th2 - data.delta_th[0] * (data.sign_th[0] ? 1 : -1);
518  if (st3)
519  ph3 = ph2 + data.delta_ph[3] * (data.sign_ph[3] ? 1 : -1);
520  // Important that phi be from adjacent station pairs (see note below)
521  if (st3 && st4)
522  ph4 = ph3 + data.delta_ph[5] * (data.sign_ph[5] ? 1 : -1);
523  else if (st4)
524  ph4 = ph2 + data.delta_ph[4] * (data.sign_ph[4] ? 1 : -1);
525  // Important that theta be from first-last station pair, not adjacent pairs: delta_th values are "best" for each pair, but
526  // thanks to duplicated CSC LCTs, are not necessarily consistent (or physical) between pairs or between delta_th and delta_ph.
527  // This is an artifact of the firmware implementation of deltas: see src/AngleCalculation.cc.
528  if (st1 && st3)
529  th3 = th1 + data.delta_th[1] * (data.sign_th[1] ? 1 : -1);
530  else if (st3)
531  th3 = th2 + data.delta_th[3] * (data.sign_th[3] ? 1 : -1);
532  if (st1 && st4)
533  th4 = th1 + data.delta_th[2] * (data.sign_th[2] ? 1 : -1);
534  else if (st4)
535  th4 = th2 + data.delta_th[4] * (data.sign_th[4] ? 1 : -1);
536  } else if (st3) {
537  ph3 = phi;
538  th3 = theta;
539  if (st1)
540  ph1 = ph3 - data.delta_ph[1] * (data.sign_ph[1] ? 1 : -1);
541  if (st1)
542  th1 = th3 - data.delta_th[1] * (data.sign_th[1] ? 1 : -1);
543  if (st4)
544  ph4 = ph3 + data.delta_ph[5] * (data.sign_ph[5] ? 1 : -1);
545  if (st1 && st4)
546  th4 = th1 + data.delta_th[2] * (data.sign_th[2] ? 1 : -1);
547  else if (st4)
548  th4 = th3 + data.delta_th[5] * (data.sign_th[5] ? 1 : -1);
549  } else if (st4) {
550  ph4 = phi;
551  th4 = theta;
552  if (st1)
553  ph1 = ph4 - data.delta_ph[2] * (data.sign_ph[2] ? 1 : -1);
554  if (st1)
555  th1 = th4 - data.delta_th[2] * (data.sign_th[2] ? 1 : -1);
556  }
557 
558  if (st1)
559  pat1 = data.cpattern[0];
560  if (st2)
561  pat2 = data.cpattern[1];
562  if (st3)
563  pat3 = data.cpattern[2];
564  if (st4)
565  pat4 = data.cpattern[3];
566 
567  // BEGIN: Identical (almost) to BDT training code in EMTFPtAssign2017/PtRegression_Apr_2017.C
568 
569  theta = aux().calcTrackTheta(th1, th2, th3, th4, St1_ring2, mode, true);
570 
571  aux().calcDeltaPhis(dPhi_12,
572  dPhi_13,
573  dPhi_14,
574  dPhi_23,
575  dPhi_24,
576  dPhi_34,
577  dPhiSign,
578  dPhiSum4,
579  dPhiSum4A,
580  dPhiSum3,
581  dPhiSum3A,
582  outStPhi,
583  ph1,
584  ph2,
585  ph3,
586  ph4,
587  mode,
588  true);
589 
590  aux().calcDeltaThetas(dTh_12, dTh_13, dTh_14, dTh_23, dTh_24, dTh_34, th1, th2, th3, th4, mode, true);
591 
592  FR_1 = (st1 ? data.fr[0] : -99);
593  FR_2 = (st2 ? data.fr[1] : -99);
594  FR_3 = (st3 ? data.fr[2] : -99);
595  FR_4 = (st4 ? data.fr[3] : -99);
596 
597  aux().calcBends(bend_1, bend_2, bend_3, bend_4, pat1, pat2, pat3, pat4, dPhiSign, endcap, mode, true);
598 
599  RPC_1 = (st1 ? (pat1 == 0) : -99);
600  RPC_2 = (st2 ? (pat2 == 0) : -99);
601  RPC_3 = (st3 ? (pat3 == 0) : -99);
602  RPC_4 = (st4 ? (pat4 == 0) : -99);
603 
604  aux().calcRPCs(RPC_1, RPC_2, RPC_3, RPC_4, mode, St1_ring2, theta, true);
605 
606  // END: Identical (almost) to BDT training code in EMTFPtAssign2017/PtRegression_Apr_2017.C
607 
608  // Fill vectors of variables for XMLs
609  // KK: sequence of variables here should exaclty match <Variables> block produced by TMVA
610  std::vector<int> predictors;
611  switch (mode) {
612  case 15: // 1-2-3-4
613  predictors = {theta, St1_ring2, dPhi_12, dPhi_23, dPhi_34, dPhi_13, dPhi_14, dPhi_24, FR_1, bend_1,
614  dPhiSum4, dPhiSum4A, dPhiSum3, dPhiSum3A, outStPhi, dTh_14, RPC_1, RPC_2, RPC_3, RPC_4};
615  break;
616  case 14: // 1-2-3
617  predictors = {theta, St1_ring2, dPhi_12, dPhi_23, dPhi_13, FR_1, FR_2, bend_1, dTh_13, RPC_1, RPC_2, RPC_3};
618  break;
619  case 13: // 1-2-4
620  predictors = {theta, St1_ring2, dPhi_12, dPhi_14, dPhi_24, FR_1, FR_2, bend_1, dTh_14, RPC_1, RPC_2, RPC_4};
621  break;
622  case 11: // 1-3-4
623  predictors = {theta, St1_ring2, dPhi_34, dPhi_13, dPhi_14, FR_1, FR_3, bend_1, dTh_14, RPC_1, RPC_3, RPC_4};
624  break;
625  case 7: // 2-3-4
626  predictors = {theta, dPhi_23, dPhi_34, dPhi_24, FR_2, bend_2, dTh_24, RPC_2, RPC_3, RPC_4};
627  break;
628  case 12: // 1-2
629  predictors = {theta, St1_ring2, dPhi_12, FR_1, FR_2, bend_1, bend_2, dTh_12, RPC_1, RPC_2};
630  break;
631  case 10: // 1-3
632  predictors = {theta, St1_ring2, dPhi_13, FR_1, FR_3, bend_1, bend_3, dTh_13, RPC_1, RPC_3};
633  break;
634  case 9: // 1-4
635  predictors = {theta, St1_ring2, dPhi_14, FR_1, FR_4, bend_1, bend_4, dTh_14, RPC_1, RPC_4};
636  break;
637  case 6: // 2-3
638  predictors = {theta, dPhi_23, FR_2, FR_3, bend_2, bend_3, dTh_23, RPC_2, RPC_3};
639  break;
640  case 5: // 2-4
641  predictors = {theta, dPhi_24, FR_2, FR_4, bend_2, bend_4, dTh_24, RPC_2, RPC_4};
642  break;
643  case 3: // 3-4
644  predictors = {theta, dPhi_34, FR_3, FR_4, bend_3, bend_4, dTh_34, RPC_3, RPC_4};
645  break;
646  }
647 
648  std::vector<double> tree_data(predictors.cbegin(), predictors.cend());
649 
650  auto tree_event = std::make_unique<emtf::Event>();
651  tree_event->predictedValue = 0;
652  tree_event->data = tree_data;
653 
654  // forests_.at(mode).predictEvent(tree_event.get(), 400);
655  emtf::Forest& forest = const_cast<emtf::Forest&>(forests_.at(mode));
656  forest.predictEvent(tree_event.get(), 400);
657 
658  // // Adjust this for different XMLs
659  // float log2_pt = tree_event->predictedValue;
660  // pt_xml = pow(2, fmax(0.0, log2_pt)); // Protect against negative values
661 
662  float inv_pt = tree_event->predictedValue;
663  pt_xml = 1.0 / fmax(0.001, inv_pt); // Protect against negative values
664 
665  return pt_xml;
666 
667 } // End function: float PtAssignmentEngine2017::calculate_pt_xml(const EMTFTrack& track)
PtAssignmentEngineAux2017::calcDeltaThetas
void calcDeltaThetas(int &dTh12, int &dTh13, int &dTh14, int &dTh23, int &dTh24, int &dTh34, const int th1, const int th2, const int th3, const int th4, const int mode, const bool BIT_COMP=false) const
Definition: PtAssignmentEngineAux2017.cc:863
PtAssignmentEngineAux2017::calcTrackTheta
int calcTrackTheta(const int th1, const int th2, const int th3, const int th4, const int ring1, const int mode, const bool BIT_COMP=false) const
Definition: PtAssignmentEngineAux2017.cc:725
PtAssignmentEngineAux2017::getNLBdPhiBin
int getNLBdPhiBin(int dPhi, int bits=7, int max=512) const
Definition: PtAssignmentEngineAux2017.cc:72
HLT_FULL_cff.track
track
Definition: HLT_FULL_cff.py:11713
PtAssignmentEngineAux2017.h
PtAssignmentEngineAux2017::calcBends
void calcBends(int &bend1, int &bend2, int &bend3, int &bend4, const int pat1, const int pat2, const int pat3, const int pat4, const int dPhSign, const int endcap, const int mode, const bool BIT_COMP=false) const
Definition: PtAssignmentEngineAux2017.cc:895
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
l1t::EMTFPtLUT
Definition: EMTFTrack.h:24
PtAssignmentEngine2017::unscale_pt
float unscale_pt(const float pt, const int mode=15) const override
Definition: PtAssignmentEngine2017.cc:36
PtAssignmentEngineAux2017::unpackSt1Ring2
int unpackSt1Ring2(int theta, int bits) const
Definition: PtAssignmentEngineAux2017.cc:473
PtAssignmentEngine::ptLUTVersion_
int ptLUTVersion_
Definition: PtAssignmentEngine.h:56
ALCARECOPromptCalibProdSiPixelAli0T_cff.mode
mode
Definition: ALCARECOPromptCalibProdSiPixelAli0T_cff.py:96
PtAssignmentEngine2017::aux
const PtAssignmentEngineAux2017 & aux() const
Definition: PtAssignmentEngine2017.cc:7
emtf::Forest::predictEvent
void predictEvent(Event *e, unsigned int trees)
Definition: Forest.cc:440
makeMuonMisalignmentScenario.endcap
endcap
Definition: makeMuonMisalignmentScenario.py:320
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
PtAssignmentEngineAux2017::calcDeltaPhis
void calcDeltaPhis(int &dPh12, int &dPh13, int &dPh14, int &dPh23, int &dPh24, int &dPh34, int &dPhSign, int &dPhSum4, int &dPhSum4A, int &dPhSum3, int &dPhSum3A, int &outStPh, const int ph1, const int ph2, const int ph3, const int ph4, const int mode, const bool BIT_COMP=false) const
Definition: PtAssignmentEngineAux2017.cc:751
PtAssignmentEngineAux2017
Definition: PtAssignmentEngineAux2017.h:6
PtAssignmentEngineAux2017::getdPhiFromBin
int getdPhiFromBin(int dPhiBin, int bits=7, int max=512) const
Definition: PtAssignmentEngineAux2017.cc:115
PtAssignmentEngineAux2017::calcRPCs
void calcRPCs(int &RPC1, int &RPC2, int &RPC3, int &RPC4, const int mode, const int st1_ring2, const int theta, const bool BIT_COMP=false) const
Definition: PtAssignmentEngineAux2017.cc:929
PtAssignmentEngineAux2017::getCLCT
int getCLCT(int clct, int endcap, int dPhiSign, int bits=3) const
Definition: PtAssignmentEngineAux2017.cc:139
b1
static constexpr float b1
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
PtAssignmentEngine2017::calculate_pt_xml
float calculate_pt_xml(const address_t &address) const override
Definition: PtAssignmentEngine2017.cc:246
PtAssignmentEngine::address_t
uint64_t address_t
Definition: PtAssignmentEngine.h:21
PtAssignmentEngineAux2017::unpack8bMode15
void unpack8bMode15(int mode15_8b, int &theta, int &st1_ring2, int endcap, int sPhiAB, int &clctA, int &rpcA, int &rpcB, int &rpcC, int &rpcD) const
Definition: PtAssignmentEngineAux2017.cc:586
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
l1t::EMTFTrack
Definition: EMTFTrack.h:42
PtAssignmentEngine2017::calculate_address
address_t calculate_address(const EMTFTrack &track) const override
Definition: PtAssignmentEngine2017.cc:50
PtAssignmentEngine::forests_
std::array< emtf::Forest, 16 > forests_
Definition: PtAssignmentEngine.h:51
PtAssignmentEngine2017.h
PtAssignmentEngineAux2017::get2bRPC
int get2bRPC(int clctA, int clctB, int clctC) const
Definition: PtAssignmentEngineAux2017.cc:492
instance
static PFTauRenderPlugin instance
Definition: PFTauRenderPlugin.cc:70
emtf::Forest
Definition: Forest.h:12
PtAssignmentEngineAux2017::unpack2bRPC
void unpack2bRPC(int rpc_2b, int &rpcA, int &rpcB, int &rpcC) const
Definition: PtAssignmentEngineAux2017.cc:508
PtAssignmentEngine2017::scale_pt
float scale_pt(const float pt, const int mode=15) const override
Definition: PtAssignmentEngine2017.cc:13
DDAxes::phi
PtAssignmentEngine::allowedModes_
std::vector< int > allowedModes_
Definition: PtAssignmentEngine.h:50
PtAssignmentEngineAux2017::calcDeltaPhiSums
void calcDeltaPhiSums(int &dPhSum4, int &dPhSum4A, int &dPhSum3, int &dPhSum3A, int &outStPh, const int dPh12, const int dPh13, const int dPh14, const int dPh23, const int dPh24, const int dPh34) const
Definition: PtAssignmentEngineAux2017.cc:1019
emtf_assert
#define emtf_assert(expr)
Definition: DebugTools.h:18
data
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:79
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
PtAssignmentEngineAux2017::getTheta
int getTheta(int theta, int ring2, int bits=5) const
Definition: PtAssignmentEngineAux2017.cc:411
PtAssignmentEngineAux2017::get8bMode15
int get8bMode15(int theta, int st1_ring2, int endcap, int sPhiAB, int clctA, int clctB, int clctC, int clctD) const
Definition: PtAssignmentEngineAux2017.cc:524
PtAssignmentEngineAux2017::getdTheta
int getdTheta(int dTheta, int bits=3) const
Definition: PtAssignmentEngineAux2017.cc:313
Generated for CMSSW Reference Manual by doxygen 1.8.16