Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-29 23:10:56

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           reco::CandidatePtr constit = *m;
0392           if (constit->pt() == 0) {
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             fjInputs.push_back(
0403                 fastjet::PseudoJet(constit->px() * w, constit->py() * w, constit->pz() * w, constit->energy() * w));
0404           } else {
0405             fjInputs.push_back(fastjet::PseudoJet(constit->px(), constit->py(), constit->pz(), constit->energy()));
0406           }
0407         }
0408       }
0409     } else {
0410       for (typename std::vector<IPTI>::const_iterator it = trackIPTagInfos->begin(); it != trackIPTagInfos->end();
0411            ++it) {
0412         std::vector<edm::Ptr<reco::Candidate> > constituents = it->jet()->getJetConstituents();
0413         std::vector<edm::Ptr<reco::Candidate> >::const_iterator m;
0414         for (m = constituents.begin(); m != constituents.end(); ++m) {
0415           reco::CandidatePtr constit = *m;
0416           if (constit->pt() == 0) {
0417             edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
0418             continue;
0419           }
0420           if (it->jet()->isWeighted()) {
0421             if (token_weights.isUninitialized())
0422               throw cms::Exception("MissingConstituentWeight")
0423                   << "TemplatedSecondaryVertexProducer: No weights (e.g. PUPPI) given for weighted jet collection"
0424                   << std::endl;
0425             float w = (*weightsHandle)[constit];
0426             fjInputs.push_back(
0427                 fastjet::PseudoJet(constit->px() * w, constit->py() * w, constit->pz() * w, constit->energy() * w));
0428           } else {
0429             fjInputs.push_back(fastjet::PseudoJet(constit->px(), constit->py(), constit->pz(), constit->energy()));
0430           }
0431         }
0432       }
0433     }
0434     // insert "ghost" SVs in the vector of constituents
0435     for (typename edm::View<VTX>::const_iterator it = extSecVertex->begin(); it != extSecVertex->end(); ++it) {
0436       const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
0437       GlobalVector dir = flightDirection(pv, *it);
0438       dir = dir.unit();
0439       fastjet::PseudoJet p(
0440           dir.x(), dir.y(), dir.z(), dir.mag());  // using SV flight direction so treating SV as massless
0441       if (useSVMomentum)
0442         p = fastjet::PseudoJet(it->p4().px(), it->p4().py(), it->p4().pz(), it->p4().energy());
0443       p *= ghostRescaling;  // rescale SV direction/momentum
0444       p.set_user_info(new VertexInfo(it - extSecVertex->begin()));
0445       fjInputs.push_back(p);
0446     }
0447 
0448     // define jet clustering sequence
0449     fjClusterSeq = std::make_shared<fastjet::ClusterSequence>(fjInputs, *fjJetDefinition);
0450     // recluster jet constituents and inserted "ghosts"
0451     std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(fjClusterSeq->inclusive_jets(jetPtMin));
0452 
0453     if (useFatJets) {
0454       if (inclusiveJets.size() < fatJetsHandle->size())
0455         edm::LogError("TooFewReclusteredJets")
0456             << "There are fewer reclustered (" << inclusiveJets.size() << ") than original fat jets ("
0457             << fatJetsHandle->size()
0458             << "). Please check that the jet algorithm and jet size match those used for the original jet collection.";
0459 
0460       // match reclustered and original fat jets
0461       std::vector<int> reclusteredIndices;
0462       matchReclusteredJets<edm::View<reco::Jet> >(fatJetsHandle, inclusiveJets, reclusteredIndices, "fat");
0463 
0464       // match groomed and original fat jets
0465       std::vector<int> groomedIndices;
0466       if (useGroomedFatJets)
0467         matchGroomedJets(fatJetsHandle, groomedFatJetsHandle, groomedIndices);
0468 
0469       // match subjets and original fat jets
0470       std::vector<std::vector<int> > subjetIndices;
0471       if (useGroomedFatJets)
0472         matchSubjets(groomedIndices, groomedFatJetsHandle, trackIPTagInfos, subjetIndices);
0473       else
0474         matchSubjets(fatJetsHandle, trackIPTagInfos, subjetIndices);
0475 
0476       // collect clustered SVs
0477       for (size_t i = 0; i < fatJetsHandle->size(); ++i) {
0478         if (reclusteredIndices.at(i) < 0)
0479           continue;  // continue if matching reclustered to original jets failed
0480 
0481         if (fatJetsHandle->at(i).pt() == 0)  // continue if the original jet has Pt=0
0482         {
0483           edm::LogWarning("NullTransverseMomentum")
0484               << "The original fat jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
0485           continue;
0486         }
0487 
0488         if (subjetIndices.at(i).empty())
0489           continue;  // continue if the original jet does not have subjets assigned
0490 
0491         // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original fat jets should in principle stay the same
0492         if ((std::abs(inclusiveJets.at(reclusteredIndices.at(i)).pt() - fatJetsHandle->at(i).pt()) /
0493              fatJetsHandle->at(i).pt()) > relPtTolerance) {
0494           if (fatJetsHandle->at(i).pt() < 10.)  // special handling for low-Pt jets (Pt<10 GeV)
0495             edm::LogWarning("JetPtMismatchAtLowPt")
0496                 << "The reclustered and original fat jet " << i << " have different Pt's ("
0497                 << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << fatJetsHandle->at(i).pt()
0498                 << " GeV, respectively).\n"
0499                 << "Please check that the jet algorithm and jet size match those used for the original fat jet "
0500                    "collection and also make sure the original fat jets are uncorrected. In addition, make sure you "
0501                    "are not using CaloJets which are presently not supported.\n"
0502                 << "Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n"
0503                 << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
0504                    "precision in which case make sure the original jet collection is produced and reclustering is "
0505                    "performed in the same job.";
0506           else
0507             edm::LogError("JetPtMismatch")
0508                 << "The reclustered and original fat jet " << i << " have different Pt's ("
0509                 << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << fatJetsHandle->at(i).pt()
0510                 << " GeV, respectively).\n"
0511                 << "Please check that the jet algorithm and jet size match those used for the original fat jet "
0512                    "collection and also make sure the original fat jets are uncorrected. In addition, make sure you "
0513                    "are not using CaloJets which are presently not supported.\n"
0514                 << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
0515                    "precision in which case make sure the original jet collection is produced and reclustering is "
0516                    "performed in the same job.";
0517         }
0518 
0519         // get jet constituents
0520         std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(i)).constituents();
0521 
0522         std::vector<int> svIndices;
0523         // loop over jet constituents and try to find "ghosts"
0524         for (std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end();
0525              ++it) {
0526           if (!it->has_user_info())
0527             continue;  // skip if not a "ghost"
0528 
0529           svIndices.push_back(it->user_info<VertexInfo>().vertexIndex());
0530         }
0531 
0532         // loop over clustered SVs and assign them to different subjets based on smallest dR
0533         for (size_t sv = 0; sv < svIndices.size(); ++sv) {
0534           const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
0535           const VTX &extSV = (*extSecVertex)[svIndices.at(sv)];
0536           GlobalVector dir = flightDirection(pv, extSV);
0537           dir = dir.unit();
0538           fastjet::PseudoJet p(
0539               dir.x(), dir.y(), dir.z(), dir.mag());  // using SV flight direction so treating SV as massless
0540           if (useSVMomentum)
0541             p = fastjet::PseudoJet(extSV.p4().px(), extSV.p4().py(), extSV.p4().pz(), extSV.p4().energy());
0542 
0543           std::vector<double> dR2toSubjets;
0544 
0545           for (size_t sj = 0; sj < subjetIndices.at(i).size(); ++sj)
0546             dR2toSubjets.push_back(Geom::deltaR2(p.rapidity(),
0547                                                  p.phi_std(),
0548                                                  trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->rapidity(),
0549                                                  trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->phi()));
0550 
0551           // find the closest subjet
0552           int closestSubjetIdx =
0553               std::distance(dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()));
0554 
0555           clusteredSVs.at(subjetIndices.at(i).at(closestSubjetIdx)).push_back(svIndices.at(sv));
0556         }
0557       }
0558     } else {
0559       if (inclusiveJets.size() < trackIPTagInfos->size())
0560         edm::LogError("TooFewReclusteredJets")
0561             << "There are fewer reclustered (" << inclusiveJets.size() << ") than original jets ("
0562             << trackIPTagInfos->size()
0563             << "). Please check that the jet algorithm and jet size match those used for the original jet collection.";
0564 
0565       // match reclustered and original jets
0566       std::vector<int> reclusteredIndices;
0567       matchReclusteredJets<std::vector<IPTI> >(trackIPTagInfos, inclusiveJets, reclusteredIndices);
0568 
0569       // collect clustered SVs
0570       for (size_t i = 0; i < trackIPTagInfos->size(); ++i) {
0571         if (reclusteredIndices.at(i) < 0)
0572           continue;  // continue if matching reclustered to original jets failed
0573 
0574         if (trackIPTagInfos->at(i).jet()->pt() == 0)  // continue if the original jet has Pt=0
0575         {
0576           edm::LogWarning("NullTransverseMomentum")
0577               << "The original jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
0578           continue;
0579         }
0580 
0581         // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original jets should in principle stay the same
0582         if ((std::abs(inclusiveJets.at(reclusteredIndices.at(i)).pt() - trackIPTagInfos->at(i).jet()->pt()) /
0583              trackIPTagInfos->at(i).jet()->pt()) > relPtTolerance) {
0584           if (trackIPTagInfos->at(i).jet()->pt() < 10.)  // special handling for low-Pt jets (Pt<10 GeV)
0585             edm::LogWarning("JetPtMismatchAtLowPt")
0586                 << "The reclustered and original jet " << i << " have different Pt's ("
0587                 << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << trackIPTagInfos->at(i).jet()->pt()
0588                 << " GeV, respectively).\n"
0589                 << "Please check that the jet algorithm and jet size match those used for the original jet collection "
0590                    "and also make sure the original jets are uncorrected. In addition, make sure you are not using "
0591                    "CaloJets which are presently not supported.\n"
0592                 << "Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n"
0593                 << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
0594                    "precision in which case make sure the original jet collection is produced and reclustering is "
0595                    "performed in the same job.";
0596           else
0597             edm::LogError("JetPtMismatch")
0598                 << "The reclustered and original jet " << i << " have different Pt's ("
0599                 << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << trackIPTagInfos->at(i).jet()->pt()
0600                 << " GeV, respectively).\n"
0601                 << "Please check that the jet algorithm and jet size match those used for the original jet collection "
0602                    "and also make sure the original jets are uncorrected. In addition, make sure you are not using "
0603                    "CaloJets which are presently not supported.\n"
0604                 << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
0605                    "precision in which case make sure the original jet collection is produced and reclustering is "
0606                    "performed in the same job.";
0607         }
0608 
0609         // get jet constituents
0610         std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(i)).constituents();
0611 
0612         // loop over jet constituents and try to find "ghosts"
0613         for (std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end();
0614              ++it) {
0615           if (!it->has_user_info())
0616             continue;  // skip if not a "ghost"
0617           // push back clustered SV indices
0618           clusteredSVs.at(i).push_back(it->user_info<VertexInfo>().vertexIndex());
0619         }
0620       }
0621     }
0622   }
0623   // case where fat jets are used to associate SVs to subjets but no SV clustering is performed
0624   else if (useExternalSV && !useSVClustering && !trackIPTagInfos->empty() && useFatJets) {
0625     // match groomed and original fat jets
0626     std::vector<int> groomedIndices;
0627     if (useGroomedFatJets)
0628       matchGroomedJets(fatJetsHandle, groomedFatJetsHandle, groomedIndices);
0629 
0630     // match subjets and original fat jets
0631     std::vector<std::vector<int> > subjetIndices;
0632     if (useGroomedFatJets)
0633       matchSubjets(groomedIndices, groomedFatJetsHandle, trackIPTagInfos, subjetIndices);
0634     else
0635       matchSubjets(fatJetsHandle, trackIPTagInfos, subjetIndices);
0636 
0637     // loop over fat jets
0638     for (size_t i = 0; i < fatJetsHandle->size(); ++i) {
0639       if (fatJetsHandle->at(i).pt() == 0)  // continue if the original jet has Pt=0
0640       {
0641         edm::LogWarning("NullTransverseMomentum")
0642             << "The original fat jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
0643         continue;
0644       }
0645 
0646       if (subjetIndices.at(i).empty())
0647         continue;  // continue if the original jet does not have subjets assigned
0648 
0649       // loop over SVs, associate them to fat jets based on dR cone and
0650       // then assign them to the closets subjet in dR
0651       for (typename edm::View<VTX>::const_iterator it = extSecVertex->begin(); it != extSecVertex->end(); ++it) {
0652         size_t sv = (it - extSecVertex->begin());
0653 
0654         const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
0655         const VTX &extSV = (*extSecVertex)[sv];
0656         GlobalVector dir = flightDirection(pv, extSV);
0657         GlobalVector jetDir(fatJetsHandle->at(i).px(), fatJetsHandle->at(i).py(), fatJetsHandle->at(i).pz());
0658         // skip SVs outside the dR cone
0659         if (Geom::deltaR2(dir, jetDir) > rParam * rParam)  // here using the jet clustering rParam as a dR cut
0660           continue;
0661 
0662         dir = dir.unit();
0663         fastjet::PseudoJet p(
0664             dir.x(), dir.y(), dir.z(), dir.mag());  // using SV flight direction so treating SV as massless
0665         if (useSVMomentum)
0666           p = fastjet::PseudoJet(extSV.p4().px(), extSV.p4().py(), extSV.p4().pz(), extSV.p4().energy());
0667 
0668         std::vector<double> dR2toSubjets;
0669 
0670         for (size_t sj = 0; sj < subjetIndices.at(i).size(); ++sj)
0671           dR2toSubjets.push_back(Geom::deltaR2(p.rapidity(),
0672                                                p.phi_std(),
0673                                                trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->rapidity(),
0674                                                trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->phi()));
0675 
0676         // find the closest subjet
0677         int closestSubjetIdx =
0678             std::distance(dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()));
0679 
0680         clusteredSVs.at(subjetIndices.at(i).at(closestSubjetIdx)).push_back(sv);
0681       }
0682     }
0683   }
0684   // ------------------------------------ SV clustering END ----------------------------------------------
0685 
0686   std::unique_ptr<ConfigurableVertexReconstructor> vertexReco;
0687   std::unique_ptr<GhostTrackVertexFinder> vertexRecoGT;
0688   if (useGhostTrack)
0689     vertexRecoGT = std::make_unique<GhostTrackVertexFinder>(
0690         vtxRecoPSet.getParameter<double>("maxFitChi2"),
0691         vtxRecoPSet.getParameter<double>("mergeThreshold"),
0692         vtxRecoPSet.getParameter<double>("primcut"),
0693         vtxRecoPSet.getParameter<double>("seccut"),
0694         getGhostTrackFitType(vtxRecoPSet.getParameter<std::string>("fitType")));
0695   else
0696     vertexReco = std::make_unique<ConfigurableVertexReconstructor>(vtxRecoPSet);
0697 
0698   TransientTrackMap primariesMap;
0699 
0700   // result secondary vertices
0701 
0702   auto tagInfos = std::make_unique<Product>();
0703 
0704   for (typename std::vector<IPTI>::const_iterator iterJets = trackIPTagInfos->begin();
0705        iterJets != trackIPTagInfos->end();
0706        ++iterJets) {
0707     TrackDataVector trackData;
0708     //            std::cout << "Jet " << iterJets-trackIPTagInfos->begin() << std::endl;
0709 
0710     const Vertex &pv = *iterJets->primaryVertex();
0711 
0712     std::set<TransientTrack> primaries;
0713     if (constraint == CONSTRAINT_PV_PRIMARIES_IN_FIT) {
0714       for (Vertex::trackRef_iterator iter = pv.tracks_begin(); iter != pv.tracks_end(); ++iter) {
0715         TransientTrackMap::iterator pos = primariesMap.lower_bound(iter->get());
0716 
0717         if (pos != primariesMap.end() && pos->first == iter->get())
0718           primaries.insert(pos->second);
0719         else {
0720           TransientTrack track = trackBuilder->build(iter->castTo<TrackRef>());
0721           primariesMap.insert(pos, std::make_pair(iter->get(), track));
0722           primaries.insert(track);
0723         }
0724       }
0725     }
0726 
0727     edm::RefToBase<Jet> jetRef = iterJets->jet();
0728 
0729     GlobalVector jetDir(jetRef->momentum().x(), jetRef->momentum().y(), jetRef->momentum().z());
0730 
0731     std::vector<std::size_t> indices = iterJets->sortedIndexes(sortCriterium);
0732 
0733     input_container trackRefs = iterJets->sortedTracks(indices);
0734 
0735     const std::vector<reco::btag::TrackIPData> &ipData = iterJets->impactParameterData();
0736 
0737     // build transient tracks used for vertex reconstruction
0738 
0739     std::vector<TransientTrack> fitTracks;
0740     std::vector<GhostTrackState> gtStates;
0741     std::unique_ptr<GhostTrackPrediction> gtPred;
0742     if (useGhostTrack)
0743       gtPred = std::make_unique<GhostTrackPrediction>(*iterJets->ghostTrack());
0744 
0745     for (unsigned int i = 0; i < indices.size(); i++) {
0746       typedef typename TemplatedSecondaryVertexTagInfo<IPTI, VTX>::IndexedTrackData IndexedTrackData;
0747 
0748       const input_item &trackRef = trackRefs[i];
0749 
0750       trackData.push_back(IndexedTrackData());
0751       trackData.back().first = indices[i];
0752 
0753       // select tracks for SV finder
0754 
0755       if (!trackSelector(
0756               *reco::btag::toTrack(trackRef), ipData[indices[i]], *jetRef, RecoVertex::convertPos(pv.position()))) {
0757         trackData.back().second.svStatus = TemplatedSecondaryVertexTagInfo<IPTI, VTX>::TrackData::trackSelected;
0758         continue;
0759       }
0760 
0761       TransientTrackMap::const_iterator pos = primariesMap.find(reco::btag::toTrack((trackRef)));
0762       TransientTrack fitTrack;
0763       if (pos != primariesMap.end()) {
0764         primaries.erase(pos->second);
0765         fitTrack = pos->second;
0766       } else
0767         fitTrack = trackBuilder->build(trackRef);
0768       fitTracks.push_back(fitTrack);
0769 
0770       trackData.back().second.svStatus = TemplatedSecondaryVertexTagInfo<IPTI, VTX>::TrackData::trackUsedForVertexFit;
0771 
0772       if (useGhostTrack) {
0773         GhostTrackState gtState(fitTrack);
0774         GlobalPoint pos = ipData[indices[i]].closestToGhostTrack;
0775         gtState.linearize(*gtPred, true, gtPred->lambda(pos));
0776         gtState.setWeight(ipData[indices[i]].ghostTrackWeight);
0777         gtStates.push_back(gtState);
0778       }
0779     }
0780 
0781     std::unique_ptr<GhostTrack> ghostTrack;
0782     if (useGhostTrack)
0783       ghostTrack = std::make_unique<GhostTrack>(
0784           GhostTrackPrediction(
0785               RecoVertex::convertPos(pv.position()),
0786               RecoVertex::convertError(pv.error()),
0787               GlobalVector(iterJets->ghostTrack()->px(), iterJets->ghostTrack()->py(), iterJets->ghostTrack()->pz()),
0788               0.05),
0789           *gtPred,
0790           gtStates,
0791           iterJets->ghostTrack()->chi2(),
0792           iterJets->ghostTrack()->ndof());
0793 
0794     // perform actual vertex finding
0795 
0796     std::vector<VTX> extAssoCollection;
0797     std::vector<TransientVertex> fittedSVs;
0798     std::vector<SecondaryVertex> SVs;
0799     if (!useExternalSV) {
0800       switch (constraint) {
0801         case CONSTRAINT_NONE:
0802           if (useGhostTrack)
0803             fittedSVs = vertexRecoGT->vertices(pv, *ghostTrack);
0804           else
0805             fittedSVs = vertexReco->vertices(fitTracks);
0806           break;
0807 
0808         case CONSTRAINT_BEAMSPOT:
0809           if (useGhostTrack)
0810             fittedSVs = vertexRecoGT->vertices(pv, *beamSpot, *ghostTrack);
0811           else
0812             fittedSVs = vertexReco->vertices(fitTracks, *beamSpot);
0813           break;
0814 
0815         case CONSTRAINT_PV_BEAMSPOT_SIZE:
0816         case CONSTRAINT_PV_BS_Z_ERRORS_SCALED:
0817         case CONSTRAINT_PV_ERROR_SCALED: {
0818           BeamSpot::CovarianceMatrix cov;
0819           for (unsigned int i = 0; i < 7; i++) {
0820             unsigned int covSrc = bsCovSrc[i];
0821             for (unsigned int j = 0; j < 7; j++) {
0822               double v = 0.0;
0823               if (!covSrc || bsCovSrc[j] != covSrc)
0824                 v = 0.0;
0825               else if (covSrc == 1)
0826                 v = beamSpot->covariance(i, j);
0827               else if (j < 3 && i < 3)
0828                 v = pv.covariance(i, j) * constraintScaling;
0829               cov(i, j) = v;
0830             }
0831           }
0832 
0833           BeamSpot bs(pv.position(),
0834                       sigmaZ,
0835                       beamSpot.isValid() ? beamSpot->dxdz() : 0.,
0836                       beamSpot.isValid() ? beamSpot->dydz() : 0.,
0837                       beamWidth,
0838                       cov,
0839                       BeamSpot::Unknown);
0840 
0841           if (useGhostTrack)
0842             fittedSVs = vertexRecoGT->vertices(pv, bs, *ghostTrack);
0843           else
0844             fittedSVs = vertexReco->vertices(fitTracks, bs);
0845         } break;
0846 
0847         case CONSTRAINT_PV_PRIMARIES_IN_FIT: {
0848           std::vector<TransientTrack> primaries_(primaries.begin(), primaries.end());
0849           if (useGhostTrack)
0850             fittedSVs = vertexRecoGT->vertices(pv, *beamSpot, primaries_, *ghostTrack);
0851           else
0852             fittedSVs = vertexReco->vertices(primaries_, fitTracks, *beamSpot);
0853         } break;
0854       }
0855       // build combined SV information and filter
0856       SVBuilder svBuilder(pv, jetDir, withPVError, minTrackWeight);
0857       std::remove_copy_if(boost::make_transform_iterator(fittedSVs.begin(), svBuilder),
0858                           boost::make_transform_iterator(fittedSVs.end(), svBuilder),
0859                           std::back_inserter(SVs),
0860                           SVFilter(vertexFilter, pv, jetDir));
0861 
0862     } else {
0863       if (useSVClustering || useFatJets) {
0864         size_t jetIdx = (iterJets - trackIPTagInfos->begin());
0865 
0866         for (size_t iExtSv = 0; iExtSv < clusteredSVs.at(jetIdx).size(); iExtSv++) {
0867           const VTX &extVertex = (*extSecVertex)[clusteredSVs.at(jetIdx).at(iExtSv)];
0868           if (extVertex.p4().M() < 0.3)
0869             continue;
0870           extAssoCollection.push_back(extVertex);
0871         }
0872       } else {
0873         for (size_t iExtSv = 0; iExtSv < extSecVertex->size(); iExtSv++) {
0874           const VTX &extVertex = (*extSecVertex)[iExtSv];
0875           if (Geom::deltaR2((position(extVertex) - pv.position()), (extSVDeltaRToJet > 0) ? jetDir : -jetDir) >
0876                   extSVDeltaRToJet * extSVDeltaRToJet ||
0877               extVertex.p4().M() < 0.3)
0878             continue;
0879           extAssoCollection.push_back(extVertex);
0880         }
0881       }
0882       // build combined SV information and filter
0883       SVBuilder svBuilder(pv, jetDir, withPVError, minTrackWeight);
0884       std::remove_copy_if(boost::make_transform_iterator(extAssoCollection.begin(), svBuilder),
0885                           boost::make_transform_iterator(extAssoCollection.end(), svBuilder),
0886                           std::back_inserter(SVs),
0887                           SVFilter(vertexFilter, pv, jetDir));
0888     }
0889     // clean up now unneeded collections
0890     gtPred.reset();
0891     ghostTrack.reset();
0892     gtStates.clear();
0893     fitTracks.clear();
0894     fittedSVs.clear();
0895     extAssoCollection.clear();
0896 
0897     // sort SVs by importance
0898 
0899     std::vector<unsigned int> vtxIndices = vertexSorting(SVs);
0900 
0901     std::vector<typename TemplatedSecondaryVertexTagInfo<IPTI, VTX>::VertexData> svData;
0902 
0903     svData.resize(vtxIndices.size());
0904     for (unsigned int idx = 0; idx < vtxIndices.size(); idx++) {
0905       const SecondaryVertex &sv = SVs[vtxIndices[idx]];
0906 
0907       svData[idx].vertex = sv;
0908       svData[idx].dist1d = sv.dist1d();
0909       svData[idx].dist2d = sv.dist2d();
0910       svData[idx].dist3d = sv.dist3d();
0911       svData[idx].direction = flightDirection(pv, sv);
0912       // mark tracks successfully used in vertex fit
0913       markUsedTracks(trackData, trackRefs, sv, idx);
0914     }
0915 
0916     // fill result into tag infos
0917 
0918     tagInfos->push_back(TemplatedSecondaryVertexTagInfo<IPTI, VTX>(
0919         trackData,
0920         svData,
0921         SVs.size(),
0922         edm::Ref<std::vector<IPTI> >(trackIPTagInfos, iterJets - trackIPTagInfos->begin())));
0923   }
0924 
0925   event.put(std::move(tagInfos));
0926 }
0927 
0928 //Need specialized template because reco::Vertex iterators are TrackBase and it is a mess to make general
0929 template <>
0930 void TemplatedSecondaryVertexProducer<TrackIPTagInfo, reco::Vertex>::markUsedTracks(TrackDataVector &trackData,
0931                                                                                     const input_container &trackRefs,
0932                                                                                     const SecondaryVertex &sv,
0933                                                                                     size_t idx) {
0934   for (Vertex::trackRef_iterator iter = sv.tracks_begin(); iter != sv.tracks_end(); ++iter) {
0935     if (sv.trackWeight(*iter) < minTrackWeight)
0936       continue;
0937 
0938     typename input_container::const_iterator pos =
0939         std::find(trackRefs.begin(), trackRefs.end(), iter->castTo<input_item>());
0940 
0941     if (pos == trackRefs.end()) {
0942       if (!useExternalSV)
0943         throw cms::Exception("TrackNotFound") << "Could not find track from secondary "
0944                                                  "vertex in original tracks."
0945                                               << std::endl;
0946     } else {
0947       unsigned int index = pos - trackRefs.begin();
0948       trackData[index].second.svStatus = (btag::TrackData::trackAssociatedToVertex + idx);
0949     }
0950   }
0951 }
0952 template <>
0953 void TemplatedSecondaryVertexProducer<CandIPTagInfo, reco::VertexCompositePtrCandidate>::markUsedTracks(
0954     TrackDataVector &trackData, const input_container &trackRefs, const SecondaryVertex &sv, size_t idx) {
0955   for (typename input_container::const_iterator iter = sv.daughterPtrVector().begin();
0956        iter != sv.daughterPtrVector().end();
0957        ++iter) {
0958     typename input_container::const_iterator pos = std::find(trackRefs.begin(), trackRefs.end(), *iter);
0959 
0960     if (pos != trackRefs.end()) {
0961       unsigned int index = pos - trackRefs.begin();
0962       trackData[index].second.svStatus = (btag::TrackData::trackAssociatedToVertex + idx);
0963     }
0964   }
0965 }
0966 
0967 template <>
0968 typename TemplatedSecondaryVertexProducer<TrackIPTagInfo, reco::Vertex>::SecondaryVertex
0969 TemplatedSecondaryVertexProducer<TrackIPTagInfo, reco::Vertex>::SVBuilder::operator()(const TransientVertex &sv) const {
0970   if (!sv.originalTracks().empty() && sv.originalTracks()[0].trackBaseRef().isNonnull())
0971     return SecondaryVertex(pv, sv, direction, withPVError);
0972   else {
0973     edm::LogError("UnexpectedInputs") << "Building from Candidates, should not happen!";
0974     return SecondaryVertex(pv, sv, direction, withPVError);
0975   }
0976 }
0977 
0978 template <>
0979 typename TemplatedSecondaryVertexProducer<CandIPTagInfo, reco::VertexCompositePtrCandidate>::SecondaryVertex
0980 TemplatedSecondaryVertexProducer<CandIPTagInfo, reco::VertexCompositePtrCandidate>::SVBuilder::operator()(
0981     const TransientVertex &sv) const {
0982   if (!sv.originalTracks().empty() && sv.originalTracks()[0].trackBaseRef().isNonnull()) {
0983     edm::LogError("UnexpectedInputs") << "Building from Tracks, should not happen!";
0984     VertexCompositePtrCandidate vtxCompPtrCand;
0985 
0986     vtxCompPtrCand.setCovariance(sv.vertexState().error().matrix());
0987     vtxCompPtrCand.setChi2AndNdof(sv.totalChiSquared(), sv.degreesOfFreedom());
0988     vtxCompPtrCand.setVertex(Candidate::Point(sv.position().x(), sv.position().y(), sv.position().z()));
0989 
0990     return SecondaryVertex(pv, vtxCompPtrCand, direction, withPVError);
0991   } else {
0992     VertexCompositePtrCandidate vtxCompPtrCand;
0993 
0994     vtxCompPtrCand.setCovariance(sv.vertexState().error().matrix());
0995     vtxCompPtrCand.setChi2AndNdof(sv.totalChiSquared(), sv.degreesOfFreedom());
0996     vtxCompPtrCand.setVertex(Candidate::Point(sv.position().x(), sv.position().y(), sv.position().z()));
0997 
0998     Candidate::LorentzVector p4;
0999     for (std::vector<reco::TransientTrack>::const_iterator tt = sv.originalTracks().begin();
1000          tt != sv.originalTracks().end();
1001          ++tt) {
1002       if (sv.trackWeight(*tt) < minTrackWeight)
1003         continue;
1004 
1005       const CandidatePtrTransientTrack *cptt =
1006           dynamic_cast<const CandidatePtrTransientTrack *>(tt->basicTransientTrack());
1007       if (cptt == nullptr)
1008         edm::LogError("DynamicCastingFailed") << "Casting of TransientTrack to CandidatePtrTransientTrack failed!";
1009       else {
1010         p4 += cptt->candidate()->p4();
1011         vtxCompPtrCand.addDaughter(cptt->candidate());
1012       }
1013     }
1014     vtxCompPtrCand.setP4(p4);
1015 
1016     return SecondaryVertex(pv, vtxCompPtrCand, direction, withPVError);
1017   }
1018 }
1019 
1020 // ------------ method that matches reclustered and original jets based on minimum dR ------------
1021 template <class IPTI, class VTX>
1022 template <class CONTAINER>
1023 void TemplatedSecondaryVertexProducer<IPTI, VTX>::matchReclusteredJets(
1024     const edm::Handle<CONTAINER> &jets,
1025     const std::vector<fastjet::PseudoJet> &reclusteredJets,
1026     std::vector<int> &matchedIndices,
1027     const std::string &jetType) {
1028   std::string type = (!jetType.empty() ? jetType + " " : jetType);
1029 
1030   std::vector<bool> matchedLocks(reclusteredJets.size(), false);
1031 
1032   for (size_t j = 0; j < jets->size(); ++j) {
1033     double matchedDR2 = 1e9;
1034     int matchedIdx = -1;
1035 
1036     for (size_t rj = 0; rj < reclusteredJets.size(); ++rj) {
1037       if (matchedLocks.at(rj))
1038         continue;  // skip jets that have already been matched
1039 
1040       double tempDR2 = Geom::deltaR2(toJet(jets->at(j))->rapidity(),
1041                                      toJet(jets->at(j))->phi(),
1042                                      reclusteredJets.at(rj).rapidity(),
1043                                      reclusteredJets.at(rj).phi_std());
1044       if (tempDR2 < matchedDR2) {
1045         matchedDR2 = tempDR2;
1046         matchedIdx = rj;
1047       }
1048     }
1049 
1050     if (matchedIdx >= 0) {
1051       if (matchedDR2 > rParam * rParam) {
1052         edm::LogError("JetMatchingFailed") << "Matched reclustered jet " << matchedIdx << " and original " << type
1053                                            << "jet " << j << " are separated by dR=" << sqrt(matchedDR2)
1054                                            << " which is greater than the jet size R=" << rParam << ".\n"
1055                                            << "This is not expected so please check that the jet algorithm and jet "
1056                                               "size match those used for the original "
1057                                            << type << "jet collection.";
1058       } else
1059         matchedLocks.at(matchedIdx) = true;
1060     } else
1061       edm::LogError("JetMatchingFailed")
1062           << "Matching reclustered to original " << type
1063           << "jets failed. Please check that the jet algorithm and jet size match those used for the original " << type
1064           << "jet collection.";
1065 
1066     matchedIndices.push_back(matchedIdx);
1067   }
1068 }
1069 
1070 // ------------ method that matches groomed and original jets based on minimum dR ------------
1071 template <class IPTI, class VTX>
1072 void TemplatedSecondaryVertexProducer<IPTI, VTX>::matchGroomedJets(const edm::Handle<edm::View<reco::Jet> > &jets,
1073                                                                    const edm::Handle<edm::View<reco::Jet> > &groomedJets,
1074                                                                    std::vector<int> &matchedIndices) {
1075   std::vector<bool> jetLocks(jets->size(), false);
1076   std::vector<int> jetIndices;
1077 
1078   for (size_t gj = 0; gj < groomedJets->size(); ++gj) {
1079     double matchedDR2 = 1e9;
1080     int matchedIdx = -1;
1081 
1082     if (groomedJets->at(gj).pt() > 0.)  // skip pathological cases of groomed jets with Pt=0
1083     {
1084       for (size_t j = 0; j < jets->size(); ++j) {
1085         if (jetLocks.at(j))
1086           continue;  // skip jets that have already been matched
1087 
1088         double tempDR2 = Geom::deltaR2(
1089             jets->at(j).rapidity(), jets->at(j).phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi());
1090         if (tempDR2 < matchedDR2) {
1091           matchedDR2 = tempDR2;
1092           matchedIdx = j;
1093         }
1094       }
1095     }
1096 
1097     if (matchedIdx >= 0) {
1098       if (matchedDR2 > rParam * rParam) {
1099         edm::LogWarning("MatchedJetsFarApart")
1100             << "Matched groomed jet " << gj << " and original jet " << matchedIdx
1101             << " are separated by dR=" << sqrt(matchedDR2) << " which is greater than the jet size R=" << rParam
1102             << ".\n"
1103             << "This is not expected so the matching of these two jets has been discarded. Please check that the two "
1104                "jet collections belong to each other.";
1105         matchedIdx = -1;
1106       } else
1107         jetLocks.at(matchedIdx) = true;
1108     }
1109     jetIndices.push_back(matchedIdx);
1110   }
1111 
1112   for (size_t j = 0; j < jets->size(); ++j) {
1113     std::vector<int>::iterator matchedIndex = std::find(jetIndices.begin(), jetIndices.end(), j);
1114 
1115     matchedIndices.push_back(matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(), matchedIndex) : -1);
1116   }
1117 }
1118 
1119 // ------------ method that matches subjets and original fat jets ------------
1120 template <class IPTI, class VTX>
1121 void TemplatedSecondaryVertexProducer<IPTI, VTX>::matchSubjets(const std::vector<int> &groomedIndices,
1122                                                                const edm::Handle<edm::View<reco::Jet> > &groomedJets,
1123                                                                const edm::Handle<std::vector<IPTI> > &subjets,
1124                                                                std::vector<std::vector<int> > &matchedIndices) {
1125   for (size_t g = 0; g < groomedIndices.size(); ++g) {
1126     std::vector<int> subjetIndices;
1127 
1128     if (groomedIndices.at(g) >= 0) {
1129       for (size_t s = 0; s < groomedJets->at(groomedIndices.at(g)).numberOfDaughters(); ++s) {
1130         const edm::Ptr<reco::Candidate> &subjet = groomedJets->at(groomedIndices.at(g)).daughterPtr(s);
1131 
1132         for (size_t sj = 0; sj < subjets->size(); ++sj) {
1133           const edm::RefToBase<reco::Jet> &subjetRef = subjets->at(sj).jet();
1134           if (subjet == edm::Ptr<reco::Candidate>(subjetRef.id(), subjetRef.get(), subjetRef.key())) {
1135             subjetIndices.push_back(sj);
1136             break;
1137           }
1138         }
1139       }
1140 
1141       if (subjetIndices.empty())
1142         edm::LogError("SubjetMatchingFailed") << "Matching subjets to original fat jets failed. Please check that the "
1143                                                  "groomed fat jet and subjet collections belong to each other.";
1144 
1145       matchedIndices.push_back(subjetIndices);
1146     } else
1147       matchedIndices.push_back(subjetIndices);
1148   }
1149 }
1150 
1151 // ------------ method that matches subjets and original fat jets ------------
1152 template <class IPTI, class VTX>
1153 void TemplatedSecondaryVertexProducer<IPTI, VTX>::matchSubjets(const edm::Handle<edm::View<reco::Jet> > &fatJets,
1154                                                                const edm::Handle<std::vector<IPTI> > &subjets,
1155                                                                std::vector<std::vector<int> > &matchedIndices) {
1156   for (size_t fj = 0; fj < fatJets->size(); ++fj) {
1157     std::vector<int> subjetIndices;
1158     size_t nSubjetCollections = 0;
1159     size_t nSubjets = 0;
1160 
1161     const pat::Jet *fatJet = dynamic_cast<const pat::Jet *>(fatJets->ptrAt(fj).get());
1162 
1163     if (!fatJet) {
1164       if (fj == 0)
1165         edm::LogError("WrongJetType")
1166             << "Wrong jet type for input fat jets. Please check that the input fat jets are of the pat::Jet type.";
1167 
1168       matchedIndices.push_back(subjetIndices);
1169       continue;
1170     } else {
1171       nSubjetCollections = fatJet->subjetCollectionNames().size();
1172 
1173       if (nSubjetCollections > 0) {
1174         for (size_t coll = 0; coll < nSubjetCollections; ++coll) {
1175           const pat::JetPtrCollection &fatJetSubjets = fatJet->subjets(coll);
1176 
1177           for (size_t fjsj = 0; fjsj < fatJetSubjets.size(); ++fjsj) {
1178             ++nSubjets;
1179 
1180             for (size_t sj = 0; sj < subjets->size(); ++sj) {
1181               const pat::Jet *subJet = dynamic_cast<const pat::Jet *>(subjets->at(sj).jet().get());
1182 
1183               if (!subJet) {
1184                 if (fj == 0 && coll == 0 && fjsj == 0 && sj == 0)
1185                   edm::LogError("WrongJetType") << "Wrong jet type for input subjets. Please check that the input "
1186                                                    "subjets are of the pat::Jet type.";
1187 
1188                 break;
1189               } else {
1190                 if (subJet->originalObjectRef() == fatJetSubjets.at(fjsj)->originalObjectRef()) {
1191                   subjetIndices.push_back(sj);
1192                   break;
1193                 }
1194               }
1195             }
1196           }
1197         }
1198 
1199         if (subjetIndices.empty() && nSubjets > 0)
1200           edm::LogError("SubjetMatchingFailed") << "Matching subjets to fat jets failed. Please check that the fat jet "
1201                                                    "and subjet collections belong to each other.";
1202 
1203         matchedIndices.push_back(subjetIndices);
1204       } else
1205         matchedIndices.push_back(subjetIndices);
1206     }
1207   }
1208 }
1209 
1210 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
1211 template <class IPTI, class VTX>
1212 void TemplatedSecondaryVertexProducer<IPTI, VTX>::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
1213   edm::ParameterSetDescription desc;
1214   desc.add<double>("extSVDeltaRToJet", 0.3);
1215   desc.add<edm::InputTag>("beamSpotTag", edm::InputTag("offlineBeamSpot"));
1216   {
1217     edm::ParameterSetDescription vertexReco;
1218     vertexReco.add<double>("primcut", 1.8);
1219     vertexReco.add<double>("seccut", 6.0);
1220     vertexReco.add<std::string>("finder", "avr");
1221     vertexReco.addOptionalNode(edm::ParameterDescription<double>("minweight", 0.5, true) and
1222                                    edm::ParameterDescription<double>("weightthreshold", 0.001, true) and
1223                                    edm::ParameterDescription<bool>("smoothing", false, true),
1224                                true);
1225     vertexReco.addOptionalNode(
1226         edm::ParameterDescription<double>("maxFitChi2", 10.0, true) and
1227             edm::ParameterDescription<double>("mergeThreshold", 3.0, true) and
1228             edm::ParameterDescription<std::string>("fitType", "RefitGhostTrackWithVertices", true),
1229         true);
1230     desc.add<edm::ParameterSetDescription>("vertexReco", vertexReco);
1231   }
1232   {
1233     edm::ParameterSetDescription vertexSelection;
1234     vertexSelection.add<std::string>("sortCriterium", "dist3dError");
1235     desc.add<edm::ParameterSetDescription>("vertexSelection", vertexSelection);
1236   }
1237   desc.add<std::string>("constraint", "BeamSpot");
1238   desc.add<edm::InputTag>("trackIPTagInfos", edm::InputTag("impactParameterTagInfos"));
1239   {
1240     edm::ParameterSetDescription vertexCuts;
1241     vertexCuts.add<double>("distSig3dMax", 99999.9);
1242     vertexCuts.add<double>("fracPV", 0.65);
1243     vertexCuts.add<double>("distVal2dMax", 2.5);
1244     vertexCuts.add<bool>("useTrackWeights", true);
1245     vertexCuts.add<double>("maxDeltaRToJetAxis", 0.4);
1246     {
1247       edm::ParameterSetDescription v0Filter;
1248       v0Filter.add<double>("k0sMassWindow", 0.05);
1249       vertexCuts.add<edm::ParameterSetDescription>("v0Filter", v0Filter);
1250     }
1251     vertexCuts.add<double>("distSig2dMin", 3.0);
1252     vertexCuts.add<unsigned int>("multiplicityMin", 2);
1253     vertexCuts.add<double>("distVal2dMin", 0.01);
1254     vertexCuts.add<double>("distSig2dMax", 99999.9);
1255     vertexCuts.add<double>("distVal3dMax", 99999.9);
1256     vertexCuts.add<double>("minimumTrackWeight", 0.5);
1257     vertexCuts.add<double>("distVal3dMin", -99999.9);
1258     vertexCuts.add<double>("massMax", 6.5);
1259     vertexCuts.add<double>("distSig3dMin", -99999.9);
1260     desc.add<edm::ParameterSetDescription>("vertexCuts", vertexCuts);
1261   }
1262   desc.add<bool>("useExternalSV", false);
1263   desc.add<double>("minimumTrackWeight", 0.5);
1264   desc.add<bool>("usePVError", true);
1265   {
1266     edm::ParameterSetDescription trackSelection;
1267     trackSelection.add<double>("b_pT", 0.3684);
1268     trackSelection.add<double>("max_pT", 500);
1269     trackSelection.add<bool>("useVariableJTA", false);
1270     trackSelection.add<double>("maxDecayLen", 99999.9);
1271     trackSelection.add<double>("sip3dValMin", -99999.9);
1272     trackSelection.add<double>("max_pT_dRcut", 0.1);
1273     trackSelection.add<double>("a_pT", 0.005263);
1274     trackSelection.add<unsigned int>("totalHitsMin", 8);
1275     trackSelection.add<double>("jetDeltaRMax", 0.3);
1276     trackSelection.add<double>("a_dR", -0.001053);
1277     trackSelection.add<double>("maxDistToAxis", 0.2);
1278     trackSelection.add<double>("ptMin", 1.0);
1279     trackSelection.add<std::string>("qualityClass", "any");
1280     trackSelection.add<unsigned int>("pixelHitsMin", 2);
1281     trackSelection.add<double>("sip2dValMax", 99999.9);
1282     trackSelection.add<double>("max_pT_trackPTcut", 3);
1283     trackSelection.add<double>("sip2dValMin", -99999.9);
1284     trackSelection.add<double>("normChi2Max", 99999.9);
1285     trackSelection.add<double>("sip3dValMax", 99999.9);
1286     trackSelection.add<double>("sip3dSigMin", -99999.9);
1287     trackSelection.add<double>("min_pT", 120);
1288     trackSelection.add<double>("min_pT_dRcut", 0.5);
1289     trackSelection.add<double>("sip2dSigMax", 99999.9);
1290     trackSelection.add<double>("sip3dSigMax", 99999.9);
1291     trackSelection.add<double>("sip2dSigMin", -99999.9);
1292     trackSelection.add<double>("b_dR", 0.6263);
1293     desc.add<edm::ParameterSetDescription>("trackSelection", trackSelection);
1294   }
1295   desc.add<std::string>("trackSort", "sip3dSig");
1296   desc.add<edm::InputTag>("extSVCollection", edm::InputTag("secondaryVertices"));
1297   desc.addOptionalNode(edm::ParameterDescription<bool>("useSVClustering", false, true) and
1298                            edm::ParameterDescription<std::string>("jetAlgorithm", true) and
1299                            edm::ParameterDescription<double>("rParam", true),
1300                        true);
1301   desc.addOptional<bool>("useSVMomentum", false);
1302   desc.addOptional<double>("ghostRescaling", 1e-18);
1303   desc.addOptional<double>("relPtTolerance", 1e-03);
1304   desc.addOptional<edm::InputTag>("fatJets");
1305   desc.addOptional<edm::InputTag>("groomedFatJets");
1306   desc.add<edm::InputTag>("weights", edm::InputTag(""));
1307   descriptions.addDefault(desc);
1308 }
1309 
1310 //define this as a plug-in
1311 typedef TemplatedSecondaryVertexProducer<TrackIPTagInfo, reco::Vertex> SecondaryVertexProducer;
1312 typedef TemplatedSecondaryVertexProducer<CandIPTagInfo, reco::VertexCompositePtrCandidate> CandSecondaryVertexProducer;
1313 
1314 DEFINE_FWK_MODULE(SecondaryVertexProducer);
1315 DEFINE_FWK_MODULE(CandSecondaryVertexProducer);