diff --git a/dis/README.md b/dis/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..ee752d7c97951bd9e5ca27ae57f0883368b50013
--- /dev/null
+++ b/dis/README.md
@@ -0,0 +1,12 @@
+# DIS Benchmarks
+
+## Compiling Pythia
+
+```
+g++ src/pythia_dis.cc -o pythia_dis \
+ -I/home/whit/stow/hepmc3/include \
+ -I../include -O2 -std=c++11 -pedantic -W -Wall -Wshadow -fPIC \
+ -L../lib -Wl,-rpath,../lib -lpythia8 -ldl \
+ -L/home/whit/stow/hepmc3/lib -Wl,-rpath,/home/whit/stow/hepmc3/lib -lHepMC3
+```
+
diff --git a/dis/dis.sh b/dis/dis.sh
index 171c792a4e6388b9942bd4df01b43114a2d6193f..977f331270988fa415132a45c5e92601ca0b28e9 100644
--- a/dis/dis.sh
+++ b/dis/dis.sh
@@ -46,6 +46,13 @@ export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
## TODO remove this
#bash util/build_detector.sh
+g++ dis/src/pythia_dis.cc -o pythia_dis \
+ -I/usr/local/include \
+ -O2 -std=c++11 -pedantic -W -Wall -Wshadow -fPIC \
+ -L/usr/local/lib -Wl,-rpath,/usr/local/lib -lpythia8 -ldl \
+ -L/usr/local/lib -Wl,-rpath,/usr/local/lib -lHepMC3
+
+
## =============================================================================
## Step 2: Run the simulation
echo "Running Geant4 simulation"
diff --git a/dis/src/pythia_dis.cc b/dis/src/pythia_dis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..71f92413c1f2790b1343487e8598e2963720d1a0
--- /dev/null
+++ b/dis/src/pythia_dis.cc
@@ -0,0 +1,189 @@
+// main45.cc is a part of the PYTHIA event generator.
+// Copyright (C) 2020 Torbjorn Sjostrand.
+// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
+// Please respect the MCnet Guidelines, see GUIDELINES for details.
+
+// Author: Stefan Prestel <stefan.prestel@thep.lu.se>.
+
+// Keywords: LHE file; hepmc;
+
+// This program (main45.cc) illustrates how a file with HepMC3 events can be
+// generated by Pythia8. See main44.cc for how to ouput HepMC2 events instead.
+// Note: both main44.cc and main45.cc can use the same main44.cmnd input card.
+
+#include "Pythia8/Pythia.h"
+#include "Pythia8Plugins/HepMC3.h"
+#include <unistd.h>
+
+using namespace Pythia8;
+
+//==========================================================================
+
+// Example main programm to illustrate merging.
+
+int main( int argc, char* argv[] ){
+
+ // Check that correct number of command-line arguments
+ if (argc != 2) {
+ cerr << " Unexpected number of command-line arguments ("<<argc<<"). \n"
+ << " You are expected to provide the arguments" << endl
+ << " 1. Output file for HepMC events" << endl
+ << " Program stopped. " << endl;
+ return 1;
+ }
+
+ // Beam energies, minimal Q2, number of events to generate.
+ double eProton = 250.;
+ double eElectron = 10.0;
+ double Q2min = 5.;
+ int nEvent = 10000;
+
+ // Generator. Shorthand for event.
+ Pythia pythia;
+ Event& event = pythia.event;
+
+ // Set up incoming beams, for frame with unequal beam energies.
+ pythia.readString("Beams:frameType = 2");
+ // BeamA = proton.
+ pythia.readString("Beams:idA = 2212");
+ pythia.settings.parm("Beams:eA", eProton);
+ // BeamB = electron.
+ pythia.readString("Beams:idB = 11");
+ pythia.settings.parm("Beams:eB", eElectron);
+
+ // Set up DIS process within some phase space.
+ // Neutral current (with gamma/Z interference).
+ pythia.readString("WeakBosonExchange:ff2ff(t:gmZ) = on");
+ // Uncomment to allow charged current.
+ //pythia.readString("WeakBosonExchange:ff2ff(t:W) = on");
+ // Phase-space cut: minimal Q2 of process.
+ pythia.settings.parm("PhaseSpace:Q2Min", Q2min);
+
+ // Set dipole recoil on. Necessary for DIS + shower.
+ pythia.readString("SpaceShower:dipoleRecoil = on");
+
+ // Allow emissions up to the kinematical limit,
+ // since rate known to match well to matrix elements everywhere.
+ pythia.readString("SpaceShower:pTmaxMatch = 2");
+
+ // QED radiation off lepton not handled yet by the new procedure.
+ pythia.readString("PDF:lepton = off");
+ pythia.readString("TimeShower:QEDshowerByL = off");
+
+ // Initialize.
+ pythia.init();
+
+ // Interface for conversion from Pythia8::Event to HepMC one.
+ HepMC3::Pythia8ToHepMC3 toHepMC;
+ // Specify file where HepMC events will be stored.
+ HepMC3::WriterAscii ascii_io(argv[1]);
+ cout << endl << endl << endl;
+
+ // Histograms.
+ double Wmax = sqrt(4.* eProton * eElectron);
+ Hist Qhist("Q [GeV]", 100, 0., 50.);
+ Hist Whist("W [GeV]", 100, 0., Wmax);
+ Hist xhist("x", 100, 0., 1.);
+ Hist yhist("y", 100, 0., 1.);
+ Hist pTehist("pT of scattered electron [GeV]", 100, 0., 50.);
+ Hist pTrhist("pT of radiated parton [GeV]", 100, 0., 50.);
+ Hist pTdhist("ratio pT_parton/pT_electron", 100, 0., 5.);
+
+ double sigmaTotal(0.), errorTotal(0.);
+ bool wroteRunInfo = false;
+ // Get the inclusive x-section by summing over all process x-sections.
+ double xs = 0.;
+ for (int i=0; i < pythia.info.nProcessesLHEF(); ++i)
+ xs += pythia.info.sigmaLHEF(i);
+
+ // Begin event loop.
+ for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
+ if (!pythia.next()) continue;
+
+ double sigmaSample = 0., errorSample = 0.;
+
+ // Four-momenta of proton, electron, virtual photon/Z^0/W^+-.
+ Vec4 pProton = event[1].p();
+ Vec4 peIn = event[4].p();
+ Vec4 peOut = event[6].p();
+ Vec4 pPhoton = peIn - peOut;
+
+ // Q2, W2, Bjorken x, y.
+ double Q2 = - pPhoton.m2Calc();
+ double W2 = (pProton + pPhoton).m2Calc();
+ double x = Q2 / (2. * pProton * pPhoton);
+ double y = (pProton * pPhoton) / (pProton * peIn);
+
+ // Fill kinematics histograms.
+ Qhist.fill( sqrt(Q2) );
+ Whist.fill( sqrt(W2) );
+ xhist.fill( x );
+ yhist.fill( y );
+ pTehist.fill( event[6].pT() );
+
+ // pT spectrum of partons being radiated in shower.
+ for (int i = 0; i < event.size(); ++i) if (event[i].statusAbs() == 43) {
+ pTrhist.fill( event[i].pT() );
+ pTdhist.fill( event[i].pT() / event[6].pT() );
+ }
+
+ // Get event weight(s).
+ double evtweight = pythia.info.weight();
+
+ // Do not print zero-weight events.
+ if ( evtweight == 0. ) continue;
+
+ // Create a GenRunInfo object with the necessary weight names and write
+ // them to the HepMC3 file only once.
+ if (!wroteRunInfo) {
+ shared_ptr<HepMC3::GenRunInfo> genRunInfo;
+ genRunInfo = make_shared<HepMC3::GenRunInfo>();
+ vector<string> weight_names = pythia.info.weightNameVector();
+ genRunInfo->set_weight_names(weight_names);
+ ascii_io.set_run_info(genRunInfo);
+ ascii_io.write_run_info();
+ wroteRunInfo = true;
+ }
+
+
+ // Construct new empty HepMC event.
+ HepMC3::GenEvent hepmcevt;
+
+ // Work with weighted (LHA strategy=-4) events.
+ double normhepmc = 1.;
+ if (abs(pythia.info.lhaStrategy()) == 4)
+ normhepmc = 1. / double(1e9*nEvent);
+ // Work with unweighted events.
+ else
+ normhepmc = xs / double(1e9*nEvent);
+
+ // Set event weight
+ //hepmcevt.weights().push_back(evtweight*normhepmc);
+ // Fill HepMC event
+ toHepMC.fill_next_event( pythia, &hepmcevt );
+ // Add the weight of the current event to the cross section.
+ sigmaTotal += evtweight*normhepmc;
+ sigmaSample += evtweight*normhepmc;
+ errorTotal += pow2(evtweight*normhepmc);
+ errorSample += pow2(evtweight*normhepmc);
+ // Report cross section to hepmc
+ shared_ptr<HepMC3::GenCrossSection> xsec;
+ xsec = make_shared<HepMC3::GenCrossSection>();
+ // First add object to event, then set cross section. This order ensures
+ // that the lengths of the cross section and the weight vector agree.
+ hepmcevt.set_cross_section( xsec );
+ xsec->set_cross_section( sigmaTotal*1e9, pythia.info.sigmaErr()*1e9 );
+ // Write the HepMC event to file. Done with it.
+ ascii_io.write_event(hepmcevt);
+
+
+
+ // End of event loop. Statistics and histograms.
+ }
+ pythia.stat();
+ cout << Qhist << Whist << xhist << yhist << pTehist << pTrhist << pTdhist;
+
+
+ return 0;
+
+}