RE: [fluka-discuss]: What phase space file format can be read into source.f

From: Joachim Vollaire <joachim.vollaire_at_cern.ch>
Date: Sun, 19 Nov 2017 20:06:16 +0000

Hi

Did you compile the source routine and create a new fluka executable ? You should also have a SOURCE card in your input file to actually call the routine (otherwise FLUKA will just use the default value or whatever you have in the BEAM/BEAMPOS... cards .... and the logic in your routine will not be considered.)

Otherwise I had a very similar case a few years ago for which I was using the position of the beam impact on a beam intercepting device. In that case the beam energy was fixed, and I was reading the X, Y coordinate and the beam direction. You can find how the source routine was written for that case in presentation providedattachment. Note that this example date from about 10 years ago, so you should adapt it to the current source routine.

Also to write user routines in general, the material from the advanced course, contains very useful tips.

https://indico.cern.ch/event/489973/contributions/2000440/attachments/1272042/1972478/09_AdvancedUserRoutines2016.pdf

Hoping this help
Joachim Vollaire


-----Original Message-----
From: owner-fluka-discuss_at_mi.infn.it [mailto:owner-fluka-discuss_at_mi.infn.it] On Behalf Of Yu-Husan Chien
Sent: 19 November 2017 14:25
To: fluka-discuss_at_fluka.org
Subject: [fluka-discuss]: What phase space file format can be read into source.f

Dear FLUKA authors and users :

I want to simulate the proton therapy and now I have a phase space file of the proton gantry resulting from TOPAS. The data format transform into IAEA-phase space file format from binary file made by TOPAS using Matlab by myself .It includes 7 rows , just like :

         Xpos Ypos CosX CosY Energy Weight Particle_type

Each row have totally 8333128 columns data representing 8333128 particles recording in the file. And it's file type is .txt .

And now I use the OAUXFI routine :

         CALL OAUXFI('filename.txt', 99, 'OLD', STAT)

try to open them but I have some problem.

Is this .txt file can be read by source.f, I mean the data rearranged format by myself can be identified by source.f ?
And how should I do to tell fluka what is the value corresponding to each variable ?

My source.f is in the bottom of the mail, it can be compiled now , but can't to run.

Can anyone help me ? Thank!!!!


Y.H. Chien
National Tsing Hua University. Taiwan R.O.C





*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
       SUBROUTINE SOURCE ( NOMORE )

       INCLUDE '(DBLPRC)'
       INCLUDE '(DIMPAR)'
       INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
* *
* Copyright (C) 1990-2010 by Alfredo Ferrari & Paola Sala *
* All Rights Reserved. *
* *
* *
* New source for FLUKA9x-FLUKA20xy: *
* *
* Created on 07 January 1990 by Alfredo Ferrari & Paola Sala *
* Infn - Milan *
* *
* Last change on 17-Oct-10 by Alfredo Ferrari *
* *
* This is just an example of a possible user written source routine. *
* note that the beam card still has some meaning - in the scoring the *
* maximum momentum used in deciding the binning is taken from the *
* beam momentum. Other beam card parameters are obsolete. *
* *
* Output variables: *
* *
* Nomore = if > 0 the run will be terminated *
* *
*----------------------------------------------------------------------*
*
       INCLUDE '(BEAMCM)'
       INCLUDE '(FHEAVY)'
       INCLUDE '(FLKSTK)'
       INCLUDE '(IOIOCM)'
       INCLUDE '(LTCLCM)'
       INCLUDE '(PAPROP)'
       INCLUDE '(SOURCM)'
       INCLUDE '(SUMCOU)'
*
       LOGICAL LFIRST
*
       SAVE LFIRST
       DATA LFIRST / .TRUE. /

           PARAMETER (MAXNUM=8333128)
          DIMENSION PTYE(MAXNUM),XPOS(MAXNUM),YPOS(MAXNUM),
      & COSX(MAXNUM),COSY(MAXNUM),ENG(MAXNUM),
      & WEG(MAXNUM)
          real i
       SAVE PTYE,XPOS,YPOS,COSX,COSY,ENG,WEG

*======================================================================*
* *
* BASIC VERSION *
* *
*======================================================================*
       NOMORE = 0
* +-------------------------------------------------------------------*
* | First call initializations:
       IF ( LFIRST ) THEN
* | *** The following 3 cards are mandatory ***
          TKESUM = ZERZER
          LFIRST = .FALSE.
          LUSSRC = .TRUE.
* | *** User initialization ***
*
          CALL OAUXFI('mysource.txt', 99, 'OLD', STAT)
           READ(99,*) (XPOS(K),YPOS(K),COSX(K),COSY(K),ENG(K),WEG(K),
      & PTYE(K),K=1,MAXNUM)

*
*
          END IF
          i=INT((FLRNDM(xxx)*MAXNUM)+1)

*
* |
* +-------------------------------------------------------------------*
* Push one source particle to the stack. Note that you could as well
* push many but this way we reserve a maximum amount of space in the
* stack for the secondaries to be generated
* Npflka is the stack counter: of course any time source is called it
* must be =0
       NPFLKA = NPFLKA + 1
* Wt is the weight of the particle
       WTFLK (NPFLKA) = ONEONE
       WEIPRI = WEIPRI + WTFLK (NPFLKA)
