Resolve "Material scan benchmark"

benchmarks/tracking_detectors/scripts/matscan_plot.py

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

benchmarks/tracking_detectors/scripts/test_matscan.cxx

+// 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);
+}
+