---
*******************************************************************
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 20-12-2018 3:47 pm, Francesc Salvat-Pujol escribió:
> 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:56:54 CET
This archive was generated by hypermail 2.3.0 : Thu Dec 20 2018 - 20:56:55 CET