Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-15 01:07:51

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/L1TCommon/interface/BitShift.h"
0006 
0007 #include <filesystem>
0008 
0009 using namespace std;
0010 using namespace trklet;
0011 
0012 TrackletLUT::TrackletLUT(const Settings& settings) : settings_(settings) {}
0013 
0014 void TrackletLUT::initmatchcut(unsigned int layerdisk, MatchType type, unsigned int region) {
0015   char cregion = 'A' + region;
0016 
0017   for (unsigned int iSeed = 0; iSeed < 12; iSeed++) {
0018     if (type == barrelphi) {
0019       table_.push_back(settings_.rphimatchcut(iSeed, layerdisk) / (settings_.kphi1() * settings_.rmean(layerdisk)));
0020     }
0021     if (type == barrelz) {
0022       table_.push_back(settings_.zmatchcut(iSeed, layerdisk) / settings_.kz());
0023     }
0024     if (type == diskPSphi) {
0025       table_.push_back(settings_.rphicutPS(iSeed, layerdisk - N_LAYER) / (settings_.kphi() * settings_.kr()));
0026     }
0027     if (type == disk2Sphi) {
0028       table_.push_back(settings_.rphicut2S(iSeed, layerdisk - N_LAYER) / (settings_.kphi() * settings_.kr()));
0029     }
0030     if (type == disk2Sr) {
0031       table_.push_back(settings_.rcut2S(iSeed, layerdisk - N_LAYER) / settings_.krprojshiftdisk());
0032     }
0033     if (type == diskPSr) {
0034       table_.push_back(settings_.rcutPS(iSeed, layerdisk - N_LAYER) / settings_.krprojshiftdisk());
0035     }
0036   }
0037 
0038   name_ = settings_.combined() ? "MP_" : "MC_";
0039 
0040   if (type == barrelphi) {
0041     name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_phicut.tab";
0042   }
0043   if (type == barrelz) {
0044     name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_zcut.tab";
0045   }
0046   if (type == diskPSphi) {
0047     name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_PSphicut.tab";
0048   }
0049   if (type == disk2Sphi) {
0050     name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_2Sphicut.tab";
0051   }
0052   if (type == disk2Sr) {
0053     name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_2Srcut.tab";
0054   }
0055   if (type == diskPSr) {
0056     name_ += TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_PSrcut.tab";
0057   }
0058 
0059   positive_ = false;
0060 
0061   writeTable();
0062 }
0063 
0064 void TrackletLUT::initTPlut(bool fillInner,
0065                             unsigned int iSeed,
0066                             unsigned int layerdisk1,
0067                             unsigned int layerdisk2,
0068                             unsigned int nbitsfinephidiff,
0069                             unsigned int iTP) {
0070   //number of fine phi bins in sector
0071   int nfinephibins = settings_.nallstubs(layerdisk2) * settings_.nvmte(1, iSeed) * (1 << settings_.nfinephi(1, iSeed));
0072   double dfinephi = settings_.dphisectorHG() / nfinephibins;
0073 
0074   int outerrbits = 3;
0075 
0076   if (iSeed == Seed::L1L2 || iSeed == Seed::L2L3 || iSeed == Seed::L3L4 || iSeed == Seed::L5L6) {
0077     outerrbits = 0;
0078   }
0079 
0080   int outerrbins = (1 << outerrbits);
0081 
0082   double dphi[2];
0083   double router[2];
0084 
0085   unsigned int nbendbitsinner = 3;
0086   unsigned int nbendbitsouter = 3;
0087   if (iSeed == Seed::L3L4) {
0088     nbendbitsouter = 4;
0089   } else if (iSeed == Seed::L5L6) {
0090     nbendbitsinner = 4;
0091     nbendbitsouter = 4;
0092   }
0093 
0094   int nbinsfinephidiff = (1 << nbitsfinephidiff);
0095 
0096   for (int iphibin = 0; iphibin < nbinsfinephidiff; iphibin++) {
0097     int iphidiff = iphibin;
0098     if (iphibin >= nbinsfinephidiff / 2) {
0099       iphidiff = iphibin - nbinsfinephidiff;
0100     }
0101     //min and max dphi
0102     dphi[0] = (iphidiff - 1.5) * dfinephi;
0103     dphi[1] = (iphidiff + 1.5) * dfinephi;
0104     for (int irouterbin = 0; irouterbin < outerrbins; irouterbin++) {
0105       if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4 || iSeed == Seed::L1D1 || iSeed == Seed::L2D1) {
0106         router[0] =
0107             settings_.rmindiskvm() + irouterbin * (settings_.rmaxdiskvm() - settings_.rmindiskvm()) / outerrbins;
0108         router[1] =
0109             settings_.rmindiskvm() + (irouterbin + 1) * (settings_.rmaxdiskvm() - settings_.rmindiskvm()) / outerrbins;
0110       } else {
0111         router[0] = settings_.rmean(layerdisk2);
0112         router[1] = settings_.rmean(layerdisk2);
0113       }
0114 
0115       double bendinnermin = 20.0;
0116       double bendinnermax = -20.0;
0117       double bendoutermin = 20.0;
0118       double bendoutermax = -20.0;
0119       double rinvmin = 1.0;
0120       for (int i2 = 0; i2 < 2; i2++) {
0121         for (int i3 = 0; i3 < 2; i3++) {
0122           double rinner = 0.0;
0123           if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4) {
0124             rinner = router[i3] * settings_.zmean(layerdisk1 - N_LAYER) / settings_.zmean(layerdisk2 - N_LAYER);
0125           } else {
0126             rinner = settings_.rmean(layerdisk1);
0127           }
0128           double rinv1 = rinv(0.0, -dphi[i2], rinner, router[i3]);
0129           double pitchinner = (rinner < settings_.rcrit()) ? settings_.stripPitch(true) : settings_.stripPitch(false);
0130           double pitchouter =
0131               (router[i3] < settings_.rcrit()) ? settings_.stripPitch(true) : settings_.stripPitch(false);
0132           double abendinner = bendstrip(rinner, rinv1, pitchinner);
0133           double abendouter = bendstrip(router[i3], rinv1, pitchouter);
0134           if (abendinner < bendinnermin)
0135             bendinnermin = abendinner;
0136           if (abendinner > bendinnermax)
0137             bendinnermax = abendinner;
0138           if (abendouter < bendoutermin)
0139             bendoutermin = abendouter;
0140           if (abendouter > bendoutermax)
0141             bendoutermax = abendouter;
0142           if (std::abs(rinv1) < rinvmin) {
0143             rinvmin = std::abs(rinv1);
0144           }
0145         }
0146       }
0147 
0148       bool passptcut = rinvmin < settings_.rinvcutte();
0149 
0150       if (fillInner) {
0151         for (int ibend = 0; ibend < (1 << nbendbitsinner); ibend++) {
0152           double bend = settings_.benddecode(ibend, layerdisk1, nbendbitsinner == 3);
0153 
0154           bool passinner = bend <= bendinnermax + settings_.bendcutte(ibend, layerdisk1, nbendbitsinner == 3) &&
0155                            bend >= bendinnermin - settings_.bendcutte(ibend, layerdisk1, nbendbitsinner == 3);
0156           table_.push_back(passinner && passptcut);
0157         }
0158       } else {
0159         for (int ibend = 0; ibend < (1 << nbendbitsouter); ibend++) {
0160           double bend = settings_.benddecode(ibend, layerdisk2, nbendbitsouter == 3);
0161 
0162           bool passouter = bend <= bendoutermax + settings_.bendcutte(ibend, layerdisk2, nbendbitsouter == 3) &&
0163                            bend >= bendoutermin - settings_.bendcutte(ibend, layerdisk2, nbendbitsouter == 3);
0164           table_.push_back(passouter && passptcut);
0165         }
0166       }
0167     }
0168   }
0169 
0170   positive_ = false;
0171   char cTP = 'A' + iTP;
0172 
0173   name_ = "TP_" + TrackletConfigBuilder::LayerName(layerdisk1) + TrackletConfigBuilder::LayerName(layerdisk2) + cTP;
0174 
0175   if (fillInner) {
0176     name_ += "_stubptinnercut.tab";
0177   } else {
0178     name_ += "_stubptoutercut.tab";
0179   }
0180 
0181   writeTable();
0182 }
0183 
0184 void TrackletLUT::initTPregionlut(unsigned int iSeed,
0185                                   unsigned int layerdisk1,
0186                                   unsigned int layerdisk2,
0187                                   unsigned int iAllStub,
0188                                   unsigned int nbitsfinephidiff,
0189                                   unsigned int nbitsfinephi,
0190                                   const TrackletLUT& tplutinner,
0191                                   unsigned int iTP) {
0192   int nirbits = 0;
0193   if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4 || iSeed == Seed::L1D1 || iSeed == Seed::L2D1) {
0194     nirbits = 3;
0195   }
0196 
0197   unsigned int nbendbitsinner = 3;
0198 
0199   if (iSeed == Seed::L5L6) {
0200     nbendbitsinner = 4;
0201   }
0202 
0203   for (int innerfinephi = 0; innerfinephi < (1 << nbitsfinephi); innerfinephi++) {
0204     for (int innerbend = 0; innerbend < (1 << nbendbitsinner); innerbend++) {
0205       for (int ir = 0; ir < (1 << nirbits); ir++) {
0206         unsigned int usereg = 0;
0207         for (unsigned int ireg = 0; ireg < settings_.nvmte(1, iSeed); ireg++) {
0208           bool match = false;
0209           for (int ifinephiouter = 0; ifinephiouter < (1 << settings_.nfinephi(1, iSeed)); ifinephiouter++) {
0210             int outerfinephi = iAllStub * (1 << (nbitsfinephi - settings_.nbitsallstubs(layerdisk2))) +
0211                                ireg * (1 << settings_.nfinephi(1, iSeed)) + ifinephiouter;
0212             int idphi = outerfinephi - innerfinephi;
0213             bool inrange = (idphi < (1 << (nbitsfinephidiff - 1))) && (idphi >= -(1 << (nbitsfinephidiff - 1)));
0214             if (idphi < 0)
0215               idphi = idphi + (1 << nbitsfinephidiff);
0216             int idphi1 = idphi;
0217             if (iSeed >= 4)
0218               idphi1 = (idphi << 3) + ir;
0219             int ptinnerindexnew = l1t::bitShift(idphi1, nbendbitsinner) + innerbend;
0220             match = match || (inrange && tplutinner.lookup(ptinnerindexnew));
0221           }
0222           if (match) {
0223             usereg = usereg | (1 << ireg);
0224           }
0225         }
0226 
0227         table_.push_back(usereg);
0228       }
0229     }
0230   }
0231 
0232   positive_ = false;
0233   char cTP = 'A' + iTP;
0234 
0235   name_ = "TP_" + TrackletConfigBuilder::LayerName(layerdisk1) + TrackletConfigBuilder::LayerName(layerdisk2) + cTP +
0236           "_usereg.tab";
0237 
0238   writeTable();
0239 }
0240 
0241 void TrackletLUT::initteptlut(bool fillInner,
0242                               bool fillTEMem,
0243                               unsigned int iSeed,
0244                               unsigned int layerdisk1,
0245                               unsigned int layerdisk2,
0246                               unsigned int innerphibits,
0247                               unsigned int outerphibits,
0248                               double innerphimin,
0249                               double innerphimax,
0250                               double outerphimin,
0251                               double outerphimax,
0252                               const std::string& innermem,
0253                               const std::string& outermem) {
0254   int outerrbits = 0;
0255   if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4 || iSeed == Seed::L1D1 || iSeed == Seed::L2D1) {
0256     outerrbits = 3;
0257   }
0258 
0259   int outerrbins = (1 << outerrbits);
0260   int innerphibins = (1 << innerphibits);
0261   int outerphibins = (1 << outerphibits);
0262 
0263   double phiinner[2];
0264   double phiouter[2];
0265   double router[2];
0266 
0267   unsigned int nbendbitsinner = 3;
0268   unsigned int nbendbitsouter = 3;
0269   if (iSeed == Seed::L3L4) {
0270     nbendbitsouter = 4;
0271   }
0272   if (iSeed == Seed::L5L6) {
0273     nbendbitsinner = 4;
0274     nbendbitsouter = 4;
0275   }
0276 
0277   if (fillTEMem) {
0278     if (fillInner) {
0279       table_.resize((1 << nbendbitsinner), false);
0280     } else {
0281       table_.resize((1 << nbendbitsouter), false);
0282     }
0283   }
0284 
0285   for (int iphiinnerbin = 0; iphiinnerbin < innerphibins; iphiinnerbin++) {
0286     phiinner[0] = innerphimin + iphiinnerbin * (innerphimax - innerphimin) / innerphibins;
0287     phiinner[1] = innerphimin + (iphiinnerbin + 1) * (innerphimax - innerphimin) / innerphibins;
0288     for (int iphiouterbin = 0; iphiouterbin < outerphibins; iphiouterbin++) {
0289       phiouter[0] = outerphimin + iphiouterbin * (outerphimax - outerphimin) / outerphibins;
0290       phiouter[1] = outerphimin + (iphiouterbin + 1) * (outerphimax - outerphimin) / outerphibins;
0291       for (int irouterbin = 0; irouterbin < outerrbins; irouterbin++) {
0292         if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4 || iSeed == Seed::L1D1 || iSeed == Seed::L2D1) {
0293           router[0] =
0294               settings_.rmindiskvm() + irouterbin * (settings_.rmaxdiskvm() - settings_.rmindiskvm()) / outerrbins;
0295           router[1] = settings_.rmindiskvm() +
0296                       (irouterbin + 1) * (settings_.rmaxdiskvm() - settings_.rmindiskvm()) / outerrbins;
0297         } else {
0298           router[0] = settings_.rmean(layerdisk2);
0299           router[1] = settings_.rmean(layerdisk2);
0300         }
0301 
0302         double bendinnermin = 20.0;
0303         double bendinnermax = -20.0;
0304         double bendoutermin = 20.0;
0305         double bendoutermax = -20.0;
0306         double rinvmin = 1.0;
0307         for (int i1 = 0; i1 < 2; i1++) {
0308           for (int i2 = 0; i2 < 2; i2++) {
0309             for (int i3 = 0; i3 < 2; i3++) {
0310               double rinner = 0.0;
0311               if (iSeed == Seed::D1D2 || iSeed == Seed::D3D4) {
0312                 rinner = router[i3] * settings_.zmean(layerdisk1 - N_LAYER) / settings_.zmean(layerdisk2 - N_LAYER);
0313               } else {
0314                 rinner = settings_.rmean(layerdisk1);
0315               }
0316               double rinv1 = -rinv(phiinner[i1], phiouter[i2], rinner, router[i3]);
0317               double pitchinner =
0318                   (rinner < settings_.rcrit()) ? settings_.stripPitch(true) : settings_.stripPitch(false);
0319               double pitchouter =
0320                   (router[i3] < settings_.rcrit()) ? settings_.stripPitch(true) : settings_.stripPitch(false);
0321               double abendinner = bendstrip(rinner, rinv1, pitchinner);
0322               double abendouter = bendstrip(router[i3], rinv1, pitchouter);
0323               if (abendinner < bendinnermin)
0324                 bendinnermin = abendinner;
0325               if (abendinner > bendinnermax)
0326                 bendinnermax = abendinner;
0327               if (abendouter < bendoutermin)
0328                 bendoutermin = abendouter;
0329               if (abendouter > bendoutermax)
0330                 bendoutermax = abendouter;
0331               if (std::abs(rinv1) < rinvmin) {
0332                 rinvmin = std::abs(rinv1);
0333               }
0334             }
0335           }
0336         }
0337 
0338         bool passptcut = rinvmin < settings_.rinvcutte();
0339 
0340         if (fillInner) {
0341           for (int ibend = 0; ibend < (1 << nbendbitsinner); ibend++) {
0342             double bend = settings_.benddecode(ibend, layerdisk1, nbendbitsinner == 3);
0343 
0344             bool passinner = bend > bendinnermin - settings_.bendcutte(ibend, layerdisk1, nbendbitsinner == 3) &&
0345                              bend < bendinnermax + settings_.bendcutte(ibend, layerdisk1, nbendbitsinner == 3);
0346 
0347             if (fillTEMem) {
0348               if (passinner) {
0349                 table_[ibend] = 1;
0350               }
0351             } else {
0352               table_.push_back(passinner && passptcut);
0353             }
0354           }
0355         } else {
0356           for (int ibend = 0; ibend < (1 << nbendbitsouter); ibend++) {
0357             double bend = settings_.benddecode(ibend, layerdisk2, nbendbitsouter == 3);
0358 
0359             bool passouter = bend > bendoutermin - settings_.bendcutte(ibend, layerdisk2, nbendbitsouter == 3) &&
0360                              bend < bendoutermax + settings_.bendcutte(ibend, layerdisk2, nbendbitsouter == 3);
0361             if (fillTEMem) {
0362               if (passouter) {
0363                 table_[ibend] = 1;
0364               }
0365             } else {
0366               table_.push_back(passouter && passptcut);
0367             }
0368           }
0369         }
0370       }
0371     }
0372   }
0373 
0374   positive_ = false;
0375 
0376   if (fillTEMem) {
0377     if (fillInner) {
0378       name_ = "VMSTE_" + innermem + "_vmbendcut.tab";
0379     } else {
0380       name_ = "VMSTE_" + outermem + "_vmbendcut.tab";
0381     }
0382   } else {
0383     name_ = "TE_" + innermem.substr(0, innermem.size() - 2) + "_" + outermem.substr(0, outermem.size() - 2);
0384     if (fillInner) {
0385       name_ += "_stubptinnercut.tab";
0386     } else {
0387       name_ += "_stubptoutercut.tab";
0388     }
0389   }
0390 
0391   writeTable();
0392 }
0393 
0394 void TrackletLUT::initProjectionBend(double k_phider,
0395                                      unsigned int idisk,
0396                                      unsigned int nrbits,
0397                                      unsigned int nphiderbits) {
0398   unsigned int nsignbins = 2;
0399   unsigned int nrbins = 1 << (nrbits);
0400   unsigned int nphiderbins = 1 << (nphiderbits);
0401 
0402   for (unsigned int isignbin = 0; isignbin < nsignbins; isignbin++) {
0403     for (unsigned int irbin = 0; irbin < nrbins; irbin++) {
0404       int ir = irbin;
0405       if (ir > (1 << (nrbits - 1)))
0406         ir -= (1 << nrbits);
0407       ir = l1t::bitShift(ir, (settings_.nrbitsstub(N_LAYER) - nrbits));
0408       for (unsigned int iphiderbin = 0; iphiderbin < nphiderbins; iphiderbin++) {
0409         int iphider = iphiderbin;
0410         if (iphider > (1 << (nphiderbits - 1)))
0411           iphider -= (1 << nphiderbits);
0412         iphider = l1t::bitShift(iphider, (settings_.nbitsphiprojderL123() - nphiderbits));
0413 
0414         double rproj = ir * settings_.krprojshiftdisk();
0415         double phider = iphider * k_phider;
0416         double t = settings_.zmean(idisk) / rproj;
0417 
0418         if (isignbin)
0419           t = -t;
0420 
0421         double rinv = -phider * (2.0 * t);
0422 
0423         double stripPitch = (rproj < settings_.rcrit()) ? settings_.stripPitch(true) : settings_.stripPitch(false);
0424         double bendproj = bendstrip(rproj, rinv, stripPitch);
0425 
0426         int ibendproj = 2.0 * bendproj + 15.5;
0427         if (ibendproj < 0)
0428           ibendproj = 0;
0429         if (ibendproj > 31)
0430           ibendproj = 31;
0431 
0432         table_.push_back(ibendproj);
0433       }
0434     }
0435   }
0436 
0437   positive_ = false;
0438   name_ = settings_.combined() ? "MP_" : "PR_";
0439   name_ += "ProjectionBend_" + TrackletConfigBuilder::LayerName(N_LAYER + idisk) + ".tab";
0440 
0441   writeTable();
0442 }
0443 
0444 void TrackletLUT::initBendMatch(unsigned int layerdisk) {
0445   unsigned int nrinv = NRINVBITS;
0446   double rinvhalf = 0.5 * ((1 << nrinv) - 1);
0447 
0448   bool barrel = layerdisk < N_LAYER;
0449   bool isPSmodule = layerdisk < N_PSLAYER;
0450   double stripPitch = settings_.stripPitch(isPSmodule);
0451 
0452   if (barrel) {
0453     unsigned int nbits = isPSmodule ? N_BENDBITS_PS : N_BENDBITS_2S;
0454 
0455     for (unsigned int irinv = 0; irinv < (1u << nrinv); irinv++) {
0456       double rinv = (irinv - rinvhalf) * (1 << (settings_.nbitsrinv() - nrinv)) * settings_.krinvpars();
0457 
0458       double projbend = bendstrip(settings_.rmean(layerdisk), rinv, stripPitch);
0459       for (unsigned int ibend = 0; ibend < (1u << nbits); ibend++) {
0460         double stubbend = settings_.benddecode(ibend, layerdisk, isPSmodule);
0461         bool pass = std::abs(stubbend - projbend) < settings_.bendcutme(ibend, layerdisk, isPSmodule);
0462         table_.push_back(pass);
0463       }
0464     }
0465   } else {
0466     for (unsigned int iprojbend = 0; iprojbend < (1u << nrinv); iprojbend++) {
0467       double projbend = 0.5 * (iprojbend - rinvhalf);
0468       for (unsigned int ibend = 0; ibend < (1 << N_BENDBITS_2S); ibend++) {
0469         double stubbend = settings_.benddecode(ibend, layerdisk, false);
0470         bool pass = std::abs(stubbend - projbend) < settings_.bendcutme(ibend, layerdisk, false);
0471         table_.push_back(pass);
0472       }
0473     }
0474     for (unsigned int iprojbend = 0; iprojbend < (1u << nrinv); iprojbend++) {
0475       double projbend = 0.5 * (iprojbend - rinvhalf);
0476       for (unsigned int ibend = 0; ibend < (1 << N_BENDBITS_PS); ibend++) {
0477         double stubbend = settings_.benddecode(ibend, layerdisk, true);
0478         bool pass = std::abs(stubbend - projbend) < settings_.bendcutme(ibend, layerdisk, true);
0479         table_.push_back(pass);
0480       }
0481     }
0482   }
0483 
0484   positive_ = false;
0485 
0486   name_ = "METable_" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
0487 
0488   writeTable();
0489 }
0490 
0491 void TrackletLUT::initVMRTable(unsigned int layerdisk, VMRTableType type, int region) {
0492   unsigned int zbits = settings_.vmrlutzbits(layerdisk);
0493   unsigned int rbits = settings_.vmrlutrbits(layerdisk);
0494 
0495   unsigned int rbins = (1 << rbits);
0496   unsigned int zbins = (1 << zbits);
0497 
0498   double zmin, zmax, rmin, rmax;
0499 
0500   if (layerdisk < N_LAYER) {
0501     zmin = -settings_.zlength();
0502     zmax = settings_.zlength();
0503     rmin = settings_.rmean(layerdisk) - settings_.drmax();
0504     rmax = settings_.rmean(layerdisk) + settings_.drmax();
0505   } else {
0506     rmin = 0;
0507     rmax = settings_.rmaxdisk();
0508     zmin = settings_.zmean(layerdisk - N_LAYER) - settings_.dzmax();
0509     zmax = settings_.zmean(layerdisk - N_LAYER) + settings_.dzmax();
0510   }
0511 
0512   double dr = (rmax - rmin) / rbins;
0513   double dz = (zmax - zmin) / zbins;
0514 
0515   int NBINS = settings_.NLONGVMBINS() * settings_.NLONGVMBINS();
0516 
0517   for (unsigned int izbin = 0; izbin < zbins; izbin++) {
0518     for (unsigned int irbin = 0; irbin < rbins; irbin++) {
0519       double r = rmin + (irbin + 0.5) * dr;
0520       double z = zmin + (izbin + 0.5) * dz;
0521 
0522       if (settings_.combined()) {
0523         int iznew = izbin - (1 << (zbits - 1));
0524         if (iznew < 0)
0525           iznew += (1 << zbits);
0526         assert(iznew >= 0);
0527         assert(iznew < (1 << zbits));
0528         z = zmin + (iznew + 0.5) * dz;
0529         if (layerdisk < N_LAYER) {
0530           int irnew = irbin - (1 << (rbits - 1));
0531           if (irnew < 0)
0532             irnew += (1 << rbits);
0533           assert(irnew >= 0);
0534           assert(irnew < (1 << rbits));
0535           r = rmin + (irnew + 0.5) * dr;
0536         }
0537       }
0538 
0539       if (layerdisk >= N_LAYER && irbin < 10)  //special case for the tabulated radii in 2S disks
0540         r = (layerdisk < N_LAYER + 2) ? settings_.rDSSinner(irbin) : settings_.rDSSouter(irbin);
0541 
0542       int bin;
0543       if (layerdisk < N_LAYER) {
0544         double zproj = z * settings_.rmean(layerdisk) / r;
0545         bin = NBINS * (zproj + settings_.zlength()) / (2 * settings_.zlength());
0546       } else {
0547         double rproj = r * settings_.zmean(layerdisk - N_LAYER) / z;
0548         bin = NBINS * (rproj - settings_.rmindiskvm()) / (settings_.rmaxdisk() - settings_.rmindiskvm());
0549       }
0550       if (bin < 0)
0551         bin = 0;
0552       if (bin >= NBINS)
0553         bin = NBINS - 1;
0554 
0555       if (type == VMRTableType::me) {
0556         table_.push_back(bin);
0557       }
0558 
0559       if (type == VMRTableType::disk) {
0560         if (layerdisk >= N_LAYER) {
0561           double rproj = r * settings_.zmean(layerdisk - N_LAYER) / z;
0562           bin = 0.5 * NBINS * (rproj - settings_.rmindiskvm()) / (settings_.rmaxdiskvm() - settings_.rmindiskvm());
0563           //bin value of zero indicates that stub is out of range
0564           if (bin < 0)
0565             bin = 0;
0566           if (bin >= NBINS / 2)
0567             bin = 0;
0568           table_.push_back(bin);
0569         }
0570       }
0571 
0572       if (type == VMRTableType::inner) {
0573         if (layerdisk == LayerDisk::L1 || layerdisk == LayerDisk::L3 || layerdisk == LayerDisk::L5 ||
0574             layerdisk == LayerDisk::D1 || layerdisk == LayerDisk::D3) {
0575           table_.push_back(getVMRLookup(layerdisk + 1, z, r, dz, dr));
0576         }
0577         if (layerdisk == LayerDisk::L2) {
0578           table_.push_back(getVMRLookup(layerdisk + 1, z, r, dz, dr, Seed::L2L3));
0579         }
0580       }
0581 
0582       if (type == VMRTableType::inneroverlap) {
0583         if (layerdisk == LayerDisk::L1 || layerdisk == LayerDisk::L2) {
0584           table_.push_back(getVMRLookup(6, z, r, dz, dr, layerdisk + 6));
0585         }
0586       }
0587 
0588       if (type == VMRTableType::innerthird) {
0589         if (layerdisk == LayerDisk::L2) {  //projection from L2 to D1 for L2L3D1 seeding
0590           table_.push_back(getVMRLookup(LayerDisk::D1, z, r, dz, dr, Seed::L2L3D1));
0591         }
0592 
0593         if (layerdisk == LayerDisk::L5) {  //projection from L5 to L4 for L5L6L4 seeding
0594           table_.push_back(getVMRLookup(LayerDisk::L4, z, r, dz, dr));
0595         }
0596 
0597         if (layerdisk == LayerDisk::L3) {  //projection from L3 to L5 for L3L4L2 seeding
0598           table_.push_back(getVMRLookup(LayerDisk::L2, z, r, dz, dr));
0599         }
0600 
0601         if (layerdisk == LayerDisk::D1) {  //projection from D1 to L2 for D1D2L2 seeding
0602           table_.push_back(getVMRLookup(LayerDisk::L2, z, r, dz, dr));
0603         }
0604       }
0605     }
0606   }
0607 
0608   if (settings_.combined()) {
0609     if (type == VMRTableType::me) {
0610       positive_ = false;
0611       name_ = "VMRME_" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
0612     }
0613     if (type == VMRTableType::disk) {
0614       positive_ = false;
0615       name_ = "VMRTE_" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
0616     }
0617     if (type == VMRTableType::inner) {
0618       positive_ = true;
0619       nbits_ = 10;
0620       name_ = "TP_" + TrackletConfigBuilder::LayerName(layerdisk) + TrackletConfigBuilder::LayerName(layerdisk + 1) +
0621               ".tab";
0622     }
0623 
0624     if (type == VMRTableType::inneroverlap) {
0625       positive_ = true;
0626       nbits_ = 10;
0627       name_ = "TP_" + TrackletConfigBuilder::LayerName(layerdisk) + TrackletConfigBuilder::LayerName(N_LAYER) + ".tab";
0628     }
0629 
0630   } else {
0631     if (type == VMRTableType::me) {
0632       //This if a hack where the same memory is used in both ME and TE modules
0633       if (layerdisk == 1 || layerdisk == 2 || layerdisk == 3 || layerdisk == 4) {
0634         positive_ = false;
0635         name_ = "VMTableOuter" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
0636         writeTable();
0637       }
0638 
0639       assert(region >= 0);
0640       char cregion = 'A' + region;
0641       name_ = "VMR_" + TrackletConfigBuilder::LayerName(layerdisk) + "PHI" + cregion + "_finebin.tab";
0642       positive_ = false;
0643     }
0644 
0645     if (type == VMRTableType::inner) {
0646       positive_ = false;
0647       name_ = "VMTableInner" + TrackletConfigBuilder::LayerName(layerdisk) +
0648               TrackletConfigBuilder::LayerName(layerdisk + 1) + ".tab";
0649     }
0650 
0651     if (type == VMRTableType::inneroverlap) {
0652       positive_ = false;
0653       name_ = "VMTableInner" + TrackletConfigBuilder::LayerName(layerdisk) + TrackletConfigBuilder::LayerName(N_LAYER) +
0654               ".tab";
0655     }
0656 
0657     if (type == VMRTableType::disk) {
0658       positive_ = false;
0659       name_ = "VMTableOuter" + TrackletConfigBuilder::LayerName(layerdisk) + ".tab";
0660     }
0661   }
0662 
0663   writeTable();
0664 }
0665 
0666 int TrackletLUT::getVMRLookup(unsigned int layerdisk, double z, double r, double dz, double dr, int iseed) const {
0667   double z0cut = settings_.z0cut();
0668 
0669   if (layerdisk < N_LAYER) {
0670     if (iseed == Seed::L2L3 && std::abs(z) < 52.0)
0671       return -1;
0672 
0673     double rmean = settings_.rmean(layerdisk);
0674 
0675     double rratio1 = rmean / (r + 0.5 * dr);
0676     double rratio2 = rmean / (r - 0.5 * dr);
0677 
0678     double z1 = (z - 0.5 * dz) * rratio1 + z0cut * (rratio1 - 1.0);
0679     double z2 = (z + 0.5 * dz) * rratio1 + z0cut * (rratio1 - 1.0);
0680     double z3 = (z - 0.5 * dz) * rratio2 + z0cut * (rratio2 - 1.0);
0681     double z4 = (z + 0.5 * dz) * rratio2 + z0cut * (rratio2 - 1.0);
0682     double z5 = (z - 0.5 * dz) * rratio1 - z0cut * (rratio1 - 1.0);
0683     double z6 = (z + 0.5 * dz) * rratio1 - z0cut * (rratio1 - 1.0);
0684     double z7 = (z - 0.5 * dz) * rratio2 - z0cut * (rratio2 - 1.0);
0685     double z8 = (z + 0.5 * dz) * rratio2 - z0cut * (rratio2 - 1.0);
0686 
0687     double zmin = std::min({z1, z2, z3, z4, z5, z6, z7, z8});
0688     double zmax = std::max({z1, z2, z3, z4, z5, z6, z7, z8});
0689 
0690     int NBINS = settings_.NLONGVMBINS() * settings_.NLONGVMBINS();
0691 
0692     int zbin1 = NBINS * (zmin + settings_.zlength()) / (2 * settings_.zlength());
0693     int zbin2 = NBINS * (zmax + settings_.zlength()) / (2 * settings_.zlength());
0694 
0695     if (zbin1 >= NBINS)
0696       return -1;
0697     if (zbin2 < 0)
0698       return -1;
0699 
0700     if (zbin2 >= NBINS)
0701       zbin2 = NBINS - 1;
0702     if (zbin1 < 0)
0703       zbin1 = 0;
0704 
0705     // This is a 10 bit word:
0706     // xxx|yyy|z|rrr
0707     // xxx is the delta z window
0708     // yyy is the z bin
0709     // z is flag to look in next bin
0710     // rrr first fine z bin
0711     // NOTE : this encoding is not efficient z is one if xxx+rrr is greater than 8
0712     //        and xxx is only 1,2, or 3
0713     //        should also reject xxx=0 as this means projection is outside range
0714 
0715     int value = zbin1 / 8;
0716     value *= 2;
0717     if (zbin2 / 8 - zbin1 / 8 > 0)
0718       value += 1;
0719     value *= 8;
0720     value += (zbin1 & 7);
0721     assert(value / 8 < 15);
0722     int deltaz = zbin2 - zbin1;
0723     if (deltaz > 7) {
0724       deltaz = 7;
0725     }
0726     assert(deltaz < 8);
0727     value += (deltaz << 7);
0728 
0729     return value;
0730 
0731   } else {
0732     if (std::abs(z) < 2.0 * z0cut)
0733       return -1;
0734 
0735     double zmean = settings_.zmean(layerdisk - N_LAYER);
0736     if (z < 0.0)
0737       zmean = -zmean;
0738 
0739     double r1 = (r + 0.5 * dr) * (zmean + z0cut) / (z + 0.5 * dz + z0cut);
0740     double r2 = (r - 0.5 * dr) * (zmean - z0cut) / (z + 0.5 * dz - z0cut);
0741     double r3 = (r + 0.5 * dr) * (zmean + z0cut) / (z - 0.5 * dz + z0cut);
0742     double r4 = (r - 0.5 * dr) * (zmean - z0cut) / (z - 0.5 * dz - z0cut);
0743     double r5 = (r + 0.5 * dr) * (zmean - z0cut) / (z + 0.5 * dz - z0cut);
0744     double r6 = (r - 0.5 * dr) * (zmean + z0cut) / (z + 0.5 * dz + z0cut);
0745     double r7 = (r + 0.5 * dr) * (zmean - z0cut) / (z - 0.5 * dz - z0cut);
0746     double r8 = (r - 0.5 * dr) * (zmean + z0cut) / (z - 0.5 * dz + z0cut);
0747 
0748     double rmin = std::min({r1, r2, r3, r4, r5, r6, r7, r8});
0749     double rmax = std::max({r1, r2, r3, r4, r5, r6, r7, r8});
0750 
0751     int NBINS = settings_.NLONGVMBINS() * settings_.NLONGVMBINS() / 2;
0752 
0753     double rmindisk = settings_.rmindiskvm();
0754     double rmaxdisk = settings_.rmaxdiskvm();
0755 
0756     if (iseed == Seed::L1D1)
0757       rmaxdisk = settings_.rmaxdiskl1overlapvm();
0758     if (iseed == Seed::L2D1)
0759       rmindisk = settings_.rmindiskl2overlapvm();
0760     if (iseed == Seed::L2L3D1)
0761       rmaxdisk = settings_.rmaxdisk();
0762 
0763     if (rmin > rmaxdisk)
0764       return -1;
0765     if (rmax > rmaxdisk)
0766       rmax = rmaxdisk;
0767 
0768     if (rmax < rmindisk)
0769       return -1;
0770     if (rmin < rmindisk)
0771       rmin = rmindisk;
0772 
0773     int rbin1 = NBINS * (rmin - settings_.rmindiskvm()) / (settings_.rmaxdiskvm() - settings_.rmindiskvm());
0774     int rbin2 = NBINS * (rmax - settings_.rmindiskvm()) / (settings_.rmaxdiskvm() - settings_.rmindiskvm());
0775 
0776     if (iseed == Seed::L2L3D1) {
0777       constexpr double rminspec = 40.0;
0778       rbin1 = NBINS * (rmin - rminspec) / (settings_.rmaxdisk() - rminspec);
0779       rbin2 = NBINS * (rmax - rminspec) / (settings_.rmaxdisk() - rminspec);
0780     }
0781 
0782     if (rbin2 >= NBINS)
0783       rbin2 = NBINS - 1;
0784     if (rbin1 < 0)
0785       rbin1 = 0;
0786 
0787     // This is a 9 bit word:
0788     // xxx|yy|z|rrr
0789     // xxx is the delta r window
0790     // yy is the r bin yy is three bits for overlaps
0791     // z is flag to look in next bin
0792     // rrr fine r bin
0793     // NOTE : this encoding is not efficient z is one if xxx+rrr is greater than 8
0794     //        and xxx is only 1,2, or 3
0795     //        should also reject xxx=0 as this means projection is outside range
0796 
0797     bool overlap = iseed == Seed::L1D1 || iseed == Seed::L2D1 || iseed == Seed::L2L3D1;
0798 
0799     int value = rbin1 / 8;
0800     if (overlap) {
0801       if (z < 0.0)
0802         value += 4;
0803     }
0804     value *= 2;
0805     if (rbin2 / 8 - rbin1 / 8 > 0)
0806       value += 1;
0807     value *= 8;
0808     value += (rbin1 & 7);
0809     assert(value / 8 < 15);
0810     int deltar = rbin2 - rbin1;
0811     if (deltar > 7)
0812       deltar = 7;
0813     if (overlap) {
0814       value += (deltar << 7);
0815     } else {
0816       value += (deltar << 6);
0817     }
0818 
0819     return value;
0820   }
0821 }
0822 
0823 void TrackletLUT::initPhiCorrTable(unsigned int layerdisk, unsigned int rbits) {
0824   bool psmodule = layerdisk < N_PSLAYER;
0825 
0826   unsigned int bendbits = psmodule ? N_BENDBITS_PS : N_BENDBITS_2S;
0827 
0828   unsigned int rbins = (1 << rbits);
0829 
0830   double rmean = settings_.rmean(layerdisk);
0831   double drmax = settings_.drmax();
0832 
0833   double dr = 2.0 * drmax / rbins;
0834 
0835   unsigned int bendbins = (1 << bendbits);
0836 
0837   for (unsigned int ibend = 0; ibend < bendbins; ibend++) {
0838     for (unsigned int irbin = 0; irbin < rbins; irbin++) {
0839       int value = getphiCorrValue(layerdisk, ibend, irbin, rmean, dr, drmax);
0840       table_.push_back(value);
0841     }
0842   }
0843 
0844   name_ = "VMPhiCorrL" + std::to_string(layerdisk + 1) + ".tab";
0845   nbits_ = 14;
0846   positive_ = false;
0847 
0848   writeTable();
0849 }
0850 
0851 int TrackletLUT::getphiCorrValue(
0852     unsigned int layerdisk, unsigned int ibend, unsigned int irbin, double rmean, double dr, double drmax) const {
0853   bool psmodule = layerdisk < N_PSLAYER;
0854 
0855   double bend = -settings_.benddecode(ibend, layerdisk, psmodule);
0856 
0857   //for the rbin - calculate the distance to the nominal layer radius
0858   double Delta = (irbin + 0.5) * dr - drmax;
0859 
0860   //calculate the phi correction - this is a somewhat approximate formula
0861   double dphi = (Delta / 0.18) * bend * settings_.stripPitch(psmodule) / rmean;
0862 
0863   double kphi = psmodule ? settings_.kphi() : settings_.kphi1();
0864 
0865   int idphi = dphi / kphi;
0866 
0867   return idphi;
0868 }
0869 
0870 // Write LUT table.
0871 void TrackletLUT::writeTable() const {
0872   if (!settings_.writeTable()) {
0873     return;
0874   }
0875 
0876   if (name_.empty()) {
0877     return;
0878   }
0879 
0880   ofstream out = openfile(settings_.tablePath(), name_, __FILE__, __LINE__);
0881 
0882   out << "{" << endl;
0883   for (unsigned int i = 0; i < table_.size(); i++) {
0884     if (i != 0) {
0885       out << "," << endl;
0886     }
0887 
0888     int itable = table_[i];
0889     if (positive_) {
0890       if (table_[i] < 0) {
0891         itable = (1 << nbits_) - 1;
0892       }
0893     }
0894 
0895     out << itable;
0896   }
0897   out << endl << "};" << endl;
0898   out.close();
0899 }
0900 
0901 int TrackletLUT::lookup(unsigned int index) const {
0902   assert(index < table_.size());
0903   return table_[index];
0904 }