# Simulations

The EIC task force has a large number of simulation tools available for investigating different types of physics processes. Unless noted otherwise, these can be accessed from /afs/rhic.bnl.gov/eic/PACKAGES.

## Event Generators

The following event generators are currently available:

- ep
- DJANGOH: (un)polarised DIS generator with QED and QCD radiative effects for NC and CC events.
- MILOU: A generator for deeply virtual Compton scattering (DVCS), the Bethe-Heitler process and their interference.
- PYTHIA: A general-purpose high energy physics event generator.
- PEPSI: A generator for polarised leptoproduction.
- RAPGAP: A generator for deeply inelastic scattering (DIS) and diffractive
*e + p*events.

- eA
- BeAGLE: Benchmark eA Generator for LEptoproduction - a generator to simulate ep/eA DIS events including nuclear shadowing effects (based on DPMJetHybrid)
- eSTARlight: Monte Carlo that simulates coherent vector meson photo- and electro-production in electron-ion collisions.

These generators are not part of the current distribution but can be installed upon request:

- ep

- eA
- DPMJet: a generator for very low Q2/real photon physics in eA
- DPMJetHybrid: an obsolete buggy beta version of BeAGLE: a generator to simulate ep/eA DIS events by employing PYTHIA in DPMJet. Not supported in any way. No physics should be published with it.
- Sartre is an event generator for exclusive diffractive vector meson production and DVCS in ep and eA collisions based on the dipole model.

There is code provided to convert the output from most of these generators into a ROOT format. It is distributed as part of eic-smear, the Monte Carlo smearing package.

## Detector simulations

The following programs are available for simulating detector geometry and response:

- eic-smear A package to apply very fast detector smearing to Monte Carlo events.

More details on detector simulations can be found here.

## Manuals

See the pages of the programmes listed above for their documentation. Other useful references are:

- BASES/SPRING v1 and v5.1: Cross section integration and Monte Carlo event generation. Used in Rapgap and MILOU.

## Helpful/Important Links

The following pages provide useful general information for Monte Carlo simulations:

- MC programs:
- A list of Monte Carlo programmes
- HepForge, high-energy physics development environment, which includes many Monte Carlo generators.
- Lecture slides from a course on QCD and Monte Carlos

- Radiative Correction Codes:
- A discussion of Radiative corrections

- Parton Distribution Function Interfaces:
- LHAPDF, the Les Houches Accord PDF Interface.
- By default, programs are linked to version 5.9.1, pointed to by
**$LHAPDF5**. - You can adapt your
**LD_LIBRARY_PATH**to point to**$LHAPDF6**instead. This is a more modern version but lacks support for some EIC-relevant PDFs. - The PDF-Grids can be accessed via
**/gpfs01/eic/data/LHAPDF59SHARE/lhapdf/PDFsets**, or more generally via**/cvmfs/sft.cern.ch/lcg/external/lhapdfsets**. These should be found automatically if using the system-default version of LHAPDF.

- By default, programs are linked to version 5.9.1, pointed to by
- The users' manual of the CERN PDFLIB

- LHAPDF, the Les Houches Accord PDF Interface.

## MC Analysis Techniques

##### How to get a cross section (for non-BeAGLE MC)

to normalize your counts to cross section you need two pieces of information

- the total number of trials, it is printed to the screen/logfile if all our MC finish
- the total integrated cross section, the unit is in general microbarn, it is printed to the screen/logfile if all our MC finish

**Counts = Luminosity x Cross Section**

==> count * total integrated cross section /total number of trials

to calculate the corresponding MC luminosity

==> total number of trials/ total integrated cross section

There are some handy ROOT functions available to get the total number of trials, the total integrated MC cross section and the total number of events in the Tree

These work on Pythia, Pepsi, Djangoh and Milou event-wise root trees

- total number of trials:

TObjString* nEventsString(NULL) file.GetObject("nEvents", nEventsString);

- total integrated MC cross section

TObjString* crossSectionString(NULL); file.GetObject("crossSection", crossSectionString);

- total number of events in the tree:

TTree* tree(NULL); file.GetObject("EICTree", tree);

##### How to get a cross-section for BeAGLE

Note: This refers to standalone BeAGLE using internal e+N event generation (Pythia or optionally RAPGAP (when available)). For BeAGLE+GCF see below.

When *BeAGLE* simulates electron-nucleus (eA) collisions, it tags each event as scattering either from a proton or from a neutron in the nucleus. The *stdout* (which can be redirected to a log file) contains the total electron-proton cross section, *σ _{ep}*, and the total electron-neutron cross section,

*σ*, for the beam energies and kinematic ranges over which we choose to generate events. As long as we simulate enough events so that the cross sections can be effectively sampled (using sequential mode), these two values (

_{en}*σ*and

_{ep}*σ*) will be independent of the number of events generated.

_{en}It should be noted that as long as *σ _{ep}* and

*σ*are equal to the desired precision, which is often the case at low x (high energy), then the simpler "non-BeAGLE" cross-section method above can be used with the exception that "number of events" should be used instead of "number of trials". For the most precise approach, we must correct for the fact that BeAGLE assumes equal ep and en cross-sections when selecting whether the struck nucleon is a proton or neutron, which is not quite right.

