Re: optical photons as primary particles

From: Giuseppe Battistoni <Giuseppe.Battistoni_at_mi.infn.it>
Date: Thu, 04 Nov 2010 09:38:11 +0100

Dear Oleg,
unfortunately in the manual it is written at 12.2) that "...light must
be generated by a user-written SOURCE routine,
   or starting from ordinary particles".
This induces people to think that optical photons can be generated in
the standard source.f user routine, but this is not true.
That routine fills the stack of ordinary particles, while optical
photons must be loaded into a separate stack.
The propagation of optical photons is driven by a different set of
routines of the core code,
the head of them being KASOPH (while for ordinary heavy particles is
KASKAD).
A more generalized source routine will be available in a next version.
In the meanwhile, if you wish, you can test my private version of a
non-standard source routine
designed to generate optical photons (IJBEAM=-1), hoping that you can
adapt it to your problem.
Notice that the stack variables have different names and that this
source makes use of
variables set in BEAM and BEAMPOS cards to assign initial position and
direction
while energy is set to a given fixed value that has no special meaning
(it was useful for my test case..).
I have used it once and I cannot assure 100% that is bug-free, and no
assistance is guaranteed.

      Good luck
              Giuseppe

On 11/3/2010 7:11 PM, Oleg Tkach wrote:
> (smtp1.mi.infn.it [192.84.138.69]); Wed, 03 Nov 2010 19:12:01 +0100 (CET)
> Sender: owner-fluka-discuss_at_mi.infn.it
>
> Dear FLUKA authors,
>
> I rewrote, compiled and linked a user-written source.f routine to
> define optical photons as primary particles (OPTIPHOT particle code -1)
> and transport them via the BEAM and SOURCE cards:
>
> BEAM -2.0E-2
> -1.OPTIPHOT
> SOURCE 1.05E-04
>
> And my modified source.f:
>
> *$ CREATE SOURCE.FOR
> *COPY SOURCE
> *
> *=== source ===========================================================*
> *
> SUBROUTINE SOURCE ( NOMORE )
>
> INCLUDE '(DBLPRC)'
> INCLUDE '(DIMPAR)'
> INCLUDE '(IOUNIT)'
> *
> *----------------------------------------------------------------------*
> * *
> * Copyright (C) 1990-2006 by Alfredo Ferrari& Paola Sala *
> * All Rights Reserved. *
> * *
> * *
> * New source for FLUKA9x-FLUKA200x: *
> * *
> * Created on 07 january 1990 by Alfredo Ferrari& Paola Sala *
> * Infn - Milan *
> * *
> * Last change on 03-mar-06 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. *
> * *
> *----------------------------------------------------------------------*
> *
> INCLUDE '(BEAMCM)'
> INCLUDE '(FHEAVY)'
> INCLUDE '(FLKSTK)'
> INCLUDE '(IOIOCM)'
> INCLUDE '(LTCLCM)'
> INCLUDE '(PAPROP)'
> INCLUDE '(SOURCM)'
> INCLUDE '(SUMCOU)'
> *
> LOGICAL LFIRST
> *
> SAVE LFIRST
> DATA LFIRST / .TRUE. /
> *======================================================================*
> * *
> * BASIC VERSION *
> * *
> *======================================================================*
> NOMORE = 0
> * +-------------------------------------------------------------------*
> * | First call initializations:
> IF ( LFIRST ) THEN
> * | *** The following 3 cards are mandatory ***
> TKESUM = ZERZER
> LFIRST = .FALSE.
> LUSSRC = .TRUE.
> * | *** User initialization ***
> write(lunout,*)
> write(lunout,*) 'modified source.f for OPTIPHOT by Tkach called'
> write(lunout,*) ' IJBEAM = ', IJBEAM, ' WHASOU(1) = ', WHASOU(1)
> END IF
> * ..............
> *
> * |
> * +-------------------------------------------------------------------*
> * 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.
> * |
> * +-------------------------------------------------------------------*
> * | Normal hadron:
> ELSE
> IONID = IJBEAM
> ILOFLK (NPFLKA) = IJBEAM
> * | Flag this is prompt radiation
> LRADDC (NPFLKA) = .FALSE.
> END IF
> * |
> * +-------------------------------------------------------------------*
> * From this point .....
> * Particle generation (1 for primaries)
> LOFLK (NPFLKA) = 1
> * User dependent flag:
> LOUSE (NPFLKA) = 0
> * 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
> * Group number for "low" energy neutrons, set to 0 anyway
> IGROUP (NPFLKA) = 0
> * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%--from Tkach-*
> * |
> IF ( IJBEAM .EQ. -1 ) THEN
> * Kinetic energy of the Optical photon (GeV)
> WV = WHASOU(1)
> ENERGY = 4.14D-24 * CLIGHT / WV
> TKEFLK (NPFLKA) = ENERGY
> * Particle momentum
> PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
> & + TWOTWO * AM (IONID) ) )
> write(lunout,*) 'IJBEAM = ', IJBEAM
> write(lunout,*) 'IONID =',IONID,'ILOFLK() =',ILOFLK (NPFLKA)
> write(lunout,*) 'WV = ', WV, ' WHASOU(1) = ', WHASOU(1)
> write(lunout,*) 'ENERGY = ',ENERGY,'TKEFLK() = ',TKEFLK(NPFLKA)
> ELSE
> * Kinetic energy of the particle (GeV)
> TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
> * Particle momentum
> PMOFLK (NPFLKA) = PBEAM
> * PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
> *& + TWOTWO * AM (IONID) ) )
> END IF
> * |
> * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%--to Tkach-*
> * Cosines (tx,ty,tz)
> TXFLK (NPFLKA) = UBEAM
> TYFLK (NPFLKA) = VBEAM
> 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) = XBEAM
> YFLK (NPFLKA) = YBEAM
> 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
>
> But the run was aborted (core dumped) and in the out file I read:
>
> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
> !!!!!!!!!!!!!!!!!!!!!!!!!!!! FULLY ANALOGUE RUN !!!!!!!!!!!!!!!!!!!!!!!!!!!!
> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>
> Total time used for initialization: 3.27 s
> +++++++++++++++++++++++++++++++++
>
> NEXT SEEDS: 0 0 0 0 0 0 181CD
> 3039 0 0
>
> modified source.f for OPTIPHOT by Tkach called
> IJBEAM = -1 WHASOU(1) = 0.000105
> IJBEAM = -1
> IONID = -1ILOFLK() = -1
> WV = 0.000105 WHASOU(1) = 0.000105
> ENERGY = 1.18203883E-09TKEFLK() = 1.18203883E-09
> Abort called from KASKAD reason ??? WHY HERE ??? Run stopped!
> STOP ??? WHY HERE ???
>
> Where is KASOPH routine? What am I doing wrong ?
> Could you please help me to solve these problems?
> Thanks in advance,
> Oleg
>
>
>
>

Received on Thu Nov 04 2010 - 10:14:07 CET

This archive was generated by hypermail 2.2.0 : Thu Nov 04 2010 - 10:14:08 CET