ex-4: χ² test and error estimation completed

ex-4/lib.c

@@ -0,0 +1,86 @@
+#include "lib.h"
+
+
+// Minimization wrapper.
+//
+double chi2(double p_max, void* params)
+  {
+  // Pars parameters.
+  struct parameters p = *((struct parameters*) params);
+  struct bin* histo = p.histo;
+  size_t n = p.n;
+  double step = p.step;
+
+  // Compute χ².
+  double chi = 0;
+  double expecto = 0; 
+  for (size_t i = 0; i < n; i++)
+    {
+    expecto = expected((i + 0.5) * step, p_max);
+    chi += pow(histo[i].sum - expecto, 2)/expecto;
+    };
+  return chi;
+  }
+
+
+// Expected function.
+//
+double expected (double x, double p_max)
+  {
+  // When p_max < x, the argument under sqrt is negative. Return great number
+  // to lead away the minimization.
+  if (p_max < x)
+    {
+    double num = 500;
+    return num;
+    }
+  return x * log(p_max/x)/atan(sqrt(fabs(pow(p_max, 2)/pow(x,2) - 1)));
+  }
+
+
+// First derivative of the expected function.
+//
+double exp1d (double x, double p_max)
+  {
+  // When p_max < x, the argument under sqrt is negative. Return great number
+  // to lead away the minimization.
+  if (p_max < x)
+    {
+    double num = 500;
+    return num;
+    }
+  double a = p_max;
+  double A = sqrt(fabs(pow(a, 2)/pow(x,2) - 1));
+  double B = atan(A);
+  double C = log(a/x);
+  return x/(a*B) - x*C/(a*A*pow(B, 2));
+  }
+
+//
+// Second derivative of the expected function.
+//
+double exp2d (double x, double p_max)
+  {
+  // When p_max < x, the argument under sqrt is negative. Return great number
+  // to lead away the minimization.
+  if (p_max < x)
+    {
+    double num = 500;
+    return num;
+    }
+  double a = p_max;
+  double a2 = pow(a, 2);
+
+  double A = sqrt(fabs(a2/pow(x,2) - 1));
+  double B = atan(A);
+  double C = log(a/x);
+  double D = 1 - a2/pow(x,2);
+
+  double E = -2/(a2*D*B);
+  double F = 1/(a2*A);
+  double G = 1/(pow(x, 2)*pow(A, 3));
+  double H = 1/a2 + 2/(a2 * A * B);
+
+  return x*((E + F + G) * C/B - H)/B;
+  }
+

ex-4/lib.h

@@ -0,0 +1,39 @@
+#pragma once
+
+#include <math.h>
+#include <stdlib.h>
+
+// Histogram struct.
+//
+struct bin
+  {
+  size_t amo;  // Amount of events in the bin.
+  double sum;  // Sum of |p_v|s of all the events in the bin.
+  };
+
+// Minimization struct.
+//
+struct parameters
+  {
+  struct bin* histo;  // Histogram
+  size_t n;           // Number of bins
+  double step;        // Bin width
+  };
+
+/////////////////////////////////////////////////////////////////
+
+// Minimization wrapper.
+//
+double chi2(double p_max, void* params);
+
+// Expected function.
+//
+double expected (double x, double p_max);
+
+// First derivative of the expected function.
+//
+double exp1d (double x, double p_max);
+
+// Second derivative of the expected function.
+//
+double exp2d (double x, double p_max);

ex-4/main.c

@@ -1,24 +1,21 @@
+#include "lib.h"
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
 #include <gsl/gsl_rng.h>
+#include <gsl/gsl_min.h>
+#include <gsl/gsl_deriv.h>
 
-int main(int argc, char **argv)
+
+// Process CLI arguments.
+//
+int parser(size_t *N, size_t *n, double *p_max, char argc, char **argv)
   {
-
-  // Set default options.
-  //
-  size_t N = 50000;          // number of events.
-  size_t n = 50;             // number of bins.
-  double p_max = 10;         // maximum value of momentum module.
-
-  // Process CLI arguments.
-  //
   for (size_t i = 1; i < argc; i++)
     {
-         if (!strcmp(argv[i], "-N")) N     = atol(argv[++i]);
+         if (!strcmp(argv[i], "-N")) *N     = atol(argv[++i]);
-    else if (!strcmp(argv[i], "-n")) n     = atol(argv[++i]);
+    else if (!strcmp(argv[i], "-n")) *n     = atol(argv[++i]);
-    else if (!strcmp(argv[i], "-p")) p_max = atof(argv[++i]);
+    else if (!strcmp(argv[i], "-p")) *p_max = atof(argv[++i]);
     else
       {
       fprintf(stderr, "Usage: %s -[hiIntp]\n", argv[0]);
@@ -29,6 +26,25 @@ int main(int argc, char **argv)
       return EXIT_FAILURE;
       }
     }
+  return EXIT_SUCCESS;
+  }
+
+
+int main(int argc, char **argv)
+  {
+
+  // Set default options.
+  //
+  size_t N = 50000;          // number of events.
+  size_t n = 50;             // number of bins.
+  double p_max = 10;         // maximum value of momentum module.
+  int res = parser(&N, &n, &p_max, argc, argv);
+  if (res == 0) printf("\nGenerating histogram with:\n"
+                       "%ld points\n"
+                       "%ld bins\n"
+                       "p_max = %.3f\n\n", N, n, p_max);
+
+  // printf("step: \t%.5f\n", step);
 
   // Initialize an RNG.
   //
@@ -51,7 +67,7 @@ int main(int argc, char **argv)
   //   p_v = p⋅cos(θ)
   //   p_h = p⋅sin(θ)
   //
-  // The histogram will be updated this way.
+  // The histogram is updated this way.
   // The j-th bin where p_h goes in is given by:
   //
   //   step = p_max / n
@@ -59,13 +75,8 @@ int main(int argc, char **argv)
   //
   // Thus an histogram was created and a structure containing the number of
   // entries in each bin and the sum of |p_v| in each of them is created and
-  // filled while generating the events.
+  // filled while generating the events (struct bin).
   //
-  struct bin
-    {
-    size_t amo;  // Amount of events in the bin.
-    double sum;  // Sum of |p_v|s of all the events in the bin.
-    };
   struct bin *histo = calloc(n, sizeof(struct bin));
 
   // Some useful variables.
@@ -94,22 +105,96 @@ int main(int argc, char **argv)
     //
     j = floor(p_h / step);
     b = &histo[j];
-    b->amo++;
+    b -> amo++;
-    b->sum += fabs(p_v);
+    b -> sum += fabs(p_v);
     }
 
   // Compute the mean value of each bin and print it to stodut
   // together with other useful things to make the histogram.
   //
