Back to home page

Project CMSSW displayed by LXR
 
 

    

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   //Returns a vector of SensorModules using T. Schuh's Setup and SensorModule classes.
0017   //Can be used 3 ways:
0018   //Default: No specified tan_range or nzbins, returns all SensorModules in specified layerdisk (unique in |z|)
0019   //tan_range: Returns modules in given tan range, where the min and max tan(theta) are measured from 0 -/+ z0 to account for displaced tracks
0020   //zbins: Returns modules in specified z bin (2 zbins = (Flat, Tilted), 13 zbins = (Flat, TR1, ..., TR12). Only for tilted barrel
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       //z0_max is swapped here so that the comparison down 5 lines is from same origin (+/- z0)
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   //Remove Duplicate Modules
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   //Given a set of modules returns a range in tan(theta), the angle is measured in the r-z(+/-z0) plane from the r-axis
0091 
0092   std::array<double, 2> tan_range = {{2147483647, 0}};  //(tan_min, tan_max)
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);  //(r, z, z0, bool z0_max), z0_max measures from +/- z0
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   //Finds range of bendstrip for given SensorModules as a function of the encoded bend. Returns in format (mid, half_range).
0120   //This uses the stub windows provided by T. Schuh's SensorModule class to determine the bend encoding. TODO test changes in stub windows
0121   //Any other change to the bend encoding requires changes here, perhaps a function that given (FEbend, isPS, stub window) and outputs an encoded bend
0122   //would be useful for consistency.
0123 
0124   unsigned int bendbits = isPS ? 3 : 4;
0125 
0126   std::vector<std::array<double, 2>> bendpars;    // mid, cut
0127   std::vector<std::array<double, 2>> bendminmax;  // min, max
0128 
0129   //Initialize array
0130   for (int i = 0; i < 1 << bendbits; i++) {
0131     bendpars.push_back({{99, 0}});
0132     bendminmax.push_back({{99, -99}});
0133   }
0134 
0135   //Loop over modules
0136   for (auto sm : sensorModules) {
0137     int window = sm->windowSize();  //Half-strip units
0138     const vector<double>& encodingBend = setup_->encodingBend(window, isPS);
0139 
0140     //Loop over FEbends
0141     for (int ibend = 0; ibend <= 2 * window; ibend++) {
0142       int FEbend = ibend - window;                                                 //Half-strip units
0143       double BEbend = setup_->stubAlgorithm()->degradeBend(isPS, window, FEbend);  //Full strip units
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);  //Encoded bend
0148 
0149       double bendmin = FEbend / 2.0 - FEbendcut;  //Full Strip units
0150       double bendmax = FEbend / 2.0 + FEbendcut;
0151 
0152       //Convert to bendstrip, calculate at module edges (z min, r max) and  (z max, r min)
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();  //z max
0157       z_mod[1] = std::abs(sm->z()) - (sm->numColumns() / 2 - 0.5) * sm->pitchCol() * sm->cosTilt();  //z min
0158 
0159       r_mod[0] = sm->r() - (sm->numColumns() / 2 - 0.5) * sm->pitchCol() * std::abs(sm->sinTilt());  //r min
0160       r_mod[1] = sm->r() + (sm->numColumns() / 2 - 0.5) * sm->pitchCol() * std::abs(sm->sinTilt());  //r max
0161 
0162       for (int i = 0; i < 2; i++) {  // 2 points to cover range in tan(theta) = z/r
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   //Convert min, max to mid, cut for ease of use
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   //number of fine phi bins in sector
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     //min and max dphi
0323     //ramge of dphi to consider due to resolution
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       //Determine bend cuts using geometry
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 {  // D1D2 D3D4
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;  //Somewhat arbitrary value, allows for better acceptance in bins with low rinv (high pt)
0438         passptcut = rinvmin < rinvcutte and rinvmax > -rinvcutte;
0439         bendfac = (rinvmin < lowrinvcutte and rinvmax > -lowrinvcutte)
0440                       ? 1.05
0441                       : 1.0;  //Somewhat arbirary value, bend cuts are 5% larger in bins with low rinv (high pt)
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     //Determine bend cuts using geometry
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 {  // D1D2 D3D4
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;  // Better acceptance for high pt
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   //When a projection to a disk is considered this offset and added and subtracted to calculate
0843   //the bin the projection is pointing to. This is to account for resolution effects such that
0844   //projections that are near a bin boundary will be assigned to both bins. The value (3 cm) should
0845   //cover the uncertanty in the resolution.
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     //Pack the data in a 8 bit word (ffffrrrd) where f is finer, r is rbin1, and d is difference
0874     int N_DIFF_FLAG = 1;  // Single bit for bool flag
0875 
0876     int word = (finer << (N_RZBITS + N_DIFF_FLAG)) + (rbin1 << N_DIFF_FLAG) + d;
0877 
0878     table_.push_back(word);
0879   }
0880 
0881   //Size of the data word from above (8 bits)
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++) {  //Should this be binned in r?
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;  //number of 2S rings in disks. This is multiplied below by two since we have two halfs of a module
1034       if (layerdisk >= N_LAYER && irbin < 2 * NRING)  //special case for the tabulated radii in 2S disks
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           //bin value of zero indicates that stub is out of range
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) {  //projection from L2 to D1 for L2L3D1 seeding
1085           table_.push_back(getVMRLookup(LayerDisk::D1, z, r, dz, dr, Seed::L2L3D1));
1086         }
1087 
1088         if (layerdisk == LayerDisk::L5) {  //projection from L5 to L4 for L5L6L4 seeding
1089           table_.push_back(getVMRLookup(LayerDisk::L4, z, r, dz, dr));
1090         }
1091 
1092         if (layerdisk == LayerDisk::L3) {  //projection from L3 to L5 for L3L4L2 seeding
1093           table_.push_back(getVMRLookup(LayerDisk::L2, z, r, dz, dr));
1094         }
1095 
1096         if (layerdisk == LayerDisk::D1) {  //projection from D1 to L2 for D1D2L2 seeding
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       //This if a hack where the same memory is used in both ME and TE modules
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;  //Stubs closer to IP in z will not be used for L2L3 seeds
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     // This is a 10 bit word:
1205     // xxx|yyy|z|rrr
1206     // xxx is the delta z window
1207     // yyy is the z bin
1208     // z is flag to look in next bin
1209     // rrr first fine z bin
1210     // NOTE : this encoding is not efficient z is one if xxx+rrr is greater than 8
1211     //        and xxx is only 1,2, or 3
1212     //        should also reject xxx=0 as this means projection is outside range
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     // This is a 9 bit word:
1287     // xxx|yy|z|rrr
1288     // xxx is the delta r window
1289     // yy is the r bin yy is three bits for overlaps
1290     // z is flag to look in next bin
1291     // rrr fine r bin
1292     // NOTE : this encoding is not efficient z is one if xxx+rrr is greater than 8
1293     //        and xxx is only 1,2, or 3
1294     //        should also reject xxx=0 as this means projection is outside range
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   //for the rbin - calculate the distance to the nominal layer radius
1366   double Delta = (irbin + 0.5) * dr - drmax;
1367 
1368   //calculate the phi correction - this is a somewhat approximate formula
1369   double drnom = 0.18;  //This is the nominal module separation for which bend is referenced
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 // Write LUT table.
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!