GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MathUtils.h
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \namespace genie::utils::math
5 
6 \brief Simple mathematical utilities not found in ROOT's TMath
7 
8 \author Costas Andreopoulos <c.andreopoulos \at cern.ch>
9  University of Liverpool
10 
11 \created May 06, 2004
12 
13 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
14  For the full text of the license visit http://copyright.genie-mc.org
15 */
16 //____________________________________________________________________________
17 
18 #ifndef _MATH_UTILS_H_
19 #define _MATH_UTILS_H_
20 
21 #include <vector>
22 #include <TMatrixD.h>
23 #include <TVectorD.h>
24 #include <TLorentzVector.h>
25 #include "Framework/Utils/Range1.h"
26 #include "cmath"
27 
28 using std::vector;
29 
30 namespace genie {
31 namespace utils {
32 
33 namespace math
34 {
35 
36  // This class has been created to perform several operations with long
37  // doubles. It is needed in HEDIS because the kinematics of the outgoing
38  // particles can be so large that the on-shell feature is not fulfilled
39  // many times due to the precission of double.
41 
42  public :
43  LongLorentzVector(double px, double py, double pz, double e) {
44  fPx = (long double) px;
45  fPy = (long double) py;
46  fPz = (long double) pz;
47  fE = (long double) e;
48  }
49  LongLorentzVector(const TLorentzVector & p4) {
50  fPx = (long double) p4.Px();
51  fPy = (long double) p4.Py();
52  fPz = (long double) p4.Pz();
53  fE = (long double) p4.E();
54  }
56 
57  long double Px (void) { return fPx; }
58  long double Py (void) { return fPy; }
59  long double Pz (void) { return fPz; }
60  long double E (void) { return fE; }
61  long double P (void) { return sqrtl(fPx*fPx+fPy*fPy+fPz*fPz); }
62  long double M (void) { return sqrtl(fE*fE-fPx*fPx-fPy*fPy-fPz*fPz); }
63  long double M2 (void) { return fE*fE-fPx*fPx-fPy*fPy-fPz*fPz; }
64  long double Dx (void) { return fPx/sqrtl(fPx*fPx+fPy*fPy+fPz*fPz); }
65  long double Dy (void) { return fPy/sqrtl(fPx*fPx+fPy*fPy+fPz*fPz); }
66  long double Dz (void) { return fPz/sqrtl(fPx*fPx+fPy*fPy+fPz*fPz); }
67 
68  void Rotate (LongLorentzVector axis) {
69  long double up = axis.Dx()*axis.Dx() + axis.Dy()*axis.Dy();
70  if (up) {
71  up = sqrtl(up);
72  long double pxaux = fPx, pyaux = fPy, pzaux = fPz;
73  fPx = (axis.Dx()*axis.Dz()*pxaux - axis.Dy()*pyaux + axis.Dx()*up*pzaux)/up;
74  fPy = (axis.Dy()*axis.Dz()*pxaux + axis.Dx()*pyaux + axis.Dy()*up*pzaux)/up;
75  fPz = (axis.Dz()*axis.Dz()*pxaux - pxaux + axis.Dz()*up*pzaux)/up;
76  }
77  else if (axis.Dz() < 0.) { // phi=0 teta=pi
78  fPx = -fPx;
79  fPz = -fPz;
80  }
81  }
82 
83  void BoostZ (long double bz) {
84  long double b2 = bz*bz;
85  long double gamma = 1.0 / sqrtl(1.0 - b2);
86  long double bp = bz*fPz;
87  long double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;
88  fPz = fPz + gamma2*bp*bz + gamma*bz*fE;
89  fE = gamma*(fE + bp);
90  }
91 
92  void BoostY (long double by) {
93  long double b2 = by*by;
94  long double gamma = 1.0 / sqrtl(1.0 - b2);
95  long double bp = by*fPy;
96  long double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;
97  fPy = fPy + gamma2*bp*by + gamma*by*fE;
98  fE = gamma*(fE + bp);
99  }
100 
101  private :
102 
103  long double fPx;
104  long double fPy;
105  long double fPz;
106  long double fE;
107  };
108 
109  // Cholesky decomposition. Returns lower triangular matrix.
110  TMatrixD CholeskyDecomposition (const TMatrixD& cov);
111  // Generates a vector of correlated parameters.
112  TVectorD CholeskyGenerateCorrelatedParams (const TMatrixD& Lch, TVectorD& mean);
113  TVectorD CholeskyGenerateCorrelatedParams (const TMatrixD& Lch, TVectorD& mean, TVectorD& g_uncorrelated);
114  // Generates a vector of correlated parameter variations.
115  TVectorD CholeskyGenerateCorrelatedParamVariations (const TMatrixD& Lch);
116  TVectorD CholeskyCalculateCorrelatedParamVariations(const TMatrixD& Lch, TVectorD& g_uncorrelated);
117 
118  double KahanSummation (double x[], unsigned int n);
119  double KahanSummation (const vector<double> & x);
120 
121  bool AreEqual (double x1, double x2);
122  bool AreEqual (float x1, float x2);
123 
124  bool IsWithinLimits (double x, Range1D_t range);
125  bool IsWithinLimits (float x, Range1F_t range);
126  bool IsWithinLimits (int i, Range1I_t range);
127 
128  double NonNegative (double x);
129  double NonNegative (float x);
130 
131 } // math namespace
132 } // utils namespace
133 } // genie namespace
134 
135 #endif // _MATH_UTILS_H_
A simple [min,max] interval for integers.
Definition: Range1.h:56
bool AreEqual(double x1, double x2)
Definition: MathUtils.cxx:236
A simple [min,max] interval for doubles.
Definition: Range1.h:42
TVectorD CholeskyGenerateCorrelatedParams(const TMatrixD &Lch, TVectorD &mean)
Definition: MathUtils.cxx:76
bool IsWithinLimits(double x, Range1D_t range)
Definition: MathUtils.cxx:258
A simple [min,max] interval for floats.
Definition: Range1.h:28
TVectorD CholeskyCalculateCorrelatedParamVariations(const TMatrixD &Lch, TVectorD &g_uncorrelated)
Definition: MathUtils.cxx:186
TMatrixD CholeskyDecomposition(const TMatrixD &cov)
Definition: MathUtils.cxx:20
double KahanSummation(double x[], unsigned int n)
Definition: MathUtils.cxx:204
LongLorentzVector(double px, double py, double pz, double e)
Definition: MathUtils.h:43
LongLorentzVector(const TLorentzVector &p4)
Definition: MathUtils.h:49
const double e
TVectorD CholeskyGenerateCorrelatedParamVariations(const TMatrixD &Lch)
Definition: MathUtils.cxx:165
double NonNegative(double x)
Definition: MathUtils.cxx:273
void Rotate(LongLorentzVector axis)
Definition: MathUtils.h:68