Merge pull request #12 from oslocyclotronlab/chi2

Fixing chi2 print of temperature fitting in counting, see #11
oslocyclotronlab · Aug 2, 2018 · 8977457 · 8977457
2 parents dd51500 + f18ce25
commit 8977457
Showing 1 changed file with 13 additions and 7 deletions.
diff --git a/prog/counting.c b/prog/counting.c
@@ -77,7 +77,7 @@ int main()
     printf("\n");
     printf("  ______________________________________________________________ \r\n");
     printf(" |                                                              |\r\n");
-    printf(" |                    C O U N T I N G   1.8.4                   |\r\n");
+    printf(" |                    C O U N T I N G   1.8.4.2                 |\r\n");
     printf(" |                                                              |\r\n");
     printf(" |Program to normalize experimental nuclear level density (NLD) |\r\n");
     printf(" |  to NLD from known low energy levels and NLD extracted from  |\r\n");
@@ -116,6 +116,7 @@ int main()
     printf(" | Modified: 15 Mar 2017 New spin cutoff opt.(5) from Alexander |\r\n");
     printf(" | Modified: 01 Apr 2017 The upbend gSF is now lin. in log.     |\r\n");
     printf(" | Modified: 17 Apr 2017 TALYS nld spin-formated outputfile     |\r\n");
+    printf(" | Modified: 02 Aug 2018 Update printing of chi2 fit T          |\r\n");
     printf(" |______________________________________________________________|\r\n");
     printf("                                                                 \r\n");
 
@@ -685,14 +686,21 @@ int main()
         i = (int)(Tlow + .5);
         Tlow = (float)i/100.;
         Thigh = 1.2*T;
-        dTx = 0.01;
+        dTx = 0.005;
         Tx = Tlow;
+        // int freeL = L2 - L1 - 1; // -1 due to one fit parameter
+        int freeH = H2 - H1 - 1; // -1 due to one fit parameter
+        printf("\n Fitting T to chosen data in the higher region");
+        printf("\n     T                E0          Chi2          Chi2/free");
         while( Tx <= Thigh){
             E0x = Bn - Tx*(log(rho0)+log(Tx));
             searchT();
-            printf("\n T = %7.3f, E0 = %7.3f, Chi2_low = %9.3f, Chi2_high = %9.3f",Tx, E0x, chi2_lowx, chi2_highx);
+            // printf("\n T = %7.3f, E0 = %7.3f, Chi2 = %9.3f, Chi2/freeH = %9.3f",Tx, E0x, chi2_highx, chi2_highx/freeH);
+            printf("\n %7.3f \t %9.3f \t %7.3f \t %7.3f",Tx, E0x, chi2_highx, chi2_highx/freeH);
             Tx = Tx + dTx;
         }
+        printf("\nNumber of free parameter in upper region: %i", freeH);
+        printf("\nNote: Uncertainty can be rerived from intersection with (Chi2_min + 1)");
 
         printf("\n\nYou may change your choice of T from the info of the Chi2 tests above");
         printf("\nTemperature parameter T (MeV) <%7.3f>:",TRobin);
@@ -1635,7 +1643,7 @@ float anchorL(){
     float corr, corrbest=0., sum, sumbest, sum0=0., cc, free, dc2, rhoL;
     corr = 0.25;
     sumbest = 1.0e+21;
-    free = (float)L2 - (float)L1;
+    free = (float)L2 - (float)L1 - 1;
     if (free <= 0)free = 1.;
     for(j = 0; j <= 3750; j++){
         corr = corr + 0.001;
@@ -1649,7 +1657,6 @@ float anchorL(){
             dc2 =sqrt(dc2*dc2 + 1.*1.);
             if(dc2 > 0) sum = sum + ((cc-rhoL)*(cc-rhoL)/dc2);
         }
-        sum = sum/free;
         if(j == 499){
             sum0 = sum;
         }
@@ -1666,7 +1673,7 @@ float anchorH(){
     float corr, corrbest=0., sum, sumbest, sum0=0., cc, free, dc2;
     corr = 0.25;
     sumbest = 1.0e+21;
-    free = (float)H2 - (float)H1;
+    free = (float)H2 - (float)H1 -1;
     if (free <= 0)free = 1.;
     for(j = 0; j <= 3750; j++){
         corr = corr + 0.001;
@@ -1679,7 +1686,6 @@ float anchorH(){
             rhofg(Amass, ex, a, Tx, E1, E0x, isig, itemp, imodel, red, b1, b2);
             if(dc2 > 0) sum=sum+(cc-eta*rhox)*(cc-eta*rhox)/dc2;
         }
-        sum = sum/free;
         if(j == 499){
             sum0 = sum;
         }