Integrating empirical equation into FLUKA routines

Hello FLUKA team,
I’ve been using comscw.f usreroutine to score LETd and RBEd by adding GETLET function and RBE formalisms and it works great. NOW, I do have an empirical equation based on the best fit to data which represents immunological response as a function of biological dose (dose+RBE). currently, I’m using MATLAB to process FLUKA output (3D USRBIN scores biological dose generated by using USRWEIGH card in FLAIR) and plug it in this empirical equation. My question here, I need to let FLUKA to score this quantity in voxel by voxel directly without using other means, how could I?
comscw.f is valid for dose weighted quantities which in this case is not helpful.

Thanks

Hello,

If I understand the question correctly, you want to score a function of the dose for each bin. Being the dose the average energy deposition per unit mass, this quantity should be not computed on the fly, but rather after the calculation of the dose.
If you could attach a working example of you current method, I could certainly help you better with some practical advices. I will include anyhow a general consideration about the way the scoring happens in FLUKA.

Assuming to work with a USRBIN scoring, during each cycle, FLUKA computes the energy deposition density, applies the correct weight (which for the dose is just 1/ρ, or whatever you want if you use comscw.f), and add it to the corresponding bin.
Then, via usbsuw.f, the results are postprocessed and merged. This means that, for each cycle, this routine compute the average of the requested quantity. It is also responsible for the computation of the statistical uncertainty of the estimator among different cycles.

In principle, still under the assumption that I understood your needs correctly, you are asking me if it possible to modify the postrprocessing procedure, including here your weighting function. Then, usbsuw.f would be the place to look at.
However, I would slightly dissuade you to proceed in such a way due to the inherent complexity and failure-proneness of that method. Indeed, post-processing with ad-hoc scripts would be much more advisable (faster) in my opinion.

Cheers,
Daniele

Thank you for responding,
I’ll give an example to clarify the problem:
If we consider scoring LETd, we could modify comscw.f routine, add GETLET function into it, and compile it to generate executable. In this scenario you need to score the dose using USRBIN card and by adding userweig card you would eventually score dose weighted LET in 3D mesh. This is a result of averaging over each voxel (fluence*LET)LET. Flair tools can postprocess this as you mentioned, but indeed I need to use other means to get LETd, this is not the issue though.
The above example directly weighed the LETv with the dose at the same voxel, but in my scenario is something else, I have a Gaussian function as follows:
immunological response = A+B * e(-((D-B)/C)^2) where ABC are constants.
D is dose
RBE is the biological dose which is already scored as dose weighted RBE the same way as LETv. I can do it in MATLAB by taking FLUKA output (dose weighted RBE) and plug it in this equation but I’m wondering if FLUKA can score it directly (if possible adding this empirical equation in any proper FLUKA’s user routines), it would make my life easier.

I hope that clarifies the issue
Thanks

Hello,

My only concern in applying such a procedure in Fluka is the following: let’s call d(i) the biological dose at the i-th step. Now, you can have fluctuation in this quantity, which eventually will average to the biological dose:

D = average( d(i) )

With FLUKA (particularly comscw.f), you can easily apply any weight during the simulation as you have correctly done. However, since you require the formula to be applied to the results (since you are satisfied with the Matlab results), the outcoming would be radically different.
Consider this example for instance.

Postprocessing: Effect = function( average( d(i) ) )
On the fly: Effect = average( function ( d(i) ) )

Unless the function is linear (i.e. a normalization factor), the two calculation will yield to different results.

Therefore, in the end, if you want to apply an empirical equation to the FLUKA output, and you want to do in in the fluka environment, you would need to do it in usbsuw.f.

As a side note, considering the work you are doing, I suggest you to read also the RAD-BIOL and TPPSCORE card documentation and the corresponding questions in the forum (i.e. About RAD-BIOL and TPSSCORE card, Question about results of RAD-BIOL and TPSSCORE card(s), Calculate with Biological model)

Hope it helped!
Cheers,
Daniele

Hi Daniele,
Sure you’re absolutely right, the result will be different because the required quantity is weighted by the dose. In FLUKA, I’m scoring the RBE (using RBE for double strand break as the endpoint(RBEdsb)) by implementing empirical equation based on the best fit to DNA damage simulation, and weight it with the dose using comscw.f routine (it’s not similar to what RAD-BIO card usually does). The FLUKA output (dose weighted RBEdsb) can now be easily plugged in the empirical equation using MATLAB, but I’m looking for a possible way to let FLUKA to do it without using MATLAB. That’s why I found comscw.f routine is not suitable to do the new task where the required effect is a function in both dose and RBEdsb (not weighting). The question remains, is there any way in FLUKA where I can use USRBIN to eventually score this quantity (i.e. immunological response = A+B * e(-((biological dose-B)/C)^2)

Thanks

Hello,

In this case, the only place you could act (at FLUKA level) is usbsuw.f.
In particular, during postprocessing, for each bin, this program collects the output data from the simulation (the one you’re correctly producing now), calculate the average value producing the requested statistics.
Let’s focus on line 301-304-305-376-377 to understand the behaviour of this program:

...
            READ (1) (GMSTOR(J), J = K0, K1)
...
               GDSTOR (J) = GDSTOR (J) + DBLE (GMSTOR(J))
     &                    * ( DBLE (NCASE) + DBLE(MCASE) * 1.D+09 )
...
         GDSTOR (J) = GDSTOR (J)
     &              / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )

This is what’s happening here: from the usrbin file you produced, you are reading the scored quantity from each bin, saving it in GMSTOR. Then, you accumulate the sum of each cycle with their (primary) weight in GDSTOR and finally you divide GDSTOR by the total (primary) particle weight to get the requested quantity.
The statistical error calculation works in a similar way, using the same data.

By default, when you process the output with flair, you are executing this program.
Since you want to plug in your function at a FLUKA level, I would suggest you to apply some small modification to this routine (copied in your working directory, to avoid any harms to the fluka installation one).

Immediately after you read the data (line 304), you could plug in your equation, i.e.:

...
            READ (1) (GMSTOR(J), J = K0, K1)
            IRECRD = IRECRD + 1
            DO J = K0, K1
               GMSTOR (J) = A + B * EXP 
     &                    ( -( GMSTOR (J) - B ) ** 2 / C )
               GDSTOR (J) = GDSTOR (J) + DBLE (GMSTOR(J))
     &                    * ( DBLE (NCASE) + DBLE(MCASE) * 1.D+09 )
               GBSTOR (J) = GBSTOR (J) + DBLE (GMSTOR(J))**2
     &                    * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
            END DO
...

This is just a snipped of how the implementation could look like. Then, you need to call the usbsuw.f via terminal, since flair will continue to call the default one. Notice that, since this process is detached from the fluka run, you don’t need to link this routine to the executable for the simulation.
Let me know if this solution could work for you.

Cheers,
Daniele

Thank you a lot for providing this solution, I’ve copied the routine, modified it, and saved it in the same simulation directory. Since I’m a beginner and only using flair to postprocess the data, Could you please provide the terminal’s command line that can call (the modified usbsuw.f), and return xxx.bnn file.

Thanks

Hello,

As a side note, before answering, notice that what flair is doing is calling the same program automatically for you. If you press on the bottom right corner, you can see the command executed:

screenshotForum

In your case, you need to compile the program, and call the executable. Then, via terminal, you just run the executable. It will ask you to provide the usrbin results as inputs and the name of the processed output.
Read also this discussion: Post-processing on the cluster

Good luck.
Cheers,
Daniele

Thank you very much, it works.