PurgeDuplicate.cc

Go to the documentation of this file.

#include "L1Trigger/TrackFindingTracklet/interface/PurgeDuplicate.h"
#include "L1Trigger/TrackFindingTracklet/interface/Tracklet.h"
#include "L1Trigger/TrackFindingTracklet/interface/Globals.h"
#include "L1Trigger/TrackFindingTracklet/interface/CleanTrackMemory.h"
#include "L1Trigger/TrackFindingTracklet/interface/Settings.h"
#include "L1Trigger/TrackFindingTracklet/interface/Stub.h"
#include "L1Trigger/TrackFindingTracklet/interface/Track.h"

#ifdef USEHYBRID
#include "DataFormats/L1TrackTrigger/interface/TTStub.h"
#include "L1Trigger/TrackFindingTMTT/interface/L1track3D.h"
#include "L1Trigger/TrackFindingTMTT/interface/KFParamsComb.h"
#include "L1Trigger/TrackFindingTracklet/interface/HybridFit.h"
#endif

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"

#include <unordered_set>
#include <algorithm>

using namespace std;
using namespace trklet;

PurgeDuplicate::PurgeDuplicate(std::string name, Settings const& settings, Globals* global)
    : ProcessBase(name, settings, global) {}

void PurgeDuplicate::addOutput(MemoryBase* memory, std::string output) {
  if (settings_.writetrace()) {
    edm::LogVerbatim("Tracklet") << "In " << name_ << " adding output to " << memory->getName() << " to output "
                                 << output;
  }
  unordered_set<string> outputs = {"trackout",
                                   "trackout1",
                                   "trackout2",
                                   "trackout3",
                                   "trackout4",
                                   "trackout5",
                                   "trackout6",
                                   "trackout7",
                                   "trackout8",
                                   "trackout9",
                                   "trackout10",
                                   "trackout11"};
  if (outputs.find(output) != outputs.end()) {
    auto* tmp = dynamic_cast<CleanTrackMemory*>(memory);
    assert(tmp != nullptr);
    outputtracklets_.push_back(tmp);
    return;
  }
  throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " could not find output: " << output;
}

void PurgeDuplicate::addInput(MemoryBase* memory, std::string input) {
  if (settings_.writetrace()) {
    edm::LogVerbatim("Tracklet") << "In " << name_ << " adding input from " << memory->getName() << " to input "
                                 << input;
  }
  unordered_set<string> inputs = {"trackin",
                                  "trackin1",
                                  "trackin2",
                                  "trackin3",
                                  "trackin4",
                                  "trackin5",
                                  "trackin6",
                                  "trackin7",
                                  "trackin8",
                                  "trackin9",
                                  "trackin10",
                                  "trackin11",
                                  "trackin12"};
  if (inputs.find(input) != inputs.end()) {
    auto* tmp = dynamic_cast<TrackFitMemory*>(memory);
    assert(tmp != nullptr);
    inputtrackfits_.push_back(tmp);
    return;
  }
  throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " could not find input: " << input;
}

