Add proper SIDIS event generation with Pythia8

.clang-format

 ---
-BasedOnStyle:  LLVM
-BreakConstructorInitializersBeforeComma: true
-ConstructorInitializerAllOnOneLineOrOnePerLine: true
+Language:        Cpp
+BasedOnStyle:  Chromium
+AccessModifierOffset: -2
+AlignAfterOpenBracket: Align
+AlignConsecutiveAssignments: true
+AlignConsecutiveDeclarations: true
+AlignEscapedNewlines: Right
+AlignOperands:   true
+AlignTrailingComments: true
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortBlocksOnASingleLine: false
+AllowShortCaseLabelsOnASingleLine: false
+AllowShortFunctionsOnASingleLine: All
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: false
+AlwaysBreakAfterDefinitionReturnType: None
+AlwaysBreakAfterReturnType: None
+AlwaysBreakBeforeMultilineStrings: false
+AlwaysBreakTemplateDeclarations: true
+BinPackArguments: true
+BinPackParameters: true
+BraceWrapping:   
+  AfterClass:      false
+  AfterControlStatement: false
+  AfterEnum:       false
+  AfterFunction:   true
+  AfterNamespace:  false
+  AfterObjCDeclaration: false
+  AfterStruct:     false
+  AfterUnion:      false
+  AfterExternBlock: false
+  BeforeCatch:     false
+  BeforeElse:      false
+  IndentBraces:    false
+  SplitEmptyFunction: true
+  SplitEmptyRecord: true
+  SplitEmptyNamespace: true
+BreakBeforeBinaryOperators: None
+BreakBeforeBraces: Custom
+BreakBeforeInheritanceComma: false
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+BreakConstructorInitializers: BeforeColon
+BreakAfterJavaFieldAnnotations: false
+BreakStringLiterals: true
+ColumnLimit:     120
+CommentPragmas:  '^ IWYU pragma:'
+CompactNamespaces: false
+ConstructorInitializerAllOnOneLineOrOnePerLine: false
+ConstructorInitializerIndentWidth: 4
+ContinuationIndentWidth: 4
 Cpp11BracedListStyle: true
-Standard: Cpp11
-#SpaceBeforeParens: ControlStatements
-SpaceAfterControlStatementKeyword: true
-PointerBindsToType: true
+DerivePointerAlignment: false
+DisableFormat:   false
+ExperimentalAutoDetectBinPacking: false
+FixNamespaceComments: true
+ForEachMacros:   
+  - foreach
+  - Q_FOREACH
+  - BOOST_FOREACH
+IncludeBlocks:   Preserve
+IncludeCategories: 
+  - Regex:           '^"(llvm|llvm-c|clang|clang-c)/'
+    Priority:        2
+  - Regex:           '^(<|"(gtest|gmock|isl|json)/)'
+    Priority:        3
+  - Regex:           '.*'
+    Priority:        1
+IncludeIsMainRegex: '(Test)?$'
+IndentCaseLabels: false
+IndentPPDirectives: None
+IndentWidth:     2
+IndentWrappedFunctionNames: false
+JavaScriptQuotes: Leave
+JavaScriptWrapImports: true
+KeepEmptyLinesAtTheStartOfBlocks: true
+MacroBlockBegin: ''
+MacroBlockEnd:   ''
+MaxEmptyLinesToKeep: 1
+NamespaceIndentation: All
+ObjCBlockIndentWidth: 2
+ObjCSpaceAfterProperty: false
+ObjCSpaceBeforeProtocolList: true
+PenaltyBreakAssignment: 2
+PenaltyBreakBeforeFirstCallParameter: 19
+PenaltyBreakComment: 300
+PenaltyBreakFirstLessLess: 120
+PenaltyBreakString: 1000
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 60
+PointerAlignment: Left
+ReflowComments:  true
+SortIncludes:    true
+SortUsingDeclarations: true
+SpaceAfterCStyleCast: false
+SpaceAfterTemplateKeyword: true
+SpaceBeforeAssignmentOperators: true
+SpaceBeforeParens: ControlStatements
+#SpaceBeforeRangeBasedForLoopColon: true
+SpaceInEmptyParentheses: false
+SpacesBeforeTrailingComments: 1
+SpacesInAngles:  false
+SpacesInContainerLiterals: true
+SpacesInCStyleCastParentheses: false
+SpacesInParentheses: false
+SpacesInSquareBrackets: false
+Standard:        Cpp11
+TabWidth:        8
+UseTab:          Never
 ...
+

dis/README.md