* Particle type (1=proton.....). Ijbeam is the type set by the BEAM
* card
* +-------------------------------------------------------------------*
* | (Radioactive) isotope:
       IF ( IJBEAM .EQ. -2 .AND. LRDBEA ) THEN
          IARES = IPROA
          IZRES = IPROZ
          IISRES = IPROM
          CALL STISBM ( IARES, IZRES, IISRES )
          IJHION = IPROZ * 1000 + IPROA
          IJHION = IJHION * 100 + KXHEAV
          IONID = IJHION
          CALL DCDION ( IONID )
          CALL SETION ( IONID )
* |
* +-------------------------------------------------------------------*
* | Heavy ion:
       ELSE IF ( IJBEAM .EQ. -2 ) THEN
          IJHION = IPROZ * 1000 + IPROA
          IJHION = IJHION * 100 + KXHEAV
          IONID = IJHION
          CALL DCDION ( IONID )
          CALL SETION ( IONID )
          ILOFLK (NPFLKA) = IJHION
* | Flag this is prompt radiation
          LRADDC (NPFLKA) = .FALSE.
* | Group number for "low" energy neutrons, set to 0 anyway
          IGROUP (NPFLKA) = 0
* |
* +-------------------------------------------------------------------*
* | Normal hadron:
       ELSE
          IJBEAM=PTYE(i)
          IONID = IJBEAM
          ILOFLK (NPFLKA) = IJBEAM
* | Flag this is prompt radiation
          LRADDC (NPFLKA) = .FALSE.
* | Group number for "low" energy neutrons, set to 0 anyway
          IGROUP (NPFLKA) = 0
       END IF
* |
* +-------------------------------------------------------------------*
* From this point .....
* Particle generation (1 for primaries)
       LOFLK (NPFLKA) = 1
* User dependent flag:
       LOUSE (NPFLKA) = 0
* No channeling:
       LCHFLK (NPFLKA) = .FALSE.
       DCHFLK (NPFLKA) = ZERZER
* User dependent spare variables:
       DO 100 ISPR = 1, MKBMX1
          SPAREK (ISPR,NPFLKA) = ZERZER
  100 CONTINUE
* User dependent spare flags:
       DO 200 ISPR = 1, MKBMX2
          ISPARK (ISPR,NPFLKA) = 0
  200 CONTINUE
* Save the track number of the stack particle:
       ISPARK (MKBMX2,NPFLKA) = NPFLKA
       NPARMA = NPARMA + 1
       NUMPAR (NPFLKA) = NPARMA
       NEVENT (NPFLKA) = 0
       DFNEAR (NPFLKA) = +ZERZER
* ... to this point: don't change anything
* Particle age (s)
       AGESTK (NPFLKA) = +ZERZER
       AKNSHR (NPFLKA) = -TWOTWO
* Kinetic energy of the particle (GeV)
* TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
       TKEFLK (NPFLKA) = ENG(i)*EMVGEV
* Particle momentum
       PMOFLK (NPFLKA) = PBEAM
* PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
* & + TWOTWO * AM (IONID) ) )
* Cosines (tx,ty,tz)
       TXFLK (NPFLKA) = COSX(i)
       TYFLK (NPFLKA) = COSY(i)
       TZFLK (NPFLKA) = WBEAM
* TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
* & - TYFLK (NPFLKA)**2 )
* Polarization cosines:
       TXPOL (NPFLKA) = -TWOTWO
       TYPOL (NPFLKA) = +ZERZER
       TZPOL (NPFLKA) = +ZERZER
* Particle coordinates
       XFLK (NPFLKA) = XPOS(i)
       YFLK (NPFLKA) = YPOS(i)
       ZFLK (NPFLKA) = ZBEAM
* Calculate the total kinetic energy of the primaries: don't change
       IF ( ILOFLK (NPFLKA) .EQ. -2 .OR. ILOFLK (NPFLKA) .GT. 100000 )
      & THEN
          TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
       ELSE IF ( ILOFLK (NPFLKA) .NE. 0 ) THEN
          TKESUM = TKESUM + ( TKEFLK (NPFLKA) + AMDISC (ILOFLK(NPFLKA)) )
      & * WTFLK (NPFLKA)
       ELSE
          TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
       END IF
       RADDLY (NPFLKA) = ZERZER
* Here we ask for the region number of the hitting point.
* NREG (NPFLKA) = ...
* The following line makes the starting region search much more
* robust if particles are starting very close to a boundary:
       CALL GEOCRS ( TXFLK (NPFLKA), TYFLK (NPFLKA), TZFLK (NPFLKA) )
       CALL GEOREG ( XFLK (NPFLKA), YFLK (NPFLKA), ZFLK (NPFLKA),
      & NRGFLK(NPFLKA), IDISC )
* Do not change these cards:
       CALL GEOHSM ( NHSPNT (NPFLKA), 1, -11, MLATTC )
       NLATTC (NPFLKA) = MLATTC
       CMPATH (NPFLKA) = ZERZER
       CALL SOEVSV
       RETURN
*=== End of subroutine Source =========================================*
       END

__________________________________________________________________________
You can manage unsubscription from this mailing list at https://www.fluka.org/fluka.php?id=acc_info




__________________________________________________________________________
You can manage unsubscription from this mailing list at https://www.fluka.org/fluka.php?id=acc_info

Received on Sun Nov 19 2017 - 23:15:09 CET

This archive was generated by hypermail 2.3.0 : Sun Nov 19 2017 - 23:15:16 CET