True photon arrival time?

[nbiederbeck]

I could not find the arrival times in the SimulatedCameraContainer. Is there a reason for this?

Maybe relevant Issue: #1341

[maxnoe]

Yes, the reason is that these are not really available in a good manner. The simtel files contain optionally the arrival time for each photon, but these times are in the CORSIKA simulation time and there is no straight forward way to convert this to a time relative to the trigger time or something else that would tell you the position inside the waveform.

So these times are not really helpful...

Also, since these are variable length lists for each pixel, it is not really straight forward to include these in the event structure / store them into hdf5 files.

[maxnoe]

Maybe it's easier if @kbernloehr could summarize again what is needed to translate CORSIKA times into the trigger time / sample times used by simtel array again instead of me trying to translate the longish email he wrote me 2.5 years ago.

But I think the important takeaway was that

a) It's only feasible for simulations performed with simtel array later than 2019-08-06 b) even then its not that easy

[kosack]

Looking into this a bit deeper, it is indeed not trivial: the photon list can be binned into a waveform, as mentioned above the absolute time is not useful as these are before the detector simulation is run, so to turn them into a detected time would basically require re-implementing the detector sim or some part of it to turn the Corsika arrival time into a detector time.

I'm not sure it even makes any sense to do that though: you still have the waveform from the simulated shower, so you know the peak time after detector sims, however that is with noise. So perhaps we could consider adding a new block to the SimTel output with the noise-free time structure, or at least the peak time.

[maxnoe]

The request here is I think for the pe list (aka photonstream) not for one value per pixel. That is #1341

[maxnoe]

I would be very easy to add the list of pes as is to SimulatedCameraContainer but we'd have problems writing it into output files (since it is variable length)

[kosack]

Converting the photon list into a sequence of images (waveform per pixel) is not hard - I just tried it. You just need to know the sampling frequency and using np.histogram2d you can get a sequence of images. It just doesn't correspond at all to the trigger time.:

times = source.file_.current_photoelectrons[1]['time']
waveforms, xed, yed = np.histogram2d(pids, times bins=[np.arange(len(geom)+1),np.linspace(0,32)])

I guess this could still be useful for checking the time slope or similar, just not the absolute peak time.

[maxnoe]

@kosack the time array is a 1d array of number of photoelectrons for the whole camera. You at least need to use also the pixel_id array at the same point. And be careful with the key of that dict, it's tel_id - 1 not tel_id.

[kosack]

@kosack the time array is a 1d array of number of photoelectrons for the whole camera. You at least need to use also the pixel_id array at the same point. And be careful with the key of that dict, it's tel_id - 1 not tel_id.

Yes this was just a quick and dirty example - I also had to subtract some offset from the times to get it to start at 0 (since often they start at e.g. -180.0). It is actually in pixel_index here, so you need to map ID to index first. But the point is you can get an image sequence. The issue is what is the scale and starting point of the time axis?

[kosack]

from ctapipe.visualization import CameraDisplay
import matplotlib.pyplot as plt
import numpy as np

tid = 10
times = source.file_.current_photoelectrons[tid-1]["time"]
pids = source.file_.current_photoelectrons[tid-1]["pixel_id"] 
geom = source.subarray.tel[tid].camera.geometry

tmin = np.int(np.round(times.min()))
tmax = np.int(np.round(times.max()))
dt = tmax-tmin

waveforms, xed, yed = np.histogram2d(
    pids,
    times,
    bins=[np.arange(len(geom) + 1), np.linspace(tmin,tmax,dt)],
)
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
disp = CameraDisplay(geom, waveforms.argmax(axis=1), ax=ax[0], title="True peak time")
disp2 = CameraDisplay(geom, event.dl1.tel[tid].peak_time, ax=ax[1], title="DL1 Peak Time")

I still seem to be off by one event or telescope or something

[maxnoe]

Here's how you can get the true image:

from ctapipe.io import SimTelEventSource
from ctapipe.visualization import CameraDisplay
import numpy as np
import matplotlib.pyplot as plt

tel_id = 1

