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

fixed AD fit, chi-squared wrong

Browse Source
This commit is contained in:
Peter DeRosa 2022-07-18 19:40:43 -04:00
parent 7efc5caa48
commit 434dd48296
3 changed files with 117 additions and 79 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
.AD.cxx.swp
View File

Binary file not shown.

112
AD.cxx
View File

@ -99,9 +99,9 @@ int main(int argc,char **argv){
}
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)
@ -226,30 +226,7 @@ int main(int argc,char **argv){
// Note, the y data must be scaled by sin(theta) of the given angle.
double a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12 = 0.;
/*
for(int i = 0; i< dangler.size(); i++){
//eq1
a1 += 1; //how many data points, not intensional but convientient.
a2 += (1.5 * pow(dangler[i],2) - .5);
a3 += (35./8. * pow(dangler[i],4) - 30./8. * pow(dangler[i],2) + 3./8. );
a4 += dydata[i];
//eq2
a5 += (1.5 * pow(dangler[i],2) - .5);
a6 += (1.5 * pow(dangler[i],2) - .5)*(1.5 * pow(dangler[i],2) - .5);
a7 += (1.5 * pow(dangler[i],2) - .5)*(35./8. * pow(dangler[i],4) - 30./8. * pow(dangler[i],2) + 3./8. );
a8 += (1.5 * pow(dangler[i],2) - .5)*dydata[i];
//eq4
a9 += (35./8. * pow(dangler[i],4) - 30./8. * pow(dangler[i],2) + 3./8. );
a10 += (1.5 * pow(dangler[i],2) - .5)*(35./8. * pow(dangler[i],4) - 30./8. * pow(dangler[i],2) + 3./8. );
a11 +=(35./8. * pow(dangler[i],4) - 30./8. * pow(dangler[i],2) + 3./8. )* (35./8. * pow(dangler[i],4) - 30./8. * pow(dangler[i],2) + 3./8. );
a12 += (35./8. * pow(dangler[i],4) - 30./8. * pow(dangler[i],2) + 3./8. )*dydata[i];
}
*/
for(int i = 0; i< dangler.size(); i++){
//eq1
@ -335,9 +312,9 @@ int main(int argc,char **argv){
}
cout<<"A0 = "<<residual[0]<<"\n";
cout<<"A2 = "<<residual[1]<<"\n";
cout<<"A4 = "<<residual[2]<<"\n";
// cout<<"A0 = "<<residual[0]<<"\n";
cout<<"a2 = "<<residual[1]/residual[0]<<"\n";
cout<<"a4 = "<<residual[2]/residual[0]<<"\n";
@ -504,9 +481,20 @@ int main(int argc,char **argv){
II = j1 - am11;
A[3] = am11;
A[4] = -am11;
if(isnan(CG2(A))){ // CG2(A) is NaN
cout << "Warning, a CG coefficient returned NaN, set to 0.\n";
cgg = 0.0;
}else{ //CG2(A) is not NaN.
cgg = CG2(A);
// cout << "cgg = " << " " << cgg << "\n";
}
// cout << "cgg = " << " " << cgg << "\n";
// cout << "am11 = " << " " << am11 << "\n";
// cout << "amsq1 = " << " " << amsq1 << "\n";
// cout << "x1 = " << " " << x1 << "\n";
@ -620,8 +608,11 @@ int main(int argc,char **argv){
double delta_min = -3.14159/2.;
double delta_max = 3.14159/2.;
double step = 0.001;
double THETA_min = 0.;
double THETA_max = 3.14159;
double step = 0.0001;
double Tangle = 0;
double delta = 0.;
double atan_delta =0.;
double A0E = residual[0];//NEED FrOM AF FIT
@ -629,6 +620,7 @@ int main(int argc,char **argv){
double A4E = residual[2];
double A2T,a2T = 0.;
double A4T,a4T = 0.;
double rd0T = 0.;
double rd2T = 0.;
double rd4T = 0.;
double X22,X24 = 0.;
@ -639,14 +631,53 @@ int main(int argc,char **argv){
int points = (delta_max - delta_min) / step ;
int t_points = (THETA_max - THETA_min) / step;
vector<double> chisqr;
vector<double> tdelta;
double YT0,YT2,YT4,YT,YE = 0.;
double Y_err = 10.;
//delta loop;
for(int i = 0; i< points; i++){
atan_delta = i*step + delta_min;
delta = tan(atan_delta);
rd2T = (rk01 + 2*delta*rk11 + pow(delta,2)*rk21)/(1+pow(delta,2));
delta = i*step + delta_min;
// delta = tan(atan_delta);
//now sum over all theta;
for(int j = 0; j < dangler.size(); j++){
Tangle = dangler[j];
//calculate Y_intensity_theory and and experiment and sum over theta
// rd0T = (rk01 + 2*delta*rk11 + pow(delta,2)*rk21)/(1+pow(delta,2));
rd2T = (rk01 + 2*delta*rk11 + pow(delta,2)*rk21)/(1+pow(delta,2));
rd4T = (rk02 + 2*delta*rk12 + pow(delta,2)*rk22)/(1+pow(delta,2));
YT0 = 1; // P0(cos(theta)) = 1
YT2 = QD2*Bk11*rd2T*((1.5*pow(cos(Tangle),2)-.5));
YT4 = QD4*Bk12*rd4T*(35./8.* pow(cos(Tangle),4) - 30./8.*pow(cos(Tangle),2) + 3./8.);
YT = YT0 + YT2 + YT4;
YE = 1 + a2E*(1.5*pow(cos(Tangle),2)-.5) + a4E*(35./8.* pow(cos(Tangle),4) - 30./8.*pow(cos(Tangle),2) + 3./8.);
X2_total += pow((YT- YE),2)/(200.);
}
//Figure out denominator.
chisqr.push_back(log(X2_total));
tdelta.push_back(delta);
X2_total = 0.;
}
/*
A2T = QD2*Bk11*rd2T;
a2T = A2T/A0E;
//using the idea that A0E is our normalizer.
@ -655,15 +686,10 @@ int main(int argc,char **argv){
A4T = QD4*Bk12*rd4T;
a4T = A4T/A0E;
X22 = pow(abs(a2E -a2T),2)/(abs(a2T));
X24 = pow(abs(a4E -a4T),2)/(abs(a4T));
X22 = pow((a2E -a2T),2)/(a2T);
X24 = pow((a4E -a4T),2)/(a4T);
X2_total = (X22 + X24)/2;
chisqr.push_back(-log(X2_total));
tdelta.push_back(atan_delta);
}
*/
vector<double> Theta;
vector<double> AD_I;
double ad_start = 0.;

84
GUI_AD.h
View File

@ -239,18 +239,22 @@ int HistoGUIad::Draw_Fit(double x_low_win, double y_low_win, double x_hi_win, do
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<x.size(); i++){
// max_x = x[0];
// max_y = y[0];
// min_x = x[0];
// min_y = y[0];
//hard set graphical boundaries
max_x = 3.1415;
max_y = 2.;
min_x = -.1;
min_y = -.1;
/* for(int i=0; i<x.size(); i++){
if(x[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);
@ -327,21 +331,28 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
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<x.size(); i++){
// max_x = x[0];
// max_y = y[0];
// min_x = x[0];
// min_y = y[0];
//hard set graphical boundaries
max_x = 3.1415;
max_y = 2.;
min_x = -.1;
min_y = -.1;
/* for(int i=0; i<x.size(); i++){
if(x[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 );
*/
// 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);
@ -352,7 +363,7 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
//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);
//printf("y_offset = %f\n", y_offst);
double axis_x = (0. + x_offset) / width_scale;
double axis_y = (0. + y_offset) / height_scale;
@ -387,15 +398,16 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
x_wid = (x[i] + x_offset) / width_scale;
y_wid = (y[i] + y_offset) / height_scale;
y_errors_wid = ((y_errors[i]/A0) / height_scale);
// 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 -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);
//draws the point
XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -3, y_wid -3, 6, 6);
//these are the vertical error bars.
// 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);
// these are the caps of the 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);
@ -425,14 +437,14 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
for(int i=0; i < (int) width; i += w_step){
double x_val = i * width_scale - x_offset;
sprintf(axis_val, "%.1f", x_val);
XDrawString(disp, wind, DefaultGC(disp, screen), i, axis_y + 10/A0, axis_val, strlen(axis_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/A0, i, axis_val, strlen(axis_val));
XDrawString(disp, wind, DefaultGC(disp, screen), axis_x + 10, i, axis_val, strlen(axis_val));
}
@ -449,15 +461,15 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
y_wid = (y[i] + y_offset) / height_scale;
y_errors_wid = ((y_errors[i]/A0) / height_scale);
// x_wid2 = (x[i + 1] + x_offset) / width_scale;
// y_wid2 = (y[i + 1] + y_offset) / height_scale;
// printf("(%f, %f), (%f,%f)\n", x_wid,y_wid,x_wid2,y_wid2);
// XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2);
//draws the point
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.
//these are the vertical error bars
// 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);
//These are the error bar caps;.
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);
}