ATLAS Offline Software
Loading...
Searching...
No Matches
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
37namespace 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
58StatusCode JetFitterVariablesFactory::initialize() {
59 ATH_MSG_INFO(" Initialization of JetFitterVariablesFactory successful");
60 return StatusCode::SUCCESS;
61}
62
63StatusCode 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
69StatusCode 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
275if ( (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
287Amg::Vector3D sumP(0.,0.,0.);
288CLHEP::HepLorentzVector massThisCluster(0.,0.,0.,0.);
289
290//in case it's a real seconday vertex track...
291
292for (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
314sumPAllVertices+=sumP;
315double ptadd=sumP.perp(flightAxis.unit());
316double 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.) {
355mass =
356 5000.+(5000./M_PI)*2.*std::atan((M_PI/2./5000.)*(mass-5000.));
357 }
358 if (mass_uncorr>5000.) {
359mass_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
701void
702JetFitterVariablesFactory::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
M_PI
#define M_PI
Definition ActiveFraction.h:11
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition AthMsgStreamMacros.h:31
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
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
Analysis::JetFitterVariablesFactory::JetFitterVariablesFactory
JetFitterVariablesFactory(const std::string &name, const std::string &n, const IInterface *p)
Definition JetFitterVariablesFactory.cxx:39
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
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
Analysis::JetFitterVariablesFactory::~JetFitterVariablesFactory
virtual ~JetFitterVariablesFactory()
Destructor - check up memory allocation delete any memory allocation on the heap.
Definition JetFitterVariablesFactory.cxx:56
AthAlgTool::AthAlgTool
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Definition AthAlgTool.cxx:16
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition AthCommonDataStore.h:145
Trk::FitQuality
Class to represent and store fit qualities from track reconstruction in terms of and number of degre...
Definition FitQuality.h:97
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::FitQuality::chiSquared
double chiSquared() const
returns the of the overall track fit
Definition FitQuality.h:56
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:162
Trk::SelectedTracksInJet::getPrimaryTrackLinks
const std::vector< const ITrackLink * > & getPrimaryTrackLinks() const
Get the priamry tracks (please do not delete the pointers)
Definition SelectedTracksInJet.cxx:153
Trk::VertexPositions::position
const Amg::VectorX & position() const
return position of vertex
Definition VertexPositions.cxx:86
Trk::VxJetCandidate::getVerticesOnJetAxis
const std::vector< VxVertexOnJetAxis * > & getVerticesOnJetAxis(void) const
Definition VxJetCandidate.cxx:533
Trk::VxJetCandidate::getPrimaryVertex
const VxVertexOnJetAxis * getPrimaryVertex(void) const
Definition VxJetCandidate.cxx:541
Trk::VxJetCandidate::getRecVertexPositions
const Trk::RecVertexPositions & getRecVertexPositions() const
Definition VxJetCandidate.cxx:509
Trk::VxJetFitterVertexInfo::getSelectedTracksInJet
const Trk::SelectedTracksInJet * getSelectedTracksInJet() const
Definition VxJetFitterVertexInfo.h:89
Trk::VxJetFitterVertexInfo::verticesJF
const std::vector< Trk::VxJetCandidate * > & verticesJF() const
Definition VxJetFitterVertexInfo.h:94
Trk::VxVertexOnJetAxis
VxVertexOnJetAxis inherits from Vertex.
Definition VxVertexOnJetAxis.h:52
Trk::VxVertexOnJetAxis::fitQuality
const Trk::FitQuality & fitQuality() const
Fit quality access method.
Definition VxVertexOnJetAxis.cxx:72
Trk::VxVertexOnJetAxis::getTracksAtVertex
const std::vector< VxTrackAtVertex * > & getTracksAtVertex(void) const
get Tracks At Vertex Method
Definition VxVertexOnJetAxis.cxx:88
xAOD::IParticle::FourMom_t
TLorentzVector FourMom_t
Definition of the 4-momentum type.
Definition Event/xAOD/xAODBase/xAODBase/IParticle.h:69
xAOD::Jet_v1::p4
virtual FourMom_t p4() const
The full 4-momentum of the particle.
Definition Jet_v1.cxx:71
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition comparitor.cxx:525
Amg::setPhi
Amg::RotationMatrix3D setPhi(Amg::RotationMatrix3D mat, double angle, int convention=0)
Definition EulerAnglesHelpers.h:102
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition EventPrimitives.h:27
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
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
Amg::VectorX
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
Definition EventPrimitives.h:30
Analysis
The namespace of all packages in PhysicsAnalysis/JetTagging.
Definition BTaggingCnvAlg.h:20
ParticleConstants::PDG2011::chargedPionMassInMeV
constexpr double chargedPionMassInMeV
the mass of the charged pion (in MeV)
Definition ParticleConstants.h:41
Trk::Perigee
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:33
Trk::TrackParameters
ParametersBase< TrackParametersDim, Charged > TrackParameters
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:27
Trk::jet_zv
@ jet_zv
position x,y,z of primary vertex
Definition JetVtxParamDefs.h:27
xAOD::Jet
Jet_v1 Jet
Definition of the current "jet version".
Definition Event/xAOD/xAODJet/xAODJet/Jet.h:17
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
xAOD::JetFitter_nSingleTracks
@ JetFitter_nSingleTracks
JetFitter : Number of single tracks.
Definition BTaggingEnums.h:46
xAOD::JetFitter_nTracksAtVtx
@ JetFitter_nTracksAtVtx
JetFitter : Number of tracks at vertex.
Definition BTaggingEnums.h:48
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
xAOD::JetFitter_significance3d
@ JetFitter_significance3d
JetFitter : 3D vertex significance.
Definition BTaggingEnums.h:54
xAOD::JetFitter_energyFraction
@ JetFitter_energyFraction
JetFitter : Jet efrac.
Definition BTaggingEnums.h:52
xAOD::JetFitter_mass
@ JetFitter_mass
JetFitter : Jet mass.
Definition BTaggingEnums.h:50
xAOD::JetFitter_nVTX
@ JetFitter_nVTX
JetFitter : Number of vertices.
Definition BTaggingEnums.h:44
xAOD::BTagging
BTagging_v1 BTagging
Definition of the current "BTagging version".
Definition BTagging.h:17
basename
std::string basename(std::string name)
Definition utils.cxx:207
Generated on for ATLAS Offline Software by doxygen 1.14.0