L2MuonSeedGeneratorFromL1T.cc

Go to the documentation of this file.

//-------------------------------------------------
//
//
//--------------------------------------------------

// Class Header
#include "RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.h"

// Framework
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"

#include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
#include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"

#include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
#include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"

using namespace std;
using namespace edm;
using namespace l1t;

//--- Functions used in output sorting
bool sortByL1Pt(L2MuonTrajectorySeed &A, L2MuonTrajectorySeed &B) {
  l1t::MuonRef Ref_L1A = A.l1tParticle();
  l1t::MuonRef Ref_L1B = B.l1tParticle();
  return (Ref_L1A->pt() > Ref_L1B->pt());
};

bool sortByL1QandPt(L2MuonTrajectorySeed &A, L2MuonTrajectorySeed &B) {
  l1t::MuonRef Ref_L1A = A.l1tParticle();
  l1t::MuonRef Ref_L1B = B.l1tParticle();

  // Compare quality first
  if (Ref_L1A->hwQual() > Ref_L1B->hwQual())
    return true;
  if (Ref_L1A->hwQual() < Ref_L1B->hwQual())
    return false;

  // For same quality L1s compare pT
  return (Ref_L1A->pt() > Ref_L1B->pt());
};

// constructors
L2MuonSeedGeneratorFromL1T::L2MuonSeedGeneratorFromL1T(const edm::ParameterSet &iConfig)
    : theSource(iConfig.getParameter<InputTag>("InputObjects")),
      theL1GMTReadoutCollection(iConfig.getParameter<InputTag>("GMTReadoutCollection")),  // to be removed
      thePropagatorName(iConfig.getParameter<string>("Propagator")),
      theL1MinPt(iConfig.getParameter<double>("L1MinPt")),
      theL1MaxEta(iConfig.getParameter<double>("L1MaxEta")),
      theL1MinQuality(iConfig.getParameter<unsigned int>("L1MinQuality")),
      theMinPtBarrel(iConfig.getParameter<double>("SetMinPtBarrelTo")),
      theMinPtEndcap(iConfig.getParameter<double>("SetMinPtEndcapTo")),
      useOfflineSeed(iConfig.getUntrackedParameter<bool>("UseOfflineSeed", false)),
      useUnassociatedL1(iConfig.getParameter<bool>("UseUnassociatedL1")),
      matchingDR(iConfig.getParameter<std::vector<double>>("MatchDR")),
      etaBins(iConfig.getParameter<std::vector<double>>("EtaMatchingBins")),
      centralBxOnly_(iConfig.getParameter<bool>("CentralBxOnly")),
      matchType(iConfig.getParameter<unsigned int>("MatchType")),
      sortType(iConfig.getParameter<unsigned int>("SortType")) {
  muCollToken_ = consumes<MuonBxCollection>(theSource);

  if (useOfflineSeed) {
    theOfflineSeedLabel = iConfig.getUntrackedParameter<InputTag>("OfflineSeedLabel");
    offlineSeedToken_ = consumes<edm::View<TrajectorySeed>>(theOfflineSeedLabel);

    // check that number of eta bins -1 matches number of dR cones
    if (matchingDR.size() != etaBins.size() - 1) {
      throw cms::Exception("Configuration") << "Size of MatchDR "
                                            << "does not match number of eta bins." << endl;
    }
  }

  if (matchType > 4 || sortType > 3) {
    throw cms::Exception("Configuration") << "Wrong MatchType or SortType" << endl;
  }

  // service parameters
  ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");

  // the services
  theService = new MuonServiceProxy(serviceParameters);

  // the estimator
  theEstimator = new Chi2MeasurementEstimator(10000.);

  produces<L2MuonTrajectorySeedCollection>();
}

// destructor
L2MuonSeedGeneratorFromL1T::~L2MuonSeedGeneratorFromL1T() {
  if (theService)
    delete theService;
  if (theEstimator)
    delete theEstimator;
}

