diff --git a/PhysicsTools/PatAlgos/plugins/PATJetUpdater.cc b/PhysicsTools/PatAlgos/plugins/PATJetUpdater.cc index 8f9d0bc09cd94..147cb01eb372e 100755 --- a/PhysicsTools/PatAlgos/plugins/PATJetUpdater.cc +++ b/PhysicsTools/PatAlgos/plugins/PATJetUpdater.cc @@ -125,9 +125,9 @@ void PATJetUpdater::produce(edm::Event & iEvent, const edm::EventSetup & iSetup) // construct the Jet from the ref -> save ref to original object unsigned int idx = itJet - jets->begin(); - edm::RefToBase jetRef = jets->refAt(idx); - edm::Ptr jetPtr = jets->ptrAt(idx); - Jet ajet( edm::RefToBase(jetRef.castTo()) ); + const edm::RefToBase jetRef = jets->refAt(idx); + const edm::RefToBase patJetRef(jetRef.castTo()); + Jet ajet( patJetRef ); if (addJetCorrFactors_) { // undo previous jet energy corrections @@ -198,6 +198,10 @@ void PATJetUpdater::produce(edm::Event & iEvent, const edm::EventSetup & iSetup) userDataHelper_.add( ajet, iEvent, iSetup ); } + // reassign the original object reference to preserve reference to the original jet the input PAT jet was derived from + // (this needs to be done at the end since cloning the input PAT jet would interfere with adding UserData) + ajet.refToOrig_ = patJetRef->originalObjectRef(); + patJets->push_back(ajet); } diff --git a/PhysicsTools/PatAlgos/python/tools/jetTools.py b/PhysicsTools/PatAlgos/python/tools/jetTools.py index 7eecd779087f1..c0f17c7ef91b6 100644 --- a/PhysicsTools/PatAlgos/python/tools/jetTools.py +++ b/PhysicsTools/PatAlgos/python/tools/jetTools.py @@ -159,14 +159,15 @@ def setupJetCorrections(process, knownModules, jetCorrections, jetSource, pvSour setattr(process, 'patMETs'+labelName+postfix, patMETs.clone(metSource = cms.InputTag(jetCorrections[0]+_labelCorrName+'Type1p2CorMet'+postfix), addMuonCorrections = False)) -def setupSVClustering(btagInfo, algo, rParam, fatJets=cms.InputTag(''), groomedFatJets=cms.InputTag('')): - btagInfo.useSVClustering = cms.bool(True) - btagInfo.jetAlgorithm = cms.string(algo) - btagInfo.rParam = cms.double(rParam) - ## if the jets is actually a subjet - if fatJets != cms.InputTag('') and groomedFatJets != cms.InputTag(''): - btagInfo.fatJets = fatJets - btagInfo.groomedFatJets = groomedFatJets +def setupSVClustering(btagInfo, svClustering, algo, rParam, fatJets=cms.InputTag(''), groomedFatJets=cms.InputTag('')): + btagInfo.useSVClustering = cms.bool(svClustering) + btagInfo.jetAlgorithm = cms.string(algo) + btagInfo.rParam = cms.double(rParam) + ## if the jet is actually a subjet + if fatJets != cms.InputTag(''): + btagInfo.fatJets = fatJets + if groomedFatJets != cms.InputTag(''): + btagInfo.groomedFatJets = groomedFatJets def setupBTagging(process, jetSource, pfCandidates, explicitJTA, pvSource, svSource, elSource, muSource, runIVF, loadStdRecoBTag, svClustering, fatJets, groomedFatJets, @@ -225,16 +226,16 @@ def setupBTagging(process, jetSource, pfCandidates, explicitJTA, pvSource, svSou setattr(process, btagInfo+labelName+postfix, btag.pfSecondaryVertexTagInfos.clone(trackIPTagInfos = cms.InputTag('pfImpactParameterTagInfos'+labelName+postfix))) if btagInfo == 'pfInclusiveSecondaryVertexFinderTagInfos': setattr(process, btagInfo+labelName+postfix, btag.pfInclusiveSecondaryVertexFinderTagInfos.clone(trackIPTagInfos = cms.InputTag('pfImpactParameterTagInfos'+labelName+postfix), extSVCollection=svSource)) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'pfInclusiveSecondaryVertexFinderAK8TagInfos': setattr(process, btagInfo+labelName+postfix, btag.pfInclusiveSecondaryVertexFinderAK8TagInfos.clone(trackIPTagInfos = cms.InputTag('pfImpactParameterAK8TagInfos'+labelName+postfix), extSVCollection=svSource)) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'pfInclusiveSecondaryVertexFinderCA15TagInfos': setattr(process, btagInfo+labelName+postfix, btag.pfInclusiveSecondaryVertexFinderCA15TagInfos.clone(trackIPTagInfos = cms.InputTag('pfImpactParameterCA15TagInfos'+labelName+postfix), extSVCollection=svSource)) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'pfInclusiveSecondaryVertexFinderCvsLTagInfos': setattr( process, @@ -244,8 +245,8 @@ def setupBTagging(process, jetSource, pfCandidates, explicitJTA, pvSource, svSou extSVCollection=svSource ) ) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'pfInclusiveSecondaryVertexFinderCvsBTagInfos': setattr( process, @@ -255,8 +256,8 @@ def setupBTagging(process, jetSource, pfCandidates, explicitJTA, pvSource, svSou extSVCollection=svSource ) ) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'pfInclusiveSecondaryVertexFinderNegativeCvsLTagInfos': setattr( process, @@ -266,8 +267,8 @@ def setupBTagging(process, jetSource, pfCandidates, explicitJTA, pvSource, svSou extSVCollection=svSource ) ) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'pfGhostTrackVertexTagInfos': setattr(process, btagInfo+labelName+postfix, btag.pfGhostTrackVertexTagInfos.clone(trackIPTagInfos = cms.InputTag('pfImpactParameterTagInfos'+labelName+postfix))) if btagInfo == 'pfGhostTrackTagVertexInfosAK8': @@ -276,30 +277,30 @@ def setupBTagging(process, jetSource, pfCandidates, explicitJTA, pvSource, svSou setattr(process, btagInfo+labelName+postfix, btag.pfSecondaryVertexNegativeTagInfos.clone(trackIPTagInfos = cms.InputTag('pfImpactParameterTagInfos'+labelName+postfix))) if btagInfo == 'pfInclusiveSecondaryVertexFinderNegativeTagInfos': setattr(process, btagInfo+labelName+postfix, btag.pfInclusiveSecondaryVertexFinderNegativeTagInfos.clone(trackIPTagInfos = cms.InputTag('pfImpactParameterTagInfos'+labelName+postfix), extSVCollection=svSource)) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'impactParameterTagInfos': setattr(process, btagInfo+labelName+postfix, btag.impactParameterTagInfos.clone(jetTracks = cms.InputTag('jetTracksAssociatorAtVertex'+labelName+postfix), primaryVertex=pvSource)) if btagInfo == 'secondaryVertexTagInfos': setattr(process, btagInfo+labelName+postfix, btag.secondaryVertexTagInfos.clone(trackIPTagInfos = cms.InputTag('impactParameterTagInfos'+labelName+postfix))) if btagInfo == 'inclusiveSecondaryVertexFinderTagInfos': setattr(process, btagInfo+labelName+postfix, btag.inclusiveSecondaryVertexFinderTagInfos.clone(trackIPTagInfos = cms.InputTag('impactParameterTagInfos'+labelName+postfix))) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'inclusiveSecondaryVertexFinderFilteredTagInfos': setattr(process, btagInfo+labelName+postfix, btag.inclusiveSecondaryVertexFinderFilteredTagInfos.clone(trackIPTagInfos = cms.InputTag('impactParameterTagInfos'+labelName+postfix))) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'secondaryVertexNegativeTagInfos': setattr(process, btagInfo+labelName+postfix, btag.secondaryVertexNegativeTagInfos.clone(trackIPTagInfos = cms.InputTag('impactParameterTagInfos'+labelName+postfix))) if btagInfo == 'inclusiveSecondaryVertexFinderNegativeTagInfos': setattr(process, btagInfo+labelName+postfix, btag.inclusiveSecondaryVertexFinderNegativeTagInfos.clone(trackIPTagInfos = cms.InputTag('impactParameterTagInfos'+labelName+postfix))) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'inclusiveSecondaryVertexFinderFilteredNegativeTagInfos': setattr(process, btagInfo+labelName+postfix, btag.inclusiveSecondaryVertexFinderFilteredNegativeTagInfos.clone(trackIPTagInfos = cms.InputTag('impactParameterTagInfos'+labelName+postfix))) - if svClustering: - setupSVClustering(getattr(process, btagInfo+labelName+postfix), algo, rParam, fatJets, groomedFatJets) + if svClustering or fatJets != cms.InputTag(''): + setupSVClustering(getattr(process, btagInfo+labelName+postfix), svClustering, algo, rParam, fatJets, groomedFatJets) if btagInfo == 'softMuonTagInfos': setattr(process, btagInfo+labelName+postfix, btag.softMuonTagInfos.clone(jets = jetSource, primaryVertex=pvSource)) if btagInfo == 'softPFMuonsTagInfos': diff --git a/PhysicsTools/PatAlgos/test/patTuple_updateJets_fromMiniAOD_cfg.py b/PhysicsTools/PatAlgos/test/patTuple_updateJets_fromMiniAOD_cfg.py index 0a56c78b755dd..7245a9a0af0f9 100644 --- a/PhysicsTools/PatAlgos/test/patTuple_updateJets_fromMiniAOD_cfg.py +++ b/PhysicsTools/PatAlgos/test/patTuple_updateJets_fromMiniAOD_cfg.py @@ -27,7 +27,7 @@ ) process.updatedPatJets.userData.userFloats.src += ['oldJetMass'] -## An example where the jet correction is undone +## An example where the jet corrections are undone updateJetCollection( process, labelName = 'UndoneJEC', @@ -36,7 +36,7 @@ ) process.updatedPatJetsUndoneJEC.userData.userFloats.src = [] -## An example where the jet correction are reapplied +## An example where the jet corrections are reapplied updateJetCollection( process, labelName = 'ReappliedJEC', @@ -45,6 +45,22 @@ ) process.updatedPatJetsReappliedJEC.userData.userFloats.src = [] +## An example where the jet energy corrections are updated to the current GlobalTag +## and specified b-tag discriminators are rerun and added to SoftDrop subjets +updateJetCollection( + process, + labelName = 'SoftDropSubjets', + jetSource = cms.InputTag('slimmedJetsAK8PFCHSSoftDropPacked:SubJets'), + jetCorrections = ('AK4PFchs', cms.vstring(['L1FastJet', 'L2Relative', 'L3Absolute']), 'None'), + btagDiscriminators = ['pfCombinedSecondaryVertexV2BJetTags', 'pfCombinedInclusiveSecondaryVertexV2BJetTags'], + explicitJTA = True, # needed for subjet b tagging + svClustering = False, # needed for subjet b tagging (IMPORTANT: Needs to be set to False to disable ghost-association which does not work with slimmed jets) + fatJets = cms.InputTag('slimmedJetsAK8'), # needed for subjet b tagging + rParam = 0.8, # needed for subjet b tagging + algo = 'ak' # has to be defined but is not used with svClustering=False +) +process.updatedPatJetsSoftDropSubjets.userData.userFloats.src = [] + ## ------------------------------------------------------ # In addition you usually want to change the following # parameters: diff --git a/RecoBTag/SecondaryVertex/plugins/TemplatedSecondaryVertexProducer.cc b/RecoBTag/SecondaryVertex/plugins/TemplatedSecondaryVertexProducer.cc index 44c98d6578360..48e255a12bcd6 100644 --- a/RecoBTag/SecondaryVertex/plugins/TemplatedSecondaryVertexProducer.cc +++ b/RecoBTag/SecondaryVertex/plugins/TemplatedSecondaryVertexProducer.cc @@ -21,6 +21,7 @@ #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/Utilities/interface/Exception.h" +#include "DataFormats/PatCandidates/interface/Jet.h" #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/TrackReco/interface/TrackFwd.h" #include "DataFormats/VertexReco/interface/Vertex.h" @@ -127,6 +128,9 @@ class TemplatedSecondaryVertexProducer : public edm::stream::EDProducer<> { const edm::Handle >& groomedJets, const edm::Handle >& subjets, std::vector >& matchedIndices); + void matchSubjets(const edm::Handle >& fatJets, + const edm::Handle >& subjets, + std::vector >& matchedIndices); const reco::Jet * toJet(const reco::Jet & j) { return &j; } const reco::Jet * toJet(const IPTI & j) { return &(*(j.jet())); } @@ -164,6 +168,7 @@ class TemplatedSecondaryVertexProducer : public edm::stream::EDProducer<> { double ghostRescaling; double relPtTolerance; bool useFatJets; + bool useGroomedFatJets; edm::EDGetTokenT > token_fatJets; edm::EDGetTokenT > token_groomedFatJets; @@ -273,15 +278,16 @@ TemplatedSecondaryVertexProducer::TemplatedSecondaryVertexProducer( constraint == CONSTRAINT_BEAMSPOT || constraint == CONSTRAINT_PV_PRIMARIES_IN_FIT ) token_BeamSpot = consumes(params.getParameter("beamSpotTag")); - useExternalSV = false; - if(params.existsAs("useExternalSV")) useExternalSV = params.getParameter ("useExternalSV"); - if(useExternalSV) { + useExternalSV = false; + if(params.existsAs("useExternalSV")) useExternalSV = params.getParameter ("useExternalSV"); + if(useExternalSV) { token_extSVCollection = consumes >(params.getParameter("extSVCollection")); - extSVDeltaRToJet = params.getParameter("extSVDeltaRToJet"); - } - useSVClustering = ( params.existsAs("useSVClustering") ? params.getParameter("useSVClustering") : false ); + extSVDeltaRToJet = params.getParameter("extSVDeltaRToJet"); + } + useSVClustering = ( params.existsAs("useSVClustering") ? params.getParameter("useSVClustering") : false ); useSVMomentum = ( params.existsAs("useSVMomentum") ? params.getParameter("useSVMomentum") : false ); - useFatJets = ( useExternalSV && useSVClustering && params.exists("fatJets") && params.exists("groomedFatJets") ); + useFatJets = ( useExternalSV && params.exists("fatJets") ); + useGroomedFatJets = ( useExternalSV && params.exists("groomedFatJets") ); if( useSVClustering ) { jetAlgorithm = params.getParameter("jetAlgorithm"); @@ -301,10 +307,11 @@ TemplatedSecondaryVertexProducer::TemplatedSecondaryVertexProducer( throw cms::Exception("InvalidJetAlgorithm") << "Jet clustering algorithm is invalid: " << jetAlgorithm << ", use CambridgeAachen | Kt | AntiKt" << std::endl; } if( useFatJets ) - { token_fatJets = consumes >(params.getParameter("fatJets")); + if( useGroomedFatJets ) token_groomedFatJets = consumes >(params.getParameter("groomedFatJets")); - } + if( useFatJets && !useSVClustering ) + rParam = params.getParameter("rParam"); // will be used later as a dR cut produces(); } @@ -339,10 +346,14 @@ void TemplatedSecondaryVertexProducer::produce(edm::Event &event, if( useFatJets ) { event.getByToken(token_fatJets, fatJetsHandle); - event.getByToken(token_groomedFatJets, groomedFatJetsHandle); - if( groomedFatJetsHandle->size() > fatJetsHandle->size() ) - edm::LogError("TooManyGroomedJets") << "There are more groomed (" << groomedFatJetsHandle->size() << ") than original fat jets (" << fatJetsHandle->size() << "). Please check that the two jet collections belong to each other."; + if( useGroomedFatJets ) + { + event.getByToken(token_groomedFatJets, groomedFatJetsHandle); + + if( groomedFatJetsHandle->size() > fatJetsHandle->size() ) + edm::LogError("TooManyGroomedJets") << "There are more groomed (" << groomedFatJetsHandle->size() << ") than original fat jets (" << fatJetsHandle->size() << "). Please check that the two jet collections belong to each other."; + } } edm::Handle beamSpot; @@ -441,7 +452,7 @@ void TemplatedSecondaryVertexProducer::produce(edm::Event &event, if( useFatJets ) { if( inclusiveJets.size() < fatJetsHandle->size() ) - edm::LogError("TooFewReclusteredJets") << "There are fewer reclustered (" << inclusiveJets.size() << ") than original fat jets (" << fatJetsHandle->size() << "). Please check that the jet algorithm and jet size match those used for the original jet collection."; + edm::LogError("TooFewReclusteredJets") << "There are fewer reclustered (" << inclusiveJets.size() << ") than original fat jets (" << fatJetsHandle->size() << "). Please check that the jet algorithm and jet size match those used for the original jet collection."; // match reclustered and original fat jets std::vector reclusteredIndices; @@ -449,11 +460,15 @@ void TemplatedSecondaryVertexProducer::produce(edm::Event &event, // match groomed and original fat jets std::vector groomedIndices; - matchGroomedJets(fatJetsHandle,groomedFatJetsHandle,groomedIndices); + if( useGroomedFatJets ) + matchGroomedJets(fatJetsHandle,groomedFatJetsHandle,groomedIndices); // match subjets and original fat jets std::vector > subjetIndices; - matchSubjets(groomedIndices,groomedFatJetsHandle,trackIPTagInfos,subjetIndices); + if( useGroomedFatJets ) + matchSubjets(groomedIndices,groomedFatJetsHandle,trackIPTagInfos,subjetIndices); + else + matchSubjets(fatJetsHandle,trackIPTagInfos,subjetIndices); // collect clustered SVs for(size_t i=0; isize(); ++i) @@ -564,6 +579,65 @@ void TemplatedSecondaryVertexProducer::produce(edm::Event &event, } } } + // case where fat jets are used to associate SVs to subjets but no SV clustering is performed + else if( useExternalSV && !useSVClustering && trackIPTagInfos->size()>0 && useFatJets ) + { + // match groomed and original fat jets + std::vector groomedIndices; + if( useGroomedFatJets ) + matchGroomedJets(fatJetsHandle,groomedFatJetsHandle,groomedIndices); + + // match subjets and original fat jets + std::vector > subjetIndices; + if( useGroomedFatJets ) + matchSubjets(groomedIndices,groomedFatJetsHandle,trackIPTagInfos,subjetIndices); + else + matchSubjets(fatJetsHandle,trackIPTagInfos,subjetIndices); + + // loop over fat jets + for(size_t i=0; isize(); ++i) + { + if( fatJetsHandle->at(i).pt() == 0 ) // continue if the original jet has Pt=0 + { + edm::LogWarning("NullTransverseMomentum") << "The original fat jet " << i << " has Pt=0. This is not expected so the jet will be skipped."; + continue; + } + + if( subjetIndices.at(i).size()==0 ) continue; // continue if the original jet does not have subjets assigned + + // loop over SVs, associate them to fat jets based on dR cone and + // then assign them to the closets subjet in dR + for(typename edm::View::const_iterator it = extSecVertex->begin(); it != extSecVertex->end(); ++it) + { + size_t sv = ( it - extSecVertex->begin() ); + + const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex()); + const VTX &extSV = (*extSecVertex)[sv]; + GlobalVector dir = flightDirection(pv, extSV); + GlobalVector jetDir(fatJetsHandle->at(i).px(), + fatJetsHandle->at(i).py(), + fatJetsHandle->at(i).pz()); + // skip SVs outside the dR cone + if( Geom::deltaR2( dir, jetDir ) > rParam*rParam ) // here using the jet clustering rParam as a dR cut + continue; + + dir = dir.unit(); + fastjet::PseudoJet p(dir.x(),dir.y(),dir.z(),dir.mag()); // using SV flight direction so treating SV as massless + if( useSVMomentum ) + p = fastjet::PseudoJet(extSV.p4().px(),extSV.p4().py(),extSV.p4().pz(),extSV.p4().energy()); + + std::vector dR2toSubjets; + + for(size_t sj=0; sjat(subjetIndices.at(i).at(sj)).jet()->rapidity(), trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->phi() ) ); + + // find the closest subjet + int closestSubjetIdx = std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) ); + + clusteredSVs.at(subjetIndices.at(i).at(closestSubjetIdx)).push_back(sv); + } + } + } // ------------------------------------ SV clustering END ---------------------------------------------- std::auto_ptr vertexReco; @@ -773,16 +847,7 @@ void TemplatedSecondaryVertexProducer::produce(edm::Event &event, SVFilter(vertexFilter, pv, jetDir)); }else{ - if( !useSVClustering ) { - for(size_t iExtSv = 0; iExtSv < extSecVertex->size(); iExtSv++){ - const VTX & extVertex = (*extSecVertex)[iExtSv]; - if( Geom::deltaR2( ( position(extVertex) - pv.position() ), jetDir ) > extSVDeltaRToJet*extSVDeltaRToJet || extVertex.p4().M() < 0.3 ) - continue; - extAssoCollection.push_back( extVertex ); - } - - } - else { + if( useSVClustering || useFatJets ) { size_t jetIdx = ( iterJets - trackIPTagInfos->begin() ); for(size_t iExtSv = 0; iExtSv < clusteredSVs.at(jetIdx).size(); iExtSv++){ @@ -792,6 +857,14 @@ void TemplatedSecondaryVertexProducer::produce(edm::Event &event, extAssoCollection.push_back( extVertex ); } } + else { + for(size_t iExtSv = 0; iExtSv < extSecVertex->size(); iExtSv++){ + const VTX & extVertex = (*extSecVertex)[iExtSv]; + if( Geom::deltaR2( ( position(extVertex) - pv.position() ), jetDir ) > extSVDeltaRToJet*extSVDeltaRToJet || extVertex.p4().M() < 0.3 ) + continue; + extAssoCollection.push_back( extVertex ); + } + } // build combined SV information and filter SVBuilder svBuilder(pv, jetDir, withPVError, minTrackWeight); std::remove_copy_if(boost::make_transform_iterator( extAssoCollection.begin(), svBuilder), @@ -1075,6 +1148,74 @@ void TemplatedSecondaryVertexProducer::matchSubjets(const std::vector< } } +// ------------ method that matches subjets and original fat jets ------------ +template +void TemplatedSecondaryVertexProducer::matchSubjets(const edm::Handle >& fatJets, + const edm::Handle >& subjets, + std::vector >& matchedIndices) +{ + for(size_t fj=0; fjsize(); ++fj) + { + std::vector subjetIndices; + size_t nSubjetCollections = 0; + size_t nSubjets = 0; + + const pat::Jet * fatJet = dynamic_cast( fatJets->ptrAt(fj).get() ); + + if( !fatJet ) + { + if( fj==0 ) edm::LogError("WrongJetType") << "Wrong jet type for input fat jets. Please check that the input fat jets are of the pat::Jet type."; + + matchedIndices.push_back(subjetIndices); + continue; + } + else + { + nSubjetCollections = fatJet->subjetCollectionNames().size(); + + if( nSubjetCollections>0 ) + { + for(size_t coll=0; collsubjets(coll); + + for(size_t fjsj=0; fjsjsize(); ++sj) + { + const pat::Jet * subJet = dynamic_cast( subjets->at(sj).jet().get() ); + + if( !subJet ) + { + if( fj==0 && coll==0 && fjsj==0 && sj==0 ) edm::LogError("WrongJetType") << "Wrong jet type for input subjets. Please check that the input subjets are of the pat::Jet type."; + + break; + } + else + { + if( subJet->originalObjectRef() == fatJetSubjets.at(fjsj)->originalObjectRef() ) + { + subjetIndices.push_back(sj); + break; + } + } + } + } + } + + if( subjetIndices.size() == 0 && nSubjets > 0) + edm::LogError("SubjetMatchingFailed") << "Matching subjets to fat jets failed. Please check that the fat jet and subjet collections belong to each other."; + + matchedIndices.push_back(subjetIndices); + } + else + matchedIndices.push_back(subjetIndices); + } + } +} + // ------------ method fills 'descriptions' with the allowed parameters for the module ------------ template void TemplatedSecondaryVertexProducer::fillDescriptions(edm::ConfigurationDescriptions & descriptions) {