TemplatedSecondaryVertexProducer.cc

Go to the documentation of this file.

#include <functional>
#include <algorithm>
#include <iterator>
#include <cstddef>
#include <string>
#include <vector>
#include <map>
#include <set>

#include <boost/iterator/transform_iterator.hpp>

#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/ParameterSet/interface/IfExistsDescription.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"

#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"

#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
#include "DataFormats/BTauReco/interface/TrackIPTagInfo.h"
#include "DataFormats/BTauReco/interface/CandIPTagInfo.h"
#include "DataFormats/BTauReco/interface/SecondaryVertexTagInfo.h"

#include "RecoVertex/VertexPrimitives/interface/VertexException.h"
#include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
#include "RecoVertex/ConfigurableVertexReco/interface/ConfigurableVertexReconstructor.h"
#include "RecoVertex/GhostTrackFitter/interface/GhostTrackVertexFinder.h"
#include "RecoVertex/GhostTrackFitter/interface/GhostTrackPrediction.h"
#include "RecoVertex/GhostTrackFitter/interface/GhostTrackState.h"
#include "RecoVertex/GhostTrackFitter/interface/GhostTrack.h"

#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/TransientTrack/interface/CandidatePtrTransientTrack.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"

#include "RecoBTag/SecondaryVertex/interface/TrackSelector.h"
#include "RecoBTag/SecondaryVertex/interface/TrackSorting.h"
#include "RecoBTag/SecondaryVertex/interface/SecondaryVertex.h"
#include "RecoBTag/SecondaryVertex/interface/VertexFilter.h"
#include "RecoBTag/SecondaryVertex/interface/VertexSorting.h"

#include "DataFormats/GeometryVector/interface/VectorUtil.h"

#include "fastjet/JetDefinition.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/PseudoJet.hh"

//
// constants, enums and typedefs
//
typedef boost::shared_ptr<fastjet::ClusterSequence>  ClusterSequencePtr;
typedef boost::shared_ptr<fastjet::JetDefinition>    JetDefPtr;


using namespace reco;

namespace {
    class VertexInfo : public fastjet::PseudoJet::UserInfoBase{
      public:
        VertexInfo(const int vertexIndex) :
          m_vertexIndex(vertexIndex) { }

        inline const int vertexIndex() const { return m_vertexIndex; }

      protected:
        int m_vertexIndex;
    };
    
    template<typename T>
    struct RefToBaseLess : public std::binary_function<edm::RefToBase<T>,
                               edm::RefToBase<T>,
                               bool> {
        inline bool operator()(const edm::RefToBase<T> &r1,
                       const edm::RefToBase<T> &r2) const
        {
            return r1.id() < r2.id() ||
                   (r1.id() == r2.id() && r1.key() < r2.key());
        }
    };
}

GlobalVector flightDirection(const reco::Vertex & pv, const reco::Vertex & sv) {
return  GlobalVector(sv.x() - pv.x(), sv.y() - pv.y(),sv.z() - pv.z());
}
GlobalVector flightDirection(const reco::Vertex & pv, const reco::VertexCompositePtrCandidate & sv) {
return  GlobalVector(sv.vertex().x() - pv.x(), sv.vertex().y() - pv.y(),sv.vertex().z() - pv.z());
}
const math::XYZPoint & position(const reco::Vertex & sv)
{return sv.position();}
const math::XYZPoint & position(const reco::VertexCompositePtrCandidate & sv)
{return sv.vertex();}


template <class IPTI,class VTX>
class TemplatedSecondaryVertexProducer : public edm::stream::EDProducer<> {
    public:
    explicit TemplatedSecondaryVertexProducer(const edm::ParameterSet &params);
    ~TemplatedSecondaryVertexProducer() override;
    static void fillDescriptions(edm::ConfigurationDescriptions & descriptions);
    typedef std::vector<TemplatedSecondaryVertexTagInfo<IPTI,VTX> > Product;
    typedef TemplatedSecondaryVertex<VTX> SecondaryVertex;
    typedef typename IPTI::input_container input_container;
    typedef typename IPTI::input_container::value_type input_item;
    typedef typename std::vector<reco::btag::IndexedTrackData> TrackDataVector;
    void produce(edm::Event &event, const edm::EventSetup &es) override;

    private:
        template<class CONTAINER>
    void matchReclusteredJets(const edm::Handle<CONTAINER>& jets,
                  const std::vector<fastjet::PseudoJet>& matchedJets,
                  std::vector<int>& matchedIndices,
                  const std::string& jetType="");
    void matchGroomedJets(const edm::Handle<edm::View<reco::Jet> >& jets,
                  const edm::Handle<edm::View<reco::Jet> >& matchedJets,
                  std::vector<int>& matchedIndices);
    void matchSubjets(const std::vector<int>& groomedIndices,
              const edm::Handle<edm::View<reco::Jet> >& groomedJets,
              const edm::Handle<std::vector<IPTI> >& subjets,
              std::vector<std::vector<int> >& matchedIndices);
    void matchSubjets(const edm::Handle<edm::View<reco::Jet> >& fatJets,
              const edm::Handle<std::vector<IPTI> >& subjets,
              std::vector<std::vector<int> >& matchedIndices);
    
    const reco::Jet * toJet(const reco::Jet & j) { return &j; }
    const reco::Jet * toJet(const IPTI & j) { return &(*(j.jet())); }

    enum ConstraintType {
        CONSTRAINT_NONE = 0,
        CONSTRAINT_BEAMSPOT,
        CONSTRAINT_PV_BEAMSPOT_SIZE,
        CONSTRAINT_PV_BS_Z_ERRORS_SCALED,
        CONSTRAINT_PV_ERROR_SCALED,
        CONSTRAINT_PV_PRIMARIES_IN_FIT
    };
    static ConstraintType getConstraintType(const std::string &name);

        edm::EDGetTokenT<reco::BeamSpot> token_BeamSpot;
        edm::EDGetTokenT<std::vector<IPTI> > token_trackIPTagInfo;
    reco::btag::SortCriteria    sortCriterium;
    TrackSelector           trackSelector;
    ConstraintType          constraint;
    double              constraintScaling;
    edm::ParameterSet       vtxRecoPSet;
    bool                useGhostTrack;
    bool                withPVError;
    double              minTrackWeight;
    VertexFilter            vertexFilter;
    VertexSorting<SecondaryVertex>  vertexSorting;
        bool                            useExternalSV;  
        double                          extSVDeltaRToJet;      
        edm::EDGetTokenT<edm::View<VTX> > token_extSVCollection;
    bool                useSVClustering;
    bool                useSVMomentum;
    std::string         jetAlgorithm;
    double              rParam;
    double              jetPtMin;
    double              ghostRescaling;
    double              relPtTolerance;
    bool                useFatJets;
    bool                useGroomedFatJets;
    edm::EDGetTokenT<edm::View<reco::Jet> > token_fatJets;
    edm::EDGetTokenT<edm::View<reco::Jet> > token_groomedFatJets;
    
    ClusterSequencePtr      fjClusterSeq;
    JetDefPtr           fjJetDefinition;

    void markUsedTracks(TrackDataVector & trackData, const input_container & trackRefs, const SecondaryVertex & sv,size_t idx);

