|
|
Figure 1: Sketch of the geometrical setup of the current example. The size of the bodies is not in scale. |
Giuseppe Battistoni; INFN, Milano
Francesco Cerutti; CERN, Geneva
G. Battistoni 1, F. Cerutti 2
1 INFN, Milano
2 CERN, Geneva
This note describes with some detail a simple case in which FLUKA is linked to the CERN library in order to provide a widely spread output file, like a basic ntupla of HBOOK. At the same time, the case presented here gives an example on how to make use of the MGDRAW routine to access track information and generate a non standard FLUKA output tailored to users's needs.
WARNING:
1) MGDRAW is to be used by advanced users, and, generally, CANNOT BE USED WHEN BIASING IS ACTIVATED.
2) When the user is interested in the multigroup low energy neutron transport, SIMPLE HISTOGRAMMING, like the standard HBOOK case here described, IS MEANIGLESS.
For the above cases, only the standard FLUKA output, followed, if needed, by offline post-processing prepared by the user, can avoid the introduction of technical mistakes.
Sometimes the use of commonly used output environment can be of help for some users. A typical example, especially for users belonging to the High Energy Physics community, is the case of HBOOK histograms or ntuple. In order to illustrate how to make use of the HBOOK interface to generate an output from FLUKA, we shall consider the case of the test of the propagation of muons in sea water, where the goal is to predict the residual energy and the scattering angle of high energy muons which have crossed a considerable depth of sea water. This kind of calculation is typically necessary in the framework of the design of large neutrino telescopes, like the case of the ANTARES experiment. In this case, a typical approach is to consider muons at different fixed energies and propagate them through a layer of water having the proper chemical composition. In our example, the quantities of interest will be recorded into a HBOOK ntupla. In order to do this, the user will prepare some specific user routine to book, fill and save the HBOOK ntupla. However, there are many different ways which can be used to arrive at the same result, and this example does not claim to represent the only efficient way to accomplish this goal.
We shall describe this application for a unix (Linux) environment, where we suppose that a standard FLUKA release has been installed. Therefore we also assume that the environmental variable FLUPRO has been defined.
In the following, all codewords addressable in the FLUKA User Manual will be typed in bold characters. File names belonging to the FLUKA package will be typed in italic style.
In summary, the problem can be faced in the following way.
In the following we shall see in detail:
A simple problem like the one chosen for this example can be managed by means of a simple data card file like the following one, that we shall name as seamu.inp. It can be noticed that the NEW-DEFA set of defaults has been chosen, the main reason being a reasonably low value of particle transport threshold, set at 10 MeV, except for neutrons (19.6 MeV), antineutrons (50 MeV), and (anti)neutrinos (0, but they are discarded by default). However other features are implied, as here summarized:
We do not need all these features, and some of them must be over-written. The critical parts of this file will be examined in detail in the next sections.
Here is the full listing of seamu.inp:
TITLE
Test of Muon Propagation in sea water
*
* Sets the standards defaults of FLUKA (see user manual)
DEFAULTS NEW-DEFA
*
* Here sets muon energy. (1 TeV momentum in this example)
BEAM 1000.0 MUON+
*
* Here sets initial coordinates of muons. They are assumed to go
* along the +z axis. Muons are injected just before the water layer
* in a vacuum layer (see description of geometry below)
BEAMPOS 0.0 0.0 -50.0
*
* Switching off Electromagnetic Showering
EMF EMF-OFF
*
*
DISCARD ELECTRON POSITRON PHOTON HAD-CHAR HAD-NEUT
*
* Neutrinos are discarded by default
*
* Starts description of geometry: a large box of water delimited
* by vacuum layers at the top and bottom. Here depth is fixed
* at 1 km. Muons are injected in the initial vacuum.
GEOBEGIN COMBNAME
0 0 water layer
SPH LARGESPH 0.0 0.0 0.0 5000000.0
RPP WORLD -1000000.0 1000000.0 -1000000.0 1000000.0 -100.0 100100.0
XYP SEATOP 0.0
XYP SEABOT 100000.0
END
* black hole
BLACKHOL 5 +LARGESPH -WORLD
* vacuum at the beginning
VACTOP 5 +WORLD +SEATOP
* water layer
SEA 5 +WORLD +SEABOT -SEATOP
* vacuum at the end
VACBOT 5 +WORLD -SEABOT
END
GEOEND
*
* Here defines materials for the sea water compound (Antares 1)
*
* 1) Hydrogen
MATERIAL 1.0 1.0079 8.99D-05 3.0 1.0HYDROGEN
* 2) Oxygen
MATERIAL 8.0 15.999 0.001429 8.0 OXYGEN
* 3) Magnesium
MATERIAL 12.0 24.305 1.738 9.0 MAGNESIU
* 4) Sodium
MATERIAL 11.0 22.99 0.038349 19.0 SODIUM
* 5) Calcium
MATERIAL 20.0 40.08 1.54 21.0 CALCIUM
* 6) Potassium
MATERIAL 19.0 39.102 0.031165 26.0 POTASSIU
* 7) Chlorine
MATERIAL 17.0 35.4529 0.0029947 27.0 CHLORINE
* 8) Sulphur
MATERIAL 16.0 32.066 2.07 28.0 SULFUR
* 9) Sea water compound (Antares 1) using atom relative contents
MATERIAL 1.0341 29.0 SEAWATER
COMPOUND 2.0 HYDROGEN 1.0088 OXYGEN 0.00943 SODIUMSEAWATER
COMPOUND 0.000209 POTASSIU 0.001087 MAGNESIU 0.000209 CALCIUMSEAWATER
COMPOUND 0.01106 CHLORINE 0.00582 SULFUR SEAWATER
*
* Here assigns materials to regions
* Sea Water in region 3
ASSIGNMA SEAWATER SEA
* External Black Hole in region 1
ASSIGNMA BLCKHOLE BLACKHOL
* Vacuum at beginning (region 2)
ASSIGNMA VACUUM VACTOP
* Vacuum at end (region 4)
ASSIGNMA VACUUM VACBOT
*
* Setting SCORE and Ouptut levels
SCORE ENERGY
*
* Activates the routine MGDRAW (what(1) = or > 100.)
* with what(3)=2 the entry BXDRAW in MGDRAW routine is activated
USERDUMP 101.0 2.0 SEAMU
*
* Activates Random Number initialization
RANDOMIZ 1.0
*
* Activates user dependent initialization
USRICALL
*
* Starts Job asking for a given number of events
START 1000.0
*
* Activates user dependent output
USROCALL
STOP
We draw the attention on the presence of a USERDUMP card with what(3) = 2. As remarked in the comments, this enables the call to the entry BXDRAW in the MGDRAW routine, where we introduce our specific code lines. The entry BXDRAW is called also for what(3) = 3,4,5,6 or what(3) ≤ 0, but with what(3) = 2 the number of calls to the other MGDRAW entries (here not necessary) is minimized.
Beam properties ::::::::::::The typical job requires the calculation of high energy muon transportation through a depth of water of the order of few km. Muons should be considered at different energies and angles. For simplicity we shall consider particles at normal incidence with respect to the water surface. The muon beam can be point-like in cross section, and the energy can vary between few hundreds of GeV to many TeV's.
It can be considered convenient, also for book-keeping purposes, to steer the primary (muon) properties by means of the BEAM and BEAMPOS cards, performing different, separated, runs according to the muon energy and angle. This way there is no need to prepare a user specific source code (namely adapting the SOURCE routine contained in the source.f file).
In view of the design of the geometry description, we can also choose the default propagation of primary particles along the positive Z axis.
Geometry description ::::::::::::Since the main purpose of this example is to describe the technicalities of interfacing the HBOOK routines, we shall not consider a sophisticated geometrical description and we shall limit ourselves to the simple approximation of a single large volume uniformly filled by water having homogeneous properties. This can be easily obtained using a Rectangular Parallelepiped body (RPP) with one edge parallel to the propagation axis (Z).
In case variation of material properties (such as the density) as a function of depth must be considered, one should prepare different bodies to be positioned one after the other, so to fill them with different materials.
We remind here that all physical regions must be embedded inside an enclosing volume, defined to be a "black hole" material, so to stop the tracking of all escaping particles. We can choose a large sphere (SPH) to describe the external boundary of this enclosing volume.
Since muons must be injected just above the water surface and since, by definition, they cannot be injected in a "black body" medium, it can be convenient to define a limited vacuum region just above the water volume. It is convenient to do the same just after the depth at which scoring is performed. All this can be easily implemented by cutting the RPP volume by two planes orthogonal to the Z axis (XYP bodies), respectively at Z = Zmin and Z = Zmax such that Zmax – Zmin = Depth, where Depth is the desired thickness of water. The extreme Z coordinates of the RPP body, Z1 < Zmin and Z2 > Zmax, will be chosen so that Zmin – Z1 and Z2 – Zmax will respectively define the desired vacuum layers above and below the water volume.
In summary, the geometrical setup will be the one depicted in Fig. 1.
|
|
|
Figure 1: Sketch of the geometrical setup of the current example. The size of the bodies is not in scale. |
The bodies can be defined in the geometry data cards (free format with names) according to the following lines:
SPH LARGESPH 0.0 0.0 0.0 5000000.0 RPP WORLD -1000000.0 1000000.0 -1000000.0 1000000.0 -100.0 100100.0 XYP SEATOP 0.0 XYP SEABOT 100000.0 END
Here Z1 = –100 cm and Z2 = 100100 cm, while Zmin = 0 and Zmax = 100000 cm. The initial and final vacuum layers are therefore 100 cm deep.
From these bodies, the following regions (to be assigned a specific material) are defined:
* black hole BLACKHOL 5 +LARGESPH -WORLD * vacuum at the beginning VACTOP 5 +WORLD +SEATOP * water layer SEA 5 +WORLD +SEABOT -SEATOP * vacuum at the end VACBOT 5 +WORLD -SEABOT END
The black hole will be the 1st region (named BLACKHOL): the space contained within the body LARGESPH (SPH) and outside the body WORLD (RPP).
The initial vacuum will be the region #2 (named VACTOP): the space within the body WORLD (RPP) and above the XYP plane named SEATOP.
The water volume will be the 3rd region (named SEA): the space within the body WORLD (RPP) and delimited by the two XYP planes (SEATOP and SEABOT).
The 4th region (named VACBOT) is at the bottom boundary of the sea-water region SEA (within the body WORLD (RPP) and below the XYP plane named SEABOT) and will be defined as a vacuum region.
Description of the material ::::::::::::The region SEA has to be filled with sea water, which has been defined as material SEAWATER (#29). In our example we have adopted one of the mixtures in use by the ANTARES experiment. In that case, sea-water is represented by the compound of 8 elements: hydrogen, oxygen, magnesium, sodium, calcium, potassium, chlorine, and sulphur. The relative proportions are given in terms of average atom content.
Adapting the parameters of the physical processes ::::::::::::In this problem we are interested only in the energy loss of the muons, with a high level of accuracy. This can be achieved only when the explicit production of secondaries is activated, but at the same time we do not need either to follow all the particles produced in muon interactions or to simulate e.m. showers. This is the reason why we have switched off the transport of e±, photons, all charged and neutral hadrons by means of the DISCARD directive (neutrinos are not transported by default) and we have asked for EMF-OFF in the EMF card (the presence of the EMF card is indeed irrelevant when e.m. particles are discarded). This way we minimize the impact in terms of computing time.
The user routines to be linked in FLUKA ::::::::::::As mentioned in the introduction, in this example the user has to prepare 4 routines: USRINI, HISTIN, MGDRAW and USROUT. Let us examine them in the following
The user routines to be linked in FLUKA: USRINI ::::::::::::Here only the call to HISTIN is added to the dummy code.
*$ CREATE USRINI.FOR
*COPY USRINI
*
*=== usrini ===========================================================*
*
SUBROUTINE USRINI ( WHAT, SDUM )
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
* *
* Copyright (C) 1991-2005 by Alfredo Ferrari & Paola Sala *
* All Rights Reserved. *
* *
* *
* USeR INItialization: this routine is called every time the *
* USRICALL card is found in the input stream *
* *
* *
* Created on 01 january 1991 by Alfredo Ferrari & Paola Sala *
* Infn - Milan *
* *
* Last change on 20-mar-05 by Alfredo Ferrari *
* *
* *
*----------------------------------------------------------------------*
*
DIMENSION WHAT (6)
CHARACTER SDUM*8
*
* Don't change the following line:
LUSRIN = .TRUE.
* *** Write from here on *** *
CALL HISTIN ! Call to HBOOK initialization
RETURN
*=== End of subroutine Usrini =========================================*
END
The user routines to be linked in FLUKA: HISTIN
::::::::::::
The main attention in the case of the use of HBOOK is that it makes use of REAL*4 variables, while the FLUKA environment is by default in double precision (REAL*8). For our example we intend to book a simple disk-resident ntupla opening the hbk.his file, having 8 variables: the muon charge, the residual energy, the 3 coordinates and the 3 direction cosines at boundary crossing. Notice that we have chosen to assign the logical unit 63 to the HBOOK file. Any number larger than 20 is good, provided that it is not already used by the user himself for other scoring purposes.
C------------------------------------------------------------
SUBROUTINE HISTIN
C------------------------------------------------------------
C
C HBOOK HISTOGRAM BOOKING
C
PARAMETER (NNTUPLE=8)
COMMON/PAWC/HMEMOR(200000)
REAL*4 HMEMOR
COMMON/QUEST/IQUEST(100)
CHARACTER*8 CHTAG(NNTUPLE)
DATA CHTAG /'MUTYP' ,'ERES' ,'XMU' ,'YMU', 'ZMU', 'THX',
$ 'THY', 'THZ'/
WRITE(*,*) 'ENTER INTO HIST INIT.'
IQUEST(10)=65000
CALL HLIMIT(200000)
C
NLUNI1 = 63
CALL HROPEN(NLUNI1, 'Muons','hbk.his','N',1024,ISTAT)
WRITE(*,*) 'HROPEN: ISTAT=',ISTAT
CALL HBOOKN(100,'Muons',NNTUPLE,'//Muons',NNTUPLE*100,CHTAG)
C
RETURN
END
The user routines to be linked in FLUKA: MGDRAW
::::::::::::
The user must define locally a REAL*4 vector to fill the ntupla, XTUP. This declaration can be placed at the beginning of MGDRAW routine, just after the inclusion of standard commons.
....
*
PARAMETER (NNTUPLE=8)
REAL*4 XTUP(NNTUPLE)
*
....
In this routine the user has to comment the WRITE fortran cards above the ENTRY BXDRAW (the other standard WRITE fortran cards are not activated since we put what(3) = 2 in the USERDUMP card), which are intended to write the file SEAMU indicated in the USERDUMP card: this is not useful here, as the relevant information is stored in a basic ntupla of HBOOK.
Then we show only the relevant section relative to ENTRY BXDRAW. The meaning of the relevant variables in BXDRAW is the following:
* Jtrack = identity number of the particle * * Etrack = total energy of the particle * * Am = mass energy of the particle * * Cx,y,ztrck = direction cosines of the current particle * * Mreg = region number before the boundary crossing * * Newreg = region number after the boundary crossing * * Xsco,Ysco,Zsco = coordinates of the particle at boundary crossing *
In order to select the desired boundary, we remember that in our geometry description, the detector depth was defined as the crossing between region 3 and 4.
....
*
*======================================================================*
* *
* Boundary-(X)crossing DRAWing: *
* *
* Icode = 1x: call from Kaskad *
* 19: boundary crossing *
* Icode = 2x: call from Emfsco *
* 29: boundary crossing *
* Icode = 3x: call from Kasneu *
* 39: boundary crossing *
* Icode = 4x: call from Kashea *
* 49: boundary crossing *
* Icode = 5x: call from Kasoph *
* 59: boundary crossing *
* *
*======================================================================*
* *
ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO )
C-STA
WRITE(*,*) 'BXDRAW ',ICODE,MREG,NEWREG,JTRACK
$ ,ETRACK,AM(JTRACK)
IF(MREG.EQ.3.AND.NEWREG.EQ.4) THEN ! Select the desired boundary
IF( JTRACK.EQ.10.OR.JTRACK.EQ.11 ) THEN ! Select muons
IF(ETRACK.GT.AM(JTRACK)) THEN ! Muon has survived
XTUP(1) = JTRACK
XTUP(2) = ETRACK-AM(JTRACK)
XTUP(3) = XSCO
XTUP(4) = YSCO
XTUP(5) = ZSCO
XTUP(6) = CXTRCK
XTUP(7) = CYTRCK
XTUP(8) = CZTRCK
CALL HFN(100,XTUP)
ENDIF
ENDIF
ENDIF
C-END
RETURN
*
....
The user routines to be linked in FLUKA: USROUT
::::::::::::
Here the ntupla and the file containing it are simply closed invoking the recommended HBOOK routines.
*$ CREATE USROUT.FOR
*COPY USROUT
*
*=== usrout ===========================================================*
*
SUBROUTINE USROUT ( WHAT, SDUM )
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
* *
* Copyright (C) 1991-2005 by Alfredo Ferrari & Paola Sala *
* All Rights Reserved. *
* *
* *
* USeR OUTput: this routine is called every time the USROCALL card *
* is found in the input stream *
* *
* *
* Created on 01 january 1991 by Alfredo Ferrari & Paola Sala *
* Infn - Milan *
* *
* Last change on 20-mar-05 by Alfredo Ferrari *
* *
* *
*----------------------------------------------------------------------*
*
DIMENSION WHAT (6)
CHARACTER SDUM*8
*
CALL HROUT(0,ICYCLE,' ')
CALL HREND('Muons')
NLUNI1 = 63
CLOSE(NLUNI1)
RETURN
*=== End of subroutine Usrout =========================================*
END
Compiling the executable module
::::::::::::
The executable module can be prepared in different ways.
The user routines must be compiled with the fff script:
$FLUPRO/flutil/fff usrini.f $FLUPRO/flutil/fff histin.f $FLUPRO/flutil/fff usrout.f $FLUPRO/flutil/fff mgdraw.f
The resulting object files can be inserted into a user library, for instance creating the file libuser.a. This is done by means of the command ar:
ar -r libuser.a usrini.o usrout.o histin.o mgdraw.o
The order of the modules is irrelevant.
Whenever one of the user routines has to be modified, it can be replaced at any time, after compilation, inside libuser.a:
ar -r libuser.a filename.o
The content of the library file, and the dates of the last update of each module, can be checked by means of the command:
ar -tv libuser.a
At this point the executable file, which we shall name for instance flukaseamu, can be linked:
$FLUPRO/flutil/lflukac -m fluka -C -O user -o flukaseamu
The option -C is necessary to invoke the link to the CERN library.
Warning:Before performing this operation, the user must check that inside the lflukac script the environmental variable CERNPATH points to the actual address of CERN library in the user's machine; for instance:
CERNPATH=/cern/2003/lib
The environmental variable X11PATH has to be checked as well.
It is important to adopt the option -O to address the external user library because it contains routines which must be substituted in the original FLUKA library (libflukahp.a).
The executable module can also be prepared avoiding the creation of a user library:
$FLUPRO/flutil/lflukac -m fluka -C usrini.f usrout.f histin.f mgdraw.f -o flukaseamu
or, if the source files have been already compiled by means of $FLUPRO/flutil/fff:
$FLUPRO/flutil/lflukac -m fluka -C usrini.o usrout.o histin.o mgdraw.o -o flukaseamu
The order in which the user's files are presented is irrelevant.
With respect to the case where a user library is used, this is a more direct way to link the executable, and for simple problems like the one described here, the difference is irrelevant. In cases in which the number of user routines is not so small, the use of a user library is much more convenient and safer, and is strongly recommended.
The number of cases to be run in a single cycle is given in the first numeric field of the directive START: 1000 in this example. Let us suppose that we want to run three independent cycles (the RANDOMIZE directive has been issued in the data card input in order to have different random number seeds for each cycle). The command line will thus be:
$FLUPRO/flutil/rfluka -e flukaseamu -N0 -M3 seamu
At the end of each cycle the HBOOK file will be named seamu00n_hbk.his where n = 1, 2, 3.
Managing the output files ::::::::::::The total number of handled cases must be checked in the seamu00n.out files. The number of survived muons will be given by the number of entries in the ntupla. In fig. 2 we show the plot of muon residual energy and in fig. 3 we plot the distribution of transverse coordinates of survived muons at detection level.
|
| Figure 2: Distribution of residual energy of survived muons (E0 = 1 TeV) at the depth of 1 km of sea water. |
|
| Figure 3: Distribution of the X, Y coordinates of survived muons (E0 = 1 TeV) at detection level (after 1 km of sea water). |
These plots can be obtained by means of a simple "kumac" like the following one:
* merging hbook ntuples nt/hmerge seamu.his seamu001_hbk.his seamu002_hbk.his seamu003_hbk.his his/fil 1 seamu.his * checking the content of the ntupla nt/print 100 * booking the 1D histogram of residual energy 1dh 10 ' ' 100 0. 1000. 0. * projecting on the histogram nt/proj 10 100.eres his/plo 10 atit 'E?[m]! GeV' * booking a 2D histogram of X vs Y muon coordinates 2dh 20 ' ' 100 -200. 200. 100 -200. 200. 0. nt/proj 20 100.ymu%xmu opt nsta h/plo 20 surf2 atit 'X (cm)' 'Y (cm)'