Transcript FLUKA Advanced Biasing

FLUKA
Advanced Biasing
FLUKA Users Meeting
CERN 31/05/2007
Outline

Sampling




FLUKA User biasing routines




Random numbers
Sampling from discrete & generic distributions
Sampling from biased distributions
USIMBS
UDCDRL
SOURCE
Application: n_TOF biasing example
CERN FLUKA Users Meeting
2
Pseudo-Random number





All computers provide a pseudo random number generator.
The FLUKA PRNG is provided by the function FLRNDM(DUMMY).
A 64-bit random number generator returning a double precision
number.
All such generators return a sequence of numbers uniformly
distributed between [0,1). To avoid overflows sometimes, you
need to change to (0,1] then use 1.0-FLRNDM(x)
Uniform and Random are independent concepts. For MC
sometimes uniformity is more important than randomness, but
not if correlations are important
It is important also that sampled histories are uncorrelated with
each other
CERN FLUKA Users Meeting
3
FLUKA Random number generators





Function FLRNDM(DUMMY)
The main FLUKA PRNG. Returns a double precision random
number uniformly distributed in the [0,1) interval
Subroutine SFECFE(SFE, CFE)
returns a random azimuthal sine and cosine
Subroutine RACO(TXX, TYY, TZZ)
returns the cosines of a randomly oriented vector
Subroutine SFLOOD(XXX, YYY, ZZZ, UXX, VXX, WXX)
returns a random position and direction on the surface of a
sphere of unit radius (generating a uniform and isotropic field
inside the sphere)
Subroutine FLNRRN(RGAUSS)
returns a random number normally distributed with unit variance
CERN FLUKA Users Meeting
4
Sampling from a discrete distribution



Suppose to have a discrete random variable x, that can assume
values x1,x2,…,xn with probability p1,p2,…,pn
Assume Sipi = 1, or normalize it
Divide the interval [0,1) in n subintervals with limits
y0=0, y1=p1, y2=p1+p2, …, yn=Sipi

Generate a uniform random number x

Find in which of the i y-intervals it is
WARNING: This can be a time consuming process:


Use optimized search techniques like Binary search (look n_TOF
example)
Use of equal-probability bins (not recommended when tails are
present)

Select the corresponding xi as sampled value

Since
x is uniformly random, we have
P(xi) = P(yi-1 < x < yi) = yi – yi-1 = pi
CERN FLUKA Users Meeting
5
Sampling from a Generic Distribution
Using one random number
 Integrate the distribution function, analytically or numerically, and
normalize to 1 to obtain the normalized cumulative distribution
x
f ( x) dx

F ( x) 
f ( x) dx

Sample a uniform pseudo-random number x
xmin
xmax

xmin

Get the desired result by finding the inverse value x=F-1(x),
analytically or most often by interpolation (table lookup)

Since
x is uniformly random, we have
b
P(a  x  b)  P[ F (a)  x  F (b)]  F (b)  F (a)   f ( x)dx
a
CERN FLUKA Users Meeting
6
Example
t = - λ log (1 – ξ)
Practical rule: a distribution can be sampled directly if and only if
its pdf can be integrated and the integral inverted
CERN FLUKA Users Meeting
7
Sampling from a Generic Distribution
Using two random numbers (Rejection technique)
 Choose a constant value C > f(x) for any x
 Sample two random numbers x1 and x2
 (x1, Cx2) represents a point on the 2D space
 If x2<f(x1)/C then X=x1
otherwise re-sample x1, x2


The probability that a x1 value is
accepted is f(x1)/C for any x1
P(x)dx = P(x1=x)dx  f(x1=x)/C
Since P(x1)dx = const,
P(x) = const  f(x)
CERN FLUKA Users Meeting
8
Sampling from a biased distribution
CERN FLUKA Users Meeting
9
Example
CERN FLUKA Users Meeting
10
Uniform sampling from a steep distribution
CERN FLUKA Users Meeting
11
User written biasing
FLUKA offers the following routines for user-written biasing
 ubsset.f: User BiaSing SETting



usimbs.f: USer defined IMportance BiaSing



called after reading in the input file and before the first event
allows to alter almost any biasing weight on a region-dependent basis
if activated, called at every particle step
allows to implement any importance biasing scheme based on region number
and/or phase space coordinates
udcdrl.f: User defined DeCay DiRection biasing and Lambda