    struct SVBuilder :
        public std::unary_function<const VTX&, SecondaryVertex> {

        SVBuilder(const reco::Vertex &pv,
                  const GlobalVector &direction,
                  bool withPVError,
              double minTrackWeight) :
            pv(pv), direction(direction),
            withPVError(withPVError),
            minTrackWeight(minTrackWeight) {}
        SecondaryVertex operator () (const TransientVertex &sv) const;

        SecondaryVertex operator () (const VTX &sv) const
        { return SecondaryVertex(pv, sv, direction, withPVError); }


        const Vertex        &pv;
        const GlobalVector  &direction;
        bool            withPVError;
        double          minTrackWeight;
    };

    struct SVFilter :
        public std::unary_function<const SecondaryVertex&, bool> {

        SVFilter(const VertexFilter &filter, const Vertex &pv,
                 const GlobalVector &direction) :
            filter(filter), pv(pv), direction(direction) {}

        inline bool operator () (const SecondaryVertex &sv) const
        { return !filter(pv, sv, direction); }

        const VertexFilter  &filter;
        const Vertex        &pv;
        const GlobalVector  &direction;
    };
};
template <class IPTI,class VTX>
typename TemplatedSecondaryVertexProducer<IPTI,VTX>::ConstraintType
TemplatedSecondaryVertexProducer<IPTI,VTX>::getConstraintType(const std::string &name)
{
    if (name == "None")
        return CONSTRAINT_NONE;
    else if (name == "BeamSpot")
        return CONSTRAINT_BEAMSPOT;
    else if (name == "BeamSpot+PVPosition")
        return CONSTRAINT_PV_BEAMSPOT_SIZE;
    else if (name == "BeamSpotZ+PVErrorScaledXY")
        return CONSTRAINT_PV_BS_Z_ERRORS_SCALED;
    else if (name == "PVErrorScaled")
        return CONSTRAINT_PV_ERROR_SCALED;
    else if (name == "BeamSpot+PVTracksInFit")
        return CONSTRAINT_PV_PRIMARIES_IN_FIT;
    else
        throw cms::Exception("InvalidArgument")
            << "TemplatedSecondaryVertexProducer: ``constraint'' parameter "
               "value \"" << name << "\" not understood."
            << std::endl;
}

static GhostTrackVertexFinder::FitType
getGhostTrackFitType(const std::string &name)
{
    if (name == "AlwaysWithGhostTrack")
        return GhostTrackVertexFinder::kAlwaysWithGhostTrack;
    else if (name == "SingleTracksWithGhostTrack")
        return GhostTrackVertexFinder::kSingleTracksWithGhostTrack;
    else if (name == "RefitGhostTrackWithVertices")
        return GhostTrackVertexFinder::kRefitGhostTrackWithVertices;
    else
        throw cms::Exception("InvalidArgument")
            << "TemplatedSecondaryVertexProducer: ``fitType'' "
               "parameter value \"" << name << "\" for "
               "GhostTrackVertexFinder settings not "
               "understood." << std::endl;
}

template <class IPTI,class VTX>
TemplatedSecondaryVertexProducer<IPTI,VTX>::TemplatedSecondaryVertexProducer(
                    const edm::ParameterSet &params) :
    sortCriterium(TrackSorting::getCriterium(params.getParameter<std::string>("trackSort"))),
    trackSelector(params.getParameter<edm::ParameterSet>("trackSelection")),
    constraint(getConstraintType(params.getParameter<std::string>("constraint"))),
    constraintScaling(1.0),
    vtxRecoPSet(params.getParameter<edm::ParameterSet>("vertexReco")),
    useGhostTrack(vtxRecoPSet.getParameter<std::string>("finder") == "gtvr"),
    withPVError(params.getParameter<bool>("usePVError")),
    minTrackWeight(params.getParameter<double>("minimumTrackWeight")),
    vertexFilter(params.getParameter<edm::ParameterSet>("vertexCuts")),
    vertexSorting(params.getParameter<edm::ParameterSet>("vertexSelection"))
{
    token_trackIPTagInfo =  consumes<std::vector<IPTI> >(params.getParameter<edm::InputTag>("trackIPTagInfos"));
    if (constraint == CONSTRAINT_PV_ERROR_SCALED ||
        constraint == CONSTRAINT_PV_BS_Z_ERRORS_SCALED)
        constraintScaling = params.getParameter<double>("pvErrorScaling");

    if (constraint == CONSTRAINT_PV_BEAMSPOT_SIZE ||
        constraint == CONSTRAINT_PV_BS_Z_ERRORS_SCALED ||
        constraint == CONSTRAINT_BEAMSPOT ||
        constraint == CONSTRAINT_PV_PRIMARIES_IN_FIT )
        token_BeamSpot = consumes<reco::BeamSpot>(params.getParameter<edm::InputTag>("beamSpotTag"));
    useExternalSV = false;
    if(params.existsAs<bool>("useExternalSV")) useExternalSV = params.getParameter<bool> ("useExternalSV");
    if(useExternalSV) {
       token_extSVCollection =  consumes<edm::View<VTX> >(params.getParameter<edm::InputTag>("extSVCollection"));
       extSVDeltaRToJet = params.getParameter<double>("extSVDeltaRToJet");
    }
    useSVClustering = ( params.existsAs<bool>("useSVClustering") ? params.getParameter<bool>("useSVClustering") : false );
    useSVMomentum = ( params.existsAs<bool>("useSVMomentum") ? params.getParameter<bool>("useSVMomentum") : false );
    useFatJets = ( useExternalSV && params.exists("fatJets") );
    useGroomedFatJets = ( useExternalSV && params.exists("groomedFatJets") );
    if( useSVClustering )
    {
      jetAlgorithm = params.getParameter<std::string>("jetAlgorithm");
      rParam = params.getParameter<double>("rParam");
      jetPtMin = 0.; // hardcoded to 0. since we simply want to recluster all input jets which already had some PtMin applied
      ghostRescaling = ( params.existsAs<double>("ghostRescaling") ? params.getParameter<double>("ghostRescaling") : 1e-18 );
      relPtTolerance = ( params.existsAs<double>("relPtTolerance") ? params.getParameter<double>("relPtTolerance") : 1e-03); // 0.1% relative difference in Pt should be sufficient to detect possible misconfigurations

      // set jet algorithm
      if (jetAlgorithm=="Kt")
        fjJetDefinition = JetDefPtr( new fastjet::JetDefinition(fastjet::kt_algorithm, rParam) );
      else if (jetAlgorithm=="CambridgeAachen")
        fjJetDefinition = JetDefPtr( new fastjet::JetDefinition(fastjet::cambridge_algorithm, rParam) );
      else if (jetAlgorithm=="AntiKt")
        fjJetDefinition = JetDefPtr( new fastjet::JetDefinition(fastjet::antikt_algorithm, rParam) );
      else
        throw cms::Exception("InvalidJetAlgorithm") << "Jet clustering algorithm is invalid: " << jetAlgorithm << ", use CambridgeAachen | Kt | AntiKt" << std::endl;
    }
    if( useFatJets )
      token_fatJets = consumes<edm::View<reco::Jet> >(params.getParameter<edm::InputTag>("fatJets"));
    if( useGroomedFatJets )
      token_groomedFatJets = consumes<edm::View<reco::Jet> >(params.getParameter<edm::InputTag>("groomedFatJets"));
    if( useFatJets && !useSVClustering )
      rParam = params.getParameter<double>("rParam"); // will be used later as a dR cut
    
    produces<Product>();
}
template <class IPTI,class VTX>
TemplatedSecondaryVertexProducer<IPTI,VTX>::~TemplatedSecondaryVertexProducer()
{
}

