Questions about the cosmic ray simulation

FLUKA: FLUKA 4-5.0
Flair: FLAIR 3.4-2
Operating system: Ubuntu 22.04.4 LTS
Compiler: gfortran 11.4.0

Description

Dear FLUKA experts
I am simulating the propagation of cosmic rays at a specific geographic location. However, I have encountered some issues here that are difficult to resolve.
1. I have set up a dipole field in the GCR-SPE section of my try.inp file, but how can I plot this magnetic field in Flair? Could you please tell me the specific steps? (This is because I usually run the simulation directly in the terminal using these commands:
$FLUPRO/bin/fff fluscw.f
$FLUPRO/bin/ldpmqmd -o myfluka -m fluka fluscw.o
$FLUPRO/bin/rfluka -e myfluka -N0 -M1 try.inp)
2. I have recorded particle information from the 100th atmospheric layer down to the ground using fluscw.f, but how do I normalize the data obtained from fluscw.f to units of counts/(m²·s·sr)?
3.In the simulation, a geomagnetic Earth model is used. So, are the direction and position information of the obtained cosmic rays also based on this geomagnetic Earth model? Additionally, I would like to know: since the geomagnetic north and south poles are opposite to the geographic poles, why is the geomagnetic position generated based on a geographic location (e.g., a certain north latitude) also at north latitude?

Dongguan_cern.zip (8.2 MB)

Input files

I have uploaded all the content and hope you can help me resolve the issues. I would be grateful.

Best regards
Yongce Gong,

Dear @chenjianqi,

  1. Note that the magnetic field specified in the GCR-SPEcard is a “special” magnetic field, which is not the same as a typical magnetic field specified with the MGNFIELDcard, and therefore cannot be plotted as the MGNFIELD. If I nonetheless find a way to plot it, I will get back to you. Generally, to plot a magnetic field defined via the MGNFIELDcard in FLAIR, you can:

    1. go to the “Plot” tab
    2. “Add” a “Geometry” plot.
    3. In the bottom right of the header, at “Type”, click and select from the drop-down menu any of the three “Field” options.
    4. Specify the coordinates for the center (e.g. z=2.2e8, x=6.02e8) and the extension (du=dv=1e6) of your plots, where your magnetic field should be located.
  2. I see in your fluscw.f routine, you are scoring each particle that arrives in your region of interest, for a total of M particles, which obviously are not normalized to anything. When you launch your simulations, you should save the number of primaries N you are using, and you would have the data as counts/primary.
    If you want to moreover normalize to counts/(m^2 s sr), please note that the <ZA>phi0465.spc files contain a distribution, where the flux is already given in particles/(m^2 s sr)/MeV. You should then compute the integrated flux from all these files, as Phi_{tot}, which will be in units ofparticles/(m^2 s sr)
    Finally, the normalization of your results should be: M * Phi_{tot} / N

  3. Please refer to the FLUKA manual, and to add more info:

    1. In FLUKA, cosmic-ray positions are generated at geographic coordinates, but their propagation is bent and filtered by the geomagnetic field model at that location. In othe words, geographic coordinates are the user-facing input, but geomagnetic coordinates are internally mapped by the model.
    2. “Why” is this so it is more a question of convention/choice, as most space weather codes use geographic coordinates as the reference grid, and then apply the geomagnetic field model on top of them.

I hope this helps,
Daniel

Dear Daniel Prelipcean
Thank you for your reply. I look forward to hearing from you about the method for drawing magnetic fields in GCR-SPE.

Additionally, I would like to ask: in fluscw.f, I recorded the weight, and the weight of each particle is not 1. So, during normalization, should it be: weight * M * Phi_{tot} / N?

I also have another question: In the SPECSOUR card, WHAT(7) can define the method for geomagnetic cutoff (WHAT(7) = 0: no geomagnetic cutoff; = 1: geomagnetic cutoff is requested (centered dipole); = 2: the vertical geomagnetic cutoff is read from WHAT(2); = 3: geomagnetic cutoff is requested (offset dipole, possibly + IGRF11)). If I am simulating the interaction of global cosmic rays with the atmosphere, why is the distribution of secondary cosmic rays reaching the ground isotropic after using the centered dipole? (As shown in the image) (Shouldn’t there be more at the poles in theory?)

Thank you for your reply,
Yongce Gong