with SimTelEventSource('tests/resources/lst_with_photons.simtel.zst', focal_length_choice="EQUIVALENT") as f:

    cam = f.subarray.tel[tel_id].camera.geometry

    for e in f:
        pe = f.file_.current_photoelectrons[tel_id - 1]

        # these two arrays belong to each other.
        # each value in time is one photon arrival time and the
        # corresponding entry in pixel_ids is the pixel the pe hit.
        pixel_indices = pe['pixel_id']
        times = pe['time']

        image = np.zeros(cam.n_pixels)

        np.add.at(image, pixel_indices, times)

        disp = CameraDisplay(cam)
        disp.image = image
        plt.show()

Strangely subtracting 1 from the pixel_id gives me the wrong answer

The off-by-one is for telescope id, not pixel id.

[maxnoe]

The issue is what is the scale and starting point of the time axis?

It's the CORSIKA system, so it's nanoseconds since the entry into the atmosphere (if TSTART is true) or since the first interaction (if TSTART is false, which shouldn't be the case for any of our simulations).

[kosack]

But that's not quite what I wanted: I wanted the waveforms, not the sum, hence needing to histogram the time axis. Could use np.digitize for that axis and still use add.at perhaps though. I updated the above example

It's the CORSIKA system, so it's nanoseconds since the entry into the atmosphere (if TSTART is true)

So need to know the difference between this and the trigger time

[kosack]

I updated the example above, and it seems to give a reasonable result: However, the absolute value of the times are not comparable

The trick I employed was just to use 1ns bins (assuming a 1ns sample width, but could be made general by looking at the actual sample frequency), and then use np.argmax() to get the sample number of the peak.

[kosack]

There the time is the time since the first photon is recorded, which is not the same as the actual trigger time, as can be seen here where the top is the true time and bottom is the r1 time:

[maxnoe]

The trick I employed was just to use 1ns bins (assuming a 1ns sample width, but could be made general by looking at the actual sample frequency), and then use np.argmax() to get the sample number of the peak.

Actually, the proper way would be to evaluate the single pe spectrum around the given arrival time, sum up and sample that with the sampling frequency, but then we are fully into resimulating.

[kosack]

Exactly - there's no real point in re-doing the instrument simulation here. The best would be just to store it somewhere in in SimTel. I suspect the simulation is done only after NSB is added, otherwise it would be easy.

[kbernloehr]

I wanted to reply to this thread a while back but got distracted after having put some effort into it. Coming back to the topic, I found a bug (now fixed, see latest 'Testing' release), and tried to complete what I had started to write up in September:

Checking out (once again) the photon bunch timing in CORSIKA IACT output and its relation to telescope readout start and trigger times in array-global time units. I chose a muon simulation as an example with a peculiarity for the zero-point of the time scale (in contrast to primaries starting 'at the top of the atmosphere)'. Unfortunately, close inspection of the output also revealed a bug in code converting the readout start and trigger offset to the values recorded in simtel files. I'll come to that in a bit ...

CORSIKA output file

Corsika run header
   Run number 10 started on 2022-09-19 with version 7.7410.
   Number of showers to be simulated: 1
   Energy range 20 to 20 GeV with slope of -2.000000.
   Simulating showers in this run with up to 4.00 m random detector offsets.
   Cherenkov light observation level ist 2.150 km a.s.l.
   ...

CORSIKA IACT positions and sizes for 1 telescopes:
   x pos.: 0.000 m
   y pos.: 0.000 m
   z pos.: 12.000 m
   radius: 8.600 m

Corsika event 1 header: primary of type 6 and energy 20 GeV at  0.00 km.
   Shower direction: theta = 0.000000, phi = 0.000000 deg (azimuth N->E = 175.467000 deg)
   Core position: 0.553 m / -0.147 m
   ...

CORSIKA IACT array offsets for 1 arrays (common time offset = 393105 ns)
   x offsets: 0.553 m
   y offsets: -0.147 m

There are 1 particles arriving at ground level (encoded like photon bunches):
   Particle of type 6 (code 6001) at -0.718560 m, 2.043629 m in direction -0.000807, 0.001927, 
   arrival time 5189.491699 ns, momentum 19.617929 GeV/c, at level 1.

MC photon or photo-electron data for array 0
Telescope no. 0 in array 0 gets 24409.097656 photons in 4920 bunches (long format)

    emission at   arrival time
    [m a.s.l.]    [ns]
    -----------   --------------
    2717.449687   -387955.156250
    ...
    2150.043281   -387955.843750

Expected time between particle and Cherenkov light

The zero points for muons at observation level and for Cherenkov light differ (with FIXCHI != 0) because the particle clock starts (in this CORSIKA version, at least) at the injection point but IACT assumes it was at the top of the atmosphere, for a zero-mass particle (at v/c=1). With that correction the expected Cherenkov arrival time, in a telescope 12 m above the observation level and based on the muon arrival time is:

gnuplot> print 5189.491699 - (120e5-2150e2)/29.9792458 - 12e2*1.0001572*1.00022/29.9792458
-387955.8372

Note: 1.0001572 herer is because of the Cherenkov opening angle resulting in non-vertical photons from a vertical muon (longer path to ground), and 1.00022 is the index of refraction (slower photon propagation).

Including the actual muon (from h ~3700 m to 2150 m) being slower makes a bit of a difference of about 0.072 ns. Multiple scattering also makes the path of the actual muon a bit longer (and the actual path length is not retained) and the muon is losing energy along the path.

While the different zero-points for particles and Cherenkov light in muon simulations are unfortunate, they are - well within nanosecond accuracy - easy to explain. For shower simulations (FIXCHI=0), they are identical. No analysis should ever depend on the overall zero point.

Sim_telarray photo-electron output

MC event 100 (shower 1):
  Core offset: 2.159924, -2.257527 m

MC photon or photo-electron data for array/event 100
Photo-electrons for telescope no. 0 in array 0 with 1768 p.e. in 1764 pixels of which 69 are non-empty (with amplitudes, no NSB).
   Pixel 12: 47 p.e. starting at offset 0.
       p.e. at time -387910.53 ns, -387910.53 ns, -387910.53 ns, -387910.53 ns, -387910.53 ns, -387910.53 ns, -387910.53 ns, -387910.53 ns, -387910.53 ns, -387910.50 ns, ...
       p.e. with ampl. 1.15 pe, 1.20 pe, 1.03 pe, 0.84 pe, 0.20 pe, 0.17 pe, 1.04 pe, 1.55 pe, 1.43 pe, 1.43 pe, ...
   Pixel 14: 24 p.e. starting at offset 47.
       p.e. at time -387909.41 ns, -387909.41 ns, -387909.41 ns, -387909.41 ns, -387909.38 ns, -387909.00 ns, -387909.00 ns, -387909.00 ns, -387909.00 ns, -387909.00 ns, ...
       p.e. with ampl. 1.04 pe, 1.92 pe, 0.98 pe, 0.88 pe, 1.33 pe, 1.38 pe, 1.05 pe, 1.26 pe, 0.75 pe, 0.79 pe, ...
   Pixel 17: 28 p.e. starting at offset 71.
       p.e. at time -387910.88 ns, -387910.88 ns, -387910.88 ns, -387910.88 ns, -387909.81 ns, -387909.81 ns, -387909.81 ns, -387909.81 ns, -387909.78 ns, -387909.78 ns, ...
       p.e. with ampl. 0.66 pe, 0.83 pe, 2.05 pe, 1.28 pe, 0.98 pe, 0.26 pe, 0.16 pe, 1.77 pe, 1.04 pe, 0.56 pe, ...
   Pixel 18: 40 p.e. starting at offset 99.
       p.e. at time -387910.88 ns, -387910.88 ns, -387910.88 ns, -387910.88 ns, -387910.88 ns, -387910.88 ns, -387910.88 ns, -387910.84 ns, -387910.53 ns, -387910.53 ns, ...
       p.e. with ampl. 1.47 pe, 1.19 pe, 1.57 pe, 2.11 pe, 1.66 pe, 0.05 pe, 1.21 pe, 1.38 pe, 0.80 pe, 0.65 pe, ...
   Pixel 36: 40 p.e. starting at offset 139.
       p.e. at time -387910.41 ns, -387910.41 ns, -387910.41 ns, -387910.41 ns, -387910.38 ns, -387910.38 ns, -387910.38 ns, -387910.38 ns, -387910.28 ns, -387910.28 ns, ...
       p.e. with ampl. 0.74 pe, 1.39 pe, 1.39 pe, 1.15 pe, 0.71 pe, 0.21 pe, 1.69 pe, 1.06 pe, 1.11 pe, 0.34 pe, ...
   Pixel 38: 57 p.e. starting at offset 179.
       p.e. at time -387910.75 ns, -387910.28 ns, -387910.28 ns, -387910.28 ns, -387910.28 ns, -387910.28 ns, -387910.28 ns, -387910.28 ns, -387910.28 ns, -387910.28 ns, ...
       p.e. with ampl. 1.19 pe, 1.59 pe, 1.87 pe, 0.09 pe, 1.22 pe, 1.51 pe, 0.35 pe, 0.06 pe, 1.21 pe, 1.81 pe, ...
   ...

The photo-electrons are registered about 45 ns later than the corresponding CORSIKA photon bunches - which depends on the placement of the mirror inside the fiducial sphere and the optical paths in the ray-tracing. For the MSTs the mirror is in front of / above the fiducial sphere center by about 1.78 m and the mean distance of mirror segments from the focal plane is between 16.0 and 16.2 m, plus 2.55 offset of pixels from that, we expect a delay of about 47.5 ns for paths near the axis and 44.5 ns for paths hitting the outer mirror segments. Note that the recorded times do not include the transit time jitter.

Sim_telarray trigger time and readout

Event 100:
  Central trigger event 100
    ...
  Telescope data event 100 for telescope 1:
    Telescope event header for telescope 1:
      Local count: 1, global count: 100
      CPU time: ...
      Trigger source: 1, flags: 0500
        43 triggered sectors:      5      7     13     15     16     17     18     19     20     21     22     23     58     64     66     72    206    211    212    213    216    217    219    252    254    255    260    261    262    263    268    397    399    402    403    404    405    406    407    409    412    450    456
                     at time:  68.00  68.00  68.00  68.00  72.00  72.00  72.00  68.00  72.00  72.00  68.00  72.00  72.00  72.00  72.00  72.00  72.00  72.00  72.00  68.00  72.00  68.00  72.00  68.00  68.00  72.00  68.00  72.00  72.00  72.00  72.00  72.00  68.00  72.00  68.00  72.00  72.00  68.00  72.00  72.00  72.00  72.00  72.00 ns
      Readout starts at -387942.969 ns, rel. trigger time: 44.000 ns.
    Sample mode data (version 3) for telescope 1:
      With 25 samples, 1 gains for 1764 pixels
      Zero suppression: 0, data reduction: 0, pixel list: 0

or the same data file illustrated with a newer hessio library version, being aware that the individual values are problematic (see below):

Event 100:
  Central trigger event 100
    ...
  Telescope data event 100 for telescope 1:
    Telescope event header for telescope 1:
      Local count: 1, global count: 100
      CPU time: 1663608494.564592000, GPS time: 0.000000000
      Trigger source: 1, flags: 0500
        43 triggered sectors:      5      7     13     15     16     17     18     19     20     21     22     23     58     64     66     72    206    211    212    213    216    217    219    252    254    255    260    261    262    263    268    397    399    402    403    404    405    406    407    409    412    450    456
                     at time:  68.00  68.00  68.00  68.00  72.00  72.00  72.00  68.00  72.00  72.00  68.00  72.00  72.00  72.00  72.00  72.00  72.00  72.00  72.00  68.00  72.00  68.00  72.00  68.00  68.00  72.00  68.00  72.00  72.00  72.00  72.00  72.00  68.00  72.00  68.00  72.00  72.00  68.00  72.00  72.00  72.00  72.00  72.00 ns
      Earliest trigger is 68.000 ns after start of simulation.
      Trigger time: -387898.969 ns (-387942.969 + 44.000 ns).

Unfortunately, that example revealed a bug in the conversion of internal values to the times stored in the simtel output file. The following patch fixes that and also adds the relative position in the simulated window where the readout begins (hess/sim_conv2hess.c):

4417a4474,4475
>       /* The bin of simulated traces where the readout starts */
>       teldata->start_readout = el->sum_start;
4419c4477
<       teldata->time_readout = el->time_offset + el->sum_offset*el->interval + array_clock_offset;
---
>       teldata->time_readout = el->time_offset + el->sum_start*el->interval + array_clock_offset;
4421c4479
<       teldata->time_trg_rel = el->trigger_time - el->sum_offset*el->interval;
---
>       teldata->time_trg_rel = el->trigger_time - el->sum_start*el->interval;

Re-running the same simulation (same CORSIKA input and same random seeds in sim_telarray) with a fixed version, writing an enhanced version of the data block, shows:

Event 100:
  Central trigger event 100
    ...
  Telescope data event 100 for telescope 1:
    Telescope event header for telescope 1:
      Local count: 1, global count: 100
      CPU time: 1666275794.906050000, GPS time: 0.000000000
      Trigger source: 1, flags: 0500
        43 triggered sectors:      5      7     13     15     16     17     18     19     20     21     22     23     58     64     66     72    206    211    212    213    216    217    219    252    254    255    260    261    262    263    268    397    399    402    403    404    405    406    407    409    412    450    456
                     at time:  68.00  68.00  68.00  68.00  72.00  72.00  72.00  68.00  72.00  72.00  68.00  72.00  72.00  72.00  72.00  72.00  72.00  72.00  72.00  68.00  72.00  68.00  72.00  68.00  68.00  72.00  68.00  72.00  72.00  72.00  72.00  72.00  68.00  72.00  68.00  72.00  72.00  68.00  72.00  72.00  72.00  72.00  72.00 ns
      Earliest trigger is 68.000 ns after start of simulation.
      Readout starts at bin 11, global time -387922.969 ns, rel. trigger time: 24.000 ns (-> trigger time: -387898.969 ns).

Note that the global trigger time remains the same but the time for the start of the readout window w.r.t. that trigger time was based on the wrong variable. FlashCam with a FADC_SUM_OFFSET=6 (*4 ns) starts readout at 24 ns before the trigger, and not 24 ns after the start of the simulation window. Data which has that fixed should come with the telescope event data block in version 4 while all data with the error comes in version 3. Although the problematic sim_telarray version would work with the newer hessio library, writing version 4 telescope event blocks with a relative start of the readout at bin zero, while the fixed sim_telarray will not work with the old hessio library, thanks to the extra variable written into the struct. Thus data with fixed time conversion can be recognized, independent of the reported sim_telarray version string.

The central trigger sub-block is of no relevance for the single telescope. The placement of the trigger time and the start of the readout is all in the event header above. In this event, the trigger is 12 ns after the typical true p.e. times and the (correct) start of the readout is 12 ns before the photon arrival. Trying to superpose the true p.e. on top of the signal peaks still requires knowledge of configuration parameters like PHOTON_DELAY (12.), FADC_MHZ (250), and FADC_PULSE_SHAPE (has a peak at relative time 19.5 ns but the shaping used in the digital sum trigger has an earlier peak).

Here is a collection of a few traces from this muon ring:

Pixel 14, ch. 0: 237 250 259 262 267 289 434 539 510 448 392 350 320 311 303 290 272 260 254 250 254 255 254 260 256
Pixel 15, ch. 0: 219 217 226 231 233 230 224 217 214 211 211 207 207 209 208 206 205 203 204 205 205 216 223 224 226
Pixel 16, ch. 0: 248 240 240 241 239 232 228 236 245 252 247 237 231 226 224 227 221 219 214 210 217 224 238 252 274
Pixel 17, ch. 0: 267 278 274 265 254 264 515 814 800 677 565 485 443 425 395 359 325 297 275 264 272 271 259 250 259
Pixel 18, ch. 0: 267 270 273 273 265 273 450 696 700 600 506 445 415 383 349 318 292 270 259 264 270 270 269 273 282
Pixel 19, ch. 0: 225 226 221 218 213 217 222 224 218 215 216 217 220 237 239 233 242 251 246 253 266 272 263 252 244

The typical peak position in/near readout sample 7 (counting from 0) can be explained by the later peak in the unshaped data than in the digital sum trigger shaping.

For verification that the above relative trigger time (that is relative to the start of the recorded signal) does not correspond to bumping into the start or end of the simulated window, the FADC_BINS, FADC_SUM_BINS, and FADC_SUM_OFFSET are also needed (here 48, 25, and 6, respectively, allowing for start positions between bins 0 and 23, inclusive). Keep in mind that the times reported for triggered sectors are w.r.t. the simulated time window (earliest here: 68 ns).

Some additional, only internally used numbers can be obtained from the log file:

Telescope 1: time_offset=-387966.97 ns, nominal_delay= 40.04 ns ( -0.00, -0.00, 40.04), als=29.972409 cm/ns,
    trigger (local) time=68.00 ns, central time=-387906.48 ns, p.e. median time=-387909.24 ns,
    fadc_delay= 57.73 ns (phase delay=0.33 samples), sum_offset=6, sum_start0=11.00 [0 to 23],
    telescope_delay=-47.55 ns, photon_delay= 12.00 ns, median_time=-387909.24

None of these numbers should be needed for the analysis though.

The global start time of the readout in the problematic telescope event data (identified as coming in version 3, also identified in the TelEvent struct as having time_readout not zero but start_readout of -1) can be recovered from the earliest trigger sector (here: 68 ns = bin 17) and the FADC_SUM_OFFSET parameter - which is not part of the data block. Extracting that parameter from the history data blocks is possible as long as no splitting of CORSIKA data and merging of sim_telarray data gets involved. Since the history data block contains the original configuration lines, where parameter names may be abbreviated, case-insensitive, and the lines may include comments, that recovery is not supported out of the box. Given the known values (FADC_SUM_OFFSET=6, FADC_MHZ=250), the recovered readout start in this event would be -387942.969 - 64.0 + (17-6)4.0 = -387922.969 ns and the recovered relative trigger time would be 44.000 - (17-6)4.0 + 64.0 = 24.000 ns. (Special treatment would be needed for events with earliest sector triggered before sample 6 or after sample 48-25+6 = 29.) The new start_readout variable, reported as an impossible -1 value when reading a version 3 televent block, would be recovered as (17-6) = 11. So, despite the unfortunate mix-up, no information on the event times is actually lost.

The updated code, producing telescope event data blocks in version 4 does, of course, not need any such recovery.

[maxnoe]

Thanks Konrad for this detailed explanation. I am still a little lost...

Let's deal first with the "simple" case of fixed files with telescope event header version 4. I produced some locally.

I think the goal of this issue is to obtain a sensible time relative to the waveform array for the single photoelectrons. Super-imposing that so it fits can come later.

The final goal could be that putting the single pe pulse template at the time of the photo electrons should produce something close to the waveform.

Checking the what I readout using pyeventio, I get this for a gamma-diffuse prod6 run:

The central trigger:

{'cpu_time': (1674033392, 781159000),
 'gps_time': (1674033392, 781159000),
 'trigger_pattern': -2084306945,
 'data_pattern': -2084306945,
 'n_triggered_telescopes': 31,
 'triggered_telescopes': array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
        18, 23, 24, 25, 26, 32, 33, 34, 35, 36, 37, 38, 39, 40],
       dtype=int16),
 'trigger_times': array([14.444741 ,  2.4879844, 11.488272 , 18.05057  ,  8.763923 ,
         6.0965447,  5.084292 ,  6.325981 ,  0.       , 13.114757 ,
        35.674904 , 35.018295 , 30.048492 , 28.853045 , 27.728024 ,
         1.4273962, 62.82173  , 56.72501  , 10.9978285, 11.129453 ,
         3.6234453, 37.94255  , 46.799103 , 69.49402  , 65.99434  ,
        29.514868 , 44.905426 , 40.86889  , 14.519685 , 11.434644 ,
        14.567381 ], dtype=float32),
 'n_telescopes_with_data': 31,
 'telescopes_with_data': array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
        18, 23, 24, 25, 26, 32, 33, 34, 35, 36, 37, 38, 39, 40],
       dtype=int16),
 'teltrg_type_mask': array([2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint8),
 'teltrg_time_by_type': {},
 'plane_wavefront_compensation': {'az': 3.1415927410125732,
  'alt': 1.2217304706573486,
  'speed_of_light': 29.972414016723633}}

