pderosa888/Angular_Distribution
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

chi not terrible

Browse Source
This commit is contained in:
Peter DeRosa 2022-07-20 16:23:21 -04:00
parent f0214213b2
commit 2ab28226d1
7 changed files with 364 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
.AD.cxx.swp Normal file
View File

Binary file not shown.

7
1147.csv Normal file
View File

@ -0,0 +1,7 @@
angle,Y,Yerr
45,2.33,.1
86,1.84,.1
90,1.73,.1
94,2.11,.1
124,2.53,.1
135,2.41,.1
1 angle Y Yerr
2 45 2.33 .1
3 86 1.84 .1
4 90 1.73 .1
5 94 2.11 .1
6 124 2.53 .1
7 135 2.41 .1

220
AD.cxx
View File

@ -31,12 +31,12 @@ int main(int argc,char **argv){
targetdistance = 4.;
detthickness = 5.;
Energy = 1062.;
Sigma = 0.5;
Energy = 1147.;
Sigma = .5;
j1 = 2.;
j2 = 1.;
j1 = 5.5;
j2 = 3.5;
//--------------------------------
//TEST MODE? y : 1 || n = 0
@ -312,9 +312,9 @@ int main(int argc,char **argv){
}
// cout<<"A0 = "<<residual[0]<<"\n";
cout<<"a2 = "<<residual[1]/residual[0]<<"\n";
cout<<"a4 = "<<residual[2]/residual[0]<<"\n";
cout<<"A0 = "<<residual[0]<<"\n";
cout<<"A2 = "<<residual[1]<<"\n";
cout<<"A4 = "<<residual[2]<<"\n";
@ -485,7 +485,7 @@ int main(int argc,char **argv){
if(isnan(CG2(A))){ // CG2(A) is NaN
cout << "Warning, a CG coefficient returned NaN, set to 0.\n";
// cout << "Warning, a CG coefficient returned NaN, set to 0.\n";
cgg = 0.0;
}else{ //CG2(A) is not NaN.
@ -590,8 +590,30 @@ int main(int argc,char **argv){
double rk02 = rk40datad[nn];
double rk12 = rk41datad[nn];
double rk22 = rk42datad[nn];
//work in progress for possible generation or racah
/*
double CC1[6] = {j1,j1,2,abs(j1-j2),abs(j1-j2),j2};
double CC2[6] = {j1,j1,2,abs(j1-j2),abs(j1-j2)+1,j2};
double CC3[6] = {j1,j1,2,abs(j1-j2)+1,abs(j1-j2)+1,j2};
double rkk1 = SixJ2(CC1);
double rkk2 = SixJ2(CC2);
double rkk3 = SixJ2(CC3);
printf("rk1 = %lf\n",rkk1);
printf("rk2 = %lf\n",rkk2);
printf("rk3 = %lf\n",rkk3);
double J1 = C[0]; // j1
double J2 = C[1]; // j1
double J3 = C[2]; // K
double J4 = C[3]; // Lp = L + 1
double J5 = C[4]; // L = |j1 - j2|
double J6 = C[5]; // j2
*/
// printf("Bk1 = %lf\n",Bk11);
// printf("Bk2 = %lf\n",Bk12);
@ -610,7 +632,7 @@ int main(int argc,char **argv){
double delta_max = 3.14159/2.;
double THETA_min = 0.;
double THETA_max = 3.14159;
double step = 0.0001;
double step = 0.01;
double Tangle = 0;
double delta = 0.;
@ -623,19 +645,20 @@ int main(int argc,char **argv){
double rd0T = 0.;
double rd2T = 0.;
double rd4T = 0.;
double X22,X24 = 0.;
double X2_total = 0.;
int points = (delta_max - delta_min) / step ;
double X2_total;
double YT_tot[points];
// vector<double> YTT;
double a2E = A2E/A0E;
double a4E = A4E/A0E;
int points = (delta_max - delta_min) / step ;
int t_points = (THETA_max - THETA_min) / step;
//int points = (delta_max - delta_min) / step ;
vector<double> chisqr;
vector<double> tdelta;
vector<double> YT_point;
vector<double> YE_point;
double YT0,YT2,YT4,YT,YE = 0.;
double Y_err = 10.;
@ -646,8 +669,50 @@ int main(int argc,char **argv){
delta = tan(atan_delta);
//now sum over all theta;
for(int j = 0; j < dangler.size(); j++){
Tangle = dangler[j];
for(int j = 0; j < points; j++){
Tangle = j*step;
//calculate Y_intensity_theory and and experiment and sum over theta
rd0T = (1 + 2*delta + pow(delta,2))/(1+pow(delta,2));
rd2T = (rk01 + 2*delta*rk11 + pow(delta,2)*rk21)/(1+pow(delta,2));
rd4T = (rk02 + 2*delta*rk12 + pow(delta,2)*rk22)/(1+pow(delta,2));
YT0 = rd0T; // k = 0;
YT2 = QD2*Bk11*rd2T*((1.5*pow(cos(Tangle),2)-.5));
YT4 = QD4*Bk12*rd4T*(35./8.* pow(cos(Tangle),4) - 30./8.*pow(cos(Tangle),2) + 3./8.);
YT = YT0 + YT2 + YT4;
// YT_tot[i] += YT;
YE = 1 + a2E*(1.5*pow(cos(Tangle),2)-.5) + a4E*(35./8.* pow(cos(Tangle),4) - 30./8.*pow(cos(Tangle),2) + 3./8.);
X2_total += pow((YT- YE),2)/((dangler.size()-1)*(pow(Y_err,2)));
// X2_total[i] = X2_total[i] + pow((YE),2)/((dangler.size()-1)*(pow(Y_err,2)));
}
// if(i%27 == 1){
// cout << "X2 = " << X2_total[i] << "\n";
// }
// YTT.push_back(YT_tot[i]);
chisqr.push_back(log(X2_total));
tdelta.push_back(atan_delta);
X2_total = 0.;
// Attempt to take 1 delta value and plot YT;
}
/*
for(int k = 0; k < dangler.size(); k++){
Tangle = dangler[k];
delta = 2.;
//calculate Y_intensity_theory and and experiment and sum over theta
@ -657,28 +722,23 @@ int main(int argc,char **argv){
rd4T = (rk02 + 2*delta*rk12 + pow(delta,2)*rk22)/(1+pow(delta,2));
YT0 = 1.; // P0(cos(theta)) = 1
YT0 = 1; // k = 0;
YT2 = QD2*Bk11*rd2T*((1.5*pow(cos(Tangle),2)-.5));
YT4 = QD4*Bk12*rd4T*(35./8.* pow(cos(Tangle),4) - 30./8.*pow(cos(Tangle),2) + 3./8.);
YT = YT0 + YT2 + YT4;
// cout << "YT = " <<YT << "\n";
// YT_point.push_back(YT);
YE = 1 + a2E*(1.5*pow(cos(Tangle),2)-.5) + a4E*(35./8.* pow(cos(Tangle),4) - 30./8.*pow(cos(Tangle),2) + 3./8.);
YE_point.push_back(YE);
X2_total += pow((YT- YE),2)/(dangler.size()-1)*(pow(Y_err,2));
// X2_total += pow((YT- YE),2)/((dangler.size()-1)*(pow(Y_err,2)));
}
//Figure out denominator.
chisqr.push_back(log(X2_total));
tdelta.push_back(atan_delta);
X2_total = 0.;
}
/*
*/
/*
A2T = QD2*Bk11*rd2T;
a2T = A2T/A0E;
//using the idea that A0E is our normalizer.
@ -690,7 +750,7 @@ int main(int argc,char **argv){
X22 = pow((a2E -a2T),2)/(a2T);
X24 = pow((a4E -a4T),2)/(a4T);
X2_total = (X22 + X24)/2;
*/
*/
vector<double> Theta;
vector<double> AD_I;
double ad_start = 0.;
@ -718,7 +778,7 @@ int main(int argc,char **argv){
double bbb = dydata[i];
dydatas.push_back(bbb/A0E);
printf("Normalized y-intensity %i = %lf\n",i+1,dydatas[i]);
// printf("Normalized y-intensity %i = %lf\n",i+1,dydatas[i]);
}
//This is to temporarily fix the issue with the AD distribution, it takes off the last point in the arrays when plotting or doesnt see them. NEED FIX
@ -739,38 +799,34 @@ int main(int argc,char **argv){
*/
//Initialize Graphics Here.
// cout<< "Dangler size = " << dangler.size() << "\n";
// cout<< "YT point size = "<< YT_point.size() << "\n";
HistoGUI gui;
HistoGUIad gui_ad;
// FIX menu input complexity.
//
// FIX menu input complexity.
//
int optnum = -1;
menu();
printf("Input option : ");
printf("Please enter Input Option : ");
scanf("%d", &optnum);
if(optnum < 0){
optnum = -1;
printf("Negative numbers are not allowed!\nInput option : ");
scanf("%d", &optnum);
}
if(optnum >= 0){
if(optnum == 0){
// menu
if(optnum == 0 and param == 1){
Readme();
optnum = -1;
printf("Input option : ");
printf("Please enter Input Option : ");
scanf("%d", &optnum);
}
//Chi-squared
}else if(optnum == 1 and param == 1){
//if all inputs are valid, plot distribution.
if(param == 1 && optnum == 1){
optnum = -1;
//plotting values here
//x will be arctan of mixing ratio.
//y will be log(X^2);
@ -780,24 +836,59 @@ int main(int argc,char **argv){
gui.Init();
gui.Loop();
gui.Close();
printf("Please enter Input Option : ");
scanf("%d", &optnum);
}else if(optnum == 2 and param == 1){
optnum = -1;
gui_ad.SetData(dangler,dydatas);
gui_ad.SetErrors(deydata);
gui_ad.SetFit(residual[0],residual[1],residual[2]);
gui_ad.Init();
gui_ad.Loop();
gui_ad.Close();
printf("Please enter Input Option : ");
scanf("%d", &optnum);
}else if(param != 1){
cout << "Invalid parameter input, something is wrong\n";
// printf("Input option : ");
// scanf("%d", &optnum);
// if(optnum < 0){
// printf("Input option : ");
// scanf("%d", &optnum);
// }
}
//Plot Angular distribution fit, with points and error bars.
if(param == 1 and optnum == 2){
/*
gui.SetData(Theta,AD_I);
}else{
do{
printf("Negative numbers are not allowed!\nRe-enter Option : ");
scanf("%d", &optnum);
}while(optnum < 0);}
if(optnum >= 0 and optnum != 3){
//repeat 3 cases.
// menu
if(optnum == 0 and param == 1){
Readme();
optnum = -1;
//Chi-squared
}else if(optnum == 1 and param == 1){
optnum = -1;
//plotting values here
//x will be arctan of mixing ratio.
//y will be log(X^2);
gui.SetData(dangler,YT_point);
// gui.SetData(tdelta,chisqr);
//gui.SetData(Theta,AD_I);
gui.Init();
gui.Loop();
gui.Close();
*/
}else if(optnum == 2 and param == 1){
optnum = -1;
gui_ad.SetData(dangler,dydatas);
gui_ad.SetErrors(deydata);
gui_ad.SetFit(residual[0],residual[1],residual[2]);
@ -806,18 +897,15 @@ int main(int argc,char **argv){
gui_ad.Loop();
gui_ad.Close();
}else if(param != 1){
cout << "Invalid parameter input, something is wrong\n";
}
}
//exit
printf("Input option : ");
scanf("%d", &optnum);
if(optnum < 0){
optnum = -1;
printf("Negative numbers are not allowed!\nInput option : ");
scanf("%d", &optnum);
}
if(param == 1 and optnum == 3){

145
Cleb.h
View File

@ -95,6 +95,151 @@ double CG2(double A[6]){
}
return A0*A1*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) * CG(J3, -m3, J1, m1, J2, m2);
}
double TJ2(double A[6]){
// ( J1 J2 J3 ) = (-1)^(J1-J2 - m3)/ sqrt(2*J3+1) * CGcoeff(J3, -m3, J1, m1, J2, m2)
// ( m1 m2 m3 )
double j1 = A[0];
double j2 = A[1];
double J = A[2];
double m1 = A[3];
double m2 = A[4];
double M = A[5];
double B[6];
B[0] = j1;
B[1] = j2;
B[2] = j1 + j2; //j3
B[3] = m1;
B[4] = m2;
B[5] = -(m1 + m2); //m3
return pow(-1, j1 - j2 - B[5])/sqrt(2*B[2]+1) * CG2(B);
}
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 j1, double j1, double j3, double Lp, double L, double j2)
// j3 = j1 + j2, L = |j1 - j2|, Lp = L + 1
//double SixJ2(double J1, double J2, double J3, double J4, double J5, double J6){
double SixJ2(double C[6]){
double J1 = C[0]; // j1
double J2 = C[1]; // j1
double J3 = C[2]; // K
double J4 = C[3]; // Lp = L + 1
double J5 = C[4]; // L = |j1 - j2|
double J6 = C[5]; // j2
// 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 C1[6] = {J2,J3,J1,-m2,-m3,-m1};
double C2[6] = {J5,J6,J1,-m5,m6,m1};
double C3[6] = {J2,J6,J4,m2,-m6,m4};
double C4[6] = {J5,J3,J4,m5,m3,-m4};
double f = (J1 - m1) + (J2 - m2) + (J3 - m3) + (J4 - m4) + (J5 - m5) + (J6 - m6);
double a1 = TJ2(C1); // J3 = j12
double a2 = TJ2(C2); // J5 = j1 + (J6 = j23)
double a3 = TJ2(C3); // J6 = j23
double a4 = TJ2(C4); // 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;
}
/*
int main(int argc, char ** argv){

6
Functlib.h
View File

@ -11,8 +11,8 @@ void menu(){
std::cout <<"==========================================\n";
std::cout<< "\t\t Menu Options \t \n";
std::cout<< "0 - Reading and Instructions \n";
std::cout<< "1 - Plot Chi Sqr \n";
std::cout<< "2 - Plot Ang.Dis Fit (Maybe)\n";
std::cout<< "1 - Plot Chi Sqr vs arc-tan delta\n";
std::cout<< "2 - Plot Angular Distribution \n";
std::cout<< "3 - Exit \n";
std::cout <<"==========================================\n";
std::cout<< " \n";
@ -29,7 +29,7 @@ void Readme(){
std::cout<< " \n";
std::cout<< " Follow the prompt in order to correctly display \n";
std::cout<< " \n";
std::cout<< " To close the gui, press any button. \n";
std::cout<< " To close the gui, press most buttons. \n";
std::cout<< " To zoom in, left click then drag and let go. To unzoom\n";
std::cout<< " press the space bar. To draw, right click\n";
std::cout<< " \n";

6
GUI_AD.h
View File

@ -397,7 +397,8 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
for(int i=0; i < x.size() - 1; i++){
x_wid = (x[i] + x_offset) / width_scale;
y_wid = (y[i] + y_offset) / height_scale;
y_errors_wid = ((y_errors[i]/A0) / height_scale);
// y_errors_wid = ((y_errors[i]/A0) / height_scale);
y_errors_wid = ((y_errors[i]) / height_scale);
//draws the point
XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -3, y_wid -3, 6, 6);
@ -460,7 +461,8 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
x_wid = (x[i] + x_offset) / width_scale;
y_wid = (y[i] + y_offset) / height_scale;
y_errors_wid = ((y_errors[i]/A0) / height_scale);
// y_errors_wid = ((y_errors[i]/A0) / height_scale);
y_errors_wid = ((y_errors[i]) / height_scale);
//draws the point
XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -3, y_wid -3, 6, 6);

8
ad.txt
View File

@ -1,7 +1,7 @@
Radius = 3 [cm]
Distance = 4 [cm]
Thickness = 5 [cm]
Atten.C = 0.294297 [cm^-1]
Gamma_E = 1062 [KeV]
QD2 = 0.801748
QD4 = 0.447498
Atten.C = 0.282206 [cm^-1]
Gamma_E = 1147 [KeV]
QD2 = 0.802565
QD4 = 0.44946