Fully functional AD program

AD.cxx

@@ -14,8 +14,15 @@ using namespace std;
 #ifndef __CINT__ 
 
 int main(int argc,char **argv){
+/*
+    for( float i = -2.; i < 6.; i += 1.0){
+        printf("%.0f! = %f \n",i, factorial(i));
+    } 
 
-
+    for( float m1 = -5.5; m1 <= 5.5; m1 += 1){
+        printf("m1: %.1f, %f \n", m1, CGcoeff(2, 0., 5.5, m1, 5.5, -m1));
+    }
+*/
     double j1, j2;
 
 
@@ -98,9 +105,9 @@ int main(int argc,char **argv){
         //then calucate QD2 and QD4, and replace the 0 with them. 
     }
 
-    double QD2 =  QK2(Energy, detradius, targetdistance, detthickness); ;
+    double QD2 =  QK2(Energy, detradius, targetdistance, detthickness);
 
-    double QD4 =  QK4(Energy, detradius, targetdistance, detthickness); ;
+    double QD4 =  QK4(Energy, detradius, targetdistance, detthickness);
 
     //	cout << "QD2 = " << " " << QD2 << "  QD4 = " << " " << QD4 << "\n";
 
@@ -199,6 +206,9 @@ int main(int argc,char **argv){
         dangler.push_back(aa*3.14159/180.);
         printf("dangle = %lf\n",dangler[i]);
     }
+    
+
+
     // if you need to scale the y data by sin(theta)	
 
 
@@ -365,7 +375,7 @@ int main(int argc,char **argv){
 
     double cg1, cg2 = 0.;	
 
-    double Bk11, Bk12 = 0.;
+    double Bk10, Bk11, Bk12 = 0.;
     if(Sigma == 0){
         if(IS == 1){
             A[0] = j1;
@@ -481,33 +491,24 @@ int main(int argc,char **argv){
                 II = j1 - am11; 
                 A[3] = am11; 
                 A[4] = -am11; 
+                
+                
 
-
+                cgg = CGcoeff(i,0,j1,am11,j1,-am11);
-                if(isnan(CG2(A))){ // CG2(A) is NaN
+               // cgg = CG2(A);
-
+              //  printf("K : %d, j1 : %.1f m1 : %.1f, %f \n", i, j1, am11, cgg);    
-                    //  cout << "Warning, a CG coefficient returned NaN, set to 0.\n";
+            //    cout << "CG =" << cgg << "\n";
-                    cgg = 0.0; 
+             //   cout << "m1 =" << am11 << "\n";
-
+               // cout << "m2 =" << -am11 << "\n";
-                }else{ //CG2(A) is not NaN.  
-
-                    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;
                 }else if(i == 4){ Bk12 = Bk12 + tTerm;} 
             }
         }
+
+//        cout << "Bk10" << " " << Bk10 << "\n";
         cout << "Bk11" << " " << Bk11 << "\n";
         cout << "Bk12" << " " << Bk12 << "\n";
 
@@ -661,7 +662,7 @@ int main(int argc,char **argv){
     vector<double> YE_point;
     double YT0,YT2,YT4,YT,YE = 0.;
 
-    double Y_err = 10.; 
+    double Y_err; 
 
     //delta loop;
     for(int i = 0; i< points; i++){
@@ -669,18 +670,19 @@ int main(int argc,char **argv){
         delta = tan(atan_delta);
         //now sum over all theta;
 
-        for(int j = 0; j < points; j++){
+        for(int j = 0; j < dangler.size(); j++){
-            Tangle = j*step;
+            Tangle = dangler[j];
+            Y_err = deydata[j];
 
             //calculate Y_intensity_theory and and experiment and sum over theta
 
-            rd0T = (1 + 2*delta + pow(delta,2))/(1+pow(delta,2));
 
             rd2T = (rk01 + 2*delta*rk11 + pow(delta,2)*rk21)/(1+pow(delta,2));
 
             rd4T = (rk02 + 2*delta*rk12 + pow(delta,2)*rk22)/(1+pow(delta,2));
 
-            YT0 = rd0T; // k = 0;
+            YT0 = 1; // k = 0; Bk10*?
+
 
             YT2 = QD2*Bk11*rd2T*((1.5*pow(cos(Tangle),2)-.5));
 
@@ -690,7 +692,7 @@ int main(int argc,char **argv){
 
             //    YT_tot[i] += YT;
 
-            YE = 1 + a2E*(1.5*pow(cos(Tangle),2)-.5) + a4E*(35./8.* pow(cos(Tangle),4) - 30./8.*pow(cos(Tangle),2) + 3./8.);
+            YE = 1. + a2E*(1.5*pow(cos(Tangle),2)-.5) + a4E*(35./8.* pow(cos(Tangle),4) - 30./8.*pow(cos(Tangle),2) + 3./8.);
 
             X2_total += pow((YT- YE),2)/((dangler.size()-1)*(pow(Y_err,2)));
             //           X2_total[i] = X2_total[i] + pow((YE),2)/((dangler.size()-1)*(pow(Y_err,2)));
@@ -709,10 +711,10 @@ int main(int argc,char **argv){
 
 
     }
-    /*
+// this is for viewing the theoretical angular distribution as a function of constant delta.    
        for(int k = 0; k < dangler.size(); k++){
        Tangle = dangler[k];
-       delta = 2.;
+       delta = .8;
 
     //calculate Y_intensity_theory and and experiment and sum over theta
 
@@ -730,14 +732,14 @@ int main(int argc,char **argv){
 
     YT = YT0 + YT2 + YT4;
     //        cout << "YT = " <<YT << "\n"; 
-    //            YT_point.push_back(YT);
+    YT_point.push_back(YT);
 
     YE = 1 + a2E*(1.5*pow(cos(Tangle),2)-.5) + a4E*(35./8.* pow(cos(Tangle),4) - 30./8.*pow(cos(Tangle),2) + 3./8.);
     YE_point.push_back(YE);
 
     //  X2_total += pow((YT- YE),2)/((dangler.size()-1)*(pow(Y_err,2)));
     }
-    */
+    
     /*
        A2T = QD2*Bk11*rd2T;
        a2T = A2T/A0E;
@@ -832,6 +834,8 @@ int main(int argc,char **argv){
             //y will be log(X^2); 
 
             gui.SetData(tdelta,chisqr);
+           // gui.SetData(dangler,YE_point);
+           // gui.SetData(dangler,YT_point);
             //gui.SetData(Theta,AD_I);
             gui.Init();
             gui.Loop();
@@ -880,8 +884,8 @@ int main(int argc,char **argv){
             //x will be arctan of mixing ratio. 
             //y will be log(X^2); 
 
-            gui.SetData(dangler,YT_point);
+            //gui.SetData(dangler,YT_point);
-            // gui.SetData(tdelta,chisqr);
+            gui.SetData(tdelta,chisqr);
             //gui.SetData(Theta,AD_I);
             gui.Init();
             gui.Loop();

Cleb.h

@@ -1,12 +1,55 @@
 
 #include "global.h"
 //To compile : g++ AD.cxx -o {Input Executable Name} -lX11
-
 //#include "Coeff.h"
 
 using namespace std;
 
+double factorial(double n){
+  if( n < 0 ) return -100.;
+  return (n == 1. || n == 0.) ? 1. : factorial(n-1) * n ;
+}
+ 
+double CGcoeff(double J, double m, double J1, double m1, double J2, double m2){
+  // (J1,m1) + (J2, m2) = (J, m)
+ 
+  if( m != m1 + m2 ) return 0;
+ 
+  double Jmin = abs(J1 - J2);
+  double Jmax = J1+J2;
+ 
+  if( J < Jmin || Jmax < J ) return 0;
+ 
+  double s0 = (2*J+1.) * factorial(J+J1-J2) * factorial(J-J1+J2) * factorial(J1+J2-J) / factorial(J+J1+J2 + 1.);
+  s0 = sqrt(s0);
+ 
+  double s = factorial(J +m ) * factorial(J -m);
+  double s1 = factorial(J1+m1) * factorial(J1-m1);
+  double s2 = factorial(J2+m2) * factorial(J2-m2);
+  s = sqrt(s * s1 * s2);
+ 
+ // printf(" s0, s = %f , %f \n", s0, s);
 
+ 
+  int kMax = min( min( J1+J2-J, J1 - m1), J2 + m2);
+ 
+  double CG = 0.;
+  for( int k = 0; k <= kMax; k++){
+    double k1 = factorial(J1+J2-J-k);
+    double k2 = factorial(J1-m1-k);
+    double k3 = factorial(J2+m2-k);
+    double k4 = factorial(J - J2 + m1 +k);
+    double k5 = factorial(J - J1 - m2 +k);
+    double temp = pow(-1, k) / (factorial((double)k) * k1 * k2 * k3 * k4 * k5);
+    if( k1 == -100. || k2 == -100. || k3 == -100. || k4 == -100. || k5 == -100. ) temp = 0.;
+ 
+  //  printf(" %d |%.12f|%.12f|%.12f|%.12f|%.12f|%.12f|%f \n", k,k1,k2,k3,k4,k5, temp,factorial(double(k)));
+    CG += temp;
+  }
+ 
+  return s0*s*CG;
+ 
+}
 
 //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){
@@ -58,7 +101,10 @@ double CG2(double A[6]){
 	double m1 = A[3];
 	double m2 = A[4];
 	double M = A[5]; 
-	if(M != m1 + m2) return 0;
+	
+    return CGcoeff(J,M,j1,m1,j2,m2);
+    
+    if(M != m1 + m2) return 0;
 
 	double Jmin  =  abs(j1 - j2);
 	double Jmax  =  j1+j2;

GUI_AD.h

@@ -370,10 +370,15 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
 		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* x_axisl = (char*)"arctan( delta ) in Radians";
+        char* y_axisl = (char*)"Log Chi-Squared";
 		char axis_val[4];
 		int w_step = width / 10;
-		for(int i=0; i < (int) width; i += w_step){
+
+		XDrawString(disp, wind, DefaultGC(disp, screen), max_x *0.9, min_y,x_axisl, strlen(x_axisl));
+		XDrawString(disp, wind, DefaultGC(disp, screen), max_y *0.9, min_x,y_axisl, strlen(y_axisl));
+        
+        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));
@@ -397,8 +402,8 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
 		for(int i=0; i < x.size() - 1; i++){
 			x_wid  = (x[i] + x_offset) / width_scale;
 			y_wid  = (y[i] + y_offset) / height_scale;
-		//	y_errors_wid = ((y_errors[i]/A0) / height_scale);
+			y_errors_wid = ((y_errors[i]/A0) / height_scale); //errors will also be scaled by A0
-			y_errors_wid = ((y_errors[i]) / height_scale);
+	    //	y_errors_wid = ((y_errors[i]) / height_scale);
             
             //draws the point
             XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -3, y_wid -3, 6, 6);
@@ -433,6 +438,8 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
 		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* x_axisl = (char*)"arctan( delta ) in Radians";
+        char* y_axisl = (char*)"Log Chi-Squared";
 		char axis_val[4];
 		int w_step = width / 10;
 		for(int i=0; i < (int) width; i += w_step){
@@ -441,7 +448,10 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
 			XDrawString(disp, wind, DefaultGC(disp, screen), i, axis_y + 10, axis_val, strlen(axis_val));
 		}
 
-		int h_step = height / 10;
+		XDrawString(disp, wind, DefaultGC(disp, screen), max_x *0.9, min_y,x_axisl, strlen(x_axisl));
+		XDrawString(disp, wind, DefaultGC(disp, screen), max_y *0.9, min_x,y_axisl, strlen(y_axisl));
+        
+        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);
@@ -461,8 +471,8 @@ int HistoGUIad::DrawData(double x_low_win, double y_low_win, double x_hi_win, do
 			x_wid  = (x[i] + x_offset) / width_scale;
 			y_wid  = (y[i] + y_offset) / height_scale;
 	
-		//	y_errors_wid = ((y_errors[i]/A0) / height_scale);
+			y_errors_wid = ((y_errors[i]/A0) / height_scale); //Error bars scaled by A0 
-			y_errors_wid = ((y_errors[i]) / height_scale);
+	    //  y_errors_wid = ((y_errors[i]) / height_scale);
             
             //draws the point
 			XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -3, y_wid -3, 6, 6);