# Difference between revisions of "BeAGLE"

BeAGLE - Benchmark eA Generator for LEptoproduction

BeAGLE, a Fortran program designed as a general purpose eA Monte-Carlo generator, was based on DPMJetHybrid, but has been upgraded to improve the nuclear handling:

1. Allowing a hard collision to occur on a neutron as well as a proton, thereby also conserving charge. In DPMJetHybrid, when a neutron was struck it was replaced with the contents of a Pythia ep collision.
2. Reinstating the Fermi momentum of the struck nucleon (lost in DPMJetHybrid), thereby conserving 4-momentum.
3. Having the target rest frame be the same for the struck nucleon and the nucleus on average
4. Optionally, we also allow for multiple nucleons to participate in the collision in the kinematic (x, Q^2 ) region where nuclear shadowing occurs.

## Program Overview

BeAGLE uses the DPMJet framework in order to handle the target (Glauber, FLUKA etc.) for ep and eA collisions, but employs PYTHIA to handle the elementary interaction as a replacement for PHOJET. The nuclear PDF is included by linking PYTHIA PDF library to LHAPDF which enables us access EPS09 nuclear PDF routines. We have included the energy loss module PyQM, written by Accardi and Dupre based on the quenching weights of Salgado & Wiedemann. This module, when turned on, applies energy loss to the partons after they have been simulated by Pythia, but before they have been hadronized.

The original DPMJet models DTUNUC, PHOJET, QNEUTRIN are NOT available in BeAGLE. If you are interested in these models for pp/pA, real photoproduction, or quasi-elastic neutrino scattering, use DPMJet directly. The history and motivation of the use of DPMJetHybrid rather than DPMJet for eRHIC can be found here.

A graphical demonstration of the program structure can be illustrated as follows:

 generator concept.

## Basic Model

The program can be run in, essentially, two modes.

Setting genShd=1 (generator Shadowing=1) is similar to DPMJetHybrid. The entire nuclear shadowing effect (where R=σ(A)/Aσ(p)<1) is interpreted as a change in the parton distribution of unspecified origin, but the DIS interaction is assumed to be point-like. One and only one nucleon participates in the interaction.

Note: This approach is in stark contrast to the original DPMJET approach which used GMVD (Generalized Vector Meson Dominance) and allowed multiple collisions.

Setting genShd>1 (specifically 2 or 3) activates a very simplistic model of nuclear shadowing loosely based on the dipole-model or aligned-jet model. Even for DIS events, the scattering process in the target rest frame is viewed as being due to a rare hadronic fluctuation of the virtual photon with a very large cross-section (few mb). In this case, the entire R(x,Q2)<1 shadowing effect is attributed to multiple nucleon collisions. Viewed in the hadronic center-of-mass frame, such collisions can still look like conventional DIS.

We implement this by making a map, outside of BeAGLE, using TGlauberMC between the "dipole" cross-section σ and the value of R. When running BeAGLE, we look up R in the nPDF (such as EPS09) and find the corresponding dipole cross-section. Using the dipole cross-section, we can then simulate the collisions with Glauber, leading to a variety of Ncoll values. If Ncoll>1, we currently pick one of the struck nucleons to undergo a Pythia interaction (for genShd=2, this is a random choice among the struck nucleons; for genShd=3 it is the first struck nucleon seen by the γ*). Any additional struck nucleons undergo an elastic scatter with the most forward parton (in the direction of the γ*) with a pt magnitude of a typical intrinsic kt.

It should be noted that for xBj>0.1 or R>1, we revert to genShd=1 behavior.

At the moment, there are many simplifying assumptions for the multinucleon mode, made for convenience, which we hope to systematically relax as we upgrade the program. These include:

   Infinite coherence length even at modest x
We apply an overall multinucleon effect based on the Rsea(x,Q2) independent of process
Note: Processes such as Photon-Gluon Fusion should depend of RG while diffraction should INCREASE as R gets smaller.
We assume that the full shadowing effect is due to multinucleon collisions with no modification of individual nucleons.


