rtang/ADPlusPlus
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

use Math/IFunction.h for the Legendre Polynormail

Browse Source
This commit is contained in:
Ryan Tang 2022-10-13 13:15:28 -04:00
parent d71f18d52e
commit 593226ab13
2 changed files with 19 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
Qk.h
View File

@ -3,6 +3,8 @@
#include <cstdio> #include <cstdio>
#include "Math/IFunction.h"
double TauCal(double Energy_MeV){ double TauCal(double Energy_MeV){
double E_log = log(Energy_MeV); double E_log = log(Energy_MeV);
@ -16,13 +18,6 @@ double TauCal(double Energy_MeV){
return exp(TT); return exp(TT);
} }
double LegendreP(int n, double theta){
if( n == 0 ) return 1;
if( n == 2 ) return (3. * cos(theta) * cos(theta) -1 )/2.;
if( n == 4 ) return (35 * pow( cos(theta), 4) - 30 * pow(cos(theta),2) + 3.) /8.;
return 0;
}
double * QK(double Energy_keV, double radius_cm, double distance_cm, double thickness_cm){ double * QK(double Energy_keV, double radius_cm, double distance_cm, double thickness_cm){
@ -61,15 +56,15 @@ double * QK(double Energy_keV, double radius_cm, double distance_cm, double thic
double ex1 = exp( - 1.0 * Tau * thickness_cm / cos(beta1) ); double ex1 = exp( - 1.0 * Tau * thickness_cm / cos(beta1) );
double term1 = A * (1-ex1) * sin(beta1) * delx1; double term1 = A * (1-ex1) * sin(beta1) * delx1;
sum1 += LegendreP(2, beta1) * term1; sum1 += std::legendre(2, cos(beta1)) * term1;
sum2 += LegendreP(4, beta1) * term1; sum2 += std::legendre(4, cos(beta1)) * term1;
sum3 += term1; sum3 += term1;
double ex2 = exp( -1 * Tau * (radius_cm / sin(beta2) - distance_cm /cos(beta2)) ); double ex2 = exp( -1 * Tau * (radius_cm / sin(beta2) - distance_cm /cos(beta2)) );
double term2 = A * (1-ex2) * sin(beta2) * delx2; double term2 = A * (1-ex2) * sin(beta2) * delx2;
sum4 += LegendreP(2, beta2) * term2; sum4 += std::legendre(2, cos(beta2)) * term2;
sum5 += LegendreP(4, beta2) * term2; sum5 += std::legendre(4, cos(beta2)) * term2;
sum6 += term2; sum6 += term2;
//if( i % 75 == 0) printf("%4d | %10.6f, %10.6f, %10.6f %d | %10.6f, %10.6f %10.6f %10.6f, %10.6f %10.6f \n", i, beta1, beta2, A, J, sum1, sum2, sum3, sum4, sum5, sum6); //if( i % 75 == 0) printf("%4d | %10.6f, %10.6f, %10.6f %d | %10.6f, %10.6f %10.6f %10.6f, %10.6f %10.6f \n", i, beta1, beta2, A, J, sum1, sum2, sum3, sum4, sum5, sum6);
@ -88,6 +83,15 @@ double * QK(double Energy_keV, double radius_cm, double distance_cm, double thic
} }
/// Relic code /// Relic code
double LegendreP(int n, double theta){
if( n == 0 ) return 1;
if( n == 2 ) return (3. * cos(theta) * cos(theta) -1 )/2.;
if( n == 4 ) return (35 * pow( cos(theta), 4) - 30 * pow(cos(theta),2) + 3.) /8.;
return 0;
}
double QK2(double Energy, double radius, double distance, double thickness){ double QK2(double Energy, double radius, double distance, double thickness){

7
ad++.h
View File

@ -13,6 +13,7 @@
#include <TApplication.h> #include <TApplication.h>
#include <TCanvas.h> #include <TCanvas.h>
#include "Math/IFunction.h"
#define PI 3.14159265358979323846 #define PI 3.14159265358979323846
@ -31,7 +32,7 @@ double YE(double * x , double *par){
/// par[1] = a2; /// par[1] = a2;
/// par[2] = a4; /// par[2] = a4;
return par[0] + par[1] * LegendreP(2, x[0]) + par[2] * LegendreP(4, x[0]); return par[0] + par[1] * std::legendre(2, cos(x[0])) + par[2] * std::legendre(4, cos(x[0]));
} }
@ -112,7 +113,7 @@ double TheoryAD(double * x, double *par){
/// par[3] = J; /// par[3] = J;
double result = 0; double result = 0;
for( int k = 0; k <= 4; k += 2){ for( int k = 0; k <= 4; k += 2){
result += Q[k]*Bk(k, par[3], par[2]) * LegendreP(k, x[0]) * ( R1[k] + 2*par[1]*R2[k] + par[1]*par[1]*R3[k] ) / (1 + par[1]*par[1] ); result += Q[k]*Bk(k, par[3], par[2]) * std::legendre(k, cos(x[0])) * ( R1[k] + 2*par[1]*R2[k] + par[1]*par[1]*R3[k] ) / (1 + par[1]*par[1] );
} }
return result * par[0]; return result * par[0];
@ -313,7 +314,7 @@ void fitOldAD(std::vector<std::vector<double>> data_deg_count_error,
double YT = 0; double YT = 0;
for( int k = 0; k <= 4; k += 2){ for( int k = 0; k <= 4; k += 2){
YT += Q[k] * B[k] * LegendreP(k, x[i] * PI/180.) * ( R1[k] + 2 * delta * R2[k] + delta*delta* R3[k] ) / (1 + delta * delta) ; YT += Q[k] * B[k] * std::legendre(k, cos(x[i] * PI/180.)) * ( R1[k] + 2 * delta * R2[k] + delta*delta* R3[k] ) / (1 + delta * delta) ;
} }
chiSq += pow( A0 * YT - y[i], 2)/ dataSize / ey[i] / ey[i]; chiSq += pow( A0 * YT - y[i], 2)/ dataSize / ey[i] / ey[i];