added reconstruction script for RICH

rich/forward_hadrons.sh

+#!/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
+

rich/options/rich_reco.py

+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
+ )
+

rich/rich_config.yml

+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
+

rich/scripts/rich_data_gen.py

+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()
+

rich/scripts/rich_samples.py

+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))
+