void PurgeDuplicate::execute(std::vector<Track>& outputtracks, unsigned int iSector) {
  inputtracklets_.clear();
  inputtracks_.clear();

  inputstubidslists_.clear();
  inputstublists_.clear();
  mergedstubidslists_.clear();

  if (settings_.removalType() != "merge") {
    for (auto& inputtrackfit : inputtrackfits_) {
      if (inputtrackfit->nTracks() == 0)
        continue;
      for (unsigned int j = 0; j < inputtrackfit->nTracks(); j++) {
        Track* aTrack = inputtrackfit->getTrack(j)->getTrack();
        aTrack->setSector(iSector);
        inputtracks_.push_back(aTrack);
      }
    }
    if (inputtracks_.empty())
      return;
  }

  unsigned int numTrk = inputtracks_.size();

  // Hybrid Removal //
#ifdef USEHYBRID

  if (settings_.removalType() == "merge") {
    // Track seed & duplicate flag
    std::vector<std::pair<int, bool>> trackInfo;
    // Flag for tracks in multiple bins that get merged but are not in the correct bin
    std::vector<bool> trackBinInfo;
    // Vector to store the relative rank of the track candidate for merging, based on seed type
    std::vector<int> seedRank;

    // Stubs on every track
    std::vector<std::vector<const Stub*>> inputstublistsall;
    // (layer, unique stub index within layer) of each stub on every track
    std::vector<std::vector<std::pair<int, int>>> mergedstubidslistsall;
    std::vector<std::vector<std::pair<int, int>>> inputstubidslistsall;
    std::vector<Tracklet*> inputtrackletsall;

    std::vector<unsigned int> prefTracks;  // Stores all the tracks that are sent to the KF from each bin
    std::vector<int> prefTrackFit;  // Stores the track seed that corresponds to the associated track in prefTracks

    for (unsigned int bin = 0; bin < settings_.rinvBins().size() - 1; bin++) {
      for (unsigned int phiBin = 0; phiBin < settings_.phiBins().size() - 1; phiBin++) {
        // Get vectors from TrackFit and save them
        // inputtracklets: Tracklet objects from the FitTrack (not actually fit yet)
        // inputstublists: L1Stubs for that track
        // inputstubidslists: Stub stubIDs for that 3rack
        // mergedstubidslists: the same as inputstubidslists, but will be used during duplicate removal
        for (unsigned int i = 0; i < inputtrackfits_.size(); i++) {
          if (inputtrackfits_[i]->nStublists() == 0)
            continue;
          if (inputtrackfits_[i]->nStublists() != inputtrackfits_[i]->nTracks())
            throw cms::Exception("LogicError")
                << __FILE__ << " " << __LINE__ << " Number of stublists and tracks don't match up! ";
          for (unsigned int j = 0; j < inputtrackfits_[i]->nStublists(); j++) {
            if (isTrackInBin(findOverlapRinvBins(inputtrackfits_[i]->getTrack(j)), bin)) {
              if (!isTrackInBin(findOverlapPhiBins(inputtrackfits_[i]->getTrack(j)), phiBin))
                continue;
              if (inputtracklets_.size() >= settings_.maxStep("DR"))
                continue;
              Tracklet* aTrack = inputtrackfits_[i]->getTrack(j);
              inputtracklets_.push_back(inputtrackfits_[i]->getTrack(j));
              std::vector<const Stub*> stublist = inputtrackfits_[i]->getStublist(j);
              inputstublists_.push_back(stublist);
              std::vector<std::pair<int, int>> stubidslist = inputtrackfits_[i]->getStubidslist(j);
              inputstubidslists_.push_back(stubidslist);
              mergedstubidslists_.push_back(stubidslist);

              // Encoding: L1L2=0, L2L3=1, L3L4=2, L5L6=3, D1D2=4, D3D4=5, L1D1=6, L2D1=7
              // Best Guess:          L1L2 > L1D1 > L2L3 > L2D1 > D1D2 > L3L4 > L5L6 > D3D4
              // Best Rank:           L1L2 > L3L4 > D3D4 > D1D2 > L2L3 > L2D1 > L5L6 > L1D1
              // Rank-Informed Guess: L1L2 > L3L4 > L1D1 > L2L3 > L2D1 > D1D2 > L5L6 > D3D4
              unsigned int curSeed = aTrack->seedIndex();
              std::vector<int> ranks{1, 5, 2, 7, 4, 3, 8, 6};
              if (settings_.extended())
                seedRank.push_back(9);
              else
                seedRank.push_back(ranks[curSeed]);

              if (stublist.size() != stubidslist.size())
                throw cms::Exception("LogicError")
                    << __FILE__ << " " << __LINE__ << " Number of stubs and stubids don't match up! ";

              trackInfo.emplace_back(i, false);
              trackBinInfo.emplace_back(false);
            } else
              continue;
          }
        }

        if (inputtracklets_.empty())
          continue;
        const unsigned int numStublists = inputstublists_.size();

        if (settings_.inventStubs()) {
          for (unsigned int itrk = 0; itrk < numStublists; itrk++) {
            inputstublists_[itrk] = getInventedSeedingStub(iSector, inputtracklets_[itrk], inputstublists_[itrk]);
          }
        }

        // Initialize all-false 2D vector of tracks being duplicates to other tracks
        vector<vector<bool>> dupMap(numStublists, vector<bool>(numStublists, false));

        // Used to check if a track is in two bins, is not a duplicate in either bin, so is sent out twice
        vector<bool> noMerge(numStublists, false);

        // Find duplicates; Fill dupMap by looping over all pairs of "tracks"
        // numStublists-1 since last track has no other to compare to
        for (unsigned int itrk = 0; itrk < numStublists - 1; itrk++) {
          for (unsigned int jtrk = itrk + 1; jtrk < numStublists; jtrk++) {
            if (itrk >= settings_.numTracksComparedPerBin())
              continue;
            // Get primary track stubids = (layer, unique stub index within layer)
            const std::vector<std::pair<int, int>>& stubsTrk1 = inputstubidslists_[itrk];

            // Get and count secondary track stubids
            const std::vector<std::pair<int, int>>& stubsTrk2 = inputstubidslists_[jtrk];

            // Count number of layers that share stubs, and the number of UR that each track hits
            unsigned int nShareLay = 0;
            if (settings_.mergeComparison() == "CompareAll") {
              bool layerArr[16];
              for (auto& i : layerArr) {
                i = false;
              };
              for (const auto& st1 : stubsTrk1) {
                for (const auto& st2 : stubsTrk2) {
                  if (st1.first == st2.first && st1.second == st2.second) {  // tracks share stub
                    // Converts layer/disk encoded in st1->first to an index in the layer array
                    int i = st1.first;  // layer/disk
                    bool barrel = (i > 0 && i < 10);
                    bool endcapA = (i > 10);
                    bool endcapB = (i < 0);
                    int lay = barrel * (i - 1) + endcapA * (i - 5) - endcapB * i;  // encode in range 0-15
                    if (!layerArr[lay]) {
                      nShareLay++;
                      layerArr[lay] = true;
                    }
                  }
                }
              }
            } else if (settings_.mergeComparison() == "CompareBest") {
              std::vector<const Stub*> fullStubslistsTrk1 = inputstublists_[itrk];
              std::vector<const Stub*> fullStubslistsTrk2 = inputstublists_[jtrk];

              // Arrays to store the index of the best stub in each layer
              int layStubidsTrk1[16];
              int layStubidsTrk2[16];
              for (int i = 0; i < 16; i++) {
                layStubidsTrk1[i] = -1;
                layStubidsTrk2[i] = -1;
              }
              // For each stub on the first track, find the stub with the best residual and store its index in the layStubidsTrk1 array
              for (unsigned int stcount = 0; stcount < stubsTrk1.size(); stcount++) {
                int i = stubsTrk1[stcount].first;  // layer/disk
                bool barrel = (i > 0 && i < 10);
                bool endcapA = (i > 10);
                bool endcapB = (i < 0);
                int lay = barrel * (i - 1) + endcapA * (i - 5) - endcapB * i;  // encode in range 0-15
                double nres = getPhiRes(inputtracklets_[itrk], fullStubslistsTrk1[stcount]);
                double ores = 0;
                if (layStubidsTrk1[lay] != -1)
                  ores = getPhiRes(inputtracklets_[itrk], fullStubslistsTrk1[layStubidsTrk1[lay]]);
                if (layStubidsTrk1[lay] == -1 || nres < ores) {
                  layStubidsTrk1[lay] = stcount;
                }
              }
              // For each stub on the second track, find the stub with the best residual and store its index in the layStubidsTrk1 array
              for (unsigned int stcount = 0; stcount < stubsTrk2.size(); stcount++) {
                int i = stubsTrk2[stcount].first;  // layer/disk
                bool barrel = (i > 0 && i < 10);
                bool endcapA = (i > 10);
                bool endcapB = (i < 0);
                int lay = barrel * (i - 1) + endcapA * (i - 5) - endcapB * i;  // encode in range 0-15
                double nres = getPhiRes(inputtracklets_[jtrk], fullStubslistsTrk2[stcount]);
                double ores = 0;
                if (layStubidsTrk2[lay] != -1)
                  ores = getPhiRes(inputtracklets_[jtrk], fullStubslistsTrk2[layStubidsTrk2[lay]]);
                if (layStubidsTrk2[lay] == -1 || nres < ores) {
                  layStubidsTrk2[lay] = stcount;
                }
              }
              // For all 16 layers (6 layers and 10 disks), count the number of layers who's best stub on both tracks are the same
              for (int i = 0; i < 16; i++) {
                int t1i = layStubidsTrk1[i];
                int t2i = layStubidsTrk2[i];
                if (t1i != -1 && t2i != -1 && stubsTrk1[t1i].first == stubsTrk2[t2i].first &&
                    stubsTrk1[t1i].second == stubsTrk2[t2i].second)
                  nShareLay++;
              }
            }

            // Fill duplicate map
            if (nShareLay >= settings_.minIndStubs()) {  // For number of shared stub merge condition
              dupMap.at(itrk).at(jtrk) = true;
              dupMap.at(jtrk).at(itrk) = true;
            }
          }
        }

        // Check to see if the track is a duplicate
        for (unsigned int itrk = 0; itrk < numStublists; itrk++) {
          for (unsigned int jtrk = 0; jtrk < numStublists; jtrk++) {
            if (dupMap.at(itrk).at(jtrk)) {
              noMerge.at(itrk) = true;
            }
          }
        }

        // 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
        for (unsigned int itrk = 0; itrk < numStublists; itrk++) {
          if (!noMerge.at(itrk)) {
            if (((findOverlapRinvBins(inputtracklets_[itrk]).size() > 1) &&
                 (findRinvBin(inputtracklets_[itrk]) != bin)) ||
                ((findOverlapPhiBins(inputtracklets_[itrk]).size() > 1) &&
                 findPhiBin(inputtracklets_[itrk]) != phiBin)) {
              trackInfo[itrk].second = true;
            }
          }
        }
        // Merge duplicate tracks
        for (unsigned int itrk = 0; itrk < numStublists - 1; itrk++) {
          for (unsigned int jtrk = itrk + 1; jtrk < numStublists; jtrk++) {
            // Merge a track with its first duplicate found.
            if (dupMap.at(itrk).at(jtrk)) {
              // Set preferred track based on seed rank
              int preftrk;
              int rejetrk;
              if (seedRank[itrk] < seedRank[jtrk]) {
                preftrk = itrk;
                rejetrk = jtrk;
              } else {
                preftrk = jtrk;
                rejetrk = itrk;
              }

              // If the preffered track is in more than one bin, but not in the proper rinv or phi bin, then mark as true
              if (((findOverlapRinvBins(inputtracklets_[preftrk]).size() > 1) &&
                   (findRinvBin(inputtracklets_[preftrk]) != bin)) ||
                  ((findOverlapPhiBins(inputtracklets_[preftrk]).size() > 1) &&
                   (findPhiBin(inputtracklets_[preftrk]) != phiBin))) {
                trackBinInfo[preftrk] = true;
                trackBinInfo[rejetrk] = true;
              } else {
                // Get a merged stub list
                std::vector<const Stub*> newStubList;
                std::vector<const Stub*> stubsTrk1 = inputstublists_[preftrk];
                std::vector<const Stub*> stubsTrk2 = inputstublists_[rejetrk];
                std::vector<unsigned int> stubsTrk1indices;
                std::vector<unsigned int> stubsTrk2indices;
                for (unsigned int stub1it = 0; stub1it < stubsTrk1.size(); stub1it++) {
                  stubsTrk1indices.push_back(stubsTrk1[stub1it]->l1tstub()->uniqueIndex());
                }
                for (unsigned int stub2it = 0; stub2it < stubsTrk2.size(); stub2it++) {
                  stubsTrk2indices.push_back(stubsTrk2[stub2it]->l1tstub()->uniqueIndex());
                }
                newStubList = stubsTrk1;
                for (unsigned int stub2it = 0; stub2it < stubsTrk2.size(); stub2it++) {
                  if (find(stubsTrk1indices.begin(), stubsTrk1indices.end(), stubsTrk2indices[stub2it]) ==
                      stubsTrk1indices.end()) {
                    newStubList.push_back(stubsTrk2[stub2it]);
                  }
                }
                //   Overwrite stublist of preferred track with merged list
                inputstublists_[preftrk] = newStubList;

                std::vector<std::pair<int, int>> newStubidsList;
                std::vector<std::pair<int, int>> stubidsTrk1 = mergedstubidslists_[preftrk];
                std::vector<std::pair<int, int>> stubidsTrk2 = mergedstubidslists_[rejetrk];
                newStubidsList = stubidsTrk1;

                for (unsigned int stub2it = 0; stub2it < stubsTrk2.size(); stub2it++) {
                  if (find(stubsTrk1indices.begin(), stubsTrk1indices.end(), stubsTrk2indices[stub2it]) ==
                      stubsTrk1indices.end()) {
                    newStubidsList.push_back(stubidsTrk2[stub2it]);
                  }
                }
                // Overwrite stubidslist of preferred track with merged list
                mergedstubidslists_[preftrk] = newStubidsList;

                // Mark that rejected track has been merged into another track
                trackInfo[rejetrk].second = true;
              }
            }
          }
        }

        for (unsigned int ktrk = 0; ktrk < numStublists; ktrk++) {
          if ((trackInfo[ktrk].second != true) && (trackBinInfo[ktrk] != true)) {
            prefTracks.push_back(ktrk);
            prefTrackFit.push_back(trackInfo[ktrk].first);
            inputtrackletsall.push_back(inputtracklets_[ktrk]);
            inputstublistsall.push_back(inputstublists_[ktrk]);
            inputstubidslistsall.push_back(inputstubidslists_[ktrk]);
            mergedstubidslistsall.push_back(mergedstubidslists_[ktrk]);
          }
        }

        // Need to clear all the vectors which will be used in the next bin
        seedRank.clear();
        trackInfo.clear();
        trackBinInfo.clear();
        inputtracklets_.clear();
        inputstublists_.clear();
        inputstubidslists_.clear();
        mergedstubidslists_.clear();
      }
    }

    // Make the final track objects, fit with KF, and send to output
    for (unsigned int itrk = 0; itrk < prefTracks.size(); itrk++) {
      Tracklet* tracklet = inputtrackletsall[itrk];
      std::vector<const Stub*> trackstublist = inputstublistsall[itrk];

      // Run KF track fit
      HybridFit hybridFitter(iSector, settings_, globals_);
      hybridFitter.Fit(tracklet, trackstublist);

      // If the track was accepted (and thus fit), add to output
      if (tracklet->fit()) {
        // Add fitted Track to output (later converted to TTTrack)
        Track* outtrack = tracklet->getTrack();
        outtrack->setSector(iSector);
        // Also store fitted track as more detailed Tracklet object.
        outputtracklets_[prefTrackFit[itrk]]->addTrack(tracklet);

        // Add all tracks to standalone root file output
        outtrack->setStubIDpremerge(inputstubidslistsall[itrk]);
        outtrack->setStubIDprefit(mergedstubidslistsall[itrk]);
        outputtracks.push_back(*outtrack);
      }
    }
  }

#endif

  // Grid removal //
  if (settings_.removalType() == "grid") {
    // Sort tracks by ichisq/DoF so that removal will keep the lower ichisq/DoF track
    std::sort(inputtracks_.begin(), inputtracks_.end(), [](const Track* lhs, const Track* rhs) {
      return lhs->ichisq() / lhs->stubID().size() < rhs->ichisq() / rhs->stubID().size();
    });
    bool grid[35][40] = {{false}};

    for (unsigned int itrk = 0; itrk < numTrk; itrk++) {
      if (inputtracks_[itrk]->duplicate())
        edm::LogPrint("Tracklet") << "WARNING: Track already tagged as duplicate!!";

      double phiBin = (inputtracks_[itrk]->phi0(settings_) - 2 * M_PI / 27 * iSector) / (2 * M_PI / 9 / 50) + 9;
      phiBin = std::max(phiBin, 0.);
      phiBin = std::min(phiBin, 34.);

      double ptBin = 1 / inputtracks_[itrk]->pt(settings_) * 40 + 20;
      ptBin = std::max(ptBin, 0.);
      ptBin = std::min(ptBin, 39.);

      if (grid[(int)phiBin][(int)ptBin])
        inputtracks_[itrk]->setDuplicate(true);
      grid[(int)phiBin][(int)ptBin] = true;

      double phiTest = inputtracks_[itrk]->phi0(settings_) - 2 * M_PI / 27 * iSector;
      if (phiTest < -2 * M_PI / 27)
        edm::LogVerbatim("Tracklet") << "track phi too small!";
      if (phiTest > 2 * 2 * M_PI / 27)
        edm::LogVerbatim("Tracklet") << "track phi too big!";
    }
  }  // end grid removal

  // ichi + nstub removal //
  if (settings_.removalType() == "ichi" || settings_.removalType() == "nstub") {
    for (unsigned int itrk = 0; itrk < numTrk - 1; itrk++) {  // numTrk-1 since last track has no other to compare to

      // If primary track is a duplicate, it cannot veto any...move on
      if (inputtracks_[itrk]->duplicate() == 1)
        continue;

      unsigned int nStubP = 0;
      vector<unsigned int> nStubS(numTrk);
      vector<unsigned int> nShare(numTrk);
      // Get and count primary stubs
      std::map<int, int> stubsTrk1 = inputtracks_[itrk]->stubID();
      nStubP = stubsTrk1.size();

      for (unsigned int jtrk = itrk + 1; jtrk < numTrk; jtrk++) {
        // Skip duplicate tracks
        if (inputtracks_[jtrk]->duplicate() == 1)
          continue;

        // Get and count secondary stubs
        std::map<int, int> stubsTrk2 = inputtracks_[jtrk]->stubID();
        nStubS[jtrk] = stubsTrk2.size();

        // Count shared stubs
        for (auto& st : stubsTrk1) {
          if (stubsTrk2.find(st.first) != stubsTrk2.end()) {
            if (st.second == stubsTrk2[st.first])
              nShare[jtrk]++;
          }
        }
      }

      // Tag duplicates
      for (unsigned int jtrk = itrk + 1; jtrk < numTrk; jtrk++) {
        // Skip duplicate tracks
        if (inputtracks_[jtrk]->duplicate() == 1)
          continue;

        // Chi2 duplicate removal
        if (settings_.removalType() == "ichi") {
          if ((nStubP - nShare[jtrk] < settings_.minIndStubs()) ||
              (nStubS[jtrk] - nShare[jtrk] < settings_.minIndStubs())) {
            if ((int)inputtracks_[itrk]->ichisq() / (2 * inputtracks_[itrk]->stubID().size() - 4) >
                (int)inputtracks_[jtrk]->ichisq() / (2 * inputtracks_[itrk]->stubID().size() - 4)) {
              inputtracks_[itrk]->setDuplicate(true);
            } else if ((int)inputtracks_[itrk]->ichisq() / (2 * inputtracks_[itrk]->stubID().size() - 4) <=
                       (int)inputtracks_[jtrk]->ichisq() / (2 * inputtracks_[itrk]->stubID().size() - 4)) {
              inputtracks_[jtrk]->setDuplicate(true);
            } else {
              edm::LogVerbatim("Tracklet") << "Error: Didn't tag either track in duplicate pair.";
            }
          }
        }  // end ichi removal

        // nStub duplicate removal
        if (settings_.removalType() == "nstub") {
          if ((nStubP - nShare[jtrk] < settings_.minIndStubs()) && (nStubP < nStubS[jtrk])) {
            inputtracks_[itrk]->setDuplicate(true);
          } else if ((nStubS[jtrk] - nShare[jtrk] < settings_.minIndStubs()) && (nStubS[jtrk] <= nStubP)) {
            inputtracks_[jtrk]->setDuplicate(true);
          } else {
            edm::LogVerbatim("Tracklet") << "Error: Didn't tag either track in duplicate pair.";
          }
        }  // end nstub removal

      }  // end tag duplicates

    }  // end loop over primary track

  }  // end ichi + nstub removal

  //Add tracks to output
  if (settings_.removalType() != "merge") {
    for (unsigned int i = 0; i < inputtrackfits_.size(); i++) {
      for (unsigned int j = 0; j < inputtrackfits_[i]->nTracks(); j++) {
        if (inputtrackfits_[i]->getTrack(j)->getTrack()->duplicate() == 0) {
          if (settings_.writeMonitorData("Seeds")) {
            ofstream fout("seeds.txt", ofstream::app);
            fout << __FILE__ << ":" << __LINE__ << " " << name_ << "_" << iSector << " "
                 << inputtrackfits_[i]->getTrack(j)->getISeed() << endl;
            fout.close();
          }
          outputtracklets_[i]->addTrack(inputtrackfits_[i]->getTrack(j));
        }
        //For root file:
        outputtracks.push_back(*inputtrackfits_[i]->getTrack(j)->getTrack());
      }
    }
  }
}

