CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CombinedSVComputer.h
Go to the documentation of this file.
1 #ifndef RecoBTag_SecondaryVertex_CombinedSVComputer_h
2 #define RecoBTag_SecondaryVertex_CombinedSVComputer_h
3 
4 #include <iostream>
5 #include <cstddef>
6 #include <string>
7 #include <cmath>
8 #include <vector>
9 
10 #include <Math/VectorUtil.h>
11 
15 
31 
36 
37 
38 #define range_for(i, x) \
39  for(int i = (x).begin; i != (x).end; i += (x).increment)
40 
41 
43  public:
44  explicit CombinedSVComputer(const edm::ParameterSet &params);
45 
47  operator () (const reco::TrackIPTagInfo &ipInfo,
48  const reco::SecondaryVertexTagInfo &svInfo) const;
50  operator () (const reco::CandIPTagInfo &ipInfo,
51  const reco::CandSecondaryVertexTagInfo &svInfo) const;
52 
53  struct IterationRange {
55  };
56  double flipValue(double value, bool vertex) const;
57  IterationRange flipIterate(int size, bool vertex) const;
58 
60 
62  threshTrack(const reco::CandIPTagInfo &trackIPTagInfo,
63  const reco::btag::SortCriteria sort,
64  const reco::Jet &jet,
65  const GlobalPoint &pv) const;
67  threshTrack(const reco::TrackIPTagInfo &trackIPTagInfo,
68  const reco::btag::SortCriteria sort,
69  const reco::Jet &jet,
70  const GlobalPoint &pv) const;
71  template <class SVTI,class IPTI>
73  const IPTI & ipInfo,const SVTI & svInfo,
74  double & vtx_track_ptSum, double & vtx_track_ESum) const;
75 
76  private:
77  bool trackFlip;
78  bool vertexFlip;
79  double charmCut;
84  unsigned int pseudoMultiplicityMin;
85  unsigned int trackMultiplicityMin;
91  std::vector<reco::btau::TaggingVariableName> taggingVariables;
92 };
93 
94 template <class SVTI,class IPTI>
96  const IPTI & ipInfo, const SVTI & svInfo,
97  double & vtx_track_ptSum, double & vtx_track_ESum) const
98 {
99  using namespace ROOT::Math;
100  using namespace reco;
101 
102  typedef typename IPTI::input_container Container;
103  typedef typename Container::value_type TrackRef;
104 
105  edm::RefToBase<Jet> jet = ipInfo.jet();
106  math::XYZVector jetDir = jet->momentum().Unit();
107  bool havePv = ipInfo.primaryVertex().isNonnull();
108  GlobalPoint pv;
109  if (havePv)
110  pv = GlobalPoint(ipInfo.primaryVertex()->x(),
111  ipInfo.primaryVertex()->y(),
112  ipInfo.primaryVertex()->z());
113 
115 
116 
117  vars.insert(btau::jetPt, jet->pt(), true);
118  vars.insert(btau::jetEta, jet->eta(), true);
119  vars.insert(btau::jetAbsEta, fabs(jet->eta()), true);
120 
121  if (ipInfo.selectedTracks().size() < trackMultiplicityMin)
122  return;
123 
124  vars.insert(btau::jetNTracks, ipInfo.selectedTracks().size(), true);
125 
126  TrackKinematics allKinematics;
127  TrackKinematics trackJetKinematics;
128 
129  double jet_track_ESum= 0.;
130 
131  int vtx = -1;
132 
133  IterationRange range = flipIterate(svInfo.nVertices(), true);
134  range_for(i , range) {
135  if (vtx < 0) vtx = i;
136  }
137 
138  if (vtx >= 0) {
139  vtxType = btag::Vertices::RecoVertex;
140 
141  vars.insert(btau::flightDistance2dVal,flipValue(svInfo.flightDistance(vtx, true).value(),true),true);
142  vars.insert(btau::flightDistance2dSig,flipValue(svInfo.flightDistance(vtx, true).significance(),true),true);
143  vars.insert(btau::flightDistance3dVal,flipValue(svInfo.flightDistance(vtx, false).value(),true),true);
144  vars.insert(btau::flightDistance3dSig,flipValue(svInfo.flightDistance(vtx, false).significance(),true),true);
145  vars.insert(btau::vertexJetDeltaR,Geom::deltaR(svInfo.flightDirection(vtx), jetDir),true);
146  vars.insert(btau::jetNSecondaryVertices, svInfo.nVertices(), true);
147  }
148 
149  std::vector<std::size_t> indices = ipInfo.sortedIndexes(sortCriterium);
150  const std::vector<reco::btag::TrackIPData> &ipData = ipInfo.impactParameterData();
151 
152  const Container &tracks = ipInfo.selectedTracks();
153  std::vector<TrackRef> pseudoVertexTracks;
154 
155  std::vector<TrackRef> trackPairV0Test(2);
156  range = flipIterate(indices.size(), false);
157  range_for(i, range) {
158  std::size_t idx = indices[i];
159  const reco::btag::TrackIPData &data = ipData[idx];
160  const TrackRef &track = tracks[idx];
161 
162  jet_track_ESum += std::sqrt(track->momentum().Mag2() + ROOT::Math::Square(ParticleMasses::piPlus));
163 
164  // add track to kinematics for all tracks in jet
165  //allKinematics.add(track); // would make more sense for some variables, e.g. vertexEnergyRatio nicely between 0 and 1, but not necessarily the best option for the discriminating power...
166 
167  // filter track -> this track selection can be tighter than the vertex track selection (used to fill the track related variables...)
168  if (!trackSelector(track, data, *jet, pv))
169  continue;
170 
171  // add track to kinematics for all tracks in jet
172  allKinematics.add(track);
173 
174  // if no vertex was reconstructed, attempt pseudo vertex
175  if (vtxType == btag::Vertices::NoVertex && trackPseudoSelector(track, data, *jet, pv)) {
176  pseudoVertexTracks.push_back(track);
177  vertexKinematics.add(track);
178  }
179 
180  // check against all other tracks for V0 track pairs
181  trackPairV0Test[0] = track;
182  bool ok = true;
183  range_for(j, range) {
184  if (i == j)
185  continue;
186 
187  std::size_t pairIdx = indices[j];
188  const reco::btag::TrackIPData &pairTrackData = ipData[pairIdx];
189  const TrackRef &pairTrack = tracks[pairIdx];
190 
191  if (!trackSelector(pairTrack, pairTrackData, *jet, pv))
192  continue;
193 
194  trackPairV0Test[1] = pairTrack;
195  if (!trackPairV0Filter(trackPairV0Test)) {
196  ok = false;
197  break;
198  }
199  }
200  if (!ok)
201  continue;
202 
203  trackJetKinematics.add(track);
204 
205  // add track variables
206  math::XYZVector trackMom = track->momentum();
207  double trackMag = std::sqrt(trackMom.Mag2());
208 
209  vars.insert(btau::trackSip3dVal, flipValue(data.ip3d.value(), false), true);
210  vars.insert(btau::trackSip3dSig, flipValue(data.ip3d.significance(), false), true);
211  vars.insert(btau::trackSip2dVal, flipValue(data.ip2d.value(), false), true);
212  vars.insert(btau::trackSip2dSig, flipValue(data.ip2d.significance(), false), true);
214 // vars.insert(btau::trackJetDistSig, data.distanceToJetAxis.significance(), true);
215 // vars.insert(btau::trackFirstTrackDist, data.distanceToFirstTrack, true);
216 // vars.insert(btau::trackGhostTrackVal, data.distanceToGhostTrack.value(), true);
217 // vars.insert(btau::trackGhostTrackSig, data.distanceToGhostTrack.significance(), true);
218  vars.insert(btau::trackDecayLenVal, havePv ? (data.closestToJetAxis - pv).mag() : -1.0, true);
219 
220  vars.insert(btau::trackMomentum, trackMag, true);
221  vars.insert(btau::trackEta, trackMom.Eta(), true);
222 
223  vars.insert(btau::trackPtRel, VectorUtil::Perp(trackMom, jetDir), true);
224  vars.insert(btau::trackPPar, jetDir.Dot(trackMom), true);
225  vars.insert(btau::trackDeltaR, VectorUtil::DeltaR(trackMom, jetDir), true);
226  vars.insert(btau::trackPtRatio, VectorUtil::Perp(trackMom, jetDir) / trackMag, true);
227  vars.insert(btau::trackPParRatio, jetDir.Dot(trackMom) / trackMag, true);
228  }
229 
230  if (vtxType == btag::Vertices::NoVertex && vertexKinematics.numberOfTracks() >= pseudoMultiplicityMin && pseudoVertexV0Filter(pseudoVertexTracks))
231  {
233  for(typename std::vector<TrackRef>::const_iterator trkIt = pseudoVertexTracks.begin(); trkIt != pseudoVertexTracks.end(); ++trkIt)
234  {
235  vars.insert(btau::trackEtaRel, reco::btau::etaRel(jetDir,(*trkIt)->momentum()), true);
236  vtx_track_ptSum += std::sqrt((*trkIt)->momentum().Perp2());
237  vtx_track_ESum += std::sqrt((*trkIt)->momentum().Mag2() + ROOT::Math::Square(ParticleMasses::piPlus));
238  }
239  }
240 
241  vars.insert(btau::vertexCategory, vtxType, true);
242 
243  vars.insert(btau::trackJetPt, trackJetKinematics.vectorSum().Pt(), true);
244  vars.insert(btau::trackSumJetDeltaR,VectorUtil::DeltaR(allKinematics.vectorSum(), jetDir), true);
245  vars.insert(btau::trackSumJetEtRatio,allKinematics.vectorSum().Et() / ipInfo.jet()->et(), true);
246 
247  vars.insert(btau::trackSip3dSigAboveCharm, flipValue(threshTrack(ipInfo, reco::btag::IP3DSig, *jet, pv).ip3d.significance(),false),true);
248  vars.insert(btau::trackSip3dValAboveCharm, flipValue(threshTrack(ipInfo, reco::btag::IP3DSig, *jet, pv).ip3d.value(),false),true);
249  vars.insert(btau::trackSip2dSigAboveCharm, flipValue(threshTrack(ipInfo, reco::btag::IP2DSig, *jet, pv).ip2d.significance(),false),true);
250  vars.insert(btau::trackSip2dValAboveCharm, flipValue(threshTrack(ipInfo, reco::btag::IP2DSig, *jet, pv).ip2d.value(),false),true);
251 
252  if (vtxType != btag::Vertices::NoVertex) {
254  ? allKinematics.weightedVectorSum()
255  : allKinematics.vectorSum();
257  ? vertexKinematics.weightedVectorSum()
258  : vertexKinematics.vectorSum();
259 
260  if (vtxType != btag::Vertices::RecoVertex) {
261  vars.insert(btau::vertexNTracks,vertexKinematics.numberOfTracks(), true);
262  vars.insert(btau::vertexJetDeltaR,VectorUtil::DeltaR(vertexSum, jetDir), true);
263  }
264 
265  double vertexMass = vertexSum.M();
266  if (vtxType == btag::Vertices::RecoVertex &&
268  GlobalVector dir = svInfo.flightDirection(vtx);
269  double vertexPt2 = math::XYZVector(dir.x(), dir.y(), dir.z()).Cross(vertexSum).Mag2() / dir.mag2();
270  vertexMass = std::sqrt(vertexMass * vertexMass + vertexPt2) + std::sqrt(vertexPt2);
271  }
272  vars.insert(btau::vertexMass, vertexMass, true);
273 
274  double varPi = (vertexMass/5.2794) * (vtx_track_ESum /jet_track_ESum); // 5.2794 should be the average B meson mass of PDG! CHECK!!!
275  vars.insert(btau::massVertexEnergyFraction, varPi / (varPi + 0.04), true);
276  double varB = (std::sqrt(5.2794) * vtx_track_ptSum) / ( vertexMass * std::sqrt(jet->pt()));
277  vars.insert(btau::vertexBoostOverSqrtJetPt,varB*varB/(varB*varB + 10.), true);
278 
279  if (allKinematics.numberOfTracks()) {
280  vars.insert(btau::vertexEnergyRatio, vertexSum.E() / allSum.E(), true);
281  }
282  else {
283  vars.insert(btau::vertexEnergyRatio, 1, true);
284  }
285  }
286 
287  reco::PFJet const * pfJet = dynamic_cast<reco::PFJet const *>( &* jet ) ;
288  pat::Jet const * patJet = dynamic_cast<pat::Jet const *>( &* jet ) ;
289  if ( pfJet != 0 )
290  {
291  vars.insert(btau::chargedHadronEnergyFraction,pfJet->chargedHadronEnergyFraction(), true);
292  vars.insert(btau::neutralHadronEnergyFraction,pfJet->neutralHadronEnergyFraction(), true);
293  vars.insert(btau::photonEnergyFraction,pfJet->photonEnergyFraction(), true);
294  vars.insert(btau::electronEnergyFraction,pfJet->electronEnergyFraction(), true);
295  vars.insert(btau::muonEnergyFraction,pfJet->muonEnergyFraction(), true);
296  vars.insert(btau::chargedHadronMultiplicity,pfJet->chargedHadronMultiplicity(), true);
297  vars.insert(btau::neutralHadronMultiplicity,pfJet->neutralHadronMultiplicity(), true);
298  vars.insert(btau::photonMultiplicity,pfJet->photonMultiplicity(), true);
299  vars.insert(btau::electronMultiplicity,pfJet->electronMultiplicity(), true);
300  vars.insert(btau::muonMultiplicity,pfJet->muonMultiplicity(), true);
301  vars.insert(btau::hadronMultiplicity,pfJet->chargedHadronMultiplicity()+pfJet->neutralHadronMultiplicity(), true);
302  vars.insert(btau::hadronPhotonMultiplicity,pfJet->chargedHadronMultiplicity()+pfJet->neutralHadronMultiplicity()+pfJet->photonMultiplicity(), true);
303  vars.insert(btau::totalMultiplicity,pfJet->chargedHadronMultiplicity()+pfJet->neutralHadronMultiplicity()+pfJet->photonMultiplicity()+pfJet->electronMultiplicity()+pfJet->muonMultiplicity(), true);
304 
305  }
306  else if( patJet != 0 && patJet->isPFJet() )
307  {
308  vars.insert(btau::chargedHadronEnergyFraction,patJet->chargedHadronEnergyFraction(), true);
309  vars.insert(btau::neutralHadronEnergyFraction,patJet->neutralHadronEnergyFraction(), true);
310  vars.insert(btau::photonEnergyFraction,patJet->photonEnergyFraction(), true);
311  vars.insert(btau::electronEnergyFraction,patJet->electronEnergyFraction(), true);
312  vars.insert(btau::muonEnergyFraction,patJet->muonEnergyFraction(), true);
313  vars.insert(btau::chargedHadronMultiplicity,patJet->chargedHadronMultiplicity(), true);
314  vars.insert(btau::neutralHadronMultiplicity,patJet->neutralHadronMultiplicity(), true);
315  vars.insert(btau::photonMultiplicity,patJet->photonMultiplicity(), true);
316  vars.insert(btau::electronMultiplicity,patJet->electronMultiplicity(), true);
317  vars.insert(btau::muonMultiplicity,patJet->muonMultiplicity(), true);
318  vars.insert(btau::hadronMultiplicity,patJet->chargedHadronMultiplicity()+patJet->neutralHadronMultiplicity(), true);
319  vars.insert(btau::hadronPhotonMultiplicity,patJet->chargedHadronMultiplicity()+patJet->neutralHadronMultiplicity()+patJet->photonMultiplicity(), true);
320  vars.insert(btau::totalMultiplicity,patJet->chargedHadronMultiplicity()+patJet->neutralHadronMultiplicity()+patJet->photonMultiplicity()+patJet->electronMultiplicity()+patJet->muonMultiplicity(), true);
321 
322  }
323  else
324  {
327  vars.insert(btau::photonEnergyFraction,0., true);
328  vars.insert(btau::electronEnergyFraction,0., true);
329  vars.insert(btau::muonEnergyFraction,0., true);
332  vars.insert(btau::photonMultiplicity,0, true);
333  vars.insert(btau::electronMultiplicity,0, true);
334  vars.insert(btau::muonMultiplicity,0, true);
335  vars.insert(btau::hadronMultiplicity,0, true);
337  vars.insert(btau::totalMultiplicity,0, true);
338  }
339 }
340 
341 
342 #endif // RecoBTag_SecondaryVertex_CombinedSVComputer_h
reco::TrackSelector trackSelector
const double piPlus
Definition: ParticleMasses.h:9
int i
Definition: DBlmapReader.cc:9
T mag2() const
Definition: PV3DBase.h:66
const math::XYZTLorentzVector & weightedVectorSum() const
unsigned int pseudoMultiplicityMin
void add(const reco::Track &track, double weight=1.0)
Base class for all types of Jets.
Definition: Jet.h:20
double etaRel(const math::XYZVector &dir, const math::XYZVector &track)
Measurement1D ip2d
Definition: IPTagInfo.h:31
reco::TrackSelector trackNoDeltaRSelector
reco::V0Filter pseudoVertexV0Filter
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:63
virtual reco::TaggingVariableList operator()(const reco::TrackIPTagInfo &ipInfo, const reco::SecondaryVertexTagInfo &svInfo) const
CombinedSVComputer(const edm::ParameterSet &params)
Jets made from PFObjects.
Definition: PFJet.h:21
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
reco::btag::SortCriteria sortCriterium
reco::TrackSelector trackPseudoSelector
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
int j
Definition: DBlmapReader.cc:9
Measurement1D ip3d
Definition: IPTagInfo.h:32
IterationRange flipIterate(int size, bool vertex) const
reco::V0Filter trackPairV0Filter
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
std::vector< reco::btau::TaggingVariableName > taggingVariables
double significance() const
Definition: Measurement1D.h:32
unsigned int numberOfTracks() const
#define range_for(i, x)
GlobalPoint closestToJetAxis
Definition: IPTagInfo.h:29
tuple tracks
Definition: testEve_cfg.py:39
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
double flipValue(double value, bool vertex) const
edm::Ref< TrackCollection > TrackRef
persistent reference to a Track
Definition: TrackFwd.h:20
Analysis-level calorimeter jet class.
Definition: Jet.h:78
edm::ParameterSet dropDeltaR(const edm::ParameterSet &pset) const
double value() const
Definition: Measurement1D.h:28
const reco::btag::TrackIPData & threshTrack(const reco::CandIPTagInfo &trackIPTagInfo, const reco::btag::SortCriteria sort, const reco::Jet &jet, const GlobalPoint &pv) const
Measurement1D distanceToJetAxis
Definition: IPTagInfo.h:33
void fillCommonVariables(reco::TaggingVariableList &vars, reco::TrackKinematics &vertexKinematics, const IPTI &ipInfo, const SVTI &svInfo, double &vtx_track_ptSum, double &vtx_track_ESum) const
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
const math::XYZTLorentzVector & vectorSum() const
dbl *** dir
Definition: mlp_gen.cc:35
T x() const
Definition: PV3DBase.h:62
edm::AssociationVector< reco::JetRefBaseProd, Values > Container
tuple size
Write out results.
void insert(const TaggingVariable &variable, bool delayed=false)
unsigned int trackMultiplicityMin