ATLAS Offline Software
JetFitterVariablesFactory.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
17 #include "VxJetVertex/VxVertexOnJetAxis.h"
18 #include "VxVertex/VxTrackAtVertex.h"
19 #include "VxJetVertex/VxJetCandidate.h"
20 #include "VxJetVertex/SelectedTracksInJet.h"
21 
22 #include "TrkParameters/TrackParameters.h"
23 #include "VxSecVertex/VxJetFitterVertexInfo.h"
24 
25 #include "TrkTrackLink/ITrackLink.h"
26 
27 #include "CLHEP/Vector/LorentzVector.h"
28 
29 #include "JetTagTools/JetFitterVariablesFactory.h"
30 
31 #include "GeoPrimitives/GeoPrimitivesHelpers.h"
32 
33 #include <vector>
34 #include <cmath>
35 
36 
37 namespace Analysis {
38 
39  JetFitterVariablesFactory::JetFitterVariablesFactory(const std::string& name,
40  const std::string& n, const IInterface* p):
41  AthAlgTool(name, n,p),
42  m_addNegativeTracksToPrimaryVertex(false),
43  m_usePtCorrectedEnergy(false),
44  m_useSingleTracksAlsoForMass(false)
45  {
46  declareProperty("addNegativeTracksToPrimaryVertex",m_addNegativeTracksToPrimaryVertex);
47  declareProperty("usePtCorrectedEnergy",m_usePtCorrectedEnergy);
48  declareProperty("useSingleTracksAlsoForMass",m_useSingleTracksAlsoForMass);
49  declareInterface<IJetFitterVariablesFactory>(this);
50  }
51 
55 
56  JetFitterVariablesFactory::~JetFitterVariablesFactory() {}
57 
58 StatusCode JetFitterVariablesFactory::initialize() {
59  ATH_MSG_INFO(" Initialization of JetFitterVariablesFactory successful");
60  return StatusCode::SUCCESS;
61 }
62 
63 StatusCode JetFitterVariablesFactory::finalize() {
64  ATH_MSG_INFO(" Finalization of JetFitterVariablesFactory successful");
65  return StatusCode::SUCCESS;
66 }
67 
68 //this is a nasty hack but doesn't disrupt BTagging workflow with fillJetFitterVariables
69 StatusCode JetFitterVariablesFactory::computeJetFitterVariables(const xAOD::Jet &myJet,
70  const Trk::VxJetFitterVertexInfo* myJetFitterInfo,
71  const std::string& basename,
72  JetFitterVariables &vars) const {
73 
74  bool nastyVsRevertPosToNeg = (basename.find("Flip") != std::string::npos);
75  int nVTX = 0, ndof = 0, nTracksAtVtx = 0, nSingleTracks = 0;
76  float energyFraction = NAN, mass = NAN, mass_uncorr = NAN, significance3d = NAN, deltaphi = NAN, deltaeta = NAN, chi2 = 0., deltaRFlightDir = NAN;
77 
78  std::vector<Trk::VxJetCandidate*> myVertices;
79  Trk::VxJetCandidate* myVxJetCandidate = nullptr;
80  if (myJetFitterInfo) myVertices = myJetFitterInfo->verticesJF();
81  if(myVertices.size() == 0){
82  ATH_MSG_DEBUG("#BTAG# Trk::VxJetCandidate not found for jet fitter ");
83  vars.massUncorr = mass_uncorr;
84  vars.chi2 = chi2;
85  vars.ndof = ndof;
86  vars.dRFlightDir = deltaRFlightDir;
87  vars.nVTX = nVTX;
88  vars.nSingleTracks = nSingleTracks;
89  vars.nTracksAtVtx = nTracksAtVtx;
90  vars.mass = mass;
91  vars.energyFraction = energyFraction;
92  vars.significance3d = significance3d;
93  vars.deltaeta = deltaeta;
94  vars.deltaphi = deltaphi;
95  return StatusCode::SUCCESS;
96  }
97  if(myVertices.size() > 0) myVxJetCandidate=dynamic_cast<Trk::VxJetCandidate*>(myVertices[0]);
98  if (myVxJetCandidate==0) {
99  ATH_MSG_WARNING("#BTAG# No correct VxJetCandidate could be retrieved." );
100  vars.massUncorr = mass_uncorr;
101  vars.chi2 = chi2;
102  vars.ndof = ndof;
103  vars.dRFlightDir = deltaRFlightDir;
104  vars.nVTX = nVTX;
105  vars.nSingleTracks = nSingleTracks;
106  vars.nTracksAtVtx = nTracksAtVtx;
107  vars.mass = mass;
108  vars.energyFraction = energyFraction;
109  vars.significance3d = significance3d;
110  vars.deltaeta = deltaeta;
111  vars.deltaphi = deltaphi;
112  return StatusCode::SUCCESS;
113  }
114 
115  const Trk::VxJetCandidate& myJetCandidate =*myVxJetCandidate;
116 
117  const Trk::SelectedTracksInJet* mySelectedTracksInJet = myJetFitterInfo->getSelectedTracksInJet();
118 
119 
120  //put all needed information inside :-)
121 
122  //const double s_massks=ParticleConstants::KZeroMassInMeV;
123  const double s_pion=ParticleConstants::chargedPionMassInMeV;
124 
125  double energyFromPrimary=0.;
126  double energyFromSecondary=0.;
127 
128  // get fit quality variables for the PV of jetfitter
129  const Trk::VxVertexOnJetAxis* pvtxjet = myVxJetCandidate->getPrimaryVertex();
130  const Trk::FitQuality& fitquality = pvtxjet->fitQuality();
131  chi2 = fitquality.chiSquared();
132  ndof = fitquality.numberDoF();
133 
134  if (mySelectedTracksInJet!=0)
135  {
136  ATH_MSG_DEBUG(" Adding the tracks from primary vertex information ");
137  const std::vector<const Trk::ITrackLink*> & myPrimaryLinks=mySelectedTracksInJet->getPrimaryTrackLinks();
138 
139  std::vector<const Trk::ITrackLink*>::const_iterator myPrimaryLinksBegin=myPrimaryLinks.begin();
140  std::vector<const Trk::ITrackLink*>::const_iterator myPrimaryLinksEnd=myPrimaryLinks.end();
141 
142  for(std::vector<const Trk::ITrackLink*>::const_iterator myPrimaryLinksIter=myPrimaryLinksBegin;
143  myPrimaryLinksIter!=myPrimaryLinksEnd;
144  ++myPrimaryLinksIter)
145  {
146  const Trk::TrackParameters* myParameters=(*myPrimaryLinksIter)->parameters();
147  if (myParameters)
148  {
149  energyFromPrimary+=std::sqrt(s_pion*s_pion+myParameters->momentum().mag2());
150  }
151  else
152  {
153  ATH_MSG_WARNING(" no perigee in track for energy computation. Skipping primary track...");
154  }
155  }
156  } else
157  {
158  ATH_MSG_DEBUG(" No information about further primary tracks available. Normal in JetFitter vs. 1");
159  }
160 
161 
162  const Trk::RecVertexPositions & recVertexPositions=myJetCandidate.getRecVertexPositions();
163  const Amg::VectorX & vertexPosition=recVertexPositions.position();
164  const Amg::MatrixX & vertexCovMatrix = recVertexPositions.covariancePosition();
165 
166  Amg::Vector3D primaryPos(vertexPosition[Trk::jet_xv],
167  vertexPosition[Trk::jet_yv],
168  vertexPosition[Trk::jet_zv]);
169 
170  Amg::Vector3D flightAxis(1,1,1);//has to be different from 0
171  Amg::setPhi(flightAxis, vertexPosition[Trk::jet_phi]);
172  Amg::setTheta(flightAxis, vertexPosition[Trk::jet_theta]);
173 
174  xAOD::IParticle::FourMom_t JetVector = myJet.p4();
175 
176  //loop over primary vertex
177  const std::vector<Trk::VxTrackAtVertex*> & TracksAtPrimary=myJetCandidate.getPrimaryVertex()->getTracksAtVertex();
178  const std::vector<Trk::VxTrackAtVertex*>::const_iterator TracksAtPrimaryBegin=TracksAtPrimary.begin();
179  const std::vector<Trk::VxTrackAtVertex*>::const_iterator TracksAtPrimaryEnd=TracksAtPrimary.end();
180 
181  for (std::vector<Trk::VxTrackAtVertex*>::const_iterator TracksAtPrimaryIter=TracksAtPrimaryBegin;
182  TracksAtPrimaryIter!=TracksAtPrimaryEnd;
183  ++TracksAtPrimaryIter) {
184 
185  //FIXME: dynamic cast necessary? neutral perigee commented out, fix when vertex supports neutral
186  if (dynamic_cast<const Trk::Perigee*>((*TracksAtPrimaryIter)->perigeeAtVertex())!=0)
187  {
188 
189  energyFromPrimary+=
190  std::sqrt(s_pion*s_pion+
191  (*TracksAtPrimaryIter)->perigeeAtVertex()->momentum().mag2());
192  }
193  else
194  {
195  ATH_MSG_ERROR(" FIXME: VERTEX DOESN'T SUPPORT NEUTRAL PERIGEE, commented out in line 163");
196  ATH_MSG_ERROR(" Track is not a normal track neither a KS. This is an ERROR (ask developer to fix it). Skipping track... ");
197  }
198  }
199 
200 
201 
202 
203  Amg::Vector3D sumPAllVertices(0.,0.,0.);
204  CLHEP::HepLorentzVector massVector(0,0,0,0);
205  double sumPtAdd(0.);
206 
207  double dist(0.);
208  double inverrordist(0.);
209 
210  //now access the vertices on the jet axis info...
211  std::vector<Trk::VxVertexOnJetAxis*> vectorOfVertices=myJetCandidate.getVerticesOnJetAxis();
212 
213  //then you have to order them...
214  std::vector<Trk::VxVertexOnJetAxis*>::iterator vectorOfClustersOfTrackBegin=vectorOfVertices.begin();
215  std::vector<Trk::VxVertexOnJetAxis*>::iterator vectorOfClustersOfTrackEnd=vectorOfVertices.end();
216 
217  for (std::vector<Trk::VxVertexOnJetAxis*>::const_iterator vectorOfClustersOfTrackIter=vectorOfClustersOfTrackBegin;
218  vectorOfClustersOfTrackIter!=vectorOfClustersOfTrackEnd;
219  ++vectorOfClustersOfTrackIter) {
220 
221  const std::vector<Trk::VxTrackAtVertex*> & tracksOfVertex=(*vectorOfClustersOfTrackIter)->getTracksAtVertex();
222 
223  int vertexSize=tracksOfVertex.size();
224  int ntrack=(*vectorOfClustersOfTrackIter)->getNumVertex()+5;//gets the right component (should add EDM method which does
225  if (!nastyVsRevertPosToNeg)
226  {
227  if (vertexPosition[ntrack]>0) {
228  if (vertexSize>1) {
229  nVTX+=1;
230  nTracksAtVtx+=vertexSize;
231  } else {
232  nSingleTracks+=1;
233  }
234  }
235  }
236  else
237  {
238  if (vertexPosition[ntrack]<=0) {
239  if (vertexSize>1) {
240  nVTX+=1;
241  nTracksAtVtx+=vertexSize;
242  } else {
243  nSingleTracks+=1;
244  }
245  }
246  }
247  }
248 
249  for (std::vector<Trk::VxVertexOnJetAxis*>::const_iterator vectorOfClustersOfTrackIter=vectorOfClustersOfTrackBegin;
250  vectorOfClustersOfTrackIter!=vectorOfClustersOfTrackEnd;
251  ++vectorOfClustersOfTrackIter) {
252 
253  const std::vector<Trk::VxTrackAtVertex*> & tracksOfVertex=(*vectorOfClustersOfTrackIter)->getTracksAtVertex();
254  std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackBegin=tracksOfVertex.begin();
255  std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackEnd=tracksOfVertex.end();
256 
257  int vertexSize=tracksOfVertex.size();
258 
259  int ntrack=(*vectorOfClustersOfTrackIter)->getNumVertex()+5;//gets the right component (should add EDM method which does
260  // this nasty little addition of 5...)
261 
262 
263  if ((vertexPosition[ntrack]<0 && (!nastyVsRevertPosToNeg))||(vertexPosition[ntrack]>=0 && nastyVsRevertPosToNeg)) {
264  if (m_addNegativeTracksToPrimaryVertex)
265  {
266  for (std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackIter=clustersOfTrackBegin;
267  clustersOfTrackIter!=clustersOfTrackEnd;++clustersOfTrackIter) {
268 
269  energyFromPrimary+=
270  std::hypot(s_pion, (*clustersOfTrackIter)->perigeeAtVertex()->momentum().mag());
271  }
272  }
273  } else {
274 
275 if ( (nVTX>0 && vertexSize>1) || nVTX==0 ) {
276  dist+=std::abs(vertexPosition[ntrack])/vertexCovMatrix(ntrack,ntrack);
277  if (vertexCovMatrix(ntrack,ntrack)>0)
278  {
279  inverrordist+=1./vertexCovMatrix(ntrack,ntrack);
280  }
281  else
282  {
283  ATH_MSG_WARNING("The diagonal element of the vertex cov matrix ("<<ntrack<<","<<ntrack<<") is "<<vertexCovMatrix(ntrack,ntrack)<<". It should be positive... Ignoring vertex when computing L/sigma(L)");
284  }
285 }
286 
287 Amg::Vector3D sumP(0.,0.,0.);
288 CLHEP::HepLorentzVector massThisCluster(0.,0.,0.,0.);
289 
290 //in case it's a real seconday vertex track...
291 
292 for (std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackIter=clustersOfTrackBegin;
293  clustersOfTrackIter!=clustersOfTrackEnd;
294  ++clustersOfTrackIter) {
295 
296 // const Trk::MeasuredPerigee* aMeasPer=static_cast<const Trk::MeasuredPerigee*>((*clustersOfTrackIter)->perigeeAtVertex());
297  const Trk::TrackParameters* aMeasPer=(*clustersOfTrackIter)->perigeeAtVertex();
298 // CLHEP::HepVector perigeeParms = aMeasPer->parameters();
299 
300  Amg::Vector3D mytrack(aMeasPer->momentum());
301  sumP+=mytrack;
302  if (dynamic_cast<const Trk::Perigee*>((*clustersOfTrackIter)->perigeeAtVertex())!=0)
303  {
304  massThisCluster+=CLHEP::HepLorentzVector(mytrack.x(), mytrack.y(), mytrack.z(), std::hypot(s_pion, mytrack.mag()));
305  }
306  else
307  {
308  ATH_MSG_ERROR("Neutral parameter has been taken out until Vertex has been rewritten to support neutral perigee again. ");
309  ATH_MSG_ERROR(" Track is not a normal track neither a KS. This is an ERROR (ask developer to fix it). Skipping track... ");
310  }
311 }
312 
313 
314 sumPAllVertices+=sumP;
315 double ptadd=sumP.perp(flightAxis.unit());
316 double masswithneutrals=std::sqrt(massThisCluster.mag2()+ptadd*ptadd)+ptadd;
317 
318  if (m_useSingleTracksAlsoForMass)
319  {
320  massVector+=massThisCluster;
321  }
322  else
323  {
324  if ( (nVTX>0 && vertexSize>1) || nVTX==0 ) {
325  massVector+=massThisCluster;
326  sumPtAdd+=ptadd;
327  }
328  }
329 
330 
331  if (m_usePtCorrectedEnergy)
332  {
333  energyFromSecondary+=std::sqrt(masswithneutrals*masswithneutrals+sumP.mag2());
334  }
335  else
336  {
337  energyFromSecondary+=std::sqrt(massThisCluster.mag2()+sumP.mag2());
338  }
339 
340 
341 
342 
343  }//end if dist<0
344  }//end vectorOfVerteces
345 
346  if (energyFromSecondary+energyFromPrimary>0) {
347  energyFraction=energyFromSecondary/(energyFromSecondary+energyFromPrimary);
348  }
349 
350  if (massVector.mag()>0) {
351  mass=std::sqrt(massVector.mag2()+sumPtAdd*sumPtAdd)+sumPtAdd;
352  mass_uncorr=massVector.mag();
353  //port range of mass to maximum 10000.
354  if (mass>5000.) {
355 mass =
356  5000.+(5000./M_PI)*2.*std::atan((M_PI/2./5000.)*(mass-5000.));
357  }
358  if (mass_uncorr>5000.) {
359 mass_uncorr =
360  5000.+(5000./M_PI)*2.*std::atan((M_PI/2./5000.)*(mass_uncorr-5000.));
361  }
362  }
363 
364  if (inverrordist!=0) {
365  significance3d=dist/std::sqrt(inverrordist);
366  //port range of significance 3d to maximum 100.
367  significance3d=100./(M_PI/2.)*std::atan((M_PI/2./100.)*significance3d);
368  }
369 
370  if (std::abs(sumPAllVertices.mag())>1e-7) {
371  deltaphi=sumPAllVertices.eta()-JetVector.Eta();
372  deltaeta=sumPAllVertices.deltaPhi(Amg::Vector3D(JetVector.Px(), JetVector.Py(), JetVector.Pz()));
373  deltaRFlightDir = std::hypot(sumPAllVertices.deltaPhi(flightAxis), sumPAllVertices.eta()-flightAxis.eta());
374  } else {
375  deltaphi=-10.;
376  deltaeta=-10.;
377  deltaRFlightDir = -10;
378  }
379 
380  // At the end, instead of writing to BTag, assign to the struct:
381  vars.massUncorr = mass_uncorr;
382  vars.chi2 = chi2;
383  vars.ndof = ndof;
384  vars.dRFlightDir = deltaRFlightDir;
385  vars.nVTX = nVTX;
386  vars.nSingleTracks = nSingleTracks;
387  vars.nTracksAtVtx = nTracksAtVtx;
388  vars.mass = mass;
389  vars.energyFraction = energyFraction;
390  vars.significance3d = significance3d;
391  vars.deltaeta = deltaeta;
392  vars.deltaphi = deltaphi;
393  return StatusCode::SUCCESS;
394 }
395 
396  StatusCode JetFitterVariablesFactory::fillJetFitterVariables(const xAOD::Jet &myJet, xAOD::BTagging* BTag, const Trk::VxJetFitterVertexInfo* myJetFitterInfo, std::string basename) const{
397 
398  //VALERIO NASTY HACK!!!!
399  bool nastyVsRevertPosToNeg = (basename.find("Flip")!=std::string::npos);
400  int nVTX(0);
401  int nTracksAtVtx(0);
402  int nSingleTracks(0);
403  float energyFraction(0);
404  float mass(0);
405  float mass_uncorr(0);
406  float significance3d(0);
407  float deltaphi(0.);
408  float deltaeta(0.);
409  float chi2(0.);
410  int ndof(0);
411  float deltaRFlightDir(0.);
412 
413  std::vector<Trk::VxJetCandidate*> myVertices;
414  Trk::VxJetCandidate* myVxJetCandidate = nullptr;
415  if (myJetFitterInfo) myVertices = myJetFitterInfo->verticesJF();
416  if(myVertices.size() == 0){
417  ATH_MSG_DEBUG("#BTAG# Trk::VxJetCandidate not found for jet fitter ");
418  fill(BTag, basename, mass_uncorr, nVTX, nSingleTracks, nTracksAtVtx, mass, energyFraction, significance3d, deltaeta, deltaphi, chi2, ndof, deltaRFlightDir);
419  return StatusCode::SUCCESS;
420  }
421  if(myVertices.size() > 0) myVxJetCandidate=dynamic_cast<Trk::VxJetCandidate*>(myVertices[0]);
422  if (myVxJetCandidate==0) {
423  ATH_MSG_WARNING("#BTAG# No correct VxJetCandidate could be retrieved." );
424  fill(BTag, basename, mass_uncorr, nVTX, nSingleTracks, nTracksAtVtx, mass, energyFraction, significance3d, deltaeta, deltaphi, chi2, ndof, deltaRFlightDir);
425  return StatusCode::SUCCESS;
426  }
427 
428  const Trk::VxJetCandidate& myJetCandidate =*myVxJetCandidate;
429 
430  const Trk::SelectedTracksInJet* mySelectedTracksInJet = myJetFitterInfo->getSelectedTracksInJet();
431 
432 
433  //put all needed information inside :-)
434 
435  //const double s_massks=ParticleConstants::KZeroMassInMeV;
436  const double s_pion=ParticleConstants::chargedPionMassInMeV;
437 
438  double energyFromPrimary=0.;
439  double energyFromSecondary=0.;
440 
441  // get fit quality variables for the PV of jetfitter
442  const Trk::VxVertexOnJetAxis* pvtxjet = myVxJetCandidate->getPrimaryVertex();
443  const Trk::FitQuality& fitquality = pvtxjet->fitQuality();
444  chi2 = fitquality.chiSquared();
445  ndof = fitquality.numberDoF();
446 
447  if (mySelectedTracksInJet!=0)
448  {
449  ATH_MSG_DEBUG(" Adding the tracks from primary vertex information ");
450  const std::vector<const Trk::ITrackLink*> & myPrimaryLinks=mySelectedTracksInJet->getPrimaryTrackLinks();
451 
452  std::vector<const Trk::ITrackLink*>::const_iterator myPrimaryLinksBegin=myPrimaryLinks.begin();
453  std::vector<const Trk::ITrackLink*>::const_iterator myPrimaryLinksEnd=myPrimaryLinks.end();
454 
455  for(std::vector<const Trk::ITrackLink*>::const_iterator myPrimaryLinksIter=myPrimaryLinksBegin;
456  myPrimaryLinksIter!=myPrimaryLinksEnd;
457  ++myPrimaryLinksIter)
458  {
459  const Trk::TrackParameters* myParameters=(*myPrimaryLinksIter)->parameters();
460  if (myParameters)
461  {
462  energyFromPrimary+=std::sqrt(s_pion*s_pion+myParameters->momentum().mag2());
463  }
464  else
465  {
466  ATH_MSG_WARNING(" no perigee in track for energy computation. Skipping primary track...");
467  }
468  }
469  } else
470  {
471  ATH_MSG_DEBUG(" No information about further primary tracks available. Normal in JetFitter vs. 1");
472  }
473 
474 
475  const Trk::RecVertexPositions & recVertexPositions=myJetCandidate.getRecVertexPositions();
476  const Amg::VectorX & vertexPosition=recVertexPositions.position();
477  const Amg::MatrixX & vertexCovMatrix = recVertexPositions.covariancePosition();
478 
479  Amg::Vector3D primaryPos(vertexPosition[Trk::jet_xv],
480  vertexPosition[Trk::jet_yv],
481  vertexPosition[Trk::jet_zv]);
482 
483  Amg::Vector3D flightAxis(1,1,1);//has to be different from 0
484  Amg::setPhi(flightAxis, vertexPosition[Trk::jet_phi]);
485  Amg::setTheta(flightAxis, vertexPosition[Trk::jet_theta]);
486 
487  xAOD::IParticle::FourMom_t JetVector = myJet.p4();
488 
489  //loop over primary vertex
490  const std::vector<Trk::VxTrackAtVertex*> & TracksAtPrimary=myJetCandidate.getPrimaryVertex()->getTracksAtVertex();
491  const std::vector<Trk::VxTrackAtVertex*>::const_iterator TracksAtPrimaryBegin=TracksAtPrimary.begin();
492  const std::vector<Trk::VxTrackAtVertex*>::const_iterator TracksAtPrimaryEnd=TracksAtPrimary.end();
493 
494  for (std::vector<Trk::VxTrackAtVertex*>::const_iterator TracksAtPrimaryIter=TracksAtPrimaryBegin;
495  TracksAtPrimaryIter!=TracksAtPrimaryEnd;
496  ++TracksAtPrimaryIter) {
497 
498  //FIXME: dynamic cast necessary? neutral perigee commented out, fix when vertex supports neutral
499  if (dynamic_cast<const Trk::Perigee*>((*TracksAtPrimaryIter)->perigeeAtVertex())!=0)
500  {
501 
502  energyFromPrimary+=
503  std::sqrt(s_pion*s_pion+
504  (*TracksAtPrimaryIter)->perigeeAtVertex()->momentum().mag2());
505  }
506  else
507  {
508  ATH_MSG_ERROR(" FIXME: VERTEX DOESN'T SUPPORT NEUTRAL PERIGEE, commented out in line 163");
509  ATH_MSG_ERROR(" Track is not a normal track neither a KS. This is an ERROR (ask developer to fix it). Skipping track... ");
510  }
511  }
512 
513 
514 
515 
516  Amg::Vector3D sumPAllVertices(0.,0.,0.);
517  CLHEP::HepLorentzVector massVector(0,0,0,0);
518  double sumPtAdd(0.);
519 
520  double dist(0.);
521  double inverrordist(0.);
522 
523  //now access the vertices on the jet axis info...
524  std::vector<Trk::VxVertexOnJetAxis*> vectorOfVertices=myJetCandidate.getVerticesOnJetAxis();
525 
526  //then you have to order them...
527  std::vector<Trk::VxVertexOnJetAxis*>::iterator vectorOfClustersOfTrackBegin=vectorOfVertices.begin();
528  std::vector<Trk::VxVertexOnJetAxis*>::iterator vectorOfClustersOfTrackEnd=vectorOfVertices.end();
529 
530  for (std::vector<Trk::VxVertexOnJetAxis*>::const_iterator vectorOfClustersOfTrackIter=vectorOfClustersOfTrackBegin;
531  vectorOfClustersOfTrackIter!=vectorOfClustersOfTrackEnd;
532  ++vectorOfClustersOfTrackIter) {
533 
534  const std::vector<Trk::VxTrackAtVertex*> & tracksOfVertex=(*vectorOfClustersOfTrackIter)->getTracksAtVertex();
535 
536  int vertexSize=tracksOfVertex.size();
537  int ntrack=(*vectorOfClustersOfTrackIter)->getNumVertex()+5;//gets the right component (should add EDM method which does
538  if (!nastyVsRevertPosToNeg)
539  {
540  if (vertexPosition[ntrack]>0) {
541  if (vertexSize>1) {
542  nVTX+=1;
543  nTracksAtVtx+=vertexSize;
544  } else {
545  nSingleTracks+=1;
546  }
547  }
548  }
549  else
550  {
551  if (vertexPosition[ntrack]<=0) {
552  if (vertexSize>1) {
553  nVTX+=1;
554  nTracksAtVtx+=vertexSize;
555  } else {
556  nSingleTracks+=1;
557  }
558  }
559  }
560  }
561 
562  for (std::vector<Trk::VxVertexOnJetAxis*>::const_iterator vectorOfClustersOfTrackIter=vectorOfClustersOfTrackBegin;
563  vectorOfClustersOfTrackIter!=vectorOfClustersOfTrackEnd;
564  ++vectorOfClustersOfTrackIter) {
565 
566  const std::vector<Trk::VxTrackAtVertex*> & tracksOfVertex=(*vectorOfClustersOfTrackIter)->getTracksAtVertex();
567  std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackBegin=tracksOfVertex.begin();
568  std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackEnd=tracksOfVertex.end();
569 
570  int vertexSize=tracksOfVertex.size();
571 
572  int ntrack=(*vectorOfClustersOfTrackIter)->getNumVertex()+5;//gets the right component (should add EDM method which does
573  // this nasty little addition of 5...)
574 
575 
576  if ((vertexPosition[ntrack]<0 && (!nastyVsRevertPosToNeg))||(vertexPosition[ntrack]>=0 && nastyVsRevertPosToNeg)) {
577  if (m_addNegativeTracksToPrimaryVertex)
578  {
579  for (std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackIter=clustersOfTrackBegin;
580  clustersOfTrackIter!=clustersOfTrackEnd;++clustersOfTrackIter) {
581 
582  energyFromPrimary+=
583  std::hypot(s_pion, (*clustersOfTrackIter)->perigeeAtVertex()->momentum().mag());
584  }
585  }
586  } else {
587 
588  if ( (nVTX>0 && vertexSize>1) || nVTX==0 ) {
589  dist+=std::abs(vertexPosition[ntrack])/vertexCovMatrix(ntrack,ntrack);
590  if (vertexCovMatrix(ntrack,ntrack)>0)
591  {
592  inverrordist+=1./vertexCovMatrix(ntrack,ntrack);
593  }
594  else
595  {
596  ATH_MSG_WARNING("The diagonal element of the vertex cov matrix ("<<ntrack<<","<<ntrack<<") is "<<vertexCovMatrix(ntrack,ntrack)<<". It should be positive... Ignoring vertex when computing L/sigma(L)");
597  }
598  }
599 
600  Amg::Vector3D sumP(0.,0.,0.);
601  CLHEP::HepLorentzVector massThisCluster(0.,0.,0.,0.);
602 
603  //in case it's a real seconday vertex track...
604 
605  for (std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackIter=clustersOfTrackBegin;
606  clustersOfTrackIter!=clustersOfTrackEnd;
607  ++clustersOfTrackIter) {
608 
609 // const Trk::MeasuredPerigee* aMeasPer=static_cast<const Trk::MeasuredPerigee*>((*clustersOfTrackIter)->perigeeAtVertex());
610  const Trk::TrackParameters* aMeasPer=(*clustersOfTrackIter)->perigeeAtVertex();
611 // CLHEP::HepVector perigeeParms = aMeasPer->parameters();
612 
613  Amg::Vector3D mytrack(aMeasPer->momentum());
614  sumP+=mytrack;
615  if (dynamic_cast<const Trk::Perigee*>((*clustersOfTrackIter)->perigeeAtVertex())!=0)
616  {
617  massThisCluster+=CLHEP::HepLorentzVector(mytrack.x(), mytrack.y(), mytrack.z(), std::hypot(s_pion, mytrack.mag()));
618  }
619  else
620  {
621  ATH_MSG_ERROR("Neutral parameter has been taken out until Vertex has been rewritten to support neutral perigee again. ");
622  ATH_MSG_ERROR(" Track is not a normal track neither a KS. This is an ERROR (ask developer to fix it). Skipping track... ");
623  }
624  }
625 
626 
627  sumPAllVertices+=sumP;
628  double ptadd=sumP.perp(flightAxis.unit());
629  double masswithneutrals=std::sqrt(massThisCluster.mag2()+ptadd*ptadd)+ptadd;
630 
631  if (m_useSingleTracksAlsoForMass)
632  {
633  massVector+=massThisCluster;
634  }
635  else
636  {
637  if ( (nVTX>0 && vertexSize>1) || nVTX==0 ) {
638  massVector+=massThisCluster;
639  sumPtAdd+=ptadd;
640  }
641  }
642 
643 
644  if (m_usePtCorrectedEnergy)
645  {
646  energyFromSecondary+=std::sqrt(masswithneutrals*masswithneutrals+sumP.mag2());
647  }
648  else
649  {
650  energyFromSecondary+=std::sqrt(massThisCluster.mag2()+sumP.mag2());
651  }
652 
653 
654 
655 
656  }//end if dist<0
657  }//end vectorOfVerteces
658 
659  if (energyFromSecondary+energyFromPrimary>0) {
660  energyFraction=energyFromSecondary/(energyFromSecondary+energyFromPrimary);
661  }
662 
663  if (massVector.mag()>0) {
664  mass=std::sqrt(massVector.mag2()+sumPtAdd*sumPtAdd)+sumPtAdd;
665  mass_uncorr=massVector.mag();
666  //port range of mass to maximum 10000.
667  if (mass>5000.) {
668  mass =
669  5000.+(5000./M_PI)*2.*std::atan((M_PI/2./5000.)*(mass-5000.));
670  }
671  if (mass_uncorr>5000.) {
672  mass_uncorr =
673  5000.+(5000./M_PI)*2.*std::atan((M_PI/2./5000.)*(mass_uncorr-5000.));
674  }
675  }
676 
677  if (inverrordist!=0) {
678  significance3d=dist/std::sqrt(inverrordist);
679  //port range of significance 3d to maximum 100.
680  significance3d=100./(M_PI/2.)*std::atan((M_PI/2./100.)*significance3d);
681  }
682 
683  if (std::abs(sumPAllVertices.mag())>1e-7) {
684  deltaphi=sumPAllVertices.eta()-JetVector.Eta();
685  deltaeta=sumPAllVertices.deltaPhi(Amg::Vector3D(JetVector.Px(), JetVector.Py(), JetVector.Pz()));
686  deltaRFlightDir = std::hypot(sumPAllVertices.deltaPhi(flightAxis), sumPAllVertices.eta()-flightAxis.eta());
687  } else {
688  deltaphi=-10.;
689  deltaeta=-10.;
690  deltaRFlightDir = -10;
691  }
692 
693 
694  fill(BTag, basename, mass_uncorr, nVTX, nSingleTracks, nTracksAtVtx, mass, energyFraction, significance3d, deltaeta, deltaphi, chi2, ndof, deltaRFlightDir);
695 
696  return StatusCode::SUCCESS;
697 
698 
699  }
700 
701 void
702 JetFitterVariablesFactory::fill(xAOD::BTagging* BTag, const std::string& basename, float mass_uncorr,
703  int nVTX, int nSingleTracks, int nTracksAtVtx, float mass, float energyFraction,
704  float significance3d, float deltaeta, float deltaphi, float chi2, int ndof, float deltaRFlightDir) const {
705 
706  BTag->setVariable<float>(basename, "massUncorr", mass_uncorr);
707  BTag->setVariable<float>(basename, "chi2", chi2);
708  BTag->setVariable<int>(basename, "ndof", ndof);
709  BTag->setVariable<float>(basename, "dRFlightDir", deltaRFlightDir);
710 
711  if (basename == "JetFitter"){
712  BTag->setTaggerInfo(nVTX, xAOD::BTagInfo::JetFitter_nVTX);
713  BTag->setTaggerInfo(nSingleTracks, xAOD::BTagInfo::JetFitter_nSingleTracks);
714  BTag->setTaggerInfo(nTracksAtVtx, xAOD::BTagInfo::JetFitter_nTracksAtVtx);
715  BTag->setTaggerInfo(mass, xAOD::BTagInfo::JetFitter_mass);
716  BTag->setTaggerInfo(energyFraction, xAOD::BTagInfo::JetFitter_energyFraction);
717  BTag->setTaggerInfo(significance3d, xAOD::BTagInfo::JetFitter_significance3d);
718  BTag->setTaggerInfo(deltaeta, xAOD::BTagInfo::JetFitter_deltaeta);
719  BTag->setTaggerInfo(deltaphi, xAOD::BTagInfo::JetFitter_deltaphi);
720  } else {
721  BTag->setVariable<int>(basename, "nVTX", nVTX);
722  BTag->setVariable<int>(basename, "nSingleTracks", nSingleTracks);
723  BTag->setVariable<int>(basename, "nTracksAtVtx", nTracksAtVtx);
724  BTag->setVariable<float>(basename, "mass", mass);
725  BTag->setVariable<float>(basename, "energyFraction", energyFraction);
726  BTag->setVariable<float>(basename, "significance3d", significance3d);
727  BTag->setVariable<float>(basename, "deltaeta", deltaeta);
728  BTag->setVariable<float>(basename, "deltaphi", deltaphi);
729  }
730 }
731 
732 }//end Analysis namespace
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
xAOD::JetFitter_deltaeta
@ JetFitter_deltaeta
JetFitter : Delta eta between jet and momentum sum of all tracks associated with displaced vertices r...
Definition: BTaggingEnums.h:56
Trk::VertexPositions::position
const Amg::VectorX & position() const
return position of vertex
Definition: VertexPositions.cxx:95
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
IDTPM::ndof
float ndof(const U &p)
Definition: TrackParametersHelper.h:134
Trk::SelectedTracksInJet
Definition: SelectedTracksInJet.h:62
Trk::VxJetCandidate::getVerticesOnJetAxis
const std::vector< VxVertexOnJetAxis * > & getVerticesOnJetAxis(void) const
Definition: VxJetCandidate.cxx:543
Amg::VectorX
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
Definition: EventPrimitives.h:30
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:27
TrackParameters.h
Analysis::JetFitterVariables::nVTX
int nVTX
Definition: IJetFitterVariablesFactory.h:36
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
Trk::VxVertexOnJetAxis::getTracksAtVertex
const std::vector< VxTrackAtVertex * > & getTracksAtVertex(void) const
get Tracks At Vertex Method
Definition: VxVertexOnJetAxis.cxx:102
Analysis::JetFitterVariablesFactory::~JetFitterVariablesFactory
virtual ~JetFitterVariablesFactory()
Destructor - check up memory allocation delete any memory allocation on the heap.
Definition: JetFitterVariablesFactory.cxx:56
Trk::jet_theta
@ jet_theta
Definition: JetVtxParamDefs.h:28
JetFitterVariablesFactory.h
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
Analysis::JetFitterVariables
Definition: IJetFitterVariablesFactory.h:31
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
xAOD::JetFitter_nSingleTracks
@ JetFitter_nSingleTracks
JetFitter : Number of single tracks.
Definition: BTaggingEnums.h:46
Analysis::JetFitterVariables::ndof
int ndof
Definition: IJetFitterVariablesFactory.h:34
VxVertexOnJetAxis.h
VxJetCandidate.h
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Trk::VxVertexOnJetAxis
VxVertexOnJetAxis inherits from Vertex.
Definition: VxVertexOnJetAxis.h:79
Trk::jet_xv
@ jet_xv
Definition: JetVtxParamDefs.h:27
SelectedTracksInJet.h
Analysis::JetFitterVariables::massUncorr
float massUncorr
Definition: IJetFitterVariablesFactory.h:32
Amg::setPhi
Amg::RotationMatrix3D setPhi(Amg::RotationMatrix3D mat, double angle, int convention=0)
Definition: EulerAnglesHelpers.h:102
drawFromPickle.atan
atan
Definition: drawFromPickle.py:36
xAOD::IParticle::FourMom_t
TLorentzVector FourMom_t
Definition of the 4-momentum type.
Definition: Event/xAOD/xAODBase/xAODBase/IParticle.h:69
Trk::VxJetCandidate::getRecVertexPositions
const Trk::RecVertexPositions & getRecVertexPositions() const
Definition: VxJetCandidate.cxx:519
xAOD::JetFitter_energyFraction
@ JetFitter_energyFraction
JetFitter : Jet efrac.
Definition: BTaggingEnums.h:52
ParticleConstants::PDG2011::chargedPionMassInMeV
constexpr double chargedPionMassInMeV
the mass of the charged pion (in MeV)
Definition: ParticleConstants.h:41
xAOD::JetFitter_nVTX
@ JetFitter_nVTX
JetFitter : Number of vertices.
Definition: BTaggingEnums.h:44
Trk::VxJetCandidate::getPrimaryVertex
const VxVertexOnJetAxis * getPrimaryVertex(void) const
Definition: VxJetCandidate.cxx:551
Analysis::JetFitterVariablesFactory::m_addNegativeTracksToPrimaryVertex
bool m_addNegativeTracksToPrimaryVertex
Definition: JetFitterVariablesFactory.h:63
Analysis::JetFitterVariables::deltaeta
float deltaeta
Definition: IJetFitterVariablesFactory.h:42
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
Analysis::JetFitterVariablesFactory::initialize
virtual StatusCode initialize()
Definition: JetFitterVariablesFactory.cxx:58
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
Analysis::JetFitterVariables::deltaphi
float deltaphi
Definition: IJetFitterVariablesFactory.h:43
Trk::SelectedTracksInJet::getPrimaryTrackLinks
const std::vector< const ITrackLink * > & getPrimaryTrackLinks() const
Get the priamry tracks (please do not delete the pointers)
Definition: SelectedTracksInJet.cxx:167
Analysis::JetFitterVariablesFactory::computeJetFitterVariables
virtual StatusCode computeJetFitterVariables(const xAOD::Jet &myJet, const Trk::VxJetFitterVertexInfo *myJetFitterInfo, const std::string &basename, JetFitterVariables &vars) const
Definition: JetFitterVariablesFactory.cxx:69
beamspotman.n
n
Definition: beamspotman.py:727
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Analysis::JetFitterVariablesFactory::fillJetFitterVariables
virtual StatusCode fillJetFitterVariables(const xAOD::Jet &, xAOD::BTagging *BTag, const Trk::VxJetFitterVertexInfo *myJetFitterInfo, std::string basename) const
Definition: JetFitterVariablesFactory.cxx:396
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition: AthCommonDataStore.h:145
Trk::jet_yv
@ jet_yv
Definition: JetVtxParamDefs.h:27
Analysis::JetFitterVariablesFactory::finalize
virtual StatusCode finalize()
Definition: JetFitterVariablesFactory.cxx:63
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:525
xAOD::JetFitter_mass
@ JetFitter_mass
JetFitter : Jet mass.
Definition: BTaggingEnums.h:50
VxTrackAtVertex.h
xAOD::JetFitter_nTracksAtVtx
@ JetFitter_nTracksAtVtx
JetFitter : Number of tracks at vertex.
Definition: BTaggingEnums.h:48
Trk::VxVertexOnJetAxis::fitQuality
const Trk::FitQuality & fitQuality() const
Fit quality access method.
Definition: VxVertexOnJetAxis.cxx:86
Analysis::JetFitterVariablesFactory::JetFitterVariablesFactory
JetFitterVariablesFactory(const std::string &name, const std::string &n, const IInterface *p)
Definition: JetFitterVariablesFactory.cxx:39
Analysis::JetFitterVariables::mass
float mass
Definition: IJetFitterVariablesFactory.h:39
Trk::FitQuality
Class to represent and store fit qualities from track reconstruction in terms of and number of degre...
Definition: FitQuality.h:97
Trk::ParametersBase
Definition: ParametersBase.h:55
Analysis::JetFitterVariables::nSingleTracks
int nSingleTracks
Definition: IJetFitterVariablesFactory.h:37
xAOD::BTagging_v1
Definition: BTagging_v1.h:39
Trk::VxJetFitterVertexInfo
Definition: VxJetFitterVertexInfo.h:58
Analysis::JetFitterVariables::nTracksAtVtx
int nTracksAtVtx
Definition: IJetFitterVariablesFactory.h:38
Analysis::JetFitterVariables::dRFlightDir
float dRFlightDir
Definition: IJetFitterVariablesFactory.h:35
xAOD::JetFitter_significance3d
@ JetFitter_significance3d
JetFitter : 3D vertex significance.
Definition: BTaggingEnums.h:54
Analysis::JetFitterVariablesFactory::m_useSingleTracksAlsoForMass
bool m_useSingleTracksAlsoForMass
Definition: JetFitterVariablesFactory.h:65
Analysis
The namespace of all packages in PhysicsAnalysis/JetTagging.
Definition: BTaggingCnvAlg.h:20
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
Amg::setTheta
void setTheta(Amg::Vector3D &v, double theta)
sets the theta of a vector without changing phi nor the magnitude
Definition: GeoPrimitivesHelpers.h:89
Analysis::JetFitterVariables::significance3d
float significance3d
Definition: IJetFitterVariablesFactory.h:41
Analysis::JetFitterVariables::energyFraction
float energyFraction
Definition: IJetFitterVariablesFactory.h:40
Trk::jet_zv
@ jet_zv
position x,y,z of primary vertex
Definition: JetVtxParamDefs.h:27
Trk::VxJetCandidate
Definition: VxJetCandidate.h:72
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Analysis::JetFitterVariablesFactory::m_usePtCorrectedEnergy
bool m_usePtCorrectedEnergy
Definition: JetFitterVariablesFactory.h:64
Trk::RecVertexPositions
Definition: RecVertexPositions.h:34
xAOD::Jet_v1
Class describing a jet.
Definition: Jet_v1.h:57
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::RecVertexPositions::covariancePosition
Amg::MatrixX const & covariancePosition() const
return the covDeltaV matrix of the vertex fit
Definition: RecVertexPositions.cxx:171
xAOD::Jet_v1::p4
virtual FourMom_t p4() const
The full 4-momentum of the particle.
Definition: Jet_v1.cxx:71
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
GeoPrimitivesHelpers.h
Analysis::JetFitterVariables::chi2
float chi2
Definition: IJetFitterVariablesFactory.h:33
xAOD::JetFitter_deltaphi
@ JetFitter_deltaphi
JetFitter : Delta phi between jet and momentum sum of all tracks associated with displaced vertices r...
Definition: BTaggingEnums.h:58
Trk::FitQuality::chiSquared
double chiSquared() const
returns the of the overall track fit
Definition: FitQuality.h:56
Trk::jet_phi
@ jet_phi
Definition: JetVtxParamDefs.h:28
Trk::FitQuality::numberDoF
int numberDoF() const
returns the number of degrees of freedom of the overall track or vertex fit as integer
Definition: FitQuality.h:60
Trk::VxJetFitterVertexInfo::verticesJF
const std::vector< Trk::VxJetCandidate * > & verticesJF() const
Definition: VxJetFitterVertexInfo.h:112
VxJetFitterVertexInfo.h
xAODType::BTag
@ BTag
The object is a b-tagging object.
Definition: ObjectType.h:60
AthAlgTool
Definition: AthAlgTool.h:26
Trk::VxJetFitterVertexInfo::getSelectedTracksInJet
const Trk::SelectedTracksInJet * getSelectedTracksInJet() const
Definition: VxJetFitterVertexInfo.h:107
Analysis::JetFitterVariablesFactory::fill
void fill(xAOD::BTagging *BTag, const std::string &basename, float mass_uncorr, int nVTX, int nSingleTracks, int nTracksAtVtx, float mass, float energyFraction, float significance3d, float deltaeta, float deltaphi, float chi2, int ndof, float deltaRFlightDir) const
Definition: JetFitterVariablesFactory.cxx:702
beamspotman.basename
basename
Definition: beamspotman.py:636
Generated on Mon Sep 29 2025 21:12:39 for ATLAS Offline Software by   doxygen 1.8.18