ATLAS Offline Software
Loading...
Searching...
No Matches
V0Tools.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5/*********************************************************************
6 V0Tools.cxx - Description in header file
7
8 authors : Evelina Bouhova-Thacker (Lancater University)
9 email : e.bouhova@cern.ch
10
11*********************************************************************/
12
13#include "TrkVertexAnalysisUtils/V0Tools.h"
14#include "TrkParametersBase/ParametersBase.h"
15#include "TrkParticleBase/LinkToTrackParticleBase.h"
16#include "TrkParticleBase/TrackParticleBase.h"
17#include "VxVertex/VxTrackAtVertex.h"
18#include "TrkExInterfaces/IExtrapolator.h"
19#include "CLHEP/GenericFunctions/CumulativeChiSquare.hh"
20#include "xAODTracking/Vertex.h"
21#include "xAODTracking/VertexContainer.h"
22#include "CLHEP/Vector/LorentzVector.h"
23namespace Trk
24{
25
26 V0Tools::V0Tools(const std::string& t, const std::string& n, const IInterface* p) :
27 AthAlgTool(t,n,p)
28 {
29 declareInterface<V0Tools>(this);
30 }
31
32 V0Tools::~V0Tools() = default;
33
34 StatusCode V0Tools::initialize()
35 {
36 ATH_MSG_DEBUG( "Initialize successful" );
37 return StatusCode::SUCCESS;
38 }
39
40 double V0Tools::invariantMass(const xAOD::Vertex * vxCandidate, double posTrackMass, double negTrackMass) const
41 {
42 std::array<double, 2> masses = {posTrackMass, negTrackMass};
43
44 return invariantMass(vxCandidate,masses);
45 }
46
47
48
49 double V0Tools::invariantMass(const xAOD::Vertex * vxCandidate, std::span<const double> masses) const
50 {
51 double px = 0., py = 0., pz = 0., e = 0.;
52 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
53 if (masses.size() != NTrk) {
54 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
55 return -999999.;
56 }
57 for( unsigned int it=0; it<NTrk; it++) {
58 if (masses[it] >= 0.) {
59 xAOD::TrackParticle::FourMom_t lorentz_trk = track4Momentum(vxCandidate,it,masses[it]);
60 px += lorentz_trk.Px();
61 py += lorentz_trk.Py();
62 pz += lorentz_trk.Pz();
63 e += lorentz_trk.E();
64 }
65 }
66 double msq = e*e - px*px - py*py - pz*pz;
67 return (msq>0.) ? sqrt(msq) : 0.;
68 }
69
70 double V0Tools::invariantMassError(const xAOD::Vertex * vxCandidate, double posTrackMass, double negTrackMass) const
71 {
72 std::array<double, 2> masses = {posTrackMass, negTrackMass};
73
74 return invariantMassError(vxCandidate,masses);
75 }
76
77 double V0Tools::invariantMassError(const xAOD::Vertex * vxCandidate, std::span<const double> masses) const
78 {
79 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
80 if (masses.size() != NTrk) {
81 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
82 return -999999.;
83 }
84 double error = -999999.;
85 auto fullCov = convertCovMatrix(vxCandidate);
86 if (fullCov.size() == 0) {
87 ATH_MSG_DEBUG("0 pointer for full covariance. Making-up one from the vertex and tracks covariances");
88 error = massErrorVxCandidate(vxCandidate,masses);
89 } else {
90 unsigned int ndim = fullCov.rows();
91 if (ndim != 0) {
92 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
93 error = massErrorV0Fitter(vxCandidate,masses);
94 } else if (ndim == 3*NTrk+3) {
95 error = massErrorVKalVrt(vxCandidate,masses);
96 }
97 }
98 }
99 return error;
100 }
101
102 double V0Tools::massErrorV0Fitter(const xAOD::Vertex * vxCandidate, double posTrackMass, double negTrackMass) const
103 {
104 std::array<double, 2> masses = {posTrackMass, negTrackMass};
105
106 return massErrorV0Fitter(vxCandidate,masses);
107 }
108
109 double V0Tools::massErrorV0Fitter(const xAOD::Vertex * vxCandidate, std::span<const double> masses) const
110 {
111 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
112 if (masses.size() != NTrk) {
113 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
114 return -999999.;
115 }
116 auto fullCov = convertCovMatrix(vxCandidate);
117 if (fullCov.size() == 0) return -999999.;
118 double E=0., Px=0., Py=0., Pz=0.;
119 std::vector<double>phi(NTrk), theta(NTrk), qOverP(NTrk), charge(NTrk), e(NTrk);
120 std::vector<double>dm2dphi(NTrk), dm2dtheta(NTrk), dm2dqOverP(NTrk);
121 for( unsigned int it=0; it<NTrk; it++) {
122 if (masses[it] >= 0.) {
123 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
124 double trkCharge = 1.;
125 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
126 charge[it] = trkCharge;
127 phi[it] = bPer->parameters()[Trk::phi];
128 theta[it] = bPer->parameters()[Trk::theta];
129 qOverP[it] = bPer->parameters()[Trk::qOverP];
130 double tmp = 1./(qOverP[it]*qOverP[it]) + masses[it]*masses[it];
131 double pe = (tmp>0.) ? sqrt(tmp) : 0.;
132 e[it] = pe;
133 E += e[it];
134 Px += bPer->momentum()[Trk::px];
135 Py += bPer->momentum()[Trk::py];
136 Pz += bPer->momentum()[Trk::pz];
137 }
138 }
139 double msq = E*E - Px*Px - Py*Py - Pz*Pz;
140 double mass = (msq>0.) ? sqrt(msq) : 0.;
141
142 for( unsigned int it=0; it<NTrk; it++) {
143 if (masses[it] >= 0.) {
144 dm2dphi[it] = 2.*(Px*sin(phi[it])-Py*cos(phi[it]))*sin(theta[it])*charge[it]/qOverP[it];
145 dm2dtheta[it] = 2.*(Pz*sin(theta[it])-(Px*cos(phi[it])+Py*sin(phi[it]))*cos(theta[it]))*charge[it]/qOverP[it];
146 dm2dqOverP[it] = 2.*(Pz*cos(theta[it])+(Px*cos(phi[it])+Py*sin(phi[it]))*sin(theta[it])-E*charge[it]/(qOverP[it]*e[it]))*charge[it]/(qOverP[it]*qOverP[it]);
147 }
148 }
149
150 Amg::MatrixX D_vec(5*NTrk,1); D_vec.setZero();
151 for( unsigned int it=0; it<NTrk; it++) {
152 D_vec(5*it+0,0) = 0.;
153 D_vec(5*it+1,0) = 0.;
154 D_vec(5*it+2,0) = dm2dphi[it];
155 D_vec(5*it+3,0) = dm2dtheta[it];
156 D_vec(5*it+4,0) = dm2dqOverP[it];
157 }
158 Amg::MatrixX V0_merr = D_vec.transpose() * fullCov.block(0,0,5*NTrk,5*NTrk) * D_vec;
159
160 double massVarsq = V0_merr(0,0);
161 if (massVarsq <= 0.) ATH_MSG_DEBUG("massError: negative sqrt massVarsq " << massVarsq);
162 double massVar = (massVarsq>0.) ? sqrt(massVarsq) : 0.;
163 double massErr = massVar/(2.*mass);
164 return massErr;
165 }
166
167 double V0Tools::massErrorVKalVrt(const xAOD::Vertex * vxCandidate, double posTrackMass, double negTrackMass) const
168 {
169 std::array<double, 2> masses = {posTrackMass, negTrackMass};
170 return massErrorVKalVrt(vxCandidate,masses);
171 }
172
173 double V0Tools::massErrorVKalVrt(const xAOD::Vertex * vxCandidate, std::span<const double> masses) const
174 {
175 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
176 if (masses.size() != NTrk) {
177 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
178 return -999999.;
179 }
180 std::vector<xAOD::TrackParticle::FourMom_t> particleMom(NTrk);
181 std::vector<AmgMatrix(3,3)> particleDeriv(NTrk);
182 xAOD::TrackParticle::FourMom_t totalMom;
183 AmgMatrix(3,3) tmpDeriv; tmpDeriv.setZero();
184 auto fullCov = convertCovMatrix(vxCandidate);
185 if (fullCov.size() == 0) return -999999.;
186
187 for( unsigned int it=0; it<NTrk; it++){
188 if (masses[it] >= 0.) {
189 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
190 double phi = bPer->parameters()[Trk::phi];
191 double theta = bPer->parameters()[Trk::theta];
192 double invP = bPer->parameters()[Trk::qOverP];
193 double px = cos(phi)*sin(theta)/fabs(invP);
194 double py = sin(phi)*sin(theta)/fabs(invP);
195 double pz = cos(theta)/fabs(invP);
196 double esq = px*px + py*py + pz*pz + masses[it]*masses[it];
197 double e = (esq>0.) ? sqrt(esq) : 0.;
198 xAOD::TrackParticle::FourMom_t tmp(0,0,0,0);
199 tmp.SetPxPyPzE(px,py,pz,e);
200 particleMom[it] = tmp;
201 totalMom += tmp;
202
203// d(Px,Py,Pz)/d(Phi,Theta,InvP)
204 tmpDeriv(0,0) = - tmp.Py();
205 tmpDeriv(1,0) = tmp.Px();
206 tmpDeriv(2,0) = 0.;
207 tmpDeriv(0,1) = cos(phi) * tmp.Pz();
208 tmpDeriv(1,1) = sin(phi) * tmp.Pz();
209 tmpDeriv(2,1) = - sin(theta)/fabs(invP);
210 tmpDeriv(0,2) = - tmp.Px()/invP;
211 tmpDeriv(1,2) = - tmp.Py()/invP;
212 tmpDeriv(2,2) = - tmp.Pz()/invP;
213 particleDeriv[it] = tmpDeriv;
214 }
215 }
216
217 std::vector<double> Deriv(3*NTrk+3, 0.);
218 for(unsigned int it=0; it<NTrk; it++){
219 if (masses[it] >= 0.) {
220 double dMdPx = ( totalMom.E() * particleMom[it].Px()/particleMom[it].E() - totalMom.Px() ) / totalMom.M();
221 double dMdPy = ( totalMom.E() * particleMom[it].Py()/particleMom[it].E() - totalMom.Py() ) / totalMom.M();
222 double dMdPz = ( totalMom.E() * particleMom[it].Pz()/particleMom[it].E() - totalMom.Pz() ) / totalMom.M();
223
224 double dMdPhi = dMdPx*particleDeriv[it](0,0) + dMdPy*particleDeriv[it](1,0) + dMdPz*particleDeriv[it](2,0);
225 double dMdTheta = dMdPx*particleDeriv[it](0,1) + dMdPy*particleDeriv[it](1,1) + dMdPz*particleDeriv[it](2,1);
226 double dMdInvP = dMdPx*particleDeriv[it](0,2) + dMdPy*particleDeriv[it](1,2) + dMdPz*particleDeriv[it](2,2);
227
228 Deriv[3*it + 3 + 0] = dMdPhi; Deriv[3*it + 3 + 1] = dMdTheta; Deriv[3*it + 3 + 2] = dMdInvP;
229 }
230 }
231
232 double err = 0;
233 for(unsigned int i=0; i<3*NTrk+3; i++){
234 for(unsigned int j=0; j<3*NTrk+3; j++){
235 err += Deriv[i]*( fullCov)(i,j)*Deriv[j];
236 }
237 }
238 if (err <= 0.) ATH_MSG_DEBUG("massError: negative sqrt err " << err);
239 return (err>0.) ? sqrt(err) : 0.;
240 }
241
242 double V0Tools::massErrorVxCandidate(const xAOD::Vertex * vxCandidate, double posTrackMass, double negTrackMass) const
243 {
244 std::array<double, 2> masses = {posTrackMass, negTrackMass};
245
246 return massErrorVxCandidate(vxCandidate,masses);
247 }
248
249 double V0Tools::massErrorVxCandidate(const xAOD::Vertex * vxCandidate, std::span<const double> masses) const
250 {
251 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
252 if (masses.size() != NTrk) {
253 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
254 return -999999.;
255 }
256 double E=0., Px=0., Py=0., Pz=0.;
257 std::vector<double>phi(NTrk), theta(NTrk), qOverP(NTrk), charge(NTrk), e(NTrk);
258 std::vector<double>dm2dphi(NTrk), dm2dtheta(NTrk), dm2dqOverP(NTrk);
259 Amg::MatrixX V0_cor(5*NTrk,5*NTrk); V0_cor.setZero();
260 for( unsigned int it=0; it<NTrk; it++) {
261 if (masses[it] >= 0.) {
262 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
263 const AmgSymMatrix(5)* cov_tmp = bPer->covariance();
264 V0_cor(5*it+2,5*it+2) = (*cov_tmp)(2,2);
265 V0_cor(5*it+2,5*it+3) = (*cov_tmp)(2,3);
266 V0_cor(5*it+2,5*it+4) = (*cov_tmp)(2,4);
267 V0_cor(5*it+3,5*it+3) = (*cov_tmp)(3,3);
268 V0_cor(5*it+3,5*it+4) = (*cov_tmp)(3,4);
269 V0_cor(5*it+4,5*it+4) = (*cov_tmp)(4,4);
270 V0_cor(5*it+3,5*it+2) = (*cov_tmp)(2,3);
271 V0_cor(5*it+4,5*it+2) = (*cov_tmp)(2,4);
272 V0_cor(5*it+4,5*it+3) = (*cov_tmp)(3,4);
273 double trkCharge = 1.;
274 if (bPer->parameters()(Trk::qOverP) < 0.) trkCharge = -1.;
275 charge[it] = trkCharge;
276 phi[it] = bPer->parameters()(Trk::phi);
277 theta[it] = bPer->parameters()(Trk::theta);
278 qOverP[it] = bPer->parameters()(Trk::qOverP);
279 double tmp = 1./(qOverP[it]*qOverP[it]) + masses[it]*masses[it];
280 double pe = (tmp>0.) ? sqrt(tmp) : 0.;
281 e[it] = pe;
282 E += e[it];
283 Px += bPer->momentum()(Trk::px);
284 Py += bPer->momentum()(Trk::py);
285 Pz += bPer->momentum()(Trk::pz);
286 }
287 }
288 double msq = E*E - Px*Px - Py*Py - Pz*Pz;
289 double mass = (msq>0.) ? sqrt(msq) : 0.;
290
291 for( unsigned int it=0; it<NTrk; it++) {
292 if (masses[it] >= 0.) {
293 dm2dphi[it] = 2.*(Px*sin(phi[it])-Py*cos(phi[it]))*sin(theta[it])*charge[it]/qOverP[it];
294 dm2dtheta[it] = 2.*(Pz*sin(theta[it])-(Px*cos(phi[it])+Py*sin(phi[it]))*cos(theta[it]))*charge[it]/qOverP[it];
295 dm2dqOverP[it] = 2.*(Pz*cos(theta[it])+(Px*cos(phi[it])+Py*sin(phi[it]))*sin(theta[it])-E*charge[it]/(qOverP[it]*e[it]))*charge[it]/(qOverP[it]*qOverP[it]);
296 }
297 }
298
299 Amg::MatrixX D_vec(5*NTrk,1); D_vec.setZero();
300 for( unsigned int it=0; it<NTrk; it++) {
301 D_vec(5*it+0,0) = 0.;
302 D_vec(5*it+1,0) = 0.;
303 D_vec(5*it+2,0) = dm2dphi[it];
304 D_vec(5*it+3,0) = dm2dtheta[it];
305 D_vec(5*it+4,0) = dm2dqOverP[it];
306 }
307
308 Amg::MatrixX V0_merr = D_vec.transpose() * V0_cor * D_vec;
309
310 double massVarsq = V0_merr(0,0);
311 if (massVarsq <= 0.) ATH_MSG_DEBUG("massError: negative sqrt massVarsq " << massVarsq);
312 double massVar = (massVarsq>0.) ? sqrt(massVarsq) : 0.;
313 double massErr = massVar/(2.*mass);
314 return massErr;
315 }
316
317 double V0Tools::invariantMassProbability(const xAOD::Vertex * vxCandidate, double V0Mass, double posTrackMass, double negTrackMass) const
318 {
319 std::array<double, 2> masses = {posTrackMass , negTrackMass};
320
321 return invariantMassProbability(vxCandidate,V0Mass,masses);
322 }
323
324 double V0Tools::invariantMassProbability(const xAOD::Vertex * vxCandidate, double V0Mass, std::span<const double> masses) const
325 {
326 double mass = invariantMass(vxCandidate, masses);
327 double massErr = invariantMassError(vxCandidate, masses);
328 if(massErr > 0.)
329 {
330 double chi2 = (V0Mass - mass)*(V0Mass - mass)/(massErr*massErr);
331 int ndf = 1;
332 Genfun::CumulativeChiSquare myCumulativeChiSquare(ndf);
333 if (chi2 > 0.) {
334 double achi2prob = 1.-myCumulativeChiSquare(chi2);
335 return achi2prob;
336 } else {
337 ATH_MSG_DEBUG("chi <= 0");
338 return -1.;
339 }
340 }
341 else {
342 return -1.;
343 }
344 }
345
346 double V0Tools::massProbability(double V0Mass, double mass, double massErr) const
347 {
348 if(massErr > 0.)
349 {
350 double chi2 = (V0Mass - mass)*(V0Mass - mass)/(massErr*massErr);
351 int ndf = 1;
352 Genfun::CumulativeChiSquare myCumulativeChiSquare(ndf);
353 if (chi2 > 0.) {
354 double achi2prob = 1.-myCumulativeChiSquare(chi2);
355 return achi2prob;
356 } else {
357 ATH_MSG_DEBUG("chi <= 0");
358 return -1.;
359 }
360 }
361 else {
362 return -1.;
363 }
364 }
365
366 Amg::Vector3D V0Tools::trackMomentum(const xAOD::Vertex * vxCandidate, unsigned int trkIndex)
367 {
368 const Trk::TrackParameters* aPerigee = vxCandidate->vxTrackAtVertex()[trkIndex].perigeeAtVertex();
369 return aPerigee->momentum();
370 }
371
372 Amg::Vector3D V0Tools::positiveTrackMomentum(const xAOD::Vertex * vxCandidate)
373 {
374 Amg::Vector3D mom; mom.setZero();
375 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
376 if (NTrk != 2) return mom;
377 for( unsigned int it=0; it<NTrk; it++) {
378 float trkCharge = vxCandidate->trackParticle(it)->charge();
379 if (trkCharge > 0) mom = trackMomentum(vxCandidate,it);
380 }
381 return mom;
382 }
383
384 Amg::Vector3D V0Tools::negativeTrackMomentum(const xAOD::Vertex * vxCandidate)
385 {
386 Amg::Vector3D mom; mom.setZero();
387 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
388 if (NTrk != 2) return mom;
389 for( unsigned int it=0; it<NTrk; it++) {
390 float trkCharge = vxCandidate->trackParticle(it)->charge();
391 if (trkCharge < 0) mom = trackMomentum(vxCandidate,it);
392 }
393 return mom;
394 }
395
396 Amg::Vector3D V0Tools::V0Momentum(const xAOD::Vertex * vxCandidate)
397 {
398 Amg::Vector3D mom; mom.setZero();
399 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
400 for( unsigned int it=0; it<NTrk; it++) {
401 mom += trackMomentum(vxCandidate,it);
402 }
403 return mom;
404 }
405
406 xAOD::TrackParticle::FourMom_t V0Tools::track4Momentum(const xAOD::Vertex * vxCandidate, unsigned int trkIndex, double mass)
407 {
408 Amg::Vector3D mom = trackMomentum(vxCandidate, trkIndex);
409 xAOD::TrackParticle::FourMom_t lorentz(0,0,0,0);
410 double tmp = mass*mass + mom.x()*mom.x() + mom.y()*mom.y() + mom.z()*mom.z();
411 double e = (tmp>0.) ? sqrt(tmp) : 0.;
412 lorentz.SetPxPyPzE(mom.x(), mom.y(), mom.z(), e);
413 return lorentz;
414 }
415
416 xAOD::TrackParticle::FourMom_t V0Tools::positiveTrack4Momentum(const xAOD::Vertex * vxCandidate, double mass)
417 {
418 Amg::Vector3D mom = positiveTrackMomentum(vxCandidate);
419 xAOD::TrackParticle::FourMom_t lorentz(0,0,0,0);
420 double tmp = mass*mass + mom.x()*mom.x() + mom.y()*mom.y() + mom.z()*mom.z();
421 double e = (tmp>0.) ? sqrt(tmp) : 0.;
422 lorentz.SetPxPyPzE(mom.x(), mom.y(), mom.z(), e);
423 return lorentz;
424 }
425
426 xAOD::TrackParticle::FourMom_t V0Tools::negativeTrack4Momentum(const xAOD::Vertex * vxCandidate, double mass)
427 {
428 Amg::Vector3D mom = negativeTrackMomentum(vxCandidate);
429 xAOD::TrackParticle::FourMom_t lorentz(0,0,0,0);
430 double tmp = mass*mass + mom.x()*mom.x() + mom.y()*mom.y() + mom.z()*mom.z();
431 double e = (tmp>0.) ? sqrt(tmp) : 0.;
432 lorentz.SetPxPyPzE(mom.x(), mom.y(), mom.z(), e);
433 return lorentz;
434 }
435
436 xAOD::TrackParticle::FourMom_t V0Tools::V04Momentum(const xAOD::Vertex * vxCandidate, double V0Mass)
437 {
438 Amg::Vector3D mom = V0Momentum(vxCandidate);
439 double tmp = V0Mass*V0Mass + mom.x()*mom.x() + mom.y()*mom.y() + mom.z()*mom.z();
440 double e = (tmp>0.) ? sqrt(tmp) : 0.;
441 xAOD::TrackParticle::FourMom_t lorentz(0,0,0,0);
442 lorentz.SetPxPyPzE(mom.x(), mom.y(), mom.z(), e);
443 return lorentz;
444 }
445
446 float V0Tools::ndof(const xAOD::Vertex * vxCandidate)
447 {
448 return vxCandidate->numberDoF();
449 }
450
451 float V0Tools::chisq(const xAOD::Vertex * vxCandidate)
452 {
453 return vxCandidate->chiSquared();
454 }
455
456 double V0Tools::vertexProbability(const xAOD::Vertex * vxCandidate) const
457 {
458 float dof = ndof(vxCandidate);
459 if (dof > 0.) {
460 Genfun::CumulativeChiSquare myCumulativeChiSquare(dof);
461 float chi = chisq(vxCandidate);
462 if (chi > 0.) {
463 double chi2prob = 1.-myCumulativeChiSquare(chi);
464 return chi2prob;
465 } else {
466 ATH_MSG_DEBUG("chi <= 0");
467 return -1.;
468 }
469 } else {
470 ATH_MSG_DEBUG("dof <= 0");
471 return -1.;
472 }
473
474 }
475
476 Amg::Vector3D V0Tools::vtx(const xAOD::Vertex * vxCandidate)
477 {
478 return vxCandidate->position();
479 }
480
481 double V0Tools::rxy(const xAOD::Vertex * vxCandidate)
482 {
483 return vxCandidate->position().perp();
484 }
485
486 double V0Tools::rxy(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex)
487 {
488 return (vxCandidate->position() - vertex->position()).perp();
489 }
490
491 double V0Tools::rxy(const xAOD::Vertex * vxCandidate, const Amg::Vector3D& vertex)
492 {
493 return (vxCandidate->position() - vertex).perp();
494 }
495
496 double V0Tools::rxy_var(double dx, double dy, const Amg::MatrixX& cov)
497 {
498 double rxysq = dx*dx + dy*dy;
499 double rxy = (rxysq>0.) ? sqrt(rxysq) : 0.;
500 double drdx = dx/rxy;
501 double drdy = dy/rxy;
502 AmgMatrix(2, 1) D_vec; D_vec.setZero();
503 D_vec(0,0) = drdx;
504 D_vec(1,0) = drdy;
505 Amg::MatrixX rxy_err = D_vec.transpose() * cov.block<2,2>(0,0) * D_vec;
506 double rxyVar = rxy_err(0,0);
507 return rxyVar;
508 }
509
510 double V0Tools::rxyError(const xAOD::Vertex * vxCandidate) const
511 {
512 const Amg::MatrixX& cov = vxCandidate->covariancePosition();
513 double dx = vxCandidate->position().x();
514 double dy = vxCandidate->position().y();
515 double rxyVar = rxy_var(dx,dy,cov);
516 if (rxyVar <= 0.) ATH_MSG_DEBUG("rxyError: negative sqrt rxyVar " << rxyVar);
517 return (rxyVar>0.) ? sqrt(rxyVar) : 0.;
518 }
519
520 double V0Tools::rxyError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex) const
521 {
522 const Amg::MatrixX cov = vxCandidate->covariancePosition() + vertex->covariancePosition();
523 auto vert = vxCandidate->position() - vertex->position();
524 double dx = vert.x();
525 double dy = vert.y();
526 double rxyVar = rxy_var(dx,dy,cov);
527 if (rxyVar <= 0.) ATH_MSG_DEBUG("rxyError: negative sqrt rxyVar " << rxyVar);
528 return (rxyVar>0.) ? sqrt(rxyVar) : 0.;
529 }
530
531 double V0Tools::rxyError(const xAOD::Vertex * vxCandidate, const Amg::Vector3D& vertex) const
532 {
533 const Amg::MatrixX& cov = vxCandidate->covariancePosition();
534 auto vert = vxCandidate->position() - vertex;
535 double dx = vert.x();
536 double dy = vert.y();
537 double rxyVar = rxy_var(dx,dy,cov);
538 if (rxyVar <= 0.) ATH_MSG_DEBUG("rxyError: negative sqrt rxyVar " << rxyVar);
539 return (rxyVar>0.) ? sqrt(rxyVar) : 0.;
540 }
541
542 double V0Tools::pT (const xAOD::Vertex * vxCandidate)
543 {
544 return V0Momentum(vxCandidate).perp();
545 }
546
547 double V0Tools::pTError(const xAOD::Vertex * vxCandidate) const
548 {
549 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
550 double Px=0., Py=0.;
551 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
552 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
553 std::vector<double>dPTdqOverP(NTrk), dPTdtheta(NTrk), dPTdphi(NTrk);
554
555 auto fullCov = convertCovMatrix(vxCandidate);
556 for( unsigned int it=0; it<NTrk; it++) {
557 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
558 double trkCharge = 1.;
559 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
560 double phi = bPer->parameters()[Trk::phi];
561 double theta = bPer->parameters()[Trk::theta];
562 double qOverP = bPer->parameters()[Trk::qOverP];
563 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
564 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
565 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
566 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
567 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
568 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
569 Px += bPer->momentum()[Trk::px];
570 Py += bPer->momentum()[Trk::py];
571 }
572 double PTsq = Px*Px+Py*Py;
573 double PT = (PTsq>0.) ? sqrt(PTsq) : 0.;
574
575 for( unsigned int it=0; it<NTrk; it++) {
576 dPTdqOverP[it] = (Px*dpxdqOverP[it]+Py*dpydqOverP[it])/PT;
577 dPTdtheta[it] = (Px*dpxdtheta[it]+Py*dpydtheta[it])/PT;
578 dPTdphi[it] = (Px*dpxdphi[it]+Py*dpydphi[it])/PT;
579 }
580
581 unsigned int ndim = 0;
582 if (fullCov.size() == 0) {
583 ndim = 5*NTrk+3;
584 } else {
585 ndim = fullCov.rows();
586 }
587
588 Amg::MatrixX PtErrSq;
589 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
590 Amg::MatrixX D_vec(5*NTrk,1); D_vec.setZero();
591 for( unsigned int it=0; it<NTrk; it++) {
592 D_vec(5*it+0,0) = 0.;
593 D_vec(5*it+1,0) = 0.;
594 D_vec(5*it+2,0) = dPTdphi[it];
595 D_vec(5*it+3,0) = dPTdtheta[it];
596 D_vec(5*it+4,0) = dPTdqOverP[it];
597 }
598 if (fullCov.size() == 0) {
599 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
600 PtErrSq = D_vec.transpose() * V0_cov * D_vec;
601 } else {
602 PtErrSq = D_vec.transpose() * fullCov.block(0,0,5*NTrk, 5*NTrk) * D_vec;
603 }
604 } else if (ndim == 3*NTrk+3) {
605 Amg::MatrixX D_vec(3*NTrk,1); D_vec.setZero();
606 for( unsigned int it=0; it<NTrk; it++) {
607 D_vec(3*it+0,0) = dPTdphi[it];
608 D_vec(3*it+1,0) = dPTdtheta[it];
609 D_vec(3*it+2,0) = dPTdqOverP[it];
610 }
611 PtErrSq = D_vec.transpose() * fullCov.block(3,3,3*NTrk,3*NTrk) * D_vec;
612 }else{
613 ATH_MSG_ERROR("This should not happen - returning zero.");
614 return 0;
615 }
616
617 double PtErrsq = PtErrSq(0,0);
618 if (PtErrsq <= 0.) ATH_MSG_DEBUG("ptError: negative sqrt PtErrsq " << PtErrsq);
619 return (PtErrsq>0.) ? sqrt(PtErrsq) : 0.;
620 }
621
622 Amg::Vector3D V0Tools::pca(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex)
623 {
624 assert(vxCandidate!=0);
625 if(nullptr == vxCandidate) {
626 Amg::Vector3D p; p.setZero();
627 return p;
628 }
629 const Amg::Vector3D& pv = vertex->position();
630 Amg::Vector3D sv = vtx(vxCandidate);
631 Amg::Vector3D P = V0Momentum(vxCandidate);
632 double p2 = P.mag2();
633 double pdr = P.dot((sv - pv));
634 return sv - P*pdr/p2;
635 }
636
637 double V0Tools::separation(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex) const
638 {
639 const Amg::MatrixX& cov = (vxCandidate->covariancePosition() + vertex->covariancePosition()).inverse().eval();
640 Amg::Vector3D D_vec; D_vec.setZero();
641 auto vert = vxCandidate->position() - vertex->position();
642 D_vec(0) = vert.x();
643 D_vec(1) = vert.y();
644 D_vec(2) = vert.z();
645 Amg::MatrixX sepVarsqMat = D_vec.transpose() * cov * D_vec;
646 double sepVarsq = sepVarsqMat(0,0);
647 if (sepVarsq <= 0.) ATH_MSG_DEBUG("separation: negative sqrt sepVarsq " << sepVarsq);
648 double sepVar = (sepVarsq>0.) ? sqrt(sepVarsq) : 0.;
649 return sepVar;
650 }
651
652 double V0Tools::separation(const xAOD::Vertex * vxCandidate, const Amg::Vector3D& vertex) const
653 {
654 const Amg::SymMatrixX& cov = vxCandidate->covariancePosition().inverse().eval();
655 Amg::Vector3D D_vec; D_vec.setZero();
656 auto vert = vxCandidate->position() - vertex;
657 D_vec(0) = vert.x();
658 D_vec(1) = vert.y();
659 D_vec(2) = vert.z();
660 Amg::MatrixX sepVarsqMat = D_vec.transpose() * cov * D_vec;
661 double sepVarsq = sepVarsqMat(0,0);
662 if (sepVarsq <= 0.) ATH_MSG_DEBUG("separation: negative sqrt sepVarsq " << sepVarsq);
663 double sepVar = (sepVarsq>0.) ? sqrt(sepVarsq) : 0.;
664 return sepVar;
665 }
666
667 double V0Tools::a0xy(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex)
668 {
669 double cosineTheta_xy = cosTheta_xy(vxCandidate,vertex);
670 double sinTheta_xy = ((1.-cosineTheta_xy*cosineTheta_xy)>0.) ? sqrt((1.-cosineTheta_xy*cosineTheta_xy)) : 0.;
671 return (vtx(vxCandidate)-vertex->position()).perp() * sinTheta_xy;
672 //return (vtx(vxCandidate)-vertex->position()).perp() * sqrt(1.-cosineTheta_xy*cosineTheta_xy);
673 }
674
675 double V0Tools::a0z(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex)
676 {
677 Amg::Vector3D pv = vertex->position();
678 Amg::Vector3D ca_point = pca(vxCandidate,vertex);
679 Amg::Vector3D a0_vec = pv - ca_point;
680 return a0_vec.z();
681 }
682
683 double V0Tools::a0(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex)
684 {
685 double cosineTheta = cosTheta(vxCandidate,vertex);
686 double sinTheta = ((1.-cosineTheta*cosineTheta)>0.) ? sqrt((1.-cosineTheta*cosineTheta)) : 0.;
687 return (vtx(vxCandidate)-vertex->position()).mag() * sinTheta;
688 //return (vtx(vxCandidate)-vertex->position()).mag() * sqrt(1.-cosineTheta*cosineTheta);
689 }
690
691 double V0Tools::a0zError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex) const
692 {
693 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
694 auto vert = vxCandidate->position() - vertex->position();
695 double dx = vert.x();
696 double dy = vert.y();
697 double dz = vert.z();
698 double Px=0., Py=0., Pz=0.;
699 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
700 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
701 std::vector<double>dpzdqOverP(NTrk), dpzdtheta(NTrk);
702 std::vector<double>da0dqOverP(NTrk), da0dtheta(NTrk), da0dphi(NTrk);
703
704 auto fullCov = convertCovMatrix(vxCandidate);
705 for( unsigned int it=0; it<NTrk; it++) {
706 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
707 double trkCharge = 1.;
708 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
709 double phi = bPer->parameters()[Trk::phi];
710 double theta = bPer->parameters()[Trk::theta];
711 double qOverP = bPer->parameters()[Trk::qOverP];
712 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
713 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
714 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
715 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
716 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
717 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
718 dpzdqOverP[it] = -(cos(theta)*trkCharge)/(qOverP*qOverP);
719 dpzdtheta[it] = -(sin(theta)*trkCharge)/qOverP;
720 Px += bPer->momentum()[Trk::px];
721 Py += bPer->momentum()[Trk::py];
722 Pz += bPer->momentum()[Trk::pz];
723 }
724 double P2 = Px*Px+Py*Py+Pz*Pz;
725 double B = Px*dx+Py*dy+Pz*dz;
726 if (P2 == 0.)[[unlikely]]{
727 ATH_MSG_ERROR("a0zError: Divisor is zero - returning zero.");
728 return 0.;
729 }
730 double da0dx = (Px*Pz)/P2;
731 double da0dy = (Py*Pz)/P2;
732 double da0dz = (Pz*Pz)/P2 - 1.;
733 double da0dx0 = -da0dx;
734 double da0dy0 = -da0dy;
735 double da0dz0 = -da0dz;
736 for( unsigned int it=0; it<NTrk; it++) {
737 double dP2dqOverP = 2.*(Px*dpxdqOverP[it]+Py*dpydqOverP[it]+Pz*dpzdqOverP[it]);
738 double dP2dtheta = 2.*(Px*dpxdtheta[it]+Py*dpydtheta[it]+Pz*dpzdtheta[it]);
739 double dP2dphi = 2.*(Px*dpxdphi[it]+Py*dpydphi[it]);
740 da0dqOverP[it] = (B*(P2*dpzdqOverP[it]-Pz*dP2dqOverP) +
741 Pz*P2*(dx*dpxdqOverP[it]+dy*dpydqOverP[it]+dz*dpzdqOverP[it]))/(P2*P2);
742 da0dtheta[it] = (B*(P2*dpzdtheta[it]-Pz*dP2dtheta) +
743 Pz*P2*(dx*dpxdtheta[it]+dy*dpydtheta[it]+dz*dpzdtheta[it]))/(P2*P2);
744 da0dphi[it] = -(B*Pz*dP2dphi -
745 Pz*P2*(dx*dpxdphi[it]+dy*dpydphi[it]))/(P2*P2);
746 }
747
748 unsigned int ndim = 0;
749 if (fullCov.size() != 0) {
750 ndim = fullCov.rows();
751 } else {
752 ndim = 5*NTrk+3;
753 }
754
755 Amg::MatrixX V0_err;
756 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
757 Amg::MatrixX D_vec(5*NTrk+6,1); D_vec.setZero();
758 for( unsigned int it=0; it<NTrk; it++) {
759 D_vec(5*it+0) = 0.;
760 D_vec(5*it+1) = 0.;
761 D_vec(5*it+2) = da0dphi[it];
762 D_vec(5*it+3) = da0dtheta[it];
763 D_vec(5*it+4) = da0dqOverP[it];
764 }
765 D_vec(5*NTrk+0) = da0dx;
766 D_vec(5*NTrk+1) = da0dy;
767 D_vec(5*NTrk+2) = da0dz;
768 D_vec(5*NTrk+3) = da0dx0;
769 D_vec(5*NTrk+4) = da0dy0;
770 D_vec(5*NTrk+5) = da0dz0;
771
772 Amg::MatrixX W_mat(5*NTrk+6,5*NTrk+6); W_mat.setZero();
773 if (fullCov.size() != 0) {
774 W_mat.block(0,0,ndim,ndim) = fullCov;
775 } else {
776 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
777 W_mat.block(0,0,V0_cov.rows(),V0_cov.rows()) = V0_cov;
778 W_mat.block<3,3>(5*NTrk,5*NTrk) = vxCandidate->covariancePosition();
779 }
780 W_mat.block<3,3>(5*NTrk+3,5*NTrk+3) = vertex->covariancePosition();
781 V0_err = D_vec.transpose() * W_mat * D_vec;
782 } else if (ndim == 3*NTrk+3) {
783 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
784 D_vec(0) = da0dx;
785 D_vec(1) = da0dy;
786 D_vec(2) = da0dz;
787 for( unsigned int it=0; it<NTrk; it++) {
788 D_vec(3*it+3) = da0dphi[it];
789 D_vec(3*it+4) = da0dtheta[it];
790 D_vec(3*it+5) = da0dqOverP[it];
791 }
792 D_vec(3*NTrk+3) = da0dx0;
793 D_vec(3*NTrk+4) = da0dy0;
794 D_vec(3*NTrk+5) = da0dz0;
795
796 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
797 W_mat.block(0,0,ndim,ndim) = fullCov;
798 W_mat.block<3,3>(3*NTrk+3,3*NTrk+3) = vertex->covariancePosition();
799 V0_err = D_vec.transpose() * W_mat * D_vec;
800 }else{
801 ATH_MSG_ERROR("This should not happen - returning zero.");
802 return 0;
803 }
804
805 double a0Errsq = V0_err(0,0);
806 if (a0Errsq <= 0.) ATH_MSG_DEBUG("a0Error: negative sqrt a0Errsq " << a0Errsq);
807 double a0Err = (a0Errsq>0.) ? sqrt(a0Errsq) : 0.;
808 return a0Err;
809 }
810
811 double V0Tools::a0xyError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex) const
812 {
813 return a0Error(vxCandidate, vertex, false);
814 }
815
816 double V0Tools::a0Error(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, bool in3D) const
817 {
818 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
819 auto vert = vxCandidate->position() - vertex->position();
820 double dx = vert.x();
821 double dy = vert.y();
822 double dz = vert.z();
823 double Px=0., Py=0., Pz=0.;
824 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
825 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
826 std::vector<double>dpzdqOverP(NTrk), dpzdtheta(NTrk);
827 std::vector<double>da0dqOverP(NTrk), da0dtheta(NTrk), da0dphi(NTrk);
828
829 auto fullCov = convertCovMatrix(vxCandidate);
830 for( unsigned int it=0; it<NTrk; it++) {
831 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
832 double trkCharge = 1.;
833 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
834 double phi = bPer->parameters()[Trk::phi];
835 double theta = bPer->parameters()[Trk::theta];
836 double qOverP = bPer->parameters()[Trk::qOverP];
837 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
838 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
839 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
840 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
841 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
842 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
843 if ( in3D ) {
844 dpzdqOverP[it] = -(cos(theta)*trkCharge)/(qOverP*qOverP);
845 dpzdtheta[it] = -(sin(theta)*trkCharge)/qOverP;
846 }
847 Px += bPer->momentum()[Trk::px];
848 Py += bPer->momentum()[Trk::py];
849 Pz += bPer->momentum()[Trk::pz];
850 }
851 double cosineTheta;
852 double a0val;
853 if ( in3D ) {
854 cosineTheta = cosTheta(vxCandidate,vertex);
855 a0val = a0(vxCandidate,vertex);
856 } else { // transforms momentum and vertex separation to transverse plane and leads also to zero dXdz derivatives
857 cosineTheta = cosTheta_xy(vxCandidate,vertex);
858 a0val = a0xy(vxCandidate,vertex);
859 Pz = 0.;
860 dz = 0.;
861 }
862 double P = sqrt(Px*Px+Py*Py+Pz*Pz);
863 double r = sqrt(dx*dx+dy*dy+dz*dz);
864
865 double da0dx = (Px/P*r*cosineTheta - dx)/a0val;
866 double da0dy = (Py/P*r*cosineTheta - dy)/a0val;
867 double da0dz = (Pz/P*r*cosineTheta - dz)/a0val;
868 double da0dx0 = -da0dx;
869 double da0dy0 = -da0dy;
870 double da0dz0 = -da0dz;
871 for( unsigned int it=0; it<NTrk; it++) {
872 da0dqOverP[it] = -(r*cosineTheta/P)*(da0dx*dpxdqOverP[it]+da0dy*dpydqOverP[it]+da0dz*dpzdqOverP[it]);
873 da0dtheta[it] = -(r*cosineTheta/P)*(da0dx*dpxdtheta[it]+da0dy*dpydtheta[it]+da0dz*dpzdtheta[it]);
874 da0dphi[it] = -(r*cosineTheta/P)*(da0dx*dpxdphi[it]+da0dy*dpydphi[it]);
875 }
876
877 unsigned int ndim = 0;
878 if (fullCov.size() != 0) {
879 ndim = fullCov.rows();
880 } else {
881 ndim = 5*NTrk+3;
882 }
883
884 Amg::MatrixX V0_err;
885 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
886 Amg::MatrixX D_vec(5*NTrk+6,1); D_vec.setZero();
887 for( unsigned int it=0; it<NTrk; it++) {
888 D_vec(5*it+0) = 0.;
889 D_vec(5*it+1) = 0.;
890 D_vec(5*it+2) = da0dphi[it];
891 D_vec(5*it+3) = da0dtheta[it];
892 D_vec(5*it+4) = da0dqOverP[it];
893 }
894 D_vec(5*NTrk+0) = da0dx;
895 D_vec(5*NTrk+1) = da0dy;
896 D_vec(5*NTrk+2) = da0dz;
897 D_vec(5*NTrk+3) = da0dx0;
898 D_vec(5*NTrk+4) = da0dy0;
899 D_vec(5*NTrk+5) = da0dz0;
900
901 Amg::MatrixX W_mat(5*NTrk+6,5*NTrk+6); W_mat.setZero();
902 if (fullCov.size() != 0) {
903 W_mat.block(0,0,ndim,ndim) = fullCov;
904 } else {
905 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
906 W_mat.block(0,0,V0_cov.rows(),V0_cov.rows()) = V0_cov;
907 W_mat.block<3,3>(5*NTrk,5*NTrk) = vxCandidate->covariancePosition();
908 }
909 W_mat.block<3,3>(5*NTrk+3,5*NTrk+3) = vertex->covariancePosition();
910 V0_err = D_vec.transpose() * W_mat * D_vec;
911 } else if (ndim == 3*NTrk+3) {
912 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
913 D_vec(0) = da0dx;
914 D_vec(1) = da0dy;
915 D_vec(2) = da0dz;
916 for( unsigned int it=0; it<NTrk; it++) {
917 D_vec(3*it+3) = da0dphi[it];
918 D_vec(3*it+4) = da0dtheta[it];
919 D_vec(3*it+5) = da0dqOverP[it];
920 }
921 D_vec(3*NTrk+3) = da0dx0;
922 D_vec(3*NTrk+4) = da0dy0;
923 D_vec(3*NTrk+5) = da0dz0;
924
925 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
926 W_mat.block(0,0,ndim,ndim) = fullCov;
927 W_mat.block<3,3>(3*NTrk+3,3*NTrk+3) = vertex->covariancePosition();
928 V0_err = D_vec.transpose() * W_mat * D_vec;
929 }else{
930 ATH_MSG_ERROR("This should not happen - returning zero.");
931 return 0;
932 }
933
934 double a0Errsq = V0_err(0,0);
935 if (a0Errsq <= 0.) ATH_MSG_DEBUG("a0Error: negative sqrt a0Errsq " << a0Errsq);
936 double a0Err = (a0Errsq>0.) ? sqrt(a0Errsq) : 0.;
937 return a0Err;
938 }
939
940 double V0Tools::lxy(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex)
941 {
942 auto vert = vxCandidate->position() - vertex->position();
943 double dx = vert.x();
944 double dy = vert.y();
945 Amg::Vector3D mom = V0Momentum(vxCandidate);
946 double dxy = (mom.x()*dx + mom.y()*dy)/mom.perp();
947 return dxy;
948 }
949
950 double V0Tools::lxyError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex) const
951 {
952 auto vert = vxCandidate->position() - vertex->position();
953 double dx = vert.x();
954 double dy = vert.y();
955 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
956 double Px=0., Py=0.;
957 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
958 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
959 std::vector<double>dLxydqOverP(NTrk), dLxydtheta(NTrk), dLxydphi(NTrk);
960
961 auto fullCov = convertCovMatrix(vxCandidate);
962 for( unsigned int it=0; it<NTrk; it++) {
963 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
964 double trkCharge = 1.;
965 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
966 double phi = bPer->parameters()[Trk::phi];
967 double theta = bPer->parameters()[Trk::theta];
968 double qOverP = bPer->parameters()[Trk::qOverP];
969 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
970 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
971 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
972 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
973 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
974 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
975 Px += bPer->momentum()[Trk::px];
976 Py += bPer->momentum()[Trk::py];
977 }
978 double PTsq = Px*Px+Py*Py;
979 double PT = (PTsq>0.) ? sqrt(PTsq) : 0.;
980 if (PTsq == 0.)[[unlikely]]{
981 ATH_MSG_ERROR("lxyError: Divisor is zero - returning zero.");
982 return 0.;
983 }
984 double LXYoverPT = (Px*dx+Py*dy)/PTsq;
985
986 for( unsigned int it=0; it<NTrk; it++) {
987 double dPTdqOverP = (Px*dpxdqOverP[it]+Py*dpydqOverP[it])/PT;
988 double dPTdtheta = (Px*dpxdtheta[it]+Py*dpydtheta[it])/PT;
989 double dPTdphi = (Px*dpxdphi[it]+Py*dpydphi[it])/PT;
990 dLxydqOverP[it] = (dx*dpxdqOverP[it]+dy*dpydqOverP[it])/PT-LXYoverPT*dPTdqOverP;
991 dLxydtheta[it] = (dx*dpxdtheta[it]+dy*dpydtheta[it])/PT-LXYoverPT*dPTdtheta;
992 dLxydphi[it] = (dx*dpxdphi[it]+dy*dpydphi[it])/PT-LXYoverPT*dPTdphi;
993 }
994 double dLxydx = Px/PT;
995 double dLxydy = Py/PT;
996 double dLxydx0 = -dLxydx;
997 double dLxydy0 = -dLxydy;
998
999 unsigned int ndim = 0;
1000 if (fullCov.size() != 0) {
1001 ndim = fullCov.rows();
1002 } else {
1003 ndim = 5*NTrk+3;
1004 }
1005
1006 Amg::MatrixX V0_err;
1007 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
1008 Amg::MatrixX D_vec(5*NTrk+6,1); D_vec.setZero();
1009 for( unsigned int it=0; it<NTrk; it++) {
1010 D_vec(5*it+0) = 0.;
1011 D_vec(5*it+1) = 0.;
1012 D_vec(5*it+2) = dLxydphi[it];
1013 D_vec(5*it+3) = dLxydtheta[it];
1014 D_vec(5*it+4) = dLxydqOverP[it];
1015 }
1016 D_vec(5*NTrk+0) = dLxydx;
1017 D_vec(5*NTrk+1) = dLxydy;
1018 D_vec(5*NTrk+2) = 0.;
1019 D_vec(5*NTrk+3) = dLxydx0;
1020 D_vec(5*NTrk+4) = dLxydy0;
1021 D_vec(5*NTrk+5) = 0.;
1022
1023 Amg::MatrixX W_mat(5*NTrk+6,5*NTrk+6); W_mat.setZero();
1024 if (fullCov.size() != 0) {
1025 W_mat.block(0,0,ndim,ndim) = fullCov;
1026 } else {
1027 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
1028 W_mat.block(0,0,V0_cov.rows(),V0_cov.rows()) = V0_cov;
1029 W_mat.block<3,3>(5*NTrk,5*NTrk) = vxCandidate->covariancePosition();
1030 }
1031 W_mat.block<3,3>(5*NTrk+3,5*NTrk+3) = vertex->covariancePosition();
1032 V0_err = D_vec.transpose() * W_mat * D_vec;
1033 } else if (ndim == 3*NTrk+3) {
1034 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
1035 D_vec(0) = dLxydx;
1036 D_vec(1) = dLxydy;
1037 D_vec(2) = 0.;
1038 for( unsigned int it=0; it<NTrk; it++) {
1039 D_vec(3*it+3) = dLxydphi[it];
1040 D_vec(3*it+4) = dLxydtheta[it];
1041 D_vec(3*it+5) = dLxydqOverP[it];
1042 }
1043 D_vec(3*NTrk+3) = dLxydx0;
1044 D_vec(3*NTrk+4) = dLxydy0;
1045 D_vec(3*NTrk+5) = 0.;
1046
1047 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
1048 W_mat.block(0,0,ndim,ndim) = fullCov;
1049 W_mat.block<3,3>(3*NTrk+3,3*NTrk+3) = vertex->covariancePosition();
1050 V0_err = D_vec.transpose() * W_mat * D_vec;
1051 }else{
1052 ATH_MSG_ERROR("This should not happen - returning zero.");
1053 return 0;
1054 }
1055
1056 double LxyErrsq = V0_err(0,0);
1057 if (LxyErrsq <= 0.) ATH_MSG_DEBUG("lxyError: negative sqrt LxyErrsq " << LxyErrsq);
1058 return (LxyErrsq>0.) ? sqrt(LxyErrsq) : 0.;
1059 }
1060
1061 double V0Tools::lxyz(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex)
1062 {
1063 auto vert = vxCandidate->position() - vertex->position();
1064 double dx = vert.x();
1065 double dy = vert.y();
1066 double dz = vert.z();
1067 Amg::Vector3D mom = V0Momentum(vxCandidate);
1068 double dxyz= (mom.x()*dx + mom.y()*dy + mom.z()*dz)/mom.mag();
1069 return dxyz;
1070 }
1071
1072 double V0Tools::lxyzError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex) const
1073 {
1074 auto vert = vxCandidate->position() - vertex->position();
1075 double dx = vert.x();
1076 double dy = vert.y();
1077 double dz = vert.z();
1078 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1079 double Px=0., Py=0., Pz=0.;
1080 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
1081 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
1082 std::vector<double>dpzdqOverP(NTrk), dpzdtheta(NTrk);
1083 std::vector<double>dLxyzdqOverP(NTrk), dLxyzdtheta(NTrk), dLxyzdphi(NTrk);
1084
1085 auto fullCov = convertCovMatrix(vxCandidate);
1086 for( unsigned int it=0; it<NTrk; it++) {
1087 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
1088 double trkCharge = 1.;
1089 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
1090 double phi = bPer->parameters()[Trk::phi];
1091 double theta = bPer->parameters()[Trk::theta];
1092 double qOverP = bPer->parameters()[Trk::qOverP];
1093 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
1094 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
1095 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
1096 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
1097 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
1098 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
1099 dpzdqOverP[it] = -(cos(theta)*trkCharge)/(qOverP*qOverP);
1100 dpzdtheta[it] = -(sin(theta)*trkCharge)/qOverP;
1101 Px += bPer->momentum()[Trk::px];
1102 Py += bPer->momentum()[Trk::py];
1103 Pz += bPer->momentum()[Trk::pz];
1104 }
1105 double Psq = Px*Px+Py*Py+Pz*Pz;
1106 double P = (Psq>0.) ? sqrt(Psq) : 0.;
1107 if (Psq == 0.)[[unlikely]]{
1108 ATH_MSG_ERROR("lxyzError: Divisor is zero - returning zero.");
1109 return 0.;
1110 }
1111 double LXYZoverP = (Px*dx+Py*dy+Pz*dz)/Psq;
1112
1113 for( unsigned int it=0; it<NTrk; it++) {
1114 double dPdqOverP = (Px*dpxdqOverP[it]+Py*dpydqOverP[it]+Pz*dpzdqOverP[it])/P;
1115 double dPdtheta = (Px*dpxdtheta[it]+Py*dpydtheta[it]+Pz*dpzdtheta[it])/P;
1116 double dPdphi = (Px*dpxdphi[it]+Py*dpydphi[it])/P;
1117 dLxyzdqOverP[it] = (dx*dpxdqOverP[it]+dy*dpydqOverP[it]+dz*dpzdqOverP[it])/P-LXYZoverP*dPdqOverP;
1118 dLxyzdtheta[it] = (dx*dpxdtheta[it]+dy*dpydtheta[it]+dz*dpzdtheta[it])/P-LXYZoverP*dPdtheta;
1119 dLxyzdphi[it] = (dx*dpxdphi[it]+dy*dpydphi[it])/P-LXYZoverP*dPdphi;
1120 }
1121 double dLxyzdx = Px/P;
1122 double dLxyzdy = Py/P;
1123 double dLxyzdz = Pz/P;
1124 double dLxyzdx0 = -dLxyzdx;
1125 double dLxyzdy0 = -dLxyzdy;
1126 double dLxyzdz0 = -dLxyzdz;
1127
1128 unsigned int ndim = 0;
1129 if (fullCov.size() != 0) {
1130 ndim = fullCov.rows();
1131 } else {
1132 ndim = 5*NTrk+3;
1133 }
1134
1135 Amg::MatrixX V0_err;
1136 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
1137 Amg::MatrixX D_vec(5*NTrk+6,1); D_vec.setZero();
1138 for( unsigned int it=0; it<NTrk; it++) {
1139 D_vec(5*it+0) = 0.;
1140 D_vec(5*it+1) = 0.;
1141 D_vec(5*it+2) = dLxyzdphi[it];
1142 D_vec(5*it+3) = dLxyzdtheta[it];
1143 D_vec(5*it+4) = dLxyzdqOverP[it];
1144 }
1145 D_vec(5*NTrk+0) = dLxyzdx;
1146 D_vec(5*NTrk+1) = dLxyzdy;
1147 D_vec(5*NTrk+2) = dLxyzdz;
1148 D_vec(5*NTrk+3) = dLxyzdx0;
1149 D_vec(5*NTrk+4) = dLxyzdy0;
1150 D_vec(5*NTrk+5) = dLxyzdz0;
1151
1152 Amg::MatrixX W_mat(5*NTrk+6,5*NTrk+6); W_mat.setZero();
1153 if (fullCov.size() != 0) {
1154 W_mat.block(0,0,ndim,ndim) = fullCov;
1155 } else {
1156 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
1157 W_mat.block(0,0,V0_cov.rows(),V0_cov.rows()) = V0_cov;
1158 W_mat.block<3,3>(5*NTrk,5*NTrk) = vxCandidate->covariancePosition();
1159 }
1160 W_mat.block<3,3>(5*NTrk+3,5*NTrk+3) = vertex->covariancePosition();
1161 V0_err = D_vec.transpose() * W_mat * D_vec;
1162 } else if (ndim == 3*NTrk+3) {
1163 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
1164 D_vec(0) = dLxyzdx;
1165 D_vec(1) = dLxyzdy;
1166 D_vec(2) = dLxyzdz;
1167 for( unsigned int it=0; it<NTrk; it++) {
1168 D_vec(3*it+3) = dLxyzdphi[it];
1169 D_vec(3*it+4) = dLxyzdtheta[it];
1170 D_vec(3*it+5) = dLxyzdqOverP[it];
1171 }
1172 D_vec(3*NTrk+3) = dLxyzdx0;
1173 D_vec(3*NTrk+4) = dLxyzdy0;
1174 D_vec(3*NTrk+5) = dLxyzdz0;
1175
1176 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
1177 W_mat.block(0,0,ndim,ndim) = fullCov;
1178 W_mat.block<3,3>(3*NTrk+3,3*NTrk+3) = vertex->covariancePosition();
1179 V0_err = D_vec.transpose() * W_mat * D_vec;
1180 }else{
1181 ATH_MSG_ERROR("This should not happen - returning zero.");
1182 return 0;
1183 }
1184
1185 double LxyzErrsq = V0_err(0,0);
1186 if (LxyzErrsq <= 0.) ATH_MSG_DEBUG("lxyzError: negative sqrt LxyzErrsq " << LxyzErrsq);
1187 return (LxyzErrsq>0.) ? sqrt(LxyzErrsq) : 0.;
1188 }
1189
1190 double V0Tools::tau(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, double posTrackMass, double negTrackMass) const
1191 {
1192 std::array<double, 2> masses = {posTrackMass, negTrackMass};
1193
1194 return tau(vxCandidate,vertex,masses);
1195 }
1196
1197 double V0Tools::tau(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, std::span<const double> masses) const
1198 {
1199 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1200 if (masses.size() != NTrk) {
1201 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
1202 return -999999.;
1203 }
1204 //double CONST = 1000./CLHEP::c_light;
1205 double CONST = 1000./299.792;
1206 double M = invariantMass(vxCandidate, masses);
1207 double LXY = lxy(vxCandidate,vertex);
1208 double PT = V0Momentum(vxCandidate).perp();
1209 return CONST*M*LXY/PT;
1210 }
1211
1212 double V0Tools::tau(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, double posTrackMass, double negTrackMass, double massV0) const
1213 {
1214 std::array<double, 2> masses = {posTrackMass, negTrackMass};
1215
1216 return tau(vxCandidate,vertex,masses,massV0);
1217 }
1218
1219 double V0Tools::tau(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, std::span<const double> masses, double massV0) const
1220 {
1221 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1222 if (masses.size() != NTrk) {
1223 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
1224 return -999999.;
1225 }
1226 Amg::MatrixX cov = tauMassCovariance(vxCandidate,vertex,masses);
1227 double Tau = tau(vxCandidate,vertex,masses);
1228 double mass = invariantMass(vxCandidate,masses);
1229 double descr = cov(0,0)*cov(1,1)-cov(0,1)*cov(0,1);
1230 double cov_i11 = cov(1,1)/descr;
1231 double cov_i12 = -cov(0,1)/descr;
1232 double deltaTau = -(massV0-mass)*cov_i12/cov_i11;
1233 return Tau + deltaTau;
1234 }
1235
1236 double V0Tools::tau(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, double M)
1237 {
1238 //double CONST = 1000./CLHEP::c_light;
1239 double CONST = 1000./299.792;
1240 double LXY = lxy(vxCandidate,vertex);
1241 double PT = V0Momentum(vxCandidate).perp();
1242 return CONST*M*LXY/PT;
1243 }
1244
1245 double V0Tools::tauError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, double posTrackMass, double negTrackMass) const
1246 {
1247 // Tau = CONST*M*(Px*dx+Py*dy)/(PT*PT)
1248 std::array<double, 2> masses = {posTrackMass, negTrackMass};
1249
1250 return tauError(vxCandidate,vertex,masses);
1251 }
1252
1253 double V0Tools::tauError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, std::span<const double> masses) const
1254 {
1255 // Tau = CONST*M*(Px*dx+Py*dy)/(PT*PT)
1256 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1257 if (masses.size() != NTrk) {
1258 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
1259 return -999999.;
1260 }
1261 //double CONST = 1000./CLHEP::c_light;
1262 double CONST = 1000./299.792;
1263 double PT = V0Momentum(vxCandidate).perp();
1264 auto vert = vxCandidate->position() - vertex->position();
1265 double dx = vert.x();
1266 double dy = vert.y();
1267 double M = invariantMass(vxCandidate, masses);
1268 double E=0., Px=0., Py=0., Pz=0.;
1269 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
1270 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
1271 std::vector<double>dpzdqOverP(NTrk), dpzdtheta(NTrk), dedqOverP(NTrk);
1272 std::vector<double>dTaudqOverP(NTrk), dTaudtheta(NTrk), dTaudphi(NTrk);
1273
1274 auto fullCov = convertCovMatrix(vxCandidate);
1275 for( unsigned int it=0; it<NTrk; it++) {
1276 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
1277 double trkCharge = 1.;
1278 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
1279 double phi = bPer->parameters()[Trk::phi];
1280 double theta = bPer->parameters()[Trk::theta];
1281 double qOverP = bPer->parameters()[Trk::qOverP];
1282 double tmp = 1./(qOverP*qOverP) + masses[it]*masses[it];
1283 double pe = (tmp>0.) ? sqrt(tmp) : 0.;
1284 dedqOverP[it] = -1./(qOverP*qOverP*qOverP*pe);
1285 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
1286 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
1287 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
1288 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
1289 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
1290 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
1291 dpzdqOverP[it] = -(cos(theta)*trkCharge)/(qOverP*qOverP);
1292 dpzdtheta[it] = -(sin(theta)*trkCharge)/qOverP;
1293 E += pe;
1294 Px += bPer->momentum()[Trk::px];
1295 Py += bPer->momentum()[Trk::py];
1296 Pz += bPer->momentum()[Trk::pz];
1297 }
1298 double LXY = Px*dx+Py*dy;
1299
1300 for( unsigned int it=0; it<NTrk; it++) {
1301 double dMdqOverP = -(Px*dpxdqOverP[it]+Py*dpydqOverP[it]+Pz*dpzdqOverP[it]-E*dedqOverP[it])/M;
1302 double dMdtheta = -(Px*dpxdtheta[it]+Py*dpydtheta[it]+Pz*dpzdtheta[it])/M;
1303 double dMdphi = -(Px*dpxdphi[it]+Py*dpydphi[it])/M;
1304 double dPTdqOverP = (Px*dpxdqOverP[it]+Py*dpydqOverP[it])/PT;
1305 double dPTdtheta = (Px*dpxdtheta[it]+Py*dpydtheta[it])/PT;
1306 double dPTdphi = (Px*dpxdphi[it]+Py*dpydphi[it])/PT;
1307 double dLXYdqOverP = dx*dpxdqOverP[it]+dy*dpydqOverP[it];
1308 double dLXYdtheta = dx*dpxdtheta[it]+dy*dpydtheta[it];
1309 double dLXYdphi = dx*dpxdphi[it]+dy*dpydphi[it];
1310 dTaudqOverP[it] = (LXY*dMdqOverP+M*dLXYdqOverP)/(PT*PT)-(2.*LXY*M*dPTdqOverP)/(PT*PT*PT);
1311 dTaudtheta[it] = (LXY*dMdtheta+M*dLXYdtheta)/(PT*PT)-(2.*LXY*M*dPTdtheta)/(PT*PT*PT);
1312 dTaudphi[it] = (LXY*dMdphi+M*dLXYdphi)/(PT*PT)-(2.*LXY*M*dPTdphi)/(PT*PT*PT);
1313 }
1314 double dTaudx = (M*Px)/(PT*PT);
1315 double dTaudy = (M*Py)/(PT*PT);
1316 double dTaudx0 = -dTaudx;
1317 double dTaudy0 = -dTaudy;
1318
1319 unsigned int ndim = 0;
1320 if (fullCov.size() != 0) {
1321 ndim = fullCov.rows();
1322 } else {
1323 ndim = 5*NTrk+3;
1324 }
1325
1326 Amg::MatrixX V0_err;
1327 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
1328 Amg::MatrixX D_vec(5*NTrk+6,1); D_vec.setZero();
1329 for( unsigned int it=0; it<NTrk; it++) {
1330 D_vec(5*it+0) = 0.;
1331 D_vec(5*it+1) = 0.;
1332 D_vec(5*it+2) = dTaudphi[it];
1333 D_vec(5*it+3) = dTaudtheta[it];
1334 D_vec(5*it+4) = dTaudqOverP[it];
1335 }
1336 D_vec(5*NTrk+0) = dTaudx;
1337 D_vec(5*NTrk+1) = dTaudy;
1338 D_vec(5*NTrk+2) = 0.;
1339 D_vec(5*NTrk+3) = dTaudx0;
1340 D_vec(5*NTrk+4) = dTaudy0;
1341 D_vec(5*NTrk+5) = 0.;
1342
1343 Amg::MatrixX W_mat(5*NTrk+6,5*NTrk+6); W_mat.setZero();
1344 if (fullCov.size() != 0) {
1345 W_mat.block(0,0,ndim,ndim) = fullCov;
1346 } else {
1347 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
1348 W_mat.block(0,0,V0_cov.rows(),V0_cov.rows()) = V0_cov;
1349 W_mat.block<3,3>(5*NTrk,5*NTrk) = vxCandidate->covariancePosition();
1350 }
1351 W_mat.block<3,3>(5*NTrk+3,5*NTrk+3) = vertex->covariancePosition();
1352 V0_err = D_vec.transpose() * W_mat * D_vec;
1353 } else if (ndim == 3*NTrk+3) {
1354 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
1355 D_vec(0) = dTaudx;
1356 D_vec(1) = dTaudy;
1357 D_vec(2) = 0.;
1358 for( unsigned int it=0; it<NTrk; it++) {
1359 D_vec(3*it+3) = dTaudphi[it];
1360 D_vec(3*it+4) = dTaudtheta[it];
1361 D_vec(3*it+5) = dTaudqOverP[it];
1362 }
1363 D_vec(3*NTrk+3) = dTaudx0;
1364 D_vec(3*NTrk+4) = dTaudy0;
1365 D_vec(3*NTrk+5) = 0.;
1366
1367 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
1368 W_mat.block(0,0,ndim,ndim) = fullCov;
1369 W_mat.block<3,3>(3*NTrk+3,3*NTrk+3) = vertex->covariancePosition();
1370 V0_err = D_vec.transpose() * W_mat * D_vec;
1371 }else{
1372 ATH_MSG_ERROR("This should not happen - returning zero.");
1373 return 0;
1374 }
1375
1376 double tauErrsq = V0_err(0,0);
1377 if (tauErrsq <= 0.) ATH_MSG_DEBUG("tauError: negative sqrt tauErrsq " << tauErrsq);
1378 double tauErr = (tauErrsq>0.) ? sqrt(tauErrsq) : 0.;
1379 return CONST*tauErr;
1380 }
1381
1382 double V0Tools::tauError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, double posTrackMass, double negTrackMass, double massV0) const
1383 {
1384 std::array<double, 2> masses = {posTrackMass, negTrackMass};
1385
1386 return tauError(vxCandidate,vertex,masses,massV0);
1387 }
1388
1389 double V0Tools::tauError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, std::span<const double> masses, double ) const
1390 {
1391 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1392 if (masses.size() != NTrk) {
1393 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
1394 return -999999.;
1395 }
1396 double error = -999999.;
1397 Amg::MatrixX cov = tauMassCovariance(vxCandidate,vertex,masses);
1398 double descr = cov(0,0)*cov(1,1)-cov(0,1)*cov(0,1);
1399 double cov_i11 = cov(1,1)/descr;
1400 if (cov_i11 > 0.) error = 1./sqrt(cov_i11);
1401 return error;
1402 }
1403
1404 double V0Tools::tauError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, double M) const
1405 {
1406 // Tau = CONST*M*(Px*dx+Py*dy)/(PT*PT)
1407 //double CONST = 1000./CLHEP::c_light;
1408 double CONST = 1000./299.792;
1409 double PT = V0Momentum(vxCandidate).perp();
1410 auto vecsub = vxCandidate->position() - vertex->position();
1411 double dx = vecsub.x();
1412 double dy = vecsub.y();
1413 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1414 double Px=0., Py=0.;
1415 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
1416 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
1417 std::vector<double>dPTdtheta(NTrk), dPTdphi(NTrk);
1418 std::vector<double>dTaudqOverP(NTrk), dTaudtheta(NTrk), dTaudphi(NTrk);
1419
1420 auto fullCov = convertCovMatrix(vxCandidate);
1421 for( unsigned int it=0; it<NTrk; it++) {
1422 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
1423 double trkCharge = 1.;
1424 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
1425 double phi = bPer->parameters()[Trk::phi];
1426 double theta = bPer->parameters()[Trk::theta];
1427 double qOverP = bPer->parameters()[Trk::qOverP];
1428 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
1429 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
1430 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
1431 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
1432 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
1433 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
1434 Px += bPer->momentum()[Trk::px];
1435 Py += bPer->momentum()[Trk::py];
1436 }
1437 double LXY = Px*dx+Py*dy;
1438
1439 for( unsigned int it=0; it<NTrk; it++) {
1440 double dPTdqOverP = (Px*dpxdqOverP[it]+Py*dpydqOverP[it])/PT;
1441 double dPTdtheta = (Px*dpxdtheta[it]+Py*dpydtheta[it])/PT;
1442 double dPTdphi = (Px*dpxdphi[it]+Py*dpydphi[it])/PT;
1443 double dLXYdqOverP = dx*dpxdqOverP[it]+dy*dpydqOverP[it];
1444 double dLXYdtheta = dx*dpxdtheta[it]+dy*dpydtheta[it];
1445 double dLXYdphi = dx*dpxdphi[it]+dy*dpydphi[it];
1446 dTaudqOverP[it] = M*dLXYdqOverP/(PT*PT)-(2.*LXY*M*dPTdqOverP)/(PT*PT*PT);
1447 dTaudtheta[it] = M*dLXYdtheta/(PT*PT)-(2.*LXY*M*dPTdtheta)/(PT*PT*PT);
1448 dTaudphi[it] = M*dLXYdphi/(PT*PT)-(2.*LXY*M*dPTdphi)/(PT*PT*PT);
1449 }
1450 double dTaudx = (M*Px)/(PT*PT);
1451 double dTaudy = (M*Py)/(PT*PT);
1452 double dTaudx0 = -dTaudx;
1453 double dTaudy0 = -dTaudy;
1454
1455 unsigned int ndim = 0;
1456 if (fullCov.size() != 0) {
1457 ndim = fullCov.rows();
1458 } else {
1459 ndim = 5*NTrk+3;
1460 }
1461
1462 Amg::MatrixX V0_err;
1463 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
1464 Amg::MatrixX D_vec(5*NTrk+6,1); D_vec.setZero();
1465 for( unsigned int it=0; it<NTrk; it++) {
1466 D_vec(5*it+0) = 0.;
1467 D_vec(5*it+1) = 0.;
1468 D_vec(5*it+2) = dTaudphi[it];
1469 D_vec(5*it+3) = dTaudtheta[it];
1470 D_vec(5*it+4) = dTaudqOverP[it];
1471 }
1472 D_vec(5*NTrk+0) = dTaudx;
1473 D_vec(5*NTrk+1) = dTaudy;
1474 D_vec(5*NTrk+2) = 0.;
1475 D_vec(5*NTrk+3) = dTaudx0;
1476 D_vec(5*NTrk+4) = dTaudy0;
1477 D_vec(5*NTrk+5) = 0.;
1478
1479 Amg::MatrixX W_mat(5*NTrk+6,5*NTrk+6); W_mat.setZero();
1480 if (fullCov.size() != 0) {
1481 W_mat.block(0,0,ndim,ndim) = fullCov;
1482 } else {
1483 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
1484 W_mat.block(0,0,V0_cov.rows(),V0_cov.rows()) = V0_cov;
1485 W_mat.block<3,3>(5*NTrk,5*NTrk) = vxCandidate->covariancePosition();
1486 }
1487 W_mat.block<3,3>(5*NTrk+3,5*NTrk+3) = vertex->covariancePosition();
1488 V0_err = D_vec.transpose() * W_mat * D_vec;
1489 } else if (ndim == 3*NTrk+3) {
1490 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
1491 D_vec(0) = dTaudx;
1492 D_vec(1) = dTaudy;
1493 D_vec(2) = 0.;
1494 for( unsigned int it=0; it<NTrk; it++) {
1495 D_vec(3*it+3) = dTaudphi[it];
1496 D_vec(3*it+4) = dTaudtheta[it];
1497 D_vec(3*it+5) = dTaudqOverP[it];
1498 }
1499 D_vec(3*NTrk+3) = dTaudx0;
1500 D_vec(3*NTrk+4) = dTaudy0;
1501 D_vec(3*NTrk+5) = 0.;
1502
1503 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
1504 W_mat.block(0,0,ndim,ndim) = fullCov;
1505 W_mat.block<3,3>(3*NTrk+3,3*NTrk+3) = vertex->covariancePosition();
1506 V0_err = D_vec.transpose() * W_mat * D_vec;
1507 }else{
1508 ATH_MSG_ERROR("This should not happen - returning zero.");
1509 return 0;
1510 }
1511
1512 double tauErrsq = V0_err(0,0);
1513 if (tauErrsq <= 0.) ATH_MSG_DEBUG("tauError: negative sqrt tauErrsq " << tauErrsq);
1514 double tauErr = (tauErrsq>0.) ? sqrt(tauErrsq) : 0.;
1515 return CONST*tauErr;
1516 }
1517
1518 double V0Tools::tau3D(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, std::span<const double> masses) const
1519 {
1520 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1521 if (masses.size() != NTrk) {
1522 ATH_MSG_DEBUG("The provided number of masses does not match the number of tracks in the vertex");
1523 return -999999.;
1524 }
1525 //double CONST = 1000./CLHEP::c_light;
1526 double CONST = 1000./299.792;
1527 double M = invariantMass(vxCandidate, masses);
1528 double LXYZ = lxyz(vxCandidate,vertex);
1529 double P = V0Momentum(vxCandidate).mag();
1530 return CONST*M*LXYZ/P;
1531 }
1532
1533 double V0Tools::tau3D(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, double M)
1534 {
1535 //double CONST = 1000./CLHEP::c_light;
1536 double CONST = 1000./299.792;
1537 double LXYZ = lxyz(vxCandidate,vertex);
1538 double P = V0Momentum(vxCandidate).mag();
1539 return CONST*M*LXYZ/P;
1540 }
1541
1542 double V0Tools::tau3DError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, std::span<const double> masses) const
1543 {
1544 // Tau = CONST*M*(Px*dx+Py*dy+Pz*dz)/(P*P)
1545 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1546 if (masses.size() != NTrk) {
1547 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
1548 return -999999.;
1549 }
1550 //double CONST = 1000./CLHEP::c_light;
1551 double CONST = 1000./299.792;
1552 double P = V0Momentum(vxCandidate).mag();
1553 auto vert = vxCandidate->position() - vertex->position();
1554 double dx = vert.x();
1555 double dy = vert.y();
1556 double dz = vert.z();
1557 double M = invariantMass(vxCandidate, masses);
1558 double E=0., Px=0., Py=0., Pz=0.;
1559 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
1560 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
1561 std::vector<double>dpzdqOverP(NTrk), dpzdtheta(NTrk), dedqOverP(NTrk);
1562 std::vector<double>dTaudqOverP(NTrk), dTaudtheta(NTrk), dTaudphi(NTrk);
1563
1564 auto fullCov = convertCovMatrix(vxCandidate);
1565 for( unsigned int it=0; it<NTrk; it++) {
1566 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
1567 double trkCharge = 1.;
1568 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
1569 double phi = bPer->parameters()[Trk::phi];
1570 double theta = bPer->parameters()[Trk::theta];
1571 double qOverP = bPer->parameters()[Trk::qOverP];
1572 double tmp = 1./(qOverP*qOverP) + masses[it]*masses[it];
1573 double pe = (tmp>0.) ? sqrt(tmp) : 0.;
1574 dedqOverP[it] = -1./(qOverP*qOverP*qOverP*pe);
1575 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
1576 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
1577 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
1578 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
1579 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
1580 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
1581 dpzdqOverP[it] = -(cos(theta)*trkCharge)/(qOverP*qOverP);
1582 dpzdtheta[it] = -(sin(theta)*trkCharge)/qOverP;
1583 E += pe;
1584 Px += bPer->momentum()[Trk::px];
1585 Py += bPer->momentum()[Trk::py];
1586 Pz += bPer->momentum()[Trk::pz];
1587 }
1588 double LXYZ = Px*dx+Py*dy+Pz*dz;
1589
1590 for( unsigned int it=0; it<NTrk; it++) {
1591 double dMdqOverP = -(Px*dpxdqOverP[it]+Py*dpydqOverP[it]+Pz*dpzdqOverP[it]-E*dedqOverP[it])/M;
1592 double dMdtheta = -(Px*dpxdtheta[it]+Py*dpydtheta[it]+Pz*dpzdtheta[it])/M;
1593 double dMdphi = -(Px*dpxdphi[it]+Py*dpydphi[it])/M;
1594 double dPdqOverP = (Px*dpxdqOverP[it]+Py*dpydqOverP[it]+Pz*dpzdqOverP[it])/P;
1595 double dPdtheta = (Px*dpxdtheta[it]+Py*dpydtheta[it]+Pz*dpzdtheta[it])/P;
1596 double dPdphi = (Px*dpxdphi[it]+Py*dpydphi[it])/P;
1597 double dLXYZdqOverP = dx*dpxdqOverP[it]+dy*dpydqOverP[it]+dz*dpzdqOverP[it];
1598 double dLXYZdtheta = dx*dpxdtheta[it]+dy*dpydtheta[it]+dz*dpzdtheta[it];
1599 double dLXYZdphi = dx*dpxdphi[it]+dy*dpydphi[it];
1600 dTaudqOverP[it] = (LXYZ*dMdqOverP+M*dLXYZdqOverP)/(P*P)-(2.*LXYZ*M*dPdqOverP)/(P*P*P);
1601 dTaudtheta[it] = (LXYZ*dMdtheta+M*dLXYZdtheta)/(P*P)-(2.*LXYZ*M*dPdtheta)/(P*P*P);
1602 dTaudphi[it] = (LXYZ*dMdphi+M*dLXYZdphi)/(P*P)-(2.*LXYZ*M*dPdphi)/(P*P*P);
1603 }
1604 double dTaudx = (M*Px)/(P*P);
1605 double dTaudy = (M*Py)/(P*P);
1606 double dTaudz = (M*Pz)/(P*P);
1607 double dTaudx0 = -dTaudx;
1608 double dTaudy0 = -dTaudy;
1609 double dTaudz0 = -dTaudz;
1610
1611 unsigned int ndim = 0;
1612 if (fullCov.size() != 0) {
1613 ndim = fullCov.rows();
1614 } else {
1615 ndim = 5*NTrk+3;
1616 }
1617
1618 Amg::MatrixX V0_err;
1619 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
1620 Amg::MatrixX D_vec(5*NTrk+6,1); D_vec.setZero();
1621 for( unsigned int it=0; it<NTrk; it++) {
1622 D_vec(5*it+0) = 0.;
1623 D_vec(5*it+1) = 0.;
1624 D_vec(5*it+2) = dTaudphi[it];
1625 D_vec(5*it+3) = dTaudtheta[it];
1626 D_vec(5*it+4) = dTaudqOverP[it];
1627 }
1628 D_vec(5*NTrk+0) = dTaudx;
1629 D_vec(5*NTrk+1) = dTaudy;
1630 D_vec(5*NTrk+2) = dTaudz;
1631 D_vec(5*NTrk+3) = dTaudx0;
1632 D_vec(5*NTrk+4) = dTaudy0;
1633 D_vec(5*NTrk+5) = dTaudz0;
1634
1635 Amg::MatrixX W_mat(5*NTrk+6,5*NTrk+6); W_mat.setZero();
1636 if (fullCov.size() != 0) {
1637 W_mat.block(0,0,ndim,ndim) = fullCov;
1638 } else {
1639 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
1640 W_mat.block(0,0,V0_cov.rows(),V0_cov.rows()) = V0_cov;
1641 W_mat.block<3,3>(5*NTrk,5*NTrk) = vxCandidate->covariancePosition();
1642 }
1643 W_mat.block<3,3>(5*NTrk+3,5*NTrk+3) = vertex->covariancePosition();
1644 V0_err = D_vec.transpose() * W_mat * D_vec;
1645 } else if (ndim == 3*NTrk+3) {
1646 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
1647 D_vec(0) = dTaudx;
1648 D_vec(1) = dTaudy;
1649 D_vec(2) = dTaudz;
1650 for( unsigned int it=0; it<NTrk; it++) {
1651 D_vec(3*it+3) = dTaudphi[it];
1652 D_vec(3*it+4) = dTaudtheta[it];
1653 D_vec(3*it+5) = dTaudqOverP[it];
1654 }
1655 D_vec(3*NTrk+3) = dTaudx0;
1656 D_vec(3*NTrk+4) = dTaudy0;
1657 D_vec(3*NTrk+5) = dTaudz0;
1658
1659 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
1660 W_mat.block(0,0,ndim,ndim) = fullCov;
1661 W_mat.block<3,3>(3*NTrk+3,3*NTrk+3) = vertex->covariancePosition();
1662 V0_err = D_vec.transpose() * W_mat * D_vec;
1663 }else{
1664 ATH_MSG_ERROR("This should not happen - returning zero.");
1665 return 0;
1666 }
1667
1668 double tauErrsq = V0_err(0,0);
1669 if (tauErrsq <= 0.) ATH_MSG_DEBUG("tauError: negative sqrt tauErrsq " << tauErrsq);
1670 double tauErr = (tauErrsq>0.) ? sqrt(tauErrsq) : 0.;
1671 return CONST*tauErr;
1672 }
1673
1674 double V0Tools::tau3DError(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, double M) const
1675 {
1676 // Tau = CONST*M*(Px*dx+Py*dy+Pz*dz)/(P*P)
1677 //double CONST = 1000./CLHEP::c_light;
1678 double CONST = 1000./299.792;
1679 double P = V0Momentum(vxCandidate).mag();
1680 auto vecsub = vxCandidate->position() - vertex->position();
1681 double dx = vecsub.x();
1682 double dy = vecsub.y();
1683 double dz = vecsub.z();
1684 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1685 double Px=0., Py=0., Pz=0.;
1686 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
1687 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
1688 std::vector<double>dpzdqOverP(NTrk), dpzdtheta(NTrk);
1689 std::vector<double>dTaudqOverP(NTrk), dTaudtheta(NTrk), dTaudphi(NTrk);
1690
1691 auto fullCov = convertCovMatrix(vxCandidate);
1692 for( unsigned int it=0; it<NTrk; it++) {
1693 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
1694 double trkCharge = 1.;
1695 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
1696 double phi = bPer->parameters()[Trk::phi];
1697 double theta = bPer->parameters()[Trk::theta];
1698 double qOverP = bPer->parameters()[Trk::qOverP];
1699 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
1700 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
1701 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
1702 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
1703 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
1704 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
1705 dpzdqOverP[it] = -(cos(theta)*trkCharge)/(qOverP*qOverP);
1706 dpzdtheta[it] = -(sin(theta)*trkCharge)/qOverP;
1707 Px += bPer->momentum()[Trk::px];
1708 Py += bPer->momentum()[Trk::py];
1709 Pz += bPer->momentum()[Trk::pz];
1710 }
1711 double LXYZ = Px*dx+Py*dy+Pz*dz;
1712
1713 for( unsigned int it=0; it<NTrk; it++) {
1714 double dPdqOverP = (Px*dpxdqOverP[it]+Py*dpydqOverP[it]+Pz*dpzdqOverP[it])/P;
1715 double dPdtheta = (Px*dpxdtheta[it]+Py*dpydtheta[it]+Pz*dpzdtheta[it])/P;
1716 double dPdphi = (Px*dpxdphi[it]+Py*dpydphi[it])/P;
1717 double dLXYZdqOverP = dx*dpxdqOverP[it]+dy*dpydqOverP[it]+dz*dpzdqOverP[it];
1718 double dLXYZdtheta = dx*dpxdtheta[it]+dy*dpydtheta[it]+dz*dpzdtheta[it];
1719 double dLXYZdphi = dx*dpxdphi[it]+dy*dpydphi[it];
1720 dTaudqOverP[it] = M*dLXYZdqOverP/(P*P)-(2.*LXYZ*M*dPdqOverP)/(P*P*P);
1721 dTaudtheta[it] = M*dLXYZdtheta/(P*P)-(2.*LXYZ*M*dPdtheta)/(P*P*P);
1722 dTaudphi[it] = M*dLXYZdphi/(P*P)-(2.*LXYZ*M*dPdphi)/(P*P*P);
1723 }
1724 double dTaudx = (M*Px)/(P*P);
1725 double dTaudy = (M*Py)/(P*P);
1726 double dTaudz = (M*Pz)/(P*P);
1727 double dTaudx0 = -dTaudx;
1728 double dTaudy0 = -dTaudy;
1729 double dTaudz0 = -dTaudz;
1730
1731 unsigned int ndim = 0;
1732 if (fullCov.size() != 0) {
1733 ndim = fullCov.rows();
1734 } else {
1735 ndim = 5*NTrk+3;
1736 }
1737
1738 Amg::MatrixX V0_err;
1739 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
1740 Amg::MatrixX D_vec(5*NTrk+6,1); D_vec.setZero();
1741 for( unsigned int it=0; it<NTrk; it++) {
1742 D_vec(5*it+0) = 0.;
1743 D_vec(5*it+1) = 0.;
1744 D_vec(5*it+2) = dTaudphi[it];
1745 D_vec(5*it+3) = dTaudtheta[it];
1746 D_vec(5*it+4) = dTaudqOverP[it];
1747 }
1748 D_vec(5*NTrk+0) = dTaudx;
1749 D_vec(5*NTrk+1) = dTaudy;
1750 D_vec(5*NTrk+2) = dTaudz;
1751 D_vec(5*NTrk+3) = dTaudx0;
1752 D_vec(5*NTrk+4) = dTaudy0;
1753 D_vec(5*NTrk+5) = dTaudz0;
1754
1755 Amg::MatrixX W_mat(5*NTrk+6,5*NTrk+6); W_mat.setZero();
1756 if (fullCov.size() != 0) {
1757 W_mat.block(0,0,ndim,ndim) = fullCov;
1758 } else {
1759 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
1760 W_mat.block(0,0,V0_cov.rows(),V0_cov.rows()) = V0_cov;
1761 W_mat.block<3,3>(5*NTrk,5*NTrk) = vxCandidate->covariancePosition();
1762 }
1763 W_mat.block<3,3>(5*NTrk+3,5*NTrk+3) = vertex->covariancePosition();
1764 V0_err = D_vec.transpose() * W_mat * D_vec;
1765 } else if (ndim == 3*NTrk+3) {
1766 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
1767 D_vec(0) = dTaudx;
1768 D_vec(1) = dTaudy;
1769 D_vec(2) = dTaudz;
1770 for( unsigned int it=0; it<NTrk; it++) {
1771 D_vec(3*it+3) = dTaudphi[it];
1772 D_vec(3*it+4) = dTaudtheta[it];
1773 D_vec(3*it+5) = dTaudqOverP[it];
1774 }
1775 D_vec(3*NTrk+3) = dTaudx0;
1776 D_vec(3*NTrk+4) = dTaudy0;
1777 D_vec(3*NTrk+5) = dTaudz0;
1778
1779 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
1780 W_mat.block(0,0,ndim,ndim) = fullCov;
1781 W_mat.block<3,3>(3*NTrk+3,3*NTrk+3) = vertex->covariancePosition();
1782 V0_err = D_vec.transpose() * W_mat * D_vec;
1783 }else{
1784 ATH_MSG_ERROR("This should not happen - returning zero.");
1785 return 0;
1786 }
1787
1788 double tauErrsq = V0_err(0,0);
1789 if (tauErrsq <= 0.) ATH_MSG_DEBUG("tauError: negative sqrt tauErrsq " << tauErrsq);
1790 double tauErr = (tauErrsq>0.) ? sqrt(tauErrsq) : 0.;
1791 return CONST*tauErr;
1792 }
1793
1794 double V0Tools::thetaStar(const xAOD::Vertex * vxCandidate, double mass1, double mass2)
1795 {
1796 double theta = 0.;
1797 xAOD::TrackParticle::FourMom_t V1 = positiveTrack4Momentum(vxCandidate, mass1);
1798 xAOD::TrackParticle::FourMom_t V2 = negativeTrack4Momentum(vxCandidate, mass2);
1799 CLHEP::HepLorentzVector v1(V1.Px(),V1.Py(),V1.Pz(),V1.E());
1800 CLHEP::HepLorentzVector v2(V2.Px(),V2.Py(),V2.Pz(),V2.E());
1801 CLHEP::HepLorentzVector v0 = v1 + v2;
1802 CLHEP::Hep3Vector boost_v0 = v0.boostVector();
1803 boost_v0 *= -1.0;
1804 v1.boost( boost_v0 );
1805 v2.boost( boost_v0 );
1806 theta = v0.angle( v1.vect() );
1807 return theta;
1808 }
1809
1810 double V0Tools::cosThetaStar(const xAOD::Vertex * vxCandidate, double posTrackMass, double negTrackMass)
1811 {
1812 xAOD::TrackParticle::FourMom_t PosMom = positiveTrack4Momentum(vxCandidate, posTrackMass);
1813 xAOD::TrackParticle::FourMom_t NegMom = negativeTrack4Momentum(vxCandidate, negTrackMass);
1814 CLHEP::HepLorentzVector posMom(PosMom.Px(),PosMom.Py(),PosMom.Pz(),PosMom.E());
1815 CLHEP::HepLorentzVector negMom(NegMom.Px(),NegMom.Py(),NegMom.Pz(),NegMom.E());
1816 return cosThetaStar(posMom, negMom);
1817 }
1818
1819 double V0Tools::cosThetaStar(const CLHEP::HepLorentzVector & posTrack, const CLHEP::HepLorentzVector & negTrack)
1820 {
1821 CLHEP::HepLorentzVector v0(posTrack + negTrack);
1822 double Mv0 = v0.m();
1823 double Mplus = posTrack.m();
1824 double Mminus= negTrack.m();
1825 double Pv0 = v0.rho();
1826 double pssq = (Mv0*Mv0-(Mplus+Mminus)*(Mplus+Mminus))*(Mv0*Mv0-(Mplus-Mminus)*(Mplus-Mminus));
1827 double ps = (pssq>0.) ? sqrt(pssq) : 0.;
1828 ps /= 2.*Mv0;
1829 double pp = v0.px()*posTrack.px() + v0.py()*posTrack.py() + v0.pz()*posTrack.pz();
1830 return ( v0.e()*pp - Pv0*Pv0*posTrack.e())/( Mv0*ps*Pv0);
1831 }
1832
1833 //double V0Tools::phiStar(xAOD::Vertex * vxCandidate) const
1834 double V0Tools::phiStar(const xAOD::Vertex * vxCandidate, double posTrackMass, double negTrackMass)
1835 {
1836 xAOD::TrackParticle::FourMom_t V_pos = positiveTrack4Momentum(vxCandidate,posTrackMass);
1837 xAOD::TrackParticle::FourMom_t V_neg = negativeTrack4Momentum(vxCandidate,negTrackMass);
1838 CLHEP::HepLorentzVector v_pos(V_pos.Px(),V_pos.Py(),V_pos.Pz(),V_pos.E());
1839 CLHEP::HepLorentzVector v_neg(V_neg.Px(),V_neg.Py(),V_neg.Pz(),V_neg.E());
1840 return phiStar(v_pos+v_neg,v_pos);
1841 }
1842
1843 double V0Tools::phiStar(const CLHEP::HepLorentzVector & v0, const CLHEP::HepLorentzVector & track)
1844 {
1845 double phiStar = -999999.;
1846 CLHEP::Hep3Vector V0 = v0.getV();
1847 CLHEP::Hep3Vector trk = track.getV();
1848 //double v0_rapidity = v0.rapidity();
1849 trk.rotateZ(-V0.phi());
1850 trk.rotateY(-V0.theta());
1851 //if (v0_rapidity < 0.) {
1852 // trk.rotateZ(-M_PI);
1853 // phiStar = -atan2(trk.y(),trk.x());
1854 //} else {
1855 // phiStar = atan2(trk.y(),trk.x());
1856 //}
1857 phiStar = atan2(trk.y(),trk.x());
1858 return phiStar;
1859 }
1860
1861 double V0Tools::cosTheta(const xAOD::Vertex * vxCandidate, const Amg::Vector3D& vertex)
1862 {
1863 auto mom = V0Momentum(vxCandidate);
1864 auto vtx1 = vtx(vxCandidate);
1865 vtx1 -= vertex;
1866 return (mom.dot(vtx1))/(mom.mag()*(vtx1).mag());
1867 }
1868
1869 double V0Tools::cosTheta(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex)
1870 {
1871 auto mom = V0Momentum(vxCandidate);
1872 auto vtx1 = vtx(vxCandidate);
1873 vtx1 -= vertex->position();
1874 return (mom.dot((vtx1)))/(mom.mag()*(vtx1).mag());
1875 }
1876
1877 double V0Tools::cosTheta_xy(const xAOD::Vertex * vxCandidate, const Amg::Vector3D& vertex)
1878 {
1879 auto mom = V0Momentum(vxCandidate);
1880 auto vtx1 = vtx(vxCandidate);
1881 vtx1 -= vertex;
1882 double pT = mom.perp();
1883 return (mom.x()*vtx1.x()+mom.y()*vtx1.y())/(pT*vtx1.perp());
1884 }
1885
1886 double V0Tools::cosTheta_xy(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex)
1887 {
1888 auto mom = V0Momentum(vxCandidate);
1889 auto vtx1 = vtx(vxCandidate);
1890 vtx1 -= vertex->position();
1891 double pT = mom.perp();
1892 return (mom.x()*vtx1.x()+mom.y()*vtx1.y())/(pT*vtx1.perp());
1893 }
1894
1895 float V0Tools::charge(const xAOD::Vertex * vxCandidate)
1896 {
1897 float ch = 0.;
1898 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1899 for( unsigned int it=0; it<NTrk; it++) {
1900 float trkCharge = vxCandidate->trackParticle(it)->charge();
1901 //const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
1902 //float trkCharge = 1.;
1903 //if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
1904 ch += trkCharge;
1905 }
1906 return ch;
1907 }
1908
1909 const xAOD::TrackParticle* V0Tools::origTrack(const xAOD::Vertex * vxCandidate, int trkIndex)
1910 {
1911 return vxCandidate->trackParticle(trkIndex);
1912 }
1913
1914 const xAOD::TrackParticle* V0Tools::positiveOrigTrack(const xAOD::Vertex * vxCandidate)
1915 {
1916 const xAOD::TrackParticle* origTrk(nullptr);
1917 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1918 if (NTrk != 2) return origTrk;
1919 for( unsigned int it=0; it<NTrk; it++) {
1920 if(vxCandidate->trackParticle(it)->charge() > 0) origTrk = vxCandidate->trackParticle(it);
1921 }
1922 return origTrk;
1923 }
1924
1925 const xAOD::TrackParticle* V0Tools::negativeOrigTrack(const xAOD::Vertex * vxCandidate)
1926 {
1927 const xAOD::TrackParticle* origTrk(nullptr);
1928 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1929 if (NTrk != 2) return origTrk;
1930 for( unsigned int it=0; it<NTrk; it++) {
1931 if(vxCandidate->trackParticle(it)->charge() < 0) origTrk = vxCandidate->trackParticle(it);
1932 }
1933 return origTrk;
1934 }
1935
1936 double V0Tools::invariantMassBeforeFitIP(const xAOD::Vertex * vxCandidate, std::span<const double> masses) const
1937 {
1938 double px = 0., py = 0., pz = 0., e = 0.;
1939 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1940 if (masses.size() != NTrk) {
1941 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
1942 return -999999.;
1943 }
1944 for( unsigned int it=0; it<NTrk; it++) {
1945 if (masses[it] >= 0.) {
1946 const xAOD::TrackParticle* TP = origTrack(vxCandidate,it);
1947 if (TP == nullptr) return -999999.;
1948 TLorentzVector Tp4 = TP->p4();
1949 px += Tp4.Px();
1950 py += Tp4.Py();
1951 pz += Tp4.Pz();
1952 double pesq = 1./(TP->qOverP()*TP->qOverP()) + masses[it]*masses[it];
1953 double pe = (pesq>0.) ? sqrt(pesq) : 0.;
1954 e += pe;
1955 }
1956 }
1957 double msq = e*e - px*px - py*py - pz*pz;
1958 return (msq>0.) ? sqrt(msq) : 0.;
1959 }
1960
1961 double V0Tools::invariantMassBeforeFit(const xAOD::Vertex * vxCandidate, std::span<const double> masses, const EventContext& ctx, const Trk::IExtrapolator* extrapolator) const
1962 {
1963 Trk::PerigeeSurface perigeeSurface(vxCandidate->position());
1964 double px = 0., py = 0., pz = 0., e = 0.;
1965 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1966 if (masses.size() != NTrk) {
1967 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
1968 return -999999.;
1969 }
1970 for( unsigned int it=0; it<NTrk; it++) {
1971 if (masses[it] >= 0.) {
1972 const xAOD::TrackParticle* TP = origTrack(vxCandidate,it);
1973 if (TP == nullptr) return -999999.;
1974 std::unique_ptr<const Trk::TrackParameters> extrPer =
1975 extrapolator->extrapolate(ctx, TP->perigeeParameters(), perigeeSurface);
1976 if (extrPer == nullptr)
1977 return -999999.;
1978 px += extrPer->momentum().x();
1979 py += extrPer->momentum().y();
1980 pz += extrPer->momentum().z();
1981 double pesq = extrPer->momentum().mag2() + masses[it]*masses[it];
1982 double pe = (pesq>0.) ? sqrt(pesq) : 0.;
1983 e += pe;
1984 }
1985 }
1986 double msq = e*e - px*px - py*py - pz*pz;
1987 return (msq>0.) ? sqrt(msq) : 0.;
1988 }
1989
1990 double V0Tools::invariantMassBeforeFit(const xAOD::Vertex * vxCandidate,
1991 std::span<const double> masses, const Amg::Vector3D& vertex, const EventContext& ctx, const Trk::IExtrapolator* extrap) const
1992 {
1993 Trk::PerigeeSurface perigeeSurface(vertex);
1994 double px = 0., py = 0., pz = 0., e = 0.;
1995 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
1996 if (masses.size() != NTrk) {
1997 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
1998 return -999999.;
1999 }
2000 for( unsigned int it=0; it<NTrk; it++) {
2001 if (masses[it] >= 0.) {
2002 const xAOD::TrackParticle* TP = origTrack(vxCandidate,it);
2003 if (TP == nullptr) return -999999.;
2004 std::unique_ptr<const Trk::TrackParameters> extrPer =
2005 extrap->extrapolate(ctx, TP->perigeeParameters(), perigeeSurface);
2006 if (extrPer == nullptr)
2007 return -999999.;
2008 px += extrPer->momentum().x();
2009 py += extrPer->momentum().y();
2010 pz += extrPer->momentum().z();
2011 double pesq = extrPer->momentum().mag2() + masses[it]*masses[it];
2012 double pe = (pesq>0.) ? sqrt(pesq) : 0.;
2013 e += pe;
2014 }
2015 }
2016 double msq = e*e - px*px - py*py - pz*pz;
2017 return (msq>0.) ? sqrt(msq) : 0.;
2018 }
2019
2020 double V0Tools::invariantMassErrorBeforeFitIP(const xAOD::Vertex * vxCandidate, std::span<const double> masses) const
2021 {
2022 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
2023 if (masses.size() != NTrk) {
2024 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
2025 return -999999.;
2026 }
2027 double mass = invariantMassBeforeFitIP(vxCandidate, masses);
2028 double E=0., Px=0., Py=0., Pz=0.;
2029 std::vector<double>phi(NTrk), theta(NTrk), qOverP(NTrk), charge(NTrk), e(NTrk);
2030 std::vector<double>dm2dphi(NTrk), dm2dtheta(NTrk), dm2dqOverP(NTrk);
2031 Amg::MatrixX V0_cor(5*NTrk,5*NTrk); V0_cor.setZero();
2032 for( unsigned int it=0; it<NTrk; it++) {
2033 if (masses[it] >= 0.) {
2034 const xAOD::TrackParticle* TP = origTrack(vxCandidate,it);
2035 if (TP == nullptr) return -999999.;
2036 const xAOD::ParametersCovMatrix_t cov_tmp = TP->definingParametersCovMatrix();
2037 V0_cor(5*it+2,5*it+2) = cov_tmp(2,2);
2038 V0_cor(5*it+2,5*it+3) = cov_tmp(2,3);
2039 V0_cor(5*it+2,5*it+4) = cov_tmp(2,4);
2040 V0_cor(5*it+3,5*it+3) = cov_tmp(3,3);
2041 V0_cor(5*it+3,5*it+4) = cov_tmp(3,4);
2042 V0_cor(5*it+4,5*it+4) = cov_tmp(4,4);
2043 V0_cor(5*it+3,5*it+2) = cov_tmp(2,3);
2044 V0_cor(5*it+4,5*it+2) = cov_tmp(2,4);
2045 V0_cor(5*it+4,5*it+3) = cov_tmp(3,4);
2046 charge[it] = TP->charge();
2047 phi[it] = TP->phi();
2048 theta[it] = TP->theta();
2049 qOverP[it] = TP->qOverP();
2050 double tmp = 1./(qOverP[it]*qOverP[it]) + masses[it]*masses[it];
2051 double pe = (tmp>0.) ? sqrt(tmp) : 0.;
2052 e[it] = pe;
2053 E += e[it];
2054 TLorentzVector p4 = TP->p4();
2055 Px += p4.Px();
2056 Py += p4.Py();
2057 Pz += p4.Pz();
2058 }
2059 }
2060
2061 for( unsigned int it=0; it<NTrk; it++) {
2062 if (masses[it] >= 0.) {
2063 dm2dphi[it] = 2.*(Px*sin(phi[it])-Py*cos(phi[it]))*sin(theta[it])*charge[it]/qOverP[it];
2064 dm2dtheta[it] = 2.*(Pz*sin(theta[it])-(Px*cos(phi[it])+Py*sin(phi[it]))*cos(theta[it]))*charge[it]/qOverP[it];
2065 dm2dqOverP[it] = 2.*(Pz*cos(theta[it])+(Px*cos(phi[it])+Py*sin(phi[it]))*sin(theta[it])-E*charge[it]/(qOverP[it]*e[it]))*charge[it]/(qOverP[it]*qOverP[it]);
2066 }
2067 }
2068
2069 Amg::MatrixX D_vec(5*NTrk,1); D_vec.setZero();
2070 for( unsigned int it=0; it<NTrk; it++) {
2071 D_vec(5*it+0,0) = 0.;
2072 D_vec(5*it+1,0) = 0.;
2073 D_vec(5*it+2,0) = dm2dphi[it];
2074 D_vec(5*it+3,0) = dm2dtheta[it];
2075 D_vec(5*it+4,0) = dm2dqOverP[it];
2076 }
2077 Amg::MatrixX V0_merr = D_vec.transpose() * V0_cor * D_vec;
2078
2079 double massVarsq = V0_merr(0,0);
2080 if (massVarsq <= 0.) ATH_MSG_DEBUG("massError: negative sqrt massVarsq " << massVarsq);
2081 double massVar = (massVarsq>0.) ? sqrt(massVarsq) : 0.;
2082 double massErr = massVar/(2.*mass);
2083 return massErr;
2084 }
2085
2086
2087 double V0Tools::invariantMassErrorBeforeFit(const xAOD::Vertex * vxCandidate, std::span<const double> masses, const EventContext& ctx, const Trk::IExtrapolator* extrap) const
2088 {
2089 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
2090 if (masses.size() != NTrk) {
2091 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
2092 return -999999.;
2093 }
2094 Amg::Vector3D vertex = vxCandidate->position();
2095 return invariantMassErrorBeforeFit(vxCandidate,masses,vertex, ctx, extrap);
2096 }
2097
2098
2099 double V0Tools::invariantMassErrorBeforeFit(const xAOD::Vertex * vxCandidate,
2100 std::span<const double> masses, const Amg::Vector3D& vertex, const EventContext& ctx, const Trk::IExtrapolator* extrap) const
2101 {
2102 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
2103 if (masses.size() != NTrk) {
2104 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
2105 return -999999.;
2106 }
2107 Trk::PerigeeSurface perigeeSurface(vertex);
2108 double E=0., Px=0., Py=0., Pz=0.;
2109 std::vector<double>phi(NTrk), theta(NTrk), qOverP(NTrk), charge(NTrk), e(NTrk);
2110 std::vector<double>dm2dphi(NTrk), dm2dtheta(NTrk), dm2dqOverP(NTrk);
2111 Amg::MatrixX V0_cor(5*NTrk,5*NTrk); V0_cor.setZero();
2112 for( unsigned int it=0; it<NTrk; it++) {
2113 if (masses[it] >= 0.) {
2114 const xAOD::TrackParticle* TP = origTrack(vxCandidate,it);
2115 if (TP == nullptr) return -999999.;
2116 std::unique_ptr<const Trk::TrackParameters> extrPer =
2117 extrap->extrapolate(ctx, TP->perigeeParameters(), perigeeSurface);
2118 if (extrPer == nullptr)
2119 return -999999.;
2120 const AmgSymMatrix(5)* cov_tmp = extrPer->covariance();
2121 V0_cor(5*it+2,5*it+2) = (*cov_tmp)(2,2);
2122 V0_cor(5*it+2,5*it+3) = (*cov_tmp)(2,3);
2123 V0_cor(5*it+2,5*it+4) = (*cov_tmp)(2,4);
2124 V0_cor(5*it+3,5*it+3) = (*cov_tmp)(3,3);
2125 V0_cor(5*it+3,5*it+4) = (*cov_tmp)(3,4);
2126 V0_cor(5*it+4,5*it+4) = (*cov_tmp)(4,4);
2127 V0_cor(5*it+3,5*it+2) = (*cov_tmp)(2,3);
2128 V0_cor(5*it+4,5*it+2) = (*cov_tmp)(2,4);
2129 V0_cor(5*it+4,5*it+3) = (*cov_tmp)(3,4);
2130 charge[it] = TP->charge();
2131 phi[it] = extrPer->parameters()[Trk::phi];
2132 theta[it] = extrPer->parameters()[Trk::theta];
2133 qOverP[it] = extrPer->parameters()[Trk::qOverP];
2134 double tmp = 1./(qOverP[it]*qOverP[it]) + masses[it]*masses[it];
2135 double pe = (tmp>0.) ? sqrt(tmp) : 0.;
2136 e[it] = pe;
2137 E += e[it];
2138 Px += extrPer->momentum().x();
2139 Py += extrPer->momentum().y();
2140 Pz += extrPer->momentum().z();
2141 }
2142 }
2143 double msq = E*E - Px*Px - Py*Py - Pz*Pz;
2144 double mass = (msq>0.) ? sqrt(msq) : 0.;
2145
2146 for( unsigned int it=0; it<NTrk; it++) {
2147 if (masses[it] >= 0.) {
2148 dm2dphi[it] = 2.*(Px*sin(phi[it])-Py*cos(phi[it]))*sin(theta[it])*charge[it]/qOverP[it];
2149 dm2dtheta[it] = 2.*(Pz*sin(theta[it])-(Px*cos(phi[it])+Py*sin(phi[it]))*cos(theta[it]))*charge[it]/qOverP[it];
2150 dm2dqOverP[it] = 2.*(Pz*cos(theta[it])+(Px*cos(phi[it])+Py*sin(phi[it]))*sin(theta[it])-E*charge[it]/(qOverP[it]*e[it]))*charge[it]/(qOverP[it]*qOverP[it]);
2151 }
2152 }
2153
2154 Amg::MatrixX D_vec(5*NTrk,1); D_vec.setZero();
2155 for( unsigned int it=0; it<NTrk; it++) {
2156 D_vec(5*it+0,0) = 0.;
2157 D_vec(5*it+1,0) = 0.;
2158 D_vec(5*it+2,0) = dm2dphi[it];
2159 D_vec(5*it+3,0) = dm2dtheta[it];
2160 D_vec(5*it+4,0) = dm2dqOverP[it];
2161 }
2162 Amg::MatrixX V0_merr = D_vec.transpose() * V0_cor * D_vec;
2163
2164 double massVarsq = V0_merr(0,0);
2165 if (massVarsq <= 0.) ATH_MSG_DEBUG("massError: negative sqrt massVarsq " << massVarsq);
2166 double massVar = (massVarsq>0.) ? sqrt(massVarsq) : 0.;
2167 double massErr = massVar/(2.*mass);
2168 return massErr;
2169 }
2170
2171
2172 double V0Tools::massTauCov(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, std::span<const double> masses) const
2173 {
2174 // Tau = CONST*M*(Px*dx+Py*dy)/(PT*PT)
2175 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
2176 if (masses.size() != NTrk) {
2177 ATH_MSG_DEBUG("The provided number of masses does not match the number of tracks in the vertex");
2178 return -999999.;
2179 }
2180 //double CONST = 1000./CLHEP::c_light;
2181 double CONST = 1000./299.792;
2182 double PT = V0Momentum(vxCandidate).perp();
2183 auto vert = vxCandidate->position() - vertex->position();
2184 double dx = vert.x();
2185 double dy = vert.y();
2186 double M = invariantMass(vxCandidate, masses);
2187 double E=0., Px=0., Py=0., Pz=0.;
2188 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
2189 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
2190 std::vector<double>dpzdqOverP(NTrk), dpzdtheta(NTrk), dedqOverP(NTrk);
2191 std::vector<double>dMdqOverP(NTrk), dMdtheta(NTrk), dMdphi(NTrk);
2192 std::vector<double>dTaudqOverP(NTrk), dTaudtheta(NTrk), dTaudphi(NTrk);
2193
2194 auto fullCov = convertCovMatrix(vxCandidate);
2195 for( unsigned int it=0; it<NTrk; it++) {
2196 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
2197 double trkCharge = 1.;
2198 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
2199 double phi = bPer->parameters()[Trk::phi];
2200 double theta = bPer->parameters()[Trk::theta];
2201 double qOverP = bPer->parameters()[Trk::qOverP];
2202 double tmp = 1./(qOverP*qOverP) + masses[it]*masses[it];
2203 double pe = (tmp>0.) ? sqrt(tmp) : 0.;
2204 dedqOverP[it] = -1./(qOverP*qOverP*qOverP*pe);
2205 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
2206 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
2207 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
2208 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
2209 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
2210 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
2211 dpzdqOverP[it] = -(cos(theta)*trkCharge)/(qOverP*qOverP);
2212 dpzdtheta[it] = -(sin(theta)*trkCharge)/qOverP;
2213 E += pe;
2214 Px += bPer->momentum()[Trk::px];
2215 Py += bPer->momentum()[Trk::py];
2216 Pz += bPer->momentum()[Trk::pz];
2217 }
2218 double LXY = Px*dx+Py*dy;
2219
2220 for( unsigned int it=0; it<NTrk; it++) {
2221 dMdqOverP[it] = -(Px*dpxdqOverP[it]+Py*dpydqOverP[it]+Pz*dpzdqOverP[it]-E*dedqOverP[it])/M;
2222 dMdtheta[it] = -(Px*dpxdtheta[it]+Py*dpydtheta[it]+Pz*dpzdtheta[it])/M;
2223 dMdphi[it] = -(Px*dpxdphi[it]+Py*dpydphi[it])/M;
2224 double dPTdqOverP = (Px*dpxdqOverP[it]+Py*dpydqOverP[it])/PT;
2225 double dPTdtheta = (Px*dpxdtheta[it]+Py*dpydtheta[it])/PT;
2226 double dPTdphi = (Px*dpxdphi[it]+Py*dpydphi[it])/PT;
2227 double dLXYdqOverP = dx*dpxdqOverP[it]+dy*dpydqOverP[it];
2228 double dLXYdtheta = dx*dpxdtheta[it]+dy*dpydtheta[it];
2229 double dLXYdphi = dx*dpxdphi[it]+dy*dpydphi[it];
2230 dTaudqOverP[it] = (LXY*dMdqOverP[it]+M*dLXYdqOverP)/(PT*PT)-(2.*LXY*M*dPTdqOverP)/(PT*PT*PT);
2231 dTaudtheta[it] = (LXY*dMdtheta[it]+M*dLXYdtheta)/(PT*PT)-(2.*LXY*M*dPTdtheta)/(PT*PT*PT);
2232 dTaudphi[it] = (LXY*dMdphi[it]+M*dLXYdphi)/(PT*PT)-(2.*LXY*M*dPTdphi)/(PT*PT*PT);
2233 }
2234 double dTaudx = (M*Px)/(PT*PT);
2235 double dTaudy = (M*Py)/(PT*PT);
2236 double dTaudx0 = -dTaudx;
2237 double dTaudy0 = -dTaudy;
2238
2239 unsigned int ndim = 0;
2240 if (fullCov.size() != 0) {
2241 ndim = fullCov.rows();
2242 } else {
2243 ndim = 5*NTrk+3;
2244 }
2245
2246 Amg::MatrixX V0_err;
2247 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
2248 Amg::MatrixX D_mat(5*NTrk+6,2); D_mat.setZero();
2249 for( unsigned int it=0; it<NTrk; it++) {
2250 D_mat(5*it+0,0) = 0.;
2251 D_mat(5*it+1,0) = 0.;
2252 D_mat(5*it+2,0) = CONST*dTaudphi[it];
2253 D_mat(5*it+3,0) = CONST*dTaudtheta[it];
2254 D_mat(5*it+4,0) = CONST*dTaudqOverP[it];
2255 D_mat(5*it+0,1) = 0.;
2256 D_mat(5*it+1,1) = 0.;
2257 D_mat(5*it+2,1) = dMdphi[it];
2258 D_mat(5*it+3,1) = dMdtheta[it];
2259 D_mat(5*it+4,1) = dMdqOverP[it];
2260 }
2261 D_mat(5*NTrk+0,0) = CONST*dTaudx;
2262 D_mat(5*NTrk+1,0) = CONST*dTaudy;
2263 D_mat(5*NTrk+2,0) = 0.;
2264 D_mat(5*NTrk+3,0) = CONST*dTaudx0;
2265 D_mat(5*NTrk+4,0) = CONST*dTaudy0;
2266 D_mat(5*NTrk+5,0) = 0.;
2267 Amg::MatrixX W_mat(5*NTrk+6,5*NTrk+6); W_mat.setZero();
2268 if (fullCov.size() != 0) {
2269 W_mat.block(0,0,ndim,ndim) = fullCov;
2270 } else {
2271 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
2272 W_mat.block(0,0,V0_cov.rows(),V0_cov.rows()) = V0_cov;
2273 W_mat.block<3,3>(5*NTrk,5*NTrk) = vxCandidate->covariancePosition();
2274 }
2275 W_mat.block<3,3>(5*NTrk+3,5*NTrk+3) = vertex->covariancePosition();
2276 V0_err = D_mat.transpose() * W_mat * D_mat;
2277 } else if (ndim == 3*NTrk+3) {
2278 Amg::MatrixX D_mat(3*NTrk+6,2); D_mat.setZero();
2279 D_mat(0,0) = CONST*dTaudx;
2280 D_mat(1,0) = CONST*dTaudy;
2281 D_mat(2,0) = 0.;
2282 for( unsigned int it=0; it<NTrk; it++) {
2283 D_mat(3*it+3,0) = CONST*dTaudphi[it];
2284 D_mat(3*it+4,0) = CONST*dTaudtheta[it];
2285 D_mat(3*it+5,0) = CONST*dTaudqOverP[it];
2286 D_mat(3*it+3,1) = dMdphi[it];
2287 D_mat(3*it+4,1) = dMdtheta[it];
2288 D_mat(3*it+5,1) = dMdqOverP[it];
2289 }
2290 D_mat(3*NTrk+3,0) = CONST*dTaudx0;
2291 D_mat(3*NTrk+4,0) = CONST*dTaudy0;
2292 D_mat(3*NTrk+5,0) = 0.;
2293 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
2294 W_mat.block(0,0,ndim,ndim) = fullCov;
2295 W_mat.block<3,3>(3*NTrk+3,3*NTrk+3) = vertex->covariancePosition();
2296 V0_err = D_mat.transpose() * W_mat * D_mat;
2297 }else{
2298 ATH_MSG_ERROR("This should not happen - returning zero.");
2299 return 0;
2300 }
2301 return V0_err(0,1);
2302 }
2303
2304 Amg::MatrixX V0Tools::makeV0Cov(const xAOD::Vertex * vxCandidate) {
2305 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
2306 Amg::MatrixX V0_cov(5*NTrk,5*NTrk); V0_cov.setZero();
2307 for( unsigned int it=0; it<NTrk; it++){
2308 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
2309 const AmgSymMatrix(5)* cov_tmp = bPer->covariance();
2310 V0_cov(5*it+2,5*it+2) = (*cov_tmp)(2,2);
2311 V0_cov(5*it+2,5*it+3) = (*cov_tmp)(2,3);
2312 V0_cov(5*it+2,5*it+4) = (*cov_tmp)(2,4);
2313 V0_cov(5*it+3,5*it+3) = (*cov_tmp)(3,3);
2314 V0_cov(5*it+3,5*it+4) = (*cov_tmp)(3,4);
2315 V0_cov(5*it+4,5*it+4) = (*cov_tmp)(4,4);
2316 V0_cov(5*it+3,5*it+2) = (*cov_tmp)(2,3);
2317 V0_cov(5*it+4,5*it+2) = (*cov_tmp)(2,4);
2318 V0_cov(5*it+4,5*it+3) = (*cov_tmp)(3,4);
2319 }
2320 return V0_cov;
2321 }
2322
2323 Amg::MatrixX V0Tools::tauMassCovariance(const xAOD::Vertex * vxCandidate, const xAOD::Vertex* vertex, std::span<const double> masses) const
2324 {
2325 // Tau = CONST*M*(Px*dx+Py*dy)/(PT*PT)
2326 Amg::MatrixX V0_err;
2327 unsigned int NTrk = vxCandidate->vxTrackAtVertex().size();
2328 if (masses.size() != NTrk) {
2329 ATH_MSG_ERROR("The provided number of masses does not match the number of tracks in the vertex");
2330 return V0_err;
2331 }
2332 //double CONST = 1000./CLHEP::c_light;
2333 double CONST = 1000./299.792;
2334 double PT = V0Momentum(vxCandidate).perp();
2335 auto vert = vxCandidate->position() - vertex->position();
2336 double dx = vert.x();
2337 double dy = vert.y();
2338 double M = invariantMass(vxCandidate, masses);
2339 double E=0., Px=0., Py=0., Pz=0.;
2340 std::vector<double>dpxdqOverP(NTrk), dpxdtheta(NTrk), dpxdphi(NTrk);
2341 std::vector<double>dpydqOverP(NTrk), dpydtheta(NTrk), dpydphi(NTrk);
2342 std::vector<double>dpzdqOverP(NTrk), dpzdtheta(NTrk), dedqOverP(NTrk);
2343 std::vector<double>dMdqOverP(NTrk), dMdtheta(NTrk), dMdphi(NTrk);
2344 std::vector<double>dTaudqOverP(NTrk), dTaudtheta(NTrk), dTaudphi(NTrk);
2345
2346 auto fullCov = convertCovMatrix(vxCandidate);
2347 for( unsigned int it=0; it<NTrk; it++) {
2348 const Trk::TrackParameters* bPer = vxCandidate->vxTrackAtVertex()[it].perigeeAtVertex();
2349 double trkCharge = 1.;
2350 if (bPer->parameters()[Trk::qOverP] < 0.) trkCharge = -1.;
2351 double phi = bPer->parameters()[Trk::phi];
2352 double theta = bPer->parameters()[Trk::theta];
2353 double qOverP = bPer->parameters()[Trk::qOverP];
2354 double tmp = 1./(qOverP*qOverP) + masses[it]*masses[it];
2355 double pe = (tmp>0.) ? sqrt(tmp) : 0.;
2356 dedqOverP[it] = -1./(qOverP*qOverP*qOverP*pe);
2357 dpxdqOverP[it] = -(sin(theta)*cos(phi)*trkCharge)/(qOverP*qOverP);
2358 dpxdtheta[it] = (cos(theta)*cos(phi)*trkCharge)/qOverP;
2359 dpxdphi[it] = -(sin(theta)*sin(phi)*trkCharge)/qOverP;
2360 dpydqOverP[it] = -(sin(theta)*sin(phi)*trkCharge)/(qOverP*qOverP);
2361 dpydtheta[it] = (cos(theta)*sin(phi)*trkCharge)/qOverP;
2362 dpydphi[it] = (sin(theta)*cos(phi)*trkCharge)/qOverP;
2363 dpzdqOverP[it] = -(cos(theta)*trkCharge)/(qOverP*qOverP);
2364 dpzdtheta[it] = -(sin(theta)*trkCharge)/qOverP;
2365 E += pe;
2366 Px += bPer->momentum()[Trk::px];
2367 Py += bPer->momentum()[Trk::py];
2368 Pz += bPer->momentum()[Trk::pz];
2369 }
2370 double LXY = Px*dx+Py*dy;
2371
2372 for( unsigned int it=0; it<NTrk; it++) {
2373 dMdqOverP[it] = -(Px*dpxdqOverP[it]+Py*dpydqOverP[it]+Pz*dpzdqOverP[it]-E*dedqOverP[it])/M;
2374 dMdtheta[it] = -(Px*dpxdtheta[it]+Py*dpydtheta[it]+Pz*dpzdtheta[it])/M;
2375 dMdphi[it] = -(Px*dpxdphi[it]+Py*dpydphi[it])/M;
2376 double dPTdqOverP = (Px*dpxdqOverP[it]+Py*dpydqOverP[it])/PT;
2377 double dPTdtheta = (Px*dpxdtheta[it]+Py*dpydtheta[it])/PT;
2378 double dPTdphi = (Px*dpxdphi[it]+Py*dpydphi[it])/PT;
2379 double dLXYdqOverP = dx*dpxdqOverP[it]+dy*dpydqOverP[it];
2380 double dLXYdtheta = dx*dpxdtheta[it]+dy*dpydtheta[it];
2381 double dLXYdphi = dx*dpxdphi[it]+dy*dpydphi[it];
2382 dTaudqOverP[it] = (LXY*dMdqOverP[it]+M*dLXYdqOverP)/(PT*PT)-(2.*LXY*M*dPTdqOverP)/(PT*PT*PT);
2383 dTaudtheta[it] = (LXY*dMdtheta[it]+M*dLXYdtheta)/(PT*PT)-(2.*LXY*M*dPTdtheta)/(PT*PT*PT);
2384 dTaudphi[it] = (LXY*dMdphi[it]+M*dLXYdphi)/(PT*PT)-(2.*LXY*M*dPTdphi)/(PT*PT*PT);
2385 }
2386 double dTaudx = (M*Px)/(PT*PT);
2387 double dTaudy = (M*Py)/(PT*PT);
2388 double dTaudx0 = -dTaudx;
2389 double dTaudy0 = -dTaudy;
2390
2391 unsigned int ndim = 0;
2392 if (fullCov.size() != 0) {
2393 ndim = fullCov.rows();
2394 } else {
2395 ndim = 5*NTrk+3;
2396 }
2397
2398 if (ndim == 5*NTrk+3 || ndim == 5*NTrk+6) {
2399 Amg::MatrixX D_mat(5*NTrk+6,2); D_mat.setZero();
2400 for( unsigned int it=0; it<NTrk; it++) {
2401 D_mat(5*it+0,0) = 0.;
2402 D_mat(5*it+1,0) = 0.;
2403 D_mat(5*it+2,0) = CONST*dTaudphi[it];
2404 D_mat(5*it+3,0) = CONST*dTaudtheta[it];
2405 D_mat(5*it+4,0) = CONST*dTaudqOverP[it];
2406 D_mat(5*it+0,1) = 0.;
2407 D_mat(5*it+1,1) = 0.;
2408 D_mat(5*it+2,1) = dMdphi[it];
2409 D_mat(5*it+3,1) = dMdtheta[it];
2410 D_mat(5*it+4,1) = dMdqOverP[it];
2411 }
2412 D_mat(5*NTrk+0,0) = CONST*dTaudx;
2413 D_mat(5*NTrk+1,0) = CONST*dTaudy;
2414 D_mat(5*NTrk+2,0) = 0.;
2415 D_mat(5*NTrk+3,0) = CONST*dTaudx0;
2416 D_mat(5*NTrk+4,0) = CONST*dTaudy0;
2417 D_mat(5*NTrk+5,0) = 0.;
2418 Amg::MatrixX W_mat(5*NTrk+6,5*NTrk+6); W_mat.setZero();
2419 if (fullCov.size() != 0) {
2420 W_mat.block(0,0,ndim,ndim) = fullCov;
2421 } else {
2422 Amg::MatrixX V0_cov = makeV0Cov(vxCandidate);
2423 W_mat.block(0,0,V0_cov.rows(),V0_cov.rows()) = V0_cov;
2424 W_mat.block<3,3>(5*NTrk,5*NTrk) = vxCandidate->covariancePosition();
2425 }
2426 W_mat.block<3,3>(5*NTrk+3,5*NTrk) = vertex->covariancePosition();
2427 V0_err = D_mat.transpose() * W_mat * D_mat;
2428 } else if (ndim == 3*NTrk+3) {
2429 Amg::MatrixX D_mat(3*NTrk+6,2); D_mat.setZero();
2430 D_mat(0,0) = CONST*dTaudx;
2431 D_mat(1,0) = CONST*dTaudy;
2432 D_mat(2,0) = 0.;
2433 for( unsigned int it=0; it<NTrk; it++) {
2434 D_mat(3*it+3,0) = CONST*dTaudphi[it];
2435 D_mat(3*it+4,0) = CONST*dTaudtheta[it];
2436 D_mat(3*it+5,0) = CONST*dTaudqOverP[it];
2437 D_mat(3*it+3,1) = dMdphi[it];
2438 D_mat(3*it+4,1) = dMdtheta[it];
2439 D_mat(3*it+5,1) = dMdqOverP[it];
2440 }
2441 D_mat(3*NTrk+3,0) = CONST*dTaudx0;
2442 D_mat(3*NTrk+4,0) = CONST*dTaudy0;
2443 D_mat(3*NTrk+5,0) = 0.;
2444 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
2445 W_mat.block(0,0,ndim,ndim) = fullCov;
2446 W_mat.block<3,3>(3*NTrk+3,3*NTrk+3) = vertex->covariancePosition();
2447 V0_err = D_mat.transpose() * W_mat * D_mat;
2448 }else{
2449 ATH_MSG_ERROR("This should not happen - returning zero.");
2450 return V0_err;
2451 }
2452 return V0_err;
2453 }
2454
2455 Amg::MatrixX V0Tools::convertCovMatrix(const xAOD::Vertex * vxCandidate)
2456 {
2457 unsigned int NTrk = vxCandidate->nTrackParticles();
2458 const std::vector<float> &matrix = vxCandidate->covariance();
2459
2460 int ndim = 0;
2461
2462 if ( matrix.size() == (3*NTrk+3)*(3*NTrk+3+1)/2) {
2463 ndim = 3*NTrk+3;
2464 } else if (matrix.size() == (5*NTrk+3)*(5*NTrk+3+1)/2) {
2465 ndim = 5*NTrk+3;
2466 } else {
2467 return Amg::MatrixX(0,0);
2468 }
2469
2470 Amg::MatrixX mtx(ndim,ndim);
2471
2472 Eigen::Index ij=0;
2473 for (int i=1; i<= ndim; i++) {
2474 for (int j=1; j<=i; j++){
2475 if (i==j) {
2476 mtx(i-1,j-1)=matrix[ij];
2477 } else {
2478 mtx.fillSymmetric(i-1,j-1,matrix[ij]);
2479 }
2480 ij++;
2481 }
2482 }
2483 // NOTE: mtx is a pointer! Take care of deleting it after you do not
2484 // need it anymore!!!!
2485 return mtx;
2486 }
2487
2488}//end of namespace definitions
mag
Scalar mag() const
mag method
Definition AmgMatrixBasePlugin.h:26
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition EventPrimitives.h:50
AmgMatrix
#define AmgMatrix(rows, cols)
Definition EventPrimitives.h:49
Tp4
static Double_t Tp4(Double_t t)
Definition LArPhysWaveHECTool.cxx:607
P
static Double_t P(Double_t *tt, Double_t *par)
Definition LArPhysWaveHECTool.cxx:280
AthAlgTool::AthAlgTool
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Definition AthAlgTool.cxx:16
Trk::IExtrapolator
Interface class for the extrapolation AlgTool, it inherits from IAlgTool Detailed information about p...
Definition IExtrapolator.h:61
Trk::IExtrapolator::extrapolate
virtual std::unique_ptr< NeutralParameters > extrapolate(const NeutralParameters &parameters, const Surface &sf, PropDirection dir=anyDirection, const BoundaryCheck &bcheck=true) const =0
Main extrapolation Interface starting from neutral parameters and aiming at surface.
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::PerigeeSurface
Class describing the Line to which the Perigee refers to.
Definition PerigeeSurface.h:43
Trk::V0Tools::makeV0Cov
static Amg::MatrixX makeV0Cov(const xAOD::Vertex *vxCandidate)
Definition V0Tools.cxx:2304
Trk::V0Tools::V0Tools
V0Tools(const std::string &t, const std::string &n, const IInterface *p)
Default constructor due to Athena interface.
Definition V0Tools.cxx:26
Trk::V0Tools::phiStar
static double phiStar(const xAOD::Vertex *vxCandidate, double posTrackMass, double negTrackMass)
phiStar
Definition V0Tools.cxx:1834
Trk::V0Tools::positiveTrack4Momentum
static xAOD::TrackParticle::FourMom_t positiveTrack4Momentum(const xAOD::Vertex *vxCandidate, double posTrackMass)
Definition V0Tools.cxx:416
Trk::V0Tools::massErrorVxCandidate
double massErrorVxCandidate(const xAOD::Vertex *vxCandidate, double posTrackMass, double negTrackMass) const
Definition V0Tools.cxx:242
Trk::V0Tools::lxyzError
double lxyzError(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex) const
error on lxyz
Definition V0Tools.cxx:1072
Trk::V0Tools::ndof
static float ndof(const xAOD::Vertex *vxCandidate)
NDoF of the vertex fit.
Definition V0Tools.cxx:446
Trk::V0Tools::massErrorVKalVrt
double massErrorVKalVrt(const xAOD::Vertex *vxCandidate, double posTrackMass, double negTrackMass) const
Definition V0Tools.cxx:167
Trk::V0Tools::negativeTrack4Momentum
static xAOD::TrackParticle::FourMom_t negativeTrack4Momentum(const xAOD::Vertex *vxCandidate, double negTrackMass)
Definition V0Tools.cxx:426
Trk::V0Tools::invariantMassProbability
double invariantMassProbability(const xAOD::Vertex *vxCandidate, double V0Mass, double posTrackMass, double negTrackMass) const
Definition V0Tools.cxx:317
Trk::V0Tools::invariantMassBeforeFitIP
double invariantMassBeforeFitIP(const xAOD::Vertex *vxCandidate, std::span< const double > masses) const
create neutral TrackParticle from vxCandidate
Definition V0Tools.cxx:1936
Trk::V0Tools::vtx
static Amg::Vector3D vtx(const xAOD::Vertex *vxCandidate)
vertex position
Definition V0Tools.cxx:476
Trk::V0Tools::cosTheta
static double cosTheta(const xAOD::Vertex *vxCandidate, const Amg::Vector3D &vertex)
cosTheta (pointing angle to an Amg::Vector3D or an xAOD::Vertex)
Definition V0Tools.cxx:1861
Trk::V0Tools::separation
double separation(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex) const
statistical separation (Mahalanobis distance)
Definition V0Tools.cxx:637
Trk::V0Tools::negativeTrackMomentum
static Amg::Vector3D negativeTrackMomentum(const xAOD::Vertex *vxCandidate)
Definition V0Tools.cxx:384
Trk::V0Tools::lxyz
static double lxyz(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex)
projection of distance in 3D of the vertex wrt an xAOD::Vertex vertex along the momentum direction (P...
Definition V0Tools.cxx:1061
Trk::V0Tools::a0z
static double a0z(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex)
Definition V0Tools.cxx:675
Trk::V0Tools::cosThetaStar
static double cosThetaStar(const xAOD::Vertex *vxCandidate, double posTrackMass, double negTrackMass)
cosThetaStar
Definition V0Tools.cxx:1810
Trk::V0Tools::massErrorV0Fitter
double massErrorV0Fitter(const xAOD::Vertex *vxCandidate, double posTrackMass, double negTrackMass) const
Definition V0Tools.cxx:102
Trk::V0Tools::invariantMassBeforeFit
double invariantMassBeforeFit(const xAOD::Vertex *vxCandidate, std::span< const double > masses, const EventContext &, const Trk::IExtrapolator *) const
Definition V0Tools.cxx:1961
Trk::V0Tools::tau3DError
double tau3DError(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex, std::span< const double > masses) const
proper time error in 3D wrt an xAOD::Vertex vertex assuming track masses
Definition V0Tools.cxx:1542
Trk::V0Tools::thetaStar
static double thetaStar(const xAOD::Vertex *vxCandidate, double mass1, double mass2)
Polarization angles in helicity frame (using positive track):
Definition V0Tools.cxx:1794
Trk::V0Tools::tau3D
double tau3D(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex, std::span< const double > masses) const
proper time in 3D wrt an xAOD::Vertex vertex assuming track masses
Definition V0Tools.cxx:1518
Trk::V0Tools::massProbability
double massProbability(double V0Mass, double mass, double massErr) const
Definition V0Tools.cxx:346
Trk::V0Tools::negativeOrigTrack
static const xAOD::TrackParticle * negativeOrigTrack(const xAOD::Vertex *vxCandidate)
Definition V0Tools.cxx:1925
Trk::V0Tools::positiveOrigTrack
static const xAOD::TrackParticle * positiveOrigTrack(const xAOD::Vertex *vxCandidate)
Definition V0Tools.cxx:1914
Trk::V0Tools::initialize
StatusCode initialize()
Standard AlgTool methods.
Definition V0Tools.cxx:34
Trk::V0Tools::tau
double tau(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex, double posTrackMass, double negTrackMass) const
proper time wrt an xAOD::Vertex vertex assuming posTrackMass and negTrackMass tau = CONST*M*lxy/pT
Definition V0Tools.cxx:1190
Trk::V0Tools::tauMassCovariance
Amg::MatrixX tauMassCovariance(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex, std::span< const double > masses) const
mass-proper time covariance
Definition V0Tools.cxx:2323
Trk::V0Tools::lxy
static double lxy(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex)
projection of distance in xy of the vertex wrt an xAOD::Vertex vertex along the momentum direction (P...
Definition V0Tools.cxx:940
Trk::V0Tools::positiveTrackMomentum
static Amg::Vector3D positiveTrackMomentum(const xAOD::Vertex *vxCandidate)
Definition V0Tools.cxx:372
Trk::V0Tools::invariantMassError
double invariantMassError(const xAOD::Vertex *vxCandidate, double posTrackMass, double negTrackMass) const
Definition V0Tools.cxx:70
Trk::V0Tools::massTauCov
double massTauCov(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex, std::span< const double > masses) const
Definition V0Tools.cxx:2172
Trk::V0Tools::a0xy
static double a0xy(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex)
closest distance in Z and the transverse plane of the momentum vector to an xAOD::Vertex a0xy using t...
Definition V0Tools.cxx:667
Trk::V0Tools::lxyError
double lxyError(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex) const
error on lxy
Definition V0Tools.cxx:950
Trk::V0Tools::a0xyError
double a0xyError(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex) const
errors on a0xy and a0z
Definition V0Tools.cxx:811
Trk::V0Tools::a0Error
double a0Error(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex, bool in3D=true) const
error on a0
Definition V0Tools.cxx:816
Trk::V0Tools::charge
static float charge(const xAOD::Vertex *vxCandidate)
sum of the charges of the tracks in the vertex
Definition V0Tools.cxx:1895
Trk::V0Tools::invariantMassErrorBeforeFitIP
double invariantMassErrorBeforeFitIP(const xAOD::Vertex *vxCandidate, std::span< const double > masses) const
Definition V0Tools.cxx:2020
Trk::V0Tools::rxy_var
static double rxy_var(double dx, double dy, const Amg::MatrixX &cov)
rxy_var
Definition V0Tools.cxx:496
Trk::V0Tools::track4Momentum
static xAOD::TrackParticle::FourMom_t track4Momentum(const xAOD::Vertex *vxCandidate, unsigned int trkIndex, double trackMass)
Methods, returning the refitted 4-momenta of the positive and negative tracks and the V0 for a given ...
Definition V0Tools.cxx:406
Trk::V0Tools::pTError
double pTError(const xAOD::Vertex *vxCandidate) const
error on the transverse momentum of the V0
Definition V0Tools.cxx:547
Trk::V0Tools::origTrack
static const xAOD::TrackParticle * origTrack(const xAOD::Vertex *vxCandidate, int trkIndex)
pointers to original tracks
Definition V0Tools.cxx:1909
Trk::V0Tools::~V0Tools
~V0Tools()
Virtual destructor.
Trk::V0Tools::vertexProbability
double vertexProbability(const xAOD::Vertex *vxCandidate) const
Probability of the vertex fit.
Definition V0Tools.cxx:456
Trk::V0Tools::invariantMass
double invariantMass(const xAOD::Vertex *vxCandidate, double posTrackMass, double negTrackMass) const
Methods, returning the invariant mass, error on the invariant mass and Chi2 probability of the invari...
Definition V0Tools.cxx:40
Trk::V0Tools::rxy
static double rxy(const xAOD::Vertex *vxCandidate)
Rxy of the vertex.
Definition V0Tools.cxx:481
Trk::V0Tools::tauError
double tauError(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex, double posTrackMass, double negTrackMass) const
proper time error wrt an xAOD::Vertex vertex assuming posTrackMass and negTrackMass tau = CONST*M*lxy...
Definition V0Tools.cxx:1245
Trk::V0Tools::pca
static Amg::Vector3D pca(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex)
point of closest approach of the momentum vector to an xAOD::Vertex
Definition V0Tools.cxx:622
Trk::V0Tools::cosTheta_xy
static double cosTheta_xy(const xAOD::Vertex *vxCandidate, const Amg::Vector3D &vertex)
cosTheta (pointing angle to an Amg::Vector3D or an xAOD::Vertex in transverse plane)
Definition V0Tools.cxx:1877
Trk::V0Tools::a0
static double a0(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex)
closest distance of the momentum vector to an xAOD::Vertex
Definition V0Tools.cxx:683
Trk::V0Tools::V04Momentum
static xAOD::TrackParticle::FourMom_t V04Momentum(const xAOD::Vertex *vxCandidate, double V0Mass)
Definition V0Tools.cxx:436
Trk::V0Tools::invariantMassErrorBeforeFit
double invariantMassErrorBeforeFit(const xAOD::Vertex *vxCandidate, std::span< const double > masses, const EventContext &ctx, const Trk::IExtrapolator *) const
Definition V0Tools.cxx:2087
Trk::V0Tools::chisq
static float chisq(const xAOD::Vertex *vxCandidate)
Chisq of the vertex fit.
Definition V0Tools.cxx:451
Trk::V0Tools::rxyError
double rxyError(const xAOD::Vertex *vxCandidate) const
error on Rxy
Definition V0Tools.cxx:510
Trk::V0Tools::a0zError
double a0zError(const xAOD::Vertex *vxCandidate, const xAOD::Vertex *vertex) const
Definition V0Tools.cxx:691
Trk::V0Tools::trackMomentum
static Amg::Vector3D trackMomentum(const xAOD::Vertex *vxCandidate, unsigned int trkIndex)
Methods, returning the refitted 3-momenta of the positive and negative tracks and the V0 (for more th...
Definition V0Tools.cxx:366
Trk::V0Tools::pT
static double pT(const xAOD::Vertex *vxCandidate)
transverse momentum of the V0
Definition V0Tools.cxx:542
Trk::V0Tools::convertCovMatrix
static Amg::MatrixX convertCovMatrix(const xAOD::Vertex *vxCandidate)
Definition V0Tools.cxx:2455
Trk::V0Tools::V0Momentum
static Amg::Vector3D V0Momentum(const xAOD::Vertex *vxCandidate)
Definition V0Tools.cxx:396
xAOD::TrackParticle_v1::theta
float theta() const
Returns the parameter, which has range 0 to .
xAOD::TrackParticle_v1::p4
virtual FourMom_t p4() const override final
The full 4-momentum of the particle.
Definition TrackParticle_v1.cxx:122
xAOD::TrackParticle_v1::perigeeParameters
const Trk::Perigee & perigeeParameters() const
Returns the Trk::MeasuredPerigee track parameters.
Definition TrackParticle_v1.cxx:437
xAOD::TrackParticle_v1::FourMom_t
IParticle::FourMom_t FourMom_t
Definition of the 4-momentum type.
Definition TrackParticle_v1.h:73
xAOD::TrackParticle_v1::phi
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .).
xAOD::TrackParticle_v1::definingParametersCovMatrix
const ParametersCovMatrix_t definingParametersCovMatrix() const
Returns the 5x5 symmetric matrix containing the defining parameters covariance matrix.
Definition TrackParticle_v1.cxx:260
xAOD::TrackParticle_v1::qOverP
float qOverP() const
Returns the parameter.
xAOD::TrackParticle_v1::charge
float charge() const
Returns the charge.
Definition TrackParticle_v1.cxx:143
xAOD::Vertex_v1::nTrackParticles
size_t nTrackParticles() const
Get the number of tracks associated with this vertex.
Definition Vertex_v1.cxx:292
xAOD::Vertex_v1::trackParticle
const TrackParticle * trackParticle(size_t i) const
Get the pointer to a given track that was used in vertex reco.
Definition Vertex_v1.cxx:271
xAOD::Vertex_v1::numberDoF
float numberDoF() const
Returns the number of degrees of freedom of the vertex fit as float.
xAOD::Vertex_v1::chiSquared
float chiSquared() const
Returns the of the vertex fit as float.
xAOD::Vertex_v1::vxTrackAtVertex
std::vector< Trk::VxTrackAtVertex > & vxTrackAtVertex()
Non-const access to the VxTrackAtVertex vector.
Definition Vertex_v1.cxx:203
xAOD::Vertex_v1::position
const Amg::Vector3D & position() const
Returns the 3-pos.
Definition Vertex_v1.cxx:106
xAOD::Vertex_v1::covariance
const std::vector< float > & covariance() const
Returns the covariance matrix as a simple vector of values.
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition comparitor.cxx:525
r
int r
Definition globals.cxx:22
a0
double a0
Definition globals.cxx:27
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition EventPrimitives.h:27
Amg::SymMatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > SymMatrixX
Definition EventPrimitives.h:28
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition FakeTrackBuilder.h:9
Trk::pz
@ pz
global momentum (cartesian)
Definition ParamDefs.h:61
Trk::theta
@ theta
Definition ParamDefs.h:66
Trk::qOverP
@ qOverP
perigee
Definition ParamDefs.h:67
Trk::phi
@ phi
Definition ParamDefs.h:75
Trk::px
@ px
Definition ParamDefs.h:59
Trk::py
@ py
Definition ParamDefs.h:60
Trk::TrackParameters
ParametersBase< TrackParametersDim, Charged > TrackParameters
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:27
xAOD::TrackParticle
TrackParticle_v1 TrackParticle
Reference the current persistent version:
Definition Event/xAOD/xAODTracking/xAODTracking/TrackParticle.h:13
xAOD::Vertex
Vertex_v1 Vertex
Define the latest version of the vertex class.
Definition Event/xAOD/xAODTracking/xAODTracking/Vertex.h:16
unlikely
#define unlikely(x)
Definition pythonic_coracool.cxx:9
Generated on for ATLAS Offline Software by doxygen 1.17.0