CMS 3D CMS Logo

ResolutionHelper.cc
Go to the documentation of this file.
3 #include <cmath>
4 #include <iostream>
5 
6 using namespace std; // faster sqrt, sqrt, pow(x,2)....
7 
9  AlgebraicSymMatrix44 &covariance,
10  const math::XYZTLorentzVector &initialP4) {
11  double inv;
12  switch (parametrization) {
15  // for us parameter[3] = mass, for KinFitter parameter[3] = mass/initialP4.mass();
16  inv = 1.0 / initialP4.mass();
17  for (int i = 0; i < 4; i++) {
18  covariance(3, i) *= inv;
19  }
20  covariance(3, 3) *= inv;
21  break;
23  // for us parameter[3] = energy, for KinFitter parameter[3] = energy/initialP4.energy();
24  inv = 1.0 / initialP4.energy();
25  for (int i = 0; i < 4; i++) {
26  covariance(3, i) *= inv;
27  }
28  covariance(3, 3) *= inv;
29  break;
30  default:; // nothing to do
31  }
32 }
33 
35  const AlgebraicSymMatrix44 &covariance,
37  switch (parametrization) {
38  // ==== CARTESIAN ====
42  // p2 = px^2 + py^2 + pz^2
43  // dp/dp_i = p_i/p ==> it is a unit vector!
44  AlgebraicVector3 derivs(p4.X(), p4.Y(), p4.Z());
45  derivs.Unit();
46  return sqrt(ROOT::Math::Similarity(derivs, covariance.Sub<AlgebraicSymMatrix33>(0, 0)));
47  }
48  // ==== SPHERICAL (P) ====
52  return sqrt(covariance(0, 0));
53  // ==== SPHERICAL (1/P) ====
55  return sqrt(covariance(0, 0)) * (p4.P2());
56  // ==== HEP CYLINDRICAL (with Pt = Et, P = E) ====
59  return getResolE(parametrization, covariance, p4);
60  // ==== OTHER ====
62  throw cms::Exception("Invalid parametrization") << parametrization;
63  default:
64  throw cms::Exception("Not Implemented")
65  << "getResolP not yet implemented for parametrization " << parametrization;
66  }
67 }
69  const AlgebraicSymMatrix44 &covariance,
71  switch (parametrization) {
72  // ==== CARTESIAN ====
76  double pti2 = 1.0 / (p4.Perp2());
77  return sqrt((covariance(0, 0) * p4.Px() * p4.Px() + covariance(1, 1) * p4.Py() * p4.Py() +
78  2 * covariance(0, 1) * p4.Px() * p4.Py()) *
79  pti2);
80  }
81  // ==== SPHERICAL (P, Theta) ====
85  // pt = p * sin(theta)
86  double a = sin(p4.Theta());
87  double b = p4.P() * cos(p4.Theta());
88  return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
89  }
90  // ==== SPHERICAL (1/P) ====
92  // pt = (1/pi) * sin(theta)
93  double p = p4.P();
94  double a = -(p * p) * sin(p4.Theta());
95  double b = p * cos(p4.Theta());
96  return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
97  }
98  // ==== HEP CYLINDRICAL (with Pt = Et) ====
101  return sqrt(covariance(0, 0));
103  throw cms::Exception("Invalid parametrization") << parametrization;
104  default:
105  throw cms::Exception("Not Implemented")
106  << "getResolPt not yet implemented for parametrization " << parametrization;
107  }
108 }
109 
111  const AlgebraicSymMatrix44 &covariance,
113  switch (parametrization) {
114  // ==== SPHERICAL (P) ====
118  return 1.0 / p4.P2() * sqrt(covariance(0, 0));
119  // ==== SPHERICAL (1/P) ====
121  return sqrt(covariance(0, 0));
122  // ==== OTHER ====
128  return 1.0 / p4.P2() * getResolP(parametrization, covariance, p4);
130  throw cms::Exception("Invalid parametrization") << parametrization;
131  default:
132  throw cms::Exception("Not Implemented")
133  << "getResolPInv not yet implemented for parametrization " << parametrization;
134  }
135 }
136 
138  const AlgebraicSymMatrix44 &covariance,
140  switch (parametrization) {
141  // ==== CARTESIAN ====
145  return sqrt(covariance(0, 0));
146  // ==== SPHERICAL (P) ====
151  // Px = P * sin(theta) * cos(phi)
152  double p = p4.P();
153  AlgebraicVector3 derivs;
154  derivs[0] = sin(p4.Theta()) * cos(p4.Phi()); // now let's hope gcc does common subexpr optimiz.
156  derivs[0] *= -(p * p);
157  }
158  derivs[1] = p * cos(p4.Theta()) * cos(p4.Phi());
159  derivs[2] = p * sin(p4.Theta()) * -sin(p4.Phi());
160  return sqrt(ROOT::Math::Similarity(derivs, covariance.Sub<AlgebraicSymMatrix33>(0, 0)));
161  }
162  // ==== HEP CYLINDRICAL (with Pt = Et) ====
165  // Px = Pt * cos(phi)
166  double a = cos(p4.Phi());
167  double b = -p4.Pt() * sin(p4.Phi());
168  return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(2, 0) + b * b * covariance(2, 2));
169  }
170  // ==== OTHERS ====
172  throw cms::Exception("Invalid parametrization") << parametrization;
173  default:
174  throw cms::Exception("Not Implemented")
175  << "getResolPx not yet implemented for parametrization " << parametrization;
176  }
177 }
179  const AlgebraicSymMatrix44 &covariance,
181  switch (parametrization) {
182  // ==== CARTESIAN ====
186  return sqrt(covariance(1, 1));
187  // ==== SPHERICAL (P) ====
192  // Py = P * sin(theta) * sin(phi)
193  double p = p4.P();
194  AlgebraicVector3 derivs;
195  derivs[0] = sin(p4.Theta()) * sin(p4.Phi()); // now let's hope gcc does common subexpr optimiz.
197  derivs[0] *= -(p * p);
198  }
199  derivs[1] = p * cos(p4.Theta()) * sin(p4.Phi());
200  derivs[2] = p * sin(p4.Theta()) * cos(p4.Phi());
201  return sqrt(ROOT::Math::Similarity(derivs, covariance.Sub<AlgebraicSymMatrix33>(0, 0)));
202  }
203  // ==== HEP CYLINDRICAL (with Pt = Et) ====
206  // Py = Pt * sin(phi)
207  double a = sin(p4.Phi());
208  double b = p4.Pt() * cos(p4.Phi());
209  return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(2, 0) + b * b * covariance(2, 2));
210  }
211  // ==== OTHERS ====
213  throw cms::Exception("Invalid parametrization") << parametrization;
214  default:
215  throw cms::Exception("Not Implemented")
216  << "getResolPy not yet implemented for parametrization " << parametrization;
217  }
218 }
220  const AlgebraicSymMatrix44 &covariance,
222  switch (parametrization) {
223  // ==== CARTESIAN ====
227  return sqrt(covariance(2, 2));
228  // ==== SPHERICAL (P) ====
232  // Pz = P * cos(theta)
233  double a = cos(p4.Theta());
234  double b = -p4.P() * sin(p4.Theta());
235  return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(1, 0) + b * b * covariance(1, 1));
236  }
238  // Pz = P * cos(theta)
239  double p = p4.P();
240  double a = -p * p * cos(p4.Theta());
241  double b = -p * sin(p4.Theta());
242  return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(1, 0) + b * b * covariance(1, 1));
243  }
244  // ==== HEP CYLINDRICAL (with Pt = Et) ====
246  // Pz = Pt * ctg(theta) d ctg(x) = -1/sin^2(x)
247  double s = sin(p4.Theta()), c = cos(p4.Theta());
248  double a = c / s;
249  double b = -p4.Pt() / (s * s);
250  return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(1, 0) + b * b * covariance(1, 1));
251  }
253  // Pz = Pt * sinh(eta)
254  double a = sinh(p4.Eta());
255  double b = p4.Et() * cosh(p4.Eta());
256  return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(1, 0) + b * b * covariance(1, 1));
257  }
258  // ==== OTHERS ====
260  throw cms::Exception("Invalid parametrization") << parametrization;
261  default:
262  throw cms::Exception("Not Implemented")
263  << "getResolPz not yet implemented for parametrization " << parametrization;
264  }
265 }
266 
268  const AlgebraicSymMatrix44 &covariance,
270  switch (parametrization) {
271  // ======= ENERGY BASED ==========
274  return sqrt(covariance(3, 3));
275  // ======= ET BASED ==========
277  // E = Et/Sin(theta)
278  double a = 1.0 / sin(p4.Theta()); // dE/dEt
279  double b = -a * a * p4.Et() * cos(p4.Theta()); // dE/dTh
280  return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
281  }
283  // E = Et/Sin(Theta(eta))
284  double th = p4.Theta();
285  double a = 1.0 / sin(th); // dE/dEt
286  double b = a * p4.Et() * cos(th); // dE/dEta: dTh/dEta = - 1.0/sin(theta) = - dE/dEt
287  return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
288  }
289  // ======= MASS BASED ==========
291  AlgebraicVector4 xoE(p4.X(), p4.Y(), p4.Z(), p4.M());
292  xoE *= 1 / p4.E();
293  return sqrt(ROOT::Math::Similarity(xoE, covariance));
294  }
296  AlgebraicVector4 xoE(p4.X(), p4.Y(), p4.Z(), 0);
297  xoE *= 1 / p4.E();
298  return sqrt(ROOT::Math::Similarity(xoE, covariance));
299  }
301  // E = sqrt(P^2 + m^2)
302  double einv = 1.0 / p4.E();
303  double a = p4.P() * einv; // dE/dP
304  double b = p4.M() * einv; // dE/dm
305  return sqrt(a * a * covariance(0, 0) + b * b * covariance(3, 3) + 2 * a * b * covariance(0, 3));
306  }
308  // E = sqrt(P^2 + m^2); |dE/dP| = |P/E| = P/E
309  return p4.P() / p4.E() * sqrt(covariance(0, 0));
310  }
312  {
313  // E = sqrt(P^2 + m^2); |dE/d(1/P)| = P^2 |dE/dP| = P^3/E
314  double p = p4.P();
315  return p * p * p / p4.E() * sqrt(covariance(0, 0));
316  }
317  // ======= OTHER ==========
319  throw cms::Exception("Invalid parametrization") << parametrization;
320  default:
321  throw cms::Exception("Not Implemented")
322  << "getResolE not yet implemented for parametrization " << parametrization;
323  }
324 }
325 
327  const AlgebraicSymMatrix44 &covariance,
329  switch (parametrization) {
330  // ======= ENERGY BASED ==========
332  // Et^2 = E^2 * (Pt^2/P^2)
333  double pt2 = p4.Perp2();
334  double pz2 = ROOT::Math::Square(p4.Pz()), p2 = pt2 + pz2;
335  double e2OverP4 = ROOT::Math::Square(p4.E() / p2);
336  AlgebraicVector4 derivs(p4.Px(), p4.Py(), p4.Pz(), p4.E());
337  derivs *= (1.0 / p4.Et());
338  derivs[0] *= pz2 * e2OverP4;
339  derivs[1] *= pz2 * e2OverP4;
340  derivs[2] *= -pt2 * e2OverP4;
341  derivs[3] *= pt2 / p2;
342  return sqrt(ROOT::Math::Similarity(derivs, covariance));
343  }
345  // Et = E * Sin(Theta)
346  double st = sin(p4.Theta()), ct = cos(p4.Theta());
347  return sqrt(st * st * covariance(3, 3) + ROOT::Math::Square(ct * p4.E()) * covariance(1, 1) +
348  2 * st * ct * p4.E() * covariance(1, 3));
349  }
350  // ======= ET BASED ==========
353  return sqrt(covariance(0, 0));
354  // ======= MASS BASED ==========
357  // Et^2 = E^2 Sin^2(th) = (p^2 + m^2) * (pt^2) / p^2
358  double pt2 = p4.Perp2();
359  double p2 = pt2 + ROOT::Math::Square(p4.Pz());
360  double e2 = p2 + p4.M2();
361  double s2 = pt2 / p2, pi2 = 1.0 / p2;
362  double et = sqrt(e2 * s2);
363  AlgebraicVector4 derivs(p4.Px(), p4.Py(), p4.Pz(), p4.M());
364  derivs *= 1.0 / et;
365  derivs[0] *= (s2 + e2 * pi2 * (1.0 - pt2 * pi2));
366  derivs[1] *= (s2 + e2 * pi2 * (1.0 - pt2 * pi2));
367  derivs[2] *= (s2 - e2 * pt2 * pi2 * pi2);
369  derivs[3] *= s2;
370  return sqrt(ROOT::Math::Similarity(derivs, covariance));
371  } else {
372  derivs[3] = 0;
373  return sqrt(ROOT::Math::Similarity(derivs, covariance)); // test if Sub<33> is faster
374  }
375  }
377  // Et = E sin(theta); dE/dM = M/E, dE/dP = P/E
378  double s = sin(p4.Theta()), c = cos(p4.Theta());
379  double e = p4.E();
380  AlgebraicVector4 derivs(p4.P() / e * s, e * c, 0, p4.M() / e * s);
381  return sqrt(ROOT::Math::Similarity(derivs, covariance));
382  }
384  // Et = E sin(theta); dE/dP = P/E
385  double s = sin(p4.Theta()), c = cos(p4.Theta());
386  double e = p4.E();
387  double a = p4.P() * s / e;
388  double b = e * c;
389  return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
390  }
392  // Et = E sin(theta); dE/dP = P/E -> dE/d(1/P) = - P^2 dE/dP = - P^3 / E
393  double s = sin(p4.Theta()), c = cos(p4.Theta());
394  double p = p4.P(), e = p4.E();
395  double a = (-p * p * p / e) * s;
396  double b = e * c;
397  return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
398  }
399  // ======= OTHER ==========
401  throw cms::Exception("Invalid parametrization") << parametrization;
402  default:
403  throw cms::Exception("Not Implemented")
404  << "getResolEt not yet implemented for parametrization " << parametrization;
405  }
406 }
407 
409  const AlgebraicSymMatrix44 &covariance,
411  switch (parametrization) {
412  // ====== MASS CONSTRAINED =====
418  return 0;
419  // ======= MASS BASED ==========
422  return sqrt(covariance(3, 3));
423  // ======= ENERGY BASED ==========
424  case pat::CandKinResolution::ESpher: { // M^2 = E^2 - P^2
425  double dMdE = p4.E() / p4.M(), dMdP = -p4.P() / p4.M();
426  return sqrt(dMdP * dMdP * covariance(0, 0) + 2 * dMdP * dMdE * covariance(0, 3) + dMdE * dMdE * covariance(3, 3));
427  }
428  case pat::CandKinResolution::ECart: { // M^2 = E^2 - sum_i P_i^2
429  AlgebraicVector4 derivs(-p4.Px(), -p4.Py(), -p4.Pz(), p4.E());
430  derivs *= 1.0 / p4.M();
431  return sqrt(ROOT::Math::Similarity(derivs, covariance));
432  }
433  throw cms::Exception("Not Implemented")
434  << "getResolM not yet implemented for parametrization " << parametrization;
436  throw cms::Exception("Invalid parametrization") << parametrization;
437  default:
438  throw cms::Exception("Not Implemented")
439  << "getResolM not yet implemented for parametrization " << parametrization;
440  }
441 }
442 
443 inline double DetaDtheta(double theta) {
444  // y = -ln(tg(x/2)) =>
445  // y' = - 1/tg(x/2) * 1/(cos(x/2))^2 * 1/2 = - 1 / (2 * sin(x/2) * cos(x/2)) = -1/sin(x)
446  return -1.0 / sin(theta);
447 }
448 inline double DthetaDeta(double eta) {
449  // y = 2 atan(exp(-x))
450  // y' = 2 * 1/(1+exp^2) * exp(-x) * (-1) = - 2 * exp/(1+exp^2) = - 2 / (exp + 1/exp)
451  double e = exp(-eta);
452  return -2.0 / (e + 1.0 / e);
453 }
454 
456  const AlgebraicSymMatrix44 &covariance,
458  switch (parametrization) {
462  // dEta = dTheta * dEta/dTheta
463  return abs(DetaDtheta(p4.Theta())) * getResolTheta(parametrization, covariance, p4);
464  case pat::CandKinResolution::ESpher: // all the ones which have
465  case pat::CandKinResolution::MCPInvSpher: // theta as parameter 1
469  return sqrt(covariance(1, 1)) * abs(DetaDtheta(p4.Theta()));
470  case pat::CandKinResolution::EtEtaPhi: // as simple as that
471  return sqrt(covariance(1, 1));
473  throw cms::Exception("Invalid parametrization") << parametrization;
474  default:
475  throw cms::Exception("Not Implemented")
476  << "getResolEta not yet implemented for parametrization " << parametrization;
477  }
478 }
480  const AlgebraicSymMatrix44 &covariance,
482  switch (parametrization) {
486  // theta = acos( pz / p ) ; d acos(x) = - 1 / sqrt( 1 - x*x) dx = - p/pt dx
487  double pt2 = p4.Perp2();
488  double p = p4.P(), pi = 1.0 / p, pi3 = pi * pi * pi;
489  double dacos = -p / sqrt(pt2);
490  AlgebraicVector3 derivs;
491  derivs[0] = -p4.Px() * p4.Pz() * dacos * pi3;
492  derivs[1] = -p4.Py() * p4.Pz() * dacos * pi3;
493  derivs[2] = pt2 * dacos * pi3;
494  return sqrt(ROOT::Math::Similarity(derivs, covariance.Sub<AlgebraicSymMatrix33>(0, 0)));
495  }
496  case pat::CandKinResolution::ESpher: // all the ones which have
497  case pat::CandKinResolution::MCPInvSpher: // theta as parameter 1
501  return sqrt(covariance(1, 1));
503  return sqrt(covariance(1, 1)) * abs(DthetaDeta(p4.Eta()));
505  throw cms::Exception("Invalid parametrization") << parametrization;
506  default:
507  throw cms::Exception("Not Implemented")
508  << "getResolTheta not yet implemented for parametrization " << parametrization;
509  }
510 }
512  const AlgebraicSymMatrix44 &covariance,
514  double pt2 = p4.Perp2();
515  switch (parametrization) {
519  return sqrt(ROOT::Math::Square(p4.Px()) * covariance(1, 1) + ROOT::Math::Square(p4.Py()) * covariance(0, 0) +
520  -2 * p4.Px() * p4.Py() * covariance(0, 1)) /
521  pt2;
522  case pat::CandKinResolution::ESpher: // all the ones which have
523  case pat::CandKinResolution::MCPInvSpher: // phi as parameter 2
528  return sqrt(covariance(2, 2));
530  throw cms::Exception("Invalid parametrization") << parametrization;
531  default:
532  throw cms::Exception("Not Implemented")
533  << "getResolPhi not yet implemented for parametrization " << parametrization;
534  }
535 }
double getResolEt(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
double DetaDtheta(double theta)
double getResolPz(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
double getResolPy(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
double getResolP(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
double getResolPx(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
double getResolPt(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
double getResolPhi(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
double getResolM(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
void rescaleForKinFitter(const pat::CandKinResolution::Parametrization parametrization, AlgebraicSymMatrix44 &covariance, const math::XYZTLorentzVector &initialP4)
const double pi2
Definition: Thrust.cc:4
double getResolPInv(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
double getResolE(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
const Double_t pi
parametrization
specify parametrization (see SWGuidePATKinematicResolutions for more details)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
T sqrt(T t)
Definition: SSEVec.h:19
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double DthetaDeta(double eta)
ROOT::Math::SMatrix< double, 4, 4, ROOT::Math::MatRepSym< double, 4 > > AlgebraicSymMatrix44
ROOT::Math::SVector< double, 4 > AlgebraicVector4
double b
Definition: hdecay.h:118
math::XYZTLorentzVector LorentzVector
double a
Definition: hdecay.h:119
ROOT::Math::SVector< double, 3 > AlgebraicVector3
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
double getResolEta(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
double getResolTheta(pat::CandKinResolution::Parametrization parametrization, const AlgebraicSymMatrix44 &covariance, const pat::CandKinResolution::LorentzVector &p4)
Geom::Theta< T > theta() const