From 7622c5072c73677d08f5898cd6b8cf1da9fee256 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Tue, 20 Oct 2020 23:05:09 -0500
Subject: [PATCH] added reconstruction script for RICH
---
rich/forward_hadrons.sh | 70 ++++++++++++++++++++++++++++
rich/options/rich_reco.py | 33 +++++++++++++
rich/rich_config.yml | 14 ++++++
rich/scripts/rich_data_gen.py | 87 +++++++++++++++++++++++++++++++++++
rich/scripts/rich_samples.py | 44 ++++++++++++++++++
5 files changed, 248 insertions(+)
create mode 100644 rich/forward_hadrons.sh
create mode 100644 rich/options/rich_reco.py
create mode 100644 rich/rich_config.yml
create mode 100644 rich/scripts/rich_data_gen.py
create mode 100644 rich/scripts/rich_samples.py
diff --git a/rich/forward_hadrons.sh b/rich/forward_hadrons.sh
new file mode 100644
index 00000000..4ef74e96
--- /dev/null
+++ b/rich/forward_hadrons.sh
@@ -0,0 +1,70 @@
+#!/bin/bash
+
+if [[ ! -n "${JUGGLER_DETECTOR}" ]] ; then
+ export JUGGLER_DETECTOR="topside"
+fi
+
+if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
+ export JUGGLER_N_EVENTS=100
+fi
+
+if [[ ! -n "${P_start}" ]] ; then
+ export P_start=5.0
+fi
+
+if [[ ! -n "${P_end}" ]] ; then
+ export P_end=100.0
+fi
+
+if [[ ! -n "${Angle_start}" ]] ; then
+ export Angle_start=3.0
+fi
+
+if [[ ! -n "${Angle_end}" ]] ; then
+ export Angle_end=8.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}"
+
+
+# Build the detector constructors.
+git clone https://eicweb.phy.anl.gov/EIC/detectors/${JUGGLER_DETECTOR}.git
+mkdir ${JUGGLER_DETECTOR}/build
+pushd ${JUGGLER_DETECTOR}/build
+cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j30 install
+popd
+
+# generate the input events
+# note datasets is now only used to develop datasets.
+#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
+python rich/scripts/rich_data_gen.py \
+ ${JUGGLER_FILE_NAME_TAG}.hepmc \
+ -n ${JUGGLER_N_EVENTS} \
+ --pmin ${P_start} \
+ --pmax ${P_end} \
+ --angmin ${Angle_start} \
+ --angmax ${Angle_end}
+
+pushd ${JUGGLER_DETECTOR}
+ls -l
+# run geant4 simulations
+python options/ForwardRICH/simu.py \
+ --compact=${JUGGLER_DETECTOR}.xml \
+ -i ../${JUGGLER_FILE_NAME_TAG} \
+ -o ${JUGGLER_SIM_FILE}
+ -n ${JUGGLER_N_EVENTS}
+
+# @TODO changeable simulation file name and detector xml file name
+xenv -x /usr/local/Juggler.xenv gaudirun.py ../rich/options/rich_reco.py
+ls -l
+popd
+
+# @TODO add analysis scripts
+
diff --git a/rich/options/rich_reco.py b/rich/options/rich_reco.py
new file mode 100644
index 00000000..d0557aeb
--- /dev/null
+++ b/rich/options/rich_reco.py
@@ -0,0 +1,33 @@
+from Gaudi.Configuration import *
+from GaudiKernel import SystemOfUnits as units
+
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+
+geo_service = GeoSvc("GeoSvc", detectors=["topside.xml"])
+podioevent = EICDataSvc("EventDataSvc", inputs=["rich_test.root"], OutputLevel=DEBUG)
+
+from Configurables import PodioInput
+from Configurables import Jug__Digi__PhotoMultiplierDigi as PhotoMultiplierDigi
+from Configurables import Jug__Reco__PhotoMultiplierReco as PhotoMultiplierReco
+from Configurables import Jug__Reco__PhotoRingClusters as PhotoRingClusters
+
+qe_data = [(1.0, 0.25), (7.5, 0.25),]
+podioinput = PodioInput("PodioReader", collections=["mcparticles", "ForwardRICHHits"], OutputLevel=DEBUG)
+pmtdigi = PhotoMultiplierDigi(inputHitCollection="ForwardRICHHits", outputHitCollection="DigiForwardRICHHits",
+ quantumEfficiency=[(a*units.eV, b) for a, b in qe_data])
+pmtreco = PhotoMultiplierReco(inputHitCollection="DigiForwardRICHHits", outputHitCollection="RecoForwardRICHHits")
+richcluster = PhotoRingClusters(inputHitCollection="RecoForwardRICHHits", inputTrackCollection="mcparticles",
+ outputClusterCollection="RICHClusters")
+
+out = PodioOutput("out", filename="rich_test_reco.root")
+out.outputCommands = ["keep *"]
+
+ApplicationMgr(
+ TopAlg = [podioinput, pmtdigi, pmtreco, richcluster, out],
+ EvtSel = 'NONE',
+ EvtMax = 100000,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+ )
+
diff --git a/rich/rich_config.yml b/rich/rich_config.yml
new file mode 100644
index 00000000..4fc4c395
--- /dev/null
+++ b/rich/rich_config.yml
@@ -0,0 +1,14 @@
+rich_job_x:
+ image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:latest
+ tags:
+ - silicon
+ needs: ["get_data"]
+ timeout: 24 hours
+ artifacts:
+ expire_in: 20 weeks
+ paths:
+ - results/
+ stage: run
+ script:
+ - bash rich/forward_hadrons.sh
+
diff --git a/rich/scripts/rich_data_gen.py b/rich/scripts/rich_data_gen.py
new file mode 100644
index 00000000..c609a888
--- /dev/null
+++ b/rich/scripts/rich_data_gen.py
@@ -0,0 +1,87 @@
+import os
+from pyHepMC3 import HepMC3 as hm
+import numpy as np
+import argparse
+
+
+PARTICLES = [
+# (111, 0.134.9766), # pi0
+ (211, 0.13957018), # pi+
+# (-211, 0.13957018), # pi-
+# (311, 0.497648), # K0
+ (321, 0.493677), # K+
+# (-321, 0.493677), # K-
+ (2212, 0.938272), # proton
+# (2112, 0.939565), # neutron
+# (11, 0.51099895e-3), # electron
+# (-11, 0.51099895e-3), # positron
+]
+
+
+# p in GeV, angle in degree, vertex in mm
+def gen_event(prange=(8, 100), arange=(0, 20)):
+ evt = hm.GenEvent(momentum_unit=hm.Units.MomentumUnit.GEV, length_unit=hm.Units.LengthUnit.MM)
+ pid, mass = PARTICLES[np.random.randint(len(PARTICLES))]
+
+ # final state
+ state = 1
+
+ # momentum, angles, energy
+ p = np.random.uniform(*prange)
+ theta = np.random.uniform(*arange)*np.pi/180.
+ phi = np.random.uniform(0., 2.*np.pi)
+ e0 = np.sqrt(p*p + mass*mass)
+
+ px = np.cos(phi)*np.sin(theta)
+ py = np.sin(phi)*np.sin(theta)
+ pz = np.cos(theta)
+
+ # beam
+ pbeam = hm.GenParticle(hm.FourVector(0, 0, 0, 0.938272), 2212, 4)
+ ebeam = hm.GenParticle(hm.FourVector(0, 0, e0, np.sqrt(e0*e0 + 0.511e-3*0.511e-3)), -11, 4)
+
+ hout = hm.GenParticle(hm.FourVector(px*p, py*p, pz*p, e0), pid, state)
+ # evt.add_particle(part)
+ vert = hm.GenVertex()
+ vert.add_particle_in(ebeam)
+ vert.add_particle_in(pbeam)
+ vert.add_particle_out(hout)
+ evt.add_vertex(vert)
+ return evt
+
+
+if __name__ == "__main__":
+ parser = argparse.ArgumentParser('RICH dataset generator')
+
+ parser.add_argument('output', help='path to the output file')
+ parser.add_argument('-n', type=int, default=1000, dest='nev', help='number of events to generate')
+ parser.add_argument('--pmin', type=float, default=8.0, dest='pmin', help='minimum momentum in GeV')
+ parser.add_argument('--pmax', type=float, default=100.0, dest='pmax', help='maximum momentum in GeV')
+ parser.add_argument('--angmin', type=float, default=0.0, dest='angmin', help='minimum angle in degree')
+ parser.add_argument('--angmax', type=float, default=20.0, dest='angmax', help='maximum angle in degree')
+ parser.add_argument('-s', type=int, default=-1, dest='seed', help='seed for random generator')
+
+ args = parser.parse_args()
+
+ # random seed
+ if args.seed < 0:
+ args.seed = os.environ.get('SEED', int.from_bytes(os.urandom(4), byteorder='big', signed=False))
+ np.random.seed(args.seed)
+
+ output = hm.WriterAscii(args.output);
+ if output.failed():
+ print("Cannot open file \"{}\"".format(args.output))
+ sys.exit(2)
+
+ count = 0
+ while count < args.nev:
+ if (count % 1000 == 0):
+ print("Generated {} events".format(count), end='\r')
+ evt = gen_event((args.pmin, args.pmax), (args.angmin, args.angmax))
+ output.write_event(evt)
+ evt.clear()
+ count += 1
+
+ print("Generated {} events".format(args.nev))
+ output.close()
+
diff --git a/rich/scripts/rich_samples.py b/rich/scripts/rich_samples.py
new file mode 100644
index 00000000..0dcc2bdb
--- /dev/null
+++ b/rich/scripts/rich_samples.py
@@ -0,0 +1,44 @@
+import ROOT
+import numpy as np
+from ROOT import gROOT
+import argparse
+from matplotlib import pyplot as plt
+from matplotlib.patches import Ellipse
+import matplotlib.transforms as transforms
+
+parser = argparse.ArgumentParser(description='rich reconstruction analysis')
+parser.add_argument('nevents', type=int, help='number of events')
+args = parser.parse_args()
+
+gROOT.Macro('rootlogon.C')
+f = ROOT.TFile("rich_test_reco.root")
+events = f.events
+
+nev = 0
+while nev < args.nevents:
+ iev = np.random.randint(0, events.GetEntries())
+ events.GetEntry(iev)
+ if events.RICHClusters.size() == 0:
+ continue
+
+ nev += 1
+ x, y = [], []
+ for hit in events.RecoForwardRICHHits:
+ x.append(hit.position.x)
+ y.append(hit.position.y)
+
+ fig, ax = plt.subplots(figsize=(9, 9), dpi=120)
+ colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+ ax.scatter(x, y, s=9, marker='s', color=colors[0])
+
+ for cluster in events.RICHClusters:
+ ellipse = Ellipse((cluster.position.x, cluster.position.y), width=cluster.radius * 2, height=cluster.radius * 2,
+ edgecolor='red', fill=False, lw=2, alpha=0.8)
+ ellipse.set_transform(ax.transData)
+ ax.add_patch(ellipse)
+ ax.set_xlim(-60, 60)
+ ax.set_ylim(-60, 60)
+ ax.set_xlabel('X (mm)')
+ ax.set_ylabel('Y (mm)')
+ fig.savefig('rich_cluster_{}.png'.format(iev))
+
--
GitLab