template <class IPTI,class VTX>
void TemplatedSecondaryVertexProducer<IPTI,VTX>::produce(edm::Event &event,
                                      const edm::EventSetup &es)
{
//  typedef std::map<TrackBaseRef, TransientTrack,
//                   RefToBaseLess<Track> > TransientTrackMap;
    //How about good old pointers?
    typedef std::map<const Track *, TransientTrack> TransientTrackMap;
    

    edm::ESHandle<TransientTrackBuilder> trackBuilder;
    es.get<TransientTrackRecord>().get("TransientTrackBuilder",
                                       trackBuilder);

    edm::Handle<std::vector<IPTI> > trackIPTagInfos;
    event.getByToken(token_trackIPTagInfo, trackIPTagInfos);

        // External Sec Vertex collection (e.g. for IVF usage)
        edm::Handle<edm::View<VTX> > extSecVertex;          
        if(useExternalSV) event.getByToken(token_extSVCollection,extSecVertex);

    edm::Handle<edm::View<reco::Jet> > fatJetsHandle;
    edm::Handle<edm::View<reco::Jet> > groomedFatJetsHandle;
    if( useFatJets )
    {
      event.getByToken(token_fatJets, fatJetsHandle);
    
      if( useGroomedFatJets )
      {
        event.getByToken(token_groomedFatJets, groomedFatJetsHandle);
    
        if( groomedFatJetsHandle->size() > fatJetsHandle->size() )
          edm::LogError("TooManyGroomedJets") << "There are more groomed (" << groomedFatJetsHandle->size() << ") than original fat jets (" << fatJetsHandle->size() << "). Please check that the two jet collections belong to each other.";
      }
    }

    edm::Handle<BeamSpot> beamSpot;
    unsigned int bsCovSrc[7] = { 0, };
    double sigmaZ = 0.0, beamWidth = 0.0;
    switch(constraint) {
        case CONSTRAINT_PV_BEAMSPOT_SIZE:
        event.getByToken(token_BeamSpot,beamSpot);
        bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = bsCovSrc[6] = 1;
        sigmaZ = beamSpot->sigmaZ();
        beamWidth = beamSpot->BeamWidthX();
        break;

        case CONSTRAINT_PV_BS_Z_ERRORS_SCALED:
        event.getByToken(token_BeamSpot,beamSpot);
        bsCovSrc[0] = bsCovSrc[1] = 2;
        bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = 1;
        sigmaZ = beamSpot->sigmaZ();
        break;

        case CONSTRAINT_PV_ERROR_SCALED:
        bsCovSrc[0] = bsCovSrc[1] = bsCovSrc[2] = 2;
        break;

        case CONSTRAINT_BEAMSPOT:
        case CONSTRAINT_PV_PRIMARIES_IN_FIT:
        event.getByToken(token_BeamSpot,beamSpot);
        break;

        default:
        /* nothing */;
    }

    // ------------------------------------ SV clustering START --------------------------------------------
    std::vector<std::vector<int> > clusteredSVs(trackIPTagInfos->size(),std::vector<int>());
    if( useExternalSV && useSVClustering && !trackIPTagInfos->empty() )
    {
      // vector of constituents for reclustering jets and "ghost" SVs
      std::vector<fastjet::PseudoJet> fjInputs;
      // loop over all input jets and collect all their constituents
      if( useFatJets )
      {
        for(edm::View<reco::Jet>::const_iterator it = fatJetsHandle->begin(); it != fatJetsHandle->end(); ++it)
        {
          std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
          std::vector<edm::Ptr<reco::Candidate> >::const_iterator m;
          for( m = constituents.begin(); m != constituents.end(); ++m )
          {
         reco::CandidatePtr constit = *m;
         if(constit->pt() == 0)
         {
           edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
           continue;
         }
         fjInputs.push_back(fastjet::PseudoJet(constit->px(),constit->py(),constit->pz(),constit->energy()));
          }
        }
      }
      else
      {
        for(typename std::vector<IPTI>::const_iterator it = trackIPTagInfos->begin(); it != trackIPTagInfos->end(); ++it)
        {
          std::vector<edm::Ptr<reco::Candidate> > constituents = it->jet()->getJetConstituents();
          std::vector<edm::Ptr<reco::Candidate> >::const_iterator m;
          for( m = constituents.begin(); m != constituents.end(); ++m )
          {
         reco::CandidatePtr constit = *m;
         if(constit->pt() == 0)
         {
           edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
           continue;
         }
         fjInputs.push_back(fastjet::PseudoJet(constit->px(),constit->py(),constit->pz(),constit->energy()));
          }
        }
      }
      // insert "ghost" SVs in the vector of constituents
      for(typename edm::View<VTX>::const_iterator it = extSecVertex->begin(); it != extSecVertex->end(); ++it)
      {
        const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
        GlobalVector dir = flightDirection(pv, *it);
        dir = dir.unit();
        fastjet::PseudoJet p(dir.x(),dir.y(),dir.z(),dir.mag()); // using SV flight direction so treating SV as massless
        if( useSVMomentum )
          p = fastjet::PseudoJet(it->p4().px(),it->p4().py(),it->p4().pz(),it->p4().energy());
        p*=ghostRescaling; // rescale SV direction/momentum
        p.set_user_info(new VertexInfo( it - extSecVertex->begin() ));
        fjInputs.push_back(p);
      }
      
      // define jet clustering sequence
      fjClusterSeq = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs, *fjJetDefinition ) );
      // recluster jet constituents and inserted "ghosts"
      std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt( fjClusterSeq->inclusive_jets(jetPtMin) );
      
      if( useFatJets )
      {
        if( inclusiveJets.size() < fatJetsHandle->size() )
          edm::LogError("TooFewReclusteredJets") << "There are fewer reclustered (" << inclusiveJets.size() << ") than original fat jets (" << fatJetsHandle->size() << "). Please check that the jet algorithm and jet size match those used for the original jet collection.";

        // match reclustered and original fat jets
        std::vector<int> reclusteredIndices;
        matchReclusteredJets<edm::View<reco::Jet> >(fatJetsHandle,inclusiveJets,reclusteredIndices,"fat");

        // match groomed and original fat jets
        std::vector<int> groomedIndices;
        if( useGroomedFatJets )
          matchGroomedJets(fatJetsHandle,groomedFatJetsHandle,groomedIndices);

        // match subjets and original fat jets
        std::vector<std::vector<int> > subjetIndices;
        if( useGroomedFatJets )
          matchSubjets(groomedIndices,groomedFatJetsHandle,trackIPTagInfos,subjetIndices);
        else
          matchSubjets(fatJetsHandle,trackIPTagInfos,subjetIndices);
    
        // collect clustered SVs
        for(size_t i=0; i<fatJetsHandle->size(); ++i)
        {
          if( reclusteredIndices.at(i) < 0 ) continue; // continue if matching reclustered to original jets failed
    
          if( fatJetsHandle->at(i).pt() == 0 ) // continue if the original jet has Pt=0
          {
            edm::LogWarning("NullTransverseMomentum") << "The original fat jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
            continue;
          }
    
          if( subjetIndices.at(i).empty() ) continue; // continue if the original jet does not have subjets assigned
        
          // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original fat jets should in principle stay the same
          if( ( std::abs( inclusiveJets.at(reclusteredIndices.at(i)).pt() - fatJetsHandle->at(i).pt() ) / fatJetsHandle->at(i).pt() ) > relPtTolerance )
          {
         if( fatJetsHandle->at(i).pt() < 10. )  // special handling for low-Pt jets (Pt<10 GeV)
           edm::LogWarning("JetPtMismatchAtLowPt") << "The reclustered and original fat jet " << i << " have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << fatJetsHandle->at(i).pt() << " GeV, respectively).\n"
                               << "Please check that the jet algorithm and jet size match those used for the original fat jet collection and also make sure the original fat jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n"
                               << "Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n"
                               << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
         else
           edm::LogError("JetPtMismatch") << "The reclustered and original fat jet " << i << " have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << fatJetsHandle->at(i).pt() << " GeV, respectively).\n"
                          << "Please check that the jet algorithm and jet size match those used for the original fat jet collection and also make sure the original fat jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n"
                          << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
          }
    
          // get jet constituents
          std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(i)).constituents();

          std::vector<int> svIndices;
          // loop over jet constituents and try to find "ghosts"
          for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
          {
         if( !it->has_user_info() ) continue; // skip if not a "ghost"
         
         svIndices.push_back( it->user_info<VertexInfo>().vertexIndex() );
          }
    
          // loop over clustered SVs and assign them to different subjets based on smallest dR
          for(size_t sv=0; sv<svIndices.size(); ++sv)
          {
        const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
        const VTX &extSV = (*extSecVertex)[ svIndices.at(sv) ];
        GlobalVector dir = flightDirection(pv, extSV);
        dir = dir.unit();
        fastjet::PseudoJet p(dir.x(),dir.y(),dir.z(),dir.mag()); // using SV flight direction so treating SV as massless
        if( useSVMomentum )
          p = fastjet::PseudoJet(extSV.p4().px(),extSV.p4().py(),extSV.p4().pz(),extSV.p4().energy());
        
        std::vector<double> dR2toSubjets;
        
        for(size_t sj=0; sj<subjetIndices.at(i).size(); ++sj)
          dR2toSubjets.push_back( Geom::deltaR2( p.rapidity(), p.phi_std(), trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->rapidity(), trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->phi() ) );

        // find the closest subjet
        int closestSubjetIdx = std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );

        clusteredSVs.at(subjetIndices.at(i).at(closestSubjetIdx)).push_back( svIndices.at(sv) );
          }
        }
      }
      else
      {
        if( inclusiveJets.size() < trackIPTagInfos->size() )
          edm::LogError("TooFewReclusteredJets") << "There are fewer reclustered (" << inclusiveJets.size() << ") than original jets (" << trackIPTagInfos->size() << "). Please check that the jet algorithm and jet size match those used for the original jet collection.";
    
        // match reclustered and original jets
        std::vector<int> reclusteredIndices;
        matchReclusteredJets<std::vector<IPTI> >(trackIPTagInfos,inclusiveJets,reclusteredIndices);
    
        // collect clustered SVs
        for(size_t i=0; i<trackIPTagInfos->size(); ++i)
        {
          if( reclusteredIndices.at(i) < 0 ) continue; // continue if matching reclustered to original jets failed
    
          if( trackIPTagInfos->at(i).jet()->pt() == 0 ) // continue if the original jet has Pt=0
          {
            edm::LogWarning("NullTransverseMomentum") << "The original jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
            continue;
          }
    
          // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original jets should in principle stay the same
          if( ( std::abs( inclusiveJets.at(reclusteredIndices.at(i)).pt() - trackIPTagInfos->at(i).jet()->pt() ) / trackIPTagInfos->at(i).jet()->pt() ) > relPtTolerance )
          {
         if( trackIPTagInfos->at(i).jet()->pt() < 10. )  // special handling for low-Pt jets (Pt<10 GeV)
           edm::LogWarning("JetPtMismatchAtLowPt") << "The reclustered and original jet " << i << " have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << trackIPTagInfos->at(i).jet()->pt() << " GeV, respectively).\n"
                               << "Please check that the jet algorithm and jet size match those used for the original jet collection and also make sure the original jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n"
                               << "Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n"
                               << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
         else
           edm::LogError("JetPtMismatch") << "The reclustered and original jet " << i << " have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << trackIPTagInfos->at(i).jet()->pt() << " GeV, respectively).\n"
                          << "Please check that the jet algorithm and jet size match those used for the original jet collection and also make sure the original jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n"
                          << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
          }
    
          // get jet constituents
          std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(i)).constituents();

          // loop over jet constituents and try to find "ghosts"
          for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
          {
         if( !it->has_user_info() ) continue; // skip if not a "ghost"
         // push back clustered SV indices
         clusteredSVs.at(i).push_back( it->user_info<VertexInfo>().vertexIndex() );
          }
        }
      }
    }
    // case where fat jets are used to associate SVs to subjets but no SV clustering is performed
    else if( useExternalSV && !useSVClustering && !trackIPTagInfos->empty() && useFatJets )
    {
      // match groomed and original fat jets
      std::vector<int> groomedIndices;
      if( useGroomedFatJets )
        matchGroomedJets(fatJetsHandle,groomedFatJetsHandle,groomedIndices);

      // match subjets and original fat jets
      std::vector<std::vector<int> > subjetIndices;
      if( useGroomedFatJets )
        matchSubjets(groomedIndices,groomedFatJetsHandle,trackIPTagInfos,subjetIndices);
      else
        matchSubjets(fatJetsHandle,trackIPTagInfos,subjetIndices);

      // loop over fat jets
      for(size_t i=0; i<fatJetsHandle->size(); ++i)
      {
        if( fatJetsHandle->at(i).pt() == 0 ) // continue if the original jet has Pt=0
        {
          edm::LogWarning("NullTransverseMomentum") << "The original fat jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
          continue;
        }

        if( subjetIndices.at(i).empty() ) continue; // continue if the original jet does not have subjets assigned

        // loop over SVs, associate them to fat jets based on dR cone and
        // then assign them to the closets subjet in dR
        for(typename edm::View<VTX>::const_iterator it = extSecVertex->begin(); it != extSecVertex->end(); ++it)
        {
          size_t sv = ( it - extSecVertex->begin() );
    
          const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
          const VTX &extSV = (*extSecVertex)[sv];
          GlobalVector dir = flightDirection(pv, extSV);
          GlobalVector jetDir(fatJetsHandle->at(i).px(),
                  fatJetsHandle->at(i).py(),
                  fatJetsHandle->at(i).pz());
          // skip SVs outside the dR cone
              if( Geom::deltaR2( dir, jetDir ) > rParam*rParam ) // here using the jet clustering rParam as a dR cut
            continue;

          dir = dir.unit();
          fastjet::PseudoJet p(dir.x(),dir.y(),dir.z(),dir.mag()); // using SV flight direction so treating SV as massless
          if( useSVMomentum )
        p = fastjet::PseudoJet(extSV.p4().px(),extSV.p4().py(),extSV.p4().pz(),extSV.p4().energy());
    
          std::vector<double> dR2toSubjets;
    
          for(size_t sj=0; sj<subjetIndices.at(i).size(); ++sj)
        dR2toSubjets.push_back( Geom::deltaR2( p.rapidity(), p.phi_std(), trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->rapidity(), trackIPTagInfos->at(subjetIndices.at(i).at(sj)).jet()->phi() ) );

          // find the closest subjet
          int closestSubjetIdx = std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );

          clusteredSVs.at(subjetIndices.at(i).at(closestSubjetIdx)).push_back(sv);
        }
      }
    }
    // ------------------------------------ SV clustering END ----------------------------------------------

    std::auto_ptr<ConfigurableVertexReconstructor> vertexReco;
    std::auto_ptr<GhostTrackVertexFinder> vertexRecoGT;
    if (useGhostTrack)
        vertexRecoGT.reset(new GhostTrackVertexFinder(
            vtxRecoPSet.getParameter<double>("maxFitChi2"),
            vtxRecoPSet.getParameter<double>("mergeThreshold"),
            vtxRecoPSet.getParameter<double>("primcut"),
            vtxRecoPSet.getParameter<double>("seccut"),
            getGhostTrackFitType(vtxRecoPSet.getParameter<std::string>("fitType"))));
    else
        vertexReco.reset(
            new ConfigurableVertexReconstructor(vtxRecoPSet));

    TransientTrackMap primariesMap;

    // result secondary vertices

    auto tagInfos = std::make_unique<Product>();

    for(typename std::vector<IPTI>::const_iterator iterJets =
        trackIPTagInfos->begin(); iterJets != trackIPTagInfos->end();
        ++iterJets) {
        TrackDataVector trackData;
//            std::cout << "Jet " << iterJets-trackIPTagInfos->begin() << std::endl; 

        const Vertex &pv = *iterJets->primaryVertex();

        std::set<TransientTrack> primaries;
        if (constraint == CONSTRAINT_PV_PRIMARIES_IN_FIT) {
            for(Vertex::trackRef_iterator iter = pv.tracks_begin();
                iter != pv.tracks_end(); ++iter) {
                TransientTrackMap::iterator pos =
                    primariesMap.lower_bound(iter->get());

                if (pos != primariesMap.end() &&
                    pos->first == iter->get())
                    primaries.insert(pos->second);
                else {
                    TransientTrack track =
                        trackBuilder->build(
                            iter->castTo<TrackRef>());
                    primariesMap.insert(pos,
                        std::make_pair(iter->get(), track));
                    primaries.insert(track);
                }
            }
        }

        edm::RefToBase<Jet> jetRef = iterJets->jet();

        GlobalVector jetDir(jetRef->momentum().x(),
                            jetRef->momentum().y(),
                            jetRef->momentum().z());

        std::vector<std::size_t> indices =
                iterJets->sortedIndexes(sortCriterium);

        input_container trackRefs = iterJets->sortedTracks(indices);

        const std::vector<reco::btag::TrackIPData> &ipData =
                    iterJets->impactParameterData();

        // build transient tracks used for vertex reconstruction

        std::vector<TransientTrack> fitTracks;
        std::vector<GhostTrackState> gtStates;
        std::auto_ptr<GhostTrackPrediction> gtPred;
        if (useGhostTrack)
            gtPred.reset(new GhostTrackPrediction(
                        *iterJets->ghostTrack()));

        for(unsigned int i = 0; i < indices.size(); i++) {
            typedef typename TemplatedSecondaryVertexTagInfo<IPTI,VTX>::IndexedTrackData IndexedTrackData;

            const input_item &trackRef = trackRefs[i];

            trackData.push_back(IndexedTrackData());
            trackData.back().first = indices[i];

            // select tracks for SV finder

            if (!trackSelector(*reco::btag::toTrack(trackRef), ipData[indices[i]], *jetRef,
                               RecoVertex::convertPos(
                                    pv.position()))) {
                trackData.back().second.svStatus =
                      TemplatedSecondaryVertexTagInfo<IPTI,VTX>::TrackData::trackSelected;
                continue;
            }

            TransientTrackMap::const_iterator pos =
                    primariesMap.find(reco::btag::toTrack((trackRef)));
            TransientTrack fitTrack;
            if (pos != primariesMap.end()) {
                primaries.erase(pos->second);
                fitTrack = pos->second;
            } else
                fitTrack = trackBuilder->build(trackRef);
            fitTracks.push_back(fitTrack);

            trackData.back().second.svStatus =
                  TemplatedSecondaryVertexTagInfo<IPTI,VTX>::TrackData::trackUsedForVertexFit;

            if (useGhostTrack) {
                GhostTrackState gtState(fitTrack);
                GlobalPoint pos =
                    ipData[indices[i]].closestToGhostTrack;
                gtState.linearize(*gtPred, true,
                                  gtPred->lambda(pos));
                gtState.setWeight(ipData[indices[i]].ghostTrackWeight);
                gtStates.push_back(gtState);
            }
        }

        std::auto_ptr<GhostTrack> ghostTrack;
        if (useGhostTrack)
            ghostTrack.reset(new GhostTrack(
                GhostTrackPrediction(
                    RecoVertex::convertPos(pv.position()),
                    RecoVertex::convertError(pv.error()),
                    GlobalVector(
                        iterJets->ghostTrack()->px(),
                        iterJets->ghostTrack()->py(),
                        iterJets->ghostTrack()->pz()),
                    0.05),
                *gtPred, gtStates,
                iterJets->ghostTrack()->chi2(),
                iterJets->ghostTrack()->ndof()));

        // perform actual vertex finding


        std::vector<VTX>       extAssoCollection;    
        std::vector<TransientVertex> fittedSVs;
        std::vector<SecondaryVertex> SVs;
        if(!useExternalSV){ 
              switch(constraint)   {
            case CONSTRAINT_NONE:
            if (useGhostTrack)
                fittedSVs = vertexRecoGT->vertices(
                        pv, *ghostTrack);
            else
                fittedSVs = vertexReco->vertices(fitTracks);
            break;

            case CONSTRAINT_BEAMSPOT:
            if (useGhostTrack)
                fittedSVs = vertexRecoGT->vertices(
                        pv, *beamSpot, *ghostTrack);
            else
                fittedSVs = vertexReco->vertices(fitTracks,
                                                 *beamSpot);
            break;

            case CONSTRAINT_PV_BEAMSPOT_SIZE:
            case CONSTRAINT_PV_BS_Z_ERRORS_SCALED:
            case CONSTRAINT_PV_ERROR_SCALED: {
            BeamSpot::CovarianceMatrix cov;
            for(unsigned int i = 0; i < 7; i++) {
                unsigned int covSrc = bsCovSrc[i];
                for(unsigned int j = 0; j < 7; j++) {
                    double v=0.0;
                    if (!covSrc || bsCovSrc[j] != covSrc)
                        v = 0.0;
                    else if (covSrc == 1)
                        v = beamSpot->covariance(i, j);
                    else if (j<3 && i<3)
                        v = pv.covariance(i, j) *
                            constraintScaling;
                    cov(i, j) = v;
                }
            }

            BeamSpot bs(pv.position(), sigmaZ,
                        beamSpot.isValid() ? beamSpot->dxdz() : 0.,
                        beamSpot.isValid() ? beamSpot->dydz() : 0.,
                        beamWidth, cov, BeamSpot::Unknown);

            if (useGhostTrack)
                fittedSVs = vertexRecoGT->vertices(
                        pv, bs, *ghostTrack);
            else
                fittedSVs = vertexReco->vertices(fitTracks, bs);
            }   break;

            case CONSTRAINT_PV_PRIMARIES_IN_FIT: {
            std::vector<TransientTrack> primaries_(
                    primaries.begin(), primaries.end());
            if (useGhostTrack)
                fittedSVs = vertexRecoGT->vertices(
                        pv, *beamSpot, primaries_,
                        *ghostTrack);
            else
                fittedSVs = vertexReco->vertices(
                        primaries_, fitTracks,
                        *beamSpot);
            }   break;
          }
          // build combined SV information and filter
          SVBuilder svBuilder(pv, jetDir, withPVError, minTrackWeight);
          std::remove_copy_if(boost::make_transform_iterator(
                      fittedSVs.begin(), svBuilder),
                      boost::make_transform_iterator(
                      fittedSVs.end(), svBuilder),
                      std::back_inserter(SVs),
                      SVFilter(vertexFilter, pv, jetDir));

        }else{
          if( useSVClustering || useFatJets ) {
              size_t jetIdx = ( iterJets - trackIPTagInfos->begin() );
        
              for(size_t iExtSv = 0; iExtSv < clusteredSVs.at(jetIdx).size(); iExtSv++){
             const VTX & extVertex = (*extSecVertex)[ clusteredSVs.at(jetIdx).at(iExtSv) ];
             if( extVertex.p4().M() < 0.3 )
               continue;
             extAssoCollection.push_back( extVertex );
              }
          }
          else {
              for(size_t iExtSv = 0; iExtSv < extSecVertex->size(); iExtSv++){
             const VTX & extVertex = (*extSecVertex)[iExtSv];
             if( Geom::deltaR2( ( position(extVertex) - pv.position() ), jetDir ) > extSVDeltaRToJet*extSVDeltaRToJet || extVertex.p4().M() < 0.3 )
               continue;
             extAssoCollection.push_back( extVertex );
              }
          }
          // build combined SV information and filter
          SVBuilder svBuilder(pv, jetDir, withPVError, minTrackWeight);
          std::remove_copy_if(boost::make_transform_iterator( extAssoCollection.begin(), svBuilder),
                    boost::make_transform_iterator(extAssoCollection.end(), svBuilder),
                    std::back_inserter(SVs),
                    SVFilter(vertexFilter, pv, jetDir));
                }
        // clean up now unneeded collections
        gtPred.reset();
        ghostTrack.reset();
        gtStates.clear();
        fitTracks.clear();
        fittedSVs.clear();
        extAssoCollection.clear();

        // sort SVs by importance
        
        std::vector<unsigned int> vtxIndices = vertexSorting(SVs);

        std::vector<typename TemplatedSecondaryVertexTagInfo<IPTI,VTX>::VertexData> svData;

        svData.resize(vtxIndices.size());
        for(unsigned int idx = 0; idx < vtxIndices.size(); idx++) {
            const SecondaryVertex &sv = SVs[vtxIndices[idx]];

            svData[idx].vertex = sv;
            svData[idx].dist1d = sv.dist1d();
                        svData[idx].dist2d = sv.dist2d();
            svData[idx].dist3d = sv.dist3d();
            svData[idx].direction = flightDirection(pv,sv);
            // mark tracks successfully used in vertex fit
            markUsedTracks(trackData,trackRefs,sv,idx);
        }

        // fill result into tag infos

        tagInfos->push_back(
            TemplatedSecondaryVertexTagInfo<IPTI,VTX>(
                trackData, svData, SVs.size(),
                edm::Ref<std::vector<IPTI> >(trackIPTagInfos,
                    iterJets - trackIPTagInfos->begin())));
    }

    event.put(std::move(tagInfos));
}

