TrackFindingTracklet/src/TrackletCalculatorDisplaced.cc

0001 #include "L1Trigger/TrackFindingTracklet/interface/TrackletCalculatorDisplaced.h"
0002 #include "L1Trigger/TrackFindingTracklet/interface/Settings.h"
0003 #include "L1Trigger/TrackFindingTracklet/interface/Globals.h"
0004 #include "L1Trigger/TrackFindingTracklet/interface/Tracklet.h"
0005 #include "L1Trigger/TrackFindingTracklet/interface/Stub.h"
0006 #include "L1Trigger/TrackFindingTracklet/interface/L1TStub.h"
0007 #include "L1Trigger/TrackFindingTracklet/interface/IMATH_TrackletCalculator.h"
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "FWCore/Utilities/interface/isFinite.h"
0012 #include "DataFormats/Math/interface/deltaPhi.h"
0013
0014 using namespace std;
0015 using namespace trklet;
0016
0017 TrackletCalculatorDisplaced::TrackletCalculatorDisplaced(string name, Settings const& settings, Globals* global)
0018     : ProcessBase(name, settings, global) {
0019   for (unsigned int ilayer = 0; ilayer < N_LAYER; ilayer++) {
0020     vector<TrackletProjectionsMemory*> tmp(settings.nallstubs(ilayer), nullptr);
0021     trackletprojlayers_.push_back(tmp);
0022   }
0023
0024   for (unsigned int idisk = 0; idisk < N_DISK; idisk++) {
0025     vector<TrackletProjectionsMemory*> tmp(settings.nallstubs(idisk + N_LAYER), nullptr);
0026     trackletprojdisks_.push_back(tmp);
0027   }
0028
0029   layer_ = 0;
0030   disk_ = 0;
0031
0032   string name1 = name.substr(1);  //this is to correct for "TCD" having one more letter then "TC"
0033   if (name1[3] == 'L')
0034     layer_ = name1[4] - '0';
0035   if (name1[3] == 'D')
0036     disk_ = name1[4] - '0';
0037
0038   // set TC index
0039   iSeed_ = 0;
0040
0041   int iTC = name1[9] - 'A';
0042
0043   if (name1.substr(3, 6) == "L3L4L2")
0044     iSeed_ = 8;
0045   else if (name1.substr(3, 6) == "L5L6L4")
0046     iSeed_ = 9;
0047   else if (name1.substr(3, 6) == "L2L3D1")
0048     iSeed_ = 10;
0049   else if (name1.substr(3, 6) == "D1D2L2")
0050     iSeed_ = 11;
0051
0052   assert(iSeed_ != 0);
0053
0054   TCIndex_ = (iSeed_ << 4) + iTC;
0055   assert(TCIndex_ >= 128 && TCIndex_ < 191);
0056
0057   assert((layer_ != 0) || (disk_ != 0));
0058
0059   toR_.clear();
0060   toZ_.clear();
0061
0062   if (iSeed_ == 8 || iSeed_ == 9) {
0063     if (layer_ == 3) {
0064       rzmeanInv_[0] = 1.0 / settings_.rmean(2 - 1);
0065       rzmeanInv_[1] = 1.0 / settings_.rmean(3 - 1);
0066       rzmeanInv_[2] = 1.0 / settings_.rmean(4 - 1);
0067
0068       rproj_[0] = settings_.rmean(0);
0069       rproj_[1] = settings_.rmean(4);
0070       rproj_[2] = settings_.rmean(5);
0071       lproj_[0] = 1;
0072       lproj_[1] = 5;
0073       lproj_[2] = 6;
0074
0075       dproj_[0] = 1;
0076       dproj_[1] = 2;
0077       dproj_[2] = 0;
0078       toZ_.push_back(settings_.zmean(0));
0079       toZ_.push_back(settings_.zmean(1));
0080     }
0081     if (layer_ == 5) {
0082       rzmeanInv_[0] = 1.0 / settings_.rmean(4 - 1);
0083       rzmeanInv_[1] = 1.0 / settings_.rmean(5 - 1);
0084       rzmeanInv_[2] = 1.0 / settings_.rmean(6 - 1);
0085
0086       rproj_[0] = settings_.rmean(0);
0087       rproj_[1] = settings_.rmean(1);
0088       rproj_[2] = settings_.rmean(2);
0089       lproj_[0] = 1;
0090       lproj_[1] = 2;
0091       lproj_[2] = 3;
0092
0093       dproj_[0] = 0;
0094       dproj_[1] = 0;
0095       dproj_[2] = 0;
0096     }
0097     for (unsigned int i = 0; i < N_LAYER - 3; ++i)
0098       toR_.push_back(rproj_[i]);
0099   }
0100
0101   if (iSeed_ == 10 || iSeed_ == 11) {
0102     if (layer_ == 2) {
0103       rzmeanInv_[0] = 1.0 / settings_.rmean(2 - 1);
0104       rzmeanInv_[1] = 1.0 / settings_.rmean(3 - 1);
0105       rzmeanInv_[2] = 1.0 / settings_.zmean(1 - 1);
0106
0107       rproj_[0] = settings_.rmean(0);
0108       lproj_[0] = 1;
0109       lproj_[1] = -1;
0110       lproj_[2] = -1;
0111
0112       zproj_[0] = settings_.zmean(1);
0113       zproj_[1] = settings_.zmean(2);
0114       zproj_[2] = settings_.zmean(3);
0115       dproj_[0] = 2;
0116       dproj_[1] = 3;
0117       dproj_[2] = 4;
0118     }
0119     if (disk_ == 1) {
0120       rzmeanInv_[0] = 1.0 / settings_.rmean(2 - 1);
0121       rzmeanInv_[1] = 1.0 / settings_.zmean(1 - 1);
0122       rzmeanInv_[2] = 1.0 / settings_.zmean(2 - 1);
0123
0124       rproj_[0] = settings_.rmean(0);
0125       lproj_[0] = 1;
0126       lproj_[1] = -1;
0127       lproj_[2] = -1;
0128
0129       zproj_[0] = settings_.zmean(2);
0130       zproj_[1] = settings_.zmean(3);
0131       zproj_[2] = settings_.zmean(4);
0132       dproj_[0] = 3;
0133       dproj_[1] = 4;
0134       dproj_[2] = 5;
0135     }
0136     toR_.push_back(settings_.rmean(0));
0137     for (unsigned int i = 0; i < N_DISK - 2; ++i)
0138       toZ_.push_back(zproj_[i]);
0139   }
0140 }
0141
0142 void TrackletCalculatorDisplaced::addOutputProjection(TrackletProjectionsMemory*& outputProj, MemoryBase* memory) {
0143   outputProj = dynamic_cast<TrackletProjectionsMemory*>(memory);
0144   assert(outputProj != nullptr);
0145 }
0146
0147 void TrackletCalculatorDisplaced::addOutput(MemoryBase* memory, string output) {
0148   if (settings_.writetrace()) {
0149     edm::LogVerbatim("Tracklet") << "In " << name_ << " adding output to " << memory->getName() << " to output "
0150                                  << output;
0151   }
0152
0153   if (output == "trackpar") {
0154     auto* tmp = dynamic_cast<TrackletParametersMemory*>(memory);
0155     assert(tmp != nullptr);
0156     trackletpars_ = tmp;
0157     return;
0158   }
0159
0160   if (output.substr(0, 7) == "projout") {
0161     //output is on the form 'projoutL2PHIC' or 'projoutD3PHIB'
0162     auto* tmp = dynamic_cast<TrackletProjectionsMemory*>(memory);
0163     assert(tmp != nullptr);
0164
0165     unsigned int layerdisk = output[8] - '1';   //layer or disk counting from 0
0166     unsigned int phiregion = output[12] - 'A';  //phiregion counting from 0
0167
0168     if (output[7] == 'L') {
0169       assert(layerdisk < N_LAYER);
0170       assert(phiregion < trackletprojlayers_[layerdisk].size());
0171       //check that phiregion not already initialized
0172       assert(trackletprojlayers_[layerdisk][phiregion] == nullptr);
0173       trackletprojlayers_[layerdisk][phiregion] = tmp;
0174       return;
0175     }
0176
0177     if (output[7] == 'D') {
0178       assert(layerdisk < N_DISK);
0179       assert(phiregion < trackletprojdisks_[layerdisk].size());
0180       //check that phiregion not already initialized
0181       assert(trackletprojdisks_[layerdisk][phiregion] == nullptr);
0182       trackletprojdisks_[layerdisk][phiregion] = tmp;
0183       return;
0184     }
0185   }
0186
0187   throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " Could not find output : " << output;
0188 }
0189
0190 void TrackletCalculatorDisplaced::addInput(MemoryBase* memory, string input) {
0191   if (settings_.writetrace()) {
0192     edm::LogVerbatim("Tracklet") << "In " << name_ << " adding input from " << memory->getName() << " to input "
0193                                  << input;
0194   }
0195
0196   if (input == "thirdallstubin") {
0197     auto* tmp = dynamic_cast<AllStubsMemory*>(memory);
0198     assert(tmp != nullptr);
0199     innerallstubs_.push_back(tmp);
0200     return;
0201   }
0202   if (input == "firstallstubin") {
0203     auto* tmp = dynamic_cast<AllStubsMemory*>(memory);
0204     assert(tmp != nullptr);
0205     middleallstubs_.push_back(tmp);
0206     return;
0207   }
0208   if (input == "secondallstubin") {
0209     auto* tmp = dynamic_cast<AllStubsMemory*>(memory);
0210     assert(tmp != nullptr);
0211     outerallstubs_.push_back(tmp);
0212     return;
0213   }
0214   if (input.find("stubtriplet") == 0) {
0215     auto* tmp = dynamic_cast<StubTripletsMemory*>(memory);
0216     assert(tmp != nullptr);
0217     stubtriplets_.push_back(tmp);
0218     return;
0219   }
0220   throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " Could not find input : " << input;
0221 }
0222
0223 void TrackletCalculatorDisplaced::execute(unsigned int iSector, double phimin, double phimax) {
0224   unsigned int countall = 0;
0225   unsigned int countsel = 0;
0226
0227   phimin_ = phimin;
0228   phimax_ = phimax;
0229   iSector_ = iSector;
0230
0231   for (auto& stubtriplet : stubtriplets_) {
0232     if (trackletpars_->nTracklets() >= settings_.ntrackletmax()) {
0233       edm::LogVerbatim("Tracklet") << "Will break on too many tracklets in " << getName();
0234       break;
0235     }
0236     for (unsigned int i = 0; i < stubtriplet->nStubTriplets(); i++) {
0237       countall++;
0238
0239       const Stub* innerFPGAStub = stubtriplet->getFPGAStub1(i);
0240       const L1TStub* innerStub = innerFPGAStub->l1tstub();
0241
0242       const Stub* middleFPGAStub = stubtriplet->getFPGAStub2(i);
0243       const L1TStub* middleStub = middleFPGAStub->l1tstub();
0244
0245       const Stub* outerFPGAStub = stubtriplet->getFPGAStub3(i);
0246       const L1TStub* outerStub = outerFPGAStub->l1tstub();
0247
0248       if (settings_.debugTracklet())
0249         edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced execute " << getName() << "[" << iSector_ << "]";
0250
0251       if (innerFPGAStub->layerdisk() < N_LAYER && middleFPGAStub->layerdisk() < N_LAYER &&
0252           outerFPGAStub->layerdisk() < N_LAYER) {
0253         //barrel+barrel seeding
0254         bool accept = LLLSeeding(innerFPGAStub, innerStub, middleFPGAStub, middleStub, outerFPGAStub, outerStub);
0255         if (accept)
0256           countsel++;
0257       } else if (innerFPGAStub->layerdisk() >= N_LAYER && middleFPGAStub->layerdisk() >= N_LAYER &&
0258                  outerFPGAStub->layerdisk() >= N_LAYER) {
0259         throw cms::Exception("LogicError") << __FILE__ << " " << __LINE__ << " Invalid seeding!";
0260       } else {
0261         //layer+disk seeding
0262         if (innerFPGAStub->layerdisk() < N_LAYER && middleFPGAStub->layerdisk() >= N_LAYER &&
0263             outerFPGAStub->layerdisk() >= N_LAYER) {  //D1D2L2
0264           bool accept = DDLSeeding(innerFPGAStub, innerStub, middleFPGAStub, middleStub, outerFPGAStub, outerStub);
0265           if (accept)
0266             countsel++;
0267         } else if (innerFPGAStub->layerdisk() >= N_LAYER && middleFPGAStub->layerdisk() < N_LAYER &&
0268                    outerFPGAStub->layerdisk() < N_LAYER) {  //L2L3D1
0269           bool accept = LLDSeeding(innerFPGAStub, innerStub, middleFPGAStub, middleStub, outerFPGAStub, outerStub);
0270           if (accept)
0271             countsel++;
0272         } else {
0273           throw cms::Exception("LogicError") << __FILE__ << " " << __LINE__ << " Invalid seeding!";
0274         }
0275       }
0276
0277       if (trackletpars_->nTracklets() >= settings_.ntrackletmax()) {
0278         edm::LogVerbatim("Tracklet") << "Will break on number of tracklets in " << getName();
0279         break;
0280       }
0281
0282       if (countall >= settings_.maxStep("TC")) {
0283         if (settings_.debugTracklet())
0284           edm::LogVerbatim("Tracklet") << "Will break on MAXTC 1";
0285         break;
0286       }
0287       if (settings_.debugTracklet())
0288         edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced execute done";
0289     }
0290     if (countall >= settings_.maxStep("TC")) {
0291       if (settings_.debugTracklet())
0292         edm::LogVerbatim("Tracklet") << "Will break on MAXTC 2";
0293       break;
0294     }
0295   }
0296
0297   if (settings_.writeMonitorData("TPD")) {
0298     globals_->ofstream("trackletcalculatordisplaced.txt") << getName() << " " << countall << " " << countsel << endl;
0299   }
0300 }
0301
0302 void TrackletCalculatorDisplaced::addDiskProj(Tracklet* tracklet, int disk) {
0303   disk = std::abs(disk);
0304   FPGAWord fpgar = tracklet->proj(N_LAYER + disk - 1).fpgarzproj();
0305
0306   if (fpgar.value() * settings_.krprojshiftdisk() < settings_.rmindiskvm())
0307     return;
0308   if (fpgar.value() * settings_.krprojshiftdisk() > settings_.rmaxdisk())
0309     return;
0310
0311   FPGAWord fpgaphi = tracklet->proj(N_LAYER + disk - 1).fpgaphiproj();
0312
0313   int iphivmRaw = fpgaphi.value() >> (fpgaphi.nbits() - 5);
0314   int iphi = iphivmRaw / (32 / settings_.nallstubs(disk + N_LAYER - 1));
0315
0316   addProjectionDisk(disk, iphi, trackletprojdisks_[disk - 1][iphi], tracklet);
0317 }
0318
0319 bool TrackletCalculatorDisplaced::addLayerProj(Tracklet* tracklet, int layer) {
0320   assert(layer > 0);
0321
0322   FPGAWord fpgaz = tracklet->proj(layer - 1).fpgarzproj();
0323   FPGAWord fpgaphi = tracklet->proj(layer - 1).fpgaphiproj();
0324
0325   if (fpgaz.atExtreme())
0326     return false;
0327
0328   if (std::abs(fpgaz.value() * settings_.kz()) > settings_.zlength())
0329     return false;
0330
0331   int iphivmRaw = fpgaphi.value() >> (fpgaphi.nbits() - 5);
0332   int iphi = iphivmRaw / (32 / settings_.nallstubs(layer - 1));
0333
0334   addProjection(layer, iphi, trackletprojlayers_[layer - 1][iphi], tracklet);
0335
0336   return true;
0337 }
0338
0339 void TrackletCalculatorDisplaced::addProjection(int layer,
0340                                                 int iphi,
0341                                                 TrackletProjectionsMemory* trackletprojs,
0342                                                 Tracklet* tracklet) {
0343   if (trackletprojs == nullptr) {
0344     if (settings_.warnNoMem()) {
0345       edm::LogVerbatim("Tracklet") << "No projection memory exists in " << getName() << " for layer = " << layer
0346                                    << " iphi = " << iphi + 1;
0347     }
0348     return;
0349   }
0350   assert(trackletprojs != nullptr);
0351   trackletprojs->addProj(tracklet);
0352 }
0353
0354 void TrackletCalculatorDisplaced::addProjectionDisk(int disk,
0355                                                     int iphi,
0356                                                     TrackletProjectionsMemory* trackletprojs,
0357                                                     Tracklet* tracklet) {
0358   if (trackletprojs == nullptr) {
0359     if (layer_ == 3 && abs(disk) == 3)
0360       return;  //L3L4 projections to D3 are not used.
0361     if (settings_.warnNoMem()) {
0362       edm::LogVerbatim("Tracklet") << "No projection memory exists in " << getName() << " for disk = " << abs(disk)
0363                                    << " iphi = " << iphi + 1;
0364     }
0365     return;
0366   }
0367   assert(trackletprojs != nullptr);
0368   if (settings_.debugTracklet())
0369     edm::LogVerbatim("Tracklet") << getName() << " adding projection to " << trackletprojs->getName();
0370   trackletprojs->addProj(tracklet);
0371 }
0372
0373 bool TrackletCalculatorDisplaced::LLLSeeding(const Stub* innerFPGAStub,
0374                                              const L1TStub* innerStub,
0375                                              const Stub* middleFPGAStub,
0376                                              const L1TStub* middleStub,
0377                                              const Stub* outerFPGAStub,
0378                                              const L1TStub* outerStub) {
0379   if (settings_.debugTracklet())
0380     edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced " << getName() << " " << layer_
0381                                  << " trying stub triplet in layer (L L L): " << innerFPGAStub->layer().value() << " "
0382                                  << middleFPGAStub->layer().value() << " " << outerFPGAStub->layer().value();
0383
0384   assert(outerFPGAStub->layerdisk() < N_LAYER);
0385
0386   double r1 = innerStub->r();
0387   double z1 = innerStub->z();
0388   double phi1 = innerStub->phi();
0389
0390   double r2 = middleStub->r();
0391   double z2 = middleStub->z();
0392   double phi2 = middleStub->phi();
0393
0394   double r3 = outerStub->r();
0395   double z3 = outerStub->z();
0396   double phi3 = outerStub->phi();
0397
0398   int take3 = 0;
0399   if (layer_ == 5)
0400     take3 = 1;
0401   unsigned ndisks = 0;
0402
0403   double rinv, phi0, d0, t, z0;
0404
0405   Projection projs[N_LAYER + N_DISK];
0406
0407   double phiproj[N_LAYER - 2], zproj[N_LAYER - 2], phider[N_LAYER - 2], zder[N_LAYER - 2];
0408   double phiprojdisk[N_DISK], rprojdisk[N_DISK], phiderdisk[N_DISK], rderdisk[N_DISK];
0409
0410   exacttracklet(r1,
0411                 z1,
0412                 phi1,
0413                 r2,
0414                 z2,
0415                 phi2,
0416                 r3,
0417                 z3,
0418                 phi3,
0419                 take3,
0420                 rinv,
0421                 phi0,
0422                 d0,
0423                 t,
0424                 z0,
0425                 phiproj,
0426                 zproj,
0427                 phiprojdisk,
0428                 rprojdisk,
0429                 phider,
0430                 zder,
0431                 phiderdisk,
0432                 rderdisk);
0433   if (settings_.debugTracklet())
0434     edm::LogVerbatim("Tracklet") << __LINE__ << ":" << __FILE__ << "LLL Exact values " << innerFPGAStub->isBarrel()
0435                                  << middleFPGAStub->isBarrel() << outerFPGAStub->isBarrel() << " " << phi1 << ", " << z1
0436                                  << ", " << r1 << ", " << phi2 << ", " << z2 << ", " << r2 << ", " << phi3 << ", " << z3
0437                                  << ", " << r3 << endl;
0438
0439   if (settings_.useapprox()) {
0440     phi1 = innerFPGAStub->phiapprox(phimin_, phimax_);
0441     z1 = innerFPGAStub->zapprox();
0442     r1 = innerFPGAStub->rapprox();
0443
0444     phi2 = middleFPGAStub->phiapprox(phimin_, phimax_);
0445     z2 = middleFPGAStub->zapprox();
0446     r2 = middleFPGAStub->rapprox();
0447
0448     phi3 = outerFPGAStub->phiapprox(phimin_, phimax_);
0449     z3 = outerFPGAStub->zapprox();
0450     r3 = outerFPGAStub->rapprox();
0451   }
0452
0453   if (settings_.debugTracklet())
0454     edm::LogVerbatim("Tracklet") << __LINE__ << ":" << __FILE__ << "LLL Approx values " << innerFPGAStub->isBarrel()
0455                                  << middleFPGAStub->isBarrel() << outerFPGAStub->isBarrel() << " " << phi1 << ", " << z1
0456                                  << ", " << r1 << ", " << phi2 << ", " << z2 << ", " << r2 << ", " << phi3 << ", " << z3
0457                                  << ", " << r3 << endl;
0458
0459   double rinvapprox, phi0approx, d0approx, tapprox, z0approx;
0460   double phiprojapprox[N_LAYER - 2], zprojapprox[N_LAYER - 2], phiderapprox[N_LAYER - 2], zderapprox[N_LAYER - 2];
0461   double phiprojdiskapprox[N_DISK], rprojdiskapprox[N_DISK];
0462   double phiderdiskapprox[N_DISK], rderdiskapprox[N_DISK];
0463
0464   //TODO: implement the actual integer calculation
0465   if (settings_.useapprox()) {
0466     approxtracklet(r1,
0467                    z1,
0468                    phi1,
0469                    r2,
0470                    z2,
0471                    phi2,
0472                    r3,
0473                    z3,
0474                    phi3,
0475                    take3,
0476                    ndisks,
0477                    rinvapprox,
0478                    phi0approx,
0479                    d0approx,
0480                    tapprox,
0481                    z0approx,
0482                    phiprojapprox,
0483                    zprojapprox,
0484                    phiderapprox,
0485                    zderapprox,
0486                    phiprojdiskapprox,
0487                    rprojdiskapprox,
0488                    phiderdiskapprox,
0489                    rderdiskapprox);
0490   } else {
0491     rinvapprox = rinv;
0492     phi0approx = phi0;
0493     d0approx = d0;
0494     tapprox = t;
0495     z0approx = z0;
0496
0497     for (unsigned int i = 0; i < toR_.size(); ++i) {
0498       phiprojapprox[i] = phiproj[i];
0499       zprojapprox[i] = zproj[i];
0500       phiderapprox[i] = phider[i];
0501       zderapprox[i] = zder[i];
0502     }
0503
0504     for (unsigned int i = 0; i < toZ_.size(); ++i) {
0505       phiprojdiskapprox[i] = phiprojdisk[i];
0506       rprojdiskapprox[i] = rprojdisk[i];
0507       phiderdiskapprox[i] = phiderdisk[i];
0508       rderdiskapprox[i] = rderdisk[i];
0509     }
0510   }
0511
0512   //store the approcximate results
0513
0514   if (settings_.debugTracklet()) {
0515     edm::LogVerbatim("Tracklet") << "rinvapprox: " << rinvapprox << " rinv: " << rinv << endl;
0516     edm::LogVerbatim("Tracklet") << "phi0approx: " << phi0approx << " phi0: " << phi0 << endl;
0517     edm::LogVerbatim("Tracklet") << "d0approx: " << d0approx << " d0: " << d0 << endl;
0518     edm::LogVerbatim("Tracklet") << "tapprox: " << tapprox << " t: " << t << endl;
0519     edm::LogVerbatim("Tracklet") << "z0approx: " << z0approx << " z0: " << z0 << endl;
0520   }
0521
0522   for (unsigned int i = 0; i < toR_.size(); ++i) {
0523     if (settings_.debugTracklet()) {
0524       edm::LogVerbatim("Tracklet") << "phiprojapprox[" << i << "]: " << phiprojapprox[i] << " phiproj[" << i
0525                                    << "]: " << phiproj[i] << endl;
0526       edm::LogVerbatim("Tracklet") << "zprojapprox[" << i << "]: " << zprojapprox[i] << " zproj[" << i
0527                                    << "]: " << zproj[i] << endl;
0528       edm::LogVerbatim("Tracklet") << "phiderapprox[" << i << "]: " << phiderapprox[i] << " phider[" << i
0529                                    << "]: " << phider[i] << endl;
0530       edm::LogVerbatim("Tracklet") << "zderapprox[" << i << "]: " << zderapprox[i] << " zder[" << i << "]: " << zder[i]
0531                                    << endl;
0532     }
0533   }
0534
0535   for (unsigned int i = 0; i < toZ_.size(); ++i) {
0536     if (settings_.debugTracklet()) {
0537       edm::LogVerbatim("Tracklet") << "phiprojdiskapprox[" << i << "]: " << phiprojdiskapprox[i] << " phiprojdisk[" << i
0538                                    << "]: " << phiprojdisk[i] << endl;
0539       edm::LogVerbatim("Tracklet") << "rprojdiskapprox[" << i << "]: " << rprojdiskapprox[i] << " rprojdisk[" << i
0540                                    << "]: " << rprojdisk[i] << endl;
0541       edm::LogVerbatim("Tracklet") << "phiderdiskapprox[" << i << "]: " << phiderdiskapprox[i] << " phiderdisk[" << i
0542                                    << "]: " << phiderdisk[i] << endl;
0543       edm::LogVerbatim("Tracklet") << "rderdiskapprox[" << i << "]: " << rderdiskapprox[i] << " rderdisk[" << i
0544                                    << "]: " << rderdisk[i] << endl;
0545     }
0546   }
0547
0548   //now binary
0549   double krinv = settings_.kphi1() / settings_.kr() * pow(2, settings_.rinv_shift()),
0550          kphi0 = settings_.kphi1() * pow(2, settings_.phi0_shift()),
0551          kt = settings_.kz() / settings_.kr() * pow(2, settings_.t_shift()),
0552          kz0 = settings_.kz() * pow(2, settings_.z0_shift()),
0553          kphiproj = settings_.kphi1() * pow(2, settings_.SS_phiL_shift()),
0554          kphider = settings_.kphi1() / settings_.kr() * pow(2, settings_.SS_phiderL_shift()),
0555          kzproj = settings_.kz() * pow(2, settings_.PS_zL_shift()),
0556          kzder = settings_.kz() / settings_.kr() * pow(2, settings_.PS_zderL_shift()),
0557          kphiprojdisk = settings_.kphi1() * pow(2, settings_.SS_phiD_shift()),
0558          kphiderdisk = settings_.kphi1() / settings_.kr() * pow(2, settings_.SS_phiderD_shift()),
0559          krprojdisk = settings_.kr() * pow(2, settings_.PS_rD_shift()),
0560          krderdisk = settings_.kr() / settings_.kz() * pow(2, settings_.PS_rderD_shift());
0561
0562   int irinv, iphi0, id0, it, iz0;
0563   int iphiproj[N_LAYER - 2], izproj[N_LAYER - 2], iphider[N_LAYER - 2], izder[N_LAYER - 2];
0564   int iphiprojdisk[N_DISK], irprojdisk[N_DISK], iphiderdisk[N_DISK], irderdisk[N_DISK];
0565
0566   //store the binary results
0567   irinv = rinvapprox / krinv;
0568   iphi0 = phi0approx / kphi0;
0569   id0 = d0approx / settings_.kd0();
0570   it = tapprox / kt;
0571   iz0 = z0approx / kz0;
0572
0573   bool success = true;
0574   if (std::abs(rinvapprox) > settings_.rinvcut()) {
0575     if (settings_.debugTracklet())
0576       edm::LogVerbatim("Tracklet") << "TrackletCalculator::LLL Seeding irinv too large: " << rinvapprox << "(" << irinv
0577                                    << ")";
0578     success = false;
0579   }
0580   if (std::abs(z0approx) > settings_.disp_z0cut()) {
0581     if (settings_.debugTracklet())
0582       edm::LogVerbatim("Tracklet") << "Failed tracklet z0 cut " << z0approx << " in layer " << layer_;
0583     success = false;
0584   }
0585   if (std::abs(d0approx) > settings_.maxd0()) {
0586     if (settings_.debugTracklet())
0587       edm::LogVerbatim("Tracklet") << "Failed tracklet approx d0 cut " << d0approx;
0588     success = false;
0589   }
0590   if (std::abs(d0) > settings_.maxd0()) {
0591     if (settings_.debugTracklet())
0592       edm::LogVerbatim("Tracklet") << "Failed tracklet exact d0 cut " << d0;
0593     success = false;
0594   }
0595
0596   if (!success) {
0597     return false;
0598   }
0599
0600   double phicritapprox = phi0approx - asin((0.5 * settings_.rcrit() * rinvapprox) + (d0approx / settings_.rcrit()));
0601   int phicrit = iphi0 - 2 * irinv - 2 * id0;
0602
0603   int iphicritmincut = settings_.phicritminmc() / globals_->ITC_L1L2()->phi0_final.K();
0604   int iphicritmaxcut = settings_.phicritmaxmc() / globals_->ITC_L1L2()->phi0_final.K();
0605
0606   bool keepapprox = (phicritapprox > settings_.phicritminmc()) && (phicritapprox < settings_.phicritmaxmc()),
0607        keep = (phicrit > iphicritmincut) && (phicrit < iphicritmaxcut);
0608
0609   if (settings_.debugTracklet())
0610     if (keep && !keepapprox)
0611       edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced::LLLSeeding tracklet kept with exact phicrit cut "
0612                                       "but not approximate, phicritapprox: "
0613                                    << phicritapprox;
0614   if (settings_.usephicritapprox()) {
0615     if (!keepapprox) {
0616       return false;
0617     }
0618   } else {
0619     if (!keep) {
0620       return false;
0621     }
0622   }
0623
0624   for (unsigned int i = 0; i < toR_.size(); ++i) {
0625     iphiproj[i] = phiprojapprox[i] / kphiproj;
0626     izproj[i] = zprojapprox[i] / kzproj;
0627
0628     iphider[i] = phiderapprox[i] / kphider;
0629     izder[i] = zderapprox[i] / kzder;
0630
0631     //check that z projection is in range
0632     if (izproj[i] < -(1 << (settings_.nzbitsstub(0) - 1)))
0633       continue;
0634     if (izproj[i] >= (1 << (settings_.nzbitsstub(0) - 1)))
0635       continue;
0636
0637     //check that phi projection is in range
0638     if (iphiproj[i] >= (1 << settings_.nphibitsstub(N_LAYER - 1)) - 1)
0639       continue;
0640     if (iphiproj[i] <= 0)
0641       continue;
0642
0643     //adjust number of bits for phi and z projection
0644     if (rproj_[i] < settings_.rPS2S()) {
0645       iphiproj[i] >>= (settings_.nphibitsstub(N_LAYER - 1) - settings_.nphibitsstub(0));
0646       if (iphiproj[i] >= (1 << settings_.nphibitsstub(0)) - 1)
0647         iphiproj[i] = (1 << settings_.nphibitsstub(0)) - 2;  //-2 not to hit atExtreme
0648     } else {
0649       izproj[i] >>= (settings_.nzbitsstub(0) - settings_.nzbitsstub(N_LAYER - 1));
0650     }
0651
0652     if (rproj_[i] < settings_.rPS2S()) {
0653       if (iphider[i] < -(1 << (settings_.nbitsphiprojderL123() - 1))) {
0654         iphider[i] = -(1 << (settings_.nbitsphiprojderL123() - 1));
0655       }
0656       if (iphider[i] >= (1 << (settings_.nbitsphiprojderL123() - 1))) {
0657         iphider[i] = (1 << (settings_.nbitsphiprojderL123() - 1)) - 1;
0658       }
0659     } else {
0660       if (iphider[i] < -(1 << (settings_.nbitsphiprojderL456() - 1))) {
0661         iphider[i] = -(1 << (settings_.nbitsphiprojderL456() - 1));
0662       }
0663       if (iphider[i] >= (1 << (settings_.nbitsphiprojderL456() - 1))) {
0664         iphider[i] = (1 << (settings_.nbitsphiprojderL456() - 1)) - 1;
0665       }
0666     }
0667
0668     projs[lproj_[i] - 1].init(settings_,
0669                               lproj_[i] - 1,
0670                               iphiproj[i],
0671                               izproj[i],
0672                               iphider[i],
0673                               izder[i],
0674                               phiproj[i],
0675                               zproj[i],
0676                               phider[i],
0677                               zder[i],
0678                               phiprojapprox[i],
0679                               zprojapprox[i],
0680                               phiderapprox[i],
0681                               zderapprox[i],
0682                               false);
0683   }
0684
0685   if (std::abs(it * kt) > 1.0) {
0686     for (unsigned int i = 0; i < toZ_.size(); ++i) {
0687       iphiprojdisk[i] = phiprojdiskapprox[i] / kphiprojdisk;
0688       irprojdisk[i] = rprojdiskapprox[i] / krprojdisk;
0689
0690       iphiderdisk[i] = phiderdiskapprox[i] / kphiderdisk;
0691       irderdisk[i] = rderdiskapprox[i] / krderdisk;
0692
0693       //check phi projection in range
0694       if (iphiprojdisk[i] <= 0)
0695         continue;
0696       if (iphiprojdisk[i] >= (1 << settings_.nphibitsstub(0)) - 1)
0697         continue;
0698
0699       //check r projection in range
0700       if (rprojdiskapprox[i] < settings_.rmindisk() || rprojdiskapprox[i] >= settings_.rmaxdisk())
0701         continue;
0702
0703       projs[N_LAYER + i].init(settings_,
0704                               N_LAYER + i,
0705                               iphiprojdisk[i],
0706                               irprojdisk[i],
0707                               iphiderdisk[i],
0708                               irderdisk[i],
0709                               phiprojdisk[i],
0710                               rprojdisk[i],
0711                               phiderdisk[i],
0712                               rderdisk[i],
0713                               phiprojdiskapprox[i],
0714                               rprojdiskapprox[i],
0715                               phiderdisk[i],
0716                               rderdisk[i],
0717                               false);
0718     }
0719   }
0720
0721   if (settings_.writeMonitorData("TrackletPars")) {
0722     globals_->ofstream("trackletpars.txt")
0723         << layer_ << " , " << rinv << " , " << rinvapprox << " , " << phi0 << " , " << phi0approx << " , " << t << " , "
0724         << tapprox << " , " << z0 << " , " << z0approx << " , " << d0 << " , " << d0approx << endl;
0725   }
0726
0727   Tracklet* tracklet = new Tracklet(settings_,
0728                                     iSeed_,
0729                                     innerFPGAStub,
0730                                     middleFPGAStub,
0731                                     outerFPGAStub,
0732                                     rinv,
0733                                     phi0,
0734                                     d0,
0735                                     z0,
0736                                     t,
0737                                     rinvapprox,
0738                                     phi0approx,
0739                                     d0approx,
0740                                     z0approx,
0741                                     tapprox,
0742                                     irinv,
0743                                     iphi0,
0744                                     id0,
0745                                     iz0,
0746                                     it,
0747                                     projs,
0748                                     false);
0749
0750   if (settings_.debugTracklet())
0751     edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced " << getName()
0752                                  << " Found LLL tracklet in sector = " << iSector_ << " phi0 = " << phi0;
0753
0754   tracklet->setTrackletIndex(trackletpars_->nTracklets());
0755   tracklet->setTCIndex(TCIndex_);
0756
0757   if (settings_.writeMonitorData("Seeds")) {
0758     ofstream fout("seeds.txt", ofstream::app);
0759     fout << __FILE__ << ":" << __LINE__ << " " << name_ << "_" << iSector_ << " " << tracklet->getISeed() << endl;
0760     fout.close();
0761   }
0762   trackletpars_->addTracklet(tracklet);
0763
0764   bool addL5 = false;
0765   bool addL6 = false;
0766   for (unsigned int j = 0; j < toR_.size(); j++) {
0767     bool added = false;
0768
0769     if (settings_.debugTracklet())
0770       edm::LogVerbatim("Tracklet") << "adding layer projection " << j << "/" << toR_.size() << " " << lproj_[j];
0771     if (tracklet->validProj(lproj_[j] - 1)) {
0772       added = addLayerProj(tracklet, lproj_[j]);
0773       if (added && lproj_[j] == 5)
0774         addL5 = true;
0775       if (added && lproj_[j] == 6)
0776         addL6 = true;
0777     }
0778   }
0779
0780   for (unsigned int j = 0; j < toZ_.size(); j++) {
0781     int disk = dproj_[j];
0782     if (disk == 0)
0783       continue;
0784     if (disk == 2 && addL5)
0785       continue;
0786     if (disk == 1 && addL6)
0787       continue;
0788     if (it < 0)
0789       disk = -disk;
0790     if (settings_.debugTracklet())
0791       edm::LogVerbatim("Tracklet") << "adding disk projection " << j << "/" << toZ_.size() << " " << disk;
0792     if (tracklet->validProj(N_LAYER + abs(disk) - 1)) {
0793       addDiskProj(tracklet, disk);
0794     }
0795   }
0796
0797   return true;
0798 }
0799
0800 bool TrackletCalculatorDisplaced::DDLSeeding(const Stub* innerFPGAStub,
0801                                              const L1TStub* innerStub,
0802                                              const Stub* middleFPGAStub,
0803                                              const L1TStub* middleStub,
0804                                              const Stub* outerFPGAStub,
0805                                              const L1TStub* outerStub) {
0806   if (settings_.debugTracklet())
0807     edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced " << getName() << " " << layer_
0808                                  << " trying stub triplet in  (L2 D1 D2): " << innerFPGAStub->layer().value() << " "
0809                                  << middleFPGAStub->disk().value() << " " << outerFPGAStub->disk().value();
0810
0811   int take3 = 1;  //D1D2L2
0812   unsigned ndisks = 2;
0813
0814   double r1 = innerStub->r();
0815   double z1 = innerStub->z();
0816   double phi1 = innerStub->phi();
0817
0818   double r2 = middleStub->r();
0819   double z2 = middleStub->z();
0820   double phi2 = middleStub->phi();
0821
0822   double r3 = outerStub->r();
0823   double z3 = outerStub->z();
0824   double phi3 = outerStub->phi();
0825
0826   double rinv, phi0, d0, t, z0;
0827
0828   double phiproj[N_LAYER - 2], zproj[N_LAYER - 2], phider[N_LAYER - 2], zder[N_LAYER - 2];
0829   double phiprojdisk[N_DISK], rprojdisk[N_DISK], phiderdisk[N_DISK], rderdisk[N_DISK];
0830
0831   exacttracklet(r1,
0832                 z1,
0833                 phi1,
0834                 r2,
0835                 z2,
0836                 phi2,
0837                 r3,
0838                 z3,
0839                 phi3,
0840                 take3,
0841                 rinv,
0842                 phi0,
0843                 d0,
0844                 t,
0845                 z0,
0846                 phiproj,
0847                 zproj,
0848                 phiprojdisk,
0849                 rprojdisk,
0850                 phider,
0851                 zder,
0852                 phiderdisk,
0853                 rderdisk);
0854
0855   if (settings_.debugTracklet())
0856     edm::LogVerbatim("Tracklet") << __LINE__ << ":" << __FILE__ << " DLL Exact values " << innerFPGAStub->isBarrel()
0857                                  << middleFPGAStub->isBarrel() << outerFPGAStub->isBarrel() << " " << phi1 << ", " << z1
0858                                  << ", " << r1 << ", " << phi2 << ", " << z2 << ", " << r2 << ", " << phi3 << ", " << z3
0859                                  << ", " << r3 << endl;
0860
0861   if (settings_.useapprox()) {
0862     phi1 = innerFPGAStub->phiapprox(phimin_, phimax_);
0863     z1 = innerFPGAStub->zapprox();
0864     r1 = innerFPGAStub->rapprox();
0865
0866     phi2 = middleFPGAStub->phiapprox(phimin_, phimax_);
0867     z2 = middleFPGAStub->zapprox();
0868     r2 = middleFPGAStub->rapprox();
0869
0870     phi3 = outerFPGAStub->phiapprox(phimin_, phimax_);
0871     z3 = outerFPGAStub->zapprox();
0872     r3 = outerFPGAStub->rapprox();
0873   }
0874
0875   if (settings_.debugTracklet())
0876     edm::LogVerbatim("Tracklet") << __LINE__ << ":" << __FILE__ << "DLL Approx values " << innerFPGAStub->isBarrel()
0877                                  << middleFPGAStub->isBarrel() << outerFPGAStub->isBarrel() << " " << phi1 << ", " << z1
0878                                  << ", " << r1 << ", " << phi2 << ", " << z2 << ", " << r2 << ", " << phi3 << ", " << z3
0879                                  << ", " << r3 << endl;
0880
0881   double rinvapprox, phi0approx, d0approx, tapprox, z0approx;
0882   double phiprojapprox[N_LAYER - 2], zprojapprox[N_LAYER - 2], phiderapprox[N_LAYER - 2], zderapprox[N_LAYER - 2];
0883   double phiprojdiskapprox[N_DISK], rprojdiskapprox[N_DISK];
0884   double phiderdiskapprox[N_DISK], rderdiskapprox[N_DISK];
0885
0886   //TODO: implement the actual integer calculation
0887   if (settings_.useapprox()) {
0888     approxtracklet(r1,
0889                    z1,
0890                    phi1,
0891                    r2,
0892                    z2,
0893                    phi2,
0894                    r3,
0895                    z3,
0896                    phi3,
0897                    take3,
0898                    ndisks,
0899                    rinvapprox,
0900                    phi0approx,
0901                    d0approx,
0902                    tapprox,
0903                    z0approx,
0904                    phiprojapprox,
0905                    zprojapprox,
0906                    phiderapprox,
0907                    zderapprox,
0908                    phiprojdiskapprox,
0909                    rprojdiskapprox,
0910                    phiderdiskapprox,
0911                    rderdiskapprox);
0912   } else {
0913     rinvapprox = rinv;
0914     phi0approx = phi0;
0915     d0approx = d0;
0916     tapprox = t;
0917     z0approx = z0;
0918
0919     for (unsigned int i = 0; i < toR_.size(); ++i) {
0920       phiprojapprox[i] = phiproj[i];
0921       zprojapprox[i] = zproj[i];
0922       phiderapprox[i] = phider[i];
0923       zderapprox[i] = zder[i];
0924     }
0925
0926     for (unsigned int i = 0; i < toZ_.size(); ++i) {
0927       phiprojdiskapprox[i] = phiprojdisk[i];
0928       rprojdiskapprox[i] = rprojdisk[i];
0929       phiderdiskapprox[i] = phiderdisk[i];
0930       rderdiskapprox[i] = rderdisk[i];
0931     }
0932   }
0933
0934   //store the approcximate results
0935   if (settings_.debugTracklet()) {
0936     edm::LogVerbatim("Tracklet") << "rinvapprox: " << rinvapprox << " rinv: " << rinv << endl;
0937     edm::LogVerbatim("Tracklet") << "phi0approx: " << phi0approx << " phi0: " << phi0 << endl;
0938     edm::LogVerbatim("Tracklet") << "d0approx: " << d0approx << " d0: " << d0 << endl;
0939     edm::LogVerbatim("Tracklet") << "tapprox: " << tapprox << " t: " << t << endl;
0940     edm::LogVerbatim("Tracklet") << "z0approx: " << z0approx << " z0: " << z0 << endl;
0941
0942     for (unsigned int i = 0; i < toR_.size(); ++i) {
0943       edm::LogVerbatim("Tracklet") << "phiprojapprox[" << i << "]: " << phiprojapprox[i] << " phiproj[" << i
0944                                    << "]: " << phiproj[i] << endl;
0945       edm::LogVerbatim("Tracklet") << "zprojapprox[" << i << "]: " << zprojapprox[i] << " zproj[" << i
0946                                    << "]: " << zproj[i] << endl;
0947       edm::LogVerbatim("Tracklet") << "phiderapprox[" << i << "]: " << phiderapprox[i] << " phider[" << i
0948                                    << "]: " << phider[i] << endl;
0949       edm::LogVerbatim("Tracklet") << "zderapprox[" << i << "]: " << zderapprox[i] << " zder[" << i << "]: " << zder[i]
0950                                    << endl;
0951     }
0952
0953     for (unsigned int i = 0; i < toZ_.size(); ++i) {
0954       edm::LogVerbatim("Tracklet") << "phiprojdiskapprox[" << i << "]: " << phiprojdiskapprox[i] << " phiprojdisk[" << i
0955                                    << "]: " << phiprojdisk[i] << endl;
0956       edm::LogVerbatim("Tracklet") << "rprojdiskapprox[" << i << "]: " << rprojdiskapprox[i] << " rprojdisk[" << i
0957                                    << "]: " << rprojdisk[i] << endl;
0958       edm::LogVerbatim("Tracklet") << "phiderdiskapprox[" << i << "]: " << phiderdiskapprox[i] << " phiderdisk[" << i
0959                                    << "]: " << phiderdisk[i] << endl;
0960       edm::LogVerbatim("Tracklet") << "rderdiskapprox[" << i << "]: " << rderdiskapprox[i] << " rderdisk[" << i
0961                                    << "]: " << rderdisk[i] << endl;
0962     }
0963   }
0964
0965   //now binary
0966   double krinv = settings_.kphi1() / settings_.kr() * pow(2, settings_.rinv_shift()),
0967          kphi0 = settings_.kphi1() * pow(2, settings_.phi0_shift()),
0968          kt = settings_.kz() / settings_.kr() * pow(2, settings_.t_shift()),
0969          kz0 = settings_.kz() * pow(2, settings_.z0_shift()),
0970          kphiproj = settings_.kphi1() * pow(2, settings_.SS_phiL_shift()),
0971          kphider = settings_.kphi1() / settings_.kr() * pow(2, settings_.SS_phiderL_shift()),
0972          kzproj = settings_.kz() * pow(2, settings_.PS_zL_shift()),
0973          kzder = settings_.kz() / settings_.kr() * pow(2, settings_.PS_zderL_shift()),
0974          kphiprojdisk = settings_.kphi1() * pow(2, settings_.SS_phiD_shift()),
0975          kphiderdisk = settings_.kphi1() / settings_.kr() * pow(2, settings_.SS_phiderD_shift()),
0976          krprojdisk = settings_.kr() * pow(2, settings_.PS_rD_shift()),
0977          krderdisk = settings_.kr() / settings_.kz() * pow(2, settings_.PS_rderD_shift());
0978
0979   int irinv, iphi0, id0, it, iz0;
0980   int iphiproj[N_LAYER - 2], izproj[N_LAYER - 2], iphider[N_LAYER - 2], izder[N_LAYER - 2];
0981   int iphiprojdisk[N_DISK], irprojdisk[N_DISK], iphiderdisk[N_DISK], irderdisk[N_DISK];
0982
0983   //store the binary results
0984   irinv = rinvapprox / krinv;
0985   iphi0 = phi0approx / kphi0;
0986   id0 = d0approx / settings_.kd0();
0987   it = tapprox / kt;
0988   iz0 = z0approx / kz0;
0989
0990   bool success = true;
0991   if ((std::abs(rinvapprox) > settings_.rinvcut()) || (!edm::isFinite(rinvapprox))) {
0992     if (settings_.debugTracklet())
0993       edm::LogVerbatim("Tracklet") << "TrackletCalculator::DDL Seeding irinv too large: " << rinvapprox << "(" << irinv
0994                                    << ")";
0995     success = false;
0996   }
0997   if (std::abs(z0approx) > settings_.disp_z0cut()) {
0998     if (settings_.debugTracklet())
0999       edm::LogVerbatim("Tracklet") << "Failed tracklet z0 cut " << z0approx;
1000     success = false;
1001   }
1002   if ((std::abs(d0approx) > settings_.maxd0()) || (!edm::isFinite(d0approx))) {
1003     if (settings_.debugTracklet())
1004       edm::LogVerbatim("Tracklet") << "Failed tracklet approx d0 cut " << d0approx;
1005     success = false;
1006   }
1007   if (std::abs(d0) > settings_.maxd0()) {
1008     if (settings_.debugTracklet())
1009       edm::LogVerbatim("Tracklet") << "Failed tracklet exact d0 cut " << d0;
1010     success = false;
1011   }
1012
1013   if (!success)
1014     return false;
1015
1016   double phicritapprox = phi0approx - asin((0.5 * settings_.rcrit() * rinvapprox) + (d0approx / settings_.rcrit()));
1017   int phicrit = iphi0 - 2 * irinv - 2 * id0;
1018
1019   int iphicritmincut = settings_.phicritminmc() / globals_->ITC_L1L2()->phi0_final.K();
1020   int iphicritmaxcut = settings_.phicritmaxmc() / globals_->ITC_L1L2()->phi0_final.K();
1021
1022   bool keepapprox = (phicritapprox > settings_.phicritminmc()) && (phicritapprox < settings_.phicritmaxmc()),
1023        keep = (phicrit > iphicritmincut) && (phicrit < iphicritmaxcut);
1024
1025   if (settings_.debugTracklet())
1026     if (keep && !keepapprox)
1027       edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced::DDLSeeding tracklet kept with exact phicrit cut "
1028                                       "but not approximate, phicritapprox: "
1029                                    << phicritapprox;
1030   if (settings_.usephicritapprox()) {
1031     if (!keepapprox)
1032       return false;
1033   } else {
1034     if (!keep)
1035       return false;
1036   }
1037
1038   Projection projs[N_LAYER + N_DISK];
1039
1040   for (unsigned int i = 0; i < toR_.size(); ++i) {
1041     iphiproj[i] = phiprojapprox[i] / kphiproj;
1042     izproj[i] = zprojapprox[i] / kzproj;
1043
1044     iphider[i] = phiderapprox[i] / kphider;
1045     izder[i] = zderapprox[i] / kzder;
1046
1047     //check that z projection in range
1048     if (izproj[i] < -(1 << (settings_.nzbitsstub(0) - 1)))
1049       continue;
1050     if (izproj[i] >= (1 << (settings_.nzbitsstub(0) - 1)))
1051       continue;
1052
1053     //check that phi projection in range
1054     if (iphiproj[i] >= (1 << settings_.nphibitsstub(5)) - 1)
1055       continue;
1056     if (iphiproj[i] <= 0)
1057       continue;
1058
1059     if (rproj_[i] < settings_.rPS2S()) {
1060       iphiproj[i] >>= (settings_.nphibitsstub(5) - settings_.nphibitsstub(0));
1061     } else {
1062       izproj[i] >>= (settings_.nzbitsstub(0) - settings_.nzbitsstub(5));
1063     }
1064
1065     if (rproj_[i] < settings_.rPS2S()) {
1066       if (iphider[i] < -(1 << (settings_.nbitsphiprojderL123() - 1)))
1067         iphider[i] = -(1 << (settings_.nbitsphiprojderL123() - 1));
1068       if (iphider[i] >= (1 << (settings_.nbitsphiprojderL123() - 1)))
1069         iphider[i] = (1 << (settings_.nbitsphiprojderL123() - 1)) - 1;
1070     } else {
1071       if (iphider[i] < -(1 << (settings_.nbitsphiprojderL456() - 1)))
1072         iphider[i] = -(1 << (settings_.nbitsphiprojderL456() - 1));
1073       if (iphider[i] >= (1 << (settings_.nbitsphiprojderL456() - 1)))
1074         iphider[i] = (1 << (settings_.nbitsphiprojderL456() - 1)) - 1;
1075     }
1076
1077     projs[lproj_[i] - 1].init(settings_,
1078                               lproj_[i] - 1,
1079                               iphiproj[i],
1080                               izproj[i],
1081                               iphider[i],
1082                               izder[i],
1083                               phiproj[i],
1084                               zproj[i],
1085                               phider[i],
1086                               zder[i],
1087                               phiprojapprox[i],
1088                               zprojapprox[i],
1089                               phiderapprox[i],
1090                               zderapprox[i],
1091                               false);
1092   }
1093
1094   if (std::abs(it * kt) > 1.0) {
1095     for (unsigned int i = 0; i < toZ_.size(); ++i) {
1096       iphiprojdisk[i] = phiprojdiskapprox[i] / kphiprojdisk;
1097       irprojdisk[i] = rprojdiskapprox[i] / krprojdisk;
1098
1099       iphiderdisk[i] = phiderdiskapprox[i] / kphiderdisk;
1100       irderdisk[i] = rderdiskapprox[i] / krderdisk;
1101
1102       if (iphiprojdisk[i] <= 0)
1103         continue;
1104       if (iphiprojdisk[i] >= (1 << settings_.nphibitsstub(0)) - 1)
1105         continue;
1106
1107       if (irprojdisk[i] < settings_.rmindisk() / krprojdisk || irprojdisk[i] >= settings_.rmaxdisk() / krprojdisk)
1108         continue;
1109
1110       projs[N_LAYER + i + 2].init(settings_,
1111                                   N_LAYER + i + 2,
1112                                   iphiprojdisk[i],
1113                                   irprojdisk[i],
1114                                   iphiderdisk[i],
1115                                   irderdisk[i],
1116                                   phiprojdisk[i],
1117                                   rprojdisk[i],
1118                                   phiderdisk[i],
1119                                   rderdisk[i],
1120                                   phiprojdiskapprox[i],
1121                                   rprojdiskapprox[i],
1122                                   phiderdisk[i],
1123                                   rderdisk[i],
1124                                   false);
1125     }
1126   }
1127
1128   if (settings_.writeMonitorData("TrackletPars")) {
1129     globals_->ofstream("trackletpars.txt")
1130         << layer_ << " , " << rinv << " , " << rinvapprox << " , " << phi0 << " , " << phi0approx << " , " << t << " , "
1131         << tapprox << " , " << z0 << " , " << z0approx << " , " << d0 << " , " << d0approx << endl;
1132   }
1133
1134   Tracklet* tracklet = new Tracklet(settings_,
1135                                     iSeed_,
1136                                     innerFPGAStub,
1137                                     middleFPGAStub,
1138                                     outerFPGAStub,
1139                                     rinv,
1140                                     phi0,
1141                                     d0,
1142                                     z0,
1143                                     t,
1144                                     rinvapprox,
1145                                     phi0approx,
1146                                     d0approx,
1147                                     z0approx,
1148                                     tapprox,
1149                                     irinv,
1150                                     iphi0,
1151                                     id0,
1152                                     iz0,
1153                                     it,
1154                                     projs,
1155                                     true);
1156
1157   if (settings_.debugTracklet())
1158     edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced " << getName()
1159                                  << " Found DDL tracklet in sector = " << iSector_ << " phi0 = " << phi0;
1160
1161   tracklet->setTrackletIndex(trackletpars_->nTracklets());
1162   tracklet->setTCIndex(TCIndex_);
1163
1164   if (settings_.writeMonitorData("Seeds")) {
1165     ofstream fout("seeds.txt", ofstream::app);
1166     fout << __FILE__ << ":" << __LINE__ << " " << name_ << "_" << iSector_ << " " << tracklet->getISeed() << endl;
1167     fout.close();
1168   }
1169   trackletpars_->addTracklet(tracklet);
1170
1171   for (unsigned int j = 0; j < toR_.size(); j++) {
1172     if (settings_.debugTracklet())
1173       edm::LogVerbatim("Tracklet") << "adding layer projection " << j << "/" << toR_.size() << " " << lproj_[j] << " "
1174                                    << tracklet->validProj(lproj_[j] - 1);
1175     if (tracklet->validProj(lproj_[j] - 1)) {
1176       addLayerProj(tracklet, lproj_[j]);
1177     }
1178   }
1179
1180   for (unsigned int j = 0; j < toZ_.size(); j++) {
1181     int disk = dproj_[j];
1182     if (disk == 0)
1183       continue;
1184     if (it < 0)
1185       disk = -disk;
1186     if (settings_.debugTracklet())
1187       edm::LogVerbatim("Tracklet") << "adding disk projection " << j << "/" << toZ_.size() << " " << disk << " "
1188                                    << tracklet->validProj(N_LAYER + abs(disk) - 1);
1189     if (tracklet->validProj(N_LAYER + abs(disk) - 1)) {
1190       addDiskProj(tracklet, disk);
1191     }
1192   }
1193
1194   return true;
1195 }
1196
1197 bool TrackletCalculatorDisplaced::LLDSeeding(const Stub* innerFPGAStub,
1198                                              const L1TStub* innerStub,
1199                                              const Stub* middleFPGAStub,
1200                                              const L1TStub* middleStub,
1201                                              const Stub* outerFPGAStub,
1202                                              const L1TStub* outerStub) {
1203   if (settings_.debugTracklet())
1204     edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced " << getName() << " " << layer_
1205                                  << " trying stub triplet in  (L2L3D1): " << middleFPGAStub->layer().value() << " "
1206                                  << outerFPGAStub->layer().value() << " " << innerFPGAStub->disk().value();
1207
1208   int take3 = 0;  //L2L3D1
1209   unsigned ndisks = 1;
1210
1211   // N.B. The names r1, r2, r3 reflect the fact that confusingly,
1212   // innerStub is actually the one furthest from the beamline ...
1213
1214   double r3 = innerStub->r();
1215   double z3 = innerStub->z();
1216   double phi3 = innerStub->phi();
1217
1218   double r1 = middleStub->r();
1219   double z1 = middleStub->z();
1220   double phi1 = middleStub->phi();
1221
1222   double r2 = outerStub->r();
1223   double z2 = outerStub->z();
1224   double phi2 = outerStub->phi();
1225
1226   double rinv, phi0, d0, t, z0;
1227
1228   double phiproj[N_LAYER - 2], zproj[N_LAYER - 2], phider[N_LAYER - 2], zder[N_LAYER - 2];
1229   double phiprojdisk[N_DISK], rprojdisk[N_DISK], phiderdisk[N_DISK], rderdisk[N_DISK];
1230
1231   exacttracklet(r1,
1232                 z1,
1233                 phi1,
1234                 r2,
1235                 z2,
1236                 phi2,
1237                 r3,
1238                 z3,
1239                 phi3,
1240                 take3,
1241                 rinv,
1242                 phi0,
1243                 d0,
1244                 t,
1245                 z0,
1246                 phiproj,
1247                 zproj,
1248                 phiprojdisk,
1249                 rprojdisk,
1250                 phider,
1251                 zder,
1252                 phiderdisk,
1253                 rderdisk);
1254
1255   if (settings_.debugTracklet())
1256     edm::LogVerbatim("Tracklet") << __LINE__ << ":" << __FILE__ << "LLD Exact values " << innerFPGAStub->isBarrel()
1257                                  << middleFPGAStub->isBarrel() << outerFPGAStub->isBarrel() << " " << phi3 << ", " << z3
1258                                  << ", " << r3 << ", " << phi1 << ", " << z1 << ", " << r1 << ", " << phi2 << ", " << z2
1259                                  << ", " << r2 << endl;
1260
1261   if (settings_.useapprox()) {
1262     phi3 = innerFPGAStub->phiapprox(phimin_, phimax_);
1263     z3 = innerFPGAStub->zapprox();
1264     r3 = innerFPGAStub->rapprox();
1265
1266     phi1 = middleFPGAStub->phiapprox(phimin_, phimax_);
1267     z1 = middleFPGAStub->zapprox();
1268     r1 = middleFPGAStub->rapprox();
1269
1270     phi2 = outerFPGAStub->phiapprox(phimin_, phimax_);
1271     z2 = outerFPGAStub->zapprox();
1272     r2 = outerFPGAStub->rapprox();
1273   }
1274
1275   if (settings_.debugTracklet())
1276     edm::LogVerbatim("Tracklet") << __LINE__ << ":" << __FILE__ << "LLD approx values " << innerFPGAStub->isBarrel()
1277                                  << middleFPGAStub->isBarrel() << outerFPGAStub->isBarrel() << " " << phi3 << ", " << z3
1278                                  << ", " << r3 << ", " << phi1 << ", " << z1 << ", " << r1 << ", " << phi2 << ", " << z2
1279                                  << ", " << r2 << endl;
1280
1281   double rinvapprox, phi0approx, d0approx, tapprox, z0approx;
1282   double phiprojapprox[N_LAYER - 2], zprojapprox[N_LAYER - 2], phiderapprox[N_LAYER - 2], zderapprox[N_LAYER - 2];
1283   double phiprojdiskapprox[N_DISK], rprojdiskapprox[N_DISK];
1284   double phiderdiskapprox[N_DISK], rderdiskapprox[N_DISK];
1285
1286   //TODO: implement the actual integer calculation
1287   if (settings_.useapprox()) {
1288     approxtracklet(r1,
1289                    z1,
1290                    phi1,
1291                    r2,
1292                    z2,
1293                    phi2,
1294                    r3,
1295                    z3,
1296                    phi3,
1297                    take3,
1298                    ndisks,
1299                    rinvapprox,
1300                    phi0approx,
1301                    d0approx,
1302                    tapprox,
1303                    z0approx,
1304                    phiprojapprox,
1305                    zprojapprox,
1306                    phiderapprox,
1307                    zderapprox,
1308                    phiprojdiskapprox,
1309                    rprojdiskapprox,
1310                    phiderdiskapprox,
1311                    rderdiskapprox);
1312   } else {
1313     rinvapprox = rinv;
1314     phi0approx = phi0;
1315     d0approx = d0;
1316     tapprox = t;
1317     z0approx = z0;
1318
1319     for (unsigned int i = 0; i < toR_.size(); ++i) {
1320       phiprojapprox[i] = phiproj[i];
1321       zprojapprox[i] = zproj[i];
1322       phiderapprox[i] = phider[i];
1323       zderapprox[i] = zder[i];
1324     }
1325
1326     for (unsigned int i = 0; i < toZ_.size(); ++i) {
1327       phiprojdiskapprox[i] = phiprojdisk[i];
1328       rprojdiskapprox[i] = rprojdisk[i];
1329       phiderdiskapprox[i] = phiderdisk[i];
1330       rderdiskapprox[i] = rderdisk[i];
1331     }
1332   }
1333
1334   //store the approcximate results
1335   if (settings_.debugTracklet()) {
1336     edm::LogVerbatim("Tracklet") << "rinvapprox: " << rinvapprox << " rinv: " << rinv << endl;
1337     edm::LogVerbatim("Tracklet") << "phi0approx: " << phi0approx << " phi0: " << phi0 << endl;
1338     edm::LogVerbatim("Tracklet") << "d0approx: " << d0approx << " d0: " << d0 << endl;
1339     edm::LogVerbatim("Tracklet") << "tapprox: " << tapprox << " t: " << t << endl;
1340     edm::LogVerbatim("Tracklet") << "z0approx: " << z0approx << " z0: " << z0 << endl;
1341   }
1342
1343   for (unsigned int i = 0; i < toR_.size(); ++i) {
1344     if (settings_.debugTracklet()) {
1345       edm::LogVerbatim("Tracklet") << "phiprojapprox[" << i << "]: " << phiprojapprox[i] << " phiproj[" << i
1346                                    << "]: " << phiproj[i] << endl;
1347       edm::LogVerbatim("Tracklet") << "zprojapprox[" << i << "]: " << zprojapprox[i] << " zproj[" << i
1348                                    << "]: " << zproj[i] << endl;
1349       edm::LogVerbatim("Tracklet") << "phiderapprox[" << i << "]: " << phiderapprox[i] << " phider[" << i
1350                                    << "]: " << phider[i] << endl;
1351       edm::LogVerbatim("Tracklet") << "zderapprox[" << i << "]: " << zderapprox[i] << " zder[" << i << "]: " << zder[i]
1352                                    << endl;
1353     }
1354   }
1355
1356   for (unsigned int i = 0; i < toZ_.size(); ++i) {
1357     if (settings_.debugTracklet()) {
1358       edm::LogVerbatim("Tracklet") << "phiprojdiskapprox[" << i << "]: " << phiprojdiskapprox[i] << " phiprojdisk[" << i
1359                                    << "]: " << phiprojdisk[i] << endl;
1360       edm::LogVerbatim("Tracklet") << "rprojdiskapprox[" << i << "]: " << rprojdiskapprox[i] << " rprojdisk[" << i
1361                                    << "]: " << rprojdisk[i] << endl;
1362       edm::LogVerbatim("Tracklet") << "phiderdiskapprox[" << i << "]: " << phiderdiskapprox[i] << " phiderdisk[" << i
1363                                    << "]: " << phiderdisk[i] << endl;
1364       edm::LogVerbatim("Tracklet") << "rderdiskapprox[" << i << "]: " << rderdiskapprox[i] << " rderdisk[" << i
1365                                    << "]: " << rderdisk[i] << endl;
1366     }
1367   }
1368
1369   //now binary
1370   double krinv = settings_.kphi1() / settings_.kr() * pow(2, settings_.rinv_shift()),
1371          kphi0 = settings_.kphi1() * pow(2, settings_.phi0_shift()),
1372          kt = settings_.kz() / settings_.kr() * pow(2, settings_.t_shift()),
1373          kz0 = settings_.kz() * pow(2, settings_.z0_shift()),
1374          kphiproj = settings_.kphi1() * pow(2, settings_.SS_phiL_shift()),
1375          kphider = settings_.kphi1() / settings_.kr() * pow(2, settings_.SS_phiderL_shift()),
1376          kzproj = settings_.kz() * pow(2, settings_.PS_zL_shift()),
1377          kzder = settings_.kz() / settings_.kr() * pow(2, settings_.PS_zderL_shift()),
1378          kphiprojdisk = settings_.kphi1() * pow(2, settings_.SS_phiD_shift()),
1379          kphiderdisk = settings_.kphi1() / settings_.kr() * pow(2, settings_.SS_phiderD_shift()),
1380          krprojdisk = settings_.kr() * pow(2, settings_.PS_rD_shift()),
1381          krderdisk = settings_.kr() / settings_.kz() * pow(2, settings_.PS_rderD_shift());
1382
1383   int irinv, iphi0, id0, it, iz0;
1384   int iphiproj[N_LAYER - 2], izproj[N_LAYER - 2], iphider[N_LAYER - 2], izder[N_LAYER - 2];
1385   int iphiprojdisk[N_DISK], irprojdisk[N_DISK], iphiderdisk[N_DISK], irderdisk[N_DISK];
1386
1387   //store the binary results
1388   irinv = rinvapprox / krinv;
1389   iphi0 = phi0approx / kphi0;
1390   id0 = d0approx / settings_.kd0();
1391   it = tapprox / kt;
1392   iz0 = z0approx / kz0;
1393
1394   bool success = true;
1395   if (std::abs(rinvapprox) > settings_.rinvcut()) {
1396     if (settings_.debugTracklet())
1397       edm::LogVerbatim("Tracklet") << "TrackletCalculator:: LLD Seeding irinv too large: " << rinvapprox << "(" << irinv
1398                                    << ")";
1399     success = false;
1400   }
1401   if (std::abs(z0approx) > settings_.disp_z0cut()) {
1402     if (settings_.debugTracklet())
1403       edm::LogVerbatim("Tracklet") << "Failed tracklet z0 cut " << z0approx;
1404     success = false;
1405   }
1406   if (std::abs(d0approx) > settings_.maxd0()) {
1407     if (settings_.debugTracklet())
1408       edm::LogVerbatim("Tracklet") << "Failed tracklet approx d0 cut " << d0approx;
1409     success = false;
1410   }
1411   if (std::abs(d0) > settings_.maxd0()) {
1412     if (settings_.debugTracklet())
1413       edm::LogVerbatim("Tracklet") << "Failed tracklet exact d0 cut " << d0;
1414     success = false;
1415   }
1416
1417   if (!success)
1418     return false;
1419
1420   double phicritapprox = phi0approx - asin((0.5 * settings_.rcrit() * rinvapprox) + (d0approx / settings_.rcrit()));
1421   int phicrit = iphi0 - 2 * irinv - 2 * id0;
1422
1423   int iphicritmincut = settings_.phicritminmc() / globals_->ITC_L1L2()->phi0_final.K();
1424   int iphicritmaxcut = settings_.phicritmaxmc() / globals_->ITC_L1L2()->phi0_final.K();
1425
1426   bool keepapprox = (phicritapprox > settings_.phicritminmc()) && (phicritapprox < settings_.phicritmaxmc()),
1427        keep = (phicrit > iphicritmincut) && (phicrit < iphicritmaxcut);
1428
1429   if (settings_.debugTracklet())
1430     if (keep && !keepapprox)
1431       edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced::LLDSeeding tracklet kept with exact phicrit cut "
1432                                       "but not approximate, phicritapprox: "
1433                                    << phicritapprox;
1434   if (settings_.usephicritapprox()) {
1435     if (!keepapprox)
1436       return false;
1437   } else {
1438     if (!keep)
1439       return false;
1440   }
1441
1442   Projection projs[N_LAYER + N_DISK];
1443
1444   for (unsigned int i = 0; i < toR_.size(); ++i) {
1445     iphiproj[i] = phiprojapprox[i] / kphiproj;
1446     izproj[i] = zprojapprox[i] / kzproj;
1447
1448     iphider[i] = phiderapprox[i] / kphider;
1449     izder[i] = zderapprox[i] / kzder;
1450
1451     if (izproj[i] < -(1 << (settings_.nzbitsstub(0) - 1)))
1452       continue;
1453     if (izproj[i] >= (1 << (settings_.nzbitsstub(0) - 1)))
1454       continue;
1455
1456     //this is left from the original....
1457     if (iphiproj[i] >= (1 << settings_.nphibitsstub(5)) - 1)
1458       continue;
1459     if (iphiproj[i] <= 0)
1460       continue;
1461
1462     if (rproj_[i] < settings_.rPS2S()) {
1463       iphiproj[i] >>= (settings_.nphibitsstub(5) - settings_.nphibitsstub(0));
1464     } else {
1465       izproj[i] >>= (settings_.nzbitsstub(0) - settings_.nzbitsstub(5));
1466     }
1467
1468     if (rproj_[i] < settings_.rPS2S()) {
1469       if (iphider[i] < -(1 << (settings_.nbitsphiprojderL123() - 1)))
1470         iphider[i] = -(1 << (settings_.nbitsphiprojderL123() - 1));
1471       if (iphider[i] >= (1 << (settings_.nbitsphiprojderL123() - 1)))
1472         iphider[i] = (1 << (settings_.nbitsphiprojderL123() - 1)) - 1;
1473     } else {
1474       if (iphider[i] < -(1 << (settings_.nbitsphiprojderL456() - 1)))
1475         iphider[i] = -(1 << (settings_.nbitsphiprojderL456() - 1));
1476       if (iphider[i] >= (1 << (settings_.nbitsphiprojderL456() - 1)))
1477         iphider[i] = (1 << (settings_.nbitsphiprojderL456() - 1)) - 1;
1478     }
1479
1480     projs[lproj_[i] - 1].init(settings_,
1481                               lproj_[i] - 1,
1482                               iphiproj[i],
1483                               izproj[i],
1484                               iphider[i],
1485                               izder[i],
1486                               phiproj[i],
1487                               zproj[i],
1488                               phider[i],
1489                               zder[i],
1490                               phiprojapprox[i],
1491                               zprojapprox[i],
1492                               phiderapprox[i],
1493                               zderapprox[i],
1494                               false);
1495   }
1496
1497   if (std::abs(it * kt) > 1.0) {
1498     for (unsigned int i = 0; i < toZ_.size(); ++i) {
1499       iphiprojdisk[i] = phiprojdiskapprox[i] / kphiprojdisk;
1500       irprojdisk[i] = rprojdiskapprox[i] / krprojdisk;
1501
1502       iphiderdisk[i] = phiderdiskapprox[i] / kphiderdisk;
1503       irderdisk[i] = rderdiskapprox[i] / krderdisk;
1504
1505       //Check phi range of projection
1506       if (iphiprojdisk[i] <= 0)
1507         continue;
1508       if (iphiprojdisk[i] >= (1 << settings_.nphibitsstub(0)) - 1)
1509         continue;
1510
1511       //Check r range of projection
1512       if (irprojdisk[i] < settings_.rmindisk() / krprojdisk || irprojdisk[i] >= settings_.rmaxdisk() / krprojdisk)
1513         continue;
1514
1515       projs[N_LAYER + i + 1].init(settings_,
1516                                   N_LAYER + i + 1,
1517                                   iphiprojdisk[i],
1518                                   irprojdisk[i],
1519                                   iphiderdisk[i],
1520                                   irderdisk[i],
1521                                   phiprojdisk[i],
1522                                   rprojdisk[i],
1523                                   phiderdisk[i],
1524                                   rderdisk[i],
1525                                   phiprojdiskapprox[i],
1526                                   rprojdiskapprox[i],
1527                                   phiderdisk[i],
1528                                   rderdisk[i],
1529                                   false);
1530     }
1531   }
1532
1533   if (settings_.writeMonitorData("TrackletPars")) {
1534     globals_->ofstream("trackletpars.txt")
1535         << layer_ << " , " << rinv << " , " << rinvapprox << " , " << phi0 << " , " << phi0approx << " , " << t << " , "
1536         << tapprox << " , " << z0 << " , " << z0approx << " , " << d0 << " , " << d0approx << endl;
1537   }
1538
1539   Tracklet* tracklet = new Tracklet(settings_,
1540                                     iSeed_,
1541                                     innerFPGAStub,
1542                                     middleFPGAStub,
1543                                     outerFPGAStub,
1544                                     rinv,
1545                                     phi0,
1546                                     d0,
1547                                     z0,
1548                                     t,
1549                                     rinvapprox,
1550                                     phi0approx,
1551                                     d0approx,
1552                                     z0approx,
1553                                     tapprox,
1554                                     irinv,
1555                                     iphi0,
1556                                     id0,
1557                                     iz0,
1558                                     it,
1559                                     projs,
1560                                     false);
1561
1562   if (settings_.debugTracklet())
1563     edm::LogVerbatim("Tracklet") << "TrackletCalculatorDisplaced " << getName()
1564                                  << " Found LLD tracklet in sector = " << iSector_ << " phi0 = " << phi0;
1565
1566   tracklet->setTrackletIndex(trackletpars_->nTracklets());
1567   tracklet->setTCIndex(TCIndex_);
1568
1569   if (settings_.writeMonitorData("Seeds")) {
1570     ofstream fout("seeds.txt", ofstream::app);
1571     fout << __FILE__ << ":" << __LINE__ << " " << name_ << "_" << iSector_ << " " << tracklet->getISeed() << endl;
1572     fout.close();
1573   }
1574   trackletpars_->addTracklet(tracklet);
1575
1576   for (unsigned int j = 0; j < toR_.size(); j++) {
1577     if (settings_.debugTracklet())
1578       edm::LogVerbatim("Tracklet") << "adding layer projection " << j << "/" << toR_.size() << " " << lproj_[j];
1579     if (tracklet->validProj(lproj_[j] - 1)) {
1580       addLayerProj(tracklet, lproj_[j]);
1581     }
1582   }
1583
1584   for (unsigned int j = 0; j < toZ_.size(); j++) {
1585     int disk = dproj_[j];
1586     if (disk == 0)
1587       continue;
1588     if (it < 0)
1589       disk = -disk;
1590     if (settings_.debugTracklet())
1591       edm::LogVerbatim("Tracklet") << "adding disk projection " << j << "/" << toZ_.size() << " " << disk;
1592     if (tracklet->validProj(N_LAYER + abs(disk) - 1)) {
1593       addDiskProj(tracklet, disk);
1594     }
1595   }
1596
1597   return true;
1598 }
1599
1600 void TrackletCalculatorDisplaced::exactproj(double rproj,
1601                                             double rinv,
1602                                             double phi0,
1603                                             double d0,
1604                                             double t,
1605                                             double z0,
1606                                             double r0,
1607                                             double& phiproj,
1608                                             double& zproj,
1609                                             double& phider,
1610                                             double& zder) {
1611   double rho = 1 / rinv;
1612   if (rho < 0) {
1613     r0 = -r0;
1614   }
1615   phiproj = phi0 - asin((rproj * rproj + r0 * r0 - rho * rho) / (2 * rproj * r0));
1616   double beta = acos((rho * rho + r0 * r0 - rproj * rproj) / (2 * r0 * rho));
1617   zproj = z0 + t * std::abs(rho * beta);
1618
1619   //not exact, but close
1620   phider = -0.5 * rinv / sqrt(1 - pow(0.5 * rproj * rinv, 2)) + d0 / (rproj * rproj);
1621   zder = t / sqrt(1 - pow(0.5 * rproj * rinv, 2));
1622
1623   if (settings_.debugTracklet())
1624     edm::LogVerbatim("Tracklet") << "exact proj layer at " << rproj << " : " << phiproj << " " << zproj;
1625 }
1626
1627 void TrackletCalculatorDisplaced::exactprojdisk(double zproj,
1628                                                 double rinv,
1629                                                 double,
1630                                                 double,  //phi0 and d0 are not used.
1631                                                 double t,
1632                                                 double z0,
1633                                                 double x0,
1634                                                 double y0,
1635                                                 double& phiproj,
1636                                                 double& rproj,
1637                                                 double& phider,
1638                                                 double& rder) {
1639   //protect against t=0
1640   if (std::abs(t) < 0.1)
1641     t = 0.1;
1642   if (t < 0)
1643     zproj = -zproj;
1644   double rho = std::abs(1 / rinv);
1645   double beta = (zproj - z0) / (t * rho);
1646   double phiV = atan2(-y0, -x0);
1647   double c = rinv > 0 ? -1 : 1;
1648
1649   double x = x0 + rho * cos(phiV + c * beta);
1650   double y = y0 + rho * sin(phiV + c * beta);
1651
1652   phiproj = atan2(y, x);
1653
1654   phiproj = reco::reducePhiRange(phiproj - phimin_);
1655
1656   rproj = sqrt(x * x + y * y);
1657
1658   phider = c / t / (x * x + y * y) * (rho + x0 * cos(phiV + c * beta) + y0 * sin(phiV + c * beta));
1659   rder = c / t / rproj * (y0 * cos(phiV + c * beta) - x0 * sin(phiV + c * beta));
1660
1661   if (settings_.debugTracklet())
1662     edm::LogVerbatim("Tracklet") << "exact proj disk at" << zproj << " : " << phiproj << " " << rproj;
1663 }
1664
1665 void TrackletCalculatorDisplaced::exacttracklet(double r1,
1666                                                 double z1,
1667                                                 double phi1,
1668                                                 double r2,
1669                                                 double z2,
1670                                                 double phi2,
1671                                                 double r3,
1672                                                 double z3,
1673                                                 double phi3,
1674                                                 int take3,
1675                                                 double& rinv,
1676                                                 double& phi0,
1677                                                 double& d0,
1678                                                 double& t,
1679                                                 double& z0,
1680                                                 double phiproj[N_LAYER - 2],
1681                                                 double zproj[N_LAYER - 2],
1682                                                 double phiprojdisk[N_DISK],
1683                                                 double rprojdisk[N_DISK],
1684                                                 double phider[N_LAYER - 2],
1685                                                 double zder[N_LAYER - 2],
1686                                                 double phiderdisk[N_DISK],
1687                                                 double rderdisk[N_DISK]) {
1688   //two lines perpendicular to the 1->2 and 2->3
1689   double x1 = r1 * cos(phi1);
1690   double x2 = r2 * cos(phi2);
1691   double x3 = r3 * cos(phi3);
1692
1693   double y1 = r1 * sin(phi1);
1694   double y2 = r2 * sin(phi2);
1695   double y3 = r3 * sin(phi3);
1696
1697   double dy21 = y2 - y1;
1698   double dy32 = y3 - y2;
1699
1700   //Hack to protect against dividing by zero
1701   //code should be rewritten to avoid this
1702   if (dy21 == 0.0)
1703     dy21 = 1e-9;
1704   if (dy32 == 0.0)
1705     dy32 = 1e-9;
1706
1707   double k1 = -(x2 - x1) / dy21;
1708   double k2 = -(x3 - x2) / dy32;
1709   double b1 = 0.5 * (y2 + y1) - 0.5 * (x1 + x2) * k1;
1710   double b2 = 0.5 * (y3 + y2) - 0.5 * (x2 + x3) * k2;
1711   //their intersection gives the center of the circle
1712   double y0 = (b1 * k2 - b2 * k1) / (k2 - k1);
1713   double x0 = (b1 - b2) / (k2 - k1);
1714   //get the radius three ways:
1715   double R1 = sqrt(pow(x1 - x0, 2) + pow(y1 - y0, 2));
1716   double R2 = sqrt(pow(x2 - x0, 2) + pow(y2 - y0, 2));
1717   double R3 = sqrt(pow(x3 - x0, 2) + pow(y3 - y0, 2));
1718   //check if the same
1719   double eps1 = std::abs(R1 / R2 - 1);
1720   double eps2 = std::abs(R3 / R2 - 1);
1721   if (eps1 > 1e-10 || eps2 > 1e-10)
1722     edm::LogVerbatim("Tracklet") << "&&&&&&&&&&&& bad circle! " << R1 << "\t" << R2 << "\t" << R3;
1723
1724   if (settings_.debugTracklet())
1725     edm::LogVerbatim("Tracklet") << "phimin_: " << phimin_ << " phimax_: " << phimax_;
1726   //results
1727   rinv = 1. / R1;
1728   phi0 = 0.5 * M_PI + atan2(y0, x0);
1729
1730   phi0 -= phimin_;
1731
1732   d0 = -R1 + sqrt(x0 * x0 + y0 * y0);
1733   //sign of rinv:
1734   double dphi = reco::reducePhiRange(phi3 - atan2(y0, x0));
1735   if (dphi < 0) {
1736     rinv = -rinv;
1737     d0 = -d0;
1738     phi0 = phi0 + M_PI;
1739   }
1740   phi0 = angle0to2pi::make0To2pi(phi0);
1741
1742   //now in RZ:
1743   //turning angle
1744   double beta1 = reco::reducePhiRange(atan2(y1 - y0, x1 - x0) - atan2(-y0, -x0));
1745   double beta2 = reco::reducePhiRange(atan2(y2 - y0, x2 - x0) - atan2(-y0, -x0));
1746   double beta3 = reco::reducePhiRange(atan2(y3 - y0, x3 - x0) - atan2(-y0, -x0));
1747
1748   double t12 = (z2 - z1) / std::abs(beta2 - beta1) / R1;
1749   double z12 = (z1 * beta2 - z2 * beta1) / (beta2 - beta1);
1750   double t13 = (z3 - z1) / std::abs(beta3 - beta1) / R1;
1751   double z13 = (z1 * beta3 - z3 * beta1) / (beta3 - beta1);
1752
1753   if (take3 > 0) {
1754     //take 13 (large lever arm)
1755     t = t13;
1756     z0 = z13;
1757   } else {
1758     //take 12 (pixel layers)
1759     t = t12;
1760     z0 = z12;
1761   }
1762
1763   for (unsigned int i = 0; i < toR_.size(); i++) {
1764     exactproj(toR_[i], rinv, phi0, d0, t, z0, sqrt(x0 * x0 + y0 * y0), phiproj[i], zproj[i], phider[i], zder[i]);
1765   }
1766
1767   for (unsigned int i = 0; i < toZ_.size(); i++) {
1768     exactprojdisk(toZ_[i], rinv, phi0, d0, t, z0, x0, y0, phiprojdisk[i], rprojdisk[i], phiderdisk[i], rderdisk[i]);
1769   }
1770
1771   if (settings_.debugTracklet())
1772     edm::LogVerbatim("Tracklet") << "exact tracklet: " << rinv << " " << phi0 << " " << t << " " << z0 << " " << d0;
1773 }
1774
1775 void TrackletCalculatorDisplaced::approxproj(double halfRinv,
1776                                              double phi0,
1777                                              double d0,
1778                                              double t,
1779                                              double z0,
1780                                              double halfRinv_0,
1781                                              double d0_0,  // zeroeth order result for higher order terms calculation
1782                                              double rmean,
1783                                              double& phiproj,
1784                                              double& phiprojder,
1785                                              double& zproj,
1786                                              double& zprojder) {
1787   if (std::abs(2.0 * halfRinv) > settings_.rinvcut() || std::abs(z0) > settings_.disp_z0cut() ||
1788       std::abs(d0) > settings_.maxd0()) {
1789     phiproj = 0.0;
1790     return;
1791   }
1792   double rmeanInv = 1.0 / rmean;
1793
1794   phiproj = phi0 + rmean * (-halfRinv + 2.0 * d0_0 * halfRinv_0 * halfRinv_0) +
1795             rmeanInv * (-d0 + halfRinv_0 * d0_0 * d0_0) + sixth * pow(-rmean * halfRinv_0 - rmeanInv * d0_0, 3);
1796   phiprojder = -halfRinv + d0 * rmeanInv * rmeanInv;  //removed all high terms
1797
1798   zproj = z0 + t * rmean - 0.5 * rmeanInv * t * d0_0 * d0_0 - t * rmean * halfRinv * d0 +
1799           sixth * pow(rmean, 3) * t * halfRinv_0 * halfRinv_0;
1800   zprojder = t;  // removed all high terms
1801
1802   phiproj = angle0to2pi::make0To2pi(phiproj);
1803
1804   if (settings_.debugTracklet())
1805     edm::LogVerbatim("Tracklet") << "approx proj layer at " << rmean << " : " << phiproj << " " << zproj << endl;
1806 }
1807
1808 void TrackletCalculatorDisplaced::approxprojdisk(double halfRinv,
1809                                                  double phi0,
1810                                                  double d0,
1811                                                  double t,
1812                                                  double z0,
1813                                                  double halfRinv_0,
1814                                                  double d0_0,  // zeroeth order result for higher order terms calculation
1815                                                  double zmean,
1816                                                  double& phiproj,
1817                                                  double& phiprojder,
1818                                                  double& rproj,
1819                                                  double& rprojder) {
1820   if (std::abs(2.0 * halfRinv) > settings_.rinvcut() || std::abs(z0) > settings_.disp_z0cut() ||
1821       std::abs(d0) > settings_.maxd0()) {
1822     phiproj = 0.0;
1823     return;
1824   }
1825
1826   if (t < 0)
1827     zmean = -zmean;
1828
1829   double zmeanInv = 1.0 / zmean, rstar = (zmean - z0) / t,
1830          epsilon = 0.5 * zmeanInv * zmeanInv * d0_0 * d0_0 * t * t + halfRinv * d0 -
1831                    sixth * rstar * rstar * halfRinv_0 * halfRinv_0;
1832
1833   rproj = rstar * (1 + epsilon);
1834   rprojder = 1 / t;
1835
1836   double A = rproj * halfRinv;
1837   double B = -d0 * t * zmeanInv * (1 + z0 * zmeanInv) * (1 - epsilon);
1838   double C = -d0 * halfRinv;
1839   double A_0 = rproj * halfRinv_0;
1840   double B_0 = -d0_0 * t * zmeanInv * (1 + z0 * zmeanInv) * (1 - epsilon);
1841   // double C_0 = -d0_0 * halfRinv_0;
1842
1843   phiproj = phi0 - A + B * (1 + C - 2 * A_0 * A_0) + sixth * pow(-A_0 + B_0, 3);
1844   phiprojder = -halfRinv / t + d0 * t * zmeanInv * zmeanInv;
1845
1846   phiproj = angle0to2pi::make0To2pi(phiproj);
1847
1848   if (settings_.debugTracklet())
1849     edm::LogVerbatim("Tracklet") << "approx proj disk at" << zmean << " : " << phiproj << " " << rproj << endl;
1850 }
1851
1852 void TrackletCalculatorDisplaced::approxtracklet(double r1,
1853                                                  double z1,
1854                                                  double phi1,
1855                                                  double r2,
1856                                                  double z2,
1857                                                  double phi2,
1858                                                  double r3,
1859                                                  double z3,
1860                                                  double phi3,
1861                                                  bool take3,
1862                                                  unsigned ndisks,
1863                                                  double& rinv,
1864                                                  double& phi0,
1865                                                  double& d0,
1866                                                  double& t,
1867                                                  double& z0,
1868                                                  double phiproj[4],
1869                                                  double zproj[4],
1870                                                  double phider[4],
1871                                                  double zder[4],
1872                                                  double phiprojdisk[5],
1873                                                  double rprojdisk[5],
1874                                                  double phiderdisk[5],
1875                                                  double rderdisk[5]) {
1876   double a = 1.0 / ((r1 - r2) * (r1 - r3));
1877   double b = 1.0 / ((r1 - r2) * (r2 - r3));
1878   double c = 1.0 / ((r1 - r3) * (r2 - r3));
1879
1880   // first iteration in r-phi plane
1881   double halfRinv_0 = -phi1 * r1 * a + phi2 * r2 * b - phi3 * r3 * c;
1882   double d0_0 = r1 * r2 * r3 * (-phi1 * a + phi2 * b - phi3 * c);
1883
1884   // corrections to phi1, phi2, and phi3
1885   double r = r2, z = z2;
1886   if (take3)
1887     r = r3, z = z3;
1888
1889   double d0OverR1 = d0_0 * rzmeanInv_[0] * (ndisks > 2 ? std::abs((z - z1) / (r - r1)) : 1.0);
1890   double d0OverR2 = d0_0 * rzmeanInv_[1] * (ndisks > 1 ? std::abs((z - z1) / (r - r1)) : 1.0);
1891   double d0OverR3 = d0_0 * rzmeanInv_[2] * (ndisks > 0 ? std::abs((z - z1) / (r - r1)) : 1.0);
1892
1893   double d0OverR = d0OverR2;
1894   if (take3)
1895     d0OverR = d0OverR3;
1896
1897   double c1 = d0_0 * halfRinv_0 * d0OverR1 + 2.0 * d0_0 * halfRinv_0 * r1 * halfRinv_0 +
1898               sixth * pow(-r1 * halfRinv_0 - d0OverR1, 3);
1899   double c2 = d0_0 * halfRinv_0 * d0OverR2 + 2.0 * d0_0 * halfRinv_0 * r2 * halfRinv_0 +
1900               sixth * pow(-r2 * halfRinv_0 - d0OverR2, 3);
1901   double c3 = d0_0 * halfRinv_0 * d0OverR3 + 2.0 * d0_0 * halfRinv_0 * r3 * halfRinv_0 +
1902               sixth * pow(-r3 * halfRinv_0 - d0OverR3, 3);
1903
1904   double phi1c = phi1 - c1;
1905   double phi2c = phi2 - c2;
1906   double phi3c = phi3 - c3;
1907
1908   // second iteration in r-phi plane
1909   double halfRinv = -phi1c * r1 * a + phi2c * r2 * b - phi3c * r3 * c;
1910   phi0 = -phi1c * r1 * (r2 + r3) * a + phi2c * r2 * (r1 + r3) * b - phi3c * r3 * (r1 + r2) * c;
1911   d0 = r1 * r2 * r3 * (-phi1c * a + phi2c * b - phi3c * c);
1912
1913   t = ((z - z1) / (r - r1)) *
1914       (1. + d0 * halfRinv - 0.5 * d0OverR1 * d0OverR - sixth * (r1 * r1 + r2 * r2 + r1 * r2) * halfRinv_0 * halfRinv_0);
1915   z0 = z1 - t * r1 * (1.0 - d0_0 * halfRinv_0 - 0.5 * d0OverR1 * d0OverR1 + sixth * r1 * r1 * halfRinv_0 * halfRinv_0);
1916
1917   rinv = 2.0 * halfRinv;
1918   phi0 += -phimin_;
1919
1920   phi0 = angle0to2pi::make0To2pi(phi0);
1921
1922   for (unsigned int i = 0; i < toR_.size(); i++) {
1923     approxproj(halfRinv,
1924                phi0,
1925                d0,
1926                t,
1927                z0,
1928                halfRinv_0,
1929                d0_0,  // added _0 version for high term calculations
1930                toR_.at(i),
1931                phiproj[i],
1932                phider[i],
1933                zproj[i],
1934                zder[i]);
1935   }
1936
1937   for (unsigned int i = 0; i < toZ_.size(); i++) {
1938     approxprojdisk(halfRinv,
1939                    phi0,
1940                    d0,
1941                    t,
1942                    z0,
1943                    halfRinv_0,
1944                    d0_0,  // added _0 version for high term calculations
1945                    toZ_.at(i),
1946                    phiprojdisk[i],
1947                    phiderdisk[i],
1948                    rprojdisk[i],
1949                    rderdisk[i]);
1950   }
1951
1952   if (std::abs(rinv) > settings_.rinvcut() || std::abs(z0) > settings_.disp_z0cut() ||
1953       std::abs(d0) > settings_.maxd0()) {
1954     phi0 = 0.0;
1955     return;
1956   }
1957
1958   if (settings_.debugTracklet())
1959     edm::LogVerbatim("Tracklet") << "TCD approx tracklet: " << rinv << " " << phi0 << " " << t << " " << z0 << " " << d0
1960                                  << endl;
1961 }