double PurgeDuplicate::getPhiRes(Tracklet* curTracklet, const Stub* curStub) const {
  double phiproj;
  double stubphi;
  double phires;
  // Get phi position of stub
  stubphi = curStub->l1tstub()->phi();
  // Get layer that the stub is in (Layer 1->6, Disk 1->5)
  int Layer = curStub->layerdisk() + 1;
  if (Layer > N_LAYER) {
    Layer = 0;
  }
  int Disk = curStub->layerdisk() - (N_LAYER - 1);
  if (Disk < 0) {
    Disk = 0;
  }
  // Get phi projection of tracklet
  int seedindex = curTracklet->seedIndex();
  // If this stub is a seed stub, set projection=phi, so that res=0
  if (isSeedingStub(seedindex, Layer, Disk)) {
    phiproj = stubphi;
    // Otherwise, get projection of tracklet
  } else {
    phiproj = curTracklet->proj(curStub->layerdisk()).phiproj();
  }
  // Calculate residual
  phires = std::abs(stubphi - phiproj);
  return phires;
}

bool PurgeDuplicate::isSeedingStub(int seedindex, int Layer, int Disk) const {
  if ((seedindex == 0 && (Layer == 1 || Layer == 2)) || (seedindex == 1 && (Layer == 2 || Layer == 3)) ||
      (seedindex == 2 && (Layer == 3 || Layer == 4)) || (seedindex == 3 && (Layer == 5 || Layer == 6)) ||
      (seedindex == 4 && (abs(Disk) == 1 || abs(Disk) == 2)) ||
      (seedindex == 5 && (abs(Disk) == 3 || abs(Disk) == 4)) || (seedindex == 6 && (Layer == 1 || abs(Disk) == 1)) ||
      (seedindex == 7 && (Layer == 2 || abs(Disk) == 1)) ||
      (seedindex == 8 && (Layer == 2 || Layer == 3 || Layer == 4)) ||
      (seedindex == 9 && (Layer == 4 || Layer == 5 || Layer == 6)) ||
      (seedindex == 10 && (Layer == 2 || Layer == 3 || abs(Disk) == 1)) ||
      (seedindex == 11 && (Layer == 2 || abs(Disk) == 1 || abs(Disk) == 2)))
    return true;

  return false;
}