-  printf("bins: \t%ld\n", n);
+  // printf("bins: \t%ld\n", n);
-  printf("step: \t%.5f\n", step);
+  // printf("step: \t%.5f\n", step);
   for (size_t i = 0; i < n; i++)
     {
-    histo[i].sum = histo[i].sum / histo[i].amo;
+    histo[i].sum = histo[i].sum / histo[i].amo;    // Average P_v
-    printf("\n%.5f", histo[i].sum);
+    //printf("\n%.5f", histo[i].sum);
     };
+  
+  // Compare the histigram with the expected function:
+  //
+  //   x * log(p_max/x)/arctan(sqrt(p_max^2/x^2 - 1))
+  //
+  // using the χ² test.
+  //
+  struct parameters params;
+  params.histo = histo;
+  params.n     = n;
+  params.step  = step;
 
-  // free memory
+  gsl_function func;
+  func.function = &chi2;
+  func.params = &params;
+
+  double min_p = 5;
+  double max_p = 15;
+
+  // Initialize minimization.
+  //
+  double x = 10;
+  int    max_iter = 100;
+  double prec = 1e-7;
+  int    status;
+  const gsl_min_fminimizer_type *T = gsl_min_fminimizer_brent;
+  gsl_min_fminimizer *s = gsl_min_fminimizer_alloc(T);
+  gsl_min_fminimizer_set(s, &func, x, min_p, max_p);
+
+  // Minimization.
+  //
+  for (int iter = 0; status == GSL_CONTINUE && iter < max_iter; iter++)
+    {  
+    status = gsl_min_fminimizer_iterate(s);
+    x      = gsl_min_fminimizer_x_minimum(s);
+    min_p  = gsl_min_fminimizer_x_lower(s);
+    max_p  = gsl_min_fminimizer_x_upper(s);
+    status = gsl_min_test_interval(min_p, max_p, 0, prec);
+    }
+
+  double result = x;
+  printf("p_max: %.7f\n", result);
+
+  // Compute the second derivative of χ² in its minimum for the result error.
+  //   
+  //   p_max = α
+  //
+  //           (Ei - Oi)²
+  //   χ² = Σi ----------
+  //               Ei
+  //
+  //              /    Oi² \
+  //   ∂αχ² = Σi | 1 - ---  | ∂αE
+  //              \    Ei² /
+  //
+  //               /         Oi²         /    Oi² \  \
+  //   ∂²αχ² = Σi | (∂αE)² 2 --- + ∂²αE | 1 - ---  | |
+  //               \         Ei³         \    Ei² /  /
+  //
+  double expecto, A, B;
+  double error = 0;
+  for (size_t i = 0; i < n; i++)
+    {
+    x = (i + 0.5) * step;
+    expecto = expected(x, result);
+    A = 2 * pow(exp1d(x, result) * histo[i].sum / expecto, 2);
+    B = exp2d(x, result) * (1 - pow((histo[i].sum / expecto), 2));
+    error = error + A + B;
+    };
+    error = 1/error;
+  printf("ΔP_max: %.7f\n\n", error);
+
+
+  // Free memory.
+  //
+  gsl_min_fminimizer_free(s);
   gsl_rng_free(r);
   free(histo);
 

makefile

@@ -21,7 +21,7 @@ ex-3/bin/main: ex-3/main.c ex-3/common.c ex-3/likelihood.c ex-3/chisquared.c
 	$(CCOMPILE)
 
 ex-4: ex-4/bin/main
-ex-4/bin/main: ex-4/main.c
+ex-4/bin/main: ex-4/main.c ex-4/lib.c
 	$(CCOMPILE)
 
 ex-5: ex-5/bin/casino ex-5/bin/manual ex-5/bin/trifecta