The telescope event header:

{'loc_count': 1,
 'glob_count': 100,
 'cpu_time': (1674033392, 781159000),
 'gps_time': (0, 0),
 'trg_source': 1,
 'list_trgsect': array([ 23,  24,  28,  89,  90,  93,  94,  96, 155, 156, 157, 158, 159,
        160, 215, 216, 217, 218, 219, 220, 271, 273, 274, 275, 276, 278,
        325, 327, 328]),
 'time_trgsect': array([40.527344, 40.283203, 40.039062, 38.57422 , 38.81836 , 38.57422 ,
        38.085938, 39.0625  , 38.57422 , 39.55078 , 38.085938, 38.085938,
        39.0625  , 38.33008 , 39.55078 , 42.23633 , 38.57422 , 39.0625  ,
        39.0625  , 38.81836 , 42.23633 , 39.79492 , 41.748047, 39.79492 ,
        39.79492 , 42.48047 , 42.96875 , 42.48047 , 43.21289 ],
       dtype=float32),
 'readout_time': -321.872802734375,
 'relative_trigger_time': 8.7890625,
 'start_readout': 30,
 'telescope_id': 1}

the photo electrons block:

{'n_pe': 2847,
 'n_pixels': 1855,
 'non_empty': 320,
 'photoelectrons': array([0, 0, 0, ..., 0, 0, 0], dtype=int32),
 'pixel_id': array([  11,   25,   25, ..., 1707, 1805, 1832], dtype=uint32),
 'time': array([-223.49152, -235.70755, -219.59154, ..., -228.85172, -217.6489 ,
        -228.97655], dtype=float32),
 'amplitude': array([1.3140458, 1.3787273, 1.1474463, ..., 1.4695556, 1.5352523,
        1.2654272], dtype=float32)}

