/******************************************************************************* * McStas, neutron ray-tracing package * Copyright (C) 1997-2008 Risoe National Laboratory, Roskilde, Denmark * * Component: NPI_tof_dhkl_detector * * %I * Written by: Jan Saroun, saroun@ujf.cas.cz * Date: July 1, 2017 * Version: $Revision$ * Origin: NPI * Release: McStas 2.5 * Last update: 30/11/2018 * * NPI detector for estimating dhkl from tof * * %D * A cylindrical detector which converts time-of-flight data (x,y,z,time) to a 1D diffractogram in dhkl. * The detector model includes finite detection depth, spatial and time resolution. Optionally, the component * exports position and time coordinates of detection events in an ASCII file. * * The component can also handle setups employing a modulation chopper for peak multiplexing at a long pulse source, * as proposed for the diffractometer BEER@ESS. In this case, the component requires chopper parameters * (modulation period, width of the primary unmodulated pulse), and a file with estimated dhkl values. * The component then estimates valid regions on the angle/tof map, excluding areas assumed to be empty * or with overlaping lines. This map is exported together with the diffractogram. * * Tips for usage: * 1) Centre the detector at the sample axis, keep z-axis parallel to the incident beam. * 2) Set the radius equal to the distance from the sample axis to the front face of the detection volume. * 3) Linst-Lc is used to determine "instrumental" wavelength as lambda = h/m_n*time/(Linst-Lc); * * %P * INPUT PARAMETERS: * * filename: [str] Name of file in which to store the detector data. * radius: [m] Radius of detector. * yheight: [m] Height of detector. * zdepth: [m] Thickness of the detection volume. * amin: [deg] minimum of scattering angle to be detected. * amax: [deg] maximum of scattering angle to be detected. * d_min: [AA] minimum of inter-planar spacing to be detected. * d_max: [AA] maximum of inter-planar spacing to be detected. * time0: [s] Time delay of the wavelength definition chopper with respect to the source pulse. * Lc: [m] Distance from the source to the pulse chopper (set 0 for a short pulse source). * Linst: [m] Distance from the source to the detector. * res_x: [m] Spatial resolution along x (Gaussian FWHM). * res_y: [m] Spatial resolution along y (Gaussian FWHM). * res_t: [s] Time resolution (Gaussian FWHM). Readout only, path to capture is accounted for by the tracing code. * mu: [1/cm/AA] Capture coefficient for the detection medium. * modulation: [1|0] Switch on/off the modulation regime (BEER multiplexing technique). * mod_dt: [s] Modulation period. * mod_twidth: [s] Width of the primary source pulse. * mod_shift: [float] Relative shift of the modulation pattern (delta_d/d, constant for all reflections). * mod_d0_table: [str] Name of a file with the list of dhkl estimates (one per line). * verbose: [1|0] Switch for extended reporting. * restore_neutron: [1|0] Switch for restoring of previous neutron state. * * DEFINITION PARAMETERS: * * nd: Number of bins for dhkl in the 1D diffractograms. * na: Number of bins for scattering angles in the 2D maps. * nt: Number of bins for time of flight in the 2D maps. * nev: Number of detection events to export in "events.dat" (only modulation mode, set 1 for no export) * * %E * Example 1, simple ToF regime * NPI_tof_dhkl_detector(nd = 1800, filename = "result.dat", yheight = 0.5, zdepth = 0.02, * radius = 1.5, amin = 75, amax = 105, d_min = 0.4, d_max = 2.4, time0 = chopper_delay, Linst = 160, Lc = 6.5, modulation = 0) * * Example 2, modulation regime * NPI_tof_dhkl_detector(nd=1800, filename="result.dat", yheight=1.0, zdepth=0.01, radius=2, * yheight=1.0, zdepth=0.01, radius=2, amin=det_th1, amax=det_th2, d_min=0.8, d_max=2.4, * time0=chopper_delay, Linst=160, Lc = 9.3, res_x=0.002, res_y=0.005, res_t=1e-6, mu=1.0, mod_shift=5e-4, * modulation=1, mod_dt=4.464e-4, mod_twidth=0.003, mod_d0_table="dhkl.dat") * ******************************************************************************/ DEFINE COMPONENT NPI_tof_dhkl_detector SETTING PARAMETERS (string filename=0, radius=2, yheight=0.3, zdepth=0.01, amin=75, amax=105, int nd=200, int na=800, int nt=800, int nev=1, d_min=1, d_max=3, time0, Linst, Lc, res_x=0, res_y=0, res_t=0, mu=1.0, mod_shift=0, modulation=0, mod_dt=0, mod_twidth=0, string mod_d0_table=0, verbose=0, restore_neutron=1) /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ /* used for reading data table from file */ %include "read_table-lib" /* Save up to nmax events stored by the detector in a text file, one per row. Each row contains detection coordinates: x, y, t, sin(theta), p. */ void saveEvents (double** EVENTS, char* file_name, long int nmax, long int iev) { long int i, j, n; double sumP; FILE* hfile; hfile = fopen (file_name, "w"); if (!hfile) { fprintf (stderr, "saveEvents: can't open file for output (%s)\n", file_name); exit (-1); } else { n = nmax; if (n > iev) { n = iev; } sumP = 0.0; for (i = 0; i < n; i++) { sumP += EVENTS[i][5]; } fprintf (hfile, "# Detection events\n"); fprintf (hfile, "# sum(p) = %g\n", sumP); fprintf (hfile, "# sum(n) = %ld\n", n); fprintf (hfile, "# Columns: x [mm], y[mm], z[mm], t[ms], sin(theta), p\n"); for (i = 0; i < n; i++) { for (j = 0; j < 6; j++) { fprintf (hfile, "%g ", EVENTS[i][j]); } fprintf (hfile, "\n"); } fclose (hfile); } } /* Used in the modulation mode: Calculates the index of the nearest line which fits within given time interval (trange). Excludes the line with given index (iex). Parameters ----------- iex: excluding this line index sinth: sin(thetaB) tof: neutron time of flight from the source (not from the chopper !) [s] trange: time range for searching overlap [s] Returns ------- The index of the nearest diffraction line from the dhkl list, or -1 if not found. */ int getNearestLine (int iex, double sinth, double tof, double trange, int n_dhkl, double* dhkl, double mod_shift, double Ltof) { int i; int j = 0; int io = -1; double dtmax, dt, tline; const double hm = 2 * PI * K2V; dtmax = trange; for (i = 0; i < n_dhkl; i++) { if (i != iex) { // get centre of the pulse chain tline = 2 * dhkl[i] * (1.0 + mod_shift) * sinth / hm * Ltof; dt = fabs (tof - tline); if (dt < dtmax) { dtmax = dt; j = i; } } } if (dtmax < 0.5 * trange) { io = j; } return io; } /* Estimates the overlap map in the (scattering angle, tof) space for plotting with DETECTOR_OUT_2D. Returns ------- TOF_N, TOF_p, TOF_p2 arrays for 2D plots with DETECTOR_OUT_2D. The values are: 0: empty region 1: valid region (diffraction line without overlap) 2: overlaping of 2 or more lines */ void calcOverlaps (double tmin, double tmax, double ami, double ama, double trange, int na, int nt, double t2lam, double mod_shift, double Ltof, double** TOF_N, double** TOF_p, double** TOF_p2, int n_dhkl, double* dhkl) { int i, j, iline, io; double tdet, a, sth; double a1, a2, drange; a1 = ami * DEG2RAD; a2 = ama * DEG2RAD; double da = (a2 - a1) / na; double dt = (tmax - tmin) / nt; for (i = 0; i < na; i++) { a = a1 + (i + 0.5) * da; sth = sin (a / 2); drange = trange * t2lam / 2 / sth; for (j = 0; j < nt; j++) { TOF_N[i][j] = 0; TOF_p[i][j] = 0; TOF_p2[i][j] = 0; tdet = tmin + (j + 0.5) * dt; iline = getNearestLine (-1, sth, tdet, trange, n_dhkl, dhkl, mod_shift, Ltof); if (iline >= 0) { io = getNearestLine (iline, sth, tdet, trange, n_dhkl, dhkl, mod_shift, Ltof); if (io < 0) { TOF_N[i][j] += 1; TOF_p[i][j] += 1; TOF_p2[i][j] += 1; } else { TOF_N[i][j] += 1; TOF_p[i][j] += 2; TOF_p2[i][j] += 4; } } } } } %} DECLARE %{ DArray1d D_N; DArray1d D_p; DArray1d D_p2; DArray2d TOF_N; DArray2d TOF_p; DArray2d TOF_p2; DArray2d EVENTS; double* dhkl; double time_min; double time_max; double grf_th1; double grf_th2; double t2lam; double L0; double Ltof; long iev; int n_dhkl; // set to 1 for some debugging info int dbg; int dbgn; %} INITIALIZE %{ /* h/m_n in [Ang*m/s] */ double hm = 2 * PI * K2V; dbgn = 0; iev = -1; D_N = create_darr1d (nd); D_p = create_darr1d (nd); D_p2 = create_darr1d (nd); TOF_N = create_darr2d (na, nt); TOF_p = create_darr2d (na, nt); TOF_p2 = create_darr2d (na, nt); EVENTS = create_darr2d (nev, 6); // read dhkl table n_dhkl = 0; if (modulation) { MPI_MASTER (printf ("%s: Modulation on, table = %s, line shift=%g\n", NAME_CURRENT_COMP, mod_d0_table, mod_shift);) if (!mod_d0_table || !strlen (mod_d0_table) || !strcmp (mod_d0_table, "NULL")) { MPI_MASTER (fprintf (stderr, "ERROR %s: Can't read table with d0 values: %s \n", NAME_CURRENT_COMP, mod_d0_table);) exit (-1); } else { t_Table sTable; MPI_MASTER (printf ("trying to read dhkl table [%s]\n", mod_d0_table)); Table_Read (&sTable, mod_d0_table, 1); double size = sTable.rows; if (size > 0) { dhkl = (double*)malloc (sizeof (double) * size); int i; for (i = 0; i < size; i++) { dhkl[i] = Table_Index (sTable, i, 0); } n_dhkl = size; if (verbose) printf ("Read %d rows from dhkl table [%s]\n", n_dhkl, mod_d0_table); } Table_Free (&sTable); } if (!n_dhkl) { MPI_MASTER (fprintf (stderr, "ERROR %s: Can't evaluate modulated data without dhkl list\n", NAME_CURRENT_COMP);); exit (-1); } if (!n_dhkl) { MPI_MASTER (fprintf (stderr, "ERROR %s: Can't evaluate modulated data without dhkl list\n", NAME_CURRENT_COMP);); exit (-1); } if (mod_dt <= 0) { MPI_MASTER (fprintf (stderr, "ERROR %s: Modulation period must be positive\n", NAME_CURRENT_COMP);); exit (-1); } if (mod_twidth <= 0) { MPI_MASTER (fprintf (stderr, "ERROR %s: Time range in modulation mode must be positive\n", NAME_CURRENT_COMP);); exit (-1); } } else { MPI_MASTER (printf ("%s: Modulation off\n", NAME_CURRENT_COMP);); } // Distance L0 determines the wavelength L0 = Linst - Lc; Ltof = Linst; t2lam = hm / L0; // range of theta-tof plot time_min = 2 * d_min * sin (amin * DEG2RAD / 2) / t2lam; time_max = 2 * d_max * sin (amax * DEG2RAD / 2) / t2lam; grf_th1 = amin; grf_th2 = amax; if (modulation) { calcOverlaps (time_min, time_max, grf_th1, grf_th2, mod_twidth, na, nt, mod_shift, Ltof, t2lam, TOF_N, TOF_p, TOF_p2, n_dhkl, dhkl); } MPI_MASTER (if (verbose) { printf ("%s: lambda/(t-t0)=%g\n", NAME_CURRENT_COMP, t2lam); printf ("Lc=%g, L0=%g, LD=%g, time0=%g\n", Lc, L0, radius, time0 * 1000); printf ("For Fe110 at 90 deg: tof=%g [ms]\n", 2 * 2.1055 * sin (45 * DEG2RAD) / hm * Ltof); if (modulation) { printf ("mod_dt=%g [ms], mod_shift=%g\n", mod_dt, mod_shift); } }); %} TRACE %{ int i, j, io, iline, iref, valid_region; double t0, t1, lam, theta2, d, cos2; double sinth, dt, p0, tau, v0, x0, y0, z0, dx, dy, sig, tof, dd; double tref, tc; double r8ln2 = 2.355; const double hm = 2 * PI * K2V; // = h/m_n double L1; double th2_min = amin * PI / 180; double th2_max = amax * PI / 180; #ifndef OPENACC if (dbgn > 20) { #pragma acc atomic write dbg = 0; } #endif /* cross-section with the detector front face. Allow height = yheight+3*res_y to account for a random displacement due to vertical resolution. */ int cross = cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius, yheight + 3 * res_y); /* don't allow intersections with top/bottom cylinder walls, only neutrons from inside of the cylinder are allowed. */ if ((cross != 1) || (t0 > 0) || (t1 < 0)) { RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); } else { // go to cylinder surface = detector entry PROP_DT (t1); // save position before propagation through the detection volume x0 = x; y0 = y; z0 = z; // add random penetration within the detector depth v0 = sqrt (vx * vx + vy * vy + vz * vz); sig = 100 * mu * hm / v0; // mu is in 1/cm p0 = 1.0 - exp (-sig * zdepth); tau = -1.0 * log (1.0 - p0 * rand01 ()) / sig / v0; PROP_DT (tau); p *= p0; // add random shift by y-resolution dy = randnorm () * res_y / r8ln2; y0 += dy; // add random shift by x-resolution (normal to radius) dx = randnorm () * res_x / r8ln2; cos2 = z0 / radius; z0 += dx * sqrt (1. - cos2 * cos2); x0 += dx * cos2; // add random time resolution shift tof = t + randnorm () * res_t / r8ln2; // NOTE: detection coordinates are x0, y0, z0, tof // L1 = detected sample to detector distance: L1 = sqrt (x0 * x0 + z0 * z0 + y0 * y0); // Clip y on +- 0.5*yheight if (fabs (y) < (0.5 * yheight)) { dd = (d_max - d_min) / nd; // get cos(2*theta), assume primary beam axis // [0 0 1] cos2 = z0 / L1; theta2 = acos (cos2); dt = 0; if (theta2 > th2_min && theta2 < th2_max) { sinth = sin (theta2 / 2); valid_region = 0; int iex = -1; /* Modulation mode: estimate the number of modulation periods and correct dt. Exclude regions with overlaping diffraction lines */ if (modulation) { iline = getNearestLine (-1, sinth, tof, mod_twidth, n_dhkl, dhkl, mod_shift, Ltof); if (iline >= 0) { io = getNearestLine (iline, sinth, tof, mod_twidth, n_dhkl, dhkl, mod_shift, Ltof); if (io < 0) { valid_region = 1; /* Correction for modulation in steps: 1) tref = reference time for given dhkl, theta and reference chopper window 2) tc = neutron time minus tref 3) iref = new reference chopper window 4) dt = correction for tof */ tref = 2 * (1.0 + mod_shift) * dhkl[iline] * sinth / hm * L0; tc = tof - time0 - tref; iref = (int)floor (tc / mod_dt + 0.5); dt = iref * mod_dt + time0; #ifndef OPENACC if (iline == 0 && dbg) { #pragma acc atomic dbgn = dbgn + 1; printf ("iref=%d, dt=%g, t=%g, tcor=%g\n", iref, iref * mod_dt * 1000, tc * 1000, (tof - tref - dt) * 1000); } #endif } } // Normal mode: get indices in the ToF - 2theta grid and accumulate plot data } else { valid_region = 1; dt = time0; i = (int)floor ((theta2 - grf_th1 * DEG2RAD) / (grf_th2 * DEG2RAD - grf_th1 * DEG2RAD) * na + 0.5); j = (int)floor ((t - time_min) / (time_max - time_min) * nt + 0.5); if (j >= 0 && j < nt && i >= 0 && i < na) { 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; } } /* Valid event: calculate wavelength and dhkl "as recorded by the instrument", i.e. including smearing by instrument resolution. Record a list of detection event coordinates in EVENTS array. */ if (valid_region) { // subtract reference chopper time from tof and calculate wavelength lam = hm * (tof - dt) / (L0 - radius + L1); d = lam / 2 / sinth; #pragma acc atomic iev = iev + 1; if (iev < nev) { double tmp; tmp = x0 * 1e3; #pragma acc atomic write EVENTS[iev][0] = tmp; tmp = y0 * 1e3; #pragma acc atomic write EVENTS[iev][1] = tmp; tmp = z0 * 1e3; #pragma acc atomic write EVENTS[iev][2] = tmp; tmp = tof * 1e3; #pragma acc atomic write EVENTS[iev][3] = tmp; #pragma acc atomic write EVENTS[iev][4] = sinth; #pragma acc atomic write EVENTS[iev][5] = p; } i = (int)floor ((d - d_min) / dd + 0.5); if (i >= 0 && i < nd) { double p2 = p * p; #pragma acc atomic D_N[i] = D_N[i] + 1; #pragma acc atomic D_p[i] = D_p[i] + p; #pragma acc atomic D_p2[i] = D_p2[i] + p2; } } } } } if (restore_neutron) { RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); } %} SAVE %{ if (nev > 1) { saveEvents (EVENTS, "events.dat", nev, iev); } DETECTOR_OUT_1D ("ToF to dhkl detector", "dhkl [AA]", "Intensity", "d", d_min, d_max, nd, &D_N[0], &D_p[0], &D_p2[0], filename); if (modulation) { DETECTOR_OUT_2D ("Map of overlapping regions", "Scattering angle [deg]", "Time-of-flight [ms]", grf_th1, grf_th2, time_min * 1000, time_max * 1000, na, nt, &TOF_N[0][0], &TOF_p[0][0], &TOF_p2[0][0], "overlap_map.dat"); } %} FINALLY %{ free (dhkl); destroy_darr1d (D_N); destroy_darr1d (D_p); destroy_darr1d (D_p2); destroy_darr2d (TOF_N); destroy_darr2d (TOF_p); destroy_darr2d (TOF_p2); destroy_darr2d (EVENTS); %} MCDISPLAY %{ magnify ("xz"); circle ("xz", 0, 0, 0, radius); %} END