From 3ca96c6ee4d116d9b213aa1a1b85b4c09ce98d1d Mon Sep 17 00:00:00 2001
From: Jihee Kim <jihee.kim@anl.gov>
Date: Tue, 2 Mar 2021 16:30:44 +0000
Subject: [PATCH] EM Barrel and Endcap
---
benchmarks/ecal/config.yml | 12 +++
benchmarks/ecal/full_emcal_electrons.sh | 80 ++++++++++++++
.../ecal/options/full_em_calorimeter_reco.py | 101 ++++++++++++++++++
.../ecal/scripts/full_emcal_electrons.cxx | 89 +++++++++++++++
.../scripts/full_emcal_electrons_analysis.cxx | 70 ++++++++++++
5 files changed, 352 insertions(+)
create mode 100644 benchmarks/ecal/full_emcal_electrons.sh
create mode 100644 benchmarks/ecal/options/full_em_calorimeter_reco.py
create mode 100644 benchmarks/ecal/scripts/full_emcal_electrons.cxx
create mode 100644 benchmarks/ecal/scripts/full_emcal_electrons_analysis.cxx
diff --git a/benchmarks/ecal/config.yml b/benchmarks/ecal/config.yml
index a1902af7..ffdfdd73 100644
--- a/benchmarks/ecal/config.yml
+++ b/benchmarks/ecal/config.yml
@@ -21,3 +21,15 @@ ecal_1_emcal_pi0s:
stage: run
script:
- bash benchmarks/ecal/emcal_pi0s.sh
+
+emcal_barrel_electrons:
+ image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
+ needs: ["detector"]
+ timeout: 48 hours
+ artifacts:
+ expire_in: 20 weeks
+ paths:
+ - results/
+ stage: run
+ script:
+ - bash benchmarks/ecal/full_emcal_electrons.sh
diff --git a/benchmarks/ecal/full_emcal_electrons.sh b/benchmarks/ecal/full_emcal_electrons.sh
new file mode 100644
index 00000000..68ad8688
--- /dev/null
+++ b/benchmarks/ecal/full_emcal_electrons.sh
@@ -0,0 +1,80 @@
+#!/bin/bash
+
+./util/print_env.sh
+
+## To run the reconstruction, we need the following global variables:
+## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
+## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
+## - JUGGLER_DETECTOR_VERSION: the detector package we want to use for this benchmark
+## - DETECTOR_PATH: full path to the detector definitions
+##
+## You can ready options/env.sh for more in-depth explanations of the variables
+## and how they can be controlled.
+source options/env.sh
+export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
+
+if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
+ export JUGGLER_N_EVENTS=100
+fi
+
+if [[ ! -n "${E_start}" ]] ; then
+ export E_start=0.5
+fi
+
+if [[ ! -n "${E_end}" ]] ; then
+ export E_end=5.0
+fi
+
+export JUGGLER_FILE_NAME_TAG="emcal_uniform_electrons"
+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
+# note datasets is now only used to develop datasets.
+#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
+root -b -q "benchmarks/ecal/scripts/full_emcal_electrons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script"
+ exit 1
+fi
+
+# run geant4 simulations
+npsim --runType batch \
+ -v WARNING \
+ --part.minimalKineticEnergy 1*GeV \
+ --numberOfEvents ${JUGGLER_N_EVENTS} \
+ --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+ --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
+ --outputFile ${JUGGLER_SIM_FILE}
+# Need to figure out how to pass file name to juggler from the commandline
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running npdet"
+ exit 1
+fi
+
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py benchmarks/ecal/options/full_em_calorimeter_reco.py
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running juggler"
+ exit 1
+fi
+
+root -b -q "benchmarks/ecal/scripts/full_emcal_electrons_analysis.cxx(\"${JUGGLER_REC_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running analysis script"
+ exit 1
+fi
+
+root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
+if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
+ # file must be less than 10 MB to upload
+ if [[ "${root_filesize}" -lt "10000000" ]] ; then
+ cp ${JUGGLER_REC_FILE} results/.
+ fi
+fi
+
diff --git a/benchmarks/ecal/options/full_em_calorimeter_reco.py b/benchmarks/ecal/options/full_em_calorimeter_reco.py
new file mode 100644
index 00000000..a61364b0
--- /dev/null
+++ b/benchmarks/ecal/options/full_em_calorimeter_reco.py
@@ -0,0 +1,101 @@
+from Gaudi.Configuration import *
+import os
+
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
+
+detector_name = "topside"
+if "JUGGLER_DETECTOR" in os.environ :
+ detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+if "JUGGLER_DETECTOR_PATH" in os.environ :
+ detector_name = str(os.environ["JUGGLER_DETECTOR_PATH"])+"/"+detector_name
+
+input_sim_file = "jug_input.root"
+if "JUGGLER_SIM_FILE" in os.environ :
+ input_sim_file = str(os.environ["JUGGLER_SIM_FILE"])
+else :
+ print(" ERROR : JUGGLER_SIM_FILE not set" )
+
+output_rec_file = "jug_rec.root"
+if "JUGGLER_REC_FILE" in os.environ :
+ output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
+else :
+ print(" ERROR : JUGGLER_REC_FILE not set" )
+
+n_events = 100
+if "JUGGLER_N_EVENTS" in os.environ :
+ n_events = str(os.environ["JUGGLER_N_EVENTS"])
+
+geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)])
+podioevent = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=DEBUG)
+
+from Configurables import PodioInput
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
+
+from Configurables import Jug__Digi__CrystalEndcapsDigi as CrystalEndcapsDigi
+from Configurables import Jug__Digi__EcalTungstenSamplingDigi as EcalTungstenSamplingDigi
+
+from Configurables import Jug__Reco__CrystalEndcapsReco as CrystalEndcapsReco
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits","EcalBarrelHits","EcalEndcapHits"], OutputLevel=DEBUG)
+
+## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
+copier = MCCopier("MCCopier",
+ inputCollection="mcparticles",
+ outputCollection="mcparticles2",
+ OutputLevel=DEBUG)
+
+calcopier = CalCopier("CalCopier",
+ inputCollection="CrystalEcalHits",
+ outputCollection="CrystalEcalHits2",
+ OutputLevel=DEBUG)
+
+embarreldigi = EcalTungstenSamplingDigi("ec_barrel_digi",
+ inputHitCollection="EcalBarrelHits",
+ outputHitCollection="RawEcalBarrelHits",
+ energyResolution=0.11,
+ OutputLevel=DEBUG)
+
+emendcapdigi = EcalTungstenSamplingDigi("ec_endcap_digi",
+ inputHitCollection="EcalEndcapHits",
+ outputHitCollection="RawEcalEndcapHits",
+ energyResolution=0.07,
+ OutputLevel=DEBUG)
+
+emcaldigi = CrystalEndcapsDigi("ecal_digi",
+ inputHitCollection="CrystalEcalHits",
+ outputHitCollection="RawDigiEcalHits")
+
+emcalreco = CrystalEndcapsReco("ecal_reco",
+ inputHitCollection="RawDigiEcalHits",
+ outputHitCollection="RecoEcalHits",
+ minModuleEdep=0.00001*units.MeV)
+
+emcalcluster = IslandCluster("emcal_cluster",
+ inputHitCollection="RecoEcalHits",
+ outputClusterCollection="EcalClusters",
+ minClusterCenterEdep=50.0*units.MeV,
+ groupRange=2.0)
+
+clusterreco = RecoCoG("cluster_reco",
+ clusterCollection="EcalClusters",
+ logWeightBase=4.2,
+ moduleDimZName="CrystalBox_z_length")
+
+out = PodioOutput("out", filename=output_rec_file)
+
+out.outputCommands = ["keep *"]
+
+ApplicationMgr(
+ TopAlg = [podioinput, copier, calcopier,
+ embarreldigi, emendcapdigi, emcaldigi, emcalreco, emcalcluster, clusterreco, out],
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+ )
diff --git a/benchmarks/ecal/scripts/full_emcal_electrons.cxx b/benchmarks/ecal/scripts/full_emcal_electrons.cxx
new file mode 100644
index 00000000..f849c076
--- /dev/null
+++ b/benchmarks/ecal/scripts/full_emcal_electrons.cxx
@@ -0,0 +1,89 @@
+//////////////////////////////////////////////////////////////
+// Tungsten Sampling EMCAL detector
+// Barrel and Endcap
+// Single Electron dataset
+// J.KIM 02/12/2021
+//
+//////////////////////////////////////////////////////////////
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include <iostream>
+#include<random>
+#include<cmath>
+#include <math.h>
+#include <TMath.h>
+
+using namespace HepMC3;
+
+void full_emcal_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_electron_0GeVto30GeV_1000kEvt.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
+ // Momentum starting at 0 GeV/c to 30 GeV/c
+ 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.000511 * 0.000511))),
+ 11, 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/ecal/scripts/full_emcal_electrons_analysis.cxx b/benchmarks/ecal/scripts/full_emcal_electrons_analysis.cxx
new file mode 100644
index 00000000..25fbecb2
--- /dev/null
+++ b/benchmarks/ecal/scripts/full_emcal_electrons_analysis.cxx
@@ -0,0 +1,70 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TH1D.h"
+
+// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
+// export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
+// export ROOT_INCLUDE_PATH=/home/jihee/stow/development/include:$ROOT_INCLUDE_PATH
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+void full_emcal_electrons_analysis(const char* input_fname = "rec_emcal_uniform_electrons.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 E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ std::vector<float> result;
+ result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass));
+ return result;
+ };
+
+ // Define variables
+ auto d1 = d0.Define("E_thr", E_thr, {"mcparticles2"})
+ ;
+
+ // Define Histograms
+ auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; E_{thr}[GeV]; Events", 100, -0.5, 30.5}, "E_thr");
+
+ // Event Counts
+ auto nevents_thrown = d1.Count();
+
+ // Print out number of events
+ std::cout << "Number of Events: " << (*nevents_thrown) << " = all thrown events \n";
+
+ // Draw Histogram
+ TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
+ hEthr->GetYaxis()->SetTitleOffset(1.4);
+ hEthr->SetLineWidth(2);
+ hEthr->SetLineColor(kBlue);
+ hEthr->DrawClone();
+ c1->SaveAs("results/emcal_electrons_Ethr.png");
+ c1->SaveAs("results/emcal_electrons_Ethr.pdf");
+}
--
GitLab