So how do we convert the photo electron time s to be relative to the waveforms (things like pulse peak not exactly at offset=0 withstanding). I.e. the true position of the photo electron in the waveform.

[kbernloehr]

Your 'photoelectrons' array in the photo-electrons data block looks suspicious. Since the data block only covers non-empty pixels, none of these numbers should be less than 1. The size of the 'photoelectrons' array should match that of the 'pixelid' array. All other pixels can be implied to have zero photo-electrons. The read{hess,cta,simtel} output (see earlier example above) also includes the cumulative values for 'photoelectrons' which were 0, 47, 71, 99, ... for reported 47, 24, 28, 40, ... photo-electrons in pixels 12, 14, 17, 18, ... Thus the 47 photo-electrons in pixel 12 have time[0] ... time[46], the 14 p.e.in pixel 14 have time[47] ... time[70], the 28 p.e. in pixel 17 have time[71] ... time[98], and so on. Same arrangement for the amplitudes. Since these are optional, you would need to set them to 1.0 if not present, as the mean amplitude in this scale is, by definition, 1.0.

The p.e. times listed don't seem to match up with the 'readout_time' + 'relative_trigger_time' but rather some 80 or 90 nanoseconds late. My example above has p.e. listed around -387910, the pulse shape has a peak at close to 20 ns, readout starting at -387922.969 and a relative trigger time of 24 ns - which happens about 7 or 8 ns before typical p.e. time plus peak position (within 1 or 2 ns of when the reference pulse reaches 50% of its peak). That was FlashCam. In a LST sample file sitting around, the p.e. arrived between 154 and 155 ns, the pulse shape reaches 50% at 7.9 ns and its peak at 9.8 ns relative, readout started at 149.353 with the trigger 9.277 ns later (for global time 158.630 ns), just a few nanoseconds before the pulses reached 50% of their peak amplitudes. In your example, you have p.e.s arriving around -220 ns but the readout already starts at -321.8728, with the trigger 8.789 ns later. Strange ... that would need a pulse shape with a peak at around -90 ns. Please take a look at the same event with 'read_cta -S', if any discrepancy shows up.