## How to run this code

### Developers:

the code can be found on the afs area for EIC at BNL

/afs/rhic.bnl.gov/eic/PACKAGES/BeAGLE/


The main program in the file called user-eA3.0-5.f. A simple README is kept in the "Documents" directory, with a short description about what files are contained in the directory.

A Makefile script has already been provided with this package. FLUKA and LHAPDF have to be preconfigured before compiling this program. The RAPGAP library situation is complicated. You have to make a RAPGAP library and then remove the duplicate routine SFECFE which is already in FLUKA.

The code can also be found in the local BNL git repository: /afs/rhic.bnl.gov/eic/gitstore/BeAGLE/. Make a local clone for yourself if you want to modify the code. The code can also be run by pointing $BEAGLESYS at this directory instead of the .../PACKAGES/BeAGLE directory. Results should be identical. At JLAB, /u/group/ldgeom/PACKAGES/BeAGLE is a git repository and also a valid$BEAGLESYS release.

### Users:

The code can be run directly from the afs area, but for general users, it is recommended to run from your private directory.

Setup (done once): Go to the (data) directory where you want to run BeAGLE from and:

   source /afs/rhic/EIC/PACKAGES/BeAGLE/setup


This will create a subdirectory called outForPythiaMode where your output .txt files will go. It will create a symbolic (soft) link to the BeAGLE exectuable It will make a copy of the FLUKA input data file nuclear.bin.

Setup (every session) in your .cshrc or when you want to run:

   source /afs/rhic/EIC/PACKAGES/BeAGLE/setup_BeAGLE


This will create an environment variable $BEAGLESYS needed by BeAGLE, pointing to the afs BeAGLE area. If$BEAGLESYS is not defined, BeAGLE will look in the working directory (.) for a variety of files and crash if they aren't there.

To run: ./BeAGLE < yourinputfile.inp > yourlogfile.log

Example input files can be found in /afs/rhic/EIC/PACKAGES/BeAGLE/

Note: This input file (originally from DPMJET) is very particular about the spacing of a lot of the numbers so make sure that you leave the starting space the same when you change numbers.

Frequently changed cards:

    START    - the first # is the # of events requested.
TARPAR   - the first # is A and the second # is Z for the target
the third # is the n/p handling mode for the Pythia subevent:
0=sequential n,p. The first ~N/A events are en, the remaining ep (binomial prob.)
Useful for getting en & ep cross-sections from Pythia
1=en only test mode (not really for physics)
2=ep only test mode (not really for physics)
3=random mix. The events are randomly en or ep.
Useful because you can analyze just a subset of the data without bias.
TAUFOR   - the first # is the formation time parameter (tau) in fm/c for the intranuclear cascade
the second # is the number of generations followed Default=25, 0 means no cascade
MOMENTUM - the first # is for the lepton beam (GeV/c), the second for the A (or p)
both #s should be entered as positive, but the lepton beam will be multiplied by -1.
L-TAG    - these #s are cuts: yMin yMax Q2Min Q2Max thetaMin thetaMax
where y & Q2 (GeV2) are the leptoproduction kinematics
and theta (radians) refers to the lepton scattering angle in the laboratory frame.
PY-INPUT - specifies the file (with an eight-character maximum name!) used as a pythia input file.
Examples (such as eAS2noq) can be found in the /afs/rhic/EIC/PACKAGES/BeAGLE/ directory.
FERMI    - First # means: -1 = no Fermi motion at all
1 = DPMJET Fermi motion, but Pythia subevent neglects it (DPMJetHybrid mode)
2 = Fermi motion added to Pythia subevent after the fact
3 = Pythia subevent used correct Fermi motion (Not Yet Implemented)
Second #: "Scale factor" for Fermi momentum distribution in GeV
Third #: Fermi momentum distribution 0 (D) = DPMJET distribution
Fourth #: Post-processing flag:  (Not Yet Implemented)
0 (D) = no post-processing
1     = Fix energy non-conservation in ion frame (for small nuclei)
OUTLEVEL  - First 4 #s are verbosity flags -1=quiet, >=1 increasing verbosity
Fifth # is the number of events to print out and in some cases to be verbose about
PYVECTORS - Allowed Pythia vector mesons for diffraction: 0(D)=all, 1=rho, 2=omega, 3=phi, 4=J/psi
USERSET   = First # specifies the meaning of the variables User1,User2,User3 (detailed below)
Second # specifies the maximum excitation energy in GeV handed to Fluka (D=9.0)
Note: This should not come into play, but it protects against infinite loops.
The last two cards should be:
START     = First # specifies the # of events to run
Second # should be 0 (or missing)
STOP


