1 / 61

FLUKA Advanced Biasing

FLUKA Advanced Biasing. FLUKA Users Meeting CERN 31/05/2007. Outline. Sampling Random numbers Sampling from discrete & generic distributions Sampling from biased distributions FLUKA User biasing routines USIMBS UDCDRL SOURCE Application: n_TOF biasing example. Pseudo-Random number.

gloriann-stevens
Télécharger la présentation

FLUKA Advanced Biasing

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. FLUKAAdvanced Biasing FLUKA Users Meeting CERN 31/05/2007

  2. Outline • Sampling • Random numbers • Sampling from discrete & generic distributions • Sampling from biased distributions • FLUKA User biasing routines • USIMBS • UDCDRL • SOURCE • Application: n_TOF biasing example CERN FLUKA Users Meeting

  3. Pseudo-Random number • All computers provide a pseudo random number generator. • The FLUKA PRNG is provided by the functionFLRNDM(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

  4. FLUKA Random number generators • Function FLRNDM(DUMMY)The main FLUKA PRNG. Returns a double precision random number uniformly distributed in the [0,1) interval • SubroutineSFECFE(SFE, CFE)returns a random azimuthal sine and cosine • SubroutineRACO(TXX, TYY, TZZ)returns the cosines of a randomly oriented vector • SubroutineSFLOOD(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) • SubroutineFLNRRN(RGAUSS)returns a random number normally distributed with unit variance CERN FLUKA Users Meeting

  5. 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 isWARNING: 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 xis uniformly random, we have P(xi) = P(yi-1 < x < yi) = yi – yi-1 = pi CERN FLUKA Users Meeting

  6. 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 • Sample a uniform pseudo-random number x • Get the desired result by finding the inverse valuex=F-1(x), analytically or most often by interpolation (table lookup) • Since x is uniformly random, we have CERN FLUKA Users Meeting

  7. 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

  8. 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=x1otherwise re-sample x1, x2 • The probability that a x1 value isaccepted is f(x1)/C for any x1P(x)dx = P(x1=x)dx  f(x1=x)/C • Since P(x1)dx = const,P(x) = const  f(x) CERN FLUKA Users Meeting

  9. Sampling from a biased distribution CERN FLUKA Users Meeting

  10. Example CERN FLUKA Users Meeting

  11. Uniform sampling from a steep distribution CERN FLUKA Users Meeting

  12. User written biasing FLUKA offers the following routines for user-written biasing • ubsset.f: User BiaSing SETting • called after reading in the input file and before the first event • allows to alter almost any biasing weight on a region-dependent basis • usimbs.f: USer defined IMportance BiaSing • 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

  13. 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

  14. 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 (Any value) • WHAT(4) Lower bound of region • WHAT(5) Upper bound of region • WHAT(6) Step • SDUM = USER CERN FLUKA Users Meeting

  15. 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 usimbs.f – Routine CERN FLUKA Users Meeting

  16. 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

  17. usimbs.f: Simple Example • Biasing a factor of 2 every 50 cm on the Z direction from 100cm to 500cm 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) * Similarly calculate the FEND from ZEND … FIMP = FEND / FSTART Initial position Final position Importance Ratio CERN FLUKA Users Meeting

  18. 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

  19. Neutrino Decay Biasing • There is a special routine udcdrl.f where one can bias the direction of the emitted neutrino in decays. DOUBLE PRECISION FUNCTIONUDCDRL( IJ, KPB, NDCY, UDCDRB, VDCDRB, WDCDRB) Input variables: IJ decaying particle KPB outgoing neutrino Output variables: U,V,W DCDRB preferential outgoing direction for the neutrino UDCDRL 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

  20. 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

  21. 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

  22. 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

  23. 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

  24. 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 • 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

  25. 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: particles per bunch • #b: number of circulating bunches • t: beam lifetime • Final normalization:Nf = LossRate x sel_Events/all_Tracked CERN FLUKA Users Meeting

  26. 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 20 CONTINUE WRITE (LUNOUT,*) NIN,' particles loaded' CLOSE(UNIT=LUNRDB) Read all coordinates to be sampled from! CERN FLUKA Users Meeting

  27. Sampling from a File - Part 2 Select a random particle • 2nd sampling it * | 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) WTFLK (NPFLKA) = WEE Set the direction cosines The coordinates + the energy and momentum and the weight CERN FLUKA Users Meeting

  28. 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

  29. Application: n_TOF Background FLUKA Users Meeting CERN 31/05/2007

  30. 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

  31. n_TOF Virtual Tour 2nd Collimator Lead Container Sweeping Magnet Exp. Area CERN FLUKA Users Meeting

  32. 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

  33. 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

  34. Topview of the Simulated Geometry CERN FLUKA Users Meeting

  35. 1st Attempt: Elements on the neutron tube • Using a two (three) step approach: • 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,cos1obias 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

  36. 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 5.0 175.0 5.0 180.0 5.0 .. 200.0 0.2 250.0 0.2 CERN FLUKA Users Meeting

  37. USIMBS: Loading the file CERN FLUKA Users Meeting

  38. USIMBS: Importance ratio CERN FLUKA Users Meeting

  39. USIMBS: Binary Search CERN FLUKA Users Meeting

  40. 1st) Background from NEL & Collimator Neutron Background from: Neutron Escape Line Collimator CERN FLUKA Users Meeting

  41. 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

  42. 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

  43. 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

  44. 3rd) Simulation for High Energy neutrons • Starting directly from the protons on the lead target • Enabling Leading particle biasing and Weight windowsLAM-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

  45. 3rd) Neutron Fluence (per fast neutron) CERN FLUKA Users Meeting

  46. 3rd) Muon Fluence Simulated muon fluence at the entrance of the exp. area per proton pulse of 7 1012 protons CERN FLUKA Users Meeting

  47. 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

  48. 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

  49. Recording particle history (stuprf.f) Extra USER information for each particle ISPARK: array of integers SPAREK: array of floating point CERN FLUKA Users Meeting

  50. Write to file history info (fluscw.f) Note the change of name SPAREK  SPAUSR ISPARK  ISPUSR CERN FLUKA Users Meeting

More Related