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