By the way, the pulse shape that you want to fold in later, is in the pixel settings data block. Unfortunately, that does not include the time where the table started - in the case of the LST pulse shapes used at -10.1 ns (8 dyn.) and -10.6 ns (7 dyn.). Perhaps we need to complement the reference pulse shapes with these values (??)

[maxnoe]

read_cta -S shows:

Photo-electrons for telescope no. 1 in array 0 with 2847 p.e. in 1855 pixels of which 320 are non-empty (with amplitudes, no NSB).
   Pixel 11: 1 p.e. starting at offset 0.
       p.e. at time -223.49 ns
       p.e. with ampl. 1.31 pe
   Pixel 25: 2 p.e. starting at offset 1.
       p.e. at time -235.71 ns, -219.59 ns
       p.e. with ampl. 1.38 pe, 1.15 pe
   Pixel 26: 1 p.e. starting at offset 3.
       p.e. at time -232.16 ns
       p.e. with ampl. 0.54 pe
   Pixel 28: 2 p.e. starting at offset 4.
       p.e. at time -235.72 ns, -235.67 ns
       p.e. with ampl. 1.08 pe, 1.41 pe
   Pixel 29: 2 p.e. starting at offset 6.
       p.e. at time -235.65 ns, -235.65 ns
       p.e. with ampl. 0.92 pe, 1.43 pe
   Pixel 30: 1 p.e. starting at offset 8.
       p.e. at time -235.18 ns
       p.e. with ampl. 1.71 pe
   Pixel 32: 2 p.e. starting at offset 9.
       p.e. at time -232.26 ns, -232.26 ns
       p.e. with ampl. 1.32 pe, 1.71 pe
   Pixel 44: 1 p.e. starting at offset 11.
       p.e. at time -237.83 ns
       p.e. with ampl. 1.34 pe
   Pixel 45: 6 p.e. starting at offset 12.
       p.e. at time -236.79 ns, -235.74 ns, -235.73 ns, -227.56 ns, -227.55 ns, -227.55 ns
       p.e. with ampl. 0.14 pe, 1.24 pe, 1.60 pe, 1.38 pe, 1.37 pe, 1.07 pe
   Pixel 46: 2 p.e. starting at offset 18.
       p.e. at time -237.26 ns, -227.51 ns
       p.e. with ampl. 1.35 pe, 0.55 pe
   ...

