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