File indexing completed on 2024-09-07 04:37:08
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
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
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
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
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();
0469 unsigned int diskmask = der.diskMask();
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);
0482 unsigned int hits = (layerhits << 5) + diskhits;
0483 assert(hits < 4096);
0484
0485
0486 int i = 0;
0487 for (const string& seed : seedings) {
0488 unsigned int iseed1 = 0;
0489 unsigned int iseed2 = 0;
0490
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
0532 int ider = 0;
0533 if (goodseed) {
0534 for (unsigned int ihit = 1; ihit < N_FITSTUB * 2; ++ihit) {
0535
0536 if (ihit == iseed1 or ihit == iseed2) {
0537 ider++;
0538 continue;
0539 }
0540
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;
0548 if (ihit == 4 or ihit == 9)
0549 inputI = 1;
0550 if (ihit == 5 or ihit == 8)
0551 inputI = 2;
0552 if (ihit == 6 or ihit == 7)
0553 inputI = 3;
0554 } else if (seed == "L3L4") {
0555 if (ihit == 1)
0556 inputI = 0;
0557 if (ihit == 2)
0558 inputI = 1;
0559 if (ihit == 5 or ihit == 8)
0560 inputI = 2;
0561 if (ihit == 6 or ihit == 7)
0562 inputI = 3;
0563 } else if (seed == "L5L6") {
0564 if (ihit == 1)
0565 inputI = 0;
0566 if (ihit == 2)
0567 inputI = 1;
0568 if (ihit == 3)
0569 inputI = 2;
0570 if (ihit == 4)
0571 inputI = 3;
0572 } else if (seed == "D1D2") {
0573 if (ihit == 1)
0574 inputI = 0;
0575 if (ihit == 9)
0576 inputI = 1;
0577 if (ihit == 10)
0578 inputI = 2;
0579 if (ihit == 2 or ihit == 11)
0580 inputI = 3;
0581 } else if (seed == "D3D4") {
0582 if (ihit == 1)
0583 inputI = 0;
0584 if (ihit == 7)
0585 inputI = 1;
0586 if (ihit == 8)
0587 inputI = 2;
0588 if (ihit == 2 or ihit == 11)
0589 inputI = 3;
0590 } else if (seed == "D1L1" or "D1L2") {
0591 if (ihit == 8)
0592 inputI = 0;
0593 if (ihit == 9)
0594 inputI = 1;
0595 if (ihit == 10)
0596 inputI = 2;
0597 if (ihit == 11)
0598 inputI = 3;
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 }
0614 }
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 }
0696
0697 }
0698
0699
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 }
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;
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
0822 for (unsigned int i = 0; i < nlayers; i++) {
0823 double ri = r[i];
0824
0825 rnew[i] = ri;
0826
0827
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
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
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
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))) {
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)))) {
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 }