The format of the PY-INPUT file is similar to the Pythia input file, which is not really documented elsewhere, but is fairly self-explanatory due to the comments.

    line 1:    output file name for the textfile
line 2:    xmin and xmax
line 3:    ignored unless radcor is being used
line 4: 0 switch for radcor. 0=off. 1=on. 2=generate tables. radcor>0 NOT YET TESTED.
line 5: 1 GVMD model choice 0=x,Q2, 1=y,Q2
line 6: 1,1  A-Tar and Z-Tar really only used if radcor is used.
line 7:   Shadowing/multinucleon switch genShd for BeAGLE. 1=1 nucleon struck.
2=Glauber. Pythia-interacting nucleon chosen at random from struck nucleons. All others elastic scatter.
3=like 2, but Pythia-interaction is the first (most negative z in TRF with +z along γ* direction) nucleon.
line 8: 101  Means LO with error set 1 in the nuclear PDF.
line 9:      switch for quenching (1:on, 0:off) NOTE: quenching is NOT WELL TESTED yet.
line 10:     Qhat value for quenching (0.0 also means no quenching even if the switch=1)
OPTIONAL LINE: FSEED ####### - Fixes seed for random # generator for debugging. Don't use for production.


other lines make changes to PYTHIA control parameters.

The usage and explanations for the parameters in the input file and what kind of parameters do we need to give for PYTHIA can be found in the following file:

/afs/rhic/eic/PACKAGES/BeAGLE/Documents/Documentation


To find the final state particles, just use the following selection rules:

ISTHKK(i)=1


which selects all the particles from hard interaction, nuclear evaporation and heavy nuclear fragments.

## Output Data Format

The BeAGLE output is a text file intended to be read in by eic-smear and eicroot at BNL and by GEMC at JLAB.

For the output file, we propose the following text format & structure:

