11 TrackDerTable::TrackDerTable(
Settings const& settings) : settings_(settings) {
26 for (
int i = 0;
i < (1 <<
Nlay_);
i++) {
30 for (
int i = 0;
i < (1 << (2 *
Ndisk_));
i++) {
40 unsigned int diskmask,
41 unsigned int alphaindex,
42 unsigned int rinvindex)
const {
58 if (diskcode < 0 || layercode < 0) {
60 edm::LogPrint(
"Tracklet") <<
"layermask diskmask : " << layermask <<
" " << diskmask;
70 int layerdiskaddress = layercode + (diskcode <<
LayerMemBits_);
72 assert(layerdiskaddress >= 0);
79 edm::LogVerbatim(
"Tracklet") <<
"layermask diskmask : " << layermask <<
" " << diskmask;
112 int layerdiskaddress = layercode + (diskcode <<
LayerMemBits_);
114 assert(layerdiskaddress >= 0);
120 edm::LogPrint(
"Tracklet") <<
"Duplicate entry: layermask=" << layermask <<
" diskmaks=" << diskmask;
131 for (
int i = 0;
i < multiplicity;
i++) {
132 for (
int irinv = 0; irinv < nrinv; irinv++) {
134 tmp.setIndex(layermask, diskmask,
i, irinv);
144 edm::LogVerbatim(
"Tracklet") <<
" flags (good/eof/fail/bad): " <<
in.good() <<
" " <<
in.eof() <<
" " <<
in.fail()
152 in >> layerstr >> diskstr >> multiplicity;
156 if (multiplicity == 8)
158 if (multiplicity == 64)
160 if (multiplicity == 512)
165 if (multiplicity == 8)
167 if (multiplicity == 64)
169 if (multiplicity == 512)
176 char** tmpptr =
nullptr;
178 int layers = strtol(layerstr.c_str(), tmpptr, 2);
179 int disks = strtol(diskstr.c_str(), tmpptr, 2);
188 for (
int i = 0;
i < nentries;
i++) {
193 int irinv = der.
irinv();
201 edm::LogVerbatim(
"Tracklet") <<
"PRINT i " <<
i <<
" " << layermask <<
" " << diskmask <<
" " << alphamask <<
" "
209 if (layermask & (1 << (
N_LAYER - 1 -
l))) {
222 if (diskmask & (3 << (2 * (
N_DISK - 1 -
d)))) {
227 if (diskmask & (1 << (2 * (
N_DISK - 1 -
d)))) {
229 int ialpha = alphamask & 7;
230 alphamask = alphamask >> 3;
236 int ialpha = alphamask & 3;
237 alphamask = alphamask >> 2;
241 int ialpha = alphamask & 1;
242 alphamask = alphamask >> 1;
261 edm::LogVerbatim(
"Tracklet") <<
"PRINT ndisks alpha[0] z[0] t: " << ndisks <<
" " <<
alpha[0] <<
" " << z[0]
263 for (
int iii = 0; iii <
nlayers; iii++) {
268 calculateDerivatives(
settings_,
nlayers,
r, ndisks, z,
alpha,
t,
rinv,
D, iD, MinvDt, iMinvDt, sigma,
kfactor);
279 settings_,
nlayers,
r, ndisks, z,
alpha,
t,
rinv,
D, iD, MinvDtDelta, iMinvDt, sigma,
kfactor);
284 double tder = (MinvDtDelta[2][2 *
ii + 1] - MinvDt[2][2 *
ii + 1]) /
delta;
287 double zder = (MinvDtDelta[3][2 *
ii + 1] - MinvDt[3][2 *
ii + 1]) /
delta;
300 edm::LogVerbatim(
"Tracklet") <<
"iMinvDt table build : " << iMinvDt[0][10] <<
" " << iMinvDt[1][10] <<
" "
301 << iMinvDt[2][10] <<
" " << iMinvDt[3][10] <<
" " <<
t <<
" " <<
nlayers <<
" "
305 for (
int iii = 0; iii < ndisks; iii++) {
307 oss += std::to_string(
alpha[iii]);
311 for (
int iii = 0; iii < ndisks; iii++) {
313 oss += std::to_string(z[iii]);
336 edm::LogVerbatim(
"Tracklet") <<
"PRINT i " <<
i <<
" " <<
j <<
" " << iMinvDt[1][2 *
j] <<
" "
368 tmp1 = (1 << nbits) - 1;
369 tmp.set(tmp1, nbits,
true, __LINE__, __FILE__);
370 outL <<
tmp.str() << endl;
379 tmp1 = (1 << nbits) - 1;
381 tmp.set(tmp1, nbits,
true, __LINE__, __FILE__);
382 outD <<
tmp.str() << endl;
391 tmp1 = (1 << nbits) - 1;
393 tmp.set(tmp1, nbits,
true, __LINE__, __FILE__);
394 outLD <<
tmp.str() << endl;
398 const std::array<string, N_TRKLSEED> seedings = {{
"L1L2",
"L3L4",
"L5L6",
"D1D2",
"D3D4",
"D1L1",
"D1L2"}};
405 const string fname =
prefix +
"Rinvdphi_" + seedings[
i] +
".tab";
406 outrinvdphi[
i].open(
fname);
407 if (outrinvdphi[
i].fail())
408 throw cms::Exception(
"BadFile") << __FILE__ <<
" " << __LINE__ <<
" could not create file " <<
fname;
413 const string fname =
prefix +
"Rinvdzordr_" + seedings[
i] +
".tab";
414 outrinvdzordr[
i].open(
fname);
415 if (outrinvdzordr[
i].fail())
416 throw cms::Exception(
"BadFile") << __FILE__ <<
" " << __LINE__ <<
" could not create file " <<
fname;
421 const string fname =
prefix +
"Phi0dphi_" + seedings[
i] +
".tab";
422 outphi0dphi[
i].open(
fname);
423 if (outphi0dphi[
i].fail())
424 throw cms::Exception(
"BadFile") << __FILE__ <<
" " << __LINE__ <<
" could not create file " <<
fname;
429 const string fname =
prefix +
"Phi0dzordr_" + seedings[
i] +
".tab";
430 outphi0dzordr[
i].open(
fname);
431 if (outphi0dzordr[
i].fail())
432 throw cms::Exception(
"BadFile") << __FILE__ <<
" " << __LINE__ <<
" could not create file " <<
fname;
437 const string fname =
prefix +
"Tdphi_" + seedings[
i] +
".tab";
439 if (outtdphi[
i].fail())
440 throw cms::Exception(
"BadFile") << __FILE__ <<
" " << __LINE__ <<
" could not create file " <<
fname;
445 const string fname =
prefix +
"Tdzordr_" + seedings[
i] +
".tab";
446 outtdzordr[
i].open(
fname);
447 if (outtdzordr[
i].fail())
448 throw cms::Exception(
"BadFile") << __FILE__ <<
" " << __LINE__ <<
" could not create file " <<
fname;
453 const string fname =
prefix +
"Z0dphi_" + seedings[
i] +
".tab";
455 if (outz0dphi[
i].fail())
456 throw cms::Exception(
"BadFile") << __FILE__ <<
" " << __LINE__ <<
" could not create file " <<
fname;
461 string fname =
prefix +
"Z0dzordr_" + seedings[
i] +
".tab";
462 outz0dzordr[
i].open(
fname);
463 if (outz0dzordr[
i].fail())
464 throw cms::Exception(
"BadFile") << __FILE__ <<
" " << __LINE__ <<
" could not create file " <<
fname;
468 unsigned int layerhits = der.layerMask();
469 unsigned int diskmask = der.diskMask();
470 unsigned int diskhits = 0;
471 if (diskmask & (3 << 8))
473 if (diskmask & (3 << 6))
475 if (diskmask & (3 << 4))
477 if (diskmask & (3 << 2))
479 if (diskmask & (3 << 0))
482 unsigned int hits = (layerhits << 5) + diskhits;
487 for (
const string&
seed : seedings) {
488 unsigned int iseed1 = 0;
489 unsigned int iseed2 = 0;
491 if (
seed ==
"L1L2") {
495 if (
seed ==
"L3L4") {
499 if (
seed ==
"L5L6") {
503 if (
seed ==
"D1D2") {
507 if (
seed ==
"D3D4") {
511 if (
seed ==
"D1L1") {
515 if (
seed ==
"D1L2") {
520 bool goodseed = (
hits & (1 << (11 - iseed1))) and (
hits & (1 << (11 - iseed2)));
522 int itmprinvdphi[
N_PROJ] = {9999999, 9999999, 9999999, 9999999};
523 int itmprinvdzordr[
N_PROJ] = {9999999, 9999999, 9999999, 9999999};
524 int itmpphi0dphi[
N_PROJ] = {9999999, 9999999, 9999999, 9999999};
525 int itmpphi0dzordr[
N_PROJ] = {9999999, 9999999, 9999999, 9999999};
526 int itmptdphi[
N_PROJ] = {9999999, 9999999, 9999999, 9999999};
527 int itmptdzordr[
N_PROJ] = {9999999, 9999999, 9999999, 9999999};
528 int itmpz0dphi[
N_PROJ] = {9999999, 9999999, 9999999, 9999999};
529 int itmpz0dzordr[
N_PROJ] = {9999999, 9999999, 9999999, 9999999};
534 for (
unsigned int ihit = 1; ihit <
N_FITSTUB * 2; ++ihit) {
536 if (ihit == iseed1
or ihit == iseed2) {
541 if (not(
hits & (1 << (11 - ihit))))
545 if (
seed ==
"L1L2") {
546 if (ihit == 3
or ihit == 10)
548 if (ihit == 4
or ihit == 9)
550 if (ihit == 5
or ihit == 8)
552 if (ihit == 6
or ihit == 7)
554 }
else if (
seed ==
"L3L4") {
559 if (ihit == 5
or ihit == 8)
561 if (ihit == 6
or ihit == 7)
563 }
else if (
seed ==
"L5L6") {
572 }
else if (
seed ==
"D1D2") {
579 if (ihit == 2
or ihit == 11)
581 }
else if (
seed ==
"D3D4") {
588 if (ihit == 2
or ihit == 11)
590 }
else if (
seed ==
"D1L1" or "D1L2") {
600 if (inputI >= 0 and inputI < (
int)
N_PROJ) {
601 itmprinvdphi[inputI] = der.irinvdphi(ider);
602 itmprinvdzordr[inputI] = der.irinvdzordr(ider);
603 itmpphi0dphi[inputI] = der.iphi0dphi(ider);
604 itmpphi0dzordr[inputI] = der.iphi0dzordr(ider);
605 itmptdphi[inputI] = der.itdphi(ider);
606 itmptdzordr[inputI] = der.itdzordr(ider);
607 itmpz0dphi[inputI] = der.iz0dphi(ider);
608 itmpz0dzordr[inputI] = der.iz0dzordr(ider);
618 for (
unsigned int j = 0;
j <
N_PROJ; ++
j) {
619 if (itmprinvdphi[
j] > (1 << nbits))
620 itmprinvdphi[
j] = (1 << nbits) - 1;
621 tmprinvdphi[
j].
set(itmprinvdphi[
j], nbits + 1,
false, __LINE__, __FILE__);
623 outrinvdphi[
i] << tmprinvdphi[0].
str() << tmprinvdphi[1].
str() << tmprinvdphi[2].
str() << tmprinvdphi[3].
str()
628 for (
unsigned int j = 0;
j <
N_PROJ; ++
j) {
629 if (itmprinvdzordr[
j] > (1 << nbits))
630 itmprinvdzordr[
j] = (1 << nbits) - 1;
631 tmprinvdzordr[
j].
set(itmprinvdzordr[
j], nbits + 1,
false, __LINE__, __FILE__);
633 outrinvdzordr[
i] << tmprinvdzordr[0].
str() << tmprinvdzordr[1].
str() << tmprinvdzordr[2].
str()
634 << tmprinvdzordr[3].
str() << endl;
638 for (
unsigned int j = 0;
j <
N_PROJ; ++
j) {
639 if (itmpphi0dphi[
j] > (1 << nbits))
640 itmpphi0dphi[
j] = (1 << nbits) - 1;
641 tmpphi0dphi[
j].
set(itmpphi0dphi[
j], nbits + 1,
false, __LINE__, __FILE__);
643 outphi0dphi[
i] << tmpphi0dphi[0].
str() << tmpphi0dphi[1].
str() << tmpphi0dphi[2].
str() << tmpphi0dphi[3].
str()
648 for (
unsigned int j = 0;
j <
N_PROJ; ++
j) {
649 if (itmpphi0dzordr[
j] > (1 << nbits))
650 itmpphi0dzordr[
j] = (1 << nbits) - 1;
651 tmpphi0dzordr[
j].
set(itmpphi0dzordr[
j], nbits + 1,
false, __LINE__, __FILE__);
653 outphi0dzordr[
i] << tmpphi0dzordr[0].
str() << tmpphi0dzordr[1].
str() << tmpphi0dzordr[2].
str()
654 << tmpphi0dzordr[3].
str() << endl;
658 for (
unsigned int j = 0;
j <
N_PROJ; ++
j) {
659 if (itmptdphi[
j] > (1 << nbits))
660 itmptdphi[
j] = (1 << nbits) - 1;
661 tmptdphi[
j].
set(itmptdphi[
j], nbits + 1,
false, __LINE__, __FILE__);
663 outtdphi[
i] << tmptdphi[0].
str() << tmptdphi[1].
str() << tmptdphi[2].
str() << tmptdphi[3].
str() << endl;
667 for (
unsigned int j = 0;
j <
N_PROJ; ++
j) {
668 if (itmptdzordr[
j] > (1 << nbits))
669 itmptdzordr[
j] = (1 << nbits) - 1;
670 tmptdzordr[
j].
set(itmptdzordr[
j], nbits + 1,
false, __LINE__, __FILE__);
672 outtdzordr[
i] << tmptdzordr[0].
str() << tmptdzordr[1].
str() << tmptdzordr[2].
str() << tmptdzordr[3].
str()
677 for (
unsigned int j = 0;
j <
N_PROJ; ++
j) {
678 if (itmpz0dphi[
j] > (1 << nbits))
679 itmpz0dphi[
j] = (1 << nbits) - 1;
680 tmpz0dphi[
j].
set(itmpz0dphi[
j], nbits + 1,
false, __LINE__, __FILE__);
682 outz0dphi[
i] << tmpz0dphi[0].
str() << tmpz0dphi[1].
str() << tmpz0dphi[2].
str() << tmpz0dphi[3].
str() << endl;
686 for (
unsigned int j = 0;
j <
N_PROJ; ++
j) {
687 if (itmpz0dzordr[
j] > (1 << nbits))
688 itmpz0dzordr[
j] = (1 << nbits) - 1;
689 tmpz0dzordr[
j].
set(itmpz0dzordr[
j], nbits + 1,
false, __LINE__, __FILE__);
691 outz0dzordr[
i] << tmpz0dzordr[0].
str() << tmpz0dzordr[1].
str() << tmpz0dzordr[2].
str() << tmpz0dzordr[3].
str()
701 outrinvdphi[
i].close();
702 outrinvdzordr[
i].close();
703 outphi0dphi[
i].close();
704 outphi0dzordr[
i].close();
706 outtdzordr[
i].close();
707 outz0dphi[
i].close();
708 outz0dzordr[
i].close();
717 unsigned int i,
j,
k;
720 for (
i = 0;
i <
n;
i++) {
721 for (
j =
n;
j < 2 *
n;
j++) {
729 for (
i = 0;
i <
n;
i++) {
730 for (
j = 0;
j <
n;
j++) {
733 for (
k = 0;
k < 2 *
n;
k++) {
740 for (
i = 0;
i <
n;
i++) {
742 for (
j = 0;
j < 2 *
n;
j++) {
752 unsigned int i,
j,
k;
755 for (
i = 0;
i <
n;
i++) {
756 for (
j =
n;
j < 2 *
n;
j++) {
764 for (
i = 0;
i <
n;
i++) {
765 for (
j = 0;
j <
n;
j++) {
768 for (
k = 0;
k < 2 *
n;
k++) {
775 for (
i = 0;
i <
n;
i++) {
777 for (
j = 0;
j < 2 *
n;
j++) {
801 double sigmazpsbarrel = sigmaz;
803 sigmazpsbarrel = sigmaz *
std::abs(
t) / 2.0;
828 D[0][
j] = -0.5 * ri * ri /
sqrt(1 - 0.25 * ri * ri *
rinv *
rinv) / sigmax;
829 D[1][
j] = ri / sigmax;
838 if (ri < settings.
rPS2S()) {
839 D[2][
j] = (2 /
rinv) * asin(0.5 * ri *
rinv) / sigmazpsbarrel;
840 D[3][
j] = 1.0 / sigmazpsbarrel;
841 sigma[
j] = sigmazpsbarrel;
844 D[2][
j] = (2 /
rinv) * asin(0.5 * ri *
rinv) / sigmaz2;
845 D[3][
j] = 1.0 / sigmaz2;
853 for (
unsigned int i = 0;
i < ndisks;
i++) {
858 double rmultiplier =
alpha[
i] * zi /
t;
860 double phimultiplier = zi /
t;
865 double drdt = -(zi -
z0) *
cos(0.5 *
rinv * (zi -
z0) /
t) / (
t *
t);
866 double drdz0 = -
cos(0.5 *
rinv * (zi -
z0) /
t) /
t;
868 double dphidrinv = -0.5 * (zi -
z0) /
t;
869 double dphidphi0 = 1.0;
870 double dphidt = 0.5 *
rinv * (zi -
z0) / (
t *
t);
871 double dphidz0 = 0.5 *
rinv /
t;
873 double r = (zi -
z0) /
t;
877 sigma[
j] = sigmax2sdisk;
879 sigma[
j] = sigmaxpsdisk;
882 D[0][
j] = (phimultiplier * dphidrinv + rmultiplier * drdrinv) / sigma[
j];
883 D[1][
j] = (phimultiplier * dphidphi0 + rmultiplier * drdphi0) / sigma[
j];
884 D[2][
j] = (phimultiplier * dphidt + rmultiplier * drdt) / sigma[
j];
885 D[3][
j] = (phimultiplier * dphidz0 + rmultiplier * drdz0) / sigma[
j];
891 D[0][
j] = drdrinv / sigmazpsdisk;
892 D[1][
j] = drdphi0 / sigmazpsdisk;
893 D[2][
j] = drdt / sigmazpsdisk;
894 D[3][
j] = drdz0 / sigmazpsdisk;
895 sigma[
j] = sigmazpsdisk;
898 D[0][
j] = drdrinv / sigmaz2sdisk;
899 D[1][
j] = drdphi0 / sigmaz2sdisk;
900 D[2][
j] = drdt / sigmaz2sdisk;
901 D[3][
j] = drdz0 / sigmaz2sdisk;
902 sigma[
j] = sigmaz2sdisk;
911 for (
unsigned int i1 = 0;
i1 < 4;
i1++) {
912 for (
unsigned int i2 = 0;
i2 < 4;
i2++) {
914 for (
unsigned int j = 0;
j < 2 *
n;
j++) {
929 for (
unsigned int j = 0;
j < 2 *
n;
j++) {
930 for (
unsigned int i1 = 0;
i1 < 4;
i1++) {
931 for (
unsigned int i2 = 0;
i2 < 4;
i2++) {
937 for (
unsigned int i = 0;
i <
n;
i++) {
958 MinvDt[0][2 *
i] *= rnew[
i] / sigmax;
959 MinvDt[1][2 *
i] *= rnew[
i] / sigmax;
960 MinvDt[2][2 *
i] *= rnew[
i] / sigmax;
961 MinvDt[3][2 *
i] *= rnew[
i] / sigmax;
970 if (rnew[
i] < settings.
rPS2S()) {
971 MinvDt[0][2 *
i + 1] /= sigmazpsbarrel;
972 MinvDt[1][2 *
i + 1] /= sigmazpsbarrel;
973 MinvDt[2][2 *
i + 1] /= sigmazpsbarrel;
974 MinvDt[3][2 *
i + 1] /= sigmazpsbarrel;
976 iMinvDt[0][2 *
i + 1] =
978 iMinvDt[1][2 *
i + 1] =
980 iMinvDt[2][2 *
i + 1] =
982 iMinvDt[3][2 *
i + 1] =
985 MinvDt[0][2 *
i + 1] /= sigmaz2;
986 MinvDt[1][2 *
i + 1] /= sigmaz2;
987 MinvDt[2][2 *
i + 1] /= sigmaz2;
988 MinvDt[3][2 *
i + 1] /= sigmaz2;
992 iMinvDt[0][2 *
i + 1] =
994 iMinvDt[1][2 *
i + 1] =
996 iMinvDt[2][2 *
i + 1] =
998 iMinvDt[3][2 *
i + 1] =
1007 MinvDt[0][2 *
i] *= (rnew[
i] /
denom);
1008 MinvDt[1][2 *
i] *= (rnew[
i] /
denom);
1009 MinvDt[2][2 *
i] *= (rnew[
i] /
denom);
1010 MinvDt[3][2 *
i] *= (rnew[
i] /
denom);
1012 assert(MinvDt[0][2 *
i] == MinvDt[0][2 *
i]);
1017 iMinvDt[3][2 *
i] = (1 << settings.
fitz0bitshift()) * MinvDt[3][2 *
i] * settings.
kphi() / settings.
kz();
1021 MinvDt[0][2 *
i + 1] /=
denom;
1022 MinvDt[1][2 *
i + 1] /=
denom;
1023 MinvDt[2][2 *
i + 1] /=
denom;
1024 MinvDt[3][2 *
i + 1] /=
denom;
1026 iMinvDt[0][2 *
i + 1] =
1028 iMinvDt[1][2 *
i + 1] =
1030 iMinvDt[2][2 *
i + 1] =
1032 iMinvDt[3][2 *
i + 1] =
1042 double tmax = 1000.0;
1046 if (diskmask & (1 << (2 * (5 -
d) + 1))) {
1047 double dmax = settings.
zmean(
d - 1) / 22.0;
1048 if (dmax > sinh(2.4))
1050 double dmin = settings.
zmean(
d - 1) / 65.0;
1057 if (diskmask & (1 << (2 * (5 -
d)))) {
1058 double dmax = settings.
zmean(
d - 1) / 65.0;
1059 double dmin = settings.
zmean(
d - 1) / 105.0;
1068 if (layermask & (1 << (6 -
l))) {