I'm doing DM->pbar->propagation simulations and have several questions on how to configure a galdef file.
Currently my setup looks like:
Code:
Title = DM->pbar
######################################################
# Dimensions and grid properties
######################################################
n_spatial_dimensions = 3
x_min = -11.0 min x
x_max = +11.0 max x
dx = 1.0 delta x
y_min = -11.0 min y
y_max = +11.0 max y
dy = 1.0 delta y
z_min = -11.0 min z
z_max = +11.0 max z
dz = 1.0 delta z
p_min = 10 min momentum (MV)
p_max = 1000000 max momentum
p_factor = 2.50 momentum factor
Ekin_min = 1.e1 min kinetic energy per nucleon (MeV)
Ekin_max = 1.e7 max kinetic energy per nucleon
Ekin_factor = 2.2 kinetic energy per nucleon factor
p_Ekin_grid = Ekin p||Ekin alignment
######################################################
# Propagation properties
######################################################
prop_x = 1 1=propagate in x (2D,3D)
prop_y = 1 1=propagate in y (3D)
prop_z = 1 1=propagate in z (3D)
prop_p = 1 1=propagate in momentum
######################################################
# Timestamp and iteration properties
######################################################
start_timestep =1.0e9
end_timestep =1.0e2
timestep_factor =0.25
timestep_repeat =20 number of repeats per timestep in timetep_mode=1
timestep_repeat2 =0 number of timesteps in timetep_mode=2
timestep_print =10000 number of timesteps between printings
timestep_diagnostics =10000 number of timesteps between diagnostics
control_diagnostics =0 control detail of diagnostics
network_iterations = 1000 number of iterations of entire network
######################################################
# Diffusion properties
######################################################
D0_xx =5.80e28 diffusion coefficient at reference rigidity
D_rigid_br =4.0e3 reference rigidity for diffusion coefficient in MV
D_g_1 = 0.33 diffusion coefficient index below reference rigidity
D_g_2 = 0.33 diffusion coefficient index above reference rigidity
diff_reacc =1 0=no reacc.; 1,2=incl.diff.reacc.; -1==beta^3 Dxx; 11=Kolmogorov+damping; 12=Kraichnan+damping
v_Alfven =30. Alfven speed in km s-1
######################################################
# Dark Matter properties
######################################################
DM_positrons = 0 1=compute DM positrons
DM_electrons = 0 1=compute DM electrons
DM_antiprotons = 1 1=compute DM antiprotons
DM_gammas = 0 1=compute DM gammas
DM_double0 = 25.0 Rc, core radius, kpc
DM_double1 = 0.3 rho0 = local DM mass density, GeV cm-3
DM_double2 = 5.0e+02 DM mass, GeV
DM_double3 = 0.
DM_double4 = 0.
DM_double5 = 1.97e+07 a parameter
DM_double6 = 2.00e+00 b parameter
DM_double7 = 1.71e+01 c parameter
DM_double8 = 1.39e-01 d parameter
DM_double9 = 3.e-25 <cross_sec*V>-thermally overaged, cm3 s-1
DM_int0 = 0 NFW profile
DM_int1 = 777 KK production
DM_int2 = 1
DM_int3 = 1
DM_int4 = 1
DM_int5 = 1
DM_int6 = 1
DM_int7 = 1
DM_int7 = 1
DM_int8 = 1
DM_int9 = 1
######################################################
# Nuclei properties
######################################################
max_Z = 0 maximum number of nucleus Z listed
use_Z_1 = 0
iso_abundance_01_001 = 1.06e+06 H
iso_abundance_01_002 = 0. 34.8
iso_abundance_02_003 = 9.033 He
iso_abundance_02_004 = 7.199e+04
iso_abundance_03_006 = 0 Li
iso_abundance_03_007 = 0
iso_abundance_04_009 = 0 Be
iso_abundance_05_010 = 0 B
iso_abundance_05_011 = 0
iso_abundance_06_012 = 2819 C
iso_abundance_06_013 = 5.268e-07
iso_abundance_07_014 = 182.8 N
iso_abundance_07_015 = 5.961e-05
iso_abundance_08_016 = 3822 O
iso_abundance_08_017 = 6.713e-07
iso_abundance_08_018 = 1.286
iso_abundance_09_019 = 2.664e-08 F
iso_abundance_10_020 = 312.5 Ne
iso_abundance_10_021 = 0.003556
iso_abundance_10_022 = 100.1
iso_abundance_11_023 = 22.84 Na
iso_abundance_12_024 = 658.1 Mg
iso_abundance_12_025 = 82.5
iso_abundance_12_026 = 104.7
iso_abundance_13_027 = 76.42 Al
iso_abundance_14_028 = 725.7 Si
iso_abundance_14_029 = 35.02
iso_abundance_14_030 = 24.68
iso_abundance_15_031 = 4.242 P
iso_abundance_16_032 = 89.12 S
iso_abundance_16_033 = 0.3056
iso_abundance_16_034 = 3.417
iso_abundance_16_036 = 0.0004281
iso_abundance_17_035 = 0.7044 Cl
iso_abundance_17_037 = 0.001167
iso_abundance_18_036 = 9.829 Ar
iso_abundance_18_038 = 0.6357
iso_abundance_18_040 = 0.001744
iso_abundance_19_039 = 1.389 K
iso_abundance_19_040 = 3.022
iso_abundance_19_041 = 0.0003339
iso_abundance_20_040 = 51.13 Ca
iso_abundance_20_041 = 1.974
iso_abundance_20_042 = 1.134e-06
iso_abundance_20_043 = 2.117e-06
iso_abundance_20_044 = 9.928e-05
iso_abundance_20_048 = 0.1099
iso_abundance_21_045 = 1.635 Sc
iso_abundance_22_046 = 5.558 Ti
iso_abundance_22_047 = 8.947e-06
iso_abundance_22_048 = 6.05e-07
iso_abundance_22_049 = 5.854e-09
iso_abundance_22_050 = 6.083e-07
iso_abundance_23_050 = 1.818e-05 V
iso_abundance_23_051 = 5.987e-09
iso_abundance_24_050 = 2.873 Cr
iso_abundance_24_052 = 8.065
iso_abundance_24_053 = 0.003014
iso_abundance_24_054 = 0.4173
iso_abundance_25_053 = 6.499 Mn
iso_abundance_25_055 = 1.273
iso_abundance_26_054 = 49.08 Fe
iso_abundance_26_056 = 697.7
iso_abundance_26_057 = 21.67
iso_abundance_26_058 = 3.335
iso_abundance_27_059 = 2.214 Co
iso_abundance_28_058 = 28.88 Ni
iso_abundance_28_060 = 11.9
iso_abundance_28_061 = 0.5992
iso_abundance_28_062 = 1.426
iso_abundance_28_064 = 0.3039
######################################################
# Other properties
######################################################
E_gamma_min = 100. min gamma-ray energy (MeV)
E_gamma_max = 1000. max gamma-ray energy (MeV)
E_gamma_factor = 2.0 gamma-ray energy factor
integration_mode = 0 integr.over part.spec.: =1-old E*logE; =0-PL analyt.
nu_synch_min = 1.0e6 min synchrotron frequency (Hz)
nu_synch_max = 1.0e10 max synchrotron frequency (Hz)
nu_synch_factor = 2.0 synchrotron frequency factor
long_min = 0.25 gamma-ray intensity skymap longitude minimum (deg); 0 -automatic binning required to get correct results!
long_max =359.75 gamma-ray intensity skymap longitude maximum (deg); 360 -automatic binning
lat_min =-89.75 gamma-ray intensity skymap latitude minimum (deg); -90 -automatic binning
lat_max =+89.75 gamma-ray intensity skymap latitude maximum (deg); +90 -automatic binning
d_long = 5.0 gamma-ray intensity skymap longitude binsize (deg)
d_lat = 5.0 gamma-ray intensity skymap latitude binsize (deg)
healpix_order = 4 order for healpix skymaps. 7 gives ~0.5 deg and it changes by an order of 2
lat_substep_number = 1 latitude bin splitting (0,1=no split, 2=split in 2...)
LoS_step = 0.01 kpc, Line of Sight (LoS) integration step
LoS_substep_number = 1 number of substeps per LoS integration step (0,1=no substeps)
damping_p0 = 1.e6 MV -some rigidity (where CR density is low)
damping_const_G = 0.02 a const derived from fitting B/C
damping_max_path_L = 3.e21 Lmax~1 kpc, max free path
convection =0 1=include convection
v0_conv =0. km s-1 v_conv=v0_conv+dvdz_conv*dz
dvdz_conv =10. km s-1 kpc-1 v_conv=v0_conv+dvdz_conv*dz
nuc_rigid_br = 9.0e3 reference rigidity for nucleus injection index in MV
nuc_g_1 =1.98 nucleus injection index below reference rigidity
nuc_g_2 =2.42 nucleus injection index index above reference rigidity
inj_spectrum_type = rigidity rigidity||beta_rig||Etot nucleon injection spectrum type
electron_g_0 =1.60 electron injection index below electron_rigid_br0
electron_rigid_br0 =4.0e3 reference rigidity0 for electron injection index in MV
electron_g_1 =2.42 2.3 2.54 electron injection index below reference rigidity
electron_rigid_br =1.0e9 reference rigidity for electron injection index in MV
electron_g_2 =5.0 electron injection index index above reference rigidity
He_H_ratio = 0.11 He/H of ISM, by number
n_X_CO = 10 an option to select functional dependence of X_CO=X_CO(R) 0=constant as below, 9=standard variation as in A&A 2004 paper 10=an exponential
X_CO = 1.9E20 conversion factor from CO integrated temperature to H2 column density
propagation_X_CO = 2 option to control H2 from CO for propagation
nHI_model = 1 selection of HI gas density model (not yet implemented)
nH2_model = 1 selection of H2 gas density model (not yet implemented)
nHII_model = 1 selection of HII gas density model 1=Cordes et al 1991 2, 3 = other models
COR_filename = rbands_co10mm_v3_2001_hdeg.fits.gz rbands_co9_qdeg.fits.gz
HIR_filename = rbands_hi12_v5_hdeg_zmax1_Ts125.fits.gz rbands_hi11_qdeg.fits.gz rbands_hi10_garbage.fits
B_field_name = galprop_original han_ran 3D models:han .... galprop_original: exponential model as in original (parameters Bo, rscale, zscale)
n_B_field_parameters = 10 number of B-field parameters (must be 10 at present)
B_field_parameters = 7.0e-6,50.0,2.00,12.26,5.00e-6,0.0,0.0,0.0,0.0,0.2 parameters for 3D models, meaning depends on the model
HII_Te = 7000 free electron temperature for free-free absorption
HII_clumping_factor = 60 free electron clumping factor for free-free absorption
#ISRF_file = ISRF_RMax20_ZMax5_DR0.5_DZ0.1_MW_BB_24092007.fits input ISRF file
ISRF_file = ISRF/Standard/Standard.dat input ISRF file THIS IS THE ONE FOR ANISOTROPIC IC
ISRF_filetype = 3 0=CMB only,1=obsolete, 2=standard FITS, 3=healpix with angular dependence
ISRF_healpixOrder = 3 for output of ISRF skymaps
ISRF_factors = 1.0,1.0,1.0 ISRF factors for IC calculation: optical, FIR, CMB
fragmentation = 1 1=include fragmentation
momentum_losses = 1 1=include momentum losses
radioactive_decay = 1 1=include radioactive decay
K_capture = 1 1=include K-capture
ionization_rate = 1 1=compute ionization rate
use_symmetry = 0 0=no symmetry, 1=optimized symmetry, 2=xyz symmetry by copying(3D)
vectorized = 0 0=unvectorized code, 1=vectorized code
SNR_events = 0 handle stochastic SNR events
SNR_interval = 1.0e4 time interval in years between SNR in 1 kpc^-3 volume
SNR_livetime = 1.0e4 CR-producing live-time in years of an SNR
SNR_electron_sdg = 0.00 delta electron source index Gaussian sigma
SNR_nuc_sdg = 0.00 delta nucleus source index Gaussian sigma
SNR_electron_dgpivot = 5.0e3 delta electron source index pivot rigidity (MeV)
SNR_nuc_dgpivot = 5.0e3 delta nucleus source index pivot rigidity (MeV)
proton_norm_Ekin = 1.00e+5 proton kinetic energy for normalization (MeV)
proton_norm_flux = 5.75e-9 6.75e-9 5.75e-9 5.00e-9 to renorm nuclei/flux of protons at norm energy (cm^-2 sr^-1 s^-1 MeV^-1)
electron_norm_Ekin = 3.45e4 electron kinetic energy for normalization (MeV)
electron_norm_flux = 0.32e-9 0.56e-9 0.32e-9 0.4e-9 0.6e-9 0.40e-9 flux of electrons at normalization energy (cm^-2 sr^-1 s^-1 MeV^-1)
total_cross_section = 2 total cross section option: 0=L83 1=WA96 2=BP01
cross_section_option = 012 100*i+j i=1: use Heinbach-Simon C,O->B j=kopt j=11=Webber, 21=ST
t_half_limit = 1.0e4 year - lower limit on radioactive half-life for explicit inclusion
primary_electrons = 0 1=compute primary electrons
secondary_positrons = 0 1=compute secondary positrons
secondary_electrons = 0 1=compute secondary electrons
knock_on_electrons = 0 1,2 1=compute knock-on electrons (p,He) 2= use factor 1.75 to scale pp,pHe
pairproduction = 0 1=compute pair production on ISRF
secondary_antiprotons= 1 1,2= calculate: 1=uses nuclear scaling; 2=uses nuclear factors (Simon et al 1998) NB parameter name changed at r548
tertiary_antiprotons = 1 1=compute tertiary antiprotons NB parameter name changed at r548
secondary_protons = 0 1=compute secondary protons
gamma_rays = 2 2 1=compute gamma rays, 2=compute HI,H2 skymaps separately
pi0_decay = 3 1 1= old formalism 2=Blattnig et al. 3=Kamae et al.
IC_isotropic = 0 1,2= compute isotropic IC: 1=compute full, 2=store skymap components
IC_anisotropic = 0 1,2,3= compute anisotropic IC: 1=full, 2=approx., 3=isotropic
synchrotron = 0 2=compute synchrotron using B_field_model and B_field_parameters 3=all Stokes (3D)
free_free_absorption = 0 >=1 free-free absorption of synchrotron
bremss = 0 1=compute bremsstrahlung
globalLuminosities = 0 1 1=compute global luminosities, 0=don't (default)
skymap_format = 3 fitsfile format: 0=old format (the default), 1=mapcube for glast science tools, 2=both, 3=healpix
output_gcr_full = 1 output full galactic cosmic ray array
warm_start = 0 read in nuclei file and continue run
verbose = 0 -456 -455 -454 -453 verbosity: 0=min,10=max <0: selected debugs
test_suite = 0 run test suite instead of normal run
Looking through the manual I've finally understood how GALPROP works, but still have lots of questions.
1) What are the most important GALPROP parameters one should concern when investigating the final antiproton spectrum from DM annihilations/decay? The ones that I in general understand and can control are in the "Dark Matter properties" section. I take care of them in the gen_DM_source.cc file.
Which ones could I safely drop?
2) I'm new to FITS format. Astropy library for Python prints out the following:
Code:
Filename: nuclei_full_54_test
No. Name Type Cards Dimensions Format
0 PRIMARY PrimaryHDU 25 (23, 23, 23, 19, 1) float32
25
SIMPLE=True
BITPIX=-32
NAXIS=5
NAXIS1=23
NAXIS2=23
NAXIS3=23
NAXIS4=19
NAXIS5=1
EXTEND=True
COMMENT= FITS (Flexible Image Transport System) format is defined in 'Astronomy
and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
COMMENT= FITS (Flexible Image Transport System) format is defined in 'Astronomy
and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
EXPOSURE=1500
CRVAL1=-11.0
CRVAL2=-11.0
CRVAL3=-11.0
CRVAL4=1.0
CRVAL5=1.0
CDELT1=1.0
CDELT2=1.0
CDELT3=1.0
CDELT4=0.342422680822206
CDELT5=1.0
NUCZ001=-1
NUCA001=1
NUCK001=0
But working with data it turns out the the last axis (#5) is actually the first and vice versa. Is that the way it should be?
3) I'm using skymap_format = 3. It was not said explicitely in the tutorial, but on the page 13 there is a hint on the units: cm-2 sr-1 s-1 MeV-1. Please tell me if I'm right or wrong on the following.
This code will give the flux in the units mentioned above of the only particle species (pbar) that I'm concerning at the center of the galaxy:
Code:
from astropy.io import fits
hdulist = fits.open('nuclei_full_54_test')
val = hdulist[0].data[0][i][10][10][10]
e = math.pow(10, hdulist[0].header["CRVAL4"] + i * hdulist[0].header["CDELT4"])
Where hdulist[0] is the first (and only) table with data, i is the index of the energy bin, val - the desired flux, e - energy value for i'th bin. There are total 23 bins in each dimension, so 10 corresponds to the center of the axis.
Am I correct in my calculations and units?
4) Let's say I set network_iterations to be equal to 10.
After running the code I have this output:
Code:
...
propagate_particles: Entry
nuclei_normalize: Entry
nuclei_normalize: CR protons not found!
nuclei_normalize: Entry
nuclei_normalize: CR protons not found!
nuclei_normalize: Entry
nuclei_normalize: CR protons not found!
nuclei_normalize: Entry
nuclei_normalize: CR protons not found!
nuclei_normalize: Entry
nuclei_normalize: CR protons not found!
nuclei_normalize: Entry
nuclei_normalize: CR protons not found!
nuclei_normalize: Entry
nuclei_normalize: CR protons not found!
nuclei_normalize: Entry
nuclei_normalize: CR protons not found!
create_transport_arrays: Entry
create_transport_arrays: Assigning primary source function
create_transport_arrays: DM_antiprotons g_0=0 rigid_br0= 0 g_1=1.98 rigid_br= 9000 g_2=2.42
create_transport_arrays: >>>>>>>>>>>>>>>>>>> norm 1 >>>>>>>>>>>>>>>>>>>
create_transport_arrays: assigning diffusion coefficient
create_transport_arrays: beta at break rigidity=0.973589 rigidity_br=4000?= galdef.D_rigid_br=4000
create_transport_arrays: Calculating diffusion coefficients (ix, iy, iz, ip): (23, 23, 23, 19)
create_transport_arrays: ======== assigning fragmentation rate ========
create_transport_arrays: ======== assigning momentum loss rate ========
create_transport_arrays: ============== completed creation of transport arrays for DM_antiprotons
create_transport_arrays: Exit
propagate_particles:
Network iteration 9 species 0 DM_antiprotons (Z,A) = (-1,1)
propel: Entry
propel_diagnostics: Entry
propel_diagnostics: initializing variables for new species
propel: Zero primary and secondary source functions. No propagation will be done
propel: Exit
nuclei_normalize: Entry
nuclei_normalize: CR protons not found!
create_transport_arrays: Entry
create_transport_arrays: Assigning primary source function
create_transport_arrays: DM_antiprotons g_0=0 rigid_br0= 0 g_1=1.98 rigid_br= 9000 g_2=2.42
create_transport_arrays: >>>>>>>>>>>>>>>>>>> norm 1 >>>>>>>>>>>>>>>>>>>
create_transport_arrays: assigning diffusion coefficient
create_transport_arrays: beta at break rigidity=0.973589 rigidity_br=4000?= galdef.D_rigid_br=4000
create_transport_arrays: Calculating diffusion coefficients (ix, iy, iz, ip): (23, 23, 23, 19)
create_transport_arrays: ======== assigning fragmentation rate ========
create_transport_arrays: ======== assigning momentum loss rate ========
create_transport_arrays: ============== completed creation of transport arrays for DM_antiprotons
create_transport_arrays: Exit
gen_secondary_source: Entry
gen_DM_source (v.2)
generating DM_antiprotons source function for n_spatial_dimensions=3
DMcs_v=4.73251e-28
<<<< gen_DM_source
gen_secondary_source: Exit
propagate_particles:
Network iteration 10 species 0 DM_antiprotons (Z,A) = (-1,1)
propel: Entry
propel_diagnostics: Entry
propel_diagnostics: initializing variables for new species
propel: Generating alpha for 3D
propel: timestep_mode=1
propel: timestep_mode=2
propel: Exit
nuclei_normalize: Entry
nuclei_normalize: CR protons not found!
propagate_particles: Cleaning up for: 1 number of species.
propagate_particles: Exit
propagate_particles: Time elapsed: 1m: 52.4s (Real); 1m: 52s (Process); 0.388s (System)
cr_luminosity: Entry
cr_luminosity: computing cr_luminosity for n_spatial_dimensions=3
cr_luminosity: CR protons not found!
store_gcr: Entry
store_gcr: iz value for output = 11, z[iz1] = 0
store_gcr: NUCZ001 -1
store_gcr: NUCA001 1
store_gcr: NUCK001 0
store_gcr: Exit
store_gcr: Time elapsed: 0.000358s (Real); 0s (Process); 0s (System)
store_gcr_full: Entry
store_gcr_full: Storing full nuclei array in !../out/nuclei_full_54_test.gz
store_gcr_full: NUCZ001 -1
store_gcr_full: NUCA001 1
store_gcr_full: NUCK001 0
store_gcr_full: Exit
store_gcr_full: Time elapsed: 0.0084s (Real); 0.008s (Process); 0s (System)
>>>> read_gas_maps HIR
reading HIR from /soft/science/_galprop/fitsdata/rbands_hi12_v5_hdeg_zmax1_Ts125.fits.gz
NAXIS = 3
NAXIS1,2,3 = 720 360 9
CRVAL1,2 = 0.25 -89.75
CDELT1,2 = 0.5 0.5
generating galaxy.HIR:
All done
<<<< read_gas_maps HIR
read_gas_maps: Time elapsed: 0.306s (Real); 0.284s (Process); 0.024s (System)
>>>> read_gas_maps COR
reading COR from /soft/science/_galprop/fitsdata/rbands_co10mm_v3_2001_hdeg.fits.gz
NAXIS = 3
NAXIS1,2,3 = 720 360 9
CRVAL1,2 = 0.25 -89.75
CDELT1,2 = 0.5 0.5
generating galaxy.COR:
All done
<<<< read_gas_maps COR
read_gas_maps: Time elapsed: 0.287s (Real); 0.276s (Process); 0.008s (System)
9
gen_pi0_decay_emiss: Entry
gen_pi0_decay_emiss: Generating pi0-decay emissivity for n_spatial_dimensions=3
CR protons not found!
gen_pi0_decay_emiss: Time elapsed: 0.00073s (Real); 0s (Process); 0s (System)
main: Exit
It looks like the first 9 iterations were done before the antiproton sources were specified (before the line gen_DM_source (v.2) was printed). And the 10'th iteration was done after that.
Also I'm wondering what all those "CR protons not found!" mean? What do I need protons for? I though everything excpet antiprotons have been desactivated.
FAQ
Search
Register
Login
