CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
L1TMuonBarrelKalmanAlgo.cc
Go to the documentation of this file.
1 #include <cmath>
3 #include "ap_int.h"
4 #include "ap_fixed.h"
5 
7  : verbose_(settings.getParameter<bool>("verbose")),
8  lutService_(new L1TMuonBarrelKalmanLUTs(settings.getParameter<std::string>("lutFile"))),
9  initK_(settings.getParameter<std::vector<double> >("initialK")),
10  initK2_(settings.getParameter<std::vector<double> >("initialK2")),
11  eLoss_(settings.getParameter<std::vector<double> >("eLoss")),
12  aPhi_(settings.getParameter<std::vector<double> >("aPhi")),
13  aPhiB_(settings.getParameter<std::vector<double> >("aPhiB")),
14  aPhiBNLO_(settings.getParameter<std::vector<double> >("aPhiBNLO")),
15  bPhi_(settings.getParameter<std::vector<double> >("bPhi")),
16  bPhiB_(settings.getParameter<std::vector<double> >("bPhiB")),
17  phiAt2_(settings.getParameter<double>("phiAt2")),
18 
19  chiSquare_(settings.getParameter<std::vector<double> >("chiSquare")),
20  chiSquareCutPattern_(settings.getParameter<std::vector<int> >("chiSquareCutPattern")),
21  chiSquareCutCurv_(settings.getParameter<std::vector<int> >("chiSquareCutCurvMax")),
22  chiSquareCut_(settings.getParameter<std::vector<int> >("chiSquareCut")),
23  trackComp_(settings.getParameter<std::vector<double> >("trackComp")),
24  trackCompErr1_(settings.getParameter<std::vector<double> >("trackCompErr1")),
25  trackCompErr2_(settings.getParameter<std::vector<double> >("trackCompErr2")),
26  trackCompPattern_(settings.getParameter<std::vector<int> >("trackCompCutPattern")),
27  trackCompCutCurv_(settings.getParameter<std::vector<int> >("trackCompCutCurvMax")),
28  trackCompCut_(settings.getParameter<std::vector<int> >("trackCompCut")),
29  chiSquareCutTight_(settings.getParameter<std::vector<int> >("chiSquareCutTight")),
30  combos4_(settings.getParameter<std::vector<int> >("combos4")),
31  combos3_(settings.getParameter<std::vector<int> >("combos3")),
32  combos2_(settings.getParameter<std::vector<int> >("combos2")),
33  combos1_(settings.getParameter<std::vector<int> >("combos1")),
34 
35  useOfflineAlgo_(settings.getParameter<bool>("useOfflineAlgo")),
36  mScatteringPhi_(settings.getParameter<std::vector<double> >("mScatteringPhi")),
37  mScatteringPhiB_(settings.getParameter<std::vector<double> >("mScatteringPhiB")),
38  pointResolutionPhi_(settings.getParameter<double>("pointResolutionPhi")),
39  pointResolutionPhiB_(settings.getParameter<double>("pointResolutionPhiB")),
40  pointResolutionPhiBH_(settings.getParameter<std::vector<double> >("pointResolutionPhiBH")),
41  pointResolutionPhiBL_(settings.getParameter<std::vector<double> >("pointResolutionPhiBL")),
42  pointResolutionVertex_(settings.getParameter<double>("pointResolutionVertex"))
43 
44 {}
45 
46 std::pair<bool, uint> L1TMuonBarrelKalmanAlgo::getByCode(const L1MuKBMTrackCollection& tracks, int mask) {
47  for (uint i = 0; i < tracks.size(); ++i) {
48  printf("Code=%d, track=%d\n", tracks[i].hitPattern(), mask);
49  if (tracks[i].hitPattern() == mask)
50  return std::make_pair(true, i);
51  }
52  return std::make_pair(false, 0);
53 }
54 
56  // int K = fabs(track.curvatureAtVertex());
57 
58  //calibration
59  int sign, signValid;
60 
61  if (track.curvatureAtVertex() == 0) {
62  sign = 0;
63  signValid = 0;
64  } else if (track.curvatureAtVertex() > 0) {
65  sign = 0;
66  signValid = 1;
67  } else {
68  sign = 1;
69  signValid = 1;
70  }
71 
72  // if (K<22)
73  // K=22;
74 
75  // if (K>4095)
76  // K=4095;
77 
78  int pt = ptLUT(track.curvatureAtVertex());
79 
80  // int K2 = fabs(track.curvatureAtMuon());
81  // if (K2<22)
82  // K2=22;
83 
84  // if (K2>4095)
85  // K2=4095;
86  int pt2 = ptLUT(track.curvatureAtMuon()) / 2;
87  int eta = track.hasFineEta() ? track.fineEta() : track.coarseEta();
88 
89  // int phi2 = track.phiAtMuon()>>2;
90  // float phi_f = float(phi2);
91  //double kPhi = 57.2958/0.625/1024.;
92  //int phi = 24+int(floor(kPhi*phi_f));
93  // if (phi > 69) phi = 69;
94  // if (phi < -8) phi = -8;
95  int phi2 = track.phiAtMuon() >> 2;
96  int tmp = fp_product(0.0895386, phi2, 14);
97  int phi = 24 + tmp;
98 
99  int processor = track.sector();
100  int HF = track.hasFineEta();
101 
102  int quality = 12 | (rank(track) >> 6);
103 
104  int dxy = abs(track.dxy()) >> 8;
105  if (dxy > 3)
106  dxy = 3;
107 
108  int trackAddr;
109  std::map<int, int> addr = trackAddress(track, trackAddr);
110 
111  l1t::RegionalMuonCand muon(pt, phi, eta, sign, signValid, quality, processor, l1t::bmtf, addr);
112  muon.setHwHF(HF);
113  muon.setHwPtUnconstrained(pt2);
114  muon.setHwDXY(dxy);
115 
116  //nw the words!
117  uint32_t word1 = pt;
118  word1 = word1 | quality << 9;
119  word1 = word1 | (twosCompToBits(eta)) << 13;
120  word1 = word1 | HF << 22;
121  word1 = word1 | (twosCompToBits(phi)) << 23;
122 
123  uint32_t word2 = sign;
124  word2 = word2 | signValid << 1;
125  word2 = word2 | dxy << 2;
126  word2 = word2 | trackAddr << 4;
127  word2 = word2 | (twosCompToBits(track.wheel())) << 20;
128  word2 = word2 | pt2 << 23;
129  muon.setDataword(word2, word1);
130  return muon;
131 }
132 
134  const L1MuKBMTrack& track,
135  std::unique_ptr<l1t::RegionalMuonCandBxCollection>& out) {
136  out->push_back(bx, convertToBMTF(track));
137 }
138 
139 // std::pair<bool,uint> L1TMuonBarrelKalmanAlgo::match(const L1MuKBMTCombinedStubRef& seed, const L1MuKBMTCombinedStubRefVector& stubs,int step) {
140 // L1MuKBMTCombinedStubRefVector selected;
141 
142 // bool found=false;
143 // uint best=0;
144 // int distance=100000;
145 // uint N=0;
146 // for (const auto& stub :stubs) {
147 // N=N+1;
148 // if (stub->stNum()!=step)
149 // continue;
150 
151 // int d = fabs(wrapAround(((correctedPhi(seed,seed->scNum())-correctedPhi(stub,seed->scNum()))>>3),1024));
152 // if (d<distance) {
153 // distance = d;
154 // best=N-1;
155 // found=true;
156 // }
157 // }
158 // return std::make_pair(found,best);
159 // }
160 
162  if (info[i] < info[j])
163  return i;
164  else
165  return j;
166 }
167 
169  const L1MuKBMTCombinedStubRefVector& stubs,
170  int step) {
172 
173  std::map<uint, uint> diffInfo;
174  for (uint i = 0; i < 12; ++i) {
175  diffInfo[i] = 60000;
176  }
177 
178  std::map<uint, uint> stubInfo;
179 
180  int sector = seed->scNum();
181  int previousSector = sector - 1;
182  int nextSector = sector + 1;
183  if (sector == 0) {
184  previousSector = 11;
185  }
186  if (sector == 11) {
187  nextSector = 0;
188  }
189 
190  int wheel = seed->whNum();
191  int innerWheel = 0;
192  if (wheel == -2)
193  innerWheel = -1;
194  if (wheel == -1)
195  innerWheel = 0;
196  if (wheel == 0)
197  innerWheel = 1982;
198  if (wheel == 1)
199  innerWheel = 0;
200  if (wheel == 2)
201  innerWheel = 1;
202 
203  //First align the data
204  uint N = 0;
205  for (const auto& stub : stubs) {
206  N = N + 1;
207 
208  if (stub->stNum() != step)
209  continue;
210 
211  uint distance =
212  fabs(wrapAround(((correctedPhi(seed, seed->scNum()) - correctedPhi(stub, seed->scNum())) >> 3), 1024));
213 
214  if (stub->scNum() == previousSector) {
215  if (stub->whNum() == wheel) {
216  if (!stub->tag()) {
217  diffInfo[0] = distance;
218  stubInfo[0] = N - 1;
219  } else {
220  diffInfo[1] = distance;
221  stubInfo[1] = N - 1;
222  }
223  } else if (stub->whNum() == innerWheel) {
224  if (!stub->tag()) {
225  diffInfo[2] = distance;
226  stubInfo[2] = N - 1;
227  } else {
228  diffInfo[3] = distance;
229  stubInfo[3] = N - 1;
230  }
231  }
232  } else if (stub->scNum() == sector) {
233  if (stub->whNum() == wheel) {
234  if (!stub->tag()) {
235  diffInfo[4] = distance;
236  stubInfo[4] = N - 1;
237  } else {
238  diffInfo[5] = distance;
239  stubInfo[5] = N - 1;
240  }
241  } else if (stub->whNum() == innerWheel) {
242  if (!stub->tag()) {
243  diffInfo[6] = distance;
244  stubInfo[6] = N - 1;
245  } else {
246  diffInfo[7] = distance;
247  stubInfo[7] = N - 1;
248  }
249  }
250  } else if (stub->scNum() == nextSector) {
251  if (stub->whNum() == wheel) {
252  if (!stub->tag()) {
253  diffInfo[8] = distance;
254  stubInfo[8] = N - 1;
255  } else {
256  diffInfo[9] = distance;
257  stubInfo[9] = N - 1;
258  }
259  } else if (stub->whNum() == innerWheel) {
260  if (!stub->tag()) {
261  diffInfo[10] = distance;
262  stubInfo[10] = N - 1;
263  } else {
264  diffInfo[11] = distance;
265  stubInfo[11] = N - 1;
266  }
267  }
268  }
269  }
270 
271  uint s1_1 = matchAbs(diffInfo, 0, 1);
272  uint s1_2 = matchAbs(diffInfo, 2, 3);
273  uint s1_3 = matchAbs(diffInfo, 4, 5);
274  uint s1_4 = matchAbs(diffInfo, 6, 7);
275  uint s1_5 = matchAbs(diffInfo, 8, 9);
276  uint s1_6 = matchAbs(diffInfo, 10, 11);
277 
278  uint s2_1 = matchAbs(diffInfo, s1_1, s1_2);
279  uint s2_2 = matchAbs(diffInfo, s1_3, s1_4);
280  uint s2_3 = matchAbs(diffInfo, s1_5, s1_6);
281 
282  uint s3_1 = matchAbs(diffInfo, s2_1, s2_2);
283 
284  uint s4 = matchAbs(diffInfo, s3_1, s2_3);
285 
286  if (diffInfo[s4] != 60000)
287  return std::make_pair(true, stubInfo[s4]);
288  else
289  return std::make_pair(false, 0);
290 }
291 
293  //Promote phiB to 12 bits
294  return 8 * stub->phiB();
295 }
296 
298  if (stub->scNum() == sector) {
299  return stub->phi();
300  } else if ((stub->scNum() == sector - 1) || (stub->scNum() == 11 && sector == 0)) {
301  return stub->phi() - 2144;
302  } else if ((stub->scNum() == sector + 1) || (stub->scNum() == 0 && sector == 11)) {
303  return stub->phi() + 2144;
304  }
305  return stub->phi();
306 }
307 
309  unsigned int mask = 0;
310  for (const auto& stub : track.stubs()) {
311  mask = mask + round(pow(2, stub->stNum() - 1));
312  }
313  return mask;
314 }
315 
316 int L1TMuonBarrelKalmanAlgo::customBitmask(unsigned int bit1, unsigned int bit2, unsigned int bit3, unsigned int bit4) {
317  return bit1 * 1 + bit2 * 2 + bit3 * 4 + bit4 * 8;
318 }
319 
320 bool L1TMuonBarrelKalmanAlgo::getBit(int bitmask, int pos) { return (bitmask & (1 << pos)) >> pos; }
321 
323  int K = track.curvature();
324  int phi = track.positionAngle();
325  int phiB = track.bendingAngle();
326  unsigned int step = track.step();
327 
328  //energy loss term only for MU->VERTEX
329  //int offset=int(charge*eLoss_[step-1]*K*K);
330  // if (fabs(offset)>4096)
331  // offset=4096*offset/fabs(offset);
332  int charge = 1;
333  if (K != 0)
334  charge = K / fabs(K);
335 
336  int KBound = K;
337 
338  if (KBound > 4095)
339  KBound = 4095;
340  if (KBound < -4095)
341  KBound = -4095;
342 
343  int deltaK = 0;
344  int KNew = 0;
345  if (step == 1) {
346  int addr = KBound / 2;
347  if (addr < 0)
348  addr = (-KBound) / 2;
349  deltaK = 2 * addr - int(2 * addr / (1 + eLoss_[step - 1] * addr));
350 
351  if (verbose_)
352  printf("propagate to vertex K=%d deltaK=%d addr=%d\n", K, deltaK, addr);
353  }
354 
355  if (K >= 0)
356  KNew = K - deltaK;
357  else
358  KNew = K + deltaK;
359 
360  //phi propagation
361  ap_fixed<BITSCURV, BITSCURV> phi11 = ap_fixed<BITSPARAM + 1, 2>(aPhi_[step - 1]) * ap_fixed<BITSCURV, BITSCURV>(K);
362  ap_fixed<BITSPHIB, BITSPHIB> phi12 =
363  ap_fixed<BITSPARAM + 1, 2>(-bPhi_[step - 1]) * ap_fixed<BITSPHIB, BITSPHIB>(phiB);
364 
365  if (verbose_) {
366  printf("phi prop = %d * %f = %d, %d * %f = %d\n",
367  K,
368  ap_fixed<BITSPARAM + 1, 2>(aPhi_[step - 1]).to_float(),
369  phi11.to_int(),
370  phiB,
371  ap_fixed<BITSPARAM + 1, 2>(-bPhi_[step - 1]).to_float(),
372  phi12.to_int());
373  }
374  int phiNew = ap_fixed<BITSPHI, BITSPHI>(phi + phi11 + phi12);
375 
376  //phiB propagation
377  ap_fixed<BITSCURV, BITSCURV> phiB11 = ap_fixed<BITSPARAM, 1>(aPhiB_[step - 1]) * ap_fixed<BITSCURV, BITSCURV>(K);
378  ap_fixed<BITSPHIB + 1, BITSPHIB + 1> phiB12 =
379  ap_ufixed<BITSPARAM + 1, 1>(bPhiB_[step - 1]) * ap_fixed<BITSPHIB, BITSPHIB>(phiB);
380  int phiBNew = ap_fixed<13, 13>(phiB11 + phiB12);
381  if (verbose_) {
382  printf("phiB prop = %d * %f = %d, %d * %f = %d\n",
383  K,
384  ap_fixed<BITSPARAM + 1, 2>(aPhiB_[step - 1]).to_float(),
385  phiB11.to_int(),
386  phiB,
387  ap_ufixed<BITSPARAM + 1, 1>(bPhiB_[step - 1]).to_float(),
388  phiB12.to_int());
389  }
390 
391  //Only for the propagation to vertex we use the LUT for better precision and the full function
392  if (step == 1) {
393  int addr = KBound / 2;
394  // Extra steps to mimic firmware for vertex prop
395  ap_ufixed<11, 11> dxyOffset = (int)fabs(aPhiB_[step - 1] * addr / (1 + charge * aPhiBNLO_[step - 1] * addr));
396  ap_fixed<12, 12> DXY;
397  if (addr > 0)
398  DXY = -dxyOffset;
399  else
400  DXY = dxyOffset;
401  phiBNew = ap_fixed<BITSPHIB, BITSPHIB>(DXY - ap_fixed<BITSPHIB, BITSPHIB>(phiB));
402  if (verbose_) {
403  printf("Vertex phiB prop = %d - %d = %d\n", DXY.to_int(), ap_fixed<BITSPHIB, BITSPHIB>(phiB).to_int(), phiBNew);
404  }
405  }
407  //Rest of the stuff is for the offline version only
408  //where we want to check what is happening in the covariaznce matrix
409 
410  //Create the transformation matrix
411  double a[9];
412  a[0] = 1.;
413  a[1] = 0.0;
414  a[2] = 0.0;
415  a[3] = aPhi_[step - 1];
416  // a[3] = 0.0;
417  a[4] = 1.0;
418  a[5] = -bPhi_[step - 1];
419  //a[6]=0.0;
420  a[6] = aPhiB_[step - 1];
421  if (step == 1)
422  a[6] = aPhiB_[step - 1] / 2.0;
423 
424  a[7] = 0.0;
425  a[8] = bPhiB_[step - 1];
426 
427  ROOT::Math::SMatrix<double, 3> P(a, 9);
428 
429  const std::vector<double>& covLine = track.covariance();
430  L1MuKBMTrack::CovarianceMatrix cov(covLine.begin(), covLine.end());
431  cov = ROOT::Math::Similarity(P, cov);
432 
433  //Add the multiple scattering
434  double phiRMS = mScatteringPhi_[step - 1] * K * K;
435  double phiBRMS = mScatteringPhiB_[step - 1] * K * K;
436 
437  std::vector<double> b(6);
438  b[0] = 0;
439  b[1] = 0;
440  b[2] = phiRMS;
441  b[3] = 0;
442  b[4] = 0;
443  b[5] = phiBRMS;
444 
445  reco::Candidate::CovarianceMatrix MS(b.begin(), b.end());
446 
447  cov = cov + MS;
448 
449  if (verbose_) {
450  printf("Covariance term for phiB = %f\n", cov(2, 2));
451  printf("Multiple scattering term for phiB = %f\n", MS(2, 2));
452  }
453 
454  track.setCovariance(cov);
455  track.setCoordinates(step - 1, KNew, phiNew, phiBNew);
456 }
457 
458 bool L1TMuonBarrelKalmanAlgo::update(L1MuKBMTrack& track, const L1MuKBMTCombinedStubRef& stub, int mask, int seedQual) {
459  updateEta(track, stub);
460  if (useOfflineAlgo_) {
461  if (mask == 3 || mask == 5 || mask == 9 || mask == 6 || mask == 10 || mask == 12)
462  return updateOffline(track, stub);
463  else
464  return updateOffline1D(track, stub);
465 
466  } else
467  return updateLUT(track, stub, mask, seedQual);
468 }
469 
471  int trackK = track.curvature();
472  int trackPhi = track.positionAngle();
473  int trackPhiB = track.bendingAngle();
474 
475  int phi = correctedPhi(stub, track.sector());
476  int phiB = correctedPhiB(stub);
477 
478  Vector2 residual;
479  residual[0] = phi - trackPhi;
480  residual[1] = phiB - trackPhiB;
481 
482  Matrix23 H;
483  H(0, 0) = 0.0;
484  H(0, 1) = 1.0;
485  H(0, 2) = 0.0;
486  H(1, 0) = 0.0;
487  H(1, 1) = 0.0;
488  H(1, 2) = 1.0;
489 
491  R(0, 0) = pointResolutionPhi_;
492  R(0, 1) = 0.0;
493  R(1, 0) = 0.0;
494  if (stub->quality() < 4)
495  R(1, 1) = pointResolutionPhiBL_[track.step() - 1];
496  else
497  R(1, 1) = pointResolutionPhiBH_[track.step() - 1];
498 
499  const std::vector<double>& covLine = track.covariance();
500  L1MuKBMTrack::CovarianceMatrix cov(covLine.begin(), covLine.end());
501 
502  CovarianceMatrix2 S = ROOT::Math::Similarity(H, cov) + R;
503  if (!S.Invert())
504  return false;
505  Matrix32 Gain = cov * ROOT::Math::Transpose(H) * S;
506 
507  track.setKalmanGain(
508  track.step(), fabs(trackK), Gain(0, 0), Gain(0, 1), Gain(1, 0), Gain(1, 1), Gain(2, 0), Gain(2, 1));
509 
510  int KNew = (trackK + int(Gain(0, 0) * residual(0) + Gain(0, 1) * residual(1)));
511  if (fabs(KNew) > 8192)
512  return false;
513 
514  int phiNew = wrapAround(trackPhi + residual(0), 8192);
515  int phiBNew = wrapAround(trackPhiB + int(Gain(2, 0) * residual(0) + Gain(2, 1) * residual(1)), 4096);
516 
517  track.setResidual(stub->stNum() - 1, fabs(phi - phiNew) + fabs(phiB - phiBNew) / 8);
518 
519  if (verbose_) {
520  printf("residual %d - %d = %d %d - %d = %d\n", phi, trackPhi, int(residual[0]), phiB, trackPhiB, int(residual[1]));
521  printf("Gains offline: %f %f %f %f\n", Gain(0, 0), Gain(0, 1), Gain(2, 0), Gain(2, 1));
522  printf(" K = %d + %f * %f + %f * %f\n", trackK, Gain(0, 0), residual(0), Gain(0, 1), residual(1));
523  printf(" phiB = %d + %f * %f + %f * %f\n", trackPhiB, Gain(2, 0), residual(0), Gain(2, 1), residual(1));
524  }
525 
526  track.setCoordinates(track.step(), KNew, phiNew, phiBNew);
527  Matrix33 covNew = cov - Gain * (H * cov);
529 
530  c(0, 0) = covNew(0, 0);
531  c(0, 1) = covNew(0, 1);
532  c(0, 2) = covNew(0, 2);
533  c(1, 0) = covNew(1, 0);
534  c(1, 1) = covNew(1, 1);
535  c(1, 2) = covNew(1, 2);
536  c(2, 0) = covNew(2, 0);
537  c(2, 1) = covNew(2, 1);
538  c(2, 2) = covNew(2, 2);
539  if (verbose_) {
540  printf("Post Fit Covariance Matrix %f %f %f\n", cov(0, 0), cov(1, 1), cov(2, 2));
541  }
542 
543  track.setCovariance(c);
544  track.addStub(stub);
545  track.setHitPattern(hitPattern(track));
546 
547  return true;
548 }
549 
551  int trackK = track.curvature();
552  int trackPhi = track.positionAngle();
553  int trackPhiB = track.bendingAngle();
554 
555  int phi = correctedPhi(stub, track.sector());
556 
557  double residual = phi - trackPhi;
558 
559  if (verbose_)
560  printf("residuals %d - %d = %d\n", phi, trackPhi, int(residual));
561 
562  Matrix13 H;
563  H(0, 0) = 0.0;
564  H(0, 1) = 1.0;
565  H(0, 2) = 0.0;
566 
567  const std::vector<double>& covLine = track.covariance();
568  L1MuKBMTrack::CovarianceMatrix cov(covLine.begin(), covLine.end());
569 
570  double S = ROOT::Math::Similarity(H, cov)(0, 0) + pointResolutionPhi_;
571 
572  if (S == 0.0)
573  return false;
574  Matrix31 Gain = cov * ROOT::Math::Transpose(H) / S;
575 
576  track.setKalmanGain(track.step(), fabs(trackK), Gain(0, 0), 0.0, Gain(1, 0), 0.0, Gain(2, 0), 0.0);
577  if (verbose_)
578  printf("Gains: %f %f\n", Gain(0, 0), Gain(2, 0));
579 
580  int KNew = wrapAround(trackK + int(Gain(0, 0) * residual), 8192);
581  int phiNew = wrapAround(trackPhi + residual, 8192);
582  int phiBNew = wrapAround(trackPhiB + int(Gain(2, 0) * residual), 4096);
583  track.setCoordinates(track.step(), KNew, phiNew, phiBNew);
584  Matrix33 covNew = cov - Gain * (H * cov);
586 
587  if (verbose_) {
588  printf("phiUpdate: %d %d\n", int(Gain(0, 0) * residual), int(Gain(2, 0) * residual));
589  }
590 
591  c(0, 0) = covNew(0, 0);
592  c(0, 1) = covNew(0, 1);
593  c(0, 2) = covNew(0, 2);
594  c(1, 0) = covNew(1, 0);
595  c(1, 1) = covNew(1, 1);
596  c(1, 2) = covNew(1, 2);
597  c(2, 0) = covNew(2, 0);
598  c(2, 1) = covNew(2, 1);
599  c(2, 2) = covNew(2, 2);
600  track.setCovariance(c);
601  track.addStub(stub);
602  track.setHitPattern(hitPattern(track));
603 
604  return true;
605 }
606 
608  const L1MuKBMTCombinedStubRef& stub,
609  int mask,
610  int seedQual) {
611  int trackK = track.curvature();
612  int trackPhi = track.positionAngle();
613  int trackPhiB = track.bendingAngle();
614 
615  int phi = correctedPhi(stub, track.sector());
616  int phiB = correctedPhiB(stub);
617 
618  Vector2 residual;
619  ap_fixed<BITSPHI + 1, BITSPHI + 1> residualPhi = phi - trackPhi;
620  ap_fixed<BITSPHIB + 1, BITSPHIB + 1> residualPhiB = phiB - trackPhiB;
621 
622  if (verbose_)
623  printf("residual %d - %d = %d %d - %d = %d\n",
624  phi,
625  trackPhi,
626  residualPhi.to_int(),
627  phiB,
628  trackPhiB,
629  residualPhiB.to_int());
630 
631  uint absK = fabs(trackK);
632  if (absK > 4095)
633  absK = 4095;
634 
635  std::vector<float> GAIN;
636  //For the three stub stuff use only gains 0 and 4
637  if (!(mask == 3 || mask == 5 || mask == 9 || mask == 6 || mask == 10 || mask == 12)) {
638  GAIN = lutService_->trackGain(track.step(), track.hitPattern(), absK / 4);
639  GAIN[1] = 0.0;
640  GAIN[3] = 0.0;
641 
642  } else {
643  GAIN = lutService_->trackGain2(track.step(), track.hitPattern(), absK / 8, seedQual, stub->quality());
644  }
645  if (verbose_) {
646  printf("Gains (fp): %f %f %f %f\n", GAIN[0], GAIN[1], GAIN[2], GAIN[3]);
647  if (!(mask == 3 || mask == 5 || mask == 9 || mask == 6 || mask == 10 || mask == 12))
648  printf("Addr=%d gain0=%f gain4=-%f\n",
649  absK / 4,
650  ap_ufixed<GAIN_0, GAIN_0INT>(GAIN[0]).to_float(),
651  ap_ufixed<GAIN_4, GAIN_4INT>(GAIN[2]).to_float());
652  else
653  printf("Addr=%d %f -%f %f %f\n",
654  absK / 4,
655  ap_fixed<GAIN2_0, GAIN2_0INT>(GAIN[0]).to_float(),
656  ap_ufixed<GAIN2_1, GAIN2_1INT>(GAIN[1]).to_float(),
657  ap_ufixed<GAIN2_4, GAIN2_4INT>(GAIN[2]).to_float(),
658  ap_ufixed<GAIN2_5, GAIN2_5INT>(GAIN[3]).to_float());
659  }
660 
661  track.setKalmanGain(track.step(), fabs(trackK), GAIN[0], GAIN[1], 1, 0, GAIN[2], GAIN[3]);
662 
663  int KNew;
664  if (!(mask == 3 || mask == 5 || mask == 9 || mask == 6 || mask == 10 || mask == 12)) {
665  KNew = ap_fixed<BITSPHI + 9, BITSPHI + 9>(ap_fixed<BITSCURV, BITSCURV>(trackK) +
666  ap_ufixed<GAIN_0, GAIN_0INT>(GAIN[0]) * residualPhi);
667  } else {
668  ap_fixed<BITSPHI + 7, BITSPHI + 7> k11 = ap_fixed<GAIN2_0, GAIN2_0INT>(GAIN[0]) * residualPhi;
669  ap_fixed<BITSPHIB + 4, BITSPHIB + 4> k12 = ap_ufixed<GAIN2_1, GAIN2_1INT>(GAIN[1]) * residualPhiB;
670  KNew = ap_fixed<BITSPHI + 9, BITSPHI + 9>(ap_fixed<BITSCURV, BITSCURV>(trackK) + k11 - k12);
671  }
672  if (fabs(KNew) >= 8191)
673  return false;
674  KNew = wrapAround(KNew, 8192);
675  int phiNew = phi;
676 
677  //different products for different firmware logic
678  ap_fixed<BITSPHI + 5, BITSPHI + 5> pbdouble_0 = ap_ufixed<GAIN2_4, GAIN2_4INT>(GAIN[2]) * residualPhi;
679  ap_fixed<BITSPHIB + 24, BITSPHIB + 4> pb_1 = ap_ufixed<GAIN2_5, GAIN2_5INT>(GAIN[3]) * residualPhiB;
680  ap_fixed<BITSPHI + 9, BITSPHI + 5> pb_0 = ap_ufixed<GAIN_4, GAIN_4INT>(GAIN[2]) * residualPhi;
681 
682  if (verbose_) {
683  printf("phiupdate %f %f %f\n", pb_0.to_float(), pb_1.to_float(), pbdouble_0.to_float());
684  }
685 
686  int phiBNew;
687  if (!(mask == 3 || mask == 5 || mask == 9 || mask == 6 || mask == 10 || mask == 12)) {
688  phiBNew = ap_fixed<BITSPHI + 8, BITSPHI + 8>(ap_fixed<BITSPHIB, BITSPHIB>(trackPhiB) -
689  ap_ufixed<GAIN_4, GAIN_4INT>(GAIN[2]) * residualPhi);
690 
691  if (fabs(phiBNew) >= 4095)
692  return false;
693  } else {
694  phiBNew = ap_fixed<BITSPHI + 7, BITSPHI + 7>(ap_fixed<BITSPHIB, BITSPHIB>(trackPhiB) + pb_1 - pbdouble_0);
695  if (fabs(phiBNew) >= 4095)
696  return false;
697  }
698  track.setCoordinates(track.step(), KNew, phiNew, phiBNew);
699  track.addStub(stub);
700  track.setHitPattern(hitPattern(track));
701  return true;
702 }
703 
705 
707  if (useOfflineAlgo_)
709  else
710  vertexConstraintLUT(track);
711 }
712 
714  double residual = -track.dxy();
715  Matrix13 H;
716  H(0, 0) = 0;
717  H(0, 1) = 0;
718  H(0, 2) = 1;
719 
720  const std::vector<double>& covLine = track.covariance();
721  L1MuKBMTrack::CovarianceMatrix cov(covLine.begin(), covLine.end());
722 
723  double S = (ROOT::Math::Similarity(H, cov))(0, 0) + pointResolutionVertex_;
724  S = 1.0 / S;
725  Matrix31 Gain = cov * (ROOT::Math::Transpose(H)) * S;
726  track.setKalmanGain(track.step(), fabs(track.curvature()), Gain(0, 0), Gain(1, 0), Gain(2, 0));
727 
728  if (verbose_) {
729  printf("sigma3=%f sigma6=%f\n", cov(0, 3), cov(3, 3));
730  printf(" K = %d + %f * %f\n", track.curvature(), Gain(0, 0), residual);
731  }
732 
733  int KNew = wrapAround(int(track.curvature() + Gain(0, 0) * residual), 8192);
734  int phiNew = wrapAround(int(track.positionAngle() + Gain(1, 0) * residual), 8192);
735  int dxyNew = wrapAround(int(track.dxy() + Gain(2, 0) * residual), 8192);
736  if (verbose_)
737  printf("Post fit impact parameter=%d\n", dxyNew);
738  track.setCoordinatesAtVertex(KNew, phiNew, -residual);
739  Matrix33 covNew = cov - Gain * (H * cov);
741  c(0, 0) = covNew(0, 0);
742  c(0, 1) = covNew(0, 1);
743  c(0, 2) = covNew(0, 2);
744  c(1, 0) = covNew(1, 0);
745  c(1, 1) = covNew(1, 1);
746  c(1, 2) = covNew(1, 2);
747  c(2, 0) = covNew(2, 0);
748  c(2, 1) = covNew(2, 1);
749  c(2, 2) = covNew(2, 2);
750  track.setCovariance(c);
751  // track.covariance = track.covariance - Gain*H*track.covariance;
752 }
753 
755  double residual = -track.dxy();
756  uint absK = fabs(track.curvature());
757  if (absK > 2047)
758  absK = 2047;
759 
760  std::pair<float, float> GAIN = lutService_->vertexGain(track.hitPattern(), absK / 2);
761  track.setKalmanGain(track.step(), fabs(track.curvature()), GAIN.first, GAIN.second, -1);
762 
763  ap_fixed<BITSCURV, BITSCURV> k_0 =
764  -(ap_ufixed<GAIN_V0, GAIN_V0INT>(fabs(GAIN.first))) * ap_fixed<BITSPHIB, BITSPHIB>(residual);
765  int KNew = ap_fixed<BITSCURV, BITSCURV>(k_0 + ap_fixed<BITSCURV, BITSCURV>(track.curvature()));
766 
767  if (verbose_) {
768  printf("VERTEX GAIN(%d)= -%f * %d = %d\n",
769  absK / 2,
770  ap_ufixed<GAIN_V0, GAIN_V0INT>(fabs(GAIN.first)).to_float(),
771  ap_fixed<BITSPHIB, BITSPHIB>(residual).to_int(),
772  k_0.to_int());
773  }
774 
775  int p_0 = fp_product(GAIN.second, int(residual), 7);
776  int phiNew = wrapAround(track.positionAngle() + p_0, 8192);
777  track.setCoordinatesAtVertex(KNew, phiNew, -residual);
778 }
779 
781  int K, etaINT;
782 
783  if (track.hasFineEta())
784  etaINT = track.fineEta();
785  else
786  etaINT = track.coarseEta();
787 
788  double lsb = 1.25 / float(1 << 13);
789  double lsbEta = 0.010875;
790 
791  if (vertex) {
792  int charge = 1;
793  if (track.curvatureAtVertex() < 0)
794  charge = -1;
795  double pt = double(ptLUT(track.curvatureAtVertex())) / 2.0;
796 
797  double phi = track.sector() * M_PI / 6.0 + track.phiAtVertex() * M_PI / (6 * 2048.) - 2 * M_PI;
798 
799  double eta = etaINT * lsbEta;
800  track.setPtEtaPhi(pt, eta, phi);
801  track.setCharge(charge);
802  } else {
803  K = track.curvatureAtMuon();
804  if (K == 0)
805  K = 1;
806 
807  if (fabs(K) < 46)
808  K = 46 * K / fabs(K);
809  double pt = 1.0 / (lsb * fabs(K));
810  if (pt < 1.6)
811  pt = 1.6;
812  track.setPtUnconstrained(pt);
813  }
814 }
815 
816 std::pair<bool, L1MuKBMTrack> L1TMuonBarrelKalmanAlgo::chain(const L1MuKBMTCombinedStubRef& seed,
817  const L1MuKBMTCombinedStubRefVector& stubs) {
818  L1MuKBMTrackCollection pretracks;
819  std::vector<int> combinatorics;
820  int seedQual;
821  switch (seed->stNum()) {
822  case 1:
823  combinatorics = combos1_;
824  break;
825  case 2:
826  combinatorics = combos2_;
827  break;
828 
829  case 3:
830  combinatorics = combos3_;
831  break;
832 
833  case 4:
834  combinatorics = combos4_;
835  break;
836 
837  default:
838  printf("Something really bad happend\n");
839  }
840 
841  L1MuKBMTrack nullTrack(seed, correctedPhi(seed, seed->scNum()), correctedPhiB(seed));
842  seedQual = seed->quality();
843  for (const auto& mask : combinatorics) {
844  L1MuKBMTrack track(seed, correctedPhi(seed, seed->scNum()), correctedPhiB(seed));
845  int phiB = correctedPhiB(seed);
846  int charge;
847  if (phiB == 0)
848  charge = 0;
849  else
850  charge = phiB / fabs(phiB);
851 
852  int address = phiB;
853  if (track.step() == 4 && (fabs(seed->phiB()) > 15))
854  address = charge * 15 * 8;
855 
856  if (track.step() == 3 && (fabs(seed->phiB()) > 30))
857  address = charge * 30 * 8;
858  if (track.step() == 2 && (fabs(seed->phiB()) > 127))
859  address = charge * 127 * 8;
860  int initialK = int(initK_[seed->stNum() - 1] * address / (1 + initK2_[seed->stNum() - 1] * charge * address));
861  if (initialK > 8191)
862  initialK = 8191;
863  if (initialK < -8191)
864  initialK = -8191;
865 
866  track.setCoordinates(seed->stNum(), initialK, correctedPhi(seed, seed->scNum()), phiB);
867  if (seed->quality() < 4) {
868  track.setCoordinates(seed->stNum(), 0, correctedPhi(seed, seed->scNum()), 0);
869  }
870 
871  track.setHitPattern(hitPattern(track));
872  //set covariance
874 
875  float DK = 512 * 512.;
876  covariance(0, 0) = DK;
877  covariance(0, 1) = 0;
878  covariance(0, 2) = 0;
879  covariance(1, 0) = 0;
880  covariance(1, 1) = float(pointResolutionPhi_);
881  covariance(1, 2) = 0;
882  covariance(2, 0) = 0;
883  covariance(2, 1) = 0;
884  if (!(mask == 3 || mask == 5 || mask == 9 || mask == 6 || mask == 10 || mask == 12))
885  covariance(2, 2) = float(pointResolutionPhiB_);
886  else {
887  if (seed->quality() < 4)
888  covariance(2, 2) = float(pointResolutionPhiBL_[seed->stNum() - 1]);
889  else
890  covariance(2, 2) = float(pointResolutionPhiBH_[seed->stNum() - 1]);
891  }
892  track.setCovariance(covariance);
893 
894  //
895  if (verbose_) {
896  printf("New Kalman fit staring at step=%d, phi=%d,phiB=%d with curvature=%d\n",
897  track.step(),
898  track.positionAngle(),
899  track.bendingAngle(),
900  track.curvature());
901  printf("BITMASK:");
902  for (unsigned int i = 0; i < 4; ++i)
903  printf("%d", getBit(mask, i));
904  printf("\n");
905  printf("------------------------------------------------------\n");
906  printf("------------------------------------------------------\n");
907  printf("------------------------------------------------------\n");
908  printf("stubs:\n");
909  for (const auto& stub : stubs)
910  printf("station=%d phi=%d phiB=%d qual=%d tag=%d sector=%d wheel=%d fineEta= %d %d\n",
911  stub->stNum(),
912  correctedPhi(stub, seed->scNum()),
913  correctedPhiB(stub),
914  stub->quality(),
915  stub->tag(),
916  stub->scNum(),
917  stub->whNum(),
918  stub->eta1(),
919  stub->eta2());
920  printf("------------------------------------------------------\n");
921  printf("------------------------------------------------------\n");
922  }
923 
924  int phiAtStation2 = 0;
925 
926  while (track.step() > 0) {
927  // muon station 1
928  if (track.step() == 1) {
929  track.setCoordinatesAtMuon(track.curvature(), track.positionAngle(), track.bendingAngle());
930  phiAtStation2 = phiAt2(track);
931  bool passed = estimateChiSquare(track);
932  if (!passed)
933  break;
936  //calculate coarse eta
938 
939  if (verbose_)
940  printf("Unconstrained PT in Muon System: pt=%f\n", track.ptUnconstrained());
941  }
942 
943  propagate(track);
944  if (verbose_)
945  printf("propagated Coordinates step:%d,phi=%d,phiB=%d,K=%d\n",
946  track.step(),
947  track.positionAngle(),
948  track.bendingAngle(),
949  track.curvature());
950 
951  if (track.step() > 0)
952  if (getBit(mask, track.step() - 1)) {
953  std::pair<bool, uint> bestStub = match(seed, stubs, track.step());
954  if ((!bestStub.first) || (!update(track, stubs[bestStub.second], mask, seedQual)))
955  break;
956  if (verbose_) {
957  printf("updated Coordinates step:%d,phi=%d,phiB=%d,K=%d\n",
958  track.step(),
959  track.positionAngle(),
960  track.bendingAngle(),
961  track.curvature());
962  }
963  }
964 
965  if (track.step() == 0) {
966  track.setCoordinatesAtVertex(track.curvature(), track.positionAngle(), track.bendingAngle());
967  if (verbose_)
968  printf(" Coordinates before vertex constraint step:%d,phi=%d,dxy=%d,K=%d\n",
969  track.step(),
970  track.phiAtVertex(),
971  track.dxy(),
972  track.curvatureAtVertex());
973  if (verbose_)
974  printf("Chi Square = %d\n", track.approxChi2());
975 
978  if (verbose_) {
979  printf(" Coordinates after vertex constraint step:%d,phi=%d,dxy=%d,K=%d maximum local chi2=%d\n",
980  track.step(),
981  track.phiAtVertex(),
982  track.dxy(),
983  track.curvatureAtVertex(),
984  track.approxChi2());
985  printf("------------------------------------------------------\n");
986  printf("------------------------------------------------------\n");
987  }
989  //set the coordinates at muon to include phi at station 2
990  track.setCoordinatesAtMuon(track.curvatureAtMuon(), phiAtStation2, track.phiBAtMuon());
991  track.setRank(rank(track));
992  if (verbose_)
993  printf("Floating point coordinates at vertex: pt=%f, eta=%f phi=%f\n", track.pt(), track.eta(), track.phi());
994  pretracks.push_back(track);
995  }
996  }
997  }
998 
999  //Now for all the pretracks we need only one
1000  L1MuKBMTrackCollection cleaned = clean(pretracks, seed->stNum());
1001 
1002  if (!cleaned.empty()) {
1003  return std::make_pair(true, cleaned[0]);
1004  }
1005  return std::make_pair(false, nullTrack);
1006 }
1007 
1009  //here we have a simplification of the algorithm for the sake of the emulator - rsult is identical
1010  // we apply cuts on the firmware as |u -u'|^2 < a+b *K^2
1011  int K = track.curvatureAtMuon();
1012 
1013  uint chi = 0;
1014 
1015  // const double PHI[4]={0.0,0.249,0.543,0.786};
1016  // const double DROR[4]={0.0,0.182,0.430,0.677};
1017 
1018  int coords = wrapAround((track.phiAtMuon() + track.phiBAtMuon()) >> 4, 512);
1019  for (const auto& stub : track.stubs()) {
1020  int AK = wrapAround(fp_product(-chiSquare_[stub->stNum() - 1], K >> 4, 8), 256);
1021  int stubCoords = wrapAround((correctedPhi(stub, track.sector()) >> 4) + (stub->phiB() >> 1), 512);
1022  int diff1 = wrapAround(stubCoords - coords, 1024);
1023  uint delta = wrapAround(abs(diff1 + AK), 2048);
1024  chi = chi + delta;
1025  if (verbose_)
1026  printf("Chi Square stub for track with pattern=%d coords=%d -> AK=%d stubCoords=%d diff=%d delta=%d\n",
1027  track.hitPattern(),
1028  coords,
1029  AK,
1030  stubCoords,
1031  diff1,
1032  delta);
1033  }
1034 
1035  if (chi > 127)
1036  chi = 127;
1037  // for (const auto& stub: track.stubs()) {
1038  // int deltaPhi = (correctedPhi(stub,track.sector())-track.phiAtMuon())>>3;
1039  // int AK = fp_product(PHI[stub->stNum()-1],K>>3,8);
1040  // int BPB = fp_product(DROR[stub->stNum()-1],track.phiBAtMuon()>>3,8);
1041  // chi=chi+abs(deltaPhi-AK-BPB);
1042  // }
1043  // // }
1044 
1045  track.setApproxChi2(chi);
1046  for (uint i = 0; i < chiSquareCutPattern_.size(); ++i) {
1047  if (track.hitPattern() == chiSquareCutPattern_[i] && fabs(K) < chiSquareCutCurv_[i] &&
1048  track.approxChi2() > chiSquareCut_[i])
1049  return false;
1050  }
1051  return true;
1052 }
1053 
1055  int K = track.curvatureAtVertex() >> 4;
1056 
1057  if (track.stubs().size() != 2) {
1058  track.setTrackCompatibility(0);
1059  return;
1060  }
1061 
1062  uint stubSel = 1;
1063  if (track.stubs()[0]->quality() > track.stubs()[1]->quality())
1064  stubSel = 0;
1065  const L1MuKBMTCombinedStubRef& stub = track.stubs()[stubSel];
1066 
1067  if (verbose_) {
1068  printf("stubsel %d phi=%d phiB=%d\n", stubSel, stub->phi(), stub->phiB());
1069  }
1070 
1071  ap_ufixed<BITSCURV - 5, BITSCURV - 5> absK;
1072  if (K < 0)
1073  absK = -K;
1074  else
1075  absK = K;
1076 
1077  ap_fixed<12, 12> diff = ap_int<10>(stub->phiB()) -
1078  ap_ufixed<5, 1>(trackComp_[stub->stNum() - 1]) * ap_fixed<BITSCURV - 4, BITSCURV - 4>(K);
1079  ap_ufixed<11, 11> delta;
1080  if (diff.is_neg())
1081  delta = -diff;
1082  else
1083  delta = diff;
1084 
1085  ap_ufixed<BITSCURV - 5, BITSCURV - 5> err =
1086  ap_uint<3>(trackCompErr1_[stub->stNum() - 1]) + ap_ufixed<5, 0>(trackCompErr2_[stub->stNum() - 1]) * absK;
1087  track.setTrackCompatibility(((int)delta) / ((int)err));
1088  for (uint i = 0; i < trackCompPattern_.size(); ++i) {
1089  int deltaMax = ap_ufixed<BITSCURV, BITSCURV>(err * trackCompCut_[i]);
1090  if (verbose_) {
1091  if (track.hitPattern() == trackCompPattern_[i]) {
1092  printf("delta = %d = abs(%d - %f*%d\n", delta.to_int(), stub->phiB(), trackComp_[stub->stNum() - 1], K);
1093  printf("err = %d = %f + %f*%d\n",
1094  err.to_int(),
1095  trackCompErr1_[stub->stNum() - 1],
1096  trackCompErr2_[stub->stNum() - 1],
1097  absK.to_int());
1098  printf("deltaMax = %d = %d*%d\n", deltaMax, err.to_int(), trackCompCut_[i]);
1099  }
1100  }
1101  if ((track.hitPattern() == trackCompPattern_[i]) && ((int)absK < trackCompCutCurv_[i]) &&
1102  ((track.approxChi2() > chiSquareCutTight_[i]) || (delta > deltaMax))) {
1103  track.setCoordinatesAtVertex(8191, track.phiAtVertex(), track.dxy());
1104  break;
1105  }
1106  }
1107 }
1108 
1110  // int offset=0;
1111  uint chi = track.approxChi2() > 127 ? 127 : track.approxChi2();
1112  if (hitPattern(track) == customBitmask(0, 0, 1, 1)) {
1113  return 60;
1114  }
1115  // return offset+(track.stubs().size()*2+track.quality())*80-track.approxChi2();
1116  return 160 + (track.stubs().size()) * 20 - chi;
1117 }
1118 
1120  if (value > maximum - 1)
1121  return wrapAround(value - 2 * maximum, maximum);
1122  if (value < -maximum)
1123  return wrapAround(value + 2 * maximum, maximum);
1124  return value;
1125 }
1126 
1127 int L1TMuonBarrelKalmanAlgo::encode(bool ownwheel, int sector, bool tag) {
1128  if (ownwheel) {
1129  if (sector == 0) {
1130  if (tag)
1131  return 9;
1132  else
1133  return 8;
1134  } else if (sector == 1) {
1135  if (tag)
1136  return 11;
1137  else
1138  return 10;
1139 
1140  } else {
1141  if (tag)
1142  return 13;
1143  else
1144  return 12;
1145  }
1146 
1147  } else {
1148  if (sector == 0) {
1149  if (tag)
1150  return 1;
1151  else
1152  return 0;
1153  } else if (sector == 1) {
1154  if (tag)
1155  return 3;
1156  else
1157  return 2;
1158 
1159  } else {
1160  if (tag)
1161  return 5;
1162  else
1163  return 4;
1164  }
1165  }
1166  return 15;
1167 }
1168 
1170  std::map<int, int> out;
1171 
1172  out[l1t::RegionalMuonCand::kWheelSide] = track.wheel() < 0;
1173  if (track.wheel() == -2)
1175  else if (track.wheel() == -1)
1177  else if (track.wheel() == 0)
1179  else if (track.wheel() == 1)
1181  else if (track.wheel() == 2)
1183  else
1193  //out[l1t::RegionalMuonCand::kNumBmtfSubAddr]=0; // This is commented out for better data/MC agreement
1194 
1195  for (const auto& stub : track.stubs()) {
1196  bool ownwheel = stub->whNum() == track.wheel();
1197  int sector = 0;
1198  if ((stub->scNum() == track.sector() + 1) || (stub->scNum() == 0 && track.sector() == 11))
1199  sector = +1;
1200  if ((stub->scNum() == track.sector() - 1) || (stub->scNum() == 11 && track.sector() == 0))
1201  sector = -1;
1202  int addr = encode(ownwheel, sector, stub->tag());
1203 
1204  if (stub->stNum() == 4) {
1205  if (stub->tag())
1206  addr = 1;
1207  else
1208  addr = 2;
1210  }
1211  if (stub->stNum() == 3) {
1213  }
1214  if (stub->stNum() == 2) {
1216  }
1217  if (stub->stNum() == 1) {
1219  }
1220  }
1221 
1222  word = 0;
1223  word = word | out[l1t::RegionalMuonCand::kStat4] << 12;
1224  word = word | out[l1t::RegionalMuonCand::kStat3] << 8;
1225  word = word | out[l1t::RegionalMuonCand::kStat2] << 4;
1226  word = word | out[l1t::RegionalMuonCand::kStat1];
1227 
1228  return out;
1229 }
1230 
1232  if (q >= 0)
1233  return q;
1234  else
1235  return (~q) + 1;
1236 }
1237 
1239  // return long(a*(1<<bits)*b)>>bits;
1240  return (long((a * (1 << bits)) * b)) >> bits;
1241 }
1242 
1244  int charge = (K >= 0) ? +1 : -1;
1245  float lsb = 1.25 / float(1 << 13);
1246  float FK = fabs(K);
1247 
1248  if (FK > 2047)
1249  FK = 2047.;
1250  if (FK < 26)
1251  FK = 26.;
1252 
1253  FK = FK * lsb;
1254 
1255  //step 1 -material and B-field
1256  FK = .8569 * FK / (1.0 + 0.1144 * FK);
1257  //step 2 - misalignment
1258  FK = FK - charge * 1.23e-03;
1259  //Get to BMTF scale
1260  FK = FK / 1.17;
1261 
1262  int pt = 0;
1263  if (FK != 0)
1264  pt = int(2.0 / FK);
1265 
1266  if (pt > 511)
1267  pt = 511;
1268 
1269  if (pt < 8)
1270  pt = 8;
1271 
1272  return pt;
1273 }
1274 
1277 
1278  std::map<uint, int> infoRank;
1279  std::map<uint, L1MuKBMTrack> infoTrack;
1280  for (uint i = 3; i <= 15; ++i) {
1281  if (i == 4 || i == 8)
1282  continue;
1283  infoRank[i] = -1;
1284  }
1285 
1286  for (const auto& track : tracks) {
1287  infoRank[track.hitPattern()] = rank(track);
1288  infoTrack[track.hitPattern()] = track;
1289  }
1290 
1291  int selected = 15;
1292  if (seed == 4) //station 4 seeded
1293  {
1294  int sel6 = infoRank[10] >= infoRank[12] ? 10 : 12;
1295  int sel5 = infoRank[14] >= infoRank[9] ? 14 : 9;
1296  int sel4 = infoRank[11] >= infoRank[13] ? 11 : 13;
1297  int sel3 = infoRank[sel6] >= infoRank[sel5] ? sel6 : sel5;
1298  int sel2 = infoRank[sel4] >= infoRank[sel3] ? sel4 : sel3;
1299  selected = infoRank[15] >= infoRank[sel2] ? 15 : sel2;
1300  }
1301  if (seed == 3) //station 3 seeded
1302  {
1303  int sel2 = infoRank[5] >= infoRank[6] ? 5 : 6;
1304  selected = infoRank[7] >= infoRank[sel2] ? 7 : sel2;
1305  }
1306  if (seed == 2) //station 2 seeded
1307  selected = 3;
1308 
1309  auto search = infoTrack.find(selected);
1310  if (search != infoTrack.end())
1311  out.push_back(search->second);
1312 
1313  return out;
1314 }
1315 
1317  if (stub->qeta1() != 0 && stub->qeta2() != 0) {
1318  return 0;
1319  }
1320  if (stub->qeta1() == 0) {
1321  return 0;
1322  }
1323  // return (stub->qeta1()*4+stub->stNum());
1324  return (stub->qeta1());
1325 }
1326 
1328  uint pattern = track.hitPattern();
1329  int wheel = track.stubs()[0]->whNum();
1330  uint awheel = fabs(wheel);
1331  int sign = 1;
1332  if (wheel < 0)
1333  sign = -1;
1334  uint nstubs = track.stubs().size();
1335  uint mask = 0;
1336  for (unsigned int i = 0; i < track.stubs().size(); ++i) {
1337  if (fabs(track.stubs()[i]->whNum()) != awheel)
1338  mask = mask | (1 << i);
1339  }
1340  mask = (awheel << nstubs) | mask;
1341  track.setCoarseEta(sign * lutService_->coarseEta(pattern, mask));
1342  int sumweights = 0;
1343  int sums = 0;
1344 
1345  for (const auto& stub : track.stubs()) {
1346  uint rank = etaStubRank(stub);
1347  if (rank == 0)
1348  continue;
1349  // printf("Stub station=%d rank=%d values=%d %d\n",stub->stNum(),rank,stub->eta1(),stub->eta2());
1350  sumweights += rank;
1351  sums += rank * stub->eta1();
1352  }
1353  //0.5 0.332031 0.25 0.199219 0.164063
1354  float factor;
1355  if (sumweights == 1)
1356  factor = 1.0;
1357  else if (sumweights == 2)
1358  factor = 0.5;
1359  else if (sumweights == 3)
1360  factor = 0.332031;
1361  else if (sumweights == 4)
1362  factor = 0.25;
1363  else if (sumweights == 5)
1364  factor = 0.199219;
1365  else if (sumweights == 6)
1366  factor = 0.164063;
1367  else
1368  factor = 0.0;
1369 
1370  int eta = 0;
1371  if (sums > 0)
1372  eta = fp_product(factor, sums, 10);
1373  else
1374  eta = -fp_product(factor, fabs(sums), 10);
1375 
1376  //int eta=int(factor*sums);
1377  // printf("Eta debug %f *%d=%d\n",factor,sums,eta);
1378 
1379  if (sumweights > 0)
1380  track.setFineEta(eta);
1381 }
1382 
1384  //If there is stub at station 2 use this else propagate from 1
1385  for (const auto& stub : track.stubs())
1386  if (stub->stNum() == 2)
1387  return correctedPhi(stub, track.sector());
1388 
1389  ap_fixed<BITSPHI, BITSPHI> phi = track.phiAtMuon();
1390  ap_fixed<BITSPHIB, BITSPHIB> phiB = track.phiBAtMuon();
1391  ap_fixed<BITSPARAM, 1> phiAt2 = phiAt2_;
1392  int phiNew = ap_fixed<BITSPHI + 1, BITSPHI + 1, AP_TRN_ZERO, AP_SAT>(phi + phiAt2 * phiB);
1393 
1394  if (verbose_)
1395  printf("Phi at second station=%d\n", phiNew);
1396  return phiNew;
1397 }
std::vector< double > bPhiB_
int correctedPhiB(const L1MuKBMTCombinedStubRef &)
int fineEta() const
Definition: L1MuKBMTrack.cc:68
l1t::RegionalMuonCand convertToBMTF(const L1MuKBMTrack &track)
std::vector< double > aPhi_
static const TGPicture * info(bool iBackgroundIsBlack)
const edm::EventSetup & c
void setCoordinatesAtVertex(int, int, int)
Definition: L1MuKBMTrack.cc:87
std::vector< double > aPhiB_
std::vector< double > aPhiBNLO_
L1TMuonBarrelKalmanAlgo(const edm::ParameterSet &settings)
std::pair< bool, uint > match(const L1MuKBMTCombinedStubRef &, const L1MuKBMTCombinedStubRefVector &, int)
std::vector< T >::const_iterator search(const cond::Time_t &val, const std::vector< T > &container)
Definition: IOVProxy.cc:21
int step() const
Definition: L1MuKBMTrack.cc:60
std::vector< double > mScatteringPhiB_
uint32_t const *__restrict__ Quality * quality
int bendingAngle() const
Definition: L1MuKBMTrack.cc:51
int coarseEta() const
Definition: L1MuKBMTrack.cc:53
double sign(double x)
void setTrackCompatibility(int)
int phiAt2(const L1MuKBMTrack &track)
bool update(L1MuKBMTrack &, const L1MuKBMTCombinedStubRef &, int, int)
int phiAtVertex() const
Definition: L1MuKBMTrack.cc:43
auto const & tracks
cannot be loose
int rank(const L1MuKBMTrack &)
void setPtUnconstrained(float)
std::vector< int > trackCompPattern_
std::vector< edm::Ref< L1MuKBMTCombinedStubCollection > > L1MuKBMTCombinedStubRefVector
void setHitPattern(int)
bool updateOffline1D(L1MuKBMTrack &, const L1MuKBMTCombinedStubRef &)
void estimateCompatibility(L1MuKBMTrack &)
void setPtEtaPhi(double, double, double)
bool updateLUT(L1MuKBMTrack &, const L1MuKBMTCombinedStubRef &, int, int)
void vertexConstraintLUT(L1MuKBMTrack &)
void setCharge(Charge q) final
set electric charge
int phiBAtMuon() const
Definition: L1MuKBMTrack.cc:39
std::pair< bool, uint > getByCode(const L1MuKBMTrackCollection &tracks, int mask)
int fp_product(float, int, uint)
std::vector< double > chiSquare_
uint64_t word
int curvatureAtMuon() const
Definition: L1MuKBMTrack.cc:37
std::vector< double > pointResolutionPhiBH_
bool hasFineEta() const
Definition: L1MuKBMTrack.cc:70
std::vector< double > initK2_
ROOT::Math::SMatrix< double, 1, 3 > Matrix13
std::pair< float, float > vertexGain(uint, uint)
void setCoordinates(int, int, int, int)
Definition: L1MuKBMTrack.cc:80
printf("params %d %f %f %f\n", minT, eps, errmax, chi2max)
int approxChi2() const
Definition: L1MuKBMTrack.cc:55
std::vector< double > trackComp_
void setCovariance(const CovarianceMatrix &)
std::vector< double > mScatteringPhi_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
L1TMuonBarrelKalmanLUTs * lutService_
void setCoarseEta(int)
Definition: L1MuKBMTrack.cc:99
void setDataword(uint32_t msbs, uint32_t lsbs)
void setHwHF(bool bit)
Set HF (halo / fine eta) bit (EMTF: halo -&gt; 1; BMTF: fine eta -&gt; 1)
void vertexConstraint(L1MuKBMTrack &)
int customBitmask(unsigned int, unsigned int, unsigned int, unsigned int)
std::vector< double > eLoss_
ROOT::Math::SMatrix< double, 3, 1 > Matrix31
bool estimateChiSquare(L1MuKBMTrack &)
void vertexConstraintOffline(L1MuKBMTrack &)
void setFineEta(int)
const L1MuKBMTCombinedStubRefVector & stubs() const
Definition: L1MuKBMTrack.cc:76
ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2 > > CovarianceMatrix2
int curvature() const
Definition: L1MuKBMTrack.cc:47
uint matchAbs(std::map< uint, uint > &, uint, uint)
std::vector< int > chiSquareCutTight_
int phiAtMuon() const
Definition: L1MuKBMTrack.cc:38
std::vector< double > bPhi_
std::vector< double > initK_
void updateEta(L1MuKBMTrack &, const L1MuKBMTCombinedStubRef &)
#define M_PI
int wheel() const
Definition: L1MuKBMTrack.cc:62
void setFloatingPointValues(L1MuKBMTrack &, bool)
void setHwDXY(int bits)
Set compressed impact parameter with respect to beamspot (4 bits)
L1MuKBMTrackCollection clean(const L1MuKBMTrackCollection &, uint)
#define N
Definition: blowfish.cc:9
int encode(bool ownwheel, int sector, bool tag)
std::vector< int > chiSquareCutPattern_
const std::vector< double > & covariance() const
int dxy() const
Definition: L1MuKBMTrack.cc:45
std::pair< bool, L1MuKBMTrack > chain(const L1MuKBMTCombinedStubRef &, const L1MuKBMTCombinedStubRefVector &)
void calculateEta(L1MuKBMTrack &track)
int curvatureAtVertex() const
Definition: L1MuKBMTrack.cc:41
double b
Definition: hdecay.h:118
std::pair< OmniClusterRef, TrackingParticleRef > P
std::vector< int > chiSquareCut_
double S(const TLorentzVector &, const TLorentzVector &)
Definition: Particle.cc:97
void addBMTFMuon(int, const L1MuKBMTrack &, std::unique_ptr< l1t::RegionalMuonCandBxCollection > &)
int hitPattern(const L1MuKBMTrack &)
ROOT::Math::SVector< double, 2 > Vector2
void setKalmanGain(unsigned int step, unsigned int K, float a1, float a2, float a3, float a4=0, float a5=0, float a6=0)
std::vector< L1MuKBMTrack > L1MuKBMTrackCollection
Definition: L1MuKBMTrack.h:15
double a
Definition: hdecay.h:119
ROOT::Math::SMatrix< double, 2, 3 > Matrix23
int correctedPhi(const L1MuKBMTCombinedStubRef &, int)
std::vector< double > trackCompErr1_
std::vector< float > trackGain2(uint, uint, uint, uint, uint)
int sector() const
Definition: L1MuKBMTrack.cc:61
ROOT::Math::SMatrix< double, 3, 2 > Matrix32
std::vector< int > chiSquareCutCurv_
ROOT::Math::SMatrix< double, 3, 3 > Matrix33
void addStub(const L1MuKBMTCombinedStubRef &)
step
Definition: StallMonitor.cc:94
int hitPattern() const
Definition: L1MuKBMTrack.cc:58
std::vector< double > trackCompErr2_
void setResidual(uint, int)
std::vector< float > trackGain(uint, uint, uint)
math::Error< dimension >::type CovarianceMatrix
covariance error matrix (3x3)
Definition: Candidate.h:46
tmp
align.sh
Definition: createJobs.py:716
void setHwPtUnconstrained(int bits)
Set compressed second displaced pT as transmitted by hardware LSB = 1.0 (8 bits)
int positionAngle() const
Definition: L1MuKBMTrack.cc:49
bool updateOffline(L1MuKBMTrack &, const L1MuKBMTCombinedStubRef &)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
uint etaStubRank(const L1MuKBMTCombinedStubRef &)
std::vector< int > trackCompCutCurv_
void setApproxChi2(int)
std::map< int, int > trackAddress(const L1MuKBMTrack &, int &)
std::vector< double > pointResolutionPhiBL_