Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-07 22:33:35

0001 #include "L1Trigger/TrackFindingTracklet/interface/TrackDerTable.h"
0002 #include "L1Trigger/TrackFindingTracklet/interface/FPGAWord.h"
0003 #include "L1Trigger/TrackFindingTracklet/interface/Globals.h"
0004 #include "L1Trigger/TrackFindingTracklet/interface/Util.h"
0005 
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 
0008 using namespace std;
0009 using namespace trklet;
0010 
0011 TrackDerTable::TrackDerTable(Settings const& settings) : settings_(settings) {
0012   Nlay_ = N_LAYER;
0013   Ndisk_ = N_DISK;
0014 
0015   LayerMemBits_ = 6;
0016   DiskMemBits_ = 7;
0017   LayerDiskMemBits_ = 18;
0018 
0019   alphaBits_ = settings_.alphaBitsTable();
0020 
0021   nextLayerValue_ = 0;
0022   nextDiskValue_ = 0;
0023   nextLayerDiskValue_ = 0;
0024   lastMultiplicity_ = (1 << (3 * alphaBits_));
0025 
0026   for (int i = 0; i < (1 << Nlay_); i++) {
0027     LayerMem_.push_back(-1);
0028   }
0029 
0030   for (int i = 0; i < (1 << (2 * Ndisk_)); i++) {
0031     DiskMem_.push_back(-1);
0032   }
0033 
0034   for (int i = 0; i < (1 << (LayerMemBits_ + DiskMemBits_)); i++) {
0035     LayerDiskMem_.push_back(-1);
0036   }
0037 }
0038 
0039 const TrackDer* TrackDerTable::getDerivatives(unsigned int layermask,
0040                                               unsigned int diskmask,
0041                                               unsigned int alphaindex,
0042                                               unsigned int rinvindex) const {
0043   int index = getIndex(layermask, diskmask);
0044   if (index < 0) {
0045     return nullptr;
0046   }
0047   return &derivatives_[index + alphaindex * (1 << settings_.nrinvBitsTable()) + rinvindex];
0048 }
0049 
0050 int TrackDerTable::getIndex(unsigned int layermask, unsigned int diskmask) const {
0051   assert(layermask < LayerMem_.size());
0052 
0053   assert(diskmask < DiskMem_.size());
0054 
0055   int layercode = LayerMem_[layermask];
0056   int diskcode = DiskMem_[diskmask];
0057 
0058   if (diskcode < 0 || layercode < 0) {
0059     if (settings_.warnNoDer()) {
0060       edm::LogPrint("Tracklet") << "layermask diskmask : " << layermask << " " << diskmask;
0061     }
0062     return -1;
0063   }
0064 
0065   assert(layercode >= 0);
0066   assert(layercode < (1 << LayerMemBits_));
0067   assert(diskcode >= 0);
0068   assert(diskcode < (1 << DiskMemBits_));
0069 
0070   int layerdiskaddress = layercode + (diskcode << LayerMemBits_);
0071 
0072   assert(layerdiskaddress >= 0);
0073   assert(layerdiskaddress < (1 << (LayerMemBits_ + DiskMemBits_)));
0074 
0075   int address = LayerDiskMem_[layerdiskaddress];
0076 
0077   if (address < 0) {
0078     if (settings_.warnNoDer()) {
0079       edm::LogVerbatim("Tracklet") << "layermask diskmask : " << layermask << " " << diskmask;
0080     }
0081     return -1;
0082   }
0083 
0084   assert(address >= 0);
0085   assert(address < (1 << LayerDiskMemBits_));
0086 
0087   return address;
0088 }
0089 
0090 void TrackDerTable::addEntry(unsigned int layermask, unsigned int diskmask, int multiplicity, int nrinv) {
0091   assert(multiplicity <= (1 << (3 * alphaBits_)));
0092 
0093   assert(layermask < (unsigned int)(1 << Nlay_));
0094 
0095   assert(diskmask < (unsigned int)(1 << (2 * Ndisk_)));
0096 
0097   if (LayerMem_[layermask] == -1) {
0098     LayerMem_[layermask] = nextLayerValue_++;
0099   }
0100   if (DiskMem_[diskmask] == -1) {
0101     DiskMem_[diskmask] = nextDiskValue_++;
0102   }
0103 
0104   int layercode = LayerMem_[layermask];
0105   int diskcode = DiskMem_[diskmask];
0106 
0107   assert(layercode >= 0);
0108   assert(layercode < (1 << LayerMemBits_));
0109   assert(diskcode >= 0);
0110   assert(diskcode < (1 << DiskMemBits_));
0111 
0112   int layerdiskaddress = layercode + (diskcode << LayerMemBits_);
0113 
0114   assert(layerdiskaddress >= 0);
0115   assert(layerdiskaddress < (1 << (LayerMemBits_ + DiskMemBits_)));
0116 
0117   int address = LayerDiskMem_[layerdiskaddress];
0118 
0119   if (address != -1) {
0120     edm::LogPrint("Tracklet") << "Duplicate entry:  layermask=" << layermask << " diskmaks=" << diskmask;
0121   }
0122 
0123   assert(address == -1);
0124 
0125   LayerDiskMem_[layerdiskaddress] = nextLayerDiskValue_;
0126 
0127   nextLayerDiskValue_ += multiplicity * nrinv;
0128 
0129   lastMultiplicity_ = multiplicity * nrinv;
0130 
0131   for (int i = 0; i < multiplicity; i++) {
0132     for (int irinv = 0; irinv < nrinv; irinv++) {
0133       TrackDer tmp;
0134       tmp.setIndex(layermask, diskmask, i, irinv);
0135       derivatives_.push_back(tmp);
0136     }
0137   }
0138 }
0139 
0140 void TrackDerTable::readPatternFile(std::string fileName) {
0141   ifstream in(fileName.c_str());
0142   if (settings_.debugTracklet()) {
0143     edm::LogVerbatim("Tracklet") << "reading fit pattern file " << fileName;
0144     edm::LogVerbatim("Tracklet") << "  flags (good/eof/fail/bad): " << in.good() << " " << in.eof() << " " << in.fail()
0145                                  << " " << in.bad();
0146   }
0147 
0148   while (in.good()) {
0149     std::string layerstr, diskstr;
0150     int multiplicity;
0151 
0152     in >> layerstr >> diskstr >> multiplicity;
0153 
0154     //correct multiplicity if you dont want 3 bits of alpha.
0155     if (alphaBits_ == 2) {
0156       if (multiplicity == 8)
0157         multiplicity = 4;
0158       if (multiplicity == 64)
0159         multiplicity = 16;
0160       if (multiplicity == 512)
0161         multiplicity = 64;
0162     }
0163 
0164     if (alphaBits_ == 1) {
0165       if (multiplicity == 8)
0166         multiplicity = 2;
0167       if (multiplicity == 64)
0168         multiplicity = 4;
0169       if (multiplicity == 512)
0170         multiplicity = 8;
0171     }
0172 
0173     if (!in.good())
0174       continue;
0175 
0176     char** tmpptr = nullptr;
0177 
0178     int layers = strtol(layerstr.c_str(), tmpptr, 2);
0179     int disks = strtol(diskstr.c_str(), tmpptr, 2);
0180 
0181     addEntry(layers, disks, multiplicity, (1 << settings_.nrinvBitsTable()));
0182   }
0183 }
0184 
0185 void TrackDerTable::fillTable() {
0186   int nentries = getEntries();
0187 
0188   for (int i = 0; i < nentries; i++) {
0189     TrackDer& der = derivatives_[i];
0190     int layermask = der.layerMask();
0191     int diskmask = der.diskMask();
0192     int alphamask = der.alphaMask();
0193     int irinv = der.irinv();
0194 
0195     double rinv = (irinv - ((1 << (settings_.nrinvBitsTable() - 1)) - 0.5)) * settings_.rinvmax() /
0196                   (1 << (settings_.nrinvBitsTable() - 1));
0197 
0198     bool print = false;
0199 
0200     if (print) {
0201       edm::LogVerbatim("Tracklet") << "PRINT i " << i << " " << layermask << " " << diskmask << " " << alphamask << " "
0202                                    << print;
0203     }
0204 
0205     int nlayers = 0;
0206     double r[N_LAYER];
0207 
0208     for (unsigned l = 0; l < N_LAYER; l++) {
0209       if (layermask & (1 << (N_LAYER - 1 - l))) {
0210         r[nlayers] = settings_.rmean(l);
0211         nlayers++;
0212       }
0213     }
0214 
0215     int ndisks = 0;
0216     double z[N_DISK];
0217     double alpha[N_DISK];
0218 
0219     double t = tpar(settings_, diskmask, layermask);
0220 
0221     for (unsigned d = 0; d < N_DISK; d++) {
0222       if (diskmask & (3 << (2 * (N_DISK - 1 - d)))) {
0223         z[ndisks] = settings_.zmean(d);
0224         alpha[ndisks] = 0.0;
0225         double r = settings_.zmean(d) / t;
0226         double r2 = r * r;
0227         if (diskmask & (1 << (2 * (N_DISK - 1 - d)))) {
0228           if (alphaBits_ == 3) {
0229             int ialpha = alphamask & 7;
0230             alphamask = alphamask >> 3;
0231             alpha[ndisks] = settings_.half2SmoduleWidth() * (ialpha - 3.5) / 4.0 / r2;
0232             if (print)
0233               edm::LogVerbatim("Tracklet") << "PRINT 3 alpha ialpha : " << alpha[ndisks] << " " << ialpha;
0234           }
0235           if (alphaBits_ == 2) {
0236             int ialpha = alphamask & 3;
0237             alphamask = alphamask >> 2;
0238             alpha[ndisks] = settings_.half2SmoduleWidth() * (ialpha - 1.5) / 2.0 / r2;
0239           }
0240           if (alphaBits_ == 1) {
0241             int ialpha = alphamask & 1;
0242             alphamask = alphamask >> 1;
0243             alpha[ndisks] = settings_.half2SmoduleWidth() * (ialpha - 0.5) / r2;
0244             if (print)
0245               edm::LogVerbatim("Tracklet") << "PRINT 1 alpha ialpha : " << alpha[ndisks] << " " << ialpha;
0246           }
0247         }
0248         ndisks++;
0249       }
0250     }
0251 
0252     double D[N_FITPARAM][N_FITSTUB * 2];
0253     int iD[N_FITPARAM][N_FITSTUB * 2];
0254     double MinvDt[N_FITPARAM][N_FITSTUB * 2];
0255     double MinvDtDelta[N_FITPARAM][N_FITSTUB * 2];
0256     int iMinvDt[N_FITPARAM][N_FITSTUB * 2];
0257     double sigma[N_FITSTUB * 2];
0258     double kfactor[N_FITSTUB * 2];
0259 
0260     if (print) {
0261       edm::LogVerbatim("Tracklet") << "PRINT ndisks alpha[0] z[0] t: " << ndisks << " " << alpha[0] << " " << z[0]
0262                                    << " " << t;
0263       for (int iii = 0; iii < nlayers; iii++) {
0264         edm::LogVerbatim("Tracklet") << "PRINT iii r: " << iii << " " << r[iii];
0265       }
0266     }
0267 
0268     calculateDerivatives(settings_, nlayers, r, ndisks, z, alpha, t, rinv, D, iD, MinvDt, iMinvDt, sigma, kfactor);
0269 
0270     double delta = 0.1;
0271 
0272     for (int i = 0; i < nlayers; i++) {
0273       if (r[i] > settings_.rPS2S())
0274         continue;
0275 
0276       r[i] += delta;
0277 
0278       calculateDerivatives(
0279           settings_, nlayers, r, ndisks, z, alpha, t, rinv, D, iD, MinvDtDelta, iMinvDt, sigma, kfactor);
0280 
0281       for (int ii = 0; ii < nlayers; ii++) {
0282         if (r[ii] > settings_.rPS2S())
0283           continue;
0284         double tder = (MinvDtDelta[2][2 * ii + 1] - MinvDt[2][2 * ii + 1]) / delta;
0285         int itder = (1 << (settings_.fittbitshift() + settings_.rcorrbits())) * tder * settings_.kr() * settings_.kz() /
0286                     settings_.ktpars();
0287         double zder = (MinvDtDelta[3][2 * ii + 1] - MinvDt[3][2 * ii + 1]) / delta;
0288         int izder = (1 << (settings_.fitz0bitshift() + settings_.rcorrbits())) * zder * settings_.kr() *
0289                     settings_.kz() / settings_.kz0pars();
0290         der.settdzcorr(i, ii, tder);
0291         der.setz0dzcorr(i, ii, zder);
0292         der.setitdzcorr(i, ii, itder);
0293         der.setiz0dzcorr(i, ii, izder);
0294       }
0295 
0296       r[i] -= delta;
0297     }
0298 
0299     if (print) {
0300       edm::LogVerbatim("Tracklet") << "iMinvDt table build : " << iMinvDt[0][10] << " " << iMinvDt[1][10] << " "
0301                                    << iMinvDt[2][10] << " " << iMinvDt[3][10] << " " << t << " " << nlayers << " "
0302                                    << ndisks;
0303 
0304       std::string oss = "alpha :";
0305       for (int iii = 0; iii < ndisks; iii++) {
0306         oss += " ";
0307         oss += std::to_string(alpha[iii]);
0308       }
0309       edm::LogVerbatim("Tracklet") << oss;
0310       oss = "z :";
0311       for (int iii = 0; iii < ndisks; iii++) {
0312         oss += " ";
0313         oss += std::to_string(z[iii]);
0314       }
0315       edm::LogVerbatim("Tracklet") << oss;
0316     }
0317 
0318     if (print) {
0319       edm::LogVerbatim("Tracklet") << "PRINT nlayers ndisks : " << nlayers << " " << ndisks;
0320     }
0321 
0322     for (int j = 0; j < nlayers + ndisks; j++) {
0323       der.settpar(t);
0324 
0325       //integer
0326       assert(std::abs(iMinvDt[0][2 * j]) < (1 << 23));
0327       assert(std::abs(iMinvDt[0][2 * j + 1]) < (1 << 23));
0328       assert(std::abs(iMinvDt[1][2 * j]) < (1 << 23));
0329       assert(std::abs(iMinvDt[1][2 * j + 1]) < (1 << 23));
0330       assert(std::abs(iMinvDt[2][2 * j]) < (1 << 19));
0331       assert(std::abs(iMinvDt[2][2 * j + 1]) < (1 << 19));
0332       assert(std::abs(iMinvDt[3][2 * j]) < (1 << 19));
0333       assert(std::abs(iMinvDt[3][2 * j + 1]) < (1 << 19));
0334 
0335       if (print) {
0336         edm::LogVerbatim("Tracklet") << "PRINT i " << i << " " << j << " " << iMinvDt[1][2 * j] << " "
0337                                      << std::abs(iMinvDt[1][2 * j]);
0338       }
0339 
0340       der.setirinvdphi(j, iMinvDt[0][2 * j]);
0341       der.setirinvdzordr(j, iMinvDt[0][2 * j + 1]);
0342       der.setiphi0dphi(j, iMinvDt[1][2 * j]);
0343       der.setiphi0dzordr(j, iMinvDt[1][2 * j + 1]);
0344       der.setitdphi(j, iMinvDt[2][2 * j]);
0345       der.setitdzordr(j, iMinvDt[2][2 * j + 1]);
0346       der.setiz0dphi(j, iMinvDt[3][2 * j]);
0347       der.setiz0dzordr(j, iMinvDt[3][2 * j + 1]);
0348       //floating point
0349       der.setrinvdphi(j, MinvDt[0][2 * j]);
0350       der.setrinvdzordr(j, MinvDt[0][2 * j + 1]);
0351       der.setphi0dphi(j, MinvDt[1][2 * j]);
0352       der.setphi0dzordr(j, MinvDt[1][2 * j + 1]);
0353       der.settdphi(j, MinvDt[2][2 * j]);
0354       der.settdzordr(j, MinvDt[2][2 * j + 1]);
0355       der.setz0dphi(j, MinvDt[3][2 * j]);
0356       der.setz0dzordr(j, MinvDt[3][2 * j + 1]);
0357     }
0358   }
0359 
0360   if (settings_.writeTable()) {
0361     ofstream outL = openfile(settings_.tablePath(), "FitDerTableNew_LayerMem.tab", __FILE__, __LINE__);
0362 
0363     int nbits = 6;
0364     for (unsigned int i = 0; i < LayerMem_.size(); i++) {
0365       FPGAWord tmp;
0366       int tmp1 = LayerMem_[i];
0367       if (tmp1 < 0)
0368         tmp1 = (1 << nbits) - 1;
0369       tmp.set(tmp1, nbits, true, __LINE__, __FILE__);
0370       outL << tmp.str() << endl;
0371     }
0372     outL.close();
0373 
0374     ofstream outD = openfile(settings_.tablePath(), "FitDerTableNew_DiskMem.tab", __FILE__, __LINE__);
0375 
0376     nbits = 7;
0377     for (int tmp1 : DiskMem_) {
0378       if (tmp1 < 0)
0379         tmp1 = (1 << nbits) - 1;
0380       FPGAWord tmp;
0381       tmp.set(tmp1, nbits, true, __LINE__, __FILE__);
0382       outD << tmp.str() << endl;
0383     }
0384     outD.close();
0385 
0386     ofstream outLD = openfile(settings_.tablePath(), "FitDerTableNew_LayerDiskMem.tab", __FILE__, __LINE__);
0387 
0388     nbits = 15;
0389     for (int tmp1 : LayerDiskMem_) {
0390       if (tmp1 < 0)
0391         tmp1 = (1 << nbits) - 1;
0392       FPGAWord tmp;
0393       tmp.set(tmp1, nbits, true, __LINE__, __FILE__);
0394       outLD << tmp.str() << endl;
0395     }
0396     outLD.close();
0397 
0398     const std::array<string, N_TRKLSEED> seedings = {{"L1L2", "L3L4", "L5L6", "D1D2", "D3D4", "D1L1", "D1L2"}};
0399     const string prefix = settings_.tablePath() + "FitDerTableNew_";
0400 
0401     // open files for derivative tables
0402 
0403     ofstream outrinvdphi[N_TRKLSEED];
0404     for (unsigned int i = 0; i < N_TRKLSEED; ++i) {
0405       const string fname = prefix + "Rinvdphi_" + seedings[i] + ".tab";
0406       outrinvdphi[i].open(fname);
0407       if (outrinvdphi[i].fail())
0408         throw cms::Exception("BadFile") << __FILE__ << " " << __LINE__ << " could not create file " << fname;
0409     }
0410 
0411     ofstream outrinvdzordr[N_TRKLSEED];
0412     for (unsigned int i = 0; i < N_TRKLSEED; ++i) {
0413       const string fname = prefix + "Rinvdzordr_" + seedings[i] + ".tab";
0414       outrinvdzordr[i].open(fname);
0415       if (outrinvdzordr[i].fail())
0416         throw cms::Exception("BadFile") << __FILE__ << " " << __LINE__ << " could not create file " << fname;
0417     }
0418 
0419     ofstream outphi0dphi[N_TRKLSEED];
0420     for (unsigned int i = 0; i < N_TRKLSEED; ++i) {
0421       const string fname = prefix + "Phi0dphi_" + seedings[i] + ".tab";
0422       outphi0dphi[i].open(fname);
0423       if (outphi0dphi[i].fail())
0424         throw cms::Exception("BadFile") << __FILE__ << " " << __LINE__ << " could not create file " << fname;
0425     }
0426 
0427     ofstream outphi0dzordr[N_TRKLSEED];
0428     for (unsigned int i = 0; i < N_TRKLSEED; ++i) {
0429       const string fname = prefix + "Phi0dzordr_" + seedings[i] + ".tab";
0430       outphi0dzordr[i].open(fname);
0431       if (outphi0dzordr[i].fail())
0432         throw cms::Exception("BadFile") << __FILE__ << " " << __LINE__ << " could not create file " << fname;
0433     }
0434 
0435     ofstream outtdphi[N_TRKLSEED];
0436     for (unsigned int i = 0; i < N_TRKLSEED; ++i) {
0437       const string fname = prefix + "Tdphi_" + seedings[i] + ".tab";
0438       outtdphi[i].open(fname);
0439       if (outtdphi[i].fail())
0440         throw cms::Exception("BadFile") << __FILE__ << " " << __LINE__ << " could not create file " << fname;
0441     }
0442 
0443     ofstream outtdzordr[N_TRKLSEED];
0444     for (unsigned int i = 0; i < N_TRKLSEED; ++i) {
0445       const string fname = prefix + "Tdzordr_" + seedings[i] + ".tab";
0446       outtdzordr[i].open(fname);
0447       if (outtdzordr[i].fail())
0448         throw cms::Exception("BadFile") << __FILE__ << " " << __LINE__ << " could not create file " << fname;
0449     }
0450 
0451     ofstream outz0dphi[N_TRKLSEED];
0452     for (unsigned int i = 0; i < N_TRKLSEED; ++i) {
0453       const string fname = prefix + "Z0dphi_" + seedings[i] + ".tab";
0454       outz0dphi[i].open(fname);
0455       if (outz0dphi[i].fail())
0456         throw cms::Exception("BadFile") << __FILE__ << " " << __LINE__ << " could not create file " << fname;
0457     }
0458 
0459     ofstream outz0dzordr[N_TRKLSEED];
0460     for (unsigned int i = 0; i < N_TRKLSEED; ++i) {
0461       string fname = prefix + "Z0dzordr_" + seedings[i] + ".tab";
0462       outz0dzordr[i].open(fname);
0463       if (outz0dzordr[i].fail())
0464         throw cms::Exception("BadFile") << __FILE__ << " " << __LINE__ << " could not create file " << fname;
0465     }
0466 
0467     for (auto& der : derivatives_) {
0468       unsigned int layerhits = der.layerMask();  // 6 bits layer hit pattern
0469       unsigned int diskmask = der.diskMask();    // 10 bits disk hit pattern
0470       unsigned int diskhits = 0;
0471       if (diskmask & (3 << 8))
0472         diskhits += 16;
0473       if (diskmask & (3 << 6))
0474         diskhits += 8;
0475       if (diskmask & (3 << 4))
0476         diskhits += 4;
0477       if (diskmask & (3 << 2))
0478         diskhits += 2;
0479       if (diskmask & (3 << 0))
0480         diskhits += 1;
0481       assert(diskhits < 32);                            // 5 bits
0482       unsigned int hits = (layerhits << 5) + diskhits;  // 11 bits hit pattern
0483       assert(hits < 4096);
0484 
0485       // loop over all seedings
0486       int i = 0;  // seeding index
0487       for (const string& seed : seedings) {
0488         unsigned int iseed1 = 0;
0489         unsigned int iseed2 = 0;
0490         // check if the seeding is good for the current hit pattern
0491         if (seed == "L1L2") {
0492           iseed1 = 1;
0493           iseed2 = 2;
0494         }
0495         if (seed == "L3L4") {
0496           iseed1 = 3;
0497           iseed2 = 4;
0498         }
0499         if (seed == "L5L6") {
0500           iseed1 = 5;
0501           iseed2 = 6;
0502         }
0503         if (seed == "D1D2") {
0504           iseed1 = 7;
0505           iseed2 = 8;
0506         }
0507         if (seed == "D3D4") {
0508           iseed1 = 9;
0509           iseed2 = 10;
0510         }
0511         if (seed == "D1L1") {
0512           iseed1 = 7;
0513           iseed2 = 1;
0514         }
0515         if (seed == "D1L2") {
0516           iseed1 = 7;
0517           iseed2 = 2;
0518         }
0519 
0520         bool goodseed = (hits & (1 << (11 - iseed1))) and (hits & (1 << (11 - iseed2)));
0521 
0522         int itmprinvdphi[N_PROJ] = {9999999, 9999999, 9999999, 9999999};
0523         int itmprinvdzordr[N_PROJ] = {9999999, 9999999, 9999999, 9999999};
0524         int itmpphi0dphi[N_PROJ] = {9999999, 9999999, 9999999, 9999999};
0525         int itmpphi0dzordr[N_PROJ] = {9999999, 9999999, 9999999, 9999999};
0526         int itmptdphi[N_PROJ] = {9999999, 9999999, 9999999, 9999999};
0527         int itmptdzordr[N_PROJ] = {9999999, 9999999, 9999999, 9999999};
0528         int itmpz0dphi[N_PROJ] = {9999999, 9999999, 9999999, 9999999};
0529         int itmpz0dzordr[N_PROJ] = {9999999, 9999999, 9999999, 9999999};
0530 
0531         // loop over bits in hit pattern
0532         int ider = 0;
0533         if (goodseed) {
0534           for (unsigned int ihit = 1; ihit < N_FITSTUB * 2; ++ihit) {
0535             // skip seeding layers
0536             if (ihit == iseed1 or ihit == iseed2) {
0537               ider++;
0538               continue;
0539             }
0540             // skip if no hit
0541             if (not(hits & (1 << (11 - ihit))))
0542               continue;
0543 
0544             int inputI = -1;
0545             if (seed == "L1L2") {
0546               if (ihit == 3 or ihit == 10)
0547                 inputI = 0;  // L3 or D4
0548               if (ihit == 4 or ihit == 9)
0549                 inputI = 1;  // L4 or D3
0550               if (ihit == 5 or ihit == 8)
0551                 inputI = 2;  // L5 or D2
0552               if (ihit == 6 or ihit == 7)
0553                 inputI = 3;  // L6 or D1
0554             } else if (seed == "L3L4") {
0555               if (ihit == 1)
0556                 inputI = 0;  // L1
0557               if (ihit == 2)
0558                 inputI = 1;  // L2
0559               if (ihit == 5 or ihit == 8)
0560                 inputI = 2;  // L5 or D2
0561               if (ihit == 6 or ihit == 7)
0562                 inputI = 3;  // L6 or D1
0563             } else if (seed == "L5L6") {
0564               if (ihit == 1)
0565                 inputI = 0;  // L1
0566               if (ihit == 2)
0567                 inputI = 1;  // L2
0568               if (ihit == 3)
0569                 inputI = 2;  // L3
0570               if (ihit == 4)
0571                 inputI = 3;  // L4
0572             } else if (seed == "D1D2") {
0573               if (ihit == 1)
0574                 inputI = 0;  // L1
0575               if (ihit == 9)
0576                 inputI = 1;  // D3
0577               if (ihit == 10)
0578                 inputI = 2;  // D4
0579               if (ihit == 2 or ihit == 11)
0580                 inputI = 3;  // L2 or D5
0581             } else if (seed == "D3D4") {
0582               if (ihit == 1)
0583                 inputI = 0;  // L1
0584               if (ihit == 7)
0585                 inputI = 1;  // D1
0586               if (ihit == 8)
0587                 inputI = 2;  // D2
0588               if (ihit == 2 or ihit == 11)
0589                 inputI = 3;  // L2 or D5
0590             } else if (seed == "D1L1" or "D1L2") {
0591               if (ihit == 8)
0592                 inputI = 0;  // D2
0593               if (ihit == 9)
0594                 inputI = 1;  // D3
0595               if (ihit == 10)
0596                 inputI = 2;  // D4
0597               if (ihit == 11)
0598                 inputI = 3;  // D5
0599             }
0600             if (inputI >= 0 and inputI < (int)N_PROJ) {
0601               itmprinvdphi[inputI] = der.irinvdphi(ider);
0602               itmprinvdzordr[inputI] = der.irinvdzordr(ider);
0603               itmpphi0dphi[inputI] = der.iphi0dphi(ider);
0604               itmpphi0dzordr[inputI] = der.iphi0dzordr(ider);
0605               itmptdphi[inputI] = der.itdphi(ider);
0606               itmptdzordr[inputI] = der.itdzordr(ider);
0607               itmpz0dphi[inputI] = der.iz0dphi(ider);
0608               itmpz0dzordr[inputI] = der.iz0dzordr(ider);
0609             }
0610 
0611             ider++;
0612 
0613           }  // for (unsigned int ihit = 1; ihit < 12; ++ihit)
0614         }    // if (goodseed)
0615 
0616         FPGAWord tmprinvdphi[N_PROJ];
0617         int nbits = 16;
0618         for (unsigned int j = 0; j < N_PROJ; ++j) {
0619           if (itmprinvdphi[j] > (1 << nbits))
0620             itmprinvdphi[j] = (1 << nbits) - 1;
0621           tmprinvdphi[j].set(itmprinvdphi[j], nbits + 1, false, __LINE__, __FILE__);
0622         }
0623         outrinvdphi[i] << tmprinvdphi[0].str() << tmprinvdphi[1].str() << tmprinvdphi[2].str() << tmprinvdphi[3].str()
0624                        << endl;
0625 
0626         FPGAWord tmprinvdzordr[N_PROJ];
0627         nbits = 15;
0628         for (unsigned int j = 0; j < N_PROJ; ++j) {
0629           if (itmprinvdzordr[j] > (1 << nbits))
0630             itmprinvdzordr[j] = (1 << nbits) - 1;
0631           tmprinvdzordr[j].set(itmprinvdzordr[j], nbits + 1, false, __LINE__, __FILE__);
0632         }
0633         outrinvdzordr[i] << tmprinvdzordr[0].str() << tmprinvdzordr[1].str() << tmprinvdzordr[2].str()
0634                          << tmprinvdzordr[3].str() << endl;
0635 
0636         FPGAWord tmpphi0dphi[N_PROJ];
0637         nbits = 13;
0638         for (unsigned int j = 0; j < N_PROJ; ++j) {
0639           if (itmpphi0dphi[j] > (1 << nbits))
0640             itmpphi0dphi[j] = (1 << nbits) - 1;
0641           tmpphi0dphi[j].set(itmpphi0dphi[j], nbits + 1, false, __LINE__, __FILE__);
0642         }
0643         outphi0dphi[i] << tmpphi0dphi[0].str() << tmpphi0dphi[1].str() << tmpphi0dphi[2].str() << tmpphi0dphi[3].str()
0644                        << endl;
0645 
0646         FPGAWord tmpphi0dzordr[N_PROJ];
0647         nbits = 15;
0648         for (unsigned int j = 0; j < N_PROJ; ++j) {
0649           if (itmpphi0dzordr[j] > (1 << nbits))
0650             itmpphi0dzordr[j] = (1 << nbits) - 1;
0651           tmpphi0dzordr[j].set(itmpphi0dzordr[j], nbits + 1, false, __LINE__, __FILE__);
0652         }
0653         outphi0dzordr[i] << tmpphi0dzordr[0].str() << tmpphi0dzordr[1].str() << tmpphi0dzordr[2].str()
0654                          << tmpphi0dzordr[3].str() << endl;
0655 
0656         FPGAWord tmptdphi[N_PROJ];
0657         nbits = 14;
0658         for (unsigned int j = 0; j < N_PROJ; ++j) {
0659           if (itmptdphi[j] > (1 << nbits))
0660             itmptdphi[j] = (1 << nbits) - 1;
0661           tmptdphi[j].set(itmptdphi[j], nbits + 1, false, __LINE__, __FILE__);
0662         }
0663         outtdphi[i] << tmptdphi[0].str() << tmptdphi[1].str() << tmptdphi[2].str() << tmptdphi[3].str() << endl;
0664 
0665         FPGAWord tmptdzordr[N_PROJ];
0666         nbits = 15;
0667         for (unsigned int j = 0; j < N_PROJ; ++j) {
0668           if (itmptdzordr[j] > (1 << nbits))
0669             itmptdzordr[j] = (1 << nbits) - 1;
0670           tmptdzordr[j].set(itmptdzordr[j], nbits + 1, false, __LINE__, __FILE__);
0671         }
0672         outtdzordr[i] << tmptdzordr[0].str() << tmptdzordr[1].str() << tmptdzordr[2].str() << tmptdzordr[3].str()
0673                       << endl;
0674 
0675         FPGAWord tmpz0dphi[N_PROJ];
0676         nbits = 13;
0677         for (unsigned int j = 0; j < N_PROJ; ++j) {
0678           if (itmpz0dphi[j] > (1 << nbits))
0679             itmpz0dphi[j] = (1 << nbits) - 1;
0680           tmpz0dphi[j].set(itmpz0dphi[j], nbits + 1, false, __LINE__, __FILE__);
0681         }
0682         outz0dphi[i] << tmpz0dphi[0].str() << tmpz0dphi[1].str() << tmpz0dphi[2].str() << tmpz0dphi[3].str() << endl;
0683 
0684         FPGAWord tmpz0dzordr[N_PROJ];
0685         nbits = 15;
0686         for (unsigned int j = 0; j < N_PROJ; ++j) {
0687           if (itmpz0dzordr[j] > (1 << nbits))
0688             itmpz0dzordr[j] = (1 << nbits) - 1;
0689           tmpz0dzordr[j].set(itmpz0dzordr[j], nbits + 1, false, __LINE__, __FILE__);
0690         }
0691         outz0dzordr[i] << tmpz0dzordr[0].str() << tmpz0dzordr[1].str() << tmpz0dzordr[2].str() << tmpz0dzordr[3].str()
0692                        << endl;
0693 
0694         i++;
0695       }  // for (const string & seed : seedings)
0696 
0697     }  // for (auto & der : derivatives_)
0698 
0699     // close files
0700     for (unsigned int i = 0; i < N_TRKLSEED; ++i) {
0701       outrinvdphi[i].close();
0702       outrinvdzordr[i].close();
0703       outphi0dphi[i].close();
0704       outphi0dzordr[i].close();
0705       outtdphi[i].close();
0706       outtdzordr[i].close();
0707       outz0dphi[i].close();
0708       outz0dzordr[i].close();
0709     }
0710 
0711   }  // if (writeFitDerTable)
0712 }
0713 
0714 void TrackDerTable::invert(double M[4][8], unsigned int n) {
0715   assert(n <= 4);
0716 
0717   unsigned int i, j, k;
0718   double ratio, a;
0719 
0720   for (i = 0; i < n; i++) {
0721     for (j = n; j < 2 * n; j++) {
0722       if (i == (j - n))
0723         M[i][j] = 1.0;
0724       else
0725         M[i][j] = 0.0;
0726     }
0727   }
0728 
0729   for (i = 0; i < n; i++) {
0730     for (j = 0; j < n; j++) {
0731       if (i != j) {
0732         ratio = M[j][i] / M[i][i];
0733         for (k = 0; k < 2 * n; k++) {
0734           M[j][k] -= ratio * M[i][k];
0735         }
0736       }
0737     }
0738   }
0739 
0740   for (i = 0; i < n; i++) {
0741     a = M[i][i];
0742     for (j = 0; j < 2 * n; j++) {
0743       M[i][j] /= a;
0744     }
0745   }
0746 }
0747 
0748 void TrackDerTable::invert(std::vector<std::vector<double> >& M, unsigned int n) {
0749   assert(M.size() == n);
0750   assert(M[0].size() == 2 * n);
0751 
0752   unsigned int i, j, k;
0753   double ratio, a;
0754 
0755   for (i = 0; i < n; i++) {
0756     for (j = n; j < 2 * n; j++) {
0757       if (i == (j - n))
0758         M[i][j] = 1.0;
0759       else
0760         M[i][j] = 0.0;
0761     }
0762   }
0763 
0764   for (i = 0; i < n; i++) {
0765     for (j = 0; j < n; j++) {
0766       if (i != j) {
0767         ratio = M[j][i] / M[i][i];
0768         for (k = 0; k < 2 * n; k++) {
0769           M[j][k] -= ratio * M[i][k];
0770         }
0771       }
0772     }
0773   }
0774 
0775   for (i = 0; i < n; i++) {
0776     a = M[i][i];
0777     for (j = 0; j < 2 * n; j++) {
0778       M[i][j] /= a;
0779     }
0780   }
0781 }
0782 
0783 void TrackDerTable::calculateDerivatives(Settings const& settings,
0784                                          unsigned int nlayers,
0785                                          double r[N_LAYER],
0786                                          unsigned int ndisks,
0787                                          double z[N_DISK],
0788                                          double alpha[N_DISK],
0789                                          double t,
0790                                          double rinv,
0791                                          double D[N_FITPARAM][N_FITSTUB * 2],
0792                                          int iD[N_FITPARAM][N_FITSTUB * 2],
0793                                          double MinvDt[N_FITPARAM][N_FITSTUB * 2],
0794                                          int iMinvDt[N_FITPARAM][N_FITSTUB * 2],
0795                                          double sigma[N_FITSTUB * 2],
0796                                          double kfactor[N_FITSTUB * 2]) {
0797   double sigmax = settings.stripPitch(true) / sqrt(12.0);
0798   double sigmaz = settings.stripLength(true) / sqrt(12.0);
0799   double sigmaz2 = settings.stripLength(false) / sqrt(12.0);
0800 
0801   double sigmazpsbarrel = sigmaz;  //This is a bit of a hack - these weights should be properly determined
0802   if (std::abs(t) > 2.0)
0803     sigmazpsbarrel = sigmaz * std::abs(t) / 2.0;
0804   if (std::abs(t) > 3.8)
0805     sigmazpsbarrel = sigmaz * std::abs(t);
0806 
0807   double sigmax2sdisk = settings.stripPitch(false) / sqrt(12.0);
0808   double sigmaz2sdisk = settings.stripLength(false) / sqrt(12.0);
0809 
0810   double sigmaxpsdisk = settings.stripPitch(true) / sqrt(12.0);
0811   double sigmazpsdisk = settings.stripLength(true) / sqrt(12.0);
0812 
0813   unsigned int n = nlayers + ndisks;
0814 
0815   assert(n <= N_FITSTUB);
0816 
0817   double rnew[N_FITSTUB];
0818 
0819   int j = 0;
0820 
0821   //here we handle a barrel hit
0822   for (unsigned int i = 0; i < nlayers; i++) {
0823     double ri = r[i];
0824 
0825     rnew[i] = ri;
0826 
0827     //first we have the phi position
0828     D[0][j] = -0.5 * ri * ri / sqrt(1 - 0.25 * ri * ri * rinv * rinv) / sigmax;
0829     D[1][j] = ri / sigmax;
0830     D[2][j] = 0.0;
0831     D[3][j] = 0.0;
0832     sigma[j] = sigmax;
0833     kfactor[j] = settings.kphi1();
0834     j++;
0835     //second the z position
0836     D[0][j] = 0.0;
0837     D[1][j] = 0.0;
0838     if (ri < settings.rPS2S()) {
0839       D[2][j] = (2 / rinv) * asin(0.5 * ri * rinv) / sigmazpsbarrel;
0840       D[3][j] = 1.0 / sigmazpsbarrel;
0841       sigma[j] = sigmazpsbarrel;
0842       kfactor[j] = settings.kz();
0843     } else {
0844       D[2][j] = (2 / rinv) * asin(0.5 * ri * rinv) / sigmaz2;
0845       D[3][j] = 1.0 / sigmaz2;
0846       sigma[j] = sigmaz2;
0847       kfactor[j] = settings.kz();
0848     }
0849 
0850     j++;
0851   }
0852 
0853   for (unsigned int i = 0; i < ndisks; i++) {
0854     double zi = z[i];
0855 
0856     double z0 = 0.0;
0857 
0858     double rmultiplier = alpha[i] * zi / t;
0859 
0860     double phimultiplier = zi / t;
0861 
0862     double drdrinv = -2.0 * sin(0.5 * rinv * (zi - z0) / t) / (rinv * rinv) +
0863                      (zi - z0) * cos(0.5 * rinv * (zi - z0) / t) / (rinv * t);
0864     double drdphi0 = 0;
0865     double drdt = -(zi - z0) * cos(0.5 * rinv * (zi - z0) / t) / (t * t);
0866     double drdz0 = -cos(0.5 * rinv * (zi - z0) / t) / t;
0867 
0868     double dphidrinv = -0.5 * (zi - z0) / t;
0869     double dphidphi0 = 1.0;
0870     double dphidt = 0.5 * rinv * (zi - z0) / (t * t);
0871     double dphidz0 = 0.5 * rinv / t;
0872 
0873     double r = (zi - z0) / t;
0874 
0875     rnew[i + nlayers] = r;
0876 
0877     sigma[j] = sigmax2sdisk;
0878     if (std::abs(alpha[i]) < 1e-10) {
0879       sigma[j] = sigmaxpsdisk;
0880     }
0881 
0882     D[0][j] = (phimultiplier * dphidrinv + rmultiplier * drdrinv) / sigma[j];
0883     D[1][j] = (phimultiplier * dphidphi0 + rmultiplier * drdphi0) / sigma[j];
0884     D[2][j] = (phimultiplier * dphidt + rmultiplier * drdt) / sigma[j];
0885     D[3][j] = (phimultiplier * dphidz0 + rmultiplier * drdz0) / sigma[j];
0886     kfactor[j] = settings.kphi();
0887 
0888     j++;
0889 
0890     if (std::abs(alpha[i]) < 1e-10) {
0891       D[0][j] = drdrinv / sigmazpsdisk;
0892       D[1][j] = drdphi0 / sigmazpsdisk;
0893       D[2][j] = drdt / sigmazpsdisk;
0894       D[3][j] = drdz0 / sigmazpsdisk;
0895       sigma[j] = sigmazpsdisk;
0896       kfactor[j] = settings.kr();
0897     } else {
0898       D[0][j] = drdrinv / sigmaz2sdisk;
0899       D[1][j] = drdphi0 / sigmaz2sdisk;
0900       D[2][j] = drdt / sigmaz2sdisk;
0901       D[3][j] = drdz0 / sigmaz2sdisk;
0902       sigma[j] = sigmaz2sdisk;
0903       kfactor[j] = settings.kr();
0904     }
0905 
0906     j++;
0907   }
0908 
0909   double M[4][8];
0910 
0911   for (unsigned int i1 = 0; i1 < 4; i1++) {
0912     for (unsigned int i2 = 0; i2 < 4; i2++) {
0913       M[i1][i2] = 0.0;
0914       for (unsigned int j = 0; j < 2 * n; j++) {
0915         M[i1][i2] += D[i1][j] * D[i2][j];
0916       }
0917     }
0918   }
0919 
0920   invert(M, 4);
0921 
0922   for (unsigned int j = 0; j < N_FITSTUB * 2; j++) {
0923     for (unsigned int i1 = 0; i1 < N_FITPARAM; i1++) {
0924       MinvDt[i1][j] = 0.0;
0925       iMinvDt[i1][j] = 0;
0926     }
0927   }
0928 
0929   for (unsigned int j = 0; j < 2 * n; j++) {
0930     for (unsigned int i1 = 0; i1 < 4; i1++) {
0931       for (unsigned int i2 = 0; i2 < 4; i2++) {
0932         MinvDt[i1][j] += M[i1][i2 + 4] * D[i2][j];
0933       }
0934     }
0935   }
0936 
0937   for (unsigned int i = 0; i < n; i++) {
0938     iD[0][2 * i] =
0939         D[0][2 * i] * (1 << settings.chisqphifactbits()) * settings.krinvpars() / (1 << settings.fitrinvbitshift());
0940     iD[1][2 * i] =
0941         D[1][2 * i] * (1 << settings.chisqphifactbits()) * settings.kphi0pars() / (1 << settings.fitphi0bitshift());
0942     iD[2][2 * i] =
0943         D[2][2 * i] * (1 << settings.chisqphifactbits()) * settings.ktpars() / (1 << settings.fittbitshift());
0944     iD[3][2 * i] =
0945         D[3][2 * i] * (1 << settings.chisqphifactbits()) * settings.kz0pars() / (1 << settings.fitz0bitshift());
0946 
0947     iD[0][2 * i + 1] =
0948         D[0][2 * i + 1] * (1 << settings.chisqzfactbits()) * settings.krinvpars() / (1 << settings.fitrinvbitshift());
0949     iD[1][2 * i + 1] =
0950         D[1][2 * i + 1] * (1 << settings.chisqzfactbits()) * settings.kphi0pars() / (1 << settings.fitphi0bitshift());
0951     iD[2][2 * i + 1] =
0952         D[2][2 * i + 1] * (1 << settings.chisqzfactbits()) * settings.ktpars() / (1 << settings.fittbitshift());
0953     iD[3][2 * i + 1] =
0954         D[3][2 * i + 1] * (1 << settings.chisqzfactbits()) * settings.kz0pars() / (1 << settings.fitz0bitshift());
0955 
0956     //First the barrel
0957     if (i < nlayers) {
0958       MinvDt[0][2 * i] *= rnew[i] / sigmax;
0959       MinvDt[1][2 * i] *= rnew[i] / sigmax;
0960       MinvDt[2][2 * i] *= rnew[i] / sigmax;
0961       MinvDt[3][2 * i] *= rnew[i] / sigmax;
0962 
0963       iMinvDt[0][2 * i] =
0964           (1 << settings.fitrinvbitshift()) * MinvDt[0][2 * i] * settings.kphi1() / settings.krinvpars();
0965       iMinvDt[1][2 * i] =
0966           (1 << settings.fitphi0bitshift()) * MinvDt[1][2 * i] * settings.kphi1() / settings.kphi0pars();
0967       iMinvDt[2][2 * i] = (1 << settings.fittbitshift()) * MinvDt[2][2 * i] * settings.kphi1() / settings.ktpars();
0968       iMinvDt[3][2 * i] = (1 << settings.fitz0bitshift()) * MinvDt[3][2 * i] * settings.kphi1() / settings.kz0pars();
0969 
0970       if (rnew[i] < settings.rPS2S()) {
0971         MinvDt[0][2 * i + 1] /= sigmazpsbarrel;
0972         MinvDt[1][2 * i + 1] /= sigmazpsbarrel;
0973         MinvDt[2][2 * i + 1] /= sigmazpsbarrel;
0974         MinvDt[3][2 * i + 1] /= sigmazpsbarrel;
0975 
0976         iMinvDt[0][2 * i + 1] =
0977             (1 << settings.fitrinvbitshift()) * MinvDt[0][2 * i + 1] * settings.kz() / settings.krinvpars();
0978         iMinvDt[1][2 * i + 1] =
0979             (1 << settings.fitphi0bitshift()) * MinvDt[1][2 * i + 1] * settings.kz() / settings.kphi0pars();
0980         iMinvDt[2][2 * i + 1] =
0981             (1 << settings.fittbitshift()) * MinvDt[2][2 * i + 1] * settings.kz() / settings.ktpars();
0982         iMinvDt[3][2 * i + 1] =
0983             (1 << settings.fitz0bitshift()) * MinvDt[3][2 * i + 1] * settings.kz() / settings.kz0pars();
0984       } else {
0985         MinvDt[0][2 * i + 1] /= sigmaz2;
0986         MinvDt[1][2 * i + 1] /= sigmaz2;
0987         MinvDt[2][2 * i + 1] /= sigmaz2;
0988         MinvDt[3][2 * i + 1] /= sigmaz2;
0989 
0990         int fact = (1 << (settings.nzbitsstub(0) - settings.nzbitsstub(5)));
0991 
0992         iMinvDt[0][2 * i + 1] =
0993             (1 << settings.fitrinvbitshift()) * MinvDt[0][2 * i + 1] * fact * settings.kz() / settings.krinvpars();
0994         iMinvDt[1][2 * i + 1] =
0995             (1 << settings.fitphi0bitshift()) * MinvDt[1][2 * i + 1] * fact * settings.kz() / settings.kphi0pars();
0996         iMinvDt[2][2 * i + 1] =
0997             (1 << settings.fittbitshift()) * MinvDt[2][2 * i + 1] * fact * settings.kz() / settings.ktpars();
0998         iMinvDt[3][2 * i + 1] =
0999             (1 << settings.fitz0bitshift()) * MinvDt[3][2 * i + 1] * fact * settings.kz() / settings.kz0pars();
1000       }
1001     }
1002 
1003     //Secondly the disks
1004     else {
1005       double denom = (std::abs(alpha[i - nlayers]) < 1e-10) ? sigmaxpsdisk : sigmax2sdisk;
1006 
1007       MinvDt[0][2 * i] *= (rnew[i] / denom);
1008       MinvDt[1][2 * i] *= (rnew[i] / denom);
1009       MinvDt[2][2 * i] *= (rnew[i] / denom);
1010       MinvDt[3][2 * i] *= (rnew[i] / denom);
1011 
1012       assert(MinvDt[0][2 * i] == MinvDt[0][2 * i]);
1013 
1014       iMinvDt[0][2 * i] = (1 << settings.fitrinvbitshift()) * MinvDt[0][2 * i] * settings.kphi() / settings.krinvpars();
1015       iMinvDt[1][2 * i] = (1 << settings.fitphi0bitshift()) * MinvDt[1][2 * i] * settings.kphi() / settings.kphi0pars();
1016       iMinvDt[2][2 * i] = (1 << settings.fittbitshift()) * MinvDt[2][2 * i] * settings.kphi() / settings.ktpars();
1017       iMinvDt[3][2 * i] = (1 << settings.fitz0bitshift()) * MinvDt[3][2 * i] * settings.kphi() / settings.kz();
1018 
1019       denom = (std::abs(alpha[i - nlayers]) < 1e-10) ? sigmazpsdisk : sigmaz2sdisk;
1020 
1021       MinvDt[0][2 * i + 1] /= denom;
1022       MinvDt[1][2 * i + 1] /= denom;
1023       MinvDt[2][2 * i + 1] /= denom;
1024       MinvDt[3][2 * i + 1] /= denom;
1025 
1026       iMinvDt[0][2 * i + 1] =
1027           (1 << settings.fitrinvbitshift()) * MinvDt[0][2 * i + 1] * settings.krprojshiftdisk() / settings.krinvpars();
1028       iMinvDt[1][2 * i + 1] =
1029           (1 << settings.fitphi0bitshift()) * MinvDt[1][2 * i + 1] * settings.krprojshiftdisk() / settings.kphi0pars();
1030       iMinvDt[2][2 * i + 1] =
1031           (1 << settings.fittbitshift()) * MinvDt[2][2 * i + 1] * settings.krprojshiftdisk() / settings.ktpars();
1032       iMinvDt[3][2 * i + 1] =
1033           (1 << settings.fitz0bitshift()) * MinvDt[3][2 * i + 1] * settings.krprojshiftdisk() / settings.kz();
1034     }
1035   }
1036 }
1037 
1038 double TrackDerTable::tpar(Settings const& settings, int diskmask, int layermask) {
1039   if (diskmask == 0)
1040     return 0.0;
1041 
1042   double tmax = 1000.0;
1043   double tmin = 0.0;
1044 
1045   for (int d = 1; d <= (int)N_DISK; d++) {
1046     if (diskmask & (1 << (2 * (5 - d) + 1))) {  //PS hit
1047       double dmax = settings.zmean(d - 1) / 22.0;
1048       if (dmax > sinh(2.4))
1049         dmax = sinh(2.4);
1050       double dmin = settings.zmean(d - 1) / 65.0;
1051       if (dmax < tmax)
1052         tmax = dmax;
1053       if (dmin > tmin)
1054         tmin = dmin;
1055     }
1056 
1057     if (diskmask & (1 << (2 * (5 - d)))) {  //2S hit
1058       double dmax = settings.zmean(d - 1) / 65.0;
1059       double dmin = settings.zmean(d - 1) / 105.0;
1060       if (dmax < tmax)
1061         tmax = dmax;
1062       if (dmin > tmin)
1063         tmin = dmin;
1064     }
1065   }
1066 
1067   for (int l = 1; l <= (int)N_LAYER; l++) {
1068     if (layermask & (1 << (6 - l))) {
1069       double lmax = settings.zlength() / settings.rmean(l - 1);
1070       if (lmax < tmax)
1071         tmax = lmax;
1072     }
1073   }
1074 
1075   return 0.5 * (tmax + tmin) * 1.07;
1076 }