I am currently simulating the simple setup of a 3.2GeV electron beam imping on a tungsten target with regards to electro-nuclear interactions (as mentioned before elsewhere). I have some technical questions regarding the simulation.
At the moment I want to dump all particle information into an output file.
What I would expect to see is - if we take 5 histories - five electrons and their respective secondaries, depending of course on their type of interactions. However my output files contain
almost randomly some of the 5 “events”. The output of one run is given in the cycles.txt file appended below.
Does this simply mean that the electrons of “events” 2, 3 and 5 did not interact? Or am I missing something in the logic of the mgdraw subroutine. I will include also my mgdraw.f file and myfluka.inp.
I have encountered the ICODEs to filter on specific interactions. Ideally I would want all particle information (including the recoiling electron) from electro-nuclear interactions. So if there is an event where the beam electron interacts electro-nuclear-ly I want all particles from that event to be printed in the manner explained above. From the list of ICODEs I am not sure if there is such an identification of electro-nuclear processes, if there was I could simply add a if-statement around the write part in the mgdraw.
Is there a way to apply some sort of event selection? E.g. I want to include only events where the recoiling electron has only 25-75 MeV of the orginal 3.2GeV left and has an angle theta less than 0.2deg. Is this done in stupre.f and if so how would I go about obtaining that information?
thank you for your question, I will go through it point by point.
entry USDRAW in mgdraw is called as soon as any kind of interaction has taken place in your simulation. Your result shows that some of the primary electrons don’t interact in the target region, therefore USDRAW is not called and no printing in the output.dat file is done. In fact, if you would score the particles exiting the target region, you would see that some primary electrons survive and don’t interact in the target.
you can find already the answer to this question in one of your previous post. Repeating what said there, the correct filtering for electro-nuclear reactions in USDRAW is ICODE = 101 .AND. JTRACK = 3. If you are then interested in the secondary particles produced in the electro-nuclear reactions, you can access their properties with the variables in common genstk.inc
continuing the previous answer, in genstk.inc you can find also variables describing the laboratory kinetic energy of the secondary particle (TKI(I)), and its direction (CXR(I), CYR(I), CZR(I) with I = 1,NP). These variables are accessible from USDRAW , so you can implement the desired logic already there before printing in the output file (no need of struprf.f).
Hello Giuseppe!
I have been looking at the output for electro-nuclear interactions from my simulation and have come across another question.
The way I understood it from your post is that when using the secondary information provided by genstk.inc I am able to reproduce each event. And indeed I get a list of secondaries for each electro-nuclear event. However when I add up the momenta in x, y and z direction of the secondaries I never obtain the initial values of my beam particle: (0,0,3.2)
The secondaries do not include the recoiled nucleus though, which could cause this issue. Do You have an idea how I can obtain the momentum information for the recoiling nucleus?
I get the momentum three-vector of my secondaries via CXR(IP)*PLR(IP).
I would like to show that momentum is conserved with the results that I get.
Thanks in advance,
Laney
as you correctly state, in genstk.inc you can find the information of the secondary particles (not “heavy”) produced in an interaction.
Being interested in the properties of the residual nucleus generated in the electro-nuclear interaction, you can access them with the common resnuc.inc.
I had a quick look and considering the nucleus too, I see that the total momentum is conserved. The difference between electron and products momenta reaches a factor ~1e-10, which I am assuming can be associated to numerical precision.
Thanks to your help I have now been able to add the momenta of all my particles in order to check momentum conservation. My idea was to simply add all momenta of secondary particles in x, y and z direction per event and thus obtain the three momentum of my beam electron, which is (0, 0, 3.2).
Now when adding up the momenta per event I do not seem to reach the same accuracy as you did. For x and y it is around 0e-3 and for the z direction I get 3.06GeV, which is of course very far off.
Maybe I am using the wrong momenta. Is the recoiling electron momentum part of the secondaries? Or do I take the information provided via CZTRCK*SQRT(ETRACK^2-0.000511^2) (added the ^ just for readability)?
Other than that I use the momenta provided by CZR(IP)*PLR(IP) (for the secondaries) and PZRES for the nucleus.
I will include the momentum distributions I have obtained from my simulation.
I would be very thankful for some help!
Kind regards,
Laney
Hi Giuseppe!
Quick update on my momentum conservation check. I have realised that the way I did it in the previous post I do not account for energy loss of the initial electron before the electro-nuclear interaction. Of course if the electron already interacted in the target before the electro-nuclear interaction the total momentum will not be the (0, 0, 3.2GeV) vector. My idea was thus to take the momentum of the electron before the interaction and subtract this information from my secondaries vector. Now the secondaries momentum takes the secondaries momenta from CZR(IP)*PLR(IP) and the nucleus PZRES.
I am not 100% sure where to obtain the electron momentum before the interaction, but I just guessed and took the TRCK information (so CZTRCK and ETRACK). With this approach I do get a somewhat symmetric distribution in x, y and z. However I still do not have the exactness that you got. Am I still missing some information?
sorry for the late reply but indeed you managed to find the problem in your calculations.
You need to take into account the momentum of the electron undergoing the electro-nuclear reaction and not simply the properties of the primary particle (since the electron can be subject to other reactions before the one you are interested in).
Now, to retrieve the information on the electron undergoing electro-nuclear reaction, you are correctly looking into the trackr.inc. The only thing I can comment is that you are considering ETRACK (= total energy of the particle) and not PTRACK (= the momentum).
For reference, to obtain the results I was reporting, I looked at:
electron reacting: PTRACK*CZTRCK from trackr.in
secondary particles: PLR(I)*CZR(I) from genstk.inc
PTRACK (when available, i.e. positive) is SQRT(ETRACK**2 - AMELCT**2) and refers to the electron right before the interaction (which indeed is not necessarily the originally aligned electron).
Momentum conservation is assured within the machine precision, namely down to negligible discrepancies.
Hello Francesco and Giuseppe!
I am still working on the momentum conservation part. Somehow I am not able to prove it with the precision you imply. The momentum I obtain is:
Which is far off from machine precision I would suppose. I am using PTRACKxCZTRCK, PLR(I)xCZR(I) and PZRES. Adding the latter two and subtracting the PTRACKxCZTRCK. Though I am only writing all variables to a file and then processing and calculating the momentum conservation via a separate c++ script. How can the distributions still be this broad?
I would be happy for some more help,
kind regards
Laney Klipphahn
With PLR(I)xCZR(I) do you actually mean Sum_I=1,NP [PLR(I)xCZR(I)] as it should be?
If so, please upload your mgdraw and c++ script in order to find your mistake.
Dear Francesco,
yes indeed I am using all secondaries from IP=1 to NP. But I will provide my mgdraw: mgdraw_empty.f (10.3 KB)
And also the script used to visualise the data: plotting.cc (17.8 KB)
I am very thankful for the help!
Kind regards,
Laney
Your mgdraw_empty looks correct (although it’s not ideal to repeat for every secondary the unchanging info concerning the projectile electron and the residual nucleus).
However, the cc script seems to suffer from a basic issue in the reading part, where you consider 7 integer values (col0-col6) instead of the 6 actually available.
Ah yes reason for that is, that I add a event ID to the output.dat file that starts from 1, such that I do not get any double event IDs from the 5 cycles. So the output.dat files for each cycle is merged into one file with an additional event ID that starts from 1 and goes up to the number of total eNuc events in all 5 cycles. This I do via: mergefiles_numbers.cc (1.5 KB)
Here is a mgdraw_empty.f (10.1 KB) version checking the momentum conservation and demonstrating its accuracy.
Note that I included secondary particles in the fheavy storage, where occasionally you may find some additional reaction products.
This works under the assumption that no LAM-BIAS was applied.
Otherwise, one should not overlook that among the secondaries in genstk the projectile electron is also found as survivor, and has to be excluded from the computation.