I’m simulating cosmic ray cascades and I’m interested in the muon production point and the muon flux at some height where a detector is set. To do so, I’m using the user routines USRDRAW and STURPF.
In the USRDRAW routine, I activate the BXDRAW entry which gives me the muons crossing from a region to a new region (line 74 in the usrdraw file attached).
In the USRDRAW entry, I get the muons creation point, energy, etc. by recording them using the muon ancestor ID’s (line 173 in the usrdraw).
I’d like to know where, muons crossing some plane, were created and reconstruct the atmosphere profiles of the muons that reached the detector. I was trying to flag them using the integer variable LOUSE in STURPF, so I can identify if the muon created is the same as reaching the detector.
Unfortunately, when I implement a counter at my USRDRAW entry I get for each muon the same LOUSE number = 1. So in each interaction, my variable LOUSE reset to 0.
How can I implement this counter correctly to identify each muon individually?
I think the problem is related to the way you are setting the value of LOUSE.
First of all, who is LOUSE? LOUSE is a variable from the FLUKA COMMON FLKSTK, which contains information about the particles in the stack (not all of them). You can find more detailed information about the different stacks in FLUKA here: Particle stacks in FLUKA.
When you are trying to redefine the value of LOUSE, you are doing so to a variable in FLKSTK that is not even related to the muon you want to follow (and it may not necessarily be a muon). When you transport the muon, the information in the FLKSTK is copied to the COMMON ‘trackr’, and the stack index in FLKSTK is decreased by 1. This means that you are trying to modify the LOUSE value of another particle.
Therefore, I recommend that you properly define the value of LOUSE in your stupr.f routine. I noticed that you were trying to define a COMMON ‘FOLLOW’ with the variable LLOUSE_COM. Did this not work for you?
In the USDRAW routine, you should be dumping the LLOUSE variable. Here, you want the information from the TRACKR COMMON, where the value LOUSE was copied when the muon was unloaded to be transported.
Also, I did not understand why you need the USDRAW routine and why you want to increase the value of the muon there. Don’t you want to flag a muon when you start to transport it?
Here is one way you can proceed (if I understand correctly):
In stupr.f, you can dump the information of the muon that was created (position, energy, etc.) and define the value of LOUSE that will be copied in LLOUSE.
With the BXDRAW routine, you can catch the moment when the muons cross the regions and dump the relevant information plus LLOUSE to identify the muon.
I hope this has been clear. Let me know if you have any further questions.
Thanks for your reply. Reflecting on your points, maybe the approach I’m trying is not the best option:
Also, I did not understand why you need the USDRAW routine and why you want to increase the value of the muon there. Don’t you want to flag a muon when you start to transport it?
I use the usrdraw routine to get the point, energy, etc. where muon was created. The idea of that is to write in the output all the positions where muons appeared for the first time. Just one point, without considering further transport of that muon. Then, the muon will reach the detector or not. The bxdraw gives this information. What I’m missing is how to relate these muons from usrdraw and bxdraw.
If I understood correctly, only using the bxdraw routine could be enough? For example, if one muon reaches the detector can I get the muon position, energy etc. where it was created? When I used stuprf.f for this purpose I only got the muon ancestor which is itself one step before during the muon transport.
Is there any example to extract this information?
When the routine STUPR is called you are loading the secondary particles from the GENSTK to the FLKSTK. You could use the information from the GENSTK. In particular the one related to the particle IS is “KPART”
Thus, you can do the following:
In your stupr.f, create a counter that will be incremented every time you load a muon. To do that, you could define an INTEGER with an initial value equal to 0 and with the command SAVE to be able to maintain the modified value after calling the routine:
INTEGER COUNTER
SAVE COUNTER
DATA COUNTER / 0 /
Then, when you try to load a muon:
IF((KPART(NPSECN).eq.10).or.(KPART(NPSECN).eq.11)) THEN
COUNTER = COUNTER + 1
LOUSE(NPFLKA) = COUNTER
WRITE(60,*) ... [dump the info. Here you can dump the initial values of the muon**]
ENDIF
I dumped kpart and tki from GENSTK with index (1) and XTRACK, YTRACK and ZTRACK from TRACKR with index (0). The position x, y and z is not avaiable in GENSTK, right? Those indexes will return the initial position and energy of the muon?
The last column is the COUNTER implemented.
But, I’m not sure if those muons are new. For example, the muon number COUNTER= 4 it seems to be propagated in the medium decreasing slightly in energy and altitude (ZTRACK). Did I declare something wrong?
The first time a muon crosses the detector ICODE 19 is the one with COUNTER= 707. This part worked as expected
The only problem I see is the way you are dumping the information in your stupr.f routine.
If you want to dump the information of the muon you are flagging, you should do the following:
I’ve implemented the way you suggested to dump the information but nothing has changed. I also included the muon parent ID to verify that it is transported ( parent ID == ID) as we see from Counter numbers from 5 to 21. The Counter = 22 it’s a new muon indeed. In the output below there are 5 muons (1,2,3,4,22).
Taking into account all this I set a new condition eliminating the lines where the parent ID is not equal to ID. The output now seems correct showing the first 5 muons without being transported.
I don’t know if this is the most direct way to get this information but is working, so any further comments are welcome. Meanwhile, I gonna mark your second answer as a solution but all your comments have been very valuable.
Be careful in the way you are filtering the muons that are “not new”.
When STUPR.F is called in one of these cases, the value of LOUSE might be reassigned as 0, losing the flag that you gave previously.
You should be able to check the value of LLOUSE (from TRACKR COMMON) in STUPR. In LLOUSE you have the flag value of the parent (the muon) that you assigned previouly. Then, you must pass in stupr the LLOUSE value to LOUSE(NPFLKA) when this muon is loaded into the FLKSTK.
I hope I am not confusing you. The only thing that I tell you is that you might be losing the track of your muons because the LOUSE value is redefined everytime a muon is copied from the GENSTK to FLKSTK. Then, you have to copy the LLOUSE value of the parent muon to the LOUSE value of this new (“not new”) muon.
Be careful in the way you are filtering the muons that are “not new”.
Do you refer to the condition parent ID .ne. ID or to use the numeric indexes e.g. TKI(1). In the indexes case, I’ve implemented TKI(NPSECN) as you suggested.
Anyway, I’ll be alert on the outputs counter numbers to check if it resets and take the action you’re commenting if it is so.