std::pair<int, int> PurgeDuplicate::findLayerDisk(const Stub* st) const {
  std::pair<int, int> layer_disk;
  layer_disk.first = st->layerdisk() + 1;
  if (layer_disk.first > N_LAYER) {
    layer_disk.first = 0;
  }
  layer_disk.second = st->layerdisk() - (N_LAYER - 1);
  if (layer_disk.second < 0) {
    layer_disk.second = 0;
  }
  return layer_disk;
}

std::string PurgeDuplicate::l1tinfo(const L1TStub* l1stub, std::string str = "") const {
  std::string thestr = Form("\t %s stub info:  r/z/phi:\t%f\t%f\t%f\t%d\t%f\t%d",
                            str.c_str(),
                            l1stub->r(),
                            l1stub->z(),
                            l1stub->phi(),
                            l1stub->iphi(),
                            l1stub->bend(),
                            l1stub->layerdisk());
  return thestr;
}

std::vector<double> PurgeDuplicate::getInventedCoords(unsigned int iSector,
                                                      const Stub* st,
                                                      const Tracklet* tracklet) const {
  int stubLayer = (findLayerDisk(st)).first;
  int stubDisk = (findLayerDisk(st)).second;

  double stub_phi = -99;
  double stub_z = -99;
  double stub_r = -99;

  double tracklet_rinv = tracklet->rinv();

  if (st->isBarrel()) {
    stub_r = settings_.rmean(stubLayer - 1);
    stub_phi = tracklet->phi0() - std::asin(stub_r * tracklet_rinv / 2);
    stub_phi = stub_phi + iSector * settings_.dphisector() - 0.5 * settings_.dphisectorHG();
    stub_phi = reco::reduceRange(stub_phi);
    stub_z = tracklet->z0() + 2 * tracklet->t() * 1 / tracklet_rinv * std::asin(stub_r * tracklet_rinv / 2);
  } else {
    stub_z = settings_.zmean(stubDisk - 1) * tracklet->disk() / abs(tracklet->disk());
    stub_phi = tracklet->phi0() - (stub_z - tracklet->z0()) * tracklet_rinv / 2 / tracklet->t();
    stub_phi = stub_phi + iSector * settings_.dphisector() - 0.5 * settings_.dphisectorHG();
    stub_phi = reco::reduceRange(stub_phi);
    stub_r = 2 / tracklet_rinv * std::sin((stub_z - tracklet->z0()) * tracklet_rinv / 2 / tracklet->t());
  }

  std::vector<double> invented_coords{stub_r, stub_z, stub_phi};
  return invented_coords;
}

