I’m a first time Fluka user and I need some guidance on scoring secondary particle information. I’m simulating a simple setup of an electron beam incident on a tungsten target. In this setup, I’m interested in scoring the position and momenta of the K+ K- pairs created from the photoproduction and subsequent decay of phi mesons.
In my attached input card tut1.inp (1.2 KB)
, I’ve enabled photonuclear interactions and applied biasing. By Scoring Kaons± with USRBIN, I’m able to see that the K+ K- pairs are being produced. I understand that I need to edit USDRAW in mgdraw.f in order to score the positions and momenta of the K+ K- pairs. Is there any material that explains how I can modify mgdraw, step-by-step, in order to obtain my desired output?
sorry for the delay. There is of course good material to teach you how to use mgdraw and other advanced FLUKA routines. The most relevant would be the advanced workshop offered a few years back, for which all the material is freely available here: .
For your specific scenario, I would pay particular attention to: background: module 2 and module 3, providing an overview of the fluka user routines, as well as some fortran basics (needed since mgdraw will need to be written in fortran) implementation: module 10 is the most relevant here, and I highly encourage you to look over the associated exercise, which you can find here under the section “User Routines II (MGDRAW): exercise”, see the attachments for that section.
I believe with this information, especially the exercise and its solution that could provide a backbone for the USRDRAW entry you want to write, you will be able to score secondary kaons as you intend.
Thank you for your helpful response. I’ll be doing a deep dive into the advanced workshop that you shared.
In the meanwhile I’d like to ask one more question that is related to my original post but perhaps not related to scoring. I checked the .out files using the viewer function in flair and I noticed that although charged kaons are being produced in the simulation, the number of charged kaon decay products per beam particle is 0, see the screenshot below:
Unstable particle decay is included by default, with no need for an additional card.
Charged kaon decay products are mostly muons, neutrinos, pions, positrons/electrons (and not charged kaons themselves), which you can actually find in the above list.
I apologize, I think my previous comment was poorly written. The process I’m interested in is phi → K+ K-, which I believe is included in Fluka. I can see the K+ K- are being produced in the simulation. I became puzzled once I realized that the K+ K- were not listed as decay products when I check the .out file because I am expecting K+ K- to be the decay products of phi.
Sorry, it was my misunderstanding, your mention of decay of phi mesons was actually clear.
In the above list, you can find only the decay products of transported particles, for which decay is an independent process in the course of their transport.
On the contrary, the decay products of very short-lived particles - not transported but decaying at production - appear as secondary particles generated in the inelastic interaction that went through the very short-live particle temporary production.
Thank you both for you help so far. I’ve managed to get fluka to output secondary particle information and I can see that the phi → K+ K- are being produced as expected! I’ll attach my files here in case they are useful to other readers: mgdraw_phiKK.f (10.5 KB) phiKK_4500.flair (2.0 KB) phiKK_4500.inp (1.1 KB)
I have just two more questions:
As you can see in mgdraw_phiKK.f, I’m able to dump information about the momentum of the K+K-, but how can I also access information about the position of the decay vertex? I could not find anything pertaining to this in genstk.inc.
Since I’m mostly interested in only the phi->K+K- process, is it possible to speed up the simulation by preventing omitting other processes / interactions?
great to hear you are making progress. Here are my thoughts on this:
the entry you are using, USDRAW, takes in as arguments XSCO, YSCO, and ZSCO. These are i fact the position of the interaction vertex (as Francesco said above, very short lived meson decaying to Kaons would not be transported, so they all appear together at the interaction vertex). You can print these variables.
I am not aware of a method to bias towards the production of a specific meson in inelastic interactions in FLUKA. However, if you know the pathways you are interested in the production of the Kaons (ignoring the non-transported phi meson, e.g.), you could set high thresholds for other particles, so you don’t waste CPU time on them.
Other notes:
You already use LAM-BIAS for PHOTON, which will increase all photonuclear reactions, which is good.
To speed up the simulation, if you are finding that the files generated by mgdraw are big, it could speed up to reduce the writing by filtering the events. If the simulation still proves to be slow, you may try to implement more advanced techniques, like importance biasing, 2-step solutions, etc.
Thanks for another helpful response. I’d like to follow up on suggestion 2, how can I set these thresholds / where can I learn more about this. In addition to increasing the threshold for particles I’m not interested in, I realized that I’d like to decrease the threshold for phi->K+K-. Below, I plot the energy distribution for the photoproduced phi mesons (really this is just the total energy of the K+ K- system). It appears there is a threshold set at 4 GeV for this process. Is it possible to lower this threshold?
The threshold I was referring to are transport threshold, not production threshold. In hindsight, given your simple geometry, changing the threshold will not speed up the simulation significantly, as very little time is spent transporting the particles. The hadron production thresholds are not something that can be changed.
Your result is not conveying the existence of a production threshold in FLUKA, but rather of an artificial threshold enforced during the mgdraw printing process. The logic you implemented checks if there are exactly 2 particles in the stack, meaning that it is enforcing that no other particles or evaporating residual nucleus can be produced. This might or might not be what you are looking for; however, you will be cutting out a large amount of events where K+K- are produced, but are accompanied by other reaction products.
If instead you are interested in the reaction
\gamma + A \rightarrow K^+ + K^- + X
where X is any other reaction product, then you should modify your mgdraw to loop through the genstk, and see if the two kaons are present. If so, to then print the results. When I implement this logic, I obtain a plot that looks like this, where you can see that the total kaon energy (the one you call phi energy) extends across the entire allowed energy range.
Thank you for another thoughtful and helpful response. I’ve modified mgdraw to allow for \gamma + A \rightarrow K^+ + K^- +X, as you have suggested. The distributions below show the total Kaon energy for two beam energies, 3.7 GeV (left) and 4.5 GeV (right).
These are the same as your histogram, except summed along the columns.
The distribution on the right is peculiar. It seems like the simulation behaves differently when the total kaon energy can be greater than 4 GeV (or, when it is possible to produce a 4 GeV phi meson which immediately decays to K^+ K^-). In this case the interaction tends to only have two secondary particles (K^+ K^-), whereas below 4 GeV the interaction will always have several more secondary particles. Do you know why this might be happening?
you have observed indeed a peculiar behavior, and you are thanked for bringing this to our attention. What you are seeing when comparing the two runs is the onset of the elastic photoproduction of the phi meson, which is only turned on above 3.5 GeV.
At the energy you are considering, photons can fluctuate into vector meson states, like \rho^0, \omega, and \phi, a regime of photonuclear interaction known as Vector Meson Dominance (VMD). The vector meson can then scatter with nucleons inside the target nucleus: if sufficiently low kinetic energy is transferred, the vector meson can leave the nucleus in an unexcited state (elastic); if more energy is transferred, other secondaries are produced (inelastic). Right now, the elastic cross section is only computed for energies E_\gamma > 3.5 GeV. Your two simulations were right around this energy, and that is the reason why you observe a sudden spike.
It is certainly a point that warrants further investigation on our side, in order to cure this too abrupt threshold. We will look into it and we will let you know as we implement relevant changes in the code.
I just want to clarify, is the elastic cross section only computed for E_\gamma > 3.5 GeV or E_\gamma > 4.0 GeV? I think my finding suggest the latter.
Also, is it the case that the elastic cross section has this threshold whereas the inelastic cross section does not? Below I plot the invariant mass distribution for the charged kaons in my 3.7 GeV simulation. It’s clear the phi mesons are being produced, but always with additional secondaries (inelastically).
indeed your intuition is correct. While the VMD elastic scattering does have a threshold at 3.5 GeV, (as you may want to check by looking at \pi^+ \pi^- outputs) there is a further rejection on the \phi meson specifically, which effectively raises the threshold of elastic production to 4 GeV.
Regarding additional thresholds, there are of course natural thresholds related to the mass of the phi meson and other secondary particles. Could you clarify the connection of the figure you showed to the inelastic thresholds?
The figure was just showing that inelastic phi meson photoproduction is still being simulated in the 3.7 GeV case, even though elastic vector meson production in not simulated due to the threshold that you mentioned. I was a little confused by this but I now understand the threshold is only for the elastic interaction.
Is it at all possible for users to modify the 3.5 GeV threshold on VMD elastic scattering and the 4 GeV threshold for phi mesons? Or will I have to wait until it is looked into on your side?
unfortunately it is not possible to modify these thresholds directly… I will have to look into that and see what other effects on the code such changes would imply. I will keep you posted,