Which seems consistend with what pyeventio has read out. File was produced with:

 NSHOW=10 EMIN=0.5 EMAX=5 ESLOPE=1.0 NSCAT=1 VIEWCONE=2 CSCAT=250 ./prod6_run LaPalma gamma-diffuse  

using a sim_telarray which I downloaded and compiled on 2022-11-25

[maxnoe]

File is here if you want to check: https://nextcloud.e5.physik.tu-dortmund.de/index.php/s/pmdanAF5paRtPG2

[kbernloehr]

Hi Max,

I see two problems: a) The 'photoelectrons' array for the data show above should read [1, 2, 1, 2, 2, ..., ] and not [0, 0, 0, ... ].

b) I know it is confusing but sticking to it for backward-compatibility: The photo-electron data block includes the 'telescope number' (starting at zero) while all of the type 2xxx data blocks show the telescope ID (usually starting at 1). Thus the p.e.s you showed are in telescope ID 2, which has

    Telescope event header for telescope 2:
      Local count: 1, global count: 100
      CPU time: 1674033392.781159000, GPS time: 0.000000000
      Trigger source: 1, flags: 0500
        42 triggered sectors:     11     12     15     16     19     20     23     24     40     44     45     48     49     52     53     73     77     78     81     82     85     86     89     90 
   143    144    147    148    151    152    156    209    210    213    215    216    217    218    271    273    274    276
                     at time:  40.28  40.53  39.06  39.55  38.82  39.06  39.55  39.31  41.99  39.06  41.75  38.82  39.31  39.55  39.55  43.46  40.53  41.99  39.55  40.28  39.79  40.04  40.77  41.02 
 41.99  43.95  41.26  42.24  34.67  34.67  34.42  43.70  44.19  34.42  34.18  33.94  36.13  34.67  34.67  34.67  34.91  36.38 ns
      Earliest trigger is 33.936 ns after start of simulation.
      Readout starts at bin 25, global time -249.721 ns, rel. trigger time: 9.521 ns (-> trigger time: -240.199 ns).

The trigger time still looks a couple of nanoseconds early but a complete p.e. list shows than many pixels have p.e. as early as -245 ns, and lots of them, thus a trigger only 5 ns after the first p.e.s got registered is perhaps not that surprising.

[maxnoe]

I know it is confusing but sticking to it for backward-compatibility: The photo-electron data block includes the 'telescope number' (starting at zero) while all of the type 2xxx data blocks show the telescope ID (usually starting at 1).

Ah, right, I though we change that in pyeventio, but we only do when reading in ctapipe