From 5f81476afe6146660ac4b7f23738fdbe8fe2832e Mon Sep 17 00:00:00 2001
From: Jihee Kim <jihee.kim@anl.gov>
Date: Sun, 4 Apr 2021 03:41:50 +0000
Subject: [PATCH] Resolve "CI Workflow for EM Barrel"
---
.gitlab-ci.yml | 120 ++++++------
calorimeters/calorimeters_config.yml | 50 +++--
calorimeters/options/emcal_barrel_reco.py | 107 +++++++++++
calorimeters/run_emcal_barrel_electrons.sh | 72 ++++++++
.../scripts/emcal_barrel_electrons.cxx | 80 ++++++++
.../emcal_barrel_electrons_analysis.cxx | 171 ++++++++++++++++++
.../scripts/emcal_barrel_electrons_reader.cxx | 125 +++++++++++++
7 files changed, 635 insertions(+), 90 deletions(-)
create mode 100644 calorimeters/options/emcal_barrel_reco.py
create mode 100755 calorimeters/run_emcal_barrel_electrons.sh
create mode 100644 calorimeters/scripts/emcal_barrel_electrons.cxx
create mode 100644 calorimeters/scripts/emcal_barrel_electrons_analysis.cxx
create mode 100644 calorimeters/scripts/emcal_barrel_electrons_reader.cxx
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index b1583677..301c4ae2 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,5 +1,4 @@
-image: eicweb.phy.anl.gov:4567/eic/npdet/npdet:latest
-
+image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
default:
artifacts:
@@ -9,21 +8,18 @@ default:
- datasets/
- sim_output/
- data
- # exclude:
- # - .git/
- # - datasets/.git/
+ reports:
+ dotenv: juggler.env
before_script:
- git clone https://eicweb.phy.anl.gov/EIC/NPDet.git
- git clone https://eicweb.phy.anl.gov/EIC/detectors/topside.git && mkdir topside/build && cd topside/build && cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j20 install && cd ../..
- # - cd NPDet/build && cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j10 && make install
+ - git clone https://eicweb.phy.anl.gov/EIC/detectors/accelerator.git && cd topside && ln -s ../accelerator/eic && cd ../.
stages:
- data_init
- # - ci_gen
- simulate
- benchmarks
- deploy
- #- others
get_data:
stage: data_init
@@ -44,15 +40,15 @@ get_data:
# - ./bin/gen_ci_config -p cal_test_ -i calorimeters > results/calorimeters_config.yml
# - ./bin/gen_ci_config -p pid_test_ -i pid > results/pid_config.yml
-cal_test_3_zdc_neutrons_reader:
- stage: benchmarks
- tags:
- - silicon
- needs:
- - ["zdc_simulation"]
- script:
- - root -b -q calorimeters/zdc_neutrons_reader.cxx
- allow_failure: true
+ #cal_test_3_zdc_neutrons_reader:
+ #stage: benchmarks
+ #tags:
+ #- silicon
+ #needs:
+ #- ["zdc_simulation"]
+ #script:
+ #- root -b -q calorimeters/zdc_neutrons_reader.cxx
+ #allow_failure: true
roman_pot_simu:
stage: simulate
@@ -61,8 +57,8 @@ roman_pot_simu:
tags:
- silicon
script:
- - cp NPDet/src/GenericDetectors/trackers/compact/elements.xml ./.
- - cp NPDet/src/GenericDetectors/trackers/compact/materials.xml ./.
+ - cp NPDet/src/detectors/trackers/compact/elements.xml ./.
+ - cp NPDet/src/detectors/trackers/compact/materials.xml ./.
- bash trackers/roman_pot_simu.sh
roman_pot_nhits:
@@ -85,43 +81,43 @@ roman_pot_eta:
- root -b -q trackers/roman_pot_hit_eta.cxx+
allow_failure: true
-zdc_simulation:
- stage: simulate
- needs:
- - ["get_data"]
- tags:
- - silicon
- script:
- - cp NPDet/src/GenericDetectors/calorimeters/compact/elements.xml ./.
- - cp NPDet/src/GenericDetectors/calorimeters/compact/materials.xml ./.
- - bash calorimeters/run_simulation_zdc.sh
-
-zdc_benchmark:
- stage: benchmarks
- tags:
- - silicon
- needs:
- - ["zdc_simulation"]
- dependencies:
- - zdc_simulation
- script:
- - ls -lrth sim_output
- - root -b -q calorimeters/simple_checking.cxx+
- allow_failure: true
-
-zdc_benchmark_info_histogram:
- stage: benchmarks
- needs:
- - ["zdc_simulation"]
- tags:
- - silicon
- dependencies:
- - zdc_simulation
- script:
- - cp NPDet/src/GenericDetectors/calorimeters/compact/elements.xml calorimeters/
- - cp NPDet/src/GenericDetectors/calorimeters/compact/materials.xml calorimeters/
- - root -b -q calorimeters/simple_info_plot_histograms.cxx+
- allow_failure: true
+ #zdc_simulation:
+ #stage: simulate
+ #needs:
+ #- ["get_data"]
+ #tags:
+ #- silicon
+ #script:
+ #- cp NPDet/src/detectors/calorimeters/compact/elements.xml ./.
+ #- cp NPDet/src/detectors/calorimeters/compact/materials.xml ./.
+ #- bash calorimeters/run_simulation_zdc.sh
+
+ #zdc_benchmark:
+ #stage: benchmarks
+ #tags:
+ #- silicon
+ #needs:
+ #- ["zdc_simulation"]
+ #dependencies:
+ #- zdc_simulation
+ #script:
+ #- ls -lrth sim_output
+ #- root -b -q calorimeters/simple_checking.cxx+
+ #allow_failure: true
+
+ #zdc_benchmark_info_histogram:
+ #stage: benchmarks
+ #needs:
+ #- ["zdc_simulation"]
+ #tags:
+ #- silicon
+ #dependencies:
+ #- zdc_simulation
+ #script:
+ #- cp NPDet/src/detectors/calorimeters/compact/elements.xml calorimeters/
+ #- cp NPDet/src/detectors/calorimeters/compact/materials.xml calorimeters/
+ #- root -b -q calorimeters/simple_info_plot_histograms.cxx+
+ #allow_failure: true
crystal_emcal_simulation:
stage: simulate
@@ -130,8 +126,8 @@ crystal_emcal_simulation:
tags:
- silicon
script:
- - cp NPDet/src/GenericDetectors/calorimeters/compact/elements.xml ./.
- - cp NPDet/src/GenericDetectors/calorimeters/compact/materials.xml ./.
+ - cp NPDet/src/detectors/calorimeters/compact/elements.xml ./.
+ - cp NPDet/src/detectors/calorimeters/compact/materials.xml ./.
- bash calorimeters/run_simulation_crystal.sh
crystal_benchmark:
@@ -155,16 +151,18 @@ crystal_pion_simulation:
- cd topside && ls -l
- npsim --runType batch --numberOfEvents 100 --compactFile topside.xml --inputFiles ../data/emcal_electrons.hepmc --outputFile ../sim_output/output_emcal_electrons.root
+include:
+ - local: 'calorimeters/calorimeters_config.yml'
+
deploy_results:
stage: deploy
- needs:
- - ["zdc_benchmark","zdc_benchmark_info_histogram"]
+ #needs:
+ #- ["zdc_benchmark","zdc_benchmark_info_histogram"]
tags:
- silicon
script:
- echo "deploy results!"
-
pages:
stage: deploy
rules:
diff --git a/calorimeters/calorimeters_config.yml b/calorimeters/calorimeters_config.yml
index a3f63a9d..013538ba 100644
--- a/calorimeters/calorimeters_config.yml
+++ b/calorimeters/calorimeters_config.yml
@@ -1,42 +1,34 @@
-cal_sim_1_dummy_test2:
- stage: benchmarks
- script:
- - echo "here we run simulation"
- #artifact:
- # paths:
- # - results/
- allow_failure: true
-
-cal_test_1_dummy_test2:
- stage: benchmarks
+sim_emcal_barrel_electrons:
+ image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
+ stage: simulate
tags:
- sodium
+ artifacts:
+ expire_in: 20 weeks
+ paths:
+ - results/
+ - sim_output/
script:
- - bash calorimeters/dummy_test2.sh
- #artifact:
- # paths:
- # - results/
+ - bash calorimeters/run_emcal_barrel_electrons.sh
allow_failure: true
-cal_test_2_dummy_test:
+ben_emcal_barrel_electrons:
stage: benchmarks
- tags:
- - sodium
+ needs:
+ - ["sim_emcal_barrel_electrons"]
script:
- - bash calorimeters/dummy_test.sh
- #artifact:
- # paths:
- # - results/
+ - ls -lrth sim_output
+ - root -b -q calorimeters/scripts/emcal_barrel_electrons_analysis.cxx+
allow_failure: true
-cal_test_3_zdc_neutrons_reader:
- stage: benchmarks
- tags:
- - sodium
- script:
- - root -b -q calorimeters/zdc_neutrons_reader.cxx
+ #cal_test_3_zdc_neutrons_reader:
+ #stage: benchmarks
+ #tags:
+ #- sodium
+ #script:
+ #- root -b -q calorimeters/zdc_neutrons_reader.cxx
#artifact:
# paths:
# - results/
- allow_failure: true
+ #allow_failure: true
diff --git a/calorimeters/options/emcal_barrel_reco.py b/calorimeters/options/emcal_barrel_reco.py
new file mode 100644
index 00000000..b93775f7
--- /dev/null
+++ b/calorimeters/options/emcal_barrel_reco.py
@@ -0,0 +1,107 @@
+#######################################
+# EMCAL Barrel detector Reconstruction
+# J.KIM 04/02/2021
+#######################################
+
+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"])+"/"+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=["sim_output/{}".format(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__EcalTungstenSamplingDigi as EcalTungstenSamplingDigi
+
+from Configurables import Jug__Reco__EcalTungstenSamplingReco as EcalTungstenSamplingReco
+from Configurables import Jug__Reco__SamplingECalHitsMerger as SamplingECalHitsMerger
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+podioinput = PodioInput("PodioReader", collections=["mcparticles","EcalBarrelHits"], OutputLevel=DEBUG)
+
+# Thrown Information
+copier = MCCopier("MCCopier",
+ inputCollection="mcparticles",
+ outputCollection="mcparticles2",
+ OutputLevel=DEBUG)
+# Geant4 Information
+embarrelcopier = CalCopier("CalBarrelCopier",
+ inputCollection="EcalBarrelHits",
+ outputCollection="EcalBarrelHits2",
+ OutputLevel=DEBUG)
+# Digitization
+embarreldigi = EcalTungstenSamplingDigi("ecal_barrel_digi",
+ inputHitCollection="EcalBarrelHits",
+ outputHitCollection="RawEcalBarrelHits",
+ inputEnergyUnit=units.GeV,
+ inputTimeUnit=units.ns,
+ OutputLevel=DEBUG)
+# Reconstruction
+embarrelreco = EcalTungstenSamplingReco("ecal_barrel_reco",
+ inputHitCollection="RawEcalBarrelHits",
+ outputHitCollection="RecoEcalBarrelHits",
+ OutputLevel=DEBUG)
+# 2D+1 Clusterings
+# readout id definition for barrel ecal
+# <id>system:8,barrel:3,module:4,layer:6,slice:5,x:32:-16,y:-16</id>
+# xy_merger sum layers/slices, masking (8+3+4, 8+3+4+5+6-1)
+embarrelxymerger = SamplingECalHitsMerger("ecal_barrel_xy_merger",
+ cellIDMaskRanges=[(15, 25)],
+ inputHitCollection="RecoEcalBarrelHits",
+ outputHitCollection="RecoEcalBarrelHitsXY")
+# xy_merger sum modules, masking (8+3+4+5+6, 8+3+4+5+6+32-1)
+embarrelzmerger = SamplingECalHitsMerger("ecal_barrel_z_merger",
+ cellIDMaskRanges=[(26, 57)],
+ inputHitCollection="RecoEcalBarrelHits",
+ outputHitCollection="RecoEcalBarrelHitsZ")
+# Clustering
+embarrelcluster = IslandCluster("ecal_barrel_cluster",
+ inputHitCollection="RecoEcalBarrelHitsXY",
+ outputClusterCollection="EcalBarrelClusters",
+ minClusterCenterEdep=5.0*units.MeV,
+ splitCluster=False,
+ groupRange=5.0)
+# Reconstruct the cluster with Center of Gravity method
+embarrelclusterreco = RecoCoG("ecal_barrel_clusterreco",
+ clusterCollection="EcalBarrelClusters",
+ logWeightBase=6.2)
+
+out = PodioOutput("out", filename=output_rec_file)
+
+out.outputCommands = ["keep *"]
+
+ApplicationMgr(
+ TopAlg = [podioinput, copier, embarrelcopier, embarreldigi,
+ embarrelreco, embarrelxymerger, embarrelzmerger, embarrelcluster, embarrelclusterreco, out],
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+ )
diff --git a/calorimeters/run_emcal_barrel_electrons.sh b/calorimeters/run_emcal_barrel_electrons.sh
new file mode 100755
index 00000000..3634a1cb
--- /dev/null
+++ b/calorimeters/run_emcal_barrel_electrons.sh
@@ -0,0 +1,72 @@
+#!/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_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
+root -b -q "calorimeters/scripts/emcal_barrel_electrons.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_electrons_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
+
+# Run geant4 simulations
+npsim --runType batch \
+ -v WARNING \
+ --part.minimalKineticEnergy 0.5*GeV \
+ --numberOfEvents ${JUGGLER_N_EVENTS} \
+ --compactFile topside/${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
+
+# Run Juggler
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py calorimeters/options/emcal_barrel_reco.py
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running juggler"
+ exit 1
+fi
+
+# Directory for plots
+mkdir -p results
+
+# Move ROOT output file
+mv ${JUGGLER_REC_FILE} sim_output/
+
diff --git a/calorimeters/scripts/emcal_barrel_electrons.cxx b/calorimeters/scripts/emcal_barrel_electrons.cxx
new file mode 100644
index 00000000..12a78da3
--- /dev/null
+++ b/calorimeters/scripts/emcal_barrel_electrons.cxx
@@ -0,0 +1,80 @@
+//////////////////////////////////////////////////////////////
+// EMCAL Barrel detector
+// Single Electron dataset
+// J.KIM 04/02/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_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_barrel_electrons.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 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/calorimeters/scripts/emcal_barrel_electrons_analysis.cxx b/calorimeters/scripts/emcal_barrel_electrons_analysis.cxx
new file mode 100644
index 00000000..1ad13946
--- /dev/null
+++ b/calorimeters/scripts/emcal_barrel_electrons_analysis.cxx
@@ -0,0 +1,171 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+#include "dd4pod/TrackerHitCollection.h"
+#include "eicd/RawCalorimeterHitCollection.h"
+#include "eicd/RawCalorimeterHitData.h"
+#include "eicd/CalorimeterHitCollection.h"
+#include "eicd/CalorimeterHitData.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.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 emcal_barrel_electrons_analysis(const char* input_fname = "sim_output/rec_emcal_barrel_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 Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ std::vector<double> 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;
+ };
+
+ // Reconstructed Energy [GeV] in XY merger
+ auto ErecXY = [] (const std::vector<eic::CalorimeterHitData> & evt) {
+ std::vector<double> result;
+ auto total_eng = 0.0;
+ for (const auto& i: evt)
+ total_eng += i.energy;
+ result.push_back(total_eng / 1.e+3);
+ return result;
+ };
+
+ // Reconstructed Energy [GeV] in Z merger
+ auto ErecZ = [] (const std::vector<eic::CalorimeterHitData> & evt) {
+ std::vector<double> result;
+ auto total_eng = 0.0;
+ for (const auto& i: evt)
+ total_eng += i.energy;
+ result.push_back(total_eng / 1.e+3);
+ return result;
+ };
+
+ // Number of Clusters
+ auto ncluster = [] (const std::vector<eic::ClusterData>& evt) {return (int) evt.size(); };
+
+ // Cluster Energy [GeV]
+ auto Ecluster = [] (const std::vector<eic::ClusterData>& evt) {
+ std::vector<double> result;
+ for (const auto& i: evt)
+ result.push_back(i.energy / 1.e+3);
+ return result;
+ };
+
+ // Sampling fraction = Esampling / Ethrown
+ auto fsam = [](const std::vector<double>& sampled, const std::vector<double>& thrown) {
+ std::vector<double> result;
+ for (const auto& E1 : thrown) {
+ for (const auto& E2 : sampled)
+ result.push_back(E2/E1);
+ }
+ return result;
+ };
+
+ // Define variables
+ auto d1 = d0.Define("Ethr", Ethr, {"mcparticles2"})
+ .Define("ErecXY", ErecXY, {"RecoEcalBarrelHitsXY"})
+ .Define("ErecZ", ErecZ, {"RecoEcalBarrelHitsZ"})
+ .Define("ncluster", ncluster, {"EcalBarrelClusters"})
+ .Define("Ecluster", Ecluster, {"EcalBarrelClusters"})
+ .Define("fsam", fsam, {"Ecluster","Ethr"})
+ ;
+
+ // Define Histograms
+ auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, -0.5, 10.5}, "Ethr");
+ auto hErecXY = d1.Histo1D({"hErecXY", "Reconstructed Energy in XY merger; Reconstructed Energy [GeV]; Events", 100, -0.5, 10.5}, "ErecXY");
+ auto hErecZ = d1.Histo1D({"hErecZ", "Reconstructed Energy in Z merger; Reconstructed Energy [GeV]; Events", 100, -0.5, 10.5}, "ErecZ");
+ auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events", 20, -0.5, 20.5}, "ncluster");
+ auto hEcluster = d1.Histo1D({"hEcluster", "Cluster Energy; Cluster Energy [GeV]; Events", 100, -0.5, 10.5}, "Ecluster");
+ auto hfsam = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 1.0}, "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);
+ 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");
+
+ TCanvas *c2 = new TCanvas("c2", "c2", 700, 500);
+ c2->SetLogy(1);
+ hErecXY->GetYaxis()->SetTitleOffset(1.4);
+ hErecXY->SetLineWidth(2);
+ hErecXY->SetLineColor(kBlue);
+ hErecXY->DrawClone();
+ c2->SaveAs("results/emcal_electrons_ErecXY.png");
+ c2->SaveAs("results/emcal_electrons_ErecXY.pdf");
+
+ TCanvas *c3 = new TCanvas("c3", "c3", 700, 500);
+ c3->SetLogy(1);
+ hErecZ->GetYaxis()->SetTitleOffset(1.4);
+ hErecZ->SetLineWidth(2);
+ hErecZ->SetLineColor(kBlue);
+ hErecZ->DrawClone();
+ c3->SaveAs("results/emal_electrons_ErecZ.png");
+ c3->SaveAs("results/emal_electrons_ErecZ.pdf");
+
+ TCanvas *c4 = new TCanvas("c4", "c4", 700, 500);
+ c4->SetLogy(1);
+ hNCluster->GetYaxis()->SetTitleOffset(1.6);
+ hNCluster->SetLineWidth(2);
+ hNCluster->SetLineColor(kBlue);
+ hNCluster->DrawClone();
+ c4->SaveAs("results/emcal_electrons_ncluster.png");
+ c4->SaveAs("results/emcal_electrons_ncluster.pdf");
+
+ TCanvas *c5 = new TCanvas("c5", "c5", 700, 500);
+ c5->SetLogy(1);
+ hEcluster->GetYaxis()->SetTitleOffset(1.4);
+ hEcluster->SetLineWidth(2);
+ hEcluster->SetLineColor(kBlue);
+ hEcluster->DrawClone();
+ c5->SaveAs("results/emcal_electrons_Ecluster.png");
+ c5->SaveAs("results/emcal_electrons_Ecluster.pdf");
+
+ TCanvas *c6 = new TCanvas("c6", "c6", 700, 500);
+ c6->SetLogy(1);
+ hfsam->GetYaxis()->SetTitleOffset(1.4);
+ hfsam->SetLineWidth(2);
+ hfsam->SetLineColor(kBlue);
+ hfsam->Fit("gaus","","",0.1,1.0);
+ hfsam->GetFunction("gaus")->SetLineWidth(2);
+ hfsam->GetFunction("gaus")->SetLineColor(kRed);
+ hfsam->DrawClone();
+ c6->SaveAs("results/emcal_electrons_fsam.png");
+ c6->SaveAs("results/emcal_electrons_fsam.pdf");
+}
diff --git a/calorimeters/scripts/emcal_barrel_electrons_reader.cxx b/calorimeters/scripts/emcal_barrel_electrons_reader.cxx
new file mode 100644
index 00000000..75c2dae9
--- /dev/null
+++ b/calorimeters/scripts/emcal_barrel_electrons_reader.cxx
@@ -0,0 +1,125 @@
+//////////////////////////
+// EMCAL Barrel detector
+// Electron dataset
+// J.KIM 04/02/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_electrons_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_barrel_electrons.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_electrons_energy = new TH1F("h_electron_energy", "electron energy;E [GeV];Events", 100, -0.5, 30.5);
+ TH1F* h_electrons_eta = new TH1F("h_electron_eta", "electron #eta;#eta;Events", 100, -10.0, 10.0);
+ TH1F* h_electrons_theta = new TH1F("h_electron_theta", "electron #theta;#theta [degree];Events", 100, -0.5, 180.5);
+ TH1F* h_electrons_phi = new TH1F("h_electron_phi", "electron #phi;#phi [degree];Events", 100, -180.5, 180.5);
+ TH2F* h_electrons_pzpt = new TH2F("h_electrons_pzpt", "electron pt vs pz;pt [GeV];pz [GeV]", 100, -0.5, 30.5, 100, -30.5, 30.5);
+ TH2F* h_electrons_pxpy = new TH2F("h_electrons_pxpy", "electron px vs py;px [GeV];py [GeV]", 100, -30.5, 30.5, 100, -30.5, 30.5);
+ TH3F* h_electrons_p = new TH3F("h_electron_p", "electron 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_electrons_energy->Fill(p->momentum().e());
+ h_electrons_eta->Fill(p->momentum().eta());
+ h_electrons_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
+ h_electrons_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
+ h_electrons_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
+ h_electrons_pxpy->Fill(p->momentum().px(), p->momentum().py());
+ h_electrons_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_electrons_energy->GetYaxis()->SetTitleOffset(1.8);
+ h_electrons_energy->SetLineWidth(2);
+ h_electrons_energy->SetLineColor(kBlue);
+ h_electrons_energy->DrawClone();
+ c->SaveAs("results/input_emcal_barrel_electrons_energy.png");
+ c->SaveAs("results/input_emcal_barrel_electrons_energy.pdf");
+
+ TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
+ h_electrons_eta->GetYaxis()->SetTitleOffset(1.9);
+ h_electrons_eta->SetLineWidth(2);
+ h_electrons_eta->SetLineColor(kBlue);
+ h_electrons_eta->DrawClone();
+ c1->SaveAs("results/input_emcal_barrel_electrons_eta.png");
+ c1->SaveAs("results/input_emcal_barrel_electrons_eta.pdf");
+
+ TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
+ h_electrons_theta->GetYaxis()->SetTitleOffset(1.8);
+ h_electrons_theta->SetLineWidth(2);
+ h_electrons_theta->SetLineColor(kBlue);
+ h_electrons_theta->DrawClone();
+ c2->SaveAs("results/input_emcal_barrel_electrons_theta.png");
+ c2->SaveAs("results/input_emcal_barrel_electrons_theta.pdf");
+
+ TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
+ h_electrons_phi->GetYaxis()->SetTitleOffset(1.8);
+ h_electrons_phi->SetLineWidth(2);
+ h_electrons_phi->GetYaxis()->SetRangeUser(0.0, h_electrons_phi->GetMaximum() + 100.0);
+ h_electrons_phi->SetLineColor(kBlue);
+ h_electrons_phi->DrawClone();
+ c3->SaveAs("results/input_emcal_barrel_electrons_phi.png");
+ c3->SaveAs("results/input_emcal_barrel_electrons_phi.pdf");
+
+ TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
+ h_electrons_pzpt->GetYaxis()->SetTitleOffset(1.4);
+ h_electrons_pzpt->SetLineWidth(2);
+ h_electrons_pzpt->SetLineColor(kBlue);
+ h_electrons_pzpt->DrawClone("COLZ");
+ c4->SaveAs("results/input_emcal_barrel_electrons_pzpt.png");
+ c4->SaveAs("results/input_emcal_barrel_electrons_pzpt.pdf");
+
+ TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
+ h_electrons_pxpy->GetYaxis()->SetTitleOffset(1.4);
+ h_electrons_pxpy->SetLineWidth(2);
+ h_electrons_pxpy->SetLineColor(kBlue);
+ h_electrons_pxpy->DrawClone("COLZ");
+ c5->SaveAs("results/input_emcal_barrel_electrons_pxpy.png");
+ c5->SaveAs("results/input_emcal_barrel_electrons_pxpy.pdf");
+
+ TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
+ h_electrons_p->GetYaxis()->SetTitleOffset(1.8);
+ h_electrons_p->GetXaxis()->SetTitleOffset(1.6);
+ h_electrons_p->GetZaxis()->SetTitleOffset(1.6);
+ h_electrons_p->SetLineWidth(2);
+ h_electrons_p->SetLineColor(kBlue);
+ h_electrons_p->DrawClone();
+ c6->SaveAs("results/input_emcal_barrel_electrons_p.png");
+ c6->SaveAs("results/input_emcal_barrel_electrons_p.pdf");
+}
--
GitLab