//Need specialized template because reco::Vertex iterators are TrackBase and it is a mess to make general
template <>
void  TemplatedSecondaryVertexProducer<TrackIPTagInfo,reco::Vertex>::markUsedTracks(TrackDataVector & trackData, const input_container & trackRefs, const SecondaryVertex & sv,size_t idx)
{
    for(Vertex::trackRef_iterator iter = sv.tracks_begin(); iter != sv.tracks_end(); ++iter) {
        if (sv.trackWeight(*iter) < minTrackWeight)
            continue;

        typename input_container::const_iterator pos =
            std::find(trackRefs.begin(), trackRefs.end(),
                    iter->castTo<input_item>());

        if (pos == trackRefs.end() ) {
            if(!useExternalSV)
                throw cms::Exception("TrackNotFound")
                    << "Could not find track from secondary "
                    "vertex in original tracks."
                    << std::endl;
        } else {
            unsigned int index = pos - trackRefs.begin();
            trackData[index].second.svStatus =
                (btag::TrackData::Status)
                ((unsigned int)btag::TrackData::trackAssociatedToVertex + idx);
        }
    }
}
template <>
void  TemplatedSecondaryVertexProducer<CandIPTagInfo,reco::VertexCompositePtrCandidate>::markUsedTracks(TrackDataVector & trackData, const input_container & trackRefs, const SecondaryVertex & sv,size_t idx)
{
    for(typename input_container::const_iterator iter = sv.daughterPtrVector().begin(); iter != sv.daughterPtrVector().end(); ++iter)
    {
        typename input_container::const_iterator pos =
            std::find(trackRefs.begin(), trackRefs.end(), *iter);

        if (pos != trackRefs.end() )
        {
            unsigned int index = pos - trackRefs.begin();
            trackData[index].second.svStatus =
                (btag::TrackData::Status)
                ((unsigned int)btag::TrackData::trackAssociatedToVertex + idx);
        }
    }
}


