diff --git a/.Chi_error.h.swp b/.Chi_error.h.swp new file mode 100644 index 0000000..c1253ac Binary files /dev/null and b/.Chi_error.h.swp differ diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6f13b0f --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +ADEX diff --git a/1062b.csv b/1062b.csv new file mode 100644 index 0000000..7056ee1 --- /dev/null +++ b/1062b.csv @@ -0,0 +1,4 @@ +angle,Y,Yerr +45,1.29,.10 +90,1.10,.10 +135,1.29,.10 diff --git a/AD.cxx b/AD.cxx index 368144d..93f726e 100644 --- a/AD.cxx +++ b/AD.cxx @@ -1,18 +1,9 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "global.h" -#include "GUI.h" -#include "Coeff.h" +//#include "GUI.h" +#include "GUI_Base.h" +#include "Functlib.h" #include "QDK.h" -#include "Racah.h" #include "Cleb.h" @@ -20,40 +11,12 @@ //For more information about angular distributions, read Rose and Brinks, and Frank Moore (1988). using namespace std; -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 \n"; - std::cout<< "3 - Clear J1 J2 Memory \n"; - std::cout<< "4 - Exit \n"; - std::cout<< "5 - Legfit to Data \n"; - std::cout<< "6 - Generate Data \n"; - std::cout<<"=========================================== \n";} -void Readme(){ - std::cout<< " The program calculates Chi-Squared values \n"; - std::cout<< " from experimental angular distributions as \n"; - std::cout<< " a function of multipole ratios using the \n"; - std::cout<< " theoretical angular distribution formulae \n"; - std::cout<< " in rose and brink \n"; - 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 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"; - std::cout<< " ad.txt is generated with geometric stats.\n"; - std::cout<< " Then at the end it generates CG, W3J & W6J\n"; - -} #ifndef __CINT__ int main(int argc,char **argv){ + HistoGUI gui; - // menu(); + double j1, j2; @@ -61,8 +24,7 @@ int main(int argc,char **argv){ double detradius, targetdistance, detthickness = -1; double Energy; //input .csv file - double Sigma; - double Feeding; + double Sigma; //-------------------------------- @@ -71,8 +33,8 @@ int main(int argc,char **argv){ detthickness = 5.; Energy = 1062.; - Sigma = .1; - Feeding = .1; + Sigma = 0.; + j1 = 2.; j2 = 1.; @@ -86,7 +48,7 @@ int main(int argc,char **argv){ int gamma_energy_token = -1; int ang_file_token = -1; int sigma_token = -1; - int feeding_token = -1; + int j1j2token = -1; if(test == 1){ @@ -94,11 +56,11 @@ int main(int argc,char **argv){ gamma_energy_token = 1; ang_file_token = -1; sigma_token = 1; - feeding_token = 1; + j1j2token = 1; } -/* +if( test == 0 ) { printf("Input detector radius, target distance, and detector thickness : "); scanf("%lf,%lf,%lf",&detradius,&targetdistance,&detthickness); @@ -132,17 +94,18 @@ int main(int argc,char **argv){ if(Energy > 0){ gamma_energy_token = 1; printf(" --- Gamma Energy loaded --- \n"); - }*/ + } - QKN(Energy, detradius, targetdistance, detthickness); //then calucate QD2 and QD4, and replace the 0 with them. - double QD2 = .7; - double QD4 = .3; +} - //input data file of Angular Data (theta, Yexp, Yerr) - string fname; - cout<<"Enter the file name : "; - cin>>fname; + 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; @@ -152,24 +115,66 @@ int main(int argc,char **argv){ 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()){ - 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]); - - } - cout<< " --- Angular Data File Loaded --- \n"; - }else cout<< "Could not open the file\n"; + 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"; @@ -177,50 +182,25 @@ int main(int argc,char **argv){ // cout << eydata[i] << "\n"; } - // Legendre Polinomial fit using adata, ydata, eydata; +// Legendre Polinomial fit using adata, ydata, eydata; - vector int_adata(adata.size()); -std::transform(adata.begin()+1, adata.end(), std::back_inserter(int_adata), StrToInt); + adata.erase(adata.begin()); + ydata.erase(ydata.begin()); + eydata.erase(eydata.begin()); - vector double_angle(int_adata.size() -4); - for(int i = 4; i< int_adata.size(); i++){ - // printf("angle = %d\n",int_adata[i]); - double_angle.push_back((double)int_adata[i]*3.14159/180); - } - - vector int_ydata(ydata.size()); -std::transform(ydata.begin()+1, ydata.end(), std::back_inserter(int_ydata), StrToInt); - - for(int i = 4; i< int_ydata.size(); i++){ - // printf("Ydata = %d\n",int_ydata[i]); - } - - vector double_ydata(int_ydata.size() -4); - for(int i = 4; i< int_ydata.size(); i++){ - // printf("angle = %d\n",int_adata[i]); - double_ydata.push_back((double)int_ydata[i]); - - } + vector dangle = string_to_double_vector(adata); + vector dydata = string_to_double_vector(ydata); + vector deydata = string_to_double_vector(eydata); - vector int_eydata(eydata.size()); -std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), StrToInt); + vector dangler; - for(int i = 4; i< int_eydata.size(); i++){ - // printf("Error data = %d\n",int_eydata[i]); + 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 = 3; i< double_angle.size();i++){ - // printf("rad angle = %lf\n",double_angle[i]); - - //displays the angles inputed as doubles in radians. - } - - for(int i = 3; i< double_ydata.size();i++){ - // printf("y_data = %lf\n",double_ydata[i]); - - //displays the angles inputed as doubles in radians. - } /* for(int i=0; i= 0 ){ sigma_token = 1; printf(" --- Sigma Loaded --- \n"); } - - printf("Please enter feeding parameter : "); - scanf("%lf", &Feeding); - if(Feeding >= 0 and Feeding < 1){ - feeding_token = 1; - }else{ - do{ - printf("Negative numbers are not allowed!\nRe-enter Feeding : "); - scanf("%lf", &Feeding); - }while(Feeding < 0);} - - if(Feeding >= 0 ){ - feeding_token = 1; printf(" --- Feeding Loaded --- \n"); - } - + printf("Please enter J1,J2 : "); scanf("%lf,%lf",&j1,&j2); if(j1 >= 0 && j2 >= 0 ){ @@ -380,32 +346,12 @@ std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), S if(j1 >= 0 && j2 >= 0 ){ j1j2token = 1; printf(" --- J1 & J2 Loaded --- \n"); } -*/ - - /* - printf("detector radius = %lf\n",detradius); - printf("detector distance = %lf\n",targetdistance); - printf("detector width = %lf\n",detthickness); +} - printf("Gamma-Ray energy = %lf\n", Energy); - printf("Sigma Distribution = %lf\n", Sigma); - printf("Feeding param = %lf\n", Feeding); - printf("J1,J2 = %lf,%lf\n",j1,j2); - */ -/* - double J = j1 + j2; - double m1 = j1; - double m2 = j2; - double M = m1 + m2; - double j3 = J; - double m3 = j3; - double j4 = j3; - double j5 = j1 + j2 + j3; - double j6 = j2 + j3; -*/ // Calculate BK(j1) for perfect allignment // + double A[6]; double js = sqrt(2*j1 +1); double j12 = 2*j1; @@ -413,10 +359,11 @@ std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), S double cg12,cg22,cg24,cg14 = 0.; double I1, I2 = 0.; - double Bk11, Bk12 = 0.; double cg1, cg2 = 0.; + double Bk11, Bk12 = 0.; + if(Sigma == 0){ if(IS == 1){ A[0] = j1; A[1] = j1; @@ -441,6 +388,15 @@ std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), S 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){ @@ -451,8 +407,8 @@ std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), S A[4] = 0.; A[5] = 0.; - cg1 = CG2(A); - + cg1 = CG2(A); + Bk11 = pow(-1,j1) * js * cg1; A[2] = 4.; @@ -460,19 +416,23 @@ std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), S 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; + 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); @@ -483,6 +443,8 @@ std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), S } double cn1 = 1./sum1; +// cout << "cn1" << " " << cn1 << "\n"; + double AL0 = 0.; double A0 = j1 - j2; double A0b = abs(A0); @@ -495,102 +457,130 @@ std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), S double am11,amsq1,x1,ex1 = 0.; //calculate Bk(j1) for gaussian W(m1) or non zero Sigma - A[0] = j1; - A[1] = j1; +//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"; - double sfact = sqrt(2*j1 +1); - double cgg = 0.; - double tTerm = 0.; - double II = 0.; - for(int i = 2; i<= 4; i = i + 2){ + + } + + 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; - 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; - cgg = CG2(A); - ex1 = exp(x1); - tTerm = cn1*ex1*pow(-1,II) * sfact*cgg; - if(i == 2) Bk11 = Bk11 * tTerm; - if(i == 4) Bk12 = Bk12 * tTerm; + + 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 = " < chisqr; vector tdelta; for(int i = 0; i< points; i++){ - tan_delta = i + step + delta_min; + 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; + X22 = pow(abs(a2E -a2T),2)/(abs(a2T)); + X24 = pow(abs(a4E -a4T),2)/(abs(a4T)); + X2_total = (X22 + X24)/2; - chisqr.push_back(X2_total); + 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); } - - - int param; - param = param_run(det_param_token, gamma_energy_token, ang_file_token, sigma_token, feeding_token, j1j2token); - - + 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(); @@ -701,32 +698,44 @@ std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), S if(param == 1 && optnum == 1){ optnum = -1; - //plotting values here - std::vector x; - std::vector y; - - //x will be arctan of mixing ratio. - //y will be log(X^2); - - - for(int i = 0; i < 50; i++){ - x.push_back( i ); - y.push_back( pow(i, 0.5) ); - } - gui.SetData(x,y); + //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(); -/* - printf("Input option : "); - scanf("%d", &optnum); - if(optnum < 0){ - printf("Input option : "); - scanf("%d", &optnum); - }*/ - } +// 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(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; + } return 1; diff --git a/Cleb.h b/Cleb.h index 4d4be1e..63f4983 100644 --- a/Cleb.h +++ b/Cleb.h @@ -1,13 +1,5 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include + +#include "global.h" //To compile : g++ AD.cxx -o {Input Executable Name} -lX11 //#include "Coeff.h" @@ -16,500 +8,8 @@ using namespace std; -double CG2(double A[6]){ - - double ii[6]; - double im[6]; - double ix[9]; - double nx[5]; - - for(int i = 0;i<6;i++){ - ii[i] = 2.0*A[i]; - } - double Cleb = 0.; - double IA = ii[0]; - double IB = ii[1]; - double IC = ii[2]; - double ID = ii[3]; - double IE = ii[4]; - double IF = ii[5]; - - if(ID + IE - IF != 0){return 0;} - double K0 = IA + IB + IC; - if((int)K0 % 2 !=0){return 0;} - - double K1 = IA + IB - IC; - double K2 = IC + IA - IB; - double K3 = IB + IC - IA; - double K4 = IA - abs(IB-IC); - double K5 = IB - abs(IC-IA); - double K6 = IC - abs(IA-IB); - - double karr[6] = {K1,K2,K3,K4,K5,K6}; - double kmn=karr[0]; - double kmx=karr[0]; - for(int i=0;i<6;i++) - { - if(kmn>karr[i]) - { - kmn=karr[i]; - } - else if(kmx= IB){ - if(IC >= IB){ - if(IB != 0){ - if(IE < 0){ - //assuming the logic passes; - FC2 = IC +1; - - - IT = K0/2 + 1; - ix[0] = IT-IC; - ix[1] = IT-IB; - ix[2] = IT-IA; - ix[3] = (IA+ID)/2 +1; - ix[4] = ix[3] - ID; - ix[5] = (IB+IE)/2 + 1; - ix[6] = ix[5] - IE; - ix[7] = (IC+IF)/2 +1; - ix[8] = ix[7] - IF; - //lgamma(arg) returns the log of the factorial of arg if arg is a natural number. - double sqfclg = log(FC2) - lgamma(IT+1); - double IXI; - for(int i = 0; i<9;i++){ - IXI = ix[i]; - sqfclg = sqfclg + lgamma(IXI); - } - sqfclg = 0.5 * sqfclg; - - double xarr[3] = {ix[0],ix[4],ix[5]}; - - double xmn=xarr[0]; - double xmx=xarr[0]; - for(int i=0;i<3;i++) - { - if(xmn>xarr[i]) - { - xmn=xarr[i]; - } - else if(xmxcarr[i]) - { - cmn=carr[i]; - } - else if(cmx NZMIN){return 0;} - double SS=0.; - double S1= pow((-1.),(NZMIN-1)); - for(int NZ = NZMIN;NZ<=NZMAX;NZ++){ - - double NZM1 = NZ -1; - nx[0] = ix[0]-NZM1; - nx[1] = ix[4]-NZM1; - nx[2] = ix[5]-NZM1; - nx[3] = NZ-NZ2; - nx[4] = NZ-NZ3; - double termlg = sqfclg - lgamma(NZ); - for(int i = 0; i<5; i++){ - IXI = nx[i]; - termlg= termlg -lgamma(IXI); - SS = SS + S1*exp(termlg); - S1 =-S1; - Cleb=sgnfac*SS; - } - } - }else{ - // if (IE >=0) _>>> - ID = -ID; - IE = -IE; - IF = -IF; - if((int)(IA+IB-IC)/2 %2 != 0) sgnfac = (-1)*sgnfac; - - IT = K0/2 + 1; - ix[0] = IT-IC; - ix[1] = IT-IB; - ix[2] = IT-IA; - ix[3] = (IA+ID)/2 +1; - ix[4] = ix[3] - ID; - ix[5] = (IB+IE)/2 + 1; - ix[6] = ix[5] - IE; - ix[7] = (IC+IF)/2 +1; - ix[8] = ix[7] - IF; - //lgamma(arg) returns the log of the factorial of arg if arg is a natural number. - double sqfclg = log(FC2) - lgamma(IT+1); - double IXI; - for(int i = 0; i<9;i++){ - IXI = ix[i]; - sqfclg = sqfclg + lgamma(IXI); - } - sqfclg = 0.5 * sqfclg; - - double xarr[3] = {ix[0],ix[4],ix[5]}; - - double xmn=xarr[0]; - double xmx=xarr[0]; - for(int i=0;i<3;i++) - { - if(xmn>xarr[i]) - { - xmn=xarr[i]; - } - else if(xmxcarr[i]) - { - cmn=carr[i]; - } - else if(cmx NZMIN){return 0;} - double SS=0.; - double S1= pow((-1.),(NZMIN-1)); - for(int NZ = NZMIN;NZ<=NZMAX;NZ++){ - - double NZM1 = NZ -1; - nx[0] = ix[0]-NZM1; - nx[1] = ix[4]-NZM1; - nx[2] = ix[5]-NZM1; - nx[3] = NZ-NZ2; - nx[4] = NZ-NZ3; - double termlg = sqfclg - lgamma(NZ); - for(int i = 0; i<5; i++){ - IXI = nx[i]; - termlg= termlg -lgamma(IXI); - SS = SS + S1*exp(termlg); - S1 =-S1; - Cleb=sgnfac*SS; - } - } - - - - } -// if (IB != 0) _>>>> - }else{ - Cleb = sgnfac; - return Cleb; - - } -// if (IC < IB) _>>> - - }else{ - - IT = IC; - IC = IB; - IB = IT; - IT = IF; - IF = -IE; - IE = -IT; - sgnfac = sqrt((2.*A[2]+1.0)/(2.*A[1]+1.0)); - if((int)(im[0]-im[3])/2 %2 != 0){sgnfac = (-1.)*sgnfac;} - - Cleb = sgnfac; - return Cleb; - } - -// if(IA < IB) _>> - }else{ - - if(IA >= IC){ - IT = IC; - IC = IB; - IB = IT; - IT = IF; - IF = -IE; - IE = -IT; - - Cleb = sgnfac; - return Cleb; -// if (IA < IC) as well as (IA < IB) - }else{ - - IT =IA; - IA =IB; - IB =IT; - IT =ID; - ID =IE; - IE =IT; - - if((int)(K1/2) %2 != 0) {sgnfac = -1.;} -//call to 135 again. - if(IB != 0){ - if(IE < 0){ - //assuming the logic passes; - FC2 = IC +1; - - - IT = K0/2 + 1; - ix[0] = IT-IC; - ix[1] = IT-IB; - ix[2] = IT-IA; - ix[3] = (IA+ID)/2 +1; - ix[4] = ix[3] - ID; - ix[5] = (IB+IE)/2 + 1; - ix[6] = ix[5] - IE; - ix[7] = (IC+IF)/2 +1; - ix[8] = ix[7] - IF; - //lgamma(arg) returns the log of the factorial of arg if arg is a natural number. - double sqfclg = log(FC2) - lgamma(IT+1); - double IXI; - for(int i = 0; i<9;i++){ - IXI = ix[i]; - sqfclg = sqfclg + lgamma(IXI); - } - sqfclg = 0.5 * sqfclg; - - double xarr[3] = {ix[0],ix[4],ix[5]}; - - double xmn=xarr[0]; - double xmx=xarr[0]; - for(int i=0;i<3;i++) - { - if(xmn>xarr[i]) - { - xmn=xarr[i]; - } - else if(xmxcarr[i]) - { - cmn=carr[i]; - } - else if(cmx NZMIN){return 0;} - double SS=0.; - double S1= pow((-1.),(NZMIN-1)); - for(int NZ = NZMIN;NZ<=NZMAX;NZ++){ - - double NZM1 = NZ -1; - nx[0] = ix[0]-NZM1; - nx[1] = ix[4]-NZM1; - nx[2] = ix[5]-NZM1; - nx[3] = NZ-NZ2; - nx[4] = NZ-NZ3; - double termlg = sqfclg - lgamma(NZ); - for(int i = 0; i<5; i++){ - IXI = nx[i]; - termlg= termlg -lgamma(IXI); - SS = SS + S1*exp(termlg); - S1 =-S1; - Cleb=sgnfac*SS; - } - } - }else{ - // if (IE >=0) _>>> - ID = -ID; - IE = -IE; - IF = -IF; - if((int)(IA+IB-IC)/2 %2 != 0) sgnfac = (-1.)*sgnfac; - - IT = K0/2 + 1; - ix[0] = IT-IC; - ix[1] = IT-IB; - ix[2] = IT-IA; - ix[3] = (IA+ID)/2 +1; - ix[4] = ix[3] - ID; - ix[5] = (IB+IE)/2 + 1; - ix[6] = ix[5] - IE; - ix[7] = (IC+IF)/2 +1; - ix[8] = ix[7] - IF; - //lgamma(arg) returns the log of the factorial of arg if arg is a natural number. - double sqfclg = log(FC2) - lgamma(IT+1); - double IXI; - for(int i = 0; i<9;i++){ - IXI = ix[i]; - sqfclg = sqfclg + lgamma(IXI); - } - sqfclg = 0.5 * sqfclg; - - double xarr[3] = {ix[0],ix[4],ix[5]}; - - double xmn=xarr[0]; - double xmx=xarr[0]; - for(int i=0;i<3;i++) - { - if(xmn>xarr[i]) - { - xmn=xarr[i]; - } - else if(xmxcarr[i]) - { - cmn=carr[i]; - } - else if(cmx NZMIN){return 0;} - double SS=0.; - double S1= pow((-1.),(NZMIN-1)); - for(int NZ = NZMIN;NZ<=NZMAX;NZ++){ - - double NZM1 = NZ -1; - nx[0] = ix[0]-NZM1; - nx[1] = ix[4]-NZM1; - nx[2] = ix[5]-NZM1; - nx[3] = NZ-NZ2; - nx[4] = NZ-NZ3; - double termlg = sqfclg - lgamma(NZ); - for(int i = 0; i<5; i++){ - IXI = nx[i]; - termlg= termlg -lgamma(IXI); - SS = SS + S1*exp(termlg); - S1 =-S1; - Cleb=sgnfac*SS; - } - } - - - - } -// if (IB != 0) _>>>> - }else{ - Cleb = sgnfac; - return Cleb; - - } - - - - - - - - - - - } - - - - - - } - - - - return Cleb; -} - - -double CG(double J, double M, double j1, double m1, double j2, double m2){ +//double CG(double J, double M, double j1, double m1, double j2, double m2){ +double CG(double j1, double j2, double J, double m1, double m2, double M){ //recall that j1,m1 + j2,m2 = J,M if(M != m1 + m2) return 0; @@ -517,137 +17,119 @@ double CG(double J, double M, double j1, double m1, double j2, double m2){ double Jmin = abs(j1 - j2); double Jmax = j1+j2; + if(J < Jmin || Jmax < J) return 0; - double a0 = (2*J+1.0)*factorial(J+j1-j2) * factorial(J-j1+j2) * factorial(j1+j2-J)/factorial(J+j1+j2+1.0); + double a0 = (2*J+1.0)*tgamma(J+j1-j2+1) * tgamma(J-j1+j2+1) * tgamma(j1+j2-J+1)/tgamma(J+j1+j2+1.0 +1); double A0 = sqrt(a0); - double a = factorial(J+M) *factorial(J-M); - double a1= factorial(j1+m1) *factorial(j1-m1); - double a2= factorial(j2+m2) *factorial(j2-m2); - double A = sqrt( a * a1 * a2); - int pmax = min( min(j1+j2-J,j1-m1),j2 + m2); - - double cg = 0.; + double a = tgamma(J+M+1) *tgamma(J-M+1); + double a1= tgamma(j1+m1+1) *tgamma(j1-m1+1); + double a2= tgamma(j2+m2+1) *tgamma(j2-m2+1); - for( int p =0; p<=pmax;p++){ - - double p1 = factorial(j1+j2-J-p); - double p2 = factorial(j1-m1-p); - double p3 = factorial(j2+m2-p); - double p4 = factorial(J -j2 +m1 +p); - double p5 = factorial(J -j1 -m2 +p); - double t = pow(-1,p)/(factorial(p) * p1 * p2 * p3 * p4 * p5); - - cg += t; - } - return A0*A*cg; -} - -double CG_a(double A[6]){ - //recall that j1,m1 + j2,m2 = J,M - //testing assingments until the output is the same. - double j1 = A[0]; - double j2 = A[1]; - double J = A[2]; - double m1 = A[3]; - double m2 = A[4]; - double M = A[5]; - - if(M != m1 + m2) return 0; - - double Jmin = abs(j1 - j2); - double Jmax = j1+j2; - - if(J < Jmin || Jmax < J) return 0; - - double a0 = (2*J+1.0)*factorial(J+j1-j2) * factorial(J-j1+j2) * factorial(j1+j2-J)/factorial(J+j1+j2+1.0); - double A0 = sqrt(a0); - - double a = factorial(J+M) *factorial(J-M); - double a1= factorial(j1+m1) *factorial(j1-m1); - double a2= factorial(j2+m2) *factorial(j2-m2); double A1 = sqrt( a * a1 * a2); - int pmax = min( min(j1+j2-J,j1-m1),j2 + m2); + double cg = 0.; for( int p =0; p<=pmax;p++){ - double p1 = factorial(j1+j2-J-p); - double p2 = factorial(j1-m1-p); - double p3 = factorial(j2+m2-p); - double p4 = factorial(J -j2 +m1 +p); - double p5 = factorial(J -j1 -m2 +p); - double t = pow(-1,p)/(factorial(p) * p1 * p2 * p3 * p4 * p5); - + + double p1 = tgamma(j1+j2-J-p+1); + double p2 = tgamma(j1-m1-p+1); + double p3 = tgamma(j2+m2-p+1); + double p4 = tgamma(J -j2 +m1 +p+1); + double p5 = tgamma(J -j1 -m2 +p+1); + double t = pow(-1,p)/(tgamma(p+1) * p1 * p2 * p3 * p4 * p5); + cg += t; } return A0*A1*cg; } -double ThreeJsym(double j1, double m1, double j2, double m2, double j3, double m3){ - //[j1 j2 j3] = (-1)^(j1-j2-m3)/ sqrt(2*j3+1) * CG(j3, -m3, j1, m1, j2, m2) - //[m1 m2 m3] + + +double CG2(double A[6]){ + //recall that j1,m1 + j2,m2 = J,M + double j1 = A[0]; + double j2 = A[1]; + double J = A[2]; + double m1 = A[3]; + double m2 = A[4]; + double M = A[5]; + if(M != m1 + m2) return 0; + + double Jmin = abs(j1 - j2); + double Jmax = j1+j2; + + if(J < Jmin || Jmax < J) return 0; + + + double a0 = (2*J+1.0)*tgamma(J+j1-j2+1) * tgamma(J-j1+j2+1) * tgamma(j1+j2-J+1)/tgamma(J+j1+j2+1.0 +1); + double A0 = sqrt(a0); + + + double a = tgamma(J+M+1) *tgamma(J-M+1); + double a1= tgamma(j1+m1+1) *tgamma(j1-m1+1); + double a2= tgamma(j2+m2+1) *tgamma(j2-m2+1); + + double A1 = sqrt( a * a1 * a2); + int pmax = min( min(j1+j2-J,j1-m1),j2 + m2); - return pow(-1,j1 -j2 -m3)/sqrt(2*j3+1)*CG(j3,-m3,j1,m1,j2,m2); -} -double SixJsym(double j1, double j2, double j3, double j4, double j5, double j6){ + double cg = 0.; -//--------------------------------------------------------------------------------// - // The six j symbol describes the coupling between j1 j2 and j3 to J - j1. - // essentially a triangle of angular momentum rules between these. - // j1 = j1 - // j2 = j2 - // j3 = j1 + j2 - // j4 = j3 - // j5 = J = j1 + j2 + j3 - // j6 = j2 + j3 -// ----------------------------------------------------------------------------- // -// the following conditions check the triangle selection rules + for( int p =0; p<=pmax;p++){ - 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; -// now that they have been checked, we can go ahead and calculate sixJ. + double p1 = tgamma(j1+j2-J-p+1); + double p2 = tgamma(j1-m1-p+1); + double p3 = tgamma(j2+m2-p+1); + double p4 = tgamma(J -j2 +m1 +p+1); + double p5 = tgamma(J -j1 -m2 +p+1); + double t = pow(-1,p)/(tgamma(p+1) * p1 * p2 * p3 * p4 * p5); - double sixj = 0.0; - float m1 = -j1; - float m2 = -j2; - float m3 = -j3; - float m4 = -j4; - float m5 = -j5; - float m6 = -j6; - - for(; m1 <= j1; m1 = m1 +1){ - for(; m2 <= j2; m2 = m2 +1){ - for(; m3 <= j3; m3 = m3 + 1){ - for(; m4 <= j4; m4 = m4 +1){ - for(; m5 <= j5; m5 = m5 + 1){ - for(; m6 <= j6; m6 = m6 +1){ - - double h = (j1 - m1) + (j2 - m2) + (j3 -m3) + (j4 - m4) + (j5 - m5) + (j6 - m6); - double b1 = ThreeJsym(j1, -m1, j2, -m2, j3, -m3); - double b2 = ThreeJsym(j1, m1, j5, -m5, j6, m6); - double b3 = ThreeJsym(j4, m4, j2, m2, j6, -m6); - double b4 = ThreeJsym(j4, -m4, j5, m5, j3, m3); - double b = b1 * b2 * b3 * b4; - - sixj += pow(-1,h)*b; - - } - } - } - } - } + cg += t; } + return A0*A1*cg; +} +/* - return sixj; +int main(int argc, char ** argv){ + + //if mod (2*J1,2) = 1 do this +// double A0[6] = {1,1,2,0,0,0}; +// double A1[6] = {1,1,4,0,0,0}; + + //if mod(2*j1,20 = 0 do this + + double j1 = 2.; + double j2 = 2.; + + + double A0[6] = {j1,j1,2,.5,-.5,0}; + double A1[6] = {j1,j1,2,-.5,.5,0}; + double A2[6] = {j1,j1,4,-.5,.5,0}; + double A3[6] = {j1,j1,4,.5,-.5,0}; + + double cgr0= CG(j1,j1,2,.5,-.5,0); + double cgr1= CG(j1,j1,2,-.5,.5,0); + double cgr2= CG(j1,j1,4,.5,-.5,0); + double cgr3= CG(j1,j1,4,-.5,.5,0); + + + printf("CGR0 = %lf\n",cgr0); + printf("CGR1 = %lf\n",cgr1); + printf("CGR2 = %lf\n",cgr2); + printf("CGR3 = %lf\n",cgr3); + + + + + return 0; } +*/ diff --git a/Coeff.h b/Coeff.h deleted file mode 100644 index de6ca97..0000000 --- a/Coeff.h +++ /dev/null @@ -1,71 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -//To compile : g++ AD.cxx -o {Input Executable Name} -lX11 - -using namespace std; - - -int param_run(int dt, int gt, int at, int st, int ft, int jt){ - int param = 0; - if( dt == 1 && gt == 1 && at == 1 && st == 1 && ft == 1 && jt == 1){ - param =0; - }else param = 1; - - return param; -} - - - -int StrToInt(std::string const& s) -{ - std::istringstream iss(s); - int value; - - if (!(iss >> value)) throw std::runtime_error("invalid int"); - - return value; -} - - -int factorial(int fact){ - for(int i=1;i<=fact;i++){ - fact=fact*i; - } - return fact; -} - - -double FACTLOG(int num){ - - double faclog[170]; - int RI = 0; - - faclog[0] = 0.; - faclog[1] = 0.; - - for( int i = 3; i < 170; i++){ - RI = i - 1; - faclog[i] = log(RI) + faclog[i-1]; - - } - - double flog = faclog[num]; - - - return flog; -} - - - - - - - diff --git a/Functlib.h b/Functlib.h new file mode 100644 index 0000000..6a02ab6 --- /dev/null +++ b/Functlib.h @@ -0,0 +1,137 @@ + +#include "global.h" +//To compile : g++ AD.cxx -o {Input Executable Name} -lX11 + +using namespace std; + + +void menu(){ + + std::cout<< " \n"; + 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<< "3 - Exit \n"; + std::cout <<"==========================================\n"; + std::cout<< " \n"; +} + +void Readme(){ + + std::cout<< " \n"; + std::cout<< " The program calculates Chi-Squared values \n"; + std::cout<< " from experimental angular distributions as \n"; + std::cout<< " a function of multipole ratios using the \n"; + std::cout<< " theoretical angular distribution formulae \n"; + std::cout<< " in rose and brink \n"; + 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 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"; + std::cout<< " ad.txt is generated with geometric stats.\n"; + std::cout<< " Then at the end it generates CG/Racah\n"; + std::cout<< " \n"; +} + +int param_run(int dt, int gt, int at, int st, int jt){ + int param = 0; + if( dt == 1 && gt == 1 && at == 1 && st == 1 && jt == 1){ + param = 1; + }else param = 0; + + return param; +} + + + +int StrToInt(std::string const& s) +{ + std::istringstream iss(s); + int value; + + if (!(iss >> value)) throw std::runtime_error("invalid int"); + + return value; +} + +double StrToDouble(string const& s){ + istringstream iss(s); + double value; + + if(!(iss >> value)) throw runtime_error("invalid double"); + + return value; +} + + +int factorial(int fact){ + for(int i=1;i<=fact;i++){ + fact=fact*i; + } + return fact; +} + +double factorial_d( double fact){ + + fact = tgamma(fact + 1); + + return fact; +} + + +double FACTLOG(int num){ + + double faclog[170]; + int RI = 0; + + faclog[0] = 0.; + faclog[1] = 0.; + + for( int i = 3; i < 170; i++){ + RI = i - 1; + faclog[i] = log(RI) + faclog[i-1]; + + } + + double flog = faclog[num]; + + + return flog; + + +} + + +vector string_to_double_vector( vector string_vec){ + + vector::iterator sit; + + vector dv; + + if(string_vec.size() < 0 ) throw overflow_error("Invalid String Vector Input\n"); + + for(sit = string_vec.begin(); sit < string_vec.end(); sit++){ + + double d; + + d = stod(*sit); + + dv.push_back(d); + + } + + return dv; +} + + + + + + + + diff --git a/GUI.h b/GUI.h index dea6673..a7e5be3 100644 --- a/GUI.h +++ b/GUI.h @@ -1,15 +1,7 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "global.h" + //To compile : g++ AD.cxx -o {Input Executable Name} -lX11 -// VERSION : 1.0 +// VERSION : 1.1 (Overlay for Angular Fit Display) // Created with thanks to Daniel Folds Holt - University of Cambridge Cavendish Lab @@ -38,6 +30,8 @@ class HistoGUI{ std::vector x; std::vector y; + std::vector y_errors; + double max_x; double min_x; double max_y; @@ -57,17 +51,35 @@ class HistoGUI{ double height_scale; double y_offset; + + double A0; + double A2E; + double A4E; + int SetFit(double a, double b, double c){ + A0 = a; + A2E = b; + A4E = c; + fit_exists = true; + return 1; + } + + double old_xl, old_xh, old_yl, old_yh; double old_mouse_x, old_mouse_y; int Init(); int SetData(std::vector a, std::vectorb); + int SetErrors(std::vector a); int Loop(); void Close(){ printf("Closing now!\n"); XCloseDisplay(disp); } int DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win); + int Draw_Fit(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win); int DrawCrosshairs(int mouse_x, int mouse_y); int Zoom(int mouse_x, int mouse_y); + double legval(double theta); + bool fit_exists = false; + // int Legendre_fit(std::vector d1, std::vector d2, double @@ -81,6 +93,12 @@ int HistoGUI::SetData(std::vector a, std::vectorb){ return a.size(); } +int HistoGUI::SetErrors(std::vector a){ + for(int i=0; i < a.size(); i++){ + y_errors.push_back(a[i]); + } + return a.size(); +} int HistoGUI::Init(){ disp = XOpenDisplay(NULL); @@ -89,6 +107,8 @@ int HistoGUI::Init(){ exit(1); } + //fit_exists = false; + screen = DefaultScreen(disp); wind = XCreateSimpleWindow(disp, RootWindow(disp, screen), 10, 10, 600, 400, 1, BlackPixel(disp, screen), WhitePixel(disp, screen)); @@ -185,11 +205,111 @@ int HistoGUI::Zoom(int mouse_x, int mouse_y){ return 1; } +double HistoGUI::legval(double theta){ + + double lg; + + double aaa = A0; + double aab = A2E*(1.5 * pow(theta,2) - 1); + double aac = A4E*(35./8. * pow(theta,4) - 30./8. * pow(theta,2) + 3./8. ); + + lg = aaa + aab + aac; + return lg; +} + +int HistoGUI::Draw_Fit(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win){ + + + int j1, j2; + unsigned int j3, j4; + Window root_return; + + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + + //printf("Width = %u, Height = %u, x = %i, y = %i\n", width, height, j1, j2); + //printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win); + + double x_step;// = (width * 0.8) / x.size(); + double y_step;// = (height * 0.8) / x.size(); + + + if(x_low_win == -1 and y_low_win == -1 and x_hi_win == -1 and y_hi_win == -1){ + // Audomatically decide data postition + max_x = x[0]; + max_y = y[0]; + min_x = x[0]; + min_y = y[0]; + + for(int i=0; i max_x) max_x = x[i]; + if(x[i] < min_x) min_x = x[i]; + if(y[i] > max_y) max_y = y[i]; + if(y[i] < min_y) min_y = y[i]; + } + + width_scale = (max_x - min_x) / (0.8 * width); + x_offset = 0.5 * ((max_x - min_x) / 0.8) - 0.5 * (min_x + max_x); + + height_scale = -1. * (max_y - min_y) / (0.8 * height); + y_offset = -0.5 * ((max_y - min_y) / 0.8) + -0.5 * (min_y + max_y); + + XSetForeground(disp, DefaultGC(disp,screen), xcolour_red.pixel); + + double x_wid ; + double y_wid ; + + for(int i=0; i < (int) width; i += 2){ + double x_val = i * width_scale - x_offset; + double y_val = legval(x_val); + + y_wid = (y_val + y_offset) / height_scale; + XFillRectangle(disp, wind, DefaultGC(disp, screen), i -2, y_wid -2, 4, 4); + } + // XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); + + } else { + //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win + + width_scale = (x_hi_win - x_low_win) / width; + x_offset = -1.0 * x_low_win; + + height_scale = (y_hi_win - y_low_win) / height; + y_offset = -1.0 * y_low_win; + + + XSetForeground(disp, DefaultGC(disp,screen), xcolour_red.pixel); + + double x_wid ; + double y_wid ; + + for(int i=0; i < (int) width; i += 2){ + double x_val = i * width_scale - x_offset; + double y_val = legval(x_val); + + y_wid = (y_val + y_offset) / height_scale; + XFillRectangle(disp, wind, DefaultGC(disp, screen), i -2, y_wid -2, 4, 4); + } + } + + old_xl = x_low_win; + old_yl = y_low_win; + old_xh = x_hi_win; + old_yh = y_hi_win; + + + + return 1; +} + int HistoGUI::DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win){ + + if(fit_exists){ + Draw_Fit(x_low_win,y_low_win,x_hi_win,y_hi_win); - + } + int j1, j2; unsigned int j3, j4; Window root_return; @@ -257,19 +377,26 @@ int HistoGUI::DrawData(double x_low_win, double y_low_win, double x_hi_win, doub XSetForeground(disp, DefaultGC(disp,screen), 0); double x_wid ; double y_wid ; - double x_wid2; - double y_wid2; + // double x_wid2; + // double y_wid2; + double y_errors_wid; - for(int i=0; i < x.size() - 2; i++){ + 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; - x_wid2 = (x[i + 1] + x_offset) / width_scale; - y_wid2 = (y[i + 1] + y_offset) / height_scale; + y_errors_wid = (y_errors[i] / 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); - XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2); - XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4); + // XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -3, y_wid -3, 6, 6); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -1, y_wid + y_errors_wid, 2, y_errors_wid); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -1, y_wid - y_errors_wid, 2, y_errors_wid); +//error bar testing + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -4, y_wid + y_errors_wid, 8, 2); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -4, y_wid - y_errors_wid, 8, 2); } - XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); + // XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); } else { //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win @@ -310,19 +437,29 @@ int HistoGUI::DrawData(double x_low_win, double y_low_win, double x_hi_win, doub XSetForeground(disp, DefaultGC(disp,screen), 0); double x_wid ; double y_wid ; - double x_wid2; - double y_wid2; + // double x_wid2; + // double y_wid2; - for(int i=0; i < x.size() - 2; i++){ + double y_errors_wid; + + 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; - x_wid2 = (x[i + 1] + x_offset) / width_scale; - y_wid2 = (y[i + 1] + y_offset) / height_scale; + + y_errors_wid = (y_errors[i] / 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); - XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2); - XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4); + // XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -3, y_wid -3, 6, 6); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -1, y_wid , 2, y_errors_wid); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -1, y_wid - y_errors_wid, 2, y_errors_wid); +//testig for error bars. + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -4, y_wid + y_errors_wid, 8, 2); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -4, y_wid - y_errors_wid, 8, 2); } - XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); + // XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); } diff --git a/GUI_Base.h b/GUI_Base.h new file mode 100644 index 0000000..255fb5e --- /dev/null +++ b/GUI_Base.h @@ -0,0 +1,418 @@ +#include +#include +#include +#include +#include +#include + + +class HistoGUI{ + + + public: + + HistoGUI(){} + + Display * disp; + Window wind; + XEvent evt; + int screen; + + XColor xcolour; + XColor xcolour_red; + XColor xcolour_veryred; + Colormap cmap; + + + std::vector x; + std::vector y; + double max_x; + double min_x; + double max_y; + double min_y; + + double pos_x; + double pos_y; + double new_pos_x; + double new_pos_y; + + unsigned int width; + unsigned int height; + double window_width; + double window_height; + double width_scale; + double x_offset; + double height_scale; + double y_offset; + + double old_xl, old_xh, old_yl, old_yh; + double old_mouse_x, old_mouse_y; + + int Init(); + int SetData(std::vector a, std::vectorb); + int Loop(); + void Close(){ printf("Closing now!\n"); XCloseDisplay(disp); } + int DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win); + int DrawCrosshairs(int mouse_x, int mouse_y); + int Zoom(int mouse_x, int mouse_y); + + +}; + +int HistoGUI::SetData(std::vector a, std::vectorb){ + for(int i=0; i < a.size(); i++){ + x.push_back(a[i]); + y.push_back(b[i]); + } + return a.size(); +} + +int HistoGUI::Init(){ + + disp = XOpenDisplay(NULL); + if (disp == NULL) { + fprintf(stderr, "Cannot open display\n"); + exit(1); + } + + screen = DefaultScreen(disp); + wind = XCreateSimpleWindow(disp, RootWindow(disp, screen), 10, 10, 600, 400, 1, + BlackPixel(disp, screen), WhitePixel(disp, screen)); + XGrabPointer(disp, wind, False, ButtonPressMask | ButtonReleaseMask | PointerMotionMask, GrabModeAsync, + GrabModeAsync, None, None, CurrentTime); + XSelectInput(disp, wind, ExposureMask | KeyPressMask | ButtonPressMask | ButtonReleaseMask | PointerMotionMask); + XMapWindow(disp, wind); + + // colours + cmap = DefaultColormap(disp, screen); + xcolour.red = 32000; + xcolour.green = 32000; + xcolour.blue = 42000; + xcolour.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour); + + xcolour_red.red = 42000; + xcolour_red.green = 32000; + xcolour_red.blue = 32000; + xcolour_red.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour_red); + + xcolour_veryred.red = 65000; + xcolour_veryred.green = 32000; + xcolour_veryred.blue = 32000; + xcolour_veryred.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour_veryred); + + + + return 1; + +} + +int HistoGUI::DrawCrosshairs(int mouse_x, int mouse_y){ + + int j1, j2; + unsigned int j3, j4; + Window root_return; + + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + XSetForeground(disp, DefaultGC(disp, screen), xcolour_red.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), mouse_x, 0.0, mouse_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, mouse_y, width, mouse_y); + + pos_y = mouse_y * height_scale - y_offset; + pos_x = mouse_x * width_scale - x_offset; + + char coord[32]; + sprintf(coord, "(%.2f, %.2f)", pos_x, pos_y); + XDrawString(disp, wind, DefaultGC(disp, screen), mouse_x + 5, mouse_y + 12, coord, strlen(coord)); + + return 1; +} + +int HistoGUI::Zoom(int mouse_x, int mouse_y){ + + new_pos_x = mouse_x * width_scale - x_offset; + new_pos_y = mouse_y * height_scale - y_offset; + + pos_x = old_mouse_x * width_scale - x_offset; + pos_y = old_mouse_y * height_scale - y_offset; + + + if(std::abs(mouse_x - old_mouse_x) < 20 or std::abs(mouse_y - old_mouse_y) < 20) return 1; + + + //printf("Zooming: %f - %f, %f - %f\n", mouse_x, old_mouse_x,mouse_y,old_mouse_y); + + double x_low_win, y_low_win, x_hi_win, y_hi_win; + if(new_pos_x < pos_x){ + x_low_win = new_pos_x; + x_hi_win = pos_x; + //printf("new low; [%f, %f]\n",x_low_win, x_hi_win); + } else { + x_low_win = pos_x; + x_hi_win = new_pos_x; + //printf("new high; [%f, %f]\n",x_low_win, x_hi_win); + } + if(new_pos_y > pos_y){ + y_low_win = new_pos_y; + y_hi_win = pos_y; + } else { + y_low_win = pos_y; + y_hi_win = new_pos_y; + } + + + //printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win); + XClearWindow(disp, wind); + DrawData(x_low_win, y_low_win, x_hi_win, y_hi_win); + //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win + + return 1; + +} + + +int HistoGUI::DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win){ + + + int j1, j2; + unsigned int j3, j4; + Window root_return; + + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + + //printf("Width = %u, Height = %u, x = %i, y = %i\n", width, height, j1, j2); + //printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win); + + double x_step;// = (width * 0.8) / x.size(); + double y_step;// = (height * 0.8) / x.size(); + + + if(x_low_win == -1 and y_low_win == -1 and x_hi_win == -1 and y_hi_win == -1){ + // Audomatically decide data postition + max_x = x[0]; + max_y = y[0]; + min_x = x[0]; + min_y = y[0]; + + for(int i=0; i max_x) max_x = x[i]; + if(x[i] < min_x) min_x = x[i]; + if(y[i] > max_y) max_y = y[i]; + if(y[i] < min_y) min_y = y[i]; + } + //printf(" max_x = %f\n", max_x ); + //printf(" max_y = %f\n", max_y ); + //printf(" min_x = %f\n", min_x ); + //printf(" min_y = %f\n", min_y ); + + width_scale = (max_x - min_x) / (0.8 * width); + x_offset = 0.5 * ((max_x - min_x) / 0.8) - 0.5 * (min_x + max_x); + + height_scale = -1. * (max_y - min_y) / (0.8 * height); + y_offset = -0.5 * ((max_y - min_y) / 0.8) + -0.5 * (min_y + max_y); + + //printf("width scale = %f\n", width_scale); + //printf("x_offset = %f\n", x_offset); + //printf("height scale = %f\n", height_scale); + //printf("y_offset = %f\n", y_offset); + double axis_x = (0. + x_offset) / width_scale; + double axis_y = (0. + y_offset) / height_scale; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y); + + char axis_val[4]; + int w_step = width / 10; + 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)); + } + + 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)); + } + + + XSetForeground(disp, DefaultGC(disp,screen), 0); + double x_wid ; + double y_wid ; + double x_wid2; + double y_wid2; + + for(int i=0; i < x.size() - 2; i++){ + x_wid = (x[i] + x_offset) / width_scale; + y_wid = (y[i] + y_offset) / 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); + XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4); + } + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); + + } else { + //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win + + width_scale = (x_hi_win - x_low_win) / width; + x_offset = -1.0 * x_low_win; + + height_scale = (y_hi_win - y_low_win) / height; + y_offset = -1.0 * y_low_win; + + //printf("width scale = %f\n", width_scale); + //printf("x_offset = %f\n", x_offset); + //printf("height scale = %f\n", height_scale); + //printf("y_offset = %f\n", y_offset); + double axis_x = (0. + x_offset) / width_scale; + double axis_y = (0. + y_offset) / height_scale; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y); + + char axis_val[4]; + int w_step = width / 10; + 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)); + } + + 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)); + } + + + XSetForeground(disp, DefaultGC(disp,screen), 0); + double x_wid ; + double y_wid ; + double x_wid2; + double y_wid2; + + for(int i=0; i < x.size() - 2; i++){ + x_wid = (x[i] + x_offset) / width_scale; + y_wid = (y[i] + y_offset) / 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); + XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4); + } + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); + + + } + + old_xl = x_low_win; + old_yl = y_low_win; + old_xh = x_hi_win; + old_yh = y_hi_win; + + return 1; +} + + +int HistoGUI::Loop(){ + + bool MousePressed = false; + bool MousePressed2 = false; + + while (1) { + XNextEvent(disp, &evt); + if (evt.type == Expose) { + // XFillRectangle(disp, wind, DefaultGC(disp, screen), 20, 20, 10, 10); +// XDrawString (disp, wind, DefaultGC(disp, screen), 10, 50, msg, strlen(msg)); + DrawData(-1,-1,-1,-1); + } else if (evt.type == ButtonPress){ + /* store the mouse button coordinates in 'int' variables. */ + /* also store the ID of the window on which the mouse was */ + /* pressed. */ + int mouse_x = evt.xbutton.x; + int mouse_y = evt.xbutton.y; + old_mouse_x = mouse_x; + old_mouse_y = mouse_y; + + /* check which mouse button was pressed, and act accordingly. */ + if(evt.xbutton.button == Button1){ + /* draw a pixel at the mouse position. */ + XClearWindow(disp, wind); + DrawCrosshairs(mouse_x, mouse_y); + DrawData(old_xl, old_yl, old_xh, old_yh); + MousePressed = true; + printf("(%f, %f)\n", pos_x, pos_y); + } if(evt.xbutton.button == Button3){ + MousePressed2 = true; + } if(evt.xbutton.button == Button2){ + XClearWindow(disp, wind); + DrawData(old_xl, old_yl, old_xh, old_yh); + } + } else if (evt.type == ButtonRelease){ + /* store the mouse button coordinates in 'int' variables. */ + /* also store the ID of the window on which the mouse was */ + /* pressed. */ + int mouse_x = evt.xbutton.x; + int mouse_y = evt.xbutton.y; + //printf("Released button: %i, %i\n", mouse_x, mouse_y); + //printf("Pressed button: %i, %i\n", mouse_x, mouse_y); + /* check which mouse button was pressed, and act accordingly. */ + if(evt.xbutton.button == Button1){ + /* draw a pixel at the mouse position. */ + Zoom(mouse_x, mouse_y); + //DrawData(); + MousePressed = false; + } if(evt.xbutton.button == Button3){ + MousePressed2 = false; + } + } else if (evt.type == MotionNotify and MousePressed){ + + int mouse_x = evt.xmotion.x; + int mouse_y = evt.xmotion.y; + XClearWindow(disp, wind); + + DrawCrosshairs(old_mouse_x, old_mouse_y); + DrawCrosshairs(mouse_x, mouse_y); + DrawData(old_xl, old_yl, old_xh, old_yh); + + } else if (evt.type == MotionNotify and MousePressed2){ + + int mouse_x = evt.xmotion.x; + int mouse_y = evt.xmotion.y; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour_veryred.pixel); + XFillRectangle(disp, wind, DefaultGC(disp, screen), mouse_x -2, mouse_y -2, 4, 4); + + } else if (evt.type == KeyPress){ + + printf("Key pressed = %x\n", evt.xkey.keycode); + if(evt.xkey.keycode == 0x41 or evt.xkey.keycode == 0x39){ + XClearWindow(disp, wind); + DrawData(-1,-1,-1,-1); + + } else { + break; + } + } + } + + return 1; +} + + + + + + + + + + diff --git a/Graphical_Testing/.GUI.h.swp b/Graphical_Testing/.GUI.h.swp new file mode 100644 index 0000000..276bbf8 Binary files /dev/null and b/Graphical_Testing/.GUI.h.swp differ diff --git a/Graphical_Testing/ADT b/Graphical_Testing/ADT new file mode 100755 index 0000000..5e039f6 Binary files /dev/null and b/Graphical_Testing/ADT differ diff --git a/Graphical_Testing/GUI.h b/Graphical_Testing/GUI.h new file mode 100644 index 0000000..7388347 --- /dev/null +++ b/Graphical_Testing/GUI.h @@ -0,0 +1,566 @@ +#include "global.h" + +//To compile : g++ AD.cxx -o {Input Executable Name} -lX11 +// VERSION : 1.0 +// Created with thanks to Daniel Folds Holt - University of Cambridge Cavendish Lab + + + + +using namespace std; + + +class HistoGUI{ + + + public: + + HistoGUI(){} + + Display * disp; + Window wind; + XEvent evt; + int screen; + + XColor xcolour; + XColor xcolour_red; + XColor xcolour_veryred; + Colormap cmap; + + + std::vector x; + std::vector y; + std::vector y_errors; + + double max_x; + double min_x; + double max_y; + double min_y; + + double pos_x; + double pos_y; + double new_pos_x; + double new_pos_y; + + unsigned int width; + unsigned int height; + double window_width; + double window_height; + double width_scale; + double x_offset; + double height_scale; + double y_offset; + + + double A0; + double A2E; + double A4E; + int SetFit(double a, double b, double c){ + A0 = a; + A2E = b; + A4E = c; + fit_exists = true; + return 1; + } + + + double old_xl, old_xh, old_yl, old_yh; + double old_mouse_x, old_mouse_y; + + int Init(); + int SetData(std::vector a, std::vectorb); + int SetErrors(std::vector a); + int Loop(); + void Close(){ printf("Closing now!\n"); XCloseDisplay(disp); } + int DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win); + int Draw_Fit(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win); + int DrawCrosshairs(int mouse_x, int mouse_y); + int Zoom(int mouse_x, int mouse_y); + + double legval(double theta); + bool fit_exists = false; + + +// int Legendre_fit(std::vector d1, std::vector d2, double + +}; + +int HistoGUI::SetData(std::vector a, std::vectorb){ + for(int i=0; i < a.size(); i++){ + x.push_back(a[i]); + y.push_back(b[i]); + } + + return a.size(); +} + +int HistoGUI::SetErrors(std::vector a){ + for(int i=0; i < a.size(); i++){ + y_errors.push_back(a[i]); + } + return a.size(); +} +int HistoGUI::Init(){ + + disp = XOpenDisplay(NULL); + if (disp == NULL) { + fprintf(stderr, "Cannot open display\n"); + exit(1); + } + + //fit_exists = false; + + screen = DefaultScreen(disp); + wind = XCreateSimpleWindow(disp, RootWindow(disp, screen), 10, 10, 600, 400, 1, + BlackPixel(disp, screen), WhitePixel(disp, screen)); + XGrabPointer(disp, wind, False, ButtonPressMask | ButtonReleaseMask | PointerMotionMask, GrabModeAsync, + GrabModeAsync, None, None, CurrentTime); + XSelectInput(disp, wind, ExposureMask | KeyPressMask | ButtonPressMask | ButtonReleaseMask | PointerMotionMask); + XMapWindow(disp, wind); + + // colours + cmap = DefaultColormap(disp, screen); + xcolour.red = 32000; + xcolour.green = 32000; + xcolour.blue = 42000; + xcolour.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour); + + xcolour_red.red = 42000; + xcolour_red.green = 32000; + xcolour_red.blue = 32000; + xcolour_red.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour_red); + + xcolour_veryred.red = 65000; + xcolour_veryred.green = 32000; + xcolour_veryred.blue = 32000; + xcolour_veryred.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour_veryred); + + + + return 1; + +} + +int HistoGUI::DrawCrosshairs(int mouse_x, int mouse_y){ + + int j1, j2; + unsigned int j3, j4; + Window root_return; + + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + XSetForeground(disp, DefaultGC(disp, screen), xcolour_red.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), mouse_x, 0.0, mouse_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, mouse_y, width, mouse_y); + + pos_y = mouse_y * height_scale - y_offset; + pos_x = mouse_x * width_scale - x_offset; + + char coord[32]; + sprintf(coord, "(%.2f, %.2f)", pos_x, pos_y); + XDrawString(disp, wind, DefaultGC(disp, screen), mouse_x + 5, mouse_y + 12, coord, strlen(coord)); + + return 1; +} + +int HistoGUI::Zoom(int mouse_x, int mouse_y){ + + new_pos_x = mouse_x * width_scale - x_offset; + new_pos_y = mouse_y * height_scale - y_offset; + + pos_x = old_mouse_x * width_scale - x_offset; + pos_y = old_mouse_y * height_scale - y_offset; + + + if(std::abs(mouse_x - old_mouse_x) < 20 or std::abs(mouse_y - old_mouse_y) < 20) return 1; + + + //printf("Zooming: %f - %f, %f - %f\n", mouse_x, old_mouse_x,mouse_y,old_mouse_y); + + double x_low_win, y_low_win, x_hi_win, y_hi_win; + if(new_pos_x < pos_x){ + x_low_win = new_pos_x; + x_hi_win = pos_x; + //printf("new low; [%f, %f]\n",x_low_win, x_hi_win); + } else { + x_low_win = pos_x; + x_hi_win = new_pos_x; + //printf("new high; [%f, %f]\n",x_low_win, x_hi_win); + } + if(new_pos_y > pos_y){ + y_low_win = new_pos_y; + y_hi_win = pos_y; + } else { + y_low_win = pos_y; + y_hi_win = new_pos_y; + } + + + //printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win); + XClearWindow(disp, wind); + DrawData(x_low_win, y_low_win, x_hi_win, y_hi_win); + //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win + + return 1; + +} +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. ); + + lg = aaa + aab + aac; + + + return lg; +} + +int HistoGUI::Draw_Fit(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win){ + + + int j1, j2; + unsigned int j3, j4; + Window root_return; + + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + + //printf("Width = %u, Height = %u, x = %i, y = %i\n", width, height, j1, j2); + //printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win); + + double x_step;// = (width * 0.8) / x.size(); + double y_step;// = (height * 0.8) / x.size(); + + + if(x_low_win == -1 and y_low_win == -1 and x_hi_win == -1 and y_hi_win == -1){ + // Audomatically decide data postition + max_x = x[0]; + max_y = y[0]; + min_x = x[0]; + min_y = y[0]; + + for(int i=0; i max_x) max_x = x[i]; + if(x[i] < min_x) min_x = x[i]; + if(y[i] > max_y) max_y = y[i]; + if(y[i] < min_y) min_y = y[i]; + } + + width_scale = (max_x - min_x) / (0.8 * width); + x_offset = 0.5 * ((max_x - min_x) / 0.8) - 0.5 * (min_x + max_x); + + height_scale = -1. * (max_y - min_y) / (0.8 * height); + y_offset = -0.5 * ((max_y - min_y) / 0.8) + -0.5 * (min_y + max_y); + + XSetForeground(disp, DefaultGC(disp,screen), xcolour_red.pixel); + + double x_wid ; + double y_wid ; + + for(int i=0; i < (int) width; i += 2){ + double x_val = i * width_scale - x_offset; + double y_val = legval(x_val); + + y_wid = (y_val + y_offset) / height_scale; + XFillRectangle(disp, wind, DefaultGC(disp, screen), i -2, y_wid -2, 4, 4); + } + // XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); + + } else { + //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win + + width_scale = (x_hi_win - x_low_win) / width; + x_offset = -1.0 * x_low_win; + + height_scale = (y_hi_win - y_low_win) / height; + y_offset = -1.0 * y_low_win; + + + XSetForeground(disp, DefaultGC(disp,screen), xcolour_red.pixel); + + double x_wid ; + double y_wid ; + + for(int i=0; i < (int) width; i += 2){ + double x_val = i * width_scale - x_offset; + double y_val = legval(x_val); + + y_wid = (y_val + y_offset) / height_scale; + XFillRectangle(disp, wind, DefaultGC(disp, screen), i -2, y_wid -2, 4, 4); + } + } + + old_xl = x_low_win; + old_yl = y_low_win; + old_xh = x_hi_win; + old_yh = y_hi_win; + + + + return 1; +} + +int HistoGUI::DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win){ + + if(fit_exists){ + Draw_Fit(x_low_win,y_low_win,x_hi_win,y_hi_win); + + } + + int j1, j2; + unsigned int j3, j4; + Window root_return; + + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + + //printf("Width = %u, Height = %u, x = %i, y = %i\n", width, height, j1, j2); + //printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win); + + double x_step;// = (width * 0.8) / x.size(); + double y_step;// = (height * 0.8) / x.size(); + + + if(x_low_win == -1 and y_low_win == -1 and x_hi_win == -1 and y_hi_win == -1){ + // Audomatically decide data postition + max_x = x[0]; + max_y = y[0]; + min_x = x[0]; + min_y = y[0]; + + for(int i=0; i max_x) max_x = x[i]; + if(x[i] < min_x) min_x = x[i]; + if(y[i] > max_y) max_y = y[i]; + if(y[i] < min_y) min_y = y[i]; + } + //printf(" max_x = %f\n", max_x ); + //printf(" max_y = %f\n", max_y ); + //printf(" min_x = %f\n", min_x ); + //printf(" min_y = %f\n", min_y ); + + width_scale = (max_x - min_x) / (0.8 * width); + x_offset = 0.5 * ((max_x - min_x) / 0.8) - 0.5 * (min_x + max_x); + + height_scale = -1. * (max_y - min_y) / (0.8 * height); + y_offset = -0.5 * ((max_y - min_y) / 0.8) + -0.5 * (min_y + max_y); + + + //printf("width scale = %f\n", width_scale); + //printf("x_offset = %f\n", x_offset); + //printf("height scale = %f\n", height_scale); + //printf("y_offset = %f\n", y_offset); + double axis_x = (0. + x_offset) / width_scale; + double axis_y = (0. + y_offset) / height_scale; + + //cout << "x-axis = %lf\n" < +#include +#include +#include +#include +#include + + +class HistoGUI{ + + + public: + + HistoGUI(){} + + Display * disp; + Window wind; + XEvent evt; + int screen; + + XColor xcolour; + XColor xcolour_red; + XColor xcolour_veryred; + Colormap cmap; + + + std::vector x; + std::vector y; + double max_x; + double min_x; + double max_y; + double min_y; + + double pos_x; + double pos_y; + double new_pos_x; + double new_pos_y; + + unsigned int width; + unsigned int height; + double window_width; + double window_height; + double width_scale; + double x_offset; + double height_scale; + double y_offset; + + double old_xl, old_xh, old_yl, old_yh; + double old_mouse_x, old_mouse_y; + + int Init(); + int SetData(std::vector a, std::vectorb); + int Loop(); + void Close(){ printf("Closing now!\n"); XCloseDisplay(disp); } + int DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win); + int DrawCrosshairs(int mouse_x, int mouse_y); + int Zoom(int mouse_x, int mouse_y); + + +}; + +int HistoGUI::SetData(std::vector a, std::vectorb){ + for(int i=0; i < a.size(); i++){ + x.push_back(a[i]); + y.push_back(b[i]); + } + return a.size(); +} + +int HistoGUI::Init(){ + + disp = XOpenDisplay(NULL); + if (disp == NULL) { + fprintf(stderr, "Cannot open display\n"); + exit(1); + } + + screen = DefaultScreen(disp); + wind = XCreateSimpleWindow(disp, RootWindow(disp, screen), 10, 10, 600, 400, 1, + BlackPixel(disp, screen), WhitePixel(disp, screen)); + XGrabPointer(disp, wind, False, ButtonPressMask | ButtonReleaseMask | PointerMotionMask, GrabModeAsync, + GrabModeAsync, None, None, CurrentTime); + XSelectInput(disp, wind, ExposureMask | KeyPressMask | ButtonPressMask | ButtonReleaseMask | PointerMotionMask); + XMapWindow(disp, wind); + + // colours + cmap = DefaultColormap(disp, screen); + xcolour.red = 32000; + xcolour.green = 32000; + xcolour.blue = 42000; + xcolour.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour); + + xcolour_red.red = 42000; + xcolour_red.green = 32000; + xcolour_red.blue = 32000; + xcolour_red.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour_red); + + xcolour_veryred.red = 65000; + xcolour_veryred.green = 32000; + xcolour_veryred.blue = 32000; + xcolour_veryred.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour_veryred); + + + + return 1; + +} + +int HistoGUI::DrawCrosshairs(int mouse_x, int mouse_y){ + + int j1, j2; + unsigned int j3, j4; + Window root_return; + + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + XSetForeground(disp, DefaultGC(disp, screen), xcolour_red.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), mouse_x, 0.0, mouse_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, mouse_y, width, mouse_y); + + pos_y = mouse_y * height_scale - y_offset; + pos_x = mouse_x * width_scale - x_offset; + + char coord[32]; + sprintf(coord, "(%.2f, %.2f)", pos_x, pos_y); + XDrawString(disp, wind, DefaultGC(disp, screen), mouse_x + 5, mouse_y + 12, coord, strlen(coord)); + + return 1; +} + +int HistoGUI::Zoom(int mouse_x, int mouse_y){ + + new_pos_x = mouse_x * width_scale - x_offset; + new_pos_y = mouse_y * height_scale - y_offset; + + pos_x = old_mouse_x * width_scale - x_offset; + pos_y = old_mouse_y * height_scale - y_offset; + + + if(std::abs(mouse_x - old_mouse_x) < 20 or std::abs(mouse_y - old_mouse_y) < 20) return 1; + + + //printf("Zooming: %f - %f, %f - %f\n", mouse_x, old_mouse_x,mouse_y,old_mouse_y); + + double x_low_win, y_low_win, x_hi_win, y_hi_win; + if(new_pos_x < pos_x){ + x_low_win = new_pos_x; + x_hi_win = pos_x; + //printf("new low; [%f, %f]\n",x_low_win, x_hi_win); + } else { + x_low_win = pos_x; + x_hi_win = new_pos_x; + //printf("new high; [%f, %f]\n",x_low_win, x_hi_win); + } + if(new_pos_y > pos_y){ + y_low_win = new_pos_y; + y_hi_win = pos_y; + } else { + y_low_win = pos_y; + y_hi_win = new_pos_y; + } + + + //printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win); + XClearWindow(disp, wind); + DrawData(x_low_win, y_low_win, x_hi_win, y_hi_win); + //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win + + return 1; + +} + + +int HistoGUI::DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win){ + + + int j1, j2; + unsigned int j3, j4; + Window root_return; + + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + + //printf("Width = %u, Height = %u, x = %i, y = %i\n", width, height, j1, j2); + //printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win); + + double x_step;// = (width * 0.8) / x.size(); + double y_step;// = (height * 0.8) / x.size(); + + + if(x_low_win == -1 and y_low_win == -1 and x_hi_win == -1 and y_hi_win == -1){ + // Audomatically decide data postition + max_x = x[0]; + max_y = y[0]; + min_x = x[0]; + min_y = y[0]; + + for(int i=0; i max_x) max_x = x[i]; + if(x[i] < min_x) min_x = x[i]; + if(y[i] > max_y) max_y = y[i]; + if(y[i] < min_y) min_y = y[i]; + } + //printf(" max_x = %f\n", max_x ); + //printf(" max_y = %f\n", max_y ); + //printf(" min_x = %f\n", min_x ); + //printf(" min_y = %f\n", min_y ); + + width_scale = (max_x - min_x) / (0.8 * width); + x_offset = 0.5 * ((max_x - min_x) / 0.8) - 0.5 * (min_x + max_x); + + height_scale = -1. * (max_y - min_y) / (0.8 * height); + y_offset = -0.5 * ((max_y - min_y) / 0.8) + -0.5 * (min_y + max_y); + + //printf("width scale = %f\n", width_scale); + //printf("x_offset = %f\n", x_offset); + //printf("height scale = %f\n", height_scale); + //printf("y_offset = %f\n", y_offset); + double axis_x = (0. + x_offset) / width_scale; + double axis_y = (0. + y_offset) / height_scale; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y); + + char axis_val[4]; + int w_step = width / 10; + 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)); + } + + 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)); + } + + + XSetForeground(disp, DefaultGC(disp,screen), 0); + double x_wid ; + double y_wid ; + double x_wid2; + double y_wid2; + + for(int i=0; i < x.size() - 2; i++){ + x_wid = (x[i] + x_offset) / width_scale; + y_wid = (y[i] + y_offset) / 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); + XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4); + } + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); + + } else { + //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win + + width_scale = (x_hi_win - x_low_win) / width; + x_offset = -1.0 * x_low_win; + + height_scale = (y_hi_win - y_low_win) / height; + y_offset = -1.0 * y_low_win; + + //printf("width scale = %f\n", width_scale); + //printf("x_offset = %f\n", x_offset); + //printf("height scale = %f\n", height_scale); + //printf("y_offset = %f\n", y_offset); + double axis_x = (0. + x_offset) / width_scale; + double axis_y = (0. + y_offset) / height_scale; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y); + + char axis_val[4]; + int w_step = width / 10; + 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)); + } + + 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)); + } + + + XSetForeground(disp, DefaultGC(disp,screen), 0); + double x_wid ; + double y_wid ; + double x_wid2; + double y_wid2; + + for(int i=0; i < x.size() - 2; i++){ + x_wid = (x[i] + x_offset) / width_scale; + y_wid = (y[i] + y_offset) / 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); + XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4); + } + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); + + + } + + old_xl = x_low_win; + old_yl = y_low_win; + old_xh = x_hi_win; + old_yh = y_hi_win; + + return 1; +} + + +int HistoGUI::Loop(){ + + bool MousePressed = false; + bool MousePressed2 = false; + + while (1) { + XNextEvent(disp, &evt); + if (evt.type == Expose) { + // XFillRectangle(disp, wind, DefaultGC(disp, screen), 20, 20, 10, 10); +// XDrawString (disp, wind, DefaultGC(disp, screen), 10, 50, msg, strlen(msg)); + DrawData(-1,-1,-1,-1); + } else if (evt.type == ButtonPress){ + /* store the mouse button coordinates in 'int' variables. */ + /* also store the ID of the window on which the mouse was */ + /* pressed. */ + int mouse_x = evt.xbutton.x; + int mouse_y = evt.xbutton.y; + old_mouse_x = mouse_x; + old_mouse_y = mouse_y; + + /* check which mouse button was pressed, and act accordingly. */ + if(evt.xbutton.button == Button1){ + /* draw a pixel at the mouse position. */ + XClearWindow(disp, wind); + DrawCrosshairs(mouse_x, mouse_y); + DrawData(old_xl, old_yl, old_xh, old_yh); + MousePressed = true; + printf("(%f, %f)\n", pos_x, pos_y); + } if(evt.xbutton.button == Button3){ + MousePressed2 = true; + } if(evt.xbutton.button == Button2){ + XClearWindow(disp, wind); + DrawData(old_xl, old_yl, old_xh, old_yh); + } + } else if (evt.type == ButtonRelease){ + /* store the mouse button coordinates in 'int' variables. */ + /* also store the ID of the window on which the mouse was */ + /* pressed. */ + int mouse_x = evt.xbutton.x; + int mouse_y = evt.xbutton.y; + //printf("Released button: %i, %i\n", mouse_x, mouse_y); + //printf("Pressed button: %i, %i\n", mouse_x, mouse_y); + /* check which mouse button was pressed, and act accordingly. */ + if(evt.xbutton.button == Button1){ + /* draw a pixel at the mouse position. */ + Zoom(mouse_x, mouse_y); + //DrawData(); + MousePressed = false; + } if(evt.xbutton.button == Button3){ + MousePressed2 = false; + } + } else if (evt.type == MotionNotify and MousePressed){ + + int mouse_x = evt.xmotion.x; + int mouse_y = evt.xmotion.y; + XClearWindow(disp, wind); + + DrawCrosshairs(old_mouse_x, old_mouse_y); + DrawCrosshairs(mouse_x, mouse_y); + DrawData(old_xl, old_yl, old_xh, old_yh); + + } else if (evt.type == MotionNotify and MousePressed2){ + + int mouse_x = evt.xmotion.x; + int mouse_y = evt.xmotion.y; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour_veryred.pixel); + XFillRectangle(disp, wind, DefaultGC(disp, screen), mouse_x -2, mouse_y -2, 4, 4); + + } else if (evt.type == KeyPress){ + + printf("Key pressed = %x\n", evt.xkey.keycode); + if(evt.xkey.keycode == 0x41 or evt.xkey.keycode == 0x39){ + XClearWindow(disp, wind); + DrawData(-1,-1,-1,-1); + + } else { + break; + } + } + } + + return 1; +} + + + + + + + + + + diff --git a/Graphical_Testing/ad_test.cxx b/Graphical_Testing/ad_test.cxx new file mode 100644 index 0000000..1cc15f8 --- /dev/null +++ b/Graphical_Testing/ad_test.cxx @@ -0,0 +1,74 @@ +#include "global.h" +#include "GUI.h" + +int main ( int argc, char** argv){ + + HistoGUI gui; + + double A0E = 134.327; + double A2E = -11.7874; + double A4E = 0.760896; + + double step = 0.0001; + + vector Theta; + vector AD_I; + double ad_start = 0.;//-3.14159/2; + int adpoints = (3.14159)/(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) - 0.5); + aac = A4E*(35./8. * pow(theta,4) - 30./8. * pow(theta,2) + 3./8. ); + + Iad = aaa + aab + aac; + + cout << theta << "\n"; + + AD_I.push_back(Iad); + Theta.push_back(theta); + + + } + + vector dangle; + //= {0.785398,1.5708,2.35619}; + dangle.push_back(0.785398); + dangle.push_back(1.5708); + dangle.push_back(2.35619); + 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.); + vector deydata; + //= {10.,10.,10.}; + deydata.push_back(10.); + deydata.push_back(10.); + deydata.push_back(10.); + deydata.push_back(10.); + + + + gui.SetData(dangle,dydata); + gui.SetErrors(deydata); + gui.SetFit(A0E,A2E,A4E); + + gui.Init(); + gui.Loop(); + gui.Close(); + + + + + + + + return 0; +} diff --git a/Graphical_Testing/global.h b/Graphical_Testing/global.h new file mode 100644 index 0000000..5c10f13 --- /dev/null +++ b/Graphical_Testing/global.h @@ -0,0 +1,19 @@ +#ifndef __iamauniqueid_h__ +#define __iamauniqueid_h__ + +//define any global data type here. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#endif diff --git a/LiveLook_2d.cxx b/LiveLook_2d.cxx new file mode 100644 index 0000000..070f7fa --- /dev/null +++ b/LiveLook_2d.cxx @@ -0,0 +1,1858 @@ +// This is a file to read the data outize +//ut from a the TDC + +// == INCLUDES == +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define TDCsignalLen 15 +#define READ_LIMIT 10000 +#define EVENT_LIMIT 100000 +#define cursup "\033[A" + +class HistoGUI{ + + + public: + + HistoGUI(){} + + Display * disp; + Window wind; + XEvent evt; + int screen; + + XColor xcolour; + XColor xcolour_red; + XColor xcolour_veryred; + Colormap cmap; + + XColor PixelColour[10]; + + std::vector x; + std::vector y; + + std::vector x_2; + std::vector y_2; + std::vector> content_2; + + std::vector x_3; + std::vector y_3; + std::vector> content_3; + + double max_x; + double min_x; + double max_y; + double min_y; + + double pos_x; + double pos_y; + double new_pos_x; + double new_pos_y; + + unsigned int width; + unsigned int height; + double window_width; + double window_height; + double width_scale; + double x_offset; + double height_scale; + double y_offset; + bool Draw2D_On; + + double old_xl, old_xh, old_yl, old_yh; + double old_mouse_x, old_mouse_y; + + int Init(); + int SetData(std::vector a, std::vector b); + int SetData2(std::vector a, std::vector b, std::vector> c); + int Loop(); + void Close(){ printf("Closing now!\n"); XCloseDisplay(disp); } + int DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win); + int DrawData2D(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win); + int SwapData(); + + int DrawCrosshairs(int mouse_x, int mouse_y); + int Zoom(int mouse_x, int mouse_y); + long int ReReadChannel(char* filename, char* folderName, long int seek_to); + int TDCLoop(long int *locs, char filenames[16][64], char * folderName); + + std::map pixelmap = { + //ChNo PIXEL + // H - Type + { 0 , 52 ,}, + { 1 , 57 ,}, + { 2 , 60 ,}, + { 3 , 58 ,}, + { 4 , 59 ,}, + { 5 , 50 ,}, + { 6 , 49 ,}, + { 7 , 51 ,}, + { 8 , 36 ,}, + { 9 , 41 ,}, + { 10, 44 ,}, + { 11, 42 ,}, + { 12, 43 ,}, + { 13, 33 ,}, + { 14, 27 ,}, + { 15, 34 ,}, + { 16, 35 ,}, + { 17, 18 ,}, + { 18, 19 ,}, + { 19, 25 ,}, + { 20, 26 ,}, + { 21, 28 ,}, + { 22, 9 ,}, + { 23, 10 ,}, + { 24, 20 ,}, + { 25, 1 ,}, + { 26, 3 ,}, + { 27, 12 ,}, + { 28, 11 ,}, + { 29, 4 ,}, + { 30, 2 ,}, + { 31, 17 ,} + }; + + +}; + +int HistoGUI::SetData(std::vector a, std::vectorb){ + for(int i=0; i < a.size(); i++){ + x.push_back(a[i]); + y.push_back(b[i]); + } + return a.size(); +} + +int HistoGUI::SetData2(std::vector a, std::vector b, std::vector> c){ + for(int i=0; i < b.size(); i++){ + y_2.push_back(b[i]); + } + for(int j=0; j < a.size(); j++){ + x_2.push_back(a[j]); + std::vector temp; + for(int i=0; i < b.size(); i++){ + temp.push_back(c[j][i]); + } + content_2.push_back(temp); + } + return a.size(); +} + +int HistoGUI::SwapData(){ + std::vector temp_x = x_2; + std::vector temp_y = y_2; + std::vector> temp_cont = content_2; + + content_2 = content_3; + x_2 = x_3; + y_2 = y_3; + content_3 = temp_cont; + x_3 = temp_x; + y_3 = temp_y; + + return x_2.size(); +} + + +int HistoGUI::Init(){ + + Draw2D_On = false; + disp = XOpenDisplay(NULL); + if (disp == NULL) { + fprintf(stderr, "Cannot open display\n"); + exit(1); + } + + screen = DefaultScreen(disp); + wind = XCreateSimpleWindow(disp, RootWindow(disp, screen), 10, 10, 600, 400, 1, + BlackPixel(disp, screen), WhitePixel(disp, screen)); + XGrabPointer(disp, wind, False, ButtonPressMask | ButtonReleaseMask | PointerMotionMask, GrabModeAsync, + GrabModeAsync, None, None, CurrentTime); + XSelectInput(disp, wind, ExposureMask | KeyPressMask | ButtonPressMask | ButtonReleaseMask | PointerMotionMask); + XMapWindow(disp, wind); + + // colours + cmap = DefaultColormap(disp, screen); + xcolour.red = 32000; + xcolour.green = 32000; + xcolour.blue = 42000; + xcolour.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour); + + xcolour_red.red = 42000; + xcolour_red.green = 32000; + xcolour_red.blue = 32000; + xcolour_red.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour_red); + + xcolour_veryred.red = 65000; + xcolour_veryred.green = 32000; + xcolour_veryred.blue = 32000; + xcolour_veryred.flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &xcolour_veryred); + + for(int i = 0; i < 10; i++){ + PixelColour[i].red = 65535; + PixelColour[i].green = 65535 * (10 - i)/ 10; + PixelColour[i].blue = 65535 * (10 - i)/ 10; + PixelColour[i].flags = DoRed | DoGreen | DoBlue; + XAllocColor(disp, cmap, &PixelColour[i]); + } + + + return 1; + +} + +int HistoGUI::DrawCrosshairs(int mouse_x, int mouse_y){ + + int j1, j2; + unsigned int j3, j4; + Window root_return; + + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + XSetForeground(disp, DefaultGC(disp, screen), xcolour_red.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), mouse_x, 0.0, mouse_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, mouse_y, width, mouse_y); + + pos_y = mouse_y * height_scale - y_offset; + pos_x = mouse_x * width_scale - x_offset; + + char coord[32]; + sprintf(coord, "(%.2f, %.2f)", pos_x, pos_y); + XDrawString(disp, wind, DefaultGC(disp, screen), mouse_x + 5, mouse_y + 12, coord, strlen(coord)); + + return 1; +} + +int HistoGUI::Zoom(int mouse_x, int mouse_y){ + + new_pos_x = mouse_x * width_scale - x_offset; + new_pos_y = mouse_y * height_scale - y_offset; + + pos_x = old_mouse_x * width_scale - x_offset; + pos_y = old_mouse_y * height_scale - y_offset; + + + if(std::abs(mouse_x - old_mouse_x) < 20 or std::abs(mouse_y - old_mouse_y) < 20) return 1; + + + //printf("Zooming: %f - %f, %f - %f\n", mouse_x, old_mouse_x,mouse_y,old_mouse_y); + + double x_low_win, y_low_win, x_hi_win, y_hi_win; + if(new_pos_x < pos_x){ + x_low_win = new_pos_x; + x_hi_win = pos_x; + //printf("new low; [%f, %f]\n",x_low_win, x_hi_win); + } else { + x_low_win = pos_x; + x_hi_win = new_pos_x; + //printf("new high; [%f, %f]\n",x_low_win, x_hi_win); + } + if(new_pos_y > pos_y){ + y_low_win = new_pos_y; + y_hi_win = pos_y; + } else { + y_low_win = pos_y; + y_hi_win = new_pos_y; + } + + + //printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win); + XClearWindow(disp, wind); + DrawData(x_low_win, y_low_win, x_hi_win, y_hi_win); + //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win + + return 1; + +} + + + +int HistoGUI::DrawData2D(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win){ + int j1, j2; + unsigned int j3, j4; + Window root_return; + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + + double x_step;// = (width * 0.8) / x.size(); + double y_step;// = (height * 0.8) / x.size(); + + if(x_low_win == -1 and y_low_win == -1 and x_hi_win == -1 and y_hi_win == -1){ + // Audomatically decide data postition + max_x = x_2[0]; + max_y = y_2[0]; + min_x = x_2[0]; + min_y = y_2[0]; + + for(int i=0; i max_x) max_x = x_2[i]; + if(x_2[i] < min_x) min_x = x_2[i]; + } + for(int i=0; i max_y) max_y = y_2[i]; + if(y_2[i] < min_y) min_y = y_2[i]; + } + + double max_cont = content_2[0][0]; + + for(int i=0; i max_cont) max_cont = content_2[i][j]; + } + } + + //printf(" max_x = %f\n", max_x ); + //printf(" max_y = %f\n", max_y ); + //printf(" min_x = %f\n", min_x ); + //printf(" min_y = %f\n", min_y ); + + width_scale = (max_x - min_x) / (0.8 * width); + x_offset = 0.5 * ((max_x - min_x) / 0.8) - 0.5 * (min_x + max_x); + + height_scale = -1. * (max_y - min_y) / (0.8 * height); + y_offset = -0.5 * ((max_y - min_y) / 0.8) + -0.5 * (min_y + max_y); + + //printf("width scale = %f\n", width_scale); + //printf("x_offset = %f\n", x_offset); + //printf("height scale = %f\n", height_scale); + //printf("y_offset = %f\n", y_offset); + + XSetForeground(disp, DefaultGC(disp,screen), 0); + double x_wid ; + double y_wid ; + double x_wid2; + double y_wid2; + + double binwidth_x = 0.8 * width_scale / x_2.size(); + double binwidth_y = 0.8 * height_scale / y_2.size(); + + + for(int i=0; i < width; i++){ + x_wid = 1+ (i * width_scale) - x_offset; + if(x_wid > max_x) x_wid = max_x; + if(x_wid < min_x) x_wid = min_x; + + for(int j=0; j < height; j++){ + y_wid = 1 + (j * height_scale) - y_offset; + if(y_wid > max_y) y_wid = max_y; + if(y_wid < min_y) y_wid = min_y; + int colindex = (int) 10 * content_2[(int)x_wid][(int)y_wid] / max_cont; + //int colindex = (int) 10 * i / width; + //if(colindex > 0){ + //printf("Colour = %i\n",colindex); + //printf("(%f, %f) = %f\n", x_wid,y_wid, content_2[(int)x_wid][(int)y_wid]); + //} + XSetForeground(disp, DefaultGC(disp,screen), PixelColour[colindex].pixel); + //XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid - 0.5* binwidth_x, y_wid -0.5*binwidth_y, binwidth_x, binwidth_y); + XDrawPoint(disp, wind, DefaultGC(disp, screen), i,j); + } + } + + double axis_x = (0. + x_offset) / width_scale; + double axis_y = (0. + y_offset) / height_scale; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y); + + char axis_val[4]; + int w_step = width / 10; + 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)); + } + + 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)); + } + + } else { + //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win + + double max_cont = content_2[0][0]; + for(int i=0; i max_cont) max_cont = content_2[i][j]; + } + } + + width_scale = (x_hi_win - x_low_win) / width; + x_offset = -1.0 * x_low_win; + + height_scale = (y_hi_win - y_low_win) / height; + y_offset = -1.0 * y_low_win; + + //printf("width scale = %f\n", width_scale); + //printf("x_offset = %f\n", x_offset); + //printf("height scale = %f\n", height_scale); + //printf("y_offset = %f\n", y_offset); + + XSetForeground(disp, DefaultGC(disp,screen), 0); + double x_wid ; + double y_wid ; + double x_wid2; + double y_wid2; + + double binwidth_x = 0.8 * width_scale / x_2.size(); + double binwidth_y = 0.8 * height_scale / y_2.size(); + + for(int i=0; i < width; i++){ + x_wid = (i * width_scale) - x_offset; + if(x_wid > max_x) x_wid = max_x; + if(x_wid < min_x) x_wid = min_x; + + for(int j=0; j < height; j++){ + y_wid = (j * height_scale) - y_offset; + if(y_wid > max_y) y_wid = max_y; + if(y_wid < min_y) y_wid = min_y; + int colindex = (int) 10 * content_2[(int)x_wid][(int)y_wid] / max_cont; + //int colindex = (int) 10 * i / width; + //if(colindex > 0){ + //printf("Colour = %i\n",colindex); + //printf("(%f, %f) = %f\n", x_wid,y_wid, content_2[(int)x_wid][(int)y_wid]); + //} + XSetForeground(disp, DefaultGC(disp,screen), PixelColour[colindex].pixel); + //XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid - 0.5* binwidth_x, y_wid -0.5*binwidth_y, binwidth_x, binwidth_y); + XDrawPoint(disp, wind, DefaultGC(disp, screen), i,j); + } + } + + double axis_x = (0. + x_offset) / width_scale; + double axis_y = (0. + y_offset) / height_scale; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y); + + char axis_val[4]; + int w_step = width / 10; + 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)); + } + + 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)); + } + + + } + + old_xl = x_low_win; + old_yl = y_low_win; + old_xh = x_hi_win; + old_yh = y_hi_win; + + return 1; +} + + +int HistoGUI::DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win){ + + if(Draw2D_On){ + DrawData2D(x_low_win, y_low_win, x_hi_win, y_hi_win); + return 1; + } + + int j1, j2; + unsigned int j3, j4; + Window root_return; + + XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4); + + //printf("Width = %u, Height = %u, x = %i, y = %i\n", width, height, j1, j2); + //printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win); + + double x_step;// = (width * 0.8) / x.size(); + double y_step;// = (height * 0.8) / x.size(); + + + if(x_low_win == -1 and y_low_win == -1 and x_hi_win == -1 and y_hi_win == -1){ + // Audomatically decide data postition + max_x = x[0]; + max_y = y[0]; + min_x = x[0]; + min_y = y[0]; + + for(int i=0; i max_x) max_x = x[i]; + if(x[i] < min_x) min_x = x[i]; + if(y[i] > max_y) max_y = y[i]; + if(y[i] < min_y) min_y = y[i]; + } + //printf(" max_x = %f\n", max_x ); + //printf(" max_y = %f\n", max_y ); + //printf(" min_x = %f\n", min_x ); + //printf(" min_y = %f\n", min_y ); + + width_scale = (max_x - min_x) / (0.8 * width); + x_offset = 0.5 * ((max_x - min_x) / 0.8) - 0.5 * (min_x + max_x); + + height_scale = -1. * (max_y - min_y) / (0.8 * height); + y_offset = -0.5 * ((max_y - min_y) / 0.8) + -0.5 * (min_y + max_y); + + //printf("width scale = %f\n", width_scale); + //printf("x_offset = %f\n", x_offset); + //printf("height scale = %f\n", height_scale); + //printf("y_offset = %f\n", y_offset); + double axis_x = (0. + x_offset) / width_scale; + double axis_y = (0. + y_offset) / height_scale; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y); + + char axis_val[4]; + int w_step = width / 10; + 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)); + } + + 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)); + } + + + XSetForeground(disp, DefaultGC(disp,screen), 0); + double x_wid ; + double y_wid ; + double x_wid2; + double y_wid2; + + for(int i=0; i < x.size() - 2; i++){ + x_wid = (x[i] + x_offset) / width_scale; + y_wid = (y[i] + y_offset) / 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); + XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4); + } + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); + + } else { + //double x_low_win, double y_low_win, double x_hi_win, double y_hi_win + + width_scale = (x_hi_win - x_low_win) / width; + x_offset = -1.0 * x_low_win; + + height_scale = (y_hi_win - y_low_win) / height; + y_offset = -1.0 * y_low_win; + + //printf("width scale = %f\n", width_scale); + //printf("x_offset = %f\n", x_offset); + //printf("height scale = %f\n", height_scale); + //printf("y_offset = %f\n", y_offset); + double axis_x = (0. + x_offset) / width_scale; + double axis_y = (0. + y_offset) / height_scale; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel); + XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height); + XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y); + + char axis_val[4]; + int w_step = width / 10; + 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)); + } + + 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)); + } + + + XSetForeground(disp, DefaultGC(disp,screen), 0); + double x_wid ; + double y_wid ; + double x_wid2; + double y_wid2; + + for(int i=0; i < x.size() - 2; i++){ + x_wid = (x[i] + x_offset) / width_scale; + y_wid = (y[i] + y_offset) / 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); + XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2); + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4); + } + XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4); + + + } + + old_xl = x_low_win; + old_yl = y_low_win; + old_xh = x_hi_win; + old_yh = y_hi_win; + + return 1; +} + + +int HistoGUI::Loop(){ + + bool MousePressed = false; + bool MousePressed2 = false; + + while (1) { + XNextEvent(disp, &evt); + if (evt.type == Expose) { + // XFillRectangle(disp, wind, DefaultGC(disp, screen), 20, 20, 10, 10); +// XDrawString (disp, wind, DefaultGC(disp, screen), 10, 50, msg, strlen(msg)); + DrawData(-1,-1,-1,-1); + } else if (evt.type == ButtonPress){ + /* store the mouse button coordinates in 'int' variables. */ + /* also store the ID of the window on which the mouse was */ + /* pressed. */ + int mouse_x = evt.xbutton.x; + int mouse_y = evt.xbutton.y; + old_mouse_x = mouse_x; + old_mouse_y = mouse_y; + + /* check which mouse button was pressed, and act accordingly. */ + if(evt.xbutton.button == Button1){ + /* draw a pixel at the mouse position. */ + XClearWindow(disp, wind); + DrawData(old_xl, old_yl, old_xh, old_yh); + DrawCrosshairs(mouse_x, mouse_y); + MousePressed = true; + printf("(%f, %f)\n", pos_x, pos_y); + } if(evt.xbutton.button == Button3){ + MousePressed2 = true; + } if(evt.xbutton.button == Button2){ + XClearWindow(disp, wind); + DrawData(old_xl, old_yl, old_xh, old_yh); + } + } else if (evt.type == ButtonRelease){ + /* store the mouse button coordinates in 'int' variables. */ + /* also store the ID of the window on which the mouse was */ + /* pressed. */ + int mouse_x = evt.xbutton.x; + int mouse_y = evt.xbutton.y; + //printf("Released button: %i, %i\n", mouse_x, mouse_y); + //printf("Pressed button: %i, %i\n", mouse_x, mouse_y); + /* check which mouse button was pressed, and act accordingly. */ + if(evt.xbutton.button == Button1){ + /* draw a pixel at the mouse position. */ + Zoom(mouse_x, mouse_y); + //DrawData(); + MousePressed = false; + } if(evt.xbutton.button == Button3){ + MousePressed2 = false; + } + } else if (evt.type == MotionNotify and MousePressed){ + + int mouse_x = evt.xmotion.x; + int mouse_y = evt.xmotion.y; + XClearWindow(disp, wind); + + DrawData(old_xl, old_yl, old_xh, old_yh); + DrawCrosshairs(old_mouse_x, old_mouse_y); + DrawCrosshairs(mouse_x, mouse_y); + + } else if (evt.type == MotionNotify and MousePressed2){ + + int mouse_x = evt.xmotion.x; + int mouse_y = evt.xmotion.y; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour_veryred.pixel); + XFillRectangle(disp, wind, DefaultGC(disp, screen), mouse_x -2, mouse_y -2, 4, 4); + + } else if (evt.type == KeyPress){ + + if(evt.xkey.keycode == 0x41){ + XClearWindow(disp, wind); + DrawData(-1,-1,-1,-1); + } else if(evt.xkey.keycode == 0x27){ + SwapData(); + XClearWindow(disp, wind); + DrawData(old_xl, old_yl, old_xh, old_yh); + } else if(evt.xkey.keycode == 0x19 and !Draw2D_On){ + XClearWindow(disp, wind); + Draw2D_On = true; + DrawData2D(-1,-1,-1,-1); + } else if(evt.xkey.keycode == 0x19 and Draw2D_On){ + XClearWindow(disp, wind); + Draw2D_On = false; + DrawData(-1,-1,-1,-1); + } else { + break; + } + } + } + + return 1; +} + +int HistoGUI::TDCLoop(long int *locs, char filenames[16][64], char * folderName){ + printf("Started loop\n"); + fflush(stdout); + + bool MousePressed = false; + bool MousePressed2 = false; + + while (1) { + XNextEvent(disp, &evt); + if (evt.type == Expose) { + // XFillRectangle(disp, wind, DefaultGC(disp, screen), 20, 20, 10, 10); +// XDrawString (disp, wind, DefaultGC(disp, screen), 10, 50, msg, strlen(msg)); + DrawData(-1,-1,-1,-1); + } else if (evt.type == ButtonPress){ + /* store the mouse button coordinates in 'int' variables. */ + /* also store the ID of the window on which the mouse was */ + /* pressed. */ + int mouse_x = evt.xbutton.x; + int mouse_y = evt.xbutton.y; + old_mouse_x = mouse_x; + old_mouse_y = mouse_y; + + /* check which mouse button was pressed, and act accordingly. */ + if(evt.xbutton.button == Button1){ + /* draw a pixel at the mouse position. */ + XClearWindow(disp, wind); + DrawData(old_xl, old_yl, old_xh, old_yh); + DrawCrosshairs(mouse_x, mouse_y); + MousePressed = true; + printf("(%f, %f)\n", pos_x, pos_y); + } if(evt.xbutton.button == Button3){ + MousePressed2 = true; + } if(evt.xbutton.button == Button2){ + XClearWindow(disp, wind); + DrawData(old_xl, old_yl, old_xh, old_yh); + } + } else if (evt.type == ButtonRelease){ + /* store the mouse button coordinates in 'int' variables. */ + /* also store the ID of the window on which the mouse was */ + /* pressed. */ + int mouse_x = evt.xbutton.x; + int mouse_y = evt.xbutton.y; + //printf("Released button: %i, %i\n", mouse_x, mouse_y); + //printf("Pressed button: %i, %i\n", mouse_x, mouse_y); + /* check which mouse button was pressed, and act accordingly. */ + if(evt.xbutton.button == Button1){ + /* draw a pixel at the mouse position. */ + Zoom(mouse_x, mouse_y); + //DrawData(); + MousePressed = false; + } if(evt.xbutton.button == Button3){ + MousePressed2 = false; + } + } else if (evt.type == MotionNotify and MousePressed){ + + int mouse_x = evt.xmotion.x; + int mouse_y = evt.xmotion.y; + XClearWindow(disp, wind); + + DrawData(old_xl, old_yl, old_xh, old_yh); + DrawCrosshairs(old_mouse_x, old_mouse_y); + DrawCrosshairs(mouse_x, mouse_y); + + } else if (evt.type == MotionNotify and MousePressed2){ + + int mouse_x = evt.xmotion.x; + int mouse_y = evt.xmotion.y; + + XSetForeground(disp, DefaultGC(disp,screen), xcolour_veryred.pixel); + XFillRectangle(disp, wind, DefaultGC(disp, screen), mouse_x -2, mouse_y -2, 4, 4); + + } else if (evt.type == KeyPress){ + printf("Key = %x\n", evt.xkey.keycode); + + if(evt.xkey.keycode == 0x41){ + XClearWindow(disp, wind); + DrawData(-1,-1,-1,-1); + } else if(evt.xkey.keycode == 0x27){ + SwapData(); + XClearWindow(disp, wind); + DrawData(old_xl, old_yl, old_xh, old_yh); + } else if(evt.xkey.keycode == 0x19 and !Draw2D_On){ + XClearWindow(disp, wind); + Draw2D_On = true; + DrawData2D(-1,-1,-1,-1); + } else if(evt.xkey.keycode == 0x19 and Draw2D_On){ + XClearWindow(disp, wind); + Draw2D_On = false; + DrawData(-1,-1,-1,-1); + + } else if(evt.xkey.keycode == 0x18) { + break; + } else { + printf("Going\n"); + fflush(stdout); + for(int file=0; file < 16; file++){ + printf("%s : %f = %li\n",filenames[file],folderName, locs[file]); + locs[file] = ReReadChannel(filenames[file], folderName, locs[file]); + } + XClearWindow(disp, wind); + DrawData(old_xl, old_yl, old_xh, old_yh); + } + } + } + + return 1; +} + + +int number_of_zero_len_payloads = 0; + +// ======== ====== ======== +// ======== HEADER ======== +// ======== ====== ======== + +class HeaderTDC{ + + public: + + // Public variables + int ChNo; + + + // Constructor & destructor events + HeaderTDC(int a){ + ChNo = a; + } +// ~HeaderTDC(); + + // Read data from a file + int Read(FILE *infile){ + + int BR = fread(&header_data, sizeof(unsigned char), 4, infile); + return BR; + } + + // Print data to screen + int Print(){ + printf("Ch %i header: %.2x%.2x %.2x%.2x\n", ChNo, header_data[0], header_data[1], header_data[2], header_data[3] ); + return 1; + } + + + private: + // Private variables + unsigned char header_data[4]; + +}; + + + + +// ======== ===== ======== +// ======== EVENT ======== +// ======== ===== ======== + +// This is the event class, it holds an event. +class EventTDC{ + + public: + + // Constructor & destructor events + EventTDC(int a){ + id = a; + verbose = 0; + payloadLength = 0; + noTriggers = 0; + } + EventTDC(int a, int b){ + id = a; + channel_id = b; + verbose = 0; + Ch1_ready = 1; + Ch2_ready = 1; + payloadLength = 0; + noTriggers = 0; + } + EventTDC(){ + id = 1111; + channel_id = 1111; + verbose = 0; + Ch1_ready = 1; + Ch2_ready = 1; + payloadLength = 0; + noTriggers = 0; + } +// ~EventTDC(); + + // Public functions + int Read(FILE *infile); + int Print(); + int SetVerbose(int ch); + int size(int ch); + int GetID(){ return id;}; + int GetChID(){ return channel_id;}; + int SetChID(int a){channel_id = a; return channel_id;}; + int IsVerbose(){ return verbose;}; + int SetChIDFromFilename(char* filename); + + // Data vectors + std::vector fineTime_1; + std::vector fineTime_2; + std::vector tot_1; + std::vector tot_2; + std::vector trigNo_1; + std::vector trigNo_2; + + int DataFrame1[24 * TDCsignalLen] = { 0 }; + int DataFrame2[24 * TDCsignalLen] = { 0 }; + int payloadLength = 0; + int noTriggers; + int tot_1_candidate = 0; + int tot_2_candidate = 0; + unsigned long int id = 0; + + + private: + + // Private data + int searchTerm_1; + int delay_1; + int searchTerm_2; + int delay_2; + int channel_id; + int verbose; + int Ch1_ready; + int Ch2_ready; + + unsigned char bin_data[16]; + unsigned char channel_data[2]; + unsigned char channel_data_A[2]; + + // Private functions + int GetTDCTest(unsigned char data[4]); + int GetTDCTime1(); + int GetTDCTime2(); + + +}; + + +int EventTDC::Print(){ + + printf("\n"); + printf("Event %.4lu Ch no : %.2i - %.2i\n", id, channel_id, channel_id+1); + printf("==========================\n"); + printf("Payload length = %i\n", payloadLength); + printf("Bytes read = %i\n", 16 * payloadLength); + printf("\n"); + + printf("Channel %i: %i\n",channel_id, size(1)); + printf("==================\n"); + printf("Evt# | Time | ToT\n"); + + // Print channel 1 data - not too much mind + if(size(1) > 5){ + for(int i=0; i < 5; i++) printf("%.4i | %.4i | %.4i\n", i, fineTime_1[i], tot_1[i]); + printf(" . . . . . .\n"); + printf("%.4i | %.4i | %.4i\n", (size(1) - 1), fineTime_1[size(1) - 1], tot_1[size(1) - 1]); + } else { + for(int i=0; i < size(1); i++) printf("%.4i | %.4i | %.4i\n", i, fineTime_1[i], tot_1[i]); + } + printf("==================\n"); + + printf("\n"); + printf("Channel %i: %i\n", channel_id + 1, size(2)); + printf("==================\n"); + printf("Evt# | Time | ToT\n"); + + // Print channel 2 data + if(size(2) > 5){ + for(int i=0; i < 5; i++) printf("%.4i | %.4i | %.4i\n", i, fineTime_2[i], tot_2[i]); + printf(" . . . . . .\n"); + printf("%.4i | %.4i | %.4i\n", (size(2) - 1), fineTime_2[size(2) - 1], tot_2[size(2) - 1]); + } else { + for(int i=0; i < size(2); i++) printf("%.4i | %.4i | %.4i\n", i, fineTime_2[i], tot_2[i]); + } + printf("==================\n"); + printf("\n==========================\n"); + printf("\n"); + + return 1; +} + + + +// Set the event to read verbose Ch = 1, 2, or 3 (both) +int EventTDC::SetVerbose(int ch){ + verbose = ch; + return 1; +} + + + + +// Return the number of signals in one trigger Ch = 1, or 2 +int EventTDC::size(int ch){ + int ret = 0; + if(ch == 1){ + if(fineTime_1.size() == 0){ + ret = 0; + }else if(tot_1.size() == 0){ + ret = 0; + } else { ret = (int) fineTime_1.size();} + } else if(ch == 2){ + if(fineTime_2.size() == 0){ + ret = 0; + }else if(tot_2.size() == 0){ + ret = 0; + } else { ret = (int) fineTime_2.size();} + } + return ret; +} + + + + +// Search within the array for the search term - Ch 1 +int EventTDC::GetTDCTime1(){ + + int count = 0; + bin_data[0] = '\0'; + + // Assign the binary data values using a bit mask + int j = 0; + for(int i=1; i >=0; i--){ + bin_data[(8 * i) + 0] = !!(channel_data[j] & 0b10000000); + bin_data[(8 * i) + 1] = !!(channel_data[j] & 0b01000000); + bin_data[(8 * i) + 2] = !!(channel_data[j] & 0b00100000); + bin_data[(8 * i) + 3] = !!(channel_data[j] & 0b00010000); + + bin_data[(8 * i) + 4] = !!(channel_data[j] & 0b00001000); + bin_data[(8 * i) + 5] = !!(channel_data[j] & 0b00000100); + bin_data[(8 * i) + 6] = !!(channel_data[j] & 0b00000010); + bin_data[(8 * i) + 7] = !!(channel_data[j] & 0b00000001); + j = j+1; + } + + // Search the binary data for the "search term" i.e. a change in value + for(int i=15; i >= 4; i--){ + if(bin_data[i] != searchTerm_1){ + count += 1; + } else { + delay_1 += count; + + if(searchTerm_1 == 1){ + if(Ch1_ready == 1){ + fineTime_1.push_back(delay_1); + Ch1_ready = 0; + } + delay_1 = 0; + count = 1; + searchTerm_1 = 0; + + // When searching for the end of the signal + } else if (searchTerm_1 == 0){ + if(Ch1_ready == 0){ + tot_1_candidate += delay_1; + } + delay_1 = 0; + count = 1; + searchTerm_1 = 1; + } + } + } + + // If event is verbose then print the data + if(verbose == 1 or verbose == 3){ + printf("Ch1 = [ "); + for(int j = 0; j < 4; j++) printf("%u", bin_data[j]); + printf(" "); + for(int j = 4; j < 8; j++) printf("%u", bin_data[j]); + printf(" "); + for(int j = 8; j < 12; j++) printf("%u", bin_data[j]); + printf(" "); + for(int j = 12; j < 16; j++) printf("%u", bin_data[j]); + printf("] :: %i :: %i\n", delay_1, (delay_1/12)); + } + + // Update the count value + count = delay_1 + count; + delay_1 = count; + + channel_data[0] = '\0'; + return 1; +} + + + +// Search within the array for the search term - Ch 2 +int EventTDC::GetTDCTime2(){ + + int count = 0; + bin_data[0] = '\0'; + + // Assign the binary data values using a bit mask + int j = 0; + for(int i=1; i >=0; i--){ + bin_data[(8 * i) + 0] = !!(channel_data[j] & 0b10000000); + bin_data[(8 * i) + 1] = !!(channel_data[j] & 0b01000000); + bin_data[(8 * i) + 2] = !!(channel_data[j] & 0b00100000); + bin_data[(8 * i) + 3] = !!(channel_data[j] & 0b00010000); + + bin_data[(8 * i) + 4] = !!(channel_data[j] & 0b00001000); + bin_data[(8 * i) + 5] = !!(channel_data[j] & 0b00000100); + bin_data[(8 * i) + 6] = !!(channel_data[j] & 0b00000010); + bin_data[(8 * i) + 7] = !!(channel_data[j] & 0b00000001); + j = j+1; + } + + + // Search the binary data for the "search term" i.e. a change in value + for(int i=15; i >= 4; i--){ + if(bin_data[i] != searchTerm_2){ + count += 1; + } else { + // If the data includes the search term save the value + delay_2 += count; + + // When searching for the initial time + if(searchTerm_2 == 1){ + if(Ch2_ready == 1){ + fineTime_2.push_back(delay_2); + Ch2_ready = 0; + } + delay_2 = 0; + count = 1; + searchTerm_2 = 0; + + // When searching for the end of the signal + } else if (searchTerm_2 == 0){ + if(Ch2_ready == 0){ + tot_2_candidate += delay_2; + } + delay_2 = 0; + count = 1; + searchTerm_2 = 1; + } + } + } + + // If event is verbose then print the data + if(verbose == 2 or verbose == 3){ + printf("Ch2 = [ "); + for(int j = 0; j < 4; j++) printf("%u", bin_data[j]); + printf(" "); + for(int j = 4; j < 8; j++) printf("%u", bin_data[j]); + printf(" "); + for(int j = 8; j < 12; j++) printf("%u", bin_data[j]); + printf(" "); + for(int j = 12; j < 16; j++) printf("%u", bin_data[j]); + printf("] :: %i :: %i\n", delay_2, (delay_2/12)); + } + + // If the data includes the search term save the value + count = delay_2 + count; + delay_2 = count; + + channel_data[0] = '\0'; + return 1; +} + + + + + + +// Read an event from a file given +int EventTDC::Read(FILE *infile){ + + int BytesRead = 0; + + // Read the payload length + payloadLength = 0; + BytesRead = fread(&payloadLength, sizeof(short int), 1, infile); + BytesRead -= 1; + + // Set starting search terms + searchTerm_1 = 1; + delay_1 = 0; + searchTerm_2 = 1; + delay_2 = 0; + Ch1_ready = 1; + Ch2_ready = 1; + noTriggers = 0; + + // Read the payload + int LoopLen = payloadLength; + for(int i = 0; i < LoopLen; i++){ + + + //printf("i TDCsignalLen = %i\n", (i%TDCsignalLen)); + if(i%TDCsignalLen == 0){ + searchTerm_1 = 1; + delay_1 = 0; + searchTerm_2 = 1; + delay_2 = 0; + Ch1_ready = 1; + Ch2_ready = 1; + noTriggers += 1; + tot_1_candidate = 0; + tot_2_candidate = 0; + } + + // Read 32 bits of data + BytesRead += fread(&channel_data_A, sizeof(unsigned char), 2, infile); + BytesRead -= 2; + BytesRead += fread(&channel_data, sizeof(unsigned char), 2, infile); + BytesRead -= 2; + + + // Swap word order. + GetTDCTime2(); + + + for(int bdc = 0; bdc < 12; bdc++){ + DataFrame2[(i%TDCsignalLen)*24 + bdc] += bin_data[(15 - bdc)]; + } + + memcpy(channel_data, channel_data_A, 2); + GetTDCTime2(); + for(int bdc = 0; bdc < 12; bdc++){ + DataFrame2[(i%TDCsignalLen)*24 + bdc + 12] += bin_data[(15 - bdc)]; + } + + channel_data_A[0] = '\0'; + + // Read 32 bits of data + BytesRead += fread(&channel_data_A, sizeof(unsigned char), 2, infile); + BytesRead -= 2; + BytesRead += fread(&channel_data, sizeof(unsigned char), 2, infile); + BytesRead -= 2; + + GetTDCTime1(); + for(int bdc = 0; bdc < 12; bdc++){ + DataFrame1[(i%TDCsignalLen)*24 + bdc] += bin_data[(15 - bdc)]; + } + memcpy(channel_data, channel_data_A, 2); + GetTDCTime1(); + for(int bdc = 0; bdc < 12; bdc++){ + DataFrame1[(i%TDCsignalLen)*24 + bdc + 12] += bin_data[(15 - bdc)]; + } + + + if(i%TDCsignalLen == TDCsignalLen-1){ + if(i > 0){ + id = 0; + BytesRead += fread(&id, sizeof(unsigned int), 1, infile); + BytesRead -= 1; + trigNo_2.push_back(id); + BytesRead += fread(&id, sizeof(unsigned int), 1, infile); + BytesRead -= 1; + trigNo_1.push_back(id); + LoopLen -= 1; + } + } + channel_data_A[0] = '\0'; + + if(i%TDCsignalLen == TDCsignalLen-1){ + + if(Ch1_ready == 0) tot_1.push_back(tot_1_candidate); + if(Ch2_ready == 0) tot_2.push_back(tot_2_candidate); + + if(Ch1_ready == 1 and delay_1 > 0 and searchTerm_1 == 1){ + fineTime_1.push_back(-1); + tot_1.push_back(-1); + } else if(Ch1_ready == 0 and fineTime_1.size() != tot_1.size()){ + tot_1.push_back(delay_1); + } + + if(Ch2_ready == 1 and delay_2 > 0 and searchTerm_2 == 1){ + fineTime_2.push_back(-1); + tot_2.push_back(-1); + } else if(Ch2_ready == 0 and fineTime_2.size() != tot_2.size()){ + tot_2.push_back(delay_2); + } + + } + + // Set latch incase of runaway + if(delay_1 > READ_LIMIT or delay_2 > READ_LIMIT){ + printf("READ_LIMIT reached!\n"); + printf("If this is too low change in file or investigate in verbose mode.\n"); + break; + } + + } + + return BytesRead; +} + + + + + + +// Read a single file +std::vector ReadChannel(char* filename, char* folderName, long int &file_loc){ + + // Define basepath for reading in a folder + char basepath[128]; + basepath[0] = '\0'; + strcat(basepath, folderName); + strcat(basepath, "/DaqData/"); + strcat(basepath, filename); + + //============================================= Open File + FILE *infile = fopen(basepath, "rb"); + printf("Opened file %s\n", basepath); + if(infile == NULL){ + printf("Can't open file '%s'!\n", filename); + } + printf("Opened file %s\n", filename); + + //============================================= Read Header + int BytesTotal = 0; + int BytesRead = 0; + + HeaderTDC header = HeaderTDC(0); // HeaderTDC( Channel number ) + BytesRead = header.Read(infile); + header.Print(); + + //============================================= Read Data + std::vector events; + + int ch_no; + char hold[2], hold2[1] = {'\0'}; + int offset = strlen(filename); + char hold3[1] = {'\0'}; + hold3[0] = filename[offset - 6]; + + if(strcmp(hold3, "0123456789") == 1){ + hold[0] = filename[offset - 6]; + hold[1] = filename[offset - 5]; + ch_no = atoi((char*)hold); + } else { + hold2[0] = filename[offset - 5]; + ch_no = atoi(hold2); + } + + ch_no = ch_no * 2; + // Iterate through events + for(int i=0; i <= EVENT_LIMIT; i++){ + + // Read events + EventTDC event = EventTDC(i, ch_no); // EventTDC( id , ch_no) (int) + event.SetChID(ch_no); + + BytesRead = event.Read(infile); + + if(BytesRead != 0){ + printf("Bytes finished at %i\n", i); + break; + } + + if(event.payloadLength == 0){ + event.noTriggers = 15; + number_of_zero_len_payloads += 1; + } + + events.push_back(event); + + // Set latch incase of runaway + if(i > EVENT_LIMIT){ + printf("EVENT_LIMIT reached!\n"); + printf("If this is too low change in file or investigate in verbose mode.\n"); + break; + } + } + + printf("Read file\n"); + + file_loc = ftell(infile); + fclose(infile); + + return events; +} + +// Read a single file +long int HistoGUI::ReReadChannel(char* filename, char* folderName, long int seek_to){ + + // Define basepath for reading in a folder + char basepath[128]; + basepath[0] = '\0'; + strcat(basepath, folderName); + strcat(basepath, "/DaqData/"); + strcat(basepath, filename); + + //============================================= Open File + FILE *infile = fopen(basepath, "rb"); +// printf("Opened file %s\n", basepath); + if(infile == NULL){ + printf("Can't open file '%s'!\n", filename); + } +// printf("Opened file %s\n", filename); + fseek(infile, seek_to, SEEK_SET); + + //============================================= Read Header - NO NEED + int BytesTotal = 0; + int BytesRead = 0; + + + //============================================= Read Data + + int ch_no; + char hold[2], hold2[1] = {'\0'}; + int offset = strlen(filename); + char hold3[1] = {'\0'}; + hold3[0] = filename[offset - 6]; + + if(strcmp(hold3, "0123456789") == 1){ + hold[0] = filename[offset - 6]; + hold[1] = filename[offset - 5]; + ch_no = atoi((char*)hold); + } else { + hold2[0] = filename[offset - 5]; + ch_no = atoi(hold2); + } + + ch_no = ch_no * 2; + // Iterate through events + for(int i=0; i <= EVENT_LIMIT; i++){ + + // Read events + EventTDC event = EventTDC(i, ch_no); // EventTDC( id , ch_no) (int) + event.SetChID(ch_no); + + BytesRead = event.Read(infile); + + if(BytesRead != 0){ + printf("Bytes finished at %i\n", i); + break; + } + + if(event.payloadLength == 0){ + event.noTriggers = 15; + number_of_zero_len_payloads += 1; + } + + + for(int j=0; j < event.size(1); j++){ + if(event.fineTime_1[j] > 0){ + y[event.fineTime_1[j]] += 1; + content_2[event.fineTime_1[j]][event.tot_1[j]] += 1; + content_3[div(64 - pixelmap[ch_no], 8).rem -4 + 1][div(64 - pixelmap[ch_no], 8).quot + 1] += 1; + } + if(event.fineTime_2[j] > 0){ + y[event.fineTime_2[j]] += 1; + content_2[event.fineTime_2[j]][event.tot_2[j]] += 1; + content_3[div(64 - pixelmap[ch_no +1], 8).rem -4 + 1][div(64 - pixelmap[ch_no +1], 8).quot + 1] += 1; + } + } + + // Set latch incase of runaway + if(i > EVENT_LIMIT){ +// printf("EVENT_LIMIT reached!\n"); +// printf("If this is too low change in file or investigate in verbose mode.\n"); + break; + } + } + + printf("Read file\n"); + + long int file_loc = ftell(infile); + fclose(infile); + + return file_loc; +} + + +// ======= ==== ======== +// ======== MAIN ======== +// ======== ==== ======== + + +#ifndef __CINT__ + +int main(int argc,char **argv){ + + + // READ FOLDER ============================================================= + DIR *folder; + char fpath[128]; + sprintf(fpath, "%s/DaqData", argv[1]); + folder = opendir(fpath); + + // Error message + if(folder == NULL){ + printf("Unable to read '%s'!\n",argv[1]); + return(1); + } else { + printf("Reading from '%s'.\n", argv[1]); + } + + struct dirent *entry; + int folderStart = telldir(folder); + + int files = 0; + while( (entry=readdir(folder))) { + if(strcmp(entry->d_name, ".") == 0) continue; + if(strcmp(entry->d_name, "..") == 0) continue; + files++; + } + + // Return to the top of the file + seekdir(folder,folderStart); + + + // MAP ============================================================= + std::map pixelmap = { + //ChNo PIXEL + // H - Type + { 0 , 52 ,}, + { 1 , 57 ,}, + { 2 , 60 ,}, + { 3 , 58 ,}, + { 4 , 59 ,}, + { 5 , 50 ,}, + { 6 , 49 ,}, + { 7 , 51 ,}, + { 8 , 36 ,}, + { 9 , 41 ,}, + { 10, 44 ,}, + { 11, 42 ,}, + { 12, 43 ,}, + { 13, 33 ,}, + { 14, 27 ,}, + { 15, 34 ,}, + { 16, 35 ,}, + { 17, 18 ,}, + { 18, 19 ,}, + { 19, 25 ,}, + { 20, 26 ,}, + { 21, 28 ,}, + { 22, 9 ,}, + { 23, 10 ,}, + { 24, 20 ,}, + { 25, 1 ,}, + { 26, 3 ,}, + { 27, 12 ,}, + { 28, 11 ,}, + { 29, 4 ,}, + { 30, 2 ,}, + { 31, 17 ,} + }; + + double RisingEdgesArray[24 * TDCsignalLen][32] = {0}; + std::vector RisingEdgesVect(24 * TDCsignalLen); + double HitMapArray[8][8] = {0}; + double ChannelOccupancy[32] = {0}; + + + // LOOP TRHOUGH FILES ======================================================== + + // Counters --------- + int currentFile = 0; + int noRisingEdges = 0; + std::vector> All_channels; //(32); + + printf("READING ======================================\n\n"); + long int file_loc[16] = {0}; + + + + HistoGUI gui; + for(int cn=0; cn < 24 * TDCsignalLen; cn++){ + gui.x.push_back(cn); + gui.y.push_back(0); + + gui.x_2.push_back(cn); + gui.y_2.push_back(cn); + std::vector temp; + for(int ca=0; ca < 24 * TDCsignalLen; ca++){ + temp.push_back(0); + } + gui.content_2.push_back(temp); + } + for(int cn=0; cn < 8; cn++){ + gui.x_3.push_back(cn); + gui.y_3.push_back(cn); + + std::vector temp; + for(int ca=0; ca < 8; ca++){ + temp.push_back(0); + } + gui.content_3.push_back(temp); + } + char all_filenames[16][64]; + + // Loop ------------- + while( (entry=readdir(folder))) { + if(strcmp(entry->d_name, ".") == 0) continue; + if(strcmp(entry->d_name, "..") == 0) continue; + + // GET CHANNEL NAME ======================================= + int ch_no = 0; + char hold[3], hold2[2]; + hold[2] = '\0'; + hold2[1] = '\0'; + int offset = strlen(entry->d_name); + char hold3[1] = {'\0'}; + hold3[0] = entry->d_name[offset - 6]; + + if(strcmp(hold3, "0123456789") == 1){ + hold[0] = entry->d_name[offset - 6]; + hold[1] = entry->d_name[offset - 5]; + ch_no = atoi((char*)hold); + // printf("Channel is double digit: %c %c\n", entry->d_name[offset - 6], entry->d_name[offset - 5]); + } else { + // printf("Channel is single digit: %c\n", entry->d_name[offset - 5]); + hold2[0] = entry->d_name[offset - 5]; + ch_no = atoi(hold2); + // printf("hold2 = %s -> %i\n", hold2, ch_no); + } + + ch_no = ch_no * 2; + + // Read File ------------ + std::vector events = ReadChannel(entry->d_name, argv[1], file_loc[(ch_no / 2)]); + sprintf(all_filenames[(ch_no / 2)], "%s", entry->d_name); + std::vector channel; + printf("%s : Loc = %i\n", entry->d_name, file_loc[(ch_no / 2)]); + + // Alert if empty ------- + if(events.size() == 0){ + printf("FILE %s is empty!\n", entry->d_name); + continue; + } + + // ITERATE TRHOUGH DATA ======================================= + unsigned long long int count_ch1 = 0; + unsigned long long int count_ch2 = 0; + int payload_sum = 0; + int noTriggers_channel = 0; + int noTriggers_channel_2 = 0; + + printf("Ch No = %i\n", ch_no); + printf("No Payloads = %lu\n", events.size()); + + // Loop -------------- + for(int i = 0; i < events.size(); i++){ + payload_sum += events[i].payloadLength; + channel.push_back(events[i].payloadLength); + noTriggers_channel += events[i].noTriggers; + noTriggers_channel_2 += events[i].size(1); + if(events.size() == 0) break; + //printf("%i : %i\n", i, events[i].payloadLength); + + // Loop through triggers -------------- + for(int j=0; j < events[i].size(1); j++){ + + // If ToA <= 0, set to -1 & fill tree + if(events[i].fineTime_1[j] <= 0){ + continue; + } else { + + count_ch1 += 1; + HitMapArray[div(64 - pixelmap[ch_no], 8).rem -4 + 1][div(64 - pixelmap[ch_no], 8).quot + 1] += 1;//, (1.0/noTriggers)); + //if(id) + RisingEdgesArray[events[i].fineTime_1[j]][ch_no] += 1; + gui.y[events[i].fineTime_1[j]] += 1; + gui.content_2[events[i].fineTime_1[j]][events[i].tot_1[j]] += 1; + gui.content_3[div(64 - pixelmap[ch_no], 8).rem -4 + 1][div(64 - pixelmap[ch_no], 8).quot + 1] += 1; + // Increment counters --------------- + noRisingEdges += 1; + + } + } + } + + + ch_no = ch_no+1; + printf("Ch No = %i\n", ch_no); + for(int i = 0; i < events.size(); i++){ + for(int j=0; j < events[i].size(2); j++){ + + if(events[i].fineTime_2[j] <= 0){ + continue; + }else{ + + count_ch2++; + noRisingEdges += 1; + + HitMapArray[div(64 - pixelmap[ch_no], 8).rem -4][div(64 - pixelmap[ch_no], 8).quot] += 1;//, (1.0/noTriggers)); + RisingEdgesArray[events[i].fineTime_2[j]][ch_no] += 1; + gui.y[events[i].fineTime_2[j]] += 1; + gui.content_2[events[i].fineTime_2[j]][events[i].tot_2[j]] += 1; + gui.content_3[div(64 - pixelmap[ch_no], 8).rem -4 + 1][div(64 - pixelmap[ch_no], 8).quot + 1] += 1; + } + + } + } + + printf("Payload sum = %i\n", payload_sum); + printf("Number of zero length payloads = %i\n", number_of_zero_len_payloads); + printf("\nDATA:\n"); + printf("Ch %.2i = %llu / %i = %.3f%% \n", ch_no-1, count_ch1, noTriggers_channel, 100. * (double)count_ch1/noTriggers_channel ); + ChannelOccupancy[ch_no-1] = (double)count_ch1 / noTriggers_channel; + printf("Ch %.2i = %llu / %i = %.3f%% \n", ch_no, count_ch2, noTriggers_channel, 100. * (double)count_ch2/noTriggers_channel ); + ChannelOccupancy[ch_no] = (double)count_ch2 / noTriggers_channel; + printf("No Triggers = %i\n", noTriggers_channel); +// printf("No Triggers 2 = %i\n", noTriggers_channel_2); + printf("Final Trig count = %lu\n", events[(events.size()-1)].id); +// printf("First Trig count = %lu\n", events[0].trigNo_1[0]); + // printf("Trig diff = %lu\n", events[(events.size()-1)].id - events[0].trigNo_1[0]); + printf("\n"); + number_of_zero_len_payloads = 0; + All_channels.push_back(channel); + printf("\n---------------------------------------------- \n\n"); + } + + gui.Init(); + for(int kl=0;kl<16; kl++) printf("%s\n", all_filenames[kl]); + gui.TDCLoop(file_loc, all_filenames, argv[1]); + gui.Close(); + + printf("Total number of rising edges = %i\n", noRisingEdges); + + int hm_sum = 0; + for(int k = 0; k < 8; k++){ + for(int m = 0; m < 8; m++){ + hm_sum += HitMapArray[k][m]; + } + } + + printf("\n\n\n"); + printf("GRAPHS =======================================\n\n"); + + + printf("Rough Hitmap - Not mapped!\n"); + printf("--------------------------\n\n"); + for(int k = 0; k < 8; k++){ + printf(" | "); + for(int m = 0; m < 8; m++){ + //printf("%.2f | ", HitMapArray[k][m] / noRisingEdges); + if(HitMapArray[m][k] / hm_sum > 0.05){ + printf("X | "); + } else if(HitMapArray[m][k] / hm_sum > 0.01){ + printf("+ | "); + } else { + printf(" | "); + } + } + printf("\n"); + } + + + printf("\n\n\n"); + + printf("Rough Timing Alignment - Not optimised!\n"); + printf("---------------------------------------\n\n"); + + int re_sum[32] = {0}; + double means[32] = {0}; + for(int k = 0; k < 32; k++){ + for(int m = 0; m < 24 * TDCsignalLen; m++){ + re_sum[k] += RisingEdgesArray[m][k]; + } + means[k] = re_sum[k] / (24 * TDCsignalLen); + } + + for(int c = 0; c < 32; c++){ + printf("Ch %.2i = %.3f%% : ", c, 100. * ChannelOccupancy[c]); + for(int k = 0; k < 24; k++){ + int junk_count = 0; + for(int m = 0; m < TDCsignalLen; m++){ + junk_count += RisingEdgesArray[m + (k * TDCsignalLen)][c]; + } + + if((double) junk_count / re_sum[c] > 0.20){ + printf("X"); + } else if((double) junk_count / re_sum[c] > 0.10){ + printf("+"); + } else { + printf("-"); + } + } + printf("\n"); + } + + printf("\n\n\n"); + + + + +// for (int i = 0; i < All_channels.size(); i++){ +// int pl_sum_2 = 0; +// printf("Channel = %i\n",i); +// for (int j = 0; j < All_channels[i].size(); j++){ +// printf("%.3i ", All_channels[i][j]); +// pl_sum_2 += All_channels[i][j]; +// if((j+1)%8 == 0) printf("\n"); +// } +// printf("payload sum = %i\n", pl_sum_2); +// printf("\n\n"); +// } + + + + + return 1; + +} + +#endif + + +/* +int main(int argc,char **argv){ + + HistoGUI gui; + + std::vector x; + std::vector y; + + for(int i = 0; i < 25; i++){ + x.push_back( i ); + y.push_back( pow(i, 0.5) ); + } + + gui.SetData(x,y); + gui.Init(); + gui.Loop(); + gui.Close(); + + return 1; + +} + +*/ + + + + + + + + + diff --git a/QDK.h b/QDK.h index b915169..e6b6820 100644 --- a/QDK.h +++ b/QDK.h @@ -1,17 +1,9 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "global.h" using namespace std; -int QKN(double Energy, double radius, double distance, double thickness){ +double QK2(double Energy, double radius, double distance, double thickness){ double Qkn = 0.; @@ -30,6 +22,9 @@ int QKN(double Energy, double radius, double distance, double thickness){ double Tau = exp(TT); +// printf("TLN = %lf\n",TT); +// printf("Tau = %lf\n",Tau); + //calulating attenuation angles @@ -45,19 +40,24 @@ int QKN(double Energy, double radius, double distance, double thickness){ 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,sum2,sum3 = 0.; double sum4,sum5,sum6 = 0.; - double cosb,sinb,secb,c2,c4,fac1,ex1 = 0.; + double cosb,sinb,secb,cscb,c2,c4,fac1,fac2,ex1,ex2 = 0.; double term1,term2,term3 = 0.; double term4,term5,term6 = 0.; int J=0; - int loop_length = 500; + 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); @@ -65,7 +65,18 @@ int QKN(double Energy, double radius, double distance, double thickness){ }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; @@ -74,6 +85,7 @@ int QKN(double Energy, double radius, double distance, double thickness){ 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; @@ -88,14 +100,13 @@ int QKN(double Energy, double radius, double distance, double thickness){ double ans1 = sum1/3; double ans2 = sum2/3; double ans3 = sum3/3; - - double J2,B,beta2,cosb2,sinb2,secb2,cscb2,c22,c44,fac2,ex2 = 0.; + 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); @@ -103,29 +114,30 @@ int QKN(double Energy, double radius, double distance, double thickness){ }else {B=4.;} beta2 = LB+i+delx2; }else{B=.1;beta2 = LB+i+delx2;} - - cosb2 = cos(beta2); - sinb2 = sin(beta2); - secb2 = 1.0/cosb2; - cscb2 = 1.0/sinb2; - c22 = pow(cosb2,2); - c44 = pow(cosb2,4); - fac2 = -1 *Tau *(radius*cscb2 -distance*secb2); +*/ + 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); -/* - printf("cosb2---%lf\n",cosb2); - printf("sinb2---%lf\n",sinb2); - printf("secb2---%lf\n",secb2); - printf("cscb2---%lf\n",cscb2); - printf("c22---%lf\n",c22); - printf("c44---%lf\n",c44); - printf("fac2---%lf\n",fac2); - printf("ex2---%lf\n",ex2); - printf("-------------------------------------\n"); -*/ - term4 = 0.5*(3*c22-1)*(1-ex2)*sinb2*B*delx2; - term5 = 0.125*B*(35*c44-30*c22+3)*(1-ex2)*sinb2*delx2; - term6 = B*(1-ex2)*sinb2*delx2; + + 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; @@ -137,12 +149,12 @@ int QKN(double Energy, double radius, double distance, double thickness){ double ans6=sum6/3; /* - printf("ans1:%lf\n",ans1); - printf("ans2:%lf\n",ans2); - printf("ans3:%lf\n",ans3); - printf("ans4:%lf\n",ans4); - printf("ans5:%lf\n",ans5); - printf("ans6:%lf\n",ans6); + 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); @@ -164,7 +176,189 @@ int QKN(double Energy, double radius, double distance, double thickness){ fileo << "QD4 = " << QD4 << "\n"; fileo.close(); - return 1; + + return QD2; +} + + + + + +double QK4(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,sum2,sum3 = 0.; + double sum4,sum5,sum6 = 0.; + + double cosb,sinb,secb,cscb,c2,c4,fac1,fac2,ex1,ex2 = 0.; + double term1,term2,term3 = 0.; + double term4,term5,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; + + } + + + double ans1 = sum1/3; + double ans2 = sum2/3; + double ans3 = sum3/3; + + 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 QD4; } @@ -184,21 +378,6 @@ int QKN(double Energy, double radius, double distance, double thickness){ - - - - - - - - - - - - - - - diff --git a/README.md b/README.md index 634e1d9..dee9a28 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,15 @@ # Angular_Correlations # This program is under construction in order to replace a fortran 77 angular distribution/ mixing ratio determination calculator used a FSU's Nuclear Accelerator laboratory. + +# Reading in detector width, radius, distance from the radiative source, and the gamma-ray energy of interest, the program calculates the QDk geometric attenuation coefficients. + +# The program also asks for the magnetic substate probability distribution - or sigma - and the j1 and j2 values of the states from the gamma-ray decay from state 2 to 1. Using this information the coupling of states via clebsh - gordon coefficients in which the quantum mechanical restrictions given the angular momentum of the system can be used to understand the states decay probabilities following the reaction of interest. + +# This program then reads in a [].csv, the format must match the test file in the directory. This file is a file of intensities as a function of angles or different detectors for a given decay energy. The data is then used for the legendre polinomial f(x) = A0 P0(x) + A2 P2(x) + A4 P4(x), gaussian elimination is used to return A0, A2, A4. These values are used to weight the chi-squared minimization comparison from theory to experiment. The minimization is key to locating the electromagnetic mixing ratio, as a function of angular momentum and coupling quantum mechanics. + +# Finally reading from Rose and Brink, the program reads the first 2 Racah coefficient sets, 6 total, depending on j2 and j1. The program can be modified to use more than 3 angles, and for higher multipoles of the angular distribution. Using these coefficients, theoretical A2 and A4 values can be compared, illuminating the mixing ratio. + +# This program has its own graphical interface. Using X11 graphics, there is a functional terminal based menu that will promt and execute the 2 plots of interest. Option 1 displays the chi-squared minimization. Option 2 displays the Angular distribution overlayed with a legendre polinomial fit and respective error bars. Option 3 closes the program. + +# The GUI can zoom in by dragging and letting go, you can draw with left click, with right click the coordinates of the point is printed to the terminal. With space-bar you can unzoom the veiw. diff --git a/Racah.h b/Racah.h index 2f7c85d..11403e8 100644 --- a/Racah.h +++ b/Racah.h @@ -1,13 +1,4 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "global.h" //To compile : g++ AD.cxx -o {Input Executable Name} -lX11 //#include "Coeff.h" @@ -300,4 +291,3 @@ double RACAH(double A[6]){ } - diff --git a/RkTable.csv b/RkTable.csv new file mode 100644 index 0000000..9ed7dd9 --- /dev/null +++ b/RkTable.csv @@ -0,0 +1,115 @@ +1.,0.,0.7071,0.,0.,0.,0.,0. +1.,1.,-.3536,1.0607,-0.3536,0.,0.,0. +1.,2.,0.0707,-.4743,0.3536,0.,0.,0. +1.,3.,-.101,-.378,.5303,0.,0.,0. +1.,4.,-.1768,-.3062,.6010,0.,0.,0. +2.,0.,-.5976,0.,0.,-1.069,0.,0. +2.,1.,.4183,.9354,-.2988,0.,0.,.7127 +2.,2.,-.4183,.6124,.1281,0.,0.,-.3054 +2.,3.,.1195,-.6547,.3415,0.,0.,0.0764 +2.,4.,-.1707,-.5051,.4482,-.0085,.0627,-.0297 +2.,5.,-.2988,-.4009,.4618,.004,.0815,-.1093 +3.,0.,-.866,0.,0.,.2132,0.,0. +3.,1.,-.4949,.4629,-.6495,-.4467,-1.0446,.0355 +3.,2.,.3464,.9487,-.1237,0.,0.,.6701 +3.,3.,-.4330,.433,.2268,0.,0.,-.4467 +3.,4.,.1443,-.3093,.3093,0.,0.,.1489 +3.,5.,-.2062,-.5455,.3608,-.0203,.1343,-.0549 +3.,6.,-.3608,-.427,.3346,.0097,.172,-.1946 +4.,1.,-.7835,.2714,-.8234,.1453.,-.7549,.3017 +4.,2.,-.4477,.5292,-.4701,-.3044,-.9004,-.0484 +4.,3.,.3134,.9402,-.0448,0.,0.,.6088 +4.,4.,-.4387,.3354,.2646,0.,0.,-.4981 +4.,5.,.1595,-.7568,.2849,0.,0.,.1937 +4.,6.,-.2279,-.5641,.2991,-.0298,.1844,-.0687 +4.,7.,-.3989,-.4369,.2466,.0142,.2337,-.2363 +5.,2.,-.736,.3291,-.6825,.1159,-.7774,.0802 +5.,3.,-.4206,.5563,-.368,-.2428,-.803,-.0773 +5.,4.,.2944,.9309,0.,0.,0.,.5666 +5.,5.,-.4416,.2739,.2831,0.,0.,-.523 +5.,6.,.1698,-.7783,.2669,0.,0.,.2241 +5.,7.,-.2426,-.5742,.2548,-.0374,.221,-.0773 +5.,8.,-.4246,-.4413,.1837,.0178,.2779,-.2592 +6.,3.,-.7051,.3575,-.5915,.0997-.7581,-.0296 +6.,4.,-.4029,.5698,-.3022,-.2088,-.7383,-.0902 +6.,5.,.282,.9231,.0288,0.,0.,.537 +6.,6.,-.4432,.2315,.2935,0.,0.,-.537 +6.,7.,.1773,-.7928,.2533,0.,0.,.2461 +6.,8.,-.2533,-.5803,.2216,-.0434,.2488,-.0829 +6.,9.,-.4432,-.4432,.137,.0207,.3109,-.2727 +7.,4.,-.6834,.3743,-.5281,.0894,-.7349,-.0929 +7.,5.,-.3905,.5776,-.2563,-.1874,-.6929,-.0969 +7.,6.,.2734,.9169,.0488,0.,0.,.5154 +7.,7.,-.4442,.2004,.3001,0.,0.,-.5457 +7.,8.,.1829,-.8033,.2426,0.,0.,.2627 +7.,9.,-.2613,-.5843,.1960,-.0484,.2706,-.0869 +7.,10.,-.4573,-.444,.101,.0231,.3364,-.281 +8.,5.,-.6673,.3853,-.4813,.0824,-.7138,-.1332 +8.,6.,-.3813,.5825,-.2224,-.1727,-.6595,-.1007 +8.,7.,.2669,.9117,.0636,0.,0.,.4989 +8.,8.,-.4449,.1768,.3044,0.,0.,-.5514 +8.,9.,.1873,-.8111,.2341,0.,0.,.2757 +8.,10.,-.2676,-.587,.1756,-.0525,.288,-.0898 +9.,6.,-.6549,.393,-.4453,.0773,-.6959,-.1606 +9.,7.,-.3742,.5858,-.1965,-.162,-.635,-.1031 +9.,8.,.262,.9075,.0748,0.,0.,.486 +9.,9.,-.4453,.1581,.3074,0.,0.,-.5555 +9.,10.,.1909,-.8172,.2272,0.,0.,.2861 +10.,7.,-.6451,.3986,-.4169,.0734,-.6808,-.1803 +10.,8.,-.3686,.5881,-.1759,-.1539,-.6138,-.1046 +10.,9.,.258,.9039,.0838,0.,0.,.4757 +10.,10.,-.4457,.143,.3096,0.,0.,-.5584 +1.5,.5,.5,.866,-.5,0.,0.,0. +1.5,1.5,-.4,.7746,0.,0.,0.,0. +1.5,2.5,.1,-.5916,.3571,0.,0.,0. +1.5,3.5,-.1429,-.4629,.5,0.,0.,0. +1.5,4.5,-.25,-.3708,.5409,0.,0.,0. +2.5,0.5,-.5345,.3780,-.8018,-.6172,-1.0911,.1543 +2.5,1.5,.3742,.9487,-.1909,0.,0.,.7054 +2.5,2.5,-.4276,.5071,.1909,0.,0.,-.3968 +2.5,3.5,.1336,-.6944,.3245,0.,0.,.1176 +2.5,4.5,-.1909,-.5297,.4009,-.0147,.1019,-.0444 +2.5,5.5,-.3341,-.4173,.3924,.007,.1314,-.1602 +3.5,.5,-.8183,.2113,-.9274,.1709,-.6621,.5128 +3.5,1.5,-.4676,.5051,-.5455,-.3582,-.9671,-.019 +3.5,2.5,.3273,.9449,-.0779,0.,0.,.6367 +3.5,3.5,-.4364,.378,.2494,0.,0.,-.4775 +3.5,4.5,.1528,-.7416,.2962,0.,0.,.1737 +3.5,5.5,-.2182,-.5563,.3273,-.0253,.1614,-.0627 +3.5,6.5,-.3819,-.433,.2867,.0121,.2056,-.2188 +4.5,1.5,-.7569,.3062,-.7444,.1281,-.7774,.1693 +4.5,2.5,-.4325,.5455,-.4129,-.2684,-.8465,-.066 +4.5,3.5,.3028,.9354,-.0197,0.,0.,.5857 +4.5,4.5,-.4404,.3015,.2752,0.,0.,-.5125 +4.5,5.5,.1651,-.7687,.2752,0.,0.,.2102 +4.5,6.5,-.2359,-.5698,.2752,-.0338,.2040,-.0735 +4.5,7.5,-.4129,-.4395,.2127,.0161,.2575,-.2493 +5.5,2.5,-.7191,.3454,-.6326,.1067,-.7692,.0171 +5.5,3.5,-.4109,.5614,-.3319,-.2237,-.7676,-.0849 +5.5,4.5,.2876,.9269,.0158,0.,0.,.5505 +5.5,5.5,-.4425,.2509,.289,0.,0.,-.5309 +5.5,6.5,.1738,-.7862,.2596,0.,0.,.2359 +5.5,7.5,-.2483,-.5776,.2371,-.0406,.2358,-.0804 +5.5,8.5,-.4346,-.4424,.1588,.0194,.2956,-.2668 +6.5,3.5,-.6934,.3669,-.5572,.094,-.7464,-.0651 +6.5,4.5,-.3962,.5742,-.2774,-.197,-.7138,-.094 +6.5,5.5,.2774,.9199,.0396,0.,0.,.5254 +6.5,6.5,-.4438,.2148,.2972,0.,0.,-.5418 +6.5,7.5,.1803,-.7984,.2476,0.,0.,.255 +6.5,8.5,-.2575,-.5825,.208,-.046,.2603,-.0851 +6.5,9.5,-.4507,-.4437,.1179,.022,.3245,-.2773 +7.5,4.5,-.6748,.3803,-.503,.0856,-.724,-.1151 +7.5,5.5,-.3856,.5803,-.2382,-.1794,-.6751,-.0991 +7.5,6.5,.2699,.9142,.0567,0.,0.,.5066 +7.5,7.5,-.4446,.1879,.3024,0.,0.,-.5488 +7.5,8.5,.1852,-.8074,.2382,0.,0.,.2696 +7.5,9.5,-.2646,-.5858,.1852,-.0506,.2797,-.0885 +8.5,5.5,.6607,.3894,-.4622,.0797,-.7045,-.1481 +8.5,6.5,.3776,.5843,-.2087,-.1670,-.6459,-.1021 +8.5,7.5,.2643,.9095,.0696,0.,0.,.4921 +8.5,8.5,-.4451,.1699,.3060,0.,0.,-.5536 +8.5,9.5,.1892,-.8143,.2305,0.,0.,.2812 +9.5,6.5,-.6497,.396,-.4303,.0753,-.688,-.1712 +9.5,7.5,-.3713,.587,-.1856,-.1577,-.6233,-.1039 +9.5,8.5,.2599,.9056,0.0796,0.,0.,.4806 +9.5,9.5,-.4455,.1502,.3086,0.,0.,-.557 diff --git a/Cleb.cxx b/Test_Scripts/Cleb.cxx similarity index 81% rename from Cleb.cxx rename to Test_Scripts/Cleb.cxx index 4c3796c..6e850c7 100644 --- a/Cleb.cxx +++ b/Test_Scripts/Cleb.cxx @@ -1,13 +1,5 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include + +#include "global.h" //To compile : g++ AD.cxx -o {Input Executable Name} -lX11 #include "Coeff.h" @@ -179,11 +171,13 @@ double CG2(double A[6]){ for(int i = 0; i<5; i++){ IXI = nx[i]; termlg= termlg -lgamma(IXI); - SS = SS + S1*exp(termlg); - S1 =-S1; - Cleb=sgnfac*SS; } - } + SS = SS + S1*exp(termlg); + S1 =-S1; + } + Cleb=sgnfac*SS; + + }else{ // if (IE >=0) _>>> ID = -ID; @@ -204,12 +198,12 @@ double CG2(double A[6]){ //lgamma(arg) returns the log of the factorial of arg if arg is a natural number. FC2 = IC + 1; double sqfclg = log(FC2) - lgamma(IT+1); - +/* printf("IT = %lf\n",IT); printf("sqfclg = %lf\n",sqfclg); printf("FC2 = %lf\n",FC2); - printf("lgamma(IT) = %lf\n",lgamma(IT+1)); + printf("lgamma(IT+1) = %lf\n",lgamma(IT+1)); printf("log(fC2) = %lf\n",log(FC2)); printf("IA = %lf\n",IA); @@ -218,15 +212,15 @@ double CG2(double A[6]){ printf("ID = %lf\n",ID); printf("IE = %lf\n",IE); printf("IF = %lf\n",IF); - +*/ double IXI; for(int i = 0; i<9;i++){ IXI = ix[i]; sqfclg = sqfclg + lgamma(IXI); - printf("sqfclg = %lf\n",sqfclg); + // printf("sqfclg = %lf\n",sqfclg); } sqfclg = 0.5 * sqfclg; - +/* printf("Here place 2\n"); printf("ix[0] = %lf\n",ix[0]); @@ -238,7 +232,7 @@ double CG2(double A[6]){ printf("ix[6] = %lf\n",ix[6]); printf("ix[7] = %lf\n",ix[7]); printf("ix[8] = %lf\n",ix[8]); - +*/ printf("sqdclg = %lf\n",sqfclg); double xarr[3] = {ix[0],ix[4],ix[5]}; @@ -282,6 +276,13 @@ double CG2(double A[6]){ if(NZMAX > NZMIN){return 0;} double SS=0.; double S1= pow((-1.),(NZMIN-1)); + + printf("NZMAX = %lf\n",NZMAX); + printf("NZMIN = %lf\n",NZMIN); + printf("NZ2 = %lf\n",NZ2); + printf("NZ3 = %lf\n",NZ3); + printf("S1 = %lf\n",S1); + for(int NZ = NZMIN;NZ<=NZMAX;NZ++){ double NZM1 = NZ -1; @@ -291,17 +292,25 @@ double CG2(double A[6]){ nx[3] = NZ-NZ2; nx[4] = NZ-NZ3; double termlg = sqfclg - lgamma(NZ); + + printf("nx[0] = %lf\n",nx[0]); + printf("nx[1] = %lf\n",nx[1]); + printf("nx[2] = %lf\n",nx[2]); + printf("nx[3] = %lf\n",nx[3]); + printf("nx[4] = %lf\n",nx[4]); + + printf("termlg = %lf\n",termlg); + for(int i = 0; i<5; i++){ IXI = nx[i]; termlg= termlg -lgamma(IXI); - SS = SS + S1*exp(termlg); - S1 =-S1; - Cleb=sgnfac*SS; - } - } + } + SS = SS + S1*exp(termlg); + S1 =-S1; + + } + Cleb=sgnfac*SS; - - } // if (IB != 0) _>>>> }else{ @@ -441,19 +450,24 @@ double CG2(double A[6]){ for(int i = 0; i<5; i++){ IXI = nx[i]; termlg= termlg -lgamma(IXI); - SS = SS + S1*exp(termlg); - S1 =-S1; - Cleb=sgnfac*SS; } - } + SS = SS + S1*exp(termlg); + S1 =-S1; + + } + Cleb=sgnfac*SS; + + }else{ // if (IE >=0) _>>> + ID = -ID; IE = -IE; IF = -IF; if((int)(IA+IB-IC)/2 %2 != 0) sgnfac = (-1.)*sgnfac; IT = K0/2 + 1; + FC2 = IT +1; ix[0] = IT-IC; ix[1] = IT-IB; ix[2] = IT-IA; @@ -537,13 +551,11 @@ double CG2(double A[6]){ for(int i = 0; i<5; i++){ IXI = nx[i]; termlg= termlg -lgamma(IXI); - SS = SS + S1*exp(termlg); - S1 =-S1; - Cleb=sgnfac*SS; } - } - - + SS = SS + S1*exp(termlg); + S1 =-S1; + } + Cleb=sgnfac*SS; } // if (IB != 0) _>>>> @@ -576,39 +588,72 @@ double CG2(double A[6]){ } -double CG(double J, double M, double j1, double m1, double j2, double m2){ +//double CG(double J, double M, double j1, double m1, double j2, double m2){ +double CG(double j1, double j2, double J, double m1, double m2, double M){ //recall that j1,m1 + j2,m2 = J,M + printf("-----------------\n"); if(M != m1 + m2) return 0; double Jmin = abs(j1 - j2); double Jmax = j1+j2; + printf(" Jmin = %lf\n", Jmin); + printf(" Jmax = %lf\n", Jmax); + if(J < Jmin || Jmax < J) return 0; - double a0 = (2*J+1.0)*factorial(J+j1-j2) * factorial(J-j1+j2) * factorial(j1+j2-J)/factorial(J+j1+j2+1.0); +// double a0 = (2*J+1.0)*factorial(J+j1-j2) * factorial(J-j1+j2) * factorial(j1+j2-J)/factorial(J+j1+j2+1.0); + + + double a0 = (2*J+1.0)*tgamma(J+j1-j2+1) * tgamma(J-j1+j2+1) * tgamma(j1+j2-J+1)/tgamma(J+j1+j2+1.0 +1); double A0 = sqrt(a0); - double a = factorial(J+M) *factorial(J-M); - double a1= factorial(j1+m1) *factorial(j1-m1); - double a2= factorial(j2+m2) *factorial(j2-m2); +// double a = factorial(J+M) *factorial(J-M); +// double a1= factorial(j1+m1) *factorial(j1-m1); +// double a2= factorial(j2+m2) *factorial(j2-m2); + + double a = tgamma(J+M+1) *tgamma(J-M+1); + double a1= tgamma(j1+m1+1) *tgamma(j1-m1+1); + double a2= tgamma(j2+m2+1) *tgamma(j2-m2+1); + double A = sqrt( a * a1 * a2); + printf(" a0 = %lf\n", a0); + printf(" A0 = %lf\n", A0); + printf(" a = %lf\n", a); + printf(" a1 = %lf\n", a1); + printf(" a2 = %lf\n", a2); + printf(" A = %lf\n", A); + int pmax = min( min(j1+j2-J,j1-m1),j2 + m2); + printf("pmax = %d\n",pmax); + double cg = 0.; for( int p =0; p<=pmax;p++){ - double p1 = factorial(j1+j2-J-p); - double p2 = factorial(j1-m1-p); - double p3 = factorial(j2+m2-p); - double p4 = factorial(J -j2 +m1 +p); - double p5 = factorial(J -j1 -m2 +p); - double t = pow(-1,p)/(factorial(p) * p1 * p2 * p3 * p4 * p5); - +// double p1 = factorial(j1+j2-J-p); +// double p2 = factorial(j1-m1-p); +// double p3 = factorial(j2+m2-p); +// double p4 = factorial(J -j2 +m1 +p); +// double p5 = factorial(J -j1 -m2 +p); +// double t = pow(-1,p)/(factorial(p) * p1 * p2 * p3 * p4 * p5); + + + double p1 = tgamma(j1+j2-J-p+1); + double p2 = tgamma(j1-m1-p+1); + double p3 = tgamma(j2+m2-p+1); + double p4 = tgamma(J -j2 +m1 +p+1); + double p5 = tgamma(J -j1 -m2 +p+1); + double t = pow(-1,p)/(tgamma(p+1) * p1 * p2 * p3 * p4 * p5); + + printf("t = %lf\n",t); cg += t; + printf("cg = %lf\n",cg); } + printf("-----------------\n"); return A0*A*cg; } @@ -722,15 +767,44 @@ double SixJsym(double j1, double j2, double j3, double j4, double j5, double j6) int main(int argc, char ** argv){ - - double A[6] = {1,1,2,0,0,0}; + //if mod (2*J1,2) = 1 do this +// double A0[6] = {1,1,2,0,0,0}; +// double A1[6] = {1,1,4,0,0,0}; + + //if mod(2*j1,20 = 0 do this + + double j1 = 1.; + double j2 = 2.; - double cg = CG2(A); - double cgr = CG(1,1,2,0,0,0); + double A0[6] = {j1,j1,2,.5,-.5,0}; + double A1[6] = {j1,j1,2,-.5,.5,0}; + double A2[6] = {j1,j1,4,-.5,.5,0}; + double A3[6] = {j1,j1,4,.5,-.5,0}; + +/* + double cg0 = CG2(A0); + double cg1 = CG2(A1); + double cg2 = CG2(A2); + double cg3 = CG2(A3); +*/ + double cgr0= CG(j1,j1,2,.5,-.5,0); + double cgr1= CG(j1,j1,2,-.5,.5,0); + double cgr2= CG(j1,j1,4,.5,-.5,0); + double cgr3= CG(j1,j1,4,-.5,.5,0); +/* printf("------\n"); - printf("CG = %lf\n",cg); - printf("CGR = %lf\n",cgr); + printf("CG0 = %lf\n",cg0); + printf("CG1 = %lf\n",cg1); + printf("CG2 = %lf\n",cg2); + printf("CG3 = %lf\n",cg3); +*/ + printf("------\n"); + printf("CGR0 = %lf\n",cgr0); + printf("CGR1 = %lf\n",cgr1); + printf("CGR2 = %lf\n",cgr2); + printf("CGR3 = %lf\n",cgr3); + return 0; diff --git a/Test_Scripts/QDK.cxx b/Test_Scripts/QDK.cxx new file mode 100644 index 0000000..c6b713e --- /dev/null +++ b/Test_Scripts/QDK.cxx @@ -0,0 +1,230 @@ +#include "global.h" + +using namespace std; + + +int QKN(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,sum2,sum3 = 0.; + double sum4,sum5,sum6 = 0.; + + double cosb,sinb,secb,cscb,c2,c4,fac1,fac2,ex1,ex2 = 0.; + double term1,term2,term3 = 0.; + double term4,term5,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; + + } + + + double ans1 = sum1/3; + double ans2 = sum2/3; + double ans3 = sum3/3; + + 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 1; +} + + +int main(int argc, char** argv){ + + double E = 1062.; + double dr = 3.; + double dd = 4.; + double dt = 5.; + + QKN(E,dr,dd,dt); + + + return 0; +} + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Test_Scripts/Racah.cxx b/Test_Scripts/Racah.cxx new file mode 100644 index 0000000..52f5dd1 --- /dev/null +++ b/Test_Scripts/Racah.cxx @@ -0,0 +1,481 @@ +#include "global.h" +//To compile : g++ AD.cxx -o {Input Executable Name} -lX11 + +//#include "Coeff.h" + +using namespace std; + + + +double RACAH(double A[6]){ + + double Racah = 0.; + + double ii[6],ix[6],kk[8]; + + double nn[4], nx[12],ny[8]; + + double K1,K2,K3,K4,K5,K6,K7,K8 = 0.; + + double NZMIN, NZMAX = .0; + + for(int i = 0;i<6;i++){ + ii[i] = 2.0*A[i]; + } + + double IA = ii[0]; + double IB = ii[1]; + double IC = ii[2]; + double ID = ii[3]; + double IE = ii[4]; + double IF = ii[5]; + + K1 = IA + IB -IE; + K3 = IC + ID -IE; + K5 = IA + IC -IF; + K7 = IB + ID -IF; + K2 = IE - abs(IA-IB); + K4 = IE - abs(IC-ID); + K6 = IF - abs(IA-IC); + K8 = IF - abs(IB-ID); + + + double karr[8] = {K1,K2,K3,K4,K5,K6,K7,K8}; + double kmn=karr[0]; + double kmx=karr[0]; + for(int i=0;i<8;i++) + { + if(kmn>karr[i]) + { + kmn=karr[i]; + } + else if(kmx 0){ + break; + } + } + + double sgn1 = 1.; + double KQ,IP,IQ,IT = 0.; + if(Kmin == 1){ + KQ = 7 - Kmin; + + IA = ix[(int)KQ-1]; + ID = ix[Kmin+4-1]; + IE = ix[Kmin-1]; + IF = ix[(int)KQ-4-1]; + + if(((int)(IA+ID-IE-IF)/2 %2) != 0) {sgn1 = -1.;} + + }else if(Kmin == 2){ + + KQ = 7 - Kmin; + + IA = ix[(int)KQ-1]; + ID = ix[Kmin+4-1]; + IE = ix[Kmin-1]; + IF = ix[(int)KQ-4-1]; + + if(((int)(IA+ID-IE-IF)/2 %2) != 0) {sgn1 = -1.;} + + }else if(Kmin == 3){ + KQ = 9 - Kmin; + + IB = ix[(int)KQ-1]; + IC = ix[Kmin+2-1]; + IE = ix[Kmin-1]; + IF = ix[(int)KQ -2-1]; + if(((int)(IB+IC-IE-IF)/2)%2 !=0) {sgn1 = -1.;} + }else if(Kmin == 4){ + + KQ = 9 - Kmin; + + IB = ix[(int)KQ-1]; + IC = ix[Kmin+2-1]; + IE = ix[Kmin-1]; + IF = ix[(int)KQ -2-1]; + if(((int)(IB+IC-IE-IF)/2)%2 !=0) {sgn1 = -1.;} + + }else if(Kmin == 5){ + + IE = ix[4]; + IF = ix[5]; + IB = ix[3]; + IC = ix[2]; + + }else if(Kmin == 6){ + + } + + IP = IA-IB; + IQ = IC-ID; + + if(abs(IP) >= abs(IQ)){ + if(IP >= 0){ + IQ = IC-ID; + double sgn2 = 1.0; + if(((int)(IB+ID-IF)/2 % 2) != 0) { sgn2 = -1.;} + if( IE > 0 ) { + + + printf("Here 177\n"); + printf("ia = %lf\n", IA); + printf("ib = %lf\n", IB); + printf("ic = %lf\n", IC); + printf("id = %lf\n", ID); + printf("ie = %lf\n", IE); + printf("if = %lf\n", IF); + printf("it = %lf\n", IT); + printf("ip = %lf\n", IP); + printf("iq = %lf\n", IQ); +// printf("bi = %lf\n", BI); + + nn[0] = (IA+IB+IE)/2 +1; + nn[1] = (IC+ID+IE)/2 +1; + nn[2] = (IA+IC+IF)/2 +1; + nn[3] = (IB+ID+IF)/2 +1; + + nx[0] = nn[0] -IE; + nx[1] = nn[0] -IB; + nx[2] = nn[0] -IA; + nx[3] = nn[1] -IE; + nx[4] = nn[1] -ID; + nx[5] = nn[1] -IC; + nx[6] = nn[2] -IF; + nx[7] = nn[2] -IC; + nx[8] = nn[2] -IA; + nx[9] = nn[3] -IF; + nx[10] = nn[3] -ID; + nx[11] = nn[3] -IB; + + IP = (IA+IB+IC+ID+4)/2; + IQ = (IE+IF-IA-ID)/2; + IT = (IE+IF-IB-IC)/2; + + double a = 1; + double b = -IQ+1; + if(a > b) { + + NZMIN = a; + }else{ + NZMIN = b; + + } + + double kzar[3] = {nx[1],nx[4],nx[10]}; + double kzmn=kzar[0]; + double kzmx=kzar[0]; + for(int i=0;i<3;i++) + { + if(kzmn>kzar[i]) + { + kzmn=kzar[i]; + } + else if(kzmx= 0){ + IQ = IC - ID; + }else{ + + IT = IB; + IB = IA; + IA = IT; + IT = ID; + ID = IC; + IC = IT; + IP = -IP; + + } + } + + double sgn2 = 1.; + + if((int)(IB+ID-IF)/2 %2 != 0) sgn2 = -1.; + if(IE >= 0){ + double BI = IB; + double DI = ID; + Racah = sgn1 * sgn2 / sqrt((BI+1.)*(DI+1.)); + }else{ + + + printf("Here 340\n"); + printf("IA = %lf\n", IA); + printf("IB = %lf\n", IB); + printf("IC = %lf\n", IC); + printf("ID = %lf\n", ID); + printf("IE = %lf\n", IE); + printf("IF = %lf\n", IF); + printf("IT = %lf\n", IT); + printf("IP = %lf\n", IP); + printf("IQ = %lf\n", IQ); +// printf("BI = %lf\n", BI); + + + nn[0] = (IA+IB+IE)/2 +1; + nn[1] = (IC+ID+IE)/2 +1; + nn[2] = (IA+IC+IF)/2 +1; + nn[3] = (IB+ID+IF)/2 +1; + + nx[0] = nn[0] -IE; + nx[1] = nn[0] -IB; + nx[2] = nn[0] -IA; + nx[3] = nn[1] -IE; + nx[4] = nn[1] -ID; + nx[5] = nn[1] -IC; + nx[6] = nn[2] -IF; + nx[7] = nn[2] -IC; + nx[8] = nn[2] -IA; + nx[9] = nn[3] -IF; + nx[10] = nn[3] -ID; + nx[11] = nn[3] -IB; + + IP = (IA+IB+IC+ID+4)/2; + IQ = (IE+IF-IA-ID)/2; + IT = (IE+IF-IB-IC)/2; + + double a = 1; + double b = -IQ+1; + if(a > b) { + + NZMIN = a; + }else{ + NZMIN = b; + + } + + double kzar[3] = {nx[1],nx[4],nx[10]}; + double kzmn=kzar[0]; + double kzmx=kzar[0]; + for(int i=0;i<3;i++) + { + if(kzmn>kzar[i]) + { + kzmn=kzar[i]; + } + else if(kzmx 0) { al0 = l0abs;} + if(l0abs <= 0){ al0 = 1.;} + + double al1 = al0 + 1; + + + double A[6] = {j1,j1,al0,al0,2,j2}; + double B[6] = {j1,j1,al0,al0,4,j2}; + + double RCHA = RACAH(A); + double RCHB = RACAH(B); + + printf("Racah A = %lf\n",RCHA); + printf("Racah B = %lf\n",RCHB); + + + return 0; +} diff --git a/Test_Scripts/Racah_Cleb.cxx b/Test_Scripts/Racah_Cleb.cxx new file mode 100644 index 0000000..62aa207 --- /dev/null +++ b/Test_Scripts/Racah_Cleb.cxx @@ -0,0 +1,182 @@ + +#include "global.h" +//To compile : g++ AD.cxx -o {Input Executable Name} -lX11 + +#include "Coeff.h" + +using namespace std; + + + +double CG(double j1, double j2, double J, double m1, double m2, double M){ + //recall that j1,m1 + j2,m2 = J,M + +// printf("-----------------\n"); + if(M != m1 + m2) return 0; + + double Jmin = abs(j1 - j2); + double Jmax = j1+j2; + + + if(J < Jmin || Jmax < J) return 0; + + + + double a0 = (2*J+1.0)*tgamma(J+j1-j2+1) * tgamma(J-j1+j2+1) * tgamma(j1+j2-J+1)/tgamma(J+j1+j2+1.0 +1); + double A0 = sqrt(a0); + + + double a = tgamma(J+M+1) *tgamma(J-M+1); + double a1= tgamma(j1+m1+1) *tgamma(j1-m1+1); + double a2= tgamma(j2+m2+1) *tgamma(j2-m2+1); + + double A = sqrt( a * a1 * a2); + + + int pmax = min( min(j1+j2-J,j1-m1),j2 + m2); + + + double cg = 0.; + + for( int p =0; p<=pmax;p++){ + + + + double p1 = tgamma(j1+j2-J-p+1); + double p2 = tgamma(j1-m1-p+1); + double p3 = tgamma(j2+m2-p+1); + double p4 = tgamma(J -j2 +m1 +p+1); + double p5 = tgamma(J -j1 -m2 +p+1); + double t = pow(-1,p)/(tgamma(p+1) * p1 * p2 * p3 * p4 * p5); + + cg += t; + } + return A0*A*cg; +} + + +double ThreeJsym(double j1, double m1, double j2, double m2, double j3, double m3){ + //[j1 j2 j3] = (-1)^(j1-j2-m3)/ sqrt(2*j3+1) * CG(j3, -m3, j1, m1, j2, m2) + //[m1 m2 m3] + + + //J,M,j1,m1,j2,m2) + //return pow(-1,j1 -j2 -m3)/sqrt(2*j3+1)*CG(j3,-m3,j1,m1,j2,m2); + double threej = pow(-1,j1 -j2 -m3)/sqrt(2*j3+1)*CG(j3,-m3,j1,m1,j2,m2); + + //double threej = pow(-1,j1 -j2 -m3)/sqrt(2*j3+1)*CG(j1,j2,j3,m1,m2,-m3); + printf("---------\n"); + printf("3J= %lf\n",threej); + + printf("---------\n"); + + return threej; +} + +double SixJsym(double j1, double j2, double j3, double j4, double j5, double j6){ + +//--------------------------------------------------------------------------------// + // The six j symbol describes the coupling between j1 j2 and j3 to J - j1. + // essentially a triangle of angular momentum rules between these. + // j1 = j1 + // j2 = j2 + // j3 = j1 + j2 + // j4 = j3 + // j5 = J = j1 + j2 + j3 + // j6 = j2 + j3 +// ----------------------------------------------------------------------------- // +// the following conditions check the triangle selection rules + + 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; +// now that they have been checked, we can go ahead and calculate sixJ. + + double sixj = 0.0; + float m1 = -j1; + float m2 = -j2; + float m3 = -j3; + float m4 = -j4; + float m5 = -j5; + float m6 = -j6; +printf("here\n"); + for(; m1 <= j1; m1 = m1 +1){ + for(; m2 <= j2; m2 = m2 +1){ + for(; m3 <= j3; m3 = m3 + 1){ + for(; m4 <= j4; m4 = m4 +1){ + for(; m5 <= j5; m5 = m5 + 1){ + for(; m6 <= j6; m6 = m6 +1){ + + double h = (j1 - m1) + (j2 - m2) + (j3 -m3) + (j4 - m4) + (j5 - m5) + (j6 - m6); + double b1 = ThreeJsym(j1, -m1, j2, -m2, j3, -m3); + double b2 = ThreeJsym(j1, m1, j5, -m5, j6, m6); + double b3 = ThreeJsym(j4, m4, j2, m2, j6, -m6); + double b4 = ThreeJsym(j4, -m4, j5, m5, j3, m3); + double b = b1 * b2 * b3 * b4; + + + printf("h = %lf\n", h); + printf("b1 = %lf\n", b1); + printf("b2 = %lf\n", b2); + printf("b3 = %lf\n", b3); + printf("b4 = %lf\n", b4); + printf("b = %lf\n", b); + + sixj += pow(-1,h)*b; + + } + } + } + } + } + } + + return sixj; +} + + +int main(int argc, char ** argv){ + + //if mod (2*J1,2) = 1 do this +// double A0[6] = {1,1,2,0,0,0}; +// double A1[6] = {1,1,4,0,0,0}; + + //if mod(2*j1,20 = 0 do this + + double j1 = 1.; + double j2 = 2.; + + + double A0[6] = {j1,j1,2,.5,-.5,0}; + double A1[6] = {j1,j1,2,-.5,.5,0}; + double A2[6] = {j1,j1,4,-.5,.5,0}; + double A3[6] = {j1,j1,4,.5,-.5,0}; + + double cgr0= CG(j1,j1,2,.5,-.5,0); + double cgr1= CG(j1,j1,2,-.5,.5,0); + double cgr2= CG(j1,j1,4,.5,-.5,0); + double cgr3= CG(j1,j1,4,-.5,.5,0); + + + printf("------\n"); +// printf("CGR0 = %lf\n",cgr0); +// printf("CGR1 = %lf\n",cgr1); +// printf("CGR2 = %lf\n",cgr2); +// printf("CGR3 = %lf\n",cgr3); + + double sixj = SixJsym(j1,j1,1.,1.,2.,2.); + double sixj2 = SixJsym(j1,j1,1.,2.,2.,2.); + double sixj3 = SixJsym(j1,j1,2.,2.,2.,2.); + + + + printf("6-J = %lf\n", sixj); + + + + return 0; +} + + + diff --git a/ad.txt b/ad.txt index de93b13..7d25f53 100644 --- a/ad.txt +++ b/ad.txt @@ -3,5 +3,5 @@ Distance = 4 [cm] Thickness = 5 [cm] Atten.C = 0.294297 [cm^-1] Gamma_E = 1062 [KeV] -QD2 = -0.499975 -QD4 = 0.374938 +QD2 = 0.801748 +QD4 = 0.447498 diff --git a/global.h b/global.h new file mode 100644 index 0000000..5c10f13 --- /dev/null +++ b/global.h @@ -0,0 +1,19 @@ +#ifndef __iamauniqueid_h__ +#define __iamauniqueid_h__ + +//define any global data type here. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#endif diff --git a/racah_table_reader_test.cxx b/racah_table_reader_test.cxx new file mode 100644 index 0000000..0cca1b0 --- /dev/null +++ b/racah_table_reader_test.cxx @@ -0,0 +1,112 @@ +#include "global.h" +#include "Functlib.h" + +using namespace std; + + +int main(int argc,char** argv){ + + + double j1 = 1.5; + double j2 = 3.5; + + + vector > content; + vector row; + string line, word; + + vector j1data; + vector j2data; + vector Rk20data; + vector Rk21data; + vector Rk22data; + vector Rk40data; + vector Rk41data; + vector Rk42data; + + + fstream file("RkTable.csv"); + if(file.is_open()){ + while(getline(file,line)){ + row.clear(); + + stringstream str(line); + + while(getline(str, word, ',')) + + row.push_back(word); + content.push_back(row); + j1data.push_back(row[0]); + j2data.push_back(row[1]); + Rk20data.push_back(row[2]); + Rk21data.push_back(row[3]); + Rk22data.push_back(row[4]); + Rk40data.push_back(row[5]); + Rk41data.push_back(row[6]); + Rk42data.push_back(row[7]); + + } + cout<< " --- Angular Data File 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 = " <