- 
                Notifications
    You must be signed in to change notification settings 
- Fork 20
Flat Field coefficients computation + the rise and decay times #221
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
base: main
Are you sure you want to change the base?
Conversation
…requered changes in containers files. And the script to compute raise and decay times. Cleaned branch.
| Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@            Coverage Diff             @@
##             main     #221      +/-   ##
==========================================
+ Coverage   51.79%   51.82%   +0.03%     
==========================================
  Files          78       78              
  Lines        6505     6510       +5     
==========================================
+ Hits         3369     3374       +5     
  Misses       3136     3136              ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
 | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi @anmikhno !
Thanks a lot for this PR !
I would like to help streamlining some parts of your main script analysis_photostat_simple.py before approving it. I'll contribute to the PR shortly.
| (self.gainLG * mask, gainLG_err, gainLG_err) | ||
| ).T | ||
| self._results.is_valid = mask | ||
| log.info("Trying to write charges") | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe we should make it clear everywhere that these 2 new fields are for the high-gain channel ? E.g. by renaming them charge_hg and charge_hg_std everywhere ?
| ] | ||
|  | ||
| print("[DEBUG] Running command:", " ".join(cmd)) | ||
| subprocess.run(cmd, check=True) | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We may try to find a more pythonic way, e.g. by invoking directly the API found in the script gain_PhotoStat_computation.py ?
| print(f"[INFO] {filename_ps} not found, running gain_PhotoStat_computation.py...") | ||
|  | ||
| gain_script = os.path.expanduser( | ||
| "~/local/src/nectarchain/src/nectarchain/" | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We should find another way, and avoid using such an absolute path.
| else: | ||
| print(f"[INFO] File {run_path} already exists, skipping computation.") | ||
|  | ||
| print(f"[INFO] Starting analysis on {filename_ps}") | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
print -> logging, to be implemented.
|  | ||
|  | ||
| def optimize_with_outlier_rejection(sigma, data): | ||
| def define_delete_out(sigma, data): | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Avoid nested methods ?
|  | ||
| # Define the least-squares function | ||
| def LSQ_wrap(a0, a1, a2, a3): | ||
| a4 = a3 # changed | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It may be good to add a comment here, explicitly stating that a4 = a3 means that we assume a symmetric 2D-Gaussian (is that correct, by the way ? Seems so, looking at https://ctapipe.readthedocs.io/en/latest/api/ctapipe.image.toymodel.Gaussian.html#ctapipe.image.toymodel.Gaussian)
| print(f"Max residual: {max_residual*100:.2f}%") | ||
|  | ||
| dict_missing_pix["rejected_outliers"] = ( | ||
| 1855 | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let's use ctapipe_io_nectarcam.constants.N_PIXELS instead.
|  | ||
| # main | ||
| if __name__ == "__main__": | ||
| camera = CameraGeometry.from_name("NectarCam-003").transform_to( | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let's discover the camera geometry from the data, e.g. using EventSource.subarray.tel[0].camera.geometry instead.
| ) | ||
| low_gain = Field(type=np.ndarray, dtype=np.float64, ndim=2, description="low gain") | ||
| pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids") | ||
| charge = Field(type=np.ndarray, dtype=np.float64, ndim=1, description="charge") | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would suggest to create a PhotostatContainer to let the gain container in a high abstraction level, the charge and charge_std fileds are not needed in the SPEfitContainer which inherits from the gaincontainer
| high_gain=np.zeros((self.npixels, 3)), | ||
| low_gain=np.zeros((self.npixels, 3)), | ||
| pixels_id=self._pixels_id, | ||
| charge=np.zeros((self.npixels, 1)), | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
and here return the PhotoStatContainer
| Hi @anmikhno, thanks a lot for your work, I just quickly let a comment about the abstraction of  | 
| # ctapipe modules | ||
| from ctapipe.visualization import CameraDisplay | ||
| from iminuit import Minuit | ||
| from jacobi import propagate | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This creates a new dependency for nectarchain. Is it strictly required ?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If I understand scipy/scipy#17059 correctly, this has been merged into scipy directly since version 1.15.0. A recent installation of nectarchain comes with at least scipy v1.16.0. Could we try to use the implementation from scipy instead of adding another dependency ?
…requered changes in containers files. And the script to compute raise and decay times. Cleaned branch.
The flat field coefficients are computed using the Gaussian fit, on top pf that the outliers are removed to improve the fit. The script generated pdf with the plots and a table with flat-field coefficients computed per pixel using the fit (model method) and flat field coefficients computed using the simple formula (n_pe_per_pix/mean(n_pe_over_camera))
Cleaned branch.