#include <IPProducer.h>

Inheritance diagram for IPProducer< Container, Base, Helper >:

Public Types
typedef std::vector< reco::IPTagInfo< Container, Base > >	Product

Public Types inherited from edm::stream::EDProducer<>
typedef CacheContexts< T... >	CacheTypes

typedef CacheTypes::GlobalCache	GlobalCache

typedef AbilityChecker< T... >	HasAbility

typedef CacheTypes::LuminosityBlockCache	LuminosityBlockCache

typedef LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >	LuminosityBlockContext

typedef CacheTypes::LuminosityBlockSummaryCache	LuminosityBlockSummaryCache

typedef CacheTypes::RunCache	RunCache

typedef RunContextT< RunCache, GlobalCache >	RunContext

typedef CacheTypes::RunSummaryCache	RunSummaryCache

Public Member Functions
template<>
void	fillDescriptions (edm::ConfigurationDescriptions &descriptions)

template<>
void	fillDescriptions (edm::ConfigurationDescriptions &descriptions)

	IPProducer (const edm::ParameterSet &)

void	produce (edm::Event &, const edm::EventSetup &) override

	~IPProducer () override

Public Member Functions inherited from edm::stream::EDProducer<>
	EDProducer ()=default

bool	hasAbilityToProduceInBeginLumis () const final

bool	hasAbilityToProduceInBeginRuns () const final

bool	hasAbilityToProduceInEndLumis () const final

bool	hasAbilityToProduceInEndRuns () const final

Static Public Member Functions
static void	fillDescriptions (edm::ConfigurationDescriptions &descriptions)

Private Member Functions
void	checkEventSetup (const edm::EventSetup &iSetup)

Private Attributes
unsigned long long	m_calibrationCacheId2D

unsigned long long	m_calibrationCacheId3D

bool	m_computeGhostTrack

bool	m_computeProbabilities

const edm::ParameterSet &	m_config

double	m_cutMaxChiSquared

double	m_cutMaxLIP

double	m_cutMaxTIP

double	m_cutMinPt

int	m_cutPixelHits

int	m_cutTotalHits

bool	m_directionWithGhostTrack

bool	m_directionWithTracks

double	m_ghostTrackPriorDeltaR

Helper	m_helper

std::unique_ptr< HistogramProbabilityEstimator >	m_probabilityEstimator

bool	m_useDB

bool	m_useTrackQuality

edm::EDGetTokenT< reco::VertexCollection >	token_primaryVertex

Detailed Description

template<class Container, class Base, class Helper>
class IPProducer< Container, Base, Helper >

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>="">

Definition at line 129 of file IPProducer.h.

Member Typedef Documentation

template<class Container , class Base , class Helper >

typedef std::vector<reco::IPTagInfo<Container, Base> > IPProducer< Container, Base, Helper >::Product

Definition at line 131 of file IPProducer.h.

Constructor & Destructor Documentation

template<class Container , class Base , class Helper >

IPProducer< Container, Base, Helper >::IPProducer ( const edm::ParameterSet & iConfig )

explicit

Definition at line 187 of file IPProducer.h.

     : m_config(iConfig), m_helper(iConfig, consumesCollector()) {
   m_calibrationCacheId3D = 0;
   m_calibrationCacheId2D = 0;
 
   token_primaryVertex = consumes<reco::VertexCollection>(m_config.getParameter<edm::InputTag>("primaryVertex"));
 
   m_computeProbabilities = m_config.getParameter<bool>("computeProbabilities");
   m_computeGhostTrack = m_config.getParameter<bool>("computeGhostTrack");
   m_ghostTrackPriorDeltaR = m_config.getParameter<double>("ghostTrackPriorDeltaR");
   m_cutPixelHits = m_config.getParameter<int>("minimumNumberOfPixelHits");
   m_cutTotalHits = m_config.getParameter<int>("minimumNumberOfHits");
   m_cutMaxTIP = m_config.getParameter<double>("maximumTransverseImpactParameter");
   m_cutMinPt = m_config.getParameter<double>("minimumTransverseMomentum");
   m_cutMaxChiSquared = m_config.getParameter<double>("maximumChiSquared");
   m_cutMaxLIP = m_config.getParameter<double>("maximumLongitudinalImpactParameter");
   m_directionWithTracks = m_config.getParameter<bool>("jetDirectionUsingTracks");
   m_directionWithGhostTrack = m_config.getParameter<bool>("jetDirectionUsingGhostTrack");
   m_useTrackQuality = m_config.getParameter<bool>("useTrackQuality");
 
   if (m_computeGhostTrack)
     produces<reco::TrackCollection>("ghostTracks");
   produces<Product>();
 }

