diff --git a/benchmarks/b0_tracker/analysis/b0_tracker_hits.cxx b/benchmarks/b0_tracker/analysis/b0_tracker_hits.cxx
index 721445ba9331c5a34487246b7d93011d46c9ccda..68cb5fa3ebc1709725e85dafc81334aa1cca5ba5 100644
--- a/benchmarks/b0_tracker/analysis/b0_tracker_hits.cxx
+++ b/benchmarks/b0_tracker/analysis/b0_tracker_hits.cxx
@@ -16,16 +16,13 @@ R__LOAD_LIBRARY(libfmt.so)
 #include "TChain.h"
 //#include "TF1.h"
 
-#include "dd4pod/TrackerHitCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
 
 #include "common_bench/particles.h"
 #include "common_bench/benchmark.h"
 #include "common_bench/mt.h"
 #include "common_bench/util.h"
 
-
-#include "dd4pod/TrackerHitCollection.h"
-
 void b0_tracker_hits(const char* fname = "./sim_output/sim_forward_protons.edm4hep.root"){
 
   ROOT::EnableImplicitMT(); // Tell ROOT you want to go parallel
@@ -36,7 +33,7 @@ void b0_tracker_hits(const char* fname = "./sim_output/sim_forward_protons.edm4h
 
   ROOT::RDataFrame d0(*t);
 
-  auto hits_eta = [&](const std::vector<dd4pod::TrackerHitData>& hits) {
+  auto hits_eta = [&](const std::vector<edm4hep::SimTrackerHitData>& hits) {
     std::vector<double> result;
     for (const auto& h : hits) {
       ROOT::Math::XYZVector vec(h.position.x,h.position.y,h.position.z);
@@ -51,8 +48,8 @@ void b0_tracker_hits(const char* fname = "./sim_output/sim_forward_protons.edm4h
   auto h1 = d1.Histo1D({"h1", "hits_eta", 100, 0,20}, "hits_eta");
   TCanvas* c = new TCanvas();
   h1->DrawCopy();
-  c->SaveAs("results/b0_tracker_hits_eta.png");
-  c->SaveAs("results/b0_tracker_hits_eta.pdf");
+  c->SaveAs("results/far_forward/b0/hits_eta.png");
+  c->SaveAs("results/far_forward/b0/hits_eta.pdf");
   auto n1 = h1->GetMean();
   std::cout << "Pseudorapidity of hits: " << n1 << std::endl;
 
diff --git a/benchmarks/b0_tracker/forward_protons.sh b/benchmarks/b0_tracker/forward_protons.sh
index a98c4c5faef12db9f49f995ae30bd5982245cfe2..f92d9db568e3faca7b1268a3890bd25ea3b81aa5 100755
--- a/benchmarks/b0_tracker/forward_protons.sh
+++ b/benchmarks/b0_tracker/forward_protons.sh
@@ -38,7 +38,7 @@ if [[ "$?" -ne "0" ]] ; then
 fi
 
 # Directory for plots
-mkdir -p results
+mkdir -p results/far_forward/b0
 
 rootls -t ${JUGGLER_SIM_FILE}
 # Plot the input events
diff --git a/benchmarks/zdc/ZDC_example.xml b/benchmarks/zdc/ZDC_example.xml
deleted file mode 100644
index e7027926310532431d2eac6db20f89460a6d9fa2..0000000000000000000000000000000000000000
--- a/benchmarks/zdc/ZDC_example.xml
+++ /dev/null
@@ -1,174 +0,0 @@
-<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="ZDC_example" title="Zero Degree Calorimeter detector example"
-        	author="Jihee Kim"
-        	url="https://eicweb.phy.anl.gov/EIC/NPDet"
-        	status="development"
-        	version="v2.0 2020-07-22">
-    	<comment>Zero Degree Calorimeter detector</comment>        
-  	</info>
-
-  	<!-- Use DD4hep elements and materials definitions -->
-  	<includes>
-    		<gdmlFile  ref="elements.xml"/>
-    		<gdmlFile  ref="materials.xml"/>
-  	</includes>
-
-  	<!-- Define the dimensions of the world volume -->
-  	<define>
-		<constant name="world_side" value="50*m"/>
-    		<constant name="world_x" value="world_side"/>
-    		<constant name="world_y" value="world_side"/>
-    		<constant name="world_z" value="world_side"/>
-    
-  	  	<constant name="CrossingAngle" value="0.020*rad"/>
-    		<constant name="tracker_region_zmax" value="5*m"/>
-    		<constant name="tracker_region_rmax" value="5*m"/>
-
-    		<constant name="offset_ZDC"    value="5.0*mm"/>
-		<constant name="st_length"     value="20.0*mm"/>
-		<constant name="lt_length"     value="40.0*mm"/>
-		<constant name="st_ZDC_x_pos"  value="0.0*m"/>   <!-- value="0.60*m"  -->
-		<constant name="st_ZDC_y_pos"  value="0.0*m"/>
-		<constant name="st_ZDC_z_pos"  value="1.0*m"/>   <!-- value="34.0*m"  -->
-		<constant name="lt_ZDC_x_pos"  value="0.0*m"/>   <!-- value="0.60*m"  -->
-		<constant name="lt_ZDC_y_pos"  value="offset_ZDC + (st_length+lt_length)/sqrt(2)"/>
-        	<constant name="lt_ZDC_z_pos"  value="1.0*m"/>  <!-- value="34.0*m"  -->
-  	</define>
-
-  <limits>
-  </limits>
-
-  <regions>
-  </regions>
-
-  	<!-- Common Generic visualization attributes -->
-  	<comment>Common Generic visualization attributes</comment>
-  	<display>
-    	<vis name="InvisibleNoDaughters"      showDaughters="false" visible="false"/>
-    	<vis name="InvisibleWithDaughters"    showDaughters="true" visible="false"/>
-    	<vis name="GreenVis"       alpha="0.5"  r= "0.0" g="1.0" b="0.0" showDaughters="true" visible="true"/>
-    	<vis name="RedVis"         alpha="0.5"  r= "1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/>
-    	<vis name="BlueVis"        alpha="0.5"  r= "0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/>
-    	<vis name="OrangeVis"      alpha="0.5"  r= "1.0" g="0.45" b="0.0" showDaughters="true" visible="true"/>
-    	<vis name="RedGreenVis"    alpha="0.5"  r= "1.0" g="1.0" b="0.0" showDaughters="true" visible="true"/>
-    	<vis name="BlueGreenVis"   alpha="0.5"  r= "0.0" g="1.0" b="1.0" showDaughters="true" visible="true"/>
-    	<vis name="PurpleVis"      alpha="0.5"  r= "1.0" g="0.0" b="1.0" showDaughters="true" visible="true"/>
-    	<vis name="DoubleRedG"     alpha="0.5"  r= "2.0" g=".10" b="0.0" showDaughters="true" visible="true"/>
-    	<vis name="RBG015"         alpha="0.5"  r= "0.0" g=".2"  b="1.0" showDaughters="true" visible="true"/>
-    	<vis name="RBG510"         alpha="0.5"  r= "1.0" g=".2"  b="0.0" showDaughters="true" visible="true"/>
-    	<vis name="RBG"            alpha="0.5"  r= "1.0" g="1.0" b="1.0" showDaughters="true" visible="true"/>
-    	<vis name="GrayVis"        alpha="0.5"  r= "0.75" g="0.75" b="0.75" showDaughters="true" visible="true"/>
-  	</display>
-
-	<!-- Define detector -->
-	<detectors>
-		<detector id="1" name="smallZDC" type="ZDC" readout="ZDCHits" vis="RedVis">
-			<position x="st_ZDC_x_pos" y="st_ZDC_y_pos" z="st_ZDC_z_pos"/>
-			<dimensions x = "st_length" y = "st_length"/>
-			<layer repeat="2">
-				<slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-				<slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-			</layer>
-			<layer repeat="1">
-				<slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-				<slice name="SciFi_belt"     material="PlasticScint"   thickness="1*mm" vis = "RedVis"  sensitive = "true"/>
-				<slice name="SciFi_belt"     material="PlasticScint"   thickness="1*mm" vis = "RedVis"  sensitive = "true"/>
-			</layer>	
-			<layer repeat="2">
-				<slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-				<slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-			</layer>	
-			<layer repeat="2">
-				<slice name="SciFi_belt"     material="PlasticScint"   thickness="1*mm" vis = "RedVis"  sensitive = "true"/>
-			</layer>	
-			<layer repeat="7">
-                                <slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-				<slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-			</layer>
-			<layer repeat="1">
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-				<slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-				<slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-				<slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-			</layer>
-			<layer repeat="2">
-                                <slice name="SciFi_belt"     material="PlasticScint"   thickness="1*mm" vis = "RedVis"  sensitive = "true"/>
-                        </layer>
-                        <layer repeat="3">
-                                <slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                        </layer>
-                        <layer repeat="2">
-                                <slice name="SciFi_belt"     material="PlasticScint"   thickness="1*mm" vis = "RedVis"  sensitive = "true"/>
-                        </layer>
-			<layer repeat="1">
-				<slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-			</layer>
-		</detector>
-                <detector id="2" name="largeZDC" type="ZDC" readout="ZDCHits" vis="RedVis">
-                        <position x="lt_ZDC_x_pos" y="lt_ZDC_y_pos" z="lt_ZDC_z_pos"/>
-                        <dimensions x = "lt_length" y = "lt_length"/>
-                        <layer repeat="2">
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                                <slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-                        </layer>
-                        <layer repeat="1">
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                                <slice name="SciFi_belt"     material="PlasticScint"   thickness="1*mm" vis = "RedVis"  sensitive = "true"/>
-                                <slice name="SciFi_belt"     material="PlasticScint"   thickness="1*mm" vis = "RedVis"  sensitive = "true"/>
-                        </layer>
-                        <layer repeat="2">
-                                <slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                        </layer>
-                        <layer repeat="2">
-                                <slice name="SciFi_belt"     material="PlasticScint"   thickness="1*mm" vis = "RedVis"  sensitive = "true"/>
-                        </layer>
-                        <layer repeat="7">
-                                <slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                        </layer>
-                        <layer repeat="1">
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                                <slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                        </layer>
-                        <layer repeat="2">
-                                <slice name="SciFi_belt"     material="PlasticScint"   thickness="1*mm" vis = "RedVis"  sensitive = "true"/>
-                        </layer>
-                        <layer repeat="3">
-                                <slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                        </layer>
-                        <layer repeat="2">
-                                <slice name="SciFi_belt"     material="PlasticScint"   thickness="1*mm" vis = "RedVis"  sensitive = "true"/>
-                        </layer>
-                        <layer repeat="1">
-                                <slice name="Scint_slice"    material="PlasticScint"   thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
-                                <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
-                        </layer>
-                </detector>			
-	</detectors>
-
-  <!--  Definition of the readout segmentation/definition  -->
-  <readouts>
-        <readout name="ZDCHits">
-		<segmentation type="CartesianGridXY" grid_size_x="1.0*mm" grid_size_y="1.0*mm" />
-		<id>system:8,layer:12,slice:12,x:48:-8,y:-8</id>  
-	</readout>
-  </readouts>
-
-  <plugins>
-  </plugins>
-
-  <fields>
-  </fields>
-</lccdd>
diff --git a/benchmarks/zdc/config.yml b/benchmarks/zdc/config.yml
index f542649595ff08ee0a4e89e13d9b27a70fc90406..eb72a215a2e70f7d6c229ebaf6dc77cfbb31b49d 100644
--- a/benchmarks/zdc/config.yml
+++ b/benchmarks/zdc/config.yml
@@ -2,7 +2,10 @@ sim:zdc:
   extends: .det_benchmark
   stage: simulate
   script:
-    - bash benchmarks/zdc/run_zdc_neutrons.sh
+    - bash benchmarks/zdc/run_zdc_particles.sh --particle $PARTICLE
+  parallel:
+    matrix:
+      - PARTICLE: ["neutron", "photon"]
 
 bench:zdc_benchmark:
   extends: .det_benchmark
@@ -10,7 +13,10 @@ bench:zdc_benchmark:
   needs: 
     - ["sim:zdc"]
   script:
-    - root -b -q benchmarks/zdc/scripts/zdc_neutrons_analysis.cxx+
+    - root -b -q benchmarks/zdc/scripts/zdc_analysis.cxx+
+  parallel:
+    matrix:
+      - PARTICLE: ["neutron", "photon"]
 
 collect_results:zdc:
   extends: .det_benchmark
diff --git a/benchmarks/zdc/elements.xml b/benchmarks/zdc/elements.xml
deleted file mode 100644
index 3fb9b9dd965f5f401c902505e82bfc478208d8a9..0000000000000000000000000000000000000000
--- a/benchmarks/zdc/elements.xml
+++ /dev/null
@@ -1,885 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<materials>
-  <element Z="89" formula="Ac" name="Ac">
-    <atom type="A" unit="g/mol" value="227.028"/>
-  </element>
-  <material formula="Ac" name="Actinium" state="solid">
-    <RL type="X0" unit="cm" value="0.601558"/>
-    <NIL type="lambda" unit="cm" value="21.2048"/>
-    <D type="density" unit="g/cm3" value="10.07"/>
-    <composite n="1" ref="Ac"/>
-  </material>
-  <element Z="47" formula="Ag" name="Ag">
-    <atom type="A" unit="g/mol" value="107.868"/>
-  </element>
-  <material formula="Ag" name="Silver" state="solid">
-    <RL type="X0" unit="cm" value="0.854292"/>
-    <NIL type="lambda" unit="cm" value="15.8546"/>
-    <D type="density" unit="g/cm3" value="10.5"/>
-    <composite n="1" ref="Ag"/>
-  </material>
-  <element Z="13" formula="Al" name="Al">
-    <atom type="A" unit="g/mol" value="26.9815"/>
-  </element>
-  <material formula="Al" name="Aluminum" state="solid">
-    <RL type="X0" unit="cm" value="8.89632"/>
-    <NIL type="lambda" unit="cm" value="38.8766"/>
-    <D type="density" unit="g/cm3" value="2.699"/>
-    <composite n="1" ref="Al"/>
-  </material>
-  <element Z="95" formula="Am" name="Am">
-    <atom type="A" unit="g/mol" value="243.061"/>
-  </element>
-  <material formula="Am" name="Americium" state="solid">
-    <RL type="X0" unit="cm" value="0.42431"/>
-    <NIL type="lambda" unit="cm" value="15.9812"/>
-    <D type="density" unit="g/cm3" value="13.67"/>
-    <composite n="1" ref="Am"/>
-  </material>
-  <element Z="18" formula="Ar" name="Ar">
-    <atom type="A" unit="g/mol" value="39.9477"/>
-  </element>
-  <material formula="Ar" name="Argon" state="gas">
-    <RL type="X0" unit="cm" value="11762.1"/>
-    <NIL type="lambda" unit="cm" value="71926"/>
-    <D type="density" unit="g/cm3" value="0.00166201"/>
-    <composite n="1" ref="Ar"/>
-  </material>
-  <element Z="33" formula="As" name="As">
-    <atom type="A" unit="g/mol" value="74.9216"/>
-  </element>
-  <material formula="As" name="Arsenic" state="solid">
-    <RL type="X0" unit="cm" value="2.0838"/>
-    <NIL type="lambda" unit="cm" value="25.7324"/>
-    <D type="density" unit="g/cm3" value="5.73"/>
-    <composite n="1" ref="As"/>
-  </material>
-  <element Z="85" formula="At" name="At">
-    <atom type="A" unit="g/mol" value="209.987"/>
-  </element>
-  <material formula="At" name="Astatine" state="solid">
-    <RL type="X0" unit="cm" value="0.650799"/>
-    <NIL type="lambda" unit="cm" value="22.3202"/>
-    <D type="density" unit="g/cm3" value="9.32"/>
-    <composite n="1" ref="At"/>
-  </material>
-  <element Z="79" formula="Au" name="Au">
-    <atom type="A" unit="g/mol" value="196.967"/>
-  </element>
-  <material formula="Au" name="Gold" state="solid">
-    <RL type="X0" unit="cm" value="0.334436"/>
-    <NIL type="lambda" unit="cm" value="10.5393"/>
-    <D type="density" unit="g/cm3" value="19.32"/>
-    <composite n="1" ref="Au"/>
-  </material>
-  <element Z="5" formula="B" name="B">
-    <atom type="A" unit="g/mol" value="10.811"/>
-  </element>
-  <material formula="B" name="Boron" state="solid">
-    <RL type="X0" unit="cm" value="22.2307"/>
-    <NIL type="lambda" unit="cm" value="32.2793"/>
-    <D type="density" unit="g/cm3" value="2.37"/>
-    <composite n="1" ref="B"/>
-  </material>
-  <element Z="56" formula="Ba" name="Ba">
-    <atom type="A" unit="g/mol" value="137.327"/>
-  </element>
-  <material formula="Ba" name="Barium" state="solid">
-    <RL type="X0" unit="cm" value="2.37332"/>
-    <NIL type="lambda" unit="cm" value="51.6743"/>
-    <D type="density" unit="g/cm3" value="3.5"/>
-    <composite n="1" ref="Ba"/>
-  </material>
-  <element Z="4" formula="Be" name="Be">
-    <atom type="A" unit="g/mol" value="9.01218"/>
-  </element>
-  <material formula="Be" name="Beryllium" state="solid">
-    <RL type="X0" unit="cm" value="35.276"/>
-    <NIL type="lambda" unit="cm" value="39.4488"/>
-    <D type="density" unit="g/cm3" value="1.848"/>
-    <composite n="1" ref="Be"/>
-  </material>
-  <element Z="83" formula="Bi" name="Bi">
-    <atom type="A" unit="g/mol" value="208.98"/>
-  </element>
-  <material formula="Bi" name="Bismuth" state="solid">
-    <RL type="X0" unit="cm" value="0.645388"/>
-    <NIL type="lambda" unit="cm" value="21.3078"/>
-    <D type="density" unit="g/cm3" value="9.747"/>
-    <composite n="1" ref="Bi"/>
-  </material>
-  <element Z="97" formula="Bk" name="Bk">
-    <atom type="A" unit="g/mol" value="247.07"/>
-  </element>
-  <material formula="Bk" name="Berkelium" state="solid">
-    <RL type="X0" unit="cm" value="0.406479"/>
-    <NIL type="lambda" unit="cm" value="15.6902"/>
-    <D type="density" unit="g/cm3" value="14"/>
-    <composite n="1" ref="Bk"/>
-  </material>
-  <element Z="35" formula="Br" name="Br">
-    <atom type="A" unit="g/mol" value="79.9035"/>
-  </element>
-  <material formula="Br" name="Bromine" state="gas">
-    <RL type="X0" unit="cm" value="1615.12"/>
-    <NIL type="lambda" unit="cm" value="21299"/>
-    <D type="density" unit="g/cm3" value="0.0070721"/>
-    <composite n="1" ref="Br"/>
-  </material>
-  <element Z="6" formula="C" name="C">
-    <atom type="A" unit="g/mol" value="12.0107"/>
-  </element>
-  <material formula="C" name="Carbon" state="solid">
-    <RL type="X0" unit="cm" value="21.3485"/>
-    <NIL type="lambda" unit="cm" value="40.1008"/>
-    <D type="density" unit="g/cm3" value="2"/>
-    <composite n="1" ref="C"/>
-  </material>
-  <element Z="20" formula="Ca" name="Ca">
-    <atom type="A" unit="g/mol" value="40.078"/>
-  </element>
-  <material formula="Ca" name="Calcium" state="solid">
-    <RL type="X0" unit="cm" value="10.4151"/>
-    <NIL type="lambda" unit="cm" value="77.3754"/>
-    <D type="density" unit="g/cm3" value="1.55"/>
-    <composite n="1" ref="Ca"/>
-  </material>
-  <element Z="48" formula="Cd" name="Cd">
-    <atom type="A" unit="g/mol" value="112.411"/>
-  </element>
-  <material formula="Cd" name="Cadmium" state="solid">
-    <RL type="X0" unit="cm" value="1.03994"/>
-    <NIL type="lambda" unit="cm" value="19.46"/>
-    <D type="density" unit="g/cm3" value="8.65"/>
-    <composite n="1" ref="Cd"/>
-  </material>
-  <element Z="58" formula="Ce" name="Ce">
-    <atom type="A" unit="g/mol" value="140.115"/>
-  </element>
-  <material formula="Ce" name="Cerium" state="solid">
-    <RL type="X0" unit="cm" value="1.19506"/>
-    <NIL type="lambda" unit="cm" value="27.3227"/>
-    <D type="density" unit="g/cm3" value="6.657"/>
-    <composite n="1" ref="Ce"/>
-  </material>
-  <element Z="98" formula="Cf" name="Cf">
-    <atom type="A" unit="g/mol" value="251.08"/>
-  </element>
-  <material formula="Cf" name="Californium" state="solid">
-    <RL type="X0" unit="cm" value="0.568328"/>
-    <NIL type="lambda" unit="cm" value="22.085"/>
-    <D type="density" unit="g/cm3" value="10"/>
-    <composite n="1" ref="Cf"/>
-  </material>
-  <element Z="17" formula="Cl" name="Cl">
-    <atom type="A" unit="g/mol" value="35.4526"/>
-  </element>
-  <material formula="Cl" name="Chlorine" state="gas">
-    <RL type="X0" unit="cm" value="6437.34"/>
-    <NIL type="lambda" unit="cm" value="38723.9"/>
-    <D type="density" unit="g/cm3" value="0.00299473"/>
-    <composite n="1" ref="Cl"/>
-  </material>
-  <element Z="96" formula="Cm" name="Cm">
-    <atom type="A" unit="g/mol" value="247.07"/>
-  </element>
-  <material formula="Cm" name="Curium" state="solid">
-    <RL type="X0" unit="cm" value="0.428706"/>
-    <NIL type="lambda" unit="cm" value="16.2593"/>
-    <D type="density" unit="g/cm3" value="13.51"/>
-    <composite n="1" ref="Cm"/>
-  </material>
-  <element Z="27" formula="Co" name="Co">
-    <atom type="A" unit="g/mol" value="58.9332"/>
-  </element>
-  <material formula="Co" name="Cobalt" state="solid">
-    <RL type="X0" unit="cm" value="1.53005"/>
-    <NIL type="lambda" unit="cm" value="15.2922"/>
-    <D type="density" unit="g/cm3" value="8.9"/>
-    <composite n="1" ref="Co"/>
-  </material>
-  <element Z="24" formula="Cr" name="Cr">
-    <atom type="A" unit="g/mol" value="51.9961"/>
-  </element>
-  <material formula="Cr" name="Chromium" state="solid">
-    <RL type="X0" unit="cm" value="2.0814"/>
-    <NIL type="lambda" unit="cm" value="18.1933"/>
-    <D type="density" unit="g/cm3" value="7.18"/>
-    <composite n="1" ref="Cr"/>
-  </material>
-  <element Z="55" formula="Cs" name="Cs">
-    <atom type="A" unit="g/mol" value="132.905"/>
-  </element>
-  <material formula="Cs" name="Cesium" state="solid">
-    <RL type="X0" unit="cm" value="4.4342"/>
-    <NIL type="lambda" unit="cm" value="95.317"/>
-    <D type="density" unit="g/cm3" value="1.873"/>
-    <composite n="1" ref="Cs"/>
-  </material>
-  <element Z="29" formula="Cu" name="Cu">
-    <atom type="A" unit="g/mol" value="63.5456"/>
-  </element>
-  <material formula="Cu" name="Copper" state="solid">
-    <RL type="X0" unit="cm" value="1.43558"/>
-    <NIL type="lambda" unit="cm" value="15.5141"/>
-    <D type="density" unit="g/cm3" value="8.96"/>
-    <composite n="1" ref="Cu"/>
-  </material>
-  <element Z="66" formula="Dy" name="Dy">
-    <atom type="A" unit="g/mol" value="162.497"/>
-  </element>
-  <material formula="Dy" name="Dysprosium" state="solid">
-    <RL type="X0" unit="cm" value="0.85614"/>
-    <NIL type="lambda" unit="cm" value="22.2923"/>
-    <D type="density" unit="g/cm3" value="8.55"/>
-    <composite n="1" ref="Dy"/>
-  </material>
-  <element Z="68" formula="Er" name="Er">
-    <atom type="A" unit="g/mol" value="167.256"/>
-  </element>
-  <material formula="Er" name="Erbium" state="solid">
-    <RL type="X0" unit="cm" value="0.788094"/>
-    <NIL type="lambda" unit="cm" value="21.2923"/>
-    <D type="density" unit="g/cm3" value="9.066"/>
-    <composite n="1" ref="Er"/>
-  </material>
-  <element Z="63" formula="Eu" name="Eu">
-    <atom type="A" unit="g/mol" value="151.964"/>
-  </element>
-  <material formula="Eu" name="Europium" state="solid">
-    <RL type="X0" unit="cm" value="1.41868"/>
-    <NIL type="lambda" unit="cm" value="35.6178"/>
-    <D type="density" unit="g/cm3" value="5.243"/>
-    <composite n="1" ref="Eu"/>
-  </material>
-  <element Z="9" formula="F" name="F">
-    <atom type="A" unit="g/mol" value="18.9984"/>
-  </element>
-  <material formula="F" name="Fluorine" state="gas">
-    <RL type="X0" unit="cm" value="20838.2"/>
-    <NIL type="lambda" unit="cm" value="59094.3"/>
-    <D type="density" unit="g/cm3" value="0.00158029"/>
-    <composite n="1" ref="F"/>
-  </material>
-  <element Z="26" formula="Fe" name="Fe">
-    <atom type="A" unit="g/mol" value="55.8451"/>
-  </element>
-  <material formula="Fe" name="Iron" state="solid">
-    <RL type="X0" unit="cm" value="1.75749"/>
-    <NIL type="lambda" unit="cm" value="16.959"/>
-    <D type="density" unit="g/cm3" value="7.874"/>
-    <composite n="1" ref="Fe"/>
-  </material>
-  <element Z="87" formula="Fr" name="Fr">
-    <atom type="A" unit="g/mol" value="223.02"/>
-  </element>
-  <material formula="Fr" name="Francium" state="solid">
-    <RL type="X0" unit="cm" value="6.18826"/>
-    <NIL type="lambda" unit="cm" value="212.263"/>
-    <D type="density" unit="g/cm3" value="1"/>
-    <composite n="1" ref="Fr"/>
-  </material>
-  <element Z="31" formula="Ga" name="Ga">
-    <atom type="A" unit="g/mol" value="69.7231"/>
-  </element>
-  <material formula="Ga" name="Gallium" state="solid">
-    <RL type="X0" unit="cm" value="2.1128"/>
-    <NIL type="lambda" unit="cm" value="24.3351"/>
-    <D type="density" unit="g/cm3" value="5.904"/>
-    <composite n="1" ref="Ga"/>
-  </material>
-  <element Z="64" formula="Gd" name="Gd">
-    <atom type="A" unit="g/mol" value="157.252"/>
-  </element>
-  <material formula="Gd" name="Gadolinium" state="solid">
-    <RL type="X0" unit="cm" value="0.947208"/>
-    <NIL type="lambda" unit="cm" value="23.9377"/>
-    <D type="density" unit="g/cm3" value="7.9004"/>
-    <composite n="1" ref="Gd"/>
-  </material>
-  <element Z="32" formula="Ge" name="Ge">
-    <atom type="A" unit="g/mol" value="72.6128"/>
-  </element>
-  <material formula="Ge" name="Germanium" state="solid">
-    <RL type="X0" unit="cm" value="2.3013"/>
-    <NIL type="lambda" unit="cm" value="27.3344"/>
-    <D type="density" unit="g/cm3" value="5.323"/>
-    <composite n="1" ref="Ge"/>
-  </material>
-  <element Z="1" formula="H" name="H">
-    <atom type="A" unit="g/mol" value="1.00794"/>
-  </element>
-  <material formula="H" name="Hydrogen" state="gas">
-    <RL type="X0" unit="cm" value="752776"/>
-    <NIL type="lambda" unit="cm" value="421239"/>
-    <D type="density" unit="g/cm3" value="8.3748e-05"/>
-    <composite n="1" ref="H"/>
-  </material>
-  <element Z="2" formula="He" name="He">
-    <atom type="A" unit="g/mol" value="4.00264"/>
-  </element>
-  <material formula="He" name="Helium" state="gas">
-    <RL type="X0" unit="cm" value="567113"/>
-    <NIL type="lambda" unit="cm" value="334266"/>
-    <D type="density" unit="g/cm3" value="0.000166322"/>
-    <composite n="1" ref="He"/>
-  </material>
-  <element Z="72" formula="Hf" name="Hf">
-    <atom type="A" unit="g/mol" value="178.485"/>
-  </element>
-  <material formula="Hf" name="Hafnium" state="solid">
-    <RL type="X0" unit="cm" value="0.517717"/>
-    <NIL type="lambda" unit="cm" value="14.7771"/>
-    <D type="density" unit="g/cm3" value="13.31"/>
-    <composite n="1" ref="Hf"/>
-  </material>
-  <element Z="80" formula="Hg" name="Hg">
-    <atom type="A" unit="g/mol" value="200.599"/>
-  </element>
-  <material formula="Hg" name="Mercury" state="solid">
-    <RL type="X0" unit="cm" value="0.475241"/>
-    <NIL type="lambda" unit="cm" value="15.105"/>
-    <D type="density" unit="g/cm3" value="13.546"/>
-    <composite n="1" ref="Hg"/>
-  </material>
-  <element Z="67" formula="Ho" name="Ho">
-    <atom type="A" unit="g/mol" value="164.93"/>
-  </element>
-  <material formula="Ho" name="Holmium" state="solid">
-    <RL type="X0" unit="cm" value="0.822447"/>
-    <NIL type="lambda" unit="cm" value="21.8177"/>
-    <D type="density" unit="g/cm3" value="8.795"/>
-    <composite n="1" ref="Ho"/>
-  </material>
-  <element Z="53" formula="I" name="I">
-    <atom type="A" unit="g/mol" value="126.904"/>
-  </element>
-  <material formula="I" name="Iodine" state="solid">
-    <RL type="X0" unit="cm" value="1.72016"/>
-    <NIL type="lambda" unit="cm" value="35.6583"/>
-    <D type="density" unit="g/cm3" value="4.93"/>
-    <composite n="1" ref="I"/>
-  </material>
-  <element Z="49" formula="In" name="In">
-    <atom type="A" unit="g/mol" value="114.818"/>
-  </element>
-  <material formula="In" name="Indium" state="solid">
-    <RL type="X0" unit="cm" value="1.21055"/>
-    <NIL type="lambda" unit="cm" value="23.2468"/>
-    <D type="density" unit="g/cm3" value="7.31"/>
-    <composite n="1" ref="In"/>
-  </material>
-  <element Z="77" formula="Ir" name="Ir">
-    <atom type="A" unit="g/mol" value="192.216"/>
-  </element>
-  <material formula="Ir" name="Iridium" state="solid">
-    <RL type="X0" unit="cm" value="0.294142"/>
-    <NIL type="lambda" unit="cm" value="9.01616"/>
-    <D type="density" unit="g/cm3" value="22.42"/>
-    <composite n="1" ref="Ir"/>
-  </material>
-  <element Z="19" formula="K" name="K">
-    <atom type="A" unit="g/mol" value="39.0983"/>
-  </element>
-  <material formula="K" name="Potassium" state="solid">
-    <RL type="X0" unit="cm" value="20.0871"/>
-    <NIL type="lambda" unit="cm" value="138.041"/>
-    <D type="density" unit="g/cm3" value="0.862"/>
-    <composite n="1" ref="K"/>
-  </material>
-  <element Z="36" formula="Kr" name="Kr">
-    <atom type="A" unit="g/mol" value="83.7993"/>
-  </element>
-  <material formula="Kr" name="Krypton" state="gas">
-    <RL type="X0" unit="cm" value="3269.44"/>
-    <NIL type="lambda" unit="cm" value="43962.9"/>
-    <D type="density" unit="g/cm3" value="0.00347832"/>
-    <composite n="1" ref="Kr"/>
-  </material>
-  <element Z="57" formula="La" name="La">
-    <atom type="A" unit="g/mol" value="138.905"/>
-  </element>
-  <material formula="La" name="Lanthanum" state="solid">
-    <RL type="X0" unit="cm" value="1.32238"/>
-    <NIL type="lambda" unit="cm" value="29.441"/>
-    <D type="density" unit="g/cm3" value="6.154"/>
-    <composite n="1" ref="La"/>
-  </material>
-  <element Z="3" formula="Li" name="Li">
-    <atom type="A" unit="g/mol" value="6.94003"/>
-  </element>
-  <material formula="Li" name="Lithium" state="solid">
-    <RL type="X0" unit="cm" value="154.997"/>
-    <NIL type="lambda" unit="cm" value="124.305"/>
-    <D type="density" unit="g/cm3" value="0.534"/>
-    <composite n="1" ref="Li"/>
-  </material>
-  <element Z="71" formula="Lu" name="Lu">
-    <atom type="A" unit="g/mol" value="174.967"/>
-  </element>
-  <material formula="Lu" name="Lutetium" state="solid">
-    <RL type="X0" unit="cm" value="0.703651"/>
-    <NIL type="lambda" unit="cm" value="19.8916"/>
-    <D type="density" unit="g/cm3" value="9.84"/>
-    <composite n="1" ref="Lu"/>
-  </material>
-  <element Z="12" formula="Mg" name="Mg">
-    <atom type="A" unit="g/mol" value="24.305"/>
-  </element>
-  <material formula="Mg" name="Magnesium" state="solid">
-    <RL type="X0" unit="cm" value="14.3859"/>
-    <NIL type="lambda" unit="cm" value="58.7589"/>
-    <D type="density" unit="g/cm3" value="1.74"/>
-    <composite n="1" ref="Mg"/>
-  </material>
-  <element Z="25" formula="Mn" name="Mn">
-    <atom type="A" unit="g/mol" value="54.938"/>
-  </element>
-  <material formula="Mn" name="Manganese" state="solid">
-    <RL type="X0" unit="cm" value="1.96772"/>
-    <NIL type="lambda" unit="cm" value="17.8701"/>
-    <D type="density" unit="g/cm3" value="7.44"/>
-    <composite n="1" ref="Mn"/>
-  </material>
-  <element Z="42" formula="Mo" name="Mo">
-    <atom type="A" unit="g/mol" value="95.9313"/>
-  </element>
-  <material formula="Mo" name="Molybdenum" state="solid">
-    <RL type="X0" unit="cm" value="0.959107"/>
-    <NIL type="lambda" unit="cm" value="15.6698"/>
-    <D type="density" unit="g/cm3" value="10.22"/>
-    <composite n="1" ref="Mo"/>
-  </material>
-  <element Z="7" formula="N" name="N">
-    <atom type="A" unit="g/mol" value="14.0068"/>
-  </element>
-  <material formula="N" name="Nitrogen" state="gas">
-    <RL type="X0" unit="cm" value="32602.2"/>
-    <NIL type="lambda" unit="cm" value="72430.3"/>
-    <D type="density" unit="g/cm3" value="0.0011652"/>
-    <composite n="1" ref="N"/>
-  </material>
-  <element Z="11" formula="Na" name="Na">
-    <atom type="A" unit="g/mol" value="22.9898"/>
-  </element>
-  <material formula="Na" name="Sodium" state="solid">
-    <RL type="X0" unit="cm" value="28.5646"/>
-    <NIL type="lambda" unit="cm" value="102.463"/>
-    <D type="density" unit="g/cm3" value="0.971"/>
-    <composite n="1" ref="Na"/>
-  </material>
-  <element Z="41" formula="Nb" name="Nb">
-    <atom type="A" unit="g/mol" value="92.9064"/>
-  </element>
-  <material formula="Nb" name="Niobium" state="solid">
-    <RL type="X0" unit="cm" value="1.15783"/>
-    <NIL type="lambda" unit="cm" value="18.4846"/>
-    <D type="density" unit="g/cm3" value="8.57"/>
-    <composite n="1" ref="Nb"/>
-  </material>
-  <element Z="60" formula="Nd" name="Nd">
-    <atom type="A" unit="g/mol" value="144.236"/>
-  </element>
-  <material formula="Nd" name="Neodymium" state="solid">
-    <RL type="X0" unit="cm" value="1.11667"/>
-    <NIL type="lambda" unit="cm" value="26.6308"/>
-    <D type="density" unit="g/cm3" value="6.9"/>
-    <composite n="1" ref="Nd"/>
-  </material>
-  <element Z="10" formula="Ne" name="Ne">
-    <atom type="A" unit="g/mol" value="20.18"/>
-  </element>
-  <material formula="Ne" name="Neon" state="gas">
-    <RL type="X0" unit="cm" value="34504.8"/>
-    <NIL type="lambda" unit="cm" value="114322"/>
-    <D type="density" unit="g/cm3" value="0.000838505"/>
-    <composite n="1" ref="Ne"/>
-  </material>
-  <element Z="28" formula="Ni" name="Ni">
-    <atom type="A" unit="g/mol" value="58.6933"/>
-  </element>
-  <material formula="Ni" name="Nickel" state="solid">
-    <RL type="X0" unit="cm" value="1.42422"/>
-    <NIL type="lambda" unit="cm" value="15.2265"/>
-    <D type="density" unit="g/cm3" value="8.902"/>
-    <composite n="1" ref="Ni"/>
-  </material>
-  <element Z="93" formula="Np" name="Np">
-    <atom type="A" unit="g/mol" value="237.048"/>
-  </element>
-  <material formula="Np" name="Neptunium" state="solid">
-    <RL type="X0" unit="cm" value="0.289676"/>
-    <NIL type="lambda" unit="cm" value="10.6983"/>
-    <D type="density" unit="g/cm3" value="20.25"/>
-    <composite n="1" ref="Np"/>
-  </material>
-  <element Z="8" formula="O" name="O">
-    <atom type="A" unit="g/mol" value="15.9994"/>
-  </element>
-  <material formula="O" name="Oxygen" state="gas">
-    <RL type="X0" unit="cm" value="25713.8"/>
-    <NIL type="lambda" unit="cm" value="66233.9"/>
-    <D type="density" unit="g/cm3" value="0.00133151"/>
-    <composite n="1" ref="O"/>
-  </material>
-  <element Z="76" formula="Os" name="Os">
-    <atom type="A" unit="g/mol" value="190.225"/>
-  </element>
-  <material formula="Os" name="Osmium" state="solid">
-    <RL type="X0" unit="cm" value="0.295861"/>
-    <NIL type="lambda" unit="cm" value="8.92553"/>
-    <D type="density" unit="g/cm3" value="22.57"/>
-    <composite n="1" ref="Os"/>
-  </material>
-  <element Z="15" formula="P" name="P">
-    <atom type="A" unit="g/mol" value="30.9738"/>
-  </element>
-  <material formula="P" name="Phosphorus" state="solid">
-    <RL type="X0" unit="cm" value="9.63879"/>
-    <NIL type="lambda" unit="cm" value="49.9343"/>
-    <D type="density" unit="g/cm3" value="2.2"/>
-    <composite n="1" ref="P"/>
-  </material>
-  <element Z="91" formula="Pa" name="Pa">
-    <atom type="A" unit="g/mol" value="231.036"/>
-  </element>
-  <material formula="Pa" name="Protactinium" state="solid">
-    <RL type="X0" unit="cm" value="0.38607"/>
-    <NIL type="lambda" unit="cm" value="13.9744"/>
-    <D type="density" unit="g/cm3" value="15.37"/>
-    <composite n="1" ref="Pa"/>
-  </material>
-  <element Z="82" formula="Pb" name="Pb">
-    <atom type="A" unit="g/mol" value="207.217"/>
-  </element>
-  <material formula="Pb" name="Lead" state="solid">
-    <RL type="X0" unit="cm" value="0.561253"/>
-    <NIL type="lambda" unit="cm" value="18.2607"/>
-    <D type="density" unit="g/cm3" value="11.35"/>
-    <composite n="1" ref="Pb"/>
-  </material>
-  <element Z="46" formula="Pd" name="Pd">
-    <atom type="A" unit="g/mol" value="106.415"/>
-  </element>
-  <material formula="Pd" name="Palladium" state="solid">
-    <RL type="X0" unit="cm" value="0.765717"/>
-    <NIL type="lambda" unit="cm" value="13.7482"/>
-    <D type="density" unit="g/cm3" value="12.02"/>
-    <composite n="1" ref="Pd"/>
-  </material>
-  <element Z="61" formula="Pm" name="Pm">
-    <atom type="A" unit="g/mol" value="144.913"/>
-  </element>
-  <material formula="Pm" name="Promethium" state="solid">
-    <RL type="X0" unit="cm" value="1.04085"/>
-    <NIL type="lambda" unit="cm" value="25.4523"/>
-    <D type="density" unit="g/cm3" value="7.22"/>
-    <composite n="1" ref="Pm"/>
-  </material>
-  <element Z="84" formula="Po" name="Po">
-    <atom type="A" unit="g/mol" value="208.982"/>
-  </element>
-  <material formula="Po" name="Polonium" state="solid">
-    <RL type="X0" unit="cm" value="0.661092"/>
-    <NIL type="lambda" unit="cm" value="22.2842"/>
-    <D type="density" unit="g/cm3" value="9.32"/>
-    <composite n="1" ref="Po"/>
-  </material>
-  <element Z="59" formula="Pr" name="Pr">
-    <atom type="A" unit="g/mol" value="140.908"/>
-  </element>
-  <material formula="Pr" name="Praseodymium" state="solid">
-    <RL type="X0" unit="cm" value="1.1562"/>
-    <NIL type="lambda" unit="cm" value="27.1312"/>
-    <D type="density" unit="g/cm3" value="6.71"/>
-    <composite n="1" ref="Pr"/>
-  </material>
-  <element Z="78" formula="Pt" name="Pt">
-    <atom type="A" unit="g/mol" value="195.078"/>
-  </element>
-  <material formula="Pt" name="Platinum" state="solid">
-    <RL type="X0" unit="cm" value="0.305053"/>
-    <NIL type="lambda" unit="cm" value="9.46584"/>
-    <D type="density" unit="g/cm3" value="21.45"/>
-    <composite n="1" ref="Pt"/>
-  </material>
-  <element Z="94" formula="Pu" name="Pu">
-    <atom type="A" unit="g/mol" value="244.064"/>
-  </element>
-  <material formula="Pu" name="Plutonium" state="solid">
-    <RL type="X0" unit="cm" value="0.298905"/>
-    <NIL type="lambda" unit="cm" value="11.0265"/>
-    <D type="density" unit="g/cm3" value="19.84"/>
-    <composite n="1" ref="Pu"/>
-  </material>
-  <element Z="88" formula="Ra" name="Ra">
-    <atom type="A" unit="g/mol" value="226.025"/>
-  </element>
-  <material formula="Ra" name="Radium" state="solid">
-    <RL type="X0" unit="cm" value="1.22987"/>
-    <NIL type="lambda" unit="cm" value="42.6431"/>
-    <D type="density" unit="g/cm3" value="5"/>
-    <composite n="1" ref="Ra"/>
-  </material>
-  <element Z="37" formula="Rb" name="Rb">
-    <atom type="A" unit="g/mol" value="85.4677"/>
-  </element>
-  <material formula="Rb" name="Rubidium" state="solid">
-    <RL type="X0" unit="cm" value="7.19774"/>
-    <NIL type="lambda" unit="cm" value="100.218"/>
-    <D type="density" unit="g/cm3" value="1.532"/>
-    <composite n="1" ref="Rb"/>
-  </material>
-  <element Z="75" formula="Re" name="Re">
-    <atom type="A" unit="g/mol" value="186.207"/>
-  </element>
-  <material formula="Re" name="Rhenium" state="solid">
-    <RL type="X0" unit="cm" value="0.318283"/>
-    <NIL type="lambda" unit="cm" value="9.5153"/>
-    <D type="density" unit="g/cm3" value="21.02"/>
-    <composite n="1" ref="Re"/>
-  </material>
-  <element Z="45" formula="Rh" name="Rh">
-    <atom type="A" unit="g/mol" value="102.906"/>
-  </element>
-  <material formula="Rh" name="Rhodium" state="solid">
-    <RL type="X0" unit="cm" value="0.746619"/>
-    <NIL type="lambda" unit="cm" value="13.2083"/>
-    <D type="density" unit="g/cm3" value="12.41"/>
-    <composite n="1" ref="Rh"/>
-  </material>
-  <element Z="86" formula="Rn" name="Rn">
-    <atom type="A" unit="g/mol" value="222.018"/>
-  </element>
-  <material formula="Rn" name="Radon" state="gas">
-    <RL type="X0" unit="cm" value="697.777"/>
-    <NIL type="lambda" unit="cm" value="23532"/>
-    <D type="density" unit="g/cm3" value="0.00900662"/>
-    <composite n="1" ref="Rn"/>
-  </material>
-  <element Z="44" formula="Ru" name="Ru">
-    <atom type="A" unit="g/mol" value="101.065"/>
-  </element>
-  <material formula="Ru" name="Ruthenium" state="solid">
-    <RL type="X0" unit="cm" value="0.764067"/>
-    <NIL type="lambda" unit="cm" value="13.1426"/>
-    <D type="density" unit="g/cm3" value="12.41"/>
-    <composite n="1" ref="Ru"/>
-  </material>
-  <element Z="16" formula="S" name="S">
-    <atom type="A" unit="g/mol" value="32.0661"/>
-  </element>
-  <material formula="S" name="Sulfur" state="solid">
-    <RL type="X0" unit="cm" value="9.74829"/>
-    <NIL type="lambda" unit="cm" value="55.6738"/>
-    <D type="density" unit="g/cm3" value="2"/>
-    <composite n="1" ref="S"/>
-  </material>
-  <element Z="51" formula="Sb" name="Sb">
-    <atom type="A" unit="g/mol" value="121.76"/>
-  </element>
-  <material formula="Sb" name="Antimony" state="solid">
-    <RL type="X0" unit="cm" value="1.30401"/>
-    <NIL type="lambda" unit="cm" value="25.8925"/>
-    <D type="density" unit="g/cm3" value="6.691"/>
-    <composite n="1" ref="Sb"/>
-  </material>
-  <element Z="21" formula="Sc" name="Sc">
-    <atom type="A" unit="g/mol" value="44.9559"/>
-  </element>
-  <material formula="Sc" name="Scandium" state="solid">
-    <RL type="X0" unit="cm" value="5.53545"/>
-    <NIL type="lambda" unit="cm" value="41.609"/>
-    <D type="density" unit="g/cm3" value="2.989"/>
-    <composite n="1" ref="Sc"/>
-  </material>
-  <element Z="34" formula="Se" name="Se">
-    <atom type="A" unit="g/mol" value="78.9594"/>
-  </element>
-  <material formula="Se" name="Selenium" state="solid">
-    <RL type="X0" unit="cm" value="2.64625"/>
-    <NIL type="lambda" unit="cm" value="33.356"/>
-    <D type="density" unit="g/cm3" value="4.5"/>
-    <composite n="1" ref="Se"/>
-  </material>
-  <element Z="14" formula="Si" name="Si">
-    <atom type="A" unit="g/mol" value="28.0854"/>
-  </element>
-  <material formula="Si" name="Silicon" state="solid">
-    <RL type="X0" unit="cm" value="9.36607"/>
-    <NIL type="lambda" unit="cm" value="45.7531"/>
-    <D type="density" unit="g/cm3" value="2.33"/>
-    <composite n="1" ref="Si"/>
-  </material>
-  <element Z="62" formula="Sm" name="Sm">
-    <atom type="A" unit="g/mol" value="150.366"/>
-  </element>
-  <material formula="Sm" name="Samarium" state="solid">
-    <RL type="X0" unit="cm" value="1.01524"/>
-    <NIL type="lambda" unit="cm" value="24.9892"/>
-    <D type="density" unit="g/cm3" value="7.46"/>
-    <composite n="1" ref="Sm"/>
-  </material>
-  <element Z="50" formula="Sn" name="Sn">
-    <atom type="A" unit="g/mol" value="118.71"/>
-  </element>
-  <material formula="Sn" name="Tin" state="solid">
-    <RL type="X0" unit="cm" value="1.20637"/>
-    <NIL type="lambda" unit="cm" value="23.4931"/>
-    <D type="density" unit="g/cm3" value="7.31"/>
-    <composite n="1" ref="Sn"/>
-  </material>
-  <element Z="38" formula="Sr" name="Sr">
-    <atom type="A" unit="g/mol" value="87.6166"/>
-  </element>
-  <material formula="Sr" name="Strontium" state="solid">
-    <RL type="X0" unit="cm" value="4.237"/>
-    <NIL type="lambda" unit="cm" value="61.0238"/>
-    <D type="density" unit="g/cm3" value="2.54"/>
-    <composite n="1" ref="Sr"/>
-  </material>
-  <element Z="73" formula="Ta" name="Ta">
-    <atom type="A" unit="g/mol" value="180.948"/>
-  </element>
-  <material formula="Ta" name="Tantalum" state="solid">
-    <RL type="X0" unit="cm" value="0.409392"/>
-    <NIL type="lambda" unit="cm" value="11.8846"/>
-    <D type="density" unit="g/cm3" value="16.654"/>
-    <composite n="1" ref="Ta"/>
-  </material>
-  <element Z="65" formula="Tb" name="Tb">
-    <atom type="A" unit="g/mol" value="158.925"/>
-  </element>
-  <material formula="Tb" name="Terbium" state="solid">
-    <RL type="X0" unit="cm" value="0.893977"/>
-    <NIL type="lambda" unit="cm" value="23.0311"/>
-    <D type="density" unit="g/cm3" value="8.229"/>
-    <composite n="1" ref="Tb"/>
-  </material>
-  <element Z="43" formula="Tc" name="Tc">
-    <atom type="A" unit="g/mol" value="97.9072"/>
-  </element>
-  <material formula="Tc" name="Technetium" state="solid">
-    <RL type="X0" unit="cm" value="0.833149"/>
-    <NIL type="lambda" unit="cm" value="14.0185"/>
-    <D type="density" unit="g/cm3" value="11.5"/>
-    <composite n="1" ref="Tc"/>
-  </material>
-  <element Z="52" formula="Te" name="Te">
-    <atom type="A" unit="g/mol" value="127.603"/>
-  </element>
-  <material formula="Te" name="Tellurium" state="solid">
-    <RL type="X0" unit="cm" value="1.41457"/>
-    <NIL type="lambda" unit="cm" value="28.1797"/>
-    <D type="density" unit="g/cm3" value="6.24"/>
-    <composite n="1" ref="Te"/>
-  </material>
-  <element Z="90" formula="Th" name="Th">
-    <atom type="A" unit="g/mol" value="232.038"/>
-  </element>
-  <material formula="Th" name="Thorium" state="solid">
-    <RL type="X0" unit="cm" value="0.51823"/>
-    <NIL type="lambda" unit="cm" value="18.353"/>
-    <D type="density" unit="g/cm3" value="11.72"/>
-    <composite n="1" ref="Th"/>
-  </material>
-  <element Z="22" formula="Ti" name="Ti">
-    <atom type="A" unit="g/mol" value="47.8667"/>
-  </element>
-  <material formula="Ti" name="Titanium" state="solid">
-    <RL type="X0" unit="cm" value="3.5602"/>
-    <NIL type="lambda" unit="cm" value="27.9395"/>
-    <D type="density" unit="g/cm3" value="4.54"/>
-    <composite n="1" ref="Ti"/>
-  </material>
-  <element Z="81" formula="Tl" name="Tl">
-    <atom type="A" unit="g/mol" value="204.383"/>
-  </element>
-  <material formula="Tl" name="Thallium" state="solid">
-    <RL type="X0" unit="cm" value="0.547665"/>
-    <NIL type="lambda" unit="cm" value="17.6129"/>
-    <D type="density" unit="g/cm3" value="11.72"/>
-    <composite n="1" ref="Tl"/>
-  </material>
-  <element Z="69" formula="Tm" name="Tm">
-    <atom type="A" unit="g/mol" value="168.934"/>
-  </element>
-  <material formula="Tm" name="Thulium" state="solid">
-    <RL type="X0" unit="cm" value="0.754428"/>
-    <NIL type="lambda" unit="cm" value="20.7522"/>
-    <D type="density" unit="g/cm3" value="9.321"/>
-    <composite n="1" ref="Tm"/>
-  </material>
-  <element Z="92" formula="U" name="U">
-    <atom type="A" unit="g/mol" value="238.029"/>
-  </element>
-  <material formula="U" name="Uranium" state="solid">
-    <RL type="X0" unit="cm" value="0.31663"/>
-    <NIL type="lambda" unit="cm" value="11.4473"/>
-    <D type="density" unit="g/cm3" value="18.95"/>
-    <composite n="1" ref="U"/>
-  </material>
-  <element Z="23" formula="V" name="V">
-    <atom type="A" unit="g/mol" value="50.9415"/>
-  </element>
-  <material formula="V" name="Vanadium" state="solid">
-    <RL type="X0" unit="cm" value="2.59285"/>
-    <NIL type="lambda" unit="cm" value="21.2187"/>
-    <D type="density" unit="g/cm3" value="6.11"/>
-    <composite n="1" ref="V"/>
-  </material>
-  <element Z="74" formula="W" name="W">
-    <atom type="A" unit="g/mol" value="183.842"/>
-  </element>
-  <material formula="W" name="Tungsten" state="solid">
-    <RL type="X0" unit="cm" value="0.350418"/>
-    <NIL type="lambda" unit="cm" value="10.3057"/>
-    <D type="density" unit="g/cm3" value="19.3"/>
-    <composite n="1" ref="W"/>
-  </material>
-  <element Z="54" formula="Xe" name="Xe">
-    <atom type="A" unit="g/mol" value="131.292"/>
-  </element>
-  <material formula="Xe" name="Xenon" state="gas">
-    <RL type="X0" unit="cm" value="1546.2"/>
-    <NIL type="lambda" unit="cm" value="32477.9"/>
-    <D type="density" unit="g/cm3" value="0.00548536"/>
-    <composite n="1" ref="Xe"/>
-  </material>
-  <element Z="39" formula="Y" name="Y">
-    <atom type="A" unit="g/mol" value="88.9058"/>
-  </element>
-  <material formula="Y" name="Yttrium" state="solid">
-    <RL type="X0" unit="cm" value="2.32943"/>
-    <NIL type="lambda" unit="cm" value="34.9297"/>
-    <D type="density" unit="g/cm3" value="4.469"/>
-    <composite n="1" ref="Y"/>
-  </material>
-  <element Z="70" formula="Yb" name="Yb">
-    <atom type="A" unit="g/mol" value="173.038"/>
-  </element>
-  <material formula="Yb" name="Ytterbium" state="solid">
-    <RL type="X0" unit="cm" value="1.04332"/>
-    <NIL type="lambda" unit="cm" value="28.9843"/>
-    <D type="density" unit="g/cm3" value="6.73"/>
-    <composite n="1" ref="Yb"/>
-  </material>
-  <element Z="30" formula="Zn" name="Zn">
-    <atom type="A" unit="g/mol" value="65.3955"/>
-  </element>
-  <material formula="Zn" name="Zinc" state="solid">
-    <RL type="X0" unit="cm" value="1.74286"/>
-    <NIL type="lambda" unit="cm" value="19.8488"/>
-    <D type="density" unit="g/cm3" value="7.133"/>
-    <composite n="1" ref="Zn"/>
-  </material>
-  <element Z="40" formula="Zr" name="Zr">
-    <atom type="A" unit="g/mol" value="91.2236"/>
-  </element>
-  <material formula="Zr" name="Zirconium" state="solid">
-    <RL type="X0" unit="cm" value="1.56707"/>
-    <NIL type="lambda" unit="cm" value="24.2568"/>
-    <D type="density" unit="g/cm3" value="6.506"/>
-    <composite n="1" ref="Zr"/>
-  </material>
-</materials>
diff --git a/benchmarks/zdc/materials.xml b/benchmarks/zdc/materials.xml
deleted file mode 100644
index f6f705f68423d130c840cd19e0d4eeb10f9bbe9d..0000000000000000000000000000000000000000
--- a/benchmarks/zdc/materials.xml
+++ /dev/null
@@ -1,189 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<materials>
-  <!--
-       Air by weight from
-
-       http://www.engineeringtoolbox.com/air-composition-24_212.html
-  -->
-  <material name="Air">
-    <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"/>
-  </material>
-  <!-- We model vakuum just as very thin air -->
-  <material name="Vacuum">
-    <D type="density" unit="g/cm3" value="0.0000000001"/>
-    <fraction n="0.754" ref="N"/>
-    <fraction n="0.234" ref="O"/>
-    <fraction n="0.012" ref="Ar"/>
-  </material>
-  <material name="Epoxy">
-    <D type="density" value="1.3" unit="g/cm3"/>
-    <composite n="44" ref="H"/>
-    <composite n="15" ref="C"/>
-    <composite n="7" ref="O"/>
-  </material>
-  <material name="Quartz">
-    <D type="density" value="2.2" unit="g/cm3"/>
-    <composite n="1" ref="Si"/>
-    <composite n="2" ref="O"/>
-  </material>
-  <material name="G10">
-    <D type="density" value="1.7" unit="g/cm3"/>
-    <fraction n="0.08" ref="Cl"/>
-    <fraction n="0.773" ref="Quartz"/>
-    <fraction n="0.147" ref="Epoxy"/>
-  </material>
-  <material name="Polystyrene">
-    <D value="1.032" unit="g/cm3"/>
-    <composite n="19" ref="C"/>
-    <composite n="21" ref="H"/>
-  </material>
-  <material name="Steel235">
-    <D value="7.85" unit="g/cm3"/>
-    <fraction n="0.998" ref="Fe"/>
-    <fraction n=".002" ref="C"/>
-  </material>
-  <material name="SiliconOxide">
-    <D type="density" value="2.65" unit="g/cm3"/>
-    <composite n="1" ref="Si"/>
-    <composite n="2" ref="O"/>
-  </material>
-  <material name="BoronOxide">
-    <D type="density" value="2.46" unit="g/cm3"/>
-    <composite n="2" ref="B"/>
-    <composite n="3" ref="O"/>
-  </material>
-  <material name="SodiumOxide">
-    <D type="density" value="2.65" unit="g/cm3"/>
-    <composite n="2" ref="Na"/>
-    <composite n="1" ref="O"/>
-  </material>
-  <material name="AluminumOxide">
-    <D type="density" value="3.89" unit="g/cm3"/>
-    <composite n="2" ref="Al"/>
-    <composite n="3" ref="O"/>
-  </material>
-  <material name="SiliconNitride">
-    <D type="density" value="3.17" unit="g/cm3"/>
-    <composite n="3" ref="Si"/>
-    <composite n="4" ref="N"/>
-  </material>
-  <material name="PyrexGlass">
-    <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"/>
-  </material>
-  <material name="CarbonFiber">
-    <D type="density" value="1.5" unit="g/cm3"/>
-    <fraction n="0.65" ref="C"/>
-    <fraction n="0.35" ref="Epoxy"/>
-  </material>
-  <material name="CarbonFiber_50D">
-    <D type="density" value="0.75" unit="g/cm3"/>
-    <fraction n="0.65" ref="C"/>
-    <fraction n="0.35" ref="Epoxy"/>
-  </material>
-  <material name="Rohacell31">
-    <D type="density" value="0.032" unit="g/cm3"/>
-    <composite n="9" ref="C"/>
-    <composite n="13" ref="H"/>
-    <composite n="2" ref="O"/>
-    <composite n="1" ref="N"/>
-  </material>
-  <material name="Rohacell31_50D">
-    <D type="density" value="0.016" unit="g/cm3"/>
-    <composite n="9" ref="C"/>
-    <composite n="13" ref="H"/>
-    <composite n="2" ref="O"/>
-    <composite n="1" ref="N"/>
-  </material>
-  <material name="RPCGasDefault" state="gas">
-    <D type="density" value="0.0037" unit="g/cm3"/>
-    <composite n="209" ref="C"/>
-    <composite n="239" ref="H"/>
-    <composite n="381" ref="F"/>
-  </material>
-  <material name="PolystyreneFoam">
-    <D type="density" value="0.0056" unit="g/cm3"/>
-    <fraction n="1.0" ref="Polystyrene"/>
-  </material>
-  <material name="Kapton">
-    <D value="1.43" unit="g/cm3"/>
-    <composite n="22" ref="C"/>
-    <composite n="10" ref="H"/>
-    <composite n="2" ref="N"/>
-    <composite n="5" ref="O"/>
-  </material>
-  <material name="PEEK">
-    <D value="1.37" unit="g/cm3"/>
-    <composite n="19" ref="C"/>
-    <composite n="12" ref="H"/>
-    <composite n="3" ref="O"/>
-  </material>
-  <material name="TungstenDens23">
-    <D value="17.7" unit="g / cm3"/>
-    <fraction n="0.925" ref="W"/>
-    <fraction n="0.066" ref="Ni"/>
-    <fraction n="0.009" ref="Fe"/>
-  </material>
-  <material name="TungstenDens24">
-    <D value="17.8" unit="g / cm3"/>
-    <fraction n="0.93" ref="W"/>
-    <fraction n="0.061" ref="Ni"/>
-    <fraction n="0.009" ref="Fe"/>
-  </material>
-  <material name="TungstenDens25">
-    <D value="18.2" unit="g / cm3"/>
-    <fraction n="0.950" ref="W"/>
-    <fraction n="0.044" ref="Ni"/>
-    <fraction n="0.006" ref="Fe"/>
-  </material>
-  <material name="CarbonFiber_25percent">
-    <D type="density" value="0.375" unit="g / cm3"/>
-    <fraction n="1.0" ref="CarbonFiber"/>
-  </material>
-  <material name="CarbonFiber_15percent">
-    <D type="density" value="0.225" unit="g / cm3"/>
-    <fraction n="1.0" ref="CarbonFiber"/>
-  </material>
-  <material name="Rohacell31_50percent">
-    <D type="density" value="0.016" unit="g / cm3"/>
-    <fraction n="1.0" ref="Rohacell31"/>
-  </material>
-  <material name="Rohacell31_15percent">
-    <D type="density" value="0.0048" unit="g / cm3"/>
-    <fraction n="1.0" ref="Rohacell31"/>
-  </material>
-  <material name="BoratedPolyethylene5">
-    <D value="0.93" unit="g / cm3"/>
-    <fraction n="0.612" ref="C"/>
-    <fraction n="0.222" ref="O"/>
-    <fraction n="0.116" ref="H"/>
-    <fraction n="0.050" ref="B"/>
-  </material>
-  <material name="SiliconCarbide">
-    <D value="3.1" unit="g / cm3"/>
-    <composite n="1" ref="Si"/>
-    <composite n="1" ref="C"/>
-  </material>
-  <material name="SiliconCarbide_6percent">
-    <D value="0.186" unit="g / cm3"/>
-    <fraction n="1.0" ref="SiliconCarbide"/>
-  </material>
-  <material name="PlasticScint">
-    <D type="density" unit="g/cm3" value="1.032"/>
-    <composite n="9" ref="C"/>
-    <composite n="10" ref="H"/>
-  </material>
-  <material name="PbWO4">
-    <D type="density" value="8.3" unit="g / cm3"/>
-    <composite n="1" ref="Pb"/>
-    <composite n="1" ref="W"/>
-    <composite n="4" ref="O"/>
-  </material>
-
-</materials>
diff --git a/benchmarks/zdc/run_simulation_zdc.sh b/benchmarks/zdc/run_simulation_zdc.sh
deleted file mode 100755
index 6e71f647cd03eba17ccb1d87f137e3e5d92d321a..0000000000000000000000000000000000000000
--- a/benchmarks/zdc/run_simulation_zdc.sh
+++ /dev/null
@@ -1,6 +0,0 @@
-#!/bin/bash
-
-ddsim --runType batch --numberOfEvents 10 \
-      --compactFile benchmarks/zdc/ZDC_example.xml \
-      --inputFiles  ./data/zdc_photons.hepmc \
-      --outputFile  ./sim_output/output_zdc_photons.edm4hep.root
diff --git a/benchmarks/zdc/run_zdc_neutrons.sh b/benchmarks/zdc/run_zdc_neutrons.sh
deleted file mode 100755
index feb0ef3aec52b84dac8909e4da6d5e382a122f64..0000000000000000000000000000000000000000
--- a/benchmarks/zdc/run_zdc_neutrons.sh
+++ /dev/null
@@ -1,50 +0,0 @@
-#!/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  "${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="zdc_uniform_neutrons"
-export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
-
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.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 "benchmarks/zdc/scripts/zdc_neutrons.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
-
-# Run geant4 simulations
-ddsim --runType batch \
-      -v WARNING \
-      --part.minimalKineticEnergy 0.5*GeV  \
-      --filter.tracker edep0 \
-      --numberOfEvents ${JUGGLER_N_EVENTS} \
-      --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR_CONFIG}.xml \
-      --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
-      --outputFile sim_output/${JUGGLER_SIM_FILE}
-
-if [[ "$?" -ne "0" ]] ; then
-  echo "ERROR running npdet"
-  exit 1
-fi
-
diff --git a/benchmarks/zdc/run_zdc_particles.sh b/benchmarks/zdc/run_zdc_particles.sh
new file mode 100755
index 0000000000000000000000000000000000000000..da5dff75dd4224575a7cee25ab1346937eda3dfe
--- /dev/null
+++ b/benchmarks/zdc/run_zdc_particles.sh
@@ -0,0 +1,92 @@
+#!/bin/bash
+
+PARTICLE="neutron"
+
+function print_the_help {
+  echo "USAGE: ${0} [--particle=$PARTICLE] "
+  echo "  OPTIONS: "
+  echo "            --particle   particle to generate and simulate"
+  exit 
+}
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    --particle)
+      PARTICLE="$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
+
+if [[ ! -n  "${JUGGLER_DETECTOR}" ]] ; then 
+  export JUGGLER_DETECTOR="topside"
+fi
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
+  export JUGGLER_N_EVENTS=1000
+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="zdc_${PARTICLE}"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_output/sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+export RESULTS_PATH="results/far_forward/zdc/"
+mkdir -p "${RESULTS_PATH}"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+# Generate the input events
+root -b -q "benchmarks/zdc/scripts/gen_zdc_particles.cxx+(${JUGGLER_N_EVENTS}, \"${PARTICLE}\", ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script: generating input events"
+  exit 1
+fi
+
+# Run geant4 simulations
+ddsim --runType batch \
+      -v WARNING \
+      --part.minimalKineticEnergy 0.5*GeV  \
+      --filter.tracker edep0 \
+      --numberOfEvents ${JUGGLER_N_EVENTS} \
+      --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR_CONFIG}.xml \
+      --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
+      --outputFile ${JUGGLER_SIM_FILE}
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running npdet"
+  exit 1
+fi
+
+echo "Running analysis root scripts"
+root -b -q "benchmarks/zdc/scripts/analysis_zdc_particles.cxx+(\"${JUGGLER_SIM_FILE}\", \"${RESULTS_PATH}\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running root analysis script"
+  exit 1
+fi
+
diff --git a/benchmarks/zdc/scripts/analysis_zdc_particles.cxx b/benchmarks/zdc/scripts/analysis_zdc_particles.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..1d5828e85cb843ef46367e545dfadd4ee67d673e
--- /dev/null
+++ b/benchmarks/zdc/scripts/analysis_zdc_particles.cxx
@@ -0,0 +1,190 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.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 analysis_zdc_particles(
+  const std::string& input_fname = "sim_output/sim_zdc_uniform_neutrons.edm4hep.root",
+  const std::string& results_path = "results/far_forward/zdc/"
+) {
+  // 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<edm4hep::MCParticleData> const& input) {
+    auto p = input[2];
+    auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
+    return energy;
+  };
+
+  // Number of hits
+  auto nhits = [] (const std::vector<edm4hep::SimCalorimeterHitData>& evt) {return (int) evt.size(); };
+
+  // Energy deposition [GeV]
+  auto Esim = [](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
+    auto total_edep = 0.0;
+    for (const auto& i: evt)
+      total_edep += i.energy;
+    return total_edep;
+  };
+
+  // Sampling fraction = Esampling / Ethrown
+  auto fsam = [](const double sampled, const double thrown) {
+    return sampled / thrown;
+  };
+
+  // Define variables
+  auto d1 = d0
+    .Define("Ethr", Ethr, {"MCParticles"})
+    .Define("nhits_Ecal", nhits, {"ZDCEcalHits"})
+    .Define("Esim_Ecal", Esim, {"ZDCEcalHits"})
+    .Define("fsam_Ecal", fsam, {"Esim_Ecal", "Ethr"})
+    .Define("nhits_Hcal", nhits, {"ZDCHcalHits"})
+    .Define("Esim_Hcal", Esim, {"ZDCHcalHits"})
+    .Define("fsam_Hcal", fsam, {"Esim_Hcal", "Ethr"})
+  ;
+
+  // Define Histograms
+  auto hEthr = d1.Histo1D(
+      {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5},
+      "Ethr");
+  auto hNhits_Ecal =
+      d1.Histo1D({"hNhits_Ecal", "Number of Ecal hits per events; Number of Ecal hits; Events",
+                  100, 0.0, 2000.0},
+                 "nhits_Ecal");
+  auto hEsim_Ecal = d1.Histo1D(
+      {"hEsim_Ecal", "Ecal Energy Deposit; Ecal Energy Deposit [GeV]; Events", 100, 0.0, 1.0},
+      "Esim_Ecal");
+  auto hfsam_Ecal = d1.Histo1D(
+      {"hfsam_Ecal", "Ecal Sampling Fraction; Ecal Sampling Fraction; Events", 100, 0.0, 0.1},
+      "fsam_Ecal");
+  auto hNhits_Hcal =
+      d1.Histo1D({"hNhits_Hcal", "Number of Hcal hits per events; Number of Hcal hits; Events",
+                  100, 0.0, 2000.0},
+                 "nhits_Hcal");
+  auto hEsim_Hcal = d1.Histo1D(
+      {"hEsim_Hcal", "Hcal Energy Deposit; Hcal Energy Deposit [GeV]; Events", 100, 0.0, 1.0},
+      "Esim_Hcal");
+  auto hfsam_Hcal = d1.Histo1D(
+      {"hfsam_Hcal", "Hcal Sampling Fraction; Hcal Sampling Fraction; Events", 100, 0.0, 0.1},
+      "fsam_Hcal");
+
+  // 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);
+    auto h = hEthr->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    c1->SaveAs(TString(results_path + "/zdc_Ethr.png"));
+    c1->SaveAs(TString(results_path + "/zdc_Ethr.pdf"));
+  }
+
+  // Ecal
+  {
+    TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
+    c2->SetLogy(1);
+    auto h = hNhits_Ecal->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    c2->SaveAs(TString(results_path + "/zdc_nhits_Ecal.png"));
+    c2->SaveAs(TString(results_path + "/zdc_nhits_Ecal.pdf"));
+  }
+
+  {
+    TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
+    c3->SetLogy(1);
+    auto h = hEsim_Ecal->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    c3->SaveAs(TString(results_path + "/zdc_Esim_Ecal.png"));
+    c3->SaveAs(TString(results_path + "/zdc_Esim_Ecal.pdf"));
+  }
+
+  std::cout << "derp4\n";
+  {
+    TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
+    c4->SetLogy(1);
+    auto h = hfsam_Ecal->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    //h->Fit("gaus", "", "", 0.01, 0.1);
+    //h->GetFunction("gaus")->SetLineWidth(2);
+    //h->GetFunction("gaus")->SetLineColor(kRed);
+    c4->SaveAs(TString(results_path + "/zdc_fsam_Ecal.png"));
+    c4->SaveAs(TString(results_path + "/zdc_fsam_Ecal.pdf"));
+  }
+
+  // Hcal
+  {
+    TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
+    c2->SetLogy(1);
+    auto h = hNhits_Hcal->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    c2->SaveAs(TString(results_path + "/zdc_nhits_Hcal.png"));
+    c2->SaveAs(TString(results_path + "/zdc_nhits_Hcal.pdf"));
+  }
+
+  {
+    TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
+    c3->SetLogy(1);
+    auto h = hEsim_Hcal->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    c3->SaveAs(TString(results_path + "/zdc_Esim_Hcal.png"));
+    c3->SaveAs(TString(results_path + "/zdc_Esim_Hcal.pdf"));
+  }
+
+  std::cout << "derp4\n";
+  {
+    TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
+    c4->SetLogy(1);
+    auto h = hfsam_Hcal->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    //h->Fit("gaus", "", "", 0.01, 0.1);
+    //h->GetFunction("gaus")->SetLineWidth(2);
+    //h->GetFunction("gaus")->SetLineColor(kRed);
+    c4->SaveAs(TString(results_path + "/zdc_fsam_Hcal.png"));
+    c4->SaveAs(TString(results_path + "/zdc_fsam_Hcal.pdf"));
+  }
+}
diff --git a/benchmarks/zdc/scripts/zdc_neutrons.cxx b/benchmarks/zdc/scripts/gen_zdc_particles.cxx
similarity index 63%
rename from benchmarks/zdc/scripts/zdc_neutrons.cxx
rename to benchmarks/zdc/scripts/gen_zdc_particles.cxx
index 4c1f90032e964d9201412e336613f3c0f8c6bbc6..21a5a2a4453496298bfb996df370e62cd99f9784 100644
--- a/benchmarks/zdc/scripts/zdc_neutrons.cxx
+++ b/benchmarks/zdc/scripts/gen_zdc_particles.cxx
@@ -16,9 +16,12 @@
 #include "TMath.h"
 #include "TRandom.h"
 
+#include <Math/Vector3D.h>
+#include <Math/RotationY.h>
+
 using namespace HepMC3;
 
-void zdc_neutrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/zdc_neutrons.hepmc") {
+void gen_zdc_particles(int n_events = 1e6, const std::string& particle = "neutron", double p_start = 0.0, double p_end = 30.0, const std::string& out_fname = "./data/zdc_neutrons.hepmc") {
   WriterAscii hepmc_output(out_fname);
   int events_parsed = 0;
   GenEvent evt(Units::GEV, Units::MM);
@@ -26,6 +29,9 @@ void zdc_neutrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0,
   // Random number generator
   TRandom* r1 = new TRandom();
 
+  // Set crossing angle [rad]
+  double crossing_angle = -0.025;
+
   // 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
@@ -44,21 +50,35 @@ void zdc_neutrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0,
     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 p_mag    = r1->Uniform(p_start, p_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);
+    Double_t p0_x     = p_mag * std::cos(phi) * std::sin(theta);
+    Double_t p0_y     = p_mag * std::sin(phi) * std::sin(theta);
+    Double_t p0_z     = p_mag * 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);
+    // Rotate the vector in Y by crossing angle when particles are being generated
+    ROOT::Math::XYZVector p0{p0_x, p0_y, p0_z};
+    ROOT::Math::RotationY r(crossing_angle);
+    auto p = r * p0;
+    auto px = p.X();
+    auto py = p.Y();
+    auto pz = p.Z();
 
     // 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.939 * 0.939))), 2112, 1);
+    // pdgid 22 - photon massless
+    // pdgid 2112 - neutron 939.565 MeV/c^2
+    GenParticlePtr p3;
+    if (particle == "neutron") {
+      p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, std::hypot(p_mag, 0.939565)), 2112, 1);
+    } else if (particle == "photon") {
+      p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, p_mag), 22, 1);
+    } else {
+      std::cout << "Particle type " << particle << " not recognized!" << std::endl;
+      exit(-1);
+    }
 
     GenVertexPtr v1 = std::make_shared<GenVertex>();
     v1->add_particle_in(p1);
