diff --git a/py/picca/bin/picca_compute_fvoigt.py b/py/picca/bin/picca_compute_fvoigt.py index 901a5777c..da81b16b6 100644 --- a/py/picca/bin/picca_compute_fvoigt.py +++ b/py/picca/bin/picca_compute_fvoigt.py @@ -1,265 +1,190 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -""" -Computes the convolution term Fvoigt(k) and saves it in an ASCII file. The inputs are a DLA and a QSO catalog (both as fits binary tables). The DLA table must contain the columns "MOCKID" matching qso "THING_ID", and "Z_DLA_RSD". The QSO table must contain the columns "THING_ID", and "Z" -""" -import sys +#!/usr/bin/env python3 import argparse + import numpy as np -import scipy.integrate as integrate -import matplotlib.pyplot as plt -import astropy.io.fits as fits +from scipy.constants import speed_of_light from scipy.special import wofz -import picca.constants as constants -from picca.utils import userprint +from picca import constants -def voigt(x, sigma=1, gamma=1): - return np.real(wofz((x + 1j*gamma)/(sigma*np.sqrt(2)))) -def tau(lamb, z, N_hi): - """ lamb = lambda in A and N_hi in log10(cm^-2) """ +def voigt_profile(x, sigma=1, gamma=1): + """Return the Voigt line shape at x with Lorentzian component HWHM gamma""" + return np.real(wofz((x + 1j*gamma) / (sigma * np.sqrt(2)))) - lamb_rf = lamb/(1+z) - e = 1.6021e-19 #C - epsilon0 = 8.8541e-12 #C^2.s^2.kg^-1.m^-3 - f = 0.4164 - mp = 1.6726e-27 #kg - me = 9.109e-31 #kg - c = 2.9979e8 #m.s^-1 - k = 1.3806e-23 #m^2.kg.s^-2.K-1 - T = 1e4 #K - gamma = 6.265e8 #s^-1 - lamb_alpha = constants.ABSORBER_IGM["LYA"] #A - Deltat_lamb = lamb_alpha/c*np.sqrt(2*k*T/mp) #A - - a = gamma/(4*np.pi*Deltat_lamb)*lamb_alpha**2/c*1e-10 - u = (lamb_rf - lamb_alpha)/Deltat_lamb - H = voigt(u, np.sqrt(1/2), a) - - absorb = np.sqrt(np.pi)*e**2*f*lamb_alpha**2*1e-10/(4*np.pi*epsilon0*me*c**2*Deltat_lamb)*H - #10^N_hi en cm^-2 et absorb en m^2 - return 10**N_hi*1e4*absorb - -def profile_voigt_lambda(x, z, N_hi): - t = tau(x, z, N_hi).astype(float) - return np.exp(-t) - -def profile_lambda_to_r(lamb, profile_lambda, fidcosmo): #pour Lyman_alpha --> sinon mettre une autre raie - z = lamb/constants.ABSORBER_IGM["LYA"] - 1 - r = fidcosmo.get_r_comov(z) - rr = np.linspace(r[0], r[-1], r.size) - profile_r = np.interp(rr,r,profile_lambda) #pour reavoir un echantillonage lineaire - return rr, profile_r - -def fft_profile(profile, dx): # non normalise - n = profile.size - tmp = (1-profile) - ft_profile = dx*np.fft.fftshift(np.fft.fft(tmp)) - k = np.fft.fftshift(np.fft.fftfreq(n, dx))*(2*np.pi) - return ft_profile, k - -def lambda_to_r(lamb, profile_lambda, fidcosmo): - """ - Converts a profile as a function of wavelength to a profile as a function of r in Mpc/h - """ - z = lamb/constants.ABSORBER_IGM["LYA"] - 1 - r = fidcosmo.get_r_comov(z) - rr = np.linspace(r[0], r[-1], r.size) - profile_lambda = profile_lambda*fidcosmo.get_hubble(z)*constants.ABSORBER_IGM["LYA"]/3e5 - profile_r = np.interp(rr,r,profile_lambda) - return rr, profile_r +def get_voigt_profile_wave(wave, z, N_hi): + """Compute voigt profile at input redshift and column density + on an observed wavelength grid -def compute_dla_prob_per_nhi(wavelength,nhi,dla,qso,dnhi): + Parameters + ---------- + wave : array + Observed wavelength + z : float + Redshift + N_hi : float + Log10 column density + + Returns + ------- + array + Flux corresponding to a voigt profile """ - Computes the probability of finding a DLA with NHI= n+-dn and wavelength lambda +- dlambda - - wavelength 1D array - - nhi value of column density bin center log10(cm-2) - - dla the DLA catalog (table) - - qso the QSO catalog (table) - - dnhi bin width log10(cm-2) + # wave = lambda in A and N_HI in log10 and 10**N_hi in cm^-2 + wave_rf = wave / (1 + z) + + e = 1.6021e-19 # C + epsilon0 = 8.8541e-12 # C^2.s^2.kg^-1.m^-3 + f = 0.4164 + mp = 1.6726e-27 # kg + me = 9.109e-31 # kg + c = speed_of_light # m.s^-1 + k = 1.3806e-23 # m^2.kg.s^-2.K-1 + T = 1e4 # K + gamma = 6.265e8 # s^-1 + + wave_lya = constants.ABSORBER_IGM["LYA"] # A + Deltat_wave = wave_lya / c * np.sqrt(2 * k * T / mp) # A + + a = gamma / (4 * np.pi * Deltat_wave) * wave_lya**2 / c * 1e-10 + u = (wave_rf - wave_lya) / Deltat_wave + H = voigt_profile(u, np.sqrt(1/2), a) + + absorbtion = np.sqrt(np.pi) * e**2 * f * wave_lya**2 * 1e-10 * H + absorbtion /= (4 * np.pi * epsilon0 * me * c**2 * Deltat_wave) + + # 10^N_hi in cm^-2 and absorb in m^2 + tau = (10**N_hi * 1e4 * absorbtion).astype(float) + return np.exp(-tau) + + +def profile_wave_to_comov_dist(wave, profile_wave, cosmo, differential=False): + """Convert profile from a function of wavelength to a function of comoving distance + + Parameters + ---------- + wave : array + Observed wavelength + profile_wave : array + Profile as a function of wavelength + cosmo : picca.constants.Cosmo + Cosmology object + differential : bool, optional + Whether we have to account for a dlambda / dr factor, by default False + + Returns + ------- + (array, array) + (comoving distance grid, profile as a function of comoving distance) """ + z = wave / constants.ABSORBER_IGM["LYA"] - 1 + comov_dist = cosmo.get_r_comov(z) + lin_spaced_comov_dist = np.linspace(comov_dist[0], comov_dist[-1], comov_dist.size) - dla_lamb = (1 + dla['Z_DLA_RSD'])*constants.ABSORBER_IGM["LYA"] - dla_nhi = dla['N_HI_DLA'] - qso_lamb = (1 + qso['Z'])*constants.ABSORBER_IGM["LYA"] - dwave=np.gradient(wavelength) + profile_comov_dist = profile_wave + if differential: + # We are in the f(lambda)dlambda = f(r)dr case, + # so need to account for dlambda / dr = H(z) * Lambda_Lya / c + profile_comov_dist = cosmo.get_hubble(z) * constants.ABSORBER_IGM["LYA"] / speed_of_light - wbins = np.zeros(wavelength.size+1) - wbins[:-1] = wavelength-dwave/2 - wbins[-1] = wavelength[-1]+dwave[-1]/2. + profile_comov_dist = np.interp(lin_spaced_comov_dist, comov_dist, profile_comov_dist) - ndla,_ = np.histogram(dla_lamb[np.abs(dla_nhi-nhi) wave) - if nqso>0 : - f[i] = ndla[i]/float(nqso)/dnhi/dwave[i] - return f +def fft_profile(profile, dx): + """Compute Fourier transform of a voigt profile -def compute_dla_prob(wavelength, NHI, dla, qso, weight): + Parameters + ---------- + profile : array + Input voigt profile in real space (function of comoving distance) + dx : float + Comoving distance bin size + Returns + ------- + (array, array) + (wavenumber grid, voigt profile in Fourier space) """ - Computes the probability of finding a DLA - - wavelength is 1 1D array - - NHI is a 1D array of N(HI) in units of log10(cm^-2) - - dla the DLA catalog (table) - - qso the QSO catalog (table) - - weights as a function of wavelength - """ - dnhi=np.gradient(NHI) - mean_density = np.zeros(NHI.size) - for i,nhi in enumerate(NHI): - userprint("compute prob(NHI={}) ({}/{})".format(nhi,(i+1),NHI.size)) - f = compute_dla_prob_per_nhi(wavelength,nhi=nhi,dla=dla,qso=qso,dnhi=dnhi[i]) - mean_density[i] = np.sum(f*weight)/np.sum(weight) - return mean_density - -def main(cmdargs) : - - parser = argparse.ArgumentParser( - formatter_class=argparse.ArgumentDefaultsHelpFormatter, - description=('Computes the convolution term Fvoigt(k) and saves it in ' - 'an ASCII file. The inputs are a DLA and a QSO catalog ' - '(both as fits binary tables). The DLA table must contain ' - 'the columns "MOCKID" matching qso "THING_ID", and ' - '"Z_DLA_RSD". The QSO table must contain the columns ' - '"THING_ID", and "Z"')) - - parser.add_argument( - '--dla-catalog', - type=str, - default=None, - required=True, - help=('DLA catalog fits file , like ' - '/project/projectdirs/desi/mocks/lya_forest/develop/saclay/v4.4/v4.4.3/master_DLA.fits')) - - parser.add_argument( - '--drq-catalog', - type=str, - default=None, - required=True, - help=('DRQ catalog fits file, like ' - '/project/projectdirs/desi/mocks/lya_forest/develop/saclay/v4.4/v4.4.3/eboss-0.2/zcat_desi_drq.fits')) - - parser.add_argument( - '--weight-vs-wavelength', - type=str, - default=None, - required=False, - help=('sum of delta weight as a function of wavelength (two columns ' - 'ASCII file)')) - - parser.add_argument( - '-o', '--output', - type=str, - default=None, - required=True, - help='FVoigt as a function of k in h/Mpc (two columns ASCII file)') - - parser.add_argument('-d','--debug', action='store_true', default=True) - - parser.add_argument('-p','--plot', action='store_true', - help="show some plots") - - - args = parser.parse_args(cmdargs) - - if args.debug: - userprint("read DLA catalog") - dla = fits.open(args.dla_catalog)[1].data - - if args.debug: - userprint("read DRQ catalog") - qso = fits.open(args.drq_catalog)[1].data - - if args.debug: - userprint("only keep DLAs in DRQ quasars LOS") - dla = dla[:][np.in1d(dla['MOCKID'], qso['THING_ID'])] - - coarse_wavelength = np.arange(3000, 8000, 100) - - if args.weight_vs_wavelength is None: - filename = "/global/common/software/desi/users/jguy/igmhub/code_stage_lbl/build_Fvoigt/data/weight_lambda.txt" - userprint("WARNING: Hardcoded weight vs wavelength file {}".format(filename)) - args.weight_vs_wavelength = filename - - if args.weight_vs_wavelength is not None: - if args.debug: - userprint("read weights vs wave") - tmp = np.loadtxt(args.weight_vs_wavelength) - weight = np.interp(coarse_wavelength, tmp[:, 0], tmp[:, 1]) - else : - weight = np.ones(coarse_wavelength.shape) - - zdla = np.mean(dla['Z_DLA_RSD']) - NHI = np.linspace(17, 22, (22 - 17)/0.1 + 1) - - # probability of finding a DLA at NHI in a QSO LOS (per A and per unit NHI) - # averaged over wavelength, using the provided wavelength weight - prob = compute_dla_prob(coarse_wavelength, NHI, dla=dla, qso=qso, weight=weight) - - ii=(prob > 0) - prob=prob[ii] - NHI=NHI[ii] - - if args.plot : - plt.figure("prob") - plt.plot(NHI,prob,"o-") - plt.xlabel("NHI (log10(cm-2))") - plt.ylabel("prob(DLA) per unit NHI per A") - - # now use a finer wavelength grid - wavelength = np.arange(3000, 8000, 1.) - # conversion A -> Mpc/h - fidcosmo = constants.Cosmo(Om=0.3) - r_wave = fidcosmo.get_r_comov(wavelength/constants.ABSORBER_IGM["LYA"] - 1) - # linear grid of Mpc/h (need to convert to linear grid for the FFT) - r_lin = np.linspace(r_wave[0],r_wave[-1],r_wave.size) - - for i in range(NHI.size): - - if args.debug : userprint("compute dF/dNHI for NHI={} ({}/{})".format(NHI[i],(i+1),NHI.size)) - profile = profile_voigt_lambda(wavelength, zdla, NHI[i]) - profile_r = np.interp(r_lin,r_wave,profile) # interpolation to linear r grid - # r is in Mpc h^-1 --> k in h*Mpc^-1 - ft_profile, k = fft_profile(profile_r, np.abs(r_lin[1]-r_lin[0])) - ft_profile = np.abs(ft_profile) - if i == 0: - df = np.array([ft_profile*prob[i]]) - else: - df = np.concatenate((df, np.array([ft_profile*prob[i]]))) - - if args.debug : userprint("compute F(k)=int dF/dNHI * dNHI") + # not normalized + size = profile.size + ft_profile = dx * np.fft.rfft(profile - 1) + k = np.fft.rfftfreq(size, dx) * (2 * np.pi) + + return k, np.abs(ft_profile) + + +def main(): + parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter, + description=('Compute FVoigt profile')) + + parser.add_argument('-o', '--output', type=str, default=None, + help=('Name and path of file to write FVoigt profile to')) + + parser.add_argument('--weights', type=str, default=None, + help=('File with total weight as a function of observed wavelength')) + + parser.add_argument('--cddf', type=str, default=None, + help=('File with column density distribution function (CDDF)')) + + parser.add_argument('--dla-prob', type=str, default=None, + help=('File with probability of finding DLAs in a forest ' + 'as a function of observed wavelength')) + + parser.add_argument('--z-dla', type=float, default=None, + help=('Mean DLA redshift')) + + parser.add_argument('--normalize', action='store_true', required=False, + help=('Whether the output FVoigt function is normalized')) + + parser.add_argument('--positive-fvoigt', action='store_true', required=False, + help=('Whether the output FVoigt function should be positive')) + + args = parser.parse_args() + + cddf_NHI, dN_NHI = np.loadtxt(args.cddf) + weights_wave = np.loadtxt(args.weights) + dla_prob_wave = np.loadtxt(args.dla_prob) + + fidcosmo = constants.Cosmo(Om=0.3147) + + comov_dist_prob, dla_prob_comov_dist = profile_wave_to_comov_dist( + dla_prob_wave[0], dla_prob_wave[1], fidcosmo, differential=True) + comov_dist_weights, weights_comov_dist = profile_wave_to_comov_dist( + weights_wave[:, 0], weights_wave[:, 1], fidcosmo) + + weight_interp = np.interp( + comov_dist_prob, comov_dist_weights, weights_comov_dist, left=0, right=0) + mean_density = np.average(dla_prob_comov_dist, weights=weight_interp) + + wave = np.arange(2000, 8000, 1) # TODO this grid may be too sparse + integrand = np.empty((dN_NHI.size, wave.size // 2 + 1)) + for index, NHI in enumerate(dN_NHI): + profile_wave = get_voigt_profile_wave(wave, args.z_dla, NHI) + profile_wave /= np.mean(profile_wave) + + # r is in Mpc h^-1 --> k (from tf) will be in (Mpc h^-1)^-1 = h Mpc^-1 :) + comov_dist, profile_comov_dist = profile_wave_to_comov_dist(wave, profile_wave, fidcosmo) + k, fourier_profile = fft_profile(profile_comov_dist, np.abs(comov_dist[1] - comov_dist[0])) + + integrand[index] = fourier_profile * mean_density * cddf_NHI[index] + Fvoigt = np.zeros(k.size) - for i in range(k.size): - Fvoigt[i] = integrate.trapz(df[:,i], NHI) - Fvoigt = Fvoigt/Fvoigt[k.size//2] - - ii=(k>=0) - k=k[ii] - Fvoigt=Fvoigt[ii] - - if args.debug : userprint("save in {}".format(args.output)) - np.savetxt(args.output, np.array([k,Fvoigt]).T) - - - if args.plot : - plt.figure("fk") - plt.plot(k,Fvoigt,label="this model") - plt.plot(k,np.exp(-k*10),"--",label="exp(-k*(10 Mpc/h))") - plt.plot(k,np.exp(-k*6),":",label="exp(-k*(6 Mpc/h))") - plt.xscale("log") - plt.xlabel("k (h/Mpc)") - plt.ylabel("F(k)") - plt.legend(loc="upper right") - plt.grid() - plt.show() - -if __name__=='__main__': - cmdargs=sys.argv[1:] - main(cmdargs) + for index in range(k.size): + Fvoigt[index] = np.trapz(integrand[:, index], dN_NHI) + + Fvoigt = Fvoigt[k > 0] + k = k[k > 0] + + if args.normalize: + Fvoigt /= Fvoigt[0] + if not args.positive_fvoigt: + Fvoigt = -Fvoigt + + np.savetxt(args.output, np.c_[k, Fvoigt], header='k[hMpc^-1] FVoigt') + + +if __name__ == '__main__': + main() diff --git a/py/picca/bin/picca_compute_fvoigt_hcd.py b/py/picca/bin/picca_compute_fvoigt_hcd.py deleted file mode 100644 index ff032f26c..000000000 --- a/py/picca/bin/picca_compute_fvoigt_hcd.py +++ /dev/null @@ -1,186 +0,0 @@ -#!/usr/bin/env python3 -import argparse - -import numpy as np -from scipy.constants import speed_of_light -from scipy.special import wofz - -from picca import constants -from picca.delta_extraction.masks.dla_mask import ( - GAMMA_LYA, - LAMBDA_LYA, - OSCILLATOR_STRENGTH_LYA, - dla_profile, -) - - -def profile_wave_to_comov_dist(wave, profile_wave, cosmo, differential=False): - """Convert profile from a function of wavelength to a function of comoving distance - - Parameters - ---------- - wave : array - Observed wavelength - profile_wave : array - Profile as a function of wavelength - cosmo : picca.constants.Cosmo - Cosmology object - differential : bool, optional - Whether we have to account for a dlambda / dr factor, by default False - - Returns - ------- - (array, array) - (comoving distance grid, profile as a function of comoving distance) - """ - z = wave / constants.ABSORBER_IGM["LYA"] - 1 - comov_dist = cosmo.get_r_comov(z) - lin_spaced_comov_dist = np.linspace(comov_dist[0], comov_dist[-1], comov_dist.size) - - profile_comov_dist = profile_wave - if differential: - # We are in the f(lambda)dlambda = f(r)dr case, - # so need to account for dlambda / dr = H(z) * Lambda_Lya / c - profile_comov_dist = ( - cosmo.get_hubble(z) * constants.ABSORBER_IGM["LYA"] / speed_of_light - ) - - profile_comov_dist = np.interp( - lin_spaced_comov_dist, comov_dist, profile_comov_dist - ) - - return lin_spaced_comov_dist, profile_comov_dist - - -def fft_profile(profile, dx): - """Compute Fourier transform of a voigt profile - - Parameters - ---------- - profile : array - Input voigt profile in real space (function of comoving distance) - dx : float - Comoving distance bin size - - Returns - ------- - (array, array) - (wavenumber grid, voigt profile in Fourier space) - """ - # not normalized - size = profile.size - ft_profile = dx * np.fft.rfft(profile - 1) - k = np.fft.rfftfreq(size, dx) * (2 * np.pi) - - return k, np.abs(ft_profile) - - -def main(): - parser = argparse.ArgumentParser( - formatter_class=argparse.ArgumentDefaultsHelpFormatter, - description=("Compute FVoigt profile"), - ) - - parser.add_argument( - "-o", - "--output", - type=str, - required=True, - help=("Name and path of file to write FVoigt profile to"), - ) - - parser.add_argument( - "--weights", - type=str, - required=True, - help=("File with total weight as a function of observed wavelength"), - ) - - parser.add_argument( - "--cddf", - type=str, - required=True, - help=("File with column density distribution function (CDDF)"), - ) - - parser.add_argument( - "--dla-prob", - type=str, - required=True, - help=( - "File with probability of finding DLAs in a forest " - "as a function of observed wavelength" - ), - ) - - parser.add_argument( - "--z-dla", type=float, required=True, help=("Mean DLA redshift") - ) - - parser.add_argument( - "--normalize", - action="store_true", - required=False, - help=("Whether the output FVoigt function is normalized"), - ) - - parser.add_argument( - "--positive-fvoigt", - action="store_true", - required=False, - help=("Whether the output FVoigt function should be positive"), - ) - - args = parser.parse_args() - - cddf_NHI, dN_NHI = np.loadtxt(args.cddf) - weights_wave = np.loadtxt(args.weights) - dla_prob_wave = np.loadtxt(args.dla_prob) - - fidcosmo = constants.Cosmo(Om=0.3147) - - comov_dist_prob, dla_prob_comov_dist = profile_wave_to_comov_dist( - dla_prob_wave[0], dla_prob_wave[1], fidcosmo, differential=True - ) - comov_dist_weights, weights_comov_dist = profile_wave_to_comov_dist( - weights_wave[:, 0], weights_wave[:, 1], fidcosmo - ) - - weight_interp = np.interp( - comov_dist_prob, comov_dist_weights, weights_comov_dist, left=0, right=0 - ) - mean_density = np.average(dla_prob_comov_dist, weights=weight_interp) - - wave = np.arange(2000, 8000, 1) # TODO this grid may be too sparse - integrand = np.empty((dN_NHI.size, wave.size // 2 + 1)) - for i, NHI in enumerate(dN_NHI): - profile_wave = dla_profile(wave, args.z_dla, NHI) - profile_wave /= np.mean(profile_wave) - - # r is in Mpc h^-1 --> k (from tf) will be in (Mpc h^-1)^-1 = h Mpc^-1 :) - comov_dist, profile_comov_dist = profile_wave_to_comov_dist( - wave, profile_wave, fidcosmo - ) - k, fourier_profile = fft_profile( - profile_comov_dist, np.abs(comov_dist[1] - comov_dist[0]) - ) - - integrand[i] = fourier_profile * mean_density * cddf_NHI[i] - - Fvoigt = np.zeros(k.size) - for i in range(k.size): - Fvoigt[i] = np.trapz(integrand[:, i], dN_NHI) - - Fvoigt = Fvoigt[k > 0] - k = k[k > 0] - - if args.normalize: - Fvoigt /= Fvoigt[0] - if not args.positive_fvoigt: - Fvoigt = -Fvoigt - - np.savetxt(args.output, np.c_[k, Fvoigt]) - - -if __name__ == "__main__": - main() diff --git a/requirements.txt b/requirements.txt index d7503b970..e7cf4c6b1 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ numpy>=2.0.0 -scipy>=1.14.0 +scipy==1.15.3 iminuit>=2.26.0 healpy>=1.17.1 fitsio>=1.2.4