Add multi-pulse support to LoadMcStas in mantid

[huegle]

I have (like others) recently worked a bit with correlation choppers, pulsed sources and event data, and realized that mantid's loadmcstas assumes only a single pulse (because that is almost always what we simulate). I have submitted a feature request to mantid that adds optional reading of a new column titled "Source_Pulse_Index" to it (mantidproject/mantid#41458). How exactly one wants to generate that column is left to the user, but the way I did it is to add a time offset to the trace at the start of the instrument, and subtract that offset again before detection, to make sure that the TOF within the pulse is correct. Below is the short form of the instrument:

COMPONENT sourceMantid(...)
  EXTEND
    %{
    pulse_index = floor(rand01() * number_of_source_pulses);
    t += pulse_index / source_frequency_hz;
    %}
[...instrument simulation things happen...]
  EXTEND
  %{
    t -= pulse_index / source_frequency_hz;
  %}
COMPONENT nD_Mantid_1 = Monitor_nD(xwidth = 0.2, yheight = 0.9, 
    user1 = "pulse_index",
    username1 = "Source_Pulse_Index",
    options="mantid square, x bins=16, y bins=256, neutron pixel min=0, t, list all neutrons, user1")

This of course messes with the statistical weight (the single pulse intensity is now distributed over a lot of pulses), I guess that is the risk (and responsibility) I'll have to take when inserting EXTENDs. Either way, it should not be on mantid to fix that.
I think this capability for loadmcstas should come in handy, and might be useful for future, not one-off applications. It could for example be expanded for other uservars.
Sorry for dropping this in both places at once, it's kind of a chicken-and-egg problem. If there are major problems you see with the idea, now would be a great time to discuss!

[willend]

@huegle thanks for starting this discussion. :-)

For the ESS_butterfly approach we took the approach to

• Select the number of pulses to simulate
• Always simulate the n/s current irrespective of number of pulses (somewhat confusing to some users indeed that we simulate the current from 14 pulses always, irrespective of time-structure but chosen to ensure our legacy unit of n/s)
• Simulating with an extra MC choice in the source to add the required time-offset for neutrons starting in pulse k

AFAIK this should work within Mantid - but I may be wrong?

[willend]

(While we’re considering suggesting changes to the Mantid-side of things, we might as well think of making the loader somewhat more flexible, e.g. wrt. the number of columns ‘required’ / ‘used’ from the list of events…)

[huegle]

If i remember correctly(?), the experimental data in Mantid is stored as pixelID-TOF-pulseID, with a separate list for the wall time stamps of each pulseID. That allows it to switch between wall time and TOF. In the approach you described, would that work? I think mantid might still interpret it as just one extra-long pulse?

[willend]

You’re probably right - and probably also right that most people simply simulate just one pulse…

It could also be that what I seem to remember is instead some level of initial manipulation of the already loaded event workspace, based on the global time present with each event and a known source frequency? Modulus and remainder like?

[huegle]

hmmm. It might be that i misremembered. The nexus manual (https://manual.nexusformat.org/classes/base_classes/NXevent_data.html) says it is in fact not the pulseID, but the "time that each pulse started with respect to the offset" (I guess that means pulse T0). That would in this case just mean pulling one line of math from the load algorithm into the EXTEND part at the end.
I think you could do it in the ESS_butterfly case by just adding the T0s afterwards and doing a bit of math as you suggested. That mightnow introduce some very realistic frame overlap concerns? I think with the pulseID approach we would actually be able to label and disentangle the frame overlap stuff, which is a neat idea now that I think about it.

Overall, do you think this is a useful/useable addition to loadmcstas (either as pulseID or as pulse_t0)? So far, nobody has picked it up, and I could just shelve it for now.

[willend]

Since ESS is a Scipp- rather than Mantid-shop @mads-bertelsen has in fact made us some Python middleware called https://github.com/mccode-dev/McStasToX. (We have been thinking to eventually implement more than X==Scipp.)

@mads-bertelsen do you remember if anyone has analysed multi-pulse simulations via McStasToX and Scipp - and is there any sort of logic built-in for pulse-identification along the lines of what @huegle is after above?

[huegle]

Added complication: I will be on vacation for a week :-). I will take down my request on mantid for now, attach it below, and we can come to a clever conclusion once I'm back. I would hate for them to waste time on this if there are better ways.

---

The LoadMcStas algorithm currently reads 6 columns per neutron (Intensity, Position, Position, Neutron_ID, Pixel_ID, Time_Of_Flight), but assumes only a single pulse. I would like to add an option to read an optional 7th column: "Source_Pulse_Index" to allow conversion between TOF and wall time.

Background information:
When simulating a pulsed source, McStas usually simulates only a single pulse. This is in most cases enough, except when simulating something that will interact with multiple pulses (e.g. a correlation chopper that runs at a different frequency than the source).
The easiest way to deal with that on the McStas side is to just pretend that some of the neutrons starting at the regular single-pulse source belong to one of the following pulses. McStas tracks the neutrons by time, so we just need to add a time offset (pulse_index/source_frequency_hz) on neutron generation. Just ahead of the detector, we'll need to subtract that offset again, so that the detector records the "real" TOF from the source to the detector, and we need to record the pulse_index.
So the McStas instrument becomes:

