diff --git a/.Chi_error.h.swp b/.Chi_error.h.swp deleted file mode 100644 index c1253ac..0000000 Binary files a/.Chi_error.h.swp and /dev/null differ diff --git a/.gitignore b/.gitignore index 6f13b0f..8abe69e 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,3 @@ ADEX +ADT +.gitignore diff --git a/1062b.csv b/1062b.csv deleted file mode 100644 index 7056ee1..0000000 --- a/1062b.csv +++ /dev/null @@ -1,4 +0,0 @@ -angle,Y,Yerr -45,1.29,.10 -90,1.10,.10 -135,1.29,.10 diff --git a/AD.cxx b/AD.cxx index 93f726e..9a1d762 100644 --- a/AD.cxx +++ b/AD.cxx @@ -1,6 +1,6 @@ #include "global.h" -//#include "GUI.h" +#include "GUI_AD.h" #include "GUI_Base.h" #include "Functlib.h" #include "QDK.h" @@ -15,730 +15,781 @@ using namespace std; int main(int argc,char **argv){ - HistoGUI gui; - - double j1, j2; + + double j1, j2; - //detector data input - double detradius, targetdistance, detthickness = -1; - double Energy; - //input .csv file - double Sigma; + //detector data input + double detradius, targetdistance, detthickness = -1; + double Energy; + //input .csv file + double Sigma; -//-------------------------------- + //-------------------------------- - detradius = 3.; - targetdistance = 4.; - detthickness = 5.; - - Energy = 1062.; - Sigma = 0.; - + detradius = 3.; + targetdistance = 4.; + detthickness = 5.; - j1 = 2.; - j2 = 1.; -//-------------------------------- - -//TEST MODE? y : 1 || n = 0 - - int test = 1; - - int det_param_token = -1; - int gamma_energy_token = -1; - int ang_file_token = -1; - int sigma_token = -1; - - int j1j2token = -1; - if(test == 1){ - - det_param_token = 1; - gamma_energy_token = 1; - ang_file_token = -1; - sigma_token = 1; - - j1j2token = 1; - -} -if( test == 0 ) { - printf("Input detector radius, target distance, and detector thickness : "); - scanf("%lf,%lf,%lf",&detradius,&targetdistance,&detthickness); + Energy = 1062.; + Sigma = 0.; - if(detradius > 0 && targetdistance > 0 && detthickness > 0){ - det_param_token = 1; - }else{ + j1 = 2.; + j2 = 1.; + //-------------------------------- - do{ - if(detradius< 0|| targetdistance < 0|| detthickness < 0){ - printf("Negative numbers are not allowed!\nRe-enter radius, distance, and thickness : "); - scanf("%lf,%lf,%lf", &detradius,&targetdistance,&detthickness); - } - }while(detradius< 0|| targetdistance < 0|| detthickness < 0);} + //TEST MODE? y : 1 || n = 0 - if(detradius > 0 && targetdistance > 0 && detthickness > 0){ - det_param_token = 1; printf(" --- Detector characteristics loaded --- \n"); - } + int test = 1; - printf("Enter the gamma-ray energy : "); - scanf("%lf", &Energy); - if(Energy > 0){ - gamma_energy_token = 1; - }else{ - do{ - if(Energy <= 0){ - printf("Negative numbers are not allowed!\nRe-enter gamma-ray energy : "); - scanf("%lf", &Energy);} + int det_param_token = -1; + int gamma_energy_token = -1; + int ang_file_token = -1; + int sigma_token = -1; - }while(Energy <= 0);} + int j1j2token = -1; + if(test == 1){ - if(Energy > 0){ - gamma_energy_token = 1; printf(" --- Gamma Energy loaded --- \n"); - } + det_param_token = 1; + gamma_energy_token = 1; + ang_file_token = -1; + sigma_token = 1; - //then calucate QD2 and QD4, and replace the 0 with them. -} + j1j2token = 1; - double QD2 = QK2(Energy, detradius, targetdistance, detthickness); ; - //.82;//NEED TO READ FROM ad.txt in the future. - double QD4 = QK4(Energy, detradius, targetdistance, detthickness); ; - //.44; -// cout << "QD2 = " << " " << QD2 << " QD4 = " << " " << QD4 << "\n"; - -//input data file of Angular Data (theta, Yexp, Yerr) - - vector > content; - vector row; - string line, word; - - vector adata; - vector ydata; - vector eydata; - int num = 0; - string fname,fname2; - cout<<"Enter the file name : "; - cin>>fname; - - int aaaa = -1; - fstream file(fname.c_str()); - - if(file.is_open()){ - cout<< " --- Angular Data File Loaded --- \n"; - ang_file_token = 1; - while(getline(file,line)){ - - row.clear(); - - stringstream str(line); - - while(getline(str, word, ',')) - - row.push_back(word); - content.push_back(row); - adata.push_back(row[0]); - ydata.push_back(row[1]); - eydata.push_back(row[2]); - - } - }else{ - do{ - cout<< "Could not open the file\n"; - cout<< "Enter the file name : "; - cin>>fname2; - - fstream file(fname2.c_str()); - - if(file.is_open()){ - aaaa = 1; - - while(getline(file,line)){ - - row.clear(); - - stringstream str(line); - - while(getline(str, word, ',')) - - row.push_back(word); - content.push_back(row); - adata.push_back(row[0]); - ydata.push_back(row[1]); - eydata.push_back(row[2]); - - } - - }else{ aaaa = -1; } - - }while(aaaa == -1); - - cout<< " --- Angular Data File Loaded --- \n"; - ang_file_token = 1; - - } - - for(int i = 0; i< adata.size();i++){ -// cout << adata[i] << "\n"; -// cout << ydata[i] << "\n"; -// cout << eydata[i] << "\n"; - } - -// Legendre Polinomial fit using adata, ydata, eydata; - - adata.erase(adata.begin()); - ydata.erase(ydata.begin()); - eydata.erase(eydata.begin()); - - vector dangle = string_to_double_vector(adata); - vector dydata = string_to_double_vector(ydata); - vector deydata = string_to_double_vector(eydata); - - vector dangler; - - for(int i = 0; i < dangle.size(); i++){ - double aa = dangle[i]; - - dangler.push_back(aa*3.14159/180.); - printf("dangle = %lf\n",dangler[i]); - } - - /* - for(int i=0; i 0 && targetdistance > 0 && detthickness > 0){ + det_param_token = 1; + }else{ + + do{ + if(detradius< 0|| targetdistance < 0|| detthickness < 0){ + printf("Negative numbers are not allowed!\nRe-enter radius, distance, and thickness : "); + scanf("%lf,%lf,%lf", &detradius,&targetdistance,&detthickness); + } + }while(detradius< 0|| targetdistance < 0|| detthickness < 0);} + + if(detradius > 0 && targetdistance > 0 && detthickness > 0){ + det_param_token = 1; printf(" --- Detector characteristics loaded --- \n"); + } + + printf("Enter the gamma-ray energy : "); + scanf("%lf", &Energy); + if(Energy > 0){ + gamma_energy_token = 1; + }else{ + do{ + if(Energy <= 0){ + printf("Negative numbers are not allowed!\nRe-enter gamma-ray energy : "); + scanf("%lf", &Energy);} + + }while(Energy <= 0);} + + if(Energy > 0){ + gamma_energy_token = 1; printf(" --- Gamma Energy loaded --- \n"); + } + + //then calucate QD2 and QD4, and replace the 0 with them. + } + + double QD2 = QK2(Energy, detradius, targetdistance, detthickness); ; + //.82;//NEED TO READ FROM ad.txt in the future. + double QD4 = QK4(Energy, detradius, targetdistance, detthickness); ; + //.44; + // cout << "QD2 = " << " " << QD2 << " QD4 = " << " " << QD4 << "\n"; + + //input data file of Angular Data (theta, Yexp, Yerr) + + vector > content; + vector row; + string line, word; + + vector adata; + vector ydata; + vector eydata; + int num = 0; + string fname,fname2; + cout<<"Enter the file name : "; + cin>>fname; + + int aaaa = -1; + fstream file(fname.c_str()); + + if(file.is_open()){ + cout<< " --- Angular Data File Loaded --- \n"; + ang_file_token = 1; + while(getline(file,line)){ + + row.clear(); + + stringstream str(line); + + while(getline(str, word, ',')) + + row.push_back(word); + content.push_back(row); + adata.push_back(row[0]); + ydata.push_back(row[1]); + eydata.push_back(row[2]); + + } + }else{ + do{ + cout<< "Could not open the file\n"; + cout<< "Enter the file name : "; + cin>>fname2; + + fstream file(fname2.c_str()); + + if(file.is_open()){ + aaaa = 1; + + while(getline(file,line)){ + + row.clear(); + + stringstream str(line); + + while(getline(str, word, ',')) + + row.push_back(word); + content.push_back(row); + adata.push_back(row[0]); + ydata.push_back(row[1]); + eydata.push_back(row[2]); + + } + + }else{ aaaa = -1; } + + }while(aaaa == -1); + + cout<< " --- Angular Data File Loaded --- \n"; + ang_file_token = 1; + + } + + for(int i = 0; i< adata.size();i++){ + // cout << adata[i] << "\n"; + // cout << ydata[i] << "\n"; + // cout << eydata[i] << "\n"; + } + + // Legendre Polinomial fit using adata, ydata, eydata; + + adata.erase(adata.begin()); + ydata.erase(ydata.begin()); + eydata.erase(eydata.begin()); + + vector dangle = string_to_double_vector(adata); + vector dydata = string_to_double_vector(ydata); + vector deydata = string_to_double_vector(eydata); + + vector dangler; + + for(int i = 0; i < dangle.size(); i++){ + double aa = dangle[i]; + + dangler.push_back(aa*3.14159/180.); + printf("dangle = %lf\n",dangler[i]); + } + // if you need to scale the y data by sin(theta) + + + /* + for(int i=0; i= 0; i--){ + + residual[i] = mat[i][n]; + + for( int j = i + 1; j= 0 ){ + sigma_token = 1; + }else{ + do{ + + printf("Negative numbers are not allowed!\nRe-enter Sigma : "); + scanf("%lf", &Sigma); + }while(Sigma <0);} + + if(Sigma >= 0 ){ + sigma_token = 1; printf(" --- Sigma Loaded --- \n"); } - } - } - } -//performing gaussian elimination. + printf("Please enter J1,J2 : "); + scanf("%lf,%lf",&j1,&j2); + if(j1 >= 0 && j2 >= 0 ){ + j1j2token = 1; + }else{ + do{ + printf("Negative numbers are not allowed!\nRe-enter J1,J2 : "); + scanf("%lf,%lf", &j1,&j2); + }while(j1 < 0 || j2 < 0);} - for(int i = 0; i< n-1; i++) { - - for(int j = i +1; j= 0; i--){ - - residual[i] = mat[i][n]; - - for( int j = i + 1; j= 0 ){ - sigma_token = 1; - }else{ - do{ - - printf("Negative numbers are not allowed!\nRe-enter Sigma : "); - scanf("%lf", &Sigma); - }while(Sigma <0);} - - if(Sigma >= 0 ){ - sigma_token = 1; printf(" --- Sigma Loaded --- \n"); - } - - printf("Please enter J1,J2 : "); - scanf("%lf,%lf",&j1,&j2); - if(j1 >= 0 && j2 >= 0 ){ - j1j2token = 1; - }else{ - do{ - printf("Negative numbers are not allowed!\nRe-enter J1,J2 : "); - scanf("%lf,%lf", &j1,&j2); - }while(j1 < 0 || j2 < 0);} - - if(j1 >= 0 && j2 >= 0 ){ - j1j2token = 1; printf(" --- J1 & J2 Loaded --- \n"); - } - -} - -// Calculate BK(j1) for perfect allignment -// - - double A[6]; - double js = sqrt(2*j1 +1); - double j12 = 2*j1; - int IS = (int)j12 % 2; - - double cg12,cg22,cg24,cg14 = 0.; - double I1, I2 = 0.; - - double cg1, cg2 = 0.; - - double Bk11, Bk12 = 0.; - if(Sigma == 0){ - if(IS == 1){ - A[0] = j1; - A[1] = j1; - A[2] = 2.; - A[3] = .5; - A[4] = -.5; - A[5] = 0.; - - cg12 = CG2(A); - - A[3] = -.5; - A[4] = .5; - - cg22 = CG2(A); - - A[2] = 4.; - - cg24 = CG2(A); - - A[3] = .5; - A[4] = -.5; - - cg14 = CG2(A); - - - printf("cg12 = %lf\n",cg12); - printf("cg22 = %lf\n",cg22); - printf("cg24 = %lf\n",cg24); - printf("cg24 = %lf\n",cg24); - - printf("Bk11 = %lf\n",Bk11); - printf("Bk12 = %lf\n",Bk12); - - Bk11 = (0.5)*js*(pow(-1,I1) * cg12 + pow(-1,I2) * cg22); - Bk12 = (0.5)*js*(pow(-1,I1) * cg14 + pow(-1,I2) * cg24); - }else if(IS == 0){ - A[0] = j1; - A[1] = j1; - A[2] = 2.; - A[3] = 0.; - A[4] = 0.; - A[5] = 0.; - - cg1 = CG2(A); - - Bk11 = pow(-1,j1) * js * cg1; - - A[2] = 4.; - - cg2 = CG2(A); - Bk12 = pow(-1,j1) * js * cg2; - - printf("Bk11 = %lf\n",Bk11); - printf("Bk12 = %lf\n",Bk12); - printf("cg1 = %lf\n",cg1); - printf("cg2 = %lf\n",cg2); - } - } -//normalize W(m1) -// - - j12 = 2*j1; - double j14 = 4*j1; - - double sigsq = pow(Sigma,2); - double sum1 = 0.; - - double am1,amsq,x,ex = 0.; -//-------------------------------------------------------SOMETHING WRONG WITH CN1. - for(int i = 0; i <= j14; i = i + 2){ - am1 = 0.5*(i - j12); - amsq = pow(am1,2); - x = - (amsq/(2*sigsq)); - ex = exp(x); - - sum1 = sum1 + ex; - } - double cn1 = 1./sum1; - -// cout << "cn1" << " " << cn1 << "\n"; - - double AL0 = 0.; - double A0 = j1 - j2; - double A0b = abs(A0); - if(A0b > 0){ AL0 = A0b;} - - if(A0b <= 0){AL0 = 1.;} - - double AL1 = AL0 +1; - - double am11,amsq1,x1,ex1 = 0.; -//calculate Bk(j1) for gaussian W(m1) or non zero Sigma - -//NEEDS VALUE FIX - if(Sigma != 0){ - double A[6]; - A[0] = j1; - A[1] = j1; - - double sfact = sqrt(2*j1 +1); - double cgg = 0.; - double tTerm = 0.; - double II = 0.; - for(int i = 2; i<= 4; i = i + 2){ - - A[2] = i; - Bk11 = 0.; - Bk12 = 0.; - for(int m = 0; m <= j14; m = m +2){ - am11 = 0.5*(m - j12); - amsq1 = pow(am11,2); - x1 = -(amsq/(2*sigsq)); - II = j1 - am11; - A[3] = am11; - A[4] = -am11; - A[5] = 0; - cgg = CG2(A); - // cout << "cgg = " << " " << cgg << "\n"; - // cout << "am11 = " << " " << am11 << "\n"; - // cout << "amsq1 = " << " " << amsq1 << "\n"; - // cout << "x1 = " << " " << x1 << "\n"; - - ex1 = exp(x1); - tTerm = cn1*ex1*pow(-1,II) * sfact*cgg; - - // cout <<"tTerm "<<" " << tTerm << "\n"; - - if(i == 2) Bk11 = Bk11 + tTerm; - if(i == 4) Bk12 = Bk12 + tTerm; - } - } - cout << "Bk11" << " " << Bk11 << "\n"; - cout << "Bk12" << " " << Bk12 << "\n"; - - - } - - vector > content_r; - vector row_r; - string line_r, word_r; - - vector j1data; - vector j2data; - vector Rk20data; - vector Rk21data; - vector Rk22data; - vector Rk40data; - vector Rk41data; - vector Rk42data; - - - fstream file1("RkTable.csv"); - if(file1.is_open()){ - while(getline(file1,line_r)){ - row_r.clear(); - - stringstream str(line_r); - - while(getline(str, word_r, ',')) - - row_r.push_back(word_r); - content_r.push_back(row_r); - j1data.push_back(row_r[0]); - j2data.push_back(row_r[1]); - Rk20data.push_back(row_r[2]); - Rk21data.push_back(row_r[3]); - Rk22data.push_back(row_r[4]); - Rk40data.push_back(row_r[5]); - Rk41data.push_back(row_r[6]); - Rk42data.push_back(row_r[7]); - - } - cout<< " --- Rk(llj1j2) Table Loaded --- \n"; - }else cout<< "Could not open the file\n"; - - vector::iterator dit; - vector::iterator dit2; - - vector j1datad = string_to_double_vector(j1data); - vector j2datad = string_to_double_vector(j2data); - vector rk20datad = string_to_double_vector(Rk20data); - vector rk21datad = string_to_double_vector(Rk21data); - vector rk22datad = string_to_double_vector(Rk22data); - vector rk40datad = string_to_double_vector(Rk40data); - vector rk41datad = string_to_double_vector(Rk41data); - vector rk42datad = string_to_double_vector(Rk42data); - - int nn = 0; - for(dit = j1datad.begin();dit < j1datad.end();dit++){ - - if(j1 == *dit){ - - for(dit2 = j2datad.begin();dit2 < j2datad.end();dit2++){ - - if(j2 == *dit2 and j1 == j1datad[distance(j2datad.begin(),dit2)]){ - // cout << *dit << " " << *dit2 << " Index Equals = " <= 0 && j2 >= 0 ){ + j1j2token = 1; printf(" --- J1 & J2 Loaded --- \n"); } - } + } - } - -// cout << j1<< " " << j2 << " "<< nn << "\n"; - - -// rk01,rk11,rk21, rk02, rk12,rk22 //READ OUT OF Rk TABLE - - double rk01 = rk20datad[nn]; - double rk11 = rk21datad[nn]; - double rk21 = rk22datad[nn]; - double rk02 = rk40datad[nn]; - double rk12 = rk41datad[nn]; - double rk22 = rk42datad[nn]; + + // Calculate BK(j1) for perfect allignment + // + + double A[6]; + double js = sqrt(2*j1 +1); + double j12 = 2*j1; + int IS = (int)j12 % 2; + + double cg12,cg22,cg24,cg14 = 0.; + double I1, I2 = 0.; + + double cg1, cg2 = 0.; + + double Bk11, Bk12 = 0.; + if(Sigma == 0){ + if(IS == 1){ + A[0] = j1; + A[1] = j1; + A[2] = 2.; + A[3] = .5; + A[4] = -.5; + A[5] = 0.; + + cg12 = CG2(A); + + A[3] = -.5; + A[4] = .5; + + cg22 = CG2(A); + + A[2] = 4.; + + cg24 = CG2(A); + + A[3] = .5; + A[4] = -.5; + + cg14 = CG2(A); -// printf("Bk1 = %lf\n",Bk11); -// printf("Bk2 = %lf\n",Bk12); + printf("cg12 = %lf\n",cg12); + printf("cg22 = %lf\n",cg22); + printf("cg24 = %lf\n",cg24); + printf("cg24 = %lf\n",cg24); - printf("rk01 = %lf\n",rk01); - printf("rk11 = %lf\n",rk11); - printf("rk21 = %lf\n",rk21); - printf("rk02 = %lf\n",rk02); - printf("rk12 = %lf\n",rk12); - printf("rk22 = %lf\n",rk22); - + printf("Bk11 = %lf\n",Bk11); + printf("Bk12 = %lf\n",Bk12); - // here the Chi-sqs from A2 and A4 needed to be added together. + Bk11 = (0.5)*js*(pow(-1,I1) * cg12 + pow(-1,I2) * cg22); + Bk12 = (0.5)*js*(pow(-1,I1) * cg14 + pow(-1,I2) * cg24); + }else if(IS == 0){ + A[0] = j1; + A[1] = j1; + A[2] = 2.; + A[3] = 0.; + A[4] = 0.; + A[5] = 0.; - - double delta_min = -3.14159/2.; - double delta_max = 3.14159/2.; - double step = 0.001; + cg1 = CG2(A); - double delta = 0.; - double tan_delta =0.; - double A0E = residual[0];//NEED FrOM AF FIT - double A2E = residual[1]; - double A4E = residual[2]; - double A2T,a2T = 0.; - double A4T,a4T = 0.; - double rd2T = 0.; - double rd4T = 0.; - double X22,X24 = 0.; - double X2_total = 0.; + Bk11 = pow(-1,j1) * js * cg1; - double a2E = A2E/A0E; - double a4E = A4E/A0E; + A[2] = 4.; + + cg2 = CG2(A); + Bk12 = pow(-1,j1) * js * cg2; + + printf("Bk11 = %lf\n",Bk11); + printf("Bk12 = %lf\n",Bk12); + printf("cg1 = %lf\n",cg1); + printf("cg2 = %lf\n",cg2); + } + } + //normalize W(m1) + // + + j12 = 2*j1; + double j14 = 4*j1; + + double sigsq = pow(Sigma,2); + double sum1 = 0.; + + double am1,amsq,x,ex = 0.; + //-------------------------------------------------------SOMETHING WRONG WITH CN1. + for(int i = 0; i <= j14; i = i + 2){ + am1 = 0.5*(i - j12); + amsq = pow(am1,2); + x = - (amsq/(2*sigsq)); + ex = exp(x); + + sum1 = sum1 + ex; + } + double cn1 = 1./sum1; + + // cout << "cn1" << " " << cn1 << "\n"; + + double AL0 = 0.; + double A0 = j1 - j2; + double A0b = abs(A0); + if(A0b > 0){ AL0 = A0b;} + + if(A0b <= 0){AL0 = 1.;} + + double AL1 = AL0 +1; + + double am11,amsq1,x1,ex1 = 0.; + //calculate Bk(j1) for gaussian W(m1) or non zero Sigma + + //NEEDS VALUE FIX + if(Sigma != 0){ + double A[6]; + A[0] = j1; + A[1] = j1; + + double sfact = sqrt(2*j1 +1); + double cgg = 0.; + double tTerm = 0.; + double II = 0.; + for(int i = 2; i<= 4; i = i + 2){ + + A[2] = i; + Bk11 = 0.; + Bk12 = 0.; + for(int m = 0; m <= j14; m = m +2){ + am11 = 0.5*(m - j12); + amsq1 = pow(am11,2); + x1 = -(amsq/(2*sigsq)); + II = j1 - am11; + A[3] = am11; + A[4] = -am11; + A[5] = 0; + cgg = CG2(A); + // cout << "cgg = " << " " << cgg << "\n"; + // cout << "am11 = " << " " << am11 << "\n"; + // cout << "amsq1 = " << " " << amsq1 << "\n"; + // cout << "x1 = " << " " << x1 << "\n"; + + ex1 = exp(x1); + tTerm = cn1*ex1*pow(-1,II) * sfact*cgg; + + // cout <<"tTerm "<<" " << tTerm << "\n"; + + if(i == 2) Bk11 = Bk11 + tTerm; + if(i == 4) Bk12 = Bk12 + tTerm; + } + } + cout << "Bk11" << " " << Bk11 << "\n"; + cout << "Bk12" << " " << Bk12 << "\n"; - int points = (delta_max - delta_min) / step ; - - vector chisqr; - vector tdelta; + } - for(int i = 0; i< points; i++){ - tan_delta = i*step + delta_min; - delta = tan(tan_delta); - rd2T = (rk01 + 2*delta*rk11 + pow(delta,2)*rk21)/(1+pow(delta,2)); - A2T = QD2*Bk11*rd2T; - a2T = A2T/A0E; - //using the idea that A0E is our normalizer. - - rd4T = (rk02 + 2*delta*rk12 + pow(delta,2)*rk22)/(1+pow(delta,2)); - A4T = QD4*Bk12*rd4T; - a4T = A4T/A0E; - - X22 = pow(abs(a2E -a2T),2)/(abs(a2T)); - X24 = pow(abs(a4E -a4T),2)/(abs(a4T)); - X2_total = (X22 + X24)/2; + vector > content_r; + vector row_r; + string line_r, word_r; - chisqr.push_back(-log(X2_total)); - tdelta.push_back(tan_delta); - - } -//NEEDS FIX - -//If gui does this how does the vector pushback here make sense. - vector Theta; - vector AD_I; - double ad_start = -3.1415/2; - int adpoints = (3.1415)/(step); - double theta,Iad = 0.; - double aaa,aab,aac = 0.; - for(int i = 0; i < adpoints; i++){ - - theta = i*step + ad_start; - aaa = A0E; - aab = A2E*(1.5 * pow(theta,2) - .5); - aac = A4E*(35./8. * pow(theta,4) - 30./8. * pow(theta,2) + 3./8. ); - - Iad = aaa + aab + aac; - - AD_I.push_back(Iad); - Theta.push_back(theta); - - } + vector j1data; + vector j2data; + vector Rk20data; + vector Rk21data; + vector Rk22data; + vector Rk40data; + vector Rk41data; + vector Rk42data; + fstream file1("RkTable.csv"); + if(file1.is_open()){ + while(getline(file1,line_r)){ + row_r.clear(); - int param; + stringstream str(line_r); - param = param_run(det_param_token, gamma_energy_token, ang_file_token, sigma_token, j1j2token); -/* - cout << det_param_token << "\n"; - cout << gamma_energy_token << "\n"; - cout << ang_file_token << "\n"; - cout << sigma_token << "\n"; - cout << j1j2token << "\n"; - - cout << param << "\n"; -*/ - int optnum = -1; - menu(); + while(getline(str, word_r, ',')) - printf("Input option : "); - scanf("%d", &optnum); - if(optnum < 0){ - optnum = -1; - printf("Negative numbers are not allowed!\nInput option : "); - scanf("%d", &optnum); - } + row_r.push_back(word_r); + content_r.push_back(row_r); + j1data.push_back(row_r[0]); + j2data.push_back(row_r[1]); + Rk20data.push_back(row_r[2]); + Rk21data.push_back(row_r[3]); + Rk22data.push_back(row_r[4]); + Rk40data.push_back(row_r[5]); + Rk41data.push_back(row_r[6]); + Rk42data.push_back(row_r[7]); - if(optnum == 0){ - Readme(); - optnum = -1; + } + cout<< " --- Rk(llj1j2) Table Loaded --- \n"; + }else cout<< "Could not open the file\n"; + + vector::iterator dit; + vector::iterator dit2; + + vector j1datad = string_to_double_vector(j1data); + vector j2datad = string_to_double_vector(j2data); + vector rk20datad = string_to_double_vector(Rk20data); + vector rk21datad = string_to_double_vector(Rk21data); + vector rk22datad = string_to_double_vector(Rk22data); + vector rk40datad = string_to_double_vector(Rk40data); + vector rk41datad = string_to_double_vector(Rk41data); + vector rk42datad = string_to_double_vector(Rk42data); + + int nn = 0; + for(dit = j1datad.begin();dit < j1datad.end();dit++){ + + if(j1 == *dit){ + + for(dit2 = j2datad.begin();dit2 < j2datad.end();dit2++){ + + if(j2 == *dit2 and j1 == j1datad[distance(j2datad.begin(),dit2)]){ + // cout << *dit << " " << *dit2 << " Index Equals = " < chisqr; + vector tdelta; + + for(int i = 0; i< points; i++){ + atan_delta = i*step + delta_min; + delta = tan(atan_delta); + rd2T = (rk01 + 2*delta*rk11 + pow(delta,2)*rk21)/(1+pow(delta,2)); + A2T = QD2*Bk11*rd2T; + a2T = A2T/A0E; + //using the idea that A0E is our normalizer. + + rd4T = (rk02 + 2*delta*rk12 + pow(delta,2)*rk22)/(1+pow(delta,2)); + A4T = QD4*Bk12*rd4T; + a4T = A4T/A0E; + + X22 = pow(abs(a2E -a2T),2)/(abs(a2T)); + X24 = pow(abs(a4E -a4T),2)/(abs(a4T)); + X2_total = (X22 + X24)/2; + + chisqr.push_back(-log(X2_total)); + tdelta.push_back(atan_delta); + + } + + vector Theta; + vector AD_I; + double ad_start = 0.; + int adpoints = (3.1415)/(step); + double theta,Iad = 0.; + double aaa,aab,aac = 0.; + for(int i = 0; i < adpoints; i++){ + + theta = i*step + ad_start; + aaa = 1.; + aab = (A2E/A0E)*(1.5 * pow(cos(theta),2) - .5); + aac = (A4E/A0E)*(35./8. * pow(cos(theta),4) - 30./8. * pow(cos(theta),2) + 3./8. ); + + Iad = aaa + aab + aac; + + AD_I.push_back(Iad); + Theta.push_back(theta); + //these arrays arent used, but they are the function (if one just wanted to the legendre fit function that goes over the data. + } + + + vector dydatas; + + for(int i = 0; i < dydata.size(); i++){ + double bbb = dydata[i]; + + dydatas.push_back(bbb/A0E); + 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 + dydatas.push_back(A0E); + dangler.push_back(3.1415/2.); + + int param; + + param = param_run(det_param_token, gamma_energy_token, ang_file_token, sigma_token, j1j2token); + /* + cout << det_param_token << "\n"; + cout << gamma_energy_token << "\n"; + cout << ang_file_token << "\n"; + cout << sigma_token << "\n"; + cout << j1j2token << "\n"; + + cout << param << "\n"; + */ + + //Initialize Graphics Here. + + + HistoGUI gui; + + HistoGUIad gui_ad; + + + int optnum = -1; + menu(); printf("Input option : "); scanf("%d", &optnum); - } + if(optnum < 0){ + optnum = -1; + printf("Negative numbers are not allowed!\nInput option : "); + scanf("%d", &optnum); + } - //if all inputs are valid, plot distribution. - if(param == 1 && optnum == 1){ - optnum = -1; + 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; - //plotting values here - //x will be arctan of mixing ratio. - //y will be log(X^2); - - gui.SetData(tdelta,chisqr); - gui.Init(); - gui.Loop(); - gui.Close(); + //plotting values here + //x will be arctan of mixing ratio. + //y will be log(X^2); -// 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(tdelta,chisqr); + //gui.SetData(Theta,AD_I); + gui.Init(); + gui.Loop(); + gui.Close(); + + // 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); + gui.Init(); + gui.Loop(); + gui.Close(); + */ + + 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.SetData(dangler,dydata); - gui.SetErrors(deydata); - gui.SetFit(residual[0],residual[1],residual[2]); - - gui.Init(); - gui.Loop(); - gui.Close(); -*/ + + } + //exit + if(param == 1 and optnum == 3){ + + return 0; + } - } -//exit - if(param == 1 and optnum == 3){ - - return 0; - } - - - return 1; + return 1; } diff --git a/GUI_AD.h b/GUI_AD.h index 8631aa6..2d03445 100644 --- a/GUI_AD.h +++ b/GUI_AD.h @@ -54,7 +54,8 @@ class HistoGUIad{ double A0; double A2E; - double A4E; + double A4E; + double Ierr; int SetFit(double a, double b, double c){ A0 = a; A2E = b; @@ -385,7 +386,7 @@ 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] / height_scale); + y_errors_wid = ((y_errors[i]/A0) / height_scale); // x_wid2 = (x[i + 1] + x_offset) / width_scale; // y_wid2 = (y[i + 1] + y_offset) / height_scale; //printf("(%f, %f), (%f,%f)\n", x_wid,y_wid,x_wid2,y_wid2); @@ -424,14 +425,14 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do for(int i=0; i < (int) width; i += w_step){ double x_val = i * width_scale - x_offset; sprintf(axis_val, "%.1f", x_val); - XDrawString(disp, wind, DefaultGC(disp, screen), i, axis_y + 10, axis_val, strlen(axis_val)); + XDrawString(disp, wind, DefaultGC(disp, screen), i, axis_y + 10/A0, axis_val, strlen(axis_val)); } int h_step = height / 10; for(int i=0; i < (int) height; i += h_step){ double y_val = i * height_scale - y_offset; sprintf(axis_val, "%.1f", y_val); - XDrawString(disp, wind, DefaultGC(disp, screen), axis_x + 10, i, axis_val, strlen(axis_val)); + XDrawString(disp, wind, DefaultGC(disp, screen), axis_x + 10/A0, i, axis_val, strlen(axis_val)); } @@ -447,7 +448,7 @@ 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] / height_scale); + y_errors_wid = ((y_errors[i]/A0) / height_scale); // x_wid2 = (x[i + 1] + x_offset) / width_scale; // y_wid2 = (y[i + 1] + y_offset) / height_scale; diff --git a/Graphical_Testing/.GUI.h.swp b/Graphical_Testing/.GUI.h.swp deleted file mode 100644 index 276bbf8..0000000 Binary files a/Graphical_Testing/.GUI.h.swp and /dev/null differ diff --git a/Graphical_Testing/ADT b/Graphical_Testing/ADT index 5e039f6..c1d62e1 100755 Binary files a/Graphical_Testing/ADT and b/Graphical_Testing/ADT differ diff --git a/Graphical_Testing/GUI.h b/Graphical_Testing/GUI.h index 7388347..f5e805f 100644 --- a/Graphical_Testing/GUI.h +++ b/Graphical_Testing/GUI.h @@ -210,9 +210,9 @@ double HistoGUI::legval(double theta){ double lg; - double aaa = A0; - double aab = A2E*(1.5 * pow(theta,2) - .5); - double aac = A4E*(35./8. * pow(theta,4) - 30./8. * pow(theta,2) + 3./8. ); + double aaa = A0/A0; + double aab = (A2E/A0)*(1.5 * pow(cos(theta),2) - .5); + double aac = (A4E/A0)*(35./8. * pow(cos(theta),4) - 30./8. * pow(cos(theta),2) + 3./8. ); lg = aaa + aab + aac; diff --git a/Graphical_Testing/ad_test.cxx b/Graphical_Testing/ad_test.cxx index 1cc15f8..ad443e3 100644 --- a/Graphical_Testing/ad_test.cxx +++ b/Graphical_Testing/ad_test.cxx @@ -5,10 +5,13 @@ int main ( int argc, char** argv){ HistoGUI gui; - double A0E = 134.327; - double A2E = -11.7874; - double A4E = 0.760896; + double A0E = 111.514; + double A2E = -54.3312; + double A4E = -76.4778; + //double A0E = 134.327; + //double A2E = -11.7874; + //double A4E = 0.760896; double step = 0.0001; vector Theta; @@ -27,7 +30,7 @@ int main ( int argc, char** argv){ Iad = aaa + aab + aac; - cout << theta << "\n"; + // cout << theta << "\n"; AD_I.push_back(Iad); Theta.push_back(theta); @@ -43,16 +46,16 @@ int main ( int argc, char** argv){ dangle.push_back(1.5619); vector dydata; //= {129.,110.,129.}; - dydata.push_back(129.); - dydata.push_back(110.); - dydata.push_back(129.); - dydata.push_back(115.); + dydata.push_back(129./A0E); + dydata.push_back(110./A0E); + dydata.push_back(129./A0E); + dydata.push_back(115./A0E); vector deydata; //= {10.,10.,10.}; - deydata.push_back(10.); - deydata.push_back(10.); - deydata.push_back(10.); - deydata.push_back(10.); + deydata.push_back(5.); + deydata.push_back(5.); + deydata.push_back(5.); + deydata.push_back(5.); diff --git a/Racah.h b/Racah.h index 11403e8..d354b95 100644 --- a/Racah.h +++ b/Racah.h @@ -5,7 +5,7 @@ using namespace std; - +//this program doesnt work and it is not used, I spent so much time transcribing it that maybe, one day some poor soul will come and try to fix it to match rose and brinks. double RACAH(double A[6]){