template <>
typename TemplatedSecondaryVertexProducer<TrackIPTagInfo,reco::Vertex>::SecondaryVertex  
TemplatedSecondaryVertexProducer<TrackIPTagInfo,reco::Vertex>::SVBuilder::operator () (const TransientVertex &sv) const
{
    if(!sv.originalTracks().empty() && sv.originalTracks()[0].trackBaseRef().isNonnull())
        return SecondaryVertex(pv, sv, direction, withPVError);
    else
    {
        edm::LogError("UnexpectedInputs") << "Building from Candidates, should not happen!";
        return SecondaryVertex(pv, sv, direction, withPVError);
    }
}

template <>
typename TemplatedSecondaryVertexProducer<CandIPTagInfo,reco::VertexCompositePtrCandidate>::SecondaryVertex
TemplatedSecondaryVertexProducer<CandIPTagInfo,reco::VertexCompositePtrCandidate>::SVBuilder::operator () (const TransientVertex &sv) const
{
    if(!sv.originalTracks().empty() && sv.originalTracks()[0].trackBaseRef().isNonnull())
    {
        edm::LogError("UnexpectedInputs") << "Building from Tracks, should not happen!";
        VertexCompositePtrCandidate vtxCompPtrCand;
        
        vtxCompPtrCand.setCovariance(sv.vertexState().error().matrix());
        vtxCompPtrCand.setChi2AndNdof(sv.totalChiSquared(), sv.degreesOfFreedom());
        vtxCompPtrCand.setVertex(Candidate::Point(sv.position().x(),sv.position().y(),sv.position().z()));
        
        return SecondaryVertex(pv, vtxCompPtrCand, direction, withPVError);
    }
    else
    {
        VertexCompositePtrCandidate vtxCompPtrCand;
        
        vtxCompPtrCand.setCovariance(sv.vertexState().error().matrix());
        vtxCompPtrCand.setChi2AndNdof(sv.totalChiSquared(), sv.degreesOfFreedom());
        vtxCompPtrCand.setVertex(Candidate::Point(sv.position().x(),sv.position().y(),sv.position().z()));
        
        Candidate::LorentzVector p4;
        for(std::vector<reco::TransientTrack>::const_iterator tt = sv.originalTracks().begin(); tt != sv.originalTracks().end(); ++tt)
        {
            if (sv.trackWeight(*tt) < minTrackWeight)
                continue;
            
            const CandidatePtrTransientTrack* cptt = dynamic_cast<const CandidatePtrTransientTrack*>(tt->basicTransientTrack());
            if ( cptt==nullptr )
                edm::LogError("DynamicCastingFailed") << "Casting of TransientTrack to CandidatePtrTransientTrack failed!";
            else
            {
                p4 += cptt->candidate()->p4();
                vtxCompPtrCand.addDaughter(cptt->candidate());
            }
        }
        vtxCompPtrCand.setP4(p4);
        
        return SecondaryVertex(pv, vtxCompPtrCand, direction, withPVError);
    }
}

