TrackFindingTracklet/src/PurgeDuplicate.cc

0001 #include "L1Trigger/TrackFindingTracklet/interface/PurgeDuplicate.h"
0002 #include "L1Trigger/TrackFindingTracklet/interface/Tracklet.h"
0003 #include "L1Trigger/TrackFindingTracklet/interface/Globals.h"
0004 #include "L1Trigger/TrackFindingTracklet/interface/CleanTrackMemory.h"
0005 #include "L1Trigger/TrackFindingTracklet/interface/Settings.h"
0006 #include "L1Trigger/TrackFindingTracklet/interface/Stub.h"
0007 #include "L1Trigger/TrackFindingTracklet/interface/Track.h"
0008
0009 #ifdef USEHYBRID
0010 #include "DataFormats/L1TrackTrigger/interface/TTStub.h"
0011 #include "L1Trigger/TrackFindingTMTT/interface/L1track3D.h"
0012 #include "L1Trigger/TrackFindingTMTT/interface/KFParamsComb.h"
0013 #include "L1Trigger/TrackFindingTracklet/interface/HybridFit.h"
0014 #endif
0015
0016 #include "DataFormats/Math/interface/deltaPhi.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/Utilities/interface/Exception.h"
0019
0020 #include "TString.h"
0021 #include <unordered_set>
0022 #include <algorithm>
0023
0024 using namespace std;
0025 using namespace trklet;
0026
0027 PurgeDuplicate::PurgeDuplicate(std::string name, Settings const& settings, Globals* global)
0028     : ProcessBase(name, settings, global) {}
0029
0030 void PurgeDuplicate::addOutput(MemoryBase* memory, std::string output) {
0031   if (settings_.writetrace()) {
0032     edm::LogVerbatim("Tracklet") << "In " << name_ << " adding output to " << memory->getName() << " to output "
0033                                  << output;
0034   }
0035   unordered_set<string> outputs = {"trackout",
0036                                    "trackout1",
0037                                    "trackout2",
0038                                    "trackout3",
0039                                    "trackout4",
0040                                    "trackout5",
0041                                    "trackout6",
0042                                    "trackout7",
0043                                    "trackout8",
0044                                    "trackout9",
0045                                    "trackout10",
0046                                    "trackout11"};
0047   if (outputs.find(output) != outputs.end()) {
0048     auto* tmp = dynamic_cast<CleanTrackMemory*>(memory);
0049     assert(tmp != nullptr);
0050     outputtracklets_.push_back(tmp);
0051     return;
0052   }
0053   throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " could not find output: " << output;
0054 }
0055
0056 void PurgeDuplicate::addInput(MemoryBase* memory, std::string input) {
0057   if (settings_.writetrace()) {
0058     edm::LogVerbatim("Tracklet") << "In " << name_ << " adding input from " << memory->getName() << " to input "
0059                                  << input;
0060   }
0061   unordered_set<string> inputs = {"trackin",
0062                                   "trackin1",
0063                                   "trackin2",
0064                                   "trackin3",
0065                                   "trackin4",
0066                                   "trackin5",
0067                                   "trackin6",
0068                                   "trackin7",
0069                                   "trackin8",
0070                                   "trackin9",
0071                                   "trackin10",
0072                                   "trackin11",
0073                                   "trackin12"};
0074   if (inputs.find(input) != inputs.end()) {
0075     auto* tmp = dynamic_cast<TrackFitMemory*>(memory);
0076     assert(tmp != nullptr);
0077     inputtrackfits_.push_back(tmp);
0078     return;
0079   }
0080   throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " could not find input: " << input;
0081 }
0082
0083 void PurgeDuplicate::execute(std::vector<Track>& outputtracks, unsigned int iSector) {
0084   inputtracklets_.clear();
0085   inputtracks_.clear();
0086
0087   inputstubidslists_.clear();
0088   inputstublists_.clear();
0089   mergedstubidslists_.clear();
0090
0091   if (settings_.removalType() != "merge") {
0092     for (auto& inputtrackfit : inputtrackfits_) {
0093       if (inputtrackfit->nTracks() == 0)
0094         continue;
0095       for (unsigned int j = 0; j < inputtrackfit->nTracks(); j++) {
0096         Track* aTrack = inputtrackfit->getTrack(j)->getTrack();
0097         aTrack->setSector(iSector);
0098         inputtracks_.push_back(aTrack);
0099       }
0100     }
0101     if (inputtracks_.empty())
0102       return;
0103   }
0104
0105   unsigned int numTrk = inputtracks_.size();
0106
0107   ////////////////////
0108   // Hybrid Removal //
0109   ////////////////////
0110 #ifdef USEHYBRID
0111
0112   if (settings_.removalType() == "merge") {
0113     // Track seed & duplicate flag
0114     std::vector<std::pair<int, bool>> trackInfo;
0115     // Flag for tracks in multiple bins that get merged but are not in the correct bin
0116     std::vector<bool> trackBinInfo;
0117     // Vector to store the relative rank of the track candidate for merging, based on seed type
0118     std::vector<int> seedRank;
0119
0120     // Stubs on every track
0121     std::vector<std::vector<const Stub*>> inputstublistsall;
0122     // (layer, unique stub index within layer) of each stub on every track
0123     std::vector<std::vector<std::pair<int, int>>> mergedstubidslistsall;
0124     std::vector<std::vector<std::pair<int, int>>> inputstubidslistsall;
0125     std::vector<Tracklet*> inputtrackletsall;
0126
0127     std::vector<unsigned int> prefTracks;  // Stores all the tracks that are sent to the KF from each bin
0128     std::vector<int> prefTrackFit;  // Stores the track seed that corresponds to the associated track in prefTracks
0129
0130     for (unsigned int bin = 0; bin < settings_.rinvBins().size() - 1; bin++) {
0131       for (unsigned int phiBin = 0; phiBin < settings_.phiBins().size() - 1; phiBin++) {
0132         // Get vectors from TrackFit and save them
0133         // inputtracklets: Tracklet objects from the FitTrack (not actually fit yet)
0134         // inputstublists: L1Stubs for that track
0135         // inputstubidslists: Stub stubIDs for that 3rack
0136         // mergedstubidslists: the same as inputstubidslists, but will be used during duplicate removal
0137         for (unsigned int i = 0; i < inputtrackfits_.size(); i++) {
0138           if (inputtrackfits_[i]->nStublists() == 0)
0139             continue;
0140           if (inputtrackfits_[i]->nStublists() != inputtrackfits_[i]->nTracks())
0141             throw cms::Exception("LogicError")
0142                 << __FILE__ << " " << __LINE__ << " Number of stublists and tracks don't match up! ";
0143           for (unsigned int j = 0; j < inputtrackfits_[i]->nStublists(); j++) {
0144             if (isTrackInBin(findOverlapRinvBins(inputtrackfits_[i]->getTrack(j)), bin)) {
0145               if (!isTrackInBin(findOverlapPhiBins(inputtrackfits_[i]->getTrack(j)), phiBin))
0146                 continue;
0147               if (inputtracklets_.size() >= settings_.maxStep("DR"))
0148                 continue;
0149               Tracklet* aTrack = inputtrackfits_[i]->getTrack(j);
0150               inputtracklets_.push_back(inputtrackfits_[i]->getTrack(j));
0151               std::vector<const Stub*> stublist = inputtrackfits_[i]->getStublist(j);
0152               inputstublists_.push_back(stublist);
0153               std::vector<std::pair<int, int>> stubidslist = inputtrackfits_[i]->getStubidslist(j);
0154               inputstubidslists_.push_back(stubidslist);
0155               mergedstubidslists_.push_back(stubidslist);
0156
0157               // Encoding: L1L2=0, L2L3=1, L3L4=2, L5L6=3, D1D2=4, D3D4=5, L1D1=6, L2D1=7
0158               // Best Guess:          L1L2 > L1D1 > L2L3 > L2D1 > D1D2 > L3L4 > L5L6 > D3D4
0159               // Best Rank:           L1L2 > L3L4 > D3D4 > D1D2 > L2L3 > L2D1 > L5L6 > L1D1
0160               // Rank-Informed Guess: L1L2 > L3L4 > L1D1 > L2L3 > L2D1 > D1D2 > L5L6 > D3D4
0161               const unsigned int curSeed = aTrack->seedIndex();
0162               static const std::vector<int> ranks{1, 5, 2, 7, 4, 3, 8, 6};
0163               if (curSeed < ranks.size()) {
0164                 seedRank.push_back(ranks[curSeed]);
0165               } else if (settings_.extended()) {
0166                 seedRank.push_back(9);
0167               } else {
0168                 throw cms::Exception("LogError") << __FILE__ << " " << __LINE__ << " Seed type " << curSeed
0169                                                  << " not found in list, and settings->extended() not set.";
0170               }
0171
0172               if (stublist.size() != stubidslist.size())
0173                 throw cms::Exception("LogicError")
0174                     << __FILE__ << " " << __LINE__ << " Number of stubs and stubids don't match up! ";
0175
0176               trackInfo.emplace_back(i, false);
0177               trackBinInfo.emplace_back(false);
0178             } else
0179               continue;
0180           }
0181         }
0182
0183         if (inputtracklets_.empty())
0184           continue;
0185         const unsigned int numStublists = inputstublists_.size();
0186
0187         if (settings_.inventStubs()) {
0188           for (unsigned int itrk = 0; itrk < numStublists; itrk++) {
0189             inputstublists_[itrk] = getInventedSeedingStub(iSector, inputtracklets_[itrk], inputstublists_[itrk]);
0190           }
0191         }
0192
0193         // Initialize all-false 2D vector of tracks being duplicates to other tracks
0194         vector<vector<bool>> dupMap(numStublists, vector<bool>(numStublists, false));
0195
0196         // Used to check if a track is in two bins, is not a duplicate in either bin, so is sent out twice
0197         vector<bool> noMerge(numStublists, false);
0198
0199         // Find duplicates; Fill dupMap by looping over all pairs of "tracks"
0200         // numStublists-1 since last track has no other to compare to
0201         for (unsigned int itrk = 0; itrk < numStublists - 1; itrk++) {
0202           for (unsigned int jtrk = itrk + 1; jtrk < numStublists; jtrk++) {
0203             if (itrk >= settings_.numTracksComparedPerBin())
0204               continue;
0205             // Get primary track stubids = (layer, unique stub index within layer)
0206             const std::vector<std::pair<int, int>>& stubsTrk1 = inputstubidslists_[itrk];
0207
0208             // Get and count secondary track stubids
0209             const std::vector<std::pair<int, int>>& stubsTrk2 = inputstubidslists_[jtrk];
0210
0211             // Count number of layers that share stubs, and the number of UR that each track hits
0212             unsigned int nShareLay = 0;
0213             if (settings_.mergeComparison() == "CompareAll") {
0214               bool layerArr[16];
0215               for (auto& i : layerArr) {
0216                 i = false;
0217               };
0218               for (const auto& st1 : stubsTrk1) {
0219                 for (const auto& st2 : stubsTrk2) {
0220                   if (st1.first == st2.first && st1.second == st2.second) {  // tracks share stub
0221                     // Converts layer/disk encoded in st1->first to an index in the layer array
0222                     int i = st1.first;  // layer/disk
0223                     bool barrel = (i > 0 && i < 10);
0224                     bool endcapA = (i > 10);
0225                     bool endcapB = (i < 0);
0226                     int lay = barrel * (i - 1) + endcapA * (i - 5) - endcapB * i;  // encode in range 0-15
0227                     if (!layerArr[lay]) {
0228                       nShareLay++;
0229                       layerArr[lay] = true;
0230                     }
0231                   }
0232                 }
0233               }
0234             } else if (settings_.mergeComparison() == "CompareBest") {
0235               std::vector<const Stub*> fullStubslistsTrk1 = inputstublists_[itrk];
0236               std::vector<const Stub*> fullStubslistsTrk2 = inputstublists_[jtrk];
0237
0238               // Arrays to store the index of the best stub in each layer
0239               int layStubidsTrk1[16];
0240               int layStubidsTrk2[16];
0241               for (int i = 0; i < 16; i++) {
0242                 layStubidsTrk1[i] = -1;
0243                 layStubidsTrk2[i] = -1;
0244               }
0245               // For each stub on the first track, find the stub with the best residual and store its index in the layStubidsTrk1 array
0246               for (unsigned int stcount = 0; stcount < stubsTrk1.size(); stcount++) {
0247                 int i = stubsTrk1[stcount].first;  // layer/disk
0248                 bool barrel = (i > 0 && i < 10);
0249                 bool endcapA = (i > 10);
0250                 bool endcapB = (i < 0);
0251                 int lay = barrel * (i - 1) + endcapA * (i - 5) - endcapB * i;  // encode in range 0-15
0252                 double nres = getPhiRes(inputtracklets_[itrk], fullStubslistsTrk1[stcount]);
0253                 double ores = 0;
0254                 if (layStubidsTrk1[lay] != -1)
0255                   ores = getPhiRes(inputtracklets_[itrk], fullStubslistsTrk1[layStubidsTrk1[lay]]);
0256                 if (layStubidsTrk1[lay] == -1 || nres < ores) {
0257                   layStubidsTrk1[lay] = stcount;
0258                 }
0259               }
0260               // For each stub on the second track, find the stub with the best residual and store its index in the layStubidsTrk1 array
0261               for (unsigned int stcount = 0; stcount < stubsTrk2.size(); stcount++) {
0262                 int i = stubsTrk2[stcount].first;  // layer/disk
0263                 bool barrel = (i > 0 && i < 10);
0264                 bool endcapA = (i > 10);
0265                 bool endcapB = (i < 0);
0266                 int lay = barrel * (i - 1) + endcapA * (i - 5) - endcapB * i;  // encode in range 0-15
0267                 double nres = getPhiRes(inputtracklets_[jtrk], fullStubslistsTrk2[stcount]);
0268                 double ores = 0;
0269                 if (layStubidsTrk2[lay] != -1)
0270                   ores = getPhiRes(inputtracklets_[jtrk], fullStubslistsTrk2[layStubidsTrk2[lay]]);
0271                 if (layStubidsTrk2[lay] == -1 || nres < ores) {
0272                   layStubidsTrk2[lay] = stcount;
0273                 }
0274               }
0275               // For all 16 layers (6 layers and 10 disks), count the number of layers who's best stub on both tracks are the same
0276               for (int i = 0; i < 16; i++) {
0277                 int t1i = layStubidsTrk1[i];
0278                 int t2i = layStubidsTrk2[i];
0279                 if (t1i != -1 && t2i != -1 && stubsTrk1[t1i].first == stubsTrk2[t2i].first &&
0280                     stubsTrk1[t1i].second == stubsTrk2[t2i].second)
0281                   nShareLay++;
0282               }
0283             }
0284
0285             // Fill duplicate map
0286             if (nShareLay >= settings_.minIndStubs()) {  // For number of shared stub merge condition
0287               dupMap.at(itrk).at(jtrk) = true;
0288               dupMap.at(jtrk).at(itrk) = true;
0289             }
0290           }
0291         }
0292
0293         // Check to see if the track is a duplicate
0294         for (unsigned int itrk = 0; itrk < numStublists; itrk++) {
0295           for (unsigned int jtrk = 0; jtrk < numStublists; jtrk++) {
0296             if (dupMap.at(itrk).at(jtrk)) {
0297               noMerge.at(itrk) = true;
0298             }
0299           }
0300         }
0301
0302         // If the track isn't a duplicate, and if it's in more than one bin, and it is not in the proper rinv or phi bin, then mark it so it won't be sent to output
0303         for (unsigned int itrk = 0; itrk < numStublists; itrk++) {
0304           if (!noMerge.at(itrk)) {
0305             if (((findOverlapRinvBins(inputtracklets_[itrk]).size() > 1) &&
0306                  (findRinvBin(inputtracklets_[itrk]) != bin)) ||
0307                 ((findOverlapPhiBins(inputtracklets_[itrk]).size() > 1) &&
0308                  findPhiBin(inputtracklets_[itrk]) != phiBin)) {
0309               trackInfo[itrk].second = true;
0310             }
0311           }
0312         }
0313         // Merge duplicate tracks
0314         for (unsigned int itrk = 0; itrk < numStublists - 1; itrk++) {
0315           for (unsigned int jtrk = itrk + 1; jtrk < numStublists; jtrk++) {
0316             // Merge a track with its first duplicate found.
0317             if (dupMap.at(itrk).at(jtrk)) {
0318               // Set preferred track based on seed rank
0319               int preftrk;
0320               int rejetrk;
0321               if (seedRank[itrk] < seedRank[jtrk]) {
0322                 preftrk = itrk;
0323                 rejetrk = jtrk;
0324               } else {
0325                 preftrk = jtrk;
0326                 rejetrk = itrk;
0327               }
0328
0329               // If the preffered track is in more than one bin, but not in the proper rinv or phi bin, then mark as true
0330               if (((findOverlapRinvBins(inputtracklets_[preftrk]).size() > 1) &&
0331                    (findRinvBin(inputtracklets_[preftrk]) != bin)) ||
0332                   ((findOverlapPhiBins(inputtracklets_[preftrk]).size() > 1) &&
0333                    (findPhiBin(inputtracklets_[preftrk]) != phiBin))) {
0334                 trackBinInfo[preftrk] = true;
0335                 trackBinInfo[rejetrk] = true;
0336               } else {
0337                 // Get a merged stub list
0338                 std::vector<const Stub*> newStubList;
0339                 std::vector<const Stub*> stubsTrk1 = inputstublists_[preftrk];
0340                 std::vector<const Stub*> stubsTrk2 = inputstublists_[rejetrk];
0341                 std::vector<unsigned int> stubsTrk1indices;
0342                 std::vector<unsigned int> stubsTrk2indices;
0343                 for (unsigned int stub1it = 0; stub1it < stubsTrk1.size(); stub1it++) {
0344                   stubsTrk1indices.push_back(stubsTrk1[stub1it]->l1tstub()->uniqueIndex());
0345                 }
0346                 for (unsigned int stub2it = 0; stub2it < stubsTrk2.size(); stub2it++) {
0347                   stubsTrk2indices.push_back(stubsTrk2[stub2it]->l1tstub()->uniqueIndex());
0348                 }
0349                 newStubList = stubsTrk1;
0350                 for (unsigned int stub2it = 0; stub2it < stubsTrk2.size(); stub2it++) {
0351                   if (find(stubsTrk1indices.begin(), stubsTrk1indices.end(), stubsTrk2indices[stub2it]) ==
0352                       stubsTrk1indices.end()) {
0353                     newStubList.push_back(stubsTrk2[stub2it]);
0354                   }
0355                 }
0356                 //   Overwrite stublist of preferred track with merged list
0357                 inputstublists_[preftrk] = newStubList;
0358
0359                 std::vector<std::pair<int, int>> newStubidsList;
0360                 std::vector<std::pair<int, int>> stubidsTrk1 = mergedstubidslists_[preftrk];
0361                 std::vector<std::pair<int, int>> stubidsTrk2 = mergedstubidslists_[rejetrk];
0362                 newStubidsList = stubidsTrk1;
0363
0364                 for (unsigned int stub2it = 0; stub2it < stubsTrk2.size(); stub2it++) {
0365                   if (find(stubsTrk1indices.begin(), stubsTrk1indices.end(), stubsTrk2indices[stub2it]) ==
0366                       stubsTrk1indices.end()) {
0367                     newStubidsList.push_back(stubidsTrk2[stub2it]);
0368                   }
0369                 }
0370                 // Overwrite stubidslist of preferred track with merged list
0371                 mergedstubidslists_[preftrk] = newStubidsList;
0372
0373                 // Mark that rejected track has been merged into another track
0374                 trackInfo[rejetrk].second = true;
0375               }
0376             }
0377           }
0378         }
0379
0380         for (unsigned int ktrk = 0; ktrk < numStublists; ktrk++) {
0381           if ((trackInfo[ktrk].second != true) && (trackBinInfo[ktrk] != true)) {
0382             prefTracks.push_back(ktrk);
0383             prefTrackFit.push_back(trackInfo[ktrk].first);
0384             inputtrackletsall.push_back(inputtracklets_[ktrk]);
0385             inputstublistsall.push_back(inputstublists_[ktrk]);
0386             inputstubidslistsall.push_back(inputstubidslists_[ktrk]);
0387             mergedstubidslistsall.push_back(mergedstubidslists_[ktrk]);
0388           }
0389         }
0390
0391         // Need to clear all the vectors which will be used in the next bin
0392         seedRank.clear();
0393         trackInfo.clear();
0394         trackBinInfo.clear();
0395         inputtracklets_.clear();
0396         inputstublists_.clear();
0397         inputstubidslists_.clear();
0398         mergedstubidslists_.clear();
0399       }
0400     }
0401
0402     // Make the final track objects, fit with KF, and send to output
0403     for (unsigned int itrk = 0; itrk < prefTracks.size(); itrk++) {
0404       Tracklet* tracklet = inputtrackletsall[itrk];
0405       std::vector<const Stub*> trackstublist = inputstublistsall[itrk];
0406
0407       // Run KF track fit
0408       HybridFit hybridFitter(iSector, settings_, globals_);
0409       hybridFitter.Fit(tracklet, trackstublist);
0410
0411       // If the track was accepted (and thus fit), add to output
0412       if (tracklet->fit()) {
0413         // Add fitted Track to output (later converted to TTTrack)
0414         Track* outtrack = tracklet->getTrack();
0415         outtrack->setSector(iSector);
0416         // Also store fitted track as more detailed Tracklet object.
0417         outputtracklets_[prefTrackFit[itrk]]->addTrack(tracklet);
0418
0419         // Add all tracks to standalone root file output
0420         outtrack->setStubIDpremerge(inputstubidslistsall[itrk]);
0421         outtrack->setStubIDprefit(mergedstubidslistsall[itrk]);
0422         outputtracks.push_back(*outtrack);
0423       }
0424     }
0425   }
0426
0427 #endif
0428
0429   //////////////////
0430   // Grid removal //
0431   //////////////////
0432   if (settings_.removalType() == "grid") {
0433     // Sort tracks by ichisq/DoF so that removal will keep the lower ichisq/DoF track
0434     std::sort(inputtracks_.begin(), inputtracks_.end(), [](const Track* lhs, const Track* rhs) {
0435       return lhs->ichisq() / lhs->stubID().size() < rhs->ichisq() / rhs->stubID().size();
0436     });
0437     bool grid[35][40] = {{false}};
0438
0439     for (unsigned int itrk = 0; itrk < numTrk; itrk++) {
0440       if (inputtracks_[itrk]->duplicate())
0441         edm::LogPrint("Tracklet") << "WARNING: Track already tagged as duplicate!!";
0442
0443       double phiBin = (inputtracks_[itrk]->phi0(settings_) - 2 * M_PI / 27 * iSector) / (2 * M_PI / 9 / 50) + 9;
0444       phiBin = std::max(phiBin, 0.);
0445       phiBin = std::min(phiBin, 34.);
0446
0447       double ptBin = 1 / inputtracks_[itrk]->pt(settings_) * 40 + 20;
0448       ptBin = std::max(ptBin, 0.);
0449       ptBin = std::min(ptBin, 39.);
0450
0451       if (grid[(int)phiBin][(int)ptBin])
0452         inputtracks_[itrk]->setDuplicate(true);
0453       grid[(int)phiBin][(int)ptBin] = true;
0454
0455       double phiTest = inputtracks_[itrk]->phi0(settings_) - 2 * M_PI / 27 * iSector;
0456       if (phiTest < -2 * M_PI / 27)
0457         edm::LogVerbatim("Tracklet") << "track phi too small!";
0458       if (phiTest > 2 * 2 * M_PI / 27)
0459         edm::LogVerbatim("Tracklet") << "track phi too big!";
0460     }
0461   }  // end grid removal
0462
0463   //////////////////////////
0464   // ichi + nstub removal //
0465   //////////////////////////
0466   if (settings_.removalType() == "ichi" || settings_.removalType() == "nstub") {
0467     for (unsigned int itrk = 0; itrk < numTrk - 1; itrk++) {  // numTrk-1 since last track has no other to compare to
0468
0469       // If primary track is a duplicate, it cannot veto any...move on
0470       if (inputtracks_[itrk]->duplicate() == 1)
0471         continue;
0472
0473       unsigned int nStubP = 0;
0474       vector<unsigned int> nStubS(numTrk);
0475       vector<unsigned int> nShare(numTrk);
0476       // Get and count primary stubs
0477       std::map<int, int> stubsTrk1 = inputtracks_[itrk]->stubID();
0478       nStubP = stubsTrk1.size();
0479
0480       for (unsigned int jtrk = itrk + 1; jtrk < numTrk; jtrk++) {
0481         // Skip duplicate tracks
0482         if (inputtracks_[jtrk]->duplicate() == 1)
0483           continue;
0484
0485         // Get and count secondary stubs
0486         std::map<int, int> stubsTrk2 = inputtracks_[jtrk]->stubID();
0487         nStubS[jtrk] = stubsTrk2.size();
0488
0489         // Count shared stubs
0490         for (auto& st : stubsTrk1) {
0491           if (stubsTrk2.find(st.first) != stubsTrk2.end()) {
0492             if (st.second == stubsTrk2[st.first])
0493               nShare[jtrk]++;
0494           }
0495         }
0496       }
0497
0498       // Tag duplicates
0499       for (unsigned int jtrk = itrk + 1; jtrk < numTrk; jtrk++) {
0500         // Skip duplicate tracks
0501         if (inputtracks_[jtrk]->duplicate() == 1)
0502           continue;
0503
0504         // Chi2 duplicate removal
0505         if (settings_.removalType() == "ichi") {
0506           if ((nStubP - nShare[jtrk] < settings_.minIndStubs()) ||
0507               (nStubS[jtrk] - nShare[jtrk] < settings_.minIndStubs())) {
0508             if ((int)inputtracks_[itrk]->ichisq() / (2 * inputtracks_[itrk]->stubID().size() - 4) >
0509                 (int)inputtracks_[jtrk]->ichisq() / (2 * inputtracks_[itrk]->stubID().size() - 4)) {
0510               inputtracks_[itrk]->setDuplicate(true);
0511             } else if ((int)inputtracks_[itrk]->ichisq() / (2 * inputtracks_[itrk]->stubID().size() - 4) <=
0512                        (int)inputtracks_[jtrk]->ichisq() / (2 * inputtracks_[itrk]->stubID().size() - 4)) {
0513               inputtracks_[jtrk]->setDuplicate(true);
0514             } else {
0515               edm::LogVerbatim("Tracklet") << "Error: Didn't tag either track in duplicate pair.";
0516             }
0517           }
0518         }  // end ichi removal
0519
0520         // nStub duplicate removal
0521         if (settings_.removalType() == "nstub") {
0522           if ((nStubP - nShare[jtrk] < settings_.minIndStubs()) && (nStubP < nStubS[jtrk])) {
0523             inputtracks_[itrk]->setDuplicate(true);
0524           } else if ((nStubS[jtrk] - nShare[jtrk] < settings_.minIndStubs()) && (nStubS[jtrk] <= nStubP)) {
0525             inputtracks_[jtrk]->setDuplicate(true);
0526           } else {
0527             edm::LogVerbatim("Tracklet") << "Error: Didn't tag either track in duplicate pair.";
0528           }
0529         }  // end nstub removal
0530
0531       }  // end tag duplicates
0532
0533     }  // end loop over primary track
0534
0535   }  // end ichi + nstub removal
0536
0537   //Add tracks to output
0538   if (settings_.removalType() != "merge") {
0539     for (unsigned int i = 0; i < inputtrackfits_.size(); i++) {
0540       for (unsigned int j = 0; j < inputtrackfits_[i]->nTracks(); j++) {
0541         if (inputtrackfits_[i]->getTrack(j)->getTrack()->duplicate() == 0) {
0542           if (settings_.writeMonitorData("Seeds")) {
0543             ofstream fout("seeds.txt", ofstream::app);
0544             fout << __FILE__ << ":" << __LINE__ << " " << name_ << "_" << iSector << " "
0545                  << inputtrackfits_[i]->getTrack(j)->getISeed() << endl;
0546             fout.close();
0547           }
0548           outputtracklets_[i]->addTrack(inputtrackfits_[i]->getTrack(j));
0549         }
0550         //For root file:
0551         outputtracks.push_back(*inputtrackfits_[i]->getTrack(j)->getTrack());
0552       }
0553     }
0554   }
0555 }
0556
0557 double PurgeDuplicate::getPhiRes(Tracklet* curTracklet, const Stub* curStub) const {
0558   double phiproj;
0559   double stubphi;
0560   double phires;
0561   // Get phi position of stub
0562   stubphi = curStub->l1tstub()->phi();
0563   // Get layer that the stub is in (Layer 1->6, Disk 1->5)
0564   int Layer = curStub->layerdisk() + 1;
0565   if (Layer > N_LAYER) {
0566     Layer = 0;
0567   }
0568   int Disk = curStub->layerdisk() - (N_LAYER - 1);
0569   if (Disk < 0) {
0570     Disk = 0;
0571   }
0572   // Get phi projection of tracklet
0573   int seedindex = curTracklet->seedIndex();
0574   // If this stub is a seed stub, set projection=phi, so that res=0
0575   if (isSeedingStub(seedindex, Layer, Disk)) {
0576     phiproj = stubphi;
0577     // Otherwise, get projection of tracklet
0578   } else {
0579     phiproj = curTracklet->proj(curStub->layerdisk()).phiproj();
0580   }
0581   // Calculate residual
0582   phires = std::abs(stubphi - phiproj);
0583   return phires;
0584 }
0585
0586 bool PurgeDuplicate::isSeedingStub(int seedindex, int Layer, int Disk) const {
0587   if ((seedindex == 0 && (Layer == 1 || Layer == 2)) || (seedindex == 1 && (Layer == 2 || Layer == 3)) ||
0588       (seedindex == 2 && (Layer == 3 || Layer == 4)) || (seedindex == 3 && (Layer == 5 || Layer == 6)) ||
0589       (seedindex == 4 && (abs(Disk) == 1 || abs(Disk) == 2)) ||
0590       (seedindex == 5 && (abs(Disk) == 3 || abs(Disk) == 4)) || (seedindex == 6 && (Layer == 1 || abs(Disk) == 1)) ||
0591       (seedindex == 7 && (Layer == 2 || abs(Disk) == 1)) ||
0592       (seedindex == 8 && (Layer == 2 || Layer == 3 || Layer == 4)) ||
0593       (seedindex == 9 && (Layer == 4 || Layer == 5 || Layer == 6)) ||
0594       (seedindex == 10 && (Layer == 2 || Layer == 3 || abs(Disk) == 1)) ||
0595       (seedindex == 11 && (Layer == 2 || abs(Disk) == 1 || abs(Disk) == 2)))
0596     return true;
0597
0598   return false;
0599 }
0600
0601 std::pair<int, int> PurgeDuplicate::findLayerDisk(const Stub* st) const {
0602   std::pair<int, int> layer_disk;
0603   layer_disk.first = st->layerdisk() + 1;
0604   if (layer_disk.first > N_LAYER) {
0605     layer_disk.first = 0;
0606   }
0607   layer_disk.second = st->layerdisk() - (N_LAYER - 1);
0608   if (layer_disk.second < 0) {
0609     layer_disk.second = 0;
0610   }
0611   return layer_disk;
0612 }
0613
0614 std::string PurgeDuplicate::l1tinfo(const L1TStub* l1stub, std::string str = "") const {
0615   // Uses ROOT::TString
0616   std::string thestr = Form("\t %s stub info:  r/z/phi:\t%f\t%f\t%f\t%d\t%f\t%d",
0617                             str.c_str(),
0618                             l1stub->r(),
0619                             l1stub->z(),
0620                             l1stub->phi(),
0621                             l1stub->iphi(),
0622                             l1stub->bend(),
0623                             l1stub->layerdisk());
0624   return thestr;
0625 }
0626
0627 std::vector<double> PurgeDuplicate::getInventedCoords(unsigned int iSector,
0628                                                       const Stub* st,
0629                                                       const Tracklet* tracklet) const {
0630   int stubLayer = (findLayerDisk(st)).first;
0631   int stubDisk = (findLayerDisk(st)).second;
0632
0633   double stub_phi = -99;
0634   double stub_z = -99;
0635   double stub_r = -99;
0636
0637   double tracklet_rinv = tracklet->rinv();
0638
0639   if (st->isBarrel()) {
0640     stub_r = settings_.rmean(stubLayer - 1);
0641     stub_phi = tracklet->phi0() - std::asin(stub_r * tracklet_rinv / 2);
0642     stub_phi = stub_phi + iSector * settings_.dphisector() - 0.5 * settings_.dphisectorHG();
0643     stub_phi = reco::reduceRange(stub_phi);
0644     stub_z = tracklet->z0() + 2 * tracklet->t() * 1 / tracklet_rinv * std::asin(stub_r * tracklet_rinv / 2);
0645   } else {
0646     stub_z = settings_.zmean(stubDisk - 1) * tracklet->disk() / abs(tracklet->disk());
0647     stub_phi = tracklet->phi0() - (stub_z - tracklet->z0()) * tracklet_rinv / 2 / tracklet->t();
0648     stub_phi = stub_phi + iSector * settings_.dphisector() - 0.5 * settings_.dphisectorHG();
0649     stub_phi = reco::reduceRange(stub_phi);
0650     stub_r = 2 / tracklet_rinv * std::sin((stub_z - tracklet->z0()) * tracklet_rinv / 2 / tracklet->t());
0651   }
0652
0653   std::vector<double> invented_coords{stub_r, stub_z, stub_phi};
0654   return invented_coords;
0655 }
0656
0657 std::vector<double> PurgeDuplicate::getInventedCoordsExtended(unsigned int iSector,
0658                                                               const Stub* st,
0659                                                               const Tracklet* tracklet) const {
0660   const int stubLayer = (findLayerDisk(st)).first;
0661   const int stubDisk = (findLayerDisk(st)).second;
0662
0663   double stub_phi = -99;
0664   double stub_z = -99;
0665   double stub_r = -99;
0666
0667   const double rho = 1 / tracklet->rinv();
0668   const double rho_minus_d0 = rho + tracklet->d0();  // should be -, but otherwise does not work
0669
0670   const int seed = tracklet->seedIndex();
0671
0672   // TMP: for displaced tracking, exclude one of the 3 seeding stubs
0673   // to be discussed
0674   if ((seed == L2L3L4 && stubLayer == 4) || (seed == L4L5L6 && stubLayer == 5) ||
0675       (seed == L2L3D1 && abs(stubDisk) == 1) || (seed == D1D2L2 && abs(stubDisk) == 1)) {
0676     stub_phi = st->l1tstub()->phi();
0677     stub_z = st->l1tstub()->z();
0678     stub_r = st->l1tstub()->r();
0679   } else {
0680     // exact helix
0681     if (st->isBarrel()) {
0682       stub_r = settings_.rmean(stubLayer - 1);
0683
0684       // The expanded version of this expression is more stable for extremely
0685       // high-pT (high-rho) tracks. But we also explicitly restrict sin_val to
0686       // the domain of asin.
0687       double sin_val =
0688           0.5 * (stub_r / rho_minus_d0) + 0.5 * (rho_minus_d0 / stub_r) - 0.5 * ((rho * rho) / (rho_minus_d0 * stub_r));
0689       sin_val = std::max(std::min(sin_val, 1.0), -1.0);
0690       stub_phi = tracklet->phi0() - std::asin(sin_val);
0691       stub_phi = stub_phi + iSector * settings_.dphisector() - 0.5 * settings_.dphisectorHG();
0692       stub_phi = reco::reduceRange(stub_phi);
0693
0694       // The expanded version of this expression is more stable for extremely
0695       // high-pT (high-rho) tracks. But we also explicitly restrict cos_val to
0696       // the domain of acos.
0697       double cos_val =
0698           0.5 * (rho / rho_minus_d0) + 0.5 * (rho_minus_d0 / rho) - 0.5 * ((stub_r * stub_r) / (rho * rho_minus_d0));
0699       cos_val = std::max(std::min(cos_val, 1.0), -1.0);
0700       double beta = std::acos(cos_val);
0701       stub_z = tracklet->z0() + tracklet->t() * std::abs(rho * beta);
0702     } else {
0703       stub_z = settings_.zmean(stubDisk - 1) * tracklet->disk() / abs(tracklet->disk());
0704
0705       double beta = (stub_z - tracklet->z0()) / (tracklet->t() * std::abs(rho));  // maybe rho should be abs value
0706       double r_square = -2 * rho * rho_minus_d0 * std::cos(beta) + rho * rho + rho_minus_d0 * rho_minus_d0;
0707       stub_r = sqrt(r_square);
0708
0709       // The expanded version of this expression is more stable for extremely
0710       // high-pT (high-rho) tracks. But we also explicitly restrict sin_val to
0711       // the domain of asin.
0712       double sin_val =
0713           0.5 * (stub_r / rho_minus_d0) + 0.5 * (rho_minus_d0 / stub_r) - 0.5 * ((rho * rho) / (rho_minus_d0 * stub_r));
0714       sin_val = std::max(std::min(sin_val, 1.0), -1.0);
0715       stub_phi = tracklet->phi0() - std::asin(sin_val);
0716       stub_phi = stub_phi + iSector * settings_.dphisector() - 0.5 * settings_.dphisectorHG();
0717       stub_phi = reco::reduceRange(stub_phi);
0718     }
0719   }
0720
0721   std::vector<double> invented_coords{stub_r, stub_z, stub_phi};
0722   return invented_coords;
0723 }
0724
0725 std::vector<const Stub*> PurgeDuplicate::getInventedSeedingStub(
0726     unsigned int iSector, const Tracklet* tracklet, const std::vector<const Stub*>& originalStubsList) const {
0727   std::vector<const Stub*> newStubList;
0728
0729   for (unsigned int stubit = 0; stubit < originalStubsList.size(); stubit++) {
0730     const Stub* thisStub = originalStubsList[stubit];
0731
0732     if (isSeedingStub(tracklet->seedIndex(), (findLayerDisk(thisStub)).first, (findLayerDisk(thisStub)).second)) {
0733       // get a vector containing r, z, phi
0734       std::vector<double> inv_r_z_phi;
0735       if (!settings_.extended())
0736         inv_r_z_phi = getInventedCoords(iSector, thisStub, tracklet);
0737       else {
0738         inv_r_z_phi = getInventedCoordsExtended(iSector, thisStub, tracklet);
0739       }
0740       double stub_x_invent = inv_r_z_phi[0] * std::cos(inv_r_z_phi[2]);
0741       double stub_y_invent = inv_r_z_phi[0] * std::sin(inv_r_z_phi[2]);
0742       double stub_z_invent = inv_r_z_phi[1];
0743
0744       Stub* invent_stub_ptr = new Stub(*thisStub);
0745       const L1TStub* l1stub = thisStub->l1tstub();
0746       L1TStub invent_l1stub = *l1stub;
0747       invent_l1stub.setCoords(stub_x_invent, stub_y_invent, stub_z_invent);
0748
0749       invent_stub_ptr->setl1tstub(new L1TStub(invent_l1stub));
0750       invent_stub_ptr->l1tstub()->setAllStubIndex(l1stub->allStubIndex());
0751       invent_stub_ptr->l1tstub()->setUniqueIndex(l1stub->uniqueIndex());
0752
0753       newStubList.push_back(invent_stub_ptr);
0754
0755     } else {
0756       newStubList.push_back(thisStub);
0757     }
0758   }
0759   return newStubList;
0760 }
0761
0762 // Tells us the variable bin to which a track would belong
0763 unsigned int PurgeDuplicate::findRinvBin(const Tracklet* trk) const {
0764   std::vector<double> rinvBins = settings_.rinvBins();
0765
0766   //Get rinverse of track
0767   double rinv = trk->rinv();
0768
0769   //Check between what 2 values in rinvbins rinv is between
0770   auto bins = std::upper_bound(rinvBins.begin(), rinvBins.end(), rinv);
0771
0772   //return integer for bin index
0773   unsigned int rIndx = std::distance(rinvBins.begin(), bins);
0774   if (rIndx == std::distance(rinvBins.end(), bins))
0775     return rinvBins.size() - 2;
0776   else if (bins == rinvBins.begin())
0777     return std::distance(rinvBins.begin(), bins);
0778   else
0779     return rIndx - 1;
0780 }
0781
0782 // Tells us the phi bin to which a track would belong
0783 unsigned int PurgeDuplicate::findPhiBin(const Tracklet* trk) const {
0784   std::vector<double> phiBins = settings_.phiBins();
0785
0786   //Get phi of track at rcrit
0787   double phi0 = trk->phi0();
0788   double rcrit = settings_.rcrit();
0789   double rinv = trk->rinv();
0790   double phi = phi0 - asin(0.5 * rinv * rcrit);
0791
0792   //Check between what 2 values in phibins phi is between
0793   auto bins = std::upper_bound(phiBins.begin(), phiBins.end(), phi);
0794
0795   //return integer for bin index
0796   unsigned int phiIndx = std::distance(phiBins.begin(), bins);
0797   if (phiIndx == std::distance(phiBins.end(), bins))
0798     return phiBins.size() - 2;
0799   else if (bins == phiBins.begin())
0800     return std::distance(phiBins.begin(), bins);
0801   else
0802     return phiIndx - 1;
0803 }
0804
0805 // Tells us the overlap rinv bin(s) to which a track belongs
0806 std::vector<unsigned int> PurgeDuplicate::findOverlapRinvBins(const Tracklet* trk) const {
0807   double rinv = trk->rinv();
0808
0809   const double rinvOverlapSize = settings_.rinvOverlapSize();
0810   const std::vector<double>& rinvBins = settings_.rinvBins();
0811
0812   std::vector<unsigned int> chosenBins;
0813   for (long unsigned int i = 0; i < rinvBins.size() - 1; i++) {
0814     if ((rinv < rinvBins[i + 1] + rinvOverlapSize) && (rinv > rinvBins[i] - rinvOverlapSize)) {
0815       chosenBins.push_back(i);
0816     }
0817   }
0818   return chosenBins;
0819 }
0820
0821 // Tells us the overlap phi bin(s) to which a track belongs
0822 std::vector<unsigned int> PurgeDuplicate::findOverlapPhiBins(const Tracklet* trk) const {
0823   double phi0 = trk->phi0();
0824   double rcrit = settings_.rcrit();
0825   double rinv = trk->rinv();
0826   double phi = phi0 - asin(0.5 * rinv * rcrit);
0827
0828   const double phiOverlapSize = settings_.phiOverlapSize();
0829   const std::vector<double>& phiBins = settings_.phiBins();
0830
0831   std::vector<unsigned int> chosenBins;
0832   for (long unsigned int i = 0; i < phiBins.size() - 1; i++) {
0833     if ((phi < phiBins[i + 1] + phiOverlapSize) && (phi > phiBins[i] - phiOverlapSize)) {
0834       chosenBins.push_back(i);
0835     }
0836   }
0837   return chosenBins;
0838 }
0839
0840 // Tells us if a track is in the current bin
0841 bool PurgeDuplicate::isTrackInBin(const std::vector<unsigned int>& vec, unsigned int num) const {
0842   auto result = std::find(vec.begin(), vec.end(), num);
0843   bool found = (result != vec.end());
0844   return found;
0845 }