template<class Container , class Base , class Helper >

IPProducer< Container, Base, Helper >::~IPProducer ( )

override

Definition at line 213 of file IPProducer.h.

213 {}

Member Function Documentation

template<class Container , class Base , class Helper >

void IPProducer< Container, Base, Helper >::checkEventSetup ( const edm::EventSetup & iSetup )

private

Definition at line 415 of file IPProducer.h.

References edm::eventsetup::EventSetupRecord::cacheIdentifier(), edm::EventSetup::get(), IPProducer< Container, Base, Helper >::m_calibrationCacheId2D, IPProducer< Container, Base, Helper >::m_calibrationCacheId3D, IPProducer< Container, Base, Helper >::m_probabilityEstimator, and edm::ESHandle< T >::product().

Referenced by IPProducer< Container, Base, Helper >::produce().

                                                                                      {
   const edm::eventsetup::EventSetupRecord& re2D = iSetup.get<BTagTrackProbability2DRcd>();
   const edm::eventsetup::EventSetupRecord& re3D = iSetup.get<BTagTrackProbability3DRcd>();
   unsigned long long cacheId2D = re2D.cacheIdentifier();
   unsigned long long cacheId3D = re3D.cacheIdentifier();
 
   if (cacheId2D != m_calibrationCacheId2D || cacheId3D != m_calibrationCacheId3D)  //Calibration changed
   {
     //iSetup.get<BTagTrackProbabilityRcd>().get(calib);
     edm::ESHandle<TrackProbabilityCalibration> calib2DHandle;
     iSetup.get<BTagTrackProbability2DRcd>().get(calib2DHandle);
     edm::ESHandle<TrackProbabilityCalibration> calib3DHandle;
     iSetup.get<BTagTrackProbability3DRcd>().get(calib3DHandle);
 
     const TrackProbabilityCalibration* ca2D = calib2DHandle.product();
     const TrackProbabilityCalibration* ca3D = calib3DHandle.product();
 
     m_probabilityEstimator.reset(new HistogramProbabilityEstimator(ca3D, ca2D));
   }
   m_calibrationCacheId3D = cacheId3D;
   m_calibrationCacheId2D = cacheId2D;
 }

template<class Container , class Base , class Helper >

static void IPProducer< Container, Base, Helper >::fillDescriptions ( edm::ConfigurationDescriptions & descriptions )

static

template<>

void IPProducer< reco::TrackRefVector, reco::JTATagInfo, IPProducerHelpers::FromJTA >::fillDescriptions ( edm::ConfigurationDescriptions & descriptions )

Definition at line 441 of file IPProducer.h.

