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