Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    TrackIPProducer
00004 // Class:      TrackIPProducer
00005 // 
00013 //
00014 // Original Author:  Andrea Rizzi
00015 //         Created:  Thu Apr  6 09:56:23 CEST 2006
00016 // $Id: TrackIPProducer.cc,v 1.24 2009/11/24 10:30:50 saout Exp $
00017 //
00018 //
00019 
00020 // system include files
00021 #include <cmath>
00022 #include <memory>
00023 #include <iostream>
00024 #include <algorithm>
00025 
00026 // user include files
00027 #include "FWCore/Framework/interface/Frameworkfwd.h"
00028 #include "FWCore/Framework/interface/Event.h"
00029 #include "FWCore/Framework/interface/MakerMacros.h"
00030 #include "FWCore/Framework/interface/EventSetup.h"
00031 #include "FWCore/Framework/interface/ESHandle.h"
00032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00033 
00034 #include "DataFormats/TrackReco/interface/Track.h"
00035 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
00036 #include "DataFormats/BTauReco/interface/JetTag.h"
00037 #include "DataFormats/BTauReco/interface/TrackIPTagInfo.h"
00038 
00039 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
00040 #include "TrackingTools/IPTools/interface/IPTools.h"
00041 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
00042 
00043 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
00044 #include "RecoVertex/VertexPrimitives/interface/ConvertError.h"
00045 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
00046 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
00047 #include "RecoVertex/GhostTrackFitter/interface/GhostTrack.h"
00048 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackState.h"
00049 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackPrediction.h"
00050 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackFitter.h"
00051 
00052 #include "RecoBTag/TrackProbability/interface/HistogramProbabilityEstimator.h"
00053 #include "RecoBTag/ImpactParameter/plugins/TrackIPProducer.h"
00054 
00055 
00056 using namespace std;
00057 using namespace reco;
00058 using namespace edm;
00059 using boost::bind;
00060 
00061 //
00062 // constructors and destructor
00063 //
00064 TrackIPProducer::TrackIPProducer(const edm::ParameterSet& iConfig) : 
00065   m_config(iConfig)
00066 {
00067   m_calibrationCacheId3D = 0;
00068   m_calibrationCacheId2D = 0;
00069 
00070   m_associator              = m_config.getParameter<InputTag>("jetTracks");
00071   m_primaryVertexProducer   = m_config.getParameter<InputTag>("primaryVertex");
00072   m_computeProbabilities    = m_config.getParameter<bool>("computeProbabilities");
00073   m_computeGhostTrack       = m_config.getParameter<bool>("computeGhostTrack");
00074   m_ghostTrackPriorDeltaR   = m_config.getParameter<double>("ghostTrackPriorDeltaR");
00075   m_cutPixelHits            = m_config.getParameter<int>("minimumNumberOfPixelHits");
00076   m_cutTotalHits            = m_config.getParameter<int>("minimumNumberOfHits");
00077   m_cutMaxTIP               = m_config.getParameter<double>("maximumTransverseImpactParameter");
00078   m_cutMinPt                = m_config.getParameter<double>("minimumTransverseMomentum");
00079   m_cutMaxChiSquared        = m_config.getParameter<double>("maximumChiSquared");
00080   m_cutMaxLIP               = m_config.getParameter<double>("maximumLongitudinalImpactParameter");
00081   m_directionWithTracks     = m_config.getParameter<bool>("jetDirectionUsingTracks");
00082   m_directionWithGhostTrack = m_config.getParameter<bool>("jetDirectionUsingGhostTrack");
00083   m_useTrackQuality         = m_config.getParameter<bool>("useTrackQuality");
00084 
00085   if (m_computeGhostTrack)
00086     produces<reco::TrackCollection>("ghostTracks");
00087   produces<reco::TrackIPTagInfoCollection>();
00088 }
00089 
00090 TrackIPProducer::~TrackIPProducer()
00091 {
00092 }
00093 
00094 //
00095 // member functions
00096 //
00097 // ------------ method called to produce the data  ------------
00098 void
00099 TrackIPProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00100 {
00101    // Update probability estimator if event setup is changed
00102    if (m_computeProbabilities)
00103      checkEventSetup(iSetup);
00104  
00105    //input objects 
00106    Handle<reco::JetTracksAssociationCollection> jetTracksAssociation;
00107    iEvent.getByLabel(m_associator, jetTracksAssociation);
00108    
00109    Handle<reco::VertexCollection> primaryVertex;
00110    iEvent.getByLabel(m_primaryVertexProducer, primaryVertex);
00111 
00112    edm::ESHandle<TransientTrackBuilder> builder;
00113    iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
00114    // m_algo.setTransientTrackBuilder(builder.product());
00115 
00116    // output collections 
00117    auto_ptr<reco::TrackIPTagInfoCollection> result(new reco::TrackIPTagInfoCollection);
00118 
00119    auto_ptr<reco::TrackCollection> ghostTracks;
00120    TrackRefProd ghostTrackRefProd;
00121    if (m_computeGhostTrack) {
00122      ghostTracks.reset(new reco::TrackCollection);
00123      ghostTrackRefProd = iEvent.getRefBeforePut<TrackCollection>("ghostTracks");
00124    }
00125 
00126    // use first pv of the collection
00127    Vertex dummy;
00128    const Vertex *pv = &dummy;
00129    edm::Ref<VertexCollection> pvRef;
00130    if (primaryVertex->size() != 0) {
00131      pv = &*primaryVertex->begin();
00132      // we always use the first vertex (at the moment)
00133      pvRef = edm::Ref<VertexCollection>(primaryVertex, 0);
00134    } else { // create a dummy PV
00135      Vertex::Error e;
00136      e(0, 0) = 0.0015 * 0.0015;
00137      e(1, 1) = 0.0015 * 0.0015;
00138      e(2, 2) = 15. * 15.;
00139      Vertex::Point p(0, 0, 0);
00140      dummy = Vertex(p, e, 0, 0, 0);
00141    }
00142 
00143    int i = 0;
00144    for(JetTracksAssociationCollection::const_iterator it =
00145                                         jetTracksAssociation->begin();
00146        it != jetTracksAssociation->end(); it++, i++) {
00147      TrackRefVector tracks = it->second;
00148      math::XYZVector jetMomentum = it->first->momentum();
00149 
00150      if (m_directionWithTracks) {
00151        jetMomentum *= 0.5;
00152        for(TrackRefVector::const_iterator itTrack = tracks.begin();
00153            itTrack != tracks.end(); ++itTrack)
00154          if ((**itTrack).numberOfValidHits() >= m_cutTotalHits)
00155            //minimal quality cuts
00156            jetMomentum += (*itTrack)->momentum();
00157      }
00158 
00159      TrackRefVector selectedTracks;
00160      vector<TransientTrack> transientTracks;
00161 
00162      for(TrackRefVector::const_iterator itTrack = tracks.begin();
00163          itTrack != tracks.end(); ++itTrack) {
00164        const Track & track = **itTrack;
00165        TransientTrack transientTrack = builder->build(*itTrack);
00166 /*     cout << " pt " <<  track.pt() <<
00167                " d0 " <<  fabs(track.d0()) <<
00168                " #hit " <<    track.hitPattern().numberOfValidHits()<<
00169                " ipZ " <<   fabs(track.dz()-pvZ)<<
00170                " chi2 " <<  track.normalizedChi2()<<
00171                " #pixel " <<    track.hitPattern().numberOfValidPixelHits()<< endl;
00172 */
00173        if (track.pt() > m_cutMinPt &&
00174            track.hitPattern().numberOfValidHits() >= m_cutTotalHits &&         // min num tracker hits
00175            track.hitPattern().numberOfValidPixelHits() >= m_cutPixelHits &&
00176            track.normalizedChi2() < m_cutMaxChiSquared &&
00177            std::abs(track.dxy(pv->position())) < m_cutMaxTIP &&
00178            std::abs(track.dz(pv->position())) < m_cutMaxLIP) {
00179          selectedTracks.push_back(*itTrack);
00180          transientTracks.push_back(transientTrack);
00181        }
00182      }
00183 
00184      GlobalVector direction(jetMomentum.x(), jetMomentum.y(), jetMomentum.z());
00185 
00186      auto_ptr<GhostTrack> ghostTrack;
00187      TrackRef ghostTrackRef;
00188      if (m_computeGhostTrack) {
00189        GhostTrackFitter fitter;
00190        GlobalPoint origin = RecoVertex::convertPos(pv->position());
00191        GlobalError error = RecoVertex::convertError(pv->error());
00192        ghostTrack.reset(new GhostTrack(fitter.fit(origin, error, direction,
00193                                                   m_ghostTrackPriorDeltaR,
00194                                                   transientTracks)));
00195 
00196 /*
00197         if (std::sqrt(jetMomentum.Perp2()) > 30) {
00198                 double offset = ghostTrack->prediction().lambda(origin);
00199                 std::cout << "------------------ jet pt " << std::sqrt(jetMomentum.Perp2()) << std::endl;
00200                 const std::vector<GhostTrackState> *states = &ghostTrack->states();
00201                 for(std::vector<GhostTrackState>::const_iterator state = states->begin();
00202                     state != states->end(); ++state) {
00203                         double dist = state->lambda() - offset;
00204                         double err = state->lambdaError(ghostTrack->prediction(), error);
00205                         double ipSig = IPTools::signedImpactParameter3D(state->track(), direction, *pv).second.significance();
00206                         double axisDist = state->axisDistance(ghostTrack->prediction());
00207                         std::cout << state->track().impactPointState().freeState()->momentum().perp()
00208                                   << ": " << dist << "/" << err << " [" << (dist / err) << "], ipsig = " << ipSig << ", dist = " << axisDist << ", w = " << state->weight() << std::endl;
00209                 }
00210         }
00211 */
00212        ghostTrackRef = TrackRef(ghostTrackRefProd, ghostTracks->size());
00213        ghostTracks->push_back(*ghostTrack);
00214 
00215        if (m_directionWithGhostTrack) { 
00216          const GhostTrackPrediction &pred = ghostTrack->prediction();
00217          double lambda = pred.lambda(origin);
00218          dummy = Vertex(RecoVertex::convertPos(pred.position(lambda)),
00219                         RecoVertex::convertError(pred.positionError(lambda)),
00220                         0, 0, 0);
00221          pv = &dummy;
00222          direction = pred.direction();
00223        }
00224      }
00225 
00226      vector<float> prob2D, prob3D;
00227      vector<TrackIPTagInfo::TrackIPData> ipData;
00228 
00229      for(unsigned int ind = 0; ind < transientTracks.size(); ind++) {
00230        const Track & track = *selectedTracks[ind];
00231        const TransientTrack &transientTrack = transientTracks[ind];
00232 
00233        TrackIPTagInfo::TrackIPData trackIP;
00234        trackIP.ip3d = IPTools::signedImpactParameter3D(transientTrack, direction, *pv).second;
00235        trackIP.ip2d = IPTools::signedTransverseImpactParameter(transientTrack, direction, *pv).second;
00236 
00237        TrajectoryStateOnSurface closest =
00238                IPTools::closestApproachToJet(transientTrack.impactPointState(),
00239                                              *pv, direction,
00240                                              transientTrack.field());
00241        if (closest.isValid())
00242          trackIP.closestToJetAxis = closest.globalPosition();
00243 
00244        // TODO: cross check if it is the same using other methods
00245        trackIP.distanceToJetAxis = IPTools::jetTrackDistance(transientTrack, direction, *pv).second;
00246 
00247        if (ghostTrack.get()) {
00248          const std::vector<GhostTrackState> &states = ghostTrack->states();
00249          std::vector<GhostTrackState>::const_iterator pos =
00250                 std::find_if(states.begin(), states.end(),
00251                              bind(equal_to<TransientTrack>(),
00252                                   bind(&GhostTrackState::track, _1),
00253                                   transientTrack));
00254 
00255          if (pos != states.end() && pos->isValid()) {
00256            VertexDistance3D dist;
00257            const GhostTrackPrediction &pred = ghostTrack->prediction();
00258            GlobalPoint p1 = pos->tsos().globalPosition();
00259            GlobalError e1 = pos->tsos().cartesianError().position();
00260            GlobalPoint p2 = pred.position(pos->lambda());
00261            GlobalError e2 = pred.positionError(pos->lambda());
00262            trackIP.closestToGhostTrack = p1;
00263            trackIP.distanceToGhostTrack = dist.distance(VertexState(p1, e1),
00264                                                         VertexState(p2, e2));
00265            trackIP.ghostTrackWeight = pos->weight();
00266          } else {
00267            trackIP.distanceToGhostTrack = Measurement1D(-1. -1.);
00268            trackIP.ghostTrackWeight = 0.;
00269          }
00270        } else {
00271          trackIP.distanceToGhostTrack = Measurement1D(-1. -1.);
00272          trackIP.ghostTrackWeight = 1.;
00273        }
00274 
00275        ipData.push_back(trackIP);
00276 
00277        if (m_computeProbabilities) {
00278          //probability with 3D ip
00279          pair<bool,double> probability = m_probabilityEstimator->probability(m_useTrackQuality, 0,ipData.back().ip3d.significance(),track,*(it->first),*pv);
00280          prob3D.push_back(probability.first ? probability.second : -1.);
00281 
00282          //probability with 2D ip
00283          probability = m_probabilityEstimator->probability(m_useTrackQuality,1,ipData.back().ip2d.significance(),track,*(it->first),*pv);
00284          prob2D.push_back(probability.first ? probability.second : -1.);
00285        } 
00286      }
00287 
00288      Ref<JetTracksAssociationCollection> jtaRef(jetTracksAssociation, i);
00289      result->push_back(
00290              TrackIPTagInfo(ipData, prob2D, prob3D, selectedTracks,
00291                             jtaRef, pvRef, direction, ghostTrackRef));
00292    }
00293  
00294    if (m_computeGhostTrack)
00295      iEvent.put(ghostTracks, "ghostTracks");
00296    iEvent.put(result);
00297 }
00298 
00299 
00300 #include "CondFormats/BTauObjects/interface/TrackProbabilityCalibration.h"
00301 #include "CondFormats/DataRecord/interface/BTagTrackProbability2DRcd.h"
00302 #include "CondFormats/DataRecord/interface/BTagTrackProbability3DRcd.h"
00303 #include "FWCore/Framework/interface/EventSetupRecord.h"
00304 #include "FWCore/Framework/interface/EventSetupRecordImplementation.h"
00305 #include "FWCore/Framework/interface/EventSetupRecordKey.h"
00306 
00307 void TrackIPProducer::checkEventSetup(const EventSetup & iSetup)
00308  {
00309   using namespace edm;
00310   using namespace edm::eventsetup;
00311 
00312    const EventSetupRecord & re2D= iSetup.get<BTagTrackProbability2DRcd>();
00313    const EventSetupRecord & re3D= iSetup.get<BTagTrackProbability3DRcd>();
00314    unsigned long long cacheId2D= re2D.cacheIdentifier();
00315    unsigned long long cacheId3D= re3D.cacheIdentifier();
00316 
00317    if(cacheId2D!=m_calibrationCacheId2D || cacheId3D!=m_calibrationCacheId3D  )  //Calibration changed
00318    {
00319      //iSetup.get<BTagTrackProbabilityRcd>().get(calib);
00320      ESHandle<TrackProbabilityCalibration> calib2DHandle;
00321      iSetup.get<BTagTrackProbability2DRcd>().get(calib2DHandle);
00322      ESHandle<TrackProbabilityCalibration> calib3DHandle;
00323      iSetup.get<BTagTrackProbability3DRcd>().get(calib3DHandle);
00324 
00325      const TrackProbabilityCalibration *  ca2D= calib2DHandle.product();
00326      const TrackProbabilityCalibration *  ca3D= calib3DHandle.product();
00327 
00328      m_probabilityEstimator.reset(new HistogramProbabilityEstimator(ca3D,ca2D));
00329 
00330    }
00331    m_calibrationCacheId3D=cacheId3D;
00332    m_calibrationCacheId2D=cacheId2D;
00333 }
/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/RecoBTag/ImpactParameter/plugins/TrackIPProducer.cc