References edm::ParameterSetDescription::add(), edm::ConfigurationDescriptions::addDefault(), and HLT_2018_cff::InputTag.

                                                 {
   edm::ParameterSetDescription desc;
   desc.add<double>("maximumTransverseImpactParameter", 0.2);
   desc.add<int>("minimumNumberOfHits", 8);
   desc.add<double>("minimumTransverseMomentum", 1.0);
   desc.add<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
   desc.add<double>("maximumLongitudinalImpactParameter", 17.0);
   desc.add<bool>("computeGhostTrack", true);
   desc.add<double>("ghostTrackPriorDeltaR", 0.03);
   desc.add<edm::InputTag>("jetTracks", edm::InputTag("ak4JetTracksAssociatorAtVertexPF"));
   desc.add<bool>("jetDirectionUsingGhostTrack", false);
   desc.add<int>("minimumNumberOfPixelHits", 2);
   desc.add<bool>("jetDirectionUsingTracks", false);
   desc.add<bool>("computeProbabilities", true);
   desc.add<bool>("useTrackQuality", false);
   desc.add<double>("maximumChiSquared", 5.0);
   descriptions.addDefault(desc);
 }

template<>

void IPProducer< std::vector< reco::CandidatePtr >, reco::JetTagInfo, IPProducerHelpers::FromJetAndCands >::fillDescriptions ( edm::ConfigurationDescriptions & descriptions )

Definition at line 462 of file IPProducer.h.

References edm::ParameterSetDescription::add(), edm::ConfigurationDescriptions::addDefault(), edm::ParameterSetDescription::addOptional(), and HLT_2018_cff::InputTag.

                                                 {
   edm::ParameterSetDescription desc;
   desc.add<double>("maximumTransverseImpactParameter", 0.2);
   desc.add<int>("minimumNumberOfHits", 8);
   desc.add<double>("minimumTransverseMomentum", 1.0);
   desc.add<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
   desc.add<double>("maximumLongitudinalImpactParameter", 17.0);
   desc.add<bool>("computeGhostTrack", true);
   desc.add<double>("maxDeltaR", 0.4);
   desc.add<edm::InputTag>("candidates", edm::InputTag("particleFlow"));
   desc.add<bool>("jetDirectionUsingGhostTrack", false);
   desc.add<int>("minimumNumberOfPixelHits", 2);
   desc.add<bool>("jetDirectionUsingTracks", false);
   desc.add<bool>("computeProbabilities", true);
   desc.add<bool>("useTrackQuality", false);
   desc.add<edm::InputTag>("jets", edm::InputTag("ak4PFJetsCHS"));
   desc.add<double>("ghostTrackPriorDeltaR", 0.03);
   desc.add<double>("maximumChiSquared", 5.0);
   desc.addOptional<bool>("explicitJTA", false);
   descriptions.addDefault(desc);
 }

template<class Container , class Base , class Helper >

void IPProducer< Container, Base, Helper >::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

override