only for neutrinos emitted in decays: bias the direction of emitted neutrino
Not biasing by itself, but it could be used for generating biased runs
 source.f: User written source

to sample primary particle properties from distribution in space, energy,
time…
CERN FLUKA Users Meeting
12
User Defined Importance Biasing


Typical problem: Spend a lot of time to write the problem input,
geometry and on the first runs you realize that statistics are not
good
First method (and safest) is to introduce region importance
biasing. In FLUKA you can introduce it with two ways:



1st Manually slice the geometry and increase the number of regions.
Modifying an existing geometry to introduce biasing can be a very
cumbersome process
2nd Introduce the “importance biasing” information with a user
fortran routine independent of the regions defined in the geometry
Routine: usimbs.f
USer defined IMportance BiaSing
Allows to implement any importance biasing scheme based on
region number and/or phase space coordinates
CERN FLUKA Users Meeting
13
User Defined Importance Biasing

Enable the call to USIMBS routine with the BIASING card:






WHAT(1) Particles to be biased
WHAT(2) and WHAT(3) ≠ 1.0
WHAT(4) Lower bound of region
WHAT(5) Upper bound of region
WHAT(6) Step
SDUM = USER
(Any value)
CERN FLUKA Users Meeting
14
usimbs.f – Routine
The routine is called on every particle
step!
WARNING: can slow down the
execution speed! Use with caution.
Input:
 Region information at the beginning
and end of the step
 X,Y,Z coordinates through the
TRACKR common.
Beginning:
X, Y, Ztrack(0)
End:
X, Y, Ztrack(Ntrack)
 Particle type and Energy could be
used for even more advanced
biasing schemes
Output:
The routine must return the
importance ratio to the new position
(end/beginning) in the variable FIMP
SUBROUTINE USIMBS ( MREG, NEWREG, FIMP )
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
*
*---------------------------------------------------------------------- *
*
*
*
USer defined IMportance BiaSing:
*
*
*
*
Created on
02 july 2001
by
Alfredo Ferrari & Paola Sala *
*
Infn - Milan
*
*
*
*
Last change on 09-jul-01
by
Alfredo Ferrari
*
*
*
*
Input variables:
*
*
Mreg = region at the beginning of the step
*
*
Newreg = region at the end of the step
*
*
(thru common TRACKR):
*
*
Jtrack = particle id. (Paprop numbering)
*
*
Etrack = particle total energy (GeV)
*
*
X,Y,Ztrack(0) = position at the beginning of the step
*
* X,Y,Ztrack(Ntrack) = position at the end of the step
*
*
*
*
Output variable:
*
*
Fimp = importance ratio (new position/original one)
*
*
*
*---------------------------------------------------------------------- *
e
INCLUDE '(TRACKR)'
*
FIMP
= ONEONE
RETURN
END
CERN FLUKA Users Meeting
15
usimbs.f – Important Notes



The routine has only a relative effect on the weight of the particle
 Beam particles can have any weight.
Importance ratio will be limited by WW-THRESh card
The Russian Roulette / Splitting will take place at the end of the
step, and not on a fixed region boundary.
The biasing position will be a little fuzzy depending on the particle step.
This has a visible effect when it is applied to low density materials (i.e.
air)



Results will similar but not be the same as with the manual
region biasing.
Is a great time saver for complex geometries, as well different
biasing schemes
Combined with the particle type and energy could become very
powerful
CERN FLUKA Users Meeting
16
usimbs.f: Simple Example

Biasing a factor of 2 every 50 cm on the Z direction from 100cm
to 500cm
Initial position
ZSTART = ZTRACK(0)
IF (ZSTART .LT. 100.0D0) THEN
FSTART = ONEONE
ELSE IF (ZSTART .GT. 500.0D0) THEN
FSTART = TWOTWO ** NINT((500.0D0-100.0D0)/50.0D0)
ELSE
FSTART = TWOTWO ** NINT((ZSTART-100.0D0)/50.0D0)
ENDIF
ZEND = ZTRACK(NTRACK)
Final position
* Similarly calculate the FEND from ZEND
…
FIMP = FEND / FSTART
Importance Ratio
CERN FLUKA Users Meeting
17
usimbs.f: Function example