+# 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
+```
+

dis/config.yml

 dis:generate:
   stage: initialize
+  image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:latest
   needs: []
   timeout: 1 hours
   artifacts:

dis/dis.sh

@@ -46,6 +46,7 @@ export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 ## TODO remove this
 #bash util/build_detector.sh
 
+
 ## =============================================================================
 ## Step 2: Run the simulation
 echo "Running Geant4 simulation"

dis/gen.sh

@@ -35,16 +35,30 @@ JUGGLER_FILE_NAME_TAG="dis"
 ## =============================================================================
 ## Step 1: Dummy event generator
 ## TODO better file name that encodes the actual configuration we're running
-root -b -q "dis/generator/gen_central_electrons.cxx(${JUGGLER_N_EVENTS}, \".local/${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+echo "Compiling   dis/src/pythia_dis.cc ..."
+g++ dis/src/pythia_dis.cc -o pythia_dis  \
+   -I/usr/local/include  -Iinclude \
+   -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
 if [[ "$?" -ne "0" ]] ; then
-  echo "ERROR running script"
+  echo "ERROR compiling pythia"
+  exit 1
+fi
+echo "done"
+pwd
+ls -lrth
+
+./pythia_dis dis.hepmc
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running pythia"
   exit 1
 fi
 
 ## =============================================================================
 ## Step 2: finalize
 echo "Moving event generator output into ${DATA_PATH}"
-mv .local/${JUGGLER_FILE_NAME_TAG}.hepmc ${DATA_PATH}/${JUGGLER_FILE_NAME_TAG}.hepmc
+#mv .local/${JUGGLER_FILE_NAME_TAG}.hepmc ${DATA_PATH}/${JUGGLER_FILE_NAME_TAG}.hepmc
 
 ## =============================================================================
 ## All done!

dis/src/pythia_dis.cc

+#include "Pythia8/Pythia.h"
+#include "Pythia8Plugins/HepMC3.h"
+#include <unistd.h>
+#include <cassert>
+
+using namespace Pythia8;
+
+#include "clipp.h"
+using namespace clipp;
+using std::string;
+//______________________________________________________________________________
+
+enum class mode { none, help, list, part };
+
+struct settings {
+  double      E_electron   = 10.0;  // GeV
+  double      E_ion        = 275.0; // GeV
+  std::string outfile      = "dis.hepmc";
+  int         ion_PID      = 2212;
+  int         electron_PID = 11;
+  bool        help         = false;
+  bool        success      = false;
+  double      Q2_min       = 4.0;
+  int         N_events     = 1000;
+  mode    selected       = mode::none;
+};
+
+template <typename T>
+void print_usage(T cli, const char* argv0)
+{
+  // used default formatting
+  std::cout << "Usage:\n" << usage_lines(cli, argv0) << "\nOptions:\n" << documentation(cli) << '\n';
+}
+//______________________________________________________________________________
+
+template <typename T>
+void print_man_page(T cli, const char* argv0)
+{
+  // all formatting options (with their default values)
+  auto fmt = clipp::doc_formatting{}
+                 .start_column(8)                    // column where usage lines and documentation starts
+                 .doc_column(20)                     // parameter docstring start col
+                 .indent_size(4)                     // indent of documentation lines for children of a
+                                                     // documented group
+                 .line_spacing(0)                    // number of empty lines after single documentation lines
+                 .paragraph_spacing(1)               // number of empty lines before and after paragraphs
+                 .flag_separator(", ")               // between flags of the same parameter
+                 .param_separator(" ")               // between parameters
+                 .group_separator(" ")               // between groups (in usage)
+                 .alternative_param_separator("|")   // between alternative flags
+                 .alternative_group_separator(" | ") // between alternative groups
+                 .surround_group("(", ")")           // surround groups with these
+                 .surround_alternatives("(", ")")    // surround group of alternatives with these
+                 .surround_alternative_flags("", "") // surround alternative flags with these
+                 .surround_joinable("(", ")")        // surround group of joinable flags with these
+                 .surround_optional("[", "]")        // surround optional parameters with these
+                 .surround_repeat("", "...");        // surround repeatable parameters with these
+
+  auto mp = make_man_page(cli, argv0, fmt);
+  mp.prepend_section("DESCRIPTION", "Simple pythia dis generator.");
+  mp.append_section("EXAMPLES", " $ pythia_dis [output_file]");
+  std::cout << mp << "\n";
+}
+//______________________________________________________________________________
+
+settings cmdline_settings(int argc, char* argv[]) {
+  settings s;
+  auto     lastOpt =
+      " options:" % (option("-h", "--help").set(s.selected, mode::help) % "show help",
+                     value("file", s.outfile).if_missing([] {
+                       std::cout << "You need to provide an output filename argument!\n";
+                     }) % "output file");
+
+  auto cli =
+      (command("help").set(s.selected, mode::help) | lastOpt);
+
+  assert(cli.flags_are_prefix_free());
+
+  auto res = parse(argc, argv, cli);
+
+  // if( res.any_error() ) {
+  //  s.success = false;
+  //  std::cout << make_man_page(cli, argv[0]).prepend_section("error: ",
+  //                                                           "    The best
+  //                                                           thing since
+  //                                                           sliced bread.");
+  //  return s;
+  //}
+  s.success = true;
+
+  if (s.selected == mode::help) {
+    print_man_page<decltype(cli)>(cli, argv[0]);
+  };
+  return s;
+}
+//______________________________________________________________________________
+
+int main(int argc, char* argv[]) {
+
+  settings s = cmdline_settings(argc, argv);
+  if (!s.success) {
+    return 1;
+  }
+
+  if (s.selected == mode::help) {
+    return 0;
+  }
+
+  // Beam energies, minimal Q2, number of events to generate.
+  double eProton   = s.E_ion;
+  double eElectron = s.E_electron;
+  double Q2min     = s.Q2_min;
+  int    nEvent    = s.N_events;
+
+  // 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 = " + std::to_string(s.ion_PID));
+  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;
+}
+