GALPROP
https://galprop.stanford.edu/forum/

DM galdef setup with gen_DM_source.cc tuning
https://galprop.stanford.edu/forum/viewtopic.php?f=13&t=237		 Page 1 of 2
Author:  StanislavIablokov [ Thu Jun 01, 2017 11:02 am ]
Post subject:  DM galdef setup with gen_DM_source.cc tuning
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.
Author:  StanislavIablokov [ Thu Jun 01, 2017 1:17 pm ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
I also have a question about gen_DM_source.cc routine.
There is a parameter that is responsible for pbar source function: particle.secondary_source_function.d3[ix][iy][iz].s[ip]

Could someone please explain:
1) What are the units and the sense this parameter has?
2) Why it is called secondary? What is primary then?
Author:  gudlaugu [ Fri Jun 02, 2017 3:05 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
Thank you for your interests in the GALPROP project.

Looking at your configuration I see you have a box for your galaxy and a bin width of 1 kpc in z. This is way too large because a lot of the important physics included in GALPROP has to do with the interstellar gas and the ISM (e.g. energy losses). I suggest you reduce this at least to .2 kpc so you resolve the structure of the disk. You also exclude all nuclei so I think you are only interested in testing propagation of DM. GALPROP was not really designed for that and much of the code expects to have protons included. I don't think it will affect the DM propagation but I am not sure.

I'll try to answer your questions, starting with the first post.

1. For the general propagation parameters I must refer you to the multitude of publications. I suggest starting with the review by Strong et al. (2007) and going from there to more specialized publications like our recent results using Bayesian analysis in GALPROP (Johannesson et al. 2016). There are also many more publications listed on this page. I will have to refer the questions regarding the gen_DM_source routine to someone else as I am not familiar with that. You should be able to figure it out from the code and comments, it is very much self contained in that single file.

2. Yes, that is the way astropy and formerly pyfits works. The last axis is the first index in the resulting array.

3.Skymap format refers to the gamma ray output and has no effect on the nuclei file. The units of the nuclei file are MeV/(cm2 sr s). The code you posted should give you the spectrum at the GC.

4. The network iterations are only useful for the generation of secondary nuclei (B in particular) and wave damping. The code is set up so that it does at most 2 iterations of the nuclei and only propagates other particles, such as secondary electrons and anti-protons once. This is because doing it multiple times has no effect at all. So what happens is that for 9 iterations it tries to propagate the nuclei, but you didn't specify any so it is just an empty loop.
The error message for protons not found is printed in the nuclei_normalize routine. GALPROP uses relative normalization on the propagated spectrum and it is based on the CR flux of protons at a specific energy specified with
Code:
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)

I think it does not affect the secondary DM propagation, but galprop was not designed to be run without protons. In fact, you won't get any secondary antiprotons unless you specify at least protons, and then you have them wrong because helium contributes significantly to the antiprotons.

Answering your questions from the second post

