#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	hasAbilityToProduceInLumis () const final

bool	hasAbilityToProduceInRuns () 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 135 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 137 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 196 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 222 of file IPProducer.h.

223 {

224 }

Member Function Documentation

template<class Container , class Base , class Helper >

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

private

Definition at line 432 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 461 of file IPProducer.h.

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

                                                                                                                                            {
 
   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 482 of file IPProducer.h.

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

                                                                                                                                                             {
 
   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 231 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(-1. -1.);
            trackIP.ghostTrackWeight = 0.;
          }
        } else {
          trackIP.distanceToGhostTrack = Measurement1D(-1. -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));
 }