Skip to content

Add comp-formatter to McStas monitors #2315

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
willend merged 1 commit into from
Feb 19, 2026
Merged

Add comp-formatter to McStas monitors #2315

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
215 changes: 95 additions & 120 deletions mcstas-comps/monitors/Brilliance_monitor.comp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,172 +110,147 @@ INITIALIZE
%{
prsec = 1e-6;

BRIL_N = create_darr2d(nt, nlam);
BRIL_p = create_darr2d(nt, nlam);
BRIL_p2 = create_darr2d(nt, nlam);
BRIL_mean = create_darr1d(nlam);
BRIL_meanN = create_darr1d(nlam);
BRIL_meanE = create_darr1d(nlam);
BRIL_peak = create_darr1d(nlam);
BRIL_peakN = create_darr1d(nlam);
BRIL_peakE = create_darr1d(nlam);
BRIL_N = create_darr2d (nt, nlam);
BRIL_p = create_darr2d (nt, nlam);
BRIL_p2 = create_darr2d (nt, nlam);
BRIL_mean = create_darr1d (nlam);
BRIL_meanN = create_darr1d (nlam);
BRIL_meanE = create_darr1d (nlam);
BRIL_peak = create_darr1d (nlam);
BRIL_peakN = create_darr1d (nlam);
BRIL_peakE = create_darr1d (nlam);

BRIL_shape = create_darr1d(nt);
BRIL_shapeN = create_darr1d(nt);
BRIL_shapeE = create_darr1d(nt);
BRIL_shape = create_darr1d (nt);
BRIL_shapeN = create_darr1d (nt);
BRIL_shapeE = create_darr1d (nt);

tt_0 = t_0 * prsec;
tt_1 = t_1 * prsec;
dt = (t_1 - t_0) * prsec / nt;
dlam = (lambda_1 - lambda_0) / (nlam - 1);

tt_0 = t_0*prsec;
tt_1 = t_1*prsec;
dt = (t_1-t_0)*prsec/nt;
dlam = (lambda_1-lambda_0)/(nlam-1);

xmax = xwidth/2.0;
ymax = yheight/2.0;
xmax = xwidth / 2.0;
ymax = yheight / 2.0;
xmin = -xmax;
ymin = -ymax;

ster = xwidth * yheight/(source_dist*source_dist);
ster = xwidth * yheight / (source_dist * source_dist);

// Use instance name for monitor output if no input was given
if (!strcmp(filename,"\0")) sprintf(filename,"%s",NAME_CURRENT_COMP);
if (!strcmp (filename, "\0"))
sprintf (filename, "%s", NAME_CURRENT_COMP);
%}

