diff --git a/benchmarks/clustering/scripts/cluster_plots.py b/benchmarks/clustering/scripts/cluster_plots.py
index 8d4ec745b0878ccf301e28319f5cf7df882691b7..e2a53cfdeb1713cf13697f231191f80edecfc7e9 100644
--- a/benchmarks/clustering/scripts/cluster_plots.py
+++ b/benchmarks/clustering/scripts/cluster_plots.py
@@ -27,15 +27,23 @@ def flatten_collection(rdf, collection, cols=None):
cols = [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}.'.format(collection))]
else:
cols = ['{}.{}'.format(collection, c) for c in cols]
+ if not cols:
+ print('cannot find any branch under collection {}'.format(collection))
+ return pd.DataFrame()
+
data = rdf.AsNumpy(cols)
# flatten the data, add an event id to identify clusters from different events
evns = []
for i, vec in enumerate(data[cols[0]]):
evns += [i]*vec.size()
for n, vals in data.items():
- data[n] = np.asarray([v for vec in vals for v in vec])
+ # make sure ints are not converted to floats
+ typename = vals[0].__class__.__name__.lower()
+ dtype = np.int64 if 'int' in typename or 'long' in typename else np.float64
+ # type safe creation
+ data[n] = np.asarray([v for vec in vals for v in vec], dtype=dtype)
# build data frame
- dfp = pd.DataFrame(columns=cols, data=np.vstack(list(data.values())).T)
+ dfp = pd.DataFrame({c: pd.Series(v) for c, v in data.items()})
dfp.loc[:, 'event'] = evns
return dfp
@@ -152,6 +160,7 @@ if __name__ == '__main__':
# calculate eta
if 'eta' not in df.columns:
df.loc[:, 'eta'] = -np.log(np.tan(df['polar.theta'].values/2.))
+ # print(df[['eta', 'polar.theta', 'position.x', 'position.y', 'position.z']])
fig, axs = plt.subplots(2, 2, figsize=(12, 8), dpi=160)
ncl = df.groupby('event')['clusterID'].nunique().values
axs[0][0].hist(ncl, weights=np.repeat(1./float(ncl.shape[0]), ncl.shape[0]),
diff --git a/benchmarks/ecal/scripts/draw_cluters.py b/benchmarks/ecal/scripts/draw_cluters.py
index a346731d9f1596a56f74b89790166c94b8390907..fd03f2dcca16e271e856c8d95c48653297a50ee0 100644
--- a/benchmarks/ecal/scripts/draw_cluters.py
+++ b/benchmarks/ecal/scripts/draw_cluters.py
@@ -28,15 +28,23 @@ def flatten_collection(rdf, collection, cols=None):
cols = [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}.'.format(collection))]
else:
cols = ['{}.{}'.format(collection, c) for c in cols]
+ if not cols:
+ print('cannot find any branch under collection {}'.format(collection))
+ return pd.DataFrame()
+
data = rdf.AsNumpy(cols)
# flatten the data, add an event id to identify clusters from different events
evns = []
for i, vec in enumerate(data[cols[0]]):
evns += [i]*vec.size()
for n, vals in data.items():
- data[n] = np.asarray([v for vec in vals for v in vec])
+ # make sure ints are not converted to floats
+ typename = vals[0].__class__.__name__.lower()
+ dtype = np.int64 if 'int' in typename or 'long' in typename else np.float64
+ # type safe creation
+ data[n] = np.asarray([v for vec in vals for v in vec], dtype=dtype)
# build data frame
- dfp = pd.DataFrame(columns=cols, data=np.vstack(list(data.values())).T)
+ dfp = pd.DataFrame({c: pd.Series(v) for c, v in data.items()})
dfp.loc[:, 'event'] = evns
return dfp
diff --git a/benchmarks/imaging_ecal/config.yml b/benchmarks/imaging_ecal/config.yml
index 13b684f893fc7d39758046dc22aac2347dd180e9..0eb3d818efca2d832cb8f42eea745a82ea9e384a 100644
--- a/benchmarks/imaging_ecal/config.yml
+++ b/benchmarks/imaging_ecal/config.yml
@@ -19,3 +19,10 @@ imaging_ecal_pions:
script:
- bash benchmarks/imaging_ecal/run_emcal_barrel.sh -t emcal_barrel_pions -p "pion-" -n 100
+imaging_ecal_pion0:
+ extends: .rec_benchmark
+ timeout: 48 hours
+ stage: run
+ script:
+ - bash benchmarks/imaging_ecal/run_imcal_pion0.sh -t imcal_barrel_pion0 -p "pion0" -n 100
+
diff --git a/benchmarks/imaging_ecal/options/hybrid_cluster.py b/benchmarks/imaging_ecal/options/hybrid_cluster.py
index 75705fc052d82adc95fc256aa84b56307dbc803e..707247aba045066d6050e404fb2f9a91a4e94e81 100644
--- a/benchmarks/imaging_ecal/options/hybrid_cluster.py
+++ b/benchmarks/imaging_ecal/options/hybrid_cluster.py
@@ -105,6 +105,8 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
outputHitCollection="EcalBarrelScFiHitsReco",
thresholdFactor=5.0,
readoutClass="EcalBarrelScFiHits",
+ layerField="layer",
+ sectorField="module",
localDetFields=["system", "module"], # use local coordinates in each module (stave)
**scfi_barrel_daq)
@@ -123,7 +125,8 @@ scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
outputHitCollection="EcalBarrelScFiClusterHits",
splitCluster=False,
minClusterCenterEdep=10.*MeV,
- localDistXZ=[30*mm, 30*mm])
+ localDistXZ=[30*mm, 30*mm],
+ sectorDist=5.*cm)
scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
inputHitCollection="EcalBarrelScFiClusterHits",
diff --git a/benchmarks/imaging_ecal/options/imaging_2dcluster.py b/benchmarks/imaging_ecal/options/imaging_2dcluster.py
new file mode 100644
index 0000000000000000000000000000000000000000..73d3def1a70531d0040d9da79b4a2039cb52ae8b
--- /dev/null
+++ b/benchmarks/imaging_ecal/options/imaging_2dcluster.py
@@ -0,0 +1,95 @@
+import os
+import ROOT
+from Gaudi.Configuration import *
+from GaudiKernel.SystemOfUnits import mm, MeV, rad, ns
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+
+from Configurables import PodioInput
+from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
+from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
+from Configurables import Jug__Reco__CalorimeterHitsEtaPhiProjector as CalHitsProj
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+
+# input arguments through environment variables
+kwargs = dict()
+kwargs['sf'] = float(os.environ.get('CB_EMCAL_SCFI, SAMP_FRAC', '0.0134'))
+kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.root')
+kwargs['output'] = os.environ.get('CB_EMCAL_REC_FILE', 'barrel_pion0_5GeV_cluster.root')
+kwargs['compact'] = os.environ.get('CB_EMCAL_COMPACT_PATH', '../topside/test.xml')
+kwargs['nev'] = int(os.environ.get('CB_EMCAL_NUMEV', 100))
+
+if kwargs['nev'] < 1:
+ f = ROOT.TFile(kwargs['input'])
+ kwargs['nev'] = f.events.GetEntries()
+
+print(kwargs)
+
+# get sampling fraction from system environment variable, 1.0 by default
+sf = float(os.environ.get('CB_EMCAL_SAMP_FRAC', '1.0'))
+geo_service = GeoSvc('GeoSvc', detectors=kwargs['compact'].split(','), OutputLevel=INFO)
+podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
+
+podin = PodioInput('PodioReader', collections=['mcparticles', 'EcalBarrelHits'], OutputLevel=DEBUG)
+podout = PodioOutput('out', filename=kwargs['output'])
+
+# use the same daq_setting for digi/reco pair
+imcal_barrel_daq = dict(
+ dynamicRangeADC=3.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=50)
+
+imcal_barrel_digi = CalHitDigi('imcal_barrel_digi',
+ inputHitCollection='EcalBarrelHits',
+ outputHitCollection='EcalBarrelHitsDigi',
+ **imcal_barrel_daq)
+
+imcal_barrel_reco = CalHitReco('imcal_barrel_reco',
+ inputHitCollection=imcal_barrel_digi.outputHitCollection,
+ outputHitCollection='EcalBarrelHitsReco',
+ thresholdFactor=5.0,
+ readoutClass='EcalBarrelHits',
+ localDetFields=['system', 'module'], # use local coordinates in each module (stave)
+ **imcal_barrel_daq)
+
+# merge hits to eta, phi bins
+imcal_barrel_merger = CalHitsProj('imcal_barrel_merger',
+ # OutputLevel=DEBUG,
+ inputHitCollection=imcal_barrel_reco.outputHitCollection,
+ outputHitCollection='EcalBarrelHitsMerge',
+ gridSizes=[0.004, 0.004*rad])
+
+imcal_barrel_cl = IslandCluster('imcal_barrel_cl',
+ # OutputLevel=DEBUG,
+ inputHitCollection=imcal_barrel_merger.outputHitCollection,
+ outputHitCollection='EcalBarrelClusterHits',
+ minClusterCenterEdep=1.*MeV,
+ minClusterHitEdep=0.1*MeV,
+ splitCluster=True,
+ globalDistEtaPhi=[0.005, 0.005*rad])
+
+imcal_barrel_clreco = RecoCoG('imcal_barrel_clreco',
+ # OutputLevel=DEBUG,
+ inputHitCollection=imcal_barrel_cl.outputHitCollection,
+ outputClusterCollection='EcalBarrelClusters',
+ logWeightBase=6.2,
+ samplingFraction=kwargs['sf'])
+
+
+podout.outputCommands = ['keep *']
+
+ApplicationMgr(
+ TopAlg=[podin,
+ imcal_barrel_digi, imcal_barrel_reco,
+ imcal_barrel_merger, imcal_barrel_cl, imcal_barrel_clreco,
+ podout],
+ EvtSel='NONE',
+ EvtMax=kwargs['nev'],
+ ExtSvc=[podioevent],
+ OutputLevel=DEBUG
+)
+
+
diff --git a/benchmarks/imaging_ecal/options/scfi_cluster.py b/benchmarks/imaging_ecal/options/scfi_cluster.py
index 6775a0aa921836fed1fdae4097f441786f8a6154..3ab0d5cc1cc427054681c385a234ab38da67dbf1 100644
--- a/benchmarks/imaging_ecal/options/scfi_cluster.py
+++ b/benchmarks/imaging_ecal/options/scfi_cluster.py
@@ -52,6 +52,8 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
outputHitCollection="EcalBarrelScFiHitsReco",
thresholdFactor=5.0,
readoutClass="EcalBarrelScFiHits",
+ layerField="layer",
+ sectorField="module",
localDetFields=["system", "module"], # use local coordinates in each module (stave)
**scfi_barrel_daq)
diff --git a/benchmarks/imaging_ecal/run_imcal_pion0.sh b/benchmarks/imaging_ecal/run_imcal_pion0.sh
new file mode 100644
index 0000000000000000000000000000000000000000..c3c23606b762be5f12b9032be57cfbe53d15248d
--- /dev/null
+++ b/benchmarks/imaging_ecal/run_imcal_pion0.sh
@@ -0,0 +1,125 @@
+#!/bin/bash
+
+print_env.sh
+
+function print_the_help {
+ echo "USAGE: ${0} -n <nevents> -t <nametag> -p <particle> "
+ echo " OPTIONS: "
+ echo " -n,--nevents Number of events"
+ echo " -t,--nametag name tag"
+ echo " -p,--particle particle type"
+ echo " allowed types: pion0, pion+, pion-, kaon0, kaon+, kaon-, proton, neutron, electron, positron, photon"
+ exit
+}
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+ key="$1"
+
+ case $key in
+ -h|--help)
+ shift # past argument
+ print_the_help
+ ;;
+ -t|--nametag)
+ nametag="$2"
+ shift # past argument
+ shift # past value
+ ;;
+ -p|--particle)
+ particle="$2"
+ shift # past argument
+ shift # past value
+ ;;
+ -n|--nevents)
+ export CB_EMCAL_NUMEV="$2"
+ shift # past argument
+ shift # past value
+ ;;
+ *) # unknown option
+ #POSITIONAL+=("$1") # save it in an array for later
+ echo "unknown option $1"
+ print_the_help
+ shift # past argument
+ ;;
+ esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+# export CB_EMCAL_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
+cp benchmarks/imaging_ecal/test_geo.xml ${DETECTOR_PATH}/.
+export CB_EMCAL_COMPACT_PATH=${DETECTOR_PATH}/test_geo.xml
+
+if [[ ! -n "${CB_EMCAL_NUMEV}" ]] ; then
+ export CB_EMCAL_NUMEV=1000
+fi
+
+if [[ ! -n "${CB_EMCAL_ENERGY}" ]] ; then
+ export CB_EMCAL_ENERGY=5.0
+fi
+
+if [[ ! -n "${CB_EMCAL_SAMP_FRAC}" ]] ; then
+ export CB_EMCAL_SAMP_FRAC=0.014
+fi
+
+export CB_EMCAL_NAME_TAG="${nametag}"
+export CB_EMCAL_GEN_FILE="${CB_EMCAL_NAME_TAG}.hepmc"
+
+export CB_EMCAL_SIM_FILE="sim_${CB_EMCAL_NAME_TAG}.root"
+export CB_EMCAL_REC_FILE="rec_${CB_EMCAL_NAME_TAG}.root"
+
+echo "CB_EMCAL_NUMEV = ${CB_EMCAL_NUMEV}"
+echo "CB_EMCAL_COMPACT_PATH = ${CB_EMCAL_COMPACT_PATH}"
+
+# Generate the input events
+python benchmarks/imaging_ecal/scripts/gen_particles.py ${CB_EMCAL_GEN_FILE} -n ${CB_EMCAL_NUMEV}\
+ --angmin 60 --angmax 120 --parray ${CB_EMCAL_ENERGY} --particles="${particle}"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script: generating input events"
+ exit 1
+fi
+
+ls -lh ${CB_EMCAL_GEN_FILE}
+
+# Run geant4 simulations
+npsim --runType batch \
+ -v WARNING \
+ --part.minimalKineticEnergy "1*TeV" \
+ --numberOfEvents ${CB_EMCAL_NUMEV} \
+ --compactFile ${CB_EMCAL_COMPACT_PATH} \
+ --inputFiles ${CB_EMCAL_GEN_FILE} \
+ --outputFile ${CB_EMCAL_SIM_FILE}
+
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running npdet"
+ exit 1
+fi
+
+rootls -t "${CB_EMCAL_SIM_FILE}"
+
+# Directory for plots
+mkdir -p results
+
+CB_EMCAL_OPTION_DIR=benchmarks/imaging_ecal/options
+# Run Juggler
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py ${CB_EMCAL_OPTION_DIR}/imaging_2dcluster.py
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running juggler"
+ exit 1
+fi
+
+# Plot clusters first
+FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts
+python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${CB_EMCAL_SIM_FILE} ${CB_EMCAL_REC_FILE} -o results \
+ --collections "EcalBarrelClusters"
+
+root_filesize=$(stat --format=%s "${CB_EMCAL_REC_FILE}")
+if [[ "${CB_EMCAL_NUMEV}" -lt "500" ]] ; then
+ # file must be less than 10 MB to upload
+ if [[ "${root_filesize}" -lt "10000000" ]] ; then
+ cp ${CB_EMCAL_REC_FILE} results/.
+ fi
+fi
+
diff --git a/benchmarks/imaging_ecal/scripts/get_layerids.py b/benchmarks/imaging_ecal/scripts/get_layerids.py
new file mode 100644
index 0000000000000000000000000000000000000000..b62871dfd8a0f5c165f27f97c0b278244ad78f56
--- /dev/null
+++ b/benchmarks/imaging_ecal/scripts/get_layerids.py
@@ -0,0 +1,136 @@
+'''
+ A simple analysis script to extract some basic info of Monte-Carlo hits
+'''
+import os
+import ROOT
+import pandas as pd
+import numpy as np
+import argparse
+from matplotlib import pyplot as plt
+import matplotlib.ticker as ticker
+from lxml import etree as ET
+
+
+class AthenaDecoder:
+ def __init__(self, compact, readout):
+ self.readouts = self.getReadouts(compact)
+ self.changeReadout(readout)
+
+ def changeReadout(self, readout):
+ self.fieldsmap = self.decomposeIDs(self.readouts[readout])
+
+ def get(self, idvals, field):
+ start, width = self.fieldsmap[field]
+ if width >= 0:
+ return np.bitwise_and(np.right_shift(idvals, start), (1 << width) - 1)
+ # first bit is sign bit
+ else:
+ width = abs(width) - 1
+ vals = np.bitwise_and(np.right_shift(idvals, start), (1 << width) - 1)
+ return np.where(np.bitwise_and(np.right_shift(idvals, start + width), 1), vals - (1 << width), vals)
+
+ def decode(self, idvals):
+ return {field: self.get(idvals, field) for field, _ in self.fieldsmap.items()}
+
+ @staticmethod
+ def getReadouts(path):
+ res = dict()
+ AthenaDecoder.__getReadoutsRecur(path, res)
+ return res
+
+ @staticmethod
+ def __getReadoutsRecur(path, res):
+ if not os.path.exists(path):
+ print('Xml file {} not exist! Ignored it.'.format(path))
+ return
+ lccdd = ET.parse(path).getroot()
+ readouts = lccdd.find('readouts')
+ if readouts is not None:
+ for readout in readouts.getchildren():
+ ids = readout.find('id')
+ if ids is not None:
+ res[readout.attrib['name']] = ids.text
+ for child in lccdd.getchildren():
+ if child.tag == 'include':
+ root_dir = os.path.dirname(os.path.realpath(path))
+ AthenaDecoder.__getReadoutsRecur(os.path.join(root_dir, child.attrib['ref']), res)
+
+ @staticmethod
+ def decomposeIDs(id_str):
+ res = dict()
+ curr_bit = 0
+ for field_bits in id_str.split(','):
+ elements = field_bits.split(':')
+ field_name = elements[0]
+ bit_width = int(elements[-1])
+ if len(elements) == 3:
+ curr_bit = int(elements[1])
+ res[field_name] = (curr_bit, bit_width)
+ curr_bit += abs(bit_width)
+ return res
+
+
+# read from RDataFrame and flatten a given collection, return pandas dataframe
+def flatten_collection(rdf, collection, cols=None):
+ if not cols:
+ cols = [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}.'.format(collection))]
+ else:
+ cols = ['{}.{}'.format(collection, c) for c in cols]
+ if not cols:
+ print('cannot find any branch under collection {}'.format(collection))
+ return pd.DataFrame()
+
+ data = rdf.AsNumpy(cols)
+ # flatten the data, add an event id to identify clusters from different events
+ evns = []
+ for i, vec in enumerate(data[cols[0]]):
+ evns += [i]*vec.size()
+ for n, vals in data.items():
+ # make sure ints are not converted to floats
+ typename = vals[0].__class__.__name__.lower()
+ dtype = np.int64 if 'int' in typename or 'long' in typename else np.float64
+ # type safe creation
+ data[n] = np.asarray([v for vec in vals for v in vec], dtype=dtype)
+ # build data frame
+ dfp = pd.DataFrame({c: pd.Series(v) for c, v in data.items()})
+ dfp.loc[:, 'event'] = evns
+ return dfp
+
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser()
+ parser.add_argument('rec_file', help='Path to reconstruction output file.')
+ parser.add_argument('-o', dest='outdir', default='.', help='Output directory.')
+ parser.add_argument('-c', '--compact', dest='compact', required=True,
+ help='Top-level xml file of the detector description')
+ parser.add_argument('--collection', dest='coll', required=True,
+ help='Hits collection name in the reconstruction file')
+ parser.add_argument('--readout', dest='readout', required=True,
+ help='Readout name for the hits collection')
+ args = parser.parse_args()
+
+ # decoder
+ decoder = AthenaDecoder(args.compact, args.readout)
+
+ # get hits
+ rdf_rec = ROOT.RDataFrame('events', args.rec_file)
+ df = flatten_collection(rdf_rec, args.coll)
+ df.rename(columns={c: c.replace(args.coll + '.', '') for c in df.columns}, inplace=True)
+
+ # initialize dd4hep detector
+ # import DDG4
+ # kernel = DDG4.Kernel()
+ # description = kernel.detectorDescription()
+ # kernel.loadGeometry("file:{}".format(args.compact))
+ # dd4hep_decoder = description.readout(args.readout).idSpec().decoder()
+ # lindex = dd4hep_decoder.index('x')
+ # get_layer_id = np.vectorize(lambda cid: dd4hep_decoder.get(cid, lindex))
+ # df.loc[:, 'layerID'] = get_layer_id(df['cellID'].astype(int).values)
+ # always terminate dd4hep kernel
+ # kernel.terminate()
+
+ # faster way to get layerids
+ df.loc[:, 'layerID'] = decoder.get(df['cellID'].values, 'layer')
+ df.loc[:, 'xID'] = decoder.get(df['cellID'].values, 'x')
+ print(df[['cellID', 'layerID', 'xID', 'position.x', 'position.y', 'position.z', 'energy']])
+
diff --git a/benchmarks/imaging_ecal/test_geo.xml b/benchmarks/imaging_ecal/test_geo.xml
new file mode 100644
index 0000000000000000000000000000000000000000..861dfcd44fe0389bd0f53561c7957f58c62f85c3
--- /dev/null
+++ b/benchmarks/imaging_ecal/test_geo.xml
@@ -0,0 +1,310 @@
+<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
+ xmlns:xs="http://www.w3.org/2001/XMLSchema"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">
+
+ <!-- Some information about detector -->
+ <info name="Athena Detector" title="Athena Detector"
+ author="Athena Collaboration"
+ url="https://eicweb.phy.anl.gov/EIC/detectors/athena.git"
+ status="development"
+ version="v1 2021-03-16">
+ <comment>Athena Detector
+ - https://eicweb.phy.anl.gov/EIC/detectors/athena.git
+ - https://eicweb.phy.anl.gov/EIC/detectors/ip6.git
+ </comment>
+ </info>
+
+ <define>
+ <include ref="ip6/ip6_defs.xml" /> <comment> IP definitions should be first</comment>
+ <include ref="compact/definitions.xml" />
+ </define>
+
+ <properties>
+ <matrix name="RINDEX__Vacuum" coldim="2" values="
+ 1.0*eV 1.0
+ 5.1*eV 1.0
+ "/>
+ <matrix name="RINDEX__Air" coldim="2" values="
+ 1.0*eV 1.00029
+ 5.1*eV 1.00029
+ "/>
+ <matrix name="RINDEX__Quartz" coldim="2" values="
+ 1.0*eV 1.46
+ 5.1*eV 1.46
+ "/>
+ <matrix name="RINDEX__N2" coldim="2" values="
+ 1.0*eV 1.00033
+ 4.0*eV 1.00033
+ 5.1*eV 1.00033
+ "/>
+ <matrix name="RINDEX__Pyrex" coldim="2" values="
+ 1.0*eV 1.5
+ 4.0*eV 1.5
+ 5.1*eV 1.5
+ "/>
+ <matrix name="ABSLENGTH__Pyrex" coldim="2" values="
+ 1.0*eV 10.0*cm
+ 4.0*eV 10.0*cm
+ 5.1*eV 10.0*cm
+ "/>
+ <matrix name= "REFLECTIVITY_mirror" coldim="2" values="
+ 1.0*eV 0.9
+ 4.0*eV 0.9
+ 5.1*eV 0.9
+ "/>
+ <matrix name="RINDEX__Aerogel" coldim="2" values="
+ 1.0*eV 1.030
+ 4.0*eV 1.030
+ 5.1*eV 1.030
+ "/>
+ <matrix name="ABSLENGTH__Aerogel" coldim="2" values="
+ 1.0*eV 4.0*cm
+ 4.0*eV 4.0*cm
+ 5.1*eV 4.0*cm
+ "/>
+ <matrix name="RINDEX__Acrylic" coldim="2" values="
+ 1240*eV/1100 1.49
+ 1240*eV/600 1.49
+ 1240*eV/400 1.49
+ "/>
+ </properties>
+
+ <includes>
+ <gdmlFile ref="compact/elements.xml"/>
+ <gdmlFile ref="compact/materials.xml"/>
+ </includes>
+
+ <materials>
+ <material name="AirOptical">
+ <D type="density" unit="g/cm3" value="0.0012"/>
+ <fraction n="0.754" ref="N"/>
+ <fraction n="0.234" ref="O"/>
+ <fraction n="0.012" ref="Ar"/>
+ <property name="RINDEX" ref="RINDEX__Air"/>
+ <property name="ABSLENGTH" coldim="2" values="1*eV 200*m 5*eV 200*m"/>
+ </material>
+ <material name="N2cherenkov">
+ <D type="density" value="0.00125" unit="g/cm3"/>
+ <composite n="1" ref="N"/>
+ <property name="RINDEX" ref="RINDEX__N2"/>
+ </material>
+ <material name="PyrexGlassOptical">
+ <D type="density" value="2.23" unit="g/cm3"/>
+ <fraction n="0.806" ref="SiliconOxide"/>
+ <fraction n="0.130" ref="BoronOxide"/>
+ <fraction n="0.040" ref="SodiumOxide"/>
+ <fraction n="0.023" ref="AluminumOxide"/>
+ <property name="RINDEX" ref="RINDEX__Pyrex"/>
+ <property name="ABSLENGTH" ref="ABSLENGTH__Pyrex"/>
+ </material>
+ <material name="AerogelOptical">
+ <D value="0.2" unit="g / cm3"/>
+ <fraction n="0.625" ref="SiliconOxide"/>
+ <fraction n="0.374" ref="SiliconOxide"/>
+ <fraction n="0.1" ref="C"/>
+ <property name="RINDEX" ref="RINDEX__Aerogel"/>
+ <property name="ABSLENGTH" ref="ABSLENGTH__Aerogel"/>
+ </material>
+ <material name="AcrylicOptical">
+ <D type="density" value="1.18" unit="g/cm3"/>
+ <composite n="5" ref="C"/>
+ <composite n="2" ref="O"/>
+ <composite n="8" ref="H"/>
+ <property name="RINDEX" ref="RINDEX__Acrylic"/>
+ </material>
+ </materials>
+
+ <surfaces>
+ <comment> For the values of "finish", model and type, see TGeoOpticalSurface.h !
+ </comment>
+
+ <opticalsurface finish="polished" model="glisur" name="MirrorOpticalSurface" type="dielectric_metal" value="0">
+ <property name="REFLECTIVITY" ref="REFLECTIVITY_mirror"/>
+ <property name="RINDEX" coldim="2" values="1.034*eV 1.5 4.136*eV 1.5"/>
+ <!--<property name="EFFICIENCY" ref="EFFICIENCY0x8b77240"/>-->
+ </opticalsurface>
+
+ <!--
+ <opticalsurface name="mirror2" finish="polished" model="glisur" type="dielectric_dielectric">
+ <property name="REFLECTIVITY" coldim="2" values="1.034*eV 0.8 4.136*eV 0.9"/>
+ <property name="EFFICIENCY" coldim="2" values="2.034*eV 0.8 4.136*eV 1.0"/>
+ <property name="RINDEX" coldim="2" values="1.034*eV 1.5 4.136*eV 1.5"/>
+ </opticalsurface>
+ -->
+
+ <opticalsurface finish="polished" model="unified" name="MRICH_MirrorOpticalSurface" type="dielectric_metal" value="0">
+ </opticalsurface>
+ <opticalsurface finish="polished" model="unified" name="MRICH_LensOpticalSurface" type="dielectric_dielectric" value="0">
+ <property name="REFLECTIVITY" coldim="2" values="1240*eV/1100 0.08 1240*eV/400 0.08"/>
+ <!--
+ <property name="RINDEX" coldim="2" values="2.034*eV 1.56 4.136*eV 1.56"/>
+ <property name="SPECULARLOBECONSTANT" coldim="2" values="2.034*eV 0.3 4.136*eV 0.3 "/>
+ <property name="SPECULARSPIKECONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/>
+ <property name="BACKSCATTERCONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/>
+ -->
+ </opticalsurface>
+
+ <opticalsurface finish="polished" model="unified" name="MRICH_PhotoSensorOpticalSurface" type="dielectric_dielectric" value="0">
+ <!--
+ <property name="RINDEX" coldim="2" values="2.034*eV 1.35 4.136*eV 1.40"/>
+ <property name="SPECULARLOBECONSTANT" coldim="2" values="2.034*eV 0.3 4.136*eV 0.3 "/>
+ <property name="SPECULARSPIKECONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/>
+ <property name="BACKSCATTERCONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/>
+ -->
+ </opticalsurface>
+
+ <opticalsurface finish="polished" model="unified" name="MRICH_AerogelOpticalSurface" type="dielectric_dielectric" value="0">
+ <!--
+ <property name="RINDEX" coldim="2" values="2.034*eV 1.010 4.136*eV 1.010"/>
+ <property name="SPECULARLOBECONSTANT" coldim="2" values="2.034*eV 0.3 4.136*eV 0.3 "/>
+ <property name="SPECULARSPIKECONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/>
+ <property name="BACKSCATTERCONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/>
+ -->
+ </opticalsurface>
+
+ </surfaces>
+ <limits>
+ <limitset name="EICBeamlineLimits">
+ <limit name="step_length_max" particles="*" value="1.0" unit="mm" />
+ <limit name="track_length_max" particles="*" value="1.0" unit="mm" />
+ <limit name="time_max" particles="*" value="0.1" unit="ns" />
+ <limit name="ekin_min" particles="*" value="0.001" unit="MeV" />
+ <limit name="range_min" particles="*" value="0.1" unit="mm" />
+ </limitset>
+ <limitset name="cal_limits">
+ <limit name="step_length_max" particles="*" value="5.0" unit="mm"/>
+ </limitset>
+ </limits>
+
+ <display>
+ <include ref="compact/display.xml" />
+ </display>
+
+ <comment> Include the IP components first </comment>
+ <include ref="ip6/forward_ion_beamline.xml"/>
+ <include ref="ip6/beampipe.xml"/>
+
+ <detectors>
+ <detector id="VertexBarrelSubAssembly_ID"
+ name="VertexBarrelSubAssembly"
+ type="DD4hep_SubdetectorAssembly"
+ vis="TrackerSubAssemblyVis">
+ <composite name="VertexBarrel" />
+ </detector>
+ <detector id="VertexEndcapSubAssembly_ID"
+ name="VertexEndcapSubAssembly"
+ type="DD4hep_SubdetectorAssembly"
+ vis="TrackerSubAssemblyVis">
+ <composite name="VertexEndcapN" />
+ <composite name="VertexEndcapP" />
+ </detector>
+
+ <detector id="TrackerBarrelSubAssembly_Inner_ID"
+ name="TrackerBarrelSubAssembly_Inner"
+ type="DD4hep_SubdetectorAssembly"
+ vis="TrackerSubAssemblyVis">
+ <composite name="TrackerBarrel_Inner"/>
+ </detector>
+ <detector id="TrackerEndcapSubAssembly_Inner_ID"
+ name="TrackerEndcapSubAssembly_Inner"
+ type="DD4hep_SubdetectorAssembly"
+ vis="TrackerSubAssemblyVis">
+ <composite name="TrackerEndcapN_Inner"/>
+ <composite name="TrackerEndcapP_Inner"/>
+ </detector>
+
+ <detector id="TrackerBarrelSubAssembly_Outer_ID"
+ name="TrackerSubAssembly_Outer"
+ type="DD4hep_SubdetectorAssembly"
+ vis="TrackerSubAssemblyVis">
+ <composite name="TrackerBarrel_Outer"/>
+ <composite name="TrackerEndcapP_Outer"/>
+ <composite name="TrackerEndcapN_Outer"/>
+ </detector>
+
+ <detector id="TOFSubAssembly_ID"
+ name="TOFSubAssembly"
+ type="DD4hep_SubdetectorAssembly"
+ vis="TOFSubAssemblyVis">
+ <composite name="BarrelTOF"/>
+ <composite name="ForwardTOF"/>
+ <composite name="BackwardTOF"/>
+ </detector>
+
+ <!--
+ <detector id="TrackerBarrelSubAssembly_Inner_ID"
+ name="TrackerBarrelSubAssembly_Inner"
+ type="DD4hep_SubdetectorAssembly"
+ vis="TrackerSubAssemblyVis">
+ <composite name="TrackerBarrel_Inner"/>
+ </detector>
+ <detector id="TrackerEndcapSubAssembly_Inner_ID"
+ name="TrackerEndcapSubAssembly_Inner"
+ type="DD4hep_SubdetectorAssembly"
+ vis="TrackerSubAssemblyVis">
+ <composite name="TrackerEndcapN_Inner"/>
+ <composite name="TrackerBarrel_Inner"/>
+ <composite name="TrackerEndcapP_Inner"/>
+ </detector>
+ <detector id="TrackerBarrelSubAssembly_Outer_ID"
+ name="TrackerBarrelSubAssembly_Outer"
+ type="DD4hep_SubdetectorAssembly"
+ vis="TrackerSubAssemblyVis">
+ <composite name="TrackerBarrel_Outer"/>
+ </detector>
+ <detector id="TrackerEndcapSubAssembly_Outer_ID"
+ name="TrackerEndcapSubAssembly_Outer"
+ type="DD4hep_SubdetectorAssembly"
+ vis="TrackerSubAssemblyVis">
+ <composite name="TrackerEndcapP_Outer"/>
+ <composite name="TrackerEndcapN_Outer"/>
+ </detector>
+ -->
+ </detectors>
+
+ <include ref="compact/vertex_tracker.xml"/>
+
+ <include ref="compact/central_tracker.xml"/>
+
+ <include ref="compact/tof_barrel.xml"/>
+
+ <!--
+ <include ref="compact/rwell_tracker_barrel.xml"/>
+ -->
+ <include ref="compact/cb_DIRC.xml"/>
+
+ <!-- When changing magnet, also select dimensions in definitions.xml. -->
+ <include ref="compact/solenoid.xml"/>
+
+ <include ref="compact/ci_ecal.xml"/>
+ <!--<include ref="compact/ci_ecal_shashlik.xml"/>-->
+ <!--<include ref="compact/ce_ecal.xml"/>-->
+ <include ref="compact/ce_ecal_crystal_glass.xml"/>
+ <!-- <include ref="compact/ecal_barrel.xml"/> -->
+ <!-- <include ref="compact/ecal_barrel_hybrid.xml"/> -->
+ <include ref="compact/ecal_barrel_interlayers.xml"/>
+
+ <include ref="compact/hcal.xml"/>
+
+ <!--include ref="compact/ce_GEM.xml"/-->
+ <!--include ref="compact/gem_tracker_endcap.xml"/-->
+ <include ref="compact/tof_endcap.xml"/>
+ <include ref="compact/forward_trd.xml"/>
+ <include ref="compact/gaseous_rich.xml"/>
+
+ <include ref="ip6/B0_tracker.xml"/>
+ <include ref="ip6/far_forward_offM_tracker.xml"/>
+ <include ref="ip6/far_forward_romanpots.xml"/>
+ <include ref="ip6/far_forward_detectors.xml"/>
+
+
+ <!--
+ <include ref="compact/mm_tracker_barrel.xml"/>
+ <include ref="compact/cb_VTX_Barrel.xml"/>
+ -->
+
+ <readouts>
+ </readouts>
+
+
+</lccdd>