Merge branch 'material' into 'master'

Add tutorial for Material Scan

See merge request cepc/CEPCSW!175

Analysis/MaterialScan/README.md

+# Howto scan the material
+
+## Introduction
+Use DD4hep tool [materialScan](https://dd4hep.web.cern.ch/dd4hep/reference/materialScan_8cpp_source.html) to get the material passing through for a geatino partical (i.e. a straight line) at a given angle.
+
+## Material Scan
+
+Select the dd4hep geometry files you wish to scan, and access the interactive interface of materialScan.
+
+```bash
+# Access to the interactive interface of materialScan
+materialScan ./Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01.xml -interactive
+
+# if you want to scan the material for a specific detector, disable/enable the detectors in the xml file
+materialScan ./Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyBeamPipe.xml -interactive
+```
+
+copy the following code to the interactive interface, and press enter.
+*bins* is the number of bins for theta and phi, *world_size* is the size of the scanning area (mm).
+
+```bash
+FILE* outFile = freopen("MaterialScanLogs.txt", "w", stdout);
+int bins = 100;
+double world_size = 10000;
+for(int thetabin=0; thetabin<=bins; thetabin++){
+    double theta = thetabin * M_PI / bins;
+    double z = world_size*cos(theta);
+    double tranverse = world_size*sin(theta);
+    for(int phibin=-bins; phibin<=bins; phibin++){
+        double phi = phibin * M_PI / bins;
+        double x = tranverse*cos(phi);
+        double y = tranverse*sin(phi);
+        gMaterialScan->print(0,0,0,x,y,z);
+    }
+}
+```
+
+The "FILE* outFile" truncates the "stdout" output stream to the file "MaterialScanLogs.txt", so you can't see any output. Be patient and wait for the file to be created, it may take 1-2 minutes depending on the detectors & the number of bins.
+
+The file "MaterialScanLogs.txt" contains all the material passing through for a straight line at the given angles with the size of bins of theta and phi.
+
+## Draw the material budget for single xml file
+
+We give an example script [options/extract_x0.py](options/extract_x0.py) to draw the material budget from the scan file generated above (e.x. MaterialScanLogs.txt).
+
+```bash
+# Draw the material budget for single xml file
+# --input_file: the scan file generated above
+# --output_file: the output file name (default: material_budget.png)
+# --detector: the detector name to be used as the title of the plot (default: all_world)
+# --bins: the bins of theta and phi (default: 100), the theta-bins is set to "bins"+1, the phi-bins is set to "bins"*2+1
+# --bias: the bias used to cut the theta range to [bias, bins-bias] (default: 0)
+python options/extract_x0.py \
+    --input_file MaterialScanLogs.txt \
+    --output_file material_budget.png \
+    --detector all_world \
+    --bins 100 \
+    --bias 0
+```
+
+Make sure the property "bins" is the same as the one used in the [Material Scan](#material-scan).
+
+The output is a png file shows the material budget at different theta and phi, and the average X0 for theta/phi.
+
+## Draw the material budget for multiple xml files
+
+For multiple xml files, we can use the script [options/extract_x0_multi.py](options/extract_x0_multi.py) to draw the material budget.
+
+```bash
+# Draw the material budget for multiple xml files
+# --input_files: the scan files generated above, files names separated with spaces
+# --labels: the labels name to be used as the legend of the plot
+# --output_file: the output file name (default: accumulated_material_budget.png)
+# --bins: the bins of theta and phi (default: 100), the theta-bins is set to "bins"+1, the phi-bins is set to "bins"*2+1
+# --bias: the bias used to cut the theta range to [bias, bins-bias] (default: 0)
+python options/extract_x0_multi.py \
+    --input_file BeamPipe.txt Lumical.txt VXD.txt FTD.txt SIT.txt \
+    --labels BeamPipe Lumical VXD FTD SIT \
+    --output_file accumulated_material_budget.png \
+    --bins 100 \
+    --bias 0
+```
+
+The output is a png file shows the average X0 for theta/phi with different labels.
\ No newline at end of file

Analysis/MaterialScan/options/extract_multi_x0.py

+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import argparse
+
+def extract_x0_values(filename):
+    """Extract x0 values from a single file"""
+    x0_values = []
+    counter = 0
+    with open(filename, 'r') as file:
+        for line in file:
+            if '|     0 Average Material' in line:
+                counter += 1
+                items = [item.strip() for item in line.split('|')]
+                x0_value = float(items[1].split()[10])
+                x0_values.append(x0_value)
+                print(f"Processing file {filename}: ", counter, " x0: ", x0_value)
+    return x0_values
+
+def process_multiple_files(file_list, bins, bias):
+    """Process multiple files and return their x0 matrices"""
+    theta_bins = bins + 1
+    phi_bins = bins*2 + 1
+    range_theta = [bias, theta_bins-bias]
+    
+    matrices = []
+    for file in file_list:
+        x0_values = extract_x0_values(file)
+        x0_matrix = np.array(x0_values).reshape(theta_bins, phi_bins)[range_theta[0]:range_theta[1], :]
+        matrices.append(x0_matrix)
+    
+    return matrices
+
+def plot_accumulated_projections(matrices, output_file, bins, bias, labels):
+    """Plot accumulated projections"""
+    theta_bins = bins + 1
+    phi_bins = bins*2 + 1
+    range_theta = [bias, theta_bins-bias]
+    
+    # Create angle arrays
+    phi = np.linspace(-180, 180, phi_bins)
+    theta = np.linspace(range_theta[0]*180/theta_bins, range_theta[1]*180/theta_bins, range_theta[1]-range_theta[0])
+    
+    # Create figure
+    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
+    
+    # Accumulation data
+    accumulated_theta = np.zeros_like(theta)
+    accumulated_phi = np.zeros_like(phi)
+    
+    colors = plt.cm.viridis(np.linspace(0, 1, len(matrices)))
+    
+    # Plot phi projection
+    for i, matrix in enumerate(matrices):
+        phi_projection = np.sum(matrix, axis=0) / (range_theta[1]-range_theta[0])
+        accumulated_phi += phi_projection
+        ax1.fill_between(phi, accumulated_phi, accumulated_phi - phi_projection, 
+                        label=labels[i], color=colors[i], alpha=0.6)
+    
+    # Plot theta projection
+    for i, matrix in enumerate(matrices):
+        theta_projection = np.sum(matrix, axis=1) / phi_bins
+        accumulated_theta += theta_projection
+        ax2.fill_between(theta, accumulated_theta, accumulated_theta - theta_projection,
+                        label=labels[i], color=colors[i], alpha=0.6)
+    
+    # Set phi projection plot
+    ax1.set_title('Accumulated Projection along Phi', fontsize=14, pad=10)
+    ax1.set_xlabel('Phi (degree)', fontsize=12)
+    ax1.set_ylabel('Accumulated X/X0', fontsize=12)
+    ax1.grid(True, linestyle='--', alpha=0.7)
+    ax1.legend(fontsize=10)
+    ax1.spines['top'].set_visible(False)
+    ax1.spines['right'].set_visible(False)
+    
+    # Set theta projection plot
+    ax2.set_title('Accumulated Projection along Theta', fontsize=14, pad=10)
+    ax2.set_xlabel('Theta (degree)', fontsize=12)
+    ax2.set_ylabel('Accumulated X/X0', fontsize=12)
+    ax2.grid(True, linestyle='--', alpha=0.7)
+    ax2.legend(fontsize=10)
+    ax2.spines['top'].set_visible(False)
+    ax2.spines['right'].set_visible(False)
+    
+    plt.suptitle('Detector Material Budget Accumulation Analysis', fontsize=16, y=0.95)
+    plt.tight_layout()
+    
+    plt.savefig(output_file, dpi=600, bbox_inches='tight')
+    plt.show()
+
+def main():
+    parser = argparse.ArgumentParser(description='Process multiple material scan data and generate accumulated distribution plots.')
+    parser.add_argument('--input_files', nargs='+', required=True,
+                        help='List of input files in order')
+    parser.add_argument('--labels', nargs='+', required=True,
+                        help='Labels for each input file')
+    parser.add_argument('--output_file', default='accumulated_material_budget.png',
+                        help='Output PNG file path')
+    parser.add_argument('--bins', type=int, default=100,
+                        help='theta bins is "bins"+1, phi bins is "bins"*2+1, (default: 100)')
+    parser.add_argument('--bias', type=int, default=0,
+                        help='Bias value for theta range (default: 0)')
+    
+    args = parser.parse_args()
+    
+    if len(args.input_files) != len(args.labels):
+        raise ValueError("Number of input files must match number of labels!")
+    
+    # Process all input files
+    matrices = process_multiple_files(args.input_files, args.bins, args.bias)
+    
+    # Plot accumulated graphs
+    plot_accumulated_projections(matrices, args.output_file, args.bins, args.bias, args.labels)
+
+if __name__ == '__main__':
+    main()
\ No newline at end of file

Analysis/MaterialScan/options/extract_x0.py

+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import argparse
+
+def extract_x0_values(filename):
+    x0_values = []
+    counter = 0
+    with open(filename, 'r') as file:
+        for line in file:
+            if '|     0 Average Material' in line:
+                counter += 1
+                items = [item.strip() for item in line.split('|')]
+                x0_value = float(items[1].split()[10])
+                x0_values.append(x0_value)
+                print("processing: ", counter, " x0: ", x0_value)
+    return x0_values
+
+def plot_all_in_one(x0_values, output_file, bins, bias, detector):
+    """Plot all charts in one figure"""
+    theta_bins = bins + 1
+    phi_bins = bins*2 + 1
+    range_theta = [bias, theta_bins-bias]
+    # Adjust figure size and ratio, increase height to accommodate title
+    fig = plt.figure(figsize=(28, 8))
+    
+    # Adjust spacing between subplots and top margin
+    gs = fig.add_gridspec(1, 3, width_ratios=[1.2, 1, 1], wspace=0.25, top=0.85)
+    
+    # Add three subplots
+    ax1 = fig.add_subplot(gs[0])  # Heat map
+    ax2 = fig.add_subplot(gs[1])  # Phi projection
+    ax3 = fig.add_subplot(gs[2])  # Theta projection
+    
+    # Draw heat map
+    X0_matrix = np.array(x0_values).reshape(theta_bins, phi_bins)[range_theta[0]:range_theta[1], :]
+    phi = np.linspace(-180, 180, phi_bins)
+    theta = np.linspace(range_theta[0]*180/theta_bins, range_theta[1]*180/theta_bins, range_theta[1]-range_theta[0])
+    THETA, PHI = np.meshgrid(theta, phi)
+    im = ax1.pcolormesh(THETA, PHI, X0_matrix.T, shading='auto', cmap='viridis')
+    cbar = plt.colorbar(im, ax=ax1)
+    cbar.set_label('X/X0', fontsize=12)
+    ax1.set_title('Material Budget Distribution', fontsize=12, pad=10)
+    ax1.set_xlabel('Theta (angle)', fontsize=10)
+    ax1.set_ylabel('Phi (angle)', fontsize=10)
+    
+    # Draw phi projection
+    phi_projection = np.sum(X0_matrix, axis=0) / (range_theta[1]-range_theta[0])
+    ax2.plot(phi, phi_projection, 'b-', linewidth=2)
+    ax2.fill_between(phi, phi_projection, alpha=0.3)
+    ax2.set_title('Projection along Phi', fontsize=12, pad=10)
+    ax2.set_xlabel('Phi (degree)', fontsize=10)
+    ax2.set_ylabel('Average X/X0', fontsize=10)
+    ax2.grid(True, linestyle='--', alpha=0.7)
+    ax2.spines['top'].set_visible(False)
+    ax2.spines['right'].set_visible(False)
+    
+    # Draw theta projection
+    theta_projection = np.sum(X0_matrix, axis=1) / phi_bins
+    ax3.plot(theta, theta_projection, 'r-', linewidth=2)
+    ax3.fill_between(theta, theta_projection, alpha=0.3)
+    ax3.set_title('Projection along Theta', fontsize=12, pad=10)
+    ax3.set_xlabel('Theta (degree)', fontsize=10)
+    ax3.set_ylabel('Average X/X0', fontsize=10)
+    ax3.grid(True, linestyle='--', alpha=0.7)
+    ax3.spines['top'].set_visible(False)
+    ax3.spines['right'].set_visible(False)
+    
+    # Adjust main title position
+    plt.suptitle('Material Budget Analysis for {}'.format(detector), 
+                fontsize=16,  # Increase font size
+                y=0.95)  # Adjust title vertical position
+    
+    plt.tight_layout()
+    plt.savefig(output_file, dpi=600, bbox_inches='tight')
+    plt.show()
+
+def main():
+    parser = argparse.ArgumentParser(description='Process material scan data and generate plots.')
+    parser.add_argument('--input_file', help='Path to the input text file')
+    parser.add_argument('--output_file', default='material_budget.png', help='Path to the output png file')
+    parser.add_argument('--detector', default='all_world', help='Detector name')
+    parser.add_argument('--bins', type=int, default=100, help='theta bins is "bins"+1, phi bins is "bins"*2+1, (default: 100)')
+    parser.add_argument('--bias', type=int, default=0, help='Bias value for theta range (default: 0)')
+    args = parser.parse_args()
+
+    # all_world  Coil  Ecal  FTD  Hcal  Lumical  Muon  OTK  ParaffinEndcap  SIT  TPC  VXD
+    # bias detectors: BeamPipe onlyTracker
+    # Extract detector name from filename
+    x0_values = extract_x0_values(args.input_file)
+    plot_all_in_one(x0_values, args.output_file, args.bins, args.bias, args.detector)
+
+if __name__ == '__main__':
+    main() 
\ No newline at end of file