• 1st line: BeAGLE EVENT FILE
• 2nd line: "============================================"
• 3rd line: Information on event wise variables stored in the file
 I: I 0 (line index) ievent: I eventnumber running from 1 to XXX genevent: I trials to generate this event ltype: I particle code for incoming lepton Atarg: I A of the target nucleus Ztarg: I Z of the target nucleus pzlep: D lab pz of the lepton beam (sign(pz)*p for non-zero crossing angle) pztarg: D lab pz/A of the ion beam (sign(pz)*p/A for non-zero crossing angle) pznucl: D lab pz of the struck nucleon (sign(pz)*p for non-zero crossing angle) crang: D crossing-angle (mr). crang = (1000*) atan(px/pz) for the beam momentum with nonzero px. Note: we assume one of the beams has px=py=0 and the other py=0. crori: I crossing angle orientation. Failed to parse (Missing texvc executable. Please see math/README to configure.): \pm 1 lepton beam defines Failed to parse (Missing texvc executable. Please see math/README to configure.): \pm\hat{z} direction. Failed to parse (Missing texvc executable. Please see math/README to configure.): \pm 2  hadron/ion beam defines Failed to parse (Missing texvc executable. Please see math/README to configure.): \pm\hat{z} direction. 0 means no crossing angle.  subprocess: I pythia subprocess (MSTI(1)), for details see PYTHIA nucleon: I hadron beam type (MSTI(12)) targetparton: I parton hit in the target (MSTI(16)) xtargparton: D x of target parton (PARI(34)) beamparton: I in case of resolved photon processes and soft VMD the virtual photon has a hadronic structure. This gives the info which parton interacted with the target parton (MSTI(15)) xbeamparton: D x of beam parton (PARI(33)) thetabeamparton: D theta of beam parton (PARI(53)) truey, trueQ2, truex, trueW2, trueNu: D are the kinematic variables of the event. If radiative corrections are turned on they are different from what is calculated from the scattered lepton. If radiative corrections are turned off they are the same as what is calculated from the scattered lepton leptonphi: D phi of the lepton (VINT(313)) in the laboratory frame s_hat: D s-hat of the process (PARI(14)) t_hat: D Mandelstam t (PARI(15)) u_hat: D Mandelstam u (PARI(16)) pt2_hat: D pT2-hat of the hard scattering (PARI(18)) Q2_hat: D Q2-hat of the hard scattering (PARI(22)), F2, F1, R, sigma_rad, SigRadCor: D information used and needed in the radiative correction code EBrems: D energy of the radiative photon in the nuclear rest frame photonflux: D flux factor from PYTHIA (VINT(319)) b: D Impact parameter (in fm), radial position of virtual photon w.r.t. the center of the nucleus in the nuclear TRF with z along γ* Phib: D azimuthal angle of the impact parameter of the virtual photon in the nuclear TRF with z along γ* and φ=0 for the scattered electron Thickness: D Nuclear thickness T(b)/ρ0 in units of fm. ThickScl: D Nuclear thickness T(b) in nucleons/fm2. Ncollt: I Number of collisions between the incoming γ* and nucleons in the nucleus (same as number of participating nucleons) Ncolli: I Number of inelastic collisions between the incoming γ* and nucleons in the nucleus (same as number of inelastically participating nucleons) Nwound: I Number of wounded nucleons, including those involved in the intranuclear cascade. WARNING: This variable is not finalized yet Nwdch: I Number of wounded protons, including those involved in the intranuclear cascade. WARNING: This variable is not finalized yet Nnevap: I Number of evaporation (and nuclear breakup) neutrons Npevap: I Number of evaporation (and nuclear breakup) protons Aremn: I A of the nuclear remnant after evaporation and breakup Ninc: I Number of stable hadrons from the Intra Nuclear Cascade Nincch: I Number of charged stable hadrons from the Intra Nuclear Cascade d1st: D density-weighted distance from first collision to the edge of the nucleus (amount of material traversed / ρ0) davg: D Average density-weighted distance from all inelastic collisions to the edge of the nucleus pxf,pyf,pzf: D Fermi-momentum of the struck nucleon (or sum Fermi momentum for all inelastic nucleon participants for Ncolli>1) in target rest frame with z along gamma* direction Eexc: D Excitation energy in the nuclear remnant before evaporation and breakup. RA D Nuclear PDF ratio for the up sea for the given event kinematics (x,Q2), but set to 1 if multi-nucleon shadowing is off (genShd=1) User1,User2,User3: D User variables to prevent/delay future format changes nrTracks: I number of tracks in this event, including event history
• 4th line: "============================================"

