CMS 3D CMS Logo

TrackVertexArbitration.h
Go to the documentation of this file.
1 #ifndef TrackVertexArbitration_H
2 #define TrackVertexArbitration_H
3 #include <memory>
4 #include <set>
5 #include <unordered_map>
6 
12 
13 
19 
26 
27 
33 
42 
44 
46 //#include "DataFormats/Math/interface/deltaR.h"
47 
48 //#define VTXDEBUG
49 
50 namespace svhelper {
51  inline double cov33(const reco::Vertex & sv) { return sv.covariance(3,3); }
52  inline double cov33(const reco::VertexCompositePtrCandidate & sv) { return sv.vertexCovariance(3,3); }
53 }
54 
55 
56 
57 template <class VTX>
59  public:
61 
62 
63  std::vector<VTX> trackVertexArbitrator(
65  const reco::Vertex &pv,
66  std::vector<reco::TransientTrack> & selectedTracks,
67  std::vector<VTX> & secondaryVertices
68  );
69 
70 
71  private:
72  bool trackFilterArbitrator(const reco::TransientTrack &track) const;
73 
78  double dRCut;
79  double distCut;
80  double sigCut;
81  double dLenFraction;
82  double fitterSigmacut; // = 3.0;
83  double fitterTini; // = 256.;
84  double fitterRatio;// = 0.25;
86  double trackMinPt;
89 };
90 
92 template <class VTX>
94  primaryVertexCollection (params.getParameter<edm::InputTag>("primaryVertices")),
95  secondaryVertexCollection (params.getParameter<edm::InputTag>("secondaryVertices")),
96  trackCollection (params.getParameter<edm::InputTag>("tracks")),
97  beamSpotCollection (params.getParameter<edm::InputTag>("beamSpot")),
98  dRCut (params.getParameter<double>("dRCut")),
99  distCut (params.getParameter<double>("distCut")),
100  sigCut (params.getParameter<double>("sigCut")),
101  dLenFraction (params.getParameter<double>("dLenFraction")),
102  fitterSigmacut (params.getParameter<double>("fitterSigmacut")),
103  fitterTini (params.getParameter<double>("fitterTini")),
104  fitterRatio (params.getParameter<double>("fitterRatio")),
105  trackMinLayers (params.getParameter<int32_t>("trackMinLayers")),
106  trackMinPt (params.getParameter<double>("trackMinPt")),
107  trackMinPixels (params.getParameter<int32_t>("trackMinPixels")),
108  maxTimeSignificance (params.getParameter<double>("maxTimeSignificance"))
109 {
110  dRCut*=dRCut;
111 }
112 template <class VTX>
114 {
115  if(!track.isValid())
116  return false;
118  return false;
119  if (track.track().pt() < trackMinPt)
120  return false;
122  return false;
123 
124  return true;
125 }
126 
128 {
129  return sv.trackWeight(track.trackBaseRef());
130 }
132 {
134  if ( cptt==nullptr )
135  edm::LogError("DynamicCastingFailed") << "Casting of TransientTrack to CandidatePtrTransientTrack failed!";
136  else
137  {
138  const reco::CandidatePtr & c=cptt->candidate();
139  if(std::find(sv.daughterPtrVector().begin(),sv.daughterPtrVector().end(),c)!= sv.daughterPtrVector().end())
140  return 1.0;
141  else
142  return 0;
143  }
144 return 0;
145 }
146 
147 
148 
149 template <class VTX>
152  const reco::Vertex &pv,
153  std::vector<reco::TransientTrack> & selectedTracks,
154  std::vector<VTX> & secondaryVertices)
155 {
156  using namespace reco;
157 
158  //std::cout << "PV: " << pv.position() << std::endl;
159  VertexDistance3D dist;
160  AdaptiveVertexFitter theAdaptiveFitter(
166 
167 
168 
169  std::vector<VTX> recoVertices;
170  VertexDistance3D vdist;
171  GlobalPoint ppv(pv.position().x(),pv.position().y(),pv.position().z());
172 
173  std::unordered_map<unsigned int, Measurement1D> cachedIP;
174  for(typename std::vector<VTX>::const_iterator sv = secondaryVertices.begin();
175  sv != secondaryVertices.end(); ++sv) {
176 
177  const double svTime(sv->t()), svTimeCov(svhelper::cov33(*sv));
178  const bool svHasTime = svTimeCov > 0.;
179 
180  GlobalPoint ssv(sv->position().x(),sv->position().y(),sv->position().z());
181  GlobalVector flightDir = ssv-ppv;
182 // std::cout << "Vertex : " << sv-secondaryVertices->begin() << " " << sv->position() << std::endl;
183  Measurement1D dlen= vdist.distance(pv,VertexState(RecoVertex::convertPos(sv->position()),RecoVertex::convertError(sv->error())));
184  std::vector<reco::TransientTrack> selTracks;
185  for(unsigned int itrack = 0; itrack < selectedTracks.size(); itrack++){
186  TransientTrack & tt = (selectedTracks)[itrack];
187  if (!trackFilterArbitrator(tt)) continue;
188  tt.setBeamSpot(*beamSpot);
189  float w = trackWeight(*sv,tt);
190  Measurement1D ipv;
191  if( cachedIP.count(itrack) ) {
192  ipv=cachedIP[itrack];
193  } else {
194  std::pair<bool,Measurement1D> ipvp = IPTools::absoluteImpactParameter3D(tt,pv);
195  cachedIP[itrack]=ipvp.second;
196  ipv=ipvp.second;
197  }
198 
199  AnalyticalImpactPointExtrapolator extrapolator(tt.field());
200  TrajectoryStateOnSurface tsos = extrapolator.extrapolate(tt.impactPointState(), RecoVertex::convertPos(sv->position()));
201  if(! tsos.isValid()) continue;
202  GlobalPoint refPoint = tsos.globalPosition();
203  GlobalError refPointErr = tsos.cartesianError().position();
204  Measurement1D isv = dist.distance(VertexState(RecoVertex::convertPos(sv->position()),RecoVertex::convertError(sv->error())),VertexState(refPoint, refPointErr));
205 
206  float dR = Geom::deltaR2(flightDir,tt.track()); //.eta(), flightDir.phi(), tt.track().eta(), tt.track().phi());
207 
208  double timeSig = 0.;
209  if( svHasTime && edm::isFinite(tt.timeExt()) ) {
210  double tError = std::sqrt( std::pow( tt.dtErrorExt(), 2 ) + svTimeCov );
211  timeSig = std::abs(tt.timeExt() - svTime)/tError;
212  }
213 
214  if( w > 0 || ( isv.significance() < sigCut && isv.value() < distCut && isv.value() < dlen.value()*dLenFraction && timeSig < maxTimeSignificance) )
215  {
216 
217  if(( isv.value() < ipv.value() ) && isv.value() < distCut && isv.value() < dlen.value()*dLenFraction
218  && dR < dRCut && timeSig < maxTimeSignificance )
219  {
220 #ifdef VTXDEBUG
221  if(w > 0.5) std::cout << " = ";
222  else std::cout << " + ";
223 #endif
224  selTracks.push_back(tt);
225  } else
226  {
227 #ifdef VTXDEBUG
228  if(w > 0.5 && isv.value() > ipv.value() ) std::cout << " - ";
229  else std::cout << " ";
230 #endif
231  //add also the tracks used in previous fitting that are still closer to Sv than Pv
232  if(w > 0.5 && isv.value() <= ipv.value() && dR < dRCut && timeSig < maxTimeSignificance) {
233  selTracks.push_back(tt);
234 #ifdef VTXDEBUG
235  std::cout << " = ";
236 #endif
237  }
238  if(w > 0.5 && isv.value() <= ipv.value() && dR >= dRCut) {
239 #ifdef VTXDEBUG
240  std::cout << " - ";
241 #endif
242 
243  }
244 
245 
246  }
247 #ifdef VTXDEBUG
248 
249  std::cout << "t : " << itrack << " ref " << ref.key() << " w: " << w
250  << " svip: " << isv.significance() << " " << isv.value()
251  << " pvip: " << ipv.significance() << " " << ipv.value() << " res " << tt.track().residualX(0) << "," << tt.track().residualY(0)
252 // << " tpvip: " << itpv.second.significance() << " " << itpv.second.value() << " dr: " << dR << std::endl;
253 #endif
254 
255  }
256  else
257  {
258 #ifdef VTXDEBUG
259 
260  std::cout << " . t : " << itrack << " ref " << ref.key()<< " w: " << w
261  << " svip: " << isv.second.significance() << " " << isv.second.value()
262  << " pvip: " << ipv.significance() << " " << ipv.value() << " dr: " << dR << std::endl;
263 #endif
264 
265  }
266  }
267 
268  if(selTracks.size() >= 2)
269  {
270  TransientVertex singleFitVertex;
271  singleFitVertex = theAdaptiveFitter.vertex(selTracks,ssv);
272 
273  if(singleFitVertex.isValid()) {
274  svhelper::updateVertexTime(singleFitVertex);
275  recoVertices.push_back(VTX(singleFitVertex));
276 
277 
278 #ifdef VTXDEBUG
279  const VTX & extVertex = recoVertices.back();
280  GlobalVector vtxDir = GlobalVector(extVertex.p4().X(),extVertex.p4().Y(),extVertex.p4().Z());
281 
282 
283  std::cout << " pointing : " << Geom::deltaR(extVertex.position() - pv.position(), vtxDir) << std::endl;
284 #endif
285  }
286  }
287  }
288  return recoVertices;
289 }
290 
291 
292 #endif
reco::Vertex::Point convertPos(const GlobalPoint &p)
const double w
Definition: UKUtility.cc:23
void setBeamSpot(const reco::BeamSpot &beamSpot)
bool isValid() const
Make the ReferenceCountingProxy method to check validity public.
TrackBaseRef trackBaseRef() const
reco::Vertex::Error convertError(const GlobalError &ge)
Definition: ConvertError.h:8
std::pair< bool, Measurement1D > absoluteImpactParameter3D(const reco::TransientTrack &transientTrack, const reco::Vertex &vertex)
Definition: IPTools.cc:37
bool trackFilterArbitrator(const reco::TransientTrack &track) const
double timeExt() const
CachingVertex< 5 > vertex(const std::vector< reco::TransientTrack > &) const override
const MagneticField * field() const
double covariance(int i, int j) const
(i, j)-th element of error matrix, i, j = 0, ... 2
Definition: Vertex.h:130
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:20
int trackerLayersWithMeasurement() const
Definition: HitPattern.cc:518
const Point & position() const
position
Definition: Vertex.h:109
virtual const daughters & daughterPtrVector() const
references to daughtes
bool isFinite(T x)
Measurement1D distance(const GlobalPoint &vtx1Position, const GlobalError &vtx1PositionError, const GlobalPoint &vtx2Position, const GlobalError &vtx2PositionError) const override
T sqrt(T t)
Definition: SSEVec.h:18
double pt() const
track transverse momentum
Definition: TrackBase.h:621
edm::InputTag secondaryVertexCollection
double vertexCovariance(int i, int j) const override
(i, j)-th element of error matrix, i, j = 0, ... 3
def pv(vc)
Definition: MetAnalyzer.py:6
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
float trackWeight(const TREF &r) const
returns the weight with which a Track has contributed to the vertex-fit.
Definition: Vertex.h:81
CandidatePtr candidate() const override
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
const BasicTransientTrack * basicTransientTrack() const
double significance() const
Definition: Measurement1D.h:32
const Track & track() const
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:446
double dtErrorExt() const
double value() const
Definition: Measurement1D.h:28
std::vector< VTX > trackVertexArbitrator(edm::Handle< reco::BeamSpot > &beamSpot, const reco::Vertex &pv, std::vector< reco::TransientTrack > &selectedTracks, std::vector< VTX > &secondaryVertices)
T1 deltaR2(T1 eta1, T2 phi1, T3 eta2, T4 phi2)
Definition: deltaR.h:36
void updateVertexTime(TransientVertex &vtx)
Definition: SVTimeHelpers.h:11
fixed size matrix
HLT enums.
double cov33(const reco::Vertex &sv)
TrackVertexArbitration(const edm::ParameterSet &params)
int numberOfValidPixelHits() const
Definition: HitPattern.h:838
TrajectoryStateOnSurface impactPointState() const
bool isValid() const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
float trackWeight(const reco::Vertex &sv, const reco::TransientTrack &track)
Global3DVector GlobalVector
Definition: GlobalVector.h:10