Definition at line 220 of file IPProducer.h.

                                                                                                {
   // Update probability estimator if event setup is changed
   if (m_computeProbabilities)
     checkEventSetup(iSetup);
 
   edm::Handle<reco::VertexCollection> primaryVertex;
   iEvent.getByToken(token_primaryVertex, primaryVertex);
 
   edm::ESHandle<TransientTrackBuilder> builder;
   iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
   // m_algo.setTransientTrackBuilder(builder.product());
 
   // output collections
   auto result = std::make_unique<Product>();
 
   std::unique_ptr<reco::TrackCollection> ghostTracks;
   reco::TrackRefProd ghostTrackRefProd;
   if (m_computeGhostTrack) {
     ghostTracks = std::make_unique<reco::TrackCollection>();
     ghostTrackRefProd = iEvent.getRefBeforePut<reco::TrackCollection>("ghostTracks");
   }
 
   // use first pv of the collection
   reco::Vertex dummy;
   const reco::Vertex* pv = &dummy;
   edm::Ref<reco::VertexCollection> pvRef;
   if (!primaryVertex->empty()) {
     pv = &*primaryVertex->begin();
     // we always use the first vertex (at the moment)
     pvRef = edm::Ref<reco::VertexCollection>(primaryVertex, 0);
   } else {  // create a dummy PV
     reco::Vertex::Error e;
     e(0, 0) = 0.0015 * 0.0015;
     e(1, 1) = 0.0015 * 0.0015;
     e(2, 2) = 15. * 15.;
     reco::Vertex::Point p(0, 0, 0);
     dummy = reco::Vertex(p, e, 0, 0, 0);
   }
 
   std::vector<Base> baseTagInfos = m_helper.makeBaseVector(iEvent);
   for (typename std::vector<Base>::const_iterator it = baseTagInfos.begin(); it != baseTagInfos.end(); it++) {
     Container tracks = m_helper.tracks(*it);
     math::XYZVector jetMomentum = it->jet()->momentum();
 
     if (m_directionWithTracks) {
       jetMomentum *= 0.5;
       for (typename Container::const_iterator itTrack = tracks.begin(); itTrack != tracks.end(); ++itTrack)
         if (reco::btag::toTrack(*itTrack)->numberOfValidHits() >= m_cutTotalHits)  //minimal quality cuts
           jetMomentum += (*itTrack)->momentum();
     }
 
     Container selectedTracks;
     std::vector<reco::TransientTrack> transientTracks;
 
     for (typename Container::const_iterator itTrack = tracks.begin(); itTrack != tracks.end(); ++itTrack) {
       reco::TransientTrack transientTrack = builder->build(*itTrack);
       const reco::Track& track = transientTrack.track();  //**itTrack;
                                                           /*    cout << " pt " <<  track.pt() <<
                " d0 " <<  fabs(track.d0()) <<
                " #hit " <<    track.hitPattern().numberOfValidHits()<<
                " ipZ " <<   fabs(track.dz()-pv->z())<<
                " chi2 " <<  track.normalizedChi2()<<
                " #pixel " <<    track.hitPattern().numberOfValidPixelHits()<< endl;
 */
       if (track.pt() > m_cutMinPt &&
           track.hitPattern().numberOfValidHits() >= m_cutTotalHits &&  // min num tracker hits
           track.hitPattern().numberOfValidPixelHits() >= m_cutPixelHits &&
           track.normalizedChi2() < m_cutMaxChiSquared && std::abs(track.dxy(pv->position())) < m_cutMaxTIP &&
           std::abs(track.dz(pv->position())) < m_cutMaxLIP) {
         //   std::cout << "selected" << std::endl;
         selectedTracks.push_back(*itTrack);
         transientTracks.push_back(transientTrack);
       }
     }
     //  std::cout <<"SIZE: " << transientTracks.size() << std::endl;
     GlobalVector direction(jetMomentum.x(), jetMomentum.y(), jetMomentum.z());
 
     std::unique_ptr<reco::GhostTrack> ghostTrack;
     reco::TrackRef ghostTrackRef;
     if (m_computeGhostTrack) {
       reco::GhostTrackFitter fitter;
       GlobalPoint origin = RecoVertex::convertPos(pv->position());
       GlobalError error = RecoVertex::convertError(pv->error());
       ghostTrack.reset(
           new reco::GhostTrack(fitter.fit(origin, error, direction, m_ghostTrackPriorDeltaR, transientTracks)));
 
       /*
     if (std::sqrt(jetMomentum.Perp2()) > 30) {
         double offset = ghostTrack->prediction().lambda(origin);
         std::cout << "------------------ jet pt " << std::sqrt(jetMomentum.Perp2()) << std::endl;
         const std::vector<GhostTrackState> *states = &ghostTrack->states();
         for(std::vector<GhostTrackState>::const_iterator state = states->begin();
             state != states->end(); ++state) {
             double dist = state->lambda() - offset;
             double err = state->lambdaError(ghostTrack->prediction(), error);
             double ipSig = IPTools::signedImpactParameter3D(state->track(), direction, *pv).second.significance();
             double axisDist = state->axisDistance(ghostTrack->prediction());
             std::cout << state->track().impactPointState().freeState()->momentum().perp()
                       << ": " << dist << "/" << err << " [" << (dist / err) << "], ipsig = " << ipSig << ", dist = " << axisDist << ", w = " << state->weight() << std::endl;
         }
     }
 */
       ghostTrackRef = reco::TrackRef(ghostTrackRefProd, ghostTracks->size());
       ghostTracks->push_back(*ghostTrack);
 
       if (m_directionWithGhostTrack) {
         const reco::GhostTrackPrediction& pred = ghostTrack->prediction();
         double lambda = pred.lambda(origin);
         dummy = reco::Vertex(RecoVertex::convertPos(pred.position(lambda)),
                              RecoVertex::convertError(pred.positionError(lambda)),
                              0,
                              0,
                              0);
         pv = &dummy;
         direction = pred.direction();
       }
     }
 
     std::vector<float> prob2D, prob3D;
     std::vector<reco::btag::TrackIPData> ipData;
 
     for (unsigned int ind = 0; ind < transientTracks.size(); ind++) {
       const reco::TransientTrack& transientTrack = transientTracks[ind];
       const reco::Track& track = transientTrack.track();
 
       reco::btag::TrackIPData trackIP;
       trackIP.ip3d = IPTools::signedImpactParameter3D(transientTrack, direction, *pv).second;
       trackIP.ip2d = IPTools::signedTransverseImpactParameter(transientTrack, direction, *pv).second;
 
       TrajectoryStateOnSurface closest =
           IPTools::closestApproachToJet(transientTrack.impactPointState(), *pv, direction, transientTrack.field());
       if (closest.isValid())
         trackIP.closestToJetAxis = closest.globalPosition();
 
       // TODO: cross check if it is the same using other methods
       trackIP.distanceToJetAxis = IPTools::jetTrackDistance(transientTrack, direction, *pv).second;
 
       if (ghostTrack.get()) {
         const std::vector<reco::GhostTrackState>& states = ghostTrack->states();
         std::vector<reco::GhostTrackState>::const_iterator pos = std::find_if(
             states.begin(),
             states.end(),
             bind(std::equal_to<reco::TransientTrack>(), bind(&reco::GhostTrackState::track, _1), transientTrack));
 
         if (pos != states.end() && pos->isValid()) {
           VertexDistance3D dist;
           const reco::GhostTrackPrediction& pred = ghostTrack->prediction();
           GlobalPoint p1 = pos->tsos().globalPosition();
           GlobalError e1 = pos->tsos().cartesianError().position();
           GlobalPoint p2 = pred.position(pos->lambda());
           GlobalError e2 = pred.positionError(pos->lambda());
           trackIP.closestToGhostTrack = p1;
           trackIP.distanceToGhostTrack = dist.distance(VertexState(p1, e1), VertexState(p2, e2));
           trackIP.ghostTrackWeight = pos->weight();
         } else {
           trackIP.distanceToGhostTrack = Measurement1D(0., -1.);
           trackIP.ghostTrackWeight = 0.;
         }
       } else {
         trackIP.distanceToGhostTrack = Measurement1D(0., -1.);
         trackIP.ghostTrackWeight = 1.;
       }
 
       ipData.push_back(trackIP);
 
       if (m_computeProbabilities) {
         //probability with 3D ip
         std::pair<bool, double> probability = m_probabilityEstimator->probability(
             m_useTrackQuality, 0, ipData.back().ip3d.significance(), track, *(it->jet()), *pv);
         prob3D.push_back(probability.first ? probability.second : -1.);
 
         //probability with 2D ip
         probability = m_probabilityEstimator->probability(
             m_useTrackQuality, 1, ipData.back().ip2d.significance(), track, *(it->jet()), *pv);
         prob2D.push_back(probability.first ? probability.second : -1.);
       }
     }
 
     result->push_back(
         typename Product::value_type(ipData, prob2D, prob3D, selectedTracks, *it, pvRef, direction, ghostTrackRef));
   }
 
   if (m_computeGhostTrack)
     iEvent.put(std::move(ghostTracks), "ghostTracks");
   iEvent.put(std::move(result));
 }