Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

Time lapse corr spike b #167

Merged
merged 17 commits into from
Nov 18, 2019

Conversation

pawel21
Copy link
Collaborator

@pawel21 pawel21 commented Sep 9, 2019

I have modified code according Seiya's Nozaki idea:

  1. Updating times from slices -1 to 38 remove spike b (small jump of signal) without interpolation (Seiya Nozaki idea).

  2. Interpolation spike A is only needed for even last capacitor (lc) in the first half of the DRS4 ring.

Copy link
Contributor

@rlopezcoto rlopezcoto left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would add some comments on the code about the nature of the spikes and why they need to be interpolated
It would also be great to have some pep8-ing of the code (spaces between operator...) to improve the readability.

@@ -341,15 +340,18 @@ def do_time_lapse_corr(waveform, expected_pixel_id, local_clock_list, fc, last_t
val = waveform[gain, pixel, k] - ped_time(time_diff_ms)
waveform[gain, pixel, k] = val

# updating times from slices -1 to 38 remove spike b (small jump of signal)
# without interpolation (Seiya Nozaki idea)
posads0 = int((0 + fc[nr_module, gain, pix]) % size4drs)
if posads0+40 < 4096:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

instead of 4096, I would use some variable name, like len_drs4

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

interpolate_spike_B(waveform, gain, spike_b_position, pixel)
# The correction is only needed for even last capacitor (lc) in the
# first half of the DRS4 ring
if (fc_old[nr_module, gain, pix] + 39) % 2 == 0 and (fc_old[nr_module, gain, pix] + 39) % 1024 <= 510:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

numbers like 38, 39 could be len_drs4_window (or something like that)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

else:
for k in prange(0, 39):
for k in prange(-1, 39):
posads = int((k + fc[nr_module, gain, pix]) % size4drs)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't this run to 38?

@pawel21
Copy link
Collaborator Author

pawel21 commented Oct 3, 2019

@rlopezcoto Thanks for comments.
Now, I am at MAGIC shift, so I am a bit busy and need some time to improve code.

@FrancaCassol
Copy link
Collaborator

Hi @pawel21,
is this PR valid only before run 1575?
I wonder, if the change of firmware implies only the change of the sample range (first and last slice to be considered), should we not simply define it using variables that are assigned as trailets so to keep only one corrector class for all the runs? (the main code will verify which value to use as function of the run). Or there are other changes necessary for the new firmware?

@jsitarek
Copy link
Collaborator

jsitarek commented Nov 8, 2019

Hi @FrancaCassol
In the script to calculate pedestal coefficients we plan to do exactly what you are saying, there is is very simple and indeed it is better not to duplicate the functions.

But the dt correction function and spike interpolation are messy, because multiple conditions change, and also things can be optimized differently (e.g. the old correction required updating times from one capacitor before the actual readout to 38th capacitor including, and the code was optimized to avoid too many loops. Now the special cases when you cross the ring look different.

@pawel21
Copy link
Collaborator Author

pawel21 commented Nov 15, 2019

Hi @rlopezcoto , @FrancaCassol , @SeiyaNozaki
I have updated code according to new DRAGON firmware, and I have created function with "old" code to allow using low level calibration on "old" data before 2019/11/05 (new firmware update date).
We (me and Julian) tested code on "old" data (Run 1532) and "new" data (Run 1579), seems work OK. Code itself check run id, then according run id choosing function to low level calibration.

I made some clean-up.
I know that the code still needs cleaning, but now before observing Crab I think that the most important thing is to have working code (for new firmware) in the main repository.

Please look at my new codes and get this into the lstchain master as soon as possible.

Copy link
Contributor

@rlopezcoto rlopezcoto left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @pawel21, this looks very good, I'd recommend it to be merged as soon as you have performed the small corrections I suggested

@@ -252,24 +315,80 @@ def interpolate_pseudo_pulses(waveform, expected_pixel_id, fc, fc_old, n_modules
for pix in prange(0, n_pix):
for k in prange(0, 4):
# looking for spike A first case
abspos = int(1024 - roisize - 2 - fc_old[nr_module, gain, pix] + k * 1024 + size4drs)
abspos = int(1025 - roi_size - 2 - fc_old[nr_module, gain, pix] + k * 1024 + size4drs)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we define 1025 as size1drs + 1, with size1drs = 1024 or something like that?
Could you please also add a small explanation about this change as a comment? we had some e-mail interchange that may be lost for the future...

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

# last capacitor (lc) in the first half of the
# DRS4 ring
if ((fc_old[nr_module, gain, pix] + (roi_size-1)) % 2 == 0 and (fc_old[nr_module, gain, pix] + (roi_size-1)) % 1024 <= 511):
pixel = expected_pixel_id[nr_module * 7 + pix]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here, 511 can be size1drs/2 - 1

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

for pix in prange(0, n_pix):
for k in prange(0, 4):
# looking for spike A first case
abspos = int(1024 - roi_size - 2 -fc_old[nr_module, gain, pix]+ k * 1024 + size4drs)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same comments as above

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

roi_size = 40
for nr_module in prange(0, number_of_modules):
time_now = local_clock_list[nr_module]
for gain in prange(0, 2):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(0, n_gain)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

for nr_module in prange(0, number_of_modules):
time_now = local_clock_list[nr_module]
for gain in prange(0, 2):
for pix in prange(0, 7):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(0, n_pix)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

@@ -170,25 +191,57 @@ def time_lapse_corr(self, event, tel_id):
event : `ctapipe` event-container
tel_id : id of the telescope
"""

run_id = event.lst.tel[1].svc.configuration_id
Copy link
Collaborator

@FrancaCassol FrancaCassol Nov 18, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @pawel21,
please, change tel[1] with tel[tel_id]

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

@FrancaCassol
Copy link
Collaborator

Hi @rlopezcoto ,

for me this code is good for merging

@rlopezcoto rlopezcoto merged commit 8bdd812 into cta-observatory:master Nov 18, 2019
@pawel21 pawel21 deleted the time_lapse_corr_spike_b branch November 18, 2019 15:55
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants