Loading [MathJax]/extensions/tex2jax.js
ATLAS Offline Software
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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 
11 namespace DerivationFramework {
12 
13 
14 CascadeTools::CascadeTools(const std::string& t, const std::string& n, const IInterface* p) :
15  AthAlgTool(t,n,p)
16 {
17 }
18 
19 CascadeTools::~CascadeTools() {}
20 
21 //Light speed constant for various calculations
22 constexpr double s_CONST = 1000./299.792;
23 
24 double 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 
33 double 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 
54 double 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 
75 double 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 
113 double 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 
120 double 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 
141 double 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 
152 double 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 
195 double 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 
204 double 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 
251 double 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 
259 double 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 
305 Amg::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 
319 double 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 
328 double 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 
338 double 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 
347 double 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 
395 double 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 
403 double 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 
448 double 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 
456 double 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 
505 Amg::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 
519 double 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 
537 double 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 
555 Amg::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 
588 Amg::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 
Matrix
Definition: Trigger/TrigT1/TrigT1RPChardware/TrigT1RPChardware/Matrix.h:15
covarianceTool.ndf
ndf
Definition: covarianceTool.py:678
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
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
DerivationFramework::CascadeTools::pca
Amg::Vector3D pca(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition: CascadeTools.cxx:305
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::CascadeTools
CascadeTools(const std::string &t, const std::string &n, const IInterface *p)
Default constructor due to Athena interface.
Definition: CascadeTools.cxx:14
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:27
xAOD::Vertex_v1::nTrackParticles
size_t nTrackParticles() const
Get the number of tracks associated with this vertex.
Definition: Vertex_v1.cxx:270
DerivationFramework::CascadeTools::vertexProbability
double vertexProbability(int ndf, double chi2) const
Definition: CascadeTools.cxx:537
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
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
hist_file_dump.d
d
Definition: hist_file_dump.py:143
DMTest::P
P_v1 P
Definition: P.h:23
DerivationFramework::CascadeTools::lxy
double lxy(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition: CascadeTools.cxx:141
skel.it
it
Definition: skel.GENtoEVGEN.py:407
plotBeamSpotVxVal.cov
cov
Definition: plotBeamSpotVxVal.py:201
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::a0xy
double a0xy(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition: CascadeTools.cxx:395
xAOD::Vertex_v1::position
const Amg::Vector3D & position() const
Returns the 3-pos.
DerivationFramework::CascadeTools::pT
double pT(const std::vector< TLorentzVector > &moms) const
Definition: CascadeTools.cxx:113
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
DerivationFramework::CascadeTools::massProbability
double massProbability(double V0Mass, double mass, double massErr) const
Definition: CascadeTools.cxx:519
DerivationFramework::CascadeTools::momentum
Amg::Vector3D momentum(const std::vector< TLorentzVector > &particleMom) const
Definition: CascadeTools.cxx:505
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::convertCovMatrix
Amg::MatrixX * convertCovMatrix(const xAOD::Vertex *vxCandidate) const
Definition: CascadeTools.cxx:555
DerivationFramework::CascadeTools::invariantMassError
double invariantMassError(const std::vector< TLorentzVector > &moms, const Amg::MatrixX &cov, const std::vector< double > &masses) const
Definition: CascadeTools.cxx:75
TRTCalib_cfilter.p2
p2
Definition: TRTCalib_cfilter.py:131
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
DerivationFramework::s_CONST
constexpr double s_CONST
Definition: CascadeTools.cxx:22
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:731
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
DerivationFramework::CascadeTools::a0
double a0(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition: CascadeTools.cxx:448
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:525
DerivationFramework
THE reconstruction tool.
Definition: ParticleSortingAlg.h:24
DerivationFramework::CascadeTools::cosTheta
double cosTheta(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition: CascadeTools.cxx:319
VP1PartSpect::E
@ E
Definition: VP1PartSpectFlags.h:21
DerivationFramework::CascadeTools::invariantMass
double invariantMass(const std::vector< TLorentzVector > &moms, const std::vector< double > &masses) const
Definition: CascadeTools.cxx:33
CascadeTools.h
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
xAOD::Vertex_v1::covariance
const std::vector< float > & covariance() const
Returns the covariance matrix as a simple vector of values.
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
python.testIfMatch.matrix
matrix
Definition: testIfMatch.py:66
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
DerivationFramework::CascadeTools::pTError
double pTError(const std::vector< TLorentzVector > &moms, const Amg::MatrixX &cov) const
Definition: CascadeTools.cxx:120
DerivationFramework::CascadeTools::a0z
double a0z(const std::vector< TLorentzVector > &particleMom, const xAOD::Vertex *SV, const xAOD::Vertex *PV) const
Definition: CascadeTools.cxx:338
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
LArCellConditions.sv
bool sv
Definition: LArCellConditions.py:45
python.changerun.pv
pv
Definition: changerun.py:81
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::SetFullMatrix
Amg::MatrixX SetFullMatrix(int NTrk, const std::vector< float > &Matrix) const
Definition: CascadeTools.cxx:588
AthAlgTool
Definition: AthAlgTool.h:26
H5Utils::internal::PT
H5::PredType PT
Definition: H5Traits.cxx:15
DerivationFramework::CascadeTools::~CascadeTools
~CascadeTools()
Virtual destructor.
Definition: CascadeTools.cxx:19
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
SUSY_SimplifiedModel_PreInclude.masses
dictionary masses
Definition: SUSY_SimplifiedModel_PreInclude.py:7
python.SystemOfUnits.L
float L
Definition: SystemOfUnits.py:91
Generated on Fri May 9 2025 21:07:35 for ATLAS Offline Software by   doxygen 1.8.18