COMPONENT sourceMantid(...)
  EXTEND
    %{
    pulse_index = floor(rand01() * number_of_source_pulses);
    t += pulse_index / source_frequency_hz;
    %}
[...instrument simulation things happen...]
  EXTEND
  %{
    t -= pulse_index / source_frequency_hz;
  %}
COMPONENT nD_Mantid_1 = Monitor_nD(xwidth = 0.2, yheight = 0.9, 
    user1 = "pulse_index",
    username1 = "Source_Pulse_Index",
    options="mantid square, x bins=16, y bins=256, neutron pixel min=0, t, list all neutrons, user1")

This will add a 7th column called "Source_Pulse_Index" to the resulting nexus file data, which will allow conversion between TOF and wall time similar to what happens in experimental event data. Caveat: nothing will be sorted, pulse_indexes and TOFs will appear rather randomly in the data because that is how it was generated in the simulation. Also, because of issues with uservars in McStas, the index might end up being a double instead of an integer.

Below are the relevant bits of the ai-generated response (ChatGPT 5.5 "thinking") to my rather crude sketch of an algorithm.

Proposed feature

Extend LoadMcStas to recognize a 7-column event dataset when the final column is a source pulse index, then convert that pulse index into the event pulse time.

1. Add a new optional property:

For an event with source pulse index n, the pulse time should be:
pulse_time = n / SourceFrequency
where SourceFrequency is a new loader property, defaulting to 60 Hz.

declareProperty("SourceFrequency", 60.0,
                "Source repetition frequency in Hz. Used to convert a "
                "McStas source pulse index event column into Mantid event "
                "pulse times.");

2. Allow either standard 6-column data or 7-column data with Source_Pulse_Index:

Replace the strict 6-column check:

if (nNeutrons && numberOfDataColumn != 6) {
  g_log.error() << "Event data in McStas nexus file expecting 6 columns\n";
  return scatteringWSNames;
}

with schema-aware logic:

const bool isStandardSixColumnMcStasEventData =
    numberOfDataColumn == 6;

const bool hasSourcePulseIndexColumn =
    numberOfDataColumn == 7 &&
    nameAttrValue.find("Source_Pulse_Index") != std::string::npos;

if (nNeutrons &&
    !isStandardSixColumnMcStasEventData &&
    !hasSourcePulseIndexColumn) {
  g_log.error() << "Event data in McStas nexus file expecting either "
                   "standard 6-column event data or 7-column event data "
                   "including Source_Pulse_Index. Found "
                << numberOfDataColumn << " columns.\n";
  return scatteringWSNames;
}

3. Populate event pulse time from Source_Pulse_Index:

Replace the current hard-coded pulse time:

int64_t pulse_time = 0;

with:

int64_t pulse_time = 0;

if (hasSourcePulseIndexColumn) {
  const double sourceFrequencyHz = getProperty("SourceFrequency");

  if (sourceFrequencyHz <= 0.0) {
    throw std::runtime_error("SourceFrequency must be positive.");
  }

  const auto sourcePulseIndex = static_cast<int64_t>(
      std::llround(data[6 + numberOfDataColumn * in])
  );

  pulse_time = static_cast<int64_t>(
      std::llround(
          static_cast<double>(sourcePulseIndex) / sourceFrequencyHz * 1.0e9
      )
  );
}

This requires:

#include <cmath>
for std::llround.

The conversion assumes that Mantid event pulse times are stored in nanoseconds relative to the workspace time origin.

4. Emit a warning when Source_Pulse_Index is used

Background: When using the above method the assign pulse ids, the intensity of the original single pulse is distributed amongst all the newly created pulses. It's not Mantid's job to fix that, and depending on how exactly the user choses to generate the pulse id in his simulation, it might not actually happen that way. Either way, it will be a good idea to add a warning.

When Source_Pulse_Index is detected, LoadMcStas should emit a warning that it assigns event pulse times but does not alter event weights.

Suggested warning text:

if (hasSourcePulseIndexColumn) {
  g_log.warning()
      << "McStas event data contain a Source_Pulse_Index column. "
      << "LoadMcStas will use this column to assign event pulse times, "
      << "but will not modify event weights. If the McStas instrument "
      << "assigned events from one simulated source pulse across multiple "
      << "source pulse indices, the user is responsible for applying any "
      << "required normalization or weight correction.\n";
}

This warning should be emitted once per load or once per event dataset, not once per event.

Non-goals

This feature does not add synthetic chopper logs such as motor speed or time-to-digital converter (TDC) timestamps.

This feature does not modify or renormalize event weights.

This feature does not define Source_Pulse_Index as a standard McStas convention. It only adds optional support in LoadMcStas for files that contain this explicitly named column.

Any chopper logs, pulse normalization, or absolute intensity correction should be handled by the simulation workflow, the McStas instrument definition, or a separate downstream algorithm.

Thank you for your consideration!

[willend]

@huegle thanks and let us discuss more once you are back.