std::vector<double> PurgeDuplicate::getInventedCoordsExtended(unsigned int iSector,
                                                              const Stub* st,
                                                              const Tracklet* tracklet) const {
  int stubLayer = (findLayerDisk(st)).first;
  int stubDisk = (findLayerDisk(st)).second;

  double stub_phi = -99;
  double stub_z = -99;
  double stub_r = -99;

  double rho = 1 / tracklet->rinv();
  double rho_minus_d0 = rho + tracklet->d0();  // should be -, but otherwise does not work

  // exact helix
  if (st->isBarrel()) {
    stub_r = settings_.rmean(stubLayer - 1);

    double sin_val = (stub_r * stub_r + rho_minus_d0 * rho_minus_d0 - rho * rho) / (2 * stub_r * rho_minus_d0);
    stub_phi = tracklet->phi0() - std::asin(sin_val);
    stub_phi = stub_phi + iSector * settings_.dphisector() - 0.5 * settings_.dphisectorHG();
    stub_phi = reco::reduceRange(stub_phi);

    double beta = std::acos((rho * rho + rho_minus_d0 * rho_minus_d0 - stub_r * stub_r) / (2 * rho * rho_minus_d0));
    stub_z = tracklet->z0() + tracklet->t() * std::abs(rho * beta);
  } else {
    stub_z = settings_.zmean(stubDisk - 1) * tracklet->disk() / abs(tracklet->disk());

    double beta = (stub_z - tracklet->z0()) / (tracklet->t() * std::abs(rho));  // maybe rho should be abs value
    double r_square = -2 * rho * rho_minus_d0 * std::cos(beta) + rho * rho + rho_minus_d0 * rho_minus_d0;
    stub_r = sqrt(r_square);

    double sin_val = (stub_r * stub_r + rho_minus_d0 * rho_minus_d0 - rho * rho) / (2 * stub_r * rho_minus_d0);
    stub_phi = tracklet->phi0() - std::asin(sin_val);
    stub_phi = stub_phi + iSector * settings_.dphisector() - 0.5 * settings_.dphisectorHG();
    stub_phi = reco::reduceRange(stub_phi);
  }

  // TMP: for displaced tracking, exclude one of the 3 seeding stubs
  // to be discussed
  int seed = tracklet->seedIndex();
  if ((seed == 8 && stubLayer == 4) || (seed == 9 && stubLayer == 5) || (seed == 10 && stubLayer == 3) ||
      (seed == 11 && abs(stubDisk) == 1)) {
    stub_phi = st->l1tstub()->phi();
    stub_z = st->l1tstub()->z();
    stub_r = st->l1tstub()->r();
  }

  std::vector<double> invented_coords{stub_r, stub_z, stub_phi};
  return invented_coords;
}

