Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:25:50

0001 //-------------------------------------------------
0002 //
0003 /**  \class L2MuonSeedGeneratorFromL1T
0004  * 
0005  *   L2 muon seed generator:
0006  *   Transform the L1 informations in seeds for the
0007  *   L2 muon reconstruction
0008  *
0009  *
0010  *
0011  *   \author  A.Everett, R.Bellan, J. Alcaraz
0012  *
0013  *    ORCA's author: N. Neumeister 
0014  *
0015  *    Modified by M.Oh
0016  */
0017 //
0018 //--------------------------------------------------
0019 
0020 // Class Header
0021 #include "RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.h"
0022 
0023 // Framework
0024 #include "FWCore/Framework/interface/ConsumesCollector.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0030 
0031 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0032 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
0033 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0034 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0035 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0036 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0037 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0038 
0039 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0040 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0041 
0042 using namespace std;
0043 using namespace edm;
0044 using namespace l1t;
0045 
0046 //--- Functions used in output sorting
0047 bool sortByL1Pt(L2MuonTrajectorySeed &A, L2MuonTrajectorySeed &B) {
0048   l1t::MuonRef Ref_L1A = A.l1tParticle();
0049   l1t::MuonRef Ref_L1B = B.l1tParticle();
0050   return (Ref_L1A->pt() > Ref_L1B->pt());
0051 };
0052 
0053 bool sortByL1QandPt(L2MuonTrajectorySeed &A, L2MuonTrajectorySeed &B) {
0054   l1t::MuonRef Ref_L1A = A.l1tParticle();
0055   l1t::MuonRef Ref_L1B = B.l1tParticle();
0056 
0057   // Compare quality first
0058   if (Ref_L1A->hwQual() > Ref_L1B->hwQual())
0059     return true;
0060   if (Ref_L1A->hwQual() < Ref_L1B->hwQual())
0061     return false;
0062 
0063   // For same quality L1s compare pT
0064   return (Ref_L1A->pt() > Ref_L1B->pt());
0065 };
0066 
0067 // constructors
0068 L2MuonSeedGeneratorFromL1T::L2MuonSeedGeneratorFromL1T(const edm::ParameterSet &iConfig)
0069     : theSource(iConfig.getParameter<InputTag>("InputObjects")),
0070       theL1GMTReadoutCollection(iConfig.getParameter<InputTag>("GMTReadoutCollection")),  // to be removed
0071       thePropagatorName(iConfig.getParameter<string>("Propagator")),
0072       theL1MinPt(iConfig.getParameter<double>("L1MinPt")),
0073       theL1MaxEta(iConfig.getParameter<double>("L1MaxEta")),
0074       theL1MinQuality(iConfig.getParameter<unsigned int>("L1MinQuality")),
0075       theMinPtBarrel(iConfig.getParameter<double>("SetMinPtBarrelTo")),
0076       theMinPtEndcap(iConfig.getParameter<double>("SetMinPtEndcapTo")),
0077       useOfflineSeed(iConfig.getUntrackedParameter<bool>("UseOfflineSeed", false)),
0078       useUnassociatedL1(iConfig.getParameter<bool>("UseUnassociatedL1")),
0079       matchingDR(iConfig.getParameter<std::vector<double>>("MatchDR")),
0080       etaBins(iConfig.getParameter<std::vector<double>>("EtaMatchingBins")),
0081       centralBxOnly_(iConfig.getParameter<bool>("CentralBxOnly")),
0082       matchType(iConfig.getParameter<unsigned int>("MatchType")),
0083       sortType(iConfig.getParameter<unsigned int>("SortType")) {
0084   muCollToken_ = consumes<MuonBxCollection>(theSource);
0085 
0086   if (useOfflineSeed) {
0087     theOfflineSeedLabel = iConfig.getUntrackedParameter<InputTag>("OfflineSeedLabel");
0088     offlineSeedToken_ = consumes<edm::View<TrajectorySeed>>(theOfflineSeedLabel);
0089 
0090     // check that number of eta bins -1 matches number of dR cones
0091     if (matchingDR.size() != etaBins.size() - 1) {
0092       throw cms::Exception("Configuration") << "Size of MatchDR "
0093                                             << "does not match number of eta bins." << endl;
0094     }
0095   }
0096 
0097   if (matchType > 4 || sortType > 3) {
0098     throw cms::Exception("Configuration") << "Wrong MatchType or SortType" << endl;
0099   }
0100 
0101   // service parameters
0102   ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
0103 
0104   // the services
0105   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0106 
0107   // the estimator
0108   theEstimator = new Chi2MeasurementEstimator(10000.);
0109 
0110   produces<L2MuonTrajectorySeedCollection>();
0111 }
0112 
0113 // destructor
0114 L2MuonSeedGeneratorFromL1T::~L2MuonSeedGeneratorFromL1T() {
0115   if (theService)
0116     delete theService;
0117   if (theEstimator)
0118     delete theEstimator;
0119 }
0120 
0121 void L2MuonSeedGeneratorFromL1T::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0122   edm::ParameterSetDescription desc;
0123   desc.add<edm::InputTag>("GMTReadoutCollection", edm::InputTag(""));  // to be removed
0124   desc.add<edm::InputTag>("InputObjects", edm::InputTag("hltGmtStage2Digis"));
0125   desc.add<string>("Propagator", "");
0126   desc.add<double>("L1MinPt", -1.);
0127   desc.add<double>("L1MaxEta", 5.0);
0128   desc.add<unsigned int>("L1MinQuality", 0);
0129   desc.add<double>("SetMinPtBarrelTo", 3.5);
0130   desc.add<double>("SetMinPtEndcapTo", 1.0);
0131   desc.addUntracked<bool>("UseOfflineSeed", false);
0132   desc.add<bool>("UseUnassociatedL1", true);
0133   desc.add<std::vector<double>>("MatchDR", {0.3});
0134   desc.add<std::vector<double>>("EtaMatchingBins", {0., 2.5});
0135   desc.add<bool>("CentralBxOnly", true);
0136   desc.add<unsigned int>("MatchType", 0)
0137       ->setComment(
0138           "MatchType : 0 Old matching,   1 L1 Order(1to1), 2 Min dR(1to1), 3 Higher Q(1to1), 4 All matched L1");
0139   desc.add<unsigned int>("SortType", 0)->setComment("SortType :  0 not sort,       1 Pt, 2 Q and Pt");
0140   desc.addUntracked<edm::InputTag>("OfflineSeedLabel", edm::InputTag(""));
0141 
0142   edm::ParameterSetDescription psd0;
0143   psd0.addUntracked<std::vector<std::string>>("Propagators", {"SteppingHelixPropagatorAny"});
0144   psd0.add<bool>("RPCLayers", true);
0145   psd0.addUntracked<bool>("UseMuonNavigation", true);
0146   desc.add<edm::ParameterSetDescription>("ServiceParameters", psd0);
0147   descriptions.add("L2MuonSeedGeneratorFromL1T", desc);
0148 }
0149 
0150 void L2MuonSeedGeneratorFromL1T::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0151   const std::string metname = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
0152   MuonPatternRecoDumper debug;
0153 
0154   auto output = std::make_unique<L2MuonTrajectorySeedCollection>();
0155 
0156   edm::Handle<MuonBxCollection> muColl;
0157   iEvent.getByToken(muCollToken_, muColl);
0158   LogDebug(metname) << "Number of muons " << muColl->size() << endl;
0159 
0160   //--- matchType 0 : Old logic
0161   if (matchType == 0) {
0162     edm::Handle<edm::View<TrajectorySeed>> offlineSeedHandle;
0163     vector<int> offlineSeedMap;
0164     if (useOfflineSeed) {
0165       iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
0166       LogDebug(metname) << "Number of offline seeds " << offlineSeedHandle->size() << endl;
0167       offlineSeedMap = vector<int>(offlineSeedHandle->size(), 0);
0168     }
0169 
0170     for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) {
0171       if (centralBxOnly_ && (ibx != 0))
0172         continue;
0173       for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++) {
0174         unsigned int quality = it->hwQual();
0175         int valid_charge = it->hwChargeValid();
0176 
0177         float pt = it->pt();
0178         float eta = it->eta();
0179         float theta = 2 * atan(exp(-eta));
0180         float phi = it->phi();
0181         int charge = it->charge();
0182         // Set charge=0 for the time being if the valid charge bit is zero
0183         if (!valid_charge)
0184           charge = 0;
0185 
0186         // link number indices of the optical fibres that connect the uGMT with the track finders
0187         //EMTF+ : 36-41, OMTF+ : 42-47, BMTF : 48-59, OMTF- : 60-65, EMTF- : 66-71
0188         int link = 36 + (int)(it->tfMuonIndex() / 3.);
0189         bool barrel = true;
0190         if ((link >= 36 && link <= 41) || (link >= 66 && link <= 71))
0191           barrel = false;
0192 
0193         if (pt < theL1MinPt || fabs(eta) > theL1MaxEta)
0194           continue;
0195 
0196         LogDebug(metname) << "New L2 Muon Seed";
0197         LogDebug(metname) << "Pt = " << pt << " GeV/c";
0198         LogDebug(metname) << "eta = " << eta;
0199         LogDebug(metname) << "theta = " << theta << " rad";
0200         LogDebug(metname) << "phi = " << phi << " rad";
0201         LogDebug(metname) << "charge = " << charge;
0202         LogDebug(metname) << "In Barrel? = " << barrel;
0203 
0204         if (quality <= theL1MinQuality)
0205           continue;
0206         LogDebug(metname) << "quality = " << quality;
0207 
0208         // Update the services
0209         theService->update(iSetup);
0210 
0211         const DetLayer *detLayer = nullptr;
0212         float radius = 0.;
0213 
0214         CLHEP::Hep3Vector vec(0., 1., 0.);
0215         vec.setTheta(theta);
0216         vec.setPhi(phi);
0217 
0218         DetId theid;
0219         // Get the det layer on which the state should be put
0220         if (barrel) {
0221           LogDebug(metname) << "The seed is in the barrel";
0222 
0223           // MB2
0224           theid = DTChamberId(0, 2, 0);
0225           detLayer = theService->detLayerGeometry()->idToLayer(theid);
0226           LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0227 
0228           const BoundSurface *sur = &(detLayer->surface());
0229           const BoundCylinder *bc = dynamic_cast<const BoundCylinder *>(sur);
0230 
0231           radius = fabs(bc->radius() / sin(theta));
0232 
0233           LogDebug(metname) << "radius " << radius;
0234 
0235           if (pt < theMinPtBarrel)
0236             pt = theMinPtBarrel;
0237         } else {
0238           LogDebug(metname) << "The seed is in the endcap";
0239 
0240           // ME2
0241           theid = theta < Geom::pi() / 2. ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
0242 
0243           detLayer = theService->detLayerGeometry()->idToLayer(theid);
0244           LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0245 
0246           radius = fabs(detLayer->position().z() / cos(theta));
0247 
0248           if (pt < theMinPtEndcap)
0249             pt = theMinPtEndcap;
0250         }
0251 
0252         vec.setMag(radius);
0253 
0254         GlobalPoint pos(vec.x(), vec.y(), vec.z());
0255 
0256         GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));
0257 
0258         GlobalTrajectoryParameters param(pos, mom, charge, &*theService->magneticField());
0259         AlgebraicSymMatrix55 mat;
0260 
0261         mat[0][0] = (0.25 / pt) * (0.25 / pt);  // sigma^2(charge/abs_momentum)
0262         if (!barrel)
0263           mat[0][0] = (0.4 / pt) * (0.4 / pt);
0264 
0265         //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
0266         if (!valid_charge)
0267           mat[0][0] = (1. / pt) * (1. / pt);
0268 
0269         mat[1][1] = 0.05 * 0.05;  // sigma^2(lambda)
0270         mat[2][2] = 0.2 * 0.2;    // sigma^2(phi)
0271         mat[3][3] = 20. * 20.;    // sigma^2(x_transverse))
0272         mat[4][4] = 20. * 20.;    // sigma^2(y_transverse))
0273 
0274         CurvilinearTrajectoryError error(mat);
0275 
0276         const FreeTrajectoryState state(param, error);
0277 
0278         LogDebug(metname) << "Free trajectory State from the parameters";
0279         LogDebug(metname) << debug.dumpFTS(state);
0280 
0281         // Propagate the state on the MB2/ME2 surface
0282         TrajectoryStateOnSurface tsos =
0283             theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());
0284 
0285         LogDebug(metname) << "State after the propagation on the layer";
0286         LogDebug(metname) << debug.dumpLayer(detLayer);
0287         LogDebug(metname) << debug.dumpTSOS(tsos);
0288 
0289         double dRcone = matchingDR[0];
0290         if (fabs(eta) < etaBins.back()) {
0291           std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(eta));
0292           dRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0293         }
0294 
0295         if (tsos.isValid()) {
0296           edm::OwnVector<TrackingRecHit> container;
0297 
0298           if (useOfflineSeed && (!valid_charge || charge == 0)) {
0299             const TrajectorySeed *assoOffseed =
0300                 associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, tsos, dRcone);
0301 
0302             if (assoOffseed != nullptr) {
0303               PTrajectoryStateOnDet const &seedTSOS = assoOffseed->startingState();
0304               for (auto const &recHit : assoOffseed->recHits()) {
0305                 container.push_back(recHit);
0306               }
0307               output->push_back(
0308                   L2MuonTrajectorySeed(seedTSOS,
0309                                        container,
0310                                        alongMomentum,
0311                                        MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0312             } else {
0313               if (useUnassociatedL1) {
0314                 // convert the TSOS into a PTSOD
0315                 PTrajectoryStateOnDet const &seedTSOS = trajectoryStateTransform::persistentState(tsos, theid.rawId());
0316                 output->push_back(
0317                     L2MuonTrajectorySeed(seedTSOS,
0318                                          container,
0319                                          alongMomentum,
0320                                          MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0321               }
0322             }
0323           } else if (useOfflineSeed && valid_charge) {
0324             // Get the compatible dets on the layer
0325             std::vector<pair<const GeomDet *, TrajectoryStateOnSurface>> detsWithStates =
0326                 detLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
0327 
0328             if (detsWithStates.empty() && barrel) {
0329               // Fallback solution using ME2, try again to propagate but using ME2 as reference
0330               DetId fallback_id;
0331               theta < Geom::pi() / 2. ? fallback_id = CSCDetId(1, 2, 0, 0, 0) : fallback_id = CSCDetId(2, 2, 0, 0, 0);
0332               const DetLayer *ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id);
0333 
0334               tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface());
0335               detsWithStates =
0336                   ME2DetLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
0337             }
0338 
0339             if (!detsWithStates.empty()) {
0340               TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
0341               const GeomDet *newTSOSDet = detsWithStates.front().first;
0342 
0343               LogDebug(metname) << "Most compatible det";
0344               LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
0345 
0346               if (newTSOS.isValid()) {
0347                 LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
0348                                   << ", phi=" << newTSOS.globalPosition().phi()
0349                                   << ", eta=" << newTSOS.globalPosition().eta() << ")";
0350                 LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
0351                                   << ", phi=" << newTSOS.globalMomentum().phi()
0352                                   << ", eta=" << newTSOS.globalMomentum().eta() << ")";
0353 
0354                 const TrajectorySeed *assoOffseed =
0355                     associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS, dRcone);
0356 
0357                 if (assoOffseed != nullptr) {
0358                   PTrajectoryStateOnDet const &seedTSOS = assoOffseed->startingState();
0359                   for (auto const &recHit : assoOffseed->recHits()) {
0360                     container.push_back(recHit);
0361                   }
0362                   output->push_back(
0363                       L2MuonTrajectorySeed(seedTSOS,
0364                                            container,
0365                                            alongMomentum,
0366                                            MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0367                 } else {
0368                   if (useUnassociatedL1) {
0369                     // convert the TSOS into a PTSOD
0370                     PTrajectoryStateOnDet const &seedTSOS =
0371                         trajectoryStateTransform::persistentState(newTSOS, newTSOSDet->geographicalId().rawId());
0372                     output->push_back(
0373                         L2MuonTrajectorySeed(seedTSOS,
0374                                              container,
0375                                              alongMomentum,
0376                                              MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0377                   }
0378                 }
0379               }
0380             }
0381           } else {
0382             // convert the TSOS into a PTSOD
0383             PTrajectoryStateOnDet const &seedTSOS = trajectoryStateTransform::persistentState(tsos, theid.rawId());
0384             output->push_back(L2MuonTrajectorySeed(
0385                 seedTSOS, container, alongMomentum, MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0386           }
0387         }
0388       }
0389     }
0390 
0391   }  // End of matchType 0
0392 
0393   //--- matchType > 0 : Updated logics
0394   else if (matchType > 0) {
0395     unsigned int nMuColl = muColl->size();
0396 
0397     vector<vector<double>> dRmtx;
0398     vector<vector<const TrajectorySeed *>> selOffseeds;
0399 
0400     edm::Handle<edm::View<TrajectorySeed>> offlineSeedHandle;
0401     if (useOfflineSeed) {
0402       iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
0403       unsigned int nOfflineSeed = offlineSeedHandle->size();
0404       LogDebug(metname) << "Number of offline seeds " << nOfflineSeed << endl;
0405 
0406       // Initialize dRmtx and selOffseeds
0407       dRmtx = vector<vector<double>>(nMuColl, vector<double>(nOfflineSeed, 999.0));
0408       selOffseeds =
0409           vector<vector<const TrajectorySeed *>>(nMuColl, vector<const TrajectorySeed *>(nOfflineSeed, nullptr));
0410     }
0411 
0412     // At least one L1 muons are associated with offSeed
0413     bool isAsso = false;
0414 
0415     //--- Fill dR matrix
0416     for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) {
0417       if (centralBxOnly_ && (ibx != 0))
0418         continue;
0419       for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++) {
0420         //Position of given L1mu
0421         unsigned int imu = distance(muColl->begin(muColl->getFirstBX()), it);
0422 
0423         unsigned int quality = it->hwQual();
0424         int valid_charge = it->hwChargeValid();
0425 
0426         float pt = it->pt();
0427         float eta = it->eta();
0428         float theta = 2 * atan(exp(-eta));
0429         float phi = it->phi();
0430         int charge = it->charge();
0431         // Set charge=0 for the time being if the valid charge bit is zero
0432         if (!valid_charge)
0433           charge = 0;
0434 
0435         // link number indices of the optical fibres that connect the uGMT with the track finders
0436         //EMTF+ : 36-41, OMTF+ : 42-47, BMTF : 48-59, OMTF- : 60-65, EMTF- : 66-71
0437         int link = 36 + (int)(it->tfMuonIndex() / 3.);
0438         bool barrel = true;
0439         if ((link >= 36 && link <= 41) || (link >= 66 && link <= 71))
0440           barrel = false;
0441 
0442         if (pt < theL1MinPt || fabs(eta) > theL1MaxEta)
0443           continue;
0444 
0445         LogDebug(metname) << "New L2 Muon Seed";
0446         LogDebug(metname) << "Pt = " << pt << " GeV/c";
0447         LogDebug(metname) << "eta = " << eta;
0448         LogDebug(metname) << "theta = " << theta << " rad";
0449         LogDebug(metname) << "phi = " << phi << " rad";
0450         LogDebug(metname) << "charge = " << charge;
0451         LogDebug(metname) << "In Barrel? = " << barrel;
0452 
0453         if (quality <= theL1MinQuality)
0454           continue;
0455         LogDebug(metname) << "quality = " << quality;
0456 
0457         // Update the services
0458         theService->update(iSetup);
0459 
0460         const DetLayer *detLayer = nullptr;
0461         float radius = 0.;
0462 
0463         CLHEP::Hep3Vector vec(0., 1., 0.);
0464         vec.setTheta(theta);
0465         vec.setPhi(phi);
0466 
0467         DetId theid;
0468         // Get the det layer on which the state should be put
0469         if (barrel) {
0470           LogDebug(metname) << "The seed is in the barrel";
0471 
0472           // MB2
0473           theid = DTChamberId(0, 2, 0);
0474           detLayer = theService->detLayerGeometry()->idToLayer(theid);
0475           LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0476 
0477           const BoundSurface *sur = &(detLayer->surface());
0478           const BoundCylinder *bc = dynamic_cast<const BoundCylinder *>(sur);
0479 
0480           radius = fabs(bc->radius() / sin(theta));
0481 
0482           LogDebug(metname) << "radius " << radius;
0483 
0484           if (pt < theMinPtBarrel)
0485             pt = theMinPtBarrel;
0486         } else {
0487           LogDebug(metname) << "The seed is in the endcap";
0488 
0489           // ME2
0490           theid = theta < Geom::pi() / 2. ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
0491 
0492           detLayer = theService->detLayerGeometry()->idToLayer(theid);
0493           LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0494 
0495           radius = fabs(detLayer->position().z() / cos(theta));
0496 
0497           if (pt < theMinPtEndcap)
0498             pt = theMinPtEndcap;
0499         }
0500 
0501         vec.setMag(radius);
0502 
0503         GlobalPoint pos(vec.x(), vec.y(), vec.z());
0504 
0505         GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));
0506 
0507         GlobalTrajectoryParameters param(pos, mom, charge, &*theService->magneticField());
0508         AlgebraicSymMatrix55 mat;
0509 
0510         mat[0][0] = (0.25 / pt) * (0.25 / pt);  // sigma^2(charge/abs_momentum)
0511         if (!barrel)
0512           mat[0][0] = (0.4 / pt) * (0.4 / pt);
0513 
0514         //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
0515         if (!valid_charge)
0516           mat[0][0] = (1. / pt) * (1. / pt);
0517 
0518         mat[1][1] = 0.05 * 0.05;  // sigma^2(lambda)
0519         mat[2][2] = 0.2 * 0.2;    // sigma^2(phi)
0520         mat[3][3] = 20. * 20.;    // sigma^2(x_transverse))
0521         mat[4][4] = 20. * 20.;    // sigma^2(y_transverse))
0522 
0523         CurvilinearTrajectoryError error(mat);
0524 
0525         const FreeTrajectoryState state(param, error);
0526 
0527         LogDebug(metname) << "Free trajectory State from the parameters";
0528         LogDebug(metname) << debug.dumpFTS(state);
0529 
0530         // Propagate the state on the MB2/ME2 surface
0531         TrajectoryStateOnSurface tsos =
0532             theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());
0533 
0534         LogDebug(metname) << "State after the propagation on the layer";
0535         LogDebug(metname) << debug.dumpLayer(detLayer);
0536         LogDebug(metname) << debug.dumpTSOS(tsos);
0537 
0538         double dRcone = matchingDR[0];
0539         if (fabs(eta) < etaBins.back()) {
0540           std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(eta));
0541           dRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0542         }
0543 
0544         if (tsos.isValid()) {
0545           edm::OwnVector<TrackingRecHit> container;
0546 
0547           if (useOfflineSeed) {
0548             if ((!valid_charge || charge == 0)) {
0549               bool isAssoOffseed = isAssociateOfflineSeedToL1(offlineSeedHandle, dRmtx, tsos, imu, selOffseeds, dRcone);
0550 
0551               if (isAssoOffseed) {
0552                 isAsso = true;
0553               }
0554 
0555               // Using old way
0556               else {
0557                 if (useUnassociatedL1) {
0558                   // convert the TSOS into a PTSOD
0559                   PTrajectoryStateOnDet const &seedTSOS =
0560                       trajectoryStateTransform::persistentState(tsos, theid.rawId());
0561                   output->push_back(
0562                       L2MuonTrajectorySeed(seedTSOS,
0563                                            container,
0564                                            alongMomentum,
0565                                            MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0566                 }
0567               }
0568 
0569             }
0570 
0571             else if (valid_charge) {
0572               // Get the compatible dets on the layer
0573               std::vector<pair<const GeomDet *, TrajectoryStateOnSurface>> detsWithStates =
0574                   detLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
0575 
0576               if (detsWithStates.empty() && barrel) {
0577                 // Fallback solution using ME2, try again to propagate but using ME2 as reference
0578                 DetId fallback_id;
0579                 theta < Geom::pi() / 2. ? fallback_id = CSCDetId(1, 2, 0, 0, 0) : fallback_id = CSCDetId(2, 2, 0, 0, 0);
0580                 const DetLayer *ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id);
0581 
0582                 tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface());
0583                 detsWithStates =
0584                     ME2DetLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
0585               }
0586 
0587               if (!detsWithStates.empty()) {
0588                 TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
0589                 const GeomDet *newTSOSDet = detsWithStates.front().first;
0590 
0591                 LogDebug(metname) << "Most compatible det";
0592                 LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
0593 
0594                 if (newTSOS.isValid()) {
0595                   LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
0596                                     << ", phi=" << newTSOS.globalPosition().phi()
0597                                     << ", eta=" << newTSOS.globalPosition().eta() << ")";
0598                   LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
0599                                     << ", phi=" << newTSOS.globalMomentum().phi()
0600                                     << ", eta=" << newTSOS.globalMomentum().eta() << ")";
0601 
0602                   bool isAssoOffseed =
0603                       isAssociateOfflineSeedToL1(offlineSeedHandle, dRmtx, newTSOS, imu, selOffseeds, dRcone);
0604 
0605                   if (isAssoOffseed) {
0606                     isAsso = true;
0607                   }
0608 
0609                   // Using old way
0610                   else {
0611                     if (useUnassociatedL1) {
0612                       // convert the TSOS into a PTSOD
0613                       PTrajectoryStateOnDet const &seedTSOS =
0614                           trajectoryStateTransform::persistentState(newTSOS, newTSOSDet->geographicalId().rawId());
0615                       output->push_back(
0616                           L2MuonTrajectorySeed(seedTSOS,
0617                                                container,
0618                                                alongMomentum,
0619                                                MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0620                     }
0621                   }
0622                 }
0623               }
0624             }
0625           }  // useOfflineSeed
0626 
0627           else {
0628             // convert the TSOS into a PTSOD
0629             PTrajectoryStateOnDet const &seedTSOS = trajectoryStateTransform::persistentState(tsos, theid.rawId());
0630             output->push_back(L2MuonTrajectorySeed(
0631                 seedTSOS, container, alongMomentum, MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0632           }
0633         }
0634       }
0635     }
0636 
0637     // MatchType : 0 Old matching,   1 L1 Order(1to1), 2 Min dR(1to1), 3 Higher Q(1to1), 4 All matched L1
0638     if (useOfflineSeed && isAsso) {
0639       unsigned int nOfflineSeed1 = offlineSeedHandle->size();
0640       unsigned int nL1;
0641       unsigned int i, j;  // for the matrix element
0642 
0643       if (matchType == 1) {
0644         vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);
0645 
0646         for (nL1 = 0; nL1 < nMuColl; ++nL1) {
0647           double tempDR = 99.;
0648           unsigned int theOffs = 0;
0649 
0650           for (j = 0; j < nOfflineSeed1; ++j) {
0651             if (removed_col[j])
0652               continue;
0653             if (tempDR > dRmtx[nL1][j]) {
0654               tempDR = dRmtx[nL1][j];
0655               theOffs = j;
0656             }
0657           }
0658 
0659           auto it = muColl->begin(muColl->getFirstBX()) + nL1;
0660 
0661           double newDRcone = matchingDR[0];
0662           if (fabs(it->eta()) < etaBins.back()) {
0663             std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
0664             newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0665           }
0666 
0667           if (!(tempDR < newDRcone))
0668             continue;
0669 
0670           // Remove column for given offSeed
0671           removed_col[theOffs] = true;
0672 
0673           if (selOffseeds[nL1][theOffs] != nullptr) {
0674             //put given L1mu and offSeed to output
0675             edm::OwnVector<TrackingRecHit> newContainer;
0676 
0677             PTrajectoryStateOnDet const &seedTSOS = selOffseeds[nL1][theOffs]->startingState();
0678             for (auto const &recHit : selOffseeds[nL1][theOffs]->recHits()) {
0679               newContainer.push_back(recHit);
0680             }
0681             output->push_back(L2MuonTrajectorySeed(seedTSOS,
0682                                                    newContainer,
0683                                                    alongMomentum,
0684                                                    MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0685           }
0686         }
0687       }
0688 
0689       else if (matchType == 2) {
0690         vector<bool> removed_row = vector<bool>(nMuColl, false);
0691         vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);
0692 
0693         for (nL1 = 0; nL1 < nMuColl; ++nL1) {
0694           double tempDR = 99.;
0695           unsigned int theL1 = 0;
0696           unsigned int theOffs = 0;
0697 
0698           for (i = 0; i < nMuColl; ++i) {
0699             if (removed_row[i])
0700               continue;
0701             for (j = 0; j < nOfflineSeed1; ++j) {
0702               if (removed_col[j])
0703                 continue;
0704               if (tempDR > dRmtx[i][j]) {
0705                 tempDR = dRmtx[i][j];
0706                 theL1 = i;
0707                 theOffs = j;
0708               }
0709             }
0710           }
0711 
0712           auto it = muColl->begin(muColl->getFirstBX()) + theL1;
0713 
0714           double newDRcone = matchingDR[0];
0715           if (fabs(it->eta()) < etaBins.back()) {
0716             std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
0717             newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0718           }
0719 
0720           if (!(tempDR < newDRcone))
0721             continue;
0722 
0723           // Remove row and column for given (L1mu, offSeed)
0724           removed_col[theOffs] = true;
0725           removed_row[theL1] = true;
0726 
0727           if (selOffseeds[theL1][theOffs] != nullptr) {
0728             //put given L1mu and offSeed to output
0729             edm::OwnVector<TrackingRecHit> newContainer;
0730 
0731             PTrajectoryStateOnDet const &seedTSOS = selOffseeds[theL1][theOffs]->startingState();
0732             for (auto const &recHit : selOffseeds[theL1][theOffs]->recHits()) {
0733               newContainer.push_back(recHit);
0734             }
0735             output->push_back(L2MuonTrajectorySeed(seedTSOS,
0736                                                    newContainer,
0737                                                    alongMomentum,
0738                                                    MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0739           }
0740         }
0741       }
0742 
0743       else if (matchType == 3) {
0744         vector<bool> removed_row = vector<bool>(nMuColl, false);
0745         vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);
0746 
0747         for (nL1 = 0; nL1 < nMuColl; ++nL1) {
0748           double tempDR = 99.;
0749           unsigned int theL1 = 0;
0750           unsigned int theOffs = 0;
0751           auto theit = muColl->begin(muColl->getFirstBX());
0752 
0753           // L1Q > 10 (L1Q = 12)
0754           for (i = 0; i < nMuColl; ++i) {
0755             if (removed_row[i])
0756               continue;
0757             theit = muColl->begin(muColl->getFirstBX()) + i;
0758             if (theit->hwQual() > 10) {
0759               for (j = 0; j < nOfflineSeed1; ++j) {
0760                 if (removed_col[j])
0761                   continue;
0762                 if (tempDR > dRmtx[i][j]) {
0763                   tempDR = dRmtx[i][j];
0764                   theL1 = i;
0765                   theOffs = j;
0766                 }
0767               }
0768             }
0769           }
0770 
0771           // 6 < L1Q <= 10 (L1Q = 8)
0772           if (tempDR == 99.) {
0773             for (i = 0; i < nMuColl; ++i) {
0774               if (removed_row[i])
0775                 continue;
0776               theit = muColl->begin(muColl->getFirstBX()) + i;
0777               if ((theit->hwQual() <= 10) && (theit->hwQual() > 6)) {
0778                 for (j = 0; j < nOfflineSeed1; ++j) {
0779                   if (removed_col[j])
0780                     continue;
0781                   if (tempDR > dRmtx[i][j]) {
0782                     tempDR = dRmtx[i][j];
0783                     theL1 = i;
0784                     theOffs = j;
0785                   }
0786                 }
0787               }
0788             }
0789           }
0790 
0791           // L1Q <= 6 (L1Q = 4)
0792           if (tempDR == 99.) {
0793             for (i = 0; i < nMuColl; ++i) {
0794               if (removed_row[i])
0795                 continue;
0796               theit = muColl->begin(muColl->getFirstBX()) + i;
0797               if (theit->hwQual() <= 6) {
0798                 for (j = 0; j < nOfflineSeed1; ++j) {
0799                   if (removed_col[j])
0800                     continue;
0801                   if (tempDR > dRmtx[i][j]) {
0802                     tempDR = dRmtx[i][j];
0803                     theL1 = i;
0804                     theOffs = j;
0805                   }
0806                 }
0807               }
0808             }
0809           }
0810 
0811           auto it = muColl->begin(muColl->getFirstBX()) + theL1;
0812 
0813           double newDRcone = matchingDR[0];
0814           if (fabs(it->eta()) < etaBins.back()) {
0815             std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
0816             newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0817           }
0818 
0819           if (!(tempDR < newDRcone))
0820             continue;
0821 
0822           // Remove row and column for given (L1mu, offSeed)
0823           removed_col[theOffs] = true;
0824           removed_row[theL1] = true;
0825 
0826           if (selOffseeds[theL1][theOffs] != nullptr) {
0827             //put given L1mu and offSeed to output
0828             edm::OwnVector<TrackingRecHit> newContainer;
0829 
0830             PTrajectoryStateOnDet const &seedTSOS = selOffseeds[theL1][theOffs]->startingState();
0831             for (auto const &recHit : selOffseeds[theL1][theOffs]->recHits()) {
0832               newContainer.push_back(recHit);
0833             }
0834             output->push_back(L2MuonTrajectorySeed(seedTSOS,
0835                                                    newContainer,
0836                                                    alongMomentum,
0837                                                    MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0838           }
0839         }
0840       }
0841 
0842       else if (matchType == 4) {
0843         for (i = 0; i < nMuColl; ++i) {
0844           auto it = muColl->begin(muColl->getFirstBX()) + i;
0845 
0846           double tempDR = 99.;
0847           unsigned int theOffs = 0;
0848 
0849           double newDRcone = matchingDR[0];
0850           if (fabs(it->eta()) < etaBins.back()) {
0851             std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
0852             newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0853           }
0854 
0855           for (j = 0; j < nOfflineSeed1; ++j) {
0856             if (tempDR > dRmtx[i][j]) {
0857               tempDR = dRmtx[i][j];
0858               theOffs = j;
0859             }
0860           }
0861 
0862           if (!(tempDR < newDRcone))
0863             continue;
0864 
0865           if (selOffseeds[i][theOffs] != nullptr) {
0866             //put given L1mu and offSeed to output
0867             edm::OwnVector<TrackingRecHit> newContainer;
0868 
0869             PTrajectoryStateOnDet const &seedTSOS = selOffseeds[i][theOffs]->startingState();
0870             for (auto const &recHit : selOffseeds[i][theOffs]->recHits()) {
0871               newContainer.push_back(recHit);
0872             }
0873             output->push_back(L2MuonTrajectorySeed(seedTSOS,
0874                                                    newContainer,
0875                                                    alongMomentum,
0876                                                    MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0877           }
0878         }
0879       }
0880 
0881     }  // useOfflineSeed && isAsso
0882 
0883   }  // matchType > 0
0884 
0885   //--- SortType 1 : by L1 pT
0886   if (sortType == 1) {
0887     sort(output->begin(), output->end(), sortByL1Pt);
0888   }
0889 
0890   //--- SortType 2 : by L1 quality and pT
0891   else if (sortType == 2) {
0892     sort(output->begin(), output->end(), sortByL1QandPt);
0893   }
0894 
0895   iEvent.put(std::move(output));
0896 }
0897 
0898 // FIXME: does not resolve ambiguities yet!
0899 const TrajectorySeed *L2MuonSeedGeneratorFromL1T::associateOfflineSeedToL1(
0900     edm::Handle<edm::View<TrajectorySeed>> &offseeds,
0901     std::vector<int> &offseedMap,
0902     TrajectoryStateOnSurface &newTsos,
0903     double dRcone) {
0904   const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
0905   MuonPatternRecoDumper debugtmp;
0906 
0907   edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
0908   const TrajectorySeed *selOffseed = nullptr;
0909   double bestDr = 99999.;
0910   unsigned int nOffseed(0);
0911   int lastOffseed(-1);
0912 
0913   for (offseed = offseeds->begin(); offseed != endOffseed; ++offseed, ++nOffseed) {
0914     if (offseedMap[nOffseed] != 0)
0915       continue;
0916     GlobalPoint glbPos = theService->trackingGeometry()
0917                              ->idToDet(offseed->startingState().detId())
0918                              ->surface()
0919                              .toGlobal(offseed->startingState().parameters().position());
0920     GlobalVector glbMom = theService->trackingGeometry()
0921                               ->idToDet(offseed->startingState().detId())
0922                               ->surface()
0923                               .toGlobal(offseed->startingState().parameters().momentum());
0924 
0925     // Preliminary check
0926     double preDr = deltaR(newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi());
0927     if (preDr > 1.0)
0928       continue;
0929 
0930     const FreeTrajectoryState offseedFTS(
0931         glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
0932     TrajectoryStateOnSurface offseedTsos =
0933         theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
0934     LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
0935     LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
0936     LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi()
0937                        << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
0938     LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge() * offseedFTS.momentum().perp()
0939                        << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")"
0940                        << std::endl
0941                        << std::endl;
0942 
0943     if (offseedTsos.isValid()) {
0944       LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
0945       LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag()
0946                          << ", phi=" << offseedTsos.globalPosition().phi()
0947                          << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
0948       LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge() * offseedTsos.globalMomentum().perp()
0949                          << ", phi=" << offseedTsos.globalMomentum().phi()
0950                          << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl
0951                          << std::endl;
0952       double newDr = deltaR(newTsos.globalPosition().eta(),
0953                             newTsos.globalPosition().phi(),
0954                             offseedTsos.globalPosition().eta(),
0955                             offseedTsos.globalPosition().phi());
0956       LogDebug(metlabel) << "   -- DR = " << newDr << std::endl;
0957       if (newDr < dRcone && newDr < bestDr) {
0958         LogDebug(metlabel) << "          --> OK! " << newDr << std::endl << std::endl;
0959         selOffseed = &*offseed;
0960         bestDr = newDr;
0961         offseedMap[nOffseed] = 1;
0962         if (lastOffseed > -1)
0963           offseedMap[lastOffseed] = 0;
0964         lastOffseed = nOffseed;
0965       } else {
0966         LogDebug(metlabel) << "          --> Rejected. " << newDr << std::endl << std::endl;
0967       }
0968     } else {
0969       LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
0970     }
0971   }
0972 
0973   return selOffseed;
0974 }
0975 
0976 bool L2MuonSeedGeneratorFromL1T::isAssociateOfflineSeedToL1(
0977     edm::Handle<edm::View<TrajectorySeed>> &offseeds,
0978     std::vector<std::vector<double>> &dRmtx,
0979     TrajectoryStateOnSurface &newTsos,
0980     unsigned int imu,
0981     std::vector<std::vector<const TrajectorySeed *>> &selOffseeds,
0982     double dRcone) {
0983   const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
0984   bool isAssociated = false;
0985   MuonPatternRecoDumper debugtmp;
0986 
0987   edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
0988   unsigned int nOffseed(0);
0989 
0990   for (offseed = offseeds->begin(); offseed != endOffseed; ++offseed, ++nOffseed) {
0991     GlobalPoint glbPos = theService->trackingGeometry()
0992                              ->idToDet(offseed->startingState().detId())
0993                              ->surface()
0994                              .toGlobal(offseed->startingState().parameters().position());
0995     GlobalVector glbMom = theService->trackingGeometry()
0996                               ->idToDet(offseed->startingState().detId())
0997                               ->surface()
0998                               .toGlobal(offseed->startingState().parameters().momentum());
0999 
1000     // Preliminary check
1001     double preDr = deltaR(newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi());
1002     if (preDr > 1.0)
1003       continue;
1004 
1005     const FreeTrajectoryState offseedFTS(
1006         glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
1007     TrajectoryStateOnSurface offseedTsos =
1008         theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
1009     LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
1010     LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
1011     LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi()
1012                        << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
1013     LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge() * offseedFTS.momentum().perp()
1014                        << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")"
1015                        << std::endl
1016                        << std::endl;
1017 
1018     if (offseedTsos.isValid()) {
1019       LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
1020       LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag()
1021                          << ", phi=" << offseedTsos.globalPosition().phi()
1022                          << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
1023       LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge() * offseedTsos.globalMomentum().perp()
1024                          << ", phi=" << offseedTsos.globalMomentum().phi()
1025                          << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl
1026                          << std::endl;
1027       double newDr = deltaR(newTsos.globalPosition().eta(),
1028                             newTsos.globalPosition().phi(),
1029                             offseedTsos.globalPosition().eta(),
1030                             offseedTsos.globalPosition().phi());
1031 
1032       LogDebug(metlabel) << "   -- DR = " << newDr << std::endl;
1033       if (newDr < dRcone) {
1034         LogDebug(metlabel) << "          --> OK! " << newDr << std::endl << std::endl;
1035 
1036         dRmtx[imu][nOffseed] = newDr;
1037         selOffseeds[imu][nOffseed] = &*offseed;
1038 
1039         isAssociated = true;
1040       } else {
1041         LogDebug(metlabel) << "          --> Rejected. " << newDr << std::endl << std::endl;
1042       }
1043     } else {
1044       LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
1045     }
1046   }
1047 
1048   return isAssociated;
1049 }