diff --git a/benchmarks/tracking_detectors/scripts/matscan_plot.py b/benchmarks/tracking_detectors/scripts/matscan_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..7cd5353dbf79aa1e9cf3c2159afb5c3980ffe210
--- /dev/null
+++ b/benchmarks/tracking_detectors/scripts/matscan_plot.py
@@ -0,0 +1,218 @@
+## plot test_matscan.cxx results
+## Shujie Li, 08, 2021
+## usage: python matscan_plot.py zrange rrange
+import pandas as pd
+import matplotlib.pyplot as plt
+import numpy as np
+import math as m
+import sys
+from matplotlib import colors
+from scipy import stats
+from matplotlib.colors import LogNorm
+
+
+plt.rcParams['figure.figsize'] = [10.0, 6.0]
+# plt.rcParams['fure.dpi'] = 80
+
+SMALL_SIZE = 14
+MEDIUM_SIZE = 16
+BIGGER_SIZE = 18
+
+plt.rc('font', size=SMALL_SIZE) # controls default text sizes
+plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
+plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
+plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
+plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
+plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
+plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title]
+
+
+def cart2sph(x,y,z):
+ x=np.array(x)
+ y=np.array(y)
+ z=np.array(z)
+ XsqPlusYsq = x**2 + y**2
+ r = np.sqrt(XsqPlusYsq + z**2) # r
+ elev = np.arctan2(np.sqrt(XsqPlusYsq),z) # theta (polar to z axis)
+ az = np.arctan2(y,x) # phi (azimuthal)
+ return r, elev, az
+
+## read a single scan table for cylinder volume
+def read_table_cylind(zmax=140,rmax=60,indir='./'):
+ header = np.array(['x', 'y', 'z', 'ind', 'material', 'Z', 'density', 'rad_length', 'thickness', 'path_length', 'sum_x0',"end_x","end_y","end_z"])
+ df = pd.read_csv(indir+"all_z%g_r%g.dat" %(zmax,rmax),names=header,delim_whitespace=True)
+ df["x0"]=np.array(df["thickness"])/np.array(df["rad_length"])
+ df=df.reset_index(drop=True).drop_duplicates()
+
+ # calcualte polar angle then eta
+ ax = df["x"]
+ ay = df["y"]
+ az = df["z"]
+ r,polar,azi=cart2sph(ax,ay,az)
+ eta = -np.log(np.tan(polar/2.0))
+ df["eta"] = eta
+ df["azi_angle"] = azi
+ df["pol_angle"] = polar
+
+ print(df.head())
+ print(df.material.unique())
+
+# df=df[np.isfinite(df).all(1)] #drop in
+ return df
+
+def add_kin(df):
+ # calcualte polar angle then eta
+ ax = df["x"]
+ ay = df["y"]
+ az = df["z"]
+ r,polar,azi=cart2sph(ax,ay,az)
+ eta = -np.log(np.tan(polar/2.0))
+ df["eta"] = eta
+ df["azi_angle"] = azi
+ df["pol_angle"] = polar
+ return df
+
+## group by xyz coordinate to get sum x0 per track with selected materials
+def get_x0(df0,mat_name=""):
+ if len(mat_name)>0:
+ df=df0[df0["material"]==mat_name]
+ else:
+ df=df0
+ # df=df[df["material"]!="CarbonFiber_25percent"]
+ dfg=df.groupby(["x","y","z"])["x0"].sum()
+ dfg=dfg.reset_index()
+ dfg=dfg[np.isfinite(dfg).all(1)] #drop inf
+ dfg=add_kin(dfg)
+ return dfg
+
+
+## plot mean,min,max X0 of each eta bin
+def plot_1d(dfgz,mat_name="",ax="eta"):
+ xlow = 0
+ values=np.array(dfgz["x0"])*100
+ xx = np.array(dfgz[ax])
+ if ax=="eta":
+ xhi=3.5
+ xlow = -xhi
+ elif ax=="pol_angle":
+ xhi=90
+ xx=180/2-xx*180/np.pi
+ elif ax=="z": # z projected at barrel radius
+ xhi=140
+ xx*=(43.23/60)
+ ## with bin stats
+ nbin = 50
+ x0mean,binedge,binnumber=stats.binned_statistic(xx,values, 'mean', bins=nbin,range=[xlow,xhi])
+ ## !!! possible bugs in this error bar calculation
+ x0max ,binedge,binnumber=stats.binned_statistic(xx,values, 'max', bins=nbin,range=[xlow,xhi])
+ x0min ,binedge,binnumber=stats.binned_statistic(xx,values, 'min', bins=nbin,range=[xlow,xhi])
+
+ bin_center = (binedge[0:-1]+binedge[1:])/2
+ plt.plot(bin_center,x0mean,label=mat_name)
+ plt.fill_between(bin_center,x0min,x0max,alpha=0.2)
+ plt.xlim(xlow,xhi)
+ plt.grid()
+ plt.suptitle("total X0")
+ plt.xlabel(ax)
+
+ return bin_center, x0mean, x0max, x0min
+
+
+
+if __name__ == "__main__":
+
+ ## corresponding to the scan output filenmae from test_matscan.cxx, in cm
+ zrange = int(sys.argv[1])
+ rrange = int(sys.argv[2])
+ outdir = './'
+ indir = './'
+ cols = np.array(["x","y","z","material","thickness","path_length"])
+ df = read_table_cylind(zrange,rrange,indir)
+
+
+ ## -----------------plot side view of geometry, z v.s. R-------------
+ plt.figure()
+ xe = df["end_x"]
+ ye = df["end_y"]
+ ze = df["end_z"]
+ re = np.sqrt(xe**2+ye**2)
+ # c5=(df["end_x"]**2+df["end_y"]**2)<11**2
+ # c6=df["material"]=="Aluminum"
+ plt.scatter(ze,re,s=0.001,marker=".")
+
+ # plt.xlabel("$\eta$")
+ plt.xlabel("z [cm]")
+ plt.ylabel("R [cm]")
+ # plt.xlim(0,4)
+ plt.grid()
+ plt.savefig(outdir+"/matscan_geo_z_z%g_r%g.png" %(zrange,rrange))
+
+ ## -----------------plot side view of geometry, eta v.s. R-------------
+ plt.figure()
+ xe = df["end_x"]
+ ye = df["end_y"]
+ ze = df["end_z"]
+ re = np.sqrt(xe**2+ye**2)
+ plt.scatter(df["eta"],re,s=0.001,marker=".")
+
+ plt.xlabel("$\eta$")
+ plt.ylabel("R [cm]")
+ # plt.xlim(0,4)
+ plt.grid()
+ plt.savefig(outdir+"/matscan_geo_eta_z%g_r%g.png" %(zrange,rrange))
+
+
+ ## -----------------plot 1 d matscan z v.s. X0-------------
+ df_all=get_x0(df)
+ plt.figure()
+ ax="eta"
+ print(df.columns,df_all.columns)
+ print(df_all['eta'])
+ _=plot_1d(df_all,"Total",ax)
+
+ mat_list=np.delete(df.material.unique(),0)
+ # mat_list.drop("Vacuum")
+ # mat_list=['NONE','Air','Beryllium','Aluminum','Silicon','CarbonFiber']
+ # mat_list=["Aluminum","Silicon","CarbonFiber"]
+ for mat_name in mat_list:
+ df_mat = get_x0(df,mat_name)
+ _=plot_1d(df_mat,mat_name,ax)
+
+ plt.legend(loc="upper right")
+ plt.xlabel("$\eta$")
+ plt.ylabel("X/X0 (%)")
+ plt.savefig(outdir+"/matscan_1d_z%g_r%g.png" %(zrange,rrange))
+
+ ## -----------------plot 2 d matscan z v.s. phi-------------
+
+ dfgz=df_all
+ values=dfgz["x0"]
+ xx = dfgz["eta"]
+ yy = dfgz["azi_angle"] # in rad
+ ## with bin stats
+ nbin = 50
+ x0mean_2d,x_edge,y_edge,binnumber_2d=stats.binned_statistic_2d(xx,yy, values, 'mean', bins=[nbin,nbin],range=[[-5,5],[-5,5]])
+ # x0max_2d ,_,_,_=stats.binned_statistic_2d(xx,yy,values, 'max', bins=[nbin,nbin],range=[[-5,5],[-5,5]])
+ # x0min_2d ,_,_,_=stats.binned_statistic_2d(xx,yy,values, 'min', bins=[nbin,nbin],range=[[-5,5],[-5,5]])
+
+ x_center = (x_edge[0:-1]+x_edge[1:])/2
+ y_center = (y_edge[0:-1]+y_edge[1:])/2
+ # plt.contour(x0mean_2d.transpose(),extent=[x_edge[0],x_edge[-1],y_edge[0],y_edge[-1]],
+ # linewidths=3, cmap = plt.cm.rainbow)
+ # plt.colorbar()
+
+
+ fig=plt.figure(figsize=(7, 6))
+ c = plt.pcolormesh(x_edge,y_edge/np.pi*180,x0mean_2d.transpose(),norm=LogNorm(vmin=0.001, vmax=10))
+ fig.colorbar(c)
+ # plt.colorbar()
+ # plt.ylim(-np.pi,np.pi)
+ plt.ylim(-180,180)
+ plt.xlim(-5,5)
+ plt.xlabel("$\eta$")
+ plt.ylabel("Azimuthal angle [degree]")
+ plt.suptitle("total X0 [%]")
+
+ plt.savefig(outdir+"/matscan_2d_z%g_r%g.png" %(zrange,rrange))
+
+
diff --git a/benchmarks/tracking_detectors/scripts/test_matscan.cxx b/benchmarks/tracking_detectors/scripts/test_matscan.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d178e50a1fbeccb20bbdf134efef52d95a5ad197
--- /dev/null
+++ b/benchmarks/tracking_detectors/scripts/test_matscan.cxx
@@ -0,0 +1,139 @@
+// material scan of a cylinder volume for athena
+// the scan starts at 0 for now
+// based on https://dd4hep.web.cern.ch/dd4hep/reference/MaterialScan_8cpp_source.html
+// dd4hep also provides command line utilities materialBudget and materialScan
+//
+// Shujie Li, Aug 2021
+
+R__LOAD_LIBRARY(libDDCore.so)
+R__LOAD_LIBRARY(libDDG4.so)
+R__LOAD_LIBRARY(libDDG4IO.so)
+#include "DD4hep/Detector.h"
+#include "DD4hep/DetElement.h"
+#include "DD4hep/Objects.h"
+#include "DD4hep/Detector.h"
+#include "DDG4/Geant4Data.h"
+#include "DDRec/CellIDPositionConverter.h"
+#include "DDRec/SurfaceManager.h"
+#include "DDRec/Surface.h"
+
+#include "TCanvas.h"
+#include "TChain.h"
+#include "TGeoMedium.h"
+#include "TGeoManager.h"
+#include "DDRec/MaterialScan.h"
+#include "DDRec/MaterialManager.h"
+#include "DD4hep/Detector.h"
+#include "DD4hep/Printout.h"
+#include "fmt/core.h"
+
+#include <iostream>
+#include <fstream>
+
+using namespace dd4hep;
+using namespace dd4hep::rec;
+
+ // for silicon tracker, the outer barrel rmax ~45, outer endcap zmax~121
+void test_matscan(const char* compact = "athena.xml",double zmax=130, double rmax = 61){
+
+ dd4hep::Detector& detector = dd4hep::Detector::getInstance();
+ detector.fromCompact(compact);
+ MaterialScan matscan(detector);
+ TString detname = "all";
+ fmt::print("\n");
+ fmt::print("All detector subsystem names:\n");
+ for(const auto& d : detector.detectors() ) {
+ fmt::print(" {}\n", d.first);
+ }
+ // to do: material scan of given dtectors.
+
+ // TString det_list[14]={"TrackerBarrel_Inner","TrackerBarrel_Outer","TrackerEndcapN_Inner","TrackerEndcapN_Outer","TrackerEndcapP_Inner","TrackerEndcapP_Outer","TrackerSubAssembly_Inner","TrackerSubAssembly_Outer","VertexBarrel","VertexBarrelSubAssembly","VertexEndcapN","VertexEndcapP","VertexEndcapSubAssembly","cb_DIRC"};
+ // for (int dd=0;dd<14;dd++){
+ // TString detname = det_list[dd];
+ // matscan.setDetector(detname.Data());
+
+ const char* fmt1 = "%8.3f %8.3f %8.3f %3d %-20s %3.0f %8.4f %11.4f %11.4f %11.4f %11.4f %8.3f %8.3f %8.3f\n";
+
+ // the beam vacuum is missing if starting from origin and ending at negative z. need to check later
+ double x0=0,y0=0,z0=0.001,x1,y1,z1; // cm
+ double epsilon=1e-4; // (mm) default 1e-4: Materials with a thickness smaller than epsilon (default 1e-4=1mu
+
+ const double DegToRad=3.141592653589793/180.0;
+ double phi,dphi=0.5; // Degree
+ double dz=0.5, dr=0.5; // cm
+
+ vector<double> lx,ly,lz;
+ // prepare coord. for barrel
+ for(z1=-zmax;z1<=zmax;z1+=dz){
+ for(phi=0.;phi<=360;phi+=dphi){
+ x1 = cos(phi*DegToRad)*rmax;
+ y1 = sin(phi*DegToRad)*rmax;
+ lx.push_back(x1);
+ ly.push_back(y1);
+ lz.push_back(z1);
+ }
+ }
+ // prepare coord. for endcaps
+ for (double r=1; r<=rmax; r+=dr){
+ for(phi=0.;phi<=360;phi+=dphi){
+ x1 = cos(phi*DegToRad)*r;
+ y1 = sin(phi*DegToRad)*r;
+ lx.push_back(x1);
+ ly.push_back(y1);
+ lz.push_back(zmax);
+ lx.push_back(x1);
+ ly.push_back(y1);
+ lz.push_back(-zmax);
+ }
+ }
+
+ // loop over coord ponits for material scan
+ long nn = lx.size();
+ TString fname = Form("%s_z%g_r%g.dat",detname.Data(),zmax,rmax);
+ FILE * pFile;
+ pFile = fopen (fname,"w");
+
+ for(int ii=0;ii<nn;ii++){
+ x1=lx[ii]; y1=ly[ii]; z1=lz[ii];
+ if(ii%5000==0)
+ cout<<ii<<"/"<<nn<<":"<<x1<<" "<<y1<<" "<<z1<<endl;
+
+ // if ((x1*x1+y1*y1)>(60*60))
+ // continue;
+ Vector3D p0(x0, y0, z0), p1(x1, y1, z1);
+ Vector3D end, direction;
+ direction = (p1-p0).unit();
+ const auto& placements = matscan.scan(x0,y0,z0,x1,y1,z1,epsilon);
+ // matscan.print(x0,y0,z0,x1,y1,z1,epsilon);
+
+ double sum_x0 = 0;
+ double sum_lambda = 0;
+ double path_length = 0, total_length = 0;
+
+
+ for (unsigned i=0;i<placements.size();i++){
+
+ TGeoMaterial* mat=placements[i].first->GetMaterial();
+ double length = placements[i].second;
+ double nx0 = length / mat->GetRadLen();
+ double nLambda = length / mat->GetIntLen();
+ sum_x0 += nx0;
+ sum_lambda += nLambda;
+ path_length += length;
+ total_length += length;
+ end = p0 + total_length * direction;
+
+
+ fprintf(pFile, fmt1,x1, y1, z1, i+1, mat->GetName(), mat->GetZ(),
+ mat->GetDensity(), mat->GetRadLen(),
+ length, path_length, sum_x0,end[0], end[1], end[2]);
+
+
+
+ }
+ // cout<<detname<<" "<<x1<<","<<y1<<","<<z1<<endl;
+ // cout<<x1<<","<<y1<<","<<z1<<": "<<placements.size()<<" "<<sum_x0<<" "<<total_length<<endl;
+ }
+ fclose (pFile);
+}
+