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

Fix dla model #1022

Merged
merged 13 commits into from
Jul 12, 2023
Merged

Fix dla model #1022

merged 13 commits into from
Jul 12, 2023

Conversation

cramirezpe
Copy link
Contributor

Fixes #1019 .

I used Garnett 2018 as a reference as it seemed more clear to me, in the issue the comparison with the old method is detailed.

@cramirezpe cramirezpe linked an issue Jul 11, 2023 that may be closed by this pull request
@cramirezpe
Copy link
Contributor Author

cramirezpe commented Jul 11, 2023

pylint fails with:
py/picca/delta_extraction/masks/dla_mask.py:8:0: E0611: No name 'voigt_profile' in module 'scipy.special' (no-name-in-module)

This makes no sense.

Also I had to change lots of variable names (e, T, mp, me) to more confusing ones because pylint complained...

@cramirezpe
Copy link
Contributor Author

I also added Naim because there are differences in the Lyman beta region between my code and his:
compare_voigt (1)

There might be an error with the constants somewhere..

@@ -28,7 +28,6 @@
GAUSSIAN_DIST = np.random.normal(size=NUM_POINTS) * np.sqrt(2)
Copy link
Collaborator

Choose a reason for hiding this comment

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

this is not used anymore, right? We should remove it

Copy link
Collaborator

Choose a reason for hiding this comment

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

we can also remove the NUM_POINTS line above

Copy link
Collaborator

Choose a reason for hiding this comment

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

and the np.random.seed(0)

@iprafols
Copy link
Collaborator

After discussing it with @cramirezpe I added a few changes to the code to further optimize. Basically, I merged the two tau functions (they were almost identical) and computed the prefactors that were independent of the DLA parameters outside of the function so that they are only computed once

I also improved the docstring adding the exact equations we use and renamed a few of the variables

From my side this PR is ok, but I think we should wait to merge until  @p-slash has review it and the discrepancies with the Lyman Beta profiles are understood

@p-slash
Copy link
Contributor

p-slash commented Jul 11, 2023

I also added Naim because there are differences in the Lyman beta region between my code and his: compare_voigt (1)

There might be an error with the constants somewhere..

The constants seem to be the same. Maybe we are seeing the thermal broadening. I think Lyb transition is weaker than Lya, so it will saturate at higher column densities. If this is the case, the choice of b (or T) will affect the profile. The other possibility is that we are seeing difference in approximations.

I will make a comment regarding what to use for oscillator strength and Einstein coefficients while replying to Michael.

z_abs: float
Redshift of the absorption
def compute_tau(lambda_, z_abs, log_nhi, lambda_t, oscillator_strength_f, gamma):
r"""Compute the optical depth for Lyman-alpha absorption.
Copy link
Contributor

Choose a reason for hiding this comment

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

This function doesn't seem specific to Lya

Copy link
Collaborator

Choose a reason for hiding this comment

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

I fixed the docstring

@Waelthus
Copy link
Contributor

I agree with @p-slash that Lyβ needs higher NHI to saturate and thus will be affected more by thermal broadening. What are we using for the b-parameter in both approaches? Are those fixed values for any DLA or something more sophisticated? How well do we know those values? And how large is the impact of the Lyβ differences?
Ideally, we should have a unified approach that is used for both analysis packages and the mocks with some characterization of its uncertainties.

@iprafols
Copy link
Collaborator

Regarding the difference in the profile, @p-slash which value are you using for the speed of light? I noticed that the reason why the tests are failing is because I changed the value of the speed of light that Cesar was using to the one stored in scipy.constants. Not sure if this would be enough to explain the differences, though

@p-slash
Copy link
Contributor

p-slash commented Jul 11, 2023

I use c=299792.458 km/s from Wikipedia

@iprafols
Copy link
Collaborator

Hi @p-slash
Cesar uses 2.9979e8, and scipy.constants.speed_of_light= 299792458.0, the same as you use. Probably it's better to use the same value so that we can rule out that the differences come from here. I'll do a commit to change this here

@Waelthus
Copy link
Contributor

Given that c=299792458 is an exact physical definition, I agree we should use that value and not some approximated one whereever possible...

@iprafols
Copy link
Collaborator

done @cramirezpe can you confirm if you are still seeing the difference you were seeing before?

@cramirezpe
Copy link
Contributor Author

All three agree (picca_old, picca_delta_extraction and qsonic) Thanks!
image

@iprafols
Copy link
Collaborator

then let's merge this and close the issue

@iprafols iprafols added this pull request to the merge queue Jul 12, 2023
Merged via the queue into master with commit 31a6c40 Jul 12, 2023
8 checks passed
@iprafols iprafols deleted the fix_dla_model branch July 12, 2023 14:07
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.

Errors in the Voigt profile functions
4 participants