You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
574 lines
25 KiB
574 lines
25 KiB
//Renata Kopecna
|
|
|
|
#include "integrals.hh"
|
|
|
|
#include <Math/Vector4D.h>
|
|
#include <Math/Boost.h>
|
|
#include <Math/GSLIntegrator.h>
|
|
|
|
#include <spdlog.h>
|
|
#include <funcs.hh>
|
|
|
|
#include <TMath.h>
|
|
|
|
double integral_x_to_n(double x, int n){
|
|
return pow(x,n+1)/(n+1.0);
|
|
}
|
|
double integrate_x_to_n(double a, double b, int n){
|
|
return integral_x_to_n(b,n) - integral_x_to_n(a,n);
|
|
}
|
|
|
|
//this methods are helpful for the phi integration
|
|
//calculates int x^n * sin(x) dx
|
|
double integral_x_to_n_times_sin_x(double x, int n){
|
|
assert(n >= 0);
|
|
if (n == 0) return -cos(x);
|
|
else return -pow(x, n)*cos(x) + n*integral_x_to_n_times_cos_x(x, n-1);
|
|
}
|
|
double integrate_x_to_n_times_sin_x(double a, double b, int n){
|
|
return integral_x_to_n_times_sin_x(b, n) - integral_x_to_n_times_sin_x(a, n);
|
|
}
|
|
|
|
//calculates int x^n * cos(x) dx (recursive using per partes)
|
|
double integral_x_to_n_times_cos_x(double x, int n){
|
|
assert(n >= 0);
|
|
if (n == 0) return sin(x);
|
|
else return pow(x,n)*sin(x) - n*integral_x_to_n_times_sin_x(x, n-1);
|
|
}
|
|
double integrate_x_to_n_times_cos_x(double a, double b, int n){
|
|
return integral_x_to_n_times_cos_x(b, n)-integral_x_to_n_times_cos_x(a, n);
|
|
}
|
|
|
|
//calculates int x^n * sin(2x) dx (recursive using per partes)
|
|
double integral_x_to_n_times_sin_2x(double x, int n){
|
|
assert(n >= 0);
|
|
if (n == 0) return -0.5*cos(2.0*x);
|
|
else return -pow(x,n)*0.5*cos(2.0*x)
|
|
+0.5*n*integral_x_to_n_times_cos_2x(x,n-1);
|
|
}
|
|
double integrate_x_to_n_times_sin_2x(double a, double b, int n){
|
|
return integral_x_to_n_times_sin_2x(b, n) - integral_x_to_n_times_sin_2x(a, n);
|
|
}
|
|
|
|
//calculates int x^n * cos(2x) dx (recursive using per partes)
|
|
double integral_x_to_n_times_cos_2x(double x, int n){
|
|
assert(n >= 0);
|
|
if (n == 0) return 0.5*sin(2.0*x);
|
|
else return +0.5*pow(x,n)*sin(2.0*x)
|
|
-0.5*n*integral_x_to_n_times_sin_2x(x,n-1);
|
|
}
|
|
double integrate_x_to_n_times_cos_2x(double a, double b, int n){
|
|
return integral_x_to_n_times_cos_2x(b, n) - integral_x_to_n_times_cos_2x(a, n);
|
|
}
|
|
|
|
//calculates int x^n * cos(x)^2 dx
|
|
double integral_x_to_n_times_cos_x_2(double x, int n){
|
|
assert(n >= 0);
|
|
return +1.0/(1.0 + n)*pow(x,n+1)*cos(x)*cos(x)
|
|
+1.0/(1.0+n)*integral_x_to_n_times_sin_2x(x, n+1);
|
|
}
|
|
|
|
//calculates int x^n * sin(x)^2 dx
|
|
double integral_x_to_n_times_sin_x_2(double x, int n){
|
|
assert(n >= 0);
|
|
return +1.0/(1.0 + n)*pow(x,n+1)*sin(x)*sin(x)
|
|
-1.0/(1.0+n)*integral_x_to_n_times_sin_2x(x, n+1);
|
|
}
|
|
|
|
//calculates int x^n * asin(x) dx
|
|
double integral_x_to_n_times_asin_x(double x, int n){
|
|
assert(n >= 0);
|
|
if (n == 0){
|
|
return x*asin(x)+sqrt(1-x*x);
|
|
}
|
|
else{
|
|
// return 1.0/(n+1.0)*pow(x,n)*(x*asin(x)+sqrt(1-x*x))
|
|
// -n*integral_x_to_n_times_sqrt_1_minus_x2(x, n-1);
|
|
//factor 1/(n+1) missing???
|
|
//fix below!
|
|
return 1.0/(n+1.0)*pow(x,n)*(x*asin(x)+sqrt(1-x*x))
|
|
-n/(n+1.0)*integral_x_to_n_times_sqrt_1_minus_x2(x, n-1);
|
|
}
|
|
}
|
|
|
|
//calculates int x^n * sqrt(1-x^2) dx (recursive using per partes)
|
|
double integral_x_to_n_times_sqrt_1_minus_x2(double x, int n){
|
|
assert(n >= 0);
|
|
if (n == 0) return 0.5*asin(x)+0.5*x*sqrt(1-x*x);
|
|
else return 2.0/(n+2.0)*pow(x, n)*(0.5*asin(x)+0.5*x*sqrt(1-x*x))-n/(n+2.0)*integral_x_to_n_times_asin_x(x, n-1);
|
|
}
|
|
|
|
double integrate_x_to_n_times_sqrt_1_minus_x2(double a, double b, int n){
|
|
return integral_x_to_n_times_sqrt_1_minus_x2(b, n)-integral_x_to_n_times_sqrt_1_minus_x2(a, n);
|
|
}
|
|
|
|
double integral_chebyshev(double x, int n){
|
|
assert(n >= 0);
|
|
if (n == 0) return x;
|
|
else if (n == 1) return 0.5*x*x;
|
|
else return 0.5*(fcnc::chebyshev(x, n+1)/(n+1) - fcnc::chebyshev(x, n-1)/(n-1));
|
|
}
|
|
|
|
double integrate_gauss(double sigma, double mean, double min, double max){
|
|
return 0.5*(TMath::Erf((mean-min)/sqrt(2.0)/sigma) - TMath::Erf((mean-max)/sqrt(2.0)/sigma));
|
|
}
|
|
|
|
double integrate_crystalball(double mean, double sigma, double alpha, double n, double min, double max){
|
|
|
|
//implementation taken from RooFit: https://root.cern.ch/doc/master/RooCBShape_8cxx_source.html#l00108
|
|
|
|
static const double sqrtPiOver2 = 1.2533141373;
|
|
static const double sqrt2 = 1.4142135624;
|
|
|
|
double result = 0.0;
|
|
bool useLog = false;
|
|
|
|
if( fabs(n-1.0) < 1.0e-05 )
|
|
useLog = true;
|
|
|
|
double sig = fabs((Double_t)sigma);
|
|
|
|
//double tmin = (m.min(rangeName)-m0)/sig;
|
|
//double tmax = (m.max(rangeName)-m0)/sig;
|
|
double tmin = (min - mean ) / sigma;
|
|
double tmax = (max - mean ) / sigma;
|
|
|
|
if(alpha < 0) {
|
|
double tmp = tmin;
|
|
tmin = -tmax;
|
|
tmax = -tmp;
|
|
}
|
|
|
|
double absAlpha = fabs(alpha);
|
|
|
|
if( tmin >= -absAlpha ) {
|
|
result += sig*sqrtPiOver2*( TMath::Erf(tmax/sqrt2)
|
|
- TMath::Erf(tmin/sqrt2) );
|
|
}
|
|
else if( tmax <= -absAlpha ) {
|
|
double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
|
|
double b = n/absAlpha - absAlpha;
|
|
|
|
if(useLog) {
|
|
result += a*sig*( log(b-tmin) - log(b-tmax) );
|
|
}
|
|
else {
|
|
result += a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
|
|
- 1.0/(TMath::Power(b-tmax,n-1.0)) );
|
|
}
|
|
}
|
|
else {
|
|
double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
|
|
double b = n/absAlpha - absAlpha;
|
|
|
|
double term1 = 0.0;
|
|
if(useLog) {
|
|
term1 = a*sig*( log(b-tmin) - log(n/absAlpha));
|
|
}
|
|
else {
|
|
term1 = a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
|
|
- 1.0/(TMath::Power(n/absAlpha,n-1.0)) );
|
|
}
|
|
|
|
double term2 = sig*sqrtPiOver2*( TMath::Erf(tmax/sqrt2)
|
|
- TMath::Erf(-absAlpha/sqrt2) );
|
|
|
|
result += term1 + term2;
|
|
}
|
|
|
|
return result;
|
|
|
|
}
|
|
|
|
double integrate_twotailedcrystalball(double mean, double sigma, double alpha1, double alpha2, double n1, double n2, double min, double max){
|
|
|
|
double central = 0.0;
|
|
double left = 0.0;
|
|
double right = 0.0;
|
|
double result = 0.0;
|
|
|
|
assert(alpha1 > 0.0);
|
|
assert(alpha2 > 0.0);
|
|
|
|
//implementation taken from RooFit: RooDoubleCB.cxx
|
|
|
|
static const double rootPiBy2 = sqrt(atan2(0.0,-1.0)/2.0);
|
|
static const double invRoot2 = 1.0/sqrt(2);
|
|
|
|
double invwidth = 1.0/sigma;
|
|
|
|
double tmin = (min-mean)*invwidth;
|
|
double tmax = (max-mean)*invwidth;
|
|
//spdlog::trace("tmin: {0:f}\t tmax: {0:f} ",tmin, tmax);
|
|
|
|
bool isfullrange = (tmin<-1000. && tmax>1000.);
|
|
|
|
//compute gaussian (central) contribution
|
|
double central_low = std::max(min, mean - alpha1*sigma);
|
|
double central_high = std::min(max, mean + alpha2*sigma);
|
|
|
|
double tcentral_low = (central_low-mean)*invwidth;
|
|
double tcentral_high = (central_high-mean)*invwidth;
|
|
|
|
//spdlog::trace("central_low: {0:f}\t central_high: {0:f} ",central_low, central_high);
|
|
//spdlog::trace("tcentral_low: {0:f}\t tcentral_high: {0:f} ",tcentral_low, tcentral_high);
|
|
|
|
if(central_low < central_high){// is the gaussian part in range?
|
|
central = rootPiBy2*sigma*(TMath::Erf(tcentral_high*invRoot2)-TMath::Erf(tcentral_low*invRoot2));
|
|
}
|
|
|
|
//compute left tail;
|
|
if(isfullrange && (n1-1.0)>1.e-3) {
|
|
left = sigma*TMath::Exp(-0.5*alpha1*alpha1)*n1/(alpha1*(n1-1.));
|
|
spdlog::trace("left: {0:f}", left);
|
|
}
|
|
else{
|
|
double left_low = min;
|
|
double left_high = std::min(max, mean - alpha1*sigma);
|
|
double thigh = (left_high-mean)*invwidth;
|
|
|
|
if(left_low < left_high){ //is the left tail in range?
|
|
double n1invalpha1 = n1/(fabs(alpha1));
|
|
if(fabs(n1-1.0)>1.e-3){
|
|
double invn1m1 = 1/(n1-1.);
|
|
double leftpow = TMath::Power(n1invalpha1,-n1*invn1m1);
|
|
assert(leftpow > 0.0);
|
|
double left0 = sigma*TMath::Exp(-0.5*alpha1*alpha1)*invn1m1;
|
|
double left1, left2;
|
|
|
|
if(max > (mean-alpha1*sigma))
|
|
left1 = n1invalpha1;
|
|
else
|
|
left1 = TMath::Power( leftpow*(n1invalpha1 - alpha1 - thigh), 1.-n1);
|
|
|
|
if(tmin<-1000.)
|
|
left2 = 0.;
|
|
else
|
|
left2 = TMath::Power( leftpow*(n1invalpha1 - alpha1 - tmin ), 1.-n1);
|
|
|
|
left = left0*(left1-left2);
|
|
}
|
|
else{
|
|
double A1 = TMath::Power(n1invalpha1,n1)*TMath::Exp(-0.5*alpha1*alpha1);
|
|
double B1 = n1invalpha1-fabs(alpha1);
|
|
left = A1*sigma*(TMath::Log(B1-(left_low-mean)*invwidth) - TMath::Log(B1-(left_high-mean)*invwidth) );
|
|
}
|
|
}
|
|
}
|
|
|
|
//compute right tail;
|
|
if(isfullrange && (n2-1.0)>1.e-3) {
|
|
right = sigma*TMath::Exp(-0.5*alpha2*alpha2)*n2/(alpha2*(n2-1.));
|
|
}
|
|
else{
|
|
double right_low = std::max(min,mean + alpha2*sigma);
|
|
double right_high = max;
|
|
double tlow = (right_low - mean)*invwidth;
|
|
|
|
if(right_low < right_high){ //is the right tail in range?
|
|
double n2invalpha2 = n2/(fabs(alpha2));
|
|
if(fabs(n2-1.0)>1.e-3) {
|
|
double invn2m2 = 1.0/(n2-1.);
|
|
double rightpow = TMath::Power(n2invalpha2,-n2*invn2m2);
|
|
if (rightpow <= 0.0){
|
|
spdlog::error("rightpow < 0.0! Returning NaN");
|
|
spdlog::error("TMath::Power(n2invalpha2,-n2*invn2m2) = {0:.30f}",powl(n2/(fabs(alpha2)),-n2*1.0/(n2-1.)));
|
|
spdlog::error("n2invalpha2: {0:f}\tn2: {1:f}\tinvn2m2: {2:f}",n2invalpha2,n2,invn2m2);
|
|
spdlog::error("mean: {0:f}\tsigma: {1:f}\talpha1: {2:f}\talpha2: {3:f}\tn1: {4:f}\tn2: {5:f}\tmin: {6:f}\tmax: {7:f}", mean, sigma, alpha1, alpha2, n1, n2, min, max);
|
|
return nan("");
|
|
}
|
|
double right0 = sigma*TMath::Exp(-0.5*alpha2*alpha2)*invn2m2;
|
|
double right1, right2;
|
|
|
|
if(min<(mean+alpha2*sigma)) right1 = n2invalpha2;
|
|
else right1 = TMath::Power( rightpow*(n2invalpha2 - alpha2 + tlow), 1.-n2);
|
|
|
|
if(tmax>1000.) right2 = 0.;
|
|
else right2 = TMath::Power( rightpow*(n2invalpha2 - alpha2 + tmax), 1.-n2);
|
|
|
|
right = right0*(right1-right2);
|
|
}
|
|
else{
|
|
double A2 = TMath::Power(n2invalpha2,n2)*TMath::Exp(-0.5*alpha2*alpha2);
|
|
double B2 = n2invalpha2-fabs(alpha2);
|
|
right = A2*sigma*(TMath::Log(B2+(right_high-mean)*invwidth) - TMath::Log(B2+(right_low-mean)*invwidth) );
|
|
}
|
|
}
|
|
}
|
|
|
|
result = left + central + right;
|
|
|
|
if(result <= 0.0){
|
|
spdlog::error("left = {0:e}",left);
|
|
spdlog::error("central = {0:e}",central);
|
|
spdlog::error("right = {0:e}", right);
|
|
spdlog::error("sigma = {0:e}", sigma);
|
|
spdlog::error("n1 = {0:e}", n1);
|
|
spdlog::error("n2 = {0:e}", n2);
|
|
spdlog::error("alpha1 = {0:e}", alpha1);
|
|
spdlog::error("alpha1 = {0:e}", alpha2);
|
|
assert(0);
|
|
}
|
|
|
|
if(std::isnan(result) || std::isinf(result)){
|
|
spdlog::critical("Determination of twotailed CB integral failed: {0:f}", result);
|
|
spdlog::critical("mean: {1:f}\tsigma: {2:f}\talpha1: {3:f}\talpha2: {4:f}\tn1: {5:f}\tn2: {6:f}", mean, sigma, alpha1, alpha2, n1, n2);
|
|
spdlog::critical("left: {0:f}\tcentral: {1:f}\tright: {2:f}", left, central, right);
|
|
|
|
assert(0);
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
double integral_x_to_n_times_exp_minus_x(double x, double tau, int n){
|
|
assert(n>=0);
|
|
if (n == 0) return -tau*exp(-x/tau);
|
|
else return pow(x, n)*(-tau*exp(-x/tau)) + n * tau * integral_x_to_n_times_exp_minus_x(x, tau, n-1);
|
|
}
|
|
|
|
|
|
|
|
//-----------------------
|
|
// P wave
|
|
//-----------------------
|
|
|
|
//THESE FUNCTIONS NEED TO BE CHECKED
|
|
//i correspons to ctl, j to ctk and k to phi
|
|
|
|
|
|
//int[ sin^2(tk) ] for J1s
|
|
double integral_f1(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin^2 d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return (integral_x_to_n(ctk,n_ctk)-integral_x_to_n(ctk,n_ctk+2))
|
|
*integral_x_to_n(ctl,n_ctl)
|
|
*integral_x_to_n(phi,n_phi);
|
|
}
|
|
double integrate_f1(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
return (integrate_x_to_n(ctk_a, ctk_b, n_ctk)-integrate_x_to_n(ctk_a,ctk_b,n_ctk+2))
|
|
*integrate_x_to_n(ctl_a, ctl_b, n_ctl)
|
|
*integrate_x_to_n(phi_a, phi_b, n_phi);
|
|
}
|
|
|
|
//int[ cos^2(tk) ] for J1c
|
|
double integral_f2(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate cos^2 thetak d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return integral_x_to_n(ctk,n_ctk+2)
|
|
*integral_x_to_n(ctl,n_ctl)
|
|
*integral_x_to_n(phi,n_phi);
|
|
}
|
|
double integrate_f2(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
return integrate_x_to_n(ctk_a, ctk_b, n_ctk+2)
|
|
*integrate_x_to_n(ctl_a, ctl_b, n_ctl)
|
|
*integrate_x_to_n(phi_a, phi_b, n_phi);
|
|
}
|
|
|
|
//int[ sin^2(tk) * cos(2tl) ] for J2s
|
|
double integral_f3(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin^2 thetak cos 2thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return (integral_x_to_n(ctk,n_ctk)-integral_x_to_n(ctk,n_ctk+2))
|
|
*(2*integral_x_to_n(ctl,n_ctl+2)-integral_x_to_n(ctl,n_ctl)) //cos(2x) = 2*cos^2(x)-1
|
|
*integral_x_to_n(phi,n_phi);
|
|
}
|
|
double integrate_f3(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
return (integrate_x_to_n(ctk_a, ctk_b, n_ctk)-integrate_x_to_n(ctk_a, ctk_b, n_ctk+2))
|
|
*(2*integrate_x_to_n(ctl_a, ctl_b, n_ctl+2)-integrate_x_to_n(ctl_a, ctl_b, n_ctl)) //cos(2x) = 2*cos^2(x)-1
|
|
*integrate_x_to_n(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
//int[ cos^2(tk) * cos(2tl) ] for J2c
|
|
double integral_f4(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate cos^2 thetak cos 2thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return integral_x_to_n(ctk,n_ctk+2)
|
|
*(2*integral_x_to_n(ctl,n_ctl+2)-integral_x_to_n(ctl,n_ctl))
|
|
*integral_x_to_n(phi,n_phi);
|
|
}
|
|
double integrate_f4(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate cos^2 thetak cos 2thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return integrate_x_to_n(ctk_a, ctk_b, n_ctk+2)
|
|
*(2*integrate_x_to_n(ctl_a, ctl_b, n_ctl+2)-integrate_x_to_n(ctl_a, ctl_b, n_ctl))
|
|
*integrate_x_to_n(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
//int[ sin^2(tk) * sin^2(tl) * cos(2phi) ] for J3
|
|
double integral_f5(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin^2 thetak sin^2 thetal cos 2phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return (integral_x_to_n(ctk,n_ctk)-integral_x_to_n(ctk,n_ctk+2))
|
|
*(integral_x_to_n(ctl,n_ctl)-integral_x_to_n(ctl,n_ctl+2))
|
|
*(integral_x_to_n_times_cos_2x(phi,n_phi));
|
|
}
|
|
double integrate_f5(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin^2 thetak sin^2 thetal cos 2phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return (integrate_x_to_n(ctk_a, ctk_b, n_ctk)-integrate_x_to_n(ctk_a, ctk_b, n_ctk+2))
|
|
*(integrate_x_to_n(ctl_a, ctl_b, n_ctl)-integrate_x_to_n(ctl_a, ctl_b, n_ctl+2))
|
|
*(integrate_x_to_n_times_cos_2x(phi_a,phi_b,n_phi));
|
|
}
|
|
|
|
//int[ sin(2tk) * sin(2tl) * cos(phi) ] for J4
|
|
double integral_f6(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin 2thetak sin 2thetal cos phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk+1)
|
|
*2*integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl+1)
|
|
*integral_x_to_n_times_cos_x(phi,n_phi);
|
|
}
|
|
double integrate_f6(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin 2thetak sin 2thetal cos phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a, ctk_b, n_ctk+1)
|
|
*2*integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a, ctl_b,n_ctl+1)
|
|
*integrate_x_to_n_times_cos_x(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
//int [ sin(2tk) * sin(tl) * cos(phi) ] for J5
|
|
double integral_f7(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin 2thetak sin thetal cos phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk+1)
|
|
*integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl)
|
|
*integral_x_to_n_times_cos_x(phi,n_phi);
|
|
}
|
|
double integrate_f7(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin 2thetak sin thetal cos phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a, ctk_b, n_ctk+1)
|
|
*integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a, ctl_b,n_ctl)
|
|
*integrate_x_to_n_times_cos_x(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
//int [ sin^2(tk) * cos(tl) ] for J6s
|
|
double integral_f8(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin^2 thetak cos thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return (integral_x_to_n(ctk,n_ctk)-integral_x_to_n(ctk,n_ctk+2))
|
|
*integral_x_to_n(ctl,n_ctl+1)
|
|
*integral_x_to_n(phi,n_phi);
|
|
}
|
|
double integrate_f8(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin^2 thetak cos thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return (integrate_x_to_n(ctk_a, ctk_b, n_ctk)-integrate_x_to_n(ctk_a, ctk_b, n_ctk+2))
|
|
*integrate_x_to_n(ctl_a, ctl_b,n_ctl+1)
|
|
*integrate_x_to_n(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
//int [ cos^2(tk) * cos(tl) ] for J6c
|
|
double integral_f9(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin^2 thetak cos thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return integral_x_to_n(ctk,n_ctk+2)
|
|
*integral_x_to_n(ctl,n_ctl+1)
|
|
*integral_x_to_n(phi,n_phi);
|
|
}
|
|
double integrate_f9(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin^2 thetak cos thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return integrate_x_to_n(ctk_a, ctk_b, n_ctk+2)
|
|
*integrate_x_to_n(ctl_a, ctl_b,n_ctl+1)
|
|
*integrate_x_to_n(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
//int [ sin(2tk) * sin(tl) * sin(phi) ] for J7
|
|
double integral_f10(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin 2thetak sin thetal sin phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk+1)
|
|
*integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl+1)
|
|
*integral_x_to_n_times_sin_x(phi,n_phi);
|
|
}
|
|
double integrate_f10(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin 2thetak sin thetal sin phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a, ctk_b, n_ctk+1)
|
|
*integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a, ctl_b,n_ctl+1)
|
|
*integrate_x_to_n_times_sin_x(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
//int [ sin(2tk) * sin(2tl) * sin(phi) ] for J8
|
|
double integral_f11(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin 2thetak sin 2thetal sin phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk+1)
|
|
*2*integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl+1)
|
|
*integral_x_to_n_times_sin_x(phi,n_phi);
|
|
}
|
|
double integrate_f11(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin 2thetak sin 2thetal sin phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a, ctk_b,n_ctk+1)
|
|
*2*integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a, ctl_b,n_ctl+1)
|
|
*integrate_x_to_n_times_sin_x(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
//int [ sin^2(tk) * sin^2(tl) * sin(2phi) ] for J9
|
|
double integral_f12(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin^2 thetak sin^2 thetal sin 2phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return (integral_x_to_n(ctk,n_ctk)-integral_x_to_n(ctk,n_ctk+2))
|
|
*(integral_x_to_n(ctl,n_ctl)-integral_x_to_n(ctl,n_ctl+2))
|
|
*(integral_x_to_n_times_sin_2x(phi,n_phi));
|
|
}
|
|
double integrate_f12(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
//Integrate sin^2 thetak sin^2 thetal sin 2phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
|
|
return (integrate_x_to_n(ctk_a, ctk_b,n_ctk)-integrate_x_to_n(ctk_a, ctk_b,n_ctk+2))
|
|
*(integrate_x_to_n(ctl_a, ctl_b,n_ctl)-integrate_x_to_n(ctl_a, ctl_b,n_ctl+2))
|
|
*(integrate_x_to_n_times_sin_2x(phi_a,phi_b,n_phi));
|
|
}
|
|
|
|
//-----------------------
|
|
// S wave
|
|
//-----------------------
|
|
|
|
double integral_s_f1(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
return (integral_x_to_n(ctl,n_ctl)-integral_x_to_n(ctl,n_ctl+2))
|
|
*integral_x_to_n(ctk,n_ctk)
|
|
*integral_x_to_n(phi, n_phi);
|
|
}
|
|
double integrate_s_f1(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
return (integrate_x_to_n(ctl_a,ctl_b,n_ctl)-integrate_x_to_n(ctl_a,ctl_b,n_ctl+2))
|
|
*integrate_x_to_n(ctk_a,ctk_b,n_ctk)
|
|
*integrate_x_to_n(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
double integral_s_f2(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
return (integral_x_to_n(ctl,n_ctl)-integral_x_to_n(ctl,n_ctl+2))
|
|
*integral_x_to_n(ctk,n_ctk+1)
|
|
*integral_x_to_n(phi, n_phi);
|
|
}
|
|
double integrate_s_f2(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
return (integrate_x_to_n(ctl_a,ctl_b,n_ctl)-integrate_x_to_n(ctl_a,ctl_b,n_ctl+2))
|
|
*integrate_x_to_n(ctk_a,ctk_b,n_ctk+1)
|
|
*integrate_x_to_n(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
double integral_s_f3(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl+1)
|
|
*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk)
|
|
*integral_x_to_n_times_cos_x(phi, n_phi);
|
|
}
|
|
double integrate_s_f3(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a,ctl_b,n_ctl+1)
|
|
*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a,ctk_b,n_ctk)
|
|
*integrate_x_to_n_times_cos_x(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
double integral_s_f4(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
return integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl)
|
|
*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk)
|
|
*integral_x_to_n_times_cos_x(phi, n_phi);
|
|
}
|
|
double integrate_s_f4(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
return integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a,ctl_b,n_ctl)
|
|
*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a,ctk_b,n_ctk)
|
|
*integrate_x_to_n_times_cos_x(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
double integral_s_f5(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
return integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl)
|
|
*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk)
|
|
*integral_x_to_n_times_sin_x(phi, n_phi);
|
|
}
|
|
double integrate_s_f5(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
return integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a,ctl_b,n_ctl)
|
|
*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a,ctk_b,n_ctk)
|
|
*integrate_x_to_n_times_sin_x(phi_a,phi_b,n_phi);
|
|
}
|
|
|
|
double integral_s_f6(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
|
|
return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl+1)
|
|
*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk)
|
|
*integral_x_to_n_times_sin_x(phi, n_phi);
|
|
}
|
|
double integrate_s_f6(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
|
|
return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a,ctl_b,n_ctl+1)
|
|
*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a,ctk_b,n_ctk)
|
|
*integrate_x_to_n_times_sin_x(phi_a,phi_b,n_phi);
|
|
}
|
|
|