add benchmark for pion0

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]),

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
 

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
+

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",

benchmarks/imaging_ecal/options/imaging_2dcluster.py

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

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)
 

benchmarks/imaging_ecal/run_imcal_pion0.sh

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

benchmarks/imaging_ecal/scripts/get_layerids.py

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

benchmarks/imaging_ecal/test_geo.xml

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