Forward the density for the backward sampling of vertices.

niess · Jul 4, 2017 · 800631f · 800631f
1 parent bce07bc
commit 800631f
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 9 deletions.
diff --git a/include/ent.h b/include/ent.h
@@ -248,7 +248,13 @@ typedef double(ent_random_cb)(struct ent_context * context);
  * @param daughter   The daughter state.
  * @return The relative density of the ancester.
  *
- * TODO: document.
+ * This callback **must** be provided for backward Monte-Carlo. It is expected
+ * to return an _a priori_ density estimate for an ancester at the daughter's
+ * location. Returning `0` or less indicates a null density, which disables
+ * the corresponding channel.
+ *
+ * **Warning** : if multiple contexts are used the user must ensure that this
+ * callback is thread safe.
  */
 typedef double(ent_ancester_cb)(struct ent_context * context,
     enum ent_pid ancester, struct ent_state * daughter);

diff --git a/src/ent.c b/src/ent.c
@@ -2527,7 +2527,8 @@ static enum ent_return ancester_draw(struct ent_context * context,
 /* Randomise the ancester at a backward vertex. */
 static enum ent_return transport_ancester_draw(struct ent_physics * physics,
     struct ent_context * context, struct ent_state * daughter,
-    struct ent_medium * medium, enum ent_pid * ancester, int * proget)
+    struct ent_medium * medium, double density, enum ent_pid * ancester,
+    int * proget)
 {
         /* Build the interpolation or extrapolation factors for cross-sections.
          */
@@ -2582,8 +2583,11 @@ static enum ent_return transport_ancester_draw(struct ent_physics * physics,
 
                 /* True decay process from a muon. */
                 if (medium->density != NULL) {
-                        double density; /* TODO: forward from the tranport. */
-                        medium->density(medium, daughter, &density);
+                        if (density <= 0.) {
+                                medium->density(medium, daughter, &density);
+                                if (density <= 0.)
+                                        return ENT_RETURN_DOMAIN_ERROR;
+                        }
                         if (apid != ENT_PID_NU_TAU) {
                                 int mother;
                                 if (apid == ENT_PID_MUON) {
@@ -2822,7 +2826,7 @@ enum ent_return vertex_forward(struct ent_physics * physics,
 /* Vertex randomisation in backward Monte-Carlo. */
 static enum ent_return vertex_backward(struct ent_physics * physics,
     struct ent_context * context, struct ent_state * state,
-    struct ent_medium * medium, enum ent_process process,
+    struct ent_medium * medium, double density, enum ent_process process,
     struct ent_state * product, enum proget_index * proget_)
 {
         /* Get the ancester, the process and the target. */
@@ -2834,7 +2838,8 @@ static enum ent_return vertex_backward(struct ent_physics * physics,
                  */
                 enum ent_return rc;
                 if ((rc = transport_ancester_draw(physics, context, state,
-                         medium, &ancester, &proget)) != ENT_RETURN_SUCCESS)
+                         medium, density, &ancester, &proget)) !=
+                    ENT_RETURN_SUCCESS)
                         return rc;
         } else if ((process == ENT_PROCESS_DIS_CC) ||
             (process == ENT_PROCESS_DIS_NC)) {
@@ -2997,8 +3002,8 @@ enum ent_return ent_vertex(struct ent_physics * physics,
                 ENT_RETURN(vertex_forward(
                     physics, context, state, medium, process, product));
         else
-                ENT_RETURN(vertex_backward(
-                    physics, context, state, medium, process, product, NULL));
+                ENT_RETURN(vertex_backward(physics, context, state, medium, -1.,
+                    process, product, NULL));
 }
 
 /* API interface for a Monte-Carlo tranport. */
@@ -3107,7 +3112,7 @@ enum ent_return ent_transport(struct ent_physics * physics,
                          */
                         enum proget_index proget;
                         rc = vertex_backward(physics, context, state, medium,
-                            ENT_PROCESS_NONE, product, &proget);
+                            density, ENT_PROCESS_NONE, product, &proget);
 
                         /* Then let us apply the effective weight for the
                          * transport.