std::vector<const Stub*> PurgeDuplicate::getInventedSeedingStub(
    unsigned int iSector, const Tracklet* tracklet, const std::vector<const Stub*>& originalStubsList) const {
  std::vector<const Stub*> newStubList;

  for (unsigned int stubit = 0; stubit < originalStubsList.size(); stubit++) {
    const Stub* thisStub = originalStubsList[stubit];

    if (isSeedingStub(tracklet->seedIndex(), (findLayerDisk(thisStub)).first, (findLayerDisk(thisStub)).second)) {
      // get a vector containing r, z, phi
      std::vector<double> inv_r_z_phi;
      if (!settings_.extended())
        inv_r_z_phi = getInventedCoords(iSector, thisStub, tracklet);
      else {
        inv_r_z_phi = getInventedCoordsExtended(iSector, thisStub, tracklet);
      }
      double stub_x_invent = inv_r_z_phi[0] * std::cos(inv_r_z_phi[2]);
      double stub_y_invent = inv_r_z_phi[0] * std::sin(inv_r_z_phi[2]);
      double stub_z_invent = inv_r_z_phi[1];

      Stub* invent_stub_ptr = new Stub(*thisStub);
      const L1TStub* l1stub = thisStub->l1tstub();
      L1TStub invent_l1stub = *l1stub;
      invent_l1stub.setCoords(stub_x_invent, stub_y_invent, stub_z_invent);

      invent_stub_ptr->setl1tstub(new L1TStub(invent_l1stub));
      invent_stub_ptr->l1tstub()->setAllStubIndex(l1stub->allStubIndex());
      invent_stub_ptr->l1tstub()->setUniqueIndex(l1stub->uniqueIndex());

      newStubList.push_back(invent_stub_ptr);

    } else {
      newStubList.push_back(thisStub);
    }
  }
  return newStubList;
}

