From 6eec429b1a0a1d2ac7dd656c1c65f31eb60cc384 Mon Sep 17 00:00:00 2001
From: Jihee Kim <jihee.kim@anl.gov>
Date: Mon, 17 May 2021 15:03:15 +0000
Subject: [PATCH] Resolve "ZDC benchmarks"
---
benchmarks/zdc/config.yml | 24 +---
benchmarks/zdc/run_zdc_neutrons.sh | 49 +++++++
benchmarks/zdc/scripts/zdc_neutrons.cxx | 81 +++++++++++
.../zdc/scripts/zdc_neutrons_analysis.cxx | 135 ++++++++++++++++++
4 files changed, 268 insertions(+), 21 deletions(-)
create mode 100755 benchmarks/zdc/run_zdc_neutrons.sh
create mode 100644 benchmarks/zdc/scripts/zdc_neutrons.cxx
create mode 100644 benchmarks/zdc/scripts/zdc_neutrons_analysis.cxx
diff --git a/benchmarks/zdc/config.yml b/benchmarks/zdc/config.yml
index 46bd6e79..f5426495 100644
--- a/benchmarks/zdc/config.yml
+++ b/benchmarks/zdc/config.yml
@@ -2,16 +2,7 @@ sim:zdc:
extends: .det_benchmark
stage: simulate
script:
- - bash benchmarks/zdc/run_simulation_zdc.sh
-
-zdc_neutrons:
- extends: .det_benchmark
- stage: benchmarks
- needs:
- - ["sim:zdc"]
- script:
- - echo " Not yet complete"
- #- root -b -q benchmarks/zdc/zdc_neutrons_reader.cxx
+ - bash benchmarks/zdc/run_zdc_neutrons.sh
bench:zdc_benchmark:
extends: .det_benchmark
@@ -19,21 +10,12 @@ bench:zdc_benchmark:
needs:
- ["sim:zdc"]
script:
- - echo " Not yet complete"
- #- root -b -q benchmarks/zdc/simple_checking.cxx+
-
-bench:zdc_benchmark_info_histogram:
- extends: .det_benchmark
- stage: benchmarks
- needs:
- - ["sim:zdc"]
- script:
- - root -b -q benchmarks/zdc/simple_info_plot_histograms.cxx+
+ - root -b -q benchmarks/zdc/scripts/zdc_neutrons_analysis.cxx+
collect_results:zdc:
extends: .det_benchmark
stage: collect
needs:
- - ["bench:zdc_benchmark","bench:zdc_benchmark_info_histogram"]
+ - ["bench:zdc_benchmark"]
script:
- ls -lrht
diff --git a/benchmarks/zdc/run_zdc_neutrons.sh b/benchmarks/zdc/run_zdc_neutrons.sh
new file mode 100755
index 00000000..981056e8
--- /dev/null
+++ b/benchmarks/zdc/run_zdc_neutrons.sh
@@ -0,0 +1,49 @@
+#!/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 "${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="zdc_uniform_neutrons"
+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 "benchmarks/zdc/scripts/zdc_neutrons.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
+
+# Run geant4 simulations
+npsim --runType batch \
+ -v WARNING \
+ --part.minimalKineticEnergy 0.5*GeV \
+ --numberOfEvents ${JUGGLER_N_EVENTS} \
+ --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+ --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
+ --outputFile sim_output/${JUGGLER_SIM_FILE}
+
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running npdet"
+ exit 1
+fi
+
diff --git a/benchmarks/zdc/scripts/zdc_neutrons.cxx b/benchmarks/zdc/scripts/zdc_neutrons.cxx
new file mode 100644
index 00000000..911d7674
--- /dev/null
+++ b/benchmarks/zdc/scripts/zdc_neutrons.cxx
@@ -0,0 +1,81 @@
+//////////////////////////////////////////////////////////////
+// ZDC detector
+// Single Neutron dataset
+// J.KIM 05/07/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 zdc_neutrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/zdc_neutrons.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 * (0.0 / 180.0));
+ double cos_theta_max = std::cos(M_PI * (5.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
+ // pdgid 2112 - neutron
+ 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 11 - electron 0.510 MeV/c^2
+ GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.939 * 0.939))), 2112, 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;
+}
diff --git a/benchmarks/zdc/scripts/zdc_neutrons_analysis.cxx b/benchmarks/zdc/scripts/zdc_neutrons_analysis.cxx
new file mode 100644
index 00000000..be218bb2
--- /dev/null
+++ b/benchmarks/zdc/scripts/zdc_neutrons_analysis.cxx
@@ -0,0 +1,135 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1D.h"
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+void zdc_neutrons_analysis(const char* input_fname = "sim_output/sim_zdc_uniform_neutrons.root")
+{
+ // Setting for graphs
+ gROOT->SetStyle("Plain");
+ gStyle->SetOptFit(1);
+ gStyle->SetLineWidth(2);
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetPadLeftMargin(0.14);
+ gStyle->SetPadRightMargin(0.14);
+
+ ROOT::EnableImplicitMT();
+ ROOT::RDataFrame d0("events", input_fname);
+
+ // Thrown Energy [GeV]
+ auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ auto p = input[2];
+ auto energy = TMath::Sqrt(p.psx * p.psx + p.psy * p.psy + p.psz * p.psz + p.mass * p.mass);
+ return energy;
+ };
+
+ // Number of hits
+ auto nhits = [] (const std::vector<dd4pod::CalorimeterHitData>& evt) {return (int) evt.size(); };
+
+ // Energy deposition [GeV]
+ auto Esim = [](const std::vector<dd4pod::CalorimeterHitData>& evt) {
+ auto total_edep = 0.0;
+ for (const auto& i: evt)
+ total_edep += i.energyDeposit;
+ return total_edep;
+ };
+
+ // Sampling fraction = Esampling / Ethrown
+ auto fsam = [](const double sampled, const double thrown) {
+ return sampled / thrown;
+ };
+
+ // Define variables
+ auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
+ .Define("nhits", nhits, {"EcalBarrelHits"})
+ .Define("Esim", Esim, {"EcalBarrelHits"})
+ .Define("fsam", fsam, {"Esim", "Ethr"});
+
+ // Define Histograms
+ auto hEthr = d1.Histo1D(
+ {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5},
+ "Ethr");
+ auto hNhits =
+ d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events",
+ 100, 0.0, 2000.0},
+ "nhits");
+ auto hEsim = d1.Histo1D(
+ {"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 100, 0.0, 1.0},
+ "Esim");
+ auto hfsam = d1.Histo1D(
+ {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1},
+ "fsam");
+
+ // Event Counts
+ auto nevents_thrown = d1.Count();
+ std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
+
+ // Draw Histograms
+ {
+ TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
+ c1->SetLogy(1);
+ auto h = hEthr->DrawCopy();
+ //h->GetYaxis()->SetTitleOffset(1.4);
+ h->SetLineWidth(2);
+ h->SetLineColor(kBlue);
+ c1->SaveAs("results/zdc_neutrons_Ethr.png");
+ c1->SaveAs("results/zdc_neutrons_Ethr.pdf");
+ }
+ std::cout << "derp1\n";
+
+ {
+ TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
+ c2->SetLogy(1);
+ auto h = hNhits->DrawCopy();
+ //h->GetYaxis()->SetTitleOffset(1.4);
+ h->SetLineWidth(2);
+ h->SetLineColor(kBlue);
+ c2->SaveAs("results/zdc_neutrons_nhits.png");
+ c2->SaveAs("results/zdc_neutrons_nhits.pdf");
+ }
+
+ {
+ TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
+ c3->SetLogy(1);
+ auto h = hEsim->DrawCopy();
+ //h->GetYaxis()->SetTitleOffset(1.4);
+ h->SetLineWidth(2);
+ h->SetLineColor(kBlue);
+ c3->SaveAs("results/zdc_neutrons_Esim.png");
+ c3->SaveAs("results/zdc_neutrons_Esim.pdf");
+ }
+
+ std::cout << "derp4\n";
+ {
+ TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
+ c4->SetLogy(1);
+ auto h = hfsam->DrawCopy();
+ //h->GetYaxis()->SetTitleOffset(1.4);
+ h->SetLineWidth(2);
+ h->SetLineColor(kBlue);
+ //h->Fit("gaus", "", "", 0.01, 0.1);
+ //h->GetFunction("gaus")->SetLineWidth(2);
+ //h->GetFunction("gaus")->SetLineColor(kRed);
+ c4->SaveAs("results/zdc_neutrons_fsam.png");
+ c4->SaveAs("results/zdc_neutrons_fsam.pdf");
+ }
+}
--
GitLab