diff --git a/.AD.cxx.swp b/.AD.cxx.swp new file mode 100644 index 0000000..350621b Binary files /dev/null and b/.AD.cxx.swp differ diff --git a/1147.csv b/1147.csv new file mode 100644 index 0000000..7616b1b --- /dev/null +++ b/1147.csv @@ -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 diff --git a/AD.cxx b/AD.cxx index 33bf8bf..8e3afaf 100644 --- a/AD.cxx +++ b/AD.cxx @@ -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 @@ -99,9 +99,9 @@ int main(int argc,char **argv){ } double QD2 = QK2(Energy, detradius, targetdistance, detthickness); ; - + double QD4 = QK4(Energy, detradius, targetdistance, detthickness); ; - + // cout << "QD2 = " << " " << QD2 << " QD4 = " << " " << QD4 << "\n"; //input data file of Angular Data (theta, Yexp, Yerr) @@ -226,7 +226,7 @@ int main(int argc,char **argv){ // Note, the y data must be scaled by sin(theta) of the given angle. double a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12 = 0.; - + for(int i = 0; i< dangler.size(); i++){ //eq1 @@ -312,9 +312,9 @@ int main(int argc,char **argv){ } - // cout<<"A0 = "< 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 chisqr; vector tdelta; - + vector YT_point; + vector YE_point; double YT0,YT2,YT4,YT,YE = 0.; double Y_err = 10.; @@ -644,53 +667,90 @@ int main(int argc,char **argv){ for(int i = 0; i< points; i++){ atan_delta = i*step + delta_min; delta = tan(atan_delta); - //now sum over all theta; - - for(int j = 0; j < dangler.size(); j++){ - Tangle = dangler[j]; + //now sum over all theta; - //calculate Y_intensity_theory and and experiment and sum over theta + 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)); - // 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 = 1.; // P0(cos(theta)) = 1 - + + 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; - 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)); - } - //Figure out denominator. + // 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; + + } -/* - A2T = QD2*Bk11*rd2T; - a2T = A2T/A0E; - //using the idea that A0E is our normalizer. + /* + for(int k = 0; k < dangler.size(); k++){ + Tangle = dangler[k]; + delta = 2.; - rd4T = (rk02 + 2*delta*rk12 + pow(delta,2)*rk22)/(1+pow(delta,2)); - A4T = QD4*Bk12*rd4T; - a4T = A4T/A0E; + //calculate Y_intensity_theory and and experiment and sum over theta - X22 = pow((a2E -a2T),2)/(a2T); - X24 = pow((a4E -a4T),2)/(a4T); - X2_total = (X22 + X24)/2; -*/ + // 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 = 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 = " < Theta; vector 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,85 +799,113 @@ 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){ - Readme(); - optnum = -1; - - printf("Input option : "); - scanf("%d", &optnum); - } - - //if all inputs are valid, plot distribution. - if(param == 1 && optnum == 1){ - optnum = -1; + if(optnum >= 0){ - //plotting values here - //x will be arctan of mixing ratio. - //y will be log(X^2); + // menu + if(optnum == 0 and param == 1){ + Readme(); + optnum = -1; - gui.SetData(tdelta,chisqr); - //gui.SetData(Theta,AD_I); - gui.Init(); - gui.Loop(); - gui.Close(); + printf("Please enter Input Option : "); + scanf("%d", &optnum); + //Chi-squared + }else if(optnum == 1 and param == 1){ - // printf("Input option : "); - // scanf("%d", &optnum); - // if(optnum < 0){ - // printf("Input option : "); - // scanf("%d", &optnum); - // } - } + optnum = -1; + //plotting values here + //x will be arctan of mixing ratio. + //y will be log(X^2); - //Plot Angular distribution fit, with points and error bars. - if(param == 1 and optnum == 2){ - /* - gui.SetData(Theta,AD_I); - gui.Init(); - gui.Loop(); - gui.Close(); - */ + gui.SetData(tdelta,chisqr); + //gui.SetData(Theta,AD_I); + gui.Init(); + gui.Loop(); + gui.Close(); + printf("Please enter Input Option : "); + scanf("%d", &optnum); + }else if(optnum == 2 and param == 1){ - gui_ad.SetData(dangler,dydatas); - gui_ad.SetErrors(deydata); - gui_ad.SetFit(residual[0],residual[1],residual[2]); + 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(); + 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"; + + } + + + }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]); + + gui_ad.Init(); + 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){ diff --git a/Cleb.h b/Cleb.h index 63f4983..9d598dd 100644 --- a/Cleb.h +++ b/Cleb.h @@ -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){ diff --git a/Functlib.h b/Functlib.h index 6a02ab6..75006d0 100644 --- a/Functlib.h +++ b/Functlib.h @@ -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"; diff --git a/GUI_AD.h b/GUI_AD.h index 3377f4d..2cff97a 100644 --- a/GUI_AD.h +++ b/GUI_AD.h @@ -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); diff --git a/ad.txt b/ad.txt index 7d25f53..5cdeadd 100644 --- a/ad.txt +++ b/ad.txt @@ -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