Skip to content

Commit 27afd2c

Browse files
niessniess
committed
Add a callback for custom stepping actions
1 parent 1a9d027 commit 27afd2c

File tree

2 files changed

+111
-61
lines changed

2 files changed

+111
-61
lines changed

‎include/ent.h

+22-2
Original file line numberDiff line numberDiff line change
@@ -238,7 +238,7 @@ typedef double ent_medium_cb(struct ent_context * context,
238238
* callback is thread safe, e.g. by using independant streams for each context
239239
* or a locking mechanism in order to share a single random stream.
240240
*/
241-
typedef double(ent_random_cb)(struct ent_context * context);
241+
typedef double ent_random_cb(struct ent_context * context);
242242

243243
/**
244244
* Callback providing an a priori relative density of ancestors.
@@ -256,9 +256,27 @@ typedef double(ent_random_cb)(struct ent_context * context);
256256
* **Warning** : if multiple contexts are used the user must ensure that this
257257
* callback is thread safe.
258258
*/
259-
typedef double(ent_ancestor_cb)(struct ent_context * context,
259+
typedef double ent_ancestor_cb(struct ent_context * context,
260260
enum ent_pid ancestor, struct ent_state * daughter);
261261

262+
/**
263+
* Callback for custom actions during the Monte-Carlo stepping.
264+
*
265+
* @param context The operating Monte-Carlo context.
266+
* @param medium The medium at the current step.
267+
* @param state The particle state at the current step.
268+
* @return On success `ENT_RETURN_SUCCESS` must be returned otherwise any error
269+
* code can be returned.
270+
*
271+
* Providing this callback is optionnal. It is called at the end of major
272+
* Monte-Carlo steps when invoking `ent_transport`.
273+
*
274+
* **Warning** : if multiple contexts are used the user must ensure that this
275+
* callback is thread safe.
276+
*/
277+
typedef enum ent_return ent_stepping_cb(struct ent_context * context,
278+
struct ent_medium * medium, struct ent_state * state);
279+
262280
/**
263281
* Data for a Monte-Carlo context.
264282
*
@@ -272,6 +290,8 @@ struct ent_context {
272290
ent_random_cb * random;
273291
/** The ancestor callback, or `NULL` for forward Monte-Carlo. */
274292
ent_ancestor_cb * ancestor;
293+
/** A custom stepping action callback, or `NULL` if none. */
294+
ent_stepping_cb * stepping_action;
275295
/** A user supplied distance limit for the transport, or `0`. */
276296
double distance_max;
277297
/** A user supplied grammage limit for the transport, or `0`. */

‎src/ent.c

+89-59
Original file line numberDiff line numberDiff line change
@@ -299,8 +299,9 @@ static enum ent_return file_get_line_(struct file_buffer ** buffer, int skip)
299299
if (ptr[size - 1] == check) break;
300300

301301
/* Get more memory then ... */
302-
void * tmp = realloc(*buffer, sizeof(**buffer) +
303-
(*buffer)->size + FILE_BLOCK_SIZE);
302+
void * tmp = realloc(*buffer,
303+
sizeof(**buffer) + (*buffer)->size +
304+
FILE_BLOCK_SIZE);
304305
if (tmp == NULL) return ENT_RETURN_MEMORY_ERROR;
305306
*buffer = tmp;
306307
ptr = (*buffer)->line + (*buffer)->size - 1;
@@ -818,10 +819,10 @@ static double dcs_integrate(struct ent_physics * physics,
818819
for (j = 0; j < nx * N_GQ; j++) {
819820
const double x = xmin *
820821
exp((0.5 + 0.5 * xGQ[j % N_GQ] + j / N_GQ) *
821-
dlx);
822+
dlx);
822823
J += wGQ[j % N_GQ] * x *
823824
dcs_compute(physics, projectile, energy, Z,
824-
A, process, x, y);
825+
A, process, x, y);
825826
}
826827
I += wGQ[i % N_GQ] * y * J;
827828
}
@@ -842,9 +843,9 @@ static double dcs_integrate(struct ent_physics * physics,
842843
for (i = 0; i < ny * N_GQ; i++) {
843844
const double y = ymin *
844845
exp((0.5 + 0.5 * xGQ[i % N_GQ] + i / N_GQ) * dly);
845-
I += wGQ[i % N_GQ] * y * dcs_compute(physics,
846-
projectile, energy, Z, A,
847-
process, 0., y);
846+
I += wGQ[i % N_GQ] * y *
847+
dcs_compute(physics, projectile, energy, Z, A,
848+
process, 0., y);
848849
}
849850

850851
const double y1min = mu / energy;
@@ -853,9 +854,9 @@ static double dcs_integrate(struct ent_physics * physics,
853854
for (i = 0; i < ny * N_GQ; i++) {
854855
const double y1 = y1min *
855856
exp((0.5 + 0.5 * xGQ[i % N_GQ] + i / N_GQ) * dly1);
856-
I1 += wGQ[i % N_GQ] * y1 * dcs_compute(physics,
857-
projectile, energy, Z, A,
858-
process, 0., 1. - y1);
857+
I1 += wGQ[i % N_GQ] * y1 *
858+
dcs_compute(physics, projectile, energy, Z, A,
859+
process, 0., 1. - y1);
859860
}
860861

861862
return 0.5 * (I * dly + I1 * dly1);
@@ -1219,29 +1220,32 @@ enum ent_return ent_physics_cross_section(struct ent_physics * physics,
12191220
if ((rc = proget_compute(projectile, ENT_PID_PROTON,
12201221
process, &proget)) != ENT_RETURN_SUCCESS)
12211222
ENT_RETURN(rc);
1222-
*cross_section += Z * cross_section_compute(mode,
1223-
proget, cs0, cs1, p1, p2);
1223+
*cross_section += Z *
1224+
cross_section_compute(
1225+
mode, proget, cs0, cs1, p1, p2);
12241226
}
12251227
const double N = A - Z;
12261228
if (N > 0.) {
12271229
if ((rc = proget_compute(projectile, ENT_PID_NEUTRON,
12281230
process, &proget)) != ENT_RETURN_SUCCESS)
12291231
ENT_RETURN(rc);
1230-
*cross_section += N * cross_section_compute(mode,
1231-
proget, cs0, cs1, p1, p2);
1232+
*cross_section += N *
1233+
cross_section_compute(
1234+
mode, proget, cs0, cs1, p1, p2);
12321235
}
12331236
} else {
12341237
int proget;
12351238
if ((rc = proget_compute(projectile, ENT_PID_ELECTRON, process,
12361239
&proget)) != ENT_RETURN_SUCCESS)
12371240
ENT_RETURN(rc);
12381241
if (proget < PROGET_N - 1) {
1239-
*cross_section = Z * cross_section_compute(mode, proget,
1240-
cs0, cs1, p1, p2);
1242+
*cross_section = Z *
1243+
cross_section_compute(
1244+
mode, proget, cs0, cs1, p1, p2);
12411245
} else {
12421246
*cross_section = Z *
1243-
cross_section_compute(mode, PROGET_N - 2, cs0, cs1,
1244-
p1, p2) *
1247+
cross_section_compute(
1248+
mode, PROGET_N - 2, cs0, cs1, p1, p2) *
12451249
(ENT_WIDTH_W / ENT_WIDTH_W_TO_MUON - 1.);
12461250
}
12471251
}
@@ -1413,11 +1417,12 @@ static enum ent_return transport_step(struct ent_context * context,
14131417
if (fabs(drho) < 1E-03) {
14141418
ds = 2. * dX / (density1 + *density) - *step;
14151419
} else {
1416-
ds = *step * ((sqrt(density1 * density1 +
1417-
2. * dX * drho / *step) -
1418-
density1) /
1419-
drho -
1420-
1.);
1420+
ds = *step *
1421+
((sqrt(density1 * density1 +
1422+
2. * dX * drho / *step) -
1423+
density1) /
1424+
drho -
1425+
1.);
14211426
}
14221427
state->grammage = grammage_max;
14231428
state->position[0] += ds * sgn * state->direction[0];
@@ -1445,6 +1450,8 @@ static enum ent_return transport_step(struct ent_context * context,
14451450
if ((step2 > 0.) && ((step1 <= 0.) || (step2 < step1))) step1 = step2;
14461451
*step = (step1 < 0) ? 0. : step1;
14471452
*density = density1;
1453+
if (*medium == NULL)
1454+
*event = ENT_EVENT_EXIT;
14481455
return ENT_RETURN_SUCCESS;
14491456

14501457
#undef STEP_MIN
@@ -2015,26 +2022,26 @@ static void transport_rotate(
20152022
const double a2 = fabs(direction[2]);
20162023
if (a0 > a1) {
20172024
if (a0 > a2) {
2018-
const double nrm =
2019-
1. / sqrt(direction[0] * direction[0] +
2020-
direction[2] * direction[2]);
2025+
const double nrm = 1. /
2026+
sqrt(direction[0] * direction[0] +
2027+
direction[2] * direction[2]);
20212028
u0x = -direction[2] * nrm, u0z = direction[0] * nrm;
20222029
} else {
2023-
const double nrm =
2024-
1. / sqrt(direction[1] * direction[1] +
2025-
direction[2] * direction[2]);
2030+
const double nrm = 1. /
2031+
sqrt(direction[1] * direction[1] +
2032+
direction[2] * direction[2]);
20262033
u0y = direction[2] * nrm, u0z = -direction[1] * nrm;
20272034
}
20282035
} else {
20292036
if (a1 > a2) {
2030-
const double nrm =
2031-
1. / sqrt(direction[0] * direction[0] +
2032-
direction[1] * direction[1]);
2037+
const double nrm = 1. /
2038+
sqrt(direction[0] * direction[0] +
2039+
direction[1] * direction[1]);
20332040
u0x = direction[1] * nrm, u0y = -direction[0] * nrm;
20342041
} else {
2035-
const double nrm =
2036-
1. / sqrt(direction[1] * direction[1] +
2037-
direction[2] * direction[2]);
2042+
const double nrm = 1. /
2043+
sqrt(direction[1] * direction[1] +
2044+
direction[2] * direction[2]);
20382045
u0y = direction[2] * nrm, u0z = -direction[1] * nrm;
20392046
}
20402047
}
@@ -2453,8 +2460,9 @@ static void ancestor_likeliness_fill(struct ent_context * context,
24532460
proget_v[*np] = proget[i];
24542461
const double p0 = (*np > 0) ? p[*np - 1] : 0.;
24552462
p[*np] = p0 +
2456-
rho0 * Z * cross_section_compute(
2457-
mode, proget[i], cs0, cs1, p1, p2);
2463+
rho0 * Z *
2464+
cross_section_compute(
2465+
mode, proget[i], cs0, cs1, p1, p2);
24582466
ancestor_v[(*np)++] = pid[i];
24592467
}
24602468
}
@@ -2584,11 +2592,13 @@ static enum ent_return transport_ancestor_draw(struct ent_physics * physics,
25842592
PROGET_NC_NU_NEUTRON :
25852593
PROGET_NC_NU_BAR_NEUTRON;
25862594
proget_v[1] = proget_v[0] + 4;
2587-
p[0] = rho0 * N * cross_section_compute(mode,
2588-
proget_v[0], cs0, cs1, p1, p2);
2595+
p[0] = rho0 * N *
2596+
cross_section_compute(
2597+
mode, proget_v[0], cs0, cs1, p1, p2);
25892598
p[1] = p[0] +
2590-
rho0 * Z * cross_section_compute(
2591-
mode, proget_v[1], cs0, cs1, p1, p2);
2599+
rho0 * Z *
2600+
cross_section_compute(
2601+
mode, proget_v[1], cs0, cs1, p1, p2);
25922602

25932603
/* Elastic event on an atomic electron. */
25942604
if (apid == ENT_PID_NU_E)
@@ -2599,8 +2609,9 @@ static enum ent_return transport_ancestor_draw(struct ent_physics * physics,
25992609
proget_v[2] = PROGET_ELASTIC_NU_TAU;
26002610
if (daughter->pid < 0) proget_v[2]++;
26012611
p[2] = p[1] +
2602-
rho0 * Z * cross_section_compute(
2603-
mode, proget_v[2], cs0, cs1, p1, p2);
2612+
rho0 * Z *
2613+
cross_section_compute(
2614+
mode, proget_v[2], cs0, cs1, p1, p2);
26042615

26052616
ancestor_v[np++] = daughter->pid;
26062617
ancestor_v[np++] = daughter->pid;
@@ -2651,9 +2662,9 @@ static enum ent_return transport_ancestor_draw(struct ent_physics * physics,
26512662
proget_v[np] = PROGET_INVERSE_NU_MU_MU;
26522663
const double p0 = (np > 0) ? p[np - 1] : 0.;
26532664
p[np] = p0 +
2654-
rho0 * Z * cross_section_compute(mode,
2655-
proget_v[np], cs0, cs1, p1,
2656-
p2);
2665+
rho0 * Z *
2666+
cross_section_compute(mode,
2667+
proget_v[np], cs0, cs1, p1, p2);
26572668
ancestor_v[np++] = ENT_PID_NU_MU;
26582669
}
26592670
const double rho1 = context->ancestor(
@@ -2662,9 +2673,9 @@ static enum ent_return transport_ancestor_draw(struct ent_physics * physics,
26622673
proget_v[np] = PROGET_INVERSE_NU_TAU_TAU;
26632674
const double p0 = (np > 0) ? p[np - 1] : 0.;
26642675
p[np] = p0 +
2665-
rho1 * Z * cross_section_compute(mode,
2666-
proget_v[np], cs0, cs1, p1,
2667-
p2);
2676+
rho1 * Z *
2677+
cross_section_compute(mode,
2678+
proget_v[np], cs0, cs1, p1, p2);
26682679
ancestor_v[np++] = ENT_PID_NU_TAU;
26692680
}
26702681
} else if ((daughter->pid == ENT_PID_NU_BAR_MU) ||
@@ -2678,9 +2689,9 @@ static enum ent_return transport_ancestor_draw(struct ent_physics * physics,
26782689
PROGET_INVERSE_NU_BAR_E_TAU;
26792690
const double p0 = (np > 0) ? p[np - 1] : 0.;
26802691
p[np] = p0 +
2681-
rho0 * Z * cross_section_compute(mode,
2682-
proget_v[np], cs0, cs1, p1,
2683-
p2);
2692+
rho0 * Z *
2693+
cross_section_compute(mode,
2694+
proget_v[np], cs0, cs1, p1, p2);
26842695
ancestor_v[np++] = ENT_PID_NU_BAR_E;
26852696
}
26862697
}
@@ -2694,11 +2705,13 @@ static enum ent_return transport_ancestor_draw(struct ent_physics * physics,
26942705
PROGET_CC_NU_NEUTRON :
26952706
PROGET_CC_NU_BAR_NEUTRON;
26962707
proget_v[1] = proget_v[0] + 4;
2697-
p[0] = rho0 * N * cross_section_compute(mode,
2698-
proget_v[0], cs0, cs1, p1, p2);
2708+
p[0] = rho0 * N *
2709+
cross_section_compute(
2710+
mode, proget_v[0], cs0, cs1, p1, p2);
26992711
p[1] = p[0] +
2700-
rho0 * Z * cross_section_compute(
2701-
mode, proget_v[1], cs0, cs1, p1, p2);
2712+
rho0 * Z *
2713+
cross_section_compute(
2714+
mode, proget_v[1], cs0, cs1, p1, p2);
27022715

27032716
ancestor_v[np++] = npid;
27042717
ancestor_v[np++] = npid;
@@ -2752,10 +2765,9 @@ static enum ent_return vertex_dis_compute(struct ent_physics * physics,
27522765
/* Compute the relevant cross-sections and randomise the
27532766
* target accordingly.
27542767
*/
2755-
*cs_p = (medium->Z > 0.) ?
2756-
medium->Z *
2768+
*cs_p = (medium->Z > 0.) ? medium->Z *
27572769
cross_section_compute(mode, proget_p, cs0, cs1, p1, p2) :
2758-
0.;
2770+
0.;
27592771
const double N = medium->A - medium->Z;
27602772
*cs_n = (N > 0.) ?
27612773
N * cross_section_compute(mode, proget_n, cs0, cs1, p1, p2) :
@@ -3096,6 +3108,7 @@ enum ent_return ent_transport(struct ent_physics * physics,
30963108
if (step)
30973109
for (;;) {
30983110
/* Step until an event occurs. */
3111+
struct ent_medium * m = medium;
30993112
if ((rc = transport_step(context, state, &medium, &step,
31003113
&density, Xlim, &event_)) !=
31013114
ENT_RETURN_SUCCESS)
@@ -3105,6 +3118,16 @@ enum ent_return ent_transport(struct ent_physics * physics,
31053118
rc = ENT_RETURN_DOMAIN_ERROR;
31063119
goto exit;
31073120
}
3121+
3122+
/* Call any custom stepping action if a medium
3123+
* change occured.
3124+
*/
3125+
if ((context->stepping_action != NULL) &&
3126+
(medium != m)) {
3127+
if ((rc = context->stepping_action(context,
3128+
medium, state)) != ENT_RETURN_SUCCESS)
3129+
goto exit;
3130+
}
31083131
}
31093132
else {
31103133
/* This is a uniform medium of infinite extension.
@@ -3113,6 +3136,13 @@ enum ent_return ent_transport(struct ent_physics * physics,
31133136
event_ = transport_straight(context, state, density, Xlim);
31143137
}
31153138

3139+
/* Call any custom stepping action. */
3140+
if (context->stepping_action != NULL) {
3141+
if ((rc = context->stepping_action(context, medium, state)) !=
3142+
ENT_RETURN_SUCCESS)
3143+
goto exit;
3144+
}
3145+
31163146
/* Process any interaction. */
31173147
if ((event_ == ENT_EVENT_LIMIT_GRAMMAGE) &&
31183148
(foreseen == ENT_EVENT_NONE)) {

0 commit comments

Comments
 (0)