• 5th line: Information on track wise variables stored in the file
 I: I line index, runs from 1 to nrTracks ISTHKK(I): I status code KS: KS=1 is the only stable final state particle code, Use NoBAM variable (below) to specify origin of particle IDHKK(I): I particle KF code (211: pion, 2112:n, ....). Code 80000 refers to a nucleus, specified in more detail by A=IDRES(I), Z=IDXRES(I) JMOHKK(2,I): I line number of second mother particle JMOHKK(1,I): I line number of first mother particle JDAHKK(1,I): I normally the line number of the first daughter. JDAHKK(2,I): I normally the line number of the last daughter. PHKK(1,I): D px of particle (GeV/c) PHKK(2,I): D py of particle (GeV/c) PHKK(3,I): D pz of particle (GeV/c) PHKK(4,I): D Energy of particle (GeV) PHKK(5,I): D mass of particle (GeV/c^2) VHKK(1,I): D x vertex information (mm) VHKK(2,I): D y vertex information (mm) VHKK(3,I): D z vertex information (mm) IDRES(I): I Baryon number, or A for a nucleus (IDHKK(I)=80000), fractional B set to 0. IDXRES(I): I Particle charge, (Z for a nucleus), 0 for fractional charge. NOBAM(I): I Flag describing the particle origin, particularly for final state particles.
• 6th line: "============================================"
• 7th line: event information for first event. Format descriptor:
    format ((I4,1x,$),(I10,1x,$),4(I4,1x,$),4(f12.6,1x,$),3(I4,1x,$), & I6,1x,$,f9.6,1x,$,I6,1x,$,2(f12.6,1x,$),7(f18.11,3x,$),
&     11(f19.9,3x,$),4(f10.6,1x,$),9(I5,1x,$),2(f10.6,1x,$),
&     3(f15.6,1x,$),2(f9.6,1x,$),3(e17.8,1x,$),I6,/)  • 8th line: "============================================" • 9th to X-1 line: trackwise info for first event. Format descriptor:  format (2(I6,1X,$),I10,1x,$,4(I8,1x,$),5(f15.6,1x,$),3(e15.6,1x,$),3(I8,1x,\$)/)

• Xth line "=============== Event finished ==============="

the information from line 7 to X repeats for each event.
Something to be mentioned here is the first 4 tracks in every event which deliver the information of beam parameter are added externally in order to get some handy variables like z, ptVsGamma... So the mother daughter relationship of these four tracks can not be taken too seriously. But the orig number and daughter number of the tracks after those 4 tracks are consistent and can be used for further purposes.

The 4-momentum for particle I, PHKK(mu,I) is reported in the lab frame for all stable final state particles. It is also reported in the lab frame for most intermediate/document particles with the following exceptions: Particle 5 is a copy of the γ* with 4-momentum reported as zero. Particles 6 through A+5 are the incoming original nucleons from the A, in the target rest frame with γ* momentum along +z. The momentum corresponds to the Fermi momentum of the nucleon in question.

The impact parameter, defined as the transverse position of the virtual photon with respect to the nucleus center in a frame with z along γ*, is given by bx=VHKK(1,5) and by=VHKK(2,5) in mm.

The status codes ISTHKK(I) have the following meaning. Note that final state particles have status 1 (although internally -1 and 1001 are also used). All other statuses refer to intermediate particles or informational lines.

 (-1): Internally used for final state particle from excited nucleus evaporation and breakup. Output as KS=1,NOBAM=3 1: stable final state particle, Note: internally, -1 and 1001 are also used. 2: unstable intermediate particle 3: Documentation line describing the collision 11: Projectile (γ*) documentation line 12: Incoming wounded nucleon from nucleus 14: Incoming spectator nucleon from nucleus 16: Intermediate baryon involved in the intranuclear cascade 18: Spectator nucleon (wrt hard collision) involved in the intranuclear cascade 19: Intermediated meson involved in the intranuclear cascade 21: documentation line describing the collision 1000: excited nucleus, which will evaporate and breakup (1001): Internally used to mean residual nucleus. Output as KS=1, NOBAM=4

The NOBAM particle origin flag for final state particles has the following meaning:

 0,2: Created in the hard collision process. 3: Created during the evaporation/breakup phase (formerly specified by KS=-1) 4: Largest remnant nucleus created during the evaporation/breakup phase (formerly specified by KS=1001). 11-34: Particles created during the Intra-nuclear cascade. Cascade generation (IDCH) #n specified by NOBAM=10+n