1. The units are given in the comments of the gen_DM_source.cc file, cm^-3 s^-1 MeV^-1
2. It is secondary because the DM is the primary particle. It has no effect in the end, both are treated equally in propagation.
Author:  StanislavIablokov [ Fri Jun 02, 2017 3:33 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
Quote:
The units of the nuclei file are MeV/(cm2 sr s)

So in order to get the flux per unit energy I have to divide the result from the nuclei file by energy (in MeV) squared:
Code:
val = hdulist[0].data[0][i][10][10][10]
e = math.pow(10, hdulist[0].header["CRVAL4"] + i * hdulist[0].header["CDELT4"])
flux = val / e**2  # <- flux in MeV-1 s-1 cm-1 s-1, right?   

Am I correct?

Quote:
For the general propagation parameters I must refer you to the multitude of publications.

Could you provide some sort of the default parameters?

Quote:
I think it does not affect the secondary DM propagation, but galprop was not designed to be run without protons. In fact, you won't get any secondary antiprotons unless you specify at least protons, and then you have them wrong because helium contributes significantly to the antiprotons.


So what I have to do is to set all "use_Z_n" to 1, so that the program would include the effects of antiproton production from nuclei. Am I correct?

Quote:
In fact, you won't get any secondary antiprotons unless you specify at least protons

So what is the default setup for protons?

Quote:
1. The units are given in the comments of the gen_DM_source.cc file, cm^-3 s^-1 MeV^-1

I saw that, but those c/4Pi confused me a lot. In my code I use:
Code:
             /////////////////////////////////////////////////////////////////////

                       if (strcmp(particle.name,"DM_antiprotons") == 0 && galdef.DM_int1 == 777)
                       {
                        if (particle.Etot[ip] * 1.e-3 <= DMmass)
                        {
                         particle.secondary_source_function.d3[ix][iy][iz].s[ip] += DMcs_v * pow(DM_profile_av(galaxy.x[ix], galaxy.y[ix], galaxy.z[iz], galaxy.dx, galaxy.dy, galaxy.dz, dzz) / DMmass, 2) *
                                                                                    (1/(4 * Pi)) * a * pow(particle.Etot[ip] * 0.001, b) * exp(-c * pow(particle.Etot[ip] * 0.001, d));
                        }
                        continue;
                       }

             /////////////////////////////////////////////////////////////////////

where 777 is a switch to interpret my distribution and
Code:
a * pow(particle.Etot[ip] * 0.001, b) * exp(-c * pow(particle.Etot[ip] * 0.001, d))

is a number (i.e. dimensionless) of antiprotons per annihilation of 2 DM particles.

I still wonder if I need to use 4 Pi in the denominator or not. As far as I could tell division by 4 Pi means that one is interested in emission per unit of solid angle. Could you please clarify that?
Author:  gudlaugu [ Fri Jun 02, 2017 3:42 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
StanislavIablokov wrote:
Quote:
The units of the nuclei file are MeV/(cm2 sr s)

So in order to get the flux per unit energy I have to divide the result from the nuclei file by energy (in MeV) squared:
Code:
val = hdulist[0].data[0][i][10][10][10]
e = math.pow(10, hdulist[0].header["CRVAL4"] + i * hdulist[0].header["CDELT4"])
flux = val / e**2  # <- flux in MeV-1 s-1 cm-1 s-1, right?   

Am I correct?

Yes.

StanislavIablokov wrote:
Quote:
For the general propagation parameters I must refer you to the multitude of publications.

Could you provide some sort of the default parameters?

There are parameters listed in the publications. Their values depend a bit upon the other setup so you will have to read up a bit on the matter.

StanislavIablokov wrote:
Quote:
I think it does not affect the secondary DM propagation, but galprop was not designed to be run without protons. In fact, you won't get any secondary antiprotons unless you specify at least protons, and then you have them wrong because helium contributes significantly to the antiprotons.


So what I have to do is to set all "use_Z_n" to 1, so that the program would include the effects of antiproton production from nuclei. Am I correct?

You will have to set max_Z to at least 2 and use_Z_n to 1

StanislavIablokov wrote:
Quote:
In fact, you won't get any secondary antiprotons unless you specify at least protons

So what is the default setup for protons?

This is again setup dependent. I suggest you take one of the galdef files from previous publications and modify them slightly to include your DM source.

StanislavIablokov wrote:
Quote:
1. The units are given in the comments of the gen_DM_source.cc file, cm^-3 s^-1 MeV^-1

I saw that, but those c/4Pi confused me a lot. In my code I use:
Code:
             /////////////////////////////////////////////////////////////////////

                       if (strcmp(particle.name,"DM_antiprotons") == 0 && galdef.DM_int1 == 777)
                       {
                        if (particle.Etot[ip] * 1.e-3 <= DMmass)
                        {
                         particle.secondary_source_function.d3[ix][iy][iz].s[ip] += DMcs_v * pow(DM_profile_av(galaxy.x[ix], galaxy.y[ix], galaxy.z[iz], galaxy.dx, galaxy.dy, galaxy.dz, dzz) / DMmass, 2) *
                                                                                    (1/(4 * Pi)) * a * pow(particle.Etot[ip] * 0.001, b) * exp(-c * pow(particle.Etot[ip] * 0.001, d));
                        }
                        continue;
                       }

             /////////////////////////////////////////////////////////////////////

where 777 is a switch to interpret my distribution and
Code:
a * pow(particle.Etot[ip] * 0.001, b) * exp(-c * pow(particle.Etot[ip] * 0.001, d))

is a number (i.e. dimensionless) of antiprotons per annihilation of 2 DM particles.

I still wonder if I need to use 4 Pi in the denominator or not. As far as I could tell division by 4 Pi means that one is interested in emission per unit of solid angle. Could you please clarify that?


Please see section 2.2 in the documentation.
Author:  StanislavIablokov [ Fri Jun 02, 2017 7:06 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
Quote:
Please see section 2.2 in the documentation.


Well, I still don't get that.

Section 2.2 wrote:
However the basic CR density used has
units of density per total momentum p since this is natural for propagation. The actual units used internally
are (c/(4(Pi)) * n(p), where n(p) = dn/dp in units of cm−3 MeV−1.


As far as I could tell from dimensions, n(p) is the amount of particles (in my case antiprotons) in 1 cm3 in the dP interval of momentum. But there is no s-1 dimension so that's still not a flux/rate_of_production. Right now it is just density per unit of momentum. Multyplying it by speed of light will give us s-1 dimention along with a cm+1 one, so it's not what we are aiming for.

As far as I could tell from comments in gen_DM_source.cc it is required to specify the rate of production per unit volume per second per interval of momentum.
gen_DM_source.cc wrote:
// The routine gen_DM_source calculates the source functions of the products of the
// dark matter (DM) particle annihilation [cm^-3 s^-1 MeV^-1].


And I still can't tell if particle.secondary_source_function.d3[ix][iy][iz].s[ip] has [cm^-3 s^-1 MeV^-1] dimensions or it should be multiplied by the speed of light. And this 4*Pi factor also is not clear for me.

I'm sorry, but from the manual I can't understand how to exactly specify particle.secondary_source_function.d3[ix][iy][iz].s[ip]. I've seen class declarations in .h files (Distribution.h, Spectrum.h). There are no comments on how the things are calculated and used. So in order to understand that I have to know how the internal code works - that's too much to grasp for me, therefore I'm asking here.

So, like in your formulas, the ingredients I have are "DMcs_v" for cross section times velocity and "pow(DM_profile_av(...)/DMmass, 2)" for concentration squared. And also I use the amount of antiprotons created in the interval of energy dE per 2 DM particles annihilations, let's call it F=dN/dE with dimensions MeV-1.

If particle.secondary_source_function.d3[ix][iy][iz].s[ip] has dimensions [cm^-3 s^-1 MeV^-1] then this combination is the right guy:
Code:
particle.secondary_source_function.d3[ix][iy][iz].s[ip] = 0.5 * DMcs_v * pow(DM_profile_av(...)/DMmass, 2) * F(particle.Etot[ip])


I couldn't see where to put speed of light. The same for 4*Pi.
Author:  gudlaugu [ Sat Jun 03, 2017 2:56 pm ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
The s-1 comes because the source term in the propagation equation has the units of the time derivative of the density per momentum. If you consider equation 2 then you can see that the internal units are half-way between density and flux. So if you specify the source term in cm^-3 MeV^-1 s^-1 you should multiply it with c/4pi for things to work properly.
Author:  StanislavIablokov [ Sat Jun 03, 2017 3:29 pm ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
gudlaugu wrote:
The s-1 comes because the source term in the propagation equation has the units of the time derivative of the density per momentum. If you consider equation 2 then you can see that the internal units are half-way between density and flux. So if you specify the source term in cm^-3 MeV^-1 s^-1 you should multiply it with c/4pi for things to work properly.

With c taken as 3 * 10^10 cm/s?



Code:
NUCZ001=-1
NUCA001=1
NUCK001=0
NUCZ002=-1
NUCA002=1
NUCK002=0
NUCZ003=-1
NUCA003=1
NUCK003=0

In my output I have three populations of antiprotons (as far as I could tell, you call them primary, secondary and tertiary).
Could you please clarify on the meaning of those terms?

In gen_DM_source.cc we have particle.secondary_source_function by which we define the source of antiprotons. Should I consider antiprotons from that source as primary or secondary (as the name suggests)? What is tertiary then?

Are primary, secondary and tertiary fluxes included in each other in some way or are they treated separately? The question that I'm concerned about is the total flux of pbar in the Solar System after my pbar source has been accounted for.

In your indexing scheme (in the log from the FITS file above) the 0th population are tertiary, the first is secondary and the 2nd stands for primary. Am I correct?

And also the fluxes of antiprotons I get differ by 25 orders of magnitude. Is that considred to be normal?

The 0th flux has it's maximum value of about 10^15, while 1st and 2nd are around 10^-10. I think it has something to do with NUCNORM=1.11099307365e-21 in the final results, but page 12 of the manual doesn't explain it in great details so that I'm able to understand. Should I multiply antiproton fluxes by NUCNORM to get the real values?
Author:  imos [ Mon Jun 05, 2017 2:49 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
Just to add to the discussion.

Protons should always be calculated because the proton array is used for the energy/rigidity scales and normalization to the local data at the end. All other isotopes and secondary particles, if calculated, are used the normalization derived from protons.

The DM particles (including antiprotons) are disconnected from other particles and are not using the normalization derived from protons, but they need protons for the energy scale.

Primary antiprotons - are probably were added by Gulli recently, I have to check.
Secondary - are those produced in hadronic interactions (pp, pA, Ap, AA),
tertiary - inelastically scattered secondaries that are treated as a separate component.
The secondary and tertiary antiprotons can be safely switched off (put to 0).

The problems with the bar fluxes you have is likely due to switching off protons, so the code results become insane.

The speed of light (constant.h):
C = 2.99792458e10, // cm/s, =c speed of light

however, there is also
c =3.0e10;

which is likely not used.

For the units, n*c/4pi is just the internal value of number density which has the flux units. If you use c/4pi * dn/dp, it is the spectrum with the differential flux units. Regarding the source functions, they have standard units.
Author:  StanislavIablokov [ Mon Jun 05, 2017 3:22 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
Quote:
The DM particles (including antiprotons) are disconnected from other particles and are not using the normalization derived from protons, but they need protons for the energy scale.

So the DM flux, which I currently get by
Code:
from astropy.io import fits
hdulist = fits.open('nuclei_full_54_test2', ignore_missing_end=True)
nucnorm = hdulist[0].header["NUCNORM"]

val = hdulist[0].data[0][i][10][10][18] # cm-2 s-1 sr-1 MeV-1 MeV+2
e = math.pow(10, hdulist[0].header["CRVAL4"] + i * hdulist[0].header["CDELT4"]) # MeV
flux = nucnorm * 1000 * 10000 * val / e**2  # GeV-1 m-2 sr-1 s-1

should be calculated without nucnorm?

Code:
For the units, n*c/4pi is just the internal value of number density which has the flux units. If you use c/4pi * dn/dp, it is the spectrum with the differential flux units. Regarding the source functions, they have standard units.

And now it all became more complicated, because this answer differs from the previous one by gudlaugu.
So in order to check that I'm doing it right I have to ask one more time about the correctness of source function calculation.
Right now I do it this way:
Code:
if (strcmp(particle.name,"DM_antiprotons") == 0 && galdef.DM_int1 == 777)
{
 if (particle.Etot[ip] * 1.e-3 <= DMmass)
 {
  particle.secondary_source_function.d3[ix][iy][iz].s[ip] += 0.5 * DMcs_v * pow(DM_profile_av(galaxy.x[ix], galaxy.y[ix], galaxy.z[iz], galaxy.dx, galaxy.dy, galaxy.dz, dzz) / DMmass, 2) *
  (SoL/(4 * Pi)) * param_a * pow(particle.Etot[ip] * 0.001, param_b) * exp(-param_c * pow(particle.Etot[ip] * 0.001, param_d)) / 1000;
 }
 continue;
}


DMcs_v has units cm+3 s-1 (cross section times velocity)
pow(DM_profile_av(...) / DMmass, 2) has units cm-6 (it is concentration of DM particles squared)
param_a * pow(particle.Etot[ip] * 0.001, param_b) * exp(-param_c * pow(particle.Etot[ip] * 0.001, param_d)) is my fragmentation function, which gives the number of antiprotons per 2 DM particles annihilation (coefficient 0.5 also accounts for that) per unit of energy (in GeV-1, but it is divided by 1000 to give MeV-1)

So just these quantities give me the resulting dimension of [cm-3 s-1 MeV-1], which is good. What about (SoL/(4 * Pi)) ? Should it be included in the calculation? Is yes then the dimension of the source function will be different from [cm-3 s-1 MeV-1]. If no then it looks wierd that the default sources in gen_DM_source.cc had those c / (4 * Pi). Which one is correct then?

P.S. If I'll survive this, I should probably write "GALPROP for dummies" sort of manual based on the experience I've got :)
Author:  imos [ Mon Jun 05, 2017 3:48 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
These are all good questions. I did not remember having nucnorm and SoL in this routine.
These are likely recent modifications made by Gulli or by Troy.

I will look into this issue when I am back from my vacations at the end of June.
Regarding this modified version, Gulli should know who did modifications and why.

I can suggest taking the older version of the routine and use it as a reference.
Author:  StanislavIablokov [ Mon Jun 05, 2017 4:19 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
imos wrote:
These are all good questions. I did not remember having nucnorm and SoL in this routine.
These are likely recent modifications made by Gulli or by Troy.

I will look into this issue when I am back from my vacations at the end of June.
Regarding this modified version, Gulli should know who did modifications and why.

I can suggest taking the older version of the routine and use it as a reference.

Okay, I'll look into the older ones.
Well, I hope Gulli and Troy will explain me this stuff before the end of June. My project has to be submitted in 10 days :)


I also have a question about the parameters of the medium. Normally people define MIN, MED and MAX scenarios. I've found how to set up some of those, but not all. (see the attachment)

Attachments:
models.png
models.png [ 28.17 KiB | Viewed 118831 times ]
Author:  imos [ Mon Jun 05, 2017 4:22 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
The scenarios listed in the table are all extreme. You are probably using one of the papers of 2001 made by people who did not use GALPROP.
Author:  StanislavIablokov [ Mon Jun 05, 2017 4:30 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
imos wrote:
The scenarios listed in the table are all extreme. You are probably using one of the papers of 2001 made by people who did not use GALPROP.

Well I have another papers with different numbers. Right now it is not a matter of the numbers, but rather of how to set them up.


Quote:
Primary antiprotons - are probably were added by Gulli recently, I have to check.

Are you referring to the astrophysical antiproton background? If yes, that would also be nice to see in action.
Author:  imos [ Mon Jun 05, 2017 4:33 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
Astrophysical antiproton background are secondaries and tertiaries produced in the interstellar medium.
Primary antiprotons are likely introduced to study the effect of secondary production in SNR shocks.

Re: min, med, max - we normally do not use anything like that. Results of the Bayesian analysis published recently (Johannesson et al. 2016)
is the proper way of putting the error bars on models.
Page 1 of 2 All times are UTC - 8 hours [ DST ]
Powered by phpBB® Forum Software © phpBB Group
https://www.phpbb.com/