#include <stdlib.h>#include <stdio.h>#include <math.h>#include <gsl/gsl_multimin.h>#include "globes/globes.h"#include "nova-globes.h"Go to the source code of this file.
Defines | |
| #define | NSAMPLEX 300 |
| #define | NSAMPLEY 20 |
| #define | SSTTLOW 0.0 |
| #define | SSTTHIGH 0.15 |
Functions | |
| double | chisqr (const gsl_vector *v, void *params) |
| double | hier_sens_best_chisqr (double theta12, double theta13, double theta23, double deltacp, double sdms, double ldms, double hier, double *thebestsstt, double *thebestdcp, double *thebestth23) |
| void | hier_sens () |
Variables | |
| glb_params | gs_test_values |
|
|
Definition at line 11 of file hier_sens.c. |
|
|
Definition at line 12 of file hier_sens.c. |
|
|
Definition at line 14 of file hier_sens.c. |
|
|
Definition at line 13 of file hier_sens.c. |
|
||||||||||||
|
Definition at line 18 of file hier_sens.c. References chisqr(), and gs_test_values. Referenced by chisqr(), and hier_sens(). 00019 {
00020 double theta13 = fabs(gsl_vector_get(v, 0));
00021 double dcp = gsl_vector_get(v, 1);
00022 double theta23 = gsl_vector_get(v, 2);
00023 double penalty = 0.0;
00024
00025 glbSetOscParams(gs_test_values, theta13, GLB_THETA_13);
00026 glbSetOscParams(gs_test_values, theta23, GLB_THETA_23);
00027 glbSetOscParams(gs_test_values, dcp, GLB_DELTA_CP);
00028
00029 double chisqr = glbChiSys(gs_test_values, GLB_ALL, GLB_ALL);
00030
00031 return chisqr;
00032 }
|
|
|
Definition at line 140 of file hier_sens.c. References chisqr(), gs_test_values, hier_sens_best_chisqr(), run_name(), SSTTHIGH, SSTTLOW, test_dcp(), test_hierarchy(), test_ldms(), test_sdms(), test_theta12(), test_theta13(), and test_theta23(). Referenced by main(). 00141 {
00142 int i, j, k;
00143 double theta12 = test_theta12();
00144 double theta23 = test_theta23();
00145 double theta13 = test_theta13();
00146 double deltacp = test_dcp();
00147 double sdms = test_sdms();
00148 double ldms = test_ldms();
00149 double sstt13;
00150 double dcp;
00151 double chisqr;
00152 double ssttmin;
00153 double dcpmin;
00154 double th23min;
00155 double chisqr_best;
00156 double sstt_best;
00157 double dcp_best;
00158
00159 glb_params true_values = glbAllocParams();
00160 if (test_hierarchy()>0.0) {
00161 glbDefineParams(true_values, theta12, theta13, theta23, deltacp,
00162 sdms, sdms+ldms);
00163 }
00164 else {
00165 glbDefineParams(true_values, theta12, theta13, theta23, deltacp,
00166 sdms, -ldms);
00167 }
00168 glbSetDensityParams(true_values, 1.0, GLB_ALL);
00169
00170 char fname[256];
00171 sprintf(fname,"%s_hier_sens.txt",run_name());
00172 FILE* fp = fopen(fname,"w");
00173
00174 /* Go to every point in sstt13-dcp space, compute new "true" rates
00175 * and then find the best solution using the opposite choice for the
00176 * mass hierarchy
00177 */
00178 gs_test_values = glbAllocParams();
00179 for (i=0; i<NSAMPLEX; ++i) {
00180 sstt13 = SSTTLOW+(SSTTHIGH-SSTTLOW)*((float)i+0.5)/NSAMPLEX;
00181 theta13 = 0.5*asin(sqrt(sstt13));
00182 for (j=0; j<NSAMPLEY; ++j) {
00183 dcp = 2.0*M_PI*((float)j+0.5)/NSAMPLEY;
00184
00185 glbSetOscParams(true_values, theta13, GLB_THETA_13);
00186 glbSetOscParams(true_values, dcp, GLB_DELTA_CP);
00187 if (test_hierarchy()>0.0) {
00188 glbSetOscParams(true_values, sdms+ldms, GLB_DM_31);
00189 }
00190 else {
00191 glbSetOscParams(true_values, -ldms, GLB_DM_31);
00192 }
00193 glbSetOscillationParameters(true_values);
00194 glbSetRates();
00195
00196 chisqr = hier_sens_best_chisqr(theta12,
00197 theta13,
00198 theta23,
00199 dcp,
00200 sdms,
00201 ldms,
00202 -test_hierarchy(),
00203 &ssttmin,
00204 &dcpmin,
00205 &th23min);
00206 fprintf(fp, "%f %f %f %f %f %f\n",
00207 sstt13, dcp/M_PI, chisqr, ssttmin, dcpmin, th23min);
00208 } /* j values of dcp */
00209 } /* i values of sstt13 */
00210 fclose(fp);
00211 }
|
|
||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 36 of file hier_sens.c. References gs_test_values. Referenced by hier_sens(). 00046 {
00047 int i;
00048 int which_dcp;
00049 int which_oct;
00050 double bestchis = 1e30;
00051 double seed_theta13;
00052 double seed_deltacp;
00053 double seed_theta23;
00054 double step_size = 0.1;
00055
00056 if (hier>0.0) {
00057 glbDefineParams(gs_test_values,
00058 theta12, theta13, theta23, deltacp,
00059 sdms, sdms+ldms);
00060 }
00061 else {
00062 glbDefineParams(gs_test_values,
00063 theta12, theta13, theta23, deltacp,
00064 sdms, -ldms);
00065 }
00066 glbSetDensityParams(gs_test_values, 1.0, GLB_ALL);
00067
00068 seed_theta23 = theta23;
00069 seed_theta13 = theta13;
00070 seed_deltacp = deltacp;
00071 for (i=0; i<9; ++i) {
00072 if (i==8) {
00073 seed_theta13 = *thebestsstt;
00074 seed_deltacp = *thebestdcp;
00075 seed_theta23 = *thebestth23;
00076 step_size = 1.0;
00077 }
00078 else {
00079 seed_theta13 = theta13;
00080 which_dcp = (i/4);
00081 which_oct = (i%2);
00082 seed_deltacp = deltacp - (float)which_dcp*0.25*M_PI;
00083 while (seed_deltacp<0) seed_deltacp += 2.0*M_PI;
00084 while (seed_deltacp>2.0*M_PI) seed_deltacp -= 2.0*M_PI;
00085 if (which_oct) seed_theta23 = theta23;
00086 else seed_theta23 = 0.5*M_PI - theta23;
00087 step_size = 0.5;
00088 }
00089 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
00090 gsl_multimin_fminimizer *s = NULL;
00091 gsl_vector *ss, *x;
00092 gsl_multimin_function minex_func;
00093
00094 size_t iter = 0;
00095 int status;
00096 double size;
00097
00098 /* Starting point */
00099 x = gsl_vector_alloc(3);
00100 gsl_vector_set(x, 0, seed_theta13);
00101 gsl_vector_set(x, 1, seed_deltacp);
00102 gsl_vector_set(x, 2, seed_theta23);
00103
00104 /* Set initial step sizes to 1 */
00105 ss = gsl_vector_alloc(3);
00106 gsl_vector_set_all(ss, 0.2);
00107
00108 /* Initialize method and iterate */
00109 minex_func.n = 3;
00110 minex_func.f = &chisqr;
00111 minex_func.params = 0;
00112
00113 s = gsl_multimin_fminimizer_alloc(T, 3);
00114 gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
00115
00116 do {
00117 iter++;
00118 status = gsl_multimin_fminimizer_iterate(s);
00119 if (status) break;
00120 size = gsl_multimin_fminimizer_size(s);
00121 status = gsl_multimin_test_size(size, 1E-4);
00122 }
00123 while (status == GSL_CONTINUE && iter < 500);
00124
00125 if (s->fval<bestchis) {
00126 *thebestsstt = gsl_vector_get(s->x, 0);
00127 *thebestdcp = gsl_vector_get(s->x, 1);
00128 *thebestth23 = gsl_vector_get(s->x, 2);
00129 bestchis = s->fval;
00130 }
00131 gsl_vector_free(x);
00132 gsl_vector_free(ss);
00133 gsl_multimin_fminimizer_free(s);
00134 }
00135 return bestchis;
00136 }
|
|
|
Definition at line 16 of file hier_sens.c. Referenced by chisqr(), hier_sens(), and hier_sens_best_chisqr(). |
1.3.9.1