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 2 of 2
Author:  StanislavIablokov [ Mon Jun 05, 2017 4:43 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
imos wrote:
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.


Let's say we have an AMS-02 or PAMELA data on the antiproton flux. It could be subdivided into two portions: known sources of antiprotons and unknown (potential DM decays/annihilations). The unknown-type flux could be simulated by injecting my antiproton source function into gen_DM_source.cc. What about flux due to known sources? Could GALPROP simulate that flux?
Author:  StanislavIablokov [ Mon Jun 05, 2017 4:46 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.

Actually they are from this paper dated 2013: https://arxiv.org/abs/1301.7079
Author:  imos [ Mon Jun 05, 2017 4:51 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
You should really look into GALPROP papers.
Here is one of the most recent:
http://adsabs.harvard.edu/abs/2017ApJ...840..115B

where AMS-02 pbar measurements are nicely reproduced without any fitting to pbars and without DM contribution.
It takes into account only secondary and tertiary pbars.

The paper that you look at and similar are obsolete because of many different reasons.
Author:  gudlaugu [ Mon Jun 05, 2017 2:39 pm ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
StanislavIablokov wrote:
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?


NUCNORM should not be used to calculate anything in the output file. It is there for information purposes only, it is the number used internally by galprop to normalize the source spectrum so the propagated proton flux agrees with what is specified in the GALDEF file.

StanislavIablokov wrote:
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.


This does not differ from my answer. GALPROP internally uses n*c/4pi and the source term in gen_DM_source must also have this factor for things to work. That is why it is there in the original gen_DM_source routine.

StanislavIablokov wrote:
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?


This is not the original gen_DM_source.cc routine distributed with the galprop code. I am afraid you are on your own regarding the units. Just remember that gen_DM_source.cc should return c/4pi * dn/dp /dt. And we use MeV-1 for momentum and energy units.
Author:  StanislavIablokov [ Mon Jun 05, 2017 3:53 pm ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
Quote:
NUCNORM should not be used to calculate anything in the output file.

So the final fluxes for ALL listed types of particles are already given with their absolute values, right?
Author:  gudlaugu [ Tue Jun 06, 2017 1:16 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
StanislavIablokov wrote:
So the final fluxes for ALL listed types of particles are already given with their absolute values, right?


Yes
Author:  StanislavIablokov [ Thu Jun 08, 2017 4:56 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
In the gen_DM_source.cc file there is a distribution function for annihilation products:
Quote:
DMbranching*exp(-pow((DMmass-particle.Etot[ip]*1.e-3)/DMwidth,2))/DMmass*1.e-3;

Where was it taken from? As far as I could tell it is the fragmentation function of the primary DM annihilation products (such as quarks) into stable particles (such as antiprotons).

The one that I've got with Pythia looks very different than the one in your code with the default parameters from galdef-file (DMwidth=50 GeV, DMbranching=40) and the DMmass=300GeV.

So Pythia gives the fragmentation function, for example, for u + ubar -> (hadronization) -> pbar + * very different in shape and in orders of magnitude (blue curve) than the default one in gen_DM_source.cc (black curve).

By fragmentation function I mean the number of particles (in my case antiprotons) created from the jet of quarks (here, u ubar) in the unit of energy interval.

Attachments:
fig.png
fig.png [ 46.8 KiB | Viewed 114670 times ]
Author:  StanislavIablokov [ Thu Jun 08, 2017 9:03 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
Now I'm pretty sure that there is something going on with GALPROP that I'm not aware about.

I'm running two simulations with the same source function: 2D and 3D. Here is the code of gen_DM_source.cc:
Code:
// assign the source function (2D)

   if(galaxy.n_spatial_dimensions==2)
     {
       for(int ir=0; ir<gcr[0].n_rgrid; ir++)
    {
      for(int iz=0; iz<gcr[0].n_zgrid; iz++)
        {
          for(int ip=0; ip<particle.n_pgrid; ip++)
       {

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

                       if (strcmp(particle.name,"DM_antiprotons") == 0 && galdef.DM_int1 == 777)
                       {
                        if (particle.Etot[ip] * 1.e-3 <= DMmass)
                        {
                         double source = 0.5 * 3 * pow(10, -26) * pow(800 / DMmass, 2) * pow(DM_profile_av(galaxy.r[ir], galaxy.z[iz], galaxy.dr, galaxy.dz, dzz) / DMmass, 2) *
                                         pbarqqFF(particle.Etot[ip] * 0.001 - 0.938, t_a, aSize, t_b, t_c, t_d) / 1000;

                         particle.secondary_source_function.d2[ir][iz].s[ip] += (SoL/(4 * Pi)) * source;
                        }
                        continue;
                       }

             /////////////////////////////////////////////////////////////////////
       } // ip
        }  //  iz
    }  //  ir
     }  //  particle.n_spatial_dimensions==2
   
// assign the source function (3D)

   if(galaxy.n_spatial_dimensions==3)
     {
       for(int ix=0; ix<gcr[0].n_xgrid; ix++)
    {
      for(int iy=0; iy<gcr[0].n_ygrid; iy++)
        {
          for(int iz=0; iz<gcr[0].n_zgrid; iz++)
       {
         for(int ip=0; ip<particle.n_pgrid; ip++)
           {
             /////////////////////////////////////////////////////////////////////

                       if (strcmp(particle.name,"DM_antiprotons") == 0 && galdef.DM_int1 == 777)
                       {
                        if (particle.Etot[ip] * 1.e-3 <= DMmass)
                        {
                         double source = 0.5 * 3 * pow(10, -26) * pow(800 / DMmass, 2) * pow(DM_profile_av(galaxy.x[ix], galaxy.y[ix], galaxy.z[iz], galaxy.dx, galaxy.dy, galaxy.dz, dzz) / DMmass, 2) *
                                         pbarqqFF(particle.Etot[ip] * 0.001 - 0.938, t_a, aSize, t_b, t_c, t_d) / 1000;

                         particle.secondary_source_function.d3[ix][iy][iz].s[ip] = (SoL/(4 * Pi)) * source;

                        }
                        continue;
                       }

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

           } //ip
       }  //  iz
        }  //  iy
    }  //  ix
     }  //  particle.n_spatial_dimensions==3


As you can see particle.secondary_source_function is the same for both 2D and 3D. Function pbarqqFF is my fragmentaion function of quarks into antiprotons. It is also the same in both cases.

But while the 2D results are acceptable, the 3D ones look crazy. I've deliberately drawn the proton's flux which is the same in both cases (as expected). I think something is going on with 3D source functions, something I'm not accounting for.

I've also attached a galdef file for 2D. The spatial grids for 3D are less detailed because it takes more time to compute in 3 dimensions than in 2, but that wouldn't change the results by 20 orders of magnitude.
The only difference between 2D and 3D is the following:
Code:
n_spatial_dimensions = 3     
z_min                = -10.0     min z
z_max                = +10.0     max z
dz                   = 1.0       delta z
x_min                = -10.0     min x
x_max                = +10.0     max x
dx                   = 1.0       delta x
y_min                = -10.0     min y
y_max                = +10.0     max y
dy                   = 1.0       delta y

prop_x               = 1     1=propagate in x (3D)
prop_y               = 1     1=propagate in y (3D)
prop_z               = 1     1=propagate in z (2D, 3D)
prop_p               = 1     1=propagate in p (2D, 3D)

p_Ekin_grid          = Ekin     p||Ekin alignment
p_min                = 100     min momentum (MV)
p_max                = 40000     max momentum (MV)
p_factor             = 1.5     momentum factor
Ekin_min             = 1.0e1     min kinetic energy per nucleon (MeV)
Ekin_max             = 1.0e6     max kinetic energy per nucleon (MeV)
Ekin_factor          = 1.5     kinetic energy per nucleon factor


Am I missing something?

Attachments:
File comment: 2D galdef file
galdef_54_test2.txt [17.23 KiB]
Downloaded 7621 times
File comment: 3D plot with crazy pbar flux
fig3.png
fig3.png [ 36.16 KiB | Viewed 114670 times ]
File comment: 2D plot with normal pbar flux
fig2.png
fig2.png [ 36.91 KiB | Viewed 114670 times ]
Author:  StanislavIablokov [ Thu Jun 08, 2017 10:22 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
It also turns out that when making 2D simulations with different grid steps the results could also differ by 20 orders of magnitude. The rest stays the same, I modify only the spatial grid step.

With 0.5-1.0 kpc step the results are consistent with what's expected, but switching to 0.1-0.2 kpc they go crazy.
As usual proton flux is the same in both cases.

Attachments:
File comment: Crazy results with 0.2 kpc step
fig2-0.2.png
fig2-0.2.png [ 33.51 KiB | Viewed 114670 times ]
File comment: Good results with 0.5 kpc step
fig2-0.5.png
fig2-0.5.png [ 38.72 KiB | Viewed 114670 times ]
Author:  gudlaugu [ Thu Jun 08, 2017 10:44 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
This is really strange. There is nothing in GALPROP that should make this much of a difference between 2D and 3D. And the grid size should not have this large of an effect. If I read the plots correctly, the shape is nearly always correct while the normalization is off by large factors. Is it possible that you are using an uninitialized variable somewhere in the code? That could result in weird random effects.
Author:  StanislavIablokov [ Thu Jun 08, 2017 10:53 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
I can't see any variable that I'm not controlling.

I think the shape also changes a bit from 2D to 3D and in 2D when I change step size.

What would you normally suggest to do? Are there any debugging tools inside the GALPROP?
Author:  gudlaugu [ Thu Jun 08, 2017 11:04 am ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
There are various degrees of debugging output that can be tuned by adjusting the verbose parameter. It is however not documented and some options produce a huge amount of information. If you grep the code for verbose you will have some indications about what number to choose to get the information you need.

Because you are obviously compiling the code yourself, one option is to add a lot of std::cout<< statements to print various variables. Adding --enable-debug may also help. I would start by making sure the source function at some specific point is not dependent on the grid size or 2D vs. 3D. Once the source function is found to be OK, then you can see if the distribution after propagation depends on these parameters. Slowly work you way from the source to the output so we can narrow down which code path is the culprit.

Also note that you are not required to have the same grid step in all dimensions and you usually need the grid to be a lot finer in the z direction to resolve the structure of the Galaxy. You are also using a very large z_max that does affect your selection for the diffusion parameter if you want to have a reasonable set of parameters that explain the B/C ratio.
Author:  emmalee [ Fri Jul 28, 2017 3:47 pm ]
Post subject:  Re: DM galdef setup with gen_DM_source.cc tuning
Recently new data with large statistics on both low and high energy antiprotonfluxes have become available which allow such tests to be performed.
Accurate calculation of the secondary antiproton flux provides a ``background'' for searches for exotic signals from the annihilation of supersymmetric particles and primordial black hole evaporation
Page 2 of 2 All times are UTC - 8 hours [ DST ]
Powered by phpBB® Forum Software © phpBB Group
https://www.phpbb.com/