gibuu is hosted by Hepforge, IPPP Durham

Changes between Initial Version and Version 1 of dileptons

Feb 6, 2017, 2:57:47 PM (3 years ago)

copy text into wiki


  • dileptons

    v1 v1  
     1= Tutorial: Dilepton simulations with GiBUU =
     3(Janus Weil, 30.11.2016)
     5== Overview and general principles ==
     7The 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-).
     9Normal 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.
     11In 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.
     13The currently implemented physical processes are:
     14* rho, omega, phi, eta → e+e-
     15* pi0, eta, eta-prime → gamma e+e-
     16* omega → pi0 e+e-
     17* Delta → N e+e-
     18* N*, Delta* → N e+e-
     19* Bremsstrahlung (NN → NN e+e-, piN → piN e+e-)
     20* Bethe-Heitler process (gamma N → N e+e-)
     22Note that the dilepton analysis has only been verified to work for the event types !RealPhoton, !HiPion, !HiLepton, !HeavyIon and !HadronInduced.
     25== Relevant program parts ==
     27All 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.
     28The 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.
     31== Namelist Input ==
     33All 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 ( 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).
     34In 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).
     35In 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.
     38== Generated Output ==
     40The are two kinds of output that the dilepton analysis can generate: histogram output (default) and event output (optional).
     41The 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.
     42The 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.
     45== Jobcards ==
     47There are several example jobcards for making dilepton runs in the directory `testRun/jobcards/`:
     49* 003_gammaDilep.job  --- gamma + Ca at 1.5 GeV
     50* 012_pionDilep.job       --- pi- + Ca at 1.3 GeV
     51* 012_protonDilep.job   --- p + Ca at 3.5 GeV
     52* 012_pp35_dilep.job    --- p + p at 3.5 GeV (with HADES acc. filtering)
     54The energies given are the kinetic energies of the beam particle in the lab frame (fixed target).
     57== Plotting ==
     59The 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.
     60It 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).
     61In 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`).
     64== Result Plot ==
     69The 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.
     70It 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.
     73##### Event Characteristics #####
     75Energy = 3.5
     76e_str = sprintf("%4.1f GeV",Energy)
     78proj = "p"
     79targ = "p"
     81##### column descriptions #####
     83t2 = "GiBUU total"
     84t3 = "{/Symbol r \256 }e^+e^-"
     85t4 = "{/Symbol w \256 }e^+e^-"
     86t5 = "{/Symbol f \256 }e^+e^-"
     87t6 = "{/Symbol w \256 p^0}e^+e^-"
     88t7 = "{/Symbol p^0 \256 }e^+e^-{/Symbol g}"
     89t8 = "{/Symbol h \256 }e^+e^-{/Symbol g}"
     90t9 = "{/Symbol D \256 }Ne^+e^-"
     91t10 = "{/Symbol h \256 }e^+e^-"
     92t14 = "pp Brems"
     94##### set default line styles #####
     96w  = 4
     97ww = 6
     99set style line 1 lt 1 lw ww
     100set style line 2 lt 2 lw w lc rgb "web-green"
     101set style line 3 lt 3 lw w
     102set style line 4 lt 4 lw w
     103set style line 5 lt 5 lw w
     104set style line 6 lt 6 lw w lc rgb "dark-violet"
     105set style line 7 lt 7 lw w
     106set style line 8 lt 8 lw w
     107set style line 9 lt 9 lw w
     109set style line 21 lt 3 lw w lc rgb "web-green"
     110set style line 22 lt 4 lw w lc rgb "web-green"
     111set style line 23 lt 5 lw w lc rgb "web-green"
     113##### function to sum over columns i to j #####
     115colsum(i,j) = ( i==j ? column(i) : column(i)+colsum(i+1,j) )
     117########################## mass spectrum ##########################
     119set terminal postscript eps enhanced color dashed "Helvetica" 20 dl 2 size 4,4
     121set title proj . " + " . targ . " \\100" . e_str
     122set xlabel 'dilepton mass m_{ee} [GeV]'
     123set ylabel 'd{/Symbol s}/dm_{ee} [{/Symbol m}b/GeV]'
     125set logscale y
     126set xrange [0:1.1]
     127set mxtics 2
     128set mytics 10
     129set yrange [5e-4:5e1]
     130set format y "10^{%L}"
     132set key samplen 3 box opaque
     134unset bars
     136set style data lines
     138FF = "DileptonMass.dat"
     140set output 'pp35_mass.eps'
     142plot \
     143     FF u 1:($4+$5+$6+$7+$8+$13+colsum(21,39)) ls 1 t t2, \
     144     FF u 1:4 ls 3 t t4, \
     145     FF u 1:5 ls 4 t t5, \
     146     FF u 1:6 ls 5 t t6, \
     147     FF u 1:7 ls 6 t t7, \
     148     FF u 1:8 ls 7 t t8 , \
     149     FF u 1:9 ls 8 t "{/Symbol D} QED", \
     150     FF u 1:(colsum(21,33)) ls 2 t "N* VMD", \
     151     FF u 1:(colsum(34,39)) ls 22 t "{/Symbol D}* VMD", \
     152     FF u 1:14 ls 9 t t14
     154!epstopdf pp35_mass.eps