Re: [fluka-discuss]: Beam definition

From: Francesc Salvat-Pujol <francesc.salvat.pujol_at_cern.ch>
Date: Thu, 20 Dec 2018 19:47:17 +0100

Hey Jilberto,

It could be that your linux distribution is using some interesting
flavor of awk. Just for curiosity: are you running Ubuntu? If so, which
version? Could you paste the output of

  ls -l $(which awk)
  ls -l /usr/bin/*awk*
  ls -l /bin/*awk*

Are you using the latest FLUKA tarball (i.e. downloaded after Nov 22
2018 from the FLUKA website)? If not, update. In the not-too-distant
past we spotted a similar situation and we protected the scripts
accordingly. Should the problem persist, let us know (and if you are
linking the gfortran version of FLUKA paste also the output of
"gfortran --version").

Season greetings,

Cesc

On Thu, Dec 20 2018, at 13:55 -0300, Jilberto Antonio Zamora Saá wrote:
>
>Dear Francesc
>
>I have tried to run examples that I found in previous fluka courses,
>however when I tryed to compile the source I got the following
>message:
>
>sirjazs_at_LENOVO:~/fluka-work/Beam-Test/ex08$ $FLUPRO/flutil/fff source.f
>awk: line 0: regular expression compile failed (missing '(')
>)
>/home/sirjazs/FLUKA/flutil/fff: 105: [: -le: unexpected operator
>/home/sirjazs/FLUKA/flutil/fff: 119: /home/sirjazs/FLUKA/flutil/fff:
>[[: not found
>sirjazs_at_LENOVO:~/fluka-work/Beam-Test/ex08$
>
>
>the files that I'm using are attached :-)
>
>
>Do you have any idea about what it is happening?
>
>
>saludos,
>Jilberto
>
>
>
>
>---
>*******************************************************************
> Dr. Jilberto Zamora Saa
> Dzhelepov Laboratory of Nuclear Problems
> Joint Institute for Nuclear Research
> 141980 Dubna, Moscow region, Russia
> phone: +(496)2162589 (Office)
> phone: +7(915)4671294 (Mobile)
> email: jzamorasaa_at_jinr.ru
>
>********************************************************************
>
>El 09-11-2018 5:48 pm, Francesc Salvat-Pujol escribió:
>>Hola Jilberto,
>>
>>For inspiration you may first want to examine a slightly simpler 1D
>>case
>>(sampling from a tabulated energy spectrum), as shown e.g. in
>>
>> http://www.fluka.org/web_archive/earchive/new-fluka-discuss/6798.html
>>
>>which contains an example spectrum and a custom source.f routine that
>>samples from it. Note also the modification of the input file mentioned
>>in the link. If you make a diff (or e.g. vim -d) with the distributed
>>source.f (under usermvax/ in your FLUKA directory) you will see the
>>modifications. In a nutshell, there is:
>>
>> - a block that loads a source spectrum from the provided file at
>> initialization and prepares a cumulative function for sampling
>> later,
>>
>> - a block which samples the kinetic energy from the histogram and
>>
>> - a block which carefully sets the sampled kinetic energy and
>> calculates the modulus of the linear momentum (of course with the
>> relativistic expression).
>>
>>You then compile your custom source.f
>>
>> $FLUPRO/flutil/fff source.f
>>
>>If there are no errors, this produces an object file which you then
>>bundle into your custom binary.
>>
>> $FLUPRO/flutil/lfluka -m fluka -o tuejecutable source.o
>>
>>Depending on what you do, link with ldpmqmd (see the manual). Run as
>>usual, but passing "-e tuejecutable" additionally.
>>
>>If I follow you correctly, you need instead to sample from a tabulated
>>2D distribution h(E,theta), given in terms of the energy E and a
>>(polar?) angle theta. So it's a slight generalization of the foregoing.
>>Stay alert for the integrals and the sampling of the polar angle (see
>>below).
>>
>>I assume that the theta dependency is not trivial, i.e., that you have
>>sufficiently different behaviors in theta at different energies
>>(otherwise the stuff below can be greatly simplified...).
>>
>>Presumably you have h(E,theta) tabulated on a reasonably dense grid of
>>energies E_i and angles theta_j. On paper, at each tabular energy you
>>would do an integration on the sphere (azimuthal angle aside), à la
>>
>> g(E_i) = \int_0^pi d theta sin(theta) h(E_i,theta) ,
>>
>>so that g(E) gives you an energy distribution regardless of angle
>>(treated below). In practice you evaluate it numerically in as
>>reasonable a way as you can, but in any event not forgetting about the
>>solid-angle element. In the same spirit as in the 1D case above, you
>>may
>>now use the tabulated g(E_i) to generate a cumulative distribution for
>>the sampling of the kinetic energy.
>>
>>Once you have a sampled kinetic energy Esampled, you look for the
>>"active" tabular interval verifying E_i <= Esampled < E_{i+1}. Use a
>>homogeneously distributed random number (search in the manual for
>>FLRNDM) to take E_i with probability (E_{i+1}-Esampled)/(E_{i+1}-E_i),
>>and E_{i+1} otherwise (I do not write down exactly how because every
>>time I get it wrong, but the idea is along these lines...).
>>
>>At each tabular energy E_i, you can keep the cumulative angular
>>distribution in memory (do not forget the solid-angle element when
>>integrating...) and use it to sample the polar angle in exactly the
>>same
>>spirit as above. It then remains to sample the azimuthal angle. Unless
>>you have something else in mind, you can sample e.g. homogeneously in
>>[0,2pi). With this you have all you need to set the initial
>>direction in
>>your source.f in terms of the director cosines using the variables
>>TXFLK, TYFLK, TZFLK.
>>
>> ===============================================
>> *************** WARNING... ********************
>> ===============================================
>>
>>Things can go wrong when setting up schemes like the above by hand...
>>Debug and test intensively before running any simulation. Make e.g.
>>a 2D
>>histogram of the sampled energies and angles and convince yourself that
>>you are indeed sampling what was intended.
>>
>>As a final note, I do not know what degrees of freedom you have in
>>preparing the input 2D source histogram, but any time there's a polar
>>angle theta involved, a way to sleep well is to tabulate, integrate,
>>and
>>sample in terms of mu=cos(theta). Doing the change of variables you
>>immediately see e.g.
>>
>> g(E_i) = \int_{-1}^1 d mu h(E_i,mu) ,
>>
>>that is, the term "sin(theta) d theta" in the solid-angle element
>>becomes trivially "d mu" and you can mostly forget about faux pas :)
>>
>>Hope this is reasonably helpful/accurate!
>>
>>Un saludo,
>>
>>Cesc
>>
>>PS: scoring-wise, there should be nothing special.
>>
>>On Fri, Nov 09 2018, at 12:17 -0300, jilberto Zamora Saa wrote:
>>>
>>>Dear FLUKA experts,
>>>
>>>How I should do in case I want to define my own beam of particles,
>>>let say, I would like to provide a spectrum of muons which depend
>>>on energy and Zenit angle and then use it to see the fluence in a
>>>detector.
>>>
>>>any help is welcome
>>>
>>>regards,
>>>Jilberto
>>
>>--
>>Francesc Salvat Pujol
>>CERN-EN/STI
>>CH-1211 Geneva 23
>>Switzerland
>>Tel: +41 22 76 64011
>>Fax: +41 22 76 69474

>* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7...
>TITLE
>FLUKA Course Exercise
>*
>* use names everywhere and free format for geometry
>DEFAULTS NEW-DEFA
>*
>* beam definitions
>BEAM -3.5 -0.082425 -1.7 0.0 0.0 1.0PROTON
>BEAMPOS 0.0 0.0 -0.1 0.0 0.0
>GEOBEGIN COMBNAME
> 0 0 Cylindrical Target
>*
>* Bodies
>* ------
>*
>* General definitions:
>* blackhole to include geometry
>SPH BLK 0.0 0.0 0.0 10000.0
>* void
>RPP VOI -1000.0 1000.0 -1000.0 1000.0 -1000.0 1000.0
>*
>* Lead target:
>* cylinder to contain target
>ZCC TARG 0.0 0.0 5.0
>* planes limiting the target
>XYP ZTlow 0.0
>XYP ZThigh 10.0
>* planes segmenting the target
>XYP T1seg 1.0
>XYP T2seg 2.0
>END
>*
>* Regions
>* -------
>*
>* Blackhole
>BLKHOLE 5 +BLK -VOI
>*
>* Target segment 1
>TARGS1 5 +TARG -ZTlow +T1seg
>* Target segment 2
>TARGS2 5 +TARG -T1seg +T2seg
>* Target segment 3
>TARGS3 5 +TARG -T2seg +ZThigh
>*
>* Air around target
>INAIR 5 +VOI -( +TARG -ZTlow +ZThigh )
>END
>GEOEND
>#if 0
>* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7...
>* switch on to debug this geometry
>GEOEND 6.0 0.0 11.0 -6.0 0.0 -6.0DEBUG
>GEOEND 120.0 1.0 170.0 &
>#endif
>*
>* material definitions
>MATERIAL 1.0 WATER
>COMPOUND 2.0 HYDROGEN 1.0 OXYGEN WATER
>MATERIAL 0.001225 AIR
>COMPOUND -0.9256 NITROGEN -0.2837 OXYGEN -0.01572 ARGONAIR
>ASSIGNMA BLCKHOLE BLKHOLE
>ASSIGNMA WATER TARGS1
>ASSIGNMA ALUMINUM TARGS2
>ASSIGNMA LEAD TARGS3
>ASSIGNMA AIR INAIR
>*
>* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7...
>*
>* scoring
>*
>* target: energy deposition and fluence
>USRBIN 11.0 ENERGY -40.0 10.0 0.0 15.0TargEne
>USRBIN 0.0 0.0 -5.0 100.0 1.0 200.0&
>USRBIN 11.0 HAD-CHAR -40.0 10.0 0.0 15.0TargChH
>USRBIN 0.0 0.0 -5.0 100.0 1.0 200.0&
>USRBIN 11.0 NEUTRON -40.0 10.0 0.0 15.0TargN
>USRBIN 0.0 0.0 -5.0 100.0 1.0 200.0&
>*
>* charged hadron fluence at boundaries between target segments
>USRBDX 99.0 HAD-CHAR -50.0 TARGS1 TARGS2 78.5398Sp1ChH
>USRBDX 10.0 0.001 40.0 &
>USRBDX 99.0 HAD-CHAR -50.0 TARGS2 TARGS3 78.5398Sp2ChH
>USRBDX 10.0 0.001 40.0 &
>* charged hadron fluence exiting lead target
>USRBDX 99.0 HAD-CHAR -50.0 TARGS3 INAIR 329.87Sp3ChH
>USRBDX 10.0 0.001 40.0 &
>* double-differential charged hadron fluence entering lead target
>USRBDX 99.0 HAD-CHAR -54.0 TARGS2 TARGS3 78.5398Sp2ChHA
>USRBDX 10.0 0.001 40.0 3.0&
>* charged hadron fluence in lead target
>USRTRACK -1.0 HAD-CHAR -55.0 TARGS3 628.31 40.0TrChH
>USRTRACK 10.0 0.001 &
>*
>* charged pion angular distribution exiting lead target
>USRYIELD 124.0 PIONS+- -57.0 TARGS3 INAIR 1.0YieAng
>USRYIELD 180.0 0.0 18.0 10.0 0.0 3.0&
>*
>* residual nuclei in lead target
>RESNUCLE 3.0 -60.0 TARGS3 activ
>*
>* Exercise: scoring
>* -----------------
>*
>* electron spectra, see difference between current and fluence scoring
>USRBDX 99.0 E+&E- -51.0 TARGS2 TARGS3 78.5398Sp2El
>USRBDX 10.0 0.001 40.0 &
>USRBDX -1.0 E+&E- -52.0 TARGS2 TARGS3 78.5398Sp2ElC
>USRBDX 10.0 0.001 40.0 &
>*
>* energy deposition by region, check with values in the output
>USRBIN 2.0 ENERGY 41.0 TARGS3 TbyReg
>USRBIN TARGS1 1.0 &
>*
>* energy deposited by electrons only
>USRBIN 2.0 ENERGY 42. TARGS3 TbyRegE
>USRBIN TARGS1 1.0 &
>AUXSCORE USRBIN ELECTRON TbyRegE
>RANDOMIZ 1.0
>START 1000.0 0.0
>STOP

>*$ 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
>*
>c defining and saving spectrum arrays
> DIMENSION ENEPOI(0:1000),ENEPRO(0:1000),ENECUM(0:1000)
> SAVE ENEPOI, ENEPRO, ENECUM
>c saving spectrum dimension
> SAVE IMAX
>*
> 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 ***
> 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
>c sampling uniformly between 30 and 70 MeV
> XYZ=FLRNDM(XYZ)
> ESAMPLE=3D-2+XYZ*(4D-2)
>* Kinetic energy of the particle (GeV)
> TKEFLK (NPFLKA) = ESAMPLE
>* Particle momentum
> PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
> & + TWOTWO * AM (IONID) ) )
>* 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


--
Francesc Salvat Pujol
CERN-EN/STI
CH-1211 Geneva 23
Switzerland
Tel: +41 22 76 64011
Fax: +41 22 76 69474
__________________________________________________________________________
You can manage unsubscription from this mailing list at https://www.fluka.org/fluka.php?id=acc_info
Received on Thu Dec 20 2018 - 20:55:43 CET

This archive was generated by hypermail 2.3.0 : Thu Dec 20 2018 - 20:55:45 CET