_{en}For a given electron-nucleon integrated luminosity, *lum ^{eic}*, the predicted number of events is

num^{eic}= [ (Z×σ_{ep}+ N×σ_{en}) / A ] × lum^{eic}

Here, for the nucleus being simulated, *A = Z + N*, *Z*, and *N* are the number of nucleons, protons, and neutrons, respectively. From the predicted total number of events, *num ^{eic}*, and the individual cross sections, we can calculate the ratio of the number originating from an electron-proton scattering to those originating from an electron-neutron scattering:

num^{eic}= num_{prot}+ num_{neut}= [ (Z×σ_{ep}+ N×σ_{en}) / A ] × lum^{eic}= (Z/A)×σ_{ep}×lum^{eic}+ (N/A)×σ_{en}×lum^{eic}

This gives

num_{prot}= (Z/A)×σ_{ep}×lum^{eic}num_{neut}= (N/A)×σ_{en}×lum^{eic}

From where we calculate the following ratio:

num_{prot}/ num_{neut}= (σ_{ep}/ σ_{en}) × (Z/N)

This means that if we simulate a certain number of events, *num ^{sim}*, the ratio of the simulated events which should come from a proton scattering to those which should come from a neutron scattering is given by the above equation. However, since the cross-section ratio depends strongly on the kinematic range we choose to generate over, using the true ratio is impractical. In practice,

*BeAGLE*simulates proton and neutron scattering events according to the ratio of protons to neutrons in the nucleus under consideration.

So if we generate a total number of events, *num ^{sim}*, and then apply cuts (and/or binning), we end up with a final number of events,

*num*. Then we can ask, for a given electron-nucleon integrated luminosity, lum

^{sim,cut}^{eic}, how many events are predicted with that cut,

*num*. In order calculate this accurately, we first define the following:

^{eic,cut}num^{sim}= num^{sim}_{prot}+ num^{sim}_{neut}

lum^{sim}_{prot}= num^{sim}_{prot}/ σ_{ep}lum^{sim}_{neut}= num^{sim}_{neut}/ σ_{en}

In the above equations, *num ^{sim}_{prot}* and

*num*are the total number of proton and neutron scattering events generated, respectively. These numbers can be found in the output log file, or in the output ROOT file by seeing how many events are tagged as coming from a nucleon==2212 (proton) or a nucleon==2112 (neutron).

^{sim}_{neut}Next, after applying the cuts (and/or binning) mentioned above, we can write

num^{sim,cut}= num^{sim,cut}_{prot}+ num^{sim,cut}_{neut}

Again, using the output ROOT file, we can determine the number of proton-tagged events and neutron-tagged events which pass the cut. Then we can calculate the expected number of events:

num^{eic,cut}= lum^{eic}× ( [ Z×(num^{sim,cut}_{prot}/lum^{sim}_{prot}) + N×(num^{sim,cut}_{neut}/lum^{sim}_{neut}) ] / A )

As a specific example, we can consider the case where the 'cut' is simply a range in *x* and *Q ^{2}* (i.e. a

*x-Q*bin). We use the above result to calculate the inclusive (per-nucleon) differential cross section as

^{2}d^{2}σ/dxdQ^{2}= num^{eic,cut}/ ( lum^{eic}× Δx × ΔQ^{2}) = ( [ Z×(num^{sim,cut}_{prot}/lum^{sim}_{prot}) + N×(num^{sim,cut}_{neut}/lum^{sim}_{neut}) ] / A ) / (Δx × ΔQ^{2})

In the above equation, *Δx × ΔQ ^{2}* gives the area of the

*x-Q*bin.

^{2}##### How to get a cross-section for BeAGLE+GCF

GCF (Generalized Contact Formalism) is a package which calculates cross-sections in e+A collisions with Short-Range Correlations (SRCs) for a variety of specific processes: Quasi-elastic, vector-meson diffraction, DIS. In many cases, BeAGLE can read these events in and handle the nuclear response. GCF generates events which are flat or gaussian in a variety of variables and then reports a weight which can be used to convert this to a cross-section.

The weight is carried in the "RAevt" variable and has units *nb* (some older versions used *cm ^{2}*).

##### How to scale to the MC luminosity to the luminosity we want for the measurement

Very often it is impossible to generate so many events that the MC luminosity would correspond to one month of EIC running.

For this case we generate so much MC events that all distributions are smooth and scale the uncertainties.

The factor needed to scale is the ratio **lumi-scale-factor = EIC-luminosity / generated MC luminosity**. If we have this factor there are 2 ways to scale.

- scaling of counts in histogram by

h11->Scale(lumi-scale-factor);

this will scale the number of counts in each bin of the histogram to what you would get for the EIC-luminosity

statistical uncertainties can then be calculated simply by sqrt(counts)

- scaling the statistical uncertainties only

sqrt(counts)/sqrt(lumi-scale-factor)

##### Example: reduced cross section

This example shows how to calculate the reduced cross section need to extract F_2 and F_L and how to scale the statistical uncertainties to a certain integrated luminosity

