Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:22:05 

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) {
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       if (settings_.combined()) {
1030         int iznew = izbin - (1 << (zbits - 1));
1031         if (iznew < 0)
1032           iznew += (1 << zbits);
1033         assert(iznew >= 0);
1034         assert(iznew < (1 << zbits));
1035         z = zmin + (iznew + 0.5) * dz;
1036         if (layerdisk < N_LAYER) {
1037           int irnew = irbin - (1 << (rbits - 1));
1038           if (irnew < 0)
1039             irnew += (1 << rbits);
1040           assert(irnew >= 0);
1041           assert(irnew < (1 << rbits));
1042           r = rmin + (irnew + 0.5) * dr;
1043         }
1044       }
1045 
1046       unsigned int NRING =
1047           5;  //number of 2S rings in disks. This is multiplied below by two since we have two halfs of a module
1048       if (layerdisk >= N_LAYER && irbin < 2 * NRING)  //special case for the tabulated radii in 2S disks
1049         r = (layerdisk < N_LAYER + 2) ? settings_.rDSSinner(irbin) : settings_.rDSSouter(irbin);
1050 
1051       int bin;
1052       if (layerdisk < N_LAYER) {
1053         double zproj = z * settings_.rmean(layerdisk) / r;
1054         bin = NBINS * (zproj + settings_.zlength()) / (2 * settings_.zlength());
1055       } else {
1056         double rproj = r * settings_.zmean(layerdisk - N_LAYER) / z;
1057         bin = NBINS * (rproj - settings_.rmindiskvm()) / (settings_.rmaxdisk() - settings_.rmindiskvm());
1058       }
1059       if (bin < 0)
1060         bin = 0;
1061       if (bin >= NBINS)
1062         bin = NBINS - 1;
1063 
1064       if (type == VMRTableType::me) {
1065         table_.push_back(bin);
1066       }
1067 
1068       if (type == VMRTableType::disk) {
1069         if (layerdisk >= N_LAYER) {
1070           double rproj = r * settings_.zmean(layerdisk - N_LAYER) / z;
1071           bin = 0.5 * NBINS * (rproj - settings_.rmindiskvm()) / (settings_.rmaxdiskvm() - settings_.rmindiskvm());
1072           //bin value of zero indicates that stub is out of range
1073           if (bin < 0)
1074             bin = 0;
1075           if (bin >= NBINS / 2)
1076             bin = 0;
1077           table_.push_back(bin);
1078         }
1079       }
1080 
1081       if (type == VMRTableType::inner) {
1082         if (layerdisk == LayerDisk::L1 || layerdisk == LayerDisk::L3 || layerdisk == LayerDisk::L5 ||
1083             layerdisk == LayerDisk::D1 || layerdisk == LayerDisk::D3) {
1084           table_.push_back(getVMRLookup(layerdisk + 1, z, r, dz, dr));
1085         }
1086         if (layerdisk == LayerDisk::L2) {
1087           table_.push_back(getVMRLookup(layerdisk + 1, z, r, dz, dr, Seed::L2L3));
1088         }
1089       }
1090 
1091       if (type == VMRTableType::inneroverlap) {
1092         if (layerdisk == LayerDisk::L1 || layerdisk == LayerDisk::L2) {
1093           table_.push_back(getVMRLookup(6, z, r, dz, dr, layerdisk + 6));
1094         }
1095       }
1096 
1097       if (type == VMRTableType::innerthird) {
1098         if (layerdisk == LayerDisk::L2) {  //projection from L2 to D1 for L2L3D1 seeding
1099           table_.push_back(getVMRLookup(LayerDisk::D1, z, r, dz, dr, Seed::L2L3D1));
1100         }
1101 
1102         if (layerdisk == LayerDisk::L5) {  //projection from L5 to L4 for L5L6L4 seeding
1103           table_.push_back(getVMRLookup(LayerDisk::L4, z, r, dz, dr));
1104         }
1105 
1106         if (layerdisk == LayerDisk::L3) {  //projection from L3 to L5 for L3L4L2 seeding
1107           table_.push_back(getVMRLookup(LayerDisk::L2, z, r, dz, dr));
1108         }
1109 
1110         if (layerdisk == LayerDisk::D1) {  //projection from D1 to L2 for D1D2L2 seeding
1111           table_.push_back(getVMRLookup(LayerDisk::L2, z, r, dz, dr));
1112         }
1113       }
1114     }
1115   }
1116 
1117   if (settings_.combined()) {
1118     if (type == VMRTableType::me) {
1119       nbits_ = 2 * settings_.NLONGVMBITS();
1120       positive_ = false;
1121       name_ = "VMRME_" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
1122     }
1123     if (type == VMRTableType::disk) {
1124       nbits_ = 2 * settings_.NLONGVMBITS();
1125       positive_ = false;
1126       name_ = "VMRTE_" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
1127     }
1128     if (type == VMRTableType::inner) {
1129       positive_ = true;
1130       nbits_ = 10;
1131       name_ = "TP_" + TrackletConfigBuilder::LayerName(layerdisk) + TrackletConfigBuilder::LayerName(layerdisk + 1) +
1132               ".tab";
1133     }
1134 
1135     if (type == VMRTableType::inneroverlap) {
1136       positive_ = true;
1137       nbits_ = 10;
1138       name_ = "TP_" + TrackletConfigBuilder::LayerName(layerdisk) + TrackletConfigBuilder::LayerName(N_LAYER) + ".tab";
1139     }
1140 
1141   } else {
1142     if (type == VMRTableType::me) {
1143       //This if a hack where the same memory is used in both ME and TE modules
1144       if (layerdisk == LayerDisk::L2 || layerdisk == LayerDisk::L3 || layerdisk == LayerDisk::L4 ||
1145           layerdisk == LayerDisk::L6) {
1146         nbits_ = 6;
1147         positive_ = false;
1148         name_ = "VMTableOuter" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
1149         writeTable();
1150       }
1151 
1152       assert(region >= 0);
1153       char cregion = 'A' + region;
1154       name_ = "VMR_" + TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_finebin.tab";
1155       nbits_ = 6;
1156       positive_ = false;
1157     }
1158 
1159     if (type == VMRTableType::inner) {
1160       positive_ = false;
1161       nbits_ = 10;
1162       name_ = "VMTableInner" + TrackletConfigBuilder::LayerName(layerdisk) +
1163               TrackletConfigBuilder::LayerName(layerdisk + 1) + ".tab";
1164     }
1165 
1166     if (type == VMRTableType::inneroverlap) {
1167       positive_ = false;
1168       nbits_ = 10;
1169       name_ = "VMTableInner" + TrackletConfigBuilder::LayerName(layerdisk) + TrackletConfigBuilder::LayerName(N_LAYER) +
1170               ".tab";
1171     }
1172 
1173     if (type == VMRTableType::disk) {
1174       positive_ = false;
1175       nbits_ = 10;
1176       name_ = "VMTableOuter" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
1177     }
1178   }
1179 
1180   writeTable();
1181 }
1182 
1183 int TrackletLUT::getVMRLookup(unsigned int layerdisk, double z, double r, double dz, double dr, int iseed) const {
1184   double z0cut = settings_.z0cut();
1185 
1186   if (layerdisk < N_LAYER) {
1187     double constexpr zcutL2L3 = 52.0;  //Stubs closer to IP in z will not be used for L2L3 seeds
1188     if (iseed == Seed::L2L3 && std::abs(z) < zcutL2L3)
1189       return -1;
1190 
1191     double rmean = settings_.rmean(layerdisk);
1192 
1193     double rratio1 = rmean / (r + 0.5 * dr);
1194     double rratio2 = rmean / (r - 0.5 * dr);
1195 
1196     double z1 = (z - 0.5 * dz) * rratio1 + z0cut * (rratio1 - 1.0);
1197     double z2 = (z + 0.5 * dz) * rratio1 + z0cut * (rratio1 - 1.0);
1198     double z3 = (z - 0.5 * dz) * rratio2 + z0cut * (rratio2 - 1.0);
1199     double z4 = (z + 0.5 * dz) * rratio2 + z0cut * (rratio2 - 1.0);
1200     double z5 = (z - 0.5 * dz) * rratio1 - z0cut * (rratio1 - 1.0);
1201     double z6 = (z + 0.5 * dz) * rratio1 - z0cut * (rratio1 - 1.0);
1202     double z7 = (z - 0.5 * dz) * rratio2 - z0cut * (rratio2 - 1.0);
1203     double z8 = (z + 0.5 * dz) * rratio2 - z0cut * (rratio2 - 1.0);
1204 
1205     double zmin = std::min({z1, z2, z3, z4, z5, z6, z7, z8});
1206     double zmax = std::max({z1, z2, z3, z4, z5, z6, z7, z8});
1207 
1208     int NBINS = settings_.NLONGVMBINS() * settings_.NLONGVMBINS();
1209 
1210     int zbin1 = NBINS * (zmin + settings_.zlength()) / (2 * settings_.zlength());
1211     int zbin2 = NBINS * (zmax + settings_.zlength()) / (2 * settings_.zlength());
1212 
1213     if (zbin1 >= NBINS)
1214       return -1;
1215     if (zbin2 < 0)
1216       return -1;
1217 
1218     if (zbin2 >= NBINS)
1219       zbin2 = NBINS - 1;
1220     if (zbin1 < 0)
1221       zbin1 = 0;
1222 
1223     // This is a 10 bit word:
1224     // xxx|yyy|z|rrr
1225     // xxx is the delta z window
1226     // yyy is the z bin
1227     // z is flag to look in next bin
1228     // rrr first fine z bin
1229     // NOTE : this encoding is not efficient z is one if xxx+rrr is greater than 8
1230     //        and xxx is only 1,2, or 3
1231     //        should also reject xxx=0 as this means projection is outside range
1232 
1233     int value = zbin1 / 8;
1234     value *= 2;
1235     if (zbin2 / 8 - zbin1 / 8 > 0)
1236       value += 1;
1237     value *= 8;
1238     value += (zbin1 & 7);
1239     assert(value / 8 < 15);
1240     int deltaz = zbin2 - zbin1;
1241     if (deltaz > 7) {
1242       deltaz = 7;
1243     }
1244     assert(deltaz < 8);
1245     value += (deltaz << 7);
1246 
1247     return value;
1248 
1249   } else {
1250     if (std::abs(z) < 2.0 * z0cut)
1251       return -1;
1252 
1253     double zmean = settings_.zmean(layerdisk - N_LAYER);
1254     if (z < 0.0)
1255       zmean = -zmean;
1256 
1257     double r1 = (r + 0.5 * dr) * (zmean + z0cut) / (z + 0.5 * dz + z0cut);
1258     double r2 = (r - 0.5 * dr) * (zmean - z0cut) / (z + 0.5 * dz - z0cut);
1259     double r3 = (r + 0.5 * dr) * (zmean + z0cut) / (z - 0.5 * dz + z0cut);
1260     double r4 = (r - 0.5 * dr) * (zmean - z0cut) / (z - 0.5 * dz - z0cut);
1261     double r5 = (r + 0.5 * dr) * (zmean - z0cut) / (z + 0.5 * dz - z0cut);
1262     double r6 = (r - 0.5 * dr) * (zmean + z0cut) / (z + 0.5 * dz + z0cut);
1263     double r7 = (r + 0.5 * dr) * (zmean - z0cut) / (z - 0.5 * dz - z0cut);
1264     double r8 = (r - 0.5 * dr) * (zmean + z0cut) / (z - 0.5 * dz + z0cut);
1265 
1266     double rmin = std::min({r1, r2, r3, r4, r5, r6, r7, r8});
1267     double rmax = std::max({r1, r2, r3, r4, r5, r6, r7, r8});
1268 
1269     int NBINS = settings_.NLONGVMBINS() * settings_.NLONGVMBINS() / 2;
1270 
1271     double rmindisk = settings_.rmindiskvm();
1272     double rmaxdisk = settings_.rmaxdiskvm();
1273 
1274     if (iseed == Seed::L1D1)
1275       rmaxdisk = settings_.rmaxdiskl1overlapvm();
1276     if (iseed == Seed::L2D1)
1277       rmindisk = settings_.rmindiskl2overlapvm();
1278     if (iseed == Seed::L2L3D1)
1279       rmaxdisk = settings_.rmaxdisk();
1280 
1281     if (rmin > rmaxdisk)
1282       return -1;
1283     if (rmax > rmaxdisk)
1284       rmax = rmaxdisk;
1285 
1286     if (rmax < rmindisk)
1287       return -1;
1288     if (rmin < rmindisk)
1289       rmin = rmindisk;
1290 
1291     int rbin1 = NBINS * (rmin - settings_.rmindiskvm()) / (settings_.rmaxdiskvm() - settings_.rmindiskvm());
1292     int rbin2 = NBINS * (rmax - settings_.rmindiskvm()) / (settings_.rmaxdiskvm() - settings_.rmindiskvm());
1293 
1294     if (iseed == Seed::L2L3D1) {
1295       constexpr double rminspec = 40.0;
1296       rbin1 = NBINS * (rmin - rminspec) / (settings_.rmaxdisk() - rminspec);
1297       rbin2 = NBINS * (rmax - rminspec) / (settings_.rmaxdisk() - rminspec);
1298     }
1299 
1300     if (rbin2 >= NBINS)
1301       rbin2 = NBINS - 1;
1302     if (rbin1 < 0)
1303       rbin1 = 0;
1304 
1305     // This is a 9 bit word:
1306     // xxx|yy|z|rrr
1307     // xxx is the delta r window
1308     // yy is the r bin yy is three bits for overlaps
1309     // z is flag to look in next bin
1310     // rrr fine r bin
1311     // NOTE : this encoding is not efficient z is one if xxx+rrr is greater than 8
1312     //        and xxx is only 1,2, or 3
1313     //        should also reject xxx=0 as this means projection is outside range
1314 
1315     bool overlap = iseed == Seed::L1D1 || iseed == Seed::L2D1 || iseed == Seed::L2L3D1;
1316 
1317     int value = rbin1 / 8;
1318     if (overlap) {
1319       if (z < 0.0)
1320         value += 4;
1321     }
1322     value *= 2;
1323     if (rbin2 / 8 - rbin1 / 8 > 0)
1324       value += 1;
1325     value *= 8;
1326     value += (rbin1 & 7);
1327     assert(value / 8 < 15);
1328     int deltar = rbin2 - rbin1;
1329     if (deltar > 7)
1330       deltar = 7;
1331     if (overlap) {
1332       value += (deltar << 7);
1333     } else {
1334       value += (deltar << 6);
1335     }
1336 
1337     return value;
1338   }
1339 }
1340 
1341 void TrackletLUT::initPhiCorrTable(unsigned int layerdisk, unsigned int rbits) {
1342   bool psmodule = layerdisk < N_PSLAYER;
1343 
1344   unsigned int bendbits = psmodule ? N_BENDBITS_PS : N_BENDBITS_2S;
1345 
1346   unsigned int rbins = (1 << rbits);
1347 
1348   double rmean = settings_.rmean(layerdisk);
1349   double drmax = settings_.drmax();
1350 
1351   double dr = 2.0 * drmax / rbins;
1352 
1353   std::vector<std::array<double, 2>> bend_vals;
1354 
1355   if (settings_.useCalcBendCuts) {
1356     std::vector<const tt::SensorModule*> sm = getSensorModules(layerdisk, psmodule);
1357     bend_vals = getBendCut(layerdisk, sm, psmodule);
1358 
1359   } else {
1360     for (int ibend = 0; ibend < 1 << bendbits; ibend++) {
1361       bend_vals.push_back({{settings_.benddecode(ibend, layerdisk, layerdisk < N_PSLAYER), 0}});
1362     }
1363   }
1364 
1365   for (int ibend = 0; ibend < 1 << bendbits; ibend++) {
1366     for (unsigned int irbin = 0; irbin < rbins; irbin++) {
1367       double bend = -bend_vals[ibend][0];
1368       int value = getphiCorrValue(layerdisk, bend, irbin, rmean, dr, drmax);
1369       table_.push_back(value);
1370     }
1371   }
1372 
1373   name_ = "VMPhiCorrL" + std::to_string(layerdisk + 1) + ".tab";
1374   nbits_ = 14;
1375   positive_ = false;
1376 
1377   writeTable();
1378 }
1379 
1380 int TrackletLUT::getphiCorrValue(
1381     unsigned int layerdisk, double bend, unsigned int irbin, double rmean, double dr, double drmax) const {
1382   bool psmodule = layerdisk < N_PSLAYER;
1383 
1384   //for the rbin - calculate the distance to the nominal layer radius
1385   double Delta = (irbin + 0.5) * dr - drmax;
1386 
1387   //calculate the phi correction - this is a somewhat approximate formula
1388   double drnom = 0.18;  //This is the nominal module separation for which bend is referenced
1389   double dphi = (Delta / drnom) * bend * settings_.stripPitch(psmodule) / rmean;
1390 
1391   double kphi = psmodule ? settings_.kphi() : settings_.kphi1();
1392 
1393   int idphi = dphi / kphi;
1394 
1395   return idphi;
1396 }
1397 
1398 // Write LUT table.
1399 void TrackletLUT::writeTable() const {
1400   if (name_.empty()) {
1401     return;
1402   }
1403 
1404   if (nbits_ == 0) {
1405     throw cms::Exception("LogicError") << "Error in " << __FILE__ << " nbits_ == 0 ";
1406   }
1407 
1408   if (!settings_.writeTable()) {
1409     return;
1410   }
1411 
1412   ofstream out = openfile(settings_.tablePath(), name_, __FILE__, __LINE__);
1413 
1414   out << "{" << endl;
1415   for (unsigned int i = 0; i < table_.size(); i++) {
1416     if (i != 0) {
1417       out << "," << endl;
1418     }
1419 
1420     int itable = table_[i];
1421     if (positive_) {
1422       if (table_[i] < 0) {
1423         itable = (1 << nbits_) - 1;
1424       }
1425     }
1426 
1427     out << itable;
1428   }
1429   out << endl << "};" << endl;
1430   out.close();
1431 
1432   string name = name_;
1433 
1434   name[name_.size() - 3] = 'd';
1435   name[name_.size() - 2] = 'a';
1436   name[name_.size() - 1] = 't';
1437 
1438   out = openfile(settings_.tablePath(), name, __FILE__, __LINE__);
1439 
1440   int width = (nbits_ + 3) / 4;
1441 
1442   for (unsigned int i = 0; i < table_.size(); i++) {
1443     int itable = table_[i];
1444     if (positive_) {
1445       if (table_[i] < 0) {
1446         itable = (1 << nbits_) - 1;
1447       }
1448     }
1449 
1450     out << uppercase << setfill('0') << setw(width) << hex << itable << dec << endl;
1451   }
1452 
1453   out.close();
1454 }
1455 
1456 int TrackletLUT::lookup(unsigned int index) const {
1457   if (index >= table_.size()) {
1458     throw cms::Exception("LogicError") << "Error in " << __FILE__ << " index >= size " << index << " " << table_.size()
1459                                        << " in " << name_;
1460   }
1461   assert(index < table_.size());
1462   return table_[index];
1463 }
[ 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!