Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-09 23:33:53 

0001 /**  \class phase2L2MuonSeedCreator
0002  *   See header file for a description of this class
0003  *   \author Luca Ferragina (INFN BO), 2024
0004  */
0005 
0006 #include "RecoMuon/L2MuonSeedGenerator/plugins/Phase2L2MuonSeedCreator.h"
0007 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0008 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0009 
0010 #include "MagneticField/Engine/interface/MagneticField.h"
0011 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0012 
0013 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0014 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0015 
0016 #include "DataFormats/Common/interface/OwnVector.h"
0017 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
0018 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0019 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
0020 #include "DataFormats/L1TMuonPhase2/interface/TrackerMuon.h"
0021 #include "DataFormats/Math/interface/deltaR.h"
0022 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeed.h"
0023 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
0024 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0025 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0026 #include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
0027 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0028 
0029 #include "FWCore/Framework/interface/EventSetup.h"
0030 #include "FWCore/Framework/interface/ESHandle.h"
0031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0032 
0033 #include <vdt/atan.h>
0034 #include <vdt/exp.h>
0035 #include <vdt/sincos.h>
0036 
0037 #include <vector>
0038 
0039 // Constructor
0040 Phase2L2MuonSeedCreator::Phase2L2MuonSeedCreator(const edm::ParameterSet& pset)
0041     : l1TkMuCollToken_{consumes(pset.getParameter<edm::InputTag>("inputObjects"))},
0042       cscSegmentCollToken_{consumes(pset.getParameter<edm::InputTag>("cscRecSegmentLabel"))},
0043       dtSegmentCollToken_{consumes(pset.getParameter<edm::InputTag>("dtRecSegmentLabel"))},
0044       cscGeometryToken_{esConsumes<CSCGeometry, MuonGeometryRecord>()},
0045       dtGeometryToken_{esConsumes<DTGeometry, MuonGeometryRecord>()},
0046       magneticFieldToken_{esConsumes<MagneticField, IdealMagneticFieldRecord>()},
0047       minMomentum_{pset.getParameter<double>("minPL1Tk")},
0048       maxMomentum_{pset.getParameter<double>("maxPL1Tk")},
0049       matchingPhiWindow_{pset.getParameter<double>("stubMatchDPhi")},
0050       matchingThetaWindow_{pset.getParameter<double>("stubMatchDTheta")},
0051       extrapolationDeltaPhiClose_{pset.getParameter<double>("extrapolationWindowClose")},
0052       extrapolationDeltaPhiFar_{pset.getParameter<double>("extrapolationWindowFar")},
0053       maxEtaBarrel_{pset.getParameter<double>("maximumEtaBarrel")},
0054       maxEtaOverlap_{pset.getParameter<double>("maximumEtaOverlap")},
0055       propagatorName_{pset.getParameter<string>("propagator")} {
0056   // Service parameters
0057   edm::ParameterSet serviceParameters = pset.getParameter<edm::ParameterSet>("serviceParameters");
0058   // Services
0059   service_ = std::make_unique<MuonServiceProxy>(serviceParameters, consumesCollector());
0060   estimator_ = std::make_unique<Chi2MeasurementEstimator>(pset.getParameter<double>("estimatorMaxChi2"));
0061   produces<L2MuonTrajectorySeedCollection>();
0062 }
0063 
0064 void Phase2L2MuonSeedCreator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0065   edm::ParameterSetDescription desc;
0066   desc.add<edm::InputTag>("inputObjects", edm::InputTag("l1tTkMuonsGmt"));
0067   desc.add<edm::InputTag>("cscRecSegmentLabel", edm::InputTag("hltCscSegments"));
0068   desc.add<edm::InputTag>("dtRecSegmentLabel", edm::InputTag("hltDt4DSegments"));
0069   desc.add<double>("minPL1Tk", 3.5);
0070   desc.add<double>("maxPL1Tk", 200);
0071   desc.add<double>("stubMatchDPhi", 0.05);
0072   desc.add<double>("stubMatchDTheta", 0.1);
0073   desc.add<double>("extrapolationWindowClose", 0.1);
0074   desc.add<double>("extrapolationWindowFar", 0.05);
0075   desc.add<double>("maximumEtaBarrel", 0.7);
0076   desc.add<double>("maximumEtaOverlap", 1.3);
0077   desc.add<string>("propagator", "SteppingHelixPropagatorAny");
0078 
0079   // Service parameters
0080   edm::ParameterSetDescription psd0;
0081   psd0.addUntracked<std::vector<std::string>>("Propagators", {"SteppingHelixPropagatorAny"});
0082   psd0.add<bool>("RPCLayers", true);
0083   psd0.addUntracked<bool>("UseMuonNavigation", true);
0084   desc.add<edm::ParameterSetDescription>("serviceParameters", psd0);
0085   desc.add<double>("estimatorMaxChi2", 1000.0);
0086   descriptions.add("Phase2L2MuonSeedCreator", desc);
0087 }
0088 
0089 void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0090   const std::string metname = "RecoMuon|Phase2L2MuonSeedCreator";
0091 
0092   auto output = std::make_unique<L2MuonTrajectorySeedCollection>();
0093 
0094   auto const l1TkMuColl = iEvent.getHandle(l1TkMuCollToken_);
0095 
0096   auto cscHandle = iEvent.getHandle(cscSegmentCollToken_);
0097   auto cscSegments = *cscHandle;
0098   auto dtHandle = iEvent.getHandle(dtSegmentCollToken_);
0099   auto dtSegments = *dtHandle;
0100 
0101   cscGeometry_ = iSetup.getHandle(cscGeometryToken_);
0102   dtGeometry_ = iSetup.getHandle(dtGeometryToken_);
0103 
0104   auto magneticFieldHandle = iSetup.getHandle(magneticFieldToken_);
0105 
0106   LogDebug(metname) << "Number of L1 Tracker Muons in the Event: " << l1TkMuColl->size();
0107 
0108   // Loop on all L1TkMu in event
0109   for (size_t l1TkMuIndex = 0; l1TkMuIndex != l1TkMuColl->size(); ++l1TkMuIndex) {
0110     l1t::TrackerMuonRef l1TkMuRef(l1TkMuColl, l1TkMuIndex);
0111 
0112     // Physical info of L1TkMu
0113     const float pt = l1TkMuRef->phPt();
0114     if (pt < minMomentum_) {
0115       continue;
0116     }
0117     const float eta = l1TkMuRef->phEta();
0118     const float phi = l1TkMuRef->phPhi();
0119     const int charge = l1TkMuRef->phCharge();
0120 
0121     // Calculate theta once to use it for stub-segment matching
0122     // theta == 2 * atan(e^(-eta))
0123     const float theta = 2 * vdt::fast_atanf(vdt::fast_expf(-eta));
0124 
0125     // Precompute trig functions for theta
0126     float sinTheta, cosTheta;
0127     vdt::fast_sincosf(theta, sinTheta, cosTheta);
0128 
0129     // Precompute trig functions for phi
0130     float sinPhi, cosPhi;
0131     vdt::fast_sincosf(phi, sinPhi, cosPhi);
0132 
0133     LogDebug(metname) << "L1TKMu pT: " << pt << ", eta: " << eta << ", phi: " << phi;
0134     Type muonType = overlap;
0135     if (std::abs(eta) < maxEtaBarrel_) {
0136       muonType = barrel;
0137       LogDebug(metname) << "L1TkMu found in the barrel";
0138     } else if (std::abs(eta) > maxEtaOverlap_) {
0139       muonType = endcap;
0140       LogDebug(metname) << "L1TkMu found in the endcap";
0141     }
0142 
0143     // Starting seed creation
0144     LogDebug(metname) << "Start seed creation";
0145 
0146     l1t::MuonStubRefVector stubRefs = l1TkMuRef->stubs();
0147 
0148     LogDebug(metname) << "Number of stubs per L1TkMu: " << stubRefs.size();
0149     LogDebug(metname) << "Number of DT segments in event: " << dtSegments.size();
0150     LogDebug(metname) << "Number of CSC segments in event: " << cscSegments.size();
0151 
0152     // Pairs segIndex, segQuality for matches in Barrel/Overlap/Endcap
0153     std::map<DTChamberId, std::pair<int, int>> matchesInBarrel;
0154     std::map<CSCDetId, std::pair<int, int>> matchesInEndcap;
0155 
0156     // Matching info
0157     bool atLeastOneMatch = false;
0158     bool bestInDt = false;
0159 
0160     // Variables for matches in Overlap
0161     int totalBarrelQuality = 0;
0162     int totalEndcapQuality = 0;
0163     unsigned int nDtHits = 0;
0164     unsigned int nCscHits = 0;
0165 
0166     // Loop on L1TkMu stubs to find best association to DT/CSC segments
0167     for (const auto& stub : stubRefs) {
0168 #ifdef EDM_ML_DEBUG
0169       stub->print();
0170 #endif
0171       // Separate barrel, endcap and overlap cases
0172       switch (muonType) {
0173         case barrel: {
0174           if (!stub->isBarrel()) {
0175             continue;  // skip all non-barrel stubs
0176           }
0177           // Create detId for stub
0178           DTChamberId stubId = DTChamberId(stub->etaRegion(),       // wheel
0179                                            stub->depthRegion(),     // station
0180                                            stub->phiRegion() + 1);  // sector, online to offline
0181           LogDebug(metname) << "Stub DT detId: " << stubId << ". RawId: " << stubId.rawId();
0182 
0183           auto& tmpMatch = matchingStubSegment(stubId, stub, dtSegments, theta);
0184 
0185           // Found a match -> update matching info
0186           if (tmpMatch.first != -1) {
0187             matchesInBarrel.emplace(stubId, tmpMatch);
0188             atLeastOneMatch = true;
0189             bestInDt = true;
0190           }
0191 
0192 #ifdef EDM_ML_DEBUG
0193           LogDebug(metname) << "BARREL best segments:";
0194           for (const auto& [detId, matchingPair] : matchesInBarrel) {
0195             LogDebug(metname) << "Station " << detId.station() << " (" << matchingPair.first << ", "
0196                               << matchingPair.second << ")";
0197           }
0198 #endif
0199           break;
0200         }  // End barrel
0201 
0202         case endcap: {
0203           if (!stub->isEndcap()) {
0204             continue;  // skip all non-endcap stubs
0205           }
0206           // Create detId for stub
0207           int endcap = (eta > 0) ? 1 : 2;  // CSC DetId endcap (1 -> Forward, 2 -> Backwards)
0208           CSCDetId stubId = CSCDetId(endcap,
0209                                      stub->depthRegion(),              // station
0210                                      6 - std::abs(stub->etaRegion()),  // ring, online to offline
0211                                      stub->phiRegion());               // chamber
0212           LogDebug(metname) << "Stub CSC detId: " << stubId << ". RawId: " << stubId.rawId();
0213 
0214           auto& tmpMatch = matchingStubSegment(stubId, stub, cscSegments, theta);
0215 
0216           // Found a match -> update matching info
0217           if (tmpMatch.first != -1) {
0218             matchesInEndcap.emplace(stubId, tmpMatch);
0219             atLeastOneMatch = true;
0220           }
0221 
0222 #ifdef EDM_ML_DEBUG
0223           LogDebug(metname) << "ENDCAP best segments:";
0224           for (const auto& [detId, matchingPair] : matchesInEndcap) {
0225             LogDebug(metname) << "Station " << detId.station() << " (" << matchingPair.first << ", "
0226                               << matchingPair.second << ")";
0227           }
0228 #endif
0229           break;
0230         }  // End endcap
0231 
0232         case overlap: {
0233           // Overlap runs on both DTs and CSCs and picks the best overall match
0234           if (stub->isBarrel()) {
0235             // Check DTs
0236             LogDebug(metname) << "OVERLAP stub in DTs, checking " << dtSegments.size() << " DT segments";
0237             // Create detId for stub
0238             DTChamberId stubId = DTChamberId(stub->etaRegion(),       // wheel
0239                                              stub->depthRegion(),     // station
0240                                              stub->phiRegion() + 1);  // sector, online to offline
0241             LogDebug(metname) << "Stub DT detId: " << stubId << ". RawId: " << stubId.rawId();
0242 
0243             auto& tmpMatch = matchingStubSegment(stubId, stub, dtSegments, theta);
0244             totalBarrelQuality += tmpMatch.second;
0245 
0246             // Found a match -> update matching info
0247             if (tmpMatch.first != -1) {
0248               matchesInBarrel.emplace(stubId, tmpMatch);
0249               atLeastOneMatch = true;
0250               auto dtSegment = dtSegments.begin() + tmpMatch.first;
0251               nDtHits += (dtSegment->hasPhi() ? dtSegment->phiSegment()->recHits().size() : 0);
0252               nDtHits += (dtSegment->hasZed() ? dtSegment->zSegment()->recHits().size() : 0);
0253             }
0254 
0255 #ifdef EDM_ML_DEBUG
0256             LogDebug(metname) << "OVERLAP best segments in DTs:";
0257             for (auto& [detId, matchingPair] : matchesInBarrel) {
0258               LogDebug(metname) << "Station " << detId.station() << " (" << matchingPair.first << ", "
0259                                 << matchingPair.second << ")";
0260             }
0261 #endif
0262           } else if (stub->isEndcap()) {
0263             // Check CSCs
0264             LogDebug(metname) << "OVERLAP stub in CSCs, checking " << cscSegments.size() << " CSC segments";
0265             int endcap = (eta > 0) ? 1 : 2;  // CSC DetId endcap (1 -> Forward, 2 -> Backwards)
0266             CSCDetId stubId = CSCDetId(endcap,
0267                                        stub->depthRegion(),              // station
0268                                        6 - std::abs(stub->etaRegion()),  // ring, online to offline
0269                                        stub->phiRegion());               // chamber
0270             LogDebug(metname) << "Stub CSC detId: " << stubId << ". RawId: " << stubId.rawId();
0271 
0272             auto& tmpMatch = matchingStubSegment(stubId, stub, cscSegments, theta);
0273             totalEndcapQuality += tmpMatch.second;
0274 
0275             // Found a match -> update matching info
0276             if (tmpMatch.first != -1) {
0277               matchesInEndcap.emplace(stubId, tmpMatch);
0278               atLeastOneMatch = true;
0279               auto cscSegment = cscSegments.begin() + tmpMatch.first;
0280               nCscHits += cscSegment->nRecHits();
0281             }
0282 
0283 #ifdef EDM_ML_DEBUG
0284             LogDebug(metname) << "OVERLAP best segments in CSCs:";
0285             for (auto& [detId, matchingPair] : matchesInEndcap) {
0286               LogDebug(metname) << "Station " << detId.station() << " (" << matchingPair.first << ", "
0287                                 << matchingPair.second << ")";
0288             }
0289 #endif
0290           }
0291 
0292           LogDebug(metname) << "OVERLAP comparing total qualities. DT: " << totalBarrelQuality
0293                             << ", CSC: " << totalEndcapQuality;
0294 
0295           // Pick segments in DTs / CSCs based on quality
0296           bestInDt = (totalBarrelQuality > totalEndcapQuality) ? true : false;
0297 
0298           // Same qualities, pick higher number of hits
0299           if (totalBarrelQuality == totalEndcapQuality and totalBarrelQuality > -1) {
0300             LogDebug(metname) << "Same quality " << totalBarrelQuality << ". Checking total number of hits";
0301             LogDebug(metname) << "DT hits: " << nDtHits << ", CSC hits: " << nCscHits;
0302             LogDebug(metname) << (nDtHits > nCscHits ? "More hits in DT segment" : "More hits in CSC segment");
0303             bestInDt = (nDtHits >= nCscHits) ? true : false;
0304           }
0305 #ifdef EDM_ML_DEBUG
0306           LogDebug(metname) << "OVERLAP best segments:";
0307           if (bestInDt) {
0308             LogDebug(metname) << "OVERLAP best match in DTs:";
0309             for (auto& [detId, matchingPair] : matchesInBarrel) {
0310               LogDebug(metname) << "Station " << detId.station() << " (" << matchingPair.first << ", "
0311                                 << matchingPair.second << ")";
0312             }
0313           } else if (!bestInDt) {
0314             LogDebug(metname) << "OVERLAP best match in CSCs:";
0315             for (auto& [detId, matchingPair] : matchesInEndcap) {
0316               LogDebug(metname) << "Station " << detId.station() << " (" << matchingPair.first << ", "
0317                                 << matchingPair.second << ")";
0318             }
0319           }
0320 #endif
0321           break;
0322         }  // End overlap
0323 
0324         default:
0325           edm::LogError("L1TkMu must be either barrel, endcap or overlap");
0326           break;
0327       }
0328     }  // End loop on stubs
0329 
0330     // Emplace seeds in output
0331     if (!atLeastOneMatch) {
0332       LogDebug(metname) << "No matching stub found, skipping seed";
0333       continue;  // skip unmatched L1TkMu
0334     } else {
0335       // Info for propagation to MB2 or ME2
0336       service_->update(iSetup);
0337       const DetLayer* detLayer = nullptr;
0338       float radius = 0.;
0339 
0340       DetId propagateToId;
0341 
0342       edm::OwnVector<TrackingRecHit> container;
0343       if (bestInDt) {
0344         // Found at least one matching segment in DT -> propagate to Muon Barrel 2 (MB2) for track finding
0345         LogDebug(metname) << "Found matching segment(s) in DTs, propagating L1TkMu info to MB2 to seed";
0346         // MB2 detId
0347         propagateToId = DTChamberId(0, 2, 0);
0348         detLayer = service_->detLayerGeometry()->idToLayer(propagateToId);
0349         const BoundSurface* sur = &(detLayer->surface());
0350         const BoundCylinder* bc = dynamic_cast<const BoundCylinder*>(sur);
0351         radius = std::abs(bc->radius() / sinTheta);
0352 
0353         // Propagate matched segments to the seed and try to extrapolate in unmatched chambers
0354         std::vector<int> matchedStations;
0355         matchedStations.reserve(4);
0356         for (auto& [detId, matchingPair] : matchesInBarrel) {
0357           // Add matched segments to the seed
0358           LogDebug(metname) << "Adding matched DT segment in station " << detId.station() << " to the seed";
0359           container.push_back(dtSegments[matchingPair.first]);
0360           matchedStations.push_back(detId.station());
0361         }
0362         // Loop over all barrel muon stations (1-4)
0363         for (int station = DTChamberId::minStationId; station <= DTChamberId::maxStationId; ++station) {
0364           if (std::find(matchedStations.begin(), matchedStations.end(), station) == matchedStations.end()) {
0365             // Try to extrapolate from stations with a match to the ones without
0366             LogDebug(metname) << "No matching DT segment found in station " << station;
0367             auto extrapolatedMatch = extrapolateToNearbyStation(station, matchesInBarrel, dtSegments);
0368             if (extrapolatedMatch.first != -1) {
0369               LogDebug(metname) << "Adding extrapolated DT segment " << extrapolatedMatch.first << " with quality "
0370                                 << extrapolatedMatch.second << " found in station " << station << " to the seed";
0371               container.push_back(dtSegments[extrapolatedMatch.first]);
0372             }
0373           }
0374         }
0375       } else if (!bestInDt) {
0376         // Found a matching segment in CSC -> propagate to Muon Endcap 2 (ME2) for track finding
0377         LogDebug(metname) << "Found matching segment(s) in CSCs, propagating L1TkMu info to ME2 to seed";
0378         // ME2 detId
0379         propagateToId = eta > 0 ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
0380         detLayer = service_->detLayerGeometry()->idToLayer(propagateToId);
0381         radius = std::abs(detLayer->position().z() / cosTheta);
0382 
0383         // Fill seed with matched segment(s)
0384         for (auto& [detId, matchingPair] : matchesInEndcap) {
0385           LogDebug(metname) << "Adding matched CSC segment in station " << detId.station() << " to the seed";
0386           container.push_back(cscSegments[matchingPair.first]);
0387         }
0388       }
0389       // Get Global point and direction
0390       GlobalPoint pos(radius * cosPhi * sinTheta, radius * sinPhi * sinTheta, radius * cosTheta);
0391       GlobalVector mom(pt * cosPhi, pt * sinPhi, pt * cosTheta / sinTheta);
0392 
0393       GlobalTrajectoryParameters param(pos, mom, charge, &*magneticFieldHandle);
0394 
0395       AlgebraicSymMatrix55 mat;
0396 
0397       mat[0][0] = bestInDt ? (0.25 / pt) * (0.25 / pt) : (0.4 / pt) * (0.4 / pt);  // sigma^2(charge/abs_momentum)
0398       mat[1][1] = 0.05 * 0.05;                                                     // sigma^2(lambda)
0399       mat[2][2] = 0.2 * 0.2;                                                       // sigma^2(phi)
0400       mat[3][3] = 20. * 20.;                                                       // sigma^2(x_transverse)
0401       mat[4][4] = 20. * 20.;                                                       // sigma^2(y_transverse)
0402 
0403       CurvilinearTrajectoryError error(mat);
0404 
0405       const FreeTrajectoryState state(param, error);
0406 
0407       // Create the TrajectoryStateOnSurface
0408       TrajectoryStateOnSurface tsos = service_->propagator(propagatorName_)->propagate(state, detLayer->surface());
0409       // Find valid detectors with states
0410       auto detsWithStates = detLayer->compatibleDets(tsos, *service_->propagator(propagatorName_), *estimator_);
0411       // Check that at least one valid detector was found
0412       if (!detsWithStates.empty()) {
0413         // Update the detId with the one from the first valid detector with measurments found
0414         propagateToId = detsWithStates.front().first->geographicalId();
0415         // Create the Trajectory State on that detector's surface
0416         tsos = detsWithStates.front().second;
0417       } else if (detsWithStates.empty() and bestInDt) {
0418         // Propagation to MB2 failed, fallback to ME2 (might be an overlap edge case)
0419         LogDebug(metname) << "Warning: detsWithStates collection is empty for a barrel collection. Falling back to ME2";
0420         // Get ME2 DetLayer
0421         DetId fallback_id = eta > 0 ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
0422         const DetLayer* ME2DetLayer = service_->detLayerGeometry()->idToLayer(fallback_id);
0423         // Trajectory state on ME2 DetLayer
0424         tsos = service_->propagator(propagatorName_)->propagate(state, ME2DetLayer->surface());
0425         // Find the detectors with states on ME2
0426         detsWithStates = ME2DetLayer->compatibleDets(tsos, *service_->propagator(propagatorName_), *estimator_);
0427       }
0428       // Use the valid detector found to produce the persistentState for the seed
0429       if (!detsWithStates.empty()) {
0430         LogDebug(metname) << "Found a compatible detWithStates";
0431         TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
0432         const GeomDet* newTSOSDet = detsWithStates.front().first;
0433         LogDebug(metname) << "Most compatible detector: " << newTSOSDet->geographicalId().rawId();
0434         if (newTSOS.isValid()) {
0435           LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
0436                             << ", phi=" << newTSOS.globalPosition().phi() << ", eta=" << newTSOS.globalPosition().eta()
0437                             << ")";
0438           LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
0439                             << ", phi=" << newTSOS.globalMomentum().phi() << ", eta=" << newTSOS.globalMomentum().eta()
0440                             << ")";
0441           // Transform the TrajectoryStateOnSurface in a Persistent TrajectoryStateOnDet
0442           const PTrajectoryStateOnDet& seedTSOS =
0443               trajectoryStateTransform::persistentState(newTSOS, newTSOSDet->geographicalId().rawId());
0444 
0445           // Emplace seed in output
0446           LogDebug(metname) << "Emplacing seed in output";
0447           output->emplace_back(L2MuonTrajectorySeed(seedTSOS, container, alongMomentum, l1TkMuRef));
0448         }
0449       }
0450     }  // End seed emplacing (one seed per L1TkMu)
0451   }  // End loop on L1TkMu
0452   LogDebug(metname) << "All L1TkMu in event processed";
0453   iEvent.put(std::move(output));
0454 }
0455 
0456 // In DT station 4 the top and bottom sectors are made of two chambers
0457 // due to material requirements. Online is not split:
0458 // Online sector 4 == offline sector 4 or 13, Online sector 10 == offline sector 10 or 14
0459 const std::vector<DTChamberId> Phase2L2MuonSeedCreator::matchingIds(const DTChamberId& stubId) const {
0460   std::vector<DTChamberId> matchingDtIds;
0461   matchingDtIds.reserve(2);
0462   matchingDtIds.push_back(stubId);
0463   if (stubId.station() == 4) {
0464     if (stubId.sector() == 4) {
0465       matchingDtIds.emplace_back(DTChamberId(stubId.wheel(), stubId.station(), 13));
0466     }
0467     if (stubId.sector() == 10) {
0468       matchingDtIds.emplace_back(DTChamberId(stubId.wheel(), stubId.station(), 14));
0469     }
0470   }
0471   return matchingDtIds;
0472 }
0473 
0474 // Pair bestSegIndex, quality for DT segments matching
0475 const std::pair<int, int> Phase2L2MuonSeedCreator::matchingStubSegment(const DTChamberId& stubId,
0476                                                                        const l1t::MuonStubRef stub,
0477                                                                        const DTRecSegment4DCollection& segments,
0478                                                                        const float l1TkMuTheta) const {
0479   const std::string metname = "RecoMuon|Phase2L2MuonSeedCreator";
0480 
0481   int bestSegIndex = -1;
0482   int quality = -1;
0483   unsigned int nHitsPhiBest = 0;
0484   unsigned int nHitsThetaBest = 0;
0485 
0486   LogDebug(metname) << "Matching stub with DT segment";
0487   int nMatchingIds = 0;
0488 
0489   for (const DTChamberId& id : matchingIds(stubId)) {
0490     DTRecSegment4DCollection::range segmentsInChamber = segments.get(id);
0491     for (DTRecSegment4DCollection::const_iterator segment = segmentsInChamber.first;
0492          segment != segmentsInChamber.second;
0493          ++segment) {
0494       ++nMatchingIds;
0495       DTChamberId segId = segment->chamberId();
0496       LogDebug(metname) << "Segment DT detId: " << segId << ". RawId: " << segId.rawId();
0497 
0498       // Global position of the segment
0499       GlobalPoint segPos = dtGeometry_->idToDet(segId)->toGlobal(segment->localPosition());
0500 
0501       // Check delta phi
0502       double deltaPhi = std::abs(segPos.phi() - stub->offline_coord1());
0503       LogDebug(metname) << "deltaPhi: " << deltaPhi;
0504 
0505       double deltaTheta = std::abs(segPos.theta() - l1TkMuTheta);
0506       LogDebug(metname) << "deltaTheta: " << deltaTheta;
0507 
0508       // Skip segments outside phi window or very far in the theta view
0509       if (deltaPhi > matchingPhiWindow_ or deltaTheta > 4 * matchingThetaWindow_) {
0510         continue;
0511       }
0512 
0513       // Inside phi window -> check hit multiplicity
0514       unsigned int nHitsPhi = (segment->hasPhi() ? segment->phiSegment()->recHits().size() : 0);
0515       unsigned int nHitsTheta = (segment->hasZed() ? segment->zSegment()->recHits().size() : 0);
0516       LogDebug(metname) << "DT found match in deltaPhi: " << std::distance(segments.begin(), segment) << " with "
0517                         << nHitsPhi << " hits in phi and " << nHitsTheta << " hits in theta";
0518 
0519       if (nHitsPhi == nHitsPhiBest and segment->hasZed()) {
0520         // Same phi hit multiplicity -> check delta theta
0521         LogDebug(metname) << "DT found segment with same hits in phi as previous best (" << nHitsPhiBest
0522                           << "), checking theta window";
0523 
0524         // More precise check in theta window
0525         if (deltaTheta > matchingThetaWindow_) {
0526           continue;  // skip segments outside theta window
0527         }
0528 
0529         LogDebug(metname) << "DT found match in deltaTheta: " << std::distance(segments.begin(), segment) << " with "
0530                           << nHitsPhi << " hits in phi and " << nHitsTheta << " hits in theta";
0531 
0532         // Inside theta window -> check hit multiplicity (theta)
0533         if (nHitsTheta > nHitsThetaBest) {
0534           // More hits in theta -> update bestSegment and quality
0535           LogDebug(metname) << "DT found segment with more hits in theta than previous best";
0536           bestSegIndex = std::distance(segments.begin(), segment);
0537           quality = 2;
0538           LogDebug(metname) << "DT updating bestSegIndex (nHitsTheta): " << bestSegIndex << " with "
0539                             << nHitsPhi + nHitsTheta << ">" << nHitsPhiBest + nHitsThetaBest
0540                             << " total hits and quality " << quality;
0541           nHitsThetaBest = nHitsTheta;
0542         }
0543       } else if (nHitsPhi > nHitsPhiBest) {
0544         // More hits in phi -> update bestSegment and quality
0545         LogDebug(metname) << "DT found segment with more hits in phi than previous best";
0546         bestSegIndex = std::distance(segments.begin(), segment);
0547         quality = 1;
0548         LogDebug(metname) << "DT updating bestSegIndex (nHitsPhi): " << bestSegIndex << " with " << nHitsPhi << ">"
0549                           << nHitsPhiBest << " hits in phi, " << nHitsTheta << " hits in theta and quality " << quality;
0550         nHitsPhiBest = nHitsPhi;
0551         nHitsThetaBest = nHitsTheta;
0552       }
0553     }  // End loop on segments
0554   }
0555 
0556   LogDebug(metname) << "DT looped over " << nMatchingIds << (nMatchingIds > 1 ? " segments" : " segment")
0557                     << " with same DT detId as stub";
0558 
0559   if (quality < 0) {
0560     LogDebug(metname) << "DT proposed match: " << bestSegIndex << " with quality " << quality << ". Not good enough!";
0561     return std::make_pair(-1, -1);
0562   } else {
0563     LogDebug(metname) << "Found DT segment match";
0564     LogDebug(metname) << "New DT segment: " << bestSegIndex << " with " << nHitsPhiBest + nHitsThetaBest
0565                       << " total hits and quality " << quality;
0566     return std::make_pair(bestSegIndex, quality);
0567   }
0568 }
0569 
0570 // Match online-level CSCDetIds to offline labels
0571 const std::vector<CSCDetId> Phase2L2MuonSeedCreator::matchingIds(const CSCDetId& stubId) const {
0572   std::vector<CSCDetId> matchingCscIds;
0573   matchingCscIds.push_back(stubId);
0574 
0575   if (stubId.station() == 1 and stubId.ring() == 1) {
0576     matchingCscIds.emplace_back(CSCDetId(stubId.endcap(), stubId.station(), 4, stubId.chamber()));
0577   }
0578 
0579   return matchingCscIds;
0580 }
0581 
0582 // Pair bestSegIndex, quality for CSC segments matching
0583 const std::pair<int, int> Phase2L2MuonSeedCreator::matchingStubSegment(const CSCDetId& stubId,
0584                                                                        const l1t::MuonStubRef stub,
0585                                                                        const CSCSegmentCollection& segments,
0586                                                                        const float l1TkMuTheta) const {
0587   const std::string metname = "RecoMuon|Phase2L2MuonSeedCreator";
0588 
0589   int bestSegIndex = -1;
0590   int quality = -1;
0591   unsigned int nHitsBest = 0;
0592 
0593   LogDebug(metname) << "Matching stub with CSC segment";
0594   int nMatchingIds = 0;
0595   for (CSCDetId id : matchingIds(stubId)) {
0596     CSCSegmentCollection::range segmentsInChamber = segments.get(id);
0597     for (CSCSegmentCollection::const_iterator segment = segmentsInChamber.first; segment != segmentsInChamber.second;
0598          ++segment) {
0599       ++nMatchingIds;
0600       CSCDetId segId = segment->cscDetId();
0601       LogDebug(metname) << "Segment CSC detId: " << segId << ". RawId: " << segId.rawId();
0602 
0603       // Global position of the segment
0604       GlobalPoint segPos = cscGeometry_->idToDet(segId)->toGlobal(segment->localPosition());
0605 
0606       // Check delta phi
0607       double deltaPhi = std::abs(segPos.phi() - stub->offline_coord1());
0608       LogDebug(metname) << "deltaPhi: " << deltaPhi;
0609 
0610       double deltaTheta = std::abs(segPos.theta() - l1TkMuTheta);
0611       LogDebug(metname) << "deltaTheta: " << deltaTheta;
0612 
0613       // Theta mainly used in cases where multiple matches are found
0614       // to keep only the best one. Still skip segments way outside
0615       // a reasonable window
0616       const double roughThetaWindow = 0.4;
0617       if (deltaPhi > matchingPhiWindow_ or deltaTheta > roughThetaWindow) {
0618         continue;  // skip segments outside phi window
0619       }
0620 
0621       // Inside phi window -> check hit multiplicity
0622       unsigned int nHits = segment->nRecHits();
0623       LogDebug(metname) << "CSC found match in deltaPhi: " << std::distance(segments.begin(), segment) << " with "
0624                         << nHits << " hits";
0625 
0626       if (nHits == nHitsBest) {
0627         // Same hit multiplicity -> check delta theta
0628         LogDebug(metname) << "Found CSC segment with same hits (" << nHitsBest
0629                           << ") as previous best, checking theta window";
0630 
0631         if (deltaTheta > matchingThetaWindow_) {
0632           continue;  // skip segments outside theta window
0633         }
0634 
0635         // Inside theta window -> update bestSegment and quality
0636         bestSegIndex = std::distance(segments.begin(), segment);
0637         quality = 1;
0638         LogDebug(metname) << "CSC found match in deltaTheta: " << bestSegIndex << " with " << nHits
0639                           << " hits and quality " << quality;
0640       } else if (nHits > nHitsBest) {
0641         // More hits -> update bestSegment and quality
0642         bestSegIndex = std::distance(segments.begin(), segment);
0643         quality = 2;
0644         LogDebug(metname) << "Found CSC segment with more hits. Index: " << bestSegIndex << " with " << nHits << ">"
0645                           << nHitsBest << " hits and quality " << quality;
0646         nHitsBest = nHits;
0647       }
0648     }  // End loop on segments
0649   }
0650 
0651   LogDebug(metname) << "CSC looped over " << nMatchingIds << (nMatchingIds != 1 ? " segments" : " segment")
0652                     << " with same CSC detId as stub";
0653 
0654   if (quality < 0) {
0655     LogDebug(metname) << "CSC proposed match: " << bestSegIndex << " with quality " << quality << ". Not good enough!";
0656     return std::make_pair(-1, -1);
0657   } else {
0658     LogDebug(metname) << "Found CSC segment match";
0659     LogDebug(metname) << "New CSC segment: " << bestSegIndex << " with " << nHitsBest << " hits and quality "
0660                       << quality;
0661     return std::make_pair(bestSegIndex, quality);
0662   }
0663 }
0664 
0665 const std::pair<int, int> Phase2L2MuonSeedCreator::extrapolateToNearbyStation(
0666     const int endingStation,
0667     const std::map<DTChamberId, std::pair<int, int>>& matchesInBarrel,
0668     const DTRecSegment4DCollection& segments) const {
0669   const std::string metname = "RecoMuon|Phase2L2MuonSeedCreator";
0670 
0671   std::pair<int, int> extrapolatedMatch = std::make_pair(-1, -1);
0672   bool foundExtrapolatedMatch = false;
0673   switch (endingStation) {
0674     case 1: {
0675       // Station 1. Extrapolate 2->1 or 3->1 (4->1)
0676       int startingStation = 2;
0677       while (startingStation < 5) {
0678         for (auto& [detId, matchingPair] : matchesInBarrel) {
0679           if (detId.station() == startingStation) {
0680             LogDebug(metname) << "Extrapolating from station " << startingStation << " to station " << endingStation;
0681             extrapolatedMatch = extrapolateMatch(matchingPair.first, endingStation, segments);
0682             if (extrapolatedMatch.first != -1) {
0683               LogDebug(metname) << "Found extrapolated match in station " << endingStation << " from station "
0684                                 << startingStation;
0685               foundExtrapolatedMatch = true;
0686               break;
0687             }
0688           }
0689         }
0690         if (foundExtrapolatedMatch) {
0691           break;
0692         }
0693         ++startingStation;
0694       }
0695       break;
0696     }
0697     case 2: {
0698       // Station 2. Extrapolate 1->2 or 3->2 (4->2)
0699       int startingStation = 1;
0700       while (startingStation < 5) {
0701         for (auto& [detId, matchingPair] : matchesInBarrel) {
0702           if (detId.station() == startingStation) {
0703             LogDebug(metname) << "Extrapolating from station " << startingStation << " to station " << endingStation;
0704             extrapolatedMatch = extrapolateMatch(matchingPair.first, endingStation, segments);
0705             if (extrapolatedMatch.first != -1) {
0706               LogDebug(metname) << "Found extrapolated match in station " << endingStation << " from station "
0707                                 << startingStation;
0708               foundExtrapolatedMatch = true;
0709               break;
0710             }
0711           }
0712         }
0713         if (foundExtrapolatedMatch) {
0714           break;
0715         }
0716         startingStation = startingStation == 1 ? startingStation + 2 : startingStation + 1;
0717       }
0718       break;
0719     }
0720     case 3: {
0721       // Station 3. Extrapolate 2->3 or 4->3 (1->3)
0722       int startingStation = 2;
0723       while (startingStation > 0) {
0724         for (auto& [detId, matchingPair] : matchesInBarrel) {
0725           if (detId.station() == startingStation) {
0726             LogDebug(metname) << "Extrapolating from station " << startingStation << " to station " << endingStation;
0727             extrapolatedMatch = extrapolateMatch(matchingPair.first, endingStation, segments);
0728             if (extrapolatedMatch.first != -1) {
0729               LogDebug(metname) << "Found extrapolated match in station " << endingStation << " from station "
0730                                 << startingStation;
0731               foundExtrapolatedMatch = true;
0732               break;
0733             }
0734           }
0735         }
0736         if (foundExtrapolatedMatch) {
0737           break;
0738         }
0739         startingStation = startingStation == 2 ? startingStation + 2 : startingStation - 3;
0740       }
0741       break;
0742     }
0743     case 4: {
0744       // Station 4. Extrapolate 2->4 or 3->4 (1->4)
0745       int startingStation = 2;
0746       while (startingStation > 0) {
0747         for (auto& [detId, matchingPair] : matchesInBarrel) {
0748           if (detId.station() == startingStation) {
0749             LogDebug(metname) << "Extrapolating from station " << startingStation << " to station " << endingStation;
0750             extrapolatedMatch = extrapolateMatch(matchingPair.first, endingStation, segments);
0751             if (extrapolatedMatch.first != -1) {
0752               LogDebug(metname) << "Found extrapolated match in station " << endingStation << " from station "
0753                                 << startingStation;
0754               foundExtrapolatedMatch = true;
0755               break;
0756             }
0757           }
0758         }
0759         if (foundExtrapolatedMatch) {
0760           break;
0761         }
0762         startingStation = startingStation == 2 ? startingStation + 1 : startingStation - 2;
0763       }
0764       break;
0765     }
0766     default:
0767       std::cerr << "Muon stations only go from 1 to 4" << std::endl;
0768       break;
0769   }  // end endingStation switch
0770   return extrapolatedMatch;
0771 }
0772 
0773 const std::pair<int, int> Phase2L2MuonSeedCreator::extrapolateMatch(const int bestStartingSegIndex,
0774                                                                     const int endingStation,
0775                                                                     const DTRecSegment4DCollection& segments) const {
0776   const std::string metname = "RecoMuon|Phase2L2MuonSeedCreator";
0777 
0778   const auto& segmentInStartingStation = segments.begin() + bestStartingSegIndex;
0779   auto matchId = segmentInStartingStation->chamberId();
0780   GlobalPoint matchPos = dtGeometry_->idToDet(matchId)->toGlobal(segmentInStartingStation->localPosition());
0781 
0782   int bestSegIndex = -1;
0783   int quality = -1;
0784   unsigned int nHitsPhiBest = 0;
0785   unsigned int nHitsThetaBest = 0;
0786 
0787   // Find possible extrapolation from startingStation to endingStation
0788   for (DTRecSegment4DCollection::const_iterator segment = segments.begin(), last = segments.end(); segment != last;
0789        ++segment) {
0790     auto segId = segment->chamberId();
0791 
0792     if (segId.station() != endingStation) {
0793       continue;  // skip segments outside of endingStation
0794     }
0795 
0796     // Global positions of the segment
0797     GlobalPoint segPos = dtGeometry_->idToDet(segId)->toGlobal(segment->localPosition());
0798 
0799     double deltaPhi = std::abs(segPos.phi() - matchPos.phi());
0800     LogDebug(metname) << "Extrapolation deltaPhi: " << deltaPhi;
0801 
0802     double deltaTheta = std::abs(segPos.theta() - matchPos.theta());
0803     LogDebug(metname) << "Extrapolation deltaTheta: " << deltaTheta;
0804 
0805     double matchingDeltaPhi =
0806         std::abs(matchId.station() - endingStation) == 1 ? extrapolationDeltaPhiClose_ : extrapolationDeltaPhiFar_;
0807 
0808     // Theta mainly used in cases where multiple matches are found
0809     // to keep only the best one. Still skip segments way outside
0810     // a reasonable window
0811     const double roughThetaWindow = 0.4;
0812     if (deltaPhi > matchingDeltaPhi or deltaTheta > roughThetaWindow) {
0813       continue;
0814     }
0815 
0816     // Inside phi window -> check hit multiplicity
0817     unsigned int nHitsPhi = (segment->hasPhi() ? segment->phiSegment()->recHits().size() : 0);
0818     unsigned int nHitsTheta = (segment->hasZed() ? segment->zSegment()->recHits().size() : 0);
0819     LogDebug(metname) << "Extrapolation found match in deltaPhi: " << std::distance(segments.begin(), segment)
0820                       << " with " << nHitsPhi << " hits in phi and " << nHitsTheta << " hits in theta";
0821 
0822     if (nHitsPhi == nHitsPhiBest and segment->hasZed()) {
0823       // Same phi hit multiplicity -> check delta theta
0824       LogDebug(metname) << "Extrapolation found segment with same hits in phi as previous best (" << nHitsPhiBest
0825                         << "), checking theta window";
0826       double deltaTheta = std::abs(segPos.theta() - matchPos.theta());
0827       LogDebug(metname) << "Extrapolation deltaTheta: " << deltaTheta;
0828 
0829       if (deltaTheta > matchingThetaWindow_) {
0830         continue;  // skip segments outside theta window
0831       }
0832 
0833       LogDebug(metname) << "Extrapolation found match in deltaTheta: " << std::distance(segments.begin(), segment)
0834                         << " with " << nHitsPhi << " hits in phi and " << nHitsTheta << " hits in theta";
0835 
0836       // Inside theta window -> check hit multiplicity (theta)
0837       if (nHitsTheta > nHitsThetaBest) {
0838         // More hits in theta -> update bestSegment and quality
0839         LogDebug(metname) << "Extrapolation found segment with more hits in theta than previous best";
0840         bestSegIndex = std::distance(segments.begin(), segment);
0841         quality = 2;
0842         LogDebug(metname) << "Extrapolation updating bestSegIndex (nHitsTheta): " << bestSegIndex << " with "
0843                           << nHitsPhi + nHitsTheta << ">" << nHitsPhiBest + nHitsThetaBest << " total hits and quality "
0844                           << quality;
0845         nHitsThetaBest = nHitsTheta;
0846       }
0847     } else if (nHitsPhi > nHitsPhiBest) {
0848       // More hits in phi -> update bestSegment and quality
0849       LogDebug(metname) << "Extrapolation found segment with more hits in phi than previous best";
0850       bestSegIndex = std::distance(segments.begin(), segment);
0851       quality = 1;
0852       LogDebug(metname) << "Extrapolation updating bestSegIndex (nHitsPhi): " << bestSegIndex << " with " << nHitsPhi
0853                         << ">" << nHitsPhiBest << " hits in phi, " << nHitsTheta << " hits in theta and quality "
0854                         << quality;
0855       nHitsPhiBest = nHitsPhi;
0856       nHitsThetaBest = nHitsTheta;
0857     }
0858   }  // end loop on segments
0859   return std::make_pair(bestSegIndex, quality);
0860 }
0861 
0862 DEFINE_FWK_MODULE(Phase2L2MuonSeedCreator);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!