EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GaussianIntegralGradHessian.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GaussianIntegralGradHessian.cpp
1 #include "GaussianIntegralGradHessian.h"
2 
3 #include <cmath>
4 #include <vector>
5 
6 using namespace std;
7 using namespace Eigen;
8 
9 namespace FitNewton
10 {
11  GaussianIntegralGradHessian::GaussianIntegralGradHessian() : FunctionGradHessian(3, 2)
12  {
13 
14  }
15 
16 
17  GaussianIntegralGradHessian::~GaussianIntegralGradHessian()
18  {
19 
20  }
21 
22 
23  FunctionGradHessian* GaussianIntegralGradHessian::Clone() const
24  {
25  GaussianIntegralGradHessian* clone = new GaussianIntegralGradHessian();
26  clone->setFixedPar(0, fixedpars[0]);
27  clone->setFixedPar(1, fixedpars[1]);
28  return clone;
29  }
30 
31 
32  bool GaussianIntegralGradHessian::calcValGradHessian(const VectorXd& x, double& val, VectorXd& grad, MatrixXd& hessian)
33  {
34  double temp0 = (fixedpars[0] - x(1));
35  double temp0_2 = temp0*temp0;
36  double temp1 = (fixedpars[1] - x(1));
37  double temp1_2 = temp1*temp1;
38  double temp2 = 1./x(2);
39  double sqrt_pi_o_2 = 0x1.40d931ff627059657ca41fae722cep0;
40  double sqrt_1_o_2 = 0xb.504f333f9de6484597d89b3754abep-4;
41 
42  val = sqrt_pi_o_2*x(0)*x(2)*( erf(sqrt_1_o_2*temp1*temp2) - erf(sqrt_1_o_2*temp0*temp2) );
43 
44  double temp2_2 = temp2*temp2;
45  double etemp1 = exp(-temp1_2*0.5*temp2_2);
46  double etemp0 = exp(-temp0_2*0.5*temp2_2);
47 
48  double temp4 = 1./x(0);
49  double temp5 = ( -etemp1*temp1 + etemp0*temp0 );
50 
51  grad(0) = val*temp4;
52  grad(1) = x(0)*(-etemp1 + etemp0);
53  grad(2) = temp2*( val + x(0)*temp5 );
54 
55  hessian(0,0) = 0.;
56  hessian(0,1) = temp4*grad(1);
57  hessian(0,2) = temp4*grad(2);
58  hessian(1,1) = x(0)*temp2_2*temp5;
59  hessian(1,2) = x(0)*temp2*temp2_2*( -etemp1*temp1_2 + etemp0*temp0_2 );
60  hessian(2,2) = x(0)*temp2_2*temp2_2*( -etemp1*temp1_2*temp1 + -etemp0*temp0_2*temp0 );
61  return true;
62  }
63 }
64 
65 
66 // f = sqrt(pi/2.)*x0*x2*(erf((1/sqrt(2))*(t1-x1)/x2) - erf((1/sqrt(2))*(t0-x1)/x2))
67 //
68 // df/dx0 = f/x0
69 // df/dx1 = x0*[ -exp(-(t1-x1)^2/(2*x2^2)) + exp(-(t0-x1)^2/(2*x2^2)) ]
70 // df/dx2 = (1/x2)*(f + x0*[ -exp(-(t1-x1)^2/(2*x2^2))*(t1-x1) + exp(-(t0-x1)^2/(2*x2^2))*(t0-x1) ])
71 //
72 // d^2f/dx0^2 = 0
73 // d^2f/dx0dx1 = (df/dx1)/x0
74 // d^2f/dx0dx2 = (df/dx2)/x0
75 // d^2f/dx1^2 = (x0/x2^2)*[ -exp(-(t1-x1)^2/(2*x2^2))*(t1-x1) + exp(-(t0-x1)^2/(2*x2^2))*(t0-x1) ]
76 // d^2f/dx1dx2 = (x0/x2^3)*[ -exp(-(t1-x1)^2/(2*x2^2))*(t1-x1)^2 + exp(-(t0-x1)^2/(2*x2^2))*(t0-x1)^2 ]
77 // d^2f/dx2^2 = (x0/x2^4)*[ -exp(-(t1-x1)^2/(2*x2^2))*(t1-x1)^3 + exp(-(t0-x1)^2/(2*x2^2))*(t0-x1)^3 ]
78