From 6f23bcd58bf1428d1c3ad9ba3efc8ab6244f815e Mon Sep 17 00:00:00 2001
From: Sylvester Joosten <sylvester.joosten@gmail.com>
Date: Wed, 6 Apr 2016 18:13:58 -0400
Subject: [PATCH] wrote actual pythia program

---
 program/example-config.json |  50 +++--
 program/pythia6m.cc         | 404 +++++++++++++++++++++++++++++++++++-
 2 files changed, 425 insertions(+), 29 deletions(-)

diff --git a/program/example-config.json b/program/example-config.json
index c1efd5a..a09304d 100644
--- a/program/example-config.json
+++ b/program/example-config.json
@@ -1,24 +1,32 @@
-"mc" : {
-  "run" : 1,
-  "events": 1000,
-  "target": {
-    "A": 1,
-    "Z": 1, 
-    "nucleon": "p"
+{
+  "mc" : {
+    "type": "MC Data",
+    "run" : 1,
+    "events": 1000,
+    "target": {
+      "A": 1,
+      "Z": 1, 
+      "nucleon": "p"
+    },
+    "beam": {
+      "energy": 11
+    },
+    "Q2": {
+      "min": 0.9,
+      "max": 60
+    },
+    "x": {
+      "min": 0.002,
+      "max": 0.99
+    },
+    "y": {
+      "min": 0.05,
+      "max": 0.95 
+    }
   },
-  "beam": {
-    "energy": 11
-  },
-  "Q2": {
-    "min": 0.9,
-    "max": 60
-  },
-  "x": {
-    "min": 0.002,
-    "max": 0.99
-  },
-  "y": {
-    "min": 0.05,
-    "max": 0.95 
+  "output": {
+    "type": "MC Output",
+    "root": true,
+    "text": true
   }
 }
diff --git a/program/pythia6m.cc b/program/pythia6m.cc
index fabfcd0..32e0f1f 100644
--- a/program/pythia6m.cc
+++ b/program/pythia6m.cc
@@ -1,16 +1,404 @@
 #include <pythia6m/util/configuration.hh>
+#include <pythia6m/util/exception.hh>
 #include <pythia6m/util/framework.hh>
+#include <pythia6m/util/logger.hh>
+#include <pythia6m/util/assert.hh>
+
+#include <pythia6m/gmc/gen_set.h>
+#include <pythia6m/gmc/gmc.h>
+#include <pythia6m/gmc/mc_set.h>
+#include <pythia6m/gmc/mcevnt.h>
+#include <pythia6m/gmc/random.h>
+
+#include <pythia6m/fmotion/fmotion.h>
+
+#include <pythia6m/pythia6/pythia6.h>
+
+#include <cstring>
+#include <fstream>
+#include <memory>
+
+#include <TFile.h>
+#include <TLorentzVector.h>
+#include <TTree.h>
 
 using namespace pythia6m;
 
-int run_pythia(const ptree& settings, const std::string& output) {}
+GeneratorSettings glb_gen;
+MCSettings glb_mc;
+EventData glb_event;
+
+// global for output configuration for ...sigh... consistency
+struct output_cfg {
+  TFile* root;
+  std::unique_ptr<std::ofstream> text;
+};
+output_cfg glb_out;
+
+// global event buffer...
+struct part_buffer {
+  int32_t n_part;
+  std::vector<int32_t> index;
+  std::vector<int32_t> ctype;
+  std::vector<int32_t> ltype;
+  std::vector<int32_t> parent;
+  std::vector<int32_t> daughter; // this is the first daughter, or zero
+  std::vector<double> mom;
+  std::vector<double> px;
+  std::vector<double> py;
+  std::vector<double> pz;
+  std::vector<double> E;
+  std::vector<double> mass;
+  std::vector<double> vx;
+  std::vector<double> vy;
+  std::vector<double> vz;
+  std::vector<bool> init;
+  std::vector<bool> scat;
+  std::vector<bool> virt;
+  std::vector<bool> lund;
+
+  void clear() {
+    n_part = 0;
+    index.clear();
+    ctype.clear();
+    ltype.clear();
+    parent.clear();
+    daughter.clear();
+    mom.clear();
+    px.clear();
+    py.clear();
+    pz.clear();
+    E.clear();
+    mass.clear();
+    vx.clear();
+    vy.clear();
+    vz.clear();
+    init.clear();
+    scat.clear();
+    virt.clear();
+    lund.clear();
+  }
+
+  void add_all() {
+    for (int idx = 0; idx < glb_pyJets.n; ++idx) {
+      add(idx);
+    }
+  }
+
+  void add(int32_t idx) {
+    n_part += 1;
+    index.push_back(idx);
+    ctype.push_back(glb_pyJets.K[0][idx]);
+    ltype.push_back(glb_pyJets.K[1][idx]);
+    parent.push_back(glb_pyJets.K[2][idx]);
+    daughter.push_back(glb_pyJets.K[3][idx]);
+    TLorentzVector ptmp{glb_pyJets.P[0][idx], glb_pyJets.P[1][idx],
+                        glb_pyJets.P[2][idx], glb_pyJets.P[3][idx]};
+    mom.push_back(ptmp.Vect().Mag());
+    px.push_back(ptmp.X());
+    py.push_back(ptmp.Y());
+    pz.push_back(ptmp.Z());
+    E.push_back(ptmp.E());
+    mass.push_back({glb_pyJets.P[4][idx]});
+    vx.push_back(glb_pyJets.V[0][idx]);
+    vy.push_back(glb_pyJets.V[1][idx]);
+    vz.push_back(glb_pyJets.V[2][idx]);
+    int ctype = glb_pyJets.K[0][idx];
+    if (glb_event.scatIdx == idx) {
+      init.push_back(false);
+      scat.push_back(true);
+      virt.push_back(false);
+      lund.push_back(false);
+    } else if (ctype == 21) {
+      init.push_back(true);
+      scat.push_back(false);
+      virt.push_back(false);
+      lund.push_back(false);
+    } else if (ctype >= 1 && ctype <= 10 && ctype != 2 && ctype != 3) {
+      init.push_back(false);
+      scat.push_back(false);
+      virt.push_back(false);
+      lund.push_back(true);
+    } else {
+      init.push_back(false);
+      scat.push_back(false);
+      virt.push_back(true);
+      lund.push_back(false);
+    }
+  }
+};
+part_buffer glb_parts;
+
+void conf_generator(const ptree& settings) {
+  LOG_INFO("conf", "Configuring generator");
+
+  // User control switch MSEL (manual section 9.2, also 8.3.1)
+  // MSEL=2 additionally enables the elastic and diffactive components of the
+  // VMD and GVMD parts
+  glb_gen.PhotoType = 2;
+
+  // Let LUND decay some of the semi-stable particles
+  glb_gen.DecayPiK = 0;
+  glb_gen.DecayLamKs = 1;
+
+  // RADGEN options, currently disabled
+  glb_gen.radgenGenLUT = 0;
+  glb_gen.radgenUseLUT = 0;
+  sprintf(glb_gen.GenRad, "NO\n");
+
+  // activate Fermi momentum? (experimental)
+  glb_gen.enableFM = 1;
+  glb_gen.FMNBins = 1000; // Number of bins for FM discretization
+  glb_gen.FMCutOff = 1.5; // Max FM (in GeV)
+
+  // Parameterizations for F2Pythia and R
+  // (I believe these are only relevant when running with RADGEN)
+  sprintf(glb_gen.FStruct, "F2PY");
+  sprintf(glb_gen.R, "1990");
+
+  // PYTHIA model to run with:
+  //  * RAD for realistic DIS
+  //  * GAM for real photon beam
+  sprintf(glb_gen.PyModel, "RAD");
+}
+void conf_mc(const ptree& settings) {
+  LOG_INFO("conf", "Configuring MC");
+  configuration conf{settings, "mc"};
+
+  glb_mc.RunNo = conf.get<int>("run");
+  glb_mc.NEvents = conf.get<int>("events");
+  glb_mc.TarA = conf.get<int>("target/A");
+  glb_mc.TarA = conf.get<int>("target/Z");
+  glb_mc.TargetNucleon[0] = conf.get<char>("target/nucleon");
+  glb_mc.EneBeam = conf.get<double>("beam/energy");
+  glb_mc.Q2Min = conf.get<double>("Q2/min");
+  glb_mc.Q2Max = conf.get<double>("Q2/max");
+  glb_mc.xMin = conf.get<double>("x/min");
+  glb_mc.xMax = conf.get<double>("x/max");
+  glb_mc.yMin = conf.get<double>("y/min");
+  glb_mc.yMax = conf.get<double>("y/max");
+
+  // hardcoded electron beam for now
+  sprintf(glb_mc.BeamParType, "ELEC");
+  // hardcoded Beam and target Polarization
+  // needed because RADGEN can also run with polarized generators
+  glb_mc.PBValue = 0;
+  glb_mc.PTValue = 0;
+  sprintf(glb_mc.PBeam, "L\n");
+  sprintf(glb_mc.PTarget, "L\n");
+
+  // reconsider fermi motion when running for A=1 options
+  if (glb_mc.TarA == 1) {
+    glb_gen.enableFM = 0;
+  }
+}
+// init the RNG
+void init_random() {
+  LOG_INFO("init",
+           "initializing the random generator with seed (run number): " +
+               std::to_string(glb_mc.RunNo));
+  int dummy1{0};
+  int dummy2{0};
+  const char* command{" "};
+  RNDMQ(dummy1, dummy2, glb_mc.RunNo, command);
+}
+// init the event structure
+void init_event() {
+  LOG_INFO("init", "Initializing the event common blocks");
+  // beam particle
+  if (!strncmp(glb_mc.BeamParType, "POSI", 4)) {
+    glb_event.beamLType = -11;
+  } else if (!strncmp(glb_mc.BeamParType, "ELEC", 4)) {
+    glb_event.beamLType = 11;
+  } else if (!strncmp(glb_mc.BeamParType, "GAM", 3)) {
+    glb_event.beamLType = 22;
+  } else {
+    throw exception("Unrecognized beam particle. This should NEVER HAPPEN!");
+  }
+
+  // Init FM engine
+  if (glb_gen.enableFM) {
+    INIT_FM(glb_mc.TarZ, glb_mc.TarA, glb_gen.FMCutOff, glb_gen.FMNBins);
+  }
+
+  glb_event.evNum = 0;
+  glb_event.nEvGen = 0;
+  glb_event.nEvGenPrev = 0;
+  glb_event.nErr = 0;
+
+  glb_event.weight = 1.f;
+  glb_event.Q2Gen = 0.f;
+  glb_event.nuGen = 0.f;
+  glb_event.xGen = 0.f;
+  glb_event.yGen = 0.f;
+  glb_event.W2Gen = 0.f;
+  glb_event.thetaGen = 0.f;
+  glb_event.phiGen = 0.f;
+  glb_event.EGen = 0.f;
+  glb_event.pGen = 0.f;
+  glb_event.pxGen = 0.f;
+  glb_event.pyGen = 0.f;
+  glb_event.pzGen = 0.f;
+  glb_event.vxGen = 0.f;
+  glb_event.vyGen = 0.f;
+  glb_event.vzGen = 0.f;
+}
+void init_radgen() {
+  // TODO
+}
+void init_gmc() {
+  LOG_INFO("init", "Initializing GMC/PYTHIA");
+  int ok;
+  GMC_INIT(ok);
+  if (!ok) {
+    throw exception("Failed to initialized GMC/PYTHIA");
+  }
+}
+void conf_output(const ptree& settings, const std::string& output) {
+  LOG_INFO("conf", "Configuring output");
+  configuration conf{settings, "output"};
+  if (conf.get<bool>("root")) {
+    glb_out.root =
+        new TFile{(output + ".root").c_str(), "recreate"};
+  }
+  if (conf.get<bool>("text")) {
+    glb_out.text = std::make_unique<std::ofstream>((output + ".lund").c_str());
+  }
+  tassert(glb_out.root || glb_out.text, "Need at least one output format");
+}
+
+void clear_event_kinematics() {
+  glb_event.scatIdx = -1;
+  glb_event.Q2Gen = 0.f;
+  glb_event.nuGen = 0.f;
+  glb_event.xGen = 0.f;
+  glb_event.yGen = 0.f;
+  glb_event.W2Gen = 0.f;
+  glb_event.pxGen = 0.f;
+  glb_event.pyGen = 0.f;
+  glb_event.pzGen = 0.f;
+  glb_event.vxGen = 0.f;
+  glb_event.vyGen = 0.f;
+  glb_event.vzGen = 0.f;
+}
+
+void events() {
+  TTree* tree;
+  if (glb_out.root) {
+    glb_out.root->cd();
+    tree = new TTree("pythia6m", "pythia6m");
+    tree->Branch("event", &glb_event.evNum);
+    tree->Branch("evgen", &glb_event.nEvGen);
+    tree->Branch("process", &glb_event.process);
+    tree->Branch("Q2", &glb_event.Q2Gen);
+    tree->Branch("nu", &glb_event.nuGen);
+    tree->Branch("x", &glb_event.xGen);
+    tree->Branch("y", &glb_event.yGen);
+    tree->Branch("W2", &glb_event.W2Gen);
+    tree->Branch("n_part", &glb_parts.n_part);
+    tree->Branch("index", &glb_parts.index);
+    tree->Branch("ctype", &glb_parts.ctype);
+    tree->Branch("ltype", &glb_parts.ltype);
+    tree->Branch("parent", &glb_parts.parent);
+    tree->Branch("daughter", &glb_parts.daughter);
+    tree->Branch("mom", &glb_parts.mom);
+    tree->Branch("px", &glb_parts.px);
+    tree->Branch("py", &glb_parts.py);
+    tree->Branch("pz", &glb_parts.pz);
+    tree->Branch("E", &glb_parts.E);
+    tree->Branch("mass", &glb_parts.mass);
+    tree->Branch("vx", &glb_parts.vx);
+    tree->Branch("vy", &glb_parts.vy);
+    tree->Branch("vz", &glb_parts.vz);
+    tree->Branch("init", &glb_parts.init);
+    tree->Branch("scat", &glb_parts.scat);
+    tree->Branch("virt", &glb_parts.virt);
+    tree->Branch("lund", &glb_parts.lund);
+  }
+
+  // event loop
+  for (int iev{0}; iev < glb_mc.NEvents; ++iev) {
+    // clear the event kinematics
+    clear_event_kinematics();
+    glb_parts.clear();
+
+    // generate one event
+    int nEvGen{0};
+    GMC_EVENT(nEvGen);
+
+    // update counters
+    glb_event.nEvGen += nEvGen;
+
+    // insert good event check here if needed
+
+    // we now have a good event
+    glb_event.nEvGenPrev = glb_event.nEvGen;
+    glb_event.evNum = iev;
+
+    // talk to the human
+    if (!((iev + 1) % 5000)) {
+      LOG_INFO(__func__, "Generated Events: " +
+                             std::to_string(glb_event.nEvGen) +
+                             ", Cross Section: " +
+                             std::to_string(glb_pyPars.pari[0]) +
+                             ", Good Events: " +
+                             std::to_string(iev + 1));
+    }
 
-int main(int argc, char* argv[]) {
-  try {
-    framework pythia6m{argc, argv, run_pythia};
-    pythia6m.run();
-    return 0;
-  } catch (...) {
-    return 1;
+    // get all the particles
+    glb_parts.add_all();
+
+    // write output
+    if (glb_out.root && tree) {
+      tree->Fill();
+    }
+    if (glb_out.text) {
+      char buf[1024];
+      int nparts = 0;
+      for (int ii = 0; ii < glb_parts.n_part; ++ii) {
+        if (glb_parts.lund[ii] || glb_parts.scat[ii]) {
+          nparts += 1;
+        }
+      }
+      sprintf(buf, "%4i %4i %4i %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n",
+              nparts, glb_mc.TarA, glb_mc.TarZ, 0., 0., glb_event.xGen,
+              glb_event.yGen, glb_event.W2Gen, glb_event.Q2Gen,
+              glb_event.nuGen);
+      *glb_out.text << buf;
+      for (int ii = 0; ii < glb_parts.n_part; ++ii) {
+        if (glb_parts.lund[ii] || glb_parts.scat[ii]) {
+          sprintf(buf,
+                  "%4i %4i %4i %4i %4i %4i %12.6e %12.6e %12.6e %12.6e %12.6e "
+                  "%12.6e %12.6e %12.6e\n",
+                  ii, 0, 1, glb_parts.ltype[ii], glb_parts.parent[ii],
+                  glb_parts.daughter[ii], glb_parts.px[ii], glb_parts.py[ii],
+                  glb_parts.pz[ii], glb_parts.E[ii], glb_parts.mass[ii],
+                  glb_parts.vx[ii], glb_parts.vy[ii], glb_parts.vz[ii]);
+          *glb_out.text << buf;
+        }
+      }
+    }
+  }
+  if (glb_out.root) {
+    tree->Write();
+    glb_out.root->Close();
   }
 }
+
+int run_pythia(const ptree& settings, const std::string& output) {
+  LOG_INFO("main", "PYTHIA6M");
+  // configuration
+  conf_generator(settings);
+  conf_mc(settings);
+  conf_output(settings, output);
+  // initialization
+  init_random();
+  init_event();
+  init_gmc();
+
+  events();
+
+  return 0;
+}
+
+MAKE_PYTHIA6M_FRAMEWORK(run_pythia)
-- 
GitLab