schiff

commit 0f08e67c22ab94fed7ce2997f6fb35d8e55a740d
parent 51ee422cf8413d3c09b157508d1c939dc3a6cd31
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 23 Mar 2016 11:08:38 +0100

Use a more general truncated gaussian distribution

The min and max bounds used to check the validity of the parsed `mu'
value are used to reject or not the sampled gaussian value.

Diffstat:
M
src/schiff_args.c
 | 
4
+++-
M
src/schiff_geometry.c
 | 
16
+++++++++++-----
M
src/schiff_geometry.h
 | 
2
+-

3 files changed, 15 insertions(+), 7 deletions(-)
diff --git a/src/schiff_args.c b/src/schiff_args.c
@@ -506,9 +506,11 @@ parse_yaml_param_distribution
     /* Gaussian distribution */
     } else if(!strcmp((char*)key->data.scalar.value, "gaussian")) {
       res = parse_yaml_param_mu_sigma(filename, doc, val, "gaussian", min_val,
-        max_val, &param->data.gaussian.mu, &param->data.lognormal.sigma);
+        max_val, &param->data.gaussian.mu, &param->data.gaussian.sigma);
       if(res != RES_OK) return res;
       param->distribution = SCHIFF_PARAM_GAUSSIAN;
+      param->data.gaussian.range[0] = min_val;
+      param->data.gaussian.range[1] = max_val;
 
     /* Histogram distribution */
     } else if(!strcmp((char*)key->data.scalar.value, "histogram")) {
diff --git a/src/schiff_geometry.c b/src/schiff_geometry.c
@@ -129,11 +129,17 @@ histogram_sample
 }
 
 static FINLINE double
-ran_positive_gaussian(struct ssp_rng* rng, const double mu, const double sigma)
+ran_gaussian_truncated
+  (struct ssp_rng* rng,
+   const double mu,
+   const double sigma,
+   const double range[2])
 {
   double val;
-  ASSERT(rng && mu > 0);
-  do { val = ssp_ran_gaussian(rng, mu, sigma); } while(val <= 0);
+  ASSERT(rng && mu > 0 && mu >= range[0] && mu <= range[1]);
+  do {
+    val = ssp_ran_gaussian(rng, mu, sigma);
+  } while(val < range[0] || val > range[1] );
   return val;
 }
 
@@ -152,8 +158,8 @@ eval_param(const struct schiff_param* param, struct ssp_rng* rng)
         log(param->data.lognormal.sigma));
       break;
     case SCHIFF_PARAM_GAUSSIAN:
-      val = ran_positive_gaussian
-        (rng, param->data.gaussian.mu, param->data.gaussian.sigma);
+      val = ran_gaussian_truncated(rng, param->data.gaussian.mu,
+        param->data.gaussian.sigma, param->data.gaussian.range);
       break;
     case SCHIFF_PARAM_HISTOGRAM:
       val = histogram_sample
diff --git a/src/schiff_geometry.h b/src/schiff_geometry.h
@@ -32,7 +32,7 @@ struct schiff_param {
   union {
     double constant;
     struct { double mu, sigma; } lognormal;
-    struct { double mu, sigma; } gaussian;
+    struct { double mu, sigma, range[2]; } gaussian;
     struct { double *entries, lower, upper; } histogram;
   } data;
 };