Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-12-12 02:57:19

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