// Tells us the variable bin to which a track would belong
unsigned int PurgeDuplicate::findRinvBin(const Tracklet* trk) const {
  std::vector<double> rinvBins = settings_.rinvBins();

  //Get rinverse of track
  double rinv = trk->rinv();

  //Check between what 2 values in rinvbins rinv is between
  auto bins = std::upper_bound(rinvBins.begin(), rinvBins.end(), rinv);

  //return integer for bin index
  unsigned int rIndx = std::distance(rinvBins.begin(), bins);
  if (rIndx == std::distance(rinvBins.end(), bins))
    return rinvBins.size() - 2;
  else if (bins == rinvBins.begin())
    return std::distance(rinvBins.begin(), bins);
  else
    return rIndx - 1;
}

// Tells us the phi bin to which a track would belong
unsigned int PurgeDuplicate::findPhiBin(const Tracklet* trk) const {
  std::vector<double> phiBins = settings_.phiBins();

  //Get phi of track at rcrit
  double phi0 = trk->phi0();
  double rcrit = settings_.rcrit();
  double rinv = trk->rinv();
  double phi = phi0 - asin(0.5 * rinv * rcrit);

  //Check between what 2 values in phibins phi is between
  auto bins = std::upper_bound(phiBins.begin(), phiBins.end(), phi);

  //return integer for bin index
  unsigned int phiIndx = std::distance(phiBins.begin(), bins);
  if (phiIndx == std::distance(phiBins.end(), bins))
    return phiBins.size() - 2;
  else if (bins == phiBins.begin())
    return std::distance(phiBins.begin(), bins);
  else
    return phiIndx - 1;
}