void L2MuonSeedGeneratorFromL1T::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("GMTReadoutCollection", edm::InputTag(""));  // to be removed
  desc.add<edm::InputTag>("InputObjects", edm::InputTag("hltGmtStage2Digis"));
  desc.add<string>("Propagator", "");
  desc.add<double>("L1MinPt", -1.);
  desc.add<double>("L1MaxEta", 5.0);
  desc.add<unsigned int>("L1MinQuality", 0);
  desc.add<double>("SetMinPtBarrelTo", 3.5);
  desc.add<double>("SetMinPtEndcapTo", 1.0);
  desc.addUntracked<bool>("UseOfflineSeed", false);
  desc.add<bool>("UseUnassociatedL1", true);
  desc.add<std::vector<double>>("MatchDR", {0.3});
  desc.add<std::vector<double>>("EtaMatchingBins", {0., 2.5});
  desc.add<bool>("CentralBxOnly", true);
  desc.add<unsigned int>("MatchType", 0)
      ->setComment(
          "MatchType : 0 Old matching,   1 L1 Order(1to1), 2 Min dR(1to1), 3 Higher Q(1to1), 4 All matched L1");
  desc.add<unsigned int>("SortType", 0)->setComment("SortType :  0 not sort,       1 Pt, 2 Q and Pt");
  desc.addUntracked<edm::InputTag>("OfflineSeedLabel", edm::InputTag(""));

  edm::ParameterSetDescription psd0;
  psd0.addUntracked<std::vector<std::string>>("Propagators", {"SteppingHelixPropagatorAny"});
  psd0.add<bool>("RPCLayers", true);
  psd0.addUntracked<bool>("UseMuonNavigation", true);
  desc.add<edm::ParameterSetDescription>("ServiceParameters", psd0);
  descriptions.add("L2MuonSeedGeneratorFromL1T", desc);
}

