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