Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-23 16:00:17 

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