diff --git a/mcstas-comps/monitors/Brilliance_monitor.comp b/mcstas-comps/monitors/Brilliance_monitor.comp index c7c7a9104..21417e226 100644 --- a/mcstas-comps/monitors/Brilliance_monitor.comp +++ b/mcstas-comps/monitors/Brilliance_monitor.comp @@ -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 && xymin && y 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=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 diff --git a/mcstas-comps/monitors/Cyl_monitor.comp b/mcstas-comps/monitors/Cyl_monitor.comp index aacc8d837..f0ddda25f 100644 --- a/mcstas-comps/monitors/Cyl_monitor.comp +++ b/mcstas-comps/monitors/Cyl_monitor.comp @@ -60,12 +60,13 @@ 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 @@ -73,60 +74,53 @@ 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 diff --git a/mcstas-comps/monitors/Cyl_monitor_PSD.comp b/mcstas-comps/monitors/Cyl_monitor_PSD.comp index 2300f6b7c..8c825912b 100644 --- a/mcstas-comps/monitors/Cyl_monitor_PSD.comp +++ b/mcstas-comps/monitors/Cyl_monitor_PSD.comp @@ -60,12 +60,13 @@ DECLARE INITIALIZE %{ - PSD_N = create_darr2d(nr,ny); - PSD_p = create_darr2d(nr,ny); - PSD_p2 = create_darr2d(nr,ny); + PSD_N = create_darr2d (nr, ny); + PSD_p = create_darr2d (nr, ny); + PSD_p2 = create_darr2d (nr, ny); // 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 @@ -73,60 +74,55 @@ 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)); - j=floor((ny) * (y-yheight/2.0)/(yheight)); - 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)); + j = floor ((ny) * (y - yheight / 2.0) / (yheight)); + double p2 = p * p; + #pragma acc atomic - PSD_N[i][j] = PSD_N[i][j]+1; + PSD_N[i][j] = PSD_N[i][j] + 1; #pragma acc atomic - PSD_p[i][j] = PSD_p[i][j]+p; - + PSD_p[i][j] = PSD_p[i][j] + p; + #pragma acc atomic - PSD_p2[i][j] = PSD_p2[i][j]+p2; - } - } + PSD_p2[i][j] = PSD_p2[i][j] + 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_2D( - "Cylindrical PSD monitor", - "radial position [deg]", - "y [m]", - thmin, thmax, -yheight/2.0, yheight/2.0, nr, ny, - &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("Cylindrical PSD monitor", "radial position [deg]", "y [m]", thmin, thmax, -yheight / 2.0, yheight / 2.0, nr, ny, &PSD_N[0][0], &PSD_p[0][0], + &PSD_p2[0][0], filename); + } %} FINALLY %{ - destroy_darr2d(PSD_N); - destroy_darr2d(PSD_p); - destroy_darr2d(PSD_p2); + destroy_darr2d (PSD_N); + destroy_darr2d (PSD_p); + destroy_darr2d (PSD_p2); %} MCDISPLAY %{ - circle("xz", 0, 0, 0, radius ); + circle ("xz", 0, 0, 0, radius); %} END diff --git a/mcstas-comps/monitors/Cyl_monitor_TOF.comp b/mcstas-comps/monitors/Cyl_monitor_TOF.comp index 6d5f0f8d5..e347e9a7b 100644 --- a/mcstas-comps/monitors/Cyl_monitor_TOF.comp +++ b/mcstas-comps/monitors/Cyl_monitor_TOF.comp @@ -63,12 +63,13 @@ DECLARE INITIALIZE %{ - PSD_N = create_darr2d(nr,nt); - PSD_p = create_darr2d(nr,nt); - PSD_p2 = create_darr2d(nr,nt); + PSD_N = create_darr2d (nr, nt); + PSD_p = create_darr2d (nr, nt); + PSD_p2 = create_darr2d (nr, nt); // 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 @@ -76,62 +77,57 @@ 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) { - if (t >= tmin && t < tmax) { - i=floor((nr) * (phi-thmin)/(thmax-thmin)); - j=floor((nt) * (t-tmin)/(tmax-tmin)); - 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) { + if (t >= tmin && t < tmax) { + i = floor ((nr) * (phi - thmin) / (thmax - thmin)); + j = floor ((nt) * (t - tmin) / (tmax - tmin)); + double p2 = p * p; + #pragma acc atomic - PSD_N[i][j] = PSD_N[i][j]+1; - + PSD_N[i][j] = PSD_N[i][j] + 1; + #pragma acc atomic - PSD_p[i][j] = PSD_p[i][j]+p; - + PSD_p[i][j] = PSD_p[i][j] + p; + #pragma acc atomic - PSD_p2[i][j] = PSD_p2[i][j]+p2; - } - } - } + PSD_p2[i][j] = PSD_p2[i][j] + 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_2D( - "Cylindrical TOF monitor", - "radial position [deg]", - "TOF [s]", - thmin, thmax, tmin, tmax, nr, nt, - &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("Cylindrical TOF monitor", "radial position [deg]", "TOF [s]", thmin, thmax, tmin, tmax, nr, nt, &PSD_N[0][0], &PSD_p[0][0], &PSD_p2[0][0], + filename); + } %} FINALLY %{ - destroy_darr2d(PSD_N); - destroy_darr2d(PSD_p); - destroy_darr2d(PSD_p2); + destroy_darr2d (PSD_N); + destroy_darr2d (PSD_p); + destroy_darr2d (PSD_p2); %} MCDISPLAY %{ - circle("xz", 0, 0, 0, radius ); + circle ("xz", 0, 0, 0, radius); %} END diff --git a/mcstas-comps/monitors/Div1D_monitor.comp b/mcstas-comps/monitors/Div1D_monitor.comp index 8b60b1227..9178da23e 100644 --- a/mcstas-comps/monitors/Div1D_monitor.comp +++ b/mcstas-comps/monitors/Div1D_monitor.comp @@ -67,27 +67,34 @@ INITIALIZE %{ int i; - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("Error(%s): Div1D_monitor: Null detection area !\n" - "(xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(-1); + printf ("Error(%s): Div1D_monitor: Null detection area !\n" + "(xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (-1); } - Div_N = create_darr1d(ndiv); - Div_p = create_darr1d(ndiv); - Div_p2 = create_darr1d(ndiv); - for (i=0; ixmin && xymin && y xmin && x < xmax && y > ymin && y < ymax) { + if (!vertical) { + div = RAD2DEG * atan2 (vx, vz); + } else { + div = RAD2DEG * atan2 (vy, vz); } - if (div < (double)maxdiv && div > -(double)maxdiv){ - i = floor((div + (double)maxdiv)*ndiv/(2.0*(double)maxdiv)); + if (div < (double)maxdiv && div > -(double)maxdiv) { + i = floor ((div + (double)maxdiv) * ndiv / (2.0 * (double)maxdiv)); - double p2 = p*p; + double p2 = p * p; #pragma acc atomic Div_N[i] = Div_N[i] + 1; @@ -120,41 +126,34 @@ TRACE } } 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){ - if (!vertical){ - DETECTOR_OUT_1D( - "Horizontal divergence monitor","horizontal divergence [deg]", - "Intensity","divergence", -maxdiv, maxdiv, ndiv, - &Div_N[0],&Div_p[0],&Div_p2[0],filename); - }else{ - DETECTOR_OUT_1D( - "Vertical divergence monitor","vertical divergence [deg]", - "Intensity","divergence", -maxdiv, maxdiv, ndiv, - &Div_N[0],&Div_p[0],&Div_p2[0],filename); + if (!nowritefile) { + if (!vertical) { + DETECTOR_OUT_1D ("Horizontal divergence monitor", "horizontal divergence [deg]", "Intensity", "divergence", -maxdiv, maxdiv, ndiv, &Div_N[0], &Div_p[0], + &Div_p2[0], filename); + } else { + DETECTOR_OUT_1D ("Vertical divergence monitor", "vertical divergence [deg]", "Intensity", "divergence", -maxdiv, maxdiv, ndiv, &Div_N[0], &Div_p[0], + &Div_p2[0], filename); } } %} FINALLY %{ - destroy_darr1d(Div_N); - destroy_darr1d(Div_p); - destroy_darr1d(Div_p2); + destroy_darr1d (Div_N); + destroy_darr1d (Div_p); + destroy_darr1d (Div_p2); %} 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); %} diff --git a/mcstas-comps/monitors/DivLambda_monitor.comp b/mcstas-comps/monitors/DivLambda_monitor.comp index ca0386cc2..d0f6c4ad1 100644 --- a/mcstas-comps/monitors/DivLambda_monitor.comp +++ b/mcstas-comps/monitors/DivLambda_monitor.comp @@ -68,93 +68,88 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("Divlambda_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("Divlambda_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - Div_N = create_darr2d(nL, nh); - Div_p = create_darr2d(nL, nh); - Div_p2 = create_darr2d(nL, nh); + Div_N = create_darr2d (nL, nh); + Div_p = create_darr2d (nL, nh); + Div_p2 = create_darr2d (nL, nh); - NORM(nx,ny,nz); + NORM (nx, ny, nz); // 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 v, vn; PROP_Z0; - lambda = (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz); - if (x>xmin && xymin && y Lmin && lambda < Lmax) - { + lambda = (2 * PI / V2K) / sqrt (vx * vx + vy * vy + vz * vz); + if (x > xmin && x < xmax && y > ymin && y < ymax && lambda > Lmin && lambda < Lmax) { /* Find length of projection onto the [nx ny nz] axis */ - vn = scalar_prod(vx, vy, vz, nx, ny, nz); - div = RAD2DEG*atan2(vx,vn); + vn = scalar_prod (vx, vy, vz, nx, ny, nz); + div = RAD2DEG * atan2 (vx, vn); - if (div < maxdiv_h && div > -maxdiv_h) - { - i = floor((lambda - Lmin)*nL/(Lmax - Lmin)); - j = floor((div + maxdiv_h)*nh/(2.0*maxdiv_h)); + if (div < maxdiv_h && div > -maxdiv_h) { + i = floor ((lambda - Lmin) * nL / (Lmax - Lmin)); + j = floor ((div + maxdiv_h) * nh / (2.0 * maxdiv_h)); - double p2 = p*p; + double p2 = p * p; #pragma acc atomic - Div_N[i][j] = Div_N[i][j]+1; + Div_N[i][j] = Div_N[i][j] + 1; #pragma acc atomic - Div_p[i][j] = Div_p[i][j]+p; + Div_p[i][j] = Div_p[i][j] + p; #pragma acc atomic - Div_p2[i][j] = Div_p2[i][j]+p2; + Div_p2[i][j] = Div_p2[i][j] + p2; SCATTER; } } 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_2D( - "Wavelength-divergence monitor", - "Wavelength [AA]", - "divergence [deg]", - Lmin, Lmax, -maxdiv_h, maxdiv_h, - nL, nh, - &Div_N[0][0],&Div_p[0][0],&Div_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("Wavelength-divergence monitor", "Wavelength [AA]", "divergence [deg]", Lmin, Lmax, -maxdiv_h, maxdiv_h, nL, nh, &Div_N[0][0], &Div_p[0][0], + &Div_p2[0][0], filename); + } %} FINALLY %{ - destroy_darr2d(Div_N); - destroy_darr2d(Div_p); - destroy_darr2d(Div_p2); + destroy_darr2d (Div_N); + destroy_darr2d (Div_p); + destroy_darr2d (Div_p2); %} 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 diff --git a/mcstas-comps/monitors/DivPos_monitor.comp b/mcstas-comps/monitors/DivPos_monitor.comp index a807a61c4..f2cd3df89 100644 --- a/mcstas-comps/monitors/DivPos_monitor.comp +++ b/mcstas-comps/monitors/DivPos_monitor.comp @@ -70,52 +70,57 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("DivPos_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("DivPos_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - Div_N = create_darr2d(nb, ndiv); - Div_p = create_darr2d(nb, ndiv); - Div_p2 = create_darr2d(nb, ndiv); + Div_N = create_darr2d (nb, ndiv); + Div_p = create_darr2d (nb, ndiv); + Div_p2 = create_darr2d (nb, ndiv); - NORM(nx,ny,nz); + NORM (nx, ny, nz); // 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 v, vn; PROP_Z0; - if (x>xmin && xymin && y xmin && x < xmax && y > ymin && y < ymax) { /* Find length of projection onto the [nx ny nz] axis */ - vn = scalar_prod(vx, vy, vz, nx, ny, nz); - if(!vertical){ - div = RAD2DEG*atan2(vx,vn); - }else{ - div = RAD2DEG*atan2(vy,vn); + vn = scalar_prod (vx, vy, vz, nx, ny, nz); + if (!vertical) { + div = RAD2DEG * atan2 (vx, vn); + } else { + div = RAD2DEG * atan2 (vy, vn); } - if (div < maxdiv && div > -maxdiv) - { - if(!vertical){ - i = floor((x - xmin)*nb/(xmax - xmin)); - }else{ - i = floor((y - ymin)*nb/(ymax - ymin)); + if (div < maxdiv && div > -maxdiv) { + if (!vertical) { + i = floor ((x - xmin) * nb / (xmax - xmin)); + } else { + i = floor ((y - ymin) * nb / (ymax - ymin)); } - j = floor((div + maxdiv)*ndiv/(2.0*maxdiv)); - double p2 = p*p; + j = floor ((div + maxdiv) * ndiv / (2.0 * maxdiv)); + double p2 = p * p; #pragma acc atomic Div_N[i][j] = Div_N[i][j] + 1; #pragma acc atomic @@ -126,36 +131,30 @@ TRACE } } 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_2D( - "Position-divergence monitor","pos [m]","divergence [deg]", - (!vertical?xmin:ymin), (!vertical?xmax:ymax), -maxdiv, maxdiv, nb, ndiv, - &Div_N[0][0],&Div_p[0][0],&Div_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("Position-divergence monitor", "pos [m]", "divergence [deg]", (!vertical ? xmin : ymin), (!vertical ? xmax : ymax), -maxdiv, maxdiv, nb, + ndiv, &Div_N[0][0], &Div_p[0][0], &Div_p2[0][0], filename); + } %} FINALLY %{ - destroy_darr2d(Div_N); - destroy_darr2d(Div_p); - destroy_darr2d(Div_p2); + destroy_darr2d (Div_N); + destroy_darr2d (Div_p); + destroy_darr2d (Div_p2); %} 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 diff --git a/mcstas-comps/monitors/Divergence_monitor.comp b/mcstas-comps/monitors/Divergence_monitor.comp index 732fdf089..06b68b031 100644 --- a/mcstas-comps/monitors/Divergence_monitor.comp +++ b/mcstas-comps/monitors/Divergence_monitor.comp @@ -70,45 +70,49 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("Divergence_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("Divergence_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - Div_N = create_darr2d(nh, nv); - Div_p = create_darr2d(nh, nv); - Div_p2 = create_darr2d(nh, nv); + Div_N = create_darr2d (nh, nv); + Div_p = create_darr2d (nh, nv); + Div_p2 = create_darr2d (nh, nv); - NORM(nx,ny,nz); + NORM (nx, ny, nz); // 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 h_div, v_div; double v, vn; PROP_Z0; - if (x>xmin && xymin && y xmin && x < xmax && y > ymin && y < ymax) { /* Find length of projection onto the [nx ny nz] axis */ - vn = scalar_prod(vx, vy, vz, nx, ny, nz); - h_div = RAD2DEG*atan2(vx,vn); - v_div = RAD2DEG*atan2(vy,vn); - if (h_div < maxdiv_h && h_div > -maxdiv_h && - v_div < maxdiv_v && v_div > -maxdiv_v) - { - i = floor((h_div + maxdiv_h)*nh/(2.0*maxdiv_h)); - j = floor((v_div + maxdiv_v)*nv/(2.0*maxdiv_v)); - double p2 = p*p; + vn = scalar_prod (vx, vy, vz, nx, ny, nz); + h_div = RAD2DEG * atan2 (vx, vn); + v_div = RAD2DEG * atan2 (vy, vn); + if (h_div < maxdiv_h && h_div > -maxdiv_h && v_div < maxdiv_v && v_div > -maxdiv_v) { + i = floor ((h_div + maxdiv_h) * nh / (2.0 * maxdiv_h)); + j = floor ((v_div + maxdiv_v) * nv / (2.0 * maxdiv_v)); + double p2 = p * p; #pragma acc atomic Div_N[i][j] = Div_N[i][j] + 1; #pragma acc atomic @@ -119,38 +123,29 @@ TRACE } } 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_2D( - "Divergence monitor", - "X divergence [deg]", - "Y divergence [deg]", - -maxdiv_h, maxdiv_h, -maxdiv_v, maxdiv_v, - nh, nv, - &Div_N[0][0],&Div_p[0][0],&Div_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("Divergence monitor", "X divergence [deg]", "Y divergence [deg]", -maxdiv_h, maxdiv_h, -maxdiv_v, maxdiv_v, nh, nv, &Div_N[0][0], + &Div_p[0][0], &Div_p2[0][0], filename); + } %} FINALLY %{ - destroy_darr2d(Div_N); - destroy_darr2d(Div_p); - destroy_darr2d(Div_p2); + destroy_darr2d (Div_N); + destroy_darr2d (Div_p); + destroy_darr2d (Div_p2); %} 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 diff --git a/mcstas-comps/monitors/EPSD_monitor.comp b/mcstas-comps/monitors/EPSD_monitor.comp index 892d34143..ea467809d 100644 --- a/mcstas-comps/monitors/EPSD_monitor.comp +++ b/mcstas-comps/monitors/EPSD_monitor.comp @@ -63,67 +63,70 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("EPSD_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("EPSD_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - PSD_N = create_darr2d(nx, nE); - PSD_p = create_darr2d(nx, nE); - PSD_p2 = create_darr2d(nx, nE); + PSD_N = create_darr2d (nx, nE); + PSD_p = create_darr2d (nx, nE); + PSD_p2 = create_darr2d (nx, nE); // 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 E,v; + int i, j; + double E, v; PROP_Z0; - v=vx*vx+vy*vy+vz*vz; - E=VS2E*(vx*vx+vy*vy+vz*vz); + v = vx * vx + vy * vy + vz * vz; + E = VS2E * (vx * vx + vy * vy + vz * vz); - if (x>xmin && xymin && yEmin && E xmin && x < xmax && y > ymin && y < ymax && E > Emin && E < Emax) { + i = floor ((x - xmin) * nx / (xmax - xmin)); + j = floor ((E - Emin) * nE / (Emax - Emin)); - double p2 = p*p; + double p2 = p * p; #pragma acc atomic - PSD_N[i][j] = PSD_N[i][j]+1; + PSD_N[i][j] = PSD_N[i][j] + 1; #pragma acc atomic - PSD_p[i][j] = PSD_p[i][j]+p; + PSD_p[i][j] = PSD_p[i][j] + p; #pragma acc atomic - PSD_p2[i][j] = PSD_p2[i][j]+p2; + PSD_p2[i][j] = PSD_p2[i][j] + 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_2D("EPSD monitor", "Position [m]", "Energy [meV]", - xmin, xmax, Emin, Emax, nx, nE, - &PSD_N[0][0], &PSD_p[0][0], &PSD_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("EPSD monitor", "Position [m]", "Energy [meV]", xmin, xmax, Emin, Emax, nx, nE, &PSD_N[0][0], &PSD_p[0][0], &PSD_p2[0][0], filename); + } %} FINALLY %{ - destroy_darr2d(PSD_N); - destroy_darr2d(PSD_p); - destroy_darr2d(PSD_p2); + destroy_darr2d (PSD_N); + destroy_darr2d (PSD_p); + destroy_darr2d (PSD_p2); %} MCDISPLAY diff --git a/mcstas-comps/monitors/E_monitor.comp b/mcstas-comps/monitors/E_monitor.comp index 2a301a109..687a00a56 100644 --- a/mcstas-comps/monitors/E_monitor.comp +++ b/mcstas-comps/monitors/E_monitor.comp @@ -67,24 +67,31 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("E_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("E_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - E_N = create_darr1d(nE); - E_p = create_darr1d(nE); - E_p2 = create_darr1d(nE); + E_N = create_darr1d (nE); + E_p = create_darr1d (nE); + E_p2 = create_darr1d (nE); S_p = S_pE = S_pE2 = 0; // 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 @@ -93,20 +100,18 @@ TRACE double E; PROP_Z0; - if (x>xmin && xymin && y xmin && x < xmax && y > ymin && y < ymax) { + E = VS2E * (vx * vx + vy * vy + vz * vz); S_p += p; - S_pE += p*E; - S_pE2 += p*E*E; + S_pE += p * E; + S_pE2 += p * E * E; - i = floor((E-Emin)*nE/(Emax-Emin)); - if(i >= 0 && i < nE) - { - double p2 = p*p; + i = floor ((E - Emin) * nE / (Emax - Emin)); + if (i >= 0 && i < nE) { + double p2 = p * p; #pragma acc atomic - E_N[i] = E_N[i] +1; + E_N[i] = E_N[i] + 1; #pragma acc atomic E_p[i] = E_p[i] + p; #pragma acc atomic @@ -115,40 +120,31 @@ TRACE } } 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( - "Energy monitor", - "Energy [meV]", - "Intensity", - "E", Emin, Emax, nE, - &E_N[0],&E_p[0],&E_p2[0], - filename); - if (S_p) printf(" : %g meV , E-width : %g meV \n", - S_pE/S_p,sqrt(S_pE2/S_p - S_pE*S_pE/(S_p*S_p)) ); -} + if (!nowritefile) { + DETECTOR_OUT_1D ("Energy monitor", "Energy [meV]", "Intensity", "E", Emin, Emax, nE, &E_N[0], &E_p[0], &E_p2[0], filename); + if (S_p) + printf (" : %g meV , E-width : %g meV \n", S_pE / S_p, sqrt (S_pE2 / S_p - S_pE * S_pE / (S_p * S_p))); + } %} FINALLY %{ - destroy_darr1d(E_N); - destroy_darr1d(E_p); - destroy_darr1d(E_p2); + destroy_darr1d (E_N); + destroy_darr1d (E_p); + destroy_darr1d (E_p2); %} 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 diff --git a/mcstas-comps/monitors/Event_monitor_simple.comp b/mcstas-comps/monitors/Event_monitor_simple.comp index 42cbad490..95a1a5f96 100644 --- a/mcstas-comps/monitors/Event_monitor_simple.comp +++ b/mcstas-comps/monitors/Event_monitor_simple.comp @@ -31,24 +31,24 @@ SETTING PARAMETERS (nevents=1e6) DECLARE %{ - DArray2d Events; + DArray2d Events; unsigned long Nevents; char fullfile[1024]; char outputdir[512]; %} INITIALIZE %{ - Nevents = ceil(nevents); - Events = create_darr2d(Nevents, 11); + Nevents = ceil (nevents); + Events = create_darr2d (Nevents, 11); if (dirname == NULL) { - sprintf(outputdir,"."); + sprintf (outputdir, "."); } else { - sprintf(outputdir,dirname); + sprintf (outputdir, dirname); } #ifndef USE_MPI - sprintf(fullfile,"%s/%s.%s",outputdir,NAME_CURRENT_COMP,"log"); + sprintf (fullfile, "%s/%s.%s", outputdir, NAME_CURRENT_COMP, "log"); #else - sprintf(fullfile,"%s/%s_%i.%s",outputdir,NAME_CURRENT_COMP,mpi_node_rank,"log"); + sprintf (fullfile, "%s/%s_%i.%s", outputdir, NAME_CURRENT_COMP, mpi_node_rank, "log"); #endif %} @@ -59,7 +59,7 @@ TRACE if (i < Nevents) { #pragma acc atomic write Events[i][0] = x; - #pragma acc atomic write + #pragma acc atomic write Events[i][1] = y; #pragma acc atomic write Events[i][2] = z; @@ -84,30 +84,30 @@ TRACE SAVE %{ - FILE * fp; - printf("Storing event logfile --> %s\n",fullfile); - fp = fopen(fullfile, "w+"); + FILE* fp; + printf ("Storing event logfile --> %s\n", fullfile); + fp = fopen (fullfile, "w+"); if (!fp) { - fprintf(stderr,"Could not open file %s in w+ mode\n",fullfile); - exit(-1); + fprintf (stderr, "Could not open file %s in w+ mode\n", fullfile); + exit (-1); } - fprintf(fp,"# Event-file %s, contains %lu events\n",fullfile,Nevents); - fprintf(fp,"########################################### START EVENT LIST ##############################################\n"); - fprintf(fp,"# id x y z vx vy vz t sx sy sz p # \n"); - fprintf(fp,"###########################################################################################################\n"); + fprintf (fp, "# Event-file %s, contains %lu events\n", fullfile, Nevents); + fprintf (fp, "########################################### START EVENT LIST ##############################################\n"); + fprintf (fp, "# id x y z vx vy vz t sx sy sz p # \n"); + fprintf (fp, "###########################################################################################################\n"); unsigned long i; int j; - for (i=0; i= 0 && i < nU) - { - double pp=particle_getvar(_particle,signal,&suc); - double p2 = pp*pp; - #pragma acc atomic - U_N[i] = U_N[i] +1; - #pragma acc atomic - U_p[i] = U_p[i] + pp; - #pragma acc atomic - U_p2[i] = U_p2[i] + p2; - SCATTER; - } - + int i = floor ((U - Umin) * nU / (Umax - Umin)); + if (!suc && i >= 0 && i < nU) { + double pp = particle_getvar (_particle, signal, &suc); + double p2 = pp * pp; + #pragma acc atomic + U_N[i] = U_N[i] + 1; + #pragma acc atomic + U_p[i] = U_p[i] + pp; + #pragma acc atomic + U_p2[i] = U_p2[i] + p2; + SCATTER; + } + 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( - "Flex monitor 1D", - username, - signalname, - username, Umin, Umax, nU, - &U_N[0],&U_p[0],&U_p2[0], - filename); + DETECTOR_OUT_1D ("Flex monitor 1D", username, signalname, username, Umin, Umax, nU, &U_N[0], &U_p[0], &U_p2[0], filename); } %} FINALLY %{ - destroy_darr1d(U_N); - destroy_darr1d(U_p); - destroy_darr1d(U_p2); + destroy_darr1d (U_N); + destroy_darr1d (U_p); + destroy_darr1d (U_p2); %} MCDISPLAY diff --git a/mcstas-comps/monitors/Flex_monitor_2D.comp b/mcstas-comps/monitors/Flex_monitor_2D.comp index 59acf041f..5aa51d689 100644 --- a/mcstas-comps/monitors/Flex_monitor_2D.comp +++ b/mcstas-comps/monitors/Flex_monitor_2D.comp @@ -65,82 +65,76 @@ DECLARE INITIALIZE %{ - U_N = create_darr2d(nU1,nU2); - U_p = create_darr2d(nU1,nU2); - U_p2 = create_darr2d(nU1,nU2); - - if(uid1!=-1){ - snprintf(username1,127,"uservar%d",uid1); - }else{ - snprintf(username1,127,"%s",ustring1); + U_N = create_darr2d (nU1, nU2); + U_p = create_darr2d (nU1, nU2); + U_p2 = create_darr2d (nU1, nU2); + + if (uid1 != -1) { + snprintf (username1, 127, "uservar%d", uid1); + } else { + snprintf (username1, 127, "%s", ustring1); } - if(uid2!=-1){ - snprintf(username2,127,"uservar%d",uid2); - }else{ - snprintf(username2,127,"%s",ustring2); + if (uid2 != -1) { + snprintf (username2, 127, "uservar%d", uid2); + } else { + snprintf (username2, 127, "%s", ustring2); } - snprintf(signalname,127,"%s",signal); + snprintf (signalname, 127, "%s", signal); // 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 %{ - double U1,U2; - int suc1,suc2,suc3; - - if(uid1!=-1){ - U1 = particle_getuservar_byid(_particle,uid1,&suc1); - }else{ - U1 = particle_getvar(_particle,ustring1,&suc1); + double U1, U2; + int suc1, suc2, suc3; + + if (uid1 != -1) { + U1 = particle_getuservar_byid (_particle, uid1, &suc1); + } else { + U1 = particle_getvar (_particle, ustring1, &suc1); } - if(uid2!=-1){ - U2 = particle_getuservar_byid(_particle,uid2,&suc2); - }else{ - U2 = particle_getvar(_particle,ustring2,&suc2); + if (uid2 != -1) { + U2 = particle_getuservar_byid (_particle, uid2, &suc2); + } else { + U2 = particle_getvar (_particle, ustring2, &suc2); } - - int i = floor((U1-Umin1)*nU1/(Umax1-Umin1)); - int j = floor((U2-Umin2)*nU2/(Umax2-Umin2)); - if(!suc1 && i >= 0 && i < nU1) - { - if(!suc2 && j >= 0 && j < nU2) - { - double pp=particle_getvar(_particle,signal,&suc3); - double p2 = pp*pp; - #pragma acc atomic - U_N[i][j] = U_N[i][j] +1; - #pragma acc atomic - U_p[i][j] = U_p[i][j] + pp; - #pragma acc atomic - U_p2[i][j] = U_p2[i][j] + p2; - SCATTER; - } + + int i = floor ((U1 - Umin1) * nU1 / (Umax1 - Umin1)); + int j = floor ((U2 - Umin2) * nU2 / (Umax2 - Umin2)); + if (!suc1 && i >= 0 && i < nU1) { + if (!suc2 && j >= 0 && j < nU2) { + double pp = particle_getvar (_particle, signal, &suc3); + double p2 = pp * pp; + #pragma acc atomic + U_N[i][j] = U_N[i][j] + 1; + #pragma acc atomic + U_p[i][j] = U_p[i][j] + pp; + #pragma acc atomic + U_p2[i][j] = U_p2[i][j] + p2; + SCATTER; } - + } + 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_2D( - "Flex monitor 2D", - username1, username2, - Umin1, Umax1, Umin2, Umax2, nU1, nU2, - &U_N[0][0],&U_p[0][0],&U_p2[0][0], - filename); + DETECTOR_OUT_2D ("Flex monitor 2D", username1, username2, Umin1, Umax1, Umin2, Umax2, nU1, nU2, &U_N[0][0], &U_p[0][0], &U_p2[0][0], filename); } %} FINALLY %{ - destroy_darr2d(U_N); - destroy_darr2d(U_p); - destroy_darr2d(U_p2); + destroy_darr2d (U_N); + destroy_darr2d (U_p); + destroy_darr2d (U_p2); %} MCDISPLAY diff --git a/mcstas-comps/monitors/Flex_monitor_3D.comp b/mcstas-comps/monitors/Flex_monitor_3D.comp index 9e39e938b..81b8101e0 100644 --- a/mcstas-comps/monitors/Flex_monitor_3D.comp +++ b/mcstas-comps/monitors/Flex_monitor_3D.comp @@ -77,75 +77,73 @@ DECLARE INITIALIZE %{ - U_N = create_darr3d(nU3,nU1,nU2); - U_p = create_darr3d(nU3,nU1,nU2); - U_p2 = create_darr3d(nU3,nU1,nU2); - - if(uid1!=-1){ - snprintf(username1,127,"uservar%d",uid1); - }else{ - snprintf(username1,127,"%s",ustring1); + U_N = create_darr3d (nU3, nU1, nU2); + U_p = create_darr3d (nU3, nU1, nU2); + U_p2 = create_darr3d (nU3, nU1, nU2); + + if (uid1 != -1) { + snprintf (username1, 127, "uservar%d", uid1); + } else { + snprintf (username1, 127, "%s", ustring1); } - if(uid2!=-1){ - snprintf(username2,127,"uservar%d",uid2); - }else{ - snprintf(username2,127,"%s",ustring2); + if (uid2 != -1) { + snprintf (username2, 127, "uservar%d", uid2); + } else { + snprintf (username2, 127, "%s", ustring2); } - if(uid3!=-1){ - snprintf(username3,127,"uservar%d",uid3); - }else{ - snprintf(username3,127,"%s",ustring3); + if (uid3 != -1) { + snprintf (username3, 127, "uservar%d", uid3); + } else { + snprintf (username3, 127, "%s", ustring3); } - snprintf(signalname,127,"%s",signal); + snprintf (signalname, 127, "%s", signal); // 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 %{ - double U1,U2,U3; - int suc1,suc2,suc3,suc4; - if(uid1!=-1){ - U1 = particle_getuservar_byid(_particle,uid1,&suc1); - }else{ - U1 = particle_getvar(_particle,ustring1,&suc1); + double U1, U2, U3; + int suc1, suc2, suc3, suc4; + if (uid1 != -1) { + U1 = particle_getuservar_byid (_particle, uid1, &suc1); + } else { + U1 = particle_getvar (_particle, ustring1, &suc1); } - if(uid2!=-1){ - U2 = particle_getuservar_byid(_particle,uid2,&suc2); - }else{ - U2 = particle_getvar(_particle,ustring2,&suc2); + if (uid2 != -1) { + U2 = particle_getuservar_byid (_particle, uid2, &suc2); + } else { + U2 = particle_getvar (_particle, ustring2, &suc2); } - if(uid3!=-1){ - U3 = particle_getuservar_byid(_particle,uid3,&suc3); - }else{ - U3 = particle_getvar(_particle,ustring3,&suc3); + if (uid3 != -1) { + U3 = particle_getuservar_byid (_particle, uid3, &suc3); + } else { + U3 = particle_getvar (_particle, ustring3, &suc3); } - - int i = floor((U1-Umin1)*nU1/(Umax1-Umin1)); - int j = floor((U2-Umin2)*nU2/(Umax2-Umin2)); - int k = floor((U3-Umin3)*nU3/(Umax3-Umin3)); - if(!suc1 && i >= 0 && i < nU1) - { - if(!suc2 && j >= 0 && j < nU2) - { - if(!suc3 && k >= 0 && k < nU3) - { - double pp=particle_getvar(_particle,signal,&suc4); - double p2 = pp*pp; - #pragma acc atomic - U_N[k][i][j] = U_N[k][i][j] +1; - #pragma acc atomic - U_p[k][i][j] = U_p[k][i][j] + pp; - #pragma acc atomic - U_p2[k][i][j] = U_p2[k][i][j] + p2; - SCATTER; - } - } + + int i = floor ((U1 - Umin1) * nU1 / (Umax1 - Umin1)); + int j = floor ((U2 - Umin2) * nU2 / (Umax2 - Umin2)); + int k = floor ((U3 - Umin3) * nU3 / (Umax3 - Umin3)); + if (!suc1 && i >= 0 && i < nU1) { + if (!suc2 && j >= 0 && j < nU2) { + if (!suc3 && k >= 0 && k < nU3) { + double pp = particle_getvar (_particle, signal, &suc4); + double p2 = pp * pp; + #pragma acc atomic + U_N[k][i][j] = U_N[k][i][j] + 1; + #pragma acc atomic + U_p[k][i][j] = U_p[k][i][j] + pp; + #pragma acc atomic + U_p2[k][i][j] = U_p2[k][i][j] + p2; + SCATTER; + } } - + } + 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); } %} @@ -155,24 +153,19 @@ SAVE int k; char filename_k[256]; char label_k[128]; - for(k=0; k 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("L_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("L_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - L_N = create_darr1d(nL); - L_p = create_darr1d(nL); - L_p2 = create_darr1d(nL); + L_N = create_darr1d (nL); + L_p = create_darr1d (nL); + L_p2 = create_darr1d (nL); // 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 %{ PROP_Z0; - if (x>xmin && xymin && y= 0 && i < nL) - { - double p2 = p*p; + if (x > xmin && x < xmax && y > ymin && y < ymax) { + double L = (2 * PI / V2K) / sqrt (vx * vx + vy * vy + vz * vz); + int i = floor ((L - Lmin) * nL / (Lmax - Lmin)); + if (i >= 0 && i < nL) { + double p2 = p * p; #pragma acc atomic - L_N[i] = L_N[i] +1; + L_N[i] = L_N[i] + 1; #pragma acc atomic L_p[i] = L_p[i] + p; #pragma acc atomic @@ -102,37 +107,28 @@ TRACE } } 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( - "Wavelength monitor", - "Wavelength [AA]", - "Intensity", - "L", Lmin, Lmax, nL, - &L_N[0],&L_p[0],&L_p2[0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_1D ("Wavelength monitor", "Wavelength [AA]", "Intensity", "L", Lmin, Lmax, nL, &L_N[0], &L_p[0], &L_p2[0], filename); + } %} FINALLY %{ - destroy_darr1d(L_N); - destroy_darr1d(L_p); - destroy_darr1d(L_p2); + destroy_darr1d (L_N); + destroy_darr1d (L_p); + destroy_darr1d (L_p2); %} 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 diff --git a/mcstas-comps/monitors/MeanPolLambda_monitor.comp b/mcstas-comps/monitors/MeanPolLambda_monitor.comp index 0f8c0776b..e549d1051 100644 --- a/mcstas-comps/monitors/MeanPolLambda_monitor.comp +++ b/mcstas-comps/monitors/MeanPolLambda_monitor.comp @@ -65,37 +65,41 @@ DECLARE INITIALIZE %{ // Check that input parameteters makes sense - if (Lmax<=Lmin) { - fprintf(stderr, "Pol_monitor: %s: l1 <= l0!\n" - "ERROR. Exiting", - NAME_CURRENT_COMP); - exit(1); + if (Lmax <= Lmin) { + fprintf (stderr, + "Pol_monitor: %s: l1 <= l0!\n" + "ERROR. Exiting", + NAME_CURRENT_COMP); + exit (1); } - if (mx==0 && my==0 && mz==0) { - fprintf(stderr, "Pol_monitor: %s: NULL vector defined!\n" - "ERROR (mx, my, mz). Exiting", - NAME_CURRENT_COMP); - exit(1); + if (mx == 0 && my == 0 && mz == 0) { + fprintf (stderr, + "Pol_monitor: %s: NULL vector defined!\n" + "ERROR (mx, my, mz). Exiting", + NAME_CURRENT_COMP); + exit (1); } - if ((xwidth<=0) || (yheight <= 0)) { - fprintf(stderr, "Pol_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight). Exiting", - NAME_CURRENT_COMP); - exit(1); + if ((xwidth <= 0) || (yheight <= 0)) { + fprintf (stderr, + "Pol_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight). Exiting", + NAME_CURRENT_COMP); + exit (1); } // Initialize variables - NORM(mx, my, mz); + NORM (mx, my, mz); - PolL_N = create_darr1d(nL); - PolL_p = create_darr1d(nL); - PolL_p2 = create_darr1d(nL); - HelpArray = create_darr1d(nL); + PolL_N = create_darr1d (nL); + PolL_p = create_darr1d (nL); + PolL_p2 = create_darr1d (nL); + HelpArray = create_darr1d (nL); // 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 @@ -105,81 +109,76 @@ TRACE double lambda; PROP_Z0; - lambda = (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz); + lambda = (2 * PI / V2K) / sqrt (vx * vx + vy * vy + vz * vz); - if (inside_rectangle(x, y, xwidth, yheight) && - lambda > Lmin && lambda < Lmax) { + if (inside_rectangle (x, y, xwidth, yheight) && lambda > Lmin && lambda < Lmax) { - pol_proj = scalar_prod(mx, my, mz, sx, sy, sz); + pol_proj = scalar_prod (mx, my, mz, sx, sy, sz); - i= floor((lambda - Lmin)*nL/(Lmax - Lmin)); + i = floor ((lambda - Lmin) * nL / (Lmax - Lmin)); - double p2 = p*p; + double p2 = p * p; #pragma acc atomic - PolL_N[i] = PolL_N[i]+1; + PolL_N[i] = PolL_N[i] + 1; #pragma acc atomic HelpArray[i] = HelpArray[i] + p; #pragma acc atomic - PolL_p[i] = PolL_p[i] + pol_proj*p; + PolL_p[i] = PolL_p[i] + pol_proj * p; #pragma acc atomic - PolL_p2[i] = PolL_p2[i] + pol_proj*pol_proj*p; /* Shouldn't this be p2?? */ + PolL_p2[i] = PolL_p2[i] + pol_proj * pol_proj * p; /* Shouldn't this be p2?? */ SCATTER; } 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) { - int i; - double mpifactor=1; - - /*1st order correction: As polarisation is not an additive signal we should downweight by the node count. - This asssumes that all nodes have equal weight sum.*/ -#ifdef USE_MPI - mpifactor=mpi_node_count; -#endif - - // Average output - for (i=0; i 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } - if ((xmin >= xmax) || (ymin >= ymax)) { - printf("Monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); - } + if ((xmin >= xmax) || (ymin >= ymax)) { + printf ("Monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); + } %} TRACE %{ - PROP_Z0; - if (x>xmin && xymin && y xmin && x < xmax && y > ymin && y < ymax) { + double p2 = p * p; + #pragma acc atomic + Nsum = Nsum + 1; + #pragma acc atomic + psum = psum + p; + #pragma acc atomic + p2sum = p2sum + p2; + SCATTER; + } + if (restore_neutron) { + RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); + } %} SAVE %{ - char title[1024]; - sprintf(title, "Single monitor %s", NAME_CURRENT_COMP); - DETECTOR_OUT_0D(title, Nsum, psum, p2sum); + char title[1024]; + sprintf (title, "Single monitor %s", NAME_CURRENT_COMP); + DETECTOR_OUT_0D (title, Nsum, psum, p2sum); %} 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 diff --git a/mcstas-comps/monitors/Monitor_4PI.comp b/mcstas-comps/monitors/Monitor_4PI.comp index 7a1f19915..f90719551 100644 --- a/mcstas-comps/monitors/Monitor_4PI.comp +++ b/mcstas-comps/monitors/Monitor_4PI.comp @@ -45,35 +45,35 @@ SETTING PARAMETERS () DECLARE %{ - double Nsum; - double psum; - double p2sum; + double Nsum; + double psum; + double p2sum; %} INITIALIZE %{ - Nsum = 0; - psum = 0; - p2sum = 0; + Nsum = 0; + psum = 0; + p2sum = 0; %} TRACE %{ - double p2 = p*p; - #pragma acc atomic - Nsum = Nsum + 1; - #pragma acc atomic - psum = psum + p; - #pragma acc atomic - p2sum = p2sum + p2; - SCATTER; + double p2 = p * p; + #pragma acc atomic + Nsum = Nsum + 1; + #pragma acc atomic + psum = psum + p; + #pragma acc atomic + p2sum = p2sum + p2; + SCATTER; %} SAVE %{ - char full_name[1024]; - sprintf(full_name, "4PI monitor %s", NAME_CURRENT_COMP); - DETECTOR_OUT_0D(full_name, Nsum, psum, p2sum); + char full_name[1024]; + sprintf (full_name, "4PI monitor %s", NAME_CURRENT_COMP); + DETECTOR_OUT_0D (full_name, Nsum, psum, p2sum); %} MCDISPLAY diff --git a/mcstas-comps/monitors/Monitor_nD.comp b/mcstas-comps/monitors/Monitor_nD.comp index ecee1c296..1a4cc7f39 100644 --- a/mcstas-comps/monitors/Monitor_nD.comp +++ b/mcstas-comps/monitors/Monitor_nD.comp @@ -260,152 +260,160 @@ DECLARE INITIALIZE %{ char tmp[CHAR_BUF_LENGTH]; - strcpy(Vars.compcurname, NAME_CURRENT_COMP); - Vars.compcurindex=INDEX_CURRENT_COMP; + strcpy (Vars.compcurname, NAME_CURRENT_COMP); + Vars.compcurindex = INDEX_CURRENT_COMP; if (options != NULL) - strncpy(Vars.option, options, CHAR_BUF_LENGTH); + strncpy (Vars.option, options, CHAR_BUF_LENGTH); else { - strcpy(Vars.option, "x y"); - printf("Monitor_nD: %s has no option specified. Setting to PSD ('x y') monitor.\n", NAME_CURRENT_COMP); + strcpy (Vars.option, "x y"); + printf ("Monitor_nD: %s has no option specified. Setting to PSD ('x y') monitor.\n", NAME_CURRENT_COMP); } Vars.compcurpos = POS_A_CURRENT_COMP; - if (strstr(Vars.option, "source")) - strcat(Vars.option, " list, x y z vx vy vz t sx sy sz "); + if (strstr (Vars.option, "source")) + strcat (Vars.option, " list, x y z vx vy vz t sx sy sz "); - if (bins) { sprintf(tmp, " all bins=%ld ", (long)bins); strcat(Vars.option, tmp); } - if (min > -FLT_MAX && max < FLT_MAX) { sprintf(tmp, " all limits=[%g %g]", min, max); strcat(Vars.option, tmp); } - else if (min > -FLT_MAX) { sprintf(tmp, " all min=%g", min); strcat(Vars.option, tmp); } - else if (max < FLT_MAX) { sprintf(tmp, " all max=%g", max); strcat(Vars.option, tmp); } + if (bins) { + sprintf (tmp, " all bins=%ld ", (long)bins); + strcat (Vars.option, tmp); + } + if (min > -FLT_MAX && max < FLT_MAX) { + sprintf (tmp, " all limits=[%g %g]", min, max); + strcat (Vars.option, tmp); + } else if (min > -FLT_MAX) { + sprintf (tmp, " all min=%g", min); + strcat (Vars.option, tmp); + } else if (max < FLT_MAX) { + sprintf (tmp, " all max=%g", max); + strcat (Vars.option, tmp); + } /* transfer, "zero", and check username- and user variable strings to Vars struct*/ - strncpy(Vars.UserName1, - username1 && strlen(username1) && strcmp(username1, "0") && strcmp(username1, "NULL") ? - username1 : "", 128); - strncpy(Vars.UserName2, - username2 && strlen(username2) && strcmp(username2, "0") && strcmp(username2, "NULL") ? - username2 : "", 128); - strncpy(Vars.UserName3, - username3 && strlen(username3) && strcmp(username3, "0") && strcmp(username3, "NULL") ? - username3 : "", 128); - if(user1 && strlen(user1) && strcmp(user1, "0") && strcmp(user1, "NULL")){ - strncpy(Vars.UserVariable1,user1,128); - int fail;_class_particle testparticle; - particle_getvar(&testparticle,Vars.UserVariable1,&fail); - if(fail){ - fprintf(stderr,"Warning (%s): user1=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n",NAME_CURRENT_COMP,user1); + strncpy (Vars.UserName1, username1&& strlen (username1) && strcmp (username1, "0") && strcmp (username1, "NULL") ? username1 : "", 128); + strncpy (Vars.UserName2, username2&& strlen (username2) && strcmp (username2, "0") && strcmp (username2, "NULL") ? username2 : "", 128); + strncpy (Vars.UserName3, username3&& strlen (username3) && strcmp (username3, "0") && strcmp (username3, "NULL") ? username3 : "", 128); + if (user1 && strlen (user1) && strcmp (user1, "0") && strcmp (user1, "NULL")) { + strncpy (Vars.UserVariable1, user1, 128); + int fail; + _class_particle testparticle; + particle_getvar (&testparticle, Vars.UserVariable1, &fail); + if (fail) { + fprintf (stderr, "Warning (%s): user1=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n", NAME_CURRENT_COMP, user1); } } - if(user2 && strlen(user2) && strcmp(user2, "0") && strcmp(user2, "NULL")){ - strncpy(Vars.UserVariable2,user2,128); - int fail;_class_particle testparticle; - particle_getvar(&testparticle,Vars.UserVariable2,&fail); - if(fail){ - fprintf(stderr,"Warning (%s): user2=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n",NAME_CURRENT_COMP,user2); + if (user2 && strlen (user2) && strcmp (user2, "0") && strcmp (user2, "NULL")) { + strncpy (Vars.UserVariable2, user2, 128); + int fail; + _class_particle testparticle; + particle_getvar (&testparticle, Vars.UserVariable2, &fail); + if (fail) { + fprintf (stderr, "Warning (%s): user2=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n", NAME_CURRENT_COMP, user2); } } - if(user3 && strlen(user3) && strcmp(user3, "0") && strcmp(user3, "NULL")){ - strncpy(Vars.UserVariable3,user3,128); - int fail;_class_particle testparticle; - particle_getvar(&testparticle,Vars.UserVariable3,&fail); - if(fail){ - fprintf(stderr,"Warning (%s): user3=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n",NAME_CURRENT_COMP,user3); + if (user3 && strlen (user3) && strcmp (user3, "0") && strcmp (user3, "NULL")) { + strncpy (Vars.UserVariable3, user3, 128); + int fail; + _class_particle testparticle; + particle_getvar (&testparticle, Vars.UserVariable3, &fail); + if (fail) { + fprintf (stderr, "Warning (%s): user3=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n", NAME_CURRENT_COMP, user3); } } - + /*sanitize parameters set for curved shapes*/ - if(strstr(Vars.option,"cylinder") || strstr(Vars.option,"banana") || strstr(Vars.option,"sphere")){ + if (strstr (Vars.option, "cylinder") || strstr (Vars.option, "banana") || strstr (Vars.option, "sphere")) { /*this _is_ an explicit curved shape. Should have a radius. Inherit from xwidth or zdepth (diameters), x has precedence.*/ - if (!radius){ - if(xwidth){ - radius=xwidth/2.0; - }else{ - radius=zdepth/2.0; + if (!radius) { + if (xwidth) { + radius = xwidth / 2.0; + } else { + radius = zdepth / 2.0; } - }else{ - xwidth=2*radius; + } else { + xwidth = 2 * radius; } - if(!yheight){ + if (!yheight) { /*if not set - use the diameter as height for the curved object. This will likely only happen for spheres*/ - yheight=2*radius; + yheight = 2 * radius; } - }else if (radius) { + } else if (radius) { /*radius is set - this must be a curved shape. Infer shape from yheight, and set remaining values (xwidth etc. They are used inside monitor_nd-lib.*/ - xwidth = zdepth = 2*radius; - if (yheight){ + xwidth = zdepth = 2 * radius; + if (yheight) { /*a height is given (and no shape explitly set - assume cylinder*/ - strcat(Vars.option, " banana"); - }else { - strcat(Vars.option, " sphere"); - yheight=2*radius; + strcat (Vars.option, " banana"); + } else { + strcat (Vars.option, " sphere"); + yheight = 2 * radius; } } - int offflag=0; - if (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) { + int offflag = 0; + if (geometry && strlen (geometry) && strcmp (geometry, "0") && strcmp (geometry, "NULL")) { #ifndef USE_OFF - fprintf(stderr,"Error: You are attempting to use an OFF geometry without -DUSE_OFF. You will need to recompile with that define set!\n"); - exit(-1); + fprintf (stderr, "Error: You are attempting to use an OFF geometry without -DUSE_OFF. You will need to recompile with that define set!\n"); + exit (-1); #else - if (!off_init( geometry, xwidth, yheight, zdepth, 1, &offdata )) { - printf("Monitor_nD: %s could not initiate the OFF geometry %s. \n" - " Defaulting to normal Monitor dimensions.\n", - NAME_CURRENT_COMP, geometry); - strcpy(geometry, ""); + if (!off_init (geometry, xwidth, yheight, zdepth, 1, &offdata)) { + printf ("Monitor_nD: %s could not initiate the OFF geometry %s. \n" + " Defaulting to normal Monitor dimensions.\n", + NAME_CURRENT_COMP, geometry); + strcpy (geometry, ""); } else { - offflag=1; + offflag = 1; } #endif } - if (!radius && !xwidth && !yheight && !zdepth && !xmin && !xmax && !ymin && !ymax && - !strstr(Vars.option, "previous") && (!geometry || !strlen(geometry))) - exit(printf("Monitor_nD: %s has no dimension specified. Aborting (radius, xwidth, yheight, zdepth, previous, geometry).\n", NAME_CURRENT_COMP)); + if (!radius && !xwidth && !yheight && !zdepth && !xmin && !xmax && !ymin && !ymax && !strstr (Vars.option, "previous") && (!geometry || !strlen (geometry))) + exit (printf ("Monitor_nD: %s has no dimension specified. Aborting (radius, xwidth, yheight, zdepth, previous, geometry).\n", NAME_CURRENT_COMP)); - Monitor_nD_Init(&DEFS, &Vars, xwidth, yheight, zdepth, xmin,xmax,ymin,ymax,zmin,zmax,offflag); + Monitor_nD_Init (&DEFS, &Vars, xwidth, yheight, zdepth, xmin, xmax, ymin, ymax, zmin, zmax, offflag); if (Vars.Flag_OFF) { - offdata.mantidflag=Vars.Flag_mantid; - offdata.mantidoffset=Vars.Coord_Min[Vars.Coord_Number-1]; + offdata.mantidflag = Vars.Flag_mantid; + offdata.mantidoffset = Vars.Coord_Min[Vars.Coord_Number - 1]; } - - if (filename && strlen(filename) && strcmp(filename,"NULL") && strcmp(filename,"0")) - strncpy(Vars.Mon_File, filename, 128); + if (filename && strlen (filename) && strcmp (filename, "NULL") && strcmp (filename, "0")) + strncpy (Vars.Mon_File, filename, 128); /* check if user given filename with ext will be used more than once */ - if ( ((Vars.Flag_Multiple && Vars.Coord_Number > 1) || Vars.Flag_List) && strchr(Vars.Mon_File,'.') ) - { char *XY; XY = strrchr(Vars.Mon_File,'.'); *XY='_'; } + if (((Vars.Flag_Multiple && Vars.Coord_Number > 1) || Vars.Flag_List) && strchr (Vars.Mon_File, '.')) { + char* XY; + XY = strrchr (Vars.Mon_File, '.'); + *XY = '_'; + } - if (restore_neutron) Vars.Flag_parallel=1; + if (restore_neutron) + Vars.Flag_parallel = 1; detector.m = 0; -#ifdef USE_MPI -MPI_MASTER( - if (strstr(Vars.option, "auto") && mpi_node_count > 1) - printf("Monitor_nD: %s is using automatic limits option 'auto' together with MPI.\n" - "WARNING this may create incorrect distributions (but integrated flux will be right).\n", NAME_CURRENT_COMP); -); -#else -#ifdef OPENACC - if (strstr(Vars.option, "auto")) - printf("Monitor_nD: %s is requesting automatic limits option 'auto' together with OpenACC.\n" - "WARNING this feature is NOT supported using OpenACC and has been disabled!\n", NAME_CURRENT_COMP); -#endif -#endif - + #ifdef USE_MPI + MPI_MASTER (if (strstr (Vars.option, "auto") && mpi_node_count > 1) + printf ("Monitor_nD: %s is using automatic limits option 'auto' together with MPI.\n" + "WARNING this may create incorrect distributions (but integrated flux will be right).\n", + NAME_CURRENT_COMP);); + #else + #ifdef OPENACC + if (strstr (Vars.option, "auto")) + printf ("Monitor_nD: %s is requesting automatic limits option 'auto' together with OpenACC.\n" + "WARNING this feature is NOT supported using OpenACC and has been disabled!\n", + NAME_CURRENT_COMP); + #endif + #endif %} TRACE %{ - double transmit_he3=1.0; - double multiplier_capture=1.0; - double t0 = 0; - double t1 = 0; - int pp; - int intersect = 0; - char Flag_Restore = 0; + double transmit_he3 = 1.0; + double multiplier_capture = 1.0; + double t0 = 0; + double t1 = 0; + int pp; + int intersect = 0; + char Flag_Restore = 0; #ifdef OPENACC #ifdef USE_OFF @@ -414,225 +422,203 @@ TRACE #else #define thread_offdata offdata #endif - + /* this is done automatically STORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); */ #ifdef USE_OFF - if (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) - { + if (geometry && strlen (geometry) && strcmp (geometry, "0") && strcmp (geometry, "NULL")) { /* determine intersections with object */ - intersect = off_intersect_all(&t0, &t1, NULL, NULL, - x,y,z, vx, vy, vz, 0, 0, 0, &thread_offdata ); + intersect = off_intersect_all (&t0, &t1, NULL, NULL, x, y, z, vx, vy, vz, 0, 0, 0, &thread_offdata); if (Vars.Flag_mantid) { - if(intersect) { - Vars.OFF_polyidx=thread_offdata.nextintersect; + if (intersect) { + Vars.OFF_polyidx = thread_offdata.nextintersect; } else { - Vars.OFF_polyidx=-1; + Vars.OFF_polyidx = -1; } } - } - else - #endif - if ( (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) - || (abs(Vars.Flag_Shape) == DEFS.SHAPE_DISK) ) /* square xy or disk xy */ - { - // propagate to xy plane and find intersection - // make sure the event is recoverable afterwards - t0 = t; - ALLOW_BACKPROP; - PROP_Z0; - if ( (t>=t0) && (z==0.0) ) // forward propagation to xy plane was successful + } else + #endif + if ((abs (Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_DISK)) /* square xy or disk xy */ { - if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) + // propagate to xy plane and find intersection + // make sure the event is recoverable afterwards + t0 = t; + ALLOW_BACKPROP; + PROP_Z0; + if ((t >= t0) && (z == 0.0)) // forward propagation to xy plane was successful { - // square xy - intersect = (x>=Vars.mxmin && x<=Vars.mxmax && y>=Vars.mymin && y<=Vars.mymax); - } - else - { - // disk xy - intersect = (SQR(x) + SQR(y)) <= SQR(Vars.Sphere_Radius); + if (abs (Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) { + // square xy + intersect = (x >= Vars.mxmin && x <= Vars.mxmax && y >= Vars.mymin && y <= Vars.mymax); + } else { + // disk xy + intersect = (SQR (x) + SQR (y)) <= SQR (Vars.Sphere_Radius); + } + } else { + intersect = 0; } - } - else + } else if (abs (Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) /* sphere */ { - intersect=0; - } - } - else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) /* sphere */ - { - intersect = sphere_intersect(&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius); - /* intersect = (intersect && t0 > 0); */ - } - else if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA)) /* cylinder */ - { - intersect = cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius, Vars.Cylinder_Height); - } - else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX) /* box */ - { - intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, - fabs(Vars.mxmax-Vars.mxmin), fabs(Vars.mymax-Vars.mymin), fabs(Vars.mzmax-Vars.mzmin)); - } - else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_PREVIOUS) /* previous comp */ - { intersect = 1; } - - if (intersect) - { - if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) - || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA) - || (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) ) + intersect = sphere_intersect (&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius); + /* intersect = (intersect && t0 > 0); */ + } else if ((abs (Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_BANANA)) /* cylinder */ { + intersect = cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius, Vars.Cylinder_Height); + } else if (abs (Vars.Flag_Shape) == DEFS.SHAPE_BOX) /* box */ + { + intersect = box_intersect (&t0, &t1, x, y, z, vx, vy, vz, fabs (Vars.mxmax - Vars.mxmin), fabs (Vars.mymax - Vars.mymin), fabs (Vars.mzmax - Vars.mzmin)); + } else if (abs (Vars.Flag_Shape) == DEFS.SHAPE_PREVIOUS) /* previous comp */ + { + intersect = 1; + } + + if (intersect) { + if ((abs (Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_BOX) + || (abs (Vars.Flag_Shape) == DEFS.SHAPE_BANANA) || (geometry && strlen (geometry) && strcmp (geometry, "0") && strcmp (geometry, "NULL"))) { /* check if we have to remove the top/bottom with BANANA shape */ - if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA) { - if (intersect == 1) { // Entered and left through sides - if (t0 < 0 && t1 > 0) { - t0 = t; /* neutron was already inside ! */ - } - if (t1 < 0 && t0 > 0) { /* neutron exit before entering !! */ - t1 = t; - } - /* t0 is now time of incoming intersection with the detection area */ - if ((Vars.Flag_Shape < 0) && (t1 > 0)) { - PROP_DT(t1); /* t1 outgoing beam */ - } else { - PROP_DT(t0); /* t0 incoming beam */ - } - } else if (intersect == 3 || intersect == 5) { // Entered from top or bottom, left through side - if ((Vars.Flag_Shape < 0) && (t1 > 0)) { - PROP_DT(t1); /* t1 outgoing beam */ - } else { - intersect=0; - Flag_Restore=1; - } - } else if (intersect == 9 || intersect == 17) { // Entered through side, left from top or bottom - if ((Vars.Flag_Shape < 0) && (t1 > 0)) { - intersect=0; - Flag_Restore=1; - } else { - PROP_DT(t0); /* t0 incoming beam */ - } - } else if (intersect == 13 || intersect == 19) { // Went through top/bottom on entry and exit - intersect=0; - Flag_Restore=1; - } else { - printf("Cylinder_intersect returned unexpected value %i\n", intersect); - } + if (abs (Vars.Flag_Shape) == DEFS.SHAPE_BANANA) { + if (intersect == 1) { // Entered and left through sides + if (t0 < 0 && t1 > 0) { + t0 = t; /* neutron was already inside ! */ + } + if (t1 < 0 && t0 > 0) { /* neutron exit before entering !! */ + t1 = t; + } + /* t0 is now time of incoming intersection with the detection area */ + if ((Vars.Flag_Shape < 0) && (t1 > 0)) { + PROP_DT (t1); /* t1 outgoing beam */ + } else { + PROP_DT (t0); /* t0 incoming beam */ + } + } else if (intersect == 3 || intersect == 5) { // Entered from top or bottom, left through side + if ((Vars.Flag_Shape < 0) && (t1 > 0)) { + PROP_DT (t1); /* t1 outgoing beam */ + } else { + intersect = 0; + Flag_Restore = 1; + } + } else if (intersect == 9 || intersect == 17) { // Entered through side, left from top or bottom + if ((Vars.Flag_Shape < 0) && (t1 > 0)) { + intersect = 0; + Flag_Restore = 1; + } else { + PROP_DT (t0); /* t0 incoming beam */ + } + } else if (intersect == 13 || intersect == 19) { // Went through top/bottom on entry and exit + intersect = 0; + Flag_Restore = 1; } else { - // All other shapes than the BANANA - if (t0 < 0 && t1 > 0) - t0 = t; /* neutron was already inside ! */ - if (t1 < 0 && t0 > 0) /* neutron exit before entering !! */ - t1 = t; - /* t0 is now time of incoming intersection with the detection area */ - if ((Vars.Flag_Shape < 0) && (t1 > 0)) - PROP_DT(t1); /* t1 outgoing beam */ - else - PROP_DT(t0); /* t0 incoming beam */ + printf ("Cylinder_intersect returned unexpected value %i\n", intersect); } - + } else { + // All other shapes than the BANANA + if (t0 < 0 && t1 > 0) + t0 = t; /* neutron was already inside ! */ + if (t1 < 0 && t0 > 0) /* neutron exit before entering !! */ + t1 = t; + /* t0 is now time of incoming intersection with the detection area */ + if ((Vars.Flag_Shape < 0) && (t1 > 0)) + PROP_DT (t1); /* t1 outgoing beam */ + else + PROP_DT (t0); /* t0 incoming beam */ + } + /* Final test if we are on lid / bottom of banana/sphere */ - if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA || abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) { - if (Vars.Cylinder_Height && fabs(y) >= Vars.Cylinder_Height/2 - FLT_EPSILON) { - intersect=0; - Flag_Restore=1; + if (abs (Vars.Flag_Shape) == DEFS.SHAPE_BANANA || abs (Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) { + if (Vars.Cylinder_Height && fabs (y) >= Vars.Cylinder_Height / 2 - FLT_EPSILON) { + intersect = 0; + Flag_Restore = 1; } } } } - if (intersect) - { + if (intersect) { /* Now get the data to monitor: current or keep from PreMonitor */ -/* if (Vars.Flag_UsePreMonitor != 1)*/ -/* {*/ -/* Vars.cp = p;*/ -/* Vars.cx = x;*/ -/* Vars.cvx = vx;*/ -/* Vars.csx = sx;*/ -/* Vars.cy = y;*/ -/* Vars.cvy = vy;*/ -/* Vars.csy = sy;*/ -/* Vars.cz = z;*/ -/* Vars.cvz = vz;*/ -/* Vars.csz = sz;*/ -/* Vars.ct = t;*/ -/* }*/ - - if ((Vars.He3_pressure > 0) && (t1 != t0) && ((abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX))) - { - transmit_he3 = exp(-7.417*Vars.He3_pressure*fabs(t1-t0)*2*PI*K2V); + /* if (Vars.Flag_UsePreMonitor != 1)*/ + /* {*/ + /* Vars.cp = p;*/ + /* Vars.cx = x;*/ + /* Vars.cvx = vx;*/ + /* Vars.csx = sx;*/ + /* Vars.cy = y;*/ + /* Vars.cvy = vy;*/ + /* Vars.csy = sy;*/ + /* Vars.cz = z;*/ + /* Vars.cvz = vz;*/ + /* Vars.csz = sz;*/ + /* Vars.ct = t;*/ + /* }*/ + + if ((Vars.He3_pressure > 0) && (t1 != t0) + && ((abs (Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_BOX))) { + transmit_he3 = exp (-7.417 * Vars.He3_pressure * fabs (t1 - t0) * 2 * PI * K2V); /* will monitor the absorbed part */ - p = p * (1-transmit_he3); + p = p * (1 - transmit_he3); } - if (Vars.Flag_capture) - { - multiplier_capture = V2K*sqrt(vx*vx+vy*vy+vz*vz); - if (multiplier_capture != 0) multiplier_capture = 2*PI/multiplier_capture; /* lambda. lambda(2200 m/2) = 1.7985 Angs */ - p = p * multiplier_capture/1.7985; + if (Vars.Flag_capture) { + multiplier_capture = V2K * sqrt (vx * vx + vy * vy + vz * vz); + if (multiplier_capture != 0) + multiplier_capture = 2 * PI / multiplier_capture; /* lambda. lambda(2200 m/2) = 1.7985 Angs */ + p = p * multiplier_capture / 1.7985; } - pp = Monitor_nD_Trace(&DEFS, &Vars, _particle); - if (pp==0.0) - { + pp = Monitor_nD_Trace (&DEFS, &Vars, _particle); + if (pp == 0.0) { ABSORB; - } - else if(pp==1) - { + } else if (pp == 1) { SCATTER; } /*set weight to undetected part if capture and/or he3_pressure*/ - if (Vars.He3_pressure > 0){ + if (Vars.He3_pressure > 0) { /* after monitor, only remains 1-p_detect */ - p = p * transmit_he3/(1.0-transmit_he3); + p = p * transmit_he3 / (1.0 - transmit_he3); } - if (Vars.Flag_capture){ - p = p / multiplier_capture*1.7985; + if (Vars.Flag_capture) { + p = p / multiplier_capture * 1.7985; } if (Vars.Flag_parallel) /* back to neutron state before detection */ Flag_Restore = 1; } /* end if intersection */ else { - if (Vars.Flag_Absorb && !Vars.Flag_parallel) - { + if (Vars.Flag_Absorb && !Vars.Flag_parallel) { // restore neutron ray before absorbing for correct mcdisplay - 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); ABSORB; - } - else Flag_Restore = 1; /* no intersection, back to previous state */ + } else + Flag_Restore = 1; /* no intersection, back to previous state */ } - if (Flag_Restore) - { - RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); + if (Flag_Restore) { + RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); } %} SAVE %{ -if (!nowritefile) { - /* save results, but do not free pointers */ - detector = Monitor_nD_Save(&DEFS, &Vars); -} + if (!nowritefile) { + /* save results, but do not free pointers */ + detector = Monitor_nD_Save (&DEFS, &Vars); + } %} FINALLY %{ /* free pointers */ - Monitor_nD_Finally(&DEFS, &Vars); + Monitor_nD_Finally (&DEFS, &Vars); %} MCDISPLAY %{ - if (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) - { - off_display(offdata); + if (geometry && strlen (geometry) && strcmp (geometry, "0") && strcmp (geometry, "NULL")) { + off_display (offdata); } else { - Monitor_nD_McDisplay(&DEFS, &Vars); + Monitor_nD_McDisplay (&DEFS, &Vars); } %} diff --git a/mcstas-comps/monitors/Monitor_nD_noacc.comp b/mcstas-comps/monitors/Monitor_nD_noacc.comp index c207af94f..0dc9361d1 100644 --- a/mcstas-comps/monitors/Monitor_nD_noacc.comp +++ b/mcstas-comps/monitors/Monitor_nD_noacc.comp @@ -261,135 +261,142 @@ DECLARE INITIALIZE %{ char tmp[CHAR_BUF_LENGTH]; - strcpy(Vars.compcurname, NAME_CURRENT_COMP); - Vars.compcurindex=INDEX_CURRENT_COMP; + strcpy (Vars.compcurname, NAME_CURRENT_COMP); + Vars.compcurindex = INDEX_CURRENT_COMP; if (options != NULL) - strncpy(Vars.option, options, CHAR_BUF_LENGTH); + strncpy (Vars.option, options, CHAR_BUF_LENGTH); else { - strcpy(Vars.option, "x y"); - printf("Monitor_nD: %s has no option specified. Setting to PSD ('x y') monitor.\n", NAME_CURRENT_COMP); + strcpy (Vars.option, "x y"); + printf ("Monitor_nD: %s has no option specified. Setting to PSD ('x y') monitor.\n", NAME_CURRENT_COMP); } Vars.compcurpos = POS_A_CURRENT_COMP; - if (strstr(Vars.option, "source")) - strcat(Vars.option, " list, x y z vx vy vz t sx sy sz "); + if (strstr (Vars.option, "source")) + strcat (Vars.option, " list, x y z vx vy vz t sx sy sz "); - if (bins) { sprintf(tmp, " all bins=%ld ", (long)bins); strcat(Vars.option, tmp); } - if (min > -FLT_MAX && max < FLT_MAX) { sprintf(tmp, " all limits=[%g %g]", min, max); strcat(Vars.option, tmp); } - else if (min > -FLT_MAX) { sprintf(tmp, " all min=%g", min); strcat(Vars.option, tmp); } - else if (max < FLT_MAX) { sprintf(tmp, " all max=%g", max); strcat(Vars.option, tmp); } + if (bins) { + sprintf (tmp, " all bins=%ld ", (long)bins); + strcat (Vars.option, tmp); + } + if (min > -FLT_MAX && max < FLT_MAX) { + sprintf (tmp, " all limits=[%g %g]", min, max); + strcat (Vars.option, tmp); + } else if (min > -FLT_MAX) { + sprintf (tmp, " all min=%g", min); + strcat (Vars.option, tmp); + } else if (max < FLT_MAX) { + sprintf (tmp, " all max=%g", max); + strcat (Vars.option, tmp); + } /* transfer, "zero", and check username- and user variable strings to Vars struct*/ - strncpy(Vars.UserName1, - username1 && strlen(username1) && strcmp(username1, "0") && strcmp(username1, "NULL") ? - username1 : "", 128); - strncpy(Vars.UserName2, - username2 && strlen(username2) && strcmp(username2, "0") && strcmp(username2, "NULL") ? - username2 : "", 128); - strncpy(Vars.UserName3, - username3 && strlen(username3) && strcmp(username3, "0") && strcmp(username3, "NULL") ? - username3 : "", 128); - if(user1 && strlen(user1) && strcmp(user1, "0") && strcmp(user1, "NULL")){ - strncpy(Vars.UserVariable1,user1,128); - int fail;_class_particle testparticle; - particle_getvar(&testparticle,Vars.UserVariable1,&fail); - if(fail){ - fprintf(stderr,"Warning (%s): user1=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n",NAME_CURRENT_COMP,user1); + strncpy (Vars.UserName1, username1&& strlen (username1) && strcmp (username1, "0") && strcmp (username1, "NULL") ? username1 : "", 128); + strncpy (Vars.UserName2, username2&& strlen (username2) && strcmp (username2, "0") && strcmp (username2, "NULL") ? username2 : "", 128); + strncpy (Vars.UserName3, username3&& strlen (username3) && strcmp (username3, "0") && strcmp (username3, "NULL") ? username3 : "", 128); + if (user1 && strlen (user1) && strcmp (user1, "0") && strcmp (user1, "NULL")) { + strncpy (Vars.UserVariable1, user1, 128); + int fail; + _class_particle testparticle; + particle_getvar (&testparticle, Vars.UserVariable1, &fail); + if (fail) { + fprintf (stderr, "Warning (%s): user1=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n", NAME_CURRENT_COMP, user1); } } - if(user2 && strlen(user2) && strcmp(user2, "0") && strcmp(user2, "NULL")){ - strncpy(Vars.UserVariable2,user2,128); - int fail;_class_particle testparticle; - particle_getvar(&testparticle,Vars.UserVariable2,&fail); - if(fail){ - fprintf(stderr,"Warning (%s): user2=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n",NAME_CURRENT_COMP,user2); + if (user2 && strlen (user2) && strcmp (user2, "0") && strcmp (user2, "NULL")) { + strncpy (Vars.UserVariable2, user2, 128); + int fail; + _class_particle testparticle; + particle_getvar (&testparticle, Vars.UserVariable2, &fail); + if (fail) { + fprintf (stderr, "Warning (%s): user2=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n", NAME_CURRENT_COMP, user2); } } - if(user3 && strlen(user3) && strcmp(user3, "0") && strcmp(user3, "NULL")){ - strncpy(Vars.UserVariable3,user3,128); - int fail;_class_particle testparticle; - particle_getvar(&testparticle,Vars.UserVariable3,&fail); - if(fail){ - fprintf(stderr,"Warning (%s): user3=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n",NAME_CURRENT_COMP,user3); + if (user3 && strlen (user3) && strcmp (user3, "0") && strcmp (user3, "NULL")) { + strncpy (Vars.UserVariable3, user3, 128); + int fail; + _class_particle testparticle; + particle_getvar (&testparticle, Vars.UserVariable3, &fail); + if (fail) { + fprintf (stderr, "Warning (%s): user3=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n", NAME_CURRENT_COMP, user3); } } - + /*sanitize parameters set for curved shapes*/ - if(strstr(Vars.option,"cylinder") || strstr(Vars.option,"banana") || strstr(Vars.option,"sphere")){ + if (strstr (Vars.option, "cylinder") || strstr (Vars.option, "banana") || strstr (Vars.option, "sphere")) { /*this _is_ an explicit curved shape. Should have a radius. Inherit from xwidth or zdepth (diameters), x has precedence.*/ - if (!radius){ - if(xwidth){ - radius=xwidth/2.0; - }else{ - radius=zdepth/2.0; + if (!radius) { + if (xwidth) { + radius = xwidth / 2.0; + } else { + radius = zdepth / 2.0; } - }else{ - xwidth=2*radius; + } else { + xwidth = 2 * radius; } - if(!yheight){ + if (!yheight) { /*if not set - use the diameter as height for the curved object. This will likely only happen for spheres*/ - yheight=2*radius; + yheight = 2 * radius; } - }else if (radius) { + } else if (radius) { /*radius is set - this must be a curved shape. Infer shape from yheight, and set remaining values (xwidth etc. They are used inside monitor_nd-lib.*/ - xwidth = zdepth = 2*radius; - if (yheight){ + xwidth = zdepth = 2 * radius; + if (yheight) { /*a height is given (and no shape explitly set - assume cylinder*/ - strcat(Vars.option, " banana"); - }else { - strcat(Vars.option, " sphere"); - yheight=2*radius; + strcat (Vars.option, " banana"); + } else { + strcat (Vars.option, " sphere"); + yheight = 2 * radius; } } - int offflag=0; - if (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) { + int offflag = 0; + if (geometry && strlen (geometry) && strcmp (geometry, "0") && strcmp (geometry, "NULL")) { #ifndef USE_OFF - fprintf(stderr,"Error: You are attempting to use an OFF geometry without -DUSE_OFF. You will need to recompile with that define set!\n"); - exit(-1); + fprintf (stderr, "Error: You are attempting to use an OFF geometry without -DUSE_OFF. You will need to recompile with that define set!\n"); + exit (-1); #else - if (!off_init( geometry, xwidth, yheight, zdepth, 1, &offdata )) { - printf("Monitor_nD: %s could not initiate the OFF geometry %s. \n" - " Defaulting to normal Monitor dimensions.\n", - NAME_CURRENT_COMP, geometry); - strcpy(geometry, ""); + if (!off_init (geometry, xwidth, yheight, zdepth, 1, &offdata)) { + printf ("Monitor_nD: %s could not initiate the OFF geometry %s. \n" + " Defaulting to normal Monitor dimensions.\n", + NAME_CURRENT_COMP, geometry); + strcpy (geometry, ""); } else { - offflag=1; + offflag = 1; } #endif } - if (!radius && !xwidth && !yheight && !zdepth && !xmin && !xmax && !ymin && !ymax && - !strstr(Vars.option, "previous") && (!geometry || !strlen(geometry))) - exit(printf("Monitor_nD: %s has no dimension specified. Aborting (radius, xwidth, yheight, zdepth, previous, geometry).\n", NAME_CURRENT_COMP)); + if (!radius && !xwidth && !yheight && !zdepth && !xmin && !xmax && !ymin && !ymax && !strstr (Vars.option, "previous") && (!geometry || !strlen (geometry))) + exit (printf ("Monitor_nD: %s has no dimension specified. Aborting (radius, xwidth, yheight, zdepth, previous, geometry).\n", NAME_CURRENT_COMP)); - Monitor_nd_noaccInit(&DEFS, &Vars, xwidth, yheight, zdepth, xmin,xmax,ymin,ymax,zmin,zmax,offflag); + Monitor_nd_noaccInit (&DEFS, &Vars, xwidth, yheight, zdepth, xmin, xmax, ymin, ymax, zmin, zmax, offflag); if (Vars.Flag_OFF) { - offdata.mantidflag=Vars.Flag_mantid; - offdata.mantidoffset=Vars.Coord_Min[Vars.Coord_Number-1]; + offdata.mantidflag = Vars.Flag_mantid; + offdata.mantidoffset = Vars.Coord_Min[Vars.Coord_Number - 1]; } - - if (filename && strlen(filename) && strcmp(filename,"NULL") && strcmp(filename,"0")) - strncpy(Vars.Mon_File, filename, 128); + if (filename && strlen (filename) && strcmp (filename, "NULL") && strcmp (filename, "0")) + strncpy (Vars.Mon_File, filename, 128); /* check if user given filename with ext will be used more than once */ - if ( ((Vars.Flag_Multiple && Vars.Coord_Number > 1) || Vars.Flag_List) && strchr(Vars.Mon_File,'.') ) - { char *XY; XY = strrchr(Vars.Mon_File,'.'); *XY='_'; } + if (((Vars.Flag_Multiple && Vars.Coord_Number > 1) || Vars.Flag_List) && strchr (Vars.Mon_File, '.')) { + char* XY; + XY = strrchr (Vars.Mon_File, '.'); + *XY = '_'; + } - if (restore_neutron) Vars.Flag_parallel=1; + if (restore_neutron) + Vars.Flag_parallel = 1; detector.m = 0; -#ifdef USE_MPI -MPI_MASTER( - if (strstr(Vars.option, "auto") && mpi_node_count > 1) - printf("Monitor_nD: %s is using automatic limits option 'auto' together with MPI.\n" - "WARNING this may create incorrect distributions (but integrated flux will be right).\n", NAME_CURRENT_COMP); -); -#endif - + #ifdef USE_MPI + MPI_MASTER (if (strstr (Vars.option, "auto") && mpi_node_count > 1) + printf ("Monitor_nD: %s is using automatic limits option 'auto' together with MPI.\n" + "WARNING this may create incorrect distributions (but integrated flux will be right).\n", + NAME_CURRENT_COMP);); + #endif %} TRACE @@ -401,211 +408,189 @@ TRACE #undef sprintf #undef exit #endif - - double transmit_he3=1.0; - double multiplier_capture=1.0; - double t0 = 0; - double t1 = 0; - int pp; - int intersect = 0; - char Flag_Restore = 0; + + double transmit_he3 = 1.0; + double multiplier_capture = 1.0; + double t0 = 0; + double t1 = 0; + int pp; + int intersect = 0; + char Flag_Restore = 0; #define thread_offdata offdata - + /* this is done automatically STORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); */ - if (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) - { + if (geometry && strlen (geometry) && strcmp (geometry, "0") && strcmp (geometry, "NULL")) { /* determine intersections with object */ - intersect = off_intersect_all(&t0, &t1, NULL, NULL, - x,y,z, vx, vy, vz, 0, 0, 0, &thread_offdata ); + intersect = off_intersect_all (&t0, &t1, NULL, NULL, x, y, z, vx, vy, vz, 0, 0, 0, &thread_offdata); if (Vars.Flag_mantid) { - if(intersect) { - Vars.OFF_polyidx=thread_offdata.nextintersect; + if (intersect) { + Vars.OFF_polyidx = thread_offdata.nextintersect; } else { - Vars.OFF_polyidx=-1; + Vars.OFF_polyidx = -1; } } - } - else - if ( (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) - || (abs(Vars.Flag_Shape) == DEFS.SHAPE_DISK) ) /* square xy or disk xy */ + } else if ((abs (Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_DISK)) /* square xy or disk xy */ { // propagate to xy plane and find intersection // make sure the event is recoverable afterwards t0 = t; ALLOW_BACKPROP; PROP_Z0; - if ( (t>=t0) && (z==0.0) ) // forward propagation to xy plane was successful + if ((t >= t0) && (z == 0.0)) // forward propagation to xy plane was successful { - if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) - { + if (abs (Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) { // square xy - intersect = (x>=Vars.mxmin && x<=Vars.mxmax && y>=Vars.mymin && y<=Vars.mymax); - } - else - { + intersect = (x >= Vars.mxmin && x <= Vars.mxmax && y >= Vars.mymin && y <= Vars.mymax); + } else { // disk xy - intersect = (SQR(x) + SQR(y)) <= SQR(Vars.Sphere_Radius); + intersect = (SQR (x) + SQR (y)) <= SQR (Vars.Sphere_Radius); } + } else { + intersect = 0; } - else - { - intersect=0; - } - } - else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) /* sphere */ + } else if (abs (Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) /* sphere */ { - intersect = sphere_intersect(&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius); - /* intersect = (intersect && t0 > 0); */ - } - else if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA)) /* cylinder */ + intersect = sphere_intersect (&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius); + /* intersect = (intersect && t0 > 0); */ + } else if ((abs (Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_BANANA)) /* cylinder */ { - intersect = cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius, Vars.Cylinder_Height); - } - else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX) /* box */ + intersect = cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius, Vars.Cylinder_Height); + } else if (abs (Vars.Flag_Shape) == DEFS.SHAPE_BOX) /* box */ + { + intersect = box_intersect (&t0, &t1, x, y, z, vx, vy, vz, fabs (Vars.mxmax - Vars.mxmin), fabs (Vars.mymax - Vars.mymin), fabs (Vars.mzmax - Vars.mzmin)); + } else if (abs (Vars.Flag_Shape) == DEFS.SHAPE_PREVIOUS) /* previous comp */ { - intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, - fabs(Vars.mxmax-Vars.mxmin), fabs(Vars.mymax-Vars.mymin), fabs(Vars.mzmax-Vars.mzmin)); + intersect = 1; } - else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_PREVIOUS) /* previous comp */ - { intersect = 1; } - if (intersect) - { - if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) - || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA) - || (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) ) - { + if (intersect) { + if ((abs (Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_BOX) + || (abs (Vars.Flag_Shape) == DEFS.SHAPE_BANANA) || (geometry && strlen (geometry) && strcmp (geometry, "0") && strcmp (geometry, "NULL"))) { /* check if we have to remove the top/bottom with BANANA shape */ - if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA) { - if (intersect == 1) { // Entered and left through sides - if (t0 < 0 && t1 > 0) { - t0 = t; /* neutron was already inside ! */ - } - if (t1 < 0 && t0 > 0) { /* neutron exit before entering !! */ - t1 = t; - } - /* t0 is now time of incoming intersection with the detection area */ - if ((Vars.Flag_Shape < 0) && (t1 > 0)) { - PROP_DT(t1); /* t1 outgoing beam */ - } else { - PROP_DT(t0); /* t0 incoming beam */ - } - } else if (intersect == 3 || intersect == 5) { // Entered from top or bottom, left through side - if ((Vars.Flag_Shape < 0) && (t1 > 0)) { - PROP_DT(t1); /* t1 outgoing beam */ - } else { - intersect=0; - Flag_Restore=1; - } - } else if (intersect == 9 || intersect == 17) { // Entered through side, left from top or bottom - if ((Vars.Flag_Shape < 0) && (t1 > 0)) { - intersect=0; - Flag_Restore=1; - } else { - PROP_DT(t0); /* t0 incoming beam */ - } - } else if (intersect == 13 || intersect == 19) { // Went through top/bottom on entry and exit - intersect=0; - Flag_Restore=1; - } else { - printf("Cylinder_intersect returned unexpected value %i\n", intersect); - } + if (abs (Vars.Flag_Shape) == DEFS.SHAPE_BANANA) { + if (intersect == 1) { // Entered and left through sides + if (t0 < 0 && t1 > 0) { + t0 = t; /* neutron was already inside ! */ + } + if (t1 < 0 && t0 > 0) { /* neutron exit before entering !! */ + t1 = t; + } + /* t0 is now time of incoming intersection with the detection area */ + if ((Vars.Flag_Shape < 0) && (t1 > 0)) { + PROP_DT (t1); /* t1 outgoing beam */ + } else { + PROP_DT (t0); /* t0 incoming beam */ + } + } else if (intersect == 3 || intersect == 5) { // Entered from top or bottom, left through side + if ((Vars.Flag_Shape < 0) && (t1 > 0)) { + PROP_DT (t1); /* t1 outgoing beam */ + } else { + intersect = 0; + Flag_Restore = 1; + } + } else if (intersect == 9 || intersect == 17) { // Entered through side, left from top or bottom + if ((Vars.Flag_Shape < 0) && (t1 > 0)) { + intersect = 0; + Flag_Restore = 1; + } else { + PROP_DT (t0); /* t0 incoming beam */ + } + } else if (intersect == 13 || intersect == 19) { // Went through top/bottom on entry and exit + intersect = 0; + Flag_Restore = 1; } else { - // All other shapes than the BANANA - if (t0 < 0 && t1 > 0) - t0 = t; /* neutron was already inside ! */ - if (t1 < 0 && t0 > 0) /* neutron exit before entering !! */ - t1 = t; - /* t0 is now time of incoming intersection with the detection area */ - if ((Vars.Flag_Shape < 0) && (t1 > 0)) - PROP_DT(t1); /* t1 outgoing beam */ - else - PROP_DT(t0); /* t0 incoming beam */ + printf ("Cylinder_intersect returned unexpected value %i\n", intersect); } - + } else { + // All other shapes than the BANANA + if (t0 < 0 && t1 > 0) + t0 = t; /* neutron was already inside ! */ + if (t1 < 0 && t0 > 0) /* neutron exit before entering !! */ + t1 = t; + /* t0 is now time of incoming intersection with the detection area */ + if ((Vars.Flag_Shape < 0) && (t1 > 0)) + PROP_DT (t1); /* t1 outgoing beam */ + else + PROP_DT (t0); /* t0 incoming beam */ + } + /* Final test if we are on lid / bottom of banana/sphere */ - if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA || abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) { - if (Vars.Cylinder_Height && fabs(y) >= Vars.Cylinder_Height/2 - FLT_EPSILON) { - intersect=0; - Flag_Restore=1; + if (abs (Vars.Flag_Shape) == DEFS.SHAPE_BANANA || abs (Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) { + if (Vars.Cylinder_Height && fabs (y) >= Vars.Cylinder_Height / 2 - FLT_EPSILON) { + intersect = 0; + Flag_Restore = 1; } } } } - if (intersect) - { + if (intersect) { /* Now get the data to monitor: current or keep from PreMonitor */ -/* if (Vars.Flag_UsePreMonitor != 1)*/ -/* {*/ -/* Vars.cp = p;*/ -/* Vars.cx = x;*/ -/* Vars.cvx = vx;*/ -/* Vars.csx = sx;*/ -/* Vars.cy = y;*/ -/* Vars.cvy = vy;*/ -/* Vars.csy = sy;*/ -/* Vars.cz = z;*/ -/* Vars.cvz = vz;*/ -/* Vars.csz = sz;*/ -/* Vars.ct = t;*/ -/* }*/ - - if ((Vars.He3_pressure > 0) && (t1 != t0) && ((abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX))) - { - transmit_he3 = exp(-7.417*Vars.He3_pressure*fabs(t1-t0)*2*PI*K2V); + /* if (Vars.Flag_UsePreMonitor != 1)*/ + /* {*/ + /* Vars.cp = p;*/ + /* Vars.cx = x;*/ + /* Vars.cvx = vx;*/ + /* Vars.csx = sx;*/ + /* Vars.cy = y;*/ + /* Vars.cvy = vy;*/ + /* Vars.csy = sy;*/ + /* Vars.cz = z;*/ + /* Vars.cvz = vz;*/ + /* Vars.csz = sz;*/ + /* Vars.ct = t;*/ + /* }*/ + + if ((Vars.He3_pressure > 0) && (t1 != t0) + && ((abs (Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs (Vars.Flag_Shape) == DEFS.SHAPE_BOX))) { + transmit_he3 = exp (-7.417 * Vars.He3_pressure * fabs (t1 - t0) * 2 * PI * K2V); /* will monitor the absorbed part */ - p = p * (1-transmit_he3); + p = p * (1 - transmit_he3); } - if (Vars.Flag_capture) - { - multiplier_capture = V2K*sqrt(vx*vx+vy*vy+vz*vz); - if (multiplier_capture != 0) multiplier_capture = 2*PI/multiplier_capture; /* lambda. lambda(2200 m/2) = 1.7985 Angs */ - p = p * multiplier_capture/1.7985; + if (Vars.Flag_capture) { + multiplier_capture = V2K * sqrt (vx * vx + vy * vy + vz * vz); + if (multiplier_capture != 0) + multiplier_capture = 2 * PI / multiplier_capture; /* lambda. lambda(2200 m/2) = 1.7985 Angs */ + p = p * multiplier_capture / 1.7985; } - pp = Monitor_nd_noaccTrace(&DEFS, &Vars, _particle); - if (pp==0.0) - { + pp = Monitor_nd_noaccTrace (&DEFS, &Vars, _particle); + if (pp == 0.0) { ABSORB; - } - else if(pp==1) - { + } else if (pp == 1) { SCATTER; } /*set weight to undetected part if capture and/or he3_pressure*/ - if (Vars.He3_pressure > 0){ + if (Vars.He3_pressure > 0) { /* after monitor, only remains 1-p_detect */ - p = p * transmit_he3/(1.0-transmit_he3); + p = p * transmit_he3 / (1.0 - transmit_he3); } - if (Vars.Flag_capture){ - p = p / multiplier_capture*1.7985; + if (Vars.Flag_capture) { + p = p / multiplier_capture * 1.7985; } if (Vars.Flag_parallel) /* back to neutron state before detection */ Flag_Restore = 1; } /* end if intersection */ else { - if (Vars.Flag_Absorb && !Vars.Flag_parallel) - { + if (Vars.Flag_Absorb && !Vars.Flag_parallel) { // restore neutron ray before absorbing for correct mcdisplay - 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); ABSORB; - } - else Flag_Restore = 1; /* no intersection, back to previous state */ + } else + Flag_Restore = 1; /* no intersection, back to previous state */ } - if (Flag_Restore) - { - RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); + if (Flag_Restore) { + RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); } - + #ifdef OPENACC_BAK #define OPENACC #define fprintf(stderr,...) printf(__VA_ARGS__) @@ -617,26 +602,25 @@ TRACE SAVE %{ -if (!nowritefile) { - - /* save results, but do not free pointers */ - detector = Monitor_nd_noaccSave(&DEFS, &Vars); -} + if (!nowritefile) { + + /* save results, but do not free pointers */ + detector = Monitor_nd_noaccSave (&DEFS, &Vars); + } %} FINALLY %{ /* free pointers */ - Monitor_nd_noaccFinally(&DEFS, &Vars); + Monitor_nd_noaccFinally (&DEFS, &Vars); %} MCDISPLAY %{ - if (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) - { - off_display(offdata); + if (geometry && strlen (geometry) && strcmp (geometry, "0") && strcmp (geometry, "NULL")) { + off_display (offdata); } else { - Monitor_nd_noaccMcDisplay(&DEFS, &Vars); + Monitor_nd_noaccMcDisplay (&DEFS, &Vars); } %} diff --git a/mcstas-comps/monitors/PSD_TOF_monitor.comp b/mcstas-comps/monitors/PSD_TOF_monitor.comp index 23c95b017..1feb294ab 100644 --- a/mcstas-comps/monitors/PSD_TOF_monitor.comp +++ b/mcstas-comps/monitors/PSD_TOF_monitor.comp @@ -68,32 +68,38 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("PSD_TOF_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("PSD_TOF_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } if (tmin >= tmax) { - printf("PSD_TOF_monitor: %s: Unmeaningful TOF definition!\n", - NAME_CURRENT_COMP); - exit(0); + printf ("PSD_TOF_monitor: %s: Unmeaningful TOF definition!\n", NAME_CURRENT_COMP); + exit (0); } - PSD_N = create_darr3d(nt, nx, ny); - PSD_p = create_darr3d(nt, nx, ny); - PSD_p2 = create_darr3d(nt, nx, ny); + PSD_N = create_darr3d (nt, nx, ny); + PSD_p = create_darr3d (nt, nx, ny); + PSD_p2 = create_darr3d (nt, nx, ny); - PSD_Ns = create_darr2d(nx, ny); - PSD_ps = create_darr2d(nx, ny); - PSD_p2s = create_darr2d(nx, ny); + PSD_Ns = create_darr2d (nx, ny); + PSD_ps = create_darr2d (nx, ny); + PSD_p2s = create_darr2d (nx, ny); // 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 @@ -101,83 +107,67 @@ TRACE int i, j, k; PROP_Z0; - if (x>xmin && xymin && ytmin && (1E6*t) xmin && x < xmax && y > ymin && y < ymax && (1E6 * t) > tmin && (1E6 * t) < tmax) { + i = floor ((x - xmin) * nx / (xmax - xmin)); + j = floor ((y - ymin) * ny / (ymax - ymin)); + k = floor ((1E6 * t - tmin) * nt / (tmax - tmin)); + + double p2 = p * p; #pragma acc atomic - PSD_p[k][i][j] = PSD_p[k][i][j]+p; - #pragma acc atomic - PSD_p2[k][i][j] = PSD_p2[k][i][j]+p2; + PSD_N[k][i][j] = PSD_N[k][i][j] + 1; + #pragma acc atomic + PSD_p[k][i][j] = PSD_p[k][i][j] + p; + #pragma acc atomic + PSD_p2[k][i][j] = PSD_p2[k][i][j] + p2; #pragma acc atomic - PSD_Ns[i][j] = PSD_Ns[i][j]+1; + PSD_Ns[i][j] = PSD_Ns[i][j] + 1; + #pragma acc atomic + PSD_ps[i][j] = PSD_ps[i][j] + p; #pragma acc atomic - PSD_ps[i][j] = PSD_ps[i][j]+p; - #pragma acc atomic - PSD_p2s[i][j] = PSD_p2s[i][j]+p2; + PSD_p2s[i][j] = PSD_p2s[i][j] + p2; SCATTER; } 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) { - int kk; - char ff[256]; - char tt[256]; - sprintf(ff, "%s_Sum",filename); - DETECTOR_OUT_2D( - "PSD monitor TOF Sum", - "X position [cm]", - "Y position [cm]", - xmin*100.0, xmax*100.0, ymin*100.0, ymax*100.0, - nx, ny, - &PSD_Ns[0][0],&PSD_ps[0][0],&PSD_p2s[0][0], - ff); - - for (kk=0; kk 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } - - if ((xmin >= xmax) || (ymin >= ymax)){ - printf("PSD_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } + + if ((xmin >= xmax) || (ymin >= ymax)) { + printf ("PSD_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - PSD_N = create_darr2d(nx, ny); - PSD_p = create_darr2d(nx, ny); - PSD_p2 = create_darr2d(nx, ny); + PSD_N = create_darr2d (nx, ny); + PSD_p = create_darr2d (nx, ny); + PSD_p2 = create_darr2d (nx, ny); // 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 %{ PROP_Z0; - if (x>xmin && xymin && y xmin && x < xmax && y > ymin && y < ymax) { + int i = floor ((x - xmin) * nx / (xmax - xmin)); + int j = floor ((y - ymin) * ny / (ymax - ymin)); - double p2 = p*p; + double p2 = p * p; #pragma acc atomic - PSD_N[i][j] = PSD_N[i][j]+1; + PSD_N[i][j] = PSD_N[i][j] + 1; #pragma acc atomic - PSD_p[i][j] = PSD_p[i][j]+p; - + PSD_p[i][j] = PSD_p[i][j] + p; + #pragma acc atomic PSD_p2[i][j] = PSD_p2[i][j] + p2; - + SCATTER; } 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_2D( - "PSD monitor", - "X position [cm]", - "Y position [cm]", - xmin*100.0, xmax*100.0, ymin*100.0, ymax*100.0, - nx, ny, - &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0], - filename); - } + if (!nowritefile) { + DETECTOR_OUT_2D ("PSD monitor", "X position [cm]", "Y position [cm]", xmin * 100.0, xmax * 100.0, ymin * 100.0, ymax * 100.0, nx, ny, &PSD_N[0][0], + &PSD_p[0][0], &PSD_p2[0][0], filename); + } %} FINALLY %{ @@ -124,12 +125,9 @@ FINALLY %{ 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 diff --git a/mcstas-comps/monitors/PSD_monitor_4PI.comp b/mcstas-comps/monitors/PSD_monitor_4PI.comp index 290955a76..a59dec020 100644 --- a/mcstas-comps/monitors/PSD_monitor_4PI.comp +++ b/mcstas-comps/monitors/PSD_monitor_4PI.comp @@ -61,81 +61,74 @@ DECLARE INITIALIZE %{ - PSD_N = create_darr2d(nx, ny); - PSD_p = create_darr2d(nx, ny); - PSD_p2 = create_darr2d(nx, ny); + PSD_N = create_darr2d (nx, ny); + PSD_p = create_darr2d (nx, ny); + PSD_p2 = create_darr2d (nx, ny); // 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 %{ double t0, t1, phi, theta; - int i,j; + int i, j; - if(sphere_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius) && t1 > 0) - { - if(t0 < 0) + if (sphere_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius) && t1 > 0) { + if (t0 < 0) t0 = t1; /* t0 is now time of intersection with the sphere. */ - mcPROP_DT(t0); - phi = atan2(x,z); - i = floor(nx*(phi/(2*PI)+0.5)); - if(i >= nx) - i = nx-1; /* Special case for phi = PI. */ - else if(i < 0) + mcPROP_DT (t0); + phi = atan2 (x, z); + i = floor (nx * (phi / (2 * PI) + 0.5)); + if (i >= nx) + i = nx - 1; /* Special case for phi = PI. */ + else if (i < 0) i = 0; - theta=asin(y/radius); - j = floor(ny*(theta+PI/2)/PI+0.5); - if(j >= ny) - j = ny-1; /* Special case for y = radius. */ - else if(j < 0) + theta = asin (y / radius); + j = floor (ny * (theta + PI / 2) / PI + 0.5); + if (j >= ny) + j = ny - 1; /* Special case for y = radius. */ + else if (j < 0) j = 0; - double p2 = p*p; + double p2 = p * p; #pragma acc atomic - PSD_N[i][j] = PSD_N[i][j]+1; - + PSD_N[i][j] = PSD_N[i][j] + 1; + #pragma acc atomic - PSD_p[i][j] = PSD_p[i][j]+p; - + PSD_p[i][j] = PSD_p[i][j] + p; + #pragma acc atomic PSD_p2[i][j] = PSD_p2[i][j] + p2; - + SCATTER; } 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_2D( - "4PI PSD monitor", - "Longitude [deg]", - "Latitude [deg]", - -180, 180, -90, 90, - nx, ny, - &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("4PI PSD monitor", "Longitude [deg]", "Latitude [deg]", -180, 180, -90, 90, nx, ny, &PSD_N[0][0], &PSD_p[0][0], &PSD_p2[0][0], filename); + } %} FINALLY %{ - destroy_darr2d(PSD_N); - destroy_darr2d(PSD_p); - destroy_darr2d(PSD_p2); + destroy_darr2d (PSD_N); + destroy_darr2d (PSD_p); + destroy_darr2d (PSD_p2); %} MCDISPLAY %{ - circle("xy",0,0,0,radius); - circle("xz",0,0,0,radius); - circle("yz",0,0,0,radius); + circle ("xy", 0, 0, 0, radius); + circle ("xz", 0, 0, 0, radius); + circle ("yz", 0, 0, 0, radius); %} END diff --git a/mcstas-comps/monitors/PSD_monitor_4PI_spin.comp b/mcstas-comps/monitors/PSD_monitor_4PI_spin.comp index 1da4273eb..bcca4e7bc 100644 --- a/mcstas-comps/monitors/PSD_monitor_4PI_spin.comp +++ b/mcstas-comps/monitors/PSD_monitor_4PI_spin.comp @@ -51,44 +51,44 @@ DECLARE %} INITIALIZE %{ - PSD_N = create_darr2d(nx, ny); - PSDu_p = create_darr2d(nx, ny); - PSDd_p = create_darr2d(nx, ny); - PSDu_p2 = create_darr2d(nx, ny); - PSDd_p2 = create_darr2d(nx, ny); + PSD_N = create_darr2d (nx, ny); + PSDu_p = create_darr2d (nx, ny); + PSDd_p = create_darr2d (nx, ny); + PSDu_p2 = create_darr2d (nx, ny); + PSDd_p2 = create_darr2d (nx, ny); // 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 %{ double t0, t1, phi, theta; - int i,j; + int i, j; - if(sphere_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius) && t1 > 0) - { - if(t0 < 0) + if (sphere_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius) && t1 > 0) { + if (t0 < 0) t0 = t1; /* t0 is now time of intersection with the sphere. */ - PROP_DT(t0); - phi = atan2(x,z); - i = floor(nx*(phi/(2*PI)+0.5)); - if(i >= nx) - i = nx-1; /* Special case for phi = PI. */ - else if(i < 0) + PROP_DT (t0); + phi = atan2 (x, z); + i = floor (nx * (phi / (2 * PI) + 0.5)); + if (i >= nx) + i = nx - 1; /* Special case for phi = PI. */ + else if (i < 0) i = 0; - theta=asin(y/radius); - j = floor(ny*(theta+PI/2)/PI+0.5); - if(j >= ny) - j = ny-1; /* Special case for y = radius. */ - else if(j < 0) + theta = asin (y / radius); + j = floor (ny * (theta + PI / 2) / PI + 0.5); + if (j >= ny) + j = ny - 1; /* Special case for y = radius. */ + else if (j < 0) j = 0; /*look at spin-down neutrons*/ - double n_up=(1+scalar_prod(sx,sy,sz,mx,my,mz))*p/2; + double n_up = (1 + scalar_prod (sx, sy, sz, mx, my, mz)) * p / 2; /*look at spin-down neutrons*/ - double n_down=(1-scalar_prod(sx,sy,sz,mx,my,mz))*p/2; + double n_down = (1 - scalar_prod (sx, sy, sz, mx, my, mz)) * p / 2; - double n_up2 = n_up*n_up; - double n_down2 = n_down*n_down; + double n_up2 = n_up * n_up; + double n_down2 = n_down * n_down; #pragma acc atomic PSD_N[i][j] = PSD_N[i][j] + 1; #pragma acc atomic @@ -102,35 +102,20 @@ TRACE SCATTER; } 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 %{ - char title[256]; - char fn[256]; - snprintf(title,256,"4PI PSD Monitor spin %s (parallel to m=%g,%g,%g).","up",mx,my,mz); - snprintf(fn,256,"%s.up",filename); - DETECTOR_OUT_2D( - title, - "Longitude [deg]", - "Lattitude [deg]", - -180, 180, -90, 90, - nx, ny, - &PSD_N[0][0],&PSDu_p[0][0],&PSDu_p2[0][0], - fn); - snprintf(title,256,"4PI PSD Monitor spin %s (antiparallel to m=%g,%g,%g).","down",mx,my,mz); - snprintf(fn,256,"%s.down",filename); - DETECTOR_OUT_2D( - title, - "Longitude [deg]", - "Lattitude [deg]", - -180, 180, -90, 90, - nx, ny, - &PSD_N[0][0],&PSDd_p[0][0],&PSDd_p2[0][0], - fn); - + char title[256]; + char fn[256]; + snprintf (title, 256, "4PI PSD Monitor spin %s (parallel to m=%g,%g,%g).", "up", mx, my, mz); + snprintf (fn, 256, "%s.up", filename); + DETECTOR_OUT_2D (title, "Longitude [deg]", "Lattitude [deg]", -180, 180, -90, 90, nx, ny, &PSD_N[0][0], &PSDu_p[0][0], &PSDu_p2[0][0], fn); + snprintf (title, 256, "4PI PSD Monitor spin %s (antiparallel to m=%g,%g,%g).", "down", mx, my, mz); + snprintf (fn, 256, "%s.down", filename); + DETECTOR_OUT_2D (title, "Longitude [deg]", "Lattitude [deg]", -180, 180, -90, 90, nx, ny, &PSD_N[0][0], &PSDd_p[0][0], &PSDd_p2[0][0], fn); %} FINALLY %{ @@ -143,10 +128,10 @@ FINALLY %{ MCDISPLAY %{ - magnify(""); - circle("xy",0,0,0,radius); - circle("xz",0,0,0,radius); - circle("yz",0,0,0,radius); + magnify (""); + circle ("xy", 0, 0, 0, radius); + circle ("xz", 0, 0, 0, radius); + circle ("yz", 0, 0, 0, radius); %} END diff --git a/mcstas-comps/monitors/PSD_monitor_TOF.comp b/mcstas-comps/monitors/PSD_monitor_TOF.comp index af916d651..69cf5c976 100644 --- a/mcstas-comps/monitors/PSD_monitor_TOF.comp +++ b/mcstas-comps/monitors/PSD_monitor_TOF.comp @@ -66,115 +66,105 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("PSD_TOF_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("PSD_TOF_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } if (tmin >= tmax) { - printf("PSD_TOF_monitor: %s: Unmeaningful TOF definition!\n", - NAME_CURRENT_COMP); - exit(0); + printf ("PSD_TOF_monitor: %s: Unmeaningful TOF definition!\n", NAME_CURRENT_COMP); + exit (0); } - PSD_N = create_darr3d(nx, ny, nt); - PSD_p = create_darr3d(nx, ny, nt); - PSD_p2 = create_darr3d(nx, ny, nt); - PSD_Ns = create_darr2d(nx, ny); - PSD_ps = create_darr2d(nx, ny); - PSD_p2s = create_darr2d(nx, ny); + PSD_N = create_darr3d (nx, ny, nt); + PSD_p = create_darr3d (nx, ny, nt); + PSD_p2 = create_darr3d (nx, ny, nt); + PSD_Ns = create_darr2d (nx, ny); + PSD_ps = create_darr2d (nx, ny); + PSD_p2s = create_darr2d (nx, ny); // 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,k; + int i, j, k; PROP_Z0; - if (x>xmin && xymin && ytmin && (1E6*t) xmin && x < xmax && y > ymin && y < ymax && (1E6 * t) > tmin && (1E6 * t) < tmax) { + i = floor ((x - xmin) * nx / (xmax - xmin)); + j = floor ((y - ymin) * ny / (ymax - ymin)); + k = floor ((1E6 * t - tmin) * nt / (tmax - tmin)); + + double p2 = p * p; + #pragma acc atomic + PSD_N[i][j][k] = PSD_N[i][j][k] + 1; + #pragma acc atomic + PSD_p[i][j][k] = PSD_p[i][j][k] + p; + #pragma acc atomic + PSD_p2[i][j][k] = PSD_p2[i][j][k] + p2; + + #pragma acc atomic + PSD_Ns[i][j] = PSD_Ns[i][j] + 1; + #pragma acc atomic + PSD_ps[i][j] = PSD_ps[i][j] + p; + #pragma acc atomic + PSD_p2s[i][j] = PSD_p2s[i][j] + p2; SCATTER; } 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) { - int kk,ll; - char ff[256]; - char tt[256]; - sprintf(ff, "%s_Sum",filename); - DETECTOR_OUT_2D( - "PSD monitor TOF Sum", - "X position [cm]", - "Y position [cm]", - xmin*100.0, xmax*100.0, ymin*100.0, ymax*100.0, - nx, ny, - &PSD_Ns[0][0],&PSD_ps[0][0],&PSD_p2s[0][0], - ff); - - for (kk=0; kk 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("PSD_monitor_psf: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("PSD_monitor_psf: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - PSD_N = create_darr2d(nx, ny); - PSD_p = create_darr2d(nx, ny); - PSD_p2 = create_darr2d(nx, ny); + PSD_N = create_darr2d (nx, ny); + PSD_p = create_darr2d (nx, ny); + PSD_p2 = create_darr2d (nx, ny); // 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 xp,yp; + int i, j; + double xp, yp; PROP_Z0; - xp = x + psf*randnorm(); - yp = y + psf*randnorm(); + xp = x + psf * randnorm (); + yp = y + psf * randnorm (); - if (xp>xmin && xpymin && yp xmin && xp < xmax && yp > ymin && yp < ymax) { + i = floor ((xp - xmin) * nx / (xmax - xmin)); + j = floor ((yp - ymin) * ny / (ymax - ymin)); - double p2 = p*p; + double p2 = p * p; #pragma acc atomic - PSD_N[i][j] = PSD_N[i][j]+1; + PSD_N[i][j] = PSD_N[i][j] + 1; #pragma acc atomic - PSD_p[i][j] = PSD_p[i][j]+p; + PSD_p[i][j] = PSD_p[i][j] + p; #pragma acc atomic - PSD_p2[i][j] = PSD_p2[i][j]+p2; - + PSD_p2[i][j] = PSD_p2[i][j] + p2; + SCATTER; } 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_2D( - "PSD monitor", - "X position [cm]", - "Y position [cm]", - xmin*100.0, xmax*100.0, ymin*100.0, ymax*100.0, - nx, ny, - &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("PSD monitor", "X position [cm]", "Y position [cm]", xmin * 100.0, xmax * 100.0, ymin * 100.0, ymax * 100.0, nx, ny, &PSD_N[0][0], + &PSD_p[0][0], &PSD_p2[0][0], filename); + } %} FINALLY %{ - destroy_darr2d(PSD_N); - destroy_darr2d(PSD_p); - destroy_darr2d(PSD_p2); + destroy_darr2d (PSD_N); + destroy_darr2d (PSD_p); + destroy_darr2d (PSD_p2); %} 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 diff --git a/mcstas-comps/monitors/PSD_monitor_psf_eff.comp b/mcstas-comps/monitors/PSD_monitor_psf_eff.comp index 56e171797..985e2908a 100644 --- a/mcstas-comps/monitors/PSD_monitor_psf_eff.comp +++ b/mcstas-comps/monitors/PSD_monitor_psf_eff.comp @@ -67,87 +67,84 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("PSD_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("PSD_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - PSD_N = create_darr2d(nx, ny); - PSD_p = create_darr2d(nx, ny); - PSD_p2 = create_darr2d(nx, ny); + PSD_N = create_darr2d (nx, ny); + PSD_p = create_darr2d (nx, ny); + PSD_p2 = create_darr2d (nx, ny); // 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 xp,yp; + int i, j; + double xp, yp; PROP_Z0; - xp = x + psf*randnorm(); - yp = y + psf*randnorm(); + xp = x + psf * randnorm (); + yp = y + psf * randnorm (); + + if (xp > xmin && xp < xmax && yp > ymin && yp < ymax) { + weight = eff * k0 / (V2K * sqrt (vx * vx + vy * vy + vz * vz)); + i = floor ((xp - xmin) * nx / (xmax - xmin)); + j = floor ((yp - ymin) * ny / (ymax - ymin)); - if (xp>xmin && xpymin && yp0) { - PROP_DT(t1); - /* Calculate pixel */ - if (fabs(y) < yheight/2.0) { - double phi, yprime; - phi=(atan2(x,z)+PI)*RAD2DEG; - yprime=y+yheight/2.0; - - if (phi >= thmin && phi < thmax) { - i=floor((nr-1) * (phi-thmin)/(thmax-thmin)); - j=floor((ny-1) * yprime/yheight); - - 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) { + double phi, yprime; + phi = (atan2 (x, z) + PI) * RAD2DEG; + yprime = y + yheight / 2.0; + + if (phi >= thmin && phi < thmax) { + i = floor ((nr - 1) * (phi - thmin) / (thmax - thmin)); + j = floor ((ny - 1) * yprime / yheight); + + double p2 = p * p; + #pragma acc atomic + PSD_N[i][j] = PSD_N[i][j] + 1; + #pragma acc atomic - PSD_N[i][j] = PSD_N[i][j]+1; + PSD_p[i][j] = PSD_p[i][j] + p; #pragma acc atomic - PSD_p[i][j] = PSD_p[i][j]+p; - - #pragma acc atomic - PSD_p2[i][j] = PSD_p2[i][j]+p2; - } - } + PSD_p2[i][j] = PSD_p2[i][j] + 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 %{ - DETECTOR_OUT_2D( - "PSD cylindrical monitor", - "radial position [deg]", - "yheight [m]", - thmin, thmax, -yheight/2.0, yheight/2.0, - nr, ny, - &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0], - filename); + DETECTOR_OUT_2D ("PSD cylindrical monitor", "radial position [deg]", "yheight [m]", thmin, thmax, -yheight / 2.0, yheight / 2.0, nr, ny, &PSD_N[0][0], + &PSD_p[0][0], &PSD_p2[0][0], filename); %} FINALLY %{ - destroy_darr2d(PSD_N); - destroy_darr2d(PSD_p); - destroy_darr2d(PSD_p2); + destroy_darr2d (PSD_N); + destroy_darr2d (PSD_p); + destroy_darr2d (PSD_p2); %} MCDISPLAY %{ - circle("xz", 0, 0, 0, radius ); + circle ("xz", 0, 0, 0, radius); %} END diff --git a/mcstas-comps/monitors/PSDlin_diff_monitor.comp b/mcstas-comps/monitors/PSDlin_diff_monitor.comp index 57a6a71de..d17df3afc 100644 --- a/mcstas-comps/monitors/PSDlin_diff_monitor.comp +++ b/mcstas-comps/monitors/PSDlin_diff_monitor.comp @@ -60,27 +60,34 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("PSDlin_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("PSDlin_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - PSDlin_N = create_darr1d(nx); - PSDlin_p = create_darr1d(nx); - PSDlin_p2 = create_darr1d(nx); - PSDdiff_N = create_darr1d(nx-1); - PSDdiff_p = create_darr1d(nx-1); - PSDdiff_p2 = create_darr1d(nx-1); + PSDlin_N = create_darr1d (nx); + PSDlin_p = create_darr1d (nx); + PSDlin_p2 = create_darr1d (nx); + PSDdiff_N = create_darr1d (nx - 1); + PSDdiff_p = create_darr1d (nx - 1); + PSDdiff_p2 = create_darr1d (nx - 1); // Use instance name for monitor output if no input was given - if (!strcmp(filename,"\0")) sprintf(filename,"%s",NAME_CURRENT_COMP); - - sprintf(difilename,"%s_diff",filename); + if (!strcmp (filename, "\0")) + sprintf (filename, "%s", NAME_CURRENT_COMP); + + sprintf (difilename, "%s_diff", filename); %} TRACE @@ -88,78 +95,63 @@ TRACE int i; PROP_Z0; - if (x>xmin && xymin && y xmin && x < xmax && y > ymin && y < ymax) { + i = floor (nx * (x - xmin) / (xmax - xmin)); /* Bin number */ - double p2 = p*p; - #pragma acc atomic - PSDlin_N[i] = PSDlin_N[i]+1; + double p2 = p * p; + #pragma acc atomic + PSDlin_N[i] = PSDlin_N[i] + 1; - #pragma acc atomic - PSDlin_p[i] = PSDlin_p[i]+p; + #pragma acc atomic + PSDlin_p[i] = PSDlin_p[i] + p; - #pragma acc atomic - PSDlin_p2[i] = PSDlin_p2[i]+p2; + #pragma acc atomic + PSDlin_p2[i] = PSDlin_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 %{ -int i; -if (!nowritefile) { - double mean; - for (i=0; i 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("PSDlin_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("PSDlin_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - PSDlin_N = create_darr1d(nbins); - PSDlin_p = create_darr1d(nbins); - PSDlin_p2 = create_darr1d(nbins); + PSDlin_N = create_darr1d (nbins); + PSDlin_p = create_darr1d (nbins); + PSDlin_p2 = create_darr1d (nbins); // 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 @@ -82,21 +89,20 @@ TRACE int i; PROP_Z0; - if (x>xmin && xymin && y xmin && x < xmax && y > ymin && y < ymax) { /* Bin number */ - if (!vertical){ - i = floor(nbins*(x-xmin)/(xmax-xmin)); - }else{ - i = floor(nbins*(y-ymin)/(ymax-ymin)); + if (!vertical) { + i = floor (nbins * (x - xmin) / (xmax - xmin)); + } else { + i = floor (nbins * (y - ymin) / (ymax - ymin)); } - if((i >= nbins) || (i<0)) - { - printf("ERROR: (%s) wrong positioning in linear PSD. i= %i \n",NAME_CURRENT_COMP,i); - exit(1); + if ((i >= nbins) || (i < 0)) { + printf ("ERROR: (%s) wrong positioning in linear PSD. i= %i \n", NAME_CURRENT_COMP, i); + exit (1); } - double p2 = p*p; + double p2 = p * p; #pragma acc atomic - PSDlin_N[i] = PSDlin_N[i] +1; + PSDlin_N[i] = PSDlin_N[i] + 1; #pragma acc atomic PSDlin_p[i] = PSDlin_p[i] + p; #pragma acc atomic @@ -104,39 +110,32 @@ TRACE SCATTER; } 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) { - if (!vertical){ - DETECTOR_OUT_1D( - "Linear PSD monitor","x-Position [m]","Intensity","x", xmin, xmax, nbins, - &PSDlin_N[0],&PSDlin_p[0],&PSDlin_p2[0],filename); - }else{ - DETECTOR_OUT_1D( - "Linear PSD monitor","y-Position [m]","Intensity","y", ymin, ymax, nbins, - &PSDlin_N[0],&PSDlin_p[0],&PSDlin_p2[0],filename); + if (!nowritefile) { + if (!vertical) { + DETECTOR_OUT_1D ("Linear PSD monitor", "x-Position [m]", "Intensity", "x", xmin, xmax, nbins, &PSDlin_N[0], &PSDlin_p[0], &PSDlin_p2[0], filename); + } else { + DETECTOR_OUT_1D ("Linear PSD monitor", "y-Position [m]", "Intensity", "y", ymin, ymax, nbins, &PSDlin_N[0], &PSDlin_p[0], &PSDlin_p2[0], filename); + } } -} %} FINALLY %{ - destroy_darr1d(PSDlin_N); - destroy_darr1d(PSDlin_p); - destroy_darr1d(PSDlin_p2); + destroy_darr1d (PSDlin_N); + destroy_darr1d (PSDlin_p); + destroy_darr1d (PSDlin_p2); %} 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 diff --git a/mcstas-comps/monitors/PolLambda_monitor.comp b/mcstas-comps/monitors/PolLambda_monitor.comp index eb3b04d39..3e41a5939 100644 --- a/mcstas-comps/monitors/PolLambda_monitor.comp +++ b/mcstas-comps/monitors/PolLambda_monitor.comp @@ -62,36 +62,40 @@ DECLARE INITIALIZE %{ // Check that input parameteters makes sense - if (Lmax<=Lmin) { - fprintf(stderr, "Pol_monitor: %s: l1 <= l0!\n" - "ERROR. Exiting", - NAME_CURRENT_COMP); - exit(1); + if (Lmax <= Lmin) { + fprintf (stderr, + "Pol_monitor: %s: l1 <= l0!\n" + "ERROR. Exiting", + NAME_CURRENT_COMP); + exit (1); } - if (mx==0 && my==0 && mz==0) { - fprintf(stderr, "Pol_monitor: %s: NULL vector defined!\n" - "ERROR (mx, my, mz). Exiting", - NAME_CURRENT_COMP); - exit(1); + if (mx == 0 && my == 0 && mz == 0) { + fprintf (stderr, + "Pol_monitor: %s: NULL vector defined!\n" + "ERROR (mx, my, mz). Exiting", + NAME_CURRENT_COMP); + exit (1); } - if ((xwidth<=0) || (yheight <= 0)) { - fprintf(stderr, "Pol_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight). Exiting", - NAME_CURRENT_COMP); - exit(1); + if ((xwidth <= 0) || (yheight <= 0)) { + fprintf (stderr, + "Pol_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight). Exiting", + NAME_CURRENT_COMP); + exit (1); } // Initialize variables - NORM(mx, my, mz); + NORM (mx, my, mz); - PolL_N = create_darr2d(nL, npol); - PolL_p = create_darr2d(nL, npol); - PolL_p2 = create_darr2d(nL, npol); + PolL_N = create_darr2d (nL, npol); + PolL_p = create_darr2d (nL, npol); + PolL_p2 = create_darr2d (nL, npol); // 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 @@ -101,60 +105,55 @@ TRACE double lambda; PROP_Z0; - lambda = (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz); + lambda = (2 * PI / V2K) / sqrt (vx * vx + vy * vy + vz * vz); - if (inside_rectangle(x, y, xwidth, yheight) && - lambda > Lmin && lambda < Lmax) { + if (inside_rectangle (x, y, xwidth, yheight) && lambda > Lmin && lambda < Lmax) { - pol_proj = scalar_prod(mx, my, mz, sx, sy, sz); + pol_proj = scalar_prod (mx, my, mz, sx, sy, sz); /*protect from rounding errors introduced by trig functions*/ - if (fabs(pol_proj-1)1) { - if (fabs(pol_proj)<1+FLT_EPSILON){ - pol_proj=1; - }else{ - ABSORB; - } + if (fabs (pol_proj) > 1) { + if (fabs (pol_proj) < 1 + FLT_EPSILON) { + pol_proj = 1; + } else { + ABSORB; + } } - double p2 = p*p; - double pp = pol_proj*p; - double p2p = pol_proj*pol_proj*p; + double p2 = p * p; + double pp = pol_proj * p; + double p2p = pol_proj * pol_proj * p; #pragma acc atomic - Pol_N = Pol_N+1; + Pol_N = Pol_N + 1; #pragma acc atomic - Psum = Psum+p; + Psum = Psum + p; #pragma acc atomic - Pol_p = Pol_p+pp; + Pol_p = Pol_p + pp; #pragma acc atomic - Pol_p2 = Pol_p2+p2p; - + Pol_p2 = Pol_p2 + p2p; + SCATTER; } 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) { - if (Psum && Pol_N){ - Pol_p /= Psum; -#ifdef USE_MPI + if (!nowritefile) { + if (Psum && Pol_N) { + Pol_p /= Psum; + #ifdef USE_MPI Pol_p /= mpi_node_count; -#endif + #endif Pol_p2 /= Psum; - Pol_p2 -= Pol_p*Pol_p; + Pol_p2 -= Pol_p * Pol_p; Pol_p2 /= Pol_N; - } - DETECTOR_OUT_0D(titlestring, Pol_N, Pol_p, Pol_p2); } + DETECTOR_OUT_0D (titlestring, Pol_N, Pol_p, Pol_p2); + } %} MCDISPLAY %{ - rectangle("xy", 0, 0, 0, xwidth, yheight); + rectangle ("xy", 0, 0, 0, xwidth, yheight); %} END diff --git a/mcstas-comps/monitors/Sqq_w_monitor.comp b/mcstas-comps/monitors/Sqq_w_monitor.comp index 5adfc6e56..dcfc5da2f 100644 --- a/mcstas-comps/monitors/Sqq_w_monitor.comp +++ b/mcstas-comps/monitors/Sqq_w_monitor.comp @@ -86,134 +86,130 @@ DECLARE INITIALIZE %{ /* Make checks for limits on qa, qb, w grid */ - M_N=create_darr3d(nE,nqa,nqb); - M_p=create_darr3d(nE,nqa,nqb); - M_p2=create_darr3d(nE,nqa,nqb); - M_Ns=create_darr2d(nqa,nqb); - M_ps=create_darr2d(nqa,nqb); - M_p2s=create_darr2d(nqa,nqb); - - dE=(Emax-Emin)/(1.0*nE-1); + M_N = create_darr3d (nE, nqa, nqb); + M_p = create_darr3d (nE, nqa, nqb); + M_p2 = create_darr3d (nE, nqa, nqb); + M_Ns = create_darr2d (nqa, nqb); + M_ps = create_darr2d (nqa, nqb); + M_p2s = create_darr2d (nqa, nqb); + + dE = (Emax - Emin) / (1.0 * nE - 1); // 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,k; - //double rx,ry,rz,rvx,rvy,rvz,rt,rsx,rsy,rsz,rp; - double rvx,rvy,rvz; - double Ei,Ef,E,Ki,Kf; - double qx,qy,qz; - double qqa, qqb; - double kix,kiy,kiz; - double kfx,kfy,kfz; - double t0,t1; - - - double nx,ny,nz; - nx=vx; - ny=vy; - nz=vz; - NORM(nx,ny,nz); - double v_div=scalar_prod(nx, ny, nz, 0, 1, 0); - /* Initial check to see if this neutron should be counted or not... */ - if (fabs(v_div)<=yheight/radius && cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight)) { - if(t0<0 && t1>0) { - /* This one hits our cylindrical band, arriving from inside */ - - int fail; - rvx = particle_getvar(_particle,vix,&fail); if(fail) rvx=0; - rvy = particle_getvar(_particle,viy,&fail); if(fail) rvy=0; - rvz = particle_getvar(_particle,viz,&fail); if(fail) rvz=0; - - /* Calculate energy transfer */ - Ei = VS2E*(rvx*rvx + rvy*rvy + rvz*rvz); - Ef = VS2E*( vx*vx + vy*vy + vz*vz); - E=Ei-Ef; - - /* calculate k vectors and momentum transfer*/ - kix=rvx; - kiy=rvy; - kiz=rvz; - kfx=vx; - kfy=vy; - kfz=vz; - NORM(kix, kiy, kiz); - NORM(kfx, kfy, kfz); - - /* K-vector lengths */ - Ki=V2K*sqrt((rvx*rvx)+(rvy*rvy)+(rvz*rvz)); - Kf=V2K*sqrt((vx*vx)+(vy*vy)+(vz*vz)); - kix=Ki*kix; kiy=Ki*kiy; kiz=Ki*kiz; - kfx=Kf*kfx; kfy=Kf*kfy; kfz=Kf*kfz; - - qx=kix-kfx; - qy=kiy-kfy; - qz=kiz-kfz; - /* Calculate projections along qa and qb */ - qqa = qx*qax + qz*qaz; - qqb = qx*qbx + qz*qbz; - - /* Check if we are within the selected q/e range */ - if (qqa <= qamax && qqa >=qamin && qqb <=qbmax && qqb>=qbmin && E<= Emax && E>=Emin) { - - i = floor((qqa - qamin)*nqa/(qamax - qamin)); - j = floor((qqb - qbmin)*nqb/(qbmax - qbmin)); - k = floor((E - Emin)*nE/(Emax - Emin)); - - double p2 = p*p; - #pragma acc atomic - M_N[k][i][j] = M_N[k][i][j] + 1; - #pragma acc atomic - M_p[k][i][j] = M_p[k][i][j] + p; - #pragma acc atomic - M_p2[k][i][j] = M_p2[k][i][j] + p2; - - #pragma acc atomic - M_Ns[i][j] = M_Ns[i][j] + 1 ; - #pragma acc atomic - M_ps[i][j] = M_ps[i][j] + p; - #pragma acc atomic - M_p2s[i][j] = M_p2s[i][j] + p2; - - SCATTER; - } - } - } - /* Stuff that don't hit cylinder properly is disregarded */ + int i, j, k; + // double rx,ry,rz,rvx,rvy,rvz,rt,rsx,rsy,rsz,rp; + double rvx, rvy, rvz; + double Ei, Ef, E, Ki, Kf; + double qx, qy, qz; + double qqa, qqb; + double kix, kiy, kiz; + double kfx, kfy, kfz; + double t0, t1; + + double nx, ny, nz; + nx = vx; + ny = vy; + nz = vz; + NORM (nx, ny, nz); + double v_div = scalar_prod (nx, ny, nz, 0, 1, 0); + /* Initial check to see if this neutron should be counted or not... */ + if (fabs (v_div) <= yheight / radius && cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius, yheight)) { + if (t0 < 0 && t1 > 0) { + /* This one hits our cylindrical band, arriving from inside */ + + int fail; + rvx = particle_getvar (_particle, vix, &fail); + if (fail) + rvx = 0; + rvy = particle_getvar (_particle, viy, &fail); + if (fail) + rvy = 0; + rvz = particle_getvar (_particle, viz, &fail); + if (fail) + rvz = 0; + + /* Calculate energy transfer */ + Ei = VS2E * (rvx * rvx + rvy * rvy + rvz * rvz); + Ef = VS2E * (vx * vx + vy * vy + vz * vz); + E = Ei - Ef; + + /* calculate k vectors and momentum transfer*/ + kix = rvx; + kiy = rvy; + kiz = rvz; + kfx = vx; + kfy = vy; + kfz = vz; + NORM (kix, kiy, kiz); + NORM (kfx, kfy, kfz); + + /* K-vector lengths */ + Ki = V2K * sqrt ((rvx * rvx) + (rvy * rvy) + (rvz * rvz)); + Kf = V2K * sqrt ((vx * vx) + (vy * vy) + (vz * vz)); + kix = Ki * kix; + kiy = Ki * kiy; + kiz = Ki * kiz; + kfx = Kf * kfx; + kfy = Kf * kfy; + kfz = Kf * kfz; + + qx = kix - kfx; + qy = kiy - kfy; + qz = kiz - kfz; + /* Calculate projections along qa and qb */ + qqa = qx * qax + qz * qaz; + qqb = qx * qbx + qz * qbz; + + /* Check if we are within the selected q/e range */ + if (qqa <= qamax && qqa >= qamin && qqb <= qbmax && qqb >= qbmin && E <= Emax && E >= Emin) { + + i = floor ((qqa - qamin) * nqa / (qamax - qamin)); + j = floor ((qqb - qbmin) * nqb / (qbmax - qbmin)); + k = floor ((E - Emin) * nE / (Emax - Emin)); + + double p2 = p * p; + #pragma acc atomic + M_N[k][i][j] = M_N[k][i][j] + 1; + #pragma acc atomic + M_p[k][i][j] = M_p[k][i][j] + p; + #pragma acc atomic + M_p2[k][i][j] = M_p2[k][i][j] + p2; + + #pragma acc atomic + M_Ns[i][j] = M_Ns[i][j] + 1; + #pragma acc atomic + M_ps[i][j] = M_ps[i][j] + p; + #pragma acc atomic + M_p2s[i][j] = M_p2s[i][j] + p2; + + SCATTER; + } + } + } + /* Stuff that don't hit cylinder properly is disregarded */ %} SAVE %{ if (!nowritefile) { - - int kk,ll; - char ff[256]; - char tt[256]; - if (!nosum) { - sprintf(ff, "%s_Sum",filename); - DETECTOR_OUT_2D( - "qa vs qb monitor E Sum", - "qa [AA^-1]", - "qb [AA^-1]", - qamin, qamax, qbmin, qbmax, - nqa, nqb, - &M_Ns[0][0],&M_ps[0][0],&M_p2s[0][0], - ff); - } - for (kk=0; kk0) { - PROP_DT(t1); - } else { - detect=0; - } - } else { - detect=0; - } - } - } - - /* Check if we are within the selected q/e range */ - if (detect && qq <= qmax && qq>=qmin && E<= Emax && E>=Emin) { - i = floor((qq - qmin)*nq/(qmax - qmin)); - j = floor((E - Emin)*nE/(Emax - Emin)); - - double p2 = p*p; - - #pragma acc atomic - M_N[i][j] = M_N[i][j] + 1 ; - #pragma acc atomic - M_p[i][j] = M_p[i][j] + p; - #pragma acc atomic - M_p2[i][j] = M_p2[i][j] + p2; - - SCATTER; + int i, j; + double rvx, rvy, rvz; + double Ei, Ef, E, Ki, Kf; + double qx, qy, qz; + double q, qq; + double kix, kiy, kiz; + double kfx, kfy, kfz; + double t0, t1; + double phi; + int detect = 1; + + int fail; + rvx = particle_getvar (_particle, vix, &fail); + if (fail) + rvx = 0; + rvy = particle_getvar (_particle, viy, &fail); + if (fail) + rvy = 0; + rvz = particle_getvar (_particle, viz, &fail); + if (fail) + rvz = 0; + + /* If initial state is with v=(0,0,0), detect nothing */ + if (!(rvx == 0 && rvy == 0 && rvz == 0)) { + Ei = VS2E * (rvx * rvx + rvy * rvy + rvz * rvz); + Ef = VS2E * (vx * vx + vy * vy + vz * vz); + E = Ef - Ei; + + /* calculate k vectors and momentum transfer*/ + kix = rvx; + kiy = rvy; + kiz = rvz; + kfx = vx; + kfy = vy; + kfz = vz; + NORM (kix, kiy, kiz); + NORM (kfx, kfy, kfz); + + /* K-vector lengths */ + Ki = V2K * sqrt ((rvx * rvx) + (rvy * rvy) + (rvz * rvz)); + Kf = V2K * sqrt ((vx * vx) + (vy * vy) + (vz * vz)); + kix = Ki * kix; + kiy = Ki * kiy; + kiz = Ki * kiz; + kfx = Kf * kfx; + kfy = Kf * kfy; + kfz = Kf * kfz; + + qx = kfx - kix; + qy = kfy - kiy; + qz = kfz - kiz; + + qq = sqrt (qx * qx + qy * qy + qz * qz); + + /* Check if we should detect or not */ + if (detect) { + if (radius && yheight) { + if (cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius, yheight) == 1) { + if (t0 < 0 && t1 > 0) { + PROP_DT (t1); + } else { + detect = 0; + } + } else { + detect = 0; + } } - RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); - } + } + + /* Check if we are within the selected q/e range */ + if (detect && qq <= qmax && qq >= qmin && E <= Emax && E >= Emin) { + i = floor ((qq - qmin) * nq / (qmax - qmin)); + j = floor ((E - Emin) * nE / (Emax - Emin)); + + double p2 = p * p; + + #pragma acc atomic + M_N[i][j] = M_N[i][j] + 1; + #pragma acc atomic + M_p[i][j] = M_p[i][j] + p; + #pragma acc atomic + M_p2[i][j] = M_p2[i][j] + p2; + + SCATTER; + } + RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); + } %} SAVE %{ if (!nowritefile) { - DETECTOR_OUT_2D( - "q vs E monitor", - "q [AA^-1]", - "E [meV]", - qmin, qmax, Emin, Emax, - nq, nE, - &M_N[0][0],&M_p[0][0],&M_p2[0][0], - filename); + DETECTOR_OUT_2D ("q vs E monitor", "q [AA^-1]", "E [meV]", qmin, qmax, Emin, Emax, nq, nE, &M_N[0][0], &M_p[0][0], &M_p2[0][0], filename); } %} diff --git a/mcstas-comps/monitors/TOF2E_monitor.comp b/mcstas-comps/monitors/TOF2E_monitor.comp index 08f21a4fc..effe6bae0 100644 --- a/mcstas-comps/monitors/TOF2E_monitor.comp +++ b/mcstas-comps/monitors/TOF2E_monitor.comp @@ -68,24 +68,31 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("E_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("E_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - E_N = create_darr1d(nE); - E_p = create_darr1d(nE); - E_p2 = create_darr1d(nE); + E_N = create_darr1d (nE); + E_p = create_darr1d (nE); + E_p2 = create_darr1d (nE); S_p = S_pE = S_pE2 = 0; // 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 @@ -94,61 +101,50 @@ TRACE double E; PROP_Z0; - if (x>xmin && xymin && y xmin && x < xmax && y > ymin && y < ymax) { + E = VS2E * (L_flight / (t - T_zero)) * (L_flight / (t - T_zero)); S_p += p; - S_pE += p*E; - S_pE2 += p*E*E; + S_pE += p * E; + S_pE2 += p * E * E; - i = floor((E-Emin)*nE/(Emax-Emin)); - if(i >= 0 && i < nE) - { - double p2 = p*p; + i = floor ((E - Emin) * nE / (Emax - Emin)); + if (i >= 0 && i < nE) { + double p2 = p * p; #pragma acc atomic - E_N[i] = E_N[i]+1; + E_N[i] = E_N[i] + 1; #pragma acc atomic - E_p[i] = E_p[i]+p; + E_p[i] = E_p[i] + p; #pragma acc atomic - E_p2[i] = E_p2[i]+p2; + E_p2[i] = E_p2[i] + p2; SCATTER; } } 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( - "TOF-to-Energy monitor", - "Energy [meV]", - "Intensity", - "E", Emin, Emax, nE, - &E_N[0],&E_p[0],&E_p2[0], - filename); - if (S_p) printf(" : %g meV , E-width : %g meV \n", - S_pE/S_p,sqrt(S_pE2/S_p - S_pE*S_pE/(S_p*S_p)) ); -} + if (!nowritefile) { + DETECTOR_OUT_1D ("TOF-to-Energy monitor", "Energy [meV]", "Intensity", "E", Emin, Emax, nE, &E_N[0], &E_p[0], &E_p2[0], filename); + if (S_p) + printf (" : %g meV , E-width : %g meV \n", S_pE / S_p, sqrt (S_pE2 / S_p - S_pE * S_pE / (S_p * S_p))); + } %} FINALLY %{ - destroy_darr1d(E_N); - destroy_darr1d(E_p); - destroy_darr1d(E_p2); + destroy_darr1d (E_N); + destroy_darr1d (E_p); + destroy_darr1d (E_p2); %} 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 diff --git a/mcstas-comps/monitors/TOF2Q_cylPSD_monitor.comp b/mcstas-comps/monitors/TOF2Q_cylPSD_monitor.comp index d5ef157b1..4d5a8ac77 100644 --- a/mcstas-comps/monitors/TOF2Q_cylPSD_monitor.comp +++ b/mcstas-comps/monitors/TOF2Q_cylPSD_monitor.comp @@ -62,16 +62,20 @@ DECLARE INITIALIZE %{ - int i,j; + int i, j; - TOF_N = create_darr2d(nQ, ny); - TOF_p = create_darr2d(nQ, ny); - TOF_p2 = create_darr2d(nQ, ny); + TOF_N = create_darr2d (nQ, ny); + TOF_p = create_darr2d (nQ, ny); + TOF_p2 = create_darr2d (nQ, ny); - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } // 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 @@ -80,63 +84,57 @@ TRACE double cyl_t0, cyl_t1, dt, phi; double E, Q; - if(!cylinder_intersect(&cyl_t0, &cyl_t1, x,y,z,vx,vy,vz, radius, yheight)) - /* No hit */ + if (!cylinder_intersect (&cyl_t0, &cyl_t1, x, y, z, vx, vy, vz, radius, yheight)) + /* No hit */ ABSORB; - if(cyl_t0>0) /* Neutron hits cylinder from the outside */ + if (cyl_t0 > 0) /* Neutron hits cylinder from the outside */ ABSORB; dt = cyl_t1; - PROP_DT(dt); - if(y>=yheight/2 || y<= -yheight/2) - ABSORB; /* Neutron hits cylinder ends; no detectors here */ + PROP_DT (dt); + if (y >= yheight / 2 || y <= -yheight / 2) + ABSORB; /* Neutron hits cylinder ends; no detectors here */ - E = VS2E*(L_flight/(t-T_zero))*(L_flight/(t-T_zero)); + E = VS2E * (L_flight / (t - T_zero)) * (L_flight / (t - T_zero)); - j = floor((y - ymin)*ny/(ymax - ymin)); + j = floor ((y - ymin) * ny / (ymax - ymin)); - Q=2*sqrt(E/2.072)*sin(theta); + Q = 2 * sqrt (E / 2.072) * sin (theta); - i = floor((Q-Qmin)*nQ/(Qmax-Qmin)); - - if (i>=0 && i=0 && j= 0 && i < nQ && j >= 0 && j < ny) { + double p2 = p * p; #pragma acc atomic - TOF_N[i][j] = TOF_N[i][j]+1; + TOF_N[i][j] = TOF_N[i][j] + 1; #pragma acc atomic - TOF_p[i][j] = TOF_p[i][j]+p; + TOF_p[i][j] = TOF_p[i][j] + p; #pragma acc atomic - TOF_p2[i][j] = TOF_p2[i][j]+p2; + TOF_p2[i][j] = TOF_p2[i][j] + 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_2D( - "Cylindrical TOF2E PSD monitor", - "TOF2Q [Ang^-1]", - "Y position [cm]", - Qmin, Qmax, ymin*100.0, ymax*100.0, - nQ, ny, - &TOF_N[0][0],&TOF_p[0][0],&TOF_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("Cylindrical TOF2E PSD monitor", "TOF2Q [Ang^-1]", "Y position [cm]", Qmin, Qmax, ymin * 100.0, ymax * 100.0, nQ, ny, &TOF_N[0][0], + &TOF_p[0][0], &TOF_p2[0][0], filename); + } %} FINALLY %{ - destroy_darr2d(TOF_N); - destroy_darr2d(TOF_p); - destroy_darr2d(TOF_p2); + destroy_darr2d (TOF_N); + destroy_darr2d (TOF_p); + destroy_darr2d (TOF_p2); %} MCDISPLAY %{ - magnify("y"); - circle("xz", 0,0,0,radius); + magnify ("y"); + circle ("xz", 0, 0, 0, radius); %} END diff --git a/mcstas-comps/monitors/TOFLambda_monitor.comp b/mcstas-comps/monitors/TOFLambda_monitor.comp index 5b02ba826..da3191ef4 100644 --- a/mcstas-comps/monitors/TOFLambda_monitor.comp +++ b/mcstas-comps/monitors/TOFLambda_monitor.comp @@ -68,82 +68,80 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("TOFlambda_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("TOFlambda_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - TOFL_N = create_darr2d(nt, nL); - TOFL_p = create_darr2d(nt, nL); - TOFL_p2 = create_darr2d(nt, nL); - tt_0 = tmin*1e-6; - tt_1 = tmax*1e-6; + TOFL_N = create_darr2d (nt, nL); + TOFL_p = create_darr2d (nt, nL); + TOFL_p2 = create_darr2d (nt, nL); + tt_0 = tmin * 1e-6; + tt_1 = tmax * 1e-6; // 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; PROP_Z0; - lambda = (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz); - if (x>xmin && xymin && y Lmin && lambda < Lmax){ - if (t < tt_1 && t > tt_0) - { - i = floor((lambda - Lmin)*nL/(Lmax - Lmin)); - j = floor((t-tt_0)*nt/(tt_1-tt_0)); + lambda = (2 * PI / V2K) / sqrt (vx * vx + vy * vy + vz * vz); + if (x > xmin && x < xmax && y > ymin && y < ymax && lambda > Lmin && lambda < Lmax) { + if (t < tt_1 && t > tt_0) { + i = floor ((lambda - Lmin) * nL / (Lmax - Lmin)); + j = floor ((t - tt_0) * nt / (tt_1 - tt_0)); /* printf("tt_0, tt_1, nt %g %g %i t j %g %i \n",tt_0,tt_1,nt,t,j); - */ - double p2 = p*p; + */ + double p2 = p * p; #pragma acc atomic - TOFL_N[j][i] = TOFL_N[j][i]+1; + TOFL_N[j][i] = TOFL_N[j][i] + 1; #pragma acc atomic - TOFL_p[j][i] = TOFL_p[j][i]+p; + TOFL_p[j][i] = TOFL_p[j][i] + p; #pragma acc atomic - TOFL_p2[j][i] = TOFL_p2[j][i]+p2; + TOFL_p2[j][i] = TOFL_p2[j][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_2D( - "TOF-wavelength monitor", - "Time-of-flight [\\gms]", "Wavelength [AA]", - tmin, tmax, Lmin, Lmax, - nt, nL, - &TOFL_N[0][0],&TOFL_p[0][0],&TOFL_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("TOF-wavelength monitor", "Time-of-flight [\\gms]", "Wavelength [AA]", tmin, tmax, Lmin, Lmax, nt, nL, &TOFL_N[0][0], &TOFL_p[0][0], + &TOFL_p2[0][0], filename); + } %} FINALLY %{ - destroy_darr2d(TOFL_N); - destroy_darr2d(TOFL_p); - destroy_darr2d(TOFL_p2); + destroy_darr2d (TOFL_N); + destroy_darr2d (TOFL_p); + destroy_darr2d (TOFL_p2); %} 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 diff --git a/mcstas-comps/monitors/TOF_PSD_monitor_rad.comp b/mcstas-comps/monitors/TOF_PSD_monitor_rad.comp index ac2040485..d32be50e6 100644 --- a/mcstas-comps/monitors/TOF_PSD_monitor_rad.comp +++ b/mcstas-comps/monitors/TOF_PSD_monitor_rad.comp @@ -60,91 +60,79 @@ DECLARE INITIALIZE %{ - TOFPSDr_N = create_darr2d(nt, nr); - TOFPSDr_p = create_darr2d(nt, nr); - TOFPSDr_p2 = create_darr2d(nt, nr); - TOFPSDr_av_p = create_darr2d(nt, nr); - TOFPSDr_av_p2 = create_darr2d(nt, nr); + TOFPSDr_N = create_darr2d (nt, nr); + TOFPSDr_p = create_darr2d (nt, nr); + TOFPSDr_p2 = create_darr2d (nt, nr); + TOFPSDr_av_p = create_darr2d (nt, nr); + TOFPSDr_av_p2 = create_darr2d (nt, nr); - tt_0=tmin*1e-6; - tt_1=tmax*1e-6; + tt_0 = tmin * 1e-6; + tt_1 = tmax * 1e-6; // 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); // Use instance name for monitor output if no input was given - if (!strcmp(filename_av,"\0")) sprintf(filename_av,"%s%s",NAME_CURRENT_COMP,"_avg"); + if (!strcmp (filename_av, "\0")) + sprintf (filename_av, "%s%s", NAME_CURRENT_COMP, "_avg"); %} TRACE %{ - int i,j; - double dt,radpos; + int i, j; + double dt, radpos; PROP_Z0; - radpos = sqrt(x*x+y*y); - i = floor((t-tt_0)*nt/(tt_1-tt_0)); /* Bin number */ - j = floor(nr*radpos/rmax); + radpos = sqrt (x * x + y * y); + i = floor ((t - tt_0) * nt / (tt_1 - tt_0)); /* Bin number */ + j = floor (nr * radpos / rmax); - if (radpos < rmax && i>=0 && i < nt) - { - double p2 = p*p; + if (radpos < rmax && i >= 0 && i < nt) { + double p2 = p * p; #pragma acc atomic - TOFPSDr_N[i][j] = TOFPSDr_N[i][j]+1; + TOFPSDr_N[i][j] = TOFPSDr_N[i][j] + 1; #pragma acc atomic - TOFPSDr_p[i][j] = TOFPSDr_p[i][j]+p; + TOFPSDr_p[i][j] = TOFPSDr_p[i][j] + p; #pragma acc atomic - TOFPSDr_p2[i][j] = TOFPSDr_p2[i][j]+p2; + TOFPSDr_p2[i][j] = TOFPSDr_p2[i][j] + p2; SCATTER; } 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) { - int i, j; - for(i=0; i0) /* Neutron hits cylinder from the outside */ + if (cyl_t0 > 0) /* Neutron hits cylinder from the outside */ ABSORB; - dt=cyl_t1; - PROP_DT(dt); - if(y>=yheight/2 || y<= -yheight/2) - ABSORB; /* Neutron hits cylinder ends; no detectors here */ + dt = cyl_t1; + PROP_DT (dt); + if (y >= yheight / 2 || y <= -yheight / 2) + ABSORB; /* Neutron hits cylinder ends; no detectors here */ - i = floor((t-tt_0)*nt/(tt_1-tt_0)); /* Bin number */ + i = floor ((t - tt_0) * nt / (tt_1 - tt_0)); /* Bin number */ - if (i < 0 || i >= nt) /* Do not detect */ + if (i < 0 || i >= nt) /* Do not detect */ { - } - else - { - phi = atan2(x,z); - j = floor((double)nphi/2.0 + RAD2DEG*phi/(double)binphi); - double p2 = p*p; - #pragma acc atomic - TOF_N[i][j] = TOF_N[i][j]+1; - #pragma acc atomic - TOF_p[i][j] = TOF_p[i][j]+p; - #pragma acc atomic - TOF_p2[i][j] = TOF_p2[i][j]+p2; + } else { + phi = atan2 (x, z); + j = floor ((double)nphi / 2.0 + RAD2DEG * phi / (double)binphi); + double p2 = p * p; + #pragma acc atomic + TOF_N[i][j] = TOF_N[i][j] + 1; + #pragma acc atomic + TOF_p[i][j] = TOF_p[i][j] + p; + #pragma acc atomic + TOF_p2[i][j] = TOF_p2[i][j] + 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_2D( - "Cylindrical Time-of-flight PSD monitor", - "Time-of-flight [\\gms]", - "Angle [deg]", - tmin, tmax, -180, 180, - nt, nphi, - &TOF_N[0][0],&TOF_p[0][0],&TOF_p2[0][0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_2D ("Cylindrical Time-of-flight PSD monitor", "Time-of-flight [\\gms]", "Angle [deg]", tmin, tmax, -180, 180, nt, nphi, &TOF_N[0][0], + &TOF_p[0][0], &TOF_p2[0][0], filename); + } %} FINALLY %{ - destroy_darr2d(TOF_N); - destroy_darr2d(TOF_p); - destroy_darr2d(TOF_p2); + destroy_darr2d (TOF_N); + destroy_darr2d (TOF_p); + destroy_darr2d (TOF_p2); %} MCDISPLAY %{ - magnify("y"); - circle("xz", 0,0,0,radius); + magnify ("y"); + circle ("xz", 0, 0, 0, radius); %} END diff --git a/mcstas-comps/monitors/TOF_monitor.comp b/mcstas-comps/monitors/TOF_monitor.comp index 7d1d0c825..0630218a1 100644 --- a/mcstas-comps/monitors/TOF_monitor.comp +++ b/mcstas-comps/monitors/TOF_monitor.comp @@ -57,42 +57,46 @@ DECLARE DArray1d TOF_N; DArray1d TOF_p; DArray1d TOF_p2; - double t_min; + double t_min; double t_max; double delta_t; %} INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("TOF_monitor: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("TOF_monitor: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - TOF_N = create_darr1d(nt); - TOF_p = create_darr1d(nt); - TOF_p2 = create_darr1d(nt); - - if (tmax!=0) - { - t_max=tmax; - t_min=tmin; - delta_t=(t_max-t_min)/nt; - } - else - { - delta_t=dt; - t_min=0; - t_max=nt*dt+tmin; + TOF_N = create_darr1d (nt); + TOF_p = create_darr1d (nt); + TOF_p2 = create_darr1d (nt); + + if (tmax != 0) { + t_max = tmax; + t_min = tmin; + delta_t = (t_max - t_min) / nt; + } else { + delta_t = dt; + t_min = 0; + t_max = nt * dt + tmin; } // 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 @@ -100,53 +104,43 @@ TRACE int i; PROP_Z0; - if (x>xmin && xymin && y= 0 && i < nt) { - double p2 = p*p; - #pragma acc atomic - TOF_N[i] = TOF_N[i]+1; + if (x > xmin && x < xmax && y > ymin && y < ymax) { + i = floor ((1E6 * t - t_min) / delta_t); /* Bin number */ + if (i >= 0 && i < nt) { + double p2 = p * p; + #pragma acc atomic + TOF_N[i] = TOF_N[i] + 1; + #pragma acc atomic + TOF_p[i] = TOF_p[i] + p; #pragma acc atomic - TOF_p[i] = TOF_p[i]+p; - #pragma acc atomic - TOF_p2[i] = TOF_p2[i]+p2; + TOF_p2[i] = TOF_p2[i] + p2; SCATTER; } } 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( - "Time-of-flight monitor", - "Time-of-flight [\\gms]", - "Intensity", - "t", t_min, t_max, nt, - &TOF_N[0],&TOF_p[0],&TOF_p2[0], - filename); + if (!nowritefile) { + DETECTOR_OUT_1D ("Time-of-flight monitor", "Time-of-flight [\\gms]", "Intensity", "t", t_min, t_max, nt, &TOF_N[0], &TOF_p[0], &TOF_p2[0], filename); } %} FINALLY %{ - destroy_darr1d(TOF_N); - destroy_darr1d(TOF_p); - destroy_darr1d(TOF_p2); + destroy_darr1d (TOF_N); + destroy_darr1d (TOF_p); + destroy_darr1d (TOF_p2); %} 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 diff --git a/mcstas-comps/monitors/TOFlog_monitor.comp b/mcstas-comps/monitors/TOFlog_monitor.comp index 0facef32e..3329f134a 100644 --- a/mcstas-comps/monitors/TOFlog_monitor.comp +++ b/mcstas-comps/monitors/TOFlog_monitor.comp @@ -65,35 +65,44 @@ DECLARE INITIALIZE %{ - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; + } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("TOFlog_mon: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("TOFlog_mon: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - + if (tmin <= 0 || tmax <= 0) { - printf("TOFlog_mon: %s: Only supports positive tmin and tmax !\n" - "ERROR (tmin, tmax). Exiting", NAME_CURRENT_COMP); - exit(0); + printf ("TOFlog_mon: %s: Only supports positive tmin and tmax !\n" + "ERROR (tmin, tmax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - + if (tmin >= tmax) { - printf("TOFlog_mon: %s: tmin should be smaller than tmax !\n" - "ERROR (tmin, tmax). Exiting", NAME_CURRENT_COMP); - exit(0); + printf ("TOFlog_mon: %s: tmin should be smaller than tmax !\n" + "ERROR (tmin, tmax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - nchan=(int)ceil(ndec*log(tmax/tmin)/log(10.0)); - TOF_N = create_darr1d(nchan); - TOF_p = create_darr1d(nchan); - TOF_p2 = create_darr1d(nchan); + nchan = (int)ceil (ndec * log (tmax / tmin) / log (10.0)); + TOF_N = create_darr1d (nchan); + TOF_p = create_darr1d (nchan); + TOF_p2 = create_darr1d (nchan); // 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 @@ -101,53 +110,44 @@ TRACE int i; PROP_Z0; - if (x>xmin && xymin && y= 0 && i < nchan) { - double p2 = p*p; - #pragma acc atomic - TOF_N[i] = TOF_N[i]+1; + if (x > xmin && x < xmax && y > ymin && y < ymax) { + i = (int)floor (ndec * log (1E6 * t / tmin) / log (10.0)); /* Bin number */ + if (i >= 0 && i < nchan) { + double p2 = p * p; #pragma acc atomic - TOF_p[i] = TOF_p[i]+p; - #pragma acc atomic - TOF_p2[i] = TOF_p2[i]+p2; + TOF_N[i] = TOF_N[i] + 1; + #pragma acc atomic + TOF_p[i] = TOF_p[i] + p; + #pragma acc atomic + TOF_p2[i] = TOF_p2[i] + p2; SCATTER; } } 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( - "Time-of-flight monitor", - "Log(Time-of-flight [\\gms])", - "Intensity", - "t", log(tmin)/log(10.0), log(tmax)/log(10.0), nchan, - &TOF_N[0],&TOF_p[0],&TOF_p2[0], - filename); -} + if (!nowritefile) { + DETECTOR_OUT_1D ("Time-of-flight monitor", "Log(Time-of-flight [\\gms])", "Intensity", "t", log (tmin) / log (10.0), log (tmax) / log (10.0), nchan, + &TOF_N[0], &TOF_p[0], &TOF_p2[0], filename); + } %} FINALLY %{ - destroy_darr1d(TOF_N); - destroy_darr1d(TOF_p); - destroy_darr1d(TOF_p2); + destroy_darr1d (TOF_N); + destroy_darr1d (TOF_p); + destroy_darr1d (TOF_p2); %} 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