ATLAS Offline Software
Loading...
Searching...
No Matches
CascadeTools.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
3*/
5// CascadeTools.cxx, (c) ATLAS Detector software
7#include "DerivationFrameworkBPhys/CascadeTools.h"
8#include "CLHEP/GenericFunctions/CumulativeChiSquare.hh"
9
10
11namespace DerivationFramework {
12
13
14CascadeTools::CascadeTools(const std::string& t, const std::string& n, const IInterface* p) :
15 AthAlgTool(t,n,p)
16{
17}
18
19CascadeTools::~CascadeTools() {}
20
21//Light speed constant for various calculations
22constexpr double s_CONST = 1000./299.792;
23
24double CascadeTools::invariantMass(const std::vector<TLorentzVector> &particleMom) const
25{
26 if(particleMom.size() == 0) return -999999.;
27 TLorentzVector totalMom;
28 unsigned int NTrk = particleMom.size();
29 for( unsigned int it=0; it<NTrk; it++) totalMom += particleMom[it];
30 return totalMom.M();
31}
32
33double CascadeTools::invariantMass(const std::vector<TLorentzVector> &particleMom2, const std::vector<double> &masses) const
34{
35 if(particleMom2.size() == 0) return -999999.;
36 TLorentzVector totalMom;
37 unsigned int NTrk = particleMom2.size();
38 if (masses.size() != NTrk) {
39 ATH_MSG_DEBUG("The provided number of masses does not match the number of tracks in the vertex");
40 return -999999.;
41 }
42 TLorentzVector temp;
43 for( unsigned int it=0; it<NTrk; it++) {
44 double esq = particleMom2[it].Px()*particleMom2[it].Px() + particleMom2[it].Py()*particleMom2[it].Py() +
45 particleMom2[it].Pz()*particleMom2[it].Pz() + masses[it]*masses[it];
46 double e = (esq>0.) ? sqrt(esq) : 0.;
47
48 temp.SetPxPyPzE(particleMom2[it].Px(),particleMom2[it].Py(),particleMom2[it].Pz(),e);
49 totalMom += temp;
50 }
51 return totalMom.M();
52}
53
54double CascadeTools::invariantMassError(const std::vector<TLorentzVector> &particleMom, const Amg::MatrixX& cov) const
55{
56 if(particleMom.size() == 0) return -999999.;
57 unsigned int NTrk = particleMom.size();
58 TLorentzVector totalMom;
59 for( unsigned int it=0; it<NTrk; it++) totalMom += particleMom[it];
60
61 Amg::MatrixX D_vec(3*NTrk+3,1); D_vec.setZero();
62 for( unsigned int it=0; it<NTrk; it++) {
63 D_vec(3*it+3) = 2.*(totalMom.E()*particleMom[it].Px()/particleMom[it].E()-totalMom.Px());
64 D_vec(3*it+4) = 2.*(totalMom.E()*particleMom[it].Py()/particleMom[it].E()-totalMom.Py());
65 D_vec(3*it+5) = 2.*(totalMom.E()*particleMom[it].Pz()/particleMom[it].E()-totalMom.Pz());
66 }
67 Amg::MatrixX merr = D_vec.transpose() * cov * D_vec;
68 double massVarsq = merr(0,0);
69 if (massVarsq <= 0.) ATH_MSG_DEBUG("massError: negative sqrt massVarsq " << massVarsq);
70 double massVar = (massVarsq>0.) ? sqrt(massVarsq) : 0.;
71 double massErr = massVar/(2.*totalMom.M());
72 return massErr;
73}
74
75double CascadeTools::invariantMassError(const std::vector<TLorentzVector> &particleMom2, const Amg::MatrixX& cov, const std::vector<double> &masses) const
76{
77 if(particleMom2.size() == 0) return -999999.;
78 unsigned int NTrk = particleMom2.size();
79 if (masses.size() != NTrk) {
80 ATH_MSG_DEBUG("The provided number of masses does not match the number of tracks in the vertex");
81 return -999999.;
82 }
83 std::vector<TLorentzVector> particleMom(NTrk);
84 for( unsigned int it=0; it<NTrk; it++) {
85 double esq = particleMom2[it].Px()*particleMom2[it].Px() + particleMom2[it].Py()*particleMom2[it].Py() +
86 particleMom2[it].Pz()*particleMom2[it].Pz() + masses[it]*masses[it];
87 double e = (esq>0.) ? sqrt(esq) : 0.;
88 particleMom[it].SetPxPyPzE(particleMom2[it].Px(),particleMom2[it].Py(),particleMom2[it].Pz(),e);
89 }
90 TLorentzVector totalMom;
91 for( unsigned int it=0; it<NTrk; it++) totalMom += particleMom[it];
92
93 std::vector<double>dm2dpx(NTrk), dm2dpy(NTrk), dm2dpz(NTrk);
94 for( unsigned int it=0; it<NTrk; it++) {
95 dm2dpx[it] = 2.*(totalMom.E()*particleMom[it].Px()/particleMom[it].E()-totalMom.Px());
96 dm2dpy[it] = 2.*(totalMom.E()*particleMom[it].Py()/particleMom[it].E()-totalMom.Py());
97 dm2dpz[it] = 2.*(totalMom.E()*particleMom[it].Pz()/particleMom[it].E()-totalMom.Pz());
98 }
99 Amg::MatrixX D_vec(3*NTrk+3,1); D_vec.setZero();
100 for( unsigned int it=0; it<NTrk; it++) {
101 D_vec(3*it+3) = dm2dpx[it];
102 D_vec(3*it+4) = dm2dpy[it];
103 D_vec(3*it+5) = dm2dpz[it];
104 }
105 Amg::MatrixX merr = D_vec.transpose() * cov * D_vec;
106 double massVarsq = merr(0,0);
107 if (massVarsq <= 0.) ATH_MSG_DEBUG("massError: negative sqrt massVarsq " << massVarsq);
108 double massVar = (massVarsq>0.) ? sqrt(massVarsq) : 0.;
109 double massErr = massVar/(2.*totalMom.M());
110 return massErr;
111}
112
113double CascadeTools::pT(const std::vector<TLorentzVector> &particleMom) const
114{
115 if(particleMom.size() == 0) return -999999.;
116 Amg::Vector3D P = momentum(particleMom);;
117 return P.perp();
118}
119
120double CascadeTools::pTError(const std::vector<TLorentzVector> &particleMom, const Amg::MatrixX& cov) const
121{
122 if(particleMom.size() == 0) return -999999.;
123 Amg::Vector3D P = momentum(particleMom);;
124 double Px = P.x();
125 double Py = P.y();
126 double PT = P.perp();
127
128 unsigned int NTrk = particleMom.size();
129 Amg::MatrixX D_vec(3*NTrk+3,1); D_vec.setZero();
130 for( unsigned int it=0; it<NTrk; it++) {
131 D_vec(3*it+3) = Px/PT;
132 D_vec(3*it+4) = Py/PT;
133 D_vec(3*it+5) = 0.;
134 }
135 Amg::MatrixX PtErrSq = D_vec.transpose() * cov * D_vec;
136 double PtErrsq = PtErrSq(0,0);
137 if (PtErrsq <= 0.) ATH_MSG_DEBUG("ptError: negative sqrt PtErrsq " << PtErrsq);
138 return (PtErrsq>0.) ? sqrt(PtErrsq) : 0.;
139}
140
141double CascadeTools::lxy(const std::vector<TLorentzVector> &particleMom, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
142{
143 if(particleMom.size() == 0) return -999999.;
144 auto vert = SV->position() - PV->position();
145 double dx = vert.x();
146 double dy = vert.y();
147 Amg::Vector3D P = momentum(particleMom);;
148 double dxy = (P.x()*dx + P.y()*dy)/P.perp();
149 return dxy;
150}
151
152double CascadeTools::lxyError(const std::vector<TLorentzVector> &particleMom, const Amg::MatrixX& cov, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
153{
154 if(particleMom.size() == 0) return -999999.;
155 auto vert = SV->position() - PV->position();
156 double dx = vert.x();
157 double dy = vert.y();
158 Amg::Vector3D P = momentum(particleMom);;
159 double Px = P.x();
160 double Py = P.y();
161 double PT = P.perp();
162 double LXYoverPT = (Px*dx+Py*dy)/(PT*PT);
163
164 unsigned int NTrk = particleMom.size();
165
166 double dLxydx = Px/PT;
167 double dLxydy = Py/PT;
168 double dLxydx0 = -dLxydx;
169 double dLxydy0 = -dLxydy;
170
171 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
172 D_vec(0) = dLxydx;
173 D_vec(1) = dLxydy;
174 D_vec(2) = 0.;
175 for( unsigned int it=0; it<NTrk; it++) {
176 D_vec(3*it+3) = (dx - LXYoverPT*Px)/PT;
177 D_vec(3*it+4) = (dy - LXYoverPT*Py)/PT;
178 D_vec(3*it+5) = 0.;
179 }
180 D_vec(3*NTrk+3) = dLxydx0;
181 D_vec(3*NTrk+4) = dLxydy0;
182 D_vec(3*NTrk+5) = 0.;
183
184 unsigned int ndim = 3*NTrk+3;
185 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
186 W_mat.block(0,0,ndim,ndim) = cov;
187 W_mat.block(3*NTrk+3,3*NTrk+3,3,3) = PV->covariancePosition();
188 Amg::MatrixX V_err = D_vec.transpose() * W_mat * D_vec;
189
190 double LxyErrsq = V_err(0,0);
191 if (LxyErrsq <= 0.) ATH_MSG_DEBUG("lxyError: negative sqrt LxyErrsq " << LxyErrsq);
192 return (LxyErrsq>0.) ? sqrt(LxyErrsq) : 0.;
193}
194
195double CascadeTools::tau(const std::vector<TLorentzVector> &particleMom, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
196{
197 if(particleMom.size() == 0) return -999999.;
198 double M = invariantMass(particleMom);
199 double LXY = lxy(particleMom,SV,PV);
200 double PT = pT(particleMom);
201 return s_CONST*M*LXY/PT;
202}
203
204double CascadeTools::tauError(const std::vector<TLorentzVector> &particleMom, const Amg::MatrixX& cov, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
205{
206 if(particleMom.size() == 0) return -999999.;
207 double M = invariantMass(particleMom);
208 auto vert = SV->position() - PV->position();
209 double dx = vert.x();
210 double dy = vert.y();
211 Amg::Vector3D P = momentum(particleMom);;
212 double Px = P.x();
213 double Py = P.y();
214 double PT = P.perp();
215 double LXY = Px*dx+Py*dy;
216
217 unsigned int NTrk = particleMom.size();
218 TLorentzVector totalMom;
219 for( unsigned int it=0; it<NTrk; it++) totalMom += particleMom[it];
220
221 double dTaudx = (M*Px)/(PT*PT);
222 double dTaudy = (M*Py)/(PT*PT);
223 double dTaudx0 = -dTaudx;
224 double dTaudy0 = -dTaudy;
225
226 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
227 D_vec(0) = dTaudx;
228 D_vec(1) = dTaudy;
229 D_vec(2) = 0.;
230 for( unsigned int it=0; it<NTrk; it++) {
231 D_vec(3*it+3) = (((totalMom.E()*particleMom[it].Px()*LXY)/(M*particleMom[it].E()))-Px*LXY/M+M*dx)/(PT*PT) - 2.*M*LXY*Px/(PT*PT*PT*PT);
232 D_vec(3*it+4) = (((totalMom.E()*particleMom[it].Py()*LXY)/(M*particleMom[it].E()))-Py*LXY/M+M*dy)/(PT*PT) - 2.*M*LXY*Py/(PT*PT*PT*PT);
233 D_vec(3*it+5) = 0.;
234 }
235 D_vec(3*NTrk+3) = dTaudx0;
236 D_vec(3*NTrk+4) = dTaudy0;
237 D_vec(3*NTrk+5) = 0.;
238
239 unsigned int ndim = 3*NTrk+3;
240 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
241 W_mat.block(0,0,ndim,ndim) = cov;
242 W_mat.block(3*NTrk+3,3*NTrk+3,3,3) = PV->covariancePosition();
243 Amg::MatrixX V_err = D_vec.transpose() * W_mat * D_vec;
244
245 double tauErrsq = V_err(0,0);
246 if (tauErrsq <= 0.) ATH_MSG_DEBUG("tauError: negative sqrt tauErrsq " << tauErrsq);
247 double tauErr = (tauErrsq>0.) ? sqrt(tauErrsq) : 0.;
248 return s_CONST*tauErr;
249}
250
251double CascadeTools::tau(const std::vector<TLorentzVector> &particleMom, const xAOD::Vertex* SV, const xAOD::Vertex* PV, double M) const
252{
253 if(particleMom.size() == 0) return -999999.;
254 double LXY = lxy(particleMom,SV,PV);
255 double PT = pT(particleMom);
256 return s_CONST*M*LXY/PT;
257}
258
259double CascadeTools::tauError(const std::vector<TLorentzVector> &particleMom, const Amg::MatrixX& cov, const xAOD::Vertex* SV, const xAOD::Vertex* PV, double M) const
260{
261 if(particleMom.size() == 0) return -999999.;
262 auto vert = SV->position() - PV->position();
263 double dx = vert.x();
264 double dy = vert.y();
265 Amg::Vector3D P = momentum(particleMom);;
266 double Px = P.x();
267 double Py = P.y();
268 double PT = P.perp();
269 double LXY = Px*dx+Py*dy;
270
271 unsigned int NTrk = particleMom.size();
272 TLorentzVector totalMom;
273 for( unsigned int it=0; it<NTrk; it++) totalMom += particleMom[it];
274
275 double dTaudx = (M*Px)/(PT*PT);
276 double dTaudy = (M*Py)/(PT*PT);
277 double dTaudx0 = -dTaudx;
278 double dTaudy0 = -dTaudy;
279
280 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
281 D_vec(0) = dTaudx;
282 D_vec(1) = dTaudy;
283 D_vec(2) = 0.;
284 for( unsigned int it=0; it<NTrk; it++) {
285 D_vec(3*it+3) = (M*dx)/(PT*PT) - 2.*M*LXY*Px/(PT*PT*PT*PT);
286 D_vec(3*it+4) = (M*dy)/(PT*PT) - 2.*M*LXY*Py/(PT*PT*PT*PT);
287 D_vec(3*it+5) = 0.;
288 }
289 D_vec(3*NTrk+3) = dTaudx0;
290 D_vec(3*NTrk+4) = dTaudy0;
291 D_vec(3*NTrk+5) = 0.;
292
293 unsigned int ndim = 3*NTrk+3;
294 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
295 W_mat.block(0,0,ndim,ndim) = cov;
296 W_mat.block(3*NTrk+3,3*NTrk+3,3,3) = PV->covariancePosition();
297 Amg::MatrixX V_err = D_vec.transpose() * W_mat * D_vec;
298
299 double tauErrsq = V_err(0,0);
300 if (tauErrsq <= 0.) ATH_MSG_DEBUG("tauError: negative sqrt tauErrsq " << tauErrsq);
301 double tauErr = (tauErrsq>0.) ? sqrt(tauErrsq) : 0.;
302 return s_CONST*tauErr;
303}
304
305Amg::Vector3D CascadeTools::pca(const std::vector<TLorentzVector> &particleMom, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
306{
307 if(particleMom.size() == 0) {
308 Amg::Vector3D p; p.setZero();
309 return p;
310 }
311 Amg::Vector3D pv = PV->position();
312 Amg::Vector3D sv = SV->position();
313 Amg::Vector3D P = momentum(particleMom);;
314 double p2 = P.mag2();
315 double pdr = P.dot((sv - pv));
316 return sv - P*pdr/p2;
317}
318
319double CascadeTools::cosTheta(const std::vector<TLorentzVector> &particleMom, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
320{
321 if(particleMom.size() == 0) return -999999.;
322 Amg::Vector3D P = momentum(particleMom);;
323 Amg::Vector3D vtx = SV->position();
324 vtx -= PV->position();
325 return (P.dot(vtx))/(P.mag()*vtx.mag());
326}
327
328double CascadeTools::cosTheta_xy(const std::vector<TLorentzVector> &particleMom, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
329{
330 if(particleMom.size() == 0) return -999999.;
331 Amg::Vector3D P = momentum(particleMom);;
332 Amg::Vector3D vtx = SV->position();
333 vtx -= PV->position();
334 double pT = P.perp();
335 return (P.x()*vtx.x()+P.y()*vtx.y())/(pT*vtx.perp());
336}
337
338double CascadeTools::a0z(const std::vector<TLorentzVector> &particleMom, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
339{
340 if(particleMom.size() == 0) return -999999.;
341 Amg::Vector3D pv = PV->position();
342 Amg::Vector3D ca_point = pca(particleMom,SV,PV);
343 Amg::Vector3D a0_vec = pv - ca_point;
344 return a0_vec.z();
345}
346
347double CascadeTools::a0zError(const std::vector<TLorentzVector> &particleMom, const Amg::MatrixX& cov, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
348{
349 if(particleMom.size() == 0) return -999999.;
350 auto vert = SV->position() - PV->position();
351 double dx = vert.x();
352 double dy = vert.y();
353 double dz = vert.z();
354 Amg::Vector3D P = momentum(particleMom);;
355 double Px = P.x();
356 double Py = P.y();
357 double Pz = P.z();
358 double P2 = P.mag2();
359 double L = Px*dx+Py*dy+Pz*dz;
360
361 unsigned int NTrk = particleMom.size();
362
363
364 double da0zdx = (Px*Pz)/P2;
365 double da0zdy = (Py*Pz)/P2;
366 double da0zdz = (Pz*Pz)/P2 - 1.;
367 double da0zdx0 = -da0zdx;
368 double da0zdy0 = -da0zdy;
369 double da0zdz0 = -da0zdz;
370
371 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
372 D_vec(0) = da0zdx;
373 D_vec(1) = da0zdy;
374 D_vec(2) = da0zdz;
375 for( unsigned int it=0; it<NTrk; it++) {
376 D_vec(3*it+3) = (Pz*(P2*dx-2.*L*Px))/(P2*P2);
377 D_vec(3*it+4) = (Pz*(P2*dy-2.*L*Py))/(P2*P2);
378 D_vec(3*it+5) = (Pz*(P2*dz-2.*L*Pz))/(P2*P2)+L/P2;
379 }
380 D_vec(3*NTrk+3) = da0zdx0;
381 D_vec(3*NTrk+4) = da0zdy0;
382 D_vec(3*NTrk+5) = da0zdz0;
383
384 unsigned int ndim = 3*NTrk+3;
385 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
386 W_mat.block(0,0,ndim,ndim) = cov;
387 W_mat.block(3*NTrk+3,3*NTrk+3,3,3) = PV->covariancePosition();
388 Amg::MatrixX V_err = D_vec.transpose() * W_mat * D_vec;
389
390 double a0zErrsq = V_err(0,0);
391 if (a0zErrsq <= 0.) ATH_MSG_DEBUG("a0zError: negative sqrt a0zErrsq " << a0zErrsq);
392 return (a0zErrsq>0.) ? sqrt(a0zErrsq) : 0.;
393}
394
395double CascadeTools::a0xy(const std::vector<TLorentzVector> &particleMom, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
396{
397 if(particleMom.size() == 0) return -999999.;
398 double cosineTheta_xy = cosTheta_xy(particleMom,SV,PV);
399 double sinTheta_xy = ((1.-cosineTheta_xy*cosineTheta_xy)>0.) ? sqrt((1.-cosineTheta_xy*cosineTheta_xy)) : 0.;
400 return (SV->position()-PV->position()).perp() * sinTheta_xy;
401}
402
403double CascadeTools::a0xyError(const std::vector<TLorentzVector> &particleMom, const Amg::MatrixX& cov, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
404{
405 if(particleMom.size() == 0) return -999999.;
406 auto vert = SV->position() - PV->position();
407 double dx = vert.x();
408 double dy = vert.y();
409 Amg::Vector3D P = momentum(particleMom);;
410 double Px = P.x();
411 double Py = P.y();
412 double P2 = P.perp()*P.perp();
413 double L = Px*dx+Py*dy;
414 double dR2 = vert.perp()*vert.perp();
415 double d = sqrt((P2*dR2-L*L)/P2);
416
417 unsigned int NTrk = particleMom.size();
418
419 double da0dx = (P2*dx-L*Px)/(P2*d);
420 double da0dy = (P2*dy-L*Py)/(P2*d);
421 double da0dx0 = -da0dx;
422 double da0dy0 = -da0dy;
423
424 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
425 D_vec(0) = da0dx;
426 D_vec(1) = da0dy;
427 D_vec(2) = 0.;
428 for( unsigned int it=0; it<NTrk; it++) {
429 D_vec(3*it+3) = (L*(L*Px-P2*dx))/(P2*P2*d);
430 D_vec(3*it+4) = (L*(L*Py-P2*dy))/(P2*P2*d);
431 D_vec(3*it+5) = 0.;
432 }
433 D_vec(3*NTrk+3) = da0dx0;
434 D_vec(3*NTrk+4) = da0dy0;
435 D_vec(3*NTrk+5) = 0.;
436
437 unsigned int ndim = 3*NTrk+3;
438 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
439 W_mat.block(0,0,ndim,ndim) = cov;
440 W_mat.block(3*NTrk+3,3*NTrk+3,3,3) = PV->covariancePosition();
441 Amg::MatrixX V_err = D_vec.transpose() * W_mat * D_vec;
442
443 double a0xyErrsq = V_err(0,0);
444 if (a0xyErrsq <= 0.) ATH_MSG_DEBUG("a0xyError: negative sqrt a0xyErrsq " << a0xyErrsq);
445 return (a0xyErrsq>0.) ? sqrt(a0xyErrsq) : 0.;
446}
447
448double CascadeTools::a0(const std::vector<TLorentzVector> &particleMom, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
449{
450 if(particleMom.size() == 0) return -999999.;
451 double cosineTheta = cosTheta(particleMom,SV,PV);
452 double sinTheta = ((1.-cosineTheta*cosineTheta)>0.) ? sqrt((1.-cosineTheta*cosineTheta)) : 0.;
453 return (SV->position()-PV->position()).mag() * sinTheta;
454}
455
456double CascadeTools::a0Error(const std::vector<TLorentzVector> &particleMom, const Amg::MatrixX& cov, const xAOD::Vertex* SV, const xAOD::Vertex* PV) const
457{
458 if(particleMom.size() == 0) return -999999.;
459 auto vert = SV->position() - PV->position();
460 double dx = vert.x();
461 double dy = vert.y();
462 double dz = vert.z();
463 Amg::Vector3D P = momentum(particleMom);;
464 double Px = P.x();
465 double Py = P.y();
466 double Pz = P.z();
467 double P2 = P.mag2();
468 double L = Px*dx+Py*dy+Pz*dz;
469 double dR2 = vert.mag2();
470 double d = sqrt((P2*dR2-L*L)/P2);
471
472 unsigned int NTrk = particleMom.size();
473
474 double da0dx = (P2*dx-L*Px)/(P2*d);
475 double da0dy = (P2*dy-L*Py)/(P2*d);
476 double da0dz = (P2*dz-L*Pz)/(P2*d);
477 double da0dx0 = -da0dx;
478 double da0dy0 = -da0dy;
479 double da0dz0 = -da0dz;
480
481 Amg::MatrixX D_vec(3*NTrk+6,1); D_vec.setZero();
482 D_vec(0) = da0dx;
483 D_vec(1) = da0dy;
484 D_vec(2) = da0dz;
485 for( unsigned int it=0; it<NTrk; it++) {
486 D_vec(3*it+3) = (L*(L*Px-P2*dx))/(P2*P2*d);
487 D_vec(3*it+4) = (L*(L*Py-P2*dy))/(P2*P2*d);
488 D_vec(3*it+5) = (L*(L*Pz-P2*dz))/(P2*P2*d);
489 }
490 D_vec(3*NTrk+3) = da0dx0;
491 D_vec(3*NTrk+4) = da0dy0;
492 D_vec(3*NTrk+5) = da0dz0;
493
494 unsigned int ndim = 3*NTrk+3;
495 Amg::MatrixX W_mat(3*NTrk+6,3*NTrk+6); W_mat.setZero();
496 W_mat.block(0,0,ndim,ndim) = cov;
497 W_mat.block(3*NTrk+3,3*NTrk+3,3,3) = PV->covariancePosition();
498 Amg::MatrixX V_err = D_vec.transpose() * W_mat * D_vec;
499
500 double a0Errsq = V_err(0,0);
501 if (a0Errsq <= 0.) ATH_MSG_DEBUG("a0Error: negative sqrt a0Errsq " << a0Errsq);
502 return (a0Errsq>0.) ? sqrt(a0Errsq) : 0.;
503}
504
505Amg::Vector3D CascadeTools::momentum(const std::vector<TLorentzVector> &particleMom) const
506{
507 if(particleMom.size() == 0) {
508 Amg::Vector3D p; p.setZero();
509 return p;
510 }
511 TLorentzVector totalMom;
512 unsigned int NTrk = particleMom.size();
513 for( unsigned int it=0; it<NTrk; it++) totalMom += particleMom[it];
514 TVector3 P3 = totalMom.Vect();
515 Amg::Vector3D mom(P3.Px(),P3.Py(),P3.Pz());
516 return mom;
517}
518
519double CascadeTools::massProbability(double V0Mass, double mass, double massErr) const
520{
521 if(massErr > 0.) {
522 double chi2 = (V0Mass - mass)*(V0Mass - mass)/(massErr*massErr);
523 int ndf = 1;
524 Genfun::CumulativeChiSquare myCumulativeChiSquare(ndf);
525 if (chi2 > 0.) {
526 double achi2prob = 1.-myCumulativeChiSquare(chi2);
527 return achi2prob;
528 } else {
529 ATH_MSG_DEBUG("chi2 <= 0");
530 return -1.;
531 }
532 } else {
533 return -1.;
534 }
535}
536
537double CascadeTools::vertexProbability(int ndf, double chi2) const
538{
539 if (ndf > 0.) {
540 Genfun::CumulativeChiSquare myCumulativeChiSquare(ndf);
541 if (chi2 > 0.) {
542 double chi2prob = 1.-myCumulativeChiSquare(chi2);
543 return chi2prob;
544 } else {
545 ATH_MSG_DEBUG("chi2 <= 0");
546 return -1.;
547 }
548 } else {
549 ATH_MSG_DEBUG("ndf <= 0");
550 return -1.;
551 }
552}
553
554
555Amg::MatrixX * CascadeTools::convertCovMatrix(const xAOD::Vertex * vxCandidate) const
556{
557 unsigned int NTrk = vxCandidate->nTrackParticles();
558 std::vector<float> matrix = vxCandidate->covariance();
559
560 int ndim = 0;
561
562 if ( matrix.size() == (3*NTrk+3)*(3*NTrk+3+1)/2) {
563 ndim = 3*NTrk+3;
564 } else if (matrix.size() == (5*NTrk+3)*(5*NTrk+3+1)/2) {
565 ndim = 5*NTrk+3;
566 } else {
567 return nullptr;
568 }
569
570 Amg::MatrixX* mtx = new Amg::MatrixX(ndim,ndim);
571
572 long int ij=0;
573 for (int i=1; i<= ndim; i++) {
574 for (int j=1; j<=i; j++){
575 if (i==j) {
576 (*mtx)(i-1,j-1)=matrix[ij];
577 } else {
578 (*mtx).fillSymmetric(i-1,j-1,matrix[ij]);
579 }
580 ij++;
581 }
582 }
583 // NOTE: mtx is a pointer! Take care of deleting it after you do not need it anymore!!!
584
585 return mtx;
586}
587
588Amg::MatrixX CascadeTools::SetFullMatrix(int NTrk, const std::vector<float> & Matrix) const
589{
590
591 Amg::MatrixX mtx(3+3*NTrk,3+3*NTrk); // Create identity matrix of needed size
592
593 int ij=0;
594
595 for(int i=0; i<(3+3*NTrk); i++){
596 for(int j=0; j<=i; j++){
597 mtx(i,j)=Matrix[ij];
598 if(i!=j) mtx(j,i)=Matrix[ij];
599 ij++;
600 }
601 }
602
603 return mtx;
604}
605
606} //end of namespace definitions
607
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
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
DerivationFramework::CascadeTools::invariantMass
double invariantMass(const std::vector< TLorentzVector > &moms, const std::vector< double > &masses) const
Definition CascadeTools.cxx:33
DerivationFramework::CascadeTools::pT
double pT(const std::vector< TLorentzVector > &moms) const
Definition CascadeTools.cxx:113
DerivationFramework::CascadeTools::invariantMassError
double invariantMassError(const std::vector< TLorentzVector > &moms, const Amg::MatrixX &cov, const std::vector< double > &masses) const
Definition CascadeTools.cxx:75
DerivationFramework::CascadeTools::a0Error
double a0Error(const std::vector< TLorentzVector > &particleMom, const Amg::MatrixX &cov, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:456
DerivationFramework::CascadeTools::lxy
double lxy(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:141
DerivationFramework::CascadeTools::momentum
Amg::Vector3D momentum(const std::vector< TLorentzVector > &particleMom) const
Definition CascadeTools.cxx:505
DerivationFramework::CascadeTools::a0xyError
double a0xyError(const std::vector< TLorentzVector > &particleMom, const Amg::MatrixX &cov, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:403
DerivationFramework::CascadeTools::pTError
double pTError(const std::vector< TLorentzVector > &moms, const Amg::MatrixX &cov) const
Definition CascadeTools.cxx:120
DerivationFramework::CascadeTools::a0xy
double a0xy(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:395
DerivationFramework::CascadeTools::~CascadeTools
~CascadeTools()
Virtual destructor.
Definition CascadeTools.cxx:19
DerivationFramework::CascadeTools::convertCovMatrix
Amg::MatrixX * convertCovMatrix(const xAOD::Vertex *vxCandidate) const
Definition CascadeTools.cxx:555
DerivationFramework::CascadeTools::tauError
double tauError(const std::vector< TLorentzVector > &particleMom, const Amg::MatrixX &cov, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:204
DerivationFramework::CascadeTools::massProbability
double massProbability(double V0Mass, double mass, double massErr) const
Definition CascadeTools.cxx:519
DerivationFramework::CascadeTools::tau
double tau(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:195
DerivationFramework::CascadeTools::SetFullMatrix
Amg::MatrixX SetFullMatrix(int NTrk, const std::vector< float > &Matrix) const
Definition CascadeTools.cxx:588
DerivationFramework::CascadeTools::vertexProbability
double vertexProbability(int ndf, double chi2) const
Definition CascadeTools.cxx:537
DerivationFramework::CascadeTools::a0
double a0(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:448
DerivationFramework::CascadeTools::cosTheta_xy
double cosTheta_xy(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:328
DerivationFramework::CascadeTools::cosTheta
double cosTheta(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:319
DerivationFramework::CascadeTools::a0zError
double a0zError(const std::vector< TLorentzVector > &particleMom, const Amg::MatrixX &cov, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:347
DerivationFramework::CascadeTools::a0z
double a0z(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:338
DerivationFramework::CascadeTools::lxyError
double lxyError(const std::vector< TLorentzVector > &particleMom, const Amg::MatrixX &cov, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:152
DerivationFramework::CascadeTools::CascadeTools
CascadeTools(const std::string &t, const std::string &n, const IInterface *p)
Default constructor due to Athena interface.
Definition CascadeTools.cxx:14
DerivationFramework::CascadeTools::pca
Amg::Vector3D pca(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition CascadeTools.cxx:305
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::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
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition EventPrimitives.h:27
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
DerivationFramework
THE reconstruction tool.
Definition ParticleSortingAlg.h:24
DerivationFramework::s_CONST
constexpr double s_CONST
Definition CascadeTools.cxx:22
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