Arbitrary distribution for a user defined source

Advanced Features and User Routines
1

Dear FLUKA experts,

I want to create a source distribution as an arbitrary function, e.g. f(x,y) =A exp(-x^2 - y^2) and then implement it in my FLAIR simulation, but I don’t know how to go about it. Any help would be appreciated.

Best Regards,
Carlos

2

Dear Carlos,

you need to use a source user routine, to code your arbitrary function.

If you didn’t use it before, I suggest to try the updated version: source_newgen.f

The following lecture introduces its usage, and it has some basic information on Fortran programming:

Cheers,
David

1 Like
3

Dear David,

Thank you for your answer and sorry for my late reply.
I used source_newgen.f as you recommended, but when I try to use an arbitrary distribution (or even the sample_flat_distribution) I get this error.
Screenshot from 2021-01-25 13-56-51

I attached my files for clarity.

basicsource.f (3.7 KB) basicsource.flair (9.4 KB) basicsource.inp (1.2 KB) source_library.inc (30.2 KB)

Any help would be appreciated.

Best Regards,
Carlos

5

Dear Carlos,

the explanation is simple. Your beam is not where you expect it. So, the all values in the USRBIN scoring are zeros, and Gnuplot can’t plot it with a logarithmic color bar.

Luckily, you are using the new source_newgen.f user routine, and this has a debugging option. If you set the debug_logical_flag variable to .true., then it will print basic parameters of the first 100 (changeable with the debug_lines variable) calls to the .log file. This will help you figure out what is going on.

In your case you are expecting the beam around x = 10, but you are overwriting this, when you are assigning a new value to x with your cosine_distribution() function. To keep the original position and sample around it, simply use:

coordinate_x = XBEAM + cosine_distribution()

With this modification, you will be able to plot the beam particles.

One more remark:

You are using the source_newgen.f which distributed at the FLUKA course. A newer version is available by default with the latest 4-1.0. You can find it in <fluka>/src/usr directory.

With this version the source_library.inc doesn’t need to be in the same folder, you don’t need to modify it. You can add your own function / subroutines to the end of the source_newgen.f file as well.

Cheers,
David

6

Dear David,

Thank you for your answer.

I did the modifications that you told me and It worked, but now I am trying to use sample_annular_distribution subroutine but I don’t think it’s working unless I am missing something. I attach my files.

basicsource.f (4.9 KB) basicsource.flair (3.2 KB) basicsource.inp (1.2 KB)

Also, I was having an error when compiling basicsource.f and I am unable to replicate it now, but it’s similar to this warning. If you could also take a look at it.

compilingerror.txt (1.2 KB)

Thank you for your time.

Best Regards,
Carlos

7

Dear Carlos,

you do not need to put your variables into “[”, “]” brackets, they only supposed to indicate which variables has to be replaced by a user defined one. Probably, this needs to be more clearly written in the guide.

There are two other issues:

  1. In the source routine you divide the radii of the distribution with 1000, which gives you values in the order of micrometer, this won’t show up in your scoring.

  2. In the plot, you use a projection along the Z axis. This will give you average values. You could

    a. set the range of the projection between Z = 10 and 10.05, to get only one layer, or

    b. set the diverge of the beam to zero.

Cheers,
David

8

Dear David,

Thank you for your input, about the issues:

  1. It worked perfectly after I realized what you said by chaging my input parameters
    2.a. It worked perfectly when getting one layer
    2.b It also worked perfectly by setting the divergence to zero.

Now, going back to my first post, If I wanted to make my distribution f(x,y) =A exp(-x^2 - y^2) and following this as reference https://flukafiles.web.cern.ch/manual/13.2.20.html
How could I go about normalising that 2 variable function to implement it in my code? If you could point me into the right direction.

Thank you for your time.

Best Regards,
Carlos

9

Dear Carlos,

could you elaborate what distribution you want to implement?

The f(x,y) distribution you mentioned, looks like a 2D Gaussian distribution, where X and Y variables are independent. In that case this distribution can be replaced with two separate 1D Gaussian ones.

Cheers,
David

10

Dear David,

Thank you for your answer and please ignore my last question, it’s badly worded.

I have been messing around with my cosine_distribution function and it seems to make sense, but if I change my input parameters starts behaving strangely. I would think since I am using f(x)=cos(x) as my function (and following the previously linked flukafiles link for it’s normalisation) I should be able to get more than one “peak” but I can’t. I attach my files.

basicsource.f (4.5 KB) basicsource.flair (3.2 KB) basicsource.inp (1.2 KB)

Please let me know if my wording is unclear.

Thank you for your time.

Best Regards,
Carlos

11

Dear Carlos,

unfortunately I can’t understand the cosine distribution in your source routine. Is there any reference what you want to code?

However, I can elaborate on how to sample from arbitrary 1D distribution.

Let’s try to sample the X random variable from the cos(x) function as a distribution.

f(x) \propto cos(x)

but we need to limit the value of X to

x \in \left[ -\frac{\pi}{2}, \frac{\pi}{2} \right]

because the probability can’t be negative, and Fortran uses radians.

First we need to calculate the Probability Density Function (PDF) f(x) of the distribution as such:

P(a \le X \le b) = \int_{a}^{b} f(x) dx

The integral of the PDF on the whole range of x, should be equal to 1.

P\left(-\frac{\pi}{2} \le X \le \frac{\pi}{2}\right) = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} f(x) dx = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} A \cdot cos(x) dx = 1

where A if the normalization factor. After doing the integration we get

A = \frac{1}{2}

So the PDF is

f(x) = \frac{1}{2} cos(x)

Now, we can write up the Cumulative Distribution Function (CDF) as the following:

F(x) = P(X \le x) = \int_{-\frac{\pi}{2}}^{x} f(x') dx' = \int_{-\frac{\pi}{2}}^{x} \frac{1}{2} cos(x') dx' = \frac{1}{2}\left(sin(x) +1 \right)

The value of the CDF is always between 0 and 1.

Then we can take the inverse function of the CDF, which is in this case:

F^{-1}(\xi) = asin(2\xi - 1)

Finally, we can sample \xi from an uniform distribution, where

\xi \in [0,1)

Doing so, the inverse CDF function will give us an x value according the f(x) distribution.

I implemented this in your source routine: basicsource.f (4.7 KB)

Since I wasn’t sure what you intended with the ini_ and fini_ variables, I used those to rescale the distribution, in such a way that the distribution ends at those values.

I hope this clarifies everything.

Cheers,
David

1 Like
12

Dear David,

Thank you for your explanation, it does help a lot.

I have been messing around with your implementation and what I am thinking is if there is way to show multiple times the same distribution, e.g.:

I want to sample from 0 to 20 and the beam is centered at 10. Without rescaling the distribution it should go from (10-π/2) to (10+π/2). Looking again to my x axis I know that from (10+π/2) to (10+3π/2) and from (10-π/2) to (10-3π/2) my distribution would be negative so I want to show this empty. Then I would have 3 “beams” one centered at (10+2π), one centered at (10) and one centered at (10-2π).

Something like this,
basicsource.f (6.6 KB) basicsource.flair (3.2 KB) basicsource.inp (1.2 KB)
But without having to define 3 beams.

Does this make sense?

Thank you for your time.

Best Regards,
Carlos

13

Dear Carlos,

you can set up a random number, which would give -1, 0, or 1 like

nrandom = int( FLRNDM(xdummy) * THRTHR ) - 1
sample = random_cosine + TWOPIP * nrandom

Cheers,
David

14

2 posts were split to a new topic: Analytical cosmic ray muon source sampling