diff --git a/RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.cc b/RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.cc index 12faaae128014..0d2f69e6175b1 100644 --- a/RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.cc +++ b/RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.cc @@ -11,6 +11,8 @@ * \author A.Everett, R.Bellan, J. Alcaraz * * ORCA's author: N. Neumeister + * + * Modified by M.Oh */ // //-------------------------------------------------- @@ -41,6 +43,25 @@ using namespace std; using namespace edm; using namespace l1t; +//--- Functions used in output sorting +bool sortByL1Pt(L2MuonTrajectorySeed &A, L2MuonTrajectorySeed &B) { + l1t::MuonRef Ref_L1A = A.l1tParticle(); + l1t::MuonRef Ref_L1B = B.l1tParticle(); + return (Ref_L1A->pt() > Ref_L1B->pt()); +}; + +bool sortByL1QandPt(L2MuonTrajectorySeed &A, L2MuonTrajectorySeed &B) { + l1t::MuonRef Ref_L1A = A.l1tParticle(); + l1t::MuonRef Ref_L1B = B.l1tParticle(); + + // Compare quality first + if (Ref_L1A->hwQual() > Ref_L1B->hwQual()) return true; + if (Ref_L1A->hwQual() < Ref_L1B->hwQual()) return false; + + // For same quality L1s compare pT + return (Ref_L1A->pt() > Ref_L1B->pt()); +}; + // constructors L2MuonSeedGeneratorFromL1T::L2MuonSeedGeneratorFromL1T(const edm::ParameterSet& iConfig) : theSource(iConfig.getParameter("InputObjects")), @@ -49,12 +70,15 @@ L2MuonSeedGeneratorFromL1T::L2MuonSeedGeneratorFromL1T(const edm::ParameterSet& theL1MinPt(iConfig.getParameter("L1MinPt")), theL1MaxEta(iConfig.getParameter("L1MaxEta")), theL1MinQuality(iConfig.getParameter("L1MinQuality")), + theMinPtBarrel(iConfig.getParameter("SetMinPtBarrelTo")), + theMinPtEndcap(iConfig.getParameter("SetMinPtEndcapTo")), useOfflineSeed(iConfig.getUntrackedParameter("UseOfflineSeed", false)), - useUnassociatedL1(iConfig.existsAs("UseUnassociatedL1") ? - iConfig.getParameter("UseUnassociatedL1") : true), + useUnassociatedL1(iConfig.getParameter("UseUnassociatedL1")), matchingDR(iConfig.getParameter >("MatchDR")), etaBins(iConfig.getParameter >("EtaMatchingBins")), - centralBxOnly_( iConfig.getParameter("CentralBxOnly") ) + centralBxOnly_( iConfig.getParameter("CentralBxOnly") ), + matchType( iConfig.getParameter("MatchType") ), + sortType( iConfig.getParameter("SortType") ) { muCollToken_ = consumes(theSource); @@ -72,6 +96,10 @@ L2MuonSeedGeneratorFromL1T::L2MuonSeedGeneratorFromL1T(const edm::ParameterSet& } + if(matchType>4 || sortType>3) { + throw cms::Exception("Configuration") << "Wrong MatchType or SortType" << endl; + } + // service parameters ParameterSet serviceParameters = iConfig.getParameter("ServiceParameters"); @@ -100,11 +128,15 @@ L2MuonSeedGeneratorFromL1T::fillDescriptions(edm::ConfigurationDescriptions& des desc.add("L1MinPt",-1.); desc.add("L1MaxEta",5.0); desc.add("L1MinQuality",0); + desc.add("SetMinPtBarrelTo",3.5); + desc.add("SetMinPtEndcapTo",1.0); desc.addUntracked("UseOfflineSeed",false); desc.add("UseUnassociatedL1", true); desc.add>("MatchDR", {0.3}); desc.add>("EtaMatchingBins", {0., 2.5}); desc.add("CentralBxOnly", true); + desc.add("MatchType", 0)->setComment("MatchType : 0 Old matching, 1 L1 Order(1to1), 2 Min dR(1to1), 3 Higher Q(1to1), 4 All matched L1"); + desc.add("SortType", 0)->setComment("SortType : 0 not sort, 1 Pt, 2 Q and Pt"); desc.addUntracked("OfflineSeedLabel", edm::InputTag("")); edm::ParameterSetDescription psd0; @@ -126,244 +158,725 @@ void L2MuonSeedGeneratorFromL1T::produce(edm::Event& iEvent, const edm::EventSet edm::Handle muColl; iEvent.getByToken(muCollToken_, muColl); - LogTrace(metname) << "Number of muons " << muColl->size() << endl; + LogDebug(metname) << "Number of muons " << muColl->size() << endl; - edm::Handle > offlineSeedHandle; - vector offlineSeedMap; - if(useOfflineSeed) { - iEvent.getByToken(offlineSeedToken_, offlineSeedHandle); - LogTrace(metname) << "Number of offline seeds " << offlineSeedHandle->size() << endl; - offlineSeedMap = vector(offlineSeedHandle->size(), 0); - } - for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) { - if (centralBxOnly_ && (ibx != 0)) continue; - for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++){ - - unsigned int quality = it->hwQual(); - int valid_charge = it->hwChargeValid(); - - float pt = it->pt(); - float eta = it->eta(); - float theta = 2*atan(exp(-eta)); - float phi = it->phi(); - int charge = it->charge(); - // Set charge=0 for the time being if the valid charge bit is zero - if (!valid_charge) charge = 0; - - int link = 36 + (int)(it -> tfMuonIndex() / 3.); - bool barrel = true; - if ( (link >= 36 && link <= 41) || (link >= 66 && link <= 71)) barrel = false; - - if ( pt < theL1MinPt || fabs(eta) > theL1MaxEta ) continue; - - LogTrace(metname) << "New L2 Muon Seed" ; - LogTrace(metname) << "Pt = " << pt << " GeV/c"; - LogTrace(metname) << "eta = " << eta; - LogTrace(metname) << "theta = " << theta << " rad"; - LogTrace(metname) << "phi = " << phi << " rad"; - LogTrace(metname) << "charge = " << charge; - LogTrace(metname) << "In Barrel? = " << barrel; - - if ( quality <= theL1MinQuality ) continue; - LogTrace(metname) << "quality = "<< quality; - - // Update the services - theService->update(iSetup); + //--- matchType 0 : Old logic + if(matchType == 0) { - const DetLayer *detLayer = nullptr; - float radius = 0.; + edm::Handle > offlineSeedHandle; + vector offlineSeedMap; + if(useOfflineSeed) { + iEvent.getByToken(offlineSeedToken_, offlineSeedHandle); + LogDebug(metname) << "Number of offline seeds " << offlineSeedHandle->size() << endl; + offlineSeedMap = vector(offlineSeedHandle->size(), 0); + } - CLHEP::Hep3Vector vec(0.,1.,0.); - vec.setTheta(theta); - vec.setPhi(phi); - - DetId theid; - // Get the det layer on which the state should be put - if ( barrel ){ - LogTrace(metname) << "The seed is in the barrel"; - - // MB2 - DetId id = DTChamberId(0,2,0); - detLayer = theService->detLayerGeometry()->idToLayer(id); - LogTrace(metname) << "L2 Layer: " << debug.dumpLayer(detLayer); - - const BoundSurface* sur = &(detLayer->surface()); - const BoundCylinder* bc = dynamic_cast(sur); - - radius = fabs(bc->radius()/sin(theta)); - theid = id; - - LogTrace(metname) << "radius "<detLayerGeometry()->idToLayer(id); - LogTrace(metname) << "L2 Layer: " << debug.dumpLayer(detLayer); - - radius = fabs(detLayer->position().z()/cos(theta)); - theid = id; - - if( pt < 1.0) pt = 1.0; + for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) { + if (centralBxOnly_ && (ibx != 0)) continue; + for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++){ + + unsigned int quality = it->hwQual(); + int valid_charge = it->hwChargeValid(); + + float pt = it->pt(); + float eta = it->eta(); + float theta = 2*atan(exp(-eta)); + float phi = it->phi(); + int charge = it->charge(); + // Set charge=0 for the time being if the valid charge bit is zero + if (!valid_charge) charge = 0; + + // link number indices of the optical fibres that connect the uGMT with the track finders + //EMTF+ : 36-41, OMTF+ : 42-47, BMTF : 48-59, OMTF- : 60-65, EMTF- : 66-71 + int link = 36 + (int)(it -> tfMuonIndex() / 3.); + bool barrel = true; + if ( (link >= 36 && link <= 41) || (link >= 66 && link <= 71)) barrel = false; + + if ( pt < theL1MinPt || fabs(eta) > theL1MaxEta ) continue; + + LogDebug(metname) << "New L2 Muon Seed" ; + LogDebug(metname) << "Pt = " << pt << " GeV/c"; + LogDebug(metname) << "eta = " << eta; + LogDebug(metname) << "theta = " << theta << " rad"; + LogDebug(metname) << "phi = " << phi << " rad"; + LogDebug(metname) << "charge = " << charge; + LogDebug(metname) << "In Barrel? = " << barrel; + + if ( quality <= theL1MinQuality ) continue; + LogDebug(metname) << "quality = "<< quality; + + // Update the services + theService->update(iSetup); + + const DetLayer *detLayer = nullptr; + float radius = 0.; + + CLHEP::Hep3Vector vec(0.,1.,0.); + vec.setTheta(theta); + vec.setPhi(phi); + + DetId theid; + // Get the det layer on which the state should be put + if ( barrel ){ + LogDebug(metname) << "The seed is in the barrel"; + + // MB2 + theid = DTChamberId(0,2,0); + detLayer = theService->detLayerGeometry()->idToLayer(theid); + LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer); + + const BoundSurface* sur = &(detLayer->surface()); + const BoundCylinder* bc = dynamic_cast(sur); + + radius = fabs(bc->radius()/sin(theta)); + + LogDebug(metname) << "radius "<detLayerGeometry()->idToLayer(theid); + LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer); + + radius = fabs(detLayer->position().z()/cos(theta)); + + if( pt < theMinPtEndcap) pt = theMinPtEndcap; + } + + vec.setMag(radius); + + GlobalPoint pos(vec.x(),vec.y(),vec.z()); + + GlobalVector mom(pt*cos(phi), pt*sin(phi), pt*cos(theta)/sin(theta)); + + GlobalTrajectoryParameters param(pos,mom,charge,&*theService->magneticField()); + AlgebraicSymMatrix55 mat; + + mat[0][0] = (0.25/pt)*(0.25/pt); // sigma^2(charge/abs_momentum) + if ( !barrel ) mat[0][0] = (0.4/pt)*(0.4/pt); + + //Assign q/pt = 0 +- 1/pt if charge has been declared invalid + if (!valid_charge) mat[0][0] = (1./pt)*(1./pt); + + mat[1][1] = 0.05*0.05; // sigma^2(lambda) + mat[2][2] = 0.2*0.2; // sigma^2(phi) + mat[3][3] = 20.*20.; // sigma^2(x_transverse)) + mat[4][4] = 20.*20.; // sigma^2(y_transverse)) + + CurvilinearTrajectoryError error(mat); + + const FreeTrajectoryState state(param,error); + + LogDebug(metname) << "Free trajectory State from the parameters"; + LogDebug(metname) << debug.dumpFTS(state); + + // Propagate the state on the MB2/ME2 surface + TrajectoryStateOnSurface tsos = theService->propagator(thePropagatorName)->propagate(state, detLayer->surface()); + + LogDebug(metname) << "State after the propagation on the layer"; + LogDebug(metname) << debug.dumpLayer(detLayer); + LogDebug(metname) << debug.dumpTSOS(tsos); + + double dRcone = matchingDR[0]; + if ( fabs(eta) < etaBins.back() ){ + std::vector::iterator lowEdge = std::upper_bound (etaBins.begin(), etaBins.end(), fabs(eta)); + dRcone = matchingDR.at( lowEdge - etaBins.begin() - 1); } - - // Fallback solution using ME2 - DetId fallback_id; - theta < Geom::pi()/2. ? fallback_id = CSCDetId(1,2,0,0,0) : fallback_id = CSCDetId(2,2,0,0,0); - const DetLayer* ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id); - vec.setMag(radius); - - GlobalPoint pos(vec.x(),vec.y(),vec.z()); - - GlobalVector mom(pt*cos(phi), pt*sin(phi), pt*cos(theta)/sin(theta)); + if (tsos.isValid()) { - GlobalTrajectoryParameters param(pos,mom,charge,&*theService->magneticField()); - AlgebraicSymMatrix55 mat; - - mat[0][0] = (0.25/pt)*(0.25/pt); // sigma^2(charge/abs_momentum) - if ( !barrel ) mat[0][0] = (0.4/pt)*(0.4/pt); + edm::OwnVector container; - //Assign q/pt = 0 +- 1/pt if charge has been declared invalid - if (!valid_charge) mat[0][0] = (1./pt)*(1./pt); - - mat[1][1] = 0.05*0.05; // sigma^2(lambda) - mat[2][2] = 0.2*0.2; // sigma^2(phi) - mat[3][3] = 20.*20.; // sigma^2(x_transverse)) - mat[4][4] = 20.*20.; // sigma^2(y_transverse)) - - CurvilinearTrajectoryError error(mat); - - const FreeTrajectoryState state(param,error); - - LogTrace(metname) << "Free trajectory State from the parameters"; - LogTrace(metname) << debug.dumpFTS(state); - - // Propagate the state on the MB2/ME2 surface - TrajectoryStateOnSurface tsos = theService->propagator(thePropagatorName)->propagate(state, detLayer->surface()); - - LogTrace(metname) << "State after the propagation on the layer"; - LogTrace(metname) << debug.dumpLayer(detLayer); - LogTrace(metname) << debug.dumpTSOS(tsos); - - double dRcone = matchingDR[0]; - if ( fabs(eta) < etaBins.back() ){ - std::vector::iterator lowEdge = std::upper_bound (etaBins.begin(), etaBins.end(), fabs(eta)); - dRcone = matchingDR.at( lowEdge - etaBins.begin() - 1); - } - - if (tsos.isValid()) { - - edm::OwnVector container; - - if(useOfflineSeed && ( !valid_charge || charge == 0) ) { - - const TrajectorySeed *assoOffseed = - associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, tsos, dRcone ); - - if(assoOffseed!=nullptr) { - PTrajectoryStateOnDet const & seedTSOS = assoOffseed->startingState(); - TrajectorySeed::const_iterator - tsci = assoOffseed->recHits().first, - tscie = assoOffseed->recHits().second; - for(; tsci!=tscie; ++tsci) { - container.push_back(*tsci); - } - output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, - MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); - } - else { - if(useUnassociatedL1) { - // convert the TSOS into a PTSOD - PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( tsos, theid.rawId()); + if(useOfflineSeed && ( !valid_charge || charge == 0) ) { + + const TrajectorySeed *assoOffseed = + associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, tsos, dRcone ); + + if(assoOffseed!=nullptr) { + PTrajectoryStateOnDet const & seedTSOS = assoOffseed->startingState(); + TrajectorySeed::const_iterator + tsci = assoOffseed->recHits().first, + tscie = assoOffseed->recHits().second; + for(; tsci!=tscie; ++tsci) { + container.push_back(*tsci); + } output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); } + else { + if(useUnassociatedL1) { + // convert the TSOS into a PTSOD + PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( tsos, theid.rawId()); + output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, + MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); + } + } } - } - else if (useOfflineSeed && valid_charge){ - // Get the compatible dets on the layer - std::vector< pair > - detsWithStates = detLayer->compatibleDets(tsos, - *theService->propagator(thePropagatorName), - *theEstimator); - - if (detsWithStates.empty() && barrel ) { - // try again to propagate but using ME2 as reference - tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface()); - detsWithStates = ME2DetLayer->compatibleDets(tsos, - *theService->propagator(thePropagatorName), - *theEstimator); - } + else if (useOfflineSeed && valid_charge){ + // Get the compatible dets on the layer + std::vector< pair > + detsWithStates = detLayer->compatibleDets(tsos, + *theService->propagator(thePropagatorName), + *theEstimator); + + if (detsWithStates.empty() && barrel ) { + // Fallback solution using ME2, try again to propagate but using ME2 as reference + DetId fallback_id; + theta < Geom::pi()/2. ? fallback_id = CSCDetId(1,2,0,0,0) : fallback_id = CSCDetId(2,2,0,0,0); + const DetLayer* ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id); + + tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface()); + detsWithStates = ME2DetLayer->compatibleDets(tsos, + *theService->propagator(thePropagatorName), + *theEstimator); + } - if (!detsWithStates.empty()){ - - TrajectoryStateOnSurface newTSOS = detsWithStates.front().second; - const GeomDet *newTSOSDet = detsWithStates.front().first; - - LogTrace(metname) << "Most compatible det"; - LogTrace(metname) << debug.dumpMuonId(newTSOSDet->geographicalId()); - - LogDebug(metname) << "L1 info: Det and State:"; - LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId()); - - if (newTSOS.isValid()){ - - //LogDebug(metname) << "(x, y, z) = (" << newTSOS.globalPosition().x() << ", " - // << newTSOS.globalPosition().y() << ", " << newTSOS.globalPosition().z() << ")"; - LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag() << ", phi=" - << newTSOS.globalPosition().phi() << ", eta=" << newTSOS.globalPosition().eta() << ")"; - LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge()*newTSOS.globalMomentum().perp() << ", phi=" - << newTSOS.globalMomentum().phi() << ", eta=" << newTSOS.globalMomentum().eta() << ")"; - - //LogDebug(metname) << "State on it"; - //LogDebug(metname) << debug.dumpTSOS(newTSOS); - - const TrajectorySeed *assoOffseed = - associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS, dRcone); - - if(assoOffseed!=nullptr) { - PTrajectoryStateOnDet const & seedTSOS = assoOffseed->startingState(); - TrajectorySeed::const_iterator - tsci = assoOffseed->recHits().first, - tscie = assoOffseed->recHits().second; - for(; tsci!=tscie; ++tsci) { - container.push_back(*tsci); + if (!detsWithStates.empty()){ + + TrajectoryStateOnSurface newTSOS = detsWithStates.front().second; + const GeomDet *newTSOSDet = detsWithStates.front().first; + + LogDebug(metname) << "Most compatible det"; + LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId()); + + if (newTSOS.isValid()){ + + LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag() << ", phi=" + << newTSOS.globalPosition().phi() << ", eta=" << newTSOS.globalPosition().eta() << ")"; + LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge()*newTSOS.globalMomentum().perp() << ", phi=" + << newTSOS.globalMomentum().phi() << ", eta=" << newTSOS.globalMomentum().eta() << ")"; + + const TrajectorySeed *assoOffseed = + associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS, dRcone); + + if(assoOffseed!=nullptr) { + PTrajectoryStateOnDet const & seedTSOS = assoOffseed->startingState(); + TrajectorySeed::const_iterator + tsci = assoOffseed->recHits().first, + tscie = assoOffseed->recHits().second; + for(; tsci!=tscie; ++tsci) { + container.push_back(*tsci); + } + output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, + MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); } - output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, - MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); + else { + if(useUnassociatedL1) { + // convert the TSOS into a PTSOD + PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( newTSOS,newTSOSDet->geographicalId().rawId()); + output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, + MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); + } + } } + } + } + else { + // convert the TSOS into a PTSOD + PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( tsos, theid.rawId()); + output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, + MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); + } + } + } + } + + } // End of matchType 0 + + //--- matchType > 0 : Updated logics + else if(matchType > 0) { + + unsigned int nMuColl = muColl->size(); + + vector< vector > dRmtx; + vector< vector > selOffseeds; + + edm::Handle > offlineSeedHandle; + if(useOfflineSeed) { + iEvent.getByToken(offlineSeedToken_, offlineSeedHandle); + unsigned int nOfflineSeed = offlineSeedHandle->size(); + LogDebug(metname) << "Number of offline seeds " << nOfflineSeed << endl; + + // Initialize dRmtx and selOffseeds + dRmtx = vector< vector >(nMuColl, vector(nOfflineSeed, 999.0)); + selOffseeds = vector< vector >(nMuColl, vector(nOfflineSeed, nullptr)); + } + + // At least one L1 muons are associated with offSeed + bool isAsso = false; + + //--- Fill dR matrix + for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) { + if (centralBxOnly_ && (ibx != 0)) continue; + for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++){ + + //Position of given L1mu + unsigned int imu = distance(muColl->begin(muColl->getFirstBX()),it); + + unsigned int quality = it->hwQual(); + int valid_charge = it->hwChargeValid(); + + float pt = it->pt(); + float eta = it->eta(); + float theta = 2*atan(exp(-eta)); + float phi = it->phi(); + int charge = it->charge(); + // Set charge=0 for the time being if the valid charge bit is zero + if (!valid_charge) charge = 0; + + // link number indices of the optical fibres that connect the uGMT with the track finders + //EMTF+ : 36-41, OMTF+ : 42-47, BMTF : 48-59, OMTF- : 60-65, EMTF- : 66-71 + int link = 36 + (int)(it -> tfMuonIndex() / 3.); + bool barrel = true; + if ( (link >= 36 && link <= 41) || (link >= 66 && link <= 71)) barrel = false; + + if ( pt < theL1MinPt || fabs(eta) > theL1MaxEta ) continue; + + LogDebug(metname) << "New L2 Muon Seed" ; + LogDebug(metname) << "Pt = " << pt << " GeV/c"; + LogDebug(metname) << "eta = " << eta; + LogDebug(metname) << "theta = " << theta << " rad"; + LogDebug(metname) << "phi = " << phi << " rad"; + LogDebug(metname) << "charge = " << charge; + LogDebug(metname) << "In Barrel? = " << barrel; + + if ( quality <= theL1MinQuality ) continue; + LogDebug(metname) << "quality = "<< quality; + + // Update the services + theService->update(iSetup); + + const DetLayer *detLayer = nullptr; + float radius = 0.; + + CLHEP::Hep3Vector vec(0.,1.,0.); + vec.setTheta(theta); + vec.setPhi(phi); + + DetId theid; + // Get the det layer on which the state should be put + if ( barrel ){ + LogDebug(metname) << "The seed is in the barrel"; + + // MB2 + theid = DTChamberId(0,2,0); + detLayer = theService->detLayerGeometry()->idToLayer(theid); + LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer); + + const BoundSurface* sur = &(detLayer->surface()); + const BoundCylinder* bc = dynamic_cast(sur); + + radius = fabs(bc->radius()/sin(theta)); + + LogDebug(metname) << "radius "<detLayerGeometry()->idToLayer(theid); + LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer); + + radius = fabs(detLayer->position().z()/cos(theta)); + + if( pt < theMinPtEndcap) pt = theMinPtEndcap; + } + + vec.setMag(radius); + + GlobalPoint pos(vec.x(),vec.y(),vec.z()); + + GlobalVector mom(pt*cos(phi), pt*sin(phi), pt*cos(theta)/sin(theta)); + + GlobalTrajectoryParameters param(pos,mom,charge,&*theService->magneticField()); + AlgebraicSymMatrix55 mat; + + mat[0][0] = (0.25/pt)*(0.25/pt); // sigma^2(charge/abs_momentum) + if ( !barrel ) mat[0][0] = (0.4/pt)*(0.4/pt); + + //Assign q/pt = 0 +- 1/pt if charge has been declared invalid + if (!valid_charge) mat[0][0] = (1./pt)*(1./pt); + + mat[1][1] = 0.05*0.05; // sigma^2(lambda) + mat[2][2] = 0.2*0.2; // sigma^2(phi) + mat[3][3] = 20.*20.; // sigma^2(x_transverse)) + mat[4][4] = 20.*20.; // sigma^2(y_transverse)) + + CurvilinearTrajectoryError error(mat); + + const FreeTrajectoryState state(param,error); + + LogDebug(metname) << "Free trajectory State from the parameters"; + LogDebug(metname) << debug.dumpFTS(state); + + // Propagate the state on the MB2/ME2 surface + TrajectoryStateOnSurface tsos = theService->propagator(thePropagatorName)->propagate(state, detLayer->surface()); + + LogDebug(metname) << "State after the propagation on the layer"; + LogDebug(metname) << debug.dumpLayer(detLayer); + LogDebug(metname) << debug.dumpTSOS(tsos); + + double dRcone = matchingDR[0]; + if ( fabs(eta) < etaBins.back() ){ + std::vector::iterator lowEdge = std::upper_bound (etaBins.begin(), etaBins.end(), fabs(eta)); + dRcone = matchingDR.at( lowEdge - etaBins.begin() - 1); + } + + if (tsos.isValid()) { + + edm::OwnVector container; + + if(useOfflineSeed) { + + if( ( !valid_charge || charge == 0 ) ) { + + bool isAssoOffseed = isAssociateOfflineSeedToL1(offlineSeedHandle, dRmtx, tsos, imu, selOffseeds, dRcone ); + + if(isAssoOffseed) { + isAsso = true; + } + + // Using old way else { if(useUnassociatedL1) { // convert the TSOS into a PTSOD - PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( newTSOS,newTSOSDet->geographicalId().rawId()); + PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( tsos, theid.rawId()); output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); } - } + } + + } + + else if( valid_charge ) { + + // Get the compatible dets on the layer + std::vector< pair > + detsWithStates = detLayer->compatibleDets(tsos, + *theService->propagator(thePropagatorName), + *theEstimator); + + if (detsWithStates.empty() && barrel ) { + // Fallback solution using ME2, try again to propagate but using ME2 as reference + DetId fallback_id; + theta < Geom::pi()/2. ? fallback_id = CSCDetId(1,2,0,0,0) : fallback_id = CSCDetId(2,2,0,0,0); + const DetLayer* ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id); + + tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface()); + detsWithStates = ME2DetLayer->compatibleDets(tsos, + *theService->propagator(thePropagatorName), + *theEstimator); + } + + if (!detsWithStates.empty()){ + + TrajectoryStateOnSurface newTSOS = detsWithStates.front().second; + const GeomDet *newTSOSDet = detsWithStates.front().first; + + LogDebug(metname) << "Most compatible det"; + LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId()); + + if (newTSOS.isValid()){ + + LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag() << ", phi=" + << newTSOS.globalPosition().phi() << ", eta=" << newTSOS.globalPosition().eta() << ")"; + LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge()*newTSOS.globalMomentum().perp() << ", phi=" + << newTSOS.globalMomentum().phi() << ", eta=" << newTSOS.globalMomentum().eta() << ")"; + + bool isAssoOffseed = isAssociateOfflineSeedToL1(offlineSeedHandle, dRmtx, newTSOS, imu, selOffseeds, dRcone ); + + if(isAssoOffseed) { + isAsso = true; + } + + // Using old way + else { + if(useUnassociatedL1) { + // convert the TSOS into a PTSOD + PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( newTSOS,newTSOSDet->geographicalId().rawId()); + output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, + MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); + } + } + + } + } } - } + } // useOfflineSeed + + else { + // convert the TSOS into a PTSOD + PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( tsos, theid.rawId()); + output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, + MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); + } + } - else { - // convert the TSOS into a PTSOD - PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( tsos, theid.rawId()); - output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum, - MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); - } - } + } } + + // MatchType : 0 Old matching, 1 L1 Order(1to1), 2 Min dR(1to1), 3 Higher Q(1to1), 4 All matched L1 + if(useOfflineSeed && isAsso) { + unsigned int nOfflineSeed1 = offlineSeedHandle->size(); + unsigned int nL1; + unsigned int i, j; // for the matrix element + + if(matchType==1) { + vector removed_col = vector(nOfflineSeed1, false); + + for(nL1=0; nL1 < nMuColl; ++nL1) { + double tempDR = 99.; + unsigned int theOffs = 0; + + for(j=0; j dRmtx[nL1][j] ) { + tempDR = dRmtx[nL1][j]; + theOffs = j; + } + } + + auto it = muColl->begin(muColl->getFirstBX()) + nL1; + + double newDRcone = matchingDR[0]; + if ( fabs(it->eta()) < etaBins.back() ){ + std::vector::iterator lowEdge = std::upper_bound (etaBins.begin(), etaBins.end(), fabs(it->eta())); + newDRcone = matchingDR.at( lowEdge - etaBins.begin() - 1); + } + + if( !(tempDR < newDRcone) ) continue; + + // Remove column for given offSeed + removed_col[theOffs] = true; + + if( selOffseeds[nL1][theOffs] != nullptr ) { + //put given L1mu and offSeed to output + edm::OwnVector newContainer; + + PTrajectoryStateOnDet const & seedTSOS = selOffseeds[nL1][theOffs] ->startingState(); + TrajectorySeed::const_iterator + tsci = selOffseeds[nL1][theOffs]->recHits().first, + tscie = selOffseeds[nL1][theOffs]->recHits().second; + for(; tsci!=tscie; ++tsci) { + newContainer.push_back(*tsci); + } + output->push_back(L2MuonTrajectorySeed(seedTSOS,newContainer,alongMomentum, + MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); + } + } + } + + else if(matchType==2) { + vector removed_row = vector(nMuColl, false); + vector removed_col = vector(nOfflineSeed1, false); + + for(nL1=0; nL1 < nMuColl; ++nL1) { + double tempDR = 99.; + unsigned int theL1 = 0; + unsigned int theOffs = 0; + + for(i=0; i dRmtx[i][j]) { + tempDR = dRmtx[i][j]; + theL1 = i; + theOffs = j; + } + } + } + + auto it = muColl->begin(muColl->getFirstBX()) + theL1; + + double newDRcone = matchingDR[0]; + if ( fabs(it->eta()) < etaBins.back() ){ + std::vector::iterator lowEdge = std::upper_bound (etaBins.begin(), etaBins.end(), fabs(it->eta())); + newDRcone = matchingDR.at( lowEdge - etaBins.begin() - 1); + } + + if( !(tempDR < newDRcone) ) continue; + + // Remove row and column for given (L1mu, offSeed) + removed_col[theOffs] = true; + removed_row[theL1] = true; + + if( selOffseeds[theL1][theOffs] != nullptr ) { + //put given L1mu and offSeed to output + edm::OwnVector newContainer; + + PTrajectoryStateOnDet const & seedTSOS = selOffseeds[theL1][theOffs] ->startingState(); + TrajectorySeed::const_iterator + tsci = selOffseeds[theL1][theOffs]->recHits().first, + tscie = selOffseeds[theL1][theOffs]->recHits().second; + for(; tsci!=tscie; ++tsci) { + newContainer.push_back(*tsci); + } + output->push_back(L2MuonTrajectorySeed(seedTSOS,newContainer,alongMomentum, + MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); + } + } + } + + else if(matchType==3) { + vector removed_row = vector(nMuColl, false); + vector removed_col = vector(nOfflineSeed1, false); + + for(nL1=0; nL1 < nMuColl; ++nL1) { + double tempDR = 99.; + unsigned int theL1 = 0; + unsigned int theOffs = 0; + auto theit = muColl->begin(muColl->getFirstBX()); + + // L1Q > 10 (L1Q = 12) + for(i=0; ibegin(muColl->getFirstBX()) + i; + if (theit->hwQual() > 10) { + for(j=0; j dRmtx[i][j]) { + tempDR = dRmtx[i][j]; + theL1 = i; + theOffs = j; + } + } + } + } + + // 6 < L1Q <= 10 (L1Q = 8) + if (tempDR == 99.) { + for(i=0; ibegin(muColl->getFirstBX()) + i; + if ( (theit->hwQual() <= 10) && (theit->hwQual() > 6) ) { + for(j=0; j dRmtx[i][j]) { + tempDR = dRmtx[i][j]; + theL1 = i; + theOffs = j; + } + } + } + } + } + + // L1Q <= 6 (L1Q = 4) + if (tempDR == 99.) { + for(i=0; ibegin(muColl->getFirstBX()) + i; + if (theit->hwQual() <= 6) { + for(j=0; j dRmtx[i][j]) { + tempDR = dRmtx[i][j]; + theL1 = i; + theOffs = j; + } + } + } + } + } + + auto it = muColl->begin(muColl->getFirstBX()) + theL1; + + double newDRcone = matchingDR[0]; + if ( fabs(it->eta()) < etaBins.back() ){ + std::vector::iterator lowEdge = std::upper_bound (etaBins.begin(), etaBins.end(), fabs(it->eta())); + newDRcone = matchingDR.at( lowEdge - etaBins.begin() - 1); + } + + if( !(tempDR < newDRcone) ) continue; + + // Remove row and column for given (L1mu, offSeed) + removed_col[theOffs] = true; + removed_row[theL1] = true; + + if( selOffseeds[theL1][theOffs] != nullptr ) { + //put given L1mu and offSeed to output + edm::OwnVector newContainer; + + PTrajectoryStateOnDet const & seedTSOS = selOffseeds[theL1][theOffs] ->startingState(); + TrajectorySeed::const_iterator + tsci = selOffseeds[theL1][theOffs]->recHits().first, + tscie = selOffseeds[theL1][theOffs]->recHits().second; + for(; tsci!=tscie; ++tsci) { + newContainer.push_back(*tsci); + } + output->push_back(L2MuonTrajectorySeed(seedTSOS,newContainer,alongMomentum, + MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); + } + } + } + + else if(matchType==4) { + + for(i=0; ibegin(muColl->getFirstBX()) + i; + + double tempDR = 99.; + unsigned int theOffs = 0; + + double newDRcone = matchingDR[0]; + if ( fabs(it->eta()) < etaBins.back() ){ + std::vector::iterator lowEdge = std::upper_bound (etaBins.begin(), etaBins.end(), fabs(it->eta())); + newDRcone = matchingDR.at( lowEdge - etaBins.begin() - 1); + } + + for(j=0; j dRmtx[i][j] ) { + tempDR = dRmtx[i][j]; + theOffs = j; + } + } + + if( !(tempDR < newDRcone) ) continue; + + if( selOffseeds[i][theOffs] != nullptr ) { + //put given L1mu and offSeed to output + edm::OwnVector newContainer; + + PTrajectoryStateOnDet const & seedTSOS = selOffseeds[i][theOffs] ->startingState(); + TrajectorySeed::const_iterator + tsci = selOffseeds[i][theOffs]->recHits().first, + tscie = selOffseeds[i][theOffs]->recHits().second; + for(; tsci!=tscie; ++tsci) { + newContainer.push_back(*tsci); + } + output->push_back(L2MuonTrajectorySeed(seedTSOS,newContainer,alongMomentum, + MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) ))); + } + } + } + + } // useOfflineSeed && isAsso + + } // matchType > 0 + + //--- SortType 1 : by L1 pT + if(sortType == 1) { + sort(output->begin(),output->end(),sortByL1Pt); } - + + //--- SortType 2 : by L1 quality and pT + else if(sortType == 2) { + sort(output->begin(),output->end(),sortByL1QandPt); + } + + iEvent.put(std::move(output)); } @@ -396,23 +909,17 @@ const TrajectorySeed* L2MuonSeedGeneratorFromL1T::associateOfflineSeedToL1( edm: TrajectoryStateOnSurface offseedTsos = theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface()); LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl; LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl; - //LogDebug(metlabel) << "(x, y, z) = (" << newTSOS.globalPosition().x() << ", " - // << newTSOS.globalPosition().y() << ", " << newTSOS.globalPosition().z() << ")" << std::endl; LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi() << ", eta=" << offseedFTS.position().eta() << ")" << std::endl; LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge()*offseedFTS.momentum().perp() << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")" << std::endl << std::endl; - //LogDebug(metlabel) << debugtmp.dumpFTS(offseedFTS) << std::endl; if(offseedTsos.isValid()) { LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl; - //LogDebug(metlabel) << "(x, y, z) = (" << offseedTsos.globalPosition().x() << ", " - // << offseedTsos.globalPosition().y() << ", " << offseedTsos.globalPosition().z() << ")" << std::endl; LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag() << ", phi=" << offseedTsos.globalPosition().phi() << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl; LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge()*offseedTsos.globalMomentum().perp() << ", phi=" << offseedTsos.globalMomentum().phi() << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl << std::endl; - //LogDebug(metlabel) << debugtmp.dumpTSOS(offseedTsos) << std::endl; double newDr = deltaR( newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), offseedTsos.globalPosition().eta(), offseedTsos.globalPosition().phi() ); LogDebug(metlabel) << " -- DR = " << newDr << std::endl; @@ -435,3 +942,70 @@ const TrajectorySeed* L2MuonSeedGeneratorFromL1T::associateOfflineSeedToL1( edm: return selOffseed; } + + +bool L2MuonSeedGeneratorFromL1T::isAssociateOfflineSeedToL1( edm::Handle > & offseeds, + std::vector< std::vector > & dRmtx, + TrajectoryStateOnSurface & newTsos, + unsigned int imu, + std::vector< std::vector > & selOffseeds, + double dRcone ) { + + const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T"; + bool isAssociated = false; + MuonPatternRecoDumper debugtmp; + + edm::View::const_iterator offseed, endOffseed = offseeds->end(); + unsigned int nOffseed(0); + + for(offseed=offseeds->begin(); offseed!=endOffseed; ++offseed, ++nOffseed) { + GlobalPoint glbPos = theService->trackingGeometry()->idToDet(offseed->startingState().detId())->surface().toGlobal(offseed->startingState().parameters().position()); + GlobalVector glbMom = theService->trackingGeometry()->idToDet(offseed->startingState().detId())->surface().toGlobal(offseed->startingState().parameters().momentum()); + + // Preliminary check + double preDr = deltaR( newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi() ); + if(preDr > 1.0) continue; + + const FreeTrajectoryState offseedFTS(glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField()); + TrajectoryStateOnSurface offseedTsos = theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface()); + LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl; + LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl; + LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << + ", phi=" << offseedFTS.position().phi() << + ", eta=" << offseedFTS.position().eta() << ")" << std::endl; + LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge()*offseedFTS.momentum().perp() << + ", phi=" << offseedFTS.momentum().phi() << + ", eta=" << offseedFTS.momentum().eta() << ")" << std::endl << std::endl; + + if(offseedTsos.isValid()) { + LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl; + LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag() << + ", phi=" << offseedTsos.globalPosition().phi() << + ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl; + LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge()*offseedTsos.globalMomentum().perp() << + ", phi=" << offseedTsos.globalMomentum().phi() << + ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl << std::endl; + double newDr = deltaR( newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), + offseedTsos.globalPosition().eta(), offseedTsos.globalPosition().phi() ); + + + LogDebug(metlabel) << " -- DR = " << newDr << std::endl; + if( newDr < dRcone ) { + LogDebug(metlabel) << " --> OK! " << newDr << std::endl << std::endl; + + dRmtx[imu][nOffseed] = newDr; + selOffseeds[imu][nOffseed] = &*offseed; + + isAssociated = true; + } + else { + LogDebug(metlabel) << " --> Rejected. " << newDr << std::endl << std::endl; + } + } + else { + LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl; + } + } + + return isAssociated; +} diff --git a/RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.h b/RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.h index 783c9cf64d808..6e151163c23fd 100644 --- a/RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.h +++ b/RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.h @@ -14,6 +14,8 @@ * \author A.Everett, R.Bellan * * ORCA's author: N. Neumeister + * + * Modified by M.Oh */ //L2MuonSeedGeneratorFromL1T //-------------------------------------------------- @@ -69,6 +71,8 @@ class L2MuonSeedGeneratorFromL1T : public edm::stream::EDProducer<> { const double theL1MinPt; const double theL1MaxEta; const unsigned theL1MinQuality; + const double theMinPtBarrel; + const double theMinPtEndcap; const bool useOfflineSeed; const bool useUnassociatedL1; std::vector matchingDR; @@ -76,6 +80,8 @@ class L2MuonSeedGeneratorFromL1T : public edm::stream::EDProducer<> { /// use central bx only muons bool centralBxOnly_; + const unsigned matchType; + const unsigned sortType; /// the event setup proxy, it takes care the services update MuonServiceProxy *theService; @@ -87,6 +93,13 @@ class L2MuonSeedGeneratorFromL1T : public edm::stream::EDProducer<> { TrajectoryStateOnSurface &, double ); + bool isAssociateOfflineSeedToL1( edm::Handle > &, + std::vector< std::vector > &, + TrajectoryStateOnSurface &, + unsigned int, + std::vector< std::vector > &, + double ); + }; #endif