GCR extended layers and scoring through mgdraw file

Dear Fluka users,

I’m running and testing GCR simulations to calculate muon fluxes underground.
To do so, I’ve started using the example provided by fluka (AllParticleExample). I’ll divide this topic into three subtopics, I hope this is not a problem.

  1. Geometry: I used the geometry of 100 layers implementing SPECSOUR and GCR-SPE. I manually increase the number of layers to extended them underground and implementing USRBDX detectors to monitor the flux in the surface and in some depth underground. The result of the flux of the latter is orders of magnitudes higher than in the surface. I realized that this is because of a necessary normalization that fluka calculates from XXXXXX.sur file for the 100 layers model. Is there a way to extend this normalization to the new layers or it should be done manually?

  2. mgdraw scoring: Following the same idea above, I want some other particle information when muons cross a defined boundary. The mgdraw it seems to do the work but when I plotted the muon position obtained by XSCO, YSCO and ZSCO I should expect the surface delimited by the three cutcones in the geometry where the muons are crossing instead of a cylinder as I’ve obtained (see attached). Am I doing anything wrong?

  3. Basic GCR inputs: The last topic is about the SPECSOUR using the GCR-IONF. I have two doubts on WHAT(5) and WHAT(8). For the first, is the spectral index of the GCR? 1.75 like in the example isn’t harder? And for the latter, when the user prepares the files generated by the atmloc.f is asked as input the geomagnetic cuttof. What exactly the WHAT(8) is used for and, as in the example the value 11.4, is not too high?

Thanks a lot for your time.

mgdraw2.f (6.7 KB)

Dear @tuneu ,

In reply to your questions:

  1. Indeed the normalization factor is applied by default only to the (100) standard atmospheric layers. If you want to add more layers, you can bypass this step and apply the normalization factor in post-process.
    Firstly, make sure you select the option NO-NORM in the GCR-SPE card. Please note that you may need to rename the *.spc files as 01NO-NORM.spc, 02NO-NORM.spc, etc, in order to run the simulation.
    Secondly, in the .out you can find the global normalization factor FLUXST (integral over energy and angle) and the equivalent flux per primary particle FLUX. The normalization factor is then given by the ratio FLUXST/FLUX. Please note that your USRTRACK/USRBDX scoring must be properly normalized. For example, if you don’t provide the volume in the USRTRACK card you must add this factor in the final post-process normalization.
    Below an example of what you can find the *.out:
  *** Fluxst:  0.651134312 ions/cm^2/s
  *** for full earth shadowing ***
  *** <Z_ion>:  1.16797233 ***
  *** Emission radius:  644899968. cm ***
  *** Flux:  3.82679652E-19 Part/pr/cmq ***
  1. Please make sure that the regions that you have specified in the mgdraw are correct. Without the *.inp is difficult to cross-check. Please note also the large curvature radius of the model.

  2. Please note that the numbers provided in the example shouldn’t be considered as universal/all-purposes values. You may change these parameters depending on your specific problem/application.

Cheers,

Angelo

Dear @ainfanti ,

  1. Thank you very much for your quick answer. I’ll run with the NO-NORM in the GCE-SPE to compare the results as you suggested! It’s been very useful.

  2. Please find attached the input file and the .geo and .cards if you or someone could take a look.

  3. I haven’t made my point correctly. Let me reformulate the question. The geomagnetic cutoff of the WHAT(5) for the GCR-IONF card will overwrite the cutoff inputted in the atmloc.f to generate the .geo and .cards files?
    And, the WHAT(8) in GCR-IONF or GCR-ALLF will overwrite the spectral index described for GCR-ALLF (16.1.2) of the Fluka manual?
    ( For the proton component at energies larger than 100 GeV, using the
    normalization obtained at 100 GeV, a spectral index gamma = -2.71 is assumed.
    A spectral index gamma = -3.11 is assumed above the knee at 3000 TeV.)

Thanks for your time and contributions,

All the best,

Jordi
mgdraw_bdx_tests.inp (4.9 KB)
36n140e.geo (30.1 KB)
atmomat.cards (28.6 KB)
atmlocmat.cards (6.2 KB)

Dear @tuneu ,

I have checked the routine and It looks fine to me. Indeed I was able to reproduce your plot.
Please note that you are requesting the boundary crossing from the water region to the Earth1 region. As I mentioned, please note also the large curvature radius of the earth model. Indeed, the points seemed scattered over a regular cylinder while in reality the shape is not regular and it reproduces the crossing surface between the two region of interest. You might check it by rotating the plot to XZ or YZ (see plots below).

Hope this can help.

Cheers,

Angelo


Dear @ainfanti

First of all, thank you for your answer!
If I’ve understood your point, I should expect the boundary between water and earth1 regions all around the Earth’s curvature, but it’s not. Please, take a look at the geometry (figures XZ and YZ) where the regions declared are delimited by the three cutcones (purple lines). The Earth curves at these delimited regions almost in a straight line.

Also, I’ve tried to reproduce the histogram given by USRBDX with the muons obtained through mgdraw getting a flux four orders of magnitude higher from the URBDX’s. Please find attached also the histogram comparison and the script for the construction of the bins.

Thank you!

All the best




bins.py (559 Bytes)

Dear @ainfanti

Before you elaborate any reply, I think I understood my mistake in both cases, the regions and the histogram.

For the former, I understand that what I’m scoring is the particles crossing the wall surface of the cone while I was expecting the surface of the sphere. What I want is to consider the particles crossing from one region to another but from the surface of the sphere. I suppose I should change the expression for the regions to get this information, is that right?

For the latter, I haven’t checked the weights for the scored particles and this, probably, is causing the difference. I’ll plot again today considering the weights and I’ll back here with the answer.

Cheers,

Jordi

Dear @tuneu ,

Indeed I think there is a bit of confusion between the bodies that define a region and the REGION (and its surface) itself.

Indeed, from the routine you have linked you have requested the boundary crossing between

IF(MREG.eq.101.and.NEWREG.eq.102)WRITE(51,*)ICODE,jtrack,
     &   emuon,eprim,idprim, xsco, ysco, zsco, Atrack

From the *.out:

101  Water      127  WATER     OFF           -1.00000E+02     9.99852E+04
102  Earth1     126  MOLASSE   OFF           -1.00000E+02     9.99852E+04

Therefore, the boundary surface between the two regions is consistent with the plot you got.

Nevertheless, if you are interested in scoring the boundary crossing over a different surface (not clear to me what you mean with “sphere” since there are multiple SPH used for region definition purpose) make sure that the regions number (and their definition) is consistent in the mgdraw.f

With regard to the histogram, I did not look yet at the details but it looks to me a possible normalization/scoring issue. Please note also that in the USRBDX you have requested the two-ways particle current (I2) while it looks to me that you are interested in the fluence scoring. Please have a look to this useful lecture of the last FLUKA online training → Scoring II

Hope this can help.

Cheers,

Angelo

1 Like

Dear @ainfanti

What I meant by sphere surface is to obtain, let’s say in a simple geometry, particles crossing the earth’s surface. For example in two regions delimited by a sphere, where region1 is the atmosphere and region2 the earth counts the hits or particles crossing from region1 to region2.

About the histogram, indeed, was a problem of normalization related to particle weights.

Thank you for your help,

Cheers,

Jordi