// Tells us the overlap rinv bin(s) to which a track belongs
std::vector<unsigned int> PurgeDuplicate::findOverlapRinvBins(const Tracklet* trk) const {
  double rinv = trk->rinv();

  const double rinvOverlapSize = settings_.rinvOverlapSize();
  const std::vector<double>& rinvBins = settings_.rinvBins();

  std::vector<unsigned int> chosenBins;
  for (long unsigned int i = 0; i < rinvBins.size() - 1; i++) {
    if ((rinv < rinvBins[i + 1] + rinvOverlapSize) && (rinv > rinvBins[i] - rinvOverlapSize)) {
      chosenBins.push_back(i);
    }
  }
  return chosenBins;
}

// Tells us the overlap phi bin(s) to which a track belongs
std::vector<unsigned int> PurgeDuplicate::findOverlapPhiBins(const Tracklet* trk) const {
  double phi0 = trk->phi0();
  double rcrit = settings_.rcrit();
  double rinv = trk->rinv();
  double phi = phi0 - asin(0.5 * rinv * rcrit);

  const double phiOverlapSize = settings_.phiOverlapSize();
  const std::vector<double>& phiBins = settings_.phiBins();

  std::vector<unsigned int> chosenBins;
  for (long unsigned int i = 0; i < phiBins.size() - 1; i++) {
    if ((phi < phiBins[i + 1] + phiOverlapSize) && (phi > phiBins[i] - phiOverlapSize)) {
      chosenBins.push_back(i);
    }
  }
  return chosenBins;
}

// Tells us if a track is in the current bin
bool PurgeDuplicate::isTrackInBin(const std::vector<unsigned int>& vec, unsigned int num) const {
  auto result = std::find(vec.begin(), vec.end(), num);
  bool found = (result != vec.end());
  return found;
}