// ------------ method that matches reclustered and original jets based on minimum dR ------------
template<class IPTI,class VTX>
template<class CONTAINER>
void TemplatedSecondaryVertexProducer<IPTI,VTX>::matchReclusteredJets(const edm::Handle<CONTAINER>& jets,
                                                                      const std::vector<fastjet::PseudoJet>& reclusteredJets,
                                                                      std::vector<int>& matchedIndices,
                                                                      const std::string& jetType)
{
   std::string type = ( jetType!="" ? jetType + " " : jetType );

   std::vector<bool> matchedLocks(reclusteredJets.size(),false);

   for(size_t j=0; j<jets->size(); ++j)
   {
     double matchedDR2 = 1e9;
     int matchedIdx = -1;

     for(size_t rj=0; rj<reclusteredJets.size(); ++rj)
     {
       if( matchedLocks.at(rj) ) continue; // skip jets that have already been matched

       double tempDR2 = Geom::deltaR2( toJet(jets->at(j))->rapidity(), toJet(jets->at(j))->phi(), reclusteredJets.at(rj).rapidity(), reclusteredJets.at(rj).phi_std() );
       if( tempDR2 < matchedDR2 )
       {
         matchedDR2 = tempDR2;
         matchedIdx = rj;
       }
     }

     if( matchedIdx>=0 )
     {
       if ( matchedDR2 > rParam*rParam )
       {
         edm::LogError("JetMatchingFailed") << "Matched reclustered jet " << matchedIdx << " and original " << type << "jet " << j <<" are separated by dR=" << sqrt(matchedDR2) << " which is greater than the jet size R=" << rParam << ".\n"
                                            << "This is not expected so please check that the jet algorithm and jet size match those used for the original " << type << "jet collection.";
       }
       else
         matchedLocks.at(matchedIdx) = true;
     }
     else
       edm::LogError("JetMatchingFailed") << "Matching reclustered to original " << type << "jets failed. Please check that the jet algorithm and jet size match those used for the original " << type << "jet collection.";

     matchedIndices.push_back(matchedIdx);
   }
}

