Generate events with D0 decay into pion and kaon

This steering macro is used to produce enriched samples of neutral-current DIS events with D0 that decay into pion and kaon within the ePIC acceptance of |eta| < 3.5. The minimum Q^2 of the DIS events are 1 and 100 GeV^2

config_10x100_Q2min1.cmnd.txt

@@ -0,0 +1,16 @@
+// Pythia Configuration File
+
+// Beam Settings
+Beams:idA = 2212
+Beams:idB = 11
+Beams:frameType = 2
+Beams:eA = 100
+Beams:eB = 10
+
+// Process Settings
+WeakBosonExchange:ff2ff(t:gmZ) = on
+SpaceShower:dipoleRecoil = on
+PDF:lepton = off
+
+// Phase Space Settings
+PhaseSpace:Q2min = 1

config_10x100_Q2min100.cmnd.txt

@@ -0,0 +1,16 @@
+// Pythia Configuration File
+
+// Beam Settings
+Beams:idA = 2212
+Beams:idB = 11
+Beams:frameType = 2
+Beams:eA = 100
+Beams:eB = 10
+
+// Process Settings
+WeakBosonExchange:ff2ff(t:gmZ) = on
+SpaceShower:dipoleRecoil = on
+PDF:lepton = off
+
+// Phase Space Settings
+PhaseSpace:Q2min = 100

main41_D0_from_piK.cc

