Added scripts to generate pion dataset

calorimeters/calorimeters_config.yml

@@ -5,6 +5,20 @@
 # - Run Juggler
 #####################
 
+sim_emcal_barrel_pions:
+  image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
+  stage: simulate
+  tags:
+    - silicon
+  artifacts:
+    expire_in: 20 weeks
+    paths:
+      - results/
+      - sim_output/
+  script:
+    - bash calorimeters/run_emcal_barrel_pions.sh
+  allow_failure: true
+
 sim_emcal_barrel_electrons:
   image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
   stage: simulate

calorimeters/run_emcal_barrel_pions.sh

+#!/bin/bash
+
+if [[ ! -n  "${JUGGLER_DETECTOR}" ]] ; then 
+  export JUGGLER_DETECTOR="topside"
+fi
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
+  export JUGGLER_N_EVENTS=1000
+fi
+
+if [[ ! -n  "${JUGGLER_INSTALL_PREFIX}" ]] ; then 
+  export JUGGLER_INSTALL_PREFIX="/usr/local"
+fi
+
+if [[ ! -n  "${E_start}" ]] ; then
+  export E_start=5.0
+fi
+
+if [[ ! -n  "${E_end}" ]] ; then
+  export E_end=5.0
+fi
+
+export JUGGLER_FILE_NAME_TAG="emcal_barrel_uniform_pions"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+# Generate the input events
+root -b -q "calorimeters/scripts/emcal_barrel_pions.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script: generating input events"
+  exit 1
+fi
+# Plot the input events
+root -b -q "calorimeters/scripts/emcal_barrel_pions_reader.cxx(${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script: plotting input events"
+  exit 1
+fi
+
+# Directory for plots
+mkdir -p results
+

calorimeters/scripts/emcal_barrel_pions.cxx

+//////////////////////////////////////////////////////////////
+// EMCAL Barrel detector
+// Single Pion dataset
+// J.KIM 04/04/2021
+//////////////////////////////////////////////////////////////
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/Print.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+
+#include <TMath.h>
+#include <cmath>
+#include <iostream>
+#include <math.h>
+#include <random>
+
+using namespace HepMC3;
+
+void emcal_barrel_pions(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_barrel_pions.hepmc") {
+  WriterAscii hepmc_output(out_fname);
+  int events_parsed = 0;
+  GenEvent evt(Units::GEV, Units::MM);
+
+  // Random number generator
+  TRandom* r1 = new TRandom();
+
+  // Constraining the solid angle, but larger than that subtended by the
+  // detector
+  // https://indico.bnl.gov/event/7449/contributions/35966/attachments/27177/41430/EIC-DWG-Calo-03192020.pdf
+  // See a figure on slide 26
+  double cos_theta_min = std::cos(M_PI * (45.0 / 180.0));
+  double cos_theta_max = std::cos(M_PI * (135.0 / 180.0));
+
+  for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
+    // FourVector(px,py,pz,e,pdgid,status)
+    // type 4 is beam
+    // pdgid 11 - electron
+    // pdgid 111 - pi0
+    // pdgid 2212 - proton
+    GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+    GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
+
+    // Define momentum
+    Double_t p        = r1->Uniform(e_start, e_end);
+    Double_t phi      = r1->Uniform(0.0, 2.0 * M_PI);
+    Double_t costheta = r1->Uniform(cos_theta_min, cos_theta_max);
+    Double_t theta    = std::acos(costheta);
+    Double_t px       = p * std::cos(phi) * std::sin(theta);
+    Double_t py       = p * std::sin(phi) * std::sin(theta);
+    Double_t pz       = p * std::cos(theta);
+
+    // Generates random vectors, uniformly distributed over the surface of a
+    // sphere of given radius, in this case momentum.
+    // r1->Sphere(px, py, pz, p);
+
+    // type 1 is final state
+    // pdgid 211 - pion+ 139.570 MeV/c^2
+    GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.139570 * 0.139570))), 211, 1);
+
+    GenVertexPtr v1 = std::make_shared<GenVertex>();
+    v1->add_particle_in(p1);
+    v1->add_particle_in(p2);
+
+    v1->add_particle_out(p3);
+    evt.add_vertex(v1);
+
+    if (events_parsed == 0) {
+      std::cout << "First event: " << std::endl;
+      Print::listing(evt);
+    }
+
+    hepmc_output.write_event(evt);
+    if (events_parsed % 10000 == 0) {
+      std::cout << "Event: " << events_parsed << std::endl;
+    }
+    evt.clear();
+  }
+  hepmc_output.close();
+  std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}

calorimeters/scripts/emcal_barrel_pions_reader.cxx