// ------------ method that matches groomed and original jets based on minimum dR ------------
template<class IPTI,class VTX>
void TemplatedSecondaryVertexProducer<IPTI,VTX>::matchGroomedJets(const edm::Handle<edm::View<reco::Jet> >& jets,
                                                                  const edm::Handle<edm::View<reco::Jet> >& groomedJets,
                                                                  std::vector<int>& matchedIndices)
{
   std::vector<bool> jetLocks(jets->size(),false);
   std::vector<int>  jetIndices;

   for(size_t gj=0; gj<groomedJets->size(); ++gj)
   {
     double matchedDR2 = 1e9;
     int matchedIdx = -1;

     if( groomedJets->at(gj).pt()>0. ) // skip pathological cases of groomed jets with Pt=0
     {
       for(size_t j=0; j<jets->size(); ++j)
       {
         if( jetLocks.at(j) ) continue; // skip jets that have already been matched

         double tempDR2 = Geom::deltaR2( jets->at(j).rapidity(), jets->at(j).phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi() );
         if( tempDR2 < matchedDR2 )
         {
           matchedDR2 = tempDR2;
           matchedIdx = j;
         }
       }
     }

     if( matchedIdx>=0 )
     {
       if ( matchedDR2 > rParam*rParam )
       {
         edm::LogWarning("MatchedJetsFarApart") << "Matched groomed jet " << gj << " and original jet " << matchedIdx <<" are separated by dR=" << sqrt(matchedDR2) << " which is greater than the jet size R=" << rParam << ".\n"
                                                << "This is not expected so the matching of these two jets has been discarded. Please check that the two jet collections belong to each other.";
         matchedIdx = -1;
       }
       else
         jetLocks.at(matchedIdx) = true;
     }
     jetIndices.push_back(matchedIdx);
   }

   for(size_t j=0; j<jets->size(); ++j)
   {
     std::vector<int>::iterator matchedIndex = std::find( jetIndices.begin(), jetIndices.end(), j );

     matchedIndices.push_back( matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(),matchedIndex) : -1 );
   }
}

// ------------ method that matches subjets and original fat jets ------------
template<class IPTI,class VTX>
void TemplatedSecondaryVertexProducer<IPTI,VTX>::matchSubjets(const std::vector<int>& groomedIndices,
                                                              const edm::Handle<edm::View<reco::Jet> >& groomedJets,
                                                              const edm::Handle<std::vector<IPTI> >& subjets,
                                                              std::vector<std::vector<int> >& matchedIndices)
{
   for(size_t g=0; g<groomedIndices.size(); ++g)
   {
     std::vector<int> subjetIndices;

     if( groomedIndices.at(g)>=0 )
     {
       for(size_t s=0; s<groomedJets->at(groomedIndices.at(g)).numberOfDaughters(); ++s)
       {
         const edm::Ptr<reco::Candidate> & subjet = groomedJets->at(groomedIndices.at(g)).daughterPtr(s);

         for(size_t sj=0; sj<subjets->size(); ++sj)
         {
           const edm::RefToBase<reco::Jet> &subjetRef = subjets->at(sj).jet();
           if( subjet == edm::Ptr<reco::Candidate>( subjetRef.id(), subjetRef.get(), subjetRef.key() ) )
           {
             subjetIndices.push_back(sj);
             break;
           }
         }
       }

       if( subjetIndices.empty() )
         edm::LogError("SubjetMatchingFailed") << "Matching subjets to original fat jets failed. Please check that the groomed fat jet and subjet collections belong to each other.";

       matchedIndices.push_back(subjetIndices);
     }
     else
       matchedIndices.push_back(subjetIndices);
   }
}

