Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:34 

0001 #include <algorithm>
0002 #include <cstddef>
0003 #include <functional>
0004 #include <iterator>
0005 #include <map>
0006 #include <memory>
0007 #include <set>
0008 #include <string>
0009 #include <vector>
0010 
0011 #include <boost/iterator/transform_iterator.hpp>
0012 
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0020 #include "FWCore/ParameterSet/interface/IfExistsDescription.h"
0021 #include "FWCore/Utilities/interface/InputTag.h"
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/Utilities/interface/Exception.h"
0024 #include "FWCore/Utilities/interface/ESGetToken.h"
0025 
0026 #include "FWCore/ParameterSet/interface/Registry.h"
0027 #include "FWCore/Common/interface/Provenance.h"
0028 #include "DataFormats/Provenance/interface/ProductProvenance.h"
0029 
0030 #include "DataFormats/PatCandidates/interface/Jet.h"
0031 #include "DataFormats/TrackReco/interface/Track.h"
0032 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0033 #include "DataFormats/VertexReco/interface/Vertex.h"
0034 
0035 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0036 #include "DataFormats/BTauReco/interface/TrackIPTagInfo.h"
0037 #include "DataFormats/BTauReco/interface/CandIPTagInfo.h"
0038 #include "DataFormats/BTauReco/interface/SecondaryVertexTagInfo.h"
0039 
0040 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0041 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0042 #include "RecoVertex/ConfigurableVertexReco/interface/ConfigurableVertexReconstructor.h"
0043 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackVertexFinder.h"
0044 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackPrediction.h"
0045 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackState.h"
0046 #include "RecoVertex/GhostTrackFitter/interface/GhostTrack.h"
0047 
0048 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0049 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0050 #include "TrackingTools/TransientTrack/interface/CandidatePtrTransientTrack.h"
0051 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0052 
0053 #include "RecoBTag/SecondaryVertex/interface/TrackSelector.h"
0054 #include "RecoBTag/SecondaryVertex/interface/TrackSorting.h"
0055 #include "RecoBTag/SecondaryVertex/interface/SecondaryVertex.h"
0056 #include "RecoBTag/SecondaryVertex/interface/VertexFilter.h"
0057 #include "RecoBTag/SecondaryVertex/interface/VertexSorting.h"
0058 
0059 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0060 
0061 #include "fastjet/JetDefinition.hh"
0062 #include "fastjet/ClusterSequence.hh"
0063 #include "fastjet/PseudoJet.hh"
0064 
0065 //
0066 // constants, enums and typedefs
0067 //
0068 typedef std::shared_ptr<fastjet::ClusterSequence> ClusterSequencePtr;
0069 typedef std::shared_ptr<fastjet::JetDefinition> JetDefPtr;
0070 
0071 using namespace reco;
0072 
0073 namespace {
0074   class VertexInfo : public fastjet::PseudoJet::UserInfoBase {
0075   public:
0076     VertexInfo(const int vertexIndex) : m_vertexIndex(vertexIndex) {}
0077 
0078     inline const int vertexIndex() const { return m_vertexIndex; }
0079 
0080   protected:
0081     int m_vertexIndex;
0082   };
0083 
0084   template <typename T>
0085   struct RefToBaseLess {
0086     inline bool operator()(const edm::RefToBase<T> &r1, const edm::RefToBase<T> &r2) const {
0087       return r1.id() < r2.id() || (r1.id() == r2.id() && r1.key() < r2.key());
0088     }
0089   };
0090 }  // namespace
0091 
0092 GlobalVector flightDirection(const reco::Vertex &pv, const reco::Vertex &sv) {
0093   return GlobalVector(sv.x() - pv.x(), sv.y() - pv.y(), sv.z() - pv.z());
0094 }
0095 GlobalVector flightDirection(const reco::Vertex &pv, const reco::VertexCompositePtrCandidate &sv) {
0096   return GlobalVector(sv.vertex().x() - pv.x(), sv.vertex().y() - pv.y(), sv.vertex().z() - pv.z());
0097 }
0098 const math::XYZPoint &position(const reco::Vertex &sv) { return sv.position(); }
0099 const math::XYZPoint &position(const reco::VertexCompositePtrCandidate &sv) { return sv.vertex(); }
0100 
0101 template <class IPTI, class VTX>
0102 class TemplatedSecondaryVertexProducer : public edm::stream::EDProducer<> {
0103 public:
0104   explicit TemplatedSecondaryVertexProducer(const edm::ParameterSet &params);
0105   ~TemplatedSecondaryVertexProducer() override;
0106   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0107   typedef std::vector<TemplatedSecondaryVertexTagInfo<IPTI, VTX> > Product;
0108   typedef TemplatedSecondaryVertex<VTX> SecondaryVertex;
0109   typedef typename IPTI::input_container input_container;
0110   typedef typename IPTI::input_container::value_type input_item;
0111   typedef typename std::vector<reco::btag::IndexedTrackData> TrackDataVector;
0112   void produce(edm::Event &event, const edm::EventSetup &es) override;
0113 
0114 private:
0115   template <class CONTAINER>
0116   void matchReclusteredJets(const edm::Handle<CONTAINER> &jets,
0117                             const std::vector<fastjet::PseudoJet> &matchedJets,
0118                             std::vector<int> &matchedIndices,
0119                             const std::string &jetType = "");
0120   void matchGroomedJets(const edm::Handle<edm::View<reco::Jet> > &jets,
0121                         const edm::Handle<edm::View<reco::Jet> > &matchedJets,
0122                         std::vector<int> &matchedIndices);
0123   void matchSubjets(const std::vector<int> &groomedIndices,
0124                     const edm::Handle<edm::View<reco::Jet> > &groomedJets,
0125                     const edm::Handle<std::vector<IPTI> > &subjets,
0126                     std::vector<std::vector<int> > &matchedIndices);
0127   void matchSubjets(const edm::Handle<edm::View<reco::Jet> > &fatJets,
0128                     const edm::Handle<std::vector<IPTI> > &subjets,
0129                     std::vector<std::vector<int> > &matchedIndices);
0130 
0131   const reco::Jet *toJet(const reco::Jet &j) { return &j; }
0132   const reco::Jet *toJet(const IPTI &j) { return &(*(j.jet())); }
0133 
0134   enum ConstraintType {
0135     CONSTRAINT_NONE = 0,
0136     CONSTRAINT_BEAMSPOT,
0137     CONSTRAINT_PV_BEAMSPOT_SIZE,
0138     CONSTRAINT_PV_BS_Z_ERRORS_SCALED,
0139     CONSTRAINT_PV_ERROR_SCALED,
0140     CONSTRAINT_PV_PRIMARIES_IN_FIT
0141   };
0142   static ConstraintType getConstraintType(const std::string &name);
0143 
0144   edm::EDGetTokenT<reco::BeamSpot> token_BeamSpot;
0145   edm::EDGetTokenT<std::vector<IPTI> > token_trackIPTagInfo;
0146   reco::btag::SortCriteria sortCriterium;
0147   TrackSelector trackSelector;
0148   ConstraintType constraint;
0149   double constraintScaling;
0150   edm::ParameterSet vtxRecoPSet;
0151   bool useGhostTrack;
0152   bool withPVError;
0153   double minTrackWeight;
0154   VertexFilter vertexFilter;
0155   VertexSorting<SecondaryVertex> vertexSorting;
0156   bool useExternalSV;
0157   double extSVDeltaRToJet;
0158   edm::EDGetTokenT<edm::View<VTX> > token_extSVCollection;
0159   bool useSVClustering;
0160   bool useSVMomentum;
0161   std::string jetAlgorithm;
0162   double rParam;
0163   double jetPtMin;
0164   double ghostRescaling;
0165   double relPtTolerance;
0166   bool useFatJets;
0167   bool useGroomedFatJets;
0168   edm::EDGetTokenT<edm::View<reco::Jet> > token_fatJets;
0169   edm::EDGetTokenT<edm::View<reco::Jet> > token_groomedFatJets;
0170   edm::EDGetTokenT<edm::ValueMap<float> > token_weights;
0171   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> token_trackBuilder;
0172 
0173   ClusterSequencePtr fjClusterSeq;
0174   JetDefPtr fjJetDefinition;
0175 
0176   void markUsedTracks(TrackDataVector &trackData,
0177                       const input_container &trackRefs,
0178                       const SecondaryVertex &sv,
0179                       size_t idx);
0180 
0181   struct SVBuilder {
0182     SVBuilder(const reco::Vertex &pv, const GlobalVector &direction, bool withPVError, double minTrackWeight)
0183         : pv(pv), direction(direction), withPVError(withPVError), minTrackWeight(minTrackWeight) {}
0184     SecondaryVertex operator()(const TransientVertex &sv) const;
0185 
0186     SecondaryVertex operator()(const VTX &sv) const { return SecondaryVertex(pv, sv, direction, withPVError); }
0187 
0188     const Vertex &pv;
0189     const GlobalVector &direction;
0190     bool withPVError;
0191     double minTrackWeight;
0192   };
0193 
0194   struct SVFilter {
0195     SVFilter(const VertexFilter &filter, const Vertex &pv, const GlobalVector &direction)
0196         : filter(filter), pv(pv), direction(direction) {}
0197 
0198     inline bool operator()(const SecondaryVertex &sv) const { return !filter(pv, sv, direction); }
0199 
0200     const VertexFilter &filter;
0201     const Vertex &pv;
0202     const GlobalVector &direction;
0203   };
0204 };
0205 template <class IPTI, class VTX>
0206 typename TemplatedSecondaryVertexProducer<IPTI, VTX>::ConstraintType
0207 TemplatedSecondaryVertexProducer<IPTI, VTX>::getConstraintType(const std::string &name) {
0208   if (name == "None")
0209     return CONSTRAINT_NONE;
0210   else if (name == "BeamSpot")
0211     return CONSTRAINT_BEAMSPOT;
0212   else if (name == "BeamSpot+PVPosition")
0213     return CONSTRAINT_PV_BEAMSPOT_SIZE;
0214   else if (name == "BeamSpotZ+PVErrorScaledXY")
0215     return CONSTRAINT_PV_BS_Z_ERRORS_SCALED;
0216   else if (name == "PVErrorScaled")
0217     return CONSTRAINT_PV_ERROR_SCALED;
0218   else if (name == "BeamSpot+PVTracksInFit")
0219     return CONSTRAINT_PV_PRIMARIES_IN_FIT;
0220   else
0221     throw cms::Exception("InvalidArgument") << "TemplatedSecondaryVertexProducer: ``constraint'' parameter "
0222                                                "value \""
0223                                             << name << "\" not understood." << std::endl;
0224 }
0225 
0226 static GhostTrackVertexFinder::FitType getGhostTrackFitType(const std::string &name) {
0227   if (name == "AlwaysWithGhostTrack")
0228     return GhostTrackVertexFinder::kAlwaysWithGhostTrack;
0229   else if (name == "SingleTracksWithGhostTrack")
0230     return GhostTrackVertexFinder::kSingleTracksWithGhostTrack;
0231   else if (name == "RefitGhostTrackWithVertices")
0232     return GhostTrackVertexFinder::kRefitGhostTrackWithVertices;
0233   else
0234     throw cms::Exception("InvalidArgument") << "TemplatedSecondaryVertexProducer: ``fitType'' "
0235                                                "parameter value \""
0236                                             << name
0237                                             << "\" for "
0238                                                "GhostTrackVertexFinder settings not "
0239                                                "understood."
0240                                             << std::endl;
0241 }
0242 
0243 template <class IPTI, class VTX>
0244 TemplatedSecondaryVertexProducer<IPTI, VTX>::TemplatedSecondaryVertexProducer(const edm::ParameterSet &params)
0245     : sortCriterium(TrackSorting::getCriterium(params.getParameter<std::string>("trackSort"))),
0246       trackSelector(params.getParameter<edm::ParameterSet>("trackSelection")),
0247       constraint(getConstraintType(params.getParameter<std::string>("constraint"))),
0248       constraintScaling(1.0),
0249       vtxRecoPSet(params.getParameter<edm::ParameterSet>("vertexReco")),
0250       useGhostTrack(vtxRecoPSet.getParameter<std::string>("finder") == "gtvr"),
0251       withPVError(params.getParameter<bool>("usePVError")),
0252       minTrackWeight(params.getParameter<double>("minimumTrackWeight")),
0253       vertexFilter(params.getParameter<edm::ParameterSet>("vertexCuts")),
0254       vertexSorting(params.getParameter<edm::ParameterSet>("vertexSelection")) {
0255   token_trackIPTagInfo = consumes<std::vector<IPTI> >(params.getParameter<edm::InputTag>("trackIPTagInfos"));
0256   if (constraint == CONSTRAINT_PV_ERROR_SCALED || constraint == CONSTRAINT_PV_BS_Z_ERRORS_SCALED)
0257     constraintScaling = params.getParameter<double>("pvErrorScaling");
0258 
0259   if (constraint == CONSTRAINT_PV_BEAMSPOT_SIZE || constraint == CONSTRAINT_PV_BS_Z_ERRORS_SCALED ||
0260       constraint == CONSTRAINT_BEAMSPOT || constraint == CONSTRAINT_PV_PRIMARIES_IN_FIT)
0261     token_BeamSpot = consumes<reco::BeamSpot>(params.getParameter<edm::InputTag>("beamSpotTag"));
0262   useExternalSV = params.getParameter<bool>("useExternalSV");
0263   if (useExternalSV) {
0264     token_extSVCollection = consumes<edm::View<VTX> >(params.getParameter<edm::InputTag>("extSVCollection"));
0265     extSVDeltaRToJet = params.getParameter<double>("extSVDeltaRToJet");
0266   }
0267   useSVClustering = (params.existsAs<bool>("useSVClustering") ? params.getParameter<bool>("useSVClustering") : false);
0268   useSVMomentum = (params.existsAs<bool>("useSVMomentum") ? params.getParameter<bool>("useSVMomentum") : false);
0269   useFatJets = (useExternalSV && params.exists("fatJets"));
0270   useGroomedFatJets = (useExternalSV && params.exists("groomedFatJets"));
0271   if (useSVClustering) {
0272     jetAlgorithm = params.getParameter<std::string>("jetAlgorithm");
0273     rParam = params.getParameter<double>("rParam");
0274     jetPtMin =
0275         0.;  // hardcoded to 0. since we simply want to recluster all input jets which already had some PtMin applied
0276     ghostRescaling =
0277         (params.existsAs<double>("ghostRescaling") ? params.getParameter<double>("ghostRescaling") : 1e-18);
0278     relPtTolerance =
0279         (params.existsAs<double>("relPtTolerance")
0280              ? params.getParameter<double>("relPtTolerance")
0281              : 1e-03);  // 0.1% relative difference in Pt should be sufficient to detect possible misconfigurations
0282 
0283     // set jet algorithm
0284     if (jetAlgorithm == "Kt")
0285       fjJetDefinition = std::make_shared<fastjet::JetDefinition>(fastjet::kt_algorithm, rParam);
0286     else if (jetAlgorithm == "CambridgeAachen")
0287       fjJetDefinition = std::make_shared<fastjet::JetDefinition>(fastjet::cambridge_algorithm, rParam);
0288     else if (jetAlgorithm == "AntiKt")
0289       fjJetDefinition = std::make_shared<fastjet::JetDefinition>(fastjet::antikt_algorithm, rParam);
0290     else
0291       throw cms::Exception("InvalidJetAlgorithm") << "Jet clustering algorithm is invalid: " << jetAlgorithm
0292                                                   << ", use CambridgeAachen | Kt | AntiKt" << std::endl;
0293   }
0294   token_trackBuilder =
0295       esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"));
0296   if (useFatJets) {
0297     token_fatJets = consumes<edm::View<reco::Jet> >(params.getParameter<edm::InputTag>("fatJets"));
0298   }
0299   edm::InputTag srcWeights = params.getParameter<edm::InputTag>("weights");
0300   if (!srcWeights.label().empty())
0301     token_weights = consumes<edm::ValueMap<float> >(srcWeights);
0302   if (useGroomedFatJets) {
0303     token_groomedFatJets = consumes<edm::View<reco::Jet> >(params.getParameter<edm::InputTag>("groomedFatJets"));
0304   }
0305   if (useFatJets && !useSVClustering)
0306     rParam = params.getParameter<double>("rParam");  // will be used later as a dR cut
0307 
0308   produces<Product>();
0309 }
0310 template <class IPTI, class VTX>
0311 TemplatedSecondaryVertexProducer<IPTI, VTX>::~TemplatedSecondaryVertexProducer() {}
0312 
0313 template <class IPTI, class VTX>
0314 void TemplatedSecondaryVertexProducer<IPTI, VTX>::produce(edm::Event &event, const edm::EventSetup &es) {
0315   //    typedef std::map<TrackBaseRef, TransientTrack,
0316   //                     RefToBaseLess<Track> > TransientTrackMap;
0317   //How about good old pointers?
0318   typedef std::map<const Track *, TransientTrack> TransientTrackMap;
0319 
0320   edm::ESHandle<TransientTrackBuilder> trackBuilder = es.getHandle(token_trackBuilder);
0321 
0322   edm::Handle<std::vector<IPTI> > trackIPTagInfos;
0323   event.getByToken(token_trackIPTagInfo, trackIPTagInfos);
0324 
0325   // External Sec Vertex collection (e.g. for IVF usage)
0326   edm::Handle<edm::View<VTX> > extSecVertex;
0327   if (useExternalSV)
0328     event.getByToken(token_extSVCollection, extSecVertex);
0329 
0330   edm::Handle<edm::View<reco::Jet> > fatJetsHandle;
0331   edm::Handle<edm::View<reco::Jet> > groomedFatJetsHandle;
0332   if (useFatJets) {
0333     event.getByToken(token_fatJets, fatJetsHandle);
0334     if (useGroomedFatJets) {
0335       event.getByToken(token_groomedFatJets, groomedFatJetsHandle);
0336 
0337       if (groomedFatJetsHandle->size() > fatJetsHandle->size())
0338         edm::LogError("TooManyGroomedJets")
0339             << "There are more groomed (" << groomedFatJetsHandle->size() << ") than original fat jets ("
0340             << fatJetsHandle->size() << "). Please check that the two jet collections belong to each other.";
0341     }
0342   }
0343   edm::Handle<edm::ValueMap<float> > weightsHandle;
0344   if (!token_weights.isUninitialized())
0345     event.getByToken(token_weights, weightsHandle);
0346 
0347   edm::Handle<BeamSpot> beamSpot;
0348   unsigned int bsCovSrc[7] = {
0349       0,
0350   };
0351   double sigmaZ = 0.0, beamWidth = 0.0;
0352   switch (constraint) {
0353     case CONSTRAINT_PV_BEAMSPOT_SIZE:
0354       event.getByToken(token_BeamSpot, beamSpot);
0355       bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = bsCovSrc[6] = 1;
0356       sigmaZ = beamSpot->sigmaZ();
0357       beamWidth = beamSpot->BeamWidthX();
0358       break;
0359 
0360     case CONSTRAINT_PV_BS_Z_ERRORS_SCALED:
0361       event.getByToken(token_BeamSpot, beamSpot);
0362       bsCovSrc[0] = bsCovSrc[1] = 2;
0363       bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = 1;
0364       sigmaZ = beamSpot->sigmaZ();
0365       break;
0366 
0367     case CONSTRAINT_PV_ERROR_SCALED:
0368       bsCovSrc[0] = bsCovSrc[1] = bsCovSrc[2] = 2;
0369       break;
0370 
0371     case CONSTRAINT_BEAMSPOT:
0372     case CONSTRAINT_PV_PRIMARIES_IN_FIT:
0373       event.getByToken(token_BeamSpot, beamSpot);
0374       break;
0375 
0376     default:
0377         /* nothing */;
0378   }
0379 
0380   // ------------------------------------ SV clustering START --------------------------------------------
0381   std::vector<std::vector<int> > clusteredSVs(trackIPTagInfos->size(), std::vector<int>());
0382   if (useExternalSV && useSVClustering && !trackIPTagInfos->empty()) {
0383     // vector of constituents for reclustering jets and "ghost" SVs
0384     std::vector<fastjet::PseudoJet> fjInputs;
0385     // loop over all input jets and collect all their constituents
0386     if (useFatJets) {
0387       for (edm::View<reco::Jet>::const_iterator it = fatJetsHandle->begin(); it != fatJetsHandle->end(); ++it) {
0388         std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
0389         std::vector<edm::Ptr<reco::Candidate> >::const_iterator m;
0390         for (m = constituents.begin(); m != constituents.end(); ++m) {
0391           const reco::CandidatePtr &constit = *m;
0392           if (constit.isNull() || constit->pt() <= std::numeric_limits<double>::epsilon()) {
0393             edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
0394             continue;
0395           }
0396           if (it->isWeighted()) {
0397             if (token_weights.isUninitialized())
0398               throw cms::Exception("MissingConstituentWeight")
0399                   << "TemplatedSecondaryVertexProducer: No weights (e.g. PUPPI) given for weighted jet collection"
0400                   << std::endl;
0401             float w = (*weightsHandle)[constit];
0402             if (w > 0) {
0403               fjInputs.push_back(
0404                   fastjet::PseudoJet(constit->px() * w, constit->py() * w, constit->pz() * w, constit->energy() * w));
0405             }
0406           } else {
0407             fjInputs.push_back(fastjet::PseudoJet(constit->px(), constit->py(), constit->pz(), constit->energy()));
0408           }
0409         }
0410       }
0411     } else {
0412       for (typename std::vector<IPTI>::const_iterator it = trackIPTagInfos->begin(); it != trackIPTagInfos->end();
0413            ++it) {
0414         std::vector<edm::Ptr<reco::Candidate> > constituents = it->jet()->getJetConstituents();
0415         std::vector<edm::Ptr<reco::Candidate> >::const_iterator m;
0416         for (m = constituents.begin(); m != constituents.end(); ++m) {
0417           const reco::CandidatePtr &constit = *m;
0418           if (constit.isNull() || constit->pt() <= std::numeric_limits<double>::epsilon()) {
0419             edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
0420             continue;
0421           }
0422           if (it->jet()->isWeighted()) {
0423             if (token_weights.isUninitialized())
0424               throw cms::Exception("MissingConstituentWeight")
0425                   << "TemplatedSecondaryVertexProducer: No weights (e.g. PUPPI) given for weighted jet collection"
0426                   << std::endl;
0427             float w = (*weightsHandle)[constit];
0428             if (w > 0) {
0429               fjInputs.push_back(
0430                   fastjet::PseudoJet(constit->px() * w, constit->py() * w, constit->pz() * w, constit->energy() * w));
0431             }
0432           } else {
0433             fjInputs.push_back(fastjet::PseudoJet(constit->px(), constit->py(), constit->pz(), constit->energy()));
0434           }
0435         }
0436       }
0437     }
0438     // insert "ghost" SVs in the vector of constituents
0439     for (typename edm::View<VTX>::const_iterator it = extSecVertex->begin(); it != extSecVertex->end(); ++it) {
0440       const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
0441       GlobalVector dir = flightDirection(pv, *it);
0442       dir = dir.unit();
0443       fastjet::PseudoJet p(
0444           dir.x(), dir.y(), dir.z(), dir.mag());  // using SV flight direction so treating SV as massless
0445       if (useSVMomentum)
0446         p = fastjet::PseudoJet(it->p4().px(), it->p4().py(), it->p4().pz(), it->p4().energy());
0447       p *= ghostRescaling;  // rescale SV direction/momentum
0448       p.set_user_info(new VertexInfo(it - extSecVertex->begin()));
0449       fjInputs.push_back(p);
0450     }
0451 
0452     // define jet clustering sequence
0453     fjClusterSeq = std::make_shared<fastjet::ClusterSequence>(fjInputs, *fjJetDefinition);
0454     // recluster jet constituents and inserted "ghosts"
0455     std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(fjClusterSeq->inclusive_jets(jetPtMin));
0456 
0457     if (useFatJets) {
0458       if (inclusiveJets.size() < fatJetsHandle->size())
0459         edm::LogError("TooFewReclusteredJets")
0460             << "There are fewer reclustered (" << inclusiveJets.size() << ") than original fat jets ("
0461             << fatJetsHandle->size()
0462             << "). Please check that the jet algorithm and jet size match those used for the original jet collection.";
0463 
0464       // match reclustered and original fat jets
0465       std::vector<int> reclusteredIndices;
0466       matchReclusteredJets<edm::View<reco::Jet> >(fatJetsHandle, inclusiveJets, reclusteredIndices, "fat");
0467 
0468       // match groomed and original fat jets
0469       std::vector<int> groomedIndices;
0470       if (useGroomedFatJets)
0471         matchGroomedJets(fatJetsHandle, groomedFatJetsHandle, groomedIndices);
0472 
0473       // match subjets and original fat jets
0474       std::vector<std::vector<int> > subjetIndices;
0475       if (useGroomedFatJets)
0476         matchSubjets(groomedIndices, groomedFatJetsHandle, trackIPTagInfos, subjetIndices);
0477       else
0478         matchSubjets(fatJetsHandle, trackIPTagInfos, subjetIndices);
0479 
0480       // collect clustered SVs
0481       for (size_t i = 0; i < fatJetsHandle->size(); ++i) {
0482         if (reclusteredIndices.at(i) < 0)
0483           continue;  // continue if matching reclustered to original jets failed
0484 
0485         if (fatJetsHandle->at(i).pt() == 0)  // continue if the original jet has Pt=0
0486         {
0487           edm::LogWarning("NullTransverseMomentum")
0488               << "The original fat jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
0489           continue;
0490         }
0491 
0492         if (subjetIndices.at(i).empty())
0493           continue;  // continue if the original jet does not have subjets assigned
0494 
0495         // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original fat jets should in principle stay the same
0496         if ((std::abs(inclusiveJets.at(reclusteredIndices.at(i)).pt() - fatJetsHandle->at(i).pt()) /
0497              fatJetsHandle->at(i).pt()) > relPtTolerance) {
0498           if (fatJetsHandle->at(i).pt() < 10.)  // special handling for low-Pt jets (Pt<10 GeV)
0499             edm::LogWarning("JetPtMismatchAtLowPt")
0500                 << "The reclustered and original fat jet " << i << " have different Pt's ("
0501                 << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << fatJetsHandle->at(i).pt()
0502                 << " GeV, respectively).\n"
0503                 << "Please check that the jet algorithm and jet size match those used for the original fat jet "
0504                    "collection and also make sure the original fat jets are uncorrected. In addition, make sure you "
0505                    "are not using CaloJets which are presently not supported.\n"
0506                 << "Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n"
0507                 << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
0508                    "precision in which case make sure the original jet collection is produced and reclustering is "
0509                    "performed in the same job.";
0510           else
0511             edm::LogError("JetPtMismatch")
0512                 << "The reclustered and original fat jet " << i << " have different Pt's ("
0513                 << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << fatJetsHandle->at(i).pt()
0514                 << " GeV, respectively).\n"
0515                 << "Please check that the jet algorithm and jet size match those used for the original fat jet "
0516                    "collection and also make sure the original fat jets are uncorrected. In addition, make sure you "
0517                    "are not using CaloJets which are presently not supported.\n"
0518                 << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
0519                    "precision in which case make sure the original jet collection is produced and reclustering is "
0520                    "performed in the same job.";
0521         }
0522 
0523         // get jet constituents
0524         std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(i)).constituents();
0525 
0526         std::vector<int> svIndices;
0527         // loop over jet constituents and try to find "ghosts"
0528         for (std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end();
0529              ++it) {
0530           if (!it->has_user_info())
0531             continue;  // skip if not a "ghost"
0532 
0533           svIndices.push_back(it->user_info<VertexInfo>().vertexIndex());
0534         }
0535 
0536         // loop over clustered SVs and assign them to different subjets based on smallest dR
0537         for (size_t sv = 0; sv < svIndices.size(); ++sv) {
0538           const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
0539           const VTX &extSV = (*extSecVertex)[svIndices.at(sv)];
0540           GlobalVector dir = flightDirection(pv, extSV);
0541           dir = dir.unit();
0542           fastjet::PseudoJet p(
0543               dir.x(), dir.y(), dir.z(), dir.mag());  // using SV flight direction so treating SV as massless
0544           if (useSVMomentum)
0545             p = fastjet::PseudoJet(extSV.p4().px(), extSV.p4().py(), extSV.p4().pz(), extSV.p4().energy());
0546 
0547           std::vector<double> dR2toSubjets;
0548 
0549           for (size_t sj = 0; sj < subjetIndices.at(i).size(); ++sj)
0550             dR2toSubjets.push_back(Geom::deltaR2(p.rapidity(),
0551                                                  p.phi_std(),
0552                                                  trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->rapidity(),
0553                                                  trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->phi()));
0554 
0555           // find the closest subjet
0556           int closestSubjetIdx =
0557               std::distance(dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()));
0558 
0559           clusteredSVs.at(subjetIndices.at(i).at(closestSubjetIdx)).push_back(svIndices.at(sv));
0560         }
0561       }
0562     } else {
0563       if (inclusiveJets.size() < trackIPTagInfos->size())
0564         edm::LogError("TooFewReclusteredJets")
0565             << "There are fewer reclustered (" << inclusiveJets.size() << ") than original jets ("
0566             << trackIPTagInfos->size()
0567             << "). Please check that the jet algorithm and jet size match those used for the original jet collection.";
0568 
0569       // match reclustered and original jets
0570       std::vector<int> reclusteredIndices;
0571       matchReclusteredJets<std::vector<IPTI> >(trackIPTagInfos, inclusiveJets, reclusteredIndices);
0572 
0573       // collect clustered SVs
0574       for (size_t i = 0; i < trackIPTagInfos->size(); ++i) {
0575         if (reclusteredIndices.at(i) < 0)
0576           continue;  // continue if matching reclustered to original jets failed
0577 
0578         if (trackIPTagInfos->at(i).jet()->pt() == 0)  // continue if the original jet has Pt=0
0579         {
0580           edm::LogWarning("NullTransverseMomentum")
0581               << "The original jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
0582           continue;
0583         }
0584 
0585         // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original jets should in principle stay the same
0586         if ((std::abs(inclusiveJets.at(reclusteredIndices.at(i)).pt() - trackIPTagInfos->at(i).jet()->pt()) /
0587              trackIPTagInfos->at(i).jet()->pt()) > relPtTolerance) {
0588           if (trackIPTagInfos->at(i).jet()->pt() < 10.)  // special handling for low-Pt jets (Pt<10 GeV)
0589             edm::LogWarning("JetPtMismatchAtLowPt")
0590                 << "The reclustered and original jet " << i << " have different Pt's ("
0591                 << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << trackIPTagInfos->at(i).jet()->pt()
0592                 << " GeV, respectively).\n"
0593                 << "Please check that the jet algorithm and jet size match those used for the original jet collection "
0594                    "and also make sure the original jets are uncorrected. In addition, make sure you are not using "
0595                    "CaloJets which are presently not supported.\n"
0596                 << "Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n"
0597                 << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
0598                    "precision in which case make sure the original jet collection is produced and reclustering is "
0599                    "performed in the same job.";
0600           else
0601             edm::LogError("JetPtMismatch")
0602                 << "The reclustered and original jet " << i << " have different Pt's ("
0603                 << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << trackIPTagInfos->at(i).jet()->pt()
0604                 << " GeV, respectively).\n"
0605                 << "Please check that the jet algorithm and jet size match those used for the original jet collection "
0606                    "and also make sure the original jets are uncorrected. In addition, make sure you are not using "
0607                    "CaloJets which are presently not supported.\n"
0608                 << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
0609                    "precision in which case make sure the original jet collection is produced and reclustering is "
0610                    "performed in the same job.";
0611         }
0612 
0613         // get jet constituents
0614         std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(i)).constituents();
0615 
0616         // loop over jet constituents and try to find "ghosts"
0617         for (std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end();
0618              ++it) {
0619           if (!it->has_user_info())
0620             continue;  // skip if not a "ghost"
0621           // push back clustered SV indices
0622           clusteredSVs.at(i).push_back(it->user_info<VertexInfo>().vertexIndex());
0623         }
0624       }
0625     }
0626   }
0627   // case where fat jets are used to associate SVs to subjets but no SV clustering is performed
0628   else if (useExternalSV && !useSVClustering && !trackIPTagInfos->empty() && useFatJets) {
0629     // match groomed and original fat jets
0630     std::vector<int> groomedIndices;
0631     if (useGroomedFatJets)
0632       matchGroomedJets(fatJetsHandle, groomedFatJetsHandle, groomedIndices);
0633 
0634     // match subjets and original fat jets
0635     std::vector<std::vector<int> > subjetIndices;
0636     if (useGroomedFatJets)
0637       matchSubjets(groomedIndices, groomedFatJetsHandle, trackIPTagInfos, subjetIndices);
0638     else
0639       matchSubjets(fatJetsHandle, trackIPTagInfos, subjetIndices);
0640 
0641     // loop over fat jets
0642     for (size_t i = 0; i < fatJetsHandle->size(); ++i) {
0643       if (fatJetsHandle->at(i).pt() == 0)  // continue if the original jet has Pt=0
0644       {
0645         edm::LogWarning("NullTransverseMomentum")
0646             << "The original fat jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
0647         continue;
0648       }
0649 
0650       if (subjetIndices.at(i).empty())
0651         continue;  // continue if the original jet does not have subjets assigned
0652 
0653       // loop over SVs, associate them to fat jets based on dR cone and
0654       // then assign them to the closets subjet in dR
0655       for (typename edm::View<VTX>::const_iterator it = extSecVertex->begin(); it != extSecVertex->end(); ++it) {
0656         size_t sv = (it - extSecVertex->begin());
0657 
0658         const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
0659         const VTX &extSV = (*extSecVertex)[sv];
0660         GlobalVector dir = flightDirection(pv, extSV);
0661         GlobalVector jetDir(fatJetsHandle->at(i).px(), fatJetsHandle->at(i).py(), fatJetsHandle->at(i).pz());
0662         // skip SVs outside the dR cone
0663         if (Geom::deltaR2(dir, jetDir) > rParam * rParam)  // here using the jet clustering rParam as a dR cut
0664           continue;
0665 
0666         dir = dir.unit();
0667         fastjet::PseudoJet p(
0668             dir.x(), dir.y(), dir.z(), dir.mag());  // using SV flight direction so treating SV as massless
0669         if (useSVMomentum)
0670           p = fastjet::PseudoJet(extSV.p4().px(), extSV.p4().py(), extSV.p4().pz(), extSV.p4().energy());
0671 
0672         std::vector<double> dR2toSubjets;
0673 
0674         for (size_t sj = 0; sj < subjetIndices.at(i).size(); ++sj)
0675           dR2toSubjets.push_back(Geom::deltaR2(p.rapidity(),
0676                                                p.phi_std(),
0677                                                trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->rapidity(),
0678                                                trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->phi()));
0679 
0680         // find the closest subjet
0681         int closestSubjetIdx =
0682             std::distance(dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()));
0683 
0684         clusteredSVs.at(subjetIndices.at(i).at(closestSubjetIdx)).push_back(sv);
0685       }
0686     }
0687   }
0688   // ------------------------------------ SV clustering END ----------------------------------------------
0689 
0690   std::unique_ptr<ConfigurableVertexReconstructor> vertexReco;
0691   std::unique_ptr<GhostTrackVertexFinder> vertexRecoGT;
0692   if (useGhostTrack)
0693     vertexRecoGT = std::make_unique<GhostTrackVertexFinder>(
0694         vtxRecoPSet.getParameter<double>("maxFitChi2"),
0695         vtxRecoPSet.getParameter<double>("mergeThreshold"),
0696         vtxRecoPSet.getParameter<double>("primcut"),
0697         vtxRecoPSet.getParameter<double>("seccut"),
0698         getGhostTrackFitType(vtxRecoPSet.getParameter<std::string>("fitType")));
0699   else
0700     vertexReco = std::make_unique<ConfigurableVertexReconstructor>(vtxRecoPSet);
0701 
0702   TransientTrackMap primariesMap;
0703 
0704   // result secondary vertices
0705 
0706   auto tagInfos = std::make_unique<Product>();
0707 
0708   for (typename std::vector<IPTI>::const_iterator iterJets = trackIPTagInfos->begin();
0709        iterJets != trackIPTagInfos->end();
0710        ++iterJets) {
0711     TrackDataVector trackData;
0712     //            std::cout << "Jet " << iterJets-trackIPTagInfos->begin() << std::endl;
0713 
0714     const Vertex &pv = *iterJets->primaryVertex();
0715 
0716     std::set<TransientTrack> primaries;
0717     if (constraint == CONSTRAINT_PV_PRIMARIES_IN_FIT) {
0718       for (Vertex::trackRef_iterator iter = pv.tracks_begin(); iter != pv.tracks_end(); ++iter) {
0719         TransientTrackMap::iterator pos = primariesMap.lower_bound(iter->get());
0720 
0721         if (pos != primariesMap.end() && pos->first == iter->get())
0722           primaries.insert(pos->second);
0723         else {
0724           TransientTrack track = trackBuilder->build(iter->castTo<TrackRef>());
0725           primariesMap.insert(pos, std::make_pair(iter->get(), track));
0726           primaries.insert(track);
0727         }
0728       }
0729     }
0730 
0731     edm::RefToBase<Jet> jetRef = iterJets->jet();
0732 
0733     GlobalVector jetDir(jetRef->momentum().x(), jetRef->momentum().y(), jetRef->momentum().z());
0734 
0735     std::vector<std::size_t> indices = iterJets->sortedIndexes(sortCriterium);
0736 
0737     input_container trackRefs = iterJets->sortedTracks(indices);
0738 
0739     const std::vector<reco::btag::TrackIPData> &ipData = iterJets->impactParameterData();
0740 
0741     // build transient tracks used for vertex reconstruction
0742 
0743     std::vector<TransientTrack> fitTracks;
0744     std::vector<GhostTrackState> gtStates;
0745     std::unique_ptr<GhostTrackPrediction> gtPred;
0746     if (useGhostTrack)
0747       gtPred = std::make_unique<GhostTrackPrediction>(*iterJets->ghostTrack());
0748 
0749     for (unsigned int i = 0; i < indices.size(); i++) {
0750       typedef typename TemplatedSecondaryVertexTagInfo<IPTI, VTX>::IndexedTrackData IndexedTrackData;
0751 
0752       const input_item &trackRef = trackRefs[i];
0753 
0754       trackData.push_back(IndexedTrackData());
0755       trackData.back().first = indices[i];
0756 
0757       // select tracks for SV finder
0758 
0759       if (!trackSelector(
0760               *reco::btag::toTrack(trackRef), ipData[indices[i]], *jetRef, RecoVertex::convertPos(pv.position()))) {
0761         trackData.back().second.svStatus = TemplatedSecondaryVertexTagInfo<IPTI, VTX>::TrackData::trackSelected;
0762         continue;
0763       }
0764 
0765       TransientTrackMap::const_iterator pos = primariesMap.find(reco::btag::toTrack((trackRef)));
0766       TransientTrack fitTrack;
0767       if (pos != primariesMap.end()) {
0768         primaries.erase(pos->second);
0769         fitTrack = pos->second;
0770       } else
0771         fitTrack = trackBuilder->build(trackRef);
0772       fitTracks.push_back(fitTrack);
0773 
0774       trackData.back().second.svStatus = TemplatedSecondaryVertexTagInfo<IPTI, VTX>::TrackData::trackUsedForVertexFit;
0775 
0776       if (useGhostTrack) {
0777         GhostTrackState gtState(fitTrack);
0778         GlobalPoint pos = ipData[indices[i]].closestToGhostTrack;
0779         gtState.linearize(*gtPred, true, gtPred->lambda(pos));
0780         gtState.setWeight(ipData[indices[i]].ghostTrackWeight);
0781         gtStates.push_back(gtState);
0782       }
0783     }
0784 
0785     std::unique_ptr<GhostTrack> ghostTrack;
0786     if (useGhostTrack)
0787       ghostTrack = std::make_unique<GhostTrack>(
0788           GhostTrackPrediction(
0789               RecoVertex::convertPos(pv.position()),
0790               RecoVertex::convertError(pv.error()),
0791               GlobalVector(iterJets->ghostTrack()->px(), iterJets->ghostTrack()->py(), iterJets->ghostTrack()->pz()),
0792               0.05),
0793           *gtPred,
0794           gtStates,
0795           iterJets->ghostTrack()->chi2(),
0796           iterJets->ghostTrack()->ndof());
0797 
0798     // perform actual vertex finding
0799 
0800     std::vector<VTX> extAssoCollection;
0801     std::vector<TransientVertex> fittedSVs;
0802     std::vector<SecondaryVertex> SVs;
0803     if (!useExternalSV) {
0804       switch (constraint) {
0805         case CONSTRAINT_NONE:
0806           if (useGhostTrack)
0807             fittedSVs = vertexRecoGT->vertices(pv, *ghostTrack);
0808           else
0809             fittedSVs = vertexReco->vertices(fitTracks);
0810           break;
0811 
0812         case CONSTRAINT_BEAMSPOT:
0813           if (useGhostTrack)
0814             fittedSVs = vertexRecoGT->vertices(pv, *beamSpot, *ghostTrack);
0815           else
0816             fittedSVs = vertexReco->vertices(fitTracks, *beamSpot);
0817           break;
0818 
0819         case CONSTRAINT_PV_BEAMSPOT_SIZE:
0820         case CONSTRAINT_PV_BS_Z_ERRORS_SCALED:
0821         case CONSTRAINT_PV_ERROR_SCALED: {
0822           BeamSpot::CovarianceMatrix cov;
0823           for (unsigned int i = 0; i < 7; i++) {
0824             unsigned int covSrc = bsCovSrc[i];
0825             for (unsigned int j = 0; j < 7; j++) {
0826               double v = 0.0;
0827               if (!covSrc || bsCovSrc[j] != covSrc)
0828                 v = 0.0;
0829               else if (covSrc == 1)
0830                 v = beamSpot->covariance(i, j);
0831               else if (j < 3 && i < 3)
0832                 v = pv.covariance(i, j) * constraintScaling;
0833               cov(i, j) = v;
0834             }
0835           }
0836 
0837           BeamSpot bs(pv.position(),
0838                       sigmaZ,
0839                       beamSpot.isValid() ? beamSpot->dxdz() : 0.,
0840                       beamSpot.isValid() ? beamSpot->dydz() : 0.,
0841                       beamWidth,
0842                       cov,
0843                       BeamSpot::Unknown);
0844 
0845           if (useGhostTrack)
0846             fittedSVs = vertexRecoGT->vertices(pv, bs, *ghostTrack);
0847           else
0848             fittedSVs = vertexReco->vertices(fitTracks, bs);
0849         } break;
0850 
0851         case CONSTRAINT_PV_PRIMARIES_IN_FIT: {
0852           std::vector<TransientTrack> primaries_(primaries.begin(), primaries.end());
0853           if (useGhostTrack)
0854             fittedSVs = vertexRecoGT->vertices(pv, *beamSpot, primaries_, *ghostTrack);
0855           else
0856             fittedSVs = vertexReco->vertices(primaries_, fitTracks, *beamSpot);
0857         } break;
0858       }
0859       // build combined SV information and filter
0860       SVBuilder svBuilder(pv, jetDir, withPVError, minTrackWeight);
0861       std::remove_copy_if(boost::make_transform_iterator(fittedSVs.begin(), svBuilder),
0862                           boost::make_transform_iterator(fittedSVs.end(), svBuilder),
0863                           std::back_inserter(SVs),
0864                           SVFilter(vertexFilter, pv, jetDir));
0865 
0866     } else {
0867       if (useSVClustering || useFatJets) {
0868         size_t jetIdx = (iterJets - trackIPTagInfos->begin());
0869 
0870         for (size_t iExtSv = 0; iExtSv < clusteredSVs.at(jetIdx).size(); iExtSv++) {
0871           const VTX &extVertex = (*extSecVertex)[clusteredSVs.at(jetIdx).at(iExtSv)];
0872           if (extVertex.p4().M() < 0.3)
0873             continue;
0874           extAssoCollection.push_back(extVertex);
0875         }
0876       } else {
0877         for (size_t iExtSv = 0; iExtSv < extSecVertex->size(); iExtSv++) {
0878           const VTX &extVertex = (*extSecVertex)[iExtSv];
0879           if (Geom::deltaR2((position(extVertex) - pv.position()), (extSVDeltaRToJet > 0) ? jetDir : -jetDir) >
0880                   extSVDeltaRToJet * extSVDeltaRToJet ||
0881               extVertex.p4().M() < 0.3)
0882             continue;
0883           extAssoCollection.push_back(extVertex);
0884         }
0885       }
0886       // build combined SV information and filter
0887       SVBuilder svBuilder(pv, jetDir, withPVError, minTrackWeight);
0888       std::remove_copy_if(boost::make_transform_iterator(extAssoCollection.begin(), svBuilder),
0889                           boost::make_transform_iterator(extAssoCollection.end(), svBuilder),
0890                           std::back_inserter(SVs),
0891                           SVFilter(vertexFilter, pv, jetDir));
0892     }
0893     // clean up now unneeded collections
0894     gtPred.reset();
0895     ghostTrack.reset();
0896     gtStates.clear();
0897     fitTracks.clear();
0898     fittedSVs.clear();
0899     extAssoCollection.clear();
0900 
0901     // sort SVs by importance
0902 
0903     std::vector<unsigned int> vtxIndices = vertexSorting(SVs);
0904 
0905     std::vector<typename TemplatedSecondaryVertexTagInfo<IPTI, VTX>::VertexData> svData;
0906 
0907     svData.resize(vtxIndices.size());
0908     for (unsigned int idx = 0; idx < vtxIndices.size(); idx++) {
0909       const SecondaryVertex &sv = SVs[vtxIndices[idx]];
0910 
0911       svData[idx].vertex = sv;
0912       svData[idx].dist1d = sv.dist1d();
0913       svData[idx].dist2d = sv.dist2d();
0914       svData[idx].dist3d = sv.dist3d();
0915       svData[idx].direction = flightDirection(pv, sv);
0916       // mark tracks successfully used in vertex fit
0917       markUsedTracks(trackData, trackRefs, sv, idx);
0918     }
0919 
0920     // fill result into tag infos
0921 
0922     tagInfos->push_back(TemplatedSecondaryVertexTagInfo<IPTI, VTX>(
0923         trackData,
0924         svData,
0925         SVs.size(),
0926         edm::Ref<std::vector<IPTI> >(trackIPTagInfos, iterJets - trackIPTagInfos->begin())));
0927   }
0928 
0929   event.put(std::move(tagInfos));
0930 }
0931 
0932 //Need specialized template because reco::Vertex iterators are TrackBase and it is a mess to make general
0933 template <>
0934 void TemplatedSecondaryVertexProducer<TrackIPTagInfo, reco::Vertex>::markUsedTracks(TrackDataVector &trackData,
0935                                                                                     const input_container &trackRefs,
0936                                                                                     const SecondaryVertex &sv,
0937                                                                                     size_t idx) {
0938   for (Vertex::trackRef_iterator iter = sv.tracks_begin(); iter != sv.tracks_end(); ++iter) {
0939     if (sv.trackWeight(*iter) < minTrackWeight)
0940       continue;
0941 
0942     typename input_container::const_iterator pos =
0943         std::find(trackRefs.begin(), trackRefs.end(), iter->castTo<input_item>());
0944 
0945     if (pos == trackRefs.end()) {
0946       if (!useExternalSV)
0947         throw cms::Exception("TrackNotFound") << "Could not find track from secondary "
0948                                                  "vertex in original tracks."
0949                                               << std::endl;
0950     } else {
0951       unsigned int index = pos - trackRefs.begin();
0952       trackData[index].second.svStatus = (btag::TrackData::trackAssociatedToVertex + idx);
0953     }
0954   }
0955 }
0956 template <>
0957 void TemplatedSecondaryVertexProducer<CandIPTagInfo, reco::VertexCompositePtrCandidate>::markUsedTracks(
0958     TrackDataVector &trackData, const input_container &trackRefs, const SecondaryVertex &sv, size_t idx) {
0959   for (typename input_container::const_iterator iter = sv.daughterPtrVector().begin();
0960        iter != sv.daughterPtrVector().end();
0961        ++iter) {
0962     typename input_container::const_iterator pos = std::find(trackRefs.begin(), trackRefs.end(), *iter);
0963 
0964     if (pos != trackRefs.end()) {
0965       unsigned int index = pos - trackRefs.begin();
0966       trackData[index].second.svStatus = (btag::TrackData::trackAssociatedToVertex + idx);
0967     }
0968   }
0969 }
0970 
0971 template <>
0972 typename TemplatedSecondaryVertexProducer<TrackIPTagInfo, reco::Vertex>::SecondaryVertex
0973 TemplatedSecondaryVertexProducer<TrackIPTagInfo, reco::Vertex>::SVBuilder::operator()(const TransientVertex &sv) const {
0974   if (!sv.originalTracks().empty() && sv.originalTracks()[0].trackBaseRef().isNonnull())
0975     return SecondaryVertex(pv, sv, direction, withPVError);
0976   else {
0977     edm::LogError("UnexpectedInputs") << "Building from Candidates, should not happen!";
0978     return SecondaryVertex(pv, sv, direction, withPVError);
0979   }
0980 }
0981 
0982 template <>
0983 typename TemplatedSecondaryVertexProducer<CandIPTagInfo, reco::VertexCompositePtrCandidate>::SecondaryVertex
0984 TemplatedSecondaryVertexProducer<CandIPTagInfo, reco::VertexCompositePtrCandidate>::SVBuilder::operator()(
0985     const TransientVertex &sv) const {
0986   if (!sv.originalTracks().empty() && sv.originalTracks()[0].trackBaseRef().isNonnull()) {
0987     edm::LogError("UnexpectedInputs") << "Building from Tracks, should not happen!";
0988     VertexCompositePtrCandidate vtxCompPtrCand;
0989 
0990     vtxCompPtrCand.setCovariance(sv.vertexState().error().matrix());
0991     vtxCompPtrCand.setChi2AndNdof(sv.totalChiSquared(), sv.degreesOfFreedom());
0992     vtxCompPtrCand.setVertex(Candidate::Point(sv.position().x(), sv.position().y(), sv.position().z()));
0993 
0994     return SecondaryVertex(pv, vtxCompPtrCand, direction, withPVError);
0995   } else {
0996     VertexCompositePtrCandidate vtxCompPtrCand;
0997 
0998     vtxCompPtrCand.setCovariance(sv.vertexState().error().matrix());
0999     vtxCompPtrCand.setChi2AndNdof(sv.totalChiSquared(), sv.degreesOfFreedom());
1000     vtxCompPtrCand.setVertex(Candidate::Point(sv.position().x(), sv.position().y(), sv.position().z()));
1001 
1002     Candidate::LorentzVector p4;
1003     for (std::vector<reco::TransientTrack>::const_iterator tt = sv.originalTracks().begin();
1004          tt != sv.originalTracks().end();
1005          ++tt) {
1006       if (sv.trackWeight(*tt) < minTrackWeight)
1007         continue;
1008 
1009       const CandidatePtrTransientTrack *cptt =
1010           dynamic_cast<const CandidatePtrTransientTrack *>(tt->basicTransientTrack());
1011       if (cptt == nullptr)
1012         edm::LogError("DynamicCastingFailed") << "Casting of TransientTrack to CandidatePtrTransientTrack failed!";
1013       else {
1014         p4 += cptt->candidate()->p4();
1015         vtxCompPtrCand.addDaughter(cptt->candidate());
1016       }
1017     }
1018     vtxCompPtrCand.setP4(p4);
1019 
1020     return SecondaryVertex(pv, vtxCompPtrCand, direction, withPVError);
1021   }
1022 }
1023 
1024 // ------------ method that matches reclustered and original jets based on minimum dR ------------
1025 template <class IPTI, class VTX>
1026 template <class CONTAINER>
1027 void TemplatedSecondaryVertexProducer<IPTI, VTX>::matchReclusteredJets(
1028     const edm::Handle<CONTAINER> &jets,
1029     const std::vector<fastjet::PseudoJet> &reclusteredJets,
1030     std::vector<int> &matchedIndices,
1031     const std::string &jetType) {
1032   std::string type = (!jetType.empty() ? jetType + " " : jetType);
1033 
1034   std::vector<bool> matchedLocks(reclusteredJets.size(), false);
1035 
1036   for (size_t j = 0; j < jets->size(); ++j) {
1037     double matchedDR2 = 1e9;
1038     int matchedIdx = -1;
1039 
1040     for (size_t rj = 0; rj < reclusteredJets.size(); ++rj) {
1041       if (matchedLocks.at(rj))
1042         continue;  // skip jets that have already been matched
1043 
1044       double tempDR2 = Geom::deltaR2(toJet(jets->at(j))->rapidity(),
1045                                      toJet(jets->at(j))->phi(),
1046                                      reclusteredJets.at(rj).rapidity(),
1047                                      reclusteredJets.at(rj).phi_std());
1048       if (tempDR2 < matchedDR2) {
1049         matchedDR2 = tempDR2;
1050         matchedIdx = rj;
1051       }
1052     }
1053 
1054     if (matchedIdx >= 0) {
1055       if (matchedDR2 > rParam * rParam) {
1056         edm::LogError("JetMatchingFailed") << "Matched reclustered jet " << matchedIdx << " and original " << type
1057                                            << "jet " << j << " are separated by dR=" << sqrt(matchedDR2)
1058                                            << " which is greater than the jet size R=" << rParam << ".\n"
1059                                            << "This is not expected so please check that the jet algorithm and jet "
1060                                               "size match those used for the original "
1061                                            << type << "jet collection.";
1062       } else
1063         matchedLocks.at(matchedIdx) = true;
1064     } else
1065       edm::LogError("JetMatchingFailed")
1066           << "Matching reclustered to original " << type
1067           << "jets failed. Please check that the jet algorithm and jet size match those used for the original " << type
1068           << "jet collection.";
1069 
1070     matchedIndices.push_back(matchedIdx);
1071   }
1072 }
1073 
1074 // ------------ method that matches groomed and original jets based on minimum dR ------------
1075 template <class IPTI, class VTX>
1076 void TemplatedSecondaryVertexProducer<IPTI, VTX>::matchGroomedJets(const edm::Handle<edm::View<reco::Jet> > &jets,
1077                                                                    const edm::Handle<edm::View<reco::Jet> > &groomedJets,
1078                                                                    std::vector<int> &matchedIndices) {
1079   std::vector<bool> jetLocks(jets->size(), false);
1080   std::vector<int> jetIndices;
1081 
1082   for (size_t gj = 0; gj < groomedJets->size(); ++gj) {
1083     double matchedDR2 = 1e9;
1084     int matchedIdx = -1;
1085 
1086     if (groomedJets->at(gj).pt() > 0.)  // skip pathological cases of groomed jets with Pt=0
1087     {
1088       for (size_t j = 0; j < jets->size(); ++j) {
1089         if (jetLocks.at(j))
1090           continue;  // skip jets that have already been matched
1091 
1092         double tempDR2 = Geom::deltaR2(
1093             jets->at(j).rapidity(), jets->at(j).phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi());
1094         if (tempDR2 < matchedDR2) {
1095           matchedDR2 = tempDR2;
1096           matchedIdx = j;
1097         }
1098       }
1099     }
1100 
1101     if (matchedIdx >= 0) {
1102       if (matchedDR2 > rParam * rParam) {
1103         edm::LogWarning("MatchedJetsFarApart")
1104             << "Matched groomed jet " << gj << " and original jet " << matchedIdx
1105             << " are separated by dR=" << sqrt(matchedDR2) << " which is greater than the jet size R=" << rParam
1106             << ".\n"
1107             << "This is not expected so the matching of these two jets has been discarded. Please check that the two "
1108                "jet collections belong to each other.";
1109         matchedIdx = -1;
1110       } else
1111         jetLocks.at(matchedIdx) = true;
1112     }
1113     jetIndices.push_back(matchedIdx);
1114   }
1115 
1116   for (size_t j = 0; j < jets->size(); ++j) {
1117     std::vector<int>::iterator matchedIndex = std::find(jetIndices.begin(), jetIndices.end(), j);
1118 
1119     matchedIndices.push_back(matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(), matchedIndex) : -1);
1120   }
1121 }
1122 
1123 // ------------ method that matches subjets and original fat jets ------------
1124 template <class IPTI, class VTX>
1125 void TemplatedSecondaryVertexProducer<IPTI, VTX>::matchSubjets(const std::vector<int> &groomedIndices,
1126                                                                const edm::Handle<edm::View<reco::Jet> > &groomedJets,
1127                                                                const edm::Handle<std::vector<IPTI> > &subjets,
1128                                                                std::vector<std::vector<int> > &matchedIndices) {
1129   for (size_t g = 0; g < groomedIndices.size(); ++g) {
1130     std::vector<int> subjetIndices;
1131 
1132     if (groomedIndices.at(g) >= 0) {
1133       for (size_t s = 0; s < groomedJets->at(groomedIndices.at(g)).numberOfDaughters(); ++s) {
1134         const edm::Ptr<reco::Candidate> &subjet = groomedJets->at(groomedIndices.at(g)).daughterPtr(s);
1135 
1136         for (size_t sj = 0; sj < subjets->size(); ++sj) {
1137           const edm::RefToBase<reco::Jet> &subjetRef = subjets->at(sj).jet();
1138           if (subjet == edm::Ptr<reco::Candidate>(subjetRef.id(), subjetRef.get(), subjetRef.key())) {
1139             subjetIndices.push_back(sj);
1140             break;
1141           }
1142         }
1143       }
1144 
1145       if (subjetIndices.empty())
1146         edm::LogError("SubjetMatchingFailed") << "Matching subjets to original fat jets failed. Please check that the "
1147                                                  "groomed fat jet and subjet collections belong to each other.";
1148 
1149       matchedIndices.push_back(subjetIndices);
1150     } else
1151       matchedIndices.push_back(subjetIndices);
1152   }
1153 }
1154 
1155 // ------------ method that matches subjets and original fat jets ------------
1156 template <class IPTI, class VTX>
1157 void TemplatedSecondaryVertexProducer<IPTI, VTX>::matchSubjets(const edm::Handle<edm::View<reco::Jet> > &fatJets,
1158                                                                const edm::Handle<std::vector<IPTI> > &subjets,
1159                                                                std::vector<std::vector<int> > &matchedIndices) {
1160   for (size_t fj = 0; fj < fatJets->size(); ++fj) {
1161     std::vector<int> subjetIndices;
1162     size_t nSubjetCollections = 0;
1163     size_t nSubjets = 0;
1164 
1165     const pat::Jet *fatJet = dynamic_cast<const pat::Jet *>(fatJets->ptrAt(fj).get());
1166 
1167     if (!fatJet) {
1168       if (fj == 0)
1169         edm::LogError("WrongJetType")
1170             << "Wrong jet type for input fat jets. Please check that the input fat jets are of the pat::Jet type.";
1171 
1172       matchedIndices.push_back(subjetIndices);
1173       continue;
1174     } else {
1175       nSubjetCollections = fatJet->subjetCollectionNames().size();
1176 
1177       if (nSubjetCollections > 0) {
1178         for (size_t coll = 0; coll < nSubjetCollections; ++coll) {
1179           const pat::JetPtrCollection &fatJetSubjets = fatJet->subjets(coll);
1180 
1181           for (size_t fjsj = 0; fjsj < fatJetSubjets.size(); ++fjsj) {
1182             ++nSubjets;
1183 
1184             for (size_t sj = 0; sj < subjets->size(); ++sj) {
1185               const pat::Jet *subJet = dynamic_cast<const pat::Jet *>(subjets->at(sj).jet().get());
1186 
1187               if (!subJet) {
1188                 if (fj == 0 && coll == 0 && fjsj == 0 && sj == 0)
1189                   edm::LogError("WrongJetType") << "Wrong jet type for input subjets. Please check that the input "
1190                                                    "subjets are of the pat::Jet type.";
1191 
1192                 break;
1193               } else {
1194                 if (subJet->originalObjectRef() == fatJetSubjets.at(fjsj)->originalObjectRef()) {
1195                   subjetIndices.push_back(sj);
1196                   break;
1197                 }
1198               }
1199             }
1200           }
1201         }
1202 
1203         if (subjetIndices.empty() && nSubjets > 0)
1204           edm::LogError("SubjetMatchingFailed") << "Matching subjets to fat jets failed. Please check that the fat jet "
1205                                                    "and subjet collections belong to each other.";
1206 
1207         matchedIndices.push_back(subjetIndices);
1208       } else
1209         matchedIndices.push_back(subjetIndices);
1210     }
1211   }
1212 }
1213 
1214 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
1215 template <class IPTI, class VTX>
1216 void TemplatedSecondaryVertexProducer<IPTI, VTX>::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
1217   edm::ParameterSetDescription desc;
1218   desc.add<double>("extSVDeltaRToJet", 0.3);
1219   desc.add<edm::InputTag>("beamSpotTag", edm::InputTag("offlineBeamSpot"));
1220   {
1221     edm::ParameterSetDescription vertexReco;
1222     vertexReco.add<double>("primcut", 1.8);
1223     vertexReco.add<double>("seccut", 6.0);
1224     vertexReco.add<std::string>("finder", "avr");
1225     vertexReco.addOptionalNode(edm::ParameterDescription<double>("minweight", 0.5, true) and
1226                                    edm::ParameterDescription<double>("weightthreshold", 0.001, true) and
1227                                    edm::ParameterDescription<bool>("smoothing", false, true),
1228                                true);
1229     vertexReco.addOptionalNode(
1230         edm::ParameterDescription<double>("maxFitChi2", 10.0, true) and
1231             edm::ParameterDescription<double>("mergeThreshold", 3.0, true) and
1232             edm::ParameterDescription<std::string>("fitType", "RefitGhostTrackWithVertices", true),
1233         true);
1234     desc.add<edm::ParameterSetDescription>("vertexReco", vertexReco);
1235   }
1236   {
1237     edm::ParameterSetDescription vertexSelection;
1238     vertexSelection.add<std::string>("sortCriterium", "dist3dError");
1239     desc.add<edm::ParameterSetDescription>("vertexSelection", vertexSelection);
1240   }
1241   desc.add<std::string>("constraint", "BeamSpot");
1242   desc.add<edm::InputTag>("trackIPTagInfos", edm::InputTag("impactParameterTagInfos"));
1243   {
1244     edm::ParameterSetDescription vertexCuts;
1245     vertexCuts.add<double>("distSig3dMax", 99999.9);
1246     vertexCuts.add<double>("fracPV", 0.65);
1247     vertexCuts.add<double>("distVal2dMax", 2.5);
1248     vertexCuts.add<bool>("useTrackWeights", true);
1249     vertexCuts.add<double>("maxDeltaRToJetAxis", 0.4);
1250     {
1251       edm::ParameterSetDescription v0Filter;
1252       v0Filter.add<double>("k0sMassWindow", 0.05);
1253       vertexCuts.add<edm::ParameterSetDescription>("v0Filter", v0Filter);
1254     }
1255     vertexCuts.add<double>("distSig2dMin", 3.0);
1256     vertexCuts.add<unsigned int>("multiplicityMin", 2);
1257     vertexCuts.add<double>("distVal2dMin", 0.01);
1258     vertexCuts.add<double>("distSig2dMax", 99999.9);
1259     vertexCuts.add<double>("distVal3dMax", 99999.9);
1260     vertexCuts.add<double>("minimumTrackWeight", 0.5);
1261     vertexCuts.add<double>("distVal3dMin", -99999.9);
1262     vertexCuts.add<double>("massMax", 6.5);
1263     vertexCuts.add<double>("distSig3dMin", -99999.9);
1264     desc.add<edm::ParameterSetDescription>("vertexCuts", vertexCuts);
1265   }
1266   desc.add<bool>("useExternalSV", false);
1267   desc.add<double>("minimumTrackWeight", 0.5);
1268   desc.add<bool>("usePVError", true);
1269   {
1270     edm::ParameterSetDescription trackSelection;
1271     trackSelection.add<double>("b_pT", 0.3684);
1272     trackSelection.add<double>("max_pT", 500);
1273     trackSelection.add<bool>("useVariableJTA", false);
1274     trackSelection.add<double>("maxDecayLen", 99999.9);
1275     trackSelection.add<double>("sip3dValMin", -99999.9);
1276     trackSelection.add<double>("max_pT_dRcut", 0.1);
1277     trackSelection.add<double>("a_pT", 0.005263);
1278     trackSelection.add<unsigned int>("totalHitsMin", 8);
1279     trackSelection.add<double>("jetDeltaRMax", 0.3);
1280     trackSelection.add<double>("a_dR", -0.001053);
1281     trackSelection.add<double>("maxDistToAxis", 0.2);
1282     trackSelection.add<double>("ptMin", 1.0);
1283     trackSelection.add<std::string>("qualityClass", "any");
1284     trackSelection.add<unsigned int>("pixelHitsMin", 2);
1285     trackSelection.add<double>("sip2dValMax", 99999.9);
1286     trackSelection.add<double>("max_pT_trackPTcut", 3);
1287     trackSelection.add<double>("sip2dValMin", -99999.9);
1288     trackSelection.add<double>("normChi2Max", 99999.9);
1289     trackSelection.add<double>("sip3dValMax", 99999.9);
1290     trackSelection.add<double>("sip3dSigMin", -99999.9);
1291     trackSelection.add<double>("min_pT", 120);
1292     trackSelection.add<double>("min_pT_dRcut", 0.5);
1293     trackSelection.add<double>("sip2dSigMax", 99999.9);
1294     trackSelection.add<double>("sip3dSigMax", 99999.9);
1295     trackSelection.add<double>("sip2dSigMin", -99999.9);
1296     trackSelection.add<double>("b_dR", 0.6263);
1297     desc.add<edm::ParameterSetDescription>("trackSelection", trackSelection);
1298   }
1299   desc.add<std::string>("trackSort", "sip3dSig");
1300   desc.add<edm::InputTag>("extSVCollection", edm::InputTag("secondaryVertices"));
1301   desc.addOptionalNode(edm::ParameterDescription<bool>("useSVClustering", false, true) and
1302                            edm::ParameterDescription<std::string>("jetAlgorithm", true) and
1303                            edm::ParameterDescription<double>("rParam", true),
1304                        true);
1305   desc.addOptional<bool>("useSVMomentum", false);
1306   desc.addOptional<double>("ghostRescaling", 1e-18);
1307   desc.addOptional<double>("relPtTolerance", 1e-03);
1308   desc.addOptional<edm::InputTag>("fatJets");
1309   desc.addOptional<edm::InputTag>("groomedFatJets");
1310   desc.add<edm::InputTag>("weights", edm::InputTag(""));
1311   descriptions.addDefault(desc);
1312 }
1313 
1314 //define this as a plug-in
1315 typedef TemplatedSecondaryVertexProducer<TrackIPTagInfo, reco::Vertex> SecondaryVertexProducer;
1316 typedef TemplatedSecondaryVertexProducer<CandIPTagInfo, reco::VertexCompositePtrCandidate> CandSecondaryVertexProducer;
1317 
1318 DEFINE_FWK_MODULE(SecondaryVertexProducer);
1319 DEFINE_FWK_MODULE(CandSecondaryVertexProducer);
[ 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!