+//////////////////////////
+// EMCAL Barrel detector
+// Pion dataset
+// J.KIM 04/04/2021
+//////////////////////////
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/Print.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+
+#include "TH1F.h"
+#include "TStyle.h"
+#include <iostream>
+
+using namespace HepMC3;
+
+void emcal_barrel_pions_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_barrel_pions.hepmc") {
+  // Setting for graphs
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptFit(1);
+  gStyle->SetLineWidth(1);
+  gStyle->SetPadTickX(1);
+  gStyle->SetPadTickY(1);
+  gStyle->SetPadGridX(1);
+  gStyle->SetPadGridY(1);
+  gStyle->SetPadLeftMargin(0.14);
+  gStyle->SetPadRightMargin(0.17);
+
+  ReaderAscii hepmc_input(in_fname);
+  int events_parsed = 0;
+  GenEvent evt(Units::GEV, Units::MM);
+
+  // Histograms
+  TH1F* h_pions_energy = new TH1F("h_pions_energy", "pion energy;E [GeV];Events",         100, -0.5, 30.5);
+  TH1F* h_pions_eta    = new TH1F("h_pions_eta",    "pion #eta;#eta;Events",              100, -10.0, 10.0);
+  TH1F* h_pions_theta  = new TH1F("h_pions_theta",  "pion #theta;#theta [degree];Events", 100, -0.5, 180.5);
+  TH1F* h_pions_phi    = new TH1F("h_pions_phi",    "pion #phi;#phi [degree];Events",     100, -180.5, 180.5);
+  TH2F* h_pions_pzpt   = new TH2F("h_pions_pzpt",   "pion pt vs pz;pt [GeV];pz [GeV]",    100, -0.5, 30.5, 100, -30.5, 30.5);
+  TH2F* h_pions_pxpy   = new TH2F("h_pions_pxpy",   "pion px vs py;px [GeV];py [GeV]",    100, -30.5, 30.5, 100, -30.5, 30.5);
+  TH3F* h_pions_p      = new TH3F("h_pions_p",      "pion p;px [GeV];py [GeV];pz [GeV]",  100, -30.5, 30.5, 100, -30.5, 30.5, 100, -30.5, 30.5);
+
+  while (!hepmc_input.failed()) {
+    // Read event from input file
+    hepmc_input.read_event(evt);
+    // If reading failed - exit loop
+    if (hepmc_input.failed())
+      break;
+
+    for (const auto& v : evt.vertices()) {
+      for (const auto& p : v->particles_out()) {
+        if (p->pid() == 11) {
+          h_pions_energy->Fill(p->momentum().e());
+          h_pions_eta->Fill(p->momentum().eta());
+          h_pions_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
+          h_pions_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
+          h_pions_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
+          h_pions_pxpy->Fill(p->momentum().px(), p->momentum().py());
+          h_pions_p->Fill(p->momentum().px(), p->momentum().py(), p->momentum().pz());
+        }
+      }
+    }
+    evt.clear();
+    events_parsed++;
+  }
+  std::cout << "Events parsed and written: " << events_parsed << std::endl;
+
+  TCanvas* c = new TCanvas("c", "c", 500, 500);
+  h_pions_energy->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_energy->SetLineWidth(2);
+  h_pions_energy->SetLineColor(kBlue);
+  h_pions_energy->DrawClone();
+  c->SaveAs("results/input_emcal_barrel_pions_energy.png");
+  c->SaveAs("results/input_emcal_barrel_pions_energy.pdf");
+
+  TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
+  h_pions_eta->GetYaxis()->SetTitleOffset(1.9);
+  h_pions_eta->SetLineWidth(2);
+  h_pions_eta->SetLineColor(kBlue);
+  h_pions_eta->DrawClone();
+  c1->SaveAs("results/input_emcal_barrel_pions_eta.png");
+  c1->SaveAs("results/input_emcal_barrel_pions_eta.pdf");
+
+  TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
+  h_pions_theta->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_theta->SetLineWidth(2);
+  h_pions_theta->SetLineColor(kBlue);
+  h_pions_theta->DrawClone();
+  c2->SaveAs("results/input_emcal_barrel_pions_theta.png");
+  c2->SaveAs("results/input_emcal_barrel_pions_theta.pdf");
+
+  TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
+  h_pions_phi->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_phi->SetLineWidth(2);
+  h_pions_phi->GetYaxis()->SetRangeUser(0.0, h_pions_phi->GetMaximum() + 100.0);
+  h_pions_phi->SetLineColor(kBlue);
+  h_pions_phi->DrawClone();
+  c3->SaveAs("results/input_emcal_barrel_pions_phi.png");
+  c3->SaveAs("results/input_emcal_barrel_pions_phi.pdf");
+
+  TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
+  h_pions_pzpt->GetYaxis()->SetTitleOffset(1.4);
+  h_pions_pzpt->SetLineWidth(2);
+  h_pions_pzpt->SetLineColor(kBlue);
+  h_pions_pzpt->DrawClone("COLZ");
+  c4->SaveAs("results/input_emcal_barrel_pions_pzpt.png");
+  c4->SaveAs("results/input_emcal_barrel_pions_pzpt.pdf");
+
+  TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
+  h_pions_pxpy->GetYaxis()->SetTitleOffset(1.4);
+  h_pions_pxpy->SetLineWidth(2);
+  h_pions_pxpy->SetLineColor(kBlue);
+  h_pions_pxpy->DrawClone("COLZ");
+  c5->SaveAs("results/input_emcal_barrel_pions_pxpy.png");
+  c5->SaveAs("results/input_emcal_barrel_pions_pxpy.pdf");
+
+  TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
+  h_pions_p->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_p->GetXaxis()->SetTitleOffset(1.6);
+  h_pions_p->GetZaxis()->SetTitleOffset(1.6);
+  h_pions_p->SetLineWidth(2);
+  h_pions_p->SetLineColor(kBlue);
+  h_pions_p->DrawClone();
+  c6->SaveAs("results/input_emcal_barrel_pions_p.png");
+  c6->SaveAs("results/input_emcal_barrel_pions_p.pdf");
+}