PYTHIA  8.312
MathTools.h
1 // MathTools.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Header file for some mathematics tools, like special functions.
7 #ifndef Pythia8_MathTools_H
8 #define Pythia8_MathTools_H
9 
10 // Header file for the MathTools methods.
11 #include "Pythia8/Basics.h"
12 #include "Pythia8/PythiaStdlib.h"
13 
14 namespace Pythia8 {
15 
16 //==========================================================================
17 
18 // The Gamma function for real argument.
19 double gammaReal(double x);
20 
21 // Modified Bessel functions of the first and second kinds.
22 double besselI0(double x);
23 double besselI1(double x);
24 double besselK0(double x);
25 double besselK1(double x);
26 
27 // Integrate f(x) dx over the specified range.
28 bool integrateGauss(double& resultOut, function<double(double)> f,
29  double xLo, double xHi, double tol=1e-6);
30 
31 // Solve f(x) = target for x in the specified range.
32 bool brent(double& solutionOut, function<double(double)> f,
33  double target, double xLo, double xHi, double tol=1e-6, int maxIter = 10000);
34 
35 // Gram determinant, invariants used in the argument = 2*pi*pj.
36 double gramDet(double s01tilde, double s12tilde, double s02tilde,
37  double m0, double m1, double m2);
38 double gramDet(Vec4 p0, Vec4 p1, Vec4 p2);
39 
40 // Dilogarithm.
41 double Li2 (const double, const double kmax = 100.0, const double xerr = 1e-9);
42 
43 // Standard factorial.
44 double factorial(const int);
45 
46 // Binomial coefficient.
47 int binomial (const int,int);
48 
49 // Lambert W function.
50 double lambertW (const double x);
51 
52 // Kallen function.
53 double kallenFunction(const double x, const double y, const double z);
54 
55 // Generate linearly or logarithmically spaced points.
56 vector<double> linSpace(int nPts, double xMin, double xMax);
57 vector<double> logSpace(int nPts, double xMin, double xMax);
58 
59 //==========================================================================
60 
61 // LinearInterpolator class.
62 // Used to interpolate between values in linearly spaced data.
63 // The interpolation uses the arithmetic mean.
64 
66 
67 public:
68 
69  LinearInterpolator() = default;
70 
71  // Constructor.
72  LinearInterpolator(double leftIn, double rightIn, vector<double> ysIn)
73  : leftSave(leftIn), rightSave(rightIn), ysSave(ysIn) { }
74 
75  // Function to get y-values of interpolation data.
76  const vector<double>& data() const { return ysSave; }
77 
78  // x-values are linearly spaced on the interpolation region.
79  double left() const { return leftSave; }
80  double right() const { return rightSave; }
81  double dx() const { return (rightSave - leftSave) / (ysSave.size() - 1); }
82  double xi(int i) const { return leftSave + i * dx(); }
83 
84  // Get interpolated value at the specified point.
85  double at(double x) const;
86  double operator()(double x) const { return at(x); }
87 
88  // Plot the data points of this LinearInterpolator in a histogram.
89  Hist plot(string title) const;
90  Hist plot(string title, double xMin, double xMax) const;
91 
92  // Sample a random number distributed as given by this LinearInterpolator.
93  double sample(Rndm& rndm) const;
94 
95 private:
96 
97  // Data members
98  double leftSave, rightSave;
99  vector<double> ysSave;
100 
101 };
102 
103 //==========================================================================
104 
105 // LogInterpolator class.
106 // Used to interpolate between values in logarithmically spaced data.
107 // The interpolation uses the geometric mean.
108 
110 
111 public:
112 
113  // Default constructor.
114  LogInterpolator() = default;
115 
116  // Constructor.
117  LogInterpolator(double leftIn, double rightIn, vector<double> ysIn)
118  : leftSave(leftIn), rightSave(rightIn), ysSave(ysIn) {
119  if (ysIn.size() <= 1)
120  rxSave = numeric_limits<double>::quiet_NaN();
121  else
122  rxSave = pow(rightSave / leftSave, 1. / (ysSave.size() - 1));
123  }
124 
125  // Function to get y-values of interpolation data.
126  const vector<double>& data() const { return ysSave; }
127 
128  // x-values are logarithmically spaced on the interpolation region.
129  double left() const { return leftSave; }
130  double right() const { return rightSave; }
131  double rx() const { return rxSave; }
132 
133  // Get interpolated value at the specified point.
134  double at(double x) const;
135  double operator()(double x) const { return at(x); }
136 
137  // Plot the data points of this LogInterpolator in a histogram.
138  Hist plot(string title, int nBins, double xMin, double xMax) const;
139 
140 private:
141 
142  // Data members.
143  double leftSave, rightSave, rxSave;
144  vector<double> ysSave;
145 
146 };
147 
148 //==========================================================================
149 
150 // Class for the "Hungarian" pairing algorithm. Adapted for Vincia
151 // from an implementation by M. Buehren and C. Ma, see notices below.
152 
153 // This is a C++ wrapper with slight modification of a hungarian algorithm
154 // implementation by Markus Buehren. The original implementation is a few
155 // mex-functions for use in MATLAB, found here:
156 // http://www.mathworks.com/matlabcentral/fileexchange/
157 // 6543-functions-for-the-rectangular-assignment-problem
158 //
159 // Both this code and the orignal code are published under the BSD license.
160 // by Cong Ma, 2016.
161 //
162 // Copyright (c) 2014, Markus Buehren
163 // All rights reserved.
164 //
165 // Redistribution and use in source and binary forms, with or without
166 // modification, are permitted provided that the following conditions are
167 // met:
168 //
169 // * Redistributions of source code must retain the above copyright
170 // notice, this list of conditions and the following disclaimer.
171 // * Redistributions in binary form must reproduce the above copyright
172 // notice, this list of conditions and the following disclaimer in
173 // the documentation and/or other materials provided with the distribution
174 //
175 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
176 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
177 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
178 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
179 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
180 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
181 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
182 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
183 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
184 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
185 // POSSIBILITY OF SUCH DAMAGE.
186 
188 
189 public:
190 
191  // Function wrapper for solving assignment problem.
192  double solve(vector <vector<double> >& distMatrix, vector<int>& assignment);
193 
194  private:
195 
196  // Solve optimal solution for assignment problem using Munkres algorithm,
197  // also known as the Hungarian algorithm.
198  void optimal(vector<int>& assignment, double& cost,
199  vector<double>& distMatrix, int nOfRows, int nOfColumns);
200  // Build the assignment vector.
201  void vect(vector<int>& assignment, vector<bool>& starMatrix, int nOfRows,
202  int nOfColumns);
203  // Calculate the assignment cost.
204  void calcCost(vector<int>& assignment, double& cost,
205  vector<double>& distMatrix, int nOfRows);
206  // Factorized step 2a of the algorithm.
207  void step2a(vector<int>& assignment, vector<double>& distMatrix,
208  vector<bool>& starMatrix, vector<bool>& newStarMatrix,
209  vector<bool>& primeMatrix, vector<bool>& coveredColumns,
210  vector<bool>& coveredRows, int nOfRows, int nOfColumns, int minDim);
211  // Factorized step 2b of the algorithm.
212  void step2b(vector<int>& assignment, vector<double>& distMatrix,
213  vector<bool>& starMatrix, vector<bool>& newStarMatrix,
214  vector<bool>& primeMatrix, vector<bool>& coveredColumns,
215  vector<bool>& coveredRows, int nOfRows, int nOfColumns, int minDim);
216  // Factorized step 3 of the algorithm.
217  void step3(vector<int>& assignment, vector<double>& distMatrix,
218  vector<bool>& starMatrix, vector<bool>& newStarMatrix,
219  vector<bool>& primeMatrix, vector<bool>& coveredColumns,
220  vector<bool>& coveredRows, int nOfRows, int nOfColumns, int minDim);
221  // Factorized step 4 of the algorithm.
222  void step4(vector<int>& assignment, vector<double>& distMatrix,
223  vector<bool>& starMatrix, vector<bool>& newStarMatrix,
224  vector<bool>& primeMatrix, vector<bool>& coveredColumns,
225  vector<bool>& coveredRows, int nOfRows, int nOfColumns, int minDim,
226  int row, int col);
227  // Factorized step 5 of the algorithm.
228  void step5(vector<int>& assignment, vector<double>& distMatrix,
229  vector<bool>& starMatrix, vector<bool>& newStarMatrix,
230  vector<bool>& primeMatrix, vector<bool>& coveredColumns,
231  vector<bool>& coveredRows, int nOfRows, int nOfColumns, int minDim);
232 
233 };
234 
235 //==========================================================================
236 
237 } // end namespace Pythia8
238 
239 #endif // Pythia8_MathTools_H
double at(double x) const
Get interpolated value at the specified point.
Definition: MathTools.cc:463
double kallenFunction(const double x, const double y, const double z)
Kallen function.
Definition: MathTools.cc:426
double gammaReal(double x)
The Gamma function for real argument.
Definition: MathTools.cc:22
double besselK0(double x)
Definition: MathTools.cc:104
Definition: MathTools.h:109
double left() const
x-values are linearly spaced on the interpolation region.
Definition: MathTools.h:79
int binomial(const int, int)
Binomial coefficient.
Definition: MathTools.cc:390
Definition: MathTools.h:187
Definition: Basics.h:386
double lambertW(const double x)
Lambert W function.
Definition: MathTools.cc:409
LogInterpolator(double leftIn, double rightIn, vector< double > ysIn)
Constructor.
Definition: MathTools.h:117
double left() const
x-values are logarithmically spaced on the interpolation region.
Definition: MathTools.h:129
const vector< double > & data() const
Function to get y-values of interpolation data.
Definition: MathTools.h:126
Hist plot(string title) const
Plot the data points of this LinearInterpolator in a histogram.
Definition: MathTools.cc:490
bool brent(double &solutionOut, function< double(double)> f, double target, double xLo, double xHi, double tol=1e-6, int maxIter=10000)
Solve f(x) = target for x in the specified range.
Definition: MathTools.cc:249
Definition: Basics.h:474
double besselK1(double x)
Definition: MathTools.cc:132
double besselI0(double x)
Modified Bessel functions of the first and second kinds.
Definition: MathTools.cc:44
LinearInterpolator(double leftIn, double rightIn, vector< double > ysIn)
Constructor.
Definition: MathTools.h:72
double sample(Rndm &rndm) const
Sample a random number distributed as given by this LinearInterpolator.
Definition: MathTools.cc:513
double m2(const Vec4 &v1)
The squared invariant mass of one or more four-vectors.
Definition: Basics.cc:605
double Li2(const double, const double kmax=100.0, const double xerr=1e-9)
Dilogarithm.
Definition: MathTools.cc:344
vector< double > linSpace(int nPts, double xMin, double xMax)
Generate linearly or logarithmically spaced points.
Definition: MathTools.cc:434
vector< double > logSpace(int nPts, double xMin, double xMax)
Generate logarithmically spaced points.
Definition: MathTools.cc:446
double factorial(const int)
Standard factorial.
Definition: MathTools.cc:380
double gramDet(double s01tilde, double s12tilde, double s02tilde, double m0, double m1, double m2)
Gram determinant, invariants used in the argument = 2*pi*pj.
Definition: MathTools.cc:328
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
const vector< double > & data() const
Function to get y-values of interpolation data.
Definition: MathTools.h:76
double besselI1(double x)
Definition: MathTools.cc:74
Definition: MathTools.h:65
bool integrateGauss(double &resultOut, function< double(double)> f, double xLo, double xHi, double tol=1e-6)
Integrate f(x) dx over the specified range.
Definition: MathTools.cc:163