sigma_reduced = prefactor * dsigma/dx/dQ2 with prefactor = Q^4 * x / (2*pi*alpha_em^2*(1+(1-y)^2)

this cross section would have the unit barn * GeV^2, to make it dimensionless you need to use a conversion factor for barn to 1/GeV^2 (h^2c^2/GeV^2 = 0.3894 millibarn)

sigma_reduced = counts(x,Q^2) * prefactor * total integrated MC cross section /total number of trials/ conversion-barn-to-GeV /x-binsize/Q2-binsize if the root function Scale was used the statistical uncertainty is delta sigma_reduced = sqrt(counts(x,Q^2)) * prefactor * total integrated MC cross section /total number of trials/ conversion-barn-to-GeV /x-binsize/Q2-binsize in the other case it is delta sigma_reduced = sqrt(counts(x,Q^2)) * prefactor * total integrated MC cross section /total number of trials/ conversion-barn-to-GeV /x-binsize/Q2-binsize/ sqrt(lumi-scale-factor)

**Attention:** all luminosities and cross section must be in the same unit (pb or fb or ...)

##### High-Statistics BeAGLE Simulation

Multiple BeAGLE simulation jobs can submitted to the BNL batch farm using the condor utility. A simple submission file to run 10 jobs simultaneously is shown below.

Universe = vanilla Notification = Never Executable = run_eAu.sh GetEnv = True Input = /dev/null Arguments = "$(Process)" Output = jobout/eAujob.out.$(Process) Error = jobout/eAujob.err.$(Process) Log = jobout/eAujob.log.$(Process) Queue 10

This now works fine. Random # seeds are now taken from /dev/urandom inside of BeAGLE automatically. If that is unavailable a more complicated calculation is done which includes the time AND the process id. There is no longer any problem with identical seeds on parallel farm jobs occurring at the same time.

**THE REST OF THIS SECTION IS OBSOLETE**

If using the standard BeAGLE random number seeding, however, many of the seeds can be identical if the jobs start at the same second. An example is given below. Running this command on the output log files:

grep -i "SEED = " eAu_*.log

displays the random number seeds:

eAu_0.log: SEED = 2411 eAu_1.log: SEED = 2421 eAu_2.log: SEED = 2296 eAu_3.log: SEED = 2296 eAu_4.log: SEED = 2421 eAu_5.log: SEED = 2136 eAu_6.log: SEED = 2321 eAu_7.log: SEED = 2136 eAu_8.log: SEED = 2391 eAu_9.log: SEED = 2391

As can be seen, several of the jobs have the same seed, and therefore the output of the simulation will be identical. A simple workaround is to use the "/dev/urandom" pseudo-random number generator in a shell script to set the seed directly in the BeAGLE input file. This is done in the following shell script, which is the executable in the condor submission file.

#!/usr/bin/bash

if [ -z "$1" ] then echo "No job number set." echo "Please run as ./run_eAu.sh jobnumber" echo "Exiting..." exit 1 fi

echo "-----------------------------------" echo "Running BeAGLE Simulation for eAu Collider!!!" echo "-----------------------------------" echo "Performing Job $1" echo "..." echo ""

VAR1=$1 PythiaCard="s/eAu.txt/eAu_${VAR1}.txt/g" BeAGLECard="s/S3ALL002/InpAu_${VAR1}/g"

sed "${PythiaCard}" S3ALL002 > InpAu_$1 sed "${BeAGLECard}" ./inputFiles/eAu.inp > ./inputFiles/eAu_$1.inp

SEED=`od -vAn -N2 -tu2 < /dev/urandom` echo "The Random SEED is ${SEED// /}" echo "" sed -i "s/1234567/${SEED// /}/g" InpAu_$1

echo "Running BeAGLE..." $BEAGLESYS/BeAGLE < inputFiles/eAu_$1.inp > logs/eAu_$1.log

echo "Completed Simulation!!!" echo ""

echo "Making Output ROOT File..." root -l -b -q "make_tree.C(\"eAu_$1.txt\")" echo "Done!!!" echo ""

echo "Cleaning up..." rm -vf ./outForPythiaMode/eAu_$1.txt rm -vf InpAu_$1 rm -vf inputFiles/eAu_$1.inp echo "Done!!!" echo ""

This script makes use of this line in the Pythia input card that is called by BeAGLE:

FSEED 1234567

Using the above script, the seeds for the simultaneously run jobs should differ. An example of the set seeds for 10 jobs is shown below.

eAu_0.log: SEED = 14473 eAu_1.log: SEED = 30550 eAu_2.log: SEED = 55545 eAu_3.log: SEED = 4022 eAu_4.log: SEED = 28238 eAu_5.log: SEED = 36662 eAu_6.log: SEED = 12502 eAu_7.log: SEED = 57292 eAu_8.log: SEED = 29721 eAu_9.log: SEED = 37451

As can be seen, all the random seeds are different. A larger test was performed where 100 jobs were submitted simultaneously, each job simulating 100k events (10 million events in total). The jobs were all found to have different seeds, and all the jobs completed within 5 hours of submission.