// ------------ method that matches subjets and original fat jets ------------
template<class IPTI,class VTX>
void TemplatedSecondaryVertexProducer<IPTI,VTX>::matchSubjets(const edm::Handle<edm::View<reco::Jet> >& fatJets,
                                                              const edm::Handle<std::vector<IPTI> >& subjets,
                                                              std::vector<std::vector<int> >& matchedIndices)
{
   for(size_t fj=0; fj<fatJets->size(); ++fj)
   {
     std::vector<int> subjetIndices;
     size_t nSubjetCollections = 0;
     size_t nSubjets = 0;

     const pat::Jet * fatJet = dynamic_cast<const pat::Jet *>( fatJets->ptrAt(fj).get() );

     if( !fatJet )
     {
       if( fj==0 ) edm::LogError("WrongJetType") << "Wrong jet type for input fat jets. Please check that the input fat jets are of the pat::Jet type.";

       matchedIndices.push_back(subjetIndices);
       continue;
     }
     else
     {
       nSubjetCollections = fatJet->subjetCollectionNames().size();

       if( nSubjetCollections>0 )
       {
         for(size_t coll=0; coll<nSubjetCollections; ++coll)
         {
           const pat::JetPtrCollection & fatJetSubjets = fatJet->subjets(coll);

           for(size_t fjsj=0; fjsj<fatJetSubjets.size(); ++fjsj)
           {
             ++nSubjets;

             for(size_t sj=0; sj<subjets->size(); ++sj)
             {
               const pat::Jet * subJet = dynamic_cast<const pat::Jet *>( subjets->at(sj).jet().get() );

               if( !subJet )
               {
                 if( fj==0 && coll==0 && fjsj==0 && sj==0 ) edm::LogError("WrongJetType") << "Wrong jet type for input subjets. Please check that the input subjets are of the pat::Jet type.";

                 break;
               }
               else
               {
                 if( subJet->originalObjectRef() == fatJetSubjets.at(fjsj)->originalObjectRef() )
                 {
                   subjetIndices.push_back(sj);
                   break;
                 }
               }
             }
           }
         }

         if( subjetIndices.empty() && nSubjets > 0)
           edm::LogError("SubjetMatchingFailed") << "Matching subjets to fat jets failed. Please check that the fat jet and subjet collections belong to each other.";

         matchedIndices.push_back(subjetIndices);
       }
       else
         matchedIndices.push_back(subjetIndices);
     }
   }
}

// ------------ method fills 'descriptions' with the allowed parameters for the module ------------
template<class IPTI,class VTX>
void TemplatedSecondaryVertexProducer<IPTI,VTX>::fillDescriptions(edm::ConfigurationDescriptions & descriptions) {

  edm::ParameterSetDescription desc;
  desc.add<double>("extSVDeltaRToJet",0.3);
  desc.add<edm::InputTag>("beamSpotTag",edm::InputTag("offlineBeamSpot"));
  {
    edm::ParameterSetDescription vertexReco;
    vertexReco.add<double>("primcut",1.8);
    vertexReco.add<double>("seccut",6.0);
    vertexReco.add<std::string>("finder","avr");
    vertexReco.addOptionalNode( edm::ParameterDescription<double>("minweight",0.5, true) and
                                edm::ParameterDescription<double>("weightthreshold",0.001, true) and
                                edm::ParameterDescription<bool>("smoothing",false, true), true );
    vertexReco.addOptionalNode( edm::ParameterDescription<double>("maxFitChi2",10.0, true) and
                                edm::ParameterDescription<double>("mergeThreshold",3.0, true) and
                                edm::ParameterDescription<std::string>("fitType","RefitGhostTrackWithVertices", true), true );
    desc.add<edm::ParameterSetDescription>("vertexReco",vertexReco);
  }
  {
    edm::ParameterSetDescription vertexSelection;
    vertexSelection.add<std::string>("sortCriterium","dist3dError");
    desc.add<edm::ParameterSetDescription>("vertexSelection",vertexSelection);
  }
  desc.add<std::string>("constraint","BeamSpot");
  desc.add<edm::InputTag>("trackIPTagInfos",edm::InputTag("impactParameterTagInfos"));
  {
    edm::ParameterSetDescription vertexCuts;
    vertexCuts.add<double>("distSig3dMax",99999.9);
    vertexCuts.add<double>("fracPV",0.65);
    vertexCuts.add<double>("distVal2dMax",2.5);
    vertexCuts.add<bool>("useTrackWeights",true);
    vertexCuts.add<double>("maxDeltaRToJetAxis",0.4);
    {
      edm::ParameterSetDescription v0Filter;
      v0Filter.add<double>("k0sMassWindow",0.05);
      vertexCuts.add<edm::ParameterSetDescription>("v0Filter",v0Filter);
    }
    vertexCuts.add<double>("distSig2dMin",3.0);
    vertexCuts.add<unsigned int>("multiplicityMin",2);
    vertexCuts.add<double>("distVal2dMin",0.01);
    vertexCuts.add<double>("distSig2dMax",99999.9);
    vertexCuts.add<double>("distVal3dMax",99999.9);
    vertexCuts.add<double>("minimumTrackWeight",0.5);
    vertexCuts.add<double>("distVal3dMin",-99999.9);
    vertexCuts.add<double>("massMax",6.5);
    vertexCuts.add<double>("distSig3dMin",-99999.9);
    desc.add<edm::ParameterSetDescription>("vertexCuts",vertexCuts);
  }
  desc.add<bool>("useExternalSV",false);
  desc.add<double>("minimumTrackWeight",0.5);
  desc.add<bool>("usePVError",true);
  {
    edm::ParameterSetDescription trackSelection;
    trackSelection.add<double>("b_pT",0.3684);
    trackSelection.add<double>("max_pT",500);
    trackSelection.add<bool>("useVariableJTA",false);
    trackSelection.add<double>("maxDecayLen",99999.9);
    trackSelection.add<double>("sip3dValMin",-99999.9);
    trackSelection.add<double>("max_pT_dRcut",0.1);
    trackSelection.add<double>("a_pT",0.005263);
    trackSelection.add<unsigned int>("totalHitsMin",8);
    trackSelection.add<double>("jetDeltaRMax",0.3);
    trackSelection.add<double>("a_dR",-0.001053);
    trackSelection.add<double>("maxDistToAxis",0.2);
    trackSelection.add<double>("ptMin",1.0);
    trackSelection.add<std::string>("qualityClass","any");
    trackSelection.add<unsigned int>("pixelHitsMin",2);
    trackSelection.add<double>("sip2dValMax",99999.9);
    trackSelection.add<double>("max_pT_trackPTcut",3);
    trackSelection.add<double>("sip2dValMin",-99999.9);
    trackSelection.add<double>("normChi2Max",99999.9);
    trackSelection.add<double>("sip3dValMax",99999.9);
    trackSelection.add<double>("sip3dSigMin",-99999.9);
    trackSelection.add<double>("min_pT",120);
    trackSelection.add<double>("min_pT_dRcut",0.5);
    trackSelection.add<double>("sip2dSigMax",99999.9);
    trackSelection.add<double>("sip3dSigMax",99999.9);
    trackSelection.add<double>("sip2dSigMin",-99999.9);
    trackSelection.add<double>("b_dR",0.6263);
    desc.add<edm::ParameterSetDescription>("trackSelection",trackSelection);
  }
  desc.add<std::string>("trackSort","sip3dSig");
  desc.add<edm::InputTag>("extSVCollection",edm::InputTag("secondaryVertices"));
  desc.addOptionalNode( edm::ParameterDescription<bool>("useSVClustering",false, true) and
                        edm::ParameterDescription<std::string>("jetAlgorithm", true) and
                        edm::ParameterDescription<double>("rParam", true), true );
  desc.addOptional<bool>("useSVMomentum",false);
  desc.addOptional<double>("ghostRescaling",1e-18);
  desc.addOptional<double>("relPtTolerance",1e-03);
  desc.addOptional<edm::InputTag>("fatJets");
  desc.addOptional<edm::InputTag>("groomedFatJets");
  descriptions.addDefault(desc);
}

//define this as a plug-in
typedef TemplatedSecondaryVertexProducer<TrackIPTagInfo,reco::Vertex> SecondaryVertexProducer;
typedef TemplatedSecondaryVertexProducer<CandIPTagInfo,reco::VertexCompositePtrCandidate> CandSecondaryVertexProducer;

DEFINE_FWK_MODULE(SecondaryVertexProducer);
DEFINE_FWK_MODULE(CandSecondaryVertexProducer);