TRACE
%{
int i,j;
int i, j;
double div;
double lambda;
double Pnorm;

PROP_Z0;
lambda = (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz);
if (x>xmin && x<xmax && y>ymin && y<ymax &&
lambda > lambda_0 && lambda < lambda_1) {
lambda = (2 * PI / V2K) / sqrt (vx * vx + vy * vy + vz * vz);
if (x > xmin && x < xmax && y > ymin && y < ymax && lambda > lambda_0 && lambda < lambda_1) {
if (t < tt_1 && t > tt_0) {
i = floor((lambda - lambda_0)*nlam/(lambda_1 - lambda_0));
j = floor((t-tt_0)*nt/(tt_1-tt_0));
Pnorm=p/dlam/ster/srcarea;
double Pnorm2 = Pnorm*Pnorm;
i = floor ((lambda - lambda_0) * nlam / (lambda_1 - lambda_0));
j = floor ((t - tt_0) * nt / (tt_1 - tt_0));
Pnorm = p / dlam / ster / srcarea;
double Pnorm2 = Pnorm * Pnorm;
#pragma acc atomic
BRIL_meanN[i] = BRIL_meanN[i] + 1;
#pragma acc atomic
BRIL_meanN[i] = BRIL_meanN[i]+1;
BRIL_mean[i] = BRIL_mean[i] + Pnorm;
#pragma acc atomic
BRIL_meanE[i] = BRIL_meanE[i] + Pnorm2;

Pnorm = Pnorm / Freq / dt;
Pnorm2 = Pnorm + Pnorm;
#pragma acc atomic
BRIL_mean[i] = BRIL_mean[i]+Pnorm;
#pragma acc atomic
BRIL_meanE[i] = BRIL_meanE[i]+Pnorm2;

Pnorm=Pnorm/Freq/dt;
Pnorm2=Pnorm+Pnorm;
#pragma acc atomic
BRIL_N[j][i] = BRIL_N[j][i]+1;
BRIL_N[j][i] = BRIL_N[j][i] + 1;
#pragma acc atomic
BRIL_p[j][i] = BRIL_p[j][i]+Pnorm;
BRIL_p[j][i] = BRIL_p[j][i] + Pnorm;
#pragma acc atomic
BRIL_p2[j][i] = BRIL_p2[j][i]+Pnorm2;
BRIL_p2[j][i] = BRIL_p2[j][i] + Pnorm2;
}
}
if (restore_neutron) {
RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
}
%}

SAVE
%{
if (!nowritefile) {
/* First, dump the 2D monitor */
if (!nowritefile) {
/* First, dump the 2D monitor */

/* For each Wavelength channel, find peak brilliance */
int i,j,jmax;
double Pnorm;
char ff[256];
char tt[256];
/* For each Wavelength channel, find peak brilliance */
int i, j, jmax;
double Pnorm;
char ff[256];
char tt[256];

for (i=0; i<nlam; i++) {
Pnorm = -1;
jmax = -1;
for (j=0; j<nt; j++) {
if (BRIL_p[j][i]>=Pnorm) {
Pnorm = BRIL_p[j][i];
jmax=j;
}
BRIL_shape[j] = BRIL_p[j][i];
BRIL_shapeN[j] = BRIL_N[j][i];
BRIL_shapeE[j] = BRIL_p2[j][i];
}
if (tofcuts == 1) {
sprintf(ff, "Shape_%s_%g",filename,lambda_0+i*dlam);
sprintf(tt, "Peak shape at %g AA",lambda_0+i*dlam);
DETECTOR_OUT_1D(
tt,
"TOF [us]",
"Peak Brilliance",
"Shape", t_0, t_1, nt,
&BRIL_shapeN[0],&BRIL_shape[0],&BRIL_shapeE[0],
ff);
for (i = 0; i < nlam; i++) {
Pnorm = -1;
jmax = -1;
for (j = 0; j < nt; j++) {
if (BRIL_p[j][i] >= Pnorm) {
Pnorm = BRIL_p[j][i];
jmax = j;
}
BRIL_shape[j] = BRIL_p[j][i];
BRIL_shapeN[j] = BRIL_N[j][i];
BRIL_shapeE[j] = BRIL_p2[j][i];
}
if (tofcuts == 1) {
sprintf (ff, "Shape_%s_%g", filename, lambda_0 + i * dlam);
sprintf (tt, "Peak shape at %g AA", lambda_0 + i * dlam);
DETECTOR_OUT_1D (tt, "TOF [us]", "Peak Brilliance", "Shape", t_0, t_1, nt, &BRIL_shapeN[0], &BRIL_shape[0], &BRIL_shapeE[0], ff);
}
BRIL_peakN[i] = BRIL_N[jmax][i];
BRIL_peak[i] = BRIL_p[jmax][i];
BRIL_peakE[i] = BRIL_p2[jmax][i];
}
BRIL_peakN[i] = BRIL_N[jmax][i];
BRIL_peak[i] = BRIL_p[jmax][i];
BRIL_peakE[i] = BRIL_p2[jmax][i];
}

sprintf(ff, "Mean_%s",filename);
DETECTOR_OUT_1D(
"Mean brilliance",
"Wavelength [AA]",
"Mean Brilliance",
"Mean", lambda_0, lambda_1, nlam,
&BRIL_meanN[0],&BRIL_mean[0],&BRIL_meanE[0],
ff);
sprintf(ff, "Peak_%s",filename);
DETECTOR_OUT_1D(
"Peak brilliance",
"Wavelength [AA]",
"Peak Brilliance",
"Peak", lambda_0, lambda_1, nlam,
&BRIL_peakN[0],&BRIL_peak[0],&BRIL_peakE[0],
ff);
sprintf (ff, "Mean_%s", filename);
DETECTOR_OUT_1D ("Mean brilliance", "Wavelength [AA]", "Mean Brilliance", "Mean", lambda_0, lambda_1, nlam, &BRIL_meanN[0], &BRIL_mean[0], &BRIL_meanE[0],
ff);
sprintf (ff, "Peak_%s", filename);
DETECTOR_OUT_1D ("Peak brilliance", "Wavelength [AA]", "Peak Brilliance", "Peak", lambda_0, lambda_1, nlam, &BRIL_peakN[0], &BRIL_peak[0], &BRIL_peakE[0],
ff);

/* MPI related NOTE: Order is important here! The 2D-data used to generate wavelength-slices and calculate
the peak brilliance should be done LAST, otherwise we will get a factor of MPI_node_count too much as
scatter/gather has been performed on the arrays... */
if (toflambda == 1) {
sprintf(ff, "TOFL_%s",filename);
DETECTOR_OUT_2D(
"TOF-wavelength brilliance",
"Time-of-flight [\\gms]", "Wavelength [AA]",
t_0, t_1, lambda_0, lambda_1,
nt, nlam,
&BRIL_N[0][0],&BRIL_p[0][0],&BRIL_p2[0][0],
filename);
/* MPI related NOTE: Order is important here! The 2D-data used to generate wavelength-slices and calculate
the peak brilliance should be done LAST, otherwise we will get a factor of MPI_node_count too much as
scatter/gather has been performed on the arrays... */
if (toflambda == 1) {
sprintf (ff, "TOFL_%s", filename);
DETECTOR_OUT_2D ("TOF-wavelength brilliance", "Time-of-flight [\\gms]", "Wavelength [AA]", t_0, t_1, lambda_0, lambda_1, nt, nlam, &BRIL_N[0][0],
&BRIL_p[0][0], &BRIL_p2[0][0], filename);
}
}
}
%}

FINALLY
%{
destroy_darr2d(BRIL_N);
destroy_darr2d(BRIL_p);
destroy_darr2d(BRIL_p2);
destroy_darr2d (BRIL_N);
destroy_darr2d (BRIL_p);
destroy_darr2d (BRIL_p2);

destroy_darr1d(BRIL_mean);
destroy_darr1d(BRIL_meanN);
destroy_darr1d(BRIL_meanE);
destroy_darr1d(BRIL_peak);
destroy_darr1d(BRIL_peakN);
destroy_darr1d(BRIL_peakE);
destroy_darr1d (BRIL_mean);
destroy_darr1d (BRIL_meanN);
destroy_darr1d (BRIL_meanE);
destroy_darr1d (BRIL_peak);
destroy_darr1d (BRIL_peakN);
destroy_darr1d (BRIL_peakE);

destroy_darr1d(BRIL_shape);
destroy_darr1d(BRIL_shapeN);
destroy_darr1d(BRIL_shapeE);
destroy_darr1d (BRIL_shape);
destroy_darr1d (BRIL_shapeN);
destroy_darr1d (BRIL_shapeE);
%}

MCDISPLAY
%{
multiline(5, (double)xmin, (double)ymin, 0.0,
(double)xmax, (double)ymin, 0.0,
(double)xmax, (double)ymax, 0.0,
(double)xmin, (double)ymax, 0.0,
(double)xmin, (double)ymin, 0.0);
multiline (5, (double)xmin, (double)ymin, 0.0, (double)xmax, (double)ymin, 0.0, (double)xmax, (double)ymax, 0.0, (double)xmin, (double)ymax, 0.0, (double)xmin,
(double)ymin, 0.0);
%}

END
72 changes: 33 additions & 39 deletions mcstas-comps/monitors/Cyl_monitor.comp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,73 +60,67 @@ DECLARE

INITIALIZE
%{
PSD_N = create_darr1d(nr);
PSD_p = create_darr1d(nr);
PSD_p2 = create_darr1d(nr);
PSD_N = create_darr1d (nr);
PSD_p = create_darr1d (nr);
PSD_p2 = create_darr1d (nr);

// Use instance name for monitor output if no input was given
if (!strcmp(filename,"\0")) sprintf(filename,"%s",NAME_CURRENT_COMP);
if (!strcmp (filename, "\0"))
sprintf (filename, "%s", NAME_CURRENT_COMP);
%}

TRACE
%{
int i, j;
double t0, t1, phi;

if(cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight) == 1) {
if(t0<0) {
if(t1>0) {
PROP_DT(t1);
/* Calculate pixel */
if (fabs(y) <= yheight/2.0) {
phi=atan2(x,z)*RAD2DEG;

if (phi >= thmin && phi <= thmax) {
i=floor((nr) * (phi-thmin)/(thmax-thmin));

double p2 = p*p;
if (cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius, yheight) == 1) {
if (t0 < 0) {
if (t1 > 0) {
PROP_DT (t1);
/* Calculate pixel */
if (fabs (y) <= yheight / 2.0) {
phi = atan2 (x, z) * RAD2DEG;

if (phi >= thmin && phi <= thmax) {
i = floor ((nr) * (phi - thmin) / (thmax - thmin));

double p2 = p * p;
#pragma acc atomic
PSD_N[i] = PSD_N[i] + 1;

#pragma acc atomic
PSD_N[i] = PSD_N[i]+1;
PSD_p[i] = PSD_p[i] + p;

#pragma acc atomic
PSD_p[i] = PSD_p[i]+p;

#pragma acc atomic
PSD_p2[i] = PSD_p2[i]+p2;
}
}
PSD_p2[i] = PSD_p2[i] + p2;
}
}
}
}
}
if (restore_neutron) {
RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
}
%}

SAVE
%{
if (!nowritefile) {
DETECTOR_OUT_1D(
"Cylindrical monitor",
"radial position [deg]",
"Intensity",
"Theta",
thmin, thmax, nr,
&PSD_N[0],&PSD_p[0],&PSD_p2[0],
filename);
}
if (!nowritefile) {
DETECTOR_OUT_1D ("Cylindrical monitor", "radial position [deg]", "Intensity", "Theta", thmin, thmax, nr, &PSD_N[0], &PSD_p[0], &PSD_p2[0], filename);
}
%}

FINALLY
%{
destroy_darr1d(PSD_N);
destroy_darr1d(PSD_p);
destroy_darr1d(PSD_p2);
destroy_darr1d (PSD_N);
destroy_darr1d (PSD_p);
destroy_darr1d (PSD_p2);
%}

MCDISPLAY
%{
circle("xz", 0, 0, 0, radius );
circle ("xz", 0, 0, 0, radius);
%}

END
Loading