void L2MuonSeedGeneratorFromL1T::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
  const std::string metname = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
  MuonPatternRecoDumper debug;

  auto output = std::make_unique<L2MuonTrajectorySeedCollection>();

  edm::Handle<MuonBxCollection> muColl;
  iEvent.getByToken(muCollToken_, muColl);
  LogDebug(metname) << "Number of muons " << muColl->size() << endl;

  //--- matchType 0 : Old logic
  if (matchType == 0) {
    edm::Handle<edm::View<TrajectorySeed>> offlineSeedHandle;
    vector<int> offlineSeedMap;
    if (useOfflineSeed) {
      iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
      LogDebug(metname) << "Number of offline seeds " << offlineSeedHandle->size() << endl;
      offlineSeedMap = vector<int>(offlineSeedHandle->size(), 0);
    }

    for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) {
      if (centralBxOnly_ && (ibx != 0))
        continue;
      for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++) {
        unsigned int quality = it->hwQual();
        int valid_charge = it->hwChargeValid();

        float pt = it->pt();
        float eta = it->eta();
        float theta = 2 * atan(exp(-eta));
        float phi = it->phi();
        int charge = it->charge();
        // Set charge=0 for the time being if the valid charge bit is zero
        if (!valid_charge)
          charge = 0;

        // link number indices of the optical fibres that connect the uGMT with the track finders
        //EMTF+ : 36-41, OMTF+ : 42-47, BMTF : 48-59, OMTF- : 60-65, EMTF- : 66-71
        int link = 36 + (int)(it->tfMuonIndex() / 3.);
        bool barrel = true;
        if ((link >= 36 && link <= 41) || (link >= 66 && link <= 71))
          barrel = false;

        if (pt < theL1MinPt || fabs(eta) > theL1MaxEta)
          continue;

        LogDebug(metname) << "New L2 Muon Seed";
        LogDebug(metname) << "Pt = " << pt << " GeV/c";
        LogDebug(metname) << "eta = " << eta;
        LogDebug(metname) << "theta = " << theta << " rad";
        LogDebug(metname) << "phi = " << phi << " rad";
        LogDebug(metname) << "charge = " << charge;
        LogDebug(metname) << "In Barrel? = " << barrel;

        if (quality <= theL1MinQuality)
          continue;
        LogDebug(metname) << "quality = " << quality;

        // Update the services
        theService->update(iSetup);

        const DetLayer *detLayer = nullptr;
        float radius = 0.;

        CLHEP::Hep3Vector vec(0., 1., 0.);
        vec.setTheta(theta);
        vec.setPhi(phi);

        DetId theid;
        // Get the det layer on which the state should be put
        if (barrel) {
          LogDebug(metname) << "The seed is in the barrel";

          // MB2
          theid = DTChamberId(0, 2, 0);
          detLayer = theService->detLayerGeometry()->idToLayer(theid);
          LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);

          const BoundSurface *sur = &(detLayer->surface());
          const BoundCylinder *bc = dynamic_cast<const BoundCylinder *>(sur);

          radius = fabs(bc->radius() / sin(theta));

          LogDebug(metname) << "radius " << radius;

          if (pt < theMinPtBarrel)
            pt = theMinPtBarrel;
        } else {
          LogDebug(metname) << "The seed is in the endcap";

          // ME2
          theid = theta < Geom::pi() / 2. ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);

          detLayer = theService->detLayerGeometry()->idToLayer(theid);
          LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);

          radius = fabs(detLayer->position().z() / cos(theta));

          if (pt < theMinPtEndcap)
            pt = theMinPtEndcap;
        }

        vec.setMag(radius);

        GlobalPoint pos(vec.x(), vec.y(), vec.z());

        GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));

        GlobalTrajectoryParameters param(pos, mom, charge, &*theService->magneticField());
        AlgebraicSymMatrix55 mat;

        mat[0][0] = (0.25 / pt) * (0.25 / pt);  // sigma^2(charge/abs_momentum)
        if (!barrel)
          mat[0][0] = (0.4 / pt) * (0.4 / pt);

        //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
        if (!valid_charge)
          mat[0][0] = (1. / pt) * (1. / pt);

        mat[1][1] = 0.05 * 0.05;  // sigma^2(lambda)
        mat[2][2] = 0.2 * 0.2;    // sigma^2(phi)
        mat[3][3] = 20. * 20.;    // sigma^2(x_transverse))
        mat[4][4] = 20. * 20.;    // sigma^2(y_transverse))

        CurvilinearTrajectoryError error(mat);

        const FreeTrajectoryState state(param, error);

        LogDebug(metname) << "Free trajectory State from the parameters";
        LogDebug(metname) << debug.dumpFTS(state);

        // Propagate the state on the MB2/ME2 surface
        TrajectoryStateOnSurface tsos =
            theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());

        LogDebug(metname) << "State after the propagation on the layer";
        LogDebug(metname) << debug.dumpLayer(detLayer);
        LogDebug(metname) << debug.dumpTSOS(tsos);

        double dRcone = matchingDR[0];
        if (fabs(eta) < etaBins.back()) {
          std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(eta));
          dRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
        }

        if (tsos.isValid()) {
          edm::OwnVector<TrackingRecHit> container;

          if (useOfflineSeed && (!valid_charge || charge == 0)) {
            const TrajectorySeed *assoOffseed =
                associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, tsos, dRcone);

            if (assoOffseed != nullptr) {
              PTrajectoryStateOnDet const &seedTSOS = assoOffseed->startingState();
              TrajectorySeed::const_iterator tsci = assoOffseed->recHits().first, tscie = assoOffseed->recHits().second;
              for (; tsci != tscie; ++tsci) {
                container.push_back(*tsci);
              }
              output->push_back(
                  L2MuonTrajectorySeed(seedTSOS,
                                       container,
                                       alongMomentum,
                                       MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
            } else {
              if (useUnassociatedL1) {
                // convert the TSOS into a PTSOD
                PTrajectoryStateOnDet const &seedTSOS = trajectoryStateTransform::persistentState(tsos, theid.rawId());
                output->push_back(
                    L2MuonTrajectorySeed(seedTSOS,
                                         container,
                                         alongMomentum,
                                         MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
              }
            }
          } else if (useOfflineSeed && valid_charge) {
            // Get the compatible dets on the layer
            std::vector<pair<const GeomDet *, TrajectoryStateOnSurface>> detsWithStates =
                detLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);

            if (detsWithStates.empty() && barrel) {
              // Fallback solution using ME2, try again to propagate but using ME2 as reference
              DetId fallback_id;
              theta < Geom::pi() / 2. ? fallback_id = CSCDetId(1, 2, 0, 0, 0) : fallback_id = CSCDetId(2, 2, 0, 0, 0);
              const DetLayer *ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id);

              tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface());
              detsWithStates =
                  ME2DetLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
            }

            if (!detsWithStates.empty()) {
              TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
              const GeomDet *newTSOSDet = detsWithStates.front().first;

              LogDebug(metname) << "Most compatible det";
              LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());

              if (newTSOS.isValid()) {
                LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
                                  << ", phi=" << newTSOS.globalPosition().phi()
                                  << ", eta=" << newTSOS.globalPosition().eta() << ")";
                LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
                                  << ", phi=" << newTSOS.globalMomentum().phi()
                                  << ", eta=" << newTSOS.globalMomentum().eta() << ")";

                const TrajectorySeed *assoOffseed =
                    associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS, dRcone);

                if (assoOffseed != nullptr) {
                  PTrajectoryStateOnDet const &seedTSOS = assoOffseed->startingState();
                  TrajectorySeed::const_iterator tsci = assoOffseed->recHits().first,
                                                 tscie = assoOffseed->recHits().second;
                  for (; tsci != tscie; ++tsci) {
                    container.push_back(*tsci);
                  }
                  output->push_back(
                      L2MuonTrajectorySeed(seedTSOS,
                                           container,
                                           alongMomentum,
                                           MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
                } else {
                  if (useUnassociatedL1) {
                    // convert the TSOS into a PTSOD
                    PTrajectoryStateOnDet const &seedTSOS =
                        trajectoryStateTransform::persistentState(newTSOS, newTSOSDet->geographicalId().rawId());
                    output->push_back(
                        L2MuonTrajectorySeed(seedTSOS,
                                             container,
                                             alongMomentum,
                                             MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
                  }
                }
              }
            }
          } else {
            // convert the TSOS into a PTSOD
            PTrajectoryStateOnDet const &seedTSOS = trajectoryStateTransform::persistentState(tsos, theid.rawId());
            output->push_back(L2MuonTrajectorySeed(
                seedTSOS, container, alongMomentum, MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
          }
        }
      }
    }

  }  // End of matchType 0

  //--- matchType > 0 : Updated logics
  else if (matchType > 0) {
    unsigned int nMuColl = muColl->size();

    vector<vector<double>> dRmtx;
    vector<vector<const TrajectorySeed *>> selOffseeds;

    edm::Handle<edm::View<TrajectorySeed>> offlineSeedHandle;
    if (useOfflineSeed) {
      iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
      unsigned int nOfflineSeed = offlineSeedHandle->size();
      LogDebug(metname) << "Number of offline seeds " << nOfflineSeed << endl;

      // Initialize dRmtx and selOffseeds
      dRmtx = vector<vector<double>>(nMuColl, vector<double>(nOfflineSeed, 999.0));
      selOffseeds =
          vector<vector<const TrajectorySeed *>>(nMuColl, vector<const TrajectorySeed *>(nOfflineSeed, nullptr));
    }

    // At least one L1 muons are associated with offSeed
    bool isAsso = false;

    //--- Fill dR matrix
    for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) {
      if (centralBxOnly_ && (ibx != 0))
        continue;
      for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++) {
        //Position of given L1mu
        unsigned int imu = distance(muColl->begin(muColl->getFirstBX()), it);

        unsigned int quality = it->hwQual();
        int valid_charge = it->hwChargeValid();

        float pt = it->pt();
        float eta = it->eta();
        float theta = 2 * atan(exp(-eta));
        float phi = it->phi();
        int charge = it->charge();
        // Set charge=0 for the time being if the valid charge bit is zero
        if (!valid_charge)
          charge = 0;

        // link number indices of the optical fibres that connect the uGMT with the track finders
        //EMTF+ : 36-41, OMTF+ : 42-47, BMTF : 48-59, OMTF- : 60-65, EMTF- : 66-71
        int link = 36 + (int)(it->tfMuonIndex() / 3.);
        bool barrel = true;
        if ((link >= 36 && link <= 41) || (link >= 66 && link <= 71))
          barrel = false;

        if (pt < theL1MinPt || fabs(eta) > theL1MaxEta)
          continue;

        LogDebug(metname) << "New L2 Muon Seed";
        LogDebug(metname) << "Pt = " << pt << " GeV/c";
        LogDebug(metname) << "eta = " << eta;
        LogDebug(metname) << "theta = " << theta << " rad";
        LogDebug(metname) << "phi = " << phi << " rad";
        LogDebug(metname) << "charge = " << charge;
        LogDebug(metname) << "In Barrel? = " << barrel;

        if (quality <= theL1MinQuality)
          continue;
        LogDebug(metname) << "quality = " << quality;

        // Update the services
        theService->update(iSetup);

        const DetLayer *detLayer = nullptr;
        float radius = 0.;

        CLHEP::Hep3Vector vec(0., 1., 0.);
        vec.setTheta(theta);
        vec.setPhi(phi);

        DetId theid;
        // Get the det layer on which the state should be put
        if (barrel) {
          LogDebug(metname) << "The seed is in the barrel";

          // MB2
          theid = DTChamberId(0, 2, 0);
          detLayer = theService->detLayerGeometry()->idToLayer(theid);
          LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);

          const BoundSurface *sur = &(detLayer->surface());
          const BoundCylinder *bc = dynamic_cast<const BoundCylinder *>(sur);

          radius = fabs(bc->radius() / sin(theta));

          LogDebug(metname) << "radius " << radius;

          if (pt < theMinPtBarrel)
            pt = theMinPtBarrel;
        } else {
          LogDebug(metname) << "The seed is in the endcap";

          // ME2
          theid = theta < Geom::pi() / 2. ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);

          detLayer = theService->detLayerGeometry()->idToLayer(theid);
          LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);

          radius = fabs(detLayer->position().z() / cos(theta));

          if (pt < theMinPtEndcap)
            pt = theMinPtEndcap;
        }

        vec.setMag(radius);

        GlobalPoint pos(vec.x(), vec.y(), vec.z());

        GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));

        GlobalTrajectoryParameters param(pos, mom, charge, &*theService->magneticField());
        AlgebraicSymMatrix55 mat;

        mat[0][0] = (0.25 / pt) * (0.25 / pt);  // sigma^2(charge/abs_momentum)
        if (!barrel)
          mat[0][0] = (0.4 / pt) * (0.4 / pt);

        //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
        if (!valid_charge)
          mat[0][0] = (1. / pt) * (1. / pt);

        mat[1][1] = 0.05 * 0.05;  // sigma^2(lambda)
        mat[2][2] = 0.2 * 0.2;    // sigma^2(phi)
        mat[3][3] = 20. * 20.;    // sigma^2(x_transverse))
        mat[4][4] = 20. * 20.;    // sigma^2(y_transverse))

        CurvilinearTrajectoryError error(mat);

        const FreeTrajectoryState state(param, error);

        LogDebug(metname) << "Free trajectory State from the parameters";
        LogDebug(metname) << debug.dumpFTS(state);

        // Propagate the state on the MB2/ME2 surface
        TrajectoryStateOnSurface tsos =
            theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());

        LogDebug(metname) << "State after the propagation on the layer";
        LogDebug(metname) << debug.dumpLayer(detLayer);
        LogDebug(metname) << debug.dumpTSOS(tsos);

        double dRcone = matchingDR[0];
        if (fabs(eta) < etaBins.back()) {
          std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(eta));
          dRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
        }

        if (tsos.isValid()) {
          edm::OwnVector<TrackingRecHit> container;

          if (useOfflineSeed) {
            if ((!valid_charge || charge == 0)) {
              bool isAssoOffseed = isAssociateOfflineSeedToL1(offlineSeedHandle, dRmtx, tsos, imu, selOffseeds, dRcone);

              if (isAssoOffseed) {
                isAsso = true;
              }

              // Using old way
              else {
                if (useUnassociatedL1) {
                  // convert the TSOS into a PTSOD
                  PTrajectoryStateOnDet const &seedTSOS =
                      trajectoryStateTransform::persistentState(tsos, theid.rawId());
                  output->push_back(
                      L2MuonTrajectorySeed(seedTSOS,
                                           container,
                                           alongMomentum,
                                           MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
                }
              }

            }

            else if (valid_charge) {
              // Get the compatible dets on the layer
              std::vector<pair<const GeomDet *, TrajectoryStateOnSurface>> detsWithStates =
                  detLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);

              if (detsWithStates.empty() && barrel) {
                // Fallback solution using ME2, try again to propagate but using ME2 as reference
                DetId fallback_id;
                theta < Geom::pi() / 2. ? fallback_id = CSCDetId(1, 2, 0, 0, 0) : fallback_id = CSCDetId(2, 2, 0, 0, 0);
                const DetLayer *ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id);

                tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface());
                detsWithStates =
                    ME2DetLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
              }

              if (!detsWithStates.empty()) {
                TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
                const GeomDet *newTSOSDet = detsWithStates.front().first;

                LogDebug(metname) << "Most compatible det";
                LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());

                if (newTSOS.isValid()) {
                  LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
                                    << ", phi=" << newTSOS.globalPosition().phi()
                                    << ", eta=" << newTSOS.globalPosition().eta() << ")";
                  LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
                                    << ", phi=" << newTSOS.globalMomentum().phi()
                                    << ", eta=" << newTSOS.globalMomentum().eta() << ")";

                  bool isAssoOffseed =
                      isAssociateOfflineSeedToL1(offlineSeedHandle, dRmtx, newTSOS, imu, selOffseeds, dRcone);

                  if (isAssoOffseed) {
                    isAsso = true;
                  }

                  // Using old way
                  else {
                    if (useUnassociatedL1) {
                      // convert the TSOS into a PTSOD
                      PTrajectoryStateOnDet const &seedTSOS =
                          trajectoryStateTransform::persistentState(newTSOS, newTSOSDet->geographicalId().rawId());
                      output->push_back(
                          L2MuonTrajectorySeed(seedTSOS,
                                               container,
                                               alongMomentum,
                                               MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
                    }
                  }
                }
              }
            }
          }  // useOfflineSeed

          else {
            // convert the TSOS into a PTSOD
            PTrajectoryStateOnDet const &seedTSOS = trajectoryStateTransform::persistentState(tsos, theid.rawId());
            output->push_back(L2MuonTrajectorySeed(
                seedTSOS, container, alongMomentum, MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
          }
        }
      }
    }

    // MatchType : 0 Old matching,   1 L1 Order(1to1), 2 Min dR(1to1), 3 Higher Q(1to1), 4 All matched L1
    if (useOfflineSeed && isAsso) {
      unsigned int nOfflineSeed1 = offlineSeedHandle->size();
      unsigned int nL1;
      unsigned int i, j;  // for the matrix element

      if (matchType == 1) {
        vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);

        for (nL1 = 0; nL1 < nMuColl; ++nL1) {
          double tempDR = 99.;
          unsigned int theOffs = 0;

          for (j = 0; j < nOfflineSeed1; ++j) {
            if (removed_col[j])
              continue;
            if (tempDR > dRmtx[nL1][j]) {
              tempDR = dRmtx[nL1][j];
              theOffs = j;
            }
          }

          auto it = muColl->begin(muColl->getFirstBX()) + nL1;

          double newDRcone = matchingDR[0];
          if (fabs(it->eta()) < etaBins.back()) {
            std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
            newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
          }

          if (!(tempDR < newDRcone))
            continue;

          // Remove column for given offSeed
          removed_col[theOffs] = true;

          if (selOffseeds[nL1][theOffs] != nullptr) {
            //put given L1mu and offSeed to output
            edm::OwnVector<TrackingRecHit> newContainer;

            PTrajectoryStateOnDet const &seedTSOS = selOffseeds[nL1][theOffs]->startingState();
            TrajectorySeed::const_iterator tsci = selOffseeds[nL1][theOffs]->recHits().first,
                                           tscie = selOffseeds[nL1][theOffs]->recHits().second;
            for (; tsci != tscie; ++tsci) {
              newContainer.push_back(*tsci);
            }
            output->push_back(L2MuonTrajectorySeed(seedTSOS,
                                                   newContainer,
                                                   alongMomentum,
                                                   MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
          }
        }
      }

      else if (matchType == 2) {
        vector<bool> removed_row = vector<bool>(nMuColl, false);
        vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);

        for (nL1 = 0; nL1 < nMuColl; ++nL1) {
          double tempDR = 99.;
          unsigned int theL1 = 0;
          unsigned int theOffs = 0;

          for (i = 0; i < nMuColl; ++i) {
            if (removed_row[i])
              continue;
            for (j = 0; j < nOfflineSeed1; ++j) {
              if (removed_col[j])
                continue;
              if (tempDR > dRmtx[i][j]) {
                tempDR = dRmtx[i][j];
                theL1 = i;
                theOffs = j;
              }
            }
          }

          auto it = muColl->begin(muColl->getFirstBX()) + theL1;

          double newDRcone = matchingDR[0];
          if (fabs(it->eta()) < etaBins.back()) {
            std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
            newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
          }

          if (!(tempDR < newDRcone))
            continue;

          // Remove row and column for given (L1mu, offSeed)
          removed_col[theOffs] = true;
          removed_row[theL1] = true;

          if (selOffseeds[theL1][theOffs] != nullptr) {
            //put given L1mu and offSeed to output
            edm::OwnVector<TrackingRecHit> newContainer;

            PTrajectoryStateOnDet const &seedTSOS = selOffseeds[theL1][theOffs]->startingState();
            TrajectorySeed::const_iterator tsci = selOffseeds[theL1][theOffs]->recHits().first,
                                           tscie = selOffseeds[theL1][theOffs]->recHits().second;
            for (; tsci != tscie; ++tsci) {
              newContainer.push_back(*tsci);
            }
            output->push_back(L2MuonTrajectorySeed(seedTSOS,
                                                   newContainer,
                                                   alongMomentum,
                                                   MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
          }
        }
      }

      else if (matchType == 3) {
        vector<bool> removed_row = vector<bool>(nMuColl, false);
        vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);

        for (nL1 = 0; nL1 < nMuColl; ++nL1) {
          double tempDR = 99.;
          unsigned int theL1 = 0;
          unsigned int theOffs = 0;
          auto theit = muColl->begin(muColl->getFirstBX());

          // L1Q > 10 (L1Q = 12)
          for (i = 0; i < nMuColl; ++i) {
            if (removed_row[i])
              continue;
            theit = muColl->begin(muColl->getFirstBX()) + i;
            if (theit->hwQual() > 10) {
              for (j = 0; j < nOfflineSeed1; ++j) {
                if (removed_col[j])
                  continue;
                if (tempDR > dRmtx[i][j]) {
                  tempDR = dRmtx[i][j];
                  theL1 = i;
                  theOffs = j;
                }
              }
            }
          }

          // 6 < L1Q <= 10 (L1Q = 8)
          if (tempDR == 99.) {
            for (i = 0; i < nMuColl; ++i) {
              if (removed_row[i])
                continue;
              theit = muColl->begin(muColl->getFirstBX()) + i;
              if ((theit->hwQual() <= 10) && (theit->hwQual() > 6)) {
                for (j = 0; j < nOfflineSeed1; ++j) {
                  if (removed_col[j])
                    continue;
                  if (tempDR > dRmtx[i][j]) {
                    tempDR = dRmtx[i][j];
                    theL1 = i;
                    theOffs = j;
                  }
                }
              }
            }
          }

          // L1Q <= 6 (L1Q = 4)
          if (tempDR == 99.) {
            for (i = 0; i < nMuColl; ++i) {
              if (removed_row[i])
                continue;
              theit = muColl->begin(muColl->getFirstBX()) + i;
              if (theit->hwQual() <= 6) {
                for (j = 0; j < nOfflineSeed1; ++j) {
                  if (removed_col[j])
                    continue;
                  if (tempDR > dRmtx[i][j]) {
                    tempDR = dRmtx[i][j];
                    theL1 = i;
                    theOffs = j;
                  }
                }
              }
            }
          }

          auto it = muColl->begin(muColl->getFirstBX()) + theL1;

          double newDRcone = matchingDR[0];
          if (fabs(it->eta()) < etaBins.back()) {
            std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
            newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
          }

          if (!(tempDR < newDRcone))
            continue;

          // Remove row and column for given (L1mu, offSeed)
          removed_col[theOffs] = true;
          removed_row[theL1] = true;

          if (selOffseeds[theL1][theOffs] != nullptr) {
            //put given L1mu and offSeed to output
            edm::OwnVector<TrackingRecHit> newContainer;

            PTrajectoryStateOnDet const &seedTSOS = selOffseeds[theL1][theOffs]->startingState();
            TrajectorySeed::const_iterator tsci = selOffseeds[theL1][theOffs]->recHits().first,
                                           tscie = selOffseeds[theL1][theOffs]->recHits().second;
            for (; tsci != tscie; ++tsci) {
              newContainer.push_back(*tsci);
            }
            output->push_back(L2MuonTrajectorySeed(seedTSOS,
                                                   newContainer,
                                                   alongMomentum,
                                                   MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
          }
        }
      }

      else if (matchType == 4) {
        for (i = 0; i < nMuColl; ++i) {
          auto it = muColl->begin(muColl->getFirstBX()) + i;

          double tempDR = 99.;
          unsigned int theOffs = 0;

          double newDRcone = matchingDR[0];
          if (fabs(it->eta()) < etaBins.back()) {
            std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
            newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
          }

          for (j = 0; j < nOfflineSeed1; ++j) {
            if (tempDR > dRmtx[i][j]) {
              tempDR = dRmtx[i][j];
              theOffs = j;
            }
          }

          if (!(tempDR < newDRcone))
            continue;

          if (selOffseeds[i][theOffs] != nullptr) {
            //put given L1mu and offSeed to output
            edm::OwnVector<TrackingRecHit> newContainer;

            PTrajectoryStateOnDet const &seedTSOS = selOffseeds[i][theOffs]->startingState();
            TrajectorySeed::const_iterator tsci = selOffseeds[i][theOffs]->recHits().first,
                                           tscie = selOffseeds[i][theOffs]->recHits().second;
            for (; tsci != tscie; ++tsci) {
              newContainer.push_back(*tsci);
            }
            output->push_back(L2MuonTrajectorySeed(seedTSOS,
                                                   newContainer,
                                                   alongMomentum,
                                                   MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
          }
        }
      }

    }  // useOfflineSeed && isAsso

  }  // matchType > 0

  //--- SortType 1 : by L1 pT
  if (sortType == 1) {
    sort(output->begin(), output->end(), sortByL1Pt);
  }

  //--- SortType 2 : by L1 quality and pT
  else if (sortType == 2) {
    sort(output->begin(), output->end(), sortByL1QandPt);
  }

  iEvent.put(std::move(output));
}

// FIXME: does not resolve ambiguities yet!
const TrajectorySeed *L2MuonSeedGeneratorFromL1T::associateOfflineSeedToL1(
    edm::Handle<edm::View<TrajectorySeed>> &offseeds,
    std::vector<int> &offseedMap,
    TrajectoryStateOnSurface &newTsos,
    double dRcone) {
  const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
  MuonPatternRecoDumper debugtmp;

  edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
  const TrajectorySeed *selOffseed = nullptr;
  double bestDr = 99999.;
  unsigned int nOffseed(0);
  int lastOffseed(-1);

  for (offseed = offseeds->begin(); offseed != endOffseed; ++offseed, ++nOffseed) {
    if (offseedMap[nOffseed] != 0)
      continue;
    GlobalPoint glbPos = theService->trackingGeometry()
                             ->idToDet(offseed->startingState().detId())
                             ->surface()
                             .toGlobal(offseed->startingState().parameters().position());
    GlobalVector glbMom = theService->trackingGeometry()
                              ->idToDet(offseed->startingState().detId())
                              ->surface()
                              .toGlobal(offseed->startingState().parameters().momentum());

    // Preliminary check
    double preDr = deltaR(newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi());
    if (preDr > 1.0)
      continue;

    const FreeTrajectoryState offseedFTS(
        glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
    TrajectoryStateOnSurface offseedTsos =
        theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
    LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
    LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
    LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi()
                       << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
    LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge() * offseedFTS.momentum().perp()
                       << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")"
                       << std::endl
                       << std::endl;

    if (offseedTsos.isValid()) {
      LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
      LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag()
                         << ", phi=" << offseedTsos.globalPosition().phi()
                         << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
      LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge() * offseedTsos.globalMomentum().perp()
                         << ", phi=" << offseedTsos.globalMomentum().phi()
                         << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl
                         << std::endl;
      double newDr = deltaR(newTsos.globalPosition().eta(),
                            newTsos.globalPosition().phi(),
                            offseedTsos.globalPosition().eta(),
                            offseedTsos.globalPosition().phi());
      LogDebug(metlabel) << "   -- DR = " << newDr << std::endl;
      if (newDr < dRcone && newDr < bestDr) {
        LogDebug(metlabel) << "          --> OK! " << newDr << std::endl << std::endl;
        selOffseed = &*offseed;
        bestDr = newDr;
        offseedMap[nOffseed] = 1;
        if (lastOffseed > -1)
          offseedMap[lastOffseed] = 0;
        lastOffseed = nOffseed;
      } else {
        LogDebug(metlabel) << "          --> Rejected. " << newDr << std::endl << std::endl;
      }
    } else {
      LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
    }
  }

  return selOffseed;
}

bool L2MuonSeedGeneratorFromL1T::isAssociateOfflineSeedToL1(
    edm::Handle<edm::View<TrajectorySeed>> &offseeds,
    std::vector<std::vector<double>> &dRmtx,
    TrajectoryStateOnSurface &newTsos,
    unsigned int imu,
    std::vector<std::vector<const TrajectorySeed *>> &selOffseeds,
    double dRcone) {
  const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
  bool isAssociated = false;
  MuonPatternRecoDumper debugtmp;

  edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
  unsigned int nOffseed(0);

  for (offseed = offseeds->begin(); offseed != endOffseed; ++offseed, ++nOffseed) {
    GlobalPoint glbPos = theService->trackingGeometry()
                             ->idToDet(offseed->startingState().detId())
                             ->surface()
                             .toGlobal(offseed->startingState().parameters().position());
    GlobalVector glbMom = theService->trackingGeometry()
                              ->idToDet(offseed->startingState().detId())
                              ->surface()
                              .toGlobal(offseed->startingState().parameters().momentum());

    // Preliminary check
    double preDr = deltaR(newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi());
    if (preDr > 1.0)
      continue;

    const FreeTrajectoryState offseedFTS(
        glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
    TrajectoryStateOnSurface offseedTsos =
        theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
    LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
    LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
    LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi()
                       << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
    LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge() * offseedFTS.momentum().perp()
                       << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")"
                       << std::endl
                       << std::endl;

    if (offseedTsos.isValid()) {
      LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
      LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag()
                         << ", phi=" << offseedTsos.globalPosition().phi()
                         << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
      LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge() * offseedTsos.globalMomentum().perp()
                         << ", phi=" << offseedTsos.globalMomentum().phi()
                         << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl
                         << std::endl;
      double newDr = deltaR(newTsos.globalPosition().eta(),
                            newTsos.globalPosition().phi(),
                            offseedTsos.globalPosition().eta(),
                            offseedTsos.globalPosition().phi());

      LogDebug(metlabel) << "   -- DR = " << newDr << std::endl;
      if (newDr < dRcone) {
        LogDebug(metlabel) << "          --> OK! " << newDr << std::endl << std::endl;

        dRmtx[imu][nOffseed] = newDr;
        selOffseeds[imu][nOffseed] = &*offseed;

        isAssociated = true;
      } else {
        LogDebug(metlabel) << "          --> Rejected. " << newDr << std::endl << std::endl;
      }
    } else {
      LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
    }
  }

  return isAssociated;
}