Introduce an importance biasing assuming an exponential law of
attenuation in the R direction exp(-l  R), for R>1cm
RSTART = SQRT(XTRACK(0)**2 + YTRACK(0)**2 + ZTRACK(0)**2)
REND = SQRT(XTRACK(NTRACK)**2 + YTRACK(NTRACK)**2 +
ZTRACK(NTRACK)**2)
IF (RSTART .LT. ONEONE) THEN
FSTART = ONEONE
ELSE
FSTART = EXP(-ALAMBDA * RSTART)
ENDIF
IF (REND .LT. ONEONE) THEN
FEND = ONEONE
ELSE
FEND = EXP(-ALAMBDA * REND)
ENDIF
FIMP = FEND / FSTART
CERN FLUKA Users Meeting
18
Neutrino Decay Biasing

There is a special routine udcdrl.f where one can bias the
direction of the emitted neutrino in decays.
DOUBLE PRECISION FUNCTION UDCDRL( IJ, KPB, NDCY, UDCDRB, VDCDRB, WDCDRB)
Input variables:
IJ
KPB
Output variables:
U,V,W DCDRB
UDCDRL
decaying particle
outgoing neutrino
preferential outgoing direction for the neutrino
Lambda for direction biasing (1-cos(q))
The biasing expression is of the form:
e-(1-cosq/l)
Useful for neutrino applications like CNGS, Beta Beams…
 For a fixed direction the LAM-BIAS card with SDUM=DCY-DIRE could be
used instead

CERN FLUKA Users Meeting
19
Direction Biasing Example (n_TOF)

Bias the direction of the neutrino in the pion decay so the
daughter muon to be directed to the exp. area @(-120,0,18500).
The direction is given in the lab-frame, and in this example the
energies we are dealing are small, so safely we can assume that
also in the lab-frame, the neutrino and muon go in opposite
directions. The lambda is ¼ wide enough to cover the whole exp
area.
Direction in the lab frame
Width [1-cos(q)]
CERN FLUKA Users Meeting
20
Source routine

The user can prepare a biased source.f routine, when needed to
sample from biased distributions



Most often encountered to sample from a biased source energy
Common use also to sample from a biased spatial and/or angular
distribution
Using a two step approach reading particle histories from another
program or another simulation
Of course in such cases it is the user’s responsibility to ensure
that weights are adjusted in statistically correct way.
Tip:
The user should assign the correct weight to the primary
particles generated by the source routine WTFLK, but it is
recommended that the WEIPRI variable which contains the total
weight of the primary particles to be increased by one (instead of
the weight of the particle WTFLK). At the end of the run the user
should renormalize everything to the correct total weight of the
primary particles

CERN FLUKA Users Meeting
21
Source: Spatial distribution sampling

In this example we want to create a biased-uniform distribution
on a disc surface giving emphasis to the center
Sampling from 1/sqrt(r)
q = 2p FLRNDM(x)
SFE=sin(q)
CFE=cos(q)
CERN FLUKA Users Meeting
22
Source: Biased energy sampling

In this example the energy is given by the tabulated distribution
SPPROT(E). The sampling is divided into two parts below the one
using the functions 1/E-g and the other as 1/E
CERN FLUKA Users Meeting
23
Sampling from an Existing Event History

Applications:

Linking different codes (e.g., tracking and cascade codes)

Sampling a particle distribution given as external data file


Two-step approach in order to better sample a certain part of the
geometry (e.g., detectors)
It is important to carefully adjust the source particle weight
corresponding to the way how the original source was
calculated
1)
2)
3)
If the first step of a calculation included biasing the weights have to be
adjusted accordingly  Weight = original_Weight
Source particle corresponds to a single original event
 Weight = 1.0
(or Weight = original_Weight)
Selected source particle is only part of the original event
 Weight = 1.0 / #Particles_per_Event
(or Weight = original_Weight / #Particles_per_Event)
CERN FLUKA Users Meeting
24
Important for the Final Normalization
Tracking Code Example:

Ingredients:




Number of particles tracked (and hypothetically absorbed in all machine
components) in the tracking simulations:
 all_Tracked
Number of particles intercepted by the element(s) of interest:
 sel_Events
Ratio of particles on the considered machine element(s)
 sel_Events / all_Tracked
This ratio has then to be used to rescale the respective normalization valid
for the required application

e.g., particle loss rate (LossRate):  Nf =




ppb:
#b:
t:
particles per bunch
number of circulating bunches
beam lifetime
Final normalization:
Nf = LossRate x sel_Events/all_Tracked
CERN FLUKA Users Meeting
25
Sampling from a File - Part 1

using e.g., coordinates generated by tracking codes; 1st loading the file
* | | Load file given on SDUM of the SOURCE card -> FileName FILEN
WRITE (LUNOUT,*) 'Read source particles from '// FILEN
CALL OAUXFI(FILEN,LUNRDB,'OLD',IERR)
IF ( IERR .GT. 0 )
+
CALL FLABRT('SOURCE','Error opening file '//FILEN)
NIN = 0
10
CONTINUE
READ( LUNRDB, '(A)', ERR=9999, END=20 ) LINE
IF ( LINE(1:1) .EQ. '#' ) GO TO 10
NIN = NIN + 1
IF (NIN.GT.NINMAX) CALL FLABRT('SOURCE','Increase NINMAX')
READ (LINE,*) XIN(NIN), YIN(NIN), ZIN(NIN),
+
UIN(NIN), VIN(NIN)
GOTO 10
Read all coordinates to be sampled from!
20
CONTINUE
WRITE (LUNOUT,*) NIN,' particles loaded'
CLOSE(UNIT=LUNRDB)
CERN FLUKA Users Meeting
26
Sampling from a File - Part 2

2nd sampling it
Select a random particle
* | Choose a random particle
RNDSIG = FLRNDM (RNDSIG)
N = NINT(NIN*RNDSIG)+1
IF (N.GT.NIN) N=NIN
* TEST BEAM: uncomment the following line
IJIJ = IJBEAM
WEE = BEAWEI
* Cosines (tx,ty,tz)
TXI = UIN(N)
TYI = VIN(N)
TZI = SQRT ( ONEONE - TXI**2 - TYI**2 )
IF (WBEAM.LT.ZERZER) TZI = -TZI
* Particle coordinates
XXX = XIN(N)
YYY = YIN(N)
ZZZ = ZIN(N)
* Energy and Momentum
POO = PBEAM
EKE = SQRT ( PBEAM**2 + AM (IJBEAM)**2 ) - AM (IJBEAM)
ELKE = LOG (EKE)
Set the direction cosines
The coordinates
+ the energy and momentum
WTFLK (NPFLKA) = WEE
and the weight
CERN FLUKA Users Meeting
27
Bibliography
A Third Monte Carlo Sampler
( A revision and Extension of Samplers I and II)
C.J. Everett and E.D.Cashwell
LA-9721-MS, UC-32, March 1983
CERN FLUKA Users Meeting
28
Application:
n_TOF Background
FLUKA Users Meeting
CERN 31/05/2007
n_TOF Facility
High flux spallation neutron source

~200 m time of flight
Measurements of neutron induced cross-sections


ADS applications
Astrophysics
Proton Beam

CERN-PS: 20 GeV/c, Ip=4 bunches  7 1012 pr / 14.4 s
Neutron Source Features




Wide energy range: thermal up to GeV
High flux: ~ 6105 n/cm2/7 1012pr @ 200 m
Resolution: DE/E 210-4 @ 1 keV, 210-3 @ 1 MeV
Low background 10-6
FOR MORE INFO...
http://cern.ch/ntof, CERN/LHC/98-02 (EET), CERN/INTC/2000-004, CERN/SPSC 99-8
CERN FLUKA Users Meeting
30
n_TOF Virtual Tour
2nd Collimator
Lead Container
Sweeping Magnet
Exp. Area
CERN FLUKA Users Meeting
31
Background Problem June 2001
Background Features
 50 larger than simulations
 Three time components



400 ns “flash” (after the g-flash)
20 ms – fast neutrons
> 16 ms – thermal neutrons
Position depended
 Strong Left-Right asymmetry
 Strong ionization signal
 Not sample related

CERN FLUKA Users Meeting
32
Background Possible Sources
There were plenty of theories trying to describe the origin:
 Elements in the neutron Tube






Collimators
Escape line
Insufficient concrete shielding in the exp. area
Charged Particles deflected from the magnet
High energy neutron leaking from the target area
…
Challenges:
 Tiny solid angle
 Huge attenuation factors
CERN FLUKA Users Meeting
33
Topview of the Simulated Geometry
CERN FLUKA Users Meeting
34
1st Attempt: Elements on the neutron tube

Using a two (three) step approach:
1.
2.
3.

Simulating only the spallation target. Record particle histories with a
modified FLUSCW.F routine when particles enter in the neutron time
of flight tube.
Filter particles using a condition of vz/v > cos10o, cos5o, cos1o
bias slightly the direction with preference the element of interest,
e.g. 2nd collimator, escape line.
Full n_TOF simulation, using as source particles from step 2
Extra Biasing used:

Importance biasing with USIMBS.F subroutine. Increasing a factor of
x5 for every 50cm of concrete on the shielding upstream of the
experimental area and decreasing downstream
CERN FLUKA Users Meeting
35
USIMBS
Using the USRICALL to set dynamically the filename from the input
usrini.f
tt2a.inp
USRICALL

1.0
usrbias
The usrbias.dat file defines the Z position of each change in importance,
and the importance ratio with the previous Z-section
170.0
175.0
180.0
..
200.0
250.0
5.0
5.0
5.0
0.2
0.2
CERN FLUKA Users Meeting
36
USIMBS: Loading the file
CERN FLUKA Users Meeting
37
USIMBS: Importance ratio
CERN FLUKA Users Meeting
38
USIMBS: Binary Search
CERN FLUKA Users Meeting
39
1st) Background from NEL & Collimator
Neutron Background from:
Neutron Escape Line
Collimator
CERN FLUKA Users Meeting
40
2nd Attempt: Fast neutron streaming from the spallation area

There are visible paths where the shielding is not enough

Find weak points in geometry with the use of the RAY particle
CERN FLUKA Users Meeting
41
2nd) Minimum mass calculation





Special “simulation” using the RAY particle of FLUKA.
RAY is a straight-line trajectory through the geometry. The
program tracks all the objects lying in the given direction
calculating a number of various quantities like distance traversed
in each material, mass, number of interaction lengths etc…
Using a modified source.f routine a scan was performed on all the
possible directions starting from a vertical plane located at
Z=+7m (X  [-120, 280 cm], Y  [-120, 230 cm]), to a vertical
plane at Z=+185 m (X  [-250, 200 cm], Y  [-200, 210 cm])
from the center of the lead target.
500,000 rays were generated scanning the geometry starting
from a grid of 1010 on the former plane to a grid of 5050 to
the later plane and vice versa.
Recording on a 2D histogram the position for the trajectory were
the surface density is minimized.
CERN FLUKA Users Meeting
42
2nd) Minimum Mass Calculation
Revealed paths with minimum mass of 1.4 kg/cm2
 Equivalent to 2 GeV/c energy losses for Minimum ionizing particle
 Minimum ionizing particles like muons could easily penetrate the
experimental area
 High energy protons could also generate a background

CERN FLUKA Users Meeting
43
3rd) Simulation for High Energy neutrons



Starting directly from the protons on the lead target
Enabling Leading particle biasing and Weight windows
LAM-BIAS, WW-FACTOr, WW-THRESh
Using the USIMBS biasing:





Factor of 2 every 100cm around the spallation area
Factor of 2 every 10m in the tunnel (and neutron tube)
Factor of 5 every 50cm in the concrete shielding
With LOW-BIAS stop neutrons with E<1MeV in the spallation area
No EMF calculation
CERN FLUKA Users Meeting
44
3rd) Neutron Fluence (per fast neutron)
CERN FLUKA Users Meeting
45
3rd) Muon Fluence
Simulated muon fluence at the entrance of the exp. area per
proton pulse of 7 1012 protons
CERN FLUKA Users Meeting
46
1st,2nd,3rd) Results
Bad news
All studied possible background sources: Fast neutrons, Neutron
Escape line, collimator, gave consistent result of 0.01n/cm2/pulse
in the experimental area. 50-100 times lower than what was
measured
Good News
Strong muon component, that can explain the ionization signal
measured with the TLD’s and the strong asymmetry left/right in
the tunnel (fast-flash signal)
Last hope
Following the disappointing results from all previous efforts, last
resource was to leave running the simulation as long as possible,
searching blindly using brute force. 8 CPU’s
CERN FLUKA Users Meeting
47
Nice Surprise

After ~few weeks x 8 CPU’s of simulation one event appeared
with ~100 times higher importance than the rest.

To identify the origin of the event we needed to record the
history of all particles reaching the exp. area. This includes
tracking information (which regions traverse), last interaction
point, last decay point (if any), parent id, energy etc.
CERN FLUKA Users Meeting
48
Recording particle history (stuprf.f)
Extra USER information
for each particle
ISPARK: array of integers
SPAREK: array of floating point
CERN FLUKA Users Meeting
49
Write to file history info (fluscw.f)
Note the change of name
SPAREK  SPAUSR
ISPARK  ISPUSR
CERN FLUKA Users Meeting
50
Identification of the event
The event was identified as a ~1 MeV neutron generated in the concrete
surrounding the experimental area, originated from a negative muon
being stopped at the experimental area, generated by a pion decay
 Neutron generation by muons is a known process, important for
underground sites or atmospheric showers at ground level where muons:



Photo-nuclear interactions via virtual photons
Negative muon capture with m- brought to rest through the weak process
The former process involves muons of (relatively) high energy and has a
mean free path of a few hundred meters in earth (not important for
n_TOF).
 Negative muon capture occurs when m- brought to rest is captured
forming muonic atom. The capture usually occurs in a relatively loosely
bound state and then the muon jumps into progressively deeper levels
emitting hard x-rays. Due to the mass being much larger than the
electron one, the muon wave functions are significantly overlapping the
nucleus, and nuclear capture becomes competitive with decay.
 The assumptions used in FLUKA simulation models predict that ~50% of
m- stopping in concrete undergo nuclear capture. Experimental spectra of
neutrons emitted following m- captures show an evaporation peak, plus a
low intensity tail extending up to several tens of MeV.

CERN FLUKA Users Meeting
51
Muons arriving at Exp. Area
Blue: Decay
Position
Red: Inelastic
reaction where
the parent of
muon was
created
CERN FLUKA Users Meeting
52
Biasing to enhance the m- capture
Knowing the origin of the problem a more sophisticated biasing was
introduced to enhance the muon fluence at the experimental area.



Keep all previous settings of the high energy neutron biasing
Forcing the decay of the parent pions every one meter, to
increase the number of muons arriving at the experimental area
Reducing (Increasing) the mean life of the muons by a factor of
0.01, to enhance the capture probability versus the muon decay.
CERN FLUKA Users Meeting
53
Muon direction (udcdrl.f)
Biasing the direction of the neutrino from the pion decay
such as the daughter muon to go in the direction of the
experimental area with the routine UDCDRL
Direction in the lab frame
Width [1-cos(q)]
CERN FLUKA Users Meeting
54
Neutron Fluence at Exp. Area
Computed neutron fluence (n/cm2) inside the experimental
area (topview) for a standard pulse of 71012 protons
CERN FLUKA Users Meeting
55
Background neutron sources
CERN FLUKA Users Meeting
56
Neutron Time Distribution
CERN FLUKA Users Meeting
57
Possible Remedies
Stop the parent pions

Iron & Concrete shielding near the target area
Stop the muons


Place iron shielding near the target area
Iron shielding after the sweeping magnet
Absorb the neutron in the experimental area

Use borated polyethylene lining
CERN FLUKA Users Meeting
58
The 3m “muon” wall

The simulation results have clearly demonstrated the
ineffectiveness of a possible wall close to the target area:
 50% of parent pions are still in the tube at the exit of the target
shielding

10% of muons/parent pions are still in the tube as far as 60 m
from the target

Therefore a suitable shielding should be located where the
fraction of muons/pions in the pipe is minimal or just after the
sweeping magnet
3 m of Iron will lower the muon energy of 3.5 GeV
CERN FLUKA Users Meeting
59
Muon wall effectiveness
CERN FLUKA Users Meeting
60
Final Reduction
CERN FLUKA Users Meeting
61