@@ -0,0 +1,179 @@
+// main41.cc is a part of the PYTHIA event generator.
+// Copyright (C) 2021 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.
+
+
+// Authors: Mikhail Kirsanov <[email protected]>.
+
+
+// Keywords: basic usage; hepmc;
+
+
+// This program illustrates how HepMC can be interfaced to Pythia8.
+// It studies the charged multiplicity distribution at the LHC.
+// HepMC events are output to the hepmcout41.dat file.
+
+
+// WARNING: typically one needs 25 MB/100 events at the LHC.
+// Therefore large event samples may be impractical.
+
+
+//#include "Pythia8Plugins/ROOT.h"
+
+
+#include "Pythia8/Pythia.h"
+#ifndef HEPMC2
+#include "Pythia8Plugins/HepMC3.h"
+#else
+#include "Pythia8Plugins/HepMC2.h"
+#endif
+
+
+#include <chrono>
+#include <iostream>
+#include <cstdio>  // For std::rename
+
+
+using namespace Pythia8;
+using namespace std::chrono;
+
+
+bool areVectorsDifferentIgnoreOrder(vector<int> vec1, vector<int> vec2) {
+  //if (vec1.size() != vec2.size()) return true;
+  sort(vec1.begin(), vec1.end());
+  sort(vec2.begin(), vec2.end());
+  return vec1 != vec2;
+}
+
+
+int main() {
+
+
+  // Start measuring time
+  auto start = high_resolution_clock::now();
+
+
+  // Interface for conversion from Pythia8::Event to HepMC
+  // event. Specify file where HepMC events will be stored.
+  Pythia8::Pythia8ToHepMC FirstFile("output/hepmc_D01.dat");
+  Pythia8::Pythia8ToHepMC SecondFile("output/hepmc_D02.dat");
+  Pythia8::Pythia8ToHepMC ThirdFile("output/hepmc_D03.dat");
+  Pythia8::Pythia8ToHepMC FourthFile("output/hepmc_D04.dat");
+  Pythia8::Pythia8ToHepMC FifthFile("output/hepmc_D05.dat");
+  Pythia8::Pythia8ToHepMC SixthFile("output/hepmc_D06.dat");
+  Pythia8::Pythia8ToHepMC SeventhFile("output/hepmc_D07.dat");
+  Pythia8::Pythia8ToHepMC EighthFile("output/hepmc_D08.dat");
+  Pythia8::Pythia8ToHepMC NinthFile("output/hepmc_D09.dat");
+  Pythia8::Pythia8ToHepMC TenthFile("output/hepmc_D10.dat");
+
+
+  // Generator. Process selection. LHC initialization. Histogram.
+  Pythia pythia;
+  pythia.readFile("config_10x100_Q2min1.cmnd");
+  pythia.readString("Print:quiet = on"); // Set Pythia to quiet mode
+  pythia.init();
+
+
+  // Begin event loop. Generate event. Skip if error.
+  int Mothercount = 0;
+  int EventsCount = 0;
+  const int Total_Count = 1000000;
+
+
+  while (Mothercount < Total_Count) {
+    if (!pythia.next()) continue;
+    Pythia8::Event &event = pythia.event;
+
+
+    bool hasMother = false;
+    EventsCount++;
+
+
+    vector<int> daughterPDGs = {211, 321};
+    const int MotherPDG = 421;
+    float EtaMax = 3.5;
+
+
+    for (int i = 0; i < event.size(); i++) {
+      Pythia8::Particle &currentParticle = event[i];
+      if (abs(currentParticle.id()) != MotherPDG) continue;
+      vector<int> daughters = currentParticle.daughterList();
+      if (daughters.size() != daughterPDGs.size()) continue;
+      vector<int> daughtersID = {abs(pythia.event[daughters[0]].id()), abs(pythia.event[daughters[1]].id())};
+      if (areVectorsDifferentIgnoreOrder(daughtersID, daughterPDGs)) continue;
+      bool DaughterOutOfRange = false;
+      for (int &daughter : daughters) {
+        if (fabs(pythia.event[daughter].eta()) > EtaMax) {
+          DaughterOutOfRange = true;
+          break;
+        }
+      }
+      if (DaughterOutOfRange) continue;
+      Mothercount++;
+      hasMother = true;
+    }
+
+
+    // Construct new empty HepMC event, fill it and write it out.
+    if (hasMother) {
+      if (Mothercount % 10 == 0) FirstFile.writeNextEvent(pythia);
+      if (Mothercount % 10 == 1) SecondFile.writeNextEvent(pythia);
+      if (Mothercount % 10 == 2) ThirdFile.writeNextEvent(pythia);
+      if (Mothercount % 10 == 3) FourthFile.writeNextEvent(pythia);
+      if (Mothercount % 10 == 4) FifthFile.writeNextEvent(pythia);
+      if (Mothercount % 10 == 5) SixthFile.writeNextEvent(pythia);
+      if (Mothercount % 10 == 6) SeventhFile.writeNextEvent(pythia);
+      if (Mothercount % 10 == 7) EighthFile.writeNextEvent(pythia);
+      if (Mothercount % 10 == 8) NinthFile.writeNextEvent(pythia);
+      if (Mothercount % 10 == 9) TenthFile.writeNextEvent(pythia);
+    }
+
+
+    if (EventsCount % 1770130 == 0) cout << "Mother count = " << Mothercount << endl;
+  }
+
+
+  // End of event loop. Statistics. Histogram.
+  pythia.stat();
+  cout << "Total events needed = " << EventsCount << endl;
+
+
+  // Stop measuring time and calculate the elapsed time
+  auto stop = high_resolution_clock::now();
+  auto duration = duration_cast<seconds>(stop - start);
+  int hours = duration.count() / 3600;
+  int minutes = (duration.count() % 3600) / 60;
+  int seconds = duration.count() % 60;
+
+
+  // Print the elapsed time
+  std::cout << "Time taken by function: "
+            << hours << " hours, "
+            << minutes << " minutes, "
+            << seconds << " seconds" << std::endl;
+
+
+  // Rename the output files
+  for (int i = 1; i <= 10; ++i) {
+    char oldFilename[50];
+    snprintf(oldFilename, sizeof(oldFilename), "output/hepmc_D%02d.dat", i);
+
+
+    char newFilename[256];
+    snprintf(newFilename, sizeof(newFilename),
+             "output/hepmc_file%02d_D0_from_piK_PYTHIA8_ep_NC_10x100_Q2min1_EtaMax3p5_%dh%dm%ds_%dEvents.dat",
+             i, hours, minutes, seconds, EventsCount);
+
+
+    if (std::rename(oldFilename, newFilename) == 0) {
+      std::cout << "File renamed to: " << newFilename << std::endl;
+    } else {
+      std::cerr << "Error renaming file: " << oldFilename << std::endl;
+    }
+  }
+
+
+  // Done.
+  return 0;
+}