From: Wang RuiGuang (wangrg@mail.ihep.ac.cn)
Date: Tue Oct 24 2006 - 10:46:58 CEST
Dear Fluka users and authors,
I am doing the simulation of secondaries of underground muons in water. Two
problems are met now. One is about option USRBIN, which of my input card is
shown as bellow,
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+....8
USRBIN 15.0 8.0 -51.0 500.0 500.0 900.0 Neutron
USRBIN -500.0 -500.0 0.0 10.0 10.0 900.0 &
USRBIN 15.0 208.0 -52.0 500.0 500.0 900.0 Edeposit
USRBIN -500.0 -500.0 0.0 10.0 10.0 900.0 &
The muon beam is symmetry with z axis
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+....8
BEAM 100.E+00 MUON-
BEAMPOS 0.0 0.0 -2000.0
When the number of X and Y bins is bellow 10, the running input file without
SOURCE is normal.
However, when the number of X and Y bins is increased, say 100, no resluts is
outputed with the same running except "No rantest2002 generated!
No test2001.err generated!
Saving additional files generated
No additional files to move
"
I want to known details of secondaries in every small steps. How to save this
problem?
Another problem is about using SOURCE. A option SOURCE is simply added in the
same input file.
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+....8
SOURCE
In SOURCE routine, first read data from a file, then calculate the beam
position and direction cosine, those values are loaded in the end. The SOURCE
routine was compiled and linked according to the instruction of page 302 of
FM.pdf. It is successful. When running Fluka with 'rfluka -e myfluka
inputfile', it is stop and showed
"../flutil/rfluka: line 311: 8408 Aborted (core dumped)
${EXE} <$INPF 2>$LOGF >$LOGF"
I don't know the reasons, hope to get your help!
Thanks and Best regards!
R.G. Wang
Modified SOURCE routine is as bellow.
*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
SUBROUTINE SOURCE ( NOMORE )
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
* *
* Copyright (C) 1990-2005 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 17-jun-05 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. /
******************************************************
*** For reading data from DYB muongenerator
real PI, TWO_PI, HPI,rcut
parameter(PI=3.1415926, TWO_PI=PI*2.0,HPI=PI/2.0)
integer N,ii
parameter(N=1e6)
Real emuon(N),theta(N),phi(N),a,b,c
real x0,y0,z0, emu0,theta0,phi0
real DUMMY,arandom(3)
real area_cross,sx,sy,sz,sx_p,sy_p,sinphi,cosphi
real lRock,wRock,hRock
real beamx,beamy,beamz
********************************************************
*======================================================================*
* *
* BASIC VERSION *
* *
*======================================================================*
NOMORE = 0
* +-------------------------------------------------------------------*
* | First call initializations:
IF ( LFIRST ) THEN
* | *** The following 3 cards are mandatory ***
TKESUM = ZERZER
LFIRST = .FALSE.
LUSSRC = .TRUE.
* | *** User initialization ***
* |
***********************************************************************
*---------------------------------------------------------------------*
* Reading data from a muongenerator file. *
* a,b and c is unuseful variables *
*---------------------------------------------------------------------*******************************************************************
open(21,file='mountaindyw.near',status='old')
do ii=1,N
read(21,*) emuon(ii),theta(ii),phi(ii),a,b,c
if(theta(ii).ge.90.)then
print*,emuon(ii),theta(ii),phi(ii),a,b,c
endif
theta(ii)=theta(ii)*PI/180.0 ! input angle is in degree
phi(ii)=phi(ii)*PI/180.0
enddo
close(21)
rcut=200.0 !cm
lRock=800
wRock=800
hRock=1210
**** return message from first call
write(LUNOUT,*) 'version 1 of sourcedyb1.f called'
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
* Lstack 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
* +-------------------------------------------------------------------*
* | Heavy ion:
IF ( IJBEAM .EQ. -2 ) THEN
IJHION = IPROZ * 1000 + IPROA
IJHION = IJHION * 100 + KXHEAV
IONID = IJHION
CALL DCDION ( IONID )
CALL SETION ( IONID )
ILOFLK (NPFLKA) = IJHION
* |
* +-------------------------------------------------------------------*
* | Normal hadron:
ELSE
IONID = IJBEAM
ILOFLK (NPFLKA) = IJBEAM
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
******************************************************************
ii=1+int(FLRNDM(DUMMY)*(N-1))
theta0=theta(ii) ! muon direction changed with different
coordinate system
if(phi(ii).le.HPI) then
phi0=HPI-phi(ii) ! phi<90 , phi+=90 - phi
else
phi0=5.0*HPI-phi(ii) ! phi>90 , phi+=450 - phi
endif
emu0=emuon(ii)
beamx=sin(theta0)*cos(phi0) !direction cosine of the beam with respect
to the x-axis
beamy=sin(theta0)*sin(phi0) !direction cosine of the beam with respect
to the y-axis
beamz= sqrt(1-beamx**2-beamy**2)
sx=4*(wRock+rcut)*(hRock+rcut) ! x direction plane area of the outer
cuboid.
sy=4*(hRock+rcut)*(lRock+rcut) ! y
sz=4*(lRock+rcut)*(wRock+rcut) ! z
sx=sx*abs(sin(theta0)*cos(phi0)) ! (theta,phi) direction projection area
sy=sy*abs(sin(theta0)*sin(phi0)) ! y
sz=sz*abs(cos(theta0)) ! z
area_cross=sx+sy+sz
sx_p=sx/area_cross ! the probability of muon on this
plane: Sx
sy_p=sy/area_cross+sx_p
arandom(1)=FLRNDM(DUMMY)
arandom(2)=FLRNDM(DUMMY)
arandom(3)=FLRNDM(DUMMY)
cosphi=cos(phi0)
sinphi=sin(phi0)
if(arandom(3).lt.sx_p) then
if(sx_p.ne.0.0) then
x0=(lRock+rcut)*(cosphi/abs(cosphi))*(-1.0) ! muon vertex on plane
in x direction
y0=(wRock+rcut)*(2.0*arandom(1)-1.0)
z0=(hRock+rcut)*(2.0*arandom(2)-1.0)
endif
else if(arandom(3).lt.sy_p) then
if(sy_p.ne.0.0) then
x0=(lRock+rcut)*(2.0*arandom(1)-1.0) ! muon vertex on plane
in y direction
y0=(wRock+rcut)*(sinphi/abs(sinphi))*(-1.0)
z0=(hRock+rcut)*(2.0*arandom(2)-1.0)
endif
else
x0=(lRock+rcut)*(2.0*arandom(1)-1.0) ! muon vertex on the top plane
z direction
y0=(wRock+rcut)*(2.0*arandom(2)-1.0)
z0= hRock+rcut
endif
z0 = z0-290 ! transform from one coordinate to another
******************************************************************
* Particle age (s)
AGESTK (NPFLKA) = +ZERZER
AKNSHR (NPFLKA) = -TWOTWO
* Group number for "low" energy neutrons, set to 0 anyway
IGROUP (NPFLKA) = 0
* Kinetic energy of the particle (GeV)
PBEAM = emu0
TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
* Particle momentum
PMOFLK (NPFLKA) = PBEAM
* PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
* & + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
* Cosines (tx,ty,tz)
TXFLK (NPFLKA) = beamx
TYFLK (NPFLKA) = beamy
TZFLK (NPFLKA) = beamz
c TXFLK (NPFLKA) = UBEAM
c TYFLK (NPFLKA) = VBEAM
c 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) = x0
YFLK (NPFLKA) = y0
ZFLK (NPFLKA) = z0
c XFLK (NPFLKA) = XBEAM
c YFLK (NPFLKA) = YBEAM
c 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
* Flag this is prompt radiation
LRADDC (NPFLKA) = .FALSE.
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
This archive was generated by hypermail 2.1.6 : Wed Oct 25 2006 - 09:38:23 CEST