Math.h
Go to the documentation of this file.
1 /*
2  * MoMEMta: a modular implementation of the Matrix Element Method
3  * Copyright (C) 2016 Universite catholique de Louvain (UCL), Belgium
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  */
18 
19 #pragma once
20 
21 #include <cmath>
22 #include <vector>
23 
25 #define SQ(x) ((x) * (x))
26 #define CB(x) ((x) * (x) * (x))
28 #define QU(x) ((x) * (x) * (x) * (x))
30 
41 template <typename T> T sign(const T& x) {
42  if (x > 0)
43  return 1;
44  else if (!x)
45  return 0;
46  else
47  return -1;
48 }
49 
53 template <typename T> inline double dP_over_dE(const T& v) {
54  const double rad = SQ(v.E()) - SQ(v.M());
55  if (rad <= 0)
56  return 0.;
57  else
58  return v.E() / std::sqrt(rad);
59 }
60 
61 // If the NWA is used for a particle, a multiplication factor has to be introduced
62 // because of the integrated-out delta function
63 inline double jacobianNWA(const double mass, const double width) {
64  return (M_PI / 2. + std::atan(mass / width)) * mass * width;
65 }
66 
67 // Compute cos(x +- 2*pi/3) in a more "analytical" way (pm = +- 1)
68 // Useful for solveCubic
69 inline double cosXpm2PI3(const double x, const double pm) {
70  return -0.5 * (std::cos(x) + pm * std::sin(x) * std::sqrt(3.));
71 }
72 
88 bool solveQuadratic(const double a, const double b, const double c, std::vector<double>& roots,
89  bool verbose = false);
90 
108 bool solveCubic(const double a, const double b, const double c, const double d,
109  std::vector<double>& roots, bool verbose = false);
110 
111 // Finds the real solutions to a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
112 // Handles special case a=0.
113 // Appends the solutions to the std::vector roots, making no attempt to check whether the vector is
114 // empty.
115 // Multiple roots appear multiple times.
116 // If verbose is true (default is false), the solutions are printed,
117 // as well as the polynomial evaluated on these solutions
118 //
119 // See https://en.wikipedia.org/wiki/Quartic_function#Solving_by_factoring_into_quadratics
120 // https://fr.wikipedia.org/wiki/Méthode_de_Descartes
121 //
122 // .
123 
147 bool solveQuartic(const double a, const double b, const double c, const double d, const double e,
148  std::vector<double>& roots, bool verbose = false);
149 
175 bool solve2Quads(const double a20, const double a02, const double a11, const double a10,
176  const double a01, const double a00, const double b20, const double b02,
177  const double b11, const double b10, const double b01, const double b00,
178  std::vector<double>& E1, std::vector<double>& E2, bool verbose = false);
179 
197 bool solve2QuadsDeg(const double a11, const double a10, const double a01, const double a00,
198  const double b11, const double b10, const double b01, const double b00,
199  std::vector<double>& E1, std::vector<double>& E2, bool verbose = false);
200 
215 bool solve2Linear(const double a10, const double a01, const double a00, const double b10,
216  const double b01, const double b00, std::vector<double>& E1,
217  std::vector<double>& E2, bool verbose = false);
218 
222 double BreitWigner(const double s, const double m, const double g);
223 
227 bool ApproxComparison(double value, double expected);
bool solveQuartic(const double a, const double b, const double c, const double d, const double e, std::vector< double > &roots, bool verbose=false)
Finds the real solutions to .
Definition: Math.cc:117
bool ApproxComparison(double value, double expected)
Compare two doubles and return true if they are approximatively equal.
Definition: Math.cc:409
bool solveQuadratic(const double a, const double b, const double c, std::vector< double > &roots, bool verbose=false)
Finds the real solutions to .
Definition: Math.cc:26
T sign(const T &x)
Sign function.
Definition: Math.h:41
bool solve2Linear(const double a10, const double a01, const double a00, const double b10, const double b01, const double b00, std::vector< double > &E1, std::vector< double > &E2, bool verbose=false)
Solve a system of two linear equations.
Definition: Math.cc:357
#define SQ(x)
Compute .
Definition: Math.h:25
double dP_over_dE(const T &v)
Compute .
Definition: Math.h:53
double BreitWigner(const double s, const double m, const double g)
A relativist Breit-Wigner distribution.
Definition: Math.cc:400
bool solve2QuadsDeg(const double a11, const double a10, const double a01, const double a00, const double b11, const double b10, const double b01, const double b00, std::vector< double > &E1, std::vector< double > &E2, bool verbose=false)
Solve a system of two degenerated quadratic equations.
bool solve2Quads(const double a20, const double a02, const double a11, const double a10, const double a01, const double a00, const double b20, const double b02, const double b11, const double b10, const double b01, const double b00, std::vector< double > &E1, std::vector< double > &E2, bool verbose=false)
Solve a system of two quadratic equations.
Definition: Math.cc:192
bool solveCubic(const double a, const double b, const double c, const double d, std::vector< double > &roots, bool verbose=false)
Finds the real solutions to .
Definition: Math.cc:69