diff --git a/benchmarks/zdc/scripts/zdc_neutrons_analysis.cxx b/benchmarks/zdc/scripts/zdc_neutrons_analysis.cxx
deleted file mode 100644
index 742a19744713872f8b4a6e183e32a6322455b1e9..0000000000000000000000000000000000000000
--- a/benchmarks/zdc/scripts/zdc_neutrons_analysis.cxx
+++ /dev/null
@@ -1,135 +0,0 @@
-////////////////////////////////////////
-// Read reconstruction ROOT output file
-// Plot variables
-////////////////////////////////////////
-
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "edm4hep/MCParticleCollection.h"
-#include "edm4hep/SimCalorimeterHitCollection.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 zdc_neutrons_analysis(const char* input_fname = "sim_output/sim_zdc_uniform_neutrons.edm4hep.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<edm4hep::MCParticleData> const& input) {
-    auto p = input[2];
-    auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
-    return energy;
-  };
-
-  // Number of hits
-  auto nhits = [] (const std::vector<edm4hep::SimCalorimeterHitData>& evt) {return (int) evt.size(); };
-
-  // Energy deposition [GeV]
-  auto Esim = [](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
-    auto total_edep = 0.0;
-    for (const auto& i: evt)
-      total_edep += i.energy;
-    return total_edep;
-  };
-
-  // Sampling fraction = Esampling / Ethrown
-  auto fsam = [](const double sampled, const double thrown) {
-    return sampled / thrown;
-  };
-
-  // Define variables
-  auto d1 = d0.Define("Ethr", Ethr, {"MCParticles"})
-                .Define("nhits", nhits, {"EcalBarrelHits"})
-                .Define("Esim", Esim, {"EcalBarrelHits"})
-                .Define("fsam", fsam, {"Esim", "Ethr"});
-
-  // Define Histograms
-  auto hEthr = d1.Histo1D(
-      {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5},
-      "Ethr");
-  auto hNhits =
-      d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events",
-                  100, 0.0, 2000.0},
-                 "nhits");
-  auto hEsim = d1.Histo1D(
-      {"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 100, 0.0, 1.0},
-      "Esim");
-  auto hfsam = d1.Histo1D(
-      {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1},
-      "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);
-    auto h = hEthr->DrawCopy();
-    //h->GetYaxis()->SetTitleOffset(1.4);
-    h->SetLineWidth(2);
-    h->SetLineColor(kBlue);
-    c1->SaveAs("results/zdc_neutrons_Ethr.png");
-    c1->SaveAs("results/zdc_neutrons_Ethr.pdf");
-  }
-  std::cout << "derp1\n";
-
-  {
-    TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
-    c2->SetLogy(1);
-    auto h = hNhits->DrawCopy();
-    //h->GetYaxis()->SetTitleOffset(1.4);
-    h->SetLineWidth(2);
-    h->SetLineColor(kBlue);
-    c2->SaveAs("results/zdc_neutrons_nhits.png");
-    c2->SaveAs("results/zdc_neutrons_nhits.pdf");
-  }
-
-  {
-    TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
-    c3->SetLogy(1);
-    auto h = hEsim->DrawCopy();
-    //h->GetYaxis()->SetTitleOffset(1.4);
-    h->SetLineWidth(2);
-    h->SetLineColor(kBlue);
-    c3->SaveAs("results/zdc_neutrons_Esim.png");
-    c3->SaveAs("results/zdc_neutrons_Esim.pdf");
-  }
-
-  std::cout << "derp4\n";
-  {
-    TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
-    c4->SetLogy(1);
-    auto h = hfsam->DrawCopy();
-    //h->GetYaxis()->SetTitleOffset(1.4);
-    h->SetLineWidth(2);
-    h->SetLineColor(kBlue);
-    //h->Fit("gaus", "", "", 0.01, 0.1);
-    //h->GetFunction("gaus")->SetLineWidth(2);
-    //h->GetFunction("gaus")->SetLineColor(kRed);
-    c4->SaveAs("results/zdc_neutrons_fsam.png");
-    c4->SaveAs("results/zdc_neutrons_fsam.pdf");
-  }
-}