AD code in C++ with ROOT compatible
This commit is contained in:
Ryan@iMac
commit
1a8a165b80
3
.gitignore
vendored
Normal file
3
.gitignore
vendored
Normal file
|
|
@ -0,0 +1,3 @@
|
|||
ad++
|
||||
testJSymbol
|
||||
testJSymbol.cpp
|
||||
267
Qk.h
Normal file
267
Qk.h
Normal file
|
|
@ -0,0 +1,267 @@
|
|||
double TauCal(double Energy_MeV){
|
||||
|
||||
double E_log = log(Energy_MeV);
|
||||
double EL1 = E_log;
|
||||
double EL2 = pow(E_log,2);
|
||||
double EL3 = pow(E_log,3);
|
||||
double EL4 = pow(E_log,4);
|
||||
double EL5 = pow(E_log,5);
|
||||
double TT = -1.1907 -0.5372*EL1 - 0.0438*EL2 + 0.0218*EL3 + 0.0765*EL4 + 0.0095*EL5;
|
||||
|
||||
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 E_MeV = Energy_keV/1000;
|
||||
double Tau = TauCal(E_MeV);
|
||||
|
||||
double alpha = atan( radius_cm / (distance_cm + thickness_cm) );
|
||||
double gamma = atan( radius_cm / distance_cm );
|
||||
|
||||
double sum1 = 0,sum2 = 0,sum3 = 0.;
|
||||
double sum4 = 0,sum5 = 0,sum6 = 0.;
|
||||
|
||||
double beta1 = 0;
|
||||
double beta2 = 0;
|
||||
double A = 0.;
|
||||
|
||||
int div = 1000;
|
||||
|
||||
double delx1 = (alpha)/div; // mrad
|
||||
double delx2 = (gamma-alpha)/div;
|
||||
|
||||
for(int i = 0; i<=div; i++){
|
||||
|
||||
int J = i % 2;
|
||||
if( J == 0 ){
|
||||
A = 2.0;
|
||||
}else{
|
||||
A = 4.0;
|
||||
beta1 = i * delx1;
|
||||
beta2 = alpha + i * delx2;
|
||||
}
|
||||
|
||||
if( i == 0 || i == div ) {
|
||||
A = 1.0;
|
||||
beta1 = i * delx1;
|
||||
beta2 = alpha + i * delx2;
|
||||
}
|
||||
|
||||
double cosb = cos(beta1);
|
||||
double sinb = sin(beta1);
|
||||
double ex1 = exp( - 1.0 * Tau * thickness_cm / cosb );
|
||||
double term3 = (1-ex1) * sinb * A * delx1;
|
||||
double term1 = LegendreP(2, beta1) * term3;
|
||||
double term2 = LegendreP(4, beta1) * term3;
|
||||
|
||||
sum1 += term1;
|
||||
sum2 += term2;
|
||||
sum3 += term3;
|
||||
|
||||
cosb = cos(beta2);
|
||||
sinb = sin(beta2);
|
||||
double ex2 = exp( -1 * Tau * (radius_cm / sinb - distance_cm /cosb) );
|
||||
double term6 = A * (1-ex2) * sinb * delx2;
|
||||
double term4 = LegendreP(2, beta2) * term6;
|
||||
double term5 = LegendreP(4, beta2) * term6;
|
||||
|
||||
sum4 = sum4 +term4;
|
||||
sum5 = sum5 +term5;
|
||||
sum6 = sum6 +term6;
|
||||
|
||||
//if( i % 75 == 0) printf("%4d | %10.6f, %10.6f %d | %10.6f, %10.6f %10.6f \n", i, beta, A, J, sum1, sum2, sum3);
|
||||
}
|
||||
|
||||
double ans1 = sum1/3;
|
||||
double ans2 = sum2/3;
|
||||
double ans3 = sum3/3;
|
||||
double ans4 = sum4/3;
|
||||
double ans5 = sum5/3;
|
||||
double ans6 = sum6/3;
|
||||
|
||||
double * Qk = new double[2];
|
||||
Qk[0] = (ans1+ans4)/(ans3+ans6);
|
||||
Qk[1] = (ans2+ans5)/(ans3+ans6);
|
||||
|
||||
printf("--------------\n");
|
||||
printf(" QD2 = %lf\n",Qk[0]);
|
||||
printf(" QD4 = %lf\n",Qk[1]);
|
||||
printf("--------------\n");
|
||||
|
||||
return Qk;
|
||||
}
|
||||
|
||||
/// Relic code
|
||||
double QK2(double Energy, double radius, double distance, double thickness){
|
||||
|
||||
|
||||
double Qkn = 0.;
|
||||
double E_mev = Energy/1000;
|
||||
|
||||
double E_log = log(E_mev);
|
||||
double EL1 = E_log;
|
||||
double EL2 = pow(E_log,2);
|
||||
double EL3 = EL1*EL2;
|
||||
double EL4 = pow(EL2,2);
|
||||
double EL5 = EL4*EL1;
|
||||
double TT = -1.1907 -0.5372*EL1 - 0.0438*EL2 + 0.0218*EL3 + 0.0765*EL4 + 0.0095*EL5;
|
||||
double Tau = exp(TT);
|
||||
// printf("TLN = %lf\n",TT);
|
||||
// printf("Tau = %lf\n",Tau);
|
||||
//calulating attenuation angles
|
||||
|
||||
double Z1 = radius / (distance + thickness);
|
||||
double Z2 = radius / distance;
|
||||
double alpha = atan(Z1);
|
||||
double gamma = atan(Z2);
|
||||
double beta;
|
||||
double BL = 0.;
|
||||
double BU = alpha;
|
||||
double A = 0.;
|
||||
double delx1 = (BU-BL)/1000;
|
||||
// printf("alpha = %lf\n",alpha);
|
||||
// printf("gamma = %lf\n",gamma);
|
||||
// printf("delx1 = %lf\n",delx1);
|
||||
double sum1 = 0,sum2 = 0,sum3 = 0.;
|
||||
double sum4 = 0,sum5 = 0,sum6 = 0.;
|
||||
double cosb,sinb,secb,cscb,c2,c4,fac1,fac2,ex1,ex2 = 0.;
|
||||
double term1 = 0,term2 = 0,term3 = 0.;
|
||||
double term4 = 0,term5 = 0,term6 = 0.;
|
||||
int J=0;
|
||||
|
||||
int loop_length = 1000;
|
||||
for(int i = 0; i<=loop_length; i++){
|
||||
/*
|
||||
if(i > 0 and i < loop_length){
|
||||
J = i%2;
|
||||
//printf("\t\ti = %d\nJ=%d\n",i,J);
|
||||
if(J==0){A=2.;
|
||||
}else {A=4.;}
|
||||
beta = BL+i+delx1;
|
||||
}else{A=.1;beta = BL+i+delx1;}
|
||||
*/
|
||||
if(i != 0){
|
||||
if(i != loop_length){
|
||||
J = i%2;
|
||||
if(J==0){
|
||||
A = 2.;
|
||||
}else{
|
||||
A=4.;
|
||||
beta = BL+i*delx1;
|
||||
}
|
||||
}else{
|
||||
A=1.0;
|
||||
beta = BL+i*delx1;}
|
||||
}else{
|
||||
A =1.0;
|
||||
beta = BL+i*delx1;
|
||||
}
|
||||
// printf("Beta = %lf\n",beta);
|
||||
|
||||
cosb = cos(beta);
|
||||
sinb = sin(beta);
|
||||
secb = 1.0/cosb;
|
||||
c2 = pow(cosb,2);
|
||||
c4 = pow(cosb,4);
|
||||
|
||||
fac1 = -1 *Tau *thickness *secb;
|
||||
|
||||
ex1 = exp(fac1);
|
||||
term1 = 0.5*(3*c2-1)*(1-ex1)*sinb*A*delx1;
|
||||
term2 = 0.125*A*(35*c4-30*c2+3)*(1-ex1)*sinb*delx1;
|
||||
term3 = A*(1-ex1)*sinb*delx1;
|
||||
|
||||
sum1 = sum1 +term1;
|
||||
sum2 = sum2 +term2;
|
||||
sum3 = sum3 +term3;
|
||||
|
||||
//if( i % 75 == 0) printf("%4d | %10.6f, %10.6f, %d | %10.6f, %10.6f %10.6f \n", i, beta, A, J, sum1, sum2, sum3);
|
||||
|
||||
}
|
||||
|
||||
double ans1 = sum1/3;
|
||||
double ans2 = sum2/3;
|
||||
double ans3 = sum3/3;
|
||||
|
||||
//printf("%10.6f, %10.6f %10.6f \n", ans1, ans2, ans3);
|
||||
|
||||
double LB=alpha;
|
||||
double UB=gamma;
|
||||
double delx2 = (UB-LB)/1000;
|
||||
for(int i = 0; i<=loop_length; i++){
|
||||
/*
|
||||
if(i > 0 and i < loop_length){
|
||||
J2 = i%2;
|
||||
//printf("\t\ti = %d\nJ=%d\n",i,J2);
|
||||
if(J2==0){B=2.;
|
||||
}else {B=4.;}
|
||||
beta2 = LB+i+delx2;
|
||||
}else{B=.1;beta2 = LB+i+delx2;}
|
||||
*/
|
||||
if(i != 0){
|
||||
if(i != loop_length){
|
||||
J = i%2;
|
||||
if(J==0){
|
||||
A = 2.;
|
||||
}else{A=4.;beta = LB+i*delx2;}
|
||||
}else{A=1.0;beta = LB+i*delx2;}
|
||||
}else{A =1.0;beta = LB+i*delx2;}
|
||||
// printf("Beta1 = %lf\n",beta);
|
||||
cosb = cos(beta);
|
||||
sinb = sin(beta);
|
||||
secb = 1.0/cosb;
|
||||
cscb = 1.0/sinb;
|
||||
c2 = pow(cosb,2);
|
||||
c4 = pow(cosb,4);
|
||||
fac2 = -1 *Tau *(radius*cscb -distance*secb);
|
||||
ex2 = exp(fac2);
|
||||
term4 = 0.5*A*(3*c2-1)*(1-ex2)*sinb*delx2;
|
||||
term5 = 0.125*A*(35*c4-30*c2+3)*(1-ex2)*sinb*delx2;
|
||||
term6 = A*(1-ex2)*sinb*delx2;
|
||||
sum4 = sum4 +term4;
|
||||
sum5 = sum5 +term5;
|
||||
sum6 = sum6 +term6;
|
||||
}
|
||||
double ans4=sum4/3;
|
||||
double ans5=sum5/3;
|
||||
double ans6=sum6/3;
|
||||
/*
|
||||
printf("ans1:%lf\n",ans1*100);
|
||||
printf("ans2:%lf\n",ans2*100);
|
||||
printf("ans3:%lf\n",ans3*100);
|
||||
printf("ans4:%lf\n",ans4*100);
|
||||
printf("ans5:%lf\n",ans5*100);
|
||||
printf("ans6:%lf\n",ans6*100);
|
||||
*/
|
||||
double QD2 = (ans1+ans4)/(ans3+ans6);
|
||||
double QD4 = (ans2+ans5)/(ans3+ans6);
|
||||
|
||||
printf("--------------\n");
|
||||
printf(" QD2 = %lf\n",QD2);
|
||||
printf(" QD4 = %lf\n",QD4);
|
||||
printf("--------------\n");
|
||||
|
||||
/*
|
||||
//Now output a file that contains R, D , T , gamma energy, attentuation coeff, q2 and q4
|
||||
ofstream fileo;
|
||||
fileo.open ("ad.txt");
|
||||
fileo << "Radius = " << radius <<" [cm]\n";
|
||||
fileo << "Distance = " << distance <<" [cm]\n";
|
||||
fileo << "Thickness = " << thickness <<" [cm]\n";
|
||||
fileo << "Atten.C = " << Tau <<" [cm^-1]\n";
|
||||
fileo << "Gamma_E = " << Energy <<" [KeV]\n";
|
||||
fileo << "QD2 = " << QD2 << "\n";
|
||||
fileo << "QD4 = " << QD4 << "\n";
|
||||
fileo.close();
|
||||
*/
|
||||
return QD2;
|
||||
}
|
||||
300
ad++.cpp
Normal file
300
ad++.cpp
Normal file
|
|
@ -0,0 +1,300 @@
|
|||
#include "jSymbol.h"
|
||||
#include "Qk.h"
|
||||
|
||||
#include <TROOT.h>
|
||||
#include <TSystem.h>
|
||||
#include <TAxis.h>
|
||||
#include <TF1.h>
|
||||
#include <TLatex.h>
|
||||
#include <TGraphErrors.h>
|
||||
#include <TApplication.h>
|
||||
#include <TCanvas.h>
|
||||
|
||||
|
||||
#define PI 3.14159265358979323846
|
||||
|
||||
double Racah(int j1, int j2, int J, int j3, int j12, int j23);
|
||||
double Rk(int k, int L1, int L2, int J1, int J2);
|
||||
void PrintRk(int k);
|
||||
double w(int M, int J);
|
||||
void Print_w_sum();
|
||||
double Bk(int k, int J);
|
||||
double YE(double * x , double *par);/// for root Fit
|
||||
|
||||
|
||||
/// This is for general fit a, delta
|
||||
double Q[5] = {0};
|
||||
double B[5] = {0};
|
||||
double R1[5] = {0};
|
||||
double R2[5] = {0};
|
||||
double R3[5] = {0};
|
||||
double Fit(double * x, double *par);
|
||||
|
||||
int main(int argc, char **argv){
|
||||
|
||||
TApplication theApp("App",&argc,argv);
|
||||
|
||||
//========================== User input
|
||||
double energy_keV = 832;
|
||||
double detRadius_cm = 2;
|
||||
double targetDistance_cm = 20;
|
||||
double detThickness_cm = 10;
|
||||
|
||||
double data[][3] = { { 150.00, 2441.44, 122},
|
||||
{ 131.75, 2580.11, 129},
|
||||
{ 90.00, 4652.08, 232},
|
||||
{ 48.75, 3023.17, 151} };
|
||||
|
||||
|
||||
const int Ji = 5;
|
||||
const int Jf = 4;
|
||||
//======================================
|
||||
|
||||
/// detector acceptance
|
||||
double * Qk = QK(energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm);
|
||||
Q[0] = 1;
|
||||
Q[2] = Qk[0];
|
||||
Q[4] = Qk[1];
|
||||
|
||||
printf("Qk2 : %f \n", Q[2]);
|
||||
printf("Qk4 : %f \n", Q[4]);
|
||||
|
||||
const int dataSize = sizeof(data)/sizeof(double)/3;
|
||||
|
||||
/// for TGraphErrors
|
||||
double x[dataSize];
|
||||
double y[dataSize];
|
||||
double ex[dataSize];
|
||||
double ey[dataSize];
|
||||
|
||||
printf("============= Data :\n");
|
||||
for( int i = 0; i < dataSize; i++){
|
||||
printf("%2d | %8.2f, %8.2f(%4.0f) \n", i, data[i][0], data[i][1], data[i][2]);
|
||||
|
||||
x[i] = data[i][0] * PI/180;
|
||||
y[i] = data[i][1];
|
||||
ey[i] = data[i][2];
|
||||
ex[i] = 0.;
|
||||
}
|
||||
printf("======================\n");
|
||||
|
||||
TGraphErrors * gExp = new TGraphErrors( dataSize, x, y, ex, ey);
|
||||
gExp->SetTitle("");
|
||||
gExp->GetXaxis()->SetTitle("Angle [rad]");
|
||||
gExp->GetYaxis()->SetTitle("Data");
|
||||
|
||||
TCanvas * c1 = new TCanvas("c1", "c1", 1000, 500);
|
||||
c1->Divide(2, 1);
|
||||
c1->cd(1);
|
||||
gExp->Draw("AP*");
|
||||
|
||||
c1->Modified();
|
||||
gSystem->ProcessEvents();
|
||||
|
||||
|
||||
///======== Fitting the experimental distribution with a0+ a2*P(2,cos(theta)) + a4 * P(4, cos(theta))
|
||||
TF1 * f1 = new TF1("f1", YE, 0, PI, 3);
|
||||
|
||||
f1->SetLineColor(4);
|
||||
f1->SetLineWidth(2);
|
||||
f1->SetNpx(1000);
|
||||
f1->SetParameter(0, 2000);
|
||||
f1->SetParameter(1, -2000);
|
||||
f1->SetParameter(2, 2000);
|
||||
|
||||
gExp->Fit("f1", "");
|
||||
|
||||
const Double_t* paraE = f1->GetParErrors();
|
||||
const Double_t* paraA = f1->GetParameters();
|
||||
|
||||
double A0 = paraA[0];
|
||||
double A2 = paraA[1];
|
||||
double A4 = paraA[2];
|
||||
|
||||
///=================================== Fit with Theritical
|
||||
|
||||
|
||||
int L = abs(Ji - Jf);
|
||||
if( L == 0 ) L = 1;
|
||||
|
||||
|
||||
for( int k = 0; k <= 4; k += 2){
|
||||
B[k] = Bk(k, Ji);
|
||||
R1[k] = Rk(k, L , L , Ji, Jf);
|
||||
R2[k] = Rk(k, L , L+1, Ji, Jf);
|
||||
R3[k] = Rk(k, L+1, L+1, Ji, Jf);
|
||||
}
|
||||
|
||||
std::vector<double> deltaDeg;
|
||||
std::vector<double> LogChiSq;
|
||||
|
||||
for( float deltaAngle = -90; deltaAngle <= 90 ; deltaAngle += 2. ){
|
||||
|
||||
double delta = tan(deltaAngle * PI/ 180.);
|
||||
|
||||
double chiSq = 0;
|
||||
for( int i = 0; i < dataSize; i++){
|
||||
|
||||
double YT = 0;
|
||||
|
||||
for( int k = 0; k <= 4; k += 2){
|
||||
YT += Q[k] * B[k] * LegendreP(k, x[i]) * ( R1[k] + 2 * delta * R2[k] + delta*delta* R3[k] ) / (1 + delta * delta) ;
|
||||
}
|
||||
|
||||
//printf(" YT : %f , YE : %f \n", A0 * YT, y[i]);
|
||||
|
||||
chiSq += pow( A0 * YT - y[i], 2)/ dataSize / ey[i] / ey[i];
|
||||
|
||||
}
|
||||
|
||||
deltaDeg.push_back(deltaAngle);
|
||||
LogChiSq.push_back(log(chiSq));
|
||||
//printf(" %6.2f deg, %8.3f \n", deltaAngle, log(chiSq));
|
||||
}
|
||||
|
||||
c1->cd(2);
|
||||
TGraph * gDelta = new TGraph((int) deltaDeg.size(), &deltaDeg[0], &LogChiSq[0]);
|
||||
gDelta->SetTitle("");
|
||||
gDelta->GetXaxis()->SetTitle("aTan(Mixing Ratio) [Deg]");
|
||||
gDelta->GetYaxis()->SetTitle("Log(chi-sq)");
|
||||
|
||||
gDelta->Draw("APL*");
|
||||
|
||||
c1->Modified();
|
||||
gSystem->ProcessEvents();
|
||||
|
||||
///=============================== Fit a, delta at once
|
||||
|
||||
TF1 * fit = new TF1("fit", Fit, 0, PI, 2);
|
||||
fit->SetLineColor(2);
|
||||
fit->SetLineWidth(2);
|
||||
fit->SetNpx(1000);
|
||||
fit->SetParameter(0, 3000);
|
||||
fit->SetParameter(1, 1);
|
||||
|
||||
gExp->Fit("fit", "q");
|
||||
|
||||
const Double_t * paraE2 = fit->GetParErrors();
|
||||
const Double_t * paraA2 = fit->GetParameters();
|
||||
|
||||
printf("===================== \n");
|
||||
printf("Best fit Amp = %f(%f)\n", paraA2[0], paraE2[0]);
|
||||
printf("Best fit delta = %f(%f) = %f(%f) deg\n", paraA2[1], paraE2[1], atan(paraA2[1]) * 180/PI, atan(paraE2[1])*180/PI);
|
||||
|
||||
c1->cd(1);
|
||||
fit->Draw("same");
|
||||
f1->Draw("same");
|
||||
|
||||
TLatex text;
|
||||
text.SetNDC();
|
||||
text.SetTextFont(82);
|
||||
text.SetTextSize(0.04);
|
||||
|
||||
text.DrawLatex(0.12, 0.85, Form("#delta : %5.1f(%5.1f) deg", atan(paraA2[1]) * 180/PI, atan(paraE2[1])*180/PI));
|
||||
text.DrawLatex(0.12, 0.80, Form("Amp: %5.1f(%5.1f)", paraA2[0], paraE2[0]));
|
||||
|
||||
text.SetTextColor(2);
|
||||
|
||||
text.DrawLatex(0.8, 0.8, Form("%d->%d", Ji, Jf));
|
||||
|
||||
printf("===================== Crtl+C to end.\n");
|
||||
|
||||
theApp.Run();
|
||||
|
||||
return 0;
|
||||
|
||||
}
|
||||
|
||||
//#########################
|
||||
///use for root fit
|
||||
double YE(double * x , double *par){
|
||||
/// x[0] = angle in radian;
|
||||
/// par[0] = a0;
|
||||
/// par[1] = a2;
|
||||
/// par[2] = a4;
|
||||
|
||||
return par[0] + par[1] * LegendreP(2, x[0]) + par[2] * LegendreP(4, x[0]);
|
||||
}
|
||||
|
||||
///use for fit a, delta
|
||||
double Fit(double * x, double *par){
|
||||
/// par[0] = a;
|
||||
/// par[1] = delta;
|
||||
|
||||
double result = 0;
|
||||
for( int k = 0; k <= 4; k += 2){
|
||||
|
||||
result += Q[k]*B[k] * LegendreP(k, x[0]) * ( R1[k] + 2*par[1]*R2[k] + par[1]*par[1]*R3[k] ) / (1 + par[1]*par[1] );
|
||||
}
|
||||
|
||||
return result * par[0];
|
||||
}
|
||||
|
||||
|
||||
double Racah(int j1, int j2, int J, int j3, int j12, int j23){
|
||||
|
||||
return pow(-1, j1+j2+j3+J) * SixJSymbol(j1, j2, j12, j3, J, j23);
|
||||
|
||||
}
|
||||
|
||||
double Rk(int k, int L1, int L2, int J1, int J2){
|
||||
|
||||
return pow(-1, 1+J1-J2+L2-L1-k) * pow((2*J1+1)*(2*L1+1)*(2*L2+1), 0.5) * CGcoeff(k, 0, L1, 1, L2, -1) * Racah(J1, J1, L1, L2, k, J2);
|
||||
|
||||
}
|
||||
|
||||
void PrintRk(int k){
|
||||
|
||||
for( int J1 = 1; J1 < 8; J1 ++){
|
||||
printf("==============================\n");
|
||||
for( int J2 = 0; J2 < 8; J2++){
|
||||
|
||||
int L = abs(J1 - J2);
|
||||
if( L == 0 ) L = 1;
|
||||
printf("%d %d | %10.6f, %10.6f, %10.6f \n", J1, J2, Rk(k, L, L, J1, J2), Rk(k, L, L+1, J1, J2), Rk(k, L+1, L+1, J1, J2));
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double w(int M, int J){
|
||||
double sigma = 1;
|
||||
switch (J) {
|
||||
case 0 : sigma = 0.3989422804014327 ; break;
|
||||
case 1 : sigma = 0.5723377817486753 ; break;
|
||||
case 2 : sigma = 0.7013915463848625 ; break;
|
||||
case 3 : sigma = 0.8091713162791643 ; break;
|
||||
case 4 : sigma = 0.9037290722944527 ; break;
|
||||
case 5 : sigma = 0.9890249035482789 ; break;
|
||||
case 6 : sigma = 1.0673592868302038 ; break;
|
||||
case 7 : sigma = 1.140210831444403 ; break;
|
||||
case 8 : sigma = 1.208599155456379 ; break;
|
||||
case 9 : sigma = 1.273313297516925 ; break;
|
||||
case 10 : sigma = 1.335665551821612 ; break;
|
||||
}
|
||||
|
||||
return 1./sqrt(2 * PI) / sigma * exp( - M*M / 2. / sigma/sigma);
|
||||
|
||||
}
|
||||
|
||||
void Print_w_sum(){ /// Check the normalization of w(m), sum(w(m)) = 1 for all J.
|
||||
|
||||
for( int J = 0; J < 11; J++){
|
||||
double w_sum = 0;
|
||||
for( int m = -J ; m < J+1; m++){
|
||||
w_sum += w(m, J);
|
||||
}
|
||||
printf("%2d %.5f\n", J, w_sum);
|
||||
}
|
||||
}
|
||||
|
||||
double Bk(int k, int J){
|
||||
|
||||
double sum = 0;
|
||||
for(int m = -J; m < J+1 ; m++){
|
||||
sum += w(m, J) * pow(-1, J-m) * pow(2*J+1,0.5) * CGcoeff(k, 0, J, m, J, -m);
|
||||
}
|
||||
|
||||
return sum;
|
||||
|
||||
}
|
||||
129
jSymbol.h
Normal file
129
jSymbol.h
Normal file
|
|
@ -0,0 +1,129 @@
|
|||
#include <stdlib.h>
|
||||
#include <cmath>
|
||||
#include <algorithm>
|
||||
|
||||
using namespace std;
|
||||
|
||||
double factorial(double n){
|
||||
if( n < 0 ) return -100.;
|
||||
return (n == 1. || n == 0.) ? 1. : factorial(n-1) * n ;
|
||||
}
|
||||
|
||||
double CGcoeff(double J, double m, double J1, double m1, double J2, double m2){
|
||||
// (J1,m1) + (J2, m2) = (J, m)
|
||||
|
||||
if( m != m1 + m2 ) return 0;
|
||||
|
||||
double Jmin = abs(J1 - J2);
|
||||
double Jmax = J1+J2;
|
||||
|
||||
if( J < Jmin || Jmax < J ) return 0;
|
||||
|
||||
double s0 = (2*J+1.) * factorial(J+J1-J2) * factorial(J-J1+J2) * factorial(J1+J2-J) / factorial(J+J1+J2 + 1.);
|
||||
s0 = sqrt(s0);
|
||||
|
||||
double s = factorial(J +m ) * factorial(J -m);
|
||||
double s1 = factorial(J1+m1) * factorial(J1-m1);
|
||||
double s2 = factorial(J2+m2) * factorial(J2-m2);
|
||||
s = sqrt(s * s1 * s2);
|
||||
|
||||
//printf(" s0, s = %f , %f \n", s0, s);
|
||||
|
||||
int kMax = min( min( J1+J2-J, J1 - m1), J2 + m2);
|
||||
|
||||
double CG = 0.;
|
||||
for( int k = 0; k <= kMax; k++){
|
||||
double k1 = factorial(J1+J2-J-k);
|
||||
double k2 = factorial(J1-m1-k);
|
||||
double k3 = factorial(J2+m2-k);
|
||||
double k4 = factorial(J - J2 + m1 +k);
|
||||
double k5 = factorial(J - J1 - m2 +k);
|
||||
double temp = pow(-1, k) / (factorial(k) * k1 * k2 * k3 * k4 * k5);
|
||||
if( k1 == -100. || k2 == -100. || k3 == -100. || k4 == -100. || k5 == -100. ) temp = 0.;
|
||||
|
||||
|
||||
//printf(" %d | %.12f | %.12f | %.12f | %.12f | %.12f | %.12f, %f \n", k, k1, k2, k3, k4, k5, temp, factorial(k) );
|
||||
CG += temp;
|
||||
}
|
||||
|
||||
return s0*s*CG;
|
||||
|
||||
}
|
||||
|
||||
double ThreeJSymbol(double J1, double m1, double J2, double m2, double J3, double m3){
|
||||
|
||||
// ( J1 J2 J3 ) = (-1)^(J1-J2 - m3)/ sqrt(2*J3+1) * CGcoeff(J3, -m3, J1, m1, J2, m2)
|
||||
// ( m1 m2 m3 )
|
||||
|
||||
return pow(-1, J1 - J2 - m3)/sqrt(2*J3+1) * CGcoeff(J3, -m3, J1, m1, J2, m2);
|
||||
|
||||
}
|
||||
|
||||
double SixJSymbol(double J1, double J2, double J3, double J4, double J5, double J6){
|
||||
|
||||
// coupling of j1, j2, j3 to J-J1
|
||||
// J1 = j1
|
||||
// J2 = j2
|
||||
// J3 = j12 = j1 + j2
|
||||
// J4 = j3
|
||||
// J5 = J = j1 + j2 + j3
|
||||
// J6 = j23 = j2 + j3
|
||||
|
||||
//check triangle condition
|
||||
if( J3 < abs(J1 - J2 ) || J1 + J2 < J3 ) return 0;
|
||||
if( J6 < abs(J2 - J4 ) || J2 + J4 < J6 ) return 0;
|
||||
if( J5 < abs(J1 - J6 ) || J1 + J6 < J5 ) return 0;
|
||||
if( J5 < abs(J3 - J4 ) || J3 + J4 < J5 ) return 0;
|
||||
|
||||
double sixJ = 0;
|
||||
|
||||
for( float m1 = -J1; m1 <= J1 ; m1 = m1 + 1){
|
||||
for( float m2 = -J2; m2 <= J2 ; m2 = m2 + 1){
|
||||
for( float m3 = -J3; m3 <= J3 ; m3 = m3 + 1){
|
||||
for( float m4 = -J4; m4 <= J4 ; m4 = m4 + 1){
|
||||
for( float m5 = -J5; m5 <= J5 ; m5 = m5 + 1){
|
||||
for( float m6 = -J6; m6 <= J6 ; m6 = m6 + 1){
|
||||
|
||||
double f = (J1 - m1) + (J2 - m2) + (J3 - m3) + (J4 - m4) + (J5 - m5) + (J6 - m6);
|
||||
|
||||
double a1 = ThreeJSymbol( J1, -m1, J2, -m2, J3, -m3); // J3 = j12
|
||||
double a2 = ThreeJSymbol( J1, m1, J5, -m5, J6, m6); // J5 = j1 + (J6 = j23)
|
||||
double a3 = ThreeJSymbol( J4, m4, J2, m2, J6, -m6); // J6 = j23
|
||||
double a4 = ThreeJSymbol( J4, -m4, J5, m5, J3, m3); // J5 = j3 + j12
|
||||
|
||||
double a = a1 * a2 * a3 * a4;
|
||||
//if( a != 0 ) printf("%4.1f %4.1f %4.1f %4.1f %4.1f %4.1f | %f \n", m1, m2, m3, m4, m5, m6, a);
|
||||
|
||||
sixJ += pow(-1, f) * a1 * a2 * a3 * a4;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return sixJ;
|
||||
}
|
||||
|
||||
double NineJSymbol( double J1, double J2, double J3, double J4, double J5, double J6, double J7, double J8, double J9){
|
||||
|
||||
double gMin = min( min (min( abs(J1 - J2 ), abs(J4 - J5)) , abs( J4 - J6 )) , abs(J7 - J8));
|
||||
double gMax = max( max (max( J1+J2, J4+J5), J3+J6), J7+J8);
|
||||
|
||||
//printf(" gMin, gMax = %f %f \n", gMin, gMax);
|
||||
|
||||
double nineJ = 0;
|
||||
for( float g = gMin; g <= gMax ; g = g + 1){
|
||||
double f = pow(-1, 2*g) * (2*g+1);
|
||||
double s1 = SixJSymbol(J1, J4, J7, J8, J9, g);
|
||||
if( s1 == 0 ) continue;
|
||||
double s2 = SixJSymbol(J2, J5, J8, J4, g, J6);
|
||||
if( s2 == 0 ) continue;
|
||||
double s3 = SixJSymbol(J3, J6, J9, g, J1, J2);
|
||||
if( s3 == 0 ) continue;
|
||||
nineJ += f* s1*s2*s3;
|
||||
}
|
||||
|
||||
return nineJ;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user