Dear @chenjianqi,

  1. Indeed, as you correctly observed, you should weight each primary with its respective weight, i.e. \sum w_i, which in your case will be less than M. If all particles were to have w_i =1, then the outcome would have been M.

  2. Regarding the magnetic field, when FLAIR plots the magnetic field, some temporary files are produced and then picked up/used, but currently not for the GCR-SPEC card. In the next release of FLAIR, this will be included and you will be able to obtain the magnetic field plot, e.g. as below:

  3. In principle, you are right to think that the an anisotropic magnetic field should also produce anisotropies at ground level, however you should also quantify the expected anisotropy. Not also that in your plot, your colorbar scale spans more than 6 orders of magnitute - maybe your “poles” are just not visible in this view.

  4. Could you (re)iterate which FLUKA version are you using? The IGRF11 is currently not shown in the FLUKA manual.

Best,
Daniel

Dear Daniel
0. So the correct normalization method is: I = (∑w_i / N) * Φ_tot ? And is Φ_tot the total flux of 28 types of ions? Since the SPECSOUR card specifies 28 types of cosmic ray particles as the source, will it automatically sample ions of different energies uniformly from the energy differential spectra of these 28 particles as incident ions? How can I check the actual information of my incident cosmic ray particles?

1.My version of flair is FLAIR 3.4-2, and I noticed that the CERN version of flair has released FLAIR 3.4-3. If I upgrade my flair version to FLAIR 3.4-3, will I be able to plot the magnetic field of GCR-SPEC? Could you please tell me exactly how to do it?

2.By dividing the neutron flux according to latitude, it was confirmed that the neutron flux indeed varies with latitude, and the previous graph did have issues.

3.My FLUKA version is FLUKA 4-5.0. For the geomagnetic cutoff, it can be set to IGRF11. I learned this from the following link (https://indico.cern.ch/event/334606/contributions/779790/attachments/653365/898411/AdvancedCosmicRays2014.pdf). As you can see in the screenshot, simply setting what(7) in the SPECSOUR card to 3 enables the IGRF11 geomagnetic cutoff screening. However, after using this setting, why is particle data only available in the Eastern Hemisphere? Please refer to the first image below. (But if what(7) is set to 0 or 1, particle data is available globally, as shown in the second image below.)
Actually, I want to know if I can implement my own geomagnetic cutoff screening. Are there any user routines that would allow me to do this?

4.Additionally, I have another question: why is my fluscw.f not recording neutrino information?

Thank you for your attention to these issues.

Best wishes,
Yongce

Dear @chenjianqi,

  1. Yes. If you want to check the sampled cosmic ray particles, you can adapt a userroutine (e.g. mgdraw) and dump the information of the source particles by filtering on BEAMPART particle id. Make sure to only count them once.
  2. No, this update will be part of the next release (not the current 3.4-3). In my first answer in this thread I explained the step-by-step procedure.
  3. Note that the data is not only in the Eastern Hemisphere. Rather, per the slides you reference and the input file, you only have a slice of the atmosphere, hence the cuts that you observe.
  4. For your own geomagnetic cutoff screening, there are no user routines readily available to do this, but you should write your own. What type of cutoffs (different dependency on location than those currently available?) would you like to use?
  5. The fluxcw.f is not recording neutrino information because FLUKA does not transport neutrinos through matter by default, however neutrino production is included in FLUKA.

Hope this helps.
Daniel

Dear Daniel
Thank you for addressing these issues; most of them I have been learnd.

However, regarding the “Eastern Hemisphere” problem, you may not have fully understood my point. The situation is as follows: I have already switched the atmospheric model to a global atmospheric model instead of just a localized one. When I set the what(7) card in SPECSOUR to 0 or 1, secondary cosmic ray particles are distributed globally. However, when I set the what(7) card in SPECSOUR to 3, secondary cosmic ray particles are only distributed in the Eastern Hemisphere. The only change I made was to the what(7) card in SPECSOUR.
Best wishes,
Yongce

Dear @chenjianqi,

Just to clarify: I assume that the two plots that you are showing are those at ground level, obtained from your fluscw routine?

Best,
Daniel

Dear Daniel,
Yes, you are right,
Best,
Yongce

Dear @chenjianqi,

I would like point out that for FLUKA 4-5.0, the IGRF11 option that you mention is not supported (as you can see in the manual). The presentation that you linked is for an experimental version of FLUKA, that is not publicly released.

That being said, could you please also upload the scripts you used for plotting your results, as in the screenshots you attached? I cannot find anything wrong so far - and in your results seem to indicate a hard border where particles accumulate - which is rather strange.

Best,
Daniel

Dear Daniel Prelipcean
Now I have decided not to use IGRF11 anymore.
And for the script, you can check in the GitHub - yongce3/3d_vis: The distribution of secondary particles from cosmic rays on Earth.

Best
Yongce,