optical photons as primary particles

From: Oleg Tkach <tkatsh65_at_mail.ru>
Date: Wed, 03 Nov 2010 21:11:55 +0300

(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 Wed Nov 03 2010 - 21:55:14 CET

This archive was generated by hypermail 2.2.0 : Wed Nov 03 2010 - 21:55:21 CET