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