Re-simulate a specific primary history

Dear colleagues,
in a previous project I was working in, I was using GEANT4 to simulate a large number of events. For each event, I wrote in the output the relevant simulated data, as well as the event number and the status (i.e. the random seeds) of the random-number generator at the very beginning of the event.

At the end of the simulation, I analyzed the event-by-event output that my simulation was producing, and I identified a set of “interesting” events. I was then able to re-simulate them, by setting the status of the random-number generator (i.e. the random seeds) to that corresponding to each “interesting event”. In particular, I find it useful to use intrinsic G4 verbosity (/tracking/verbose 1) to investigate the event itself.
I was even able to modify my G4 application to add verbosity (via cout) and still re-simulate the “interesting events” - specifically, I learned that as soon as the “sequence” of random number extractions in the history was not altered, the same event could be reproduced in this way.
I found this procedure very useful to study specific “rare” events.

Is there any way to do so in FLUKA? I know that, by setting the RANDOMIZ card appropriately, I can reproduce all histories for a given run, but I’d be interested in checking the same for a specific event.

Thanks,
Bests,
Andrea

Hi Andrea,
enable through the START card the generation of the random seed file at the beginning of each event (WHAT(5)) and link this usreou.f user routine coded by @cesc that saves it at the end of the each event (but the first) embedding the event number in the random seed file name.

1 Like

Dear @ceruttif, sorry for late reply.
I tried the code you suggested, and indeed I see that, at the end of the short run I tried, for each event I have a file whose content is similar to the following:
10C24FA 0 4D2 3039B66B018E 55 1596DBC4BC3FEC4DAB 92FD02983FD023A53455B8003FDF5D60E2707B003F7566C524FAF3733FE400B4 F3AE19CC3FC3FB1821F509AC3FC0D428B6E2BB903FE5AFDFBC1C76913FED288B D0DEB4803F902BE06D3327663FD9E9653B14B8893FEFCDE0653CD0383FC375C2 E6906AEC3FE07FD669C8430D3FEAAB12D0B3637C3FE39C787F3EEDFC3FEABE76 44D857DD3FE1D190E3554CDE3FEEE977E1F4E2743FD128773CED9BFC3FE080DF 29A7E143FCE2960E31CCC293FE6B854B141551F3FE87A106A405B503FE406BF FD3ED4803FBA2EF5DB60520E3FEFBD9F9D3DA8EF3FE0D9E0AF2CC93F3FE0890D D9875BF03FE6E917212D64003FA5858130EFB0D43FDBB31D4D5B66AE3FEF0802 A3310BEA3FE2005A53B4009D3FEE2C66FFE458763FEC2CFA44EF6BD43FD4820C 81547173FE66F07420044A13FE548E1 1F4E2243FD12D2FB420E4403FC336FB D0C5A1FF3FE1B180680587BB3FEA5163CF111ADC3FC0E5193DDD09543FEEDCA5 40AD12943FC51496E207A4A63FE5F7118CED86173FEF68B5F352EB783FB0CF8E 5AAA31BB3FE1988E1772DBF03FEB8B86D044D9CA3FE016AA9926C2A03FD0CB76 10CA898D3FE4FFD94B9FFFB03FBA8DCEAC0FF9423FDBB5E7474C91A03FA1EDB6 F98C43793FE9CAE6365499AA3FE72BE6C6C5351A3FE841481D80999C3FCD1489 2338A4B83FD7DC12B240D2943FCC9D61AF3D88C03F8F2846661DDC1B3FEB712C 5949B383FD992879DE4A4403FDB4FA1230F4C3C3FDCD17176A8F9043FC8EC8C 6C65E0B3FEB988A9BD8F4F03FB9AA5F8E2C03CD3FECEAE5 F63DBD23FD95C47 EC33D02E3FE79718371459603FB689DDE5DEF85E3FDE7981C02EE0423FDBD45A D5DF35D23FDF1709286E792C3FEC163080D584403F9EB3F18571E2713FEC7BAE 7E41F6A83FDC539A 64CCF203FE9035252A33E7C3FE9B532D262E72F3FE6FFD1 7DF3A03E3FD67E25D42AECCA3FDDEACC5C4FD9A83FB77F04BC7629853FE69D81 7143629C3FDCD377B2D719923FDAAC2853EA74ED3FE793A41FCF9D693FEB9520 11BDAB103FCF9848494D00B43FE4BBF7840B71573FED1C09465DBBB23FDD4152 AC42E5F83FB57A01

Is it correct that, in order to reproduce this event, I should replace the content of the file ran* with the above lines, and re-launch FLUKA?

I also have another question regarding this - since in my simulation I use an external system to launch on LXBATCH, and the WHAT(2) field of RANDOMIZ is automatically set for each job to a different value, is it possible to get this value inside any of the users routines?

Thanks,
Andrea

In order to reproduce the event, don’t touch the RANDOMIZ card (leave it as it was in the initial run), rename the random file of interest as ran[INPUT_NAME]002 and launch from the directory containing the latter a simulation (of 1 single primary) starting from cycle 2 (and come back in case of failure :wink:).

The WHAT(2) of RANDOMIZ has to be made available by a core routine. I can pre-compile and upload it based on your gfortran version.

Hi @ceruttif, thanks, I confirm that I can properly re-simulate the event without touching the RANDOMIZ and renaming appropriately the ran* file.

It would be great if you can send me a pre-compiled version of the code where I can get the WHAT(2) of RANDOMIZ. The reason why I’d like to do so is that, in this way, I can save in the output file at the event-level the random-number seeds, and at cycle level WHAT(2). My output file will be self-consistent regarding the possibility of reproducing events.

My gfortran version: 11.3.0

Thanks!
Andrea

Rename this object file to rndcrd.o and link it into your FLUKA executable. In your user routines, you will be able to access the WHAT(2) value of RANDOMIZ via the statement:

 COMMON / RUNIDE / ID_RAN

to be put right after INCLUDE statements. The value is contained in the integer variable ID_RAN.

Dear @ceruttif, thanks, it works! I have two suggestions for the next FLUKA release:

  • Have the WHAT(2) of RANDOMIZ available by default without linking to the object file you provided me :slight_smile:
  • Be able to access the content of the RANDOM seed file directly within the FLUKA users routines - at this moment, I read it by writing to the disk via START at every event, and then reading it back.

Thanks!
Andrea

Dear Andrea,

  • Accessing WHAT(2) of RANDOMIZ within the FLUKA users routines

You can directly printout the WHAT(2) of RANDOMIZE from any user routine (e.g. usreou.f) by adding the following line (no extra COMMON needed):

call flrnoc (iseed1, idum1, idum2, idum3)

flrnoc allows to access seeds values and counters as used by the random engine (without altering its state). The first integer (iseed1) is the value of WHAT(2).

NB 1: WHAT(2) is also printed in any ran[INPUT_NAME]* file: it is the third value (expressed in hexadecimal). 4D2 in the example you shared earlier in the thread.
NB 2: If ever WHAT(2) > 100 000, you will get its value [100 000], but the trick will still work, because that is the value actually used within the random engine and being checked against.
Off-topic 1: RANDOMIZE WHAT(2) / iseed1 is part of a compacted seed, in turn used for setting the 97 initial actual seeds of the random engine.

  • Accessing RANDOM seed file directly within the FLUKA users routines

Here, why not directly use what was shared earlier in this thread (in usreou.f)?

Hope this helps,
Best regards,
Gabrielle

Dear Gabrielle,
thanks for your message. The proposed method to get the WHAT(2) is very clear.

The reason why I think it can be useful to access the RANDOM seed file directly in FLUKA is the following. In my application, I save in the output event-by-event data, for selected events. In order to be able to ri-simulate a given event, for further analysis(*), this content is needed. What I do now is to use START with WHAT(5)=1, to force FLUKA to write the file on disk for each event, and then I read the content with userouf.f (slightly modified). This probably introduces a penalty on the simulation (each event requires to read and write a file to the disk), hence my request.

Thanks,
Andrea

(*) What I mean by ‘re-simulating an event for further analysis’ is that, when re-simulating, I use the same machine where the original run was executed, with the same executable and same input file. I just add many PRINT statement in user routines to be able to investigate the event, taking care to not introduce any new random number extraction that would unavoidably change the event history.

Dear Andrea,

Ok, thanks for the details.

So for point (2), let me first answer with a question, to be sure I fully grasped what you need.

Need for I/O: The analysis step does not have access to the memory from the first run, so you must already have some I/O? How do you identify in the second attempt (analysis step), which event you would like to look at?

I/O granularity: In the first run, do you dump info on a per event basis, or are you able to already identify interesting events (with conditional statements in user routines) and dump info for those interesting events only?
If so: what could potentially be done, is to dump the random engine state for interesting events only?

Thanks in advance,
Cheers,
Gabrielle

Note 1: On the need for I/O
In FLUKA v4, the random engine state at the beginning of an event depends on the status of the random engine as left at the end of the previous event. Hence you need to have I/O (shall it be the random state of the interesting event only) to be able to then directly reproduce that event in a second run.
Of course, a different design choice could have been to associate the random engine state at the beginning of an event to the eventId: in that case, it would have been possible to run events in any order (hence making a run reproducibile in multi-threaded runs), and … allowing to restart an event for analysis by simply knowing its eventId, without any I/O.

Note 2: Dumping the random engine state
In FLUKA v4, it is possible to dump the status of the random engine at any stage to a file with:

call rnwrit ( iounit )

For event reproducibility purposes, this line needs to be placed at the very beginning (source.f user routine) or end (usreou.f user routine) of an event. Note that in that case, you do not need to activate the START card: since these user routines are called for each event, you directly get the status of the random engine at the start/end of each event.

You are then free to put conditional statements to activate those dumps as you which. However, to be able to only dump seeds for an interesting event, one could think of the following:

  • Save the status of the random engine in source.f (at the beginning of the event) to a cache.
  • Identify (in any user routine, or at the end of the event) that the event is an interesting one, and in that case trigger the dump of the cache to a file (random engine state I/O).
    This is not yet possible to do in FLUKA v4.

Dear Gabrielle,
thanks for your detailed reply.
I try to reply to your points.

In my FLUKA based application, I write to in a ROOT TTree relevant data for each event - for example, the energy deposited in each cell of a calorimeter, or the momentum of the generated particle. This is done only for events that satisfy a certain condition (total energy deposited in the calorimeter less than a certain value).
I also write to this ROOT TTree the content of the random file. To do so, at this moment I have the START card with WHAT(2)>1, so that the content the random file is always written to the disk, at every event. Then, for events of interest, I read from the disk the content of the random file, and I save it as a TString inside an appropriate branch of the TTree.

What I do not like in this approach is that for all events I introduce an I/O operation: write to the disk the content of random file, then eventually read it back and save to the TTree
What I was wondering is that if, instead of writing to the disk, the content of the random file at the beginning of each event can just be saved to the memory, and then eventually saved to the TTree.

Just to give the idea, in pseudo-code, what I do now is:

void beginning_of_every_event(){
  write_content_of_random_file(name_of_output)
};
void end_of_every_event(){
  if (interesting){
    string s=read_content_of_random_file(name_of_output);
    tStringOut=s;  //defined somewhere and connected to the TTree
    tout->Write();
  }
  //somewhere
}

What I’d like to do:

void beginning_of_every_event(){
  get_content_of_random_file(s); //s is a string defined somewhere with appropriate scope
};
void end_of_every_event(){
  if (interesting){
    tStringOut=s;  //defined somewhere and connected to the TTree
    tout->Write();
  }

Thanks again,
Bests,
Andrea

Dear Andrea,

Thanks, very clear.
We have implemented this feature in the next FLUKA generation. As long as it’s not a showstopper for your simulations, would you mind waiting for the first v5 public release? It should be in the few years range.

Best regards,
Gabrielle

Dear Gabrielle,
thanks, I can definitively wait for the next major release!

Cheers
Andrea

1 Like