File indexing completed on 2023-10-25 09:56:30
0001 #include "L1Trigger/TrackFindingTracklet/interface/TrackletLUT.h"
0002 #include "L1Trigger/TrackFindingTracklet/interface/Util.h"
0003 #include "L1Trigger/TrackFindingTracklet/interface/Settings.h"
0004 #include "L1Trigger/TrackFindingTracklet/interface/TrackletConfigBuilder.h"
0005 #include "L1Trigger/TrackTrigger/interface/Setup.h"
0006
0007 #include <filesystem>
0008
0009 using namespace std;
0010 using namespace trklet;
0011
0012 TrackletLUT::TrackletLUT(const Settings& settings) : settings_(settings), setup_(settings.setup()) {}
0013
0014 std::vector<const tt::SensorModule*> TrackletLUT::getSensorModules(
0015 unsigned int layerdisk, bool isPS, std::array<double, 2> tan_range, unsigned int nzbins, unsigned int zbin) {
0016
0017
0018
0019
0020
0021
0022 bool use_tan_range = !(tan_range[0] == -1 and tan_range[1] == -1);
0023 bool use_zbins = (nzbins > 1);
0024
0025 bool barrel = layerdisk < N_LAYER;
0026
0027 int layerId = barrel ? layerdisk + 1 : layerdisk + N_LAYER - 1;
0028
0029 std::vector<const tt::SensorModule*> sensorModules;
0030
0031 double z0 = settings_.z0cut();
0032
0033 for (auto& sm : setup_->sensorModules()) {
0034 if (sm.layerId() != layerId || sm.z() < 0 || sm.psModule() != isPS) {
0035 continue;
0036 }
0037
0038 if (use_tan_range) {
0039 const double term = (sm.numColumns() / 2 - 0.5) * sm.pitchCol();
0040 double rmin = sm.r() - term * std::abs(sm.sinTilt());
0041 double rmax = sm.r() + term * std::abs(sm.sinTilt());
0042
0043 double zmin = std::abs(sm.z()) - term * std::abs(sm.cosTilt());
0044 double zmax = std::abs(sm.z()) + term * std::abs(sm.cosTilt());
0045
0046
0047 double mod_tan_max = tan_theta(rmin, zmax, z0, false);
0048 double mod_tan_min = tan_theta(rmax, zmin, z0, true);
0049
0050 if (mod_tan_max >= tan_range[0] && mod_tan_min <= tan_range[1]) {
0051 sensorModules.push_back(&sm);
0052 }
0053 } else if (use_zbins) {
0054 assert(layerdisk < 3);
0055
0056 if (nzbins == 2) {
0057 bool useFlat = (zbin == 0);
0058 bool isFlat = (sm.tilt() == 0);
0059
0060 if (useFlat and isFlat)
0061 sensorModules.push_back(&sm);
0062 else if (!useFlat and !isFlat)
0063 sensorModules.push_back(&sm);
0064 } else if (nzbins == 13) {
0065 if (sm.ringId(setup_) == zbin)
0066 sensorModules.push_back(&sm);
0067 } else {
0068 throw cms::Exception("Unspecified number of z bins");
0069 }
0070 } else {
0071 sensorModules.push_back(&sm);
0072 }
0073 }
0074
0075
0076 static constexpr double delta = 1.e-3;
0077 auto smallerR = [](const tt::SensorModule* lhs, const tt::SensorModule* rhs) { return lhs->r() < rhs->r(); };
0078 auto smallerZ = [](const tt::SensorModule* lhs, const tt::SensorModule* rhs) { return lhs->z() < rhs->z(); };
0079 auto equalRZ = [](const tt::SensorModule* lhs, const tt::SensorModule* rhs) {
0080 return abs(lhs->r() - rhs->r()) < delta && abs(lhs->z() - rhs->z()) < delta;
0081 };
0082 stable_sort(sensorModules.begin(), sensorModules.end(), smallerR);
0083 stable_sort(sensorModules.begin(), sensorModules.end(), smallerZ);
0084 sensorModules.erase(unique(sensorModules.begin(), sensorModules.end(), equalRZ), sensorModules.end());
0085
0086 return sensorModules;
0087 }
0088
0089 std::array<double, 2> TrackletLUT::getTanRange(const std::vector<const tt::SensorModule*>& sensorModules) {
0090
0091
0092 std::array<double, 2> tan_range = {{2147483647, 0}};
0093
0094 double z0 = settings_.z0cut();
0095
0096 for (auto sm : sensorModules) {
0097 const double term = (sm->numColumns() / 2 - 0.5) * sm->pitchCol();
0098 double rmin = sm->r() - term * std::abs(sm->sinTilt());
0099 double rmax = sm->r() + term * std::abs(sm->sinTilt());
0100
0101 double zmin = std::abs(sm->z()) - term * sm->cosTilt();
0102 double zmax = std::abs(sm->z()) + term * sm->cosTilt();
0103
0104 double mod_tan_max = tan_theta(rmin, zmax, z0, true);
0105 double mod_tan_min = tan_theta(rmax, zmin, z0, false);
0106
0107 if (mod_tan_min < tan_range[0])
0108 tan_range[0] = mod_tan_min;
0109 if (mod_tan_max > tan_range[1])
0110 tan_range[1] = mod_tan_max;
0111 }
0112 return tan_range;
0113 }
0114
0115 std::vector<std::array<double, 2>> TrackletLUT::getBendCut(unsigned int layerdisk,
0116 const std::vector<const tt::SensorModule*>& sensorModules,
0117 bool isPS,
0118 double FEbendcut) {
0119
0120
0121
0122
0123
0124 unsigned int bendbits = isPS ? 3 : 4;
0125
0126 std::vector<std::array<double, 2>> bendpars;
0127 std::vector<std::array<double, 2>> bendminmax;
0128
0129
0130 for (int i = 0; i < 1 << bendbits; i++) {
0131 bendpars.push_back({{99, 0}});
0132 bendminmax.push_back({{99, -99}});
0133 }
0134
0135
0136 for (auto sm : sensorModules) {
0137 int window = sm->windowSize();
0138 const vector<double>& encodingBend = setup_->encodingBend(window, isPS);
0139
0140
0141 for (int ibend = 0; ibend <= 2 * window; ibend++) {
0142 int FEbend = ibend - window;
0143 double BEbend = setup_->stubAlgorithm()->degradeBend(isPS, window, FEbend);
0144
0145 const auto pos = std::find(encodingBend.begin(), encodingBend.end(), std::abs(BEbend));
0146 int bend = std::signbit(BEbend) ? (1 << bendbits) - distance(encodingBend.begin(), pos)
0147 : distance(encodingBend.begin(), pos);
0148
0149 double bendmin = FEbend / 2.0 - FEbendcut;
0150 double bendmax = FEbend / 2.0 + FEbendcut;
0151
0152
0153 double z_mod[2];
0154 double r_mod[2];
0155
0156 z_mod[0] = std::abs(sm->z()) + (sm->numColumns() / 2 - 0.5) * sm->pitchCol() * sm->cosTilt();
0157 z_mod[1] = std::abs(sm->z()) - (sm->numColumns() / 2 - 0.5) * sm->pitchCol() * sm->cosTilt();
0158
0159 r_mod[0] = sm->r() - (sm->numColumns() / 2 - 0.5) * sm->pitchCol() * std::abs(sm->sinTilt());
0160 r_mod[1] = sm->r() + (sm->numColumns() / 2 - 0.5) * sm->pitchCol() * std::abs(sm->sinTilt());
0161
0162 for (int i = 0; i < 2; i++) {
0163 double CF = std::abs(sm->sinTilt()) * (z_mod[i] / r_mod[i]) + sm->cosTilt();
0164
0165 double cbendmin =
0166 convertFEBend(bendmin, sm->sep(), settings_.sensorSpacing2S(), CF, (layerdisk < N_LAYER), r_mod[i]);
0167 double cbendmax =
0168 convertFEBend(bendmax, sm->sep(), settings_.sensorSpacing2S(), CF, (layerdisk < N_LAYER), r_mod[i]);
0169
0170 if (cbendmin < bendminmax[bend][0])
0171 bendminmax.at(bend)[0] = cbendmin;
0172 if (cbendmax > bendminmax[bend][1])
0173 bendminmax.at(bend)[1] = cbendmax;
0174 }
0175 }
0176 }
0177
0178 for (int i = 0; i < 1 << bendbits; i++) {
0179 double mid = (bendminmax[i][1] + bendminmax[i][0]) / 2;
0180 double cut = (bendminmax[i][1] - bendminmax[i][0]) / 2;
0181
0182 bendpars[i][0] = mid;
0183 bendpars[i][1] = cut;
0184 }
0185
0186 return bendpars;
0187 }
0188
0189 void TrackletLUT::initmatchcut(unsigned int layerdisk, MatchType type, unsigned int region) {
0190 char cregion = 'A' + region;
0191
0192 for (unsigned int iSeed = 0; iSeed < N_SEED; iSeed++) {
0193 if (type == barrelphi) {
0194 table_.push_back(settings_.rphimatchcut(iSeed, layerdisk) / (settings_.kphi1() * settings_.rmean(layerdisk)));
0195 }
0196 if (type == barrelz) {
0197 table_.push_back(settings_.zmatchcut(iSeed, layerdisk) / settings_.kz());
0198 }
0199 if (type == diskPSphi) {
0200 table_.push_back(settings_.rphicutPS(iSeed, layerdisk - N_LAYER) / (settings_.kphi() * settings_.kr()));
0201 }
0202 if (type == disk2Sphi) {
0203 table_.push_back(settings_.rphicut2S(iSeed, layerdisk - N_LAYER) / (settings_.kphi() * settings_.kr()));
0204 }
0205 if (type == disk2Sr) {
0206 table_.push_back(settings_.rcut2S(iSeed, layerdisk - N_LAYER) / settings_.krprojshiftdisk());
0207 }
0208 if (type == diskPSr) {
0209 table_.push_back(settings_.rcutPS(iSeed, layerdisk - N_LAYER) / settings_.krprojshiftdisk());
0210 }
0211 }
0212 if (type == alphainner) {
0213 for (unsigned int i = 0; i < N_DSS_MOD * 2; i++) {
0214 table_.push_back((1 << settings_.alphashift()) * settings_.krprojshiftdisk() * settings_.half2SmoduleWidth() /
0215 (1 << (settings_.nbitsalpha() - 1)) / (settings_.rDSSinner(i) * settings_.rDSSinner(i)) /
0216 settings_.kphi());
0217 }
0218 }
0219 if (type == alphaouter) {
0220 for (unsigned int i = 0; i < N_DSS_MOD * 2; i++) {
0221 table_.push_back((1 << settings_.alphashift()) * settings_.krprojshiftdisk() * settings_.half2SmoduleWidth() /
0222 (1 << (settings_.nbitsalpha() - 1)) / (settings_.rDSSouter(i) * settings_.rDSSouter(i)) /
0223 settings_.kphi());
0224 }
0225 }
0226 if (type == rSSinner) {
0227 for (unsigned int i = 0; i < N_DSS_MOD * 2; i++) {
0228 table_.push_back(settings_.rDSSinner(i) / settings_.kr());
0229 }
0230 }
0231 if (type == rSSouter) {
0232 for (unsigned int i = 0; i < N_DSS_MOD * 2; i++) {
0233 table_.push_back(settings_.rDSSouter(i) / settings_.kr());
0234 }
0235 }
0236
0237 name_ = settings_.combined() ? "MP_" : "MC_";
0238
0239 if (type == barrelphi) {
0240 name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_phicut.tab";
0241 }
0242 if (type == barrelz) {
0243 name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_zcut.tab";
0244 }
0245 if (type == diskPSphi) {
0246 name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_PSphicut.tab";
0247 }
0248 if (type == disk2Sphi) {
0249 name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_2Sphicut.tab";
0250 }
0251 if (type == disk2Sr) {
0252 name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_2Srcut.tab";
0253 }
0254 if (type == diskPSr) {
0255 name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_PSrcut.tab";
0256 }
0257 if (type == alphainner) {
0258 name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_alphainner.tab";
0259 }
0260 if (type == alphaouter) {
0261 name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_alphaouter.tab";
0262 }
0263 if (type == rSSinner) {
0264 name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_rDSSinner.tab";
0265 }
0266 if (type == rSSouter) {
0267 name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_rDSSouter.tab";
0268 }
0269
0270 positive_ = false;
0271
0272 writeTable();
0273 }
0274
0275 void TrackletLUT::initTPlut(bool fillInner,
0276 unsigned int iSeed,
0277 unsigned int layerdisk1,
0278 unsigned int layerdisk2,
0279 unsigned int nbitsfinephidiff,
0280 unsigned int iTP) {
0281
0282 int nfinephibins = settings_.nallstubs(layerdisk2) * settings_.nvmte(1, iSeed) * (1 << settings_.nfinephi(1, iSeed));
0283 double dfinephi = settings_.dphisectorHG() / nfinephibins;
0284
0285 int outerrbits = 3;
0286
0287 if (iSeed == Seed::L1L2 || iSeed == Seed::L2L3 || iSeed == Seed::L3L4 || iSeed == Seed::L5L6) {
0288 outerrbits = 0;
0289 }
0290
0291 int outerrbins = (1 << outerrbits);
0292
0293 double dphi[2];
0294 double router[2];
0295
0296 bool isPSinner;
0297 bool isPSouter;
0298
0299 if (iSeed == Seed::L3L4) {
0300 isPSinner = true;
0301 isPSouter = false;
0302 } else if (iSeed == Seed::L5L6) {
0303 isPSinner = false;
0304 isPSouter = false;
0305 } else {
0306 isPSinner = true;
0307 isPSouter = true;
0308 }
0309
0310 unsigned int nbendbitsinner = isPSinner ? N_BENDBITS_PS : N_BENDBITS_2S;
0311 unsigned int nbendbitsouter = isPSouter ? N_BENDBITS_PS : N_BENDBITS_2S;
0312
0313 double z0 = settings_.z0cut();
0314
0315 int nbinsfinephidiff = (1 << nbitsfinephidiff);
0316
0317 for (int iphibin = 0; iphibin < nbinsfinephidiff; iphibin++) {
0318 int iphidiff = iphibin;
0319 if (iphibin >= nbinsfinephidiff / 2) {
0320 iphidiff = iphibin - nbinsfinephidiff;
0321 }
0322
0323
0324 double deltaphi = 1.5;
0325 dphi[0] = (iphidiff - deltaphi) * dfinephi;
0326 dphi[1] = (iphidiff + deltaphi) * dfinephi;
0327 for (int irouterbin = 0; irouterbin < outerrbins; irouterbin++) {
0328 if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4 || iSeed == Seed::L1D1 || iSeed == Seed::L2D1) {
0329 router[0] =
0330 settings_.rmindiskvm() + irouterbin * (settings_.rmaxdiskvm() - settings_.rmindiskvm()) / outerrbins;
0331 router[1] =
0332 settings_.rmindiskvm() + (irouterbin + 1) * (settings_.rmaxdiskvm() - settings_.rmindiskvm()) / outerrbins;
0333 } else {
0334 router[0] = settings_.rmean(layerdisk2);
0335 router[1] = settings_.rmean(layerdisk2);
0336 }
0337
0338
0339 std::vector<std::array<double, 2>> bend_cuts_inner;
0340 std::vector<std::array<double, 2>> bend_cuts_outer;
0341
0342 if (settings_.useCalcBendCuts) {
0343 std::vector<const tt::SensorModule*> sminner;
0344 std::vector<const tt::SensorModule*> smouter;
0345
0346 if (iSeed == Seed::L1L2 || iSeed == Seed::L2L3 || iSeed == Seed::L3L4 || iSeed == Seed::L5L6) {
0347 double outer_tan_max = tan_theta(settings_.rmean(layerdisk2), settings_.zlength(), z0, true);
0348 std::array<double, 2> tan_range = {{0, outer_tan_max}};
0349
0350 smouter = getSensorModules(layerdisk2, isPSouter, tan_range);
0351 sminner = getSensorModules(layerdisk1, isPSinner, tan_range);
0352
0353 } else if (iSeed == Seed::L1D1 || iSeed == Seed::L2D1) {
0354 double outer_tan_min = tan_theta(router[1], settings_.zmindisk(layerdisk2 - N_LAYER), z0, false);
0355 double outer_tan_max = tan_theta(router[0], settings_.zmaxdisk(layerdisk2 - N_LAYER), z0, true);
0356
0357 smouter = getSensorModules(layerdisk2, isPSouter, {{outer_tan_min, outer_tan_max}});
0358 std::array<double, 2> tan_range = getTanRange(smouter);
0359 sminner = getSensorModules(layerdisk1, isPSinner, tan_range);
0360
0361 } else {
0362
0363 double outer_tan_min = tan_theta(router[1], settings_.zmindisk(layerdisk2 - N_LAYER), z0, false);
0364 double outer_tan_max = tan_theta(router[0], settings_.zmaxdisk(layerdisk2 - N_LAYER), z0, true);
0365
0366 smouter = getSensorModules(layerdisk2, isPSouter, {{outer_tan_min, outer_tan_max}});
0367
0368 std::array<double, 2> tan_range = getTanRange(smouter);
0369 sminner = getSensorModules(layerdisk1, isPSinner, tan_range);
0370 }
0371
0372 bend_cuts_inner = getBendCut(layerdisk1, sminner, isPSinner, settings_.bendcutTE(iSeed, true));
0373 bend_cuts_outer = getBendCut(layerdisk2, smouter, isPSouter, settings_.bendcutTE(iSeed, false));
0374
0375 } else {
0376 for (int ibend = 0; ibend < (1 << nbendbitsinner); ibend++) {
0377 double mid = settings_.benddecode(ibend, layerdisk1, isPSinner);
0378 double cut = settings_.bendcutte(ibend, layerdisk1, isPSinner);
0379 bend_cuts_inner.push_back({{mid, cut}});
0380 }
0381 for (int ibend = 0; ibend < (1 << nbendbitsouter); ibend++) {
0382 double mid = settings_.benddecode(ibend, layerdisk2, isPSouter);
0383 double cut = settings_.bendcutte(ibend, layerdisk2, isPSouter);
0384 bend_cuts_outer.push_back({{mid, cut}});
0385 }
0386 }
0387
0388 double bendinnermin = 20.0;
0389 double bendinnermax = -20.0;
0390 double bendoutermin = 20.0;
0391 double bendoutermax = -20.0;
0392 double rinvmin = 1.0;
0393 double rinvmax = -1.0;
0394 double absrinvmin = 1.0;
0395
0396 for (int i2 = 0; i2 < 2; i2++) {
0397 for (int i3 = 0; i3 < 2; i3++) {
0398 double rinner = 0.0;
0399 if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4) {
0400 rinner = router[i3] * settings_.zmean(layerdisk1 - N_LAYER) / settings_.zmean(layerdisk2 - N_LAYER);
0401 } else {
0402 rinner = settings_.rmean(layerdisk1);
0403 }
0404 if (settings_.useCalcBendCuts) {
0405 if (rinner >= router[i3])
0406 continue;
0407 }
0408 double rinv1 = (rinner < router[i3]) ? rinv(0.0, -dphi[i2], rinner, router[i3]) : 20.0;
0409 double pitchinner = (rinner < settings_.rcrit()) ? settings_.stripPitch(true) : settings_.stripPitch(false);
0410 double pitchouter =
0411 (router[i3] < settings_.rcrit()) ? settings_.stripPitch(true) : settings_.stripPitch(false);
0412 double abendinner = bendstrip(rinner, rinv1, pitchinner, settings_.sensorSpacing2S());
0413 double abendouter = bendstrip(router[i3], rinv1, pitchouter, settings_.sensorSpacing2S());
0414 if (abendinner < bendinnermin)
0415 bendinnermin = abendinner;
0416 if (abendinner > bendinnermax)
0417 bendinnermax = abendinner;
0418 if (abendouter < bendoutermin)
0419 bendoutermin = abendouter;
0420 if (abendouter > bendoutermax)
0421 bendoutermax = abendouter;
0422 if (std::abs(rinv1) < absrinvmin)
0423 absrinvmin = std::abs(rinv1);
0424 if (rinv1 > rinvmax)
0425 rinvmax = rinv1;
0426 if (rinv1 < rinvmin)
0427 rinvmin = rinv1;
0428 }
0429 }
0430
0431 bool passptcut;
0432 double bendfac;
0433 double rinvcutte = settings_.rinvcutte();
0434
0435 if (settings_.useCalcBendCuts) {
0436 double lowrinvcutte =
0437 rinvcutte / 3;
0438 passptcut = rinvmin < rinvcutte and rinvmax > -rinvcutte;
0439 bendfac = (rinvmin < lowrinvcutte and rinvmax > -lowrinvcutte)
0440 ? 1.05
0441 : 1.0;
0442 } else {
0443 passptcut = absrinvmin < rinvcutte;
0444 bendfac = 1.0;
0445 }
0446
0447 if (fillInner) {
0448 for (int ibend = 0; ibend < (1 << nbendbitsinner); ibend++) {
0449 double bendminfac = (isPSinner and (ibend == 2 or ibend == 3)) ? bendfac : 1.0;
0450 double bendmaxfac = (isPSinner and (ibend == 6 or ibend == 5)) ? bendfac : 1.0;
0451
0452 double mid = bend_cuts_inner.at(ibend)[0];
0453 double cut = bend_cuts_inner.at(ibend)[1];
0454
0455 bool passinner = mid + cut * bendmaxfac > bendinnermin && mid - cut * bendminfac < bendinnermax;
0456
0457 table_.push_back(passinner && passptcut);
0458 }
0459 } else {
0460 for (int ibend = 0; ibend < (1 << nbendbitsouter); ibend++) {
0461 double bendminfac = (isPSouter and (ibend == 2 or ibend == 3)) ? bendfac : 1.0;
0462 double bendmaxfac = (isPSouter and (ibend == 6 or ibend == 5)) ? bendfac : 1.0;
0463
0464 double mid = bend_cuts_outer.at(ibend)[0];
0465 double cut = bend_cuts_outer.at(ibend)[1];
0466
0467 bool passouter = mid + cut * bendmaxfac > bendoutermin && mid - cut * bendminfac < bendoutermax;
0468
0469 table_.push_back(passouter && passptcut);
0470 }
0471 }
0472 }
0473 }
0474
0475 nbits_ = 8;
0476
0477 positive_ = false;
0478 char cTP = 'A' + iTP;
0479
0480 name_ = "TP_" + TrackletConfigBuilder::LayerName(layerdisk1) + TrackletConfigBuilder::LayerName(layerdisk2) + cTP;
0481
0482 if (fillInner) {
0483 name_ += "_stubptinnercut.tab";
0484 } else {
0485 name_ += "_stubptoutercut.tab";
0486 }
0487
0488 writeTable();
0489 }
0490
0491 void TrackletLUT::initTPregionlut(unsigned int iSeed,
0492 unsigned int layerdisk1,
0493 unsigned int layerdisk2,
0494 unsigned int iAllStub,
0495 unsigned int nbitsfinephidiff,
0496 unsigned int nbitsfinephi,
0497 const TrackletLUT& tplutinner,
0498 unsigned int iTP) {
0499 int nirbits = 0;
0500 if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4 || iSeed == Seed::L1D1 || iSeed == Seed::L2D1) {
0501 nirbits = 3;
0502 }
0503
0504 unsigned int nbendbitsinner = 3;
0505
0506 if (iSeed == Seed::L5L6) {
0507 nbendbitsinner = 4;
0508 }
0509
0510 for (int innerfinephi = 0; innerfinephi < (1 << nbitsfinephi); innerfinephi++) {
0511 for (int innerbend = 0; innerbend < (1 << nbendbitsinner); innerbend++) {
0512 for (int ir = 0; ir < (1 << nirbits); ir++) {
0513 unsigned int usereg = 0;
0514 for (unsigned int ireg = 0; ireg < settings_.nvmte(1, iSeed); ireg++) {
0515 bool match = false;
0516 for (int ifinephiouter = 0; ifinephiouter < (1 << settings_.nfinephi(1, iSeed)); ifinephiouter++) {
0517 int outerfinephi = iAllStub * (1 << (nbitsfinephi - settings_.nbitsallstubs(layerdisk2))) +
0518 ireg * (1 << settings_.nfinephi(1, iSeed)) + ifinephiouter;
0519 int idphi = outerfinephi - innerfinephi;
0520 bool inrange = (idphi < (1 << (nbitsfinephidiff - 1))) && (idphi >= -(1 << (nbitsfinephidiff - 1)));
0521 if (idphi < 0)
0522 idphi = idphi + (1 << nbitsfinephidiff);
0523 int idphi1 = idphi;
0524 if (iSeed >= 4)
0525 idphi1 = (idphi << 3) + ir;
0526 int ptinnerindexnew = (idphi1 << nbendbitsinner) + innerbend;
0527 match = match || (inrange && tplutinner.lookup(ptinnerindexnew));
0528 }
0529 if (match) {
0530 usereg = usereg | (1 << ireg);
0531 }
0532 }
0533
0534 table_.push_back(usereg);
0535 }
0536 }
0537 }
0538
0539 positive_ = false;
0540 char cTP = 'A' + iTP;
0541
0542 name_ = "TP_" + TrackletConfigBuilder::LayerName(layerdisk1) + TrackletConfigBuilder::LayerName(layerdisk2) + cTP +
0543 "_usereg.tab";
0544
0545 writeTable();
0546 }
0547
0548 void TrackletLUT::initteptlut(bool fillInner,
0549 bool fillTEMem,
0550 unsigned int iSeed,
0551 unsigned int layerdisk1,
0552 unsigned int layerdisk2,
0553 unsigned int innerphibits,
0554 unsigned int outerphibits,
0555 double innerphimin,
0556 double innerphimax,
0557 double outerphimin,
0558 double outerphimax,
0559 const std::string& innermem,
0560 const std::string& outermem) {
0561 int outerrbits = 0;
0562 if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4 || iSeed == Seed::L1D1 || iSeed == Seed::L2D1) {
0563 outerrbits = 3;
0564 }
0565
0566 int outerrbins = (1 << outerrbits);
0567 int innerphibins = (1 << innerphibits);
0568 int outerphibins = (1 << outerphibits);
0569
0570 double phiinner[2];
0571 double phiouter[2];
0572 double router[2];
0573
0574 bool isPSinner;
0575 bool isPSouter;
0576
0577 if (iSeed == Seed::L3L4) {
0578 isPSinner = true;
0579 isPSouter = false;
0580 } else if (iSeed == Seed::L5L6) {
0581 isPSinner = false;
0582 isPSouter = false;
0583 } else {
0584 isPSinner = true;
0585 isPSouter = true;
0586 }
0587
0588 unsigned int nbendbitsinner = isPSinner ? N_BENDBITS_PS : N_BENDBITS_2S;
0589 unsigned int nbendbitsouter = isPSouter ? N_BENDBITS_PS : N_BENDBITS_2S;
0590
0591 if (fillTEMem) {
0592 if (fillInner) {
0593 table_.resize((1 << nbendbitsinner), false);
0594 } else {
0595 table_.resize((1 << nbendbitsouter), false);
0596 }
0597 }
0598
0599 double z0 = settings_.z0cut();
0600
0601 for (int irouterbin = 0; irouterbin < outerrbins; irouterbin++) {
0602 if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4 || iSeed == Seed::L1D1 || iSeed == Seed::L2D1) {
0603 router[0] = settings_.rmindiskvm() + irouterbin * (settings_.rmaxdiskvm() - settings_.rmindiskvm()) / outerrbins;
0604 router[1] =
0605 settings_.rmindiskvm() + (irouterbin + 1) * (settings_.rmaxdiskvm() - settings_.rmindiskvm()) / outerrbins;
0606 } else {
0607 router[0] = settings_.rmean(layerdisk2);
0608 router[1] = settings_.rmean(layerdisk2);
0609 }
0610
0611
0612 std::vector<std::array<double, 2>> bend_cuts_inner;
0613 std::vector<std::array<double, 2>> bend_cuts_outer;
0614
0615 if (settings_.useCalcBendCuts) {
0616 std::vector<const tt::SensorModule*> sminner;
0617 std::vector<const tt::SensorModule*> smouter;
0618
0619 if (iSeed == Seed::L1L2 || iSeed == Seed::L2L3 || iSeed == Seed::L3L4 || iSeed == Seed::L5L6) {
0620 double outer_tan_max = tan_theta(settings_.rmean(layerdisk2), settings_.zlength(), z0, true);
0621 std::array<double, 2> tan_range = {{0, outer_tan_max}};
0622
0623 smouter = getSensorModules(layerdisk2, isPSouter, tan_range);
0624 sminner = getSensorModules(layerdisk1, isPSinner, tan_range);
0625
0626 } else if (iSeed == Seed::L1D1 || iSeed == Seed::L2D1) {
0627 double outer_tan_min = tan_theta(router[1], settings_.zmindisk(layerdisk2 - N_LAYER), z0, false);
0628 double outer_tan_max = tan_theta(router[0], settings_.zmaxdisk(layerdisk2 - N_LAYER), z0, true);
0629
0630 smouter = getSensorModules(layerdisk2, isPSouter, {{outer_tan_min, outer_tan_max}});
0631 std::array<double, 2> tan_range = getTanRange(smouter);
0632 sminner = getSensorModules(layerdisk1, isPSinner, tan_range);
0633
0634 } else {
0635
0636 double outer_tan_min = tan_theta(router[1], settings_.zmindisk(layerdisk2 - N_LAYER), z0, false);
0637 double outer_tan_max = tan_theta(router[0], settings_.zmaxdisk(layerdisk2 - N_LAYER), z0, true);
0638
0639 smouter = getSensorModules(layerdisk2, isPSouter, {{outer_tan_min, outer_tan_max}});
0640
0641 std::array<double, 2> tan_range = getTanRange(smouter);
0642 sminner = getSensorModules(layerdisk1, isPSinner, tan_range);
0643 }
0644
0645 bend_cuts_inner = getBendCut(layerdisk1, sminner, isPSinner, settings_.bendcutTE(iSeed, true));
0646 bend_cuts_outer = getBendCut(layerdisk2, smouter, isPSouter, settings_.bendcutTE(iSeed, false));
0647
0648 } else {
0649 for (int ibend = 0; ibend < (1 << nbendbitsinner); ibend++) {
0650 double mid = settings_.benddecode(ibend, layerdisk1, nbendbitsinner == 3);
0651 double cut = settings_.bendcutte(ibend, layerdisk1, nbendbitsinner == 3);
0652 bend_cuts_inner.push_back({{mid, cut}});
0653 }
0654 for (int ibend = 0; ibend < (1 << nbendbitsouter); ibend++) {
0655 double mid = settings_.benddecode(ibend, layerdisk2, nbendbitsouter == 3);
0656 double cut = settings_.bendcutte(ibend, layerdisk2, nbendbitsouter == 3);
0657 bend_cuts_outer.push_back({{mid, cut}});
0658 }
0659 }
0660
0661 for (int iphiinnerbin = 0; iphiinnerbin < innerphibins; iphiinnerbin++) {
0662 phiinner[0] = innerphimin + iphiinnerbin * (innerphimax - innerphimin) / innerphibins;
0663 phiinner[1] = innerphimin + (iphiinnerbin + 1) * (innerphimax - innerphimin) / innerphibins;
0664 for (int iphiouterbin = 0; iphiouterbin < outerphibins; iphiouterbin++) {
0665 phiouter[0] = outerphimin + iphiouterbin * (outerphimax - outerphimin) / outerphibins;
0666 phiouter[1] = outerphimin + (iphiouterbin + 1) * (outerphimax - outerphimin) / outerphibins;
0667
0668 double bendinnermin = 20.0;
0669 double bendinnermax = -20.0;
0670 double bendoutermin = 20.0;
0671 double bendoutermax = -20.0;
0672 double rinvmin = 1.0;
0673 double rinvmax = -1.0;
0674 double absrinvmin = 1.0;
0675
0676 for (int i1 = 0; i1 < 2; i1++) {
0677 for (int i2 = 0; i2 < 2; i2++) {
0678 for (int i3 = 0; i3 < 2; i3++) {
0679 double rinner = 0.0;
0680 if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4) {
0681 rinner = router[i3] * settings_.zmean(layerdisk1 - N_LAYER) / settings_.zmean(layerdisk2 - N_LAYER);
0682 } else {
0683 rinner = settings_.rmean(layerdisk1);
0684 }
0685
0686 if (settings_.useCalcBendCuts) {
0687 if (rinner >= router[i3])
0688 continue;
0689 }
0690
0691 double rinv1 = (rinner < router[i3]) ? -rinv(phiinner[i1], phiouter[i2], rinner, router[i3]) : -20.0;
0692 double pitchinner =
0693 (rinner < settings_.rcrit()) ? settings_.stripPitch(true) : settings_.stripPitch(false);
0694 double pitchouter =
0695 (router[i3] < settings_.rcrit()) ? settings_.stripPitch(true) : settings_.stripPitch(false);
0696
0697 double abendinner = bendstrip(rinner, rinv1, pitchinner, settings_.sensorSpacing2S());
0698 double abendouter = bendstrip(router[i3], rinv1, pitchouter, settings_.sensorSpacing2S());
0699
0700 if (abendinner < bendinnermin)
0701 bendinnermin = abendinner;
0702 if (abendinner > bendinnermax)
0703 bendinnermax = abendinner;
0704 if (abendouter < bendoutermin)
0705 bendoutermin = abendouter;
0706 if (abendouter > bendoutermax)
0707 bendoutermax = abendouter;
0708 if (std::abs(rinv1) < absrinvmin)
0709 absrinvmin = std::abs(rinv1);
0710 if (rinv1 > rinvmax)
0711 rinvmax = rinv1;
0712 if (rinv1 < rinvmin)
0713 rinvmin = rinv1;
0714 }
0715 }
0716 }
0717
0718 double lowrinvcutte = 0.002;
0719
0720 bool passptcut;
0721 double bendfac;
0722
0723 if (settings_.useCalcBendCuts) {
0724 passptcut = rinvmin < settings_.rinvcutte() and rinvmax > -settings_.rinvcutte();
0725 bendfac = (rinvmin < lowrinvcutte and rinvmax > -lowrinvcutte) ? 1.05 : 1.0;
0726 } else {
0727 passptcut = absrinvmin < settings_.rinvcutte();
0728 bendfac = 1.0;
0729 }
0730
0731 if (fillInner) {
0732 for (int ibend = 0; ibend < (1 << nbendbitsinner); ibend++) {
0733 double bendminfac = (isPSinner and (ibend == 2 or ibend == 3)) ? bendfac : 1.0;
0734 double bendmaxfac = (isPSinner and (ibend == 6 or ibend == 5)) ? bendfac : 1.0;
0735
0736 double mid = bend_cuts_inner.at(ibend)[0];
0737 double cut = bend_cuts_inner.at(ibend)[1];
0738
0739 bool passinner = mid + cut * bendmaxfac > bendinnermin && mid - cut * bendminfac < bendinnermax;
0740
0741 if (fillTEMem) {
0742 if (passinner)
0743 table_[ibend] = 1;
0744 } else {
0745 table_.push_back(passinner && passptcut);
0746 }
0747 }
0748 } else {
0749 for (int ibend = 0; ibend < (1 << nbendbitsouter); ibend++) {
0750 double bendminfac = (isPSouter and (ibend == 2 or ibend == 3)) ? bendfac : 1.0;
0751 double bendmaxfac = (isPSouter and (ibend == 6 or ibend == 5)) ? bendfac : 1.0;
0752
0753 double mid = bend_cuts_outer.at(ibend)[0];
0754 double cut = bend_cuts_outer.at(ibend)[1];
0755
0756 bool passouter = mid + cut * bendmaxfac > bendoutermin && mid - cut * bendminfac < bendoutermax;
0757
0758 if (fillTEMem) {
0759 if (passouter)
0760 table_[ibend] = 1;
0761 } else {
0762 table_.push_back(passouter && passptcut);
0763 }
0764 }
0765 }
0766 }
0767 }
0768 }
0769
0770 positive_ = false;
0771
0772 if (fillTEMem) {
0773 if (fillInner) {
0774 name_ = "VMSTE_" + innermem + "_vmbendcut.tab";
0775 } else {
0776 name_ = "VMSTE_" + outermem + "_vmbendcut.tab";
0777 }
0778 } else {
0779 name_ = "TE_" + innermem.substr(0, innermem.size() - 2) + "_" + outermem.substr(0, outermem.size() - 2);
0780 if (fillInner) {
0781 name_ += "_stubptinnercut.tab";
0782 } else {
0783 name_ += "_stubptoutercut.tab";
0784 }
0785 }
0786 writeTable();
0787 }
0788
0789 void TrackletLUT::initProjectionBend(double k_phider,
0790 unsigned int idisk,
0791 unsigned int nrbits,
0792 unsigned int nphiderbits) {
0793 unsigned int nsignbins = 2;
0794 unsigned int nrbins = 1 << (nrbits);
0795 unsigned int nphiderbins = 1 << (nphiderbits);
0796
0797 for (unsigned int isignbin = 0; isignbin < nsignbins; isignbin++) {
0798 for (unsigned int irbin = 0; irbin < nrbins; irbin++) {
0799 int ir = irbin;
0800 if (ir > (1 << (nrbits - 1)))
0801 ir -= (1 << nrbits);
0802 ir = ir << (settings_.nrbitsstub(N_LAYER) - nrbits);
0803 for (unsigned int iphiderbin = 0; iphiderbin < nphiderbins; iphiderbin++) {
0804 int iphider = iphiderbin;
0805 if (iphider > (1 << (nphiderbits - 1)))
0806 iphider -= (1 << nphiderbits);
0807 iphider = iphider << (settings_.nbitsphiprojderL123() - nphiderbits);
0808
0809 double rproj = ir * settings_.krprojshiftdisk();
0810 double phider = iphider * k_phider;
0811 double t = settings_.zmean(idisk) / rproj;
0812
0813 if (isignbin)
0814 t = -t;
0815
0816 double rinv = -phider * (2.0 * t);
0817
0818 double stripPitch = (rproj < settings_.rcrit()) ? settings_.stripPitch(true) : settings_.stripPitch(false);
0819 double bendproj = bendstrip(rproj, rinv, stripPitch, settings_.sensorSpacing2S());
0820
0821 constexpr double maxbend = (1 << NRINVBITS) - 1;
0822
0823 int ibendproj = 2.0 * bendproj + 0.5 * maxbend;
0824 if (ibendproj < 0)
0825 ibendproj = 0;
0826 if (ibendproj > maxbend)
0827 ibendproj = maxbend;
0828
0829 table_.push_back(ibendproj);
0830 }
0831 }
0832 }
0833
0834 positive_ = false;
0835 name_ = settings_.combined() ? "MP_" : "PR_";
0836 name_ += "ProjectionBend_" + TrackletConfigBuilder::LayerName(N_LAYER + idisk) + ".tab";
0837
0838 writeTable();
0839 }
0840
0841 void TrackletLUT::initProjectionDiskRadius(int nrbits) {
0842
0843
0844
0845
0846 double roffset = 3.0;
0847
0848 for (unsigned int ir = 0; ir < (1u << nrbits); ir++) {
0849 double r = ir * settings_.rmaxdisk() / (1u << nrbits);
0850
0851 int rbin1 =
0852 (1 << N_RZBITS) * (r - roffset - settings_.rmindiskvm()) / (settings_.rmaxdisk() - settings_.rmindiskvm());
0853 int rbin2 =
0854 (1 << N_RZBITS) * (r + roffset - settings_.rmindiskvm()) / (settings_.rmaxdisk() - settings_.rmindiskvm());
0855
0856 if (rbin1 < 0) {
0857 rbin1 = 0;
0858 }
0859 rbin2 = clamp(rbin2, 0, ((1 << N_RZBITS) - 1));
0860
0861 assert(rbin1 <= rbin2);
0862 assert(rbin2 - rbin1 <= 1);
0863
0864 int d = rbin1 != rbin2;
0865
0866 int finer =
0867 (1 << (N_RZBITS + NFINERZBITS)) *
0868 ((r - settings_.rmindiskvm()) - rbin1 * (settings_.rmaxdisk() - settings_.rmindiskvm()) / (1 << N_RZBITS)) /
0869 (settings_.rmaxdisk() - settings_.rmindiskvm());
0870
0871 finer = clamp(finer, 0, ((1 << (NFINERZBITS + 1)) - 1));
0872
0873
0874 int N_DIFF_FLAG = 1;
0875
0876 int word = (finer << (N_RZBITS + N_DIFF_FLAG)) + (rbin1 << N_DIFF_FLAG) + d;
0877
0878 table_.push_back(word);
0879 }
0880
0881
0882 nbits_ = NFINERZBITS + 1 + N_RZBITS + 1;
0883 positive_ = true;
0884 name_ = "ProjectionDiskRadius.tab";
0885 writeTable();
0886 }
0887
0888 void TrackletLUT::initBendMatch(unsigned int layerdisk) {
0889 unsigned int nrinv = NRINVBITS;
0890 double rinvhalf = 0.5 * ((1 << nrinv) - 1);
0891
0892 bool barrel = layerdisk < N_LAYER;
0893
0894 if (barrel) {
0895 bool isPSmodule = layerdisk < N_PSLAYER;
0896 double stripPitch = settings_.stripPitch(isPSmodule);
0897 unsigned int nbits = isPSmodule ? N_BENDBITS_PS : N_BENDBITS_2S;
0898
0899 std::vector<std::array<double, 2>> bend_cuts;
0900
0901 if (settings_.useCalcBendCuts) {
0902 double bendcutFE = settings_.bendcutME(layerdisk, isPSmodule);
0903 std::vector<const tt::SensorModule*> sm = getSensorModules(layerdisk, isPSmodule);
0904 bend_cuts = getBendCut(layerdisk, sm, isPSmodule, bendcutFE);
0905
0906 } else {
0907 for (unsigned int ibend = 0; ibend < (1u << nbits); ibend++) {
0908 double mid = settings_.benddecode(ibend, layerdisk, isPSmodule);
0909 double cut = settings_.bendcutte(ibend, layerdisk, isPSmodule);
0910 bend_cuts.push_back({{mid, cut}});
0911 }
0912 }
0913
0914 for (unsigned int irinv = 0; irinv < (1u << nrinv); irinv++) {
0915 double rinv = (irinv - rinvhalf) * (1 << (settings_.nbitsrinv() - nrinv)) * settings_.krinvpars();
0916
0917 double projbend = bendstrip(settings_.rmean(layerdisk), rinv, stripPitch, settings_.sensorSpacing2S());
0918 for (unsigned int ibend = 0; ibend < (1u << nbits); ibend++) {
0919 double mid = bend_cuts[ibend][0];
0920 double cut = bend_cuts[ibend][1];
0921
0922 double pass = mid + cut > projbend && mid - cut < projbend;
0923
0924 table_.push_back(pass);
0925 }
0926 }
0927 } else {
0928 std::vector<std::array<double, 2>> bend_cuts_2S;
0929 std::vector<std::array<double, 2>> bend_cuts_PS;
0930
0931 if (settings_.useCalcBendCuts) {
0932 double bendcutFE2S = settings_.bendcutME(layerdisk, false);
0933 std::vector<const tt::SensorModule*> sm2S = getSensorModules(layerdisk, false);
0934 bend_cuts_2S = getBendCut(layerdisk, sm2S, false, bendcutFE2S);
0935
0936 double bendcutFEPS = settings_.bendcutME(layerdisk, true);
0937 std::vector<const tt::SensorModule*> smPS = getSensorModules(layerdisk, true);
0938 bend_cuts_PS = getBendCut(layerdisk, smPS, true, bendcutFEPS);
0939
0940 } else {
0941 for (unsigned int ibend = 0; ibend < (1 << N_BENDBITS_2S); ibend++) {
0942 double mid = settings_.benddecode(ibend, layerdisk, false);
0943 double cut = settings_.bendcutme(ibend, layerdisk, false);
0944 bend_cuts_2S.push_back({{mid, cut}});
0945 }
0946 for (unsigned int ibend = 0; ibend < (1 << N_BENDBITS_PS); ibend++) {
0947 double mid = settings_.benddecode(ibend, layerdisk, true);
0948 double cut = settings_.bendcutme(ibend, layerdisk, true);
0949 bend_cuts_PS.push_back({{mid, cut}});
0950 }
0951 }
0952
0953 for (unsigned int iprojbend = 0; iprojbend < (1u << nrinv); iprojbend++) {
0954 double projbend = 0.5 * (iprojbend - rinvhalf);
0955 for (unsigned int ibend = 0; ibend < (1 << N_BENDBITS_2S); ibend++) {
0956 double mid = bend_cuts_2S[ibend][0];
0957 double cut = bend_cuts_2S[ibend][1];
0958
0959 double pass = mid + cut > projbend && mid - cut < projbend;
0960
0961 table_.push_back(pass);
0962 }
0963 }
0964 for (unsigned int iprojbend = 0; iprojbend < (1u << nrinv); iprojbend++) {
0965 double projbend = 0.5 * (iprojbend - rinvhalf);
0966 for (unsigned int ibend = 0; ibend < (1 << N_BENDBITS_PS); ibend++) {
0967 double mid = bend_cuts_PS[ibend][0];
0968 double cut = bend_cuts_PS[ibend][1];
0969
0970 double pass = mid + cut > projbend && mid - cut < projbend;
0971
0972 table_.push_back(pass);
0973 }
0974 }
0975 }
0976
0977 positive_ = false;
0978
0979 name_ = "METable_" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
0980
0981 writeTable();
0982 }
0983
0984 void TrackletLUT::initVMRTable(unsigned int layerdisk, VMRTableType type, int region) {
0985 unsigned int zbits = settings_.vmrlutzbits(layerdisk);
0986 unsigned int rbits = settings_.vmrlutrbits(layerdisk);
0987
0988 unsigned int rbins = (1 << rbits);
0989 unsigned int zbins = (1 << zbits);
0990
0991 double zmin, zmax, rmin, rmax;
0992
0993 if (layerdisk < N_LAYER) {
0994 zmin = -settings_.zlength();
0995 zmax = settings_.zlength();
0996 rmin = settings_.rmean(layerdisk) - settings_.drmax();
0997 rmax = settings_.rmean(layerdisk) + settings_.drmax();
0998 } else {
0999 rmin = 0;
1000 rmax = settings_.rmaxdisk();
1001 zmin = settings_.zmean(layerdisk - N_LAYER) - settings_.dzmax();
1002 zmax = settings_.zmean(layerdisk - N_LAYER) + settings_.dzmax();
1003 }
1004
1005 double dr = (rmax - rmin) / rbins;
1006 double dz = (zmax - zmin) / zbins;
1007
1008 int NBINS = settings_.NLONGVMBINS() * settings_.NLONGVMBINS();
1009
1010 for (unsigned int izbin = 0; izbin < zbins; izbin++) {
1011 for (unsigned int irbin = 0; irbin < rbins; irbin++) {
1012 double r = rmin + (irbin + 0.5) * dr;
1013 double z = zmin + (izbin + 0.5) * dz;
1014
1015 if (settings_.combined()) {
1016 int iznew = izbin - (1 << (zbits - 1));
1017 if (iznew < 0)
1018 iznew += (1 << zbits);
1019 assert(iznew >= 0);
1020 assert(iznew < (1 << zbits));
1021 z = zmin + (iznew + 0.5) * dz;
1022 if (layerdisk < N_LAYER) {
1023 int irnew = irbin - (1 << (rbits - 1));
1024 if (irnew < 0)
1025 irnew += (1 << rbits);
1026 assert(irnew >= 0);
1027 assert(irnew < (1 << rbits));
1028 r = rmin + (irnew + 0.5) * dr;
1029 }
1030 }
1031
1032 unsigned int NRING =
1033 5;
1034 if (layerdisk >= N_LAYER && irbin < 2 * NRING)
1035 r = (layerdisk < N_LAYER + 2) ? settings_.rDSSinner(irbin) : settings_.rDSSouter(irbin);
1036
1037 int bin;
1038 if (layerdisk < N_LAYER) {
1039 double zproj = z * settings_.rmean(layerdisk) / r;
1040 bin = NBINS * (zproj + settings_.zlength()) / (2 * settings_.zlength());
1041 } else {
1042 double rproj = r * settings_.zmean(layerdisk - N_LAYER) / z;
1043 bin = NBINS * (rproj - settings_.rmindiskvm()) / (settings_.rmaxdisk() - settings_.rmindiskvm());
1044 }
1045 if (bin < 0)
1046 bin = 0;
1047 if (bin >= NBINS)
1048 bin = NBINS - 1;
1049
1050 if (type == VMRTableType::me) {
1051 table_.push_back(bin);
1052 }
1053
1054 if (type == VMRTableType::disk) {
1055 if (layerdisk >= N_LAYER) {
1056 double rproj = r * settings_.zmean(layerdisk - N_LAYER) / z;
1057 bin = 0.5 * NBINS * (rproj - settings_.rmindiskvm()) / (settings_.rmaxdiskvm() - settings_.rmindiskvm());
1058
1059 if (bin < 0)
1060 bin = 0;
1061 if (bin >= NBINS / 2)
1062 bin = 0;
1063 table_.push_back(bin);
1064 }
1065 }
1066
1067 if (type == VMRTableType::inner) {
1068 if (layerdisk == LayerDisk::L1 || layerdisk == LayerDisk::L3 || layerdisk == LayerDisk::L5 ||
1069 layerdisk == LayerDisk::D1 || layerdisk == LayerDisk::D3) {
1070 table_.push_back(getVMRLookup(layerdisk + 1, z, r, dz, dr));
1071 }
1072 if (layerdisk == LayerDisk::L2) {
1073 table_.push_back(getVMRLookup(layerdisk + 1, z, r, dz, dr, Seed::L2L3));
1074 }
1075 }
1076
1077 if (type == VMRTableType::inneroverlap) {
1078 if (layerdisk == LayerDisk::L1 || layerdisk == LayerDisk::L2) {
1079 table_.push_back(getVMRLookup(6, z, r, dz, dr, layerdisk + 6));
1080 }
1081 }
1082
1083 if (type == VMRTableType::innerthird) {
1084 if (layerdisk == LayerDisk::L2) {
1085 table_.push_back(getVMRLookup(LayerDisk::D1, z, r, dz, dr, Seed::L2L3D1));
1086 }
1087
1088 if (layerdisk == LayerDisk::L5) {
1089 table_.push_back(getVMRLookup(LayerDisk::L4, z, r, dz, dr));
1090 }
1091
1092 if (layerdisk == LayerDisk::L3) {
1093 table_.push_back(getVMRLookup(LayerDisk::L2, z, r, dz, dr));
1094 }
1095
1096 if (layerdisk == LayerDisk::D1) {
1097 table_.push_back(getVMRLookup(LayerDisk::L2, z, r, dz, dr));
1098 }
1099 }
1100 }
1101 }
1102
1103 if (settings_.combined()) {
1104 if (type == VMRTableType::me) {
1105 nbits_ = 2 * settings_.NLONGVMBITS();
1106 positive_ = false;
1107 name_ = "VMRME_" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
1108 }
1109 if (type == VMRTableType::disk) {
1110 nbits_ = 2 * settings_.NLONGVMBITS();
1111 positive_ = false;
1112 name_ = "VMRTE_" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
1113 }
1114 if (type == VMRTableType::inner) {
1115 positive_ = true;
1116 nbits_ = 10;
1117 name_ = "TP_" + TrackletConfigBuilder::LayerName(layerdisk) + TrackletConfigBuilder::LayerName(layerdisk + 1) +
1118 ".tab";
1119 }
1120
1121 if (type == VMRTableType::inneroverlap) {
1122 positive_ = true;
1123 nbits_ = 10;
1124 name_ = "TP_" + TrackletConfigBuilder::LayerName(layerdisk) + TrackletConfigBuilder::LayerName(N_LAYER) + ".tab";
1125 }
1126
1127 } else {
1128 if (type == VMRTableType::me) {
1129
1130 if (layerdisk == LayerDisk::L2 || layerdisk == LayerDisk::L3 || layerdisk == LayerDisk::L4 ||
1131 layerdisk == LayerDisk::L6) {
1132 positive_ = false;
1133 name_ = "VMTableOuter" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
1134 writeTable();
1135 }
1136
1137 assert(region >= 0);
1138 char cregion = 'A' + region;
1139 name_ = "VMR_" + TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_finebin.tab";
1140 positive_ = false;
1141 }
1142
1143 if (type == VMRTableType::inner) {
1144 positive_ = false;
1145 name_ = "VMTableInner" + TrackletConfigBuilder::LayerName(layerdisk) +
1146 TrackletConfigBuilder::LayerName(layerdisk + 1) + ".tab";
1147 }
1148
1149 if (type == VMRTableType::inneroverlap) {
1150 positive_ = false;
1151 name_ = "VMTableInner" + TrackletConfigBuilder::LayerName(layerdisk) + TrackletConfigBuilder::LayerName(N_LAYER) +
1152 ".tab";
1153 }
1154
1155 if (type == VMRTableType::disk) {
1156 positive_ = false;
1157 name_ = "VMTableOuter" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
1158 }
1159 }
1160
1161 writeTable();
1162 }
1163
1164 int TrackletLUT::getVMRLookup(unsigned int layerdisk, double z, double r, double dz, double dr, int iseed) const {
1165 double z0cut = settings_.z0cut();
1166
1167 if (layerdisk < N_LAYER) {
1168 double constexpr zcutL2L3 = 52.0;
1169 if (iseed == Seed::L2L3 && std::abs(z) < zcutL2L3)
1170 return -1;
1171
1172 double rmean = settings_.rmean(layerdisk);
1173
1174 double rratio1 = rmean / (r + 0.5 * dr);
1175 double rratio2 = rmean / (r - 0.5 * dr);
1176
1177 double z1 = (z - 0.5 * dz) * rratio1 + z0cut * (rratio1 - 1.0);
1178 double z2 = (z + 0.5 * dz) * rratio1 + z0cut * (rratio1 - 1.0);
1179 double z3 = (z - 0.5 * dz) * rratio2 + z0cut * (rratio2 - 1.0);
1180 double z4 = (z + 0.5 * dz) * rratio2 + z0cut * (rratio2 - 1.0);
1181 double z5 = (z - 0.5 * dz) * rratio1 - z0cut * (rratio1 - 1.0);
1182 double z6 = (z + 0.5 * dz) * rratio1 - z0cut * (rratio1 - 1.0);
1183 double z7 = (z - 0.5 * dz) * rratio2 - z0cut * (rratio2 - 1.0);
1184 double z8 = (z + 0.5 * dz) * rratio2 - z0cut * (rratio2 - 1.0);
1185
1186 double zmin = std::min({z1, z2, z3, z4, z5, z6, z7, z8});
1187 double zmax = std::max({z1, z2, z3, z4, z5, z6, z7, z8});
1188
1189 int NBINS = settings_.NLONGVMBINS() * settings_.NLONGVMBINS();
1190
1191 int zbin1 = NBINS * (zmin + settings_.zlength()) / (2 * settings_.zlength());
1192 int zbin2 = NBINS * (zmax + settings_.zlength()) / (2 * settings_.zlength());
1193
1194 if (zbin1 >= NBINS)
1195 return -1;
1196 if (zbin2 < 0)
1197 return -1;
1198
1199 if (zbin2 >= NBINS)
1200 zbin2 = NBINS - 1;
1201 if (zbin1 < 0)
1202 zbin1 = 0;
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214 int value = zbin1 / 8;
1215 value *= 2;
1216 if (zbin2 / 8 - zbin1 / 8 > 0)
1217 value += 1;
1218 value *= 8;
1219 value += (zbin1 & 7);
1220 assert(value / 8 < 15);
1221 int deltaz = zbin2 - zbin1;
1222 if (deltaz > 7) {
1223 deltaz = 7;
1224 }
1225 assert(deltaz < 8);
1226 value += (deltaz << 7);
1227
1228 return value;
1229
1230 } else {
1231 if (std::abs(z) < 2.0 * z0cut)
1232 return -1;
1233
1234 double zmean = settings_.zmean(layerdisk - N_LAYER);
1235 if (z < 0.0)
1236 zmean = -zmean;
1237
1238 double r1 = (r + 0.5 * dr) * (zmean + z0cut) / (z + 0.5 * dz + z0cut);
1239 double r2 = (r - 0.5 * dr) * (zmean - z0cut) / (z + 0.5 * dz - z0cut);
1240 double r3 = (r + 0.5 * dr) * (zmean + z0cut) / (z - 0.5 * dz + z0cut);
1241 double r4 = (r - 0.5 * dr) * (zmean - z0cut) / (z - 0.5 * dz - z0cut);
1242 double r5 = (r + 0.5 * dr) * (zmean - z0cut) / (z + 0.5 * dz - z0cut);
1243 double r6 = (r - 0.5 * dr) * (zmean + z0cut) / (z + 0.5 * dz + z0cut);
1244 double r7 = (r + 0.5 * dr) * (zmean - z0cut) / (z - 0.5 * dz - z0cut);
1245 double r8 = (r - 0.5 * dr) * (zmean + z0cut) / (z - 0.5 * dz + z0cut);
1246
1247 double rmin = std::min({r1, r2, r3, r4, r5, r6, r7, r8});
1248 double rmax = std::max({r1, r2, r3, r4, r5, r6, r7, r8});
1249
1250 int NBINS = settings_.NLONGVMBINS() * settings_.NLONGVMBINS() / 2;
1251
1252 double rmindisk = settings_.rmindiskvm();
1253 double rmaxdisk = settings_.rmaxdiskvm();
1254
1255 if (iseed == Seed::L1D1)
1256 rmaxdisk = settings_.rmaxdiskl1overlapvm();
1257 if (iseed == Seed::L2D1)
1258 rmindisk = settings_.rmindiskl2overlapvm();
1259 if (iseed == Seed::L2L3D1)
1260 rmaxdisk = settings_.rmaxdisk();
1261
1262 if (rmin > rmaxdisk)
1263 return -1;
1264 if (rmax > rmaxdisk)
1265 rmax = rmaxdisk;
1266
1267 if (rmax < rmindisk)
1268 return -1;
1269 if (rmin < rmindisk)
1270 rmin = rmindisk;
1271
1272 int rbin1 = NBINS * (rmin - settings_.rmindiskvm()) / (settings_.rmaxdiskvm() - settings_.rmindiskvm());
1273 int rbin2 = NBINS * (rmax - settings_.rmindiskvm()) / (settings_.rmaxdiskvm() - settings_.rmindiskvm());
1274
1275 if (iseed == Seed::L2L3D1) {
1276 constexpr double rminspec = 40.0;
1277 rbin1 = NBINS * (rmin - rminspec) / (settings_.rmaxdisk() - rminspec);
1278 rbin2 = NBINS * (rmax - rminspec) / (settings_.rmaxdisk() - rminspec);
1279 }
1280
1281 if (rbin2 >= NBINS)
1282 rbin2 = NBINS - 1;
1283 if (rbin1 < 0)
1284 rbin1 = 0;
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296 bool overlap = iseed == Seed::L1D1 || iseed == Seed::L2D1 || iseed == Seed::L2L3D1;
1297
1298 int value = rbin1 / 8;
1299 if (overlap) {
1300 if (z < 0.0)
1301 value += 4;
1302 }
1303 value *= 2;
1304 if (rbin2 / 8 - rbin1 / 8 > 0)
1305 value += 1;
1306 value *= 8;
1307 value += (rbin1 & 7);
1308 assert(value / 8 < 15);
1309 int deltar = rbin2 - rbin1;
1310 if (deltar > 7)
1311 deltar = 7;
1312 if (overlap) {
1313 value += (deltar << 7);
1314 } else {
1315 value += (deltar << 6);
1316 }
1317
1318 return value;
1319 }
1320 }
1321
1322 void TrackletLUT::initPhiCorrTable(unsigned int layerdisk, unsigned int rbits) {
1323 bool psmodule = layerdisk < N_PSLAYER;
1324
1325 unsigned int bendbits = psmodule ? N_BENDBITS_PS : N_BENDBITS_2S;
1326
1327 unsigned int rbins = (1 << rbits);
1328
1329 double rmean = settings_.rmean(layerdisk);
1330 double drmax = settings_.drmax();
1331
1332 double dr = 2.0 * drmax / rbins;
1333
1334 std::vector<std::array<double, 2>> bend_vals;
1335
1336 if (settings_.useCalcBendCuts) {
1337 std::vector<const tt::SensorModule*> sm = getSensorModules(layerdisk, psmodule);
1338 bend_vals = getBendCut(layerdisk, sm, psmodule);
1339
1340 } else {
1341 for (int ibend = 0; ibend < 1 << bendbits; ibend++) {
1342 bend_vals.push_back({{settings_.benddecode(ibend, layerdisk, layerdisk < N_PSLAYER), 0}});
1343 }
1344 }
1345
1346 for (int ibend = 0; ibend < 1 << bendbits; ibend++) {
1347 for (unsigned int irbin = 0; irbin < rbins; irbin++) {
1348 double bend = -bend_vals[ibend][0];
1349 int value = getphiCorrValue(layerdisk, bend, irbin, rmean, dr, drmax);
1350 table_.push_back(value);
1351 }
1352 }
1353
1354 name_ = "VMPhiCorrL" + std::to_string(layerdisk + 1) + ".tab";
1355 nbits_ = 14;
1356 positive_ = false;
1357
1358 writeTable();
1359 }
1360
1361 int TrackletLUT::getphiCorrValue(
1362 unsigned int layerdisk, double bend, unsigned int irbin, double rmean, double dr, double drmax) const {
1363 bool psmodule = layerdisk < N_PSLAYER;
1364
1365
1366 double Delta = (irbin + 0.5) * dr - drmax;
1367
1368
1369 double drnom = 0.18;
1370 double dphi = (Delta / drnom) * bend * settings_.stripPitch(psmodule) / rmean;
1371
1372 double kphi = psmodule ? settings_.kphi() : settings_.kphi1();
1373
1374 int idphi = dphi / kphi;
1375
1376 return idphi;
1377 }
1378
1379
1380 void TrackletLUT::writeTable() const {
1381 if (!settings_.writeTable()) {
1382 return;
1383 }
1384
1385 if (name_.empty()) {
1386 return;
1387 }
1388
1389 ofstream out = openfile(settings_.tablePath(), name_, __FILE__, __LINE__);
1390
1391 out << "{" << endl;
1392 for (unsigned int i = 0; i < table_.size(); i++) {
1393 if (i != 0) {
1394 out << "," << endl;
1395 }
1396
1397 int itable = table_[i];
1398 if (positive_) {
1399 if (table_[i] < 0) {
1400 itable = (1 << nbits_) - 1;
1401 }
1402 }
1403
1404 out << itable;
1405 }
1406 out << endl << "};" << endl;
1407 out.close();
1408
1409 string name = name_;
1410
1411 name[name_.size() - 3] = 'd';
1412 name[name_.size() - 2] = 'a';
1413 name[name_.size() - 1] = 't';
1414
1415 out = openfile(settings_.tablePath(), name, __FILE__, __LINE__);
1416
1417 int width = (nbits_ + 3) / 4;
1418
1419 for (unsigned int i = 0; i < table_.size(); i++) {
1420 int itable = table_[i];
1421 if (positive_) {
1422 if (table_[i] < 0) {
1423 itable = (1 << nbits_) - 1;
1424 }
1425 }
1426
1427 out << uppercase << setfill('0') << setw(width) << hex << itable << dec << endl;
1428 }
1429
1430 out.close();
1431 }
1432
1433 int TrackletLUT::lookup(unsigned int index) const {
1434 assert(index < table_.size());
1435 return table_[index];
1436 }