Tutorial: Dilepton simulations with GiBUU

(Janus Weil, 30.11.2016)

Overview and general principles

The dilepton analysis is a special component of GiBUU that is not enabled by default and needs to be turned on manually if desired. Since dilepton decays are particularly rare, they are not treated like normal hadronic decays in GiBUU, but are handled with a ‘perturbative’ treatment, the so-called shining method. Note that currently only di-electron decays (e+e-) are supported, but no di-muons (mu+mu-).

Normal hadronic decays work like this: In each timestep the probability for each particle to decay within this timestep is computed, then a Monte-Carlo decision is done according to this probability whether the decay actually happens (and if yes, into which channel). Then the decayed particle is removed from the particle vector and the decay products are inserted.

In contrast, the dilepton decays do not work on the particle vector. Leptons are not inserted there, so that they can not rescatter (this is a rather good approximation since they only interact electromagnetically). Instead the dilepton-analysis routines loop over the particle vector and check for each particle if it has any dilepton decay channels. All of these decays are always(!) performed (for each particle in each timestep), in the sense that a dilepton pair is created and written directly into the output files (either a histogram or an event file). The mother particle is not removed from the particle vector, so that it can and will radiate off more dileptons in the following timesteps. All of this can only work by associating a certain weight factor with each lepton pair, which encodes the probability with which this decay happened. Dilepton spectra are generated by summing up the weights of all involved lepton pairs.

The currently implemented physical processes are:

• rho, omega, phi, eta → e+e-
• pi0, eta, eta-prime → gamma e+e-
• omega → pi0 e+e-
• Delta → N e+e-
• N*, Delta* → N e+e-
• Bremsstrahlung (NN → NN e+e-, piN → piN e+e-)
• Bethe-Heitler process (gamma N → N e+e-)

Note that the dilepton analysis has only been verified to work for the event types RealPhoton, HiPion, HiLepton, HeavyIon and HadronInduced.

Relevant program parts

All the code that is relevant for the dilepton analysis lives in the file code/analysis/DileptonAnalysis.f90, apart from a few routines related to filtering and form factors, which are in separate files in the analysis directory. The most important routines are Dilep_Init, which has to be called in the beginning to read the input and initialize everything, and Dilep_Decays, which is the ‘work horse’ that performs all implemented decays (and calls various other routines). Furthermore there are the additional routines Dilep_Brems, which produces the Bremsstrahlung contribution that comes not from decays but from 2-body collisions, and Dilep_BH, which implements the Bethe-Heitler contribution for photon-induced reactions. And finally there is Dilep_write_CS, which writes out all the output in the end.

Namelist Input

All the input parameters that influence the dilepton analysis are collected in the namelist /DileptonAnalysis/. All of the parameters are documented and explained in the Robodoc documentation on the website (​https://gibuu.hepforge.org/Documentation/code/robo_namelist.html). The Enable parameter is the ‘master switch’. It is used to turn on the dilepton analysis at all. When set to .true., the dilepton analysis will be performed and basic output is written to several histogram files (that contain different dilepton spectra). In order to compare the generated spectra to experimental data, it is usually necessary to run the lepton pairs through an acceptance filter, which models the geometry and acceptance of a given detector. By default the spectra in GiBUU are unfiltered, but the switch filter can be used to generate filtered spectra that are tailored to a specific detector (possible choices are: DLS, HADES, g7/CLAS, KEK E325, JPARC E16). Note that the implemented filters vary strongly in complexity, reaching from a few simple cuts (KEK/JPARC) to a full-blown acceptance and smearing machinery that models many detector details (HADES). In addition to the standard histogram output, a detailed output of dilepton events (including full momenta of all pairs) can be generated via the switch WriteEvents? if necessary.

Generated Output

The are two kinds of output that the dilepton analysis can generate: histogram output (default) and event output (optional). The histogram output consists of several histogram files that contains various spectra, e.g. spectra of the mass, transverse momentum or rapidity of the lepton pairs. Those files are named DileptonMass.dat, DileptonPt.dat, DileptonY.dat, etc. Each histogram contains in the first column the corresponding quantity (m, pT, y) and in the second column the total dilepton yield. There are several columns following after that with all the different contributions from the various channels. Each histogram file has a header which explains the contents. The event output is more detailed and much larger in size. It contains all the lepton pairs with their full momenta. It is useful if one wants to use GiBUU as an event generator for dilepton events, to filter the events through an external filter or to produce spectra that are not included in the default histograms. There are different variants of event output, which print either just the lepton pair, or full events including all hadronic particles.

Jobcards

There are several example jobcards for making dilepton runs in the directory testRun/jobcards/:

• 003_gammaDilep.job --- gamma + Ca at 1.5 GeV
• 012_pionDilep.job --- pi- + Ca at 1.3 GeV
• 012_protonDilep.job --- p + Ca at 3.5 GeV
• 012_pp35_dilep.job --- p + p at 3.5 GeV (with HADES acc. filtering)

The energies given are the kinetic energies of the beam particle in the lab frame (fixed target).

Plotting

The typical way of producing result plots is to just plot the contents of the histogram files (DileptonMass.dat etc.), e.g. via gnuplot or python. It is recommended to not plot the 2nd column (“total spectrum”) directly, because this sum of all the channels can contain unwanted/unreasonable contributions or even double-count some channels. Instead the sum should be taken manually of all the channels plotted (which is simple in most plotting programs). In order to improve statistics, it can be useful to make several runs with the same jobcard (possibly in parallel on a cluster), and then combine the results. For summing up the contents of several histogram files from different runs, one can use the auxiliary program sumHistMC (which can be found in the directory code/storage/testRun, together with the corresponding jobcard jobSumMC).

Result Plot

The shown plot was generated from the jobcard 012_pp35_dilep.job, using the sum of 24 runs from this jobcard, plotted with a gnuplot script. It contains the spectrum for a proton-proton collision at a kinetic energy of 3.5 GeV, filtered through the HADES acceptance and smeared with the detector resolution.

##### Event Characteristics #####

Energy = 3.5
e_str = sprintf("%4.1f GeV",Energy)

proj = "p"
targ = "p"

##### column descriptions #####

t2 = "GiBUU total"
t3 = "{/Symbol r \256 }e^+e^-"
t4 = "{/Symbol w \256 }e^+e^-"
t5 = "{/Symbol f \256 }e^+e^-"
t6 = "{/Symbol w \256 p^0}e^+e^-"
t7 = "{/Symbol p^0 \256 }e^+e^-{/Symbol g}"
t8 = "{/Symbol h \256 }e^+e^-{/Symbol g}"
t9 = "{/Symbol D \256 }Ne^+e^-"
t10 = "{/Symbol h \256 }e^+e^-"
t14 = "pp Brems"

##### set default line styles #####

w  = 4
ww = 6

set style line 1 lt 1 lw ww
set style line 2 lt 2 lw w lc rgb "web-green"
set style line 3 lt 3 lw w
set style line 4 lt 4 lw w
set style line 5 lt 5 lw w
set style line 6 lt 6 lw w lc rgb "dark-violet"
set style line 7 lt 7 lw w
set style line 8 lt 8 lw w
set style line 9 lt 9 lw w

set style line 21 lt 3 lw w lc rgb "web-green"
set style line 22 lt 4 lw w lc rgb "web-green"
set style line 23 lt 5 lw w lc rgb "web-green"

##### function to sum over columns i to j #####

colsum(i,j) = ( i==j ? column(i) : column(i)+colsum(i+1,j) )

########################## mass spectrum ##########################

set terminal postscript eps enhanced color dashed "Helvetica" 20 dl 2 size 4,4

set title proj . " + " . targ . " \\100" . e_str
set xlabel 'dilepton mass m_{ee} [GeV]'
set ylabel 'd{/Symbol s}/dm_{ee} [{/Symbol m}b/GeV]'

set logscale y
set xrange [0:1.1]
set mxtics 2
set mytics 10
set yrange [5e-4:5e1]
set format y "10^{%L}"

set key samplen 3 box opaque

unset bars

set style data lines

FF = "DileptonMass.dat"

set output 'pp35_mass.eps'

plot \
     FF u 1:($4+$5+$6+$7+$8+$13+colsum(21,39)) ls 1 t t2, \
     FF u 1:4 ls 3 t t4, \
     FF u 1:5 ls 4 t t5, \
     FF u 1:6 ls 5 t t6, \
     FF u 1:7 ls 6 t t7, \
     FF u 1:8 ls 7 t t8 , \
     FF u 1:9 ls 8 t "{/Symbol D} QED", \
     FF u 1:(colsum(21,33)) ls 2 t "N* VMD", \
     FF u 1:(colsum(34,39)) ls 22 t "{/Symbol D}* VMD", \
     FF